CMS 3D CMS Logo

SiStripGainsPCLHarvester.cc
Go to the documentation of this file.
1 // system include files
2 #include <memory>
3 
4 // CMSSW includes
27 
28 // user include files
33 #include <iostream>
34 #include <sstream>
35 
36 //********************************************************************************//
38  : doStoreOnDB(false),
39  GOOD(0),
40  BAD(0),
41  MASKED(0),
42  NStripAPVs(0),
43  NPixelDets(0),
44  bareTkGeomPtr_(nullptr),
45  tTopo_(nullptr) {
46  m_Record = ps.getUntrackedParameter<std::string>("Record", "SiStripApvGainRcd");
47  CalibrationLevel = ps.getUntrackedParameter<int>("CalibrationLevel", 0);
48  MinNrEntries = ps.getUntrackedParameter<double>("minNrEntries", 20);
49  m_DQMdir = ps.getUntrackedParameter<std::string>("DQMdir", "AlCaReco/SiStripGains");
50  m_calibrationMode = ps.getUntrackedParameter<std::string>("calibrationMode", "StdBunch");
51  tagCondition_NClusters = ps.getUntrackedParameter<double>("NClustersForTagProd", 2E8);
52  tagCondition_GoodFrac = ps.getUntrackedParameter<double>("GoodFracForTagProd", 0.95);
53  doChargeMonitorPerPlane = ps.getUntrackedParameter<bool>("doChargeMonitorPerPlane", false);
54  VChargeHisto = ps.getUntrackedParameter<std::vector<std::string>>("ChargeHisto");
55 
56  //Set the monitoring element tag and store
57  dqm_tag_.reserve(7);
58  dqm_tag_.clear();
59  dqm_tag_.push_back("StdBunch"); // statistic collection from Standard Collision Bunch @ 3.8 T
60  dqm_tag_.push_back("StdBunch0T"); // statistic collection from Standard Collision Bunch @ 0 T
61  dqm_tag_.push_back("AagBunch"); // statistic collection from First Collision After Abort Gap @ 3.8 T
62  dqm_tag_.push_back("AagBunch0T"); // statistic collection from First Collision After Abort Gap @ 0 T
63  dqm_tag_.push_back("IsoMuon"); // statistic collection from Isolated Muon @ 3.8 T
64  dqm_tag_.push_back("IsoMuon0T"); // statistic collection from Isolated Muon @ 0 T
65  dqm_tag_.push_back("Harvest"); // statistic collection: Harvest
66 }
67 
68 //********************************************************************************//
69 // ------------ method called for each event ------------
71  using namespace edm;
72  static constexpr float defaultGainTick = 690. / 640.;
73 
74  this->checkBookAPVColls(iSetup); // check whether APV colls are booked and do so if not yet done
75  this->checkAndRetrieveTopology(iSetup);
76 
77  edm::ESHandle<SiStripGain> gainHandle;
78  iSetup.get<SiStripGainRcd>().get(gainHandle);
79  if (!gainHandle.isValid()) {
80  edm::LogError("SiStripGainPCLHarvester") << "gainHandle is not valid\n";
81  exit(0);
82  }
83 
84  edm::ESHandle<SiStripQuality> SiStripQuality_;
85  iSetup.get<SiStripQualityRcd>().get(SiStripQuality_);
86 
87  for (unsigned int a = 0; a < APVsCollOrdered.size(); a++) {
88  std::shared_ptr<stAPVGain> APV = APVsCollOrdered[a];
89 
90  if (APV->SubDet == PixelSubdetector::PixelBarrel || APV->SubDet == PixelSubdetector::PixelEndcap)
91  continue;
92 
93  APV->isMasked = SiStripQuality_->IsApvBad(APV->DetId, APV->APVId);
94 
95  if (gainHandle->getNumberOfTags() != 2) {
96  edm::LogError("SiStripGainPCLHarvester") << "NUMBER OF GAIN TAG IS EXPECTED TO BE 2\n";
97  fflush(stdout);
98  exit(0);
99  };
100  float newPreviousGain = gainHandle->getApvGain(APV->APVId, gainHandle->getRange(APV->DetId, 1), 1);
101  if (APV->PreviousGain != 1 and newPreviousGain != APV->PreviousGain)
102  edm::LogWarning("SiStripGainPCLHarvester") << "WARNING: ParticleGain in the global tag changed\n";
103  APV->PreviousGain = newPreviousGain;
104 
105  float newPreviousGainTick =
106  APV->isMasked ? defaultGainTick : gainHandle->getApvGain(APV->APVId, gainHandle->getRange(APV->DetId, 0), 0);
107  if (APV->PreviousGainTick != 1 and newPreviousGainTick != APV->PreviousGainTick) {
108  edm::LogWarning("SiStripGainPCLHarvester")
109  << "WARNING: TickMarkGain in the global tag changed\n"
110  << std::endl
111  << " APV->SubDet: " << APV->SubDet << " APV->APVId:" << APV->APVId << std::endl
112  << " APV->PreviousGainTick: " << APV->PreviousGainTick << " newPreviousGainTick: " << newPreviousGainTick
113  << std::endl;
114  }
115  APV->PreviousGainTick = newPreviousGainTick;
116  }
117 }
118 
119 //********************************************************************************//
121  edm::LogInfo("SiStripGainsPCLHarvester") << "Starting harvesting statistics" << std::endl;
122 
123  std::string DQM_dir = m_DQMdir;
124 
125  std::string stag = *(std::find(dqm_tag_.begin(), dqm_tag_.end(), m_calibrationMode));
126  if (!stag.empty() && stag[0] != '_')
127  stag.insert(0, 1, '_');
128 
129  std::string cvi = DQM_dir + std::string("/Charge_Vs_Index") + stag;
130 
131  MonitorElement* Charge_Vs_Index = igetter_.get(cvi);
132 
133  if (Charge_Vs_Index == nullptr) {
134  edm::LogError("SiStripGainsPCLHarvester")
135  << "Harvesting: could not retrieve " << cvi.c_str() << ", statistics will not be summed!" << std::endl;
136  } else {
137  edm::LogInfo("SiStripGainsPCLHarvester")
138  << "Harvesting " << (Charge_Vs_Index)->getTH2S()->GetEntries() << " more clusters" << std::endl;
139  }
140 
141  algoComputeMPVandGain(Charge_Vs_Index);
142  std::unique_ptr<SiStripApvGain> theAPVGains = this->getNewObject(Charge_Vs_Index);
143 
144  // write out the APVGains record
146 
147  if (doStoreOnDB) {
148  if (poolDbService.isAvailable())
149  poolDbService->writeOne(theAPVGains.get(), poolDbService->currentTime(), m_Record);
150  else
151  throw std::runtime_error("PoolDBService required.");
152  } else {
153  edm::LogInfo("SiStripGainsPCLHarvester") << "Will not produce payload!" << std::endl;
154  }
155 
156  //Collect the statistics for monitoring and validation
157  gainQualityMonitor(ibooker_, Charge_Vs_Index);
158 }
159 
160 //********************************************************************************//
162  const MonitorElement* Charge_Vs_Index) const {
163  ibooker_.setCurrentFolder("AlCaReco/SiStripGainsHarvesting/");
164 
165  std::vector<APVGain::APVmon> new_charge_histos;
166  std::vector<std::pair<std::string, std::string>> cnames =
168  for (unsigned int i = 0; i < cnames.size(); i++) {
169  MonitorElement* monitor = ibooker_.book1DD((cnames[i]).first, (cnames[i]).second.c_str(), 100, 0., 1000.);
170  int id = APVGain::subdetectorId((cnames[i]).first);
171  int side = APVGain::subdetectorSide((cnames[i]).first);
172  int plane = APVGain::subdetectorPlane((cnames[i]).first);
173  new_charge_histos.push_back(APVGain::APVmon(id, side, plane, monitor));
174  }
175 
176  int MPVbin = 300;
177  float MPVmin = 0.;
178  float MPVmax = 600.;
179 
180  MonitorElement* MPV_Vs_EtaTIB =
181  ibooker_.book2DD("MPVvsEtaTIB", "MPV vs Eta TIB", 50, -3.0, 3.0, MPVbin, MPVmin, MPVmax);
182  MonitorElement* MPV_Vs_EtaTID =
183  ibooker_.book2DD("MPVvsEtaTID", "MPV vs Eta TID", 50, -3.0, 3.0, MPVbin, MPVmin, MPVmax);
184  MonitorElement* MPV_Vs_EtaTOB =
185  ibooker_.book2DD("MPVvsEtaTOB", "MPV vs Eta TOB", 50, -3.0, 3.0, MPVbin, MPVmin, MPVmax);
186  MonitorElement* MPV_Vs_EtaTEC =
187  ibooker_.book2DD("MPVvsEtaTEC", "MPV vs Eta TEC", 50, -3.0, 3.0, MPVbin, MPVmin, MPVmax);
188  MonitorElement* MPV_Vs_EtaTECthin =
189  ibooker_.book2DD("MPVvsEtaTEC1", "MPV vs Eta TEC-thin", 50, -3.0, 3.0, MPVbin, MPVmin, MPVmax);
190  MonitorElement* MPV_Vs_EtaTECthick =
191  ibooker_.book2DD("MPVvsEtaTEC2", "MPV vs Eta TEC-thick", 50, -3.0, 3.0, MPVbin, MPVmin, MPVmax);
192 
193  MonitorElement* MPV_Vs_PhiTIB =
194  ibooker_.book2DD("MPVvsPhiTIB", "MPV vs Phi TIB", 50, -3.4, 3.4, MPVbin, MPVmin, MPVmax);
195  MonitorElement* MPV_Vs_PhiTID =
196  ibooker_.book2DD("MPVvsPhiTID", "MPV vs Phi TID", 50, -3.4, 3.4, MPVbin, MPVmin, MPVmax);
197  MonitorElement* MPV_Vs_PhiTOB =
198  ibooker_.book2DD("MPVvsPhiTOB", "MPV vs Phi TOB", 50, -3.4, 3.4, MPVbin, MPVmin, MPVmax);
199  MonitorElement* MPV_Vs_PhiTEC =
200  ibooker_.book2DD("MPVvsPhiTEC", "MPV vs Phi TEC", 50, -3.4, 3.4, MPVbin, MPVmin, MPVmax);
201  MonitorElement* MPV_Vs_PhiTECthin =
202  ibooker_.book2DD("MPVvsPhiTEC1", "MPV vs Phi TEC-thin ", 50, -3.4, 3.4, MPVbin, MPVmin, MPVmax);
203  MonitorElement* MPV_Vs_PhiTECthick =
204  ibooker_.book2DD("MPVvsPhiTEC2", "MPV vs Phi TEC-thick", 50, -3.4, 3.4, MPVbin, MPVmin, MPVmax);
205 
206  MonitorElement* NoMPVfit = ibooker_.book2DD("NoMPVfit", "Modules with bad Landau Fit", 350, -350, 350, 240, 0, 120);
207  MonitorElement* NoMPVmasked = ibooker_.book2DD("NoMPVmasked", "Masked Modules", 350, -350, 350, 240, 0, 120);
208 
209  MonitorElement* Gains = ibooker_.book1DD("Gains", "Gains", 300, 0, 2);
210  MonitorElement* MPVs = ibooker_.book1DD("MPVs", "MPVs", MPVbin, MPVmin, MPVmax);
211  MonitorElement* MPVs320 = ibooker_.book1DD("MPV_320", "MPV 320 thickness", MPVbin, MPVmin, MPVmax);
212  MonitorElement* MPVs500 = ibooker_.book1DD("MPV_500", "MPV 500 thickness", MPVbin, MPVmin, MPVmax);
213  MonitorElement* MPVsTIB = ibooker_.book1DD("MPV_TIB", "MPV TIB", MPVbin, MPVmin, MPVmax);
214  MonitorElement* MPVsTID = ibooker_.book1DD("MPV_TID", "MPV TID", MPVbin, MPVmin, MPVmax);
215  MonitorElement* MPVsTIDP = ibooker_.book1DD("MPV_TIDP", "MPV TIDP", MPVbin, MPVmin, MPVmax);
216  MonitorElement* MPVsTIDM = ibooker_.book1DD("MPV_TIDM", "MPV TIDM", MPVbin, MPVmin, MPVmax);
217  MonitorElement* MPVsTOB = ibooker_.book1DD("MPV_TOB", "MPV TOB", MPVbin, MPVmin, MPVmax);
218  MonitorElement* MPVsTEC = ibooker_.book1DD("MPV_TEC", "MPV TEC", MPVbin, MPVmin, MPVmax);
219  MonitorElement* MPVsTECP = ibooker_.book1DD("MPV_TECP", "MPV TECP", MPVbin, MPVmin, MPVmax);
220  MonitorElement* MPVsTECM = ibooker_.book1DD("MPV_TECM", "MPV TECM", MPVbin, MPVmin, MPVmax);
221  MonitorElement* MPVsTECthin = ibooker_.book1DD("MPV_TEC1", "MPV TEC thin", MPVbin, MPVmin, MPVmax);
222  MonitorElement* MPVsTECthick = ibooker_.book1DD("MPV_TEC2", "MPV TEC thick", MPVbin, MPVmin, MPVmax);
223  MonitorElement* MPVsTECP1 = ibooker_.book1DD("MPV_TECP1", "MPV TECP thin ", MPVbin, MPVmin, MPVmax);
224  MonitorElement* MPVsTECP2 = ibooker_.book1DD("MPV_TECP2", "MPV TECP thick", MPVbin, MPVmin, MPVmax);
225  MonitorElement* MPVsTECM1 = ibooker_.book1DD("MPV_TECM1", "MPV TECM thin", MPVbin, MPVmin, MPVmax);
226  MonitorElement* MPVsTECM2 = ibooker_.book1DD("MPV_TECM2", "MPV TECM thick", MPVbin, MPVmin, MPVmax);
227 
228  MonitorElement* MPVError = ibooker_.book1DD("MPVError", "MPV Error", 150, 0, 150);
229  MonitorElement* MPVErrorVsMPV = ibooker_.book2DD("MPVErrorVsMPV", "MPV Error vs MPV", 300, 0, 600, 150, 0, 150);
230  MonitorElement* MPVErrorVsEta = ibooker_.book2DD("MPVErrorVsEta", "MPV Error vs Eta", 50, -3.0, 3.0, 150, 0, 150);
231  MonitorElement* MPVErrorVsPhi = ibooker_.book2DD("MPVErrorVsPhi", "MPV Error vs Phi", 50, -3.4, 3.4, 150, 0, 150);
232  MonitorElement* MPVErrorVsN = ibooker_.book2DD("MPVErrorVsN", "MPV Error vs N", 500, 0, 1000, 150, 0, 150);
233 
234  MonitorElement* DiffWRTPrevGainTIB = ibooker_.book1DD("DiffWRTPrevGainTIB", "Diff w.r.t. PrevGain TIB", 250, 0, 2);
235  MonitorElement* DiffWRTPrevGainTID = ibooker_.book1DD("DiffWRTPrevGainTID", "Diff w.r.t. PrevGain TID", 250, 0, 2);
236  MonitorElement* DiffWRTPrevGainTOB = ibooker_.book1DD("DiffWRTPrevGainTOB", "Diff w.r.t. PrevGain TOB", 250, 0, 2);
237  MonitorElement* DiffWRTPrevGainTEC = ibooker_.book1DD("DiffWRTPrevGainTEC", "Diff w.r.t. PrevGain TEC", 250, 0, 2);
238 
239  MonitorElement* GainVsPrevGainTIB =
240  ibooker_.book2DD("GainVsPrevGainTIB", "Gain vs PrevGain TIB", 100, 0, 2, 100, 0, 2);
241  MonitorElement* GainVsPrevGainTID =
242  ibooker_.book2DD("GainVsPrevGainTID", "Gain vs PrevGain TID", 100, 0, 2, 100, 0, 2);
243  MonitorElement* GainVsPrevGainTOB =
244  ibooker_.book2DD("GainVsPrevGainTOB", "Gain vs PrevGain TOB", 100, 0, 2, 100, 0, 2);
245  MonitorElement* GainVsPrevGainTEC =
246  ibooker_.book2DD("GainVsPrevGainTEC", "Gain vs PrevGain TEC", 100, 0, 2, 100, 0, 2);
247 
248  for (unsigned int a = 0; a < APVsCollOrdered.size(); a++) {
249  std::shared_ptr<stAPVGain> APV = APVsCollOrdered[a];
250  if (APV == nullptr)
251  continue;
252 
253  unsigned int Index = APV->Index;
254  unsigned int SubDet = APV->SubDet;
255  unsigned int DetId = APV->DetId;
256  float z = APV->z;
257  float Eta = APV->Eta;
258  float R = APV->R;
259  float Phi = APV->Phi;
260  float Thickness = APV->Thickness;
261  double FitMPV = APV->FitMPV;
262  double FitMPVErr = APV->FitMPVErr;
263  double Gain = APV->Gain;
264  double NEntries = APV->NEntries;
265  double PreviousGain = APV->PreviousGain;
266 
267  if (SubDet < 3)
268  continue; // avoid to loop over Pixel det id
269 
270  if (Gain != 1.) {
271  std::vector<MonitorElement*> charge_histos = APVGain::FetchMonitor(new_charge_histos, DetId, tTopo_);
272  TH2S* chvsidx = (Charge_Vs_Index)->getTH2S();
273  int bin = chvsidx->GetXaxis()->FindBin(Index);
274  TH1D* Proj = chvsidx->ProjectionY("proj", bin, bin);
275  for (int binId = 0; binId < Proj->GetXaxis()->GetNbins(); binId++) {
276  double new_charge = Proj->GetXaxis()->GetBinCenter(binId) / Gain;
277  if (Proj->GetBinContent(binId) != 0.) {
278  for (unsigned int h = 0; h < charge_histos.size(); h++) {
279  TH1D* chisto = (charge_histos[h])->getTH1D();
280  for (int e = 0; e < Proj->GetBinContent(binId); e++)
281  chisto->Fill(new_charge);
282  }
283  }
284  }
285  }
286 
287  if (FitMPV <= 0.) { // No fit of MPV
288  if (APV->isMasked)
289  NoMPVmasked->Fill(z, R);
290  else
291  NoMPVfit->Fill(z, R);
292 
293  } else { // Fit of MPV
294  if (FitMPV > 0.)
295  Gains->Fill(Gain);
296 
297  MPVs->Fill(FitMPV);
298  if (Thickness < 0.04)
299  MPVs320->Fill(FitMPV);
300  if (Thickness > 0.04)
301  MPVs500->Fill(FitMPV);
302 
303  MPVError->Fill(FitMPVErr);
304  MPVErrorVsMPV->Fill(FitMPV, FitMPVErr);
305  MPVErrorVsEta->Fill(Eta, FitMPVErr);
306  MPVErrorVsPhi->Fill(Phi, FitMPVErr);
307  MPVErrorVsN->Fill(NEntries, FitMPVErr);
308 
309  if (SubDet == 3) {
310  MPV_Vs_EtaTIB->Fill(Eta, FitMPV);
311  MPV_Vs_PhiTIB->Fill(Phi, FitMPV);
312  MPVsTIB->Fill(FitMPV);
313 
314  } else if (SubDet == 4) {
315  MPV_Vs_EtaTID->Fill(Eta, FitMPV);
316  MPV_Vs_PhiTID->Fill(Phi, FitMPV);
317  MPVsTID->Fill(FitMPV);
318  if (Eta < 0.)
319  MPVsTIDM->Fill(FitMPV);
320  if (Eta > 0.)
321  MPVsTIDP->Fill(FitMPV);
322 
323  } else if (SubDet == 5) {
324  MPV_Vs_EtaTOB->Fill(Eta, FitMPV);
325  MPV_Vs_PhiTOB->Fill(Phi, FitMPV);
326  MPVsTOB->Fill(FitMPV);
327 
328  } else if (SubDet == 6) {
329  MPV_Vs_EtaTEC->Fill(Eta, FitMPV);
330  MPV_Vs_PhiTEC->Fill(Phi, FitMPV);
331  MPVsTEC->Fill(FitMPV);
332  if (Eta < 0.)
333  MPVsTECM->Fill(FitMPV);
334  if (Eta > 0.)
335  MPVsTECP->Fill(FitMPV);
336  if (Thickness < 0.04) {
337  MPV_Vs_EtaTECthin->Fill(Eta, FitMPV);
338  MPV_Vs_PhiTECthin->Fill(Phi, FitMPV);
339  MPVsTECthin->Fill(FitMPV);
340  if (Eta > 0.)
341  MPVsTECP1->Fill(FitMPV);
342  if (Eta < 0.)
343  MPVsTECM1->Fill(FitMPV);
344  }
345  if (Thickness > 0.04) {
346  MPV_Vs_EtaTECthick->Fill(Eta, FitMPV);
347  MPV_Vs_PhiTECthick->Fill(Phi, FitMPV);
348  MPVsTECthick->Fill(FitMPV);
349  if (Eta > 0.)
350  MPVsTECP2->Fill(FitMPV);
351  if (Eta < 0.)
352  MPVsTECM2->Fill(FitMPV);
353  }
354  }
355  }
356 
357  if (SubDet == 3 && PreviousGain != 0.)
358  DiffWRTPrevGainTIB->Fill(Gain / PreviousGain);
359  else if (SubDet == 4 && PreviousGain != 0.)
360  DiffWRTPrevGainTID->Fill(Gain / PreviousGain);
361  else if (SubDet == 5 && PreviousGain != 0.)
362  DiffWRTPrevGainTOB->Fill(Gain / PreviousGain);
363  else if (SubDet == 6 && PreviousGain != 0.)
364  DiffWRTPrevGainTEC->Fill(Gain / PreviousGain);
365 
366  if (SubDet == 3)
367  GainVsPrevGainTIB->Fill(PreviousGain, Gain);
368  else if (SubDet == 4)
369  GainVsPrevGainTID->Fill(PreviousGain, Gain);
370  else if (SubDet == 5)
371  GainVsPrevGainTOB->Fill(PreviousGain, Gain);
372  else if (SubDet == 6)
373  GainVsPrevGainTEC->Fill(PreviousGain, Gain);
374  }
375 }
376 
377 //********************************************************************************//
379  unsigned int I = 0;
380  TH1F* Proj = nullptr;
381  double FitResults[6];
382  double MPVmean = 300;
383 
384  if (Charge_Vs_Index == nullptr) {
385  edm::LogError("SiStripGainsPCLHarvester")
386  << "Harvesting: could not execute algoComputeMPVandGain method because " << m_calibrationMode
387  << " statistics cannot be retrieved.\n"
388  << "Please check if input contains " << m_calibrationMode << " data." << std::endl;
389  return;
390  }
391 
392  TH2S* chvsidx = (Charge_Vs_Index)->getTH2S();
393 
394  printf("Progressing Bar :0%% 20%% 40%% 60%% 80%% 100%%\n");
395  printf("Fitting Charge Distribution :");
396  int TreeStep = APVsColl.size() / 50;
397 
398  for (auto it = APVsColl.begin(); it != APVsColl.end(); it++, I++) {
399  if (I % TreeStep == 0) {
400  printf(".");
401  fflush(stdout);
402  }
403  std::shared_ptr<stAPVGain> APV = it->second;
404  if (APV->Bin < 0)
405  APV->Bin = chvsidx->GetXaxis()->FindBin(APV->Index);
406 
407  if (APV->isMasked) {
408  APV->Gain = APV->PreviousGain;
409  MASKED++;
410  continue;
411  }
412 
413  Proj = (TH1F*)(chvsidx->ProjectionY(
414  "", chvsidx->GetXaxis()->FindBin(APV->Index), chvsidx->GetXaxis()->FindBin(APV->Index), "e"));
415  if (!Proj)
416  continue;
417 
418  if (CalibrationLevel == 0) {
419  } else if (CalibrationLevel == 1) {
420  int SecondAPVId = APV->APVId;
421  if (SecondAPVId % 2 == 0) {
422  SecondAPVId = SecondAPVId + 1;
423  } else {
424  SecondAPVId = SecondAPVId - 1;
425  }
426  std::shared_ptr<stAPVGain> APV2 = APVsColl[(APV->DetId << 4) | SecondAPVId];
427  if (APV2->Bin < 0)
428  APV2->Bin = chvsidx->GetXaxis()->FindBin(APV2->Index);
429  TH1F* Proj2 = (TH1F*)(chvsidx->ProjectionY("", APV2->Bin, APV2->Bin, "e"));
430  if (Proj2) {
431  Proj->Add(Proj2, 1);
432  delete Proj2;
433  }
434  } else if (CalibrationLevel == 2) {
435  for (unsigned int i = 0; i < 16; i++) { //loop up to 6APV for Strip and up to 16 for Pixels
436  auto tmpit = APVsColl.find((APV->DetId << 4) | i);
437  if (tmpit == APVsColl.end())
438  continue;
439  std::shared_ptr<stAPVGain> APV2 = tmpit->second;
440  if (APV2->DetId != APV->DetId || APV2->APVId == APV->APVId)
441  continue;
442  if (APV2->Bin < 0)
443  APV2->Bin = chvsidx->GetXaxis()->FindBin(APV2->Index);
444  TH1F* Proj2 = (TH1F*)(chvsidx->ProjectionY("", APV2->Bin, APV2->Bin, "e"));
445  if (Proj2) {
446  Proj->Add(Proj2, 1);
447  delete Proj2;
448  }
449  }
450  } else {
451  CalibrationLevel = 0;
452  printf("Unknown Calibration Level, will assume %i\n", CalibrationLevel);
453  }
454 
455  getPeakOfLandau(Proj, FitResults);
456  APV->FitMPV = FitResults[0];
457  APV->FitMPVErr = FitResults[1];
458  APV->FitWidth = FitResults[2];
459  APV->FitWidthErr = FitResults[3];
460  APV->FitChi2 = FitResults[4];
461  APV->FitNorm = FitResults[5];
462  APV->NEntries = Proj->GetEntries();
463 
464  if (IsGoodLandauFit(FitResults)) {
465  APV->Gain = APV->FitMPV / MPVmean;
466  if (APV->SubDet > 2)
467  GOOD++;
468  } else {
469  APV->Gain = APV->PreviousGain;
470  if (APV->SubDet > 2)
471  BAD++;
472  }
473  if (APV->Gain <= 0)
474  APV->Gain = 1;
475 
476  delete Proj;
477  }
478  printf("\n");
479 }
480 
481 //********************************************************************************//
482 void SiStripGainsPCLHarvester::getPeakOfLandau(TH1* InputHisto, double* FitResults, double LowRange, double HighRange) {
483  FitResults[0] = -0.5; //MPV
484  FitResults[1] = 0; //MPV error
485  FitResults[2] = -0.5; //Width
486  FitResults[3] = 0; //Width error
487  FitResults[4] = -0.5; //Fit Chi2/NDF
488  FitResults[5] = 0; //Normalization
489 
490  if (InputHisto->GetEntries() < MinNrEntries)
491  return;
492 
493  // perform fit with standard landau
494  TF1 MyLandau("MyLandau", "landau", LowRange, HighRange);
495  MyLandau.SetParameter(1, 300);
496  InputHisto->Fit(&MyLandau, "0QR WW");
497 
498  // MPV is parameter 1 (0=constant, 1=MPV, 2=Sigma)
499  FitResults[0] = MyLandau.GetParameter(1); //MPV
500  FitResults[1] = MyLandau.GetParError(1); //MPV error
501  FitResults[2] = MyLandau.GetParameter(2); //Width
502  FitResults[3] = MyLandau.GetParError(2); //Width error
503  FitResults[4] = MyLandau.GetChisquare() / MyLandau.GetNDF(); //Fit Chi2/NDF
504  FitResults[5] = MyLandau.GetParameter(0);
505 }
506 
507 //********************************************************************************//
509  if (FitResults[0] <= 0)
510  return false;
511  // if(FitResults[1] > MaxMPVError )return false;
512  // if(FitResults[4] > MaxChi2OverNDF)return false;
513  return true;
514 }
515 
516 //********************************************************************************//
517 // ------------ method called once each job just before starting event loop ------------
520  const TrackerGeometry* newBareTkGeomPtr = &(*tkGeom_);
521  if (newBareTkGeomPtr == bareTkGeomPtr_)
522  return; // already filled APVColls, nothing changed
523 
524  if (!bareTkGeomPtr_) { // pointer not yet set: called the first time => fill the APVColls
525  auto const& Det = newBareTkGeomPtr->dets();
526 
527  unsigned int Index = 0;
528 
529  for (unsigned int i = 0; i < Det.size(); i++) {
530  DetId Detid = Det[i]->geographicalId();
531  int SubDet = Detid.subdetId();
532 
533  if (SubDet == StripSubdetector::TIB || SubDet == StripSubdetector::TID || SubDet == StripSubdetector::TOB ||
534  SubDet == StripSubdetector::TEC) {
535  auto DetUnit = dynamic_cast<const StripGeomDetUnit*>(Det[i]);
536  if (!DetUnit)
537  continue;
538 
539  const StripTopology& Topo = DetUnit->specificTopology();
540  unsigned int NAPV = Topo.nstrips() / 128;
541 
542  for (unsigned int j = 0; j < NAPV; j++) {
543  auto APV = std::make_shared<stAPVGain>();
544  APV->Index = Index;
545  APV->Bin = -1;
546  APV->DetId = Detid.rawId();
547  APV->APVId = j;
548  APV->SubDet = SubDet;
549  APV->FitMPV = -1;
550  APV->FitMPVErr = -1;
551  APV->FitWidth = -1;
552  APV->FitWidthErr = -1;
553  APV->FitChi2 = -1;
554  APV->FitNorm = -1;
555  APV->Gain = -1;
556  APV->PreviousGain = 1;
557  APV->PreviousGainTick = 1;
558  APV->x = DetUnit->position().basicVector().x();
559  APV->y = DetUnit->position().basicVector().y();
560  APV->z = DetUnit->position().basicVector().z();
561  APV->Eta = DetUnit->position().basicVector().eta();
562  APV->Phi = DetUnit->position().basicVector().phi();
563  APV->R = DetUnit->position().basicVector().transverse();
564  APV->Thickness = DetUnit->surface().bounds().thickness();
565  APV->NEntries = 0;
566  APV->isMasked = false;
567 
568  APVsCollOrdered.push_back(APV);
569  APVsColl[(APV->DetId << 4) | APV->APVId] = APV;
570  Index++;
571  NStripAPVs++;
572  } // loop on APVs
573  } // if is Strips
574  } // loop on dets
575 
576  for (unsigned int i = 0; i < Det.size();
577  i++) { //Make two loop such that the Pixel information is added at the end --> make transition simpler
578  DetId Detid = Det[i]->geographicalId();
579  int SubDet = Detid.subdetId();
581  auto DetUnit = dynamic_cast<const PixelGeomDetUnit*>(Det[i]);
582  if (!DetUnit)
583  continue;
584 
585  const PixelTopology& Topo = DetUnit->specificTopology();
586  unsigned int NROCRow = Topo.nrows() / (80.);
587  unsigned int NROCCol = Topo.ncolumns() / (52.);
588 
589  for (unsigned int j = 0; j < NROCRow; j++) {
590  for (unsigned int i = 0; i < NROCCol; i++) {
591  auto APV = std::make_shared<stAPVGain>();
592  APV->Index = Index;
593  APV->Bin = -1;
594  APV->DetId = Detid.rawId();
595  APV->APVId = (j << 3 | i);
596  APV->SubDet = SubDet;
597  APV->FitMPV = -1;
598  APV->FitMPVErr = -1;
599  APV->FitWidth = -1;
600  APV->FitWidthErr = -1;
601  APV->FitChi2 = -1;
602  APV->Gain = -1;
603  APV->PreviousGain = 1;
604  APV->PreviousGainTick = 1;
605  APV->x = DetUnit->position().basicVector().x();
606  APV->y = DetUnit->position().basicVector().y();
607  APV->z = DetUnit->position().basicVector().z();
608  APV->Eta = DetUnit->position().basicVector().eta();
609  APV->Phi = DetUnit->position().basicVector().phi();
610  APV->R = DetUnit->position().basicVector().transverse();
611  APV->Thickness = DetUnit->surface().bounds().thickness();
612  APV->isMasked = false; //SiPixelQuality_->IsModuleBad(Detid.rawId());
613  APV->NEntries = 0;
614 
615  APVsCollOrdered.push_back(APV);
616  APVsColl[(APV->DetId << 4) | APV->APVId] = APV;
617  Index++;
618  NPixelDets++;
619 
620  } // loop on ROC cols
621  } // loop on ROC rows
622  } // if Pixel
623  } // loop on Dets
624  } //if (!bareTkGeomPtr_) ...
625  bareTkGeomPtr_ = newBareTkGeomPtr;
626 }
627 
629  if (!tTopo_) {
631  setup.get<TrackerTopologyRcd>().get(TopoHandle);
632  tTopo_ = TopoHandle.product();
633  }
634 }
635 
636 //********************************************************************************//
638  // The goal of this function is to check wether or not there is enough statistics
639  // to produce a meaningful tag for the DB
640 
641  if (Charge_Vs_Index == nullptr) {
642  edm::LogError("SiStripGainsPCLHarvester")
643  << "produceTagFilter -> Return false: could not retrieve the " << m_calibrationMode << " statistics.\n"
644  << "Please check if input contains " << m_calibrationMode << " data." << std::endl;
645  return false;
646  }
647 
648  float integral = (Charge_Vs_Index)->getTH2S()->Integral();
649  if ((Charge_Vs_Index)->getTH2S()->Integral(0, NStripAPVs + 1, 0, 99999) < tagCondition_NClusters) {
650  edm::LogWarning("SiStripGainsPCLHarvester")
651  << "calibrationMode -> " << m_calibrationMode << "\n"
652  << "produceTagFilter -> Return false: Statistics is too low : " << integral << std::endl;
653  return false;
654  }
655  if ((1.0 * GOOD) / (GOOD + BAD) < tagCondition_GoodFrac) {
656  edm::LogWarning("SiStripGainsPCLHarvester")
657  << "calibrationMode -> " << m_calibrationMode << "\n"
658  << "produceTagFilter -> Return false: ratio of GOOD/TOTAL is too low: " << (1.0 * GOOD) / (GOOD + BAD)
659  << std::endl;
660  return false;
661  }
662  return true;
663 }
664 
665 //********************************************************************************//
666 std::unique_ptr<SiStripApvGain> SiStripGainsPCLHarvester::getNewObject(const MonitorElement* Charge_Vs_Index) {
667  std::unique_ptr<SiStripApvGain> obj = std::unique_ptr<SiStripApvGain>(new SiStripApvGain());
668 
669  if (!produceTagFilter(Charge_Vs_Index)) {
670  edm::LogWarning("SiStripGainsPCLHarvester")
671  << "getNewObject -> will not produce a paylaod because produceTagFilter returned false " << std::endl;
672  return obj;
673  } else {
674  doStoreOnDB = true;
675  }
676 
677  std::vector<float> theSiStripVector;
678  unsigned int PreviousDetId = 0;
679  for (unsigned int a = 0; a < APVsCollOrdered.size(); a++) {
680  std::shared_ptr<stAPVGain> APV = APVsCollOrdered[a];
681  if (APV == nullptr) {
682  printf("Bug\n");
683  continue;
684  }
685  if (APV->SubDet <= 2)
686  continue;
687  if (APV->DetId != PreviousDetId) {
688  if (!theSiStripVector.empty()) {
689  SiStripApvGain::Range range(theSiStripVector.begin(), theSiStripVector.end());
690  if (!obj->put(PreviousDetId, range))
691  printf("Bug to put detId = %i\n", PreviousDetId);
692  }
693  theSiStripVector.clear();
694  PreviousDetId = APV->DetId;
695  }
696  theSiStripVector.push_back(APV->Gain);
697 
698  LogDebug("SiStripGainsPCLHarvester") << " DetId: " << APV->DetId << " APV: " << APV->APVId
699  << " Gain: " << APV->Gain << std::endl;
700  }
701  if (!theSiStripVector.empty()) {
702  SiStripApvGain::Range range(theSiStripVector.begin(), theSiStripVector.end());
703  if (!obj->put(PreviousDetId, range))
704  printf("Bug to put detId = %i\n", PreviousDetId);
705  }
706 
707  return obj;
708 }
709 
710 //********************************************************************************//
711 // ------------ method fills 'descriptions' with the allowed parameters for the module ------------
714  desc.setUnknown();
715  descriptions.addDefault(desc);
716 }
717 
718 //********************************************************************************//
#define LogDebug(id)
std::unique_ptr< SiStripApvGain > getNewObject(const MonitorElement *Charge_Vs_Index)
bool IsApvBad(const uint32_t &detid, const short &apvNb) const
static constexpr auto TEC
T getUntrackedParameter(std::string const &, T const &) const
FWCore Framework interface EventSetupRecordImplementation h
Helper function to determine trigger accepts.
virtual int nrows() const =0
edm::ESHandle< TrackerGeometry > tkGeom_
void dqmEndJob(DQMStore::IBooker &ibooker_, DQMStore::IGetter &igetter_) override
size_t getNumberOfTags() const
Definition: SiStripGain.h:95
void setCurrentFolder(std::string const &fullpath)
Definition: DQMStore.cc:418
MonitorElement * book1DD(TString const &name, TString const &title, int nchX, double lowX, double highX)
Definition: DQMStore.cc:257
std::vector< MonitorElement * > FetchMonitor(std::vector< APVmon >, uint32_t, const TrackerTopology *topo=0)
#define nullptr
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:57
const TrackerTopology * tTopo_
int subdetectorPlane(uint32_t, const TrackerTopology *)
std::vector< std::string > dqm_tag_
void beginRun(edm::Run const &run, edm::EventSetup const &isetup) override
bool produceTagFilter(const MonitorElement *Charge_Vs_Index)
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:19
void getPeakOfLandau(TH1 *InputHisto, double *FitResults, double LowRange=50, double HighRange=5400)
void algoComputeMPVandGain(const MonitorElement *Charge_Vs_Index)
static float getApvGain(const uint16_t &apv, const SiStripApvGain::Range &range)
Definition: SiStripGain.h:76
void Fill(long long x)
U second(std::pair< T, U > const &p)
SiStripGainsPCLHarvester(const edm::ParameterSet &ps)
const DetContainer & dets() const override
Returm a vector of all GeomDet (including all GeomDetUnits)
Data available (reporting)
void addDefault(ParameterSetDescription const &psetDescription)
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
bool isAvailable() const
Definition: Service.h:40
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:48
std::pair< ContainerIterator, ContainerIterator > Range
const std::complex< double > I
Definition: I.h:8
void writeOne(T *payload, Time_t time, const std::string &recordName, bool withlogging=false)
double const GOOD
Definition: Constants.h:13
static constexpr auto TOB
std::vector< std::string > VChargeHisto
Integral< F, X >::type integral(const F &f)
Definition: Integral.h:70
virtual void checkAndRetrieveTopology(const edm::EventSetup &setup)
Definition: DetId.h:17
bool IsGoodLandauFit(double *FitResults)
double const BAD
Definition: Constants.h:15
static constexpr auto TIB
void gainQualityMonitor(DQMStore::IBooker &ibooker_, const MonitorElement *Charge_Vs_Index) const
virtual int nstrips() const =0
int subdetectorSide(uint32_t, const TrackerTopology *)
const TrackerGeometry * bareTkGeomPtr_
void endRun(edm::Run const &run, edm::EventSetup const &isetup) override
virtual void checkBookAPVColls(const edm::EventSetup &setup)
std::vector< std::shared_ptr< stAPVGain > > APVsCollOrdered
HLT enums.
virtual int ncolumns() const =0
MonitorElement * book2DD(TString const &name, TString const &title, int nchX, double lowX, double highX, int nchY, double lowY, double highY)
Definition: DQMStore.cc:302
double a
Definition: hdecay.h:119
T get() const
Definition: EventSetup.h:73
std::vector< std::pair< std::string, std::string > > monHnames(std::vector< std::string >, bool, const char *tag)
bool isValid() const
Definition: ESHandle.h:44
MonitorElement * get(std::string const &path)
Definition: DQMStore.cc:437
static constexpr auto TID
T const * product() const
Definition: ESHandle.h:86
int subdetectorId(uint32_t)
#define constexpr
Definition: Run.h:45
constexpr DetId()
Create an empty or null id (also for persistence)
Definition: DetId.h:38
const SiStripApvGain::Range getRange(uint32_t detID) const
Definition: SiStripGain.h:71
def exit(msg="")
std::unordered_map< unsigned int, std::shared_ptr< stAPVGain > > APVsColl