11 #include <pxd/calibration/PXDHotPixelMaskCalibrationAlgorithm.h>
12 #include <pxd/dbobjects/PXDMaskedPixelPar.h>
13 #include <pxd/dbobjects/PXDDeadPixelPar.h>
14 #include <pxd/dbobjects/PXDOccupancyInfoPar.h>
16 #include <boost/format.hpp>
27 PXDHotPixelMaskCalibrationAlgorithm::PXDHotPixelMaskCalibrationAlgorithm():
CalibrationAlgorithm(
"PXDHotPixelMaskCollector"),
28 forceContinueMasking(false), minEvents(10000), minHits(20), pixelMultiplier(7), maskDrains(true),
29 drainMultiplier(7), maskRows(true), rowMultiplier(7)
32 " -------------------------- PXDHotPixelMak Calibration Algorithm ------------------------\n"
34 " Algorithm which masks hot pixels with too large occupancy and dead pixels w/o no hits. \n"
35 " ----------------------------------------------------------------------------------------\n"
46 TDirectory* currentDir = gDirectory;
50 std::string fileName = (this->
getPrefix()) +
"debug.root";
51 B2INFO(
"Creating file " << fileName);
52 m_file = std::make_shared<TFile>(fileName.c_str(),
"RECREATE");
56 string iov_str = str(format(
"E%1%_R%2%_E%3%_R%4%") % expRuns.front().first
57 % expRuns.front().second % expRuns.back().first % expRuns.back().second);
58 m_file->mkdir(iov_str.c_str());
59 m_file->cd(iov_str.c_str());
62 auto collector_pxdhits = getObjectPtr<TH1I>(
"PXDHits");
63 auto collector_pxdhitcounts = getObjectPtr<TH1I>(
"PXDHitCounts");
65 auto nevents = collector_pxdhits->GetEntries();
67 for (
auto sensBin = 1; sensBin <= collector_pxdhitcounts->GetXaxis()->GetNbins(); sensBin++) {
68 string sensorDescr = collector_pxdhitcounts->GetXaxis()->GetBinLabel(sensBin);
69 VxdID id(sensorDescr);
72 string name = str(format(
"PXD_%1%_PixelHitmap") %
id.
getID());
73 auto collector_pxdhitmap = getObjectPtr<TH1I>(name.c_str());
74 if (collector_pxdhitmap ==
nullptr) {
77 TH1I* debugHisto =
new TH1I(*collector_pxdhitmap);
79 debugHisto->SetBinContent(debugHisto->GetNbinsX() + 1, nevents);
85 if (expRuns.back() == expRunsAll.back()) {
86 B2INFO(
"Reached Final ExpRun: (" << expRuns.back().first <<
", " << expRuns.back().second <<
")");
96 auto collector_pxdhits = getObjectPtr<TH1I>(
"PXDHits");
97 auto collector_pxdhitcounts = getObjectPtr<TH1I>(
"PXDHitCounts");
100 auto nevents = collector_pxdhits->GetEntries();
103 B2INFO(
"Not enough data: Only " << nevents <<
" events were collected!");
106 B2WARNING(
"Not enough data: Only " << nevents <<
" events were collected! The masking continous but the mask may be empty.");
110 B2RESULT(
"Found total of " << nevents <<
" events in collected data.");
113 unsigned long long int nPXDHits = 0;
115 for (
auto sensBin = 1; sensBin <= collector_pxdhitcounts->GetXaxis()->GetNbins(); sensBin++) {
117 string sensorDescr = collector_pxdhitcounts->GetXaxis()->GetBinLabel(sensBin);
118 VxdID id(sensorDescr);
122 unsigned long long int nSensorHits = collector_pxdhitcounts->GetBinContent(sensBin);
123 nPXDHits += nSensorHits;
125 B2RESULT(
"Number of hits for sensor sensor=" <<
id <<
" is " << nSensorHits);
130 if (nPXDHits < minPXDHits) {
132 B2INFO(
"Not enough data: Only " << nPXDHits <<
" raw hits were collected!");
135 B2WARNING(
"Not enough data: Only " << nPXDHits <<
" raw hits were collected! The masking continous but the mask may be empty.");
139 B2RESULT(
"Found total of " << nPXDHits <<
" raw hits in collected data.");
142 for (
auto sensBin = 1; sensBin <= collector_pxdhitcounts->GetXaxis()->GetNbins(); sensBin++) {
144 string sensorDescr = collector_pxdhitcounts->GetXaxis()->GetBinLabel(sensBin);
145 VxdID id(sensorDescr);
150 string name = str(format(
"PXD_%1%_PixelHitmap") %
id.
getID());
151 auto collector_pxdhitmap = getObjectPtr<TH1I>(name.c_str());
154 if (collector_pxdhitmap ==
nullptr)
continue;
157 int nBins = collector_pxdhitmap->GetXaxis()->GetNbins();
159 vector<double> hitVec(nBins);
161 for (
auto bin = 1; bin <= nBins; bin++) {
162 hitVec[bin - 1] = (double) collector_pxdhitmap->GetBinContent(bin);
164 double medianNumberOfHits;
165 TMath::Quantiles(nBins, 1, &hitVec[0], &medianNumberOfHits, &prob, kFALSE);
166 if (medianNumberOfHits <= 0) {
167 B2WARNING(
"Median number of hits per senor is smaller <1. Raise median to 1 instead.");
168 medianNumberOfHits = 1;
170 B2RESULT(
"Median of hits is " << medianNumberOfHits <<
" for sensor " <<
id);
177 if (medianNumberOfHits <
minHits) {
179 B2INFO(
"Not enough data: Median number of his is only " << medianNumberOfHits <<
"!");
182 B2WARNING(
"Not enough data: Median number of hits is only " << medianNumberOfHits <<
183 "! The masking continous but the mask may be empty.");
201 for (
auto sensBin = 1; sensBin <= collector_pxdhitcounts->GetXaxis()->GetNbins(); sensBin++) {
203 string sensorDescr = collector_pxdhitcounts->GetXaxis()->GetBinLabel(sensBin);
204 VxdID id(sensorDescr);
207 int nSensorHits = collector_pxdhitcounts->GetBinContent(sensBin);
208 if (nSensorHits == 0) {
214 string name = str(format(
"PXD_%1%_PixelHitmap") %
id.
getID());
215 auto collector_pxdhitmap = getObjectPtr<TH1I>(name.c_str());
216 if (collector_pxdhitmap ==
nullptr) {
217 B2WARNING(
"Cannot find PixelHitmap although there should be hits. This is strange!");
222 int nBins = collector_pxdhitmap->GetXaxis()->GetNbins();
225 if (medianNumberOfHits >=
minHits) {
229 vector<int> hitsAlongDrain(
c_nDrains, 0);
232 for (
auto bin = 1; bin <= nBins; bin++) {
237 int drainID = uCell * 4 + vCell % 4;
238 int nhits = collector_pxdhitmap->GetBinContent(bin);
239 hitsAlongDrain[drainID] += nhits;
240 hitsAlongRow[vCell] += nhits;
244 B2INFO(
"Entering masking of dead rows ...");
245 for (
auto vCell = 0; vCell <
c_nVCells; vCell++) {
247 int nhits = hitsAlongRow[vCell];
251 B2RESULT(
"Dead row with vCell=" << vCell <<
" on sensor " <<
id);
256 B2INFO(
"Entering masking of dead drains ...");
257 for (
auto drainID = 0; drainID <
c_nDrains; drainID++) {
259 int nhits = hitsAlongDrain[drainID];
263 B2RESULT(
"Dead drain line at drainID=" << drainID <<
" on sensor " <<
id);
268 B2INFO(
"Entering masking of single dead pixels ...");
269 for (
auto bin = 1; bin <= nBins; bin++) {
271 int nhits = collector_pxdhitmap->GetBinContent(bin);
275 int drainID = uCell * 4 + vCell % 4;
277 if (nhits == 0 && !deadPixelsPar->
isDeadRow(
id.getID(), vCell) && !deadPixelsPar->
isDeadDrain(
id.getID(), drainID)) {
280 B2RESULT(
"Dead single pixel with ucell=" << uCell <<
", vcell=" << vCell <<
" on sensor " <<
id);
288 vector<float> unmaskedHitsAlongDrain(
c_nDrains, 0);
289 vector<int> unmaskedCellsAlongDrain(
c_nDrains, 0);
292 vector<float> unmaskedHitsAlongRow(
c_nVCells, 0);
293 vector<int> unmaskedCellsAlongRow(
c_nVCells, 0);
297 B2RESULT(
"Pixel hit threshold is " << pixelHitThr <<
" for sensor " <<
id);
300 for (
auto bin = 1; bin <= nBins; bin++) {
305 int drainID = uCell * 4 + vCell % 4;
308 float nhits = collector_pxdhitmap->GetBinContent(bin);
311 if (nhits > pixelHitThr) {
315 B2RESULT(
"Masking single pixel with ucell=" << uCell <<
", vcell=" << vCell <<
" on sensor " <<
id);
321 ++unmaskedCellsAlongDrain[drainID];
322 unmaskedHitsAlongDrain[drainID] += nhits;
323 ++unmaskedCellsAlongRow[vCell];
324 unmaskedHitsAlongRow[vCell] += nhits;
330 B2RESULT(
"Drain hit threshold is " << drainHitThr <<
" for sensor " <<
id);
332 for (
auto drainID = 0; drainID <
c_nDrains; drainID++) {
333 if (unmaskedCellsAlongDrain[drainID] > 0) {
335 float nhits = unmaskedHitsAlongDrain[drainID] / unmaskedCellsAlongDrain[drainID];
337 if (nhits > drainHitThr) {
338 for (
auto iGate = 0; iGate < 192; iGate++) {
339 int uCell = drainID / 4;
340 int vCell = drainID % 4 + iGate * 4;
343 B2RESULT(
"Masking drain line with drainID=" << drainID <<
" on sensor " <<
id);
351 B2RESULT(
"Row hit threshold is " << rowHitThr <<
" for sensor " <<
id);
353 for (
auto vCell = 0; vCell <
c_nVCells; vCell++) {
354 if (unmaskedCellsAlongRow[vCell] > 0) {
356 float nhits = unmaskedHitsAlongRow[vCell] / unmaskedCellsAlongRow[vCell];
358 if (nhits > rowHitThr) {
359 for (
auto uCell = 0; uCell <
c_nUCells; uCell++)
362 B2RESULT(
"Masking complete row with vCell=" << vCell <<
" on sensor " <<
id);
372 int numberOfUnmaskedHits = 0;
374 for (
auto bin = 1; bin <= nBins; bin++) {
378 if (maskedPixelsPar->
pixelOK(
id.getID(), pixID)) {
379 numberOfUnmaskedHits += collector_pxdhitmap->GetBinContent(bin);
384 float meanOccupancyAfterMasking = (float)numberOfUnmaskedHits / nevents /
c_nVCells /
c_nUCells;
385 B2RESULT(
"Hotpixel filtered occupancy sensor=" <<
id <<
" is " << meanOccupancyAfterMasking);
391 auto id = elem.first;
392 auto singles = elem.second;
393 B2RESULT(
"SensorID " <<
VxdID(
id) <<
" has filtered occupancy of " << occupancyInfoPar->
getOccupancy(
id));
394 B2RESULT(
"SensorID " <<
VxdID(
id) <<
" has fraction of masked pixels of " << (
float)singles.size() / (
c_nVCells *
c_nUCells));
406 B2INFO(
"PXDHotPixelMask Calibration Successful");