CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy.cc
Go to the documentation of this file.
1 #include "CalibTracker/SiStripQuality/interface/SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy.h"
2 
3 
4 #include "FWCore/MessageLogger/interface/MessageLogger.h"
5 #include "DataFormats/DetId/interface/DetId.h"
6 #include "DataFormats/SiStripDetId/interface/TIBDetId.h"
7 #include "DataFormats/SiStripDetId/interface/TIDDetId.h"
8 #include "DataFormats/SiStripDetId/interface/TOBDetId.h"
9 #include "DataFormats/SiStripDetId/interface/TECDetId.h"
10 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
11 #include "Geometry/CommonTopologies/interface/StripTopology.h"
12 #include "CalibFormats/SiStripObjects/interface/SiStripQuality.h"
13 
14 
15 SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy::SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy(const edm::ParameterSet& iConfig):
16  lowoccupancy_(0),
17  highoccupancy_(100),
18  absolutelow_(0),
19  numberiterations_(2),
20  Nevents_(0),
21  absolute_occupancy_(0),
22  OutFileName_("Occupancy.root"),
23  UseInputDB_(iConfig.getUntrackedParameter<bool>("UseInputDB",false))
24  {
25  minNevents_=Nevents_*absolute_occupancy_;
26  }
27 
28 SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy::~SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy(){
29  LogTrace("SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy")<<"[SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy::~SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy] "<<std::endl;
30 }
31 
32 void SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy::extractBadAPVSandStrips(SiStripQuality* siStripQuality,HistoMap& DM, edm::ESHandle<SiStripQuality>& inSiStripQuality){
33 
34  LogTrace("SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy")<<"[SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy::extractBadAPVs] "<<std::endl;
35 
36  if (WriteOutputFile_== true)
37  {
38  f = new TFile(OutFileName_.c_str(),"RECREATE");
39  f->cd();
40 
41  apvtree = new TTree("moduleOccupancy","tree");
42 
43  apvtree->Branch("DetRawId", &detrawid, "DetRawId/I");
44  apvtree->Branch("SubDetId", &subdetid, "SubDetId/I");
45  apvtree->Branch("Layer_Ring", &layer_ring, "Layer_Ring/I");
46  apvtree->Branch("Disc", &disc, "Disc/I");
47  apvtree->Branch("IsBack", &isback, "IsBack/I");
48  apvtree->Branch("IsExternalString", &isexternalstring, "IsExternalString/I");
49  apvtree->Branch("IsZMinusSide", &iszminusside, "IsZMinusSide/I");
50  apvtree->Branch("RodStringPetal", &rodstringpetal, "RodStringPetal/I");
51  apvtree->Branch("IsStereo", &isstereo, "IsStereo/I");
52  apvtree->Branch("ModuleNumber", &module_number, "ModuleNumber/I");
53  apvtree->Branch("NumberOfStrips", &number_strips, "NumberOfStrips/I");
54  apvtree->Branch("APVGlobalPositionX", &global_position_x, "APVGlobalPositionX/F");
55  apvtree->Branch("APVGlobalPositionY", &global_position_y, "APVGlobalPositionY/F");
56  apvtree->Branch("APVGlobalPositionZ", &global_position_z, "APVGlobalPositionZ/F");
57  apvtree->Branch("APVNumber", &apv_number, "APVNumber/I");
58  apvtree->Branch("APVAbsoluteOccupancy", &apvAbsoluteOccupancy, "apvAbsoluteOccupancy/I");
59  apvtree->Branch("APVMedianOccupancy", &apvMedianOccupancy, "apvMedianOccupancy/D");
60  apvtree->Branch("IsBad", &isBad, "IsBad/I");
61 
62  striptree = new TTree("stripOccupancy","tree");
63 
64  striptree->Branch("DetRawId", &detrawid, "DetRawId/I");
65  striptree->Branch("SubDetId", &subdetid, "SubDetId/I");
66  striptree->Branch("Layer_Ring", &layer_ring, "Layer_Ring/I");
67  striptree->Branch("Disc", &disc, "Disc/I");
68  striptree->Branch("IsBack", &isback, "IsBack/I");
69  striptree->Branch("IsExternalString", &isexternalstring, "IsExternalString/I");
70  striptree->Branch("IsZMinusSide", &iszminusside, "IsZMinusSide/I");
71  striptree->Branch("RodStringPetal", &rodstringpetal, "RodStringPetal/I");
72  striptree->Branch("IsStereo", &isstereo, "IsStereo/I");
73  striptree->Branch("ModulePosition", &module_number, "ModulePosition/I");
74  striptree->Branch("NumberOfStrips", &number_strips, "NumberOfStrips/I");
75  striptree->Branch("StripNumber", &strip_number, "StripNumber/I");
76  striptree->Branch("APVChannel", &apv_channel, "APVChannel/I");
77  striptree->Branch("StripGlobalPositionX", &strip_global_position_x, "StripGlobalPositionX/F");
78  striptree->Branch("StripGlobalPositionY", &strip_global_position_y, "StripGlobalPositionY/F");
79  striptree->Branch("StripGlobalPositionZ", &strip_global_position_z, "StripGlobalPositionZ/F");
80  striptree->Branch("IsHot", &isHot, "IsHot/I");
81  striptree->Branch("HotStripsPerAPV", &hotStripsPerAPV, "HotStripsPerAPV/I");
82  striptree->Branch("HotStripsPerModule", &hotStripsPerModule,"HotStripsPerModule/I");
83  striptree->Branch("StripOccupancy", &singleStripOccupancy, "StripOccupancy/D");
84  striptree->Branch("StripHits", &stripHits, "StripHits/I");
85  striptree->Branch("PoissonProb", &poissonProb, "PoissonProb/D");
86  }
87 
88  HistoMap::iterator it=DM.begin();
89  HistoMap::iterator itEnd=DM.end();
90  std::vector<unsigned int> badStripList;
91  uint32_t detid;
92  for (;it!=itEnd;++it){
93 
94  Apv APV;
95 
96  for (int apv=0; apv<6; apv++)
97  {
98  APV.apvMedian[apv] = 0;
99  APV.apvabsoluteOccupancy[apv] = 0;
100  APV.NEntries[apv] = 0;
101  APV.NEmptyBins[apv] = 0;
102 
103  for (int strip=0; strip<128; strip++)
104  {
105  stripOccupancy[apv][strip] = 0;
106  stripWeight[apv][strip] = 0;
107  }
108  }
109 
110  number_strips = (int)((it->second.get())->GetNbinsX());
111  number_apvs = number_strips/128;
112  APV.numberApvs = number_apvs;
113 
114  for (int apv=0; apv<number_apvs; apv++)
115  {
116  APV.th1f[apv] = new TH1F("tmp","tmp",128,0.5,128.5);
117  int NumberEntriesPerAPV=0;
118 
119  for (int strip=0; strip<128; strip++)
120  {
121  stripOccupancy[apv][strip] = (it->second.get())->GetBinContent((apv*128)+strip+1); // Remember: Bin=0 is underflow bin!
122  stripWeight[apv][strip] = 1;
123  APV.apvabsoluteOccupancy[apv] += (it->second.get())->GetBinContent((apv*128)+strip+1); // Remember: Bin=0 is underflow bin!
124  APV.th1f[apv]->SetBinContent(strip+1,(it->second.get())->GetBinContent((apv*128)+strip+1));
125  NumberEntriesPerAPV += (int)(it->second.get())->GetBinContent((apv*128)+strip+1);
126  }
127 
128  APV.th1f[apv]->SetEntries(NumberEntriesPerAPV);
129  APV.NEntries[apv]=(int)APV.th1f[apv]->GetEntries();
130  }
131 
132  for (int apv=0; apv<number_apvs; apv++)
133  {
134  APV.apvMedian[apv] = TMath::Median(128,stripOccupancy[apv],stripWeight[apv]);
135  }
136 
137  detid=it->first;
138  DetId detectorId=DetId(detid);
139 
140  if (edm::isDebugEnabled())
141  LogTrace("SiStripBadAPV") << "Analyzing detid " << detid<< std::endl;
142 
143  detrawid = detid;
144  APV.detrawId = detrawid;
145  subdetid = detectorId.subdetId();
146 
147  switch (detectorId.subdetId())
148  {
149  case StripSubdetector::TIB :
150  layer_ring = TIBDetId(detrawid).layer();
151  module_number = TIBDetId(detrawid).moduleNumber();
152  APV.modulePosition = module_number;
153 
154  if (layer_ring == 1) medianValues_TIB_Layer1.push_back(APV);
155  else if (layer_ring == 2) medianValues_TIB_Layer2.push_back(APV);
156  else if (layer_ring == 3) medianValues_TIB_Layer3.push_back(APV);
157  else if (layer_ring == 4) medianValues_TIB_Layer4.push_back(APV);
158  break;
159 
160  case StripSubdetector::TID :
161  layer_ring = TIDDetId(detrawid).ring();
162  disc = TIDDetId(detrawid).wheel();
163  APV.modulePosition = layer_ring;
164 
165  if (TIDDetId(detrawid).isZMinusSide()) iszminusside = 1;
166  else iszminusside = 0;
167 
168  if (iszminusside==0)
169  {
170  if (disc==1) medianValues_TIDPlus_Disc1.push_back(APV);
171  else if (disc==2) medianValues_TIDPlus_Disc2.push_back(APV);
172  else if (disc==3) medianValues_TIDPlus_Disc3.push_back(APV);
173  }
174  else if (iszminusside==1)
175  {
176  if (disc==1) medianValues_TIDMinus_Disc1.push_back(APV);
177  else if (disc==2) medianValues_TIDMinus_Disc2.push_back(APV);
178  else if (disc==3) medianValues_TIDMinus_Disc3.push_back(APV);
179  }
180  break;
181 
182  case StripSubdetector::TOB :
183  layer_ring = TOBDetId(detrawid).layer();
184  module_number = TOBDetId(detrawid).moduleNumber();
185  APV.modulePosition = module_number;
186 
187  if (layer_ring == 1) medianValues_TOB_Layer1.push_back(APV);
188  else if (layer_ring == 2) medianValues_TOB_Layer2.push_back(APV);
189  else if (layer_ring == 3) medianValues_TOB_Layer3.push_back(APV);
190  else if (layer_ring == 4) medianValues_TOB_Layer4.push_back(APV);
191  else if (layer_ring == 5) medianValues_TOB_Layer5.push_back(APV);
192  else if (layer_ring == 6) medianValues_TOB_Layer6.push_back(APV);
193  break;
194 
195  case StripSubdetector::TEC :
196  layer_ring = TECDetId(detrawid).ring();
197  disc = TECDetId(detrawid).wheel();
198  APV.modulePosition = layer_ring;
199 
200  if (TECDetId(detrawid).isZMinusSide()) iszminusside = 1;
201  else iszminusside = 0;
202 
203  if (iszminusside==0)
204  {
205  if (disc==1) medianValues_TECPlus_Disc1.push_back(APV);
206  else if (disc==2) medianValues_TECPlus_Disc2.push_back(APV);
207  else if (disc==3) medianValues_TECPlus_Disc3.push_back(APV);
208  else if (disc==4) medianValues_TECPlus_Disc4.push_back(APV);
209  else if (disc==5) medianValues_TECPlus_Disc5.push_back(APV);
210  else if (disc==6) medianValues_TECPlus_Disc6.push_back(APV);
211  else if (disc==7) medianValues_TECPlus_Disc7.push_back(APV);
212  else if (disc==8) medianValues_TECPlus_Disc8.push_back(APV);
213  else if (disc==9) medianValues_TECPlus_Disc9.push_back(APV);
214  }
215  else if (iszminusside==1)
216  {
217  if (disc==1) medianValues_TECMinus_Disc1.push_back(APV);
218  else if (disc==2) medianValues_TECMinus_Disc2.push_back(APV);
219  else if (disc==3) medianValues_TECMinus_Disc3.push_back(APV);
220  else if (disc==4) medianValues_TECMinus_Disc4.push_back(APV);
221  else if (disc==5) medianValues_TECMinus_Disc5.push_back(APV);
222  else if (disc==6) medianValues_TECMinus_Disc6.push_back(APV);
223  else if (disc==7) medianValues_TECMinus_Disc7.push_back(APV);
224  else if (disc==8) medianValues_TECMinus_Disc8.push_back(APV);
225  else if (disc==9) medianValues_TECMinus_Disc9.push_back(APV);
226  }
227  break;
228 
229  default :
230  std::cout << "### Detector does not belong to TIB, TID, TOB or TEC !? ###" << std::endl;
231  std::cout << "### DetRawId: " << detrawid << " ###" << std::endl;
232  }
233 
234  } // end loop on modules
235 
236  // Calculate Mean and RMS for each Layer
237  CalculateMeanAndRMS(medianValues_TIB_Layer1,MeanAndRms_TIB_Layer1,numberiterations_);
238  CalculateMeanAndRMS(medianValues_TIB_Layer2,MeanAndRms_TIB_Layer2,numberiterations_);
239  CalculateMeanAndRMS(medianValues_TIB_Layer3,MeanAndRms_TIB_Layer3,numberiterations_);
240  CalculateMeanAndRMS(medianValues_TIB_Layer4,MeanAndRms_TIB_Layer4,numberiterations_);
241 
242  CalculateMeanAndRMS(medianValues_TOB_Layer1,MeanAndRms_TOB_Layer1,numberiterations_);
243  CalculateMeanAndRMS(medianValues_TOB_Layer2,MeanAndRms_TOB_Layer2,numberiterations_);
244  CalculateMeanAndRMS(medianValues_TOB_Layer3,MeanAndRms_TOB_Layer3,numberiterations_);
245  CalculateMeanAndRMS(medianValues_TOB_Layer4,MeanAndRms_TOB_Layer4,numberiterations_);
246  CalculateMeanAndRMS(medianValues_TOB_Layer5,MeanAndRms_TOB_Layer5,numberiterations_);
247  CalculateMeanAndRMS(medianValues_TOB_Layer6,MeanAndRms_TOB_Layer6,numberiterations_);
248 
249  CalculateMeanAndRMS(medianValues_TIDPlus_Disc1,MeanAndRms_TIDPlus_Disc1,numberiterations_);
250  CalculateMeanAndRMS(medianValues_TIDPlus_Disc2,MeanAndRms_TIDPlus_Disc2,numberiterations_);
251  CalculateMeanAndRMS(medianValues_TIDPlus_Disc3,MeanAndRms_TIDPlus_Disc3,numberiterations_);
252  CalculateMeanAndRMS(medianValues_TIDMinus_Disc1,MeanAndRms_TIDMinus_Disc1,numberiterations_);
253  CalculateMeanAndRMS(medianValues_TIDMinus_Disc2,MeanAndRms_TIDMinus_Disc2,numberiterations_);
254  CalculateMeanAndRMS(medianValues_TIDMinus_Disc3,MeanAndRms_TIDMinus_Disc3,numberiterations_);
255 
256  CalculateMeanAndRMS(medianValues_TECPlus_Disc1,MeanAndRms_TECPlus_Disc1,numberiterations_);
257  CalculateMeanAndRMS(medianValues_TECPlus_Disc2,MeanAndRms_TECPlus_Disc2,numberiterations_);
258  CalculateMeanAndRMS(medianValues_TECPlus_Disc3,MeanAndRms_TECPlus_Disc3,numberiterations_);
259  CalculateMeanAndRMS(medianValues_TECPlus_Disc4,MeanAndRms_TECPlus_Disc4,numberiterations_);
260  CalculateMeanAndRMS(medianValues_TECPlus_Disc5,MeanAndRms_TECPlus_Disc5,numberiterations_);
261  CalculateMeanAndRMS(medianValues_TECPlus_Disc6,MeanAndRms_TECPlus_Disc6,numberiterations_);
262  CalculateMeanAndRMS(medianValues_TECPlus_Disc7,MeanAndRms_TECPlus_Disc7,numberiterations_);
263  CalculateMeanAndRMS(medianValues_TECPlus_Disc8,MeanAndRms_TECPlus_Disc8,numberiterations_);
264  CalculateMeanAndRMS(medianValues_TECPlus_Disc9,MeanAndRms_TECPlus_Disc9,numberiterations_);
265 
266  CalculateMeanAndRMS(medianValues_TECMinus_Disc1,MeanAndRms_TECMinus_Disc1,numberiterations_);
267  CalculateMeanAndRMS(medianValues_TECMinus_Disc2,MeanAndRms_TECMinus_Disc2,numberiterations_);
268  CalculateMeanAndRMS(medianValues_TECMinus_Disc3,MeanAndRms_TECMinus_Disc3,numberiterations_);
269  CalculateMeanAndRMS(medianValues_TECMinus_Disc4,MeanAndRms_TECMinus_Disc4,numberiterations_);
270  CalculateMeanAndRMS(medianValues_TECMinus_Disc5,MeanAndRms_TECMinus_Disc5,numberiterations_);
271  CalculateMeanAndRMS(medianValues_TECMinus_Disc6,MeanAndRms_TECMinus_Disc6,numberiterations_);
272  CalculateMeanAndRMS(medianValues_TECMinus_Disc7,MeanAndRms_TECMinus_Disc7,numberiterations_);
273  CalculateMeanAndRMS(medianValues_TECMinus_Disc8,MeanAndRms_TECMinus_Disc8,numberiterations_);
274  CalculateMeanAndRMS(medianValues_TECMinus_Disc9,MeanAndRms_TECMinus_Disc9,numberiterations_);
275 
276  pQuality=siStripQuality;
277  badStripList.clear();
278 
279  // Analyze the Occupancy for both APVs and Strips
280  AnalyzeOccupancy(siStripQuality,medianValues_TIB_Layer1,MeanAndRms_TIB_Layer1,badStripList,inSiStripQuality);
281  AnalyzeOccupancy(siStripQuality,medianValues_TIB_Layer2,MeanAndRms_TIB_Layer2,badStripList,inSiStripQuality);
282  AnalyzeOccupancy(siStripQuality,medianValues_TIB_Layer3,MeanAndRms_TIB_Layer3,badStripList,inSiStripQuality);
283  AnalyzeOccupancy(siStripQuality,medianValues_TIB_Layer4,MeanAndRms_TIB_Layer4,badStripList,inSiStripQuality);
284 
285  AnalyzeOccupancy(siStripQuality,medianValues_TOB_Layer1,MeanAndRms_TOB_Layer1,badStripList,inSiStripQuality);
286  AnalyzeOccupancy(siStripQuality,medianValues_TOB_Layer2,MeanAndRms_TOB_Layer2,badStripList,inSiStripQuality);
287  AnalyzeOccupancy(siStripQuality,medianValues_TOB_Layer3,MeanAndRms_TOB_Layer3,badStripList,inSiStripQuality);
288  AnalyzeOccupancy(siStripQuality,medianValues_TOB_Layer4,MeanAndRms_TOB_Layer4,badStripList,inSiStripQuality);
289  AnalyzeOccupancy(siStripQuality,medianValues_TOB_Layer5,MeanAndRms_TOB_Layer5,badStripList,inSiStripQuality);
290  AnalyzeOccupancy(siStripQuality,medianValues_TOB_Layer6,MeanAndRms_TOB_Layer6,badStripList,inSiStripQuality);
291 
292  AnalyzeOccupancy(siStripQuality,medianValues_TIDPlus_Disc1,MeanAndRms_TIDPlus_Disc1,badStripList,inSiStripQuality);
293  AnalyzeOccupancy(siStripQuality,medianValues_TIDPlus_Disc2,MeanAndRms_TIDPlus_Disc2,badStripList,inSiStripQuality);
294  AnalyzeOccupancy(siStripQuality,medianValues_TIDPlus_Disc3,MeanAndRms_TIDPlus_Disc3,badStripList,inSiStripQuality);
295  AnalyzeOccupancy(siStripQuality,medianValues_TIDMinus_Disc1,MeanAndRms_TIDMinus_Disc1,badStripList,inSiStripQuality);
296  AnalyzeOccupancy(siStripQuality,medianValues_TIDMinus_Disc2,MeanAndRms_TIDMinus_Disc2,badStripList,inSiStripQuality);
297  AnalyzeOccupancy(siStripQuality,medianValues_TIDMinus_Disc3,MeanAndRms_TIDMinus_Disc3,badStripList,inSiStripQuality);
298 
299  AnalyzeOccupancy(siStripQuality,medianValues_TECPlus_Disc1,MeanAndRms_TECPlus_Disc1,badStripList,inSiStripQuality);
300  AnalyzeOccupancy(siStripQuality,medianValues_TECPlus_Disc2,MeanAndRms_TECPlus_Disc2,badStripList,inSiStripQuality);
301  AnalyzeOccupancy(siStripQuality,medianValues_TECPlus_Disc3,MeanAndRms_TECPlus_Disc3,badStripList,inSiStripQuality);
302  AnalyzeOccupancy(siStripQuality,medianValues_TECPlus_Disc4,MeanAndRms_TECPlus_Disc4,badStripList,inSiStripQuality);
303  AnalyzeOccupancy(siStripQuality,medianValues_TECPlus_Disc5,MeanAndRms_TECPlus_Disc5,badStripList,inSiStripQuality);
304  AnalyzeOccupancy(siStripQuality,medianValues_TECPlus_Disc6,MeanAndRms_TECPlus_Disc6,badStripList,inSiStripQuality);
305  AnalyzeOccupancy(siStripQuality,medianValues_TECPlus_Disc7,MeanAndRms_TECPlus_Disc7,badStripList,inSiStripQuality);
306  AnalyzeOccupancy(siStripQuality,medianValues_TECPlus_Disc8,MeanAndRms_TECPlus_Disc8,badStripList,inSiStripQuality);
307  AnalyzeOccupancy(siStripQuality,medianValues_TECPlus_Disc9,MeanAndRms_TECPlus_Disc9,badStripList,inSiStripQuality);
308 
309  AnalyzeOccupancy(siStripQuality,medianValues_TECMinus_Disc1,MeanAndRms_TECMinus_Disc1,badStripList,inSiStripQuality);
310  AnalyzeOccupancy(siStripQuality,medianValues_TECMinus_Disc2,MeanAndRms_TECMinus_Disc2,badStripList,inSiStripQuality);
311  AnalyzeOccupancy(siStripQuality,medianValues_TECMinus_Disc3,MeanAndRms_TECMinus_Disc3,badStripList,inSiStripQuality);
312  AnalyzeOccupancy(siStripQuality,medianValues_TECMinus_Disc4,MeanAndRms_TECMinus_Disc4,badStripList,inSiStripQuality);
313  AnalyzeOccupancy(siStripQuality,medianValues_TECMinus_Disc5,MeanAndRms_TECMinus_Disc5,badStripList,inSiStripQuality);
314  AnalyzeOccupancy(siStripQuality,medianValues_TECMinus_Disc6,MeanAndRms_TECMinus_Disc6,badStripList,inSiStripQuality);
315  AnalyzeOccupancy(siStripQuality,medianValues_TECMinus_Disc7,MeanAndRms_TECMinus_Disc7,badStripList,inSiStripQuality);
316  AnalyzeOccupancy(siStripQuality,medianValues_TECMinus_Disc8,MeanAndRms_TECMinus_Disc8,badStripList,inSiStripQuality);
317  AnalyzeOccupancy(siStripQuality,medianValues_TECMinus_Disc9,MeanAndRms_TECMinus_Disc9,badStripList,inSiStripQuality);
318 
319  siStripQuality->fillBadComponents();
320 
321  if (WriteOutputFile_==true){
322  f->cd();
323  apvtree->Write();
324  striptree->Write();
325  f->Close();
326  }
327 
328  LogTrace("SiStripBadAPV") << ss.str() << std::endl;
329 }
330 
331 
332 void SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy::CalculateMeanAndRMS(std::vector<Apv> a, std::pair<double,double>* MeanRMS, int number_iterations)
333 {
334  Double_t tot[7], tot2[7];
335  Double_t n[7];
336 
337  Double_t Mean[7] = {0};
338  Double_t Rms[7] = {1000,1000,1000,1000,1000,1000,1000};
339 
340  int Moduleposition;
341 
342  for (int i=0; i<number_iterations; i++)
343  {
344  for (int j=0; j<7; j++)
345  {
346  n[j] = 0;
347  tot[j] = 0;
348  tot2[j] = 0;
349  }
350 
351  for (uint32_t it=0; it<a.size(); it++)
352  {
353  Moduleposition = (a[it].modulePosition)-1;
354 
355  for (int apv=0; apv<a[it].numberApvs; apv++)
356  {
357  if (i>0)
358  {
359  if (a[it].apvMedian[apv]<(Mean[Moduleposition]-3*Rms[Moduleposition]) || (a[it].apvMedian[apv]>(Mean[Moduleposition]+5*Rms[Moduleposition])))
360  {
361  continue;
362  }
363  }
364  tot[Moduleposition] += a[it].apvMedian[apv];
365  tot2[Moduleposition] += (a[it].apvMedian[apv])*(a[it].apvMedian[apv]);
366  n[Moduleposition]++;
367  }
368  }
369 
370  for (int j=0; j<7; j++)
371  {
372  if (n[j]!=0)
373  {
374  Mean[j] = tot[j]/n[j];
375  Rms[j] = TMath::Sqrt(TMath::Abs(tot2[j]/n[j] -Mean[j]*Mean[j]));
376  }
377  }
378  }
379 
380  for (int j=0; j<7; j++)
381  {
382  MeanRMS[j] = std::make_pair(Mean[j],Rms[j]);
383  }
384 
385 }
386 
387 void SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy::AnalyzeOccupancy(SiStripQuality* quality, std::vector<Apv>& medianValues, std::pair<double,double>* MeanAndRms, std::vector<unsigned int>& BadStripList, edm::ESHandle<SiStripQuality>& InSiStripQuality)
388 {
389  int Moduleposition;
390  uint32_t Detid;
391 
392  for (uint32_t it=0; it<medianValues.size(); it++)
393  {
394  Moduleposition = (medianValues[it].modulePosition)-1;
395  Detid = medianValues[it].detrawId;
396 
397  setBasicTreeParameters(Detid);
398 
399  DetId DetectorId=DetId(Detid);
400  const StripGeomDetUnit* TheStripDet = dynamic_cast<const StripGeomDetUnit*>( (TkGeom->idToDet(DetectorId)) );
401  const StripTopology* TheStripTopol = dynamic_cast<const StripTopology*>( &(TheStripDet->specificTopology()) );
402 
403  //Analyze the occupancies
404  hotstripspermodule = 0;
405 
406  for (int apv=0; apv<medianValues[it].numberApvs; apv++)
407  {
408 
409  for (int i=0; i<128; i++)
410  {
411  ishot[i] = 0;
412  stripoccupancy[i] = 0;
413  striphits[i] = 0;
414  poissonprob[i] = 0;
415  }
416 
417  apv_number = apv+1;
418  apvMedianOccupancy = medianValues[it].apvMedian[apv];
419  apvAbsoluteOccupancy = medianValues[it].apvabsoluteOccupancy[apv];
420  isBad = 0;
421  hotstripsperapv[apv] = 0;
422 
423  LocalPoint pos_apv_local = TheStripTopol->localPosition((apv*128));
424  GlobalPoint pos_apv_global = (TkGeom->idToDet(DetectorId))->surface().toGlobal(pos_apv_local);
425 
426  global_position_x = pos_apv_global.x();
427  global_position_y = pos_apv_global.y();
428  global_position_z = pos_apv_global.z();
429 
430  if(UseInputDB_)
431  {
432  if(InSiStripQuality->IsApvBad(Detid,apv) )
433  {
434  if (WriteOutputFile_==true)
435  {
436  apvtree->Fill();
437  for (int strip=0; strip<128; strip++)
438  {
439  strip_number = (apv*128)+strip+1;
440  apv_channel = apv+1;
441  isHot = ishot[strip];
442  singleStripOccupancy = stripoccupancy[strip];
443  stripHits = striphits[strip];
444  poissonProb = poissonprob[strip];
445 
446  hotStripsPerModule = hotstripspermodule;
447  hotStripsPerAPV = hotstripsperapv[apv];
448 
449  LocalPoint pos_strip_local = TheStripTopol->localPosition(strip);
450  GlobalPoint pos_strip_global = (TkGeom->idToDet(DetectorId))->surface().toGlobal(pos_strip_local);
451 
452  strip_global_position_x = pos_strip_global.x();
453  strip_global_position_y = pos_strip_global.y();
454  strip_global_position_z = pos_strip_global.z();
455  striptree->Fill();
456  }
457  }
458  continue;//if the apv is already flagged as bad, continue.
459  }
460  }
461  if (medianValues[it].apvMedian[apv] > minNevents_)
462  {
463  if ((medianValues[it].apvMedian[apv]>(MeanAndRms[Moduleposition].first+highoccupancy_*MeanAndRms[Moduleposition].second)) && (medianValues[it].apvMedian[apv]>absolutelow_))
464  {
465  BadStripList.push_back(pQuality->encode((apv*128),128,0));
466  isBad = 1;
467  }
468  }
469  else if (medianValues[it].apvMedian[apv]<(MeanAndRms[Moduleposition].first-lowoccupancy_*MeanAndRms[Moduleposition].second) && (MeanAndRms[Moduleposition].first>2 || medianValues[it].apvabsoluteOccupancy[apv]==0))
470  {
471  BadStripList.push_back(pQuality->encode((apv*128),128,0));
472  isBad = 1;
473  }
474 
475  if (isBad!=1)
476  {
477  iterativeSearch(medianValues[it],BadStripList,apv);
478  }
479 
480  if (WriteOutputFile_==true)
481  {
482  apvtree->Fill();
483  for (int strip=0; strip<128; strip++)
484  {
485  strip_number = (apv*128)+strip+1;
486  apv_channel = apv+1;
487  isHot = ishot[strip];
488  singleStripOccupancy = stripoccupancy[strip];
489  stripHits = striphits[strip];
490  poissonProb = poissonprob[strip];
491 
492  hotStripsPerModule = hotstripspermodule;
493  hotStripsPerAPV = hotstripsperapv[apv];
494 
495  LocalPoint pos_strip_local = TheStripTopol->localPosition(strip);
496  GlobalPoint pos_strip_global = (TkGeom->idToDet(DetectorId))->surface().toGlobal(pos_strip_local);
497 
498  strip_global_position_x = pos_strip_global.x();
499  strip_global_position_y = pos_strip_global.y();
500  strip_global_position_z = pos_strip_global.z();
501  striptree->Fill();
502  }
503  }
504  }
505 
506  if (BadStripList.begin()!=BadStripList.end())
507  {
508  quality->compact(Detid,BadStripList);
509  SiStripQuality::Range range(BadStripList.begin(),BadStripList.end());
510  quality->put(Detid,range);
511  }
512  BadStripList.clear();
513  }
514 }
515 
516 void SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy::iterativeSearch(Apv& histo,std::vector<unsigned int>& vect, int apv){
517  if (!histo.NEntries[apv] || histo.NEntries[apv] <=MinNumEntries_ || histo.NEntries[apv] <= minNevents_)
518  return;
519 
520  size_t startingSize=vect.size();
521  long double diff=1.-prob_;
522 
523  int Nbins = histo.th1f[apv]->GetNbinsX();
524  int ibinStart = 1;
525  int ibinStop = Nbins+1;
526  int MaxEntry = (int)histo.th1f[apv]->GetMaximum();
527 
528  std::vector<long double> vPoissonProbs(MaxEntry+1,0);
529  long double meanVal=1.*histo.NEntries[apv]/(1.*Nbins-histo.NEmptyBins[apv]);
530  evaluatePoissonian(vPoissonProbs,meanVal);
531 
532  for (Int_t i=ibinStart; i<ibinStop; ++i){
533  unsigned int entries= (unsigned int)histo.th1f[apv]->GetBinContent(i);
534 
535  if (ishot[i-1]==0){
536  stripoccupancy[i-1] = entries/(double) Nevents_;
537  striphits[i-1] = entries;
538  poissonprob[i-1] = 1-vPoissonProbs[entries];
539  }
540 
541  if (entries<=MinNumEntriesPerStrip_ || entries <= minNevents_)
542  continue;
543 
544  if(diff<vPoissonProbs[entries]){
545  ishot[i-1] = 1;
546  hotstripspermodule++;
547  hotstripsperapv[apv]++;
548  histo.th1f[apv]->SetBinContent(i,0.);
549  histo.NEntries[apv]-=entries;
550  histo.NEmptyBins[apv]++;
551  if (edm::isDebugEnabled())
552  LogTrace("SiStripHotStrip")<< " rejecting strip " << (apv*128)+i-1 << " value " << entries << " diff " << diff << " prob " << vPoissonProbs[entries]<< std::endl;
553  vect.push_back(pQuality->encode((apv*128)+i-1,1,0));
554  }
555 
556  }
557  if (edm::isDebugEnabled())
558  LogTrace("SiStripHotStrip") << " [SiStripHotStripAlgorithmFromClusterOccupancy::iterativeSearch] Nbins="<< Nbins << " MaxEntry="<<MaxEntry << " meanVal=" << meanVal << " NEmptyBins="<<histo.NEmptyBins[apv]<< " NEntries=" << histo.NEntries[apv] << " thEntries " << histo.th1f[apv]->GetEntries()<< " startingSize " << startingSize << " vector.size " << vect.size() << std::endl;
559 
560  if (vect.size()!=startingSize)
561  iterativeSearch(histo,vect,apv);
562 }
563 
564 void SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy::evaluatePoissonian(std::vector<long double>& vPoissonProbs, long double& meanVal){
565  for(size_t i=0;i<vPoissonProbs.size();++i){
566  vPoissonProbs[i]= (i==0)?TMath::Poisson(i,meanVal):vPoissonProbs[i-1]+TMath::Poisson(i,meanVal);
567  }
568 }
569 
570 void SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy::setBasicTreeParameters(int detid){
571  DetId DetectorID=DetId(detid);
572 
573  if (SiStripDetId(detid).stereo() !=0 ) isstereo = 1; // It's a stereo module
574  else isstereo = 0; // It's an rphi module
575  switch (DetectorID.subdetId())
576  {
577  case StripSubdetector::TIB :
578  layer_ring = TIBDetId(detid).layer();
579  disc = -1;
580  isback = -1;
581  if (TIBDetId(detid).isExternalString()) isexternalstring = 1;
582  else isexternalstring = 0;
583  if (TIBDetId(detid).isZMinusSide()) iszminusside = 1;
584  else iszminusside = 0;
585  rodstringpetal = TIBDetId(detid).stringNumber();
586  module_number = TIBDetId(detid).moduleNumber();
587 
588  break;
589 
590  case StripSubdetector::TID :
591  layer_ring = TIDDetId(detid).ring();
592  disc = TIDDetId(detid).wheel();
593  if (TIDDetId(detid).isBackRing()) isback = 1;
594  else isback = 0;
595  if (TIDDetId(detid).isZMinusSide()) iszminusside = 1;
596  else iszminusside = 0;
597  isexternalstring = -1;
598  rodstringpetal = -1;
599  module_number = TIDDetId(detid).moduleNumber();
600 
601  break;
602 
603  case StripSubdetector::TOB :
604  layer_ring = TOBDetId(detid).layer();
605  disc = -1;
606  isback = -1;
607  if (TOBDetId(detid).isZMinusSide()) iszminusside = 1;
608  else iszminusside = 0;
609  isexternalstring = -1;
610  rodstringpetal = TOBDetId(detid).rodNumber();
611  module_number = TOBDetId(detid).moduleNumber();
612 
613  break;
614 
615  case StripSubdetector::TEC :
616  layer_ring = TECDetId(detid).ring();
617  disc = TECDetId(detid).wheel();
618  if (TECDetId(detid).isBackPetal()) isback = 1;
619  else isback = 0;
620  if (TECDetId(detid).isZMinusSide()) iszminusside = 1;
621  else iszminusside = 0;
622  isexternalstring = -1;
623  rodstringpetal = TECDetId(detid).petalNumber();
624  module_number = TECDetId(detid).moduleNumber();
625 
626  break;
627 
628  default :
629  std::cout << "### Detector does not belong to TIB, TID, TOB or TEC !? ###" << std::endl;
630  std::cout << "### DetRawId: " << detid << " ###" << std::endl;
631  }
632 }
633 
634 void SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy::setMinNumOfEvents()
635 {
636  minNevents_=absolute_occupancy_*Nevents_;
637 }
TOBDetId::rodNumber
unsigned int rodNumber() const
Definition: TOBDetId.h:77
edm::isDebugEnabled
bool isDebugEnabled()
Definition: MessageLogger.cc:68
i
int i
Definition: DBlmapReader.cc:9
TECDetId::petalNumber
unsigned int petalNumber() const
Definition: TECDetId.h:94
TIBDetId::stringNumber
unsigned int stringNumber() const
Definition: TIBDetId.h:87
TOBDetId::layer
unsigned int layer() const
layer id
Definition: TOBDetId.h:39
strip
void strip(std::string &input, const std::string &blanks=" \n\t")
Definition: stringTools.cc:16
PV3DBase::y
T y() const
Definition: PV3DBase.h:57
trackerHits.histo
tuple histo
Definition: trackerHits.py:12
SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy::CalculateMeanAndRMS
void CalculateMeanAndRMS(std::vector< Apv >, std::pair< double, double > *, int)
Definition: SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy.cc:332
StripGeomDetUnit::specificTopology
virtual const StripTopology & specificTopology() const
Returns a reference to the strip proxy topology.
Definition: StripGeomDetUnit.cc:19
edm::second
U second(std::pair< T, U > const &p)
Definition: ParameterSet.cc:185
PV3DBase::z
T z() const
Definition: PV3DBase.h:58
TIBDetId::moduleNumber
unsigned int moduleNumber() const
Definition: TIBDetId.h:91
SiStripQuality::compact
void compact(unsigned int &, std::vector< unsigned int > &)
Definition: SiStripQuality.cc:326
TIDDetId::ring
unsigned int ring() const
ring id
Definition: TIDDetId.h:55
j
int j
Definition: DBlmapReader.cc:9
first
bool first
Definition: L1TdeRCT.cc:79
SiStripQuality::fillBadComponents
void fillBadComponents()
Definition: SiStripQuality.cc:520
DetId::subdetId
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:39
TrackerGeometry::idToDet
virtual const GeomDet * idToDet(DetId) const
Definition: TrackerGeometry.cc:135
LogTrace
#define LogTrace(id)
Definition: MessageLogger.h:511
TECDetId::moduleNumber
unsigned int moduleNumber() const
Definition: TECDetId.h:102
SiStripDetId
Detector identifier class for the strip tracker.
Definition: SiStripDetId.h:17
DetId
Definition: DetId.h:20
SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy::AnalyzeOccupancy
void AnalyzeOccupancy(SiStripQuality *, std::vector< Apv > &, std::pair< double, double > *, std::vector< unsigned int > &, edm::ESHandle< SiStripQuality > &)
Definition: SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy.cc:387
TECDetId::wheel
unsigned int wheel() const
wheel id
Definition: TECDetId.h:52
TIBDetId::layer
unsigned int layer() const
layer id
Definition: TIBDetId.h:41
a
double a
Definition: hdecay.h:121
SiStripBadStrip::Range
std::pair< ContainerIterator, ContainerIterator > Range
Definition: SiStripBadStrip.h:50
TECDetId::ring
unsigned int ring() const
ring id
Definition: TECDetId.h:71
StripTopology::localPosition
virtual LocalPoint localPosition(float strip) const =0
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:41
SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy::extractBadAPVSandStrips
void extractBadAPVSandStrips(SiStripQuality *, HistoMap &, edm::ESHandle< SiStripQuality > &)
Definition: SiStripBadAPVandHotStripAlgorithmFromClusterOccupancy.cc:32
SiStripBadStrip::put
bool put(const uint32_t &detID, const InputVector &vect)
Definition: SiStripBadStrip.h:62
PV3DBase::x
T x() const
Definition: PV3DBase.h:56
TIDDetId::moduleNumber
unsigned int moduleNumber() const
Definition: TIDDetId.h:101
SiStripBadStrip::encode
unsigned int encode(const unsigned short &first, const unsigned short &NconsecutiveBadStrips, const unsigned short &flag=0)
Definition: SiStripBadStrip.h:83
TOBDetId::moduleNumber
unsigned int moduleNumber() const
Definition: TOBDetId.h:81
TIDDetId::wheel
unsigned int wheel() const
wheel id
Definition: TIDDetId.h:50
Generated for CMSSW Reference Manual by doxygen 1.8.5