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