Belle II Software development
PXDHotPixelMaskCalibrationAlgorithm.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include <pxd/calibration/PXDHotPixelMaskCalibrationAlgorithm.h>
10#include <pxd/dbobjects/PXDMaskedPixelPar.h>
11#include <pxd/dbobjects/PXDDeadPixelPar.h>
12#include <pxd/dbobjects/PXDOccupancyInfoPar.h>
13#include <vxd/dataobjects/VxdID.h>
14
15#include <boost/format.hpp>
16#include <string>
17#include <vector>
18#include "TH1I.h"
19#include "TMath.h"
20
21using namespace std;
22using boost::format;
23using namespace Belle2;
24
25
29 deadPixelPayloadName("PXDDeadPixelPar"),
31{
33 " -------------------------- PXDHotPixelMak Calibration Algorithm ------------------------\n"
34 " \n"
35 " Algorithm which masks hot pixels with too large occupancy and dead pixels w/o no hits. \n"
36 " ----------------------------------------------------------------------------------------\n"
37 );
38}
39
41{
42 // Run list of current IoV and all runs to be calibrated
43 auto expRuns = getRunList();
44 auto expRunsAll = getRunListFromAllData();
45
46 // Save the current directory to change back later
47 TDirectory* currentDir = gDirectory;
48
49 // Create TFile if not exist
50 if (!m_file) {
51 std::string fileName = (this->getPrefix()) + "debug.root";
52 B2INFO("Creating file " << fileName);
53 m_file = std::make_shared<TFile>(fileName.c_str(), "RECREATE");
54 }
55
56 //m_file->cd();
57 string iov_str = str(format("E%1%_R%2%_E%3%_R%4%") % expRuns.front().first
58 % expRuns.front().second % expRuns.back().first % expRuns.back().second);
59 m_file->mkdir(iov_str.c_str());
60 m_file->cd(iov_str.c_str());
61
62 // Collector info
63 auto collector_pxdhits = getObjectPtr<TH1I>("PXDHits");
64 auto collector_pxdhitcounts = getObjectPtr<TH1I>("PXDHitCounts");
65
66 auto nevents = collector_pxdhits->GetEntries();
67
68 for (auto sensBin = 1; sensBin <= collector_pxdhitcounts->GetXaxis()->GetNbins(); sensBin++) {
69 string sensorDescr = collector_pxdhitcounts->GetXaxis()->GetBinLabel(sensBin);
70 VxdID id(sensorDescr);
71
72 // Get hitmap from collector
73 string name = str(format("PXD_%1%_PixelHitmap") % id.getID());
74 auto collector_pxdhitmap = getObjectPtr<TH1I>(name.c_str());
75 if (collector_pxdhitmap == nullptr) {
76 continue;
77 }
78 TH1I* debugHisto = new TH1I(*collector_pxdhitmap);
79 // Set overflow to total number of events
80 debugHisto->SetBinContent(debugHisto->GetNbinsX() + 1, nevents);
81 debugHisto->Write();
82 delete debugHisto;
83 }
84
85 // Close TFile
86 if (expRuns.back() == expRunsAll.back()) {
87 B2INFO("Reached Final ExpRun: (" << expRuns.back().first << ", " << expRuns.back().second << ")");
88 m_file->Close();
89
90 }
91 currentDir->cd();
92}
93
95{
96
97 auto collector_pxdhits = getObjectPtr<TH1I>("PXDHits");
98 auto collector_pxdhitcounts = getObjectPtr<TH1I>("PXDHitCounts");
99
100 // Check if there is any PXD hit
101 if (!collector_pxdhits) {
103 return c_OK;
104 else
105 return c_NotEnoughData;
106 }
107
108 // We should have some minimum number of events
109 auto nevents = collector_pxdhits->GetEntries();
110 if (nevents < minEvents) {
111 if (not forceContinueMasking) {
112 B2INFO("Not enough data: Only " << nevents << " events were collected!");
113 return c_NotEnoughData;
114 } else {
115 B2WARNING("Not enough data: Only " << nevents << " events were collected! The masking continuous but the mask may be empty.");
116 }
117 }
118
119 B2RESULT("Found total of " << nevents << " events in collected data.");
120
121 // Get the total number of PXD hits and sensors
122 unsigned long long int nPXDHits = 0;
123 int nPXDSensors = 0;
124 for (auto sensBin = 1; sensBin <= collector_pxdhitcounts->GetXaxis()->GetNbins(); sensBin++) {
125 // The bin label is assumed to be a string representation of VxdID
126 string sensorDescr = collector_pxdhitcounts->GetXaxis()->GetBinLabel(sensBin);
127 VxdID id(sensorDescr);
128 //Increment number of sensors
129 nPXDSensors += 1;
130 // Increment number of of collected hits
131 unsigned long long int nSensorHits = collector_pxdhitcounts->GetBinContent(sensBin);
132 nPXDHits += nSensorHits;
133
134 B2RESULT("Number of hits for sensor sensor=" << id << " is " << nSensorHits);
135 }
136
137 // We should have enough hits in the PXD before we decide a single sensor is dead
138 unsigned long long int minPXDHits = minHits * nPXDSensors * c_nUCells * c_nVCells;
139 if (nPXDHits < minPXDHits) {
140 if (not forceContinueMasking) {
141 B2INFO("Not enough data: Only " << nPXDHits << " raw hits were collected!");
142 return c_NotEnoughData;
143 } else {
144 B2WARNING("Not enough data: Only " << nPXDHits << " raw hits were collected! The masking continuous but the mask may be empty.");
145 }
146 }
147
148 B2RESULT("Found total of " << nPXDHits << " raw hits in collected data.");
149
150 // Check that the median number of hits is large enough
151 for (auto sensBin = 1; sensBin <= collector_pxdhitcounts->GetXaxis()->GetNbins(); sensBin++) {
152 // The bin label is assumed to be a string representation of VxdID
153 string sensorDescr = collector_pxdhitcounts->GetXaxis()->GetBinLabel(sensBin);
154 VxdID id(sensorDescr);
155
156 m_medianOfHitsMap[id] = 0.0;
157
158 // Get hitmap from collector
159 string name = str(format("PXD_%1%_PixelHitmap") % id.getID());
160 auto collector_pxdhitmap = getObjectPtr<TH1I>(name.c_str());
161
162 // Check if there was data collected for this sensor
163 if (collector_pxdhitmap == nullptr) continue;
164
165 // Compute the median number of hits to define a robust baseline for judging a channel fires too often
166 int nBins = collector_pxdhitmap->GetXaxis()->GetNbins();
167 double prob = 0.5;
168 vector<double> hitVec(nBins);
169
170 for (auto bin = 1; bin <= nBins; bin++) {
171 hitVec[bin - 1] = (double) collector_pxdhitmap->GetBinContent(bin);
172 }
173 double medianNumberOfHits;
174 TMath::Quantiles(nBins, 1, &hitVec[0], &medianNumberOfHits, &prob, kFALSE);
175 if (medianNumberOfHits <= 0) {
176 B2WARNING("Median number of hits per senor is smaller <1. Raise median to 1 instead.");
177 medianNumberOfHits = 1;
178 } else {
179 B2RESULT("Median of hits is " << medianNumberOfHits << " for sensor " << id);
180 }
181
182 // Keep the median of later use
183 m_medianOfHitsMap[id] = medianNumberOfHits;
184
185 // Check median number of hits is large enough
186 if (medianNumberOfHits < minHits) {
187 if (not forceContinueMasking) {
188 B2INFO("Not enough data: Median number of his is only " << medianNumberOfHits << "!");
189 return c_NotEnoughData;
190 } else {
191 B2WARNING("Not enough data: Median number of hits is only " << medianNumberOfHits <<
192 "! The masking continuous but the mask may be empty.");
193 }
194 }
195 }
196
197 // This is the occupancy info payload for conditions DB
198 PXDOccupancyInfoPar* occupancyInfoPar = new PXDOccupancyInfoPar();
199
200 // This is the dead pixels payload for conditions DB
201 PXDDeadPixelPar* deadPixelsPar = new PXDDeadPixelPar();
202
203 // This is the hot pixel masking payload for conditions DB
204 PXDMaskedPixelPar* maskedPixelsPar = new PXDMaskedPixelPar();
205
206 // Remember the number of events, helps to judge the reliability of calibrations.
207 occupancyInfoPar->setNumberOfEvents(nevents);
208
209 // Compute the masks for all sensors
210 for (auto sensBin = 1; sensBin <= collector_pxdhitcounts->GetXaxis()->GetNbins(); sensBin++) {
211 // The bin label is assumed to be a string representation of VxdID
212 string sensorDescr = collector_pxdhitcounts->GetXaxis()->GetBinLabel(sensBin);
213 VxdID id(sensorDescr);
214
215 // If we reach here, a sensor with no hits is deemed dead
216 int nSensorHits = collector_pxdhitcounts->GetBinContent(sensBin);
217 if (nSensorHits == 0) {
218 deadPixelsPar->maskSensor(id.getID());
219 continue;
220 }
221
222 // Get hitmap from collector
223 string name = str(format("PXD_%1%_PixelHitmap") % id.getID());
224 auto collector_pxdhitmap = getObjectPtr<TH1I>(name.c_str());
225 if (collector_pxdhitmap == nullptr) {
226 B2WARNING("Cannot find PixelHitmap although there should be hits. This is strange!");
227 continue;
228 }
229
230 double medianNumberOfHits = m_medianOfHitsMap[id];
231 int nBins = collector_pxdhitmap->GetXaxis()->GetNbins();
232
233 // Dead pixel masking
234 if (medianNumberOfHits >= minHits) {
235
236 // Bookkeeping for masking of drains and rows
237 vector<int> hitsAlongRow(c_nVCells, 0);
238 vector<int> hitsAlongDrain(c_nDrains, 0);
239 vector<int> nDeadAlongRow(c_nVCells, 0);
240
241 // Mask all single pixels exceeding medianNumberOfHits x multiplier
242 double inefficientPixelHitThr = inefficientPixelMultiplier * medianNumberOfHits;
243 B2RESULT("Pixel hit threshold for dead rows is " << inefficientPixelHitThr << " for sensor " << id);
244 // Accumulate hits along drains and rows
245 for (auto bin = 1; bin <= nBins; bin++) {
246 // Find the current pixel cell
247 int pixID = bin - 1;
248 int uCell = pixID / c_nVCells;
249 int vCell = pixID % c_nVCells;
250 int drainID = uCell * 4 + vCell % 4;
251 int nhits = collector_pxdhitmap->GetBinContent(bin);
252 hitsAlongDrain[drainID] += nhits;
253 hitsAlongRow[vCell] += nhits;
254 if (inefficientPixelMultiplier > 0 && nhits < inefficientPixelHitThr)
255 nDeadAlongRow[vCell] += 1;
256 }
257
258 // Dead row masking
259 B2INFO("Entering masking of dead rows ...");
260 for (auto vCell = 0; vCell < c_nVCells; vCell++) {
261 // Get number of hits per row
262 int nhits = hitsAlongRow[vCell];
263 int nDeadPixels = nDeadAlongRow[vCell];
264 // Mask dead row
265 if (nhits == 0 || nDeadPixels >= minInefficientPixels) {
266 deadPixelsPar->maskRow(id.getID(), vCell);
267 B2RESULT("Dead row with vCell=" << vCell << " on sensor " << id);
268 }
269
270 }
271
272 // Dead drain masking
273 B2INFO("Entering masking of dead drains ...");
274 for (auto drainID = 0; drainID < c_nDrains; drainID++) {
275 // Get number of hits per drain
276 int nhits = hitsAlongDrain[drainID];
277 // Mask dead drain
278 if (nhits == 0) {
279 deadPixelsPar->maskDrain(id.getID(), drainID);
280 B2RESULT("Dead drain line at drainID=" << drainID << " on sensor " << id);
281 }
282 }
283
284 // Dead pixel masking
285 B2INFO("Entering masking of single dead pixels ...");
286 for (auto bin = 1; bin <= nBins; bin++) {
287 // First, we mask single pixels exceeding hit threshold
288 int nhits = collector_pxdhitmap->GetBinContent(bin);
289 int pixID = bin - 1;
290 int uCell = pixID / c_nVCells;
291 int vCell = pixID % c_nVCells;
292 int drainID = uCell * 4 + vCell % 4;
293 // Mask dead pixel
294 if (nhits == 0 && !deadPixelsPar->isDeadRow(id.getID(), vCell) && !deadPixelsPar->isDeadDrain(id.getID(), drainID)) {
295 // This pixel is dead, we have to mask it
296 deadPixelsPar->maskSinglePixel(id.getID(), pixID);
297 B2RESULT("Dead single pixel with ucell=" << uCell << ", vcell=" << vCell << " on sensor " << id);
298 }
299 }
300 }
301
302 // Hot pixel masking
303
304 // Bookkeeping for masking hot drains
305 vector<float> unmaskedHitsAlongDrain(c_nDrains, 0);
306 vector<int> unmaskedCellsAlongDrain(c_nDrains, 0);
307
308 // Bookkeeping for maskign hot rows
309 vector<float> unmaskedHitsAlongRow(c_nVCells, 0);
310 vector<int> unmaskedCellsAlongRow(c_nVCells, 0);
311
312 // Mask all single pixels exceeding medianNumberOfHits x multiplier
313 double pixelHitThr = pixelMultiplier * medianNumberOfHits;
314 B2RESULT("Pixel hit threshold is " << pixelHitThr << " for sensor " << id);
315
316 // Mask all hot pixel for this sensor
317 for (auto bin = 1; bin <= nBins; bin++) {
318 // Find the current pixel cell
319 int pixID = bin - 1;
320 int uCell = pixID / c_nVCells;
321 int vCell = pixID % c_nVCells;
322 int drainID = uCell * 4 + vCell % 4;
323
324 // First, we mask single pixels exceeding hit threshold
325 float nhits = collector_pxdhitmap->GetBinContent(bin);
326 bool masked = false;
327
328 if (nhits > pixelHitThr) {
329 // This pixel is hot, we have to mask it
330 maskedPixelsPar->maskSinglePixel(id.getID(), pixID);
331 masked = true;
332 B2RESULT("Masking single pixel with ucell=" << uCell << ", vcell=" << vCell << " on sensor " << id);
333 }
334
335 // Then we accumulate hits along u and v direction for unmasked
336 // pixels
337 if (not masked) {
338 ++unmaskedCellsAlongDrain[drainID];
339 unmaskedHitsAlongDrain[drainID] += nhits;
340 ++unmaskedCellsAlongRow[vCell];
341 unmaskedHitsAlongRow[vCell] += nhits;
342 }
343 }
344
345 if (maskDrains) {
346 double drainHitThr = drainMultiplier * medianNumberOfHits;
347 B2RESULT("Drain hit threshold is " << drainHitThr << " for sensor " << id);
348
349 for (auto drainID = 0; drainID < c_nDrains; drainID++) {
350 if (unmaskedCellsAlongDrain[drainID] > 0) {
351 // Compute average number of hits per drain
352 float nhits = unmaskedHitsAlongDrain[drainID] / unmaskedCellsAlongDrain[drainID];
353 // Mask residual hot drain
354 if (nhits > drainHitThr) {
355 for (auto iGate = 0; iGate < 192; iGate++) {
356 int uCell = drainID / 4;
357 int vCell = drainID % 4 + iGate * 4;
358 maskedPixelsPar->maskSinglePixel(id.getID(), uCell * c_nVCells + vCell);
359 }
360 B2RESULT("Masking drain line with drainID=" << drainID << " on sensor " << id);
361 }
362 }
363 }
364 }
365
366 if (maskRows) {
367 double rowHitThr = rowMultiplier * medianNumberOfHits;
368 B2RESULT("Row hit threshold is " << rowHitThr << " for sensor " << id);
369
370 for (auto vCell = 0; vCell < c_nVCells; vCell++) {
371 if (unmaskedCellsAlongRow[vCell] > 0) {
372 // Compute average number of hits per row
373 float nhits = unmaskedHitsAlongRow[vCell] / unmaskedCellsAlongRow[vCell];
374 // Mask residual hot row
375 if (nhits > rowHitThr) {
376 for (auto uCell = 0; uCell < c_nUCells; uCell++)
377 maskedPixelsPar->maskSinglePixel(id.getID(), uCell * c_nVCells + vCell);
378
379 B2RESULT("Masking complete row with vCell=" << vCell << " on sensor " << id);
380 }
381 }
382 }
383 }
384
385 // After the masking is done, we compute the average sensor occupancy after
386 // hot pixel masking.
387
388 // Count all unmasked hits
389 int numberOfUnmaskedHits = 0;
390
391 for (auto bin = 1; bin <= nBins; bin++) {
392 // Find the current pixel cell
393 int pixID = bin - 1;
394
395 if (maskedPixelsPar->pixelOK(id.getID(), pixID)) {
396 numberOfUnmaskedHits += collector_pxdhitmap->GetBinContent(bin);
397 }
398 }
399
400 // Compute mean occupancy before masking
401 float meanOccupancyAfterMasking = (float)numberOfUnmaskedHits / nevents / c_nVCells / c_nUCells;
402 B2RESULT("Hotpixel filtered occupancy sensor=" << id << " is " << meanOccupancyAfterMasking);
403
404 occupancyInfoPar->setOccupancy(id.getID(), meanOccupancyAfterMasking);
405 }
406
407 for (auto elem : maskedPixelsPar->getMaskedPixelMap()) {
408 auto id = elem.first;
409 auto singles = elem.second;
410 B2RESULT("SensorID " << VxdID(id) << " has filtered occupancy of " << occupancyInfoPar->getOccupancy(id));
411 B2RESULT("SensorID " << VxdID(id) << " has fraction of masked pixels of " << (float)singles.size() / (c_nVCells * c_nUCells));
412 }
413
414 // Save the hot pixel mask to database. Note that this will set the database object name to the same as the collector but you
415 // are free to change it.
416 saveCalibration(maskedPixelsPar, "PXDMaskedPixelPar");
418 saveCalibration(occupancyInfoPar, "PXDOccupancyInfoPar");
419
420 // Create debugging histo if we asked for it
422
423 B2INFO("PXDHotPixelMask Calibration Successful");
424 return c_OK;
425}
void saveCalibration(TClonesArray *data, const std::string &name)
Store DBArray payload with given name with default IOV.
std::vector< Calibration::ExpRun > getRunListFromAllData() const
Get the complete list of runs from inspection of collected data.
void setDescription(const std::string &description)
Set algorithm description (in constructor)
const std::vector< Calibration::ExpRun > & getRunList() const
Get the list of runs for which calibration is called.
std::string getPrefix() const
Get the prefix used for getting calibration data.
EResult
The result of calibration.
@ c_OK
Finished successfully =0 in Python.
@ c_NotEnoughData
Needs more data =2 in Python.
CalibrationAlgorithm(const std::string &collectorModuleName)
Constructor - sets the prefix for collected objects (won't be accesses until execute(....
The payload telling which PXD pixel is dead (=Readout system does not receive signals)
void maskSensor(unsigned short sensorID)
Mask sensor.
void maskRow(unsigned short sensorID, unsigned int vCellID)
Mask single row.
void maskSinglePixel(unsigned short sensorID, unsigned int pixID)
Mask single pixel.
void maskDrain(unsigned short sensorID, unsigned int drainID)
Mask single drain.
bool isDeadDrain(unsigned short sensorID, unsigned int drainID) const
Check whether a drain is dead.
bool isDeadRow(unsigned short sensorID, unsigned int vCellID) const
Check whether a row is dead.
float inefficientPixelMultiplier
The occupancy threshold for inefficient (or dead) pixels is the median occupancy x inefficientPixelMu...
bool maskRows
Mask hot rows with too high average occupancy after single pixel masking.
std::map< VxdID, double > m_medianOfHitsMap
map of VxdID to median hits of each sensor
float drainMultiplier
The occupancy threshold for masking hot drains is the median occupancy x drainMultiplier.
int minInefficientPixels
The minimum number of inefficient (or dead) pixels per row to define an inefficient row.
const unsigned short c_nDrains
Number of drain lines of Belle II PXD sensors.
bool maskDrains
Mask hot drain lines with too high average occupancy after single pixel masking.
const unsigned short c_nVCells
Number of vCells of Belle II PXD sensors.
int m_debugHisto
Set if a debugging histogram should be created in the algorithm output directory.
float pixelMultiplier
The occupancy threshold for masking hot single pixels is the median occupancy x pixelMultiplier.
const unsigned short c_nUCells
Number of uCells of Belle II PXD sensors.
PXDHotPixelMaskCalibrationAlgorithm()
Constructor set the prefix to PXDHotPixelMaskCalibrationAlgorithm.
float rowMultiplier
The occupancy threshold for masking hot rows is the median occupancy x rowMultiplier.
virtual EResult calibrate() override
Run algo on data.
bool forceContinueMasking
Force continue masking in almost empty runs instead of returning c_NotEnoughData.
int minHits
Minimum median number of hits per pixel needed for dead pixel masking.
std::string deadPixelPayloadName
Payload name for PXDDeadPixelPar used for more defective pixels from damaged gates.
void createDebugHistogram()
Perform debug histogram file creation.
The payload telling which PXD pixel to mask (ignore)
bool pixelOK(unsigned short sensorID, unsigned int pixID) const
Check whether a pixel on a given sensor is OK or not.
const std::unordered_map< unsigned short, MaskedSinglePixelsSet > & getMaskedPixelMap() const
Return unordered_map with all masked single pixels in PXD.
void maskSinglePixel(unsigned short sensorID, unsigned int pixID)
Mask single pixel.
The payload collecting some meta information from running the masking algorithm.
void setOccupancy(unsigned short sensorID, float occupancy)
Set occupancy.
void setNumberOfEvents(int nEvents)
Set number of events used for occupancy estimation.
float getOccupancy(unsigned short sensorID) const
Get occupancy.
Class to uniquely identify a any structure of the PXD and SVD.
Definition VxdID.h:32
std::shared_ptr< T > getObjectPtr(const std::string &name, const std::vector< Calibration::ExpRun > &requestedRuns)
Get calibration data object by name and list of runs, the Merge function will be called to generate t...
int getID(const std::vector< double > &breaks, double t)
get id of the time point t
Definition calibTools.h:60
Abstract base class for different kinds of events.
STL namespace.