Belle II Software development
SVDHotStripFinderModule.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#include <geometry/GeometryManager.h>
9#include <time.h>
10#include <vxd/geometry/GeoCache.h>
11
12// framework - Database
13#include <framework/database/IntervalOfValidity.h>
14#include <framework/database/DBImportObjPtr.h>
15
16#include <svd/calibration/SVDHotStripsCalibrations.h>
17#include <svd/calibration/SVDOccupancyCalibrations.h>
18#include <svd/modules/svdCalibration/SVDHotStripFinderModule.h>
19
20
21using namespace std;
22
23using namespace Belle2;
24
25REG_MODULE(SVDHotStripFinder);
26SVDHotStripFinderModule::SVDHotStripFinderModule() : Module()
27{
28 setDescription("The svdHotStripFinder module finds hot strips in SVD data SVDShaperDigit");
29
30 addParam("outputFileName", m_rootFileName, "Name of output root file.", std::string("SVDHotStripFinder.root"));
31 addParam("threshold", m_thr, "Threshold cut for Hot strip finder in percent", float(4.0));
32 addParam("searchBase", m_base,
33 "number of strips used to compute the average occupancy, possible choices = 32, 64, 128. Default = -1, use all sensor strips.",
34 int(-1));
35 //additional parameters to import on the DB:
36 addParam("zeroSuppression", m_zs, "ZeroSuppression cut of the input SVDShaperDigits", float(5));
37 addParam("firstExp", m_firstExp, "experiment number", int(-1));
38 addParam("firstRun", m_firstRun, "run number", int(-1));
39 addParam("lastExp", m_lastExp, "open iov", int(-1));
40 addParam("lastRun", m_lastRun, "open iov", int(-1));
41 addParam("ShaperDigits", m_ShaperDigitName, "shaper digit name", std::string(""));
42 //new parameters for HSFinderV2:
43 addParam("useHSFinderV1", m_useHSFinderV1, "Set to false only if you want to test the second version of the algorithm",
44 bool(false));
45 addParam("absOccThreshold", m_absThr,
46 "Absolute occupancy threshold: at a first loop, flag as Hot Strip (HS) all those whose occupancy > absOccThreshold", float(0.2));
47 addParam("relOccPrec", m_relOccPrec,
48 "Number of times the average sensor occupancy considered to fix the sensor dependent threshold, as for example occ_threshold = relOccPrec x occ_average",
49 float(5));
50 addParam("verbose", m_verbose, " True by default, it allows to switch off the printing of all found HS.", bool(true));
51
52}
53
54
58
59
60void SVDHotStripFinderModule::initialize()
61{
62
63 m_eventMetaData.isRequired();
64 m_storeDigits.isRequired(m_ShaperDigitName);
65
66 B2DEBUG(25, " ShaperDigits: " << m_ShaperDigitName);
67
68 m_histoList_occu = new TList;
69
70 m_rootFilePtr = new TFile(m_rootFileName.c_str(), "RECREATE");
71 if (m_useHSFinderV1)
72 B2RESULT("You are using the first version of the HSFinder algorithm (see SVDHotStripFinder::terminate in the module)");
73 else B2RESULT("You are using the modified version of the HSFinder algorithm (see SVDHotStripFinder::endRun in the module)");
74
75}
76void SVDHotStripFinderModule::beginRun()
77{
78
79 //create histograms
80
81 TH1F hOccupancy768("Occupancy768_L@layerL@ladderS@sensor@view", "Strip Occupancy of @layer.@ladder.@sensor @view/@side side", 768,
82 0,
83 768);
84 hOccupancy768.GetXaxis()->SetTitle("cellID");
85 TH1F hOccupancy512("Occupancy512_L@layerL@ladderS@sensor@view", "Strip Occupancy of @layer.@ladder.@sensor @view/@side side", 512,
86 0,
87 512);
88 hOccupancy512.GetXaxis()->SetTitle("cellID");
89 hm_occupancy = new SVDHistograms<TH1F>(hOccupancy768, hOccupancy768, hOccupancy768, hOccupancy512);
90
91 TH1F hHotStrips768("HotStrips768_L@layerL@ladderS@sensor@view", "Hot Strips of @layer.@ladder.@sensor @view/@side side", 768, 0,
92 768);
93 hHotStrips768.GetXaxis()->SetTitle("cellID");
94 TH1F hHotStrips512("HotStrips512_L@layerL@ladderS@sensor@view", "Hot Strips of @layer.@ladder.@sensor @view/@side side", 512, 0,
95 512);
96 hHotStrips512.GetXaxis()->SetTitle("cellID");
97 hm_hot_strips = new SVDHistograms<TH1F>(hHotStrips768, hHotStrips768, hHotStrips768, hHotStrips512);
98
99 TH1F hOccupancy_after768("OccupancyAfter768_L@layerL@ladderS@sensor@view",
100 "Non-Hot Strip Occupancy after HSF of @layer.@ladder.@sensor @view/@side side", 768, 0, 768);
101 hOccupancy_after768.GetXaxis()->SetTitle("cellID");
102 TH1F hOccupancy_after512("OccupancyAfter512_L@layerL@ladderS@sensor@view",
103 "Non-Hot Strip Occupancy after HSF of @layer.@ladder.@sensor @view/@side side", 512, 0, 512);
104 hOccupancy_after512.GetXaxis()->SetTitle("cellID");
105 hm_occupancy_after = new SVDHistograms<TH1F>(hOccupancy_after768, hOccupancy_after768, hOccupancy_after768, hOccupancy_after512);
106
107 TH1F hOccAll("occAll_L@layerL@ladderS@sensor@view", "Strip Occupancy Distribution of @layer.@ladder.@sensor @view/@side side", 1000,
108 0, 1);
109 hOccAll.GetXaxis()->SetTitle("occupancy");
110 hm_occAll = new SVDHistograms<TH1F>(hOccAll);
111
112 TH1F hOccHot("occHot_L@layerL@ladderS@sensor@view", "Hot Strip Occupancy Distribution of @layer.@ladder.@sensor @view/@side side",
113 10000, 0, 1);
114 hOccHot.GetXaxis()->SetTitle("occupancy");
115 hm_occHot = new SVDHistograms<TH1F>(hOccHot);
116
117 TH1F hOccAfter("occAfter_L@layerL@ladderS@sensor@view",
118 "Non-Hot Strip Occupancy Distribution of @layer.@ladder.@sensor @view/@side side", 1000, 0, 0.05);
119 hOccAfter.GetXaxis()->SetTitle("occupancy");
120 hm_occAfter = new SVDHistograms<TH1F>(hOccAfter);
121
122 //
123 TH1F hDist("dist_L@layerL@ladderS@sensor@view", "DSSD occupancy distribution of @layer.@ladder.@sensor @view/@side side", 100, 0,
124 0.05);
125 hDist.GetXaxis()->SetTitle("occupancy");
126 hm_dist = new SVDHistograms<TH1F>(hDist);
127
128 TH1F hDist1("dist1_L@layerL@ladderS@sensor@view", "DSSD true occupancy distribution of @layer.@ladder.@sensor @view/@side side",
129 100, 0, 0.05);
130 hm_dist1 = new SVDHistograms<TH1F>(hDist1);
131 hDist.GetXaxis()->SetTitle("occupancy");
132
133 TH2F hDist12("dist2d_L@layerL@ladderS@sensor@view",
134 "DSSD true vs sensor occupancy distribution of @layer.@ladder.@sensor @view/@side side", 1000, 0, 0.05, 1000, 0, 0.05);
135 hDist12.GetXaxis()->SetTitle("sensor occupancy");
136 hDist12.GetYaxis()->SetTitle("occupancy");
137 hm_dist12 = new SVDHistograms<TH2F>(hDist12);
138
139 //summary plot of the hot strips per sensor
140 m_hHotStripsSummary = new SVDSummaryPlots("hotStripsSummary@view", "Number of HotStrips on @view/@side Side");
141
142
143 // DQM style histogram number of hs vs sensor plane
144 h_tot_dqm = createHistogram1D("htodqm", "HS per sensor", 28, 0, 28.0, "HS per sensor", m_histoList_occu);
145 h_tot_dqm1 = createHistogram1D("htodqm1", "HS per sensor1", 350, 0, 350.0, "HS per sensor ", m_histoList_occu);
146
147 h_tot_dist = createHistogram1D("htotdist", "Occupancy distribution", 1000, 0, 0.05, "Relative occupancy", m_histoList_occu);
148 h_tot_dist1 = createHistogram1D("htotdist1", "True occupancy distribution", 1000, 0, 0.05, "occupancy", m_histoList_occu);
149 h_tot_dist12 = createHistogram2D("htotdist2d", "True vs sensor occupancy distribution", 1000, 0, 0.05, "sensor occupancy", 1000, 0,
150 0.05, "occupancy", m_histoList_occu);
151 h_nevents = createHistogram1D("hnevents", "Number of events", 1, 0, 1, "", m_histoList_occu);
152
153}
154
155void SVDHotStripFinderModule::event()
156{
157
158 int nDigits = m_storeDigits.getEntries();
159 h_nevents->Fill(0.0); // number of events count
160
161 if (nDigits == 0)
162 return;
163
164 //loop over the SVDShaperDigits
165 int i = 0;
166 while (i < nDigits) {
167 VxdID theVxdID = m_storeDigits[i]->getSensorID();
168 int side = m_storeDigits[i]->isUStrip();
169 int CellID = m_storeDigits[i]->getCellID();
170
171 hm_occupancy->fill(theVxdID, side, CellID);
172
173 i++;
174 }
175
176
177
178}
179
180void SVDHotStripFinderModule::endRun()
181{
182
183 int exp = m_eventMetaData->getExperiment();
184 int run = m_eventMetaData->getRun();
185
186 if (!m_useHSFinderV1) {
187
188 TDirectory* oldDir = nullptr;
189 TDirectory* dir_occuL[4] = {nullptr, nullptr, nullptr, nullptr};
190
191 //prepare ROOT FILE
192 if (m_rootFilePtr != nullptr) {
193 m_rootFilePtr->cd();
194 oldDir = gDirectory;
195 dir_occuL[0] = oldDir->mkdir("layer3");
196 dir_occuL[1] = oldDir->mkdir("layer4");
197 dir_occuL[2] = oldDir->mkdir("layer5");
198 dir_occuL[3] = oldDir->mkdir("layer6");
199 }
200
201 //Scale strip occupancy plots (per each sensor side) by the number of events. Fill SVDOccupancyCalibrations payload with the measured strip occupancy.
202
203 int hsflag[768]; //found hot strips list; hsflag[i]==1 for indetified Hot Strip
204 int nevents = h_nevents->GetEntries(); //number of events processed in events
205
206 //Define the DBObj pointers to create the needed payloads
207
208 DBImportObjPtr< SVDOccupancyCalibrations::t_payload> occDBObjPtr(SVDOccupancyCalibrations::name);
209 occDBObjPtr.construct(-99., Form("SVDOccupancy_exp%d_run%d_zs%1.1f", exp, run, m_zs));
210
211 DBImportObjPtr< SVDHotStripsCalibrations::t_payload> hotStripsDBObjPtr(SVDHotStripsCalibrations::name);
212 hotStripsDBObjPtr.construct(0, Form("SVDHotStrips_exp%d_run%d_zs%1.1f_absThr%f_relOccPrec%f", exp, run, m_zs, m_absThr,
213 m_relOccPrec));
214
215 B2RESULT("number of events " << nevents);
216
217 VXD::GeoCache& aGeometry = VXD::GeoCache::getInstance();
218 std::set<Belle2::VxdID> svdLayers = aGeometry.getLayers(VXD::SensorInfoBase::SVD);
219 std::set<Belle2::VxdID>::iterator itSvdLayers = svdLayers.begin();
220 while ((itSvdLayers != svdLayers.end())
221 && (itSvdLayers->getLayerNumber() != 7)) { //loop on Layers
222
223 std::set<Belle2::VxdID> svdLadders = aGeometry.getLadders(*itSvdLayers);
224 std::set<Belle2::VxdID>::iterator itSvdLadders = svdLadders.begin();
225
226 while (itSvdLadders != svdLadders.end()) { //loop on Ladders
227
228 std::set<Belle2::VxdID> svdSensors = aGeometry.getSensors(*itSvdLadders);
229 std::set<Belle2::VxdID>::iterator itSvdSensors = svdSensors.begin();
230
231 while (itSvdSensors != svdSensors.end()) { //loop on sensors
232
233 for (int k = 0; k < m_nSides; k ++) { //loop on Sides , k = isU(), k=0 is v-side, k=1 is u-side
234
235 /* we start indexing from 0 to avoid empty histograms */
236 int layer = itSvdSensors->getLayerNumber();
237 int ladder = itSvdSensors->getLadderNumber();
238 int sensor = itSvdSensors->getSensorNumber();
239
240 // int nafter =0; //number of good strips after first preselection cut
241 int nstrips = 768;
242 if (!k && layer != 3) nstrips = 512;
243
244 double stripOcc[768];
245 for (int i = 0; i < nstrips; i++) {stripOcc[i] = 0; hsflag[i] = 0;} //initialize vector to zero
246 double stripOccAfterAbsCut[768]; // vector of strip occupancy after first preselection based on absOccupThres cut
247 (hm_occupancy->getHistogram(*itSvdSensors, k))->Scale(1. / nevents);
248 for (int l = 0; l < nstrips; l++) {
249
250 //normalized to the total number of events to have the correct occupancy per strip and fill the corresponding dbobject
251
252 stripOcc[l] = (double)(hm_occupancy->getHistogram(*itSvdSensors, k)->GetBinContent(l + 1));
253
254 //0. Fill SVDOccupancyCalibrations Payload with the measured strip occupancy
255 occDBObjPtr->set(layer, ladder, sensor, k, l, stripOcc[l]);
256 hm_occAll->fill(*itSvdSensors, k, stripOcc[l]);
257
258 //1. Cut based on absOccupancyThreshold
259 if (stripOcc[l] > m_absThr) {
260 stripOccAfterAbsCut[l] = 0;
261 hsflag[l] = 1;
262 } else {
263 stripOccAfterAbsCut[l] = stripOcc[l];
264 hsflag[l] = 0;
265 }
266 B2DEBUG(1, "Measured strip occupancy for strip " << l << ":" << stripOccAfterAbsCut[l]);
267 }
268
269 // 2. flag hot strips that has occ_Strip > sensor average Occupancy * relOccPrec
270 bool moreHS = true;
271
272 while (moreHS && theHSFinder(stripOccAfterAbsCut, hsflag, nstrips)) {
273 moreHS = theHSFinder(stripOccAfterAbsCut, hsflag, nstrips);
274 }
275
276 //3. after second step: fill HS histograms and occupancy histograms of survived strips; fill HS payload, SVDHotStripsCalibrations
277 for (int l = 0; l < nstrips; l++) {
278 hotStripsDBObjPtr->set(layer, ladder, sensor, k, l, (int)hsflag[l]);
279 if (hsflag[l] == 0) {
280 hm_occupancy_after->getHistogram(*itSvdSensors, k)->SetBinContent(l + 1, stripOccAfterAbsCut[l]);
281 hm_occAfter->fill(*itSvdSensors, k, stripOccAfterAbsCut[l]);
282 } else {
283 hm_hot_strips->getHistogram(*itSvdSensors, k)->SetBinContent(l + 1, 1);
284 hm_occHot->fill(*itSvdSensors, k, stripOcc[l]);
285
286 TString aux_side = "V/N";
287 if (k) aux_side = "U/P";
288 if (m_verbose) B2RESULT("HS found, occupancy = " << stripOcc[l] << ", Layer: " << layer << " Ladder: " << ladder << " Sensor: "
289 << sensor <<
290 " Side: " << k << " channel: " << l);
291
292 }
293
294 }
295
296 for (int s = 0; s < hm_hot_strips->getHistogram(*itSvdSensors, k)->GetEntries(); s++)
297 m_hHotStripsSummary->fill(*itSvdSensors, k, 1);
298
299 if (m_rootFilePtr != nullptr) {
300 dir_occuL[layer - 3]->cd();
301 hm_occupancy->getHistogram(*itSvdSensors, k)->Write();
302 hm_hot_strips->getHistogram(*itSvdSensors, k)->SetLineColor(kBlack);
303 hm_hot_strips->getHistogram(*itSvdSensors, k)->SetMarkerColor(kBlack);
304 hm_hot_strips->getHistogram(*itSvdSensors, k)->SetFillStyle(3001);
305 hm_hot_strips->getHistogram(*itSvdSensors, k)->SetFillColor(kBlack);
306 hm_hot_strips->getHistogram(*itSvdSensors, k)->Write();
307 hm_occupancy_after->getHistogram(*itSvdSensors, k)->SetLineColor(kRed);
308 hm_occupancy_after->getHistogram(*itSvdSensors, k)->SetMarkerColor(kRed);
309 hm_occupancy_after->getHistogram(*itSvdSensors, k)->Write();
310 hm_occAll->getHistogram(*itSvdSensors, k)->Write();
311 hm_occHot->getHistogram(*itSvdSensors, k)->SetLineColor(kBlack);
312 hm_occHot->getHistogram(*itSvdSensors, k)->SetFillStyle(3001);
313 hm_occHot->getHistogram(*itSvdSensors, k)->SetFillColor(kBlack);
314 hm_occHot->getHistogram(*itSvdSensors, k)->SetMarkerColor(kBlack);
315 hm_occHot->getHistogram(*itSvdSensors, k)->Write();
316 hm_occAfter->getHistogram(*itSvdSensors, k)->SetLineColor(kRed);
317 hm_occAfter->getHistogram(*itSvdSensors, k)->SetMarkerColor(kRed);
318 hm_occAfter->getHistogram(*itSvdSensors, k)->Write();
319 }
320
321 B2DEBUG(1, " L" << layer << "." << ladder << "." << sensor << ".isU=" << k);
322
323 }
324 ++itSvdSensors;
325 }
326 ++itSvdLadders;
327 }
328 ++itSvdLayers;
329 }
330
331 if (m_rootFilePtr != nullptr) {
332 oldDir->cd();
333 m_hHotStripsSummary->getHistogram(0)->Write();
334 m_hHotStripsSummary->getHistogram(1)->Write();
335 m_rootFilePtr->Close();
336 }
337
338 //import the filled dbobjects to the ConditionDB
339 if (m_firstExp == -1)
340 m_firstExp = exp;
341 if (m_lastExp == -1)
342 m_lastExp = exp;
343 if (m_firstRun == -1)
344 m_firstRun = run;
345 if (m_lastRun == -1)
346 m_lastRun = run;
347
348 IntervalOfValidity iov(m_firstExp, m_firstRun, m_lastExp, m_lastRun);
349 occDBObjPtr.import(iov);
350 hotStripsDBObjPtr.import(iov);
351 B2RESULT("Imported to database.");
352 }
353}
354
355void SVDHotStripFinderModule::terminate()
356{
357 if (m_useHSFinderV1) {
358 TDirectory* oldDir = nullptr;
359
360 TDirectory* dir_occuL[4] = {nullptr, nullptr, nullptr, nullptr};
361
362 //prepare ROOT FILE
363 if (m_rootFilePtr != nullptr) {
364 m_rootFilePtr->cd();
365 oldDir = gDirectory;
366 dir_occuL[0] = oldDir->mkdir("layer3");
367 dir_occuL[1] = oldDir->mkdir("layer4");
368 dir_occuL[2] = oldDir->mkdir("layer5");
369 dir_occuL[3] = oldDir->mkdir("layer6");
370 }
371
372 /**************************************************************************
373 * Hotstrips finding algorithm *
374 ***************************************************************************/
375
376
377 //Find low charged clusters with high occupancy.
378 int flag[768]; // list of working (non zero) strips
379 int hsflag[768]; //found hot strips list; hsflag[i]==1 for indetified Hot Strip
380 int nevents = h_nevents->GetEntries(); //number of events processed in events
381 int ibase = 768; // interval used for the hot strip finding
382 if (m_base != -1)
383 ibase = m_base;
384
385 B2DEBUG(1, "number of events " << nevents);
386
387 VXD::GeoCache& aGeometry = VXD::GeoCache::getInstance();
388 std::set<Belle2::VxdID> svdLayers = aGeometry.getLayers(VXD::SensorInfoBase::SVD);
389 std::set<Belle2::VxdID>::iterator itSvdLayers = svdLayers.begin();
390 int itsensor = 0; //sensor numbering
391 while ((itSvdLayers != svdLayers.end())
392 && (itSvdLayers->getLayerNumber() != 7)) { //loop on Layers
393
394 std::set<Belle2::VxdID> svdLadders = aGeometry.getLadders(*itSvdLayers);
395 std::set<Belle2::VxdID>::iterator itSvdLadders = svdLadders.begin();
396
397 while (itSvdLadders != svdLadders.end()) { //loop on Ladders
398
399 std::set<Belle2::VxdID> svdSensors = aGeometry.getSensors(*itSvdLadders);
400 std::set<Belle2::VxdID>::iterator itSvdSensors = svdSensors.begin();
401
402 while (itSvdSensors != svdSensors.end()) { //loop on sensors
403
404 for (int k = 0; k < m_nSides; k ++) { //loop on Sides
405
406 // we start indexing from 0 to avoid empty histograms
407 int i = itSvdSensors->getLayerNumber() - 3;
408 int m = itSvdSensors->getLadderNumber() - 1;
409 int j = itSvdSensors->getSensorNumber() - 1;
410 float position1[768]; // vector of hits in the sensor
411 float nCltrk[24]; // index to interval if we search in smaller intervals then full sensor
412 int it = 0;
413 int iths = 0;
414 // it is safer to initialize the nCltrk vector
415 for (int l = 0; l < 24; l++) {
416 nCltrk[l] = 0.0;
417 }
418
419 for (int l = 0; l < 768; l++) {
420
421 position1[l] = hm_occupancy->getHistogram(*itSvdSensors, k)->GetBinContent(l + 1);
422 //if no hits in strip, mark the strip as bad
423
424 if (position1[l] == 0) { flag[l] = 0;}
425 else {
426 flag[l] = 1;
427 it++; //number of good (non zero ) strips
428 // find in which interval the strip lays
429 div_t test = div(l, ibase);
430 nCltrk[test.quot] = nCltrk[test.quot] + position1[l]; // number of entries in given interval
431 }
432
433 }
434
435 for (int l = 0; l < 768; l++) {
436 div_t test = div(l, ibase);
437
438 // tmp_occ - relative occupancy for Hot Strip search interval, for channel l
439 // tmp_occ1 - occupancy
440
441 float tmp_occ = position1[l] / (float)nCltrk[test.quot]; //for hot strip search
442 float tmp_occ1 = position1[l] / (float)nevents; //for SVDOccupancyCalibration
443 position1[l] = tmp_occ; //vector used for Hot strip search<
444 if (tmp_occ > 0.0) {
445 hm_dist->fill(*itSvdSensors, k, tmp_occ); // ..
446 h_tot_dist->Fill(tmp_occ);
447 hm_dist1->fill(*itSvdSensors, k, tmp_occ1); //occupancy as probability to fire the strip
448 h_tot_dist1->Fill(tmp_occ1);
449 hm_dist12->fill(*itSvdSensors, k, tmp_occ, tmp_occ1); // 2D distribution
450 h_tot_dist12->Fill(tmp_occ, tmp_occ1);
451 }
452 }
453 float occupancy[24]; //occupancy for second pass
454 for (int l = 0; l < 24; l++) {
455 occupancy[l] = 0.0;
456 }
457 int it1st = it; //first pass:number of good strips
458 it = 0;
459 // first pass
460 for (int l = 0; l < 768; l++) {
461 div_t test = div(l, ibase);
462 float threshold_corrections = 1.0;
463 /*
464 threshold is corrected by the real number of alive strips
465 */
466 threshold_corrections = threshold_corrections * sqrt(768.0 / (float)it1st);
467 if (ibase == 32) threshold_corrections = 24.0;
468 if (ibase == 64) threshold_corrections = 12.0;
469 if (ibase == 128) threshold_corrections = 6.0;
470
471 if (position1[l] > 0.01 * m_thr * threshold_corrections) { // if probability is larger then threshold mark as Hot strip
472 hsflag[l] = 1; // HS vector
473 flag[l] = 0; // mark strip as bad for second pass
474 iths++;
475 B2RESULT("1st pass HS found! Layer: " << i + 3 << " Ladder: " << m << " Sensor: " << j << " Side: " << k << " channel: " << l);
476 } else {
477 hsflag[l] = 0; // not a HS
478 //recalculate the occupancy in DSSD only for good strip after first pass
479 occupancy[test.quot] = occupancy[test.quot] + hm_occupancy->getHistogram(*itSvdSensors, k)->GetBinContent(l + 1);
480 it++; //number of good strips after first pass
481 }
482
483 }
484 /*
485 Second pass of Hot strip finder, After the first pass we remove already found Host strips we do the second pass with the same threshold
486 */
487 // second pass
488 for (int l = 0; l < 768; l++) {
489 div_t test = div(l, ibase);
490 position1[l] = position1[l] * nCltrk[test.quot] / (float)occupancy[test.quot];
491 float threshold_corrections = 1.0;
492 threshold_corrections = threshold_corrections * sqrt(768.0 / (float)it);
493 if (ibase == 32) threshold_corrections = 24.0;
494 if (ibase == 64) threshold_corrections = 12.0;
495 if (ibase == 128) threshold_corrections = 6.0;
496
497 if ((flag[l]) && (position1[l] > 0.01 * m_thr * threshold_corrections)) { //HS
498 hsflag[l] = 1;// HS vector
499 flag[l] = 0; // mark strip as bad after second pass
500 iths++;
501 B2RESULT("2nd pass HS FOUND! Layer: " << i + 3 << " Ladder: " << m << " Sensor: " << j << " Side: " << k << " channel: " << l);
502 }
503 }
504 // for Laura .. HS flags, place interface for DB
505 // outputs : hsflag[l] 1- HS 0- non HS, flag[l] 0- bad strip, 1 working strip, h_tot_dist1 occupancy as probability of firing the strip
506
507
508 for (int l = 0; l < 768; l++) {
509
510 B2DEBUG(1, hsflag[l]);
511
512 float tmpOcc = hm_occupancy->getHistogram(*itSvdSensors, k)->GetBinContent(l + 1) / (double)nevents;
513 hm_occAll->fill(*itSvdSensors, k, tmpOcc);
514
515 if (hsflag[l] == 0) {
516 hm_occupancy_after->getHistogram(*itSvdSensors, k)->SetBinContent(l + 1, tmpOcc); //alive strips without identified HS
517 hm_occAfter->fill(*itSvdSensors, k, tmpOcc);
518 } else {
519 hm_hot_strips->getHistogram(*itSvdSensors, k)->SetBinContent(l + 1, 1);
520 hm_occHot->fill(*itSvdSensors, k, tmpOcc);
521 }
522 }
523
524 if (m_rootFilePtr != nullptr) {
525 hm_occupancy->getHistogram(*itSvdSensors, k)->Scale(1.0 / (double)nevents);
526
527 dir_occuL[i]->cd();
528 hm_occupancy->getHistogram(*itSvdSensors, k)->Write();
529 hm_hot_strips->getHistogram(*itSvdSensors, k)->SetLineColor(kBlack);
530 hm_hot_strips->getHistogram(*itSvdSensors, k)->SetFillStyle(3001);
531 hm_hot_strips->getHistogram(*itSvdSensors, k)->SetFillColor(kBlack);
532 hm_hot_strips->getHistogram(*itSvdSensors, k)->SetMarkerColor(kBlack);
533 hm_hot_strips->getHistogram(*itSvdSensors, k)->Write();
534 hm_occupancy_after->getHistogram(*itSvdSensors, k)->SetLineColor(kRed);
535 hm_occupancy_after->getHistogram(*itSvdSensors, k)->SetMarkerColor(kRed);
536 hm_occupancy_after->getHistogram(*itSvdSensors, k)->Write();
537 hm_occAll->getHistogram(*itSvdSensors, k)->Write();
538 hm_occHot->getHistogram(*itSvdSensors, k)->SetLineColor(kBlack);
539 hm_occHot->getHistogram(*itSvdSensors, k)->SetFillStyle(3001);
540 hm_occHot->getHistogram(*itSvdSensors, k)->SetFillColor(kBlack);
541 hm_occHot->getHistogram(*itSvdSensors, k)->SetMarkerColor(kBlack);
542 hm_occHot->getHistogram(*itSvdSensors, k)->Write();
543 hm_occAfter->getHistogram(*itSvdSensors, k)->SetLineColor(kRed);
544 hm_occAfter->getHistogram(*itSvdSensors, k)->SetMarkerColor(kRed);
545 hm_occAfter->getHistogram(*itSvdSensors, k)->Write();
546 // hm_dist12->getHistogram(*itSvdSensors, k)->Write();
547 }
548
549 B2DEBUG(1, " side " << i << " " << j << " " << m << " " << k);
550 /* end */
551
552
553 // store number hot strips per sensor
554
555
556 for (int iy = 0; iy < iths; iy++) {
557 h_tot_dqm->Fill(float(itsensor));
558 h_tot_dqm1->Fill(float(itsensor));
559 }
560
561 for (int s = 0; s < hm_hot_strips->getHistogram(*itSvdSensors, k)->GetEntries(); s++)
562 m_hHotStripsSummary->fill(*itSvdSensors, k, 1);
563
564 itsensor++;
565 }
566 ++itSvdSensors;
567 }
568 ++itSvdLadders;
569 }
570 ++itSvdLayers;
571
572 }
573
574 if (m_rootFilePtr != nullptr) {
575 oldDir->cd();
576
577 m_hHotStripsSummary->getHistogram(0)->Write();
578 m_hHotStripsSummary->getHistogram(1)->Write();
579
580 TObject* obj;
581 TIter nextH_occu(m_histoList_occu);
582 while ((obj = nextH_occu()))
583 obj->Write();
584
585 m_rootFilePtr->Close();
586 }
587
588 }
589}
590
591TH1F* SVDHotStripFinderModule::createHistogram1D(const char* name, const char* title,
592 Int_t nbins, Double_t min, Double_t max,
593 const char* xtitle, TList* histoList)
594{
595
596 TH1F* h = new TH1F(name, title, nbins, min, max);
597
598 h->GetXaxis()->SetTitle(xtitle);
599
600 if (histoList)
601 histoList->Add(h);
602
603
604 return h;
605}
606
607
608TH2F* SVDHotStripFinderModule::createHistogram2D(const char* name, const char* title,
609 Int_t nbinsX, Double_t minX, Double_t maxX,
610 const char* titleX,
611 Int_t nbinsY, Double_t minY, Double_t maxY,
612 const char* titleY, TList* histoList)
613{
614
615 TH2F* h = new TH2F(name, title, nbinsX, minX, maxX, nbinsY, minY, maxY);
616 h->GetXaxis()->SetTitle(titleX);
617 h->GetYaxis()->SetTitle(titleY);
618
619 if (histoList)
620 histoList->Add(h);
621
622 return h;
623}
624
625
626bool SVDHotStripFinderModule::theHSFinder(double* stripOccAfterAbsCut, int* hsflag, int nstrips)
627{
628 bool found = false;
629
630 if (m_base == -1)
631 m_base = nstrips;
632
633 int N = nstrips / m_base;
634
635 for (int sector = 0; sector < N; sector++) {
636
637 int nafter = 0;
638 double sensorOccAverage = 0;
639
640 for (int l = sector * m_base; l < sector * m_base + m_base; l++) {
641 sensorOccAverage = sensorOccAverage + stripOccAfterAbsCut[l];
642 if (stripOccAfterAbsCut[l] > 0) nafter++;
643 }
644 sensorOccAverage = sensorOccAverage / nafter;
645
646 B2DEBUG(1, "Average occupancy: " << sensorOccAverage);
647
648 for (int l = sector * m_base; l < sector * m_base + m_base; l++) {
649
650 // flag additional HS by comparing each strip occupancy with the sensor-based average occupancy
651
652 if (stripOccAfterAbsCut[l] > sensorOccAverage * m_relOccPrec) {
653 hsflag[l] = 1;
654 found = true;
655 stripOccAfterAbsCut[l] = 0;
656 }
657 // else hsflag[l]=0;
658 }
659 }
660
661 return found;
662}
Belle2::DBImportBase::import
bool import(const IntervalOfValidity &iov)
Import the object to database.
Definition DBImportBase.cc:36
Belle2::DBImportObjPtr
Class for importing a single object to the database.
Definition DBImportObjPtr.h:23
Belle2::DBImportObjPtr::construct
void construct(Args &&... params)
Construct an object of type T in this DBImportObjPtr using the provided constructor arguments.
Definition DBImportObjPtr.h:47
Belle2::IntervalOfValidity
A class that describes the interval of experiments/runs for which an object in the database is valid.
Definition IntervalOfValidity.h:25
Belle2::Module::setDescription
void setDescription(const std::string &description)
Sets the description of the module.
Definition Module.cc:214
Belle2::Module::Module
Module()
Constructor.
Definition Module.cc:30
Belle2::SVDHistograms
template class for SVd histograms
Definition SVDHistograms.h:24
Belle2::SVDHotStripFinderModule::m_rootFileName
std::string m_rootFileName
root file name
Definition SVDHotStripFinderModule.h:78
Belle2::SVDHotStripFinderModule::createHistogram1D
TH1F * createHistogram1D(const char *name, const char *title, Int_t nbins, Double_t min, Double_t max, const char *xtitle, TList *histoList=nullptr)
create 1D histograms
Definition SVDHotStripFinderModule.cc:591
Belle2::SVDHotStripFinderModule::m_storeDigits
StoreArray< SVDShaperDigit > m_storeDigits
shaper digits store array
Definition SVDHotStripFinderModule.h:76
Belle2::SVDHotStripFinderModule::m_absThr
float m_absThr
Absolute Occupancy Threshold cut for Hot strip finder.
Definition SVDHotStripFinderModule.h:83
Belle2::SVDHotStripFinderModule::hm_occupancy
SVDHistograms< TH1F > * hm_occupancy
strip occupancy per sensor
Definition SVDHotStripFinderModule.h:105
Belle2::SVDHotStripFinderModule::SVDHotStripFinderModule
SVDHotStripFinderModule()
Constructor, for setting module description and parameters.
Definition SVDHotStripFinderModule.cc:26
Belle2::SVDHotStripFinderModule::m_ShaperDigitName
std::string m_ShaperDigitName
shaper digits name
Definition SVDHotStripFinderModule.h:79
Belle2::SVDHotStripFinderModule::hm_hot_strips
SVDHistograms< TH1F > * hm_hot_strips
hot strips per sensor
Definition SVDHotStripFinderModule.h:106
Belle2::SVDHotStripFinderModule::initialize
virtual void initialize() override
Use this to initialize resources or memory your module needs.
Definition SVDHotStripFinderModule.cc:60
Belle2::SVDHotStripFinderModule::hm_dist12
SVDHistograms< TH2F > * hm_dist12
true occupancy VS sensor histograms
Definition SVDHotStripFinderModule.h:114
Belle2::SVDHotStripFinderModule::event
virtual void event() override
Called once for each event.
Definition SVDHotStripFinderModule.cc:155
Belle2::SVDHotStripFinderModule::m_histoList_occu
TList * m_histoList_occu
occupancy for low charge clusters
Definition SVDHotStripFinderModule.h:103
Belle2::SVDHotStripFinderModule::m_relOccPrec
float m_relOccPrec
Relative precision on occupancy which is defined to be negligible for the hit background rate estimat...
Definition SVDHotStripFinderModule.h:84
Belle2::SVDHotStripFinderModule::endRun
virtual void endRun() override
Called once when a run ends.
Definition SVDHotStripFinderModule.cc:180
Belle2::SVDHotStripFinderModule::h_tot_dist1
TH1F * h_tot_dist1
absolute occupany histogram
Definition SVDHotStripFinderModule.h:121
Belle2::SVDHotStripFinderModule::hm_dist1
SVDHistograms< TH1F > * hm_dist1
true occupancy histograms
Definition SVDHotStripFinderModule.h:113
Belle2::SVDHotStripFinderModule::terminate
virtual void terminate() override
Clean up anything you created in initialize().
Definition SVDHotStripFinderModule.cc:355
Belle2::SVDHotStripFinderModule::m_hHotStripsSummary
SVDSummaryPlots * m_hHotStripsSummary
hot strip summary histo
Definition SVDHotStripFinderModule.h:116
Belle2::SVDHotStripFinderModule::theHSFinder
bool theHSFinder(double *stripOccAfterAbsCut, int *hsflag, int nstrips)
return true if the strip is hot
Definition SVDHotStripFinderModule.cc:626
Belle2::SVDHotStripFinderModule::m_nSides
static const int m_nSides
number of sides
Definition SVDHotStripFinderModule.h:101
Belle2::SVDHotStripFinderModule::m_eventMetaData
StoreObjPtr< EventMetaData > m_eventMetaData
event meta data store array
Definition SVDHotStripFinderModule.h:77
Belle2::SVDHotStripFinderModule::hm_dist
SVDHistograms< TH1F > * hm_dist
occupancy histograms
Definition SVDHotStripFinderModule.h:112
Belle2::SVDHotStripFinderModule::beginRun
virtual void beginRun() override
Called once before a new run begins.
Definition SVDHotStripFinderModule.cc:76
Belle2::SVDHotStripFinderModule::hm_occAll
SVDHistograms< TH1F > * hm_occAll
occupancy distribution - all strips
Definition SVDHotStripFinderModule.h:108
Belle2::SVDHotStripFinderModule::m_zs
float m_zs
zero suppression cut for the input shaper digits
Definition SVDHotStripFinderModule.h:91
Belle2::SVDHotStripFinderModule::~SVDHotStripFinderModule
virtual ~SVDHotStripFinderModule()
Use to clean up anything you created in the constructor.
Definition SVDHotStripFinderModule.cc:55
Belle2::SVDHotStripFinderModule::h_nevents
TH1F * h_nevents
number of events counting
Definition SVDHotStripFinderModule.h:123
Belle2::SVDHotStripFinderModule::h_tot_dist
TH1F * h_tot_dist
relative occupancy histogram
Definition SVDHotStripFinderModule.h:120
Belle2::SVDHotStripFinderModule::hm_occAfter
SVDHistograms< TH1F > * hm_occAfter
occupancy distribution - not hot strips
Definition SVDHotStripFinderModule.h:110
Belle2::SVDHotStripFinderModule::h_tot_dist12
TH2F * h_tot_dist12
2d distributiuons of occupancies
Definition SVDHotStripFinderModule.h:122
Belle2::SVDHotStripFinderModule::hm_occHot
SVDHistograms< TH1F > * hm_occHot
occupancy distribution - hot strips
Definition SVDHotStripFinderModule.h:109
Belle2::SVDHotStripFinderModule::h_tot_dqm
TH1F * h_tot_dqm
number of hot strips per sensor
Definition SVDHotStripFinderModule.h:118
Belle2::SVDHotStripFinderModule::createHistogram2D
TH2F * createHistogram2D(const char *name, const char *title, Int_t nbinsX, Double_t minX, Double_t maxX, const char *titleX, Int_t nbinsY, Double_t minY, Double_t maxY, const char *titleY, TList *histoList=nullptr)
create 2D histograms
Definition SVDHotStripFinderModule.cc:608
Belle2::SVDHotStripFinderModule::hm_occupancy_after
SVDHistograms< TH1F > * hm_occupancy_after
strip occupancy after removal of hot strips, per sensor
Definition SVDHotStripFinderModule.h:107
Belle2::SVDHotStripFinderModule::m_rootFilePtr
TFile * m_rootFilePtr
pointer at root file used for storing histograms
Definition SVDHotStripFinderModule.h:88
Belle2::SVDHotStripFinderModule::h_tot_dqm1
TH1F * h_tot_dqm1
number of hot strips per sensor for layer 3
Definition SVDHotStripFinderModule.h:119
Belle2::SVDHotStripFinderModule::m_verbose
bool m_verbose
False by default, it allows to switch on the printing of all found HS.
Definition SVDHotStripFinderModule.h:85
Belle2::SVDHotStripsCalibrations::name
static std::string name
name of the SVDHotStripsCalibrations payload
Definition SVDHotStripsCalibrations.h:34
Belle2::SVDOccupancyCalibrations::name
static std::string name
name of the ccupancy payload
Definition SVDOccupancyCalibrations.h:32
Belle2::SVDSummaryPlots
class to summarize SVD quantities per sensor and side
Definition SVDSummaryPlots.h:26
Belle2::VXD::GeoCache
Class to facilitate easy access to sensor information of the VXD like coordinate transformations or p...
Definition GeoCache.h:39
Belle2::VXD::GeoCache::getLayers
const std::set< Belle2::VxdID > getLayers(SensorInfoBase::SensorType sensortype=SensorInfoBase::VXD)
Return a set of all known Layers.
Definition GeoCache.cc:176
Belle2::VXD::GeoCache::getSensors
const std::set< Belle2::VxdID > & getSensors(Belle2::VxdID ladder) const
Return a set of all sensor IDs belonging to a given ladder.
Definition GeoCache.cc:204
Belle2::VXD::GeoCache::getInstance
static GeoCache & getInstance()
Return a reference to the singleton instance.
Definition GeoCache.cc:214
Belle2::VXD::GeoCache::getLadders
const std::set< Belle2::VxdID > & getLadders(Belle2::VxdID layer) const
Return a set of all ladder IDs belonging to a given layer.
Definition GeoCache.cc:193
Belle2::VxdID
Class to uniquely identify a any structure of the PXD and SVD.
Definition VxdID.h:33
Belle2::Module::addParam
void addParam(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
Definition Module.h:559
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition Module.h:649
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition beamHelpers.h:28
Belle2
Abstract base class for different kinds of events.
Definition MillepedeAlgorithm.h:17
std
STL namespace.