CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
SiStripBadAPVAlgorithmFromClusterOccupancy.cc
Go to the documentation of this file.
2 
3 
13 
14 
16  lowoccupancy_(0),
17  highoccupancy_(100),
18  absolutelow_(0),
19  numberiterations_(2),
20  Nevents_(0),
21  occupancy_(0),
22  OutFileName_("Occupancy.root"),
23  UseInputDB_(iConfig.getUntrackedParameter<bool>("UseInputDB",false))
24  {
26  }
27 
29  LogTrace("SiStripBadAPVAlgorithmFromClusterOccupancy")<<"[SiStripBadAPVAlgorithmFromClusterOccupancy::~SiStripBadAPVAlgorithmFromClusterOccupancy] "<<std::endl;
30 }
31 
33 
34  LogTrace("SiStripBadAPVAlgorithmFromClusterOccupancy")<<"[SiStripBadAPVAlgorithmFromClusterOccupancy::extractBadAPVs] "<<std::endl;
35 
36  if (WriteOutputFile_==true){
37  f = new TFile(OutFileName_.c_str(),"RECREATE");
38  f->cd();
39 
40  apvtree = new TTree("moduleOccupancy","tree");
41 
42  apvtree->Branch("DetRawId", &detrawid, "DetRawId/I");
43  apvtree->Branch("SubDetId", &subdetid, "SubDetId/I");
44  apvtree->Branch("Layer_Ring", &layer_ring, "Layer_Ring/I");
45  apvtree->Branch("Disc", &disc, "Disc/I");
46  apvtree->Branch("IsBack", &isback, "IsBack/I");
47  apvtree->Branch("IsExternalString", &isexternalstring, "IsExternalString/I");
48  apvtree->Branch("IsZMinusSide", &iszminusside, "IsZMinusSide/I");
49  apvtree->Branch("RodStringPetal", &rodstringpetal, "RodStringPetal/I");
50  apvtree->Branch("IsStereo", &isstereo, "IsStereo/I");
51  apvtree->Branch("ModuleNumber", &module_number, "ModuleNumber/I");
52  apvtree->Branch("NumberOfStrips", &number_strips, "NumberOfStrips/I");
53  apvtree->Branch("APVGlobalPositionX", &global_position_x, "APVGlobalPositionX/F");
54  apvtree->Branch("APVGlobalPositionY", &global_position_y, "APVGlobalPositionY/F");
55  apvtree->Branch("APVGlobalPositionZ", &global_position_z, "APVGlobalPositionZ/F");
56  apvtree->Branch("APVNumber", &apv_number, "APVNumber/I");
57  apvtree->Branch("APVAbsoluteOccupancy", &apvAbsoluteOccupancy, "apvAbsoluteOccupancy/D");
58  apvtree->Branch("APVMedianOccupancy", &apvMedianOccupancy, "apvMedianOccupancy/D");
59  }
60 
61  HistoMap::iterator it=DM.begin();
62  HistoMap::iterator itEnd=DM.end();
63  std::vector<unsigned int> badStripList;
64  uint32_t detid;
65  for (;it!=itEnd;++it){
66 
67  Apv APV;
68 
69  for (int apv=0; apv<6; apv++)
70  {
71  APV.apvMedian[apv] = 0;
72  APV.apvabsoluteOccupancy[apv] = 0;
73 
74  for (int strip=0; strip<128; strip++)
75  {
76  stripOccupancy[apv][strip] = 0;
77  stripWeight[apv][strip] = 0;
78  }
79  }
80 
81  pHisto phisto;
82  phisto._th1f = it->second.get();
83  phisto._NEntries = (int)phisto._th1f->GetEntries();
84  phisto._NBins = phisto._th1f->GetNbinsX();
85 
86  number_strips = (int)phisto._NBins;
88  APV.numberApvs = number_apvs;
89 
90  for (int apv=0; apv<number_apvs; apv++)
91  {
92  for (int strip=0; strip<128; strip++)
93  {
94  stripOccupancy[apv][strip] = phisto._th1f->GetBinContent((apv*128)+strip+1); // Remember: Bin=0 is underflow bin!
95  stripWeight[apv][strip] = 1;
96  APV.apvabsoluteOccupancy[apv] += phisto._th1f->GetBinContent((apv*128)+strip+1); // Remember: Bin=0 is underflow bin!
97  }
98  }
99 
100  for (int apv=0; apv<number_apvs; apv++)
101  {
102  APV.apvMedian[apv] = TMath::Median(128,stripOccupancy[apv],stripWeight[apv]);
103  }
104 
105  detid=it->first;
106  DetId detectorId=DetId(detid);
107 
108  if (edm::isDebugEnabled())
109  LogTrace("SiStripBadAPV") << "Analyzing detid " << detid<< std::endl;
110 
111  detrawid = detid;
112  APV.detrawId = detrawid;
113  subdetid = detectorId.subdetId();
114  if (SiStripDetId(detrawid).stereo() !=0 ) isstereo = 1; // It's a stereo module
115  else isstereo = 0; // It's an rphi module
116  switch (detectorId.subdetId())
117  {
118  case StripSubdetector::TIB :
120  disc = -1;
121  isback = -1;
122  if (TIBDetId(detrawid).isExternalString()) isexternalstring = 1;
123  else isexternalstring = 0;
124  if (TIBDetId(detrawid).isZMinusSide()) iszminusside = 1;
125  else iszminusside = 0;
129 
130  if (layer_ring == 1) medianValues_TIB_Layer1.push_back(APV);
131  else if (layer_ring == 2) medianValues_TIB_Layer2.push_back(APV);
132  else if (layer_ring == 3) medianValues_TIB_Layer3.push_back(APV);
133  else if (layer_ring == 4) medianValues_TIB_Layer4.push_back(APV);
134  break;
135 
136  case StripSubdetector::TID :
139  if (TIDDetId(detrawid).isBackRing()) isback = 1;
140  else isback = 0;
141  if (TIDDetId(detrawid).isZMinusSide()) iszminusside = 1;
142  else iszminusside = 0;
143  isexternalstring = -1;
144  rodstringpetal = -1;
147 
148  if (iszminusside==0)
149  {
150  if (disc==1) medianValues_TIDPlus_Disc1.push_back(APV);
151  else if (disc==2) medianValues_TIDPlus_Disc2.push_back(APV);
152  else if (disc==3) medianValues_TIDPlus_Disc3.push_back(APV);
153  }
154  else if (iszminusside==1)
155  {
156  if (disc==1) medianValues_TIDMinus_Disc1.push_back(APV);
157  else if (disc==2) medianValues_TIDMinus_Disc2.push_back(APV);
158  else if (disc==3) medianValues_TIDMinus_Disc3.push_back(APV);
159  }
160  break;
161 
162  case StripSubdetector::TOB :
164  disc = -1;
165  isback = -1;
166  if (TOBDetId(detrawid).isZMinusSide()) iszminusside = 1;
167  else iszminusside = 0;
168  isexternalstring = -1;
172 
173  if (layer_ring == 1) medianValues_TOB_Layer1.push_back(APV);
174  else if (layer_ring == 2) medianValues_TOB_Layer2.push_back(APV);
175  else if (layer_ring == 3) medianValues_TOB_Layer3.push_back(APV);
176  else if (layer_ring == 4) medianValues_TOB_Layer4.push_back(APV);
177  else if (layer_ring == 5) medianValues_TOB_Layer5.push_back(APV);
178  else if (layer_ring == 6) medianValues_TOB_Layer6.push_back(APV);
179  break;
180 
181  case StripSubdetector::TEC :
184  if (TECDetId(detrawid).isBackPetal()) isback = 1;
185  else isback = 0;
186  if (TECDetId(detrawid).isZMinusSide()) iszminusside = 1;
187  else iszminusside = 0;
188  isexternalstring = -1;
192 
193  if (iszminusside==0)
194  {
195  if (disc==1) medianValues_TECPlus_Disc1.push_back(APV);
196  else if (disc==2) medianValues_TECPlus_Disc2.push_back(APV);
197  else if (disc==3) medianValues_TECPlus_Disc3.push_back(APV);
198  else if (disc==4) medianValues_TECPlus_Disc4.push_back(APV);
199  else if (disc==5) medianValues_TECPlus_Disc5.push_back(APV);
200  else if (disc==6) medianValues_TECPlus_Disc6.push_back(APV);
201  else if (disc==7) medianValues_TECPlus_Disc7.push_back(APV);
202  else if (disc==8) medianValues_TECPlus_Disc8.push_back(APV);
203  else if (disc==9) medianValues_TECPlus_Disc9.push_back(APV);
204  }
205  else if (iszminusside==1)
206  {
207  if (disc==1) medianValues_TECMinus_Disc1.push_back(APV);
208  else if (disc==2) medianValues_TECMinus_Disc2.push_back(APV);
209  else if (disc==3) medianValues_TECMinus_Disc3.push_back(APV);
210  else if (disc==4) medianValues_TECMinus_Disc4.push_back(APV);
211  else if (disc==5) medianValues_TECMinus_Disc5.push_back(APV);
212  else if (disc==6) medianValues_TECMinus_Disc6.push_back(APV);
213  else if (disc==7) medianValues_TECMinus_Disc7.push_back(APV);
214  else if (disc==8) medianValues_TECMinus_Disc8.push_back(APV);
215  else if (disc==9) medianValues_TECMinus_Disc9.push_back(APV);
216  }
217  break;
218 
219  default :
220  std::cout << "### Detector does not belong to TIB, TID, TOB or TEC !? ###" << std::endl;
221  std::cout << "### DetRawId: " << detrawid << " ###" << std::endl;
222  }
223 
224  const StripGeomDetUnit* theStripDet = dynamic_cast<const StripGeomDetUnit*>( (TkGeom->idToDet(detectorId)) );
225  const StripTopology* theStripTopol = dynamic_cast<const StripTopology*>( &(theStripDet->specificTopology()) );
226 
227  for (int apv=0; apv<number_apvs; apv++)
228  {
229  apv_number = apv+1;
230  apvMedianOccupancy = APV.apvMedian[apv];
232 
233  LocalPoint pos_strip_local = theStripTopol->localPosition((apv*128));
234  GlobalPoint pos_strip_global = (TkGeom->idToDet(detectorId))->surface().toGlobal(pos_strip_local);
235 
236  global_position_x = pos_strip_global.x();
237  global_position_y = pos_strip_global.y();
238  global_position_z = pos_strip_global.z();
239 
240  if (WriteOutputFile_==true) apvtree->Fill();
241  }
242 
243  } // end loop on modules
244 
245  // Calculate Mean and RMS for each Layer
250 
257 
264 
274 
284 
285  pQuality=siStripQuality;
286  badStripList.clear();
287 
288  // Analyze the APV Occupancy for hot APVs
289  AnalyzeOccupancy(siStripQuality,medianValues_TIB_Layer1,MeanAndRms_TIB_Layer1,badStripList,inSiStripQuality);
290  AnalyzeOccupancy(siStripQuality,medianValues_TIB_Layer2,MeanAndRms_TIB_Layer2,badStripList,inSiStripQuality);
291  AnalyzeOccupancy(siStripQuality,medianValues_TIB_Layer3,MeanAndRms_TIB_Layer3,badStripList,inSiStripQuality);
292  AnalyzeOccupancy(siStripQuality,medianValues_TIB_Layer4,MeanAndRms_TIB_Layer4,badStripList,inSiStripQuality);
293 
294  AnalyzeOccupancy(siStripQuality,medianValues_TOB_Layer1,MeanAndRms_TOB_Layer1,badStripList,inSiStripQuality);
295  AnalyzeOccupancy(siStripQuality,medianValues_TOB_Layer2,MeanAndRms_TOB_Layer2,badStripList,inSiStripQuality);
296  AnalyzeOccupancy(siStripQuality,medianValues_TOB_Layer3,MeanAndRms_TOB_Layer3,badStripList,inSiStripQuality);
297  AnalyzeOccupancy(siStripQuality,medianValues_TOB_Layer4,MeanAndRms_TOB_Layer4,badStripList,inSiStripQuality);
298  AnalyzeOccupancy(siStripQuality,medianValues_TOB_Layer5,MeanAndRms_TOB_Layer5,badStripList,inSiStripQuality);
299  AnalyzeOccupancy(siStripQuality,medianValues_TOB_Layer6,MeanAndRms_TOB_Layer6,badStripList,inSiStripQuality);
300 
301  AnalyzeOccupancy(siStripQuality,medianValues_TIDPlus_Disc1,MeanAndRms_TIDPlus_Disc1,badStripList,inSiStripQuality);
302  AnalyzeOccupancy(siStripQuality,medianValues_TIDPlus_Disc2,MeanAndRms_TIDPlus_Disc2,badStripList,inSiStripQuality);
303  AnalyzeOccupancy(siStripQuality,medianValues_TIDPlus_Disc3,MeanAndRms_TIDPlus_Disc3,badStripList,inSiStripQuality);
304  AnalyzeOccupancy(siStripQuality,medianValues_TIDMinus_Disc1,MeanAndRms_TIDMinus_Disc1,badStripList,inSiStripQuality);
305  AnalyzeOccupancy(siStripQuality,medianValues_TIDMinus_Disc2,MeanAndRms_TIDMinus_Disc2,badStripList,inSiStripQuality);
306  AnalyzeOccupancy(siStripQuality,medianValues_TIDMinus_Disc3,MeanAndRms_TIDMinus_Disc3,badStripList,inSiStripQuality);
307 
308  AnalyzeOccupancy(siStripQuality,medianValues_TECPlus_Disc1,MeanAndRms_TECPlus_Disc1,badStripList,inSiStripQuality);
309  AnalyzeOccupancy(siStripQuality,medianValues_TECPlus_Disc2,MeanAndRms_TECPlus_Disc2,badStripList,inSiStripQuality);
310  AnalyzeOccupancy(siStripQuality,medianValues_TECPlus_Disc3,MeanAndRms_TECPlus_Disc3,badStripList,inSiStripQuality);
311  AnalyzeOccupancy(siStripQuality,medianValues_TECPlus_Disc4,MeanAndRms_TECPlus_Disc4,badStripList,inSiStripQuality);
312  AnalyzeOccupancy(siStripQuality,medianValues_TECPlus_Disc5,MeanAndRms_TECPlus_Disc5,badStripList,inSiStripQuality);
313  AnalyzeOccupancy(siStripQuality,medianValues_TECPlus_Disc6,MeanAndRms_TECPlus_Disc6,badStripList,inSiStripQuality);
314  AnalyzeOccupancy(siStripQuality,medianValues_TECPlus_Disc7,MeanAndRms_TECPlus_Disc7,badStripList,inSiStripQuality);
315  AnalyzeOccupancy(siStripQuality,medianValues_TECPlus_Disc8,MeanAndRms_TECPlus_Disc8,badStripList,inSiStripQuality);
316  AnalyzeOccupancy(siStripQuality,medianValues_TECPlus_Disc9,MeanAndRms_TECPlus_Disc9,badStripList,inSiStripQuality);
317 
318  AnalyzeOccupancy(siStripQuality,medianValues_TECMinus_Disc1,MeanAndRms_TECMinus_Disc1,badStripList,inSiStripQuality);
319  AnalyzeOccupancy(siStripQuality,medianValues_TECMinus_Disc2,MeanAndRms_TECMinus_Disc2,badStripList,inSiStripQuality);
320  AnalyzeOccupancy(siStripQuality,medianValues_TECMinus_Disc3,MeanAndRms_TECMinus_Disc3,badStripList,inSiStripQuality);
321  AnalyzeOccupancy(siStripQuality,medianValues_TECMinus_Disc4,MeanAndRms_TECMinus_Disc4,badStripList,inSiStripQuality);
322  AnalyzeOccupancy(siStripQuality,medianValues_TECMinus_Disc5,MeanAndRms_TECMinus_Disc5,badStripList,inSiStripQuality);
323  AnalyzeOccupancy(siStripQuality,medianValues_TECMinus_Disc6,MeanAndRms_TECMinus_Disc6,badStripList,inSiStripQuality);
324  AnalyzeOccupancy(siStripQuality,medianValues_TECMinus_Disc7,MeanAndRms_TECMinus_Disc7,badStripList,inSiStripQuality);
325  AnalyzeOccupancy(siStripQuality,medianValues_TECMinus_Disc8,MeanAndRms_TECMinus_Disc8,badStripList,inSiStripQuality);
326  AnalyzeOccupancy(siStripQuality,medianValues_TECMinus_Disc9,MeanAndRms_TECMinus_Disc9,badStripList,inSiStripQuality);
327 
328  siStripQuality->fillBadComponents();
329 
330  if (WriteOutputFile_==true){
331  f->cd();
332  apvtree->Write();
333  f->Close();
334  }
335 
336  LogTrace("SiStripBadAPV") << ss.str() << std::endl;
337 }
338 
339 
340 void SiStripBadAPVAlgorithmFromClusterOccupancy::CalculateMeanAndRMS(std::vector<Apv> a, std::pair<double,double>* MeanRMS, int number_iterations)
341 {
342  Double_t tot[7], tot2[7];
343  Double_t n[7];
344 
345  Double_t Mean[7] = {0};
346  Double_t Rms[7] = {1000,1000,1000,1000,1000,1000,1000};
347 
348  int Moduleposition;
349 
350  for (int i=0; i<number_iterations; i++)
351  {
352  for (int j=0; j<7; j++)
353  {
354  n[j] = 0;
355  tot[j] = 0;
356  tot2[j] = 0;
357  }
358 
359  for (uint32_t it=0; it<a.size(); it++)
360  {
361  Moduleposition = (a[it].modulePosition)-1;
362 
363  for (int apv=0; apv<a[it].numberApvs; apv++)
364  {
365  if (i>0)
366  {
367  if (a[it].apvMedian[apv]<(Mean[Moduleposition]-3*Rms[Moduleposition]) || (a[it].apvMedian[apv]>(Mean[Moduleposition]+5*Rms[Moduleposition])))
368  {
369  continue;
370  }
371  }
372  tot[Moduleposition] += a[it].apvMedian[apv];
373  tot2[Moduleposition] += (a[it].apvMedian[apv])*(a[it].apvMedian[apv]);
374  n[Moduleposition]++;
375  }
376  }
377 
378  for (int j=0; j<7; j++)
379  {
380  if (n[j]!=0)
381  {
382  Mean[j] = tot[j]/n[j];
383  Rms[j] = TMath::Sqrt(TMath::Abs(tot2[j]/n[j] -Mean[j]*Mean[j]));
384  }
385  }
386  }
387 
388  for (int j=0; j<7; j++)
389  {
390  MeanRMS[j] = std::make_pair(Mean[j],Rms[j]);
391  }
392 
393 }
394 
395 void SiStripBadAPVAlgorithmFromClusterOccupancy::AnalyzeOccupancy(SiStripQuality* quality, std::vector<Apv>& medianValues, std::pair<double,double>* MeanAndRms, std::vector<unsigned int>& BadStripList, edm::ESHandle<SiStripQuality>& InSiStripQuality)
396 {
397  int Moduleposition;
398  uint32_t Detid;
399 
400  for (uint32_t it=0; it<medianValues.size(); it++)
401  {
402  Moduleposition = (medianValues[it].modulePosition)-1;
403  Detid = medianValues[it].detrawId;
404 
405  for (int apv=0; apv<medianValues[it].numberApvs; apv++)
406  {
407  if(UseInputDB_)
408  {
409  if(InSiStripQuality->IsApvBad(Detid,apv) )
410  {
411  continue;//if the apv is already flagged as bad, continue.
412  }
413  }
414  if (medianValues[it].apvMedian[apv] > minNevents_)
415  {
416  if ((medianValues[it].apvMedian[apv]>(MeanAndRms[Moduleposition].first+highoccupancy_*MeanAndRms[Moduleposition].second)) && (medianValues[it].apvMedian[apv]>absolutelow_))
417  BadStripList.push_back(pQuality->encode((apv*128),128,0));
418  }
419  else if (medianValues[it].apvMedian[apv]<(MeanAndRms[Moduleposition].first-lowoccupancy_*MeanAndRms[Moduleposition].second) && (MeanAndRms[Moduleposition].first>2 || medianValues[it].apvabsoluteOccupancy[apv]==0))
420  {
421  BadStripList.push_back(pQuality->encode((apv*128),128,0));
422  std::cout << "Dead APV! DetId: " << medianValues[it].detrawId << ", APV number: " << apv+1 << ", APVMedian: " << medianValues[it].apvMedian[apv] << ", Mean: " << MeanAndRms[Moduleposition].first << ", RMS: " << MeanAndRms[Moduleposition].second << ", LowThreshold: " << lowoccupancy_ << ", Mean-Low*RMS: " << (MeanAndRms[Moduleposition].first-lowoccupancy_*MeanAndRms[Moduleposition].second) << std::endl;
423  }
424  }
425  if (BadStripList.begin()!=BadStripList.end())
426  {
427  quality->compact(Detid,BadStripList);
428  SiStripQuality::Range range(BadStripList.begin(),BadStripList.end());
429  quality->put(Detid,range);
430  }
431  BadStripList.clear();
432  }
433 }
434 
436 {
438 }
unsigned int rodNumber() const
Definition: TOBDetId.h:77
bool isDebugEnabled()
int i
Definition: DBlmapReader.cc:9
unsigned int petalNumber() const
Definition: TECDetId.h:94
unsigned int stringNumber() const
Definition: TIBDetId.h:87
unsigned int layer() const
layer id
Definition: TOBDetId.h:39
void strip(std::string &input, const std::string &blanks=" \n\t")
Definition: stringTools.cc:16
void AnalyzeOccupancy(SiStripQuality *, std::vector< Apv > &, std::pair< double, double > *, std::vector< unsigned int > &, edm::ESHandle< SiStripQuality > &)
T y() const
Definition: PV3DBase.h:62
virtual const StripTopology & specificTopology() const
Returns a reference to the strip proxy topology.
static const char occupancy_[]
U second(std::pair< T, U > const &p)
T z() const
Definition: PV3DBase.h:63
unsigned int moduleNumber() const
Definition: TIBDetId.h:91
void compact(unsigned int &, std::vector< unsigned int > &)
unsigned int ring() const
ring id
Definition: TIDDetId.h:55
int j
Definition: DBlmapReader.cc:9
bool first
Definition: L1TdeRCT.cc:94
void fillBadComponents()
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:39
virtual const GeomDet * idToDet(DetId) const
#define LogTrace(id)
void CalculateMeanAndRMS(std::vector< Apv >, std::pair< double, double > *, int)
unsigned int moduleNumber() const
Definition: TECDetId.h:102
Detector identifier class for the strip tracker.
Definition: SiStripDetId.h:17
Definition: DetId.h:20
unsigned int wheel() const
wheel id
Definition: TECDetId.h:52
unsigned int layer() const
layer id
Definition: TIBDetId.h:41
double a
Definition: hdecay.h:121
std::pair< ContainerIterator, ContainerIterator > Range
unsigned int ring() const
ring id
Definition: TECDetId.h:71
void extractBadAPVs(SiStripQuality *, HistoMap &, edm::ESHandle< SiStripQuality > &)
virtual LocalPoint localPosition(float strip) const =0
tuple cout
Definition: gather_cfg.py:121
bool put(const uint32_t &detID, const InputVector &vect)
T x() const
Definition: PV3DBase.h:61
unsigned int moduleNumber() const
Definition: TIDDetId.h:101
unsigned int encode(const unsigned short &first, const unsigned short &NconsecutiveBadStrips, const unsigned short &flag=0)
unsigned int moduleNumber() const
Definition: TOBDetId.h:81
unsigned int wheel() const
wheel id
Definition: TIDDetId.h:50