Belle II Software  release-06-01-15
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 
21 using namespace std;
22 
23 using namespace Belle2;
24 
25 REG_MODULE(SVDHotStripFinder)
26 SVDHotStripFinderModule::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 paramters 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 
55 SVDHotStripFinderModule::~SVDHotStripFinderModule()
56 {
57 }
58 
59 
60 void 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 }
76 void 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 historgram 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 
155 void 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 
180 void 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  int itsensor = 0; //sensor numbering
221  while ((itSvdLayers != svdLayers.end())
222  && (itSvdLayers->getLayerNumber() != 7)) { //loop on Layers
223 
224  std::set<Belle2::VxdID> svdLadders = aGeometry.getLadders(*itSvdLayers);
225  std::set<Belle2::VxdID>::iterator itSvdLadders = svdLadders.begin();
226 
227  while (itSvdLadders != svdLadders.end()) { //loop on Ladders
228 
229  std::set<Belle2::VxdID> svdSensors = aGeometry.getSensors(*itSvdLadders);
230  std::set<Belle2::VxdID>::iterator itSvdSensors = svdSensors.begin();
231 
232  while (itSvdSensors != svdSensors.end()) { //loop on sensors
233 
234  for (int k = 0; k < m_nSides; k ++) { //loop on Sides , k = isU(), k=0 is v-side, k=1 is u-side
235 
236  /* we start indexing from 0 to avoid empty histograms */
237  int layer = itSvdSensors->getLayerNumber();
238  int ladder = itSvdSensors->getLadderNumber();
239  int sensor = itSvdSensors->getSensorNumber();
240 
241  // int nafter =0; //number of good strips after first preselection cut
242  int nstrips = 768;
243  if (!k && layer != 3) nstrips = 512;
244 
245  double stripOcc[768];
246  for (int i = 0; i < nstrips; i++) {stripOcc[i] = 0; hsflag[i] = 0;} //initialize vector to zero
247  double stripOccAfterAbsCut[768]; // vector of strip occupancy after first preselection based on absOccupThres cut
248  (hm_occupancy->getHistogram(*itSvdSensors, k))->Scale(1. / nevents);
249  for (int l = 0; l < nstrips; l++) {
250 
251  //normalized to the total number of events to have the correct occupancy per strip and fill the corresponding dbobject
252 
253  stripOcc[l] = (double)(hm_occupancy->getHistogram(*itSvdSensors, k)->GetBinContent(l + 1));
254 
255  //0. Fill SVDOccupancyCalibrations Payload with the measured strip occupancy
256  occDBObjPtr->set(layer, ladder, sensor, k, l, stripOcc[l]);
257  hm_occAll->fill(*itSvdSensors, k, stripOcc[l]);
258 
259  //1. Cut based on absOccupancyThreshold
260  if (stripOcc[l] > m_absThr) {
261  stripOccAfterAbsCut[l] = 0;
262  hsflag[l] = 1;
263  } else {
264  stripOccAfterAbsCut[l] = stripOcc[l];
265  hsflag[l] = 0;
266  }
267  B2DEBUG(1, "Measured strip occupancy for strip " << l << ":" << stripOccAfterAbsCut[l]);
268  }
269 
270  // 2. flag hot strips that has occ_Strip > sensor average Occupancy * relOccPrec
271  bool moreHS = true;
272 
273  while (moreHS && theHSFinder(stripOccAfterAbsCut, hsflag, nstrips)) {
274  moreHS = theHSFinder(stripOccAfterAbsCut, hsflag, nstrips);
275  }
276 
277  //3. after second step: fill HS histograms and occupancy histograms of survived strips; fill HS payload, SVDHotStripsCalibrations
278  for (int l = 0; l < nstrips; l++) {
279  hotStripsDBObjPtr->set(layer, ladder, sensor, k, l, (int)hsflag[l]);
280  if (hsflag[l] == 0) {
281  hm_occupancy_after->getHistogram(*itSvdSensors, k)->SetBinContent(l + 1 , stripOccAfterAbsCut[l]);
282  hm_occAfter->fill(*itSvdSensors, k, stripOccAfterAbsCut[l]);
283  } else {
284  hm_hot_strips->getHistogram(*itSvdSensors, k)->SetBinContent(l + 1, 1);
285  hm_occHot->fill(*itSvdSensors, k, stripOcc[l]);
286 
287  TString aux_side = "V/N";
288  if (k) aux_side = "U/P";
289  if (m_verbose) B2RESULT("HS found, occupancy = " << stripOcc[l] << ", Layer: " << layer << " Ladder: " << ladder << " Sensor: "
290  << sensor <<
291  " Side: " << k << " channel: " << l);
292 
293  }
294 
295  }
296 
297  for (int s = 0; s < hm_hot_strips->getHistogram(*itSvdSensors, k)->GetEntries(); s++)
298  m_hHotStripsSummary->fill(*itSvdSensors, k, 1);
299 
300  if (m_rootFilePtr != nullptr) {
301  dir_occuL[layer - 3]->cd();
302  hm_occupancy->getHistogram(*itSvdSensors, k)->Write();
303  hm_hot_strips->getHistogram(*itSvdSensors, k)->SetLineColor(kBlack);
304  hm_hot_strips->getHistogram(*itSvdSensors, k)->SetMarkerColor(kBlack);
305  hm_hot_strips->getHistogram(*itSvdSensors, k)->SetFillStyle(3001);
306  hm_hot_strips->getHistogram(*itSvdSensors, k)->SetFillColor(kBlack);
307  hm_hot_strips->getHistogram(*itSvdSensors, k)->Write();
308  hm_occupancy_after->getHistogram(*itSvdSensors, k)->SetLineColor(kRed);
309  hm_occupancy_after->getHistogram(*itSvdSensors, k)->SetMarkerColor(kRed);
310  hm_occupancy_after->getHistogram(*itSvdSensors, k)->Write();
311  hm_occAll->getHistogram(*itSvdSensors, k)->Write();
312  hm_occHot->getHistogram(*itSvdSensors, k)->SetLineColor(kBlack);
313  hm_occHot->getHistogram(*itSvdSensors, k)->SetFillStyle(3001);
314  hm_occHot->getHistogram(*itSvdSensors, k)->SetFillColor(kBlack);
315  hm_occHot->getHistogram(*itSvdSensors, k)->SetMarkerColor(kBlack);
316  hm_occHot->getHistogram(*itSvdSensors, k)->Write();
317  hm_occAfter->getHistogram(*itSvdSensors, k)->SetLineColor(kRed);
318  hm_occAfter->getHistogram(*itSvdSensors, k)->SetMarkerColor(kRed);
319  hm_occAfter->getHistogram(*itSvdSensors, k)->Write();
320  }
321 
322  B2DEBUG(1, " L" << layer << "." << ladder << "." << sensor << ".isU=" << k);
323 
324  itsensor++;
325  }
326  ++itSvdSensors;
327  }
328  ++itSvdLadders;
329  }
330  ++itSvdLayers;
331  }
332 
333  if (m_rootFilePtr != nullptr) {
334  oldDir->cd();
335  m_hHotStripsSummary->getHistogram(0)->Write();
336  m_hHotStripsSummary->getHistogram(1)->Write();
337  }
338 
339  m_rootFilePtr->Close();
340  //import the filled dbobjects to the ConditionDB
341  if (m_firstExp == -1)
342  m_firstExp = exp;
343  if (m_lastExp == -1)
344  m_lastExp = exp;
345  if (m_firstRun == -1)
346  m_firstRun = run;
347  if (m_lastRun == -1)
348  m_lastRun = run;
349 
350  IntervalOfValidity iov(m_firstExp, m_firstRun, m_lastExp, m_lastRun);
351  occDBObjPtr.import(iov);
352  hotStripsDBObjPtr.import(iov);
353  B2RESULT("Imported to database.");
354  }
355 }
356 
357 void SVDHotStripFinderModule::terminate()
358 {
359  if (m_useHSFinderV1) {
360  TDirectory* oldDir = nullptr;
361 
362  TDirectory* dir_occuL[4] = {nullptr, nullptr, nullptr, nullptr};
363 
364  //prepare ROOT FILE
365  if (m_rootFilePtr != nullptr) {
366  m_rootFilePtr->cd();
367  oldDir = gDirectory;
368  dir_occuL[0] = oldDir->mkdir("layer3");
369  dir_occuL[1] = oldDir->mkdir("layer4");
370  dir_occuL[2] = oldDir->mkdir("layer5");
371  dir_occuL[3] = oldDir->mkdir("layer6");
372  }
373 
374  /**************************************************************************
375  * Hotstrips finding algorithm *
376  ***************************************************************************/
377 
378 
379  //Find low charged clusters with high occupancy.
380  int flag[768]; // list of working (non zero) strips
381  int hsflag[768]; //found hot strips list; hsflag[i]==1 for indetified Hot Strip
382  int nevents = h_nevents->GetEntries(); //number of events processed in events
383  int ibase = 768; // interval used for the hot strip finding
384  if (m_base != -1)
385  ibase = m_base;
386 
387  B2DEBUG(1, "number of events " << nevents);
388 
389  VXD::GeoCache& aGeometry = VXD::GeoCache::getInstance();
390  std::set<Belle2::VxdID> svdLayers = aGeometry.getLayers(VXD::SensorInfoBase::SVD);
391  std::set<Belle2::VxdID>::iterator itSvdLayers = svdLayers.begin();
392  int itsensor = 0; //sensor numbering
393  while ((itSvdLayers != svdLayers.end())
394  && (itSvdLayers->getLayerNumber() != 7)) { //loop on Layers
395 
396  std::set<Belle2::VxdID> svdLadders = aGeometry.getLadders(*itSvdLayers);
397  std::set<Belle2::VxdID>::iterator itSvdLadders = svdLadders.begin();
398 
399  while (itSvdLadders != svdLadders.end()) { //loop on Ladders
400 
401  std::set<Belle2::VxdID> svdSensors = aGeometry.getSensors(*itSvdLadders);
402  std::set<Belle2::VxdID>::iterator itSvdSensors = svdSensors.begin();
403 
404  while (itSvdSensors != svdSensors.end()) { //loop on sensors
405 
406  for (int k = 0; k < m_nSides; k ++) { //loop on Sides
407 
408  // we start indexing from 0 to avoid empty histograms
409  int i = itSvdSensors->getLayerNumber() - 3;
410  int m = itSvdSensors->getLadderNumber() - 1;
411  int j = itSvdSensors->getSensorNumber() - 1;
412  float position1[768]; // vector of hits in the sensor
413  float nCltrk[24]; // index to interval if we search in smaller intervals then full sensor
414  int it = 0;
415  int iths = 0;
416  // it is safer to initialize the nCltrk vector
417  for (int l = 0; l < 24; l++) {
418  nCltrk[l] = 0.0;
419  }
420 
421  for (int l = 0; l < 768; l++) {
422 
423  position1[l] = hm_occupancy->getHistogram(*itSvdSensors, k)->GetBinContent(l + 1);
424  //if no hits in strip, mark the strip as bad
425 
426  if (position1[l] == 0) { flag[l] = 0;}
427  else {
428  flag[l] = 1;
429  it++; //number of good (non zero ) strips
430  // find in which interval the strip lays
431  div_t test = div(l, ibase);
432  nCltrk[test.quot] = nCltrk[test.quot] + position1[l]; // number of entries in given interval
433  }
434 
435  }
436 
437  for (int l = 0; l < 768; l++) {
438  div_t test = div(l, ibase);
439 
440  // tmp_occ - relative occupancy for Hot Strip search interval, for channel l
441  // tmp_occ1 - occupancy
442 
443  float tmp_occ = position1[l] / (float)nCltrk[test.quot]; //for hot strip search
444  float tmp_occ1 = position1[l] / (float)nevents; //for SVDOccupancyCalibration
445  position1[l] = tmp_occ; //vector used for Hot strip search<
446  if (tmp_occ > 0.0) {
447  hm_dist->fill(*itSvdSensors, k, tmp_occ); // ..
448  h_tot_dist->Fill(tmp_occ);
449  hm_dist1->fill(*itSvdSensors, k, tmp_occ1); //occupancy as probablity to fire the strip
450  h_tot_dist1->Fill(tmp_occ1);
451  hm_dist12->fill(*itSvdSensors, k, tmp_occ, tmp_occ1); // 2D distribution
452  h_tot_dist12->Fill(tmp_occ, tmp_occ1);
453  }
454  }
455  float occupancy[24]; //occupancy for second pass
456  for (int l = 0; l < 24; l++) {
457  occupancy[l] = 0.0;
458  }
459  int it1st = it; //first pass:number of good strips
460  it = 0;
461  // first pass
462  for (int l = 0; l < 768; l++) {
463  div_t test = div(l, ibase);
464  float threshold_corrections = 1.0;
465  /*
466  threshold is corrected by the real number of alive strips
467  */
468  threshold_corrections = threshold_corrections * sqrt(768.0 / (float)it1st);
469  if (ibase == 32) threshold_corrections = 24.0;
470  if (ibase == 64) threshold_corrections = 12.0;
471  if (ibase == 128) threshold_corrections = 6.0;
472 
473  if (position1[l] > 0.01 * m_thr * threshold_corrections) { // if probablity is larger then threshold mark as Hot strip
474  hsflag[l] = 1; // HS vector
475  flag[l] = 0; // mark strip as bad for second pass
476  iths++;
477  B2RESULT("1st pass HS found! Layer: " << i + 3 << " Ladder: " << m << " Sensor: " << j << " Side: " << k << " channel: " << l);
478  } else {
479  hsflag[l] = 0; // not a HS
480  //recalculate the occupancy in DSSD only for good strip after first pass
481  occupancy[test.quot] = occupancy[test.quot] + hm_occupancy->getHistogram(*itSvdSensors, k)->GetBinContent(l + 1);
482  it++; //number of good strips after first pass
483  }
484 
485  }
486  /*
487  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
488  */
489  // second pass
490  for (int l = 0; l < 768; l++) {
491  div_t test = div(l, ibase);
492  position1[l] = position1[l] * nCltrk[test.quot] / (float)occupancy[test.quot];
493  float threshold_corrections = 1.0;
494  threshold_corrections = threshold_corrections * sqrt(768.0 / (float)it);
495  if (ibase == 32) threshold_corrections = 24.0;
496  if (ibase == 64) threshold_corrections = 12.0;
497  if (ibase == 128) threshold_corrections = 6.0;
498 
499  if ((flag[l]) && (position1[l] > 0.01 * m_thr * threshold_corrections)) { //HS
500  hsflag[l] = 1;// HS vector
501  flag[l] = 0; // mark strip as bad after second pass
502  iths++;
503  B2RESULT("2nd pass HS FOUND! Layer: " << i + 3 << " Ladder: " << m << " Sensor: " << j << " Side: " << k << " channel: " << l);
504  }
505  }
506  // for Laura .. HS flags, place interface for DB
507  // 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
508 
509 
510  for (int l = 0; l < 768; l++) {
511 
512  B2DEBUG(1, hsflag[l]);
513 
514  float tmpOcc = hm_occupancy->getHistogram(*itSvdSensors, k)->GetBinContent(l + 1) / (double)nevents;
515  hm_occAll->fill(*itSvdSensors, k, tmpOcc);
516 
517  if (hsflag[l] == 0) {
518  hm_occupancy_after->getHistogram(*itSvdSensors, k)->SetBinContent(l + 1, tmpOcc); //alive strips without identified HS
519  hm_occAfter->fill(*itSvdSensors, k, tmpOcc);
520  } else {
521  hm_hot_strips->getHistogram(*itSvdSensors, k)->SetBinContent(l + 1, 1);
522  hm_occHot->fill(*itSvdSensors, k, tmpOcc);
523  }
524  }
525 
526  if (m_rootFilePtr != nullptr) {
527  hm_occupancy->getHistogram(*itSvdSensors, k)->Scale(1.0 / (double)nevents);
528 
529  dir_occuL[i]->cd();
530  hm_occupancy->getHistogram(*itSvdSensors, k)->Write();
531  hm_hot_strips->getHistogram(*itSvdSensors, k)->SetLineColor(kBlack);
532  hm_hot_strips->getHistogram(*itSvdSensors, k)->SetFillStyle(3001);
533  hm_hot_strips->getHistogram(*itSvdSensors, k)->SetFillColor(kBlack);
534  hm_hot_strips->getHistogram(*itSvdSensors, k)->SetMarkerColor(kBlack);
535  hm_hot_strips->getHistogram(*itSvdSensors, k)->Write();
536  hm_occupancy_after->getHistogram(*itSvdSensors, k)->SetLineColor(kRed);
537  hm_occupancy_after->getHistogram(*itSvdSensors, k)->SetMarkerColor(kRed);
538  hm_occupancy_after->getHistogram(*itSvdSensors, k)->Write();
539  hm_occAll->getHistogram(*itSvdSensors, k)->Write();
540  hm_occHot->getHistogram(*itSvdSensors, k)->SetLineColor(kBlack);
541  hm_occHot->getHistogram(*itSvdSensors, k)->SetFillStyle(3001);
542  hm_occHot->getHistogram(*itSvdSensors, k)->SetFillColor(kBlack);
543  hm_occHot->getHistogram(*itSvdSensors, k)->SetMarkerColor(kBlack);
544  hm_occHot->getHistogram(*itSvdSensors, k)->Write();
545  hm_occAfter->getHistogram(*itSvdSensors, k)->SetLineColor(kRed);
546  hm_occAfter->getHistogram(*itSvdSensors, k)->SetMarkerColor(kRed);
547  hm_occAfter->getHistogram(*itSvdSensors, k)->Write();
548  // hm_dist12->getHistogram(*itSvdSensors, k)->Write();
549  }
550 
551  B2DEBUG(1, " side " << i << " " << j << " " << m << " " << k);
552  /* end */
553 
554 
555  // store number hot strips per sensor
556 
557 
558  for (int iy = 0; iy < iths; iy++) {
559  h_tot_dqm->Fill(float(itsensor));
560  h_tot_dqm1->Fill(float(itsensor));
561  }
562 
563  for (int s = 0; s < hm_hot_strips->getHistogram(*itSvdSensors, k)->GetEntries(); s++)
564  m_hHotStripsSummary->fill(*itSvdSensors, k, 1);
565 
566  itsensor++;
567  }
568  ++itSvdSensors;
569  }
570  ++itSvdLadders;
571  }
572  ++itSvdLayers;
573 
574  }
575 
576  if (m_rootFilePtr != nullptr) {
577  oldDir->cd();
578 
579  m_hHotStripsSummary->getHistogram(0)->Write();
580  m_hHotStripsSummary->getHistogram(1)->Write();
581 
582  TObject* obj;
583  TIter nextH_occu(m_histoList_occu);
584  while ((obj = nextH_occu()))
585  obj->Write();
586 
587  m_rootFilePtr->Close();
588  }
589 
590  }
591 }
592 
593 TH1F* SVDHotStripFinderModule::createHistogram1D(const char* name, const char* title,
594  Int_t nbins, Double_t min, Double_t max,
595  const char* xtitle, TList* histoList)
596 {
597 
598  TH1F* h = new TH1F(name, title, nbins, min, max);
599 
600  h->GetXaxis()->SetTitle(xtitle);
601 
602  if (histoList)
603  histoList->Add(h);
604 
605 
606  return h;
607 }
608 
609 
610 TH2F* SVDHotStripFinderModule::createHistogram2D(const char* name, const char* title,
611  Int_t nbinsX, Double_t minX, Double_t maxX,
612  const char* titleX,
613  Int_t nbinsY, Double_t minY, Double_t maxY,
614  const char* titleY, TList* histoList)
615 {
616 
617  TH2F* h = new TH2F(name, title, nbinsX, minX, maxX, nbinsY, minY, maxY);
618  h->GetXaxis()->SetTitle(titleX);
619  h->GetYaxis()->SetTitle(titleY);
620 
621  if (histoList)
622  histoList->Add(h);
623 
624  return h;
625 }
626 
627 
628 bool SVDHotStripFinderModule::theHSFinder(double* stripOccAfterAbsCut, int* hsflag, int nstrips)
629 {
630  bool found = false;
631 
632  if (m_base == -1)
633  m_base = nstrips;
634 
635  int N = nstrips / m_base;
636 
637  for (int sector = 0; sector < N; sector++) {
638 
639  int nafter = 0;
640  double sensorOccAverage = 0;
641 
642  for (int l = sector * m_base; l < sector * m_base + m_base; l++) {
643  sensorOccAverage = sensorOccAverage + stripOccAfterAbsCut[l];
644  if (stripOccAfterAbsCut[l] > 0) nafter++;
645  }
646  sensorOccAverage = sensorOccAverage / nafter;
647 
648  B2DEBUG(1, "Average occupancy: " << sensorOccAverage);
649 
650  for (int l = sector * m_base; l < sector * m_base + m_base; l++) {
651 
652  // flag additional HS by comparing each strip occupancy with the sensor-based average occupancy
653 
654  if (stripOccAfterAbsCut[l] > sensorOccAverage * m_relOccPrec) {
655  hsflag[l] = 1;
656  found = true;
657  stripOccAfterAbsCut[l] = 0;
658  }
659  // else hsflag[l]=0;
660  }
661  }
662 
663  return found;
664 }
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
Base class for Modules.
Definition: Module.h:72
Belle2::SVDSummaryPlots
class to summarize SVD quantities per sensor and side
Definition: SVDSummaryPlots.h:23
Belle2::VXD::GeoCache
Class to faciliate easy access to sensor information of the VXD like coordinate transformations or pi...
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:175
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:203
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:192
Belle2::VxdID
Class to uniquely identify a any structure of the PXD and SVD.
Definition: VxdID.h:33
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:650
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17