Belle II Software  release-05-01-25
PXDHotPixelMaskCalibrationAlgorithm.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2013 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Benjamin Schwenker *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <vxd/dataobjects/VxdID.h>
12 #include <pxd/calibration/PXDHotPixelMaskCalibrationAlgorithm.h>
13 #include <pxd/dbobjects/PXDMaskedPixelPar.h>
14 #include <pxd/dbobjects/PXDDeadPixelPar.h>
15 #include <pxd/dbobjects/PXDOccupancyInfoPar.h>
16 
17 #include <boost/format.hpp>
18 #include <string>
19 #include <vector>
20 #include <map>
21 #include "TH1I.h"
22 #include "TMath.h"
23 
24 using namespace std;
25 using boost::format;
26 using namespace Belle2;
27 
28 
29 PXDHotPixelMaskCalibrationAlgorithm::PXDHotPixelMaskCalibrationAlgorithm(): CalibrationAlgorithm("PXDHotPixelMaskCollector"),
30  forceContinueMasking(false), minEvents(10000), minHits(20), pixelMultiplier(7), maskDrains(true),
31  drainMultiplier(7), maskRows(true), rowMultiplier(7)
32 {
34  " -------------------------- PXDHotPixelMak Calibration Algorithm ------------------------\n"
35  " \n"
36  " Algorithm which masks hot pixels with too large occupancy and dead pixels w/o no hits. \n"
37  " ----------------------------------------------------------------------------------------\n"
38  );
39 }
40 
42 {
43 
44  auto collector_pxdhits = getObjectPtr<TH1I>("PXDHits");
45  auto collector_pxdhitcounts = getObjectPtr<TH1I>("PXDHitCounts");
46 
47  // We should have some minimum number of events
48  auto nevents = collector_pxdhits->GetEntries();
49  if (nevents < minEvents) {
50  if (not forceContinueMasking) {
51  B2INFO("Not enough data: Only " << nevents << " events were collected!");
52  return c_NotEnoughData;
53  } else {
54  B2WARNING("Not enough data: Only " << nevents << " events were collected! The masking continous but the mask may be empty.");
55  }
56  }
57 
58  B2RESULT("Found total of " << nevents << " events in collected data.");
59 
60  // Get the total number of PXD hits and sensors
61  unsigned long long int nPXDHits = 0;
62  int nPXDSensors = 0;
63  for (auto sensBin = 1; sensBin <= collector_pxdhitcounts->GetXaxis()->GetNbins(); sensBin++) {
64  // The bin label is assumed to be a string representation of VxdID
65  string sensorDescr = collector_pxdhitcounts->GetXaxis()->GetBinLabel(sensBin);
66  VxdID id(sensorDescr);
67  //Increment number of sensors
68  nPXDSensors += 1;
69  // Increment number of of collected hits
70  unsigned long long int nSensorHits = collector_pxdhitcounts->GetBinContent(sensBin);
71  nPXDHits += nSensorHits;
72 
73  B2RESULT("Number of hits for sensor sensor=" << id << " is " << nSensorHits);
74  }
75 
76  // We should have enough hits in the PXD before we decide a single sensor is dead
77  unsigned long long int minPXDHits = minHits * nPXDSensors * c_nUCells * c_nVCells;
78  if (nPXDHits < minPXDHits) {
79  if (not forceContinueMasking) {
80  B2INFO("Not enough data: Only " << nPXDHits << " raw hits were collected!");
81  return c_NotEnoughData;
82  } else {
83  B2WARNING("Not enough data: Only " << nPXDHits << " raw hits were collected! The masking continous but the mask may be empty.");
84  }
85  }
86 
87  B2RESULT("Found total of " << nPXDHits << " raw hits in collected data.");
88 
89  // Check that the median number of hits is large enough
90  map<VxdID, double> medianOfHitsMap;
91  for (auto sensBin = 1; sensBin <= collector_pxdhitcounts->GetXaxis()->GetNbins(); sensBin++) {
92  // The bin label is assumed to be a string representation of VxdID
93  string sensorDescr = collector_pxdhitcounts->GetXaxis()->GetBinLabel(sensBin);
94  VxdID id(sensorDescr);
95 
96  medianOfHitsMap[id] = 0.0;
97 
98  // Get hitmap from collector
99  string name = str(format("PXD_%1%_PixelHitmap") % id.getID());
100  auto collector_pxdhitmap = getObjectPtr<TH1I>(name.c_str());
101 
102  // Check if there was data collected for this sensor
103  if (collector_pxdhitmap == nullptr) continue;
104 
105  // Compute the median number of hits to define a robust baseline for judging a channel fires too often
106  int nBins = collector_pxdhitmap->GetXaxis()->GetNbins();
107  double prob = 0.5;
108  vector<double> hitVec(nBins);
109 
110  for (auto bin = 1; bin <= nBins; bin++) {
111  hitVec[bin - 1] = (double) collector_pxdhitmap->GetBinContent(bin);
112  }
113  double medianNumberOfHits;
114  TMath::Quantiles(nBins, 1, &hitVec[0], &medianNumberOfHits, &prob, kFALSE);
115  if (medianNumberOfHits <= 0) {
116  B2WARNING("Median number of hits per senor is smaller <1. Raise median to 1 instead.");
117  medianNumberOfHits = 1;
118  } else {
119  B2RESULT("Median of hits is " << medianNumberOfHits << " for sensor " << id);
120  }
121 
122  // Keep the median of later use
123  medianOfHitsMap[id] = medianNumberOfHits;
124 
125  // Check median number of hits is large enough
126  if (medianNumberOfHits < minHits) {
127  if (not forceContinueMasking) {
128  B2INFO("Not enough data: Median number of his is only " << medianNumberOfHits << "!");
129  return c_NotEnoughData;
130  } else {
131  B2WARNING("Not enough data: Median number of hits is only " << medianNumberOfHits <<
132  "! The masking continous but the mask may be empty.");
133  }
134  }
135  }
136 
137  // This is the occupancy info payload for conditions DB
138  PXDOccupancyInfoPar* occupancyInfoPar = new PXDOccupancyInfoPar();
139 
140  // This is the dead pixels payload for conditions DB
141  PXDDeadPixelPar* deadPixelsPar = new PXDDeadPixelPar();
142 
143  // This is the hot pixel masking payload for conditions DB
144  PXDMaskedPixelPar* maskedPixelsPar = new PXDMaskedPixelPar();
145 
146  // Remember the number of events, helps to judge the reliability of calibrations.
147  occupancyInfoPar->setNumberOfEvents(nevents);
148 
149  // Compute the masks for all sensors
150  for (auto sensBin = 1; sensBin <= collector_pxdhitcounts->GetXaxis()->GetNbins(); sensBin++) {
151  // The bin label is assumed to be a string representation of VxdID
152  string sensorDescr = collector_pxdhitcounts->GetXaxis()->GetBinLabel(sensBin);
153  VxdID id(sensorDescr);
154 
155  // If we reach here, a sensor with no hits is deemed dead
156  int nSensorHits = collector_pxdhitcounts->GetBinContent(sensBin);
157  if (nSensorHits == 0) {
158  deadPixelsPar->maskSensor(id.getID());
159  continue;
160  }
161 
162  // Get hitmap from collector
163  string name = str(format("PXD_%1%_PixelHitmap") % id.getID());
164  auto collector_pxdhitmap = getObjectPtr<TH1I>(name.c_str());
165  if (collector_pxdhitmap == nullptr) {
166  B2WARNING("Cannot find PixelHitmap although there should be hits. This is strange!");
167  continue;
168  }
169 
170  double medianNumberOfHits = medianOfHitsMap[id];
171  int nBins = collector_pxdhitmap->GetXaxis()->GetNbins();
172 
173  // Dead pixel masking
174  if (medianNumberOfHits >= minHits) {
175 
176  // Bookkeeping for masking of drains and rows
177  vector<int> hitsAlongRow(c_nVCells, 0);
178  vector<int> hitsAlongDrain(c_nDrains, 0);
179 
180  // Accumulate hits along drains and rows
181  for (auto bin = 1; bin <= nBins; bin++) {
182  // Find the current pixel cell
183  int pixID = bin - 1;
184  int uCell = pixID / c_nVCells;
185  int vCell = pixID % c_nVCells;
186  int drainID = uCell * 4 + vCell % 4;
187  int nhits = collector_pxdhitmap->GetBinContent(bin);
188  hitsAlongDrain[drainID] += nhits;
189  hitsAlongRow[vCell] += nhits;
190  }
191 
192  // Dead row masking
193  B2INFO("Entering masking of dead rows ...");
194  for (auto vCell = 0; vCell < c_nVCells; vCell++) {
195  // Get number of hits per row
196  int nhits = hitsAlongRow[vCell];
197  // Mask dead row
198  if (nhits == 0) {
199  deadPixelsPar->maskRow(id.getID(), vCell);
200  B2RESULT("Dead row with vCell=" << vCell << " on sensor " << id);
201  }
202  }
203 
204  // Dead drain masking
205  B2INFO("Entering masking of dead drains ...");
206  for (auto drainID = 0; drainID < c_nDrains; drainID++) {
207  // Get number of hits per drain
208  int nhits = hitsAlongDrain[drainID];
209  // Mask dead drain
210  if (nhits == 0) {
211  deadPixelsPar->maskDrain(id.getID(), drainID);
212  B2RESULT("Dead drain line at drainID=" << drainID << " on sensor " << id);
213  }
214  }
215 
216  // Dead pixel masking
217  B2INFO("Entering masking of single dead pixels ...");
218  for (auto bin = 1; bin <= nBins; bin++) {
219  // First, we mask single pixels exceeding hit threshold
220  int nhits = collector_pxdhitmap->GetBinContent(bin);
221  int pixID = bin - 1;
222  int uCell = pixID / c_nVCells;
223  int vCell = pixID % c_nVCells;
224  int drainID = uCell * 4 + vCell % 4;
225  // Mask dead pixel
226  if (nhits == 0 && !deadPixelsPar->isDeadRow(id.getID(), vCell) && !deadPixelsPar->isDeadDrain(id.getID(), drainID)) {
227  // This pixel is dead, we have to mask it
228  deadPixelsPar->maskSinglePixel(id.getID(), pixID);
229  B2RESULT("Dead single pixel with ucell=" << uCell << ", vcell=" << vCell << " on sensor " << id);
230  }
231  }
232  }
233 
234  // Hot pixel masking
235 
236  // Bookkeeping for masking hot drains
237  vector<float> unmaskedHitsAlongDrain(c_nDrains, 0);
238  vector<int> unmaskedCellsAlongDrain(c_nDrains, 0);
239 
240  // Bookkeeping for maskign hot rows
241  vector<float> unmaskedHitsAlongRow(c_nVCells, 0);
242  vector<int> unmaskedCellsAlongRow(c_nVCells, 0);
243 
244  // Mask all single pixels exceeding medianNumberOfHits x multiplier
245  double pixelHitThr = pixelMultiplier * medianNumberOfHits;
246  B2RESULT("Pixel hit threshold is " << pixelHitThr << " for sensor " << id);
247 
248  // Mask all hot pixel for this sensor
249  for (auto bin = 1; bin <= nBins; bin++) {
250  // Find the current pixel cell
251  int pixID = bin - 1;
252  int uCell = pixID / c_nVCells;
253  int vCell = pixID % c_nVCells;
254  int drainID = uCell * 4 + vCell % 4;
255 
256  // First, we mask single pixels exceeding hit threshold
257  float nhits = collector_pxdhitmap->GetBinContent(bin);
258  bool masked = false;
259 
260  if (nhits > pixelHitThr) {
261  // This pixel is hot, we have to mask it
262  maskedPixelsPar->maskSinglePixel(id.getID(), pixID);
263  masked = true;
264  B2RESULT("Masking single pixel with ucell=" << uCell << ", vcell=" << vCell << " on sensor " << id);
265  }
266 
267  // Then we accumulate hits along u and v direction for unmasked
268  // pixels
269  if (not masked) {
270  ++unmaskedCellsAlongDrain[drainID];
271  unmaskedHitsAlongDrain[drainID] += nhits;
272  ++unmaskedCellsAlongRow[vCell];
273  unmaskedHitsAlongRow[vCell] += nhits;
274  }
275  }
276 
277  if (maskDrains) {
278  double drainHitThr = drainMultiplier * medianNumberOfHits;
279  B2RESULT("Drain hit threshold is " << drainHitThr << " for sensor " << id);
280 
281  for (auto drainID = 0; drainID < c_nDrains; drainID++) {
282  if (unmaskedCellsAlongDrain[drainID] > 0) {
283  // Compute average number of hits per drain
284  float nhits = unmaskedHitsAlongDrain[drainID] / unmaskedCellsAlongDrain[drainID];
285  // Mask residual hot drain
286  if (nhits > drainHitThr) {
287  for (auto iGate = 0; iGate < 192; iGate++) {
288  int uCell = drainID / 4;
289  int vCell = drainID % 4 + iGate * 4;
290  maskedPixelsPar->maskSinglePixel(id.getID(), uCell * c_nVCells + vCell);
291  }
292  B2RESULT("Masking drain line with drainID=" << drainID << " on sensor " << id);
293  }
294  }
295  }
296  }
297 
298  if (maskRows) {
299  double rowHitThr = rowMultiplier * medianNumberOfHits;
300  B2RESULT("Row hit threshold is " << rowHitThr << " for sensor " << id);
301 
302  for (auto vCell = 0; vCell < c_nVCells; vCell++) {
303  if (unmaskedCellsAlongRow[vCell] > 0) {
304  // Compute average number of hits per row
305  float nhits = unmaskedHitsAlongRow[vCell] / unmaskedCellsAlongRow[vCell];
306  // Mask residual hot row
307  if (nhits > rowHitThr) {
308  for (auto uCell = 0; uCell < c_nUCells; uCell++)
309  maskedPixelsPar->maskSinglePixel(id.getID(), uCell * c_nVCells + vCell);
310 
311  B2RESULT("Masking complete row with vCell=" << vCell << " on sensor " << id);
312  }
313  }
314  }
315  }
316 
317  // After the masking is done, we compute the average sensor occupancy after
318  // hot pixel masking.
319 
320  // Count all unmasked hits
321  int numberOfUnmaskedHits = 0;
322 
323  for (auto bin = 1; bin <= nBins; bin++) {
324  // Find the current pixel cell
325  int pixID = bin - 1;
326 
327  if (maskedPixelsPar->pixelOK(id.getID(), pixID)) {
328  numberOfUnmaskedHits += collector_pxdhitmap->GetBinContent(bin);
329  }
330  }
331 
332  // Compute mean occupancy before masking
333  float meanOccupancyAfterMasking = (float)numberOfUnmaskedHits / nevents / c_nVCells / c_nUCells;
334  B2RESULT("Hotpixel filtered occupancy sensor=" << id << " is " << meanOccupancyAfterMasking);
335 
336  occupancyInfoPar->setOccupancy(id.getID(), meanOccupancyAfterMasking);
337  }
338 
339  for (auto elem : maskedPixelsPar->getMaskedPixelMap()) {
340  auto id = elem.first;
341  auto singles = elem.second;
342  B2RESULT("SensorID " << VxdID(id) << " has filtered occupancy of " << occupancyInfoPar->getOccupancy(id));
343  B2RESULT("SensorID " << VxdID(id) << " has fraction of masked pixels of " << (float)singles.size() / (c_nVCells * c_nUCells));
344  }
345 
346  // Save the hot pixel mask to database. Note that this will set the database object name to the same as the collector but you
347  // are free to change it.
348  saveCalibration(maskedPixelsPar, "PXDMaskedPixelPar");
349  saveCalibration(deadPixelsPar, "PXDDeadPixelPar");
350  saveCalibration(occupancyInfoPar, "PXDOccupancyInfoPar");
351 
352  B2INFO("PXDHotPixelMask Calibration Successful");
353  return c_OK;
354 }
Belle2::PXDHotPixelMaskCalibrationAlgorithm::forceContinueMasking
bool forceContinueMasking
Force continue masking in almost empty runs instead of returning c_NotEnoughData.
Definition: PXDHotPixelMaskCalibrationAlgorithm.h:37
Belle2::PXDDeadPixelPar::maskDrain
void maskDrain(unsigned short sensorID, unsigned int drainID)
Mask single drain.
Definition: PXDDeadPixelPar.h:69
Belle2::PXDHotPixelMaskCalibrationAlgorithm::minHits
int minHits
Minimum median number of hits per pixel needed for dead pixel masking.
Definition: PXDHotPixelMaskCalibrationAlgorithm.h:43
Belle2::VxdID
Class to uniquely identify a any structure of the PXD and SVD.
Definition: VxdID.h:43
Belle2::PXDHotPixelMaskCalibrationAlgorithm::drainMultiplier
float drainMultiplier
The occupancy threshold for masking hot drains is the median occupancy x drainMultiplier.
Definition: PXDHotPixelMaskCalibrationAlgorithm.h:52
Belle2::PXDDeadPixelPar::isDeadRow
bool isDeadRow(unsigned short sensorID, unsigned int vCellID) const
Check whether a row is dead.
Definition: PXDDeadPixelPar.h:150
Belle2::PXDHotPixelMaskCalibrationAlgorithm::pixelMultiplier
float pixelMultiplier
The occupancy threshold for masking hot single pixels is the median occupancy x pixelMultiplier.
Definition: PXDHotPixelMaskCalibrationAlgorithm.h:46
Belle2::CalibrationAlgorithm::saveCalibration
void saveCalibration(TClonesArray *data, const std::string &name)
Store DBArray payload with given name with default IOV.
Definition: CalibrationAlgorithm.cc:290
Belle2::getID
int getID(const std::vector< double > &breaks, double t)
get id of the time point t
Definition: calibTools.h:71
Belle2::PXDHotPixelMaskCalibrationAlgorithm::c_nDrains
const unsigned short c_nDrains
Number of drain lines of Belle II PXD sensors.
Definition: PXDHotPixelMaskCalibrationAlgorithm.h:71
Belle2::PXDHotPixelMaskCalibrationAlgorithm::calibrate
virtual EResult calibrate() override
Run algo on data.
Definition: PXDHotPixelMaskCalibrationAlgorithm.cc:41
Belle2::PXDHotPixelMaskCalibrationAlgorithm::minEvents
int minEvents
Minimum number of collected events.
Definition: PXDHotPixelMaskCalibrationAlgorithm.h:40
Belle2::PXDDeadPixelPar::maskSensor
void maskSensor(unsigned short sensorID)
Mask sensor.
Definition: PXDDeadPixelPar.h:55
Belle2::CalibrationAlgorithm::c_OK
@ c_OK
Finished successfuly =0 in Python.
Definition: CalibrationAlgorithm.h:51
Belle2::CalibrationAlgorithm::setDescription
void setDescription(const std::string &description)
Set algorithm description (in constructor)
Definition: CalibrationAlgorithm.h:331
Belle2::PXDOccupancyInfoPar::setOccupancy
void setOccupancy(unsigned short sensorID, float occupancy)
Set occupancy.
Definition: PXDOccupancyInfoPar.h:57
Belle2::PXDDeadPixelPar
The payload telling which PXD pixel is dead (=Readout system does not receive signals)
Definition: PXDDeadPixelPar.h:43
Belle2::PXDHotPixelMaskCalibrationAlgorithm::maskDrains
bool maskDrains
Mask hot drain lines with too high average occupancy after single pixel masking.
Definition: PXDHotPixelMaskCalibrationAlgorithm.h:49
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::PXDOccupancyInfoPar::setNumberOfEvents
void setNumberOfEvents(int nEvents)
Set number of events used for occupancy estimation.
Definition: PXDOccupancyInfoPar.h:46
Belle2::PXDHotPixelMaskCalibrationAlgorithm::rowMultiplier
float rowMultiplier
The occupancy threshold for masking hot rows is the median occupancy x rowMultiplier.
Definition: PXDHotPixelMaskCalibrationAlgorithm.h:58
Belle2::PXDDeadPixelPar::maskSinglePixel
void maskSinglePixel(unsigned short sensorID, unsigned int pixID)
Mask single pixel.
Definition: PXDDeadPixelPar.h:115
Belle2::PXDMaskedPixelPar::getMaskedPixelMap
const std::unordered_map< unsigned short, MaskedSinglePixelsSet > & getMaskedPixelMap() const
Return unordered_map with all masked single pixels in PXD.
Definition: PXDMaskedPixelPar.h:87
Belle2::CalibrationAlgorithm::EResult
EResult
The result of calibration.
Definition: CalibrationAlgorithm.h:50
Belle2::CalibrationAlgorithm::c_NotEnoughData
@ c_NotEnoughData
Needs more data =2 in Python.
Definition: CalibrationAlgorithm.h:53
Belle2::PXDOccupancyInfoPar::getOccupancy
float getOccupancy(unsigned short sensorID) const
Get occupancy.
Definition: PXDOccupancyInfoPar.h:67
Belle2::PXDMaskedPixelPar::maskSinglePixel
void maskSinglePixel(unsigned short sensorID, unsigned int pixID)
Mask single pixel.
Definition: PXDMaskedPixelPar.h:49
Belle2::CalibrationAlgorithm
Base class for calibration algorithms.
Definition: CalibrationAlgorithm.h:47
Belle2::PXDHotPixelMaskCalibrationAlgorithm::c_nUCells
const unsigned short c_nUCells
Number of uCells of Belle II PXD sensors.
Definition: PXDHotPixelMaskCalibrationAlgorithm.h:69
Belle2::PXDDeadPixelPar::isDeadDrain
bool isDeadDrain(unsigned short sensorID, unsigned int drainID) const
Check whether a drain is dead.
Definition: PXDDeadPixelPar.h:168
Belle2::PXDOccupancyInfoPar
The payload collecting some meta information from running the masking algorithm.
Definition: PXDOccupancyInfoPar.h:38
Belle2::PXDMaskedPixelPar::pixelOK
bool pixelOK(unsigned short sensorID, unsigned int pixID) const
Check whether a pixel on a given sensor is OK or not.
Definition: PXDMaskedPixelPar.h:72
Belle2::PXDDeadPixelPar::maskRow
void maskRow(unsigned short sensorID, unsigned int vCellID)
Mask single row.
Definition: PXDDeadPixelPar.h:92
Belle2::PXDHotPixelMaskCalibrationAlgorithm::c_nVCells
const unsigned short c_nVCells
Number of vCells of Belle II PXD sensors.
Definition: PXDHotPixelMaskCalibrationAlgorithm.h:67
Belle2::PXDHotPixelMaskCalibrationAlgorithm::maskRows
bool maskRows
Mask hot rows with too high average occupancy after single pixel masking.
Definition: PXDHotPixelMaskCalibrationAlgorithm.h:55
Belle2::PXDMaskedPixelPar
The payload telling which PXD pixel to mask (ignore)
Definition: PXDMaskedPixelPar.h:34