CMS 3D CMS Logo

ZeeCalibration.cc
Go to the documentation of this file.
1 #include <fstream>
2 #include <iostream>
3 #include <utility>
4 
5 #include <TBranch.h>
6 #include <TCanvas.h>
7 #include <TF1.h>
8 #include <TProfile.h>
9 #include <TRandom.h>
10 
11 #include <CLHEP/Vector/LorentzVector.h>
12 
35 
36 #define MZ 91.1876
37 
38 #define DEBUG 1
39 
41  : hlTriggerResults_(iConfig.getParameter<edm::InputTag>("HLTriggerResults")),
42  mcProducer_(iConfig.getUntrackedParameter<std::string>("mcProducer", "")),
43  rechitProducer_(iConfig.getParameter<std::string>("rechitProducer")),
44  rechitCollection_(iConfig.getParameter<std::string>("rechitCollection")),
45  erechitProducer_(iConfig.getParameter<std::string>("erechitProducer")),
46  erechitCollection_(iConfig.getParameter<std::string>("erechitCollection")),
47  scProducer_(iConfig.getParameter<std::string>("scProducer")),
48  scCollection_(iConfig.getParameter<std::string>("scCollection")),
49  scIslandProducer_(iConfig.getParameter<std::string>("scIslandProducer")),
50  scIslandCollection_(iConfig.getParameter<std::string>("scIslandCollection")),
51  electronProducer_(iConfig.getParameter<std::string>("electronProducer")),
52  electronCollection_(iConfig.getParameter<std::string>("electronCollection")),
53  trigResultsToken_(consumes<edm::TriggerResults>(hlTriggerResults_)),
54  hepMCToken_(consumes<edm::HepMCProduct>(edm::InputTag(mcProducer_))),
55  ebRecHitToken_(consumes<EBRecHitCollection>(edm::InputTag(rechitProducer_, rechitCollection_))),
56  eeRecHitToken_(consumes<EERecHitCollection>(edm::InputTag(erechitProducer_, erechitCollection_))),
57  scToken_(consumes<reco::SuperClusterCollection>(edm::InputTag(scProducer_, scCollection_))),
58  islandSCToken_(consumes<reco::SuperClusterCollection>(edm::InputTag(scIslandProducer_, scIslandCollection_))),
59  gsfElectronToken_(consumes<reco::GsfElectronCollection>(edm::InputTag(electronProducer_, electronCollection_))),
60  geometryToken_(esConsumes()) {
61 #ifdef DEBUG
62  std::cout << "[ZeeCalibration] Starting the ctor" << std::endl;
63 #endif
64 
65  theMaxLoops = iConfig.getUntrackedParameter<unsigned int>("maxLoops", 0);
66 
67  wantEtaCorrection_ = iConfig.getUntrackedParameter<bool>("wantEtaCorrection", true);
68 
69  outputFileName_ = iConfig.getParameter<std::string>("outputFile");
70 
71  minInvMassCut_ = iConfig.getUntrackedParameter<double>("minInvMassCut", 70.);
72  maxInvMassCut_ = iConfig.getUntrackedParameter<double>("maxInvMassCut", 110.);
73 
74  calibMode_ = iConfig.getUntrackedParameter<std::string>("ZCalib_CalibType");
75 
76  outputFile_ = TFile::Open(outputFileName_.c_str(), "RECREATE"); // open output file to store histograms
77 
78  myTree = new TTree("myTree", "myTree");
79  // myTree->Branch("zMass","zMass", &mass);
80  myTree->Branch("zMass", &mass4tree, "mass/F");
81  myTree->Branch("zMassDiff", &massDiff4tree, "massDiff/F");
82 
83  barrelfile_ = iConfig.getUntrackedParameter<std::string>("initialMiscalibrationBarrel", "");
84  endcapfile_ = iConfig.getUntrackedParameter<std::string>("initialMiscalibrationEndcap", "");
85 
87  iConfig.getUntrackedParameter<unsigned int>("electronSelection", 0); //option for electron selection
88 
89  etaBins_ = iConfig.getUntrackedParameter<unsigned int>("etaBins", 10);
90  etBins_ = iConfig.getUntrackedParameter<unsigned int>("etBins", 10);
91 
92  etaMin_ = iConfig.getUntrackedParameter<double>("etaMin", 0.);
93  etMin_ = iConfig.getUntrackedParameter<double>("etMin", 0.);
94  etaMax_ = iConfig.getUntrackedParameter<double>("etaMax", 3.);
95  etMax_ = iConfig.getUntrackedParameter<double>("etMax", 100.);
96 
97  // new ZeePlots("zeePlots.root");
98  // ZeePlots->bookHistos();
99 
100  //ZeeCalibrationPLots("zeeCalibPlots");
101  //ZeecaPlots->bookHistos(maxsIter);
102 
103  theParameterSet = iConfig;
104  EcalIndexingTools* myIndexTool = nullptr;
105 
106  myIndexTool = EcalIndexingTools::getInstance();
107 
109 
110  //creating the algorithm
112 
113  // Tell the framework what data is being produced
114  //setWhatProduced(this);
116  findingRecord<EcalIntercalibConstantsRcd>();
117 
118  for (int i = 0; i < 50; i++) {
120  loopArray[i] = -1.;
121  sigmaArray[i] = -1.;
122  sigmaErrorArray[i] = -1.;
123  }
124 
125 #ifdef DEBUG
126  std::cout << "[ZeeCalibration] Done with the ctor" << std::endl;
127 #endif
128 }
129 
131  // if (theAlgorithm_)
132  // delete theAlgorithm_;
133 }
134 
135 //_____________________________________________________________________________
136 // Produce EcalIntercalibConstants
137 std::shared_ptr<EcalIntercalibConstants> ZeeCalibration::produceEcalIntercalibConstants(
138  const EcalIntercalibConstantsRcd& iRecord) {
139  std::cout << "@SUB=ZeeCalibration::produceEcalIntercalibConstants" << std::endl;
140  return ical;
141 }
142 
144 
145 //========================================================================
147  printStatistics();
148 
149  if (calibMode_ != "ETA_ET_MODE") {
151 
152  //Writing out calibration coefficients
153  calibXMLwriter* barrelWriter = new calibXMLwriter(EcalBarrel);
154  for (int ieta = -EBDetId::MAX_IETA; ieta <= EBDetId::MAX_IETA; ++ieta) {
155  if (ieta == 0)
156  continue;
157  for (int iphi = EBDetId::MIN_IPHI; iphi <= EBDetId::MAX_IPHI; ++iphi) {
158  // make an EBDetId since we need EBDetId::rawId() to be used as the key for the pedestals
159  if (EBDetId::validDetId(ieta, iphi)) {
160  EBDetId ebid(ieta, iphi);
161  barrelWriter->writeLine(ebid, *(ical->getMap().find(ebid.rawId())));
162  }
163  }
164  }
165 
166  calibXMLwriter* endcapWriter = new calibXMLwriter(EcalEndcap);
167  for (int iX = EEDetId::IX_MIN; iX <= EEDetId::IX_MAX; ++iX) {
168  for (int iY = EEDetId::IY_MIN; iY <= EEDetId::IY_MAX; ++iY) {
169  // make an EEDetId since we need EEDetId::rawId() to be used as the key for the pedestals
170  if (EEDetId::validDetId(iX, iY, 1)) {
171  EEDetId eeid(iX, iY, 1);
172  endcapWriter->writeLine(eeid, *(ical->getMap().find(eeid.rawId())));
173  }
174  if (EEDetId::validDetId(iX, iY, -1)) {
175  EEDetId eeid(iX, iY, -1);
176  endcapWriter->writeLine(eeid, *(ical->getMap().find(eeid.rawId())));
177  }
178  }
179  }
180 
181  }
182 
183  std::cout << "Writing histos..." << std::endl;
184  outputFile_->cd();
185 
186  // zeeplts->Write();
187 
190 
193 
195 
199 
200  h1_electronCosTheta_SC_->Write();
201  h1_electronCosTheta_TK_->Write();
203 
204  h1_zMassResol_->Write();
205  h1_zEtaResol_->Write();
206  h1_zPhiResol_->Write();
207  h1_eleEtaResol_->Write();
208  h1_elePhiResol_->Write();
209  h1_seedOverSC_->Write();
210  h1_preshowerOverSC_->Write();
211 
212  for (unsigned int i = 0; i < 25; i++) {
213  if (i < theMaxLoops) {
214  h_ESCEtrueVsEta_[i]->Write();
215  h_ESCEtrue_[i]->Write();
216 
217  h_ESCcorrEtrueVsEta_[i]->Write();
218  h_ESCcorrEtrue_[i]->Write();
219 
220  h2_chi2_[i]->Write();
221  h2_iterations_[i]->Write();
222 
223  // h_DiffZMassDistr_[i]->Write();
224 
225  //h_ZMassDistr_[i]->Write();
226  }
227  }
228 
229  h2_fEtaBarrelGood_->Write();
230  h2_fEtaBarrelBad_->Write();
231  h2_fEtaEndcapGood_->Write();
232  h2_fEtaEndcapBad_->Write();
233  h1_eleClasses_->Write();
234 
235  h_eleEffEta_[0]->Write();
236  h_eleEffPhi_[0]->Write();
237  h_eleEffPt_[0]->Write();
238 
239  h_eleEffEta_[1]->Write();
240  h_eleEffPhi_[1]->Write();
241  h_eleEffPt_[1]->Write();
242 
243  int j = 0;
244 
245  int flag = 0;
246 
247  Double_t mean[25] = {0.};
248  Double_t num[25] = {0.};
249  Double_t meanErr[25] = {0.};
250  Float_t rms[25] = {0.};
251  Float_t tempRms[10][25];
252 
253  for (int ia = 0; ia < 10; ia++) {
254  for (int ib = 0; ib < 25; ib++) {
255  tempRms[ia][ib] = 0.;
256  }
257  }
258 
259  int aa = 0;
260 
261  for (int k = 0; k < theAlgorithm_->getNumberOfChannels(); k++) {
263  bool isNearCrack = false;
264 
265  if (calibMode_ == "RING") {
266  isNearCrack = (abs(ringNumberCorrector(k)) == 1 || abs(ringNumberCorrector(k)) == 25 ||
267  abs(ringNumberCorrector(k)) == 26 || abs(ringNumberCorrector(k)) == 45 ||
268  abs(ringNumberCorrector(k)) == 46 || abs(ringNumberCorrector(k)) == 65 ||
269  abs(ringNumberCorrector(k)) == 66 || abs(ringNumberCorrector(k)) == 85 ||
270  abs(ringNumberCorrector(k)) == 86 || abs(ringNumberCorrector(k)) == 124);
271  }
272 
273  if (k < 85) {
274  if ((k + 1) % 5 != 0) {
275  if (!isNearCrack) {
276  mean[j] += calibCoeff[k];
277  mean[j] += calibCoeff[169 - k];
278 
279  num[j] += 2.;
280 
281  //meanErr[j]+= calibCoeffError[k];
282  //meanErr[j]+= calibCoeffError[169 - k];
283 
284  meanErr[j] += 1. / pow(calibCoeffError[k], 2);
285  meanErr[j] += 1. / pow(calibCoeffError[169 - k], 2);
286 
287  tempRms[aa][j] += calibCoeff[k];
288  aa++;
289  tempRms[aa][j] += calibCoeff[169 - k];
290  aa++;
291  }
292  }
293 
294  else {
295  if (!isNearCrack) {
296  mean[j] += calibCoeff[k];
297  mean[j] += calibCoeff[169 - k];
298 
299  num[j] += 2.;
300 
301  //meanErr[j]+= calibCoeffError[k];
302  //meanErr[j]+= calibCoeffError[169 - k];
303 
304  meanErr[j] += 1. / pow(calibCoeffError[k], 2);
305  meanErr[j] += 1. / pow(calibCoeffError[169 - k], 2);
306 
307  tempRms[aa][j] += calibCoeff[k];
308  aa++;
309  tempRms[aa][j] += calibCoeff[169 - k];
310  aa++;
311  }
312  j++;
313  aa = 0;
314  }
315  }
316  //EE begin
317 
318  if (k >= 170 && k <= 204) {
319  if (flag < 4) {
320  //make groups of 5 Xtals in #eta
321  mean[j] += calibCoeff[k] / 10.;
322  mean[j] += calibCoeff[k + 39] / 10.;
323 
324  meanErr[j] += calibCoeffError[k] / 30.;
325  meanErr[j] += calibCoeffError[k + 39] / 30.;
326 
327  tempRms[aa][j] += calibCoeff[k];
328  aa++;
329  tempRms[aa][j] += calibCoeff[k + 39];
330  aa++;
331 
332  flag++;
333  }
334 
335  else if (flag == 4) {
336  //make groups of 5 Xtals in #eta
337  mean[j] += calibCoeff[k] / 10.;
338  mean[j] += calibCoeff[k + 39] / 10.;
339 
340  meanErr[j] += calibCoeffError[k] / 30.;
341  meanErr[j] += calibCoeffError[k + 39] / 30.;
342 
343  tempRms[aa][j] += calibCoeff[k];
344  aa++;
345  tempRms[aa][j] += calibCoeff[k + 39];
346  aa++;
347 
348  flag = 0;
349  // std::cout<<" index(>85) "<<k<<" j is "<<j<<" mean[j] is "<<mean[j]<<std::endl;
350  j++;
351  aa = 0;
352  }
353  }
354  if (k >= 205 && k <= 208) {
355  mean[j] += calibCoeff[k] / 8.;
356  mean[j] += calibCoeff[k + 39] / 8.;
357 
358  meanErr[j] += calibCoeffError[k] / 30.;
359  meanErr[j] += calibCoeffError[k + 39] / 30.;
360 
361  tempRms[aa][j] += calibCoeff[k];
362  aa++;
363  tempRms[aa][j] += calibCoeff[k + 39];
364  aa++;
365  }
366  //EE end
367 
368  /*
369  for(int jj =0; jj< 25; jj++){
370  if(meanErr[jj] > 0.)
371  std::cout<<" meanErr[jj] before sqrt: "<<meanErr[jj]<<std::endl;
372 
373  meanErr[jj] = 1./sqrt( meanErr[jj] );
374 
375  std::cout<<" meanErr[jj] after sqrt: "<<meanErr[jj]<<std::endl;
376  }
377  */
378 
379  if (!isNearCrack) {
382  h1_mc_->Fill(initCalibCoeff[k] * calibCoeff[k] - 1.);
383 
384  if (k < 170) {
386  h1_mcEB_->Fill(initCalibCoeff[k] * calibCoeff[k] - 1.);
387  }
388 
389  if (k >= 170) {
391  h1_mcEE_->Fill(initCalibCoeff[k] * calibCoeff[k] - 1.);
392  }
393  }
394  }
395 
396  for (int ic = 0; ic < 17; ic++) {
397  mean[ic] = mean[ic] / num[ic]; //find mean of recalib coeff on group of rings
398  //meanErr[ic] = meanErr[ic] / ( sqrt( num[ic] ) * num[ic] ); //find mean of recalib coeff on group of rings
399  meanErr[ic] = 1. / sqrt(meanErr[ic]); //find mean of recalib coeff on group of rings
400  }
401 
402  //build array of RMS
403  for (int ic = 0; ic < 25; ic++) {
404  for (int id = 0; id < 10; id++) {
405  if (tempRms[id][ic] > 0.) {
406  rms[ic] += (tempRms[id][ic] - mean[j]) * (tempRms[id][ic] - mean[j]);
407  }
408  }
409  rms[ic] /= 10.; //this is approximate
410  rms[ic] = sqrt(rms[ic]);
411  }
412 
413  //build array of RMS
414 
415  Double_t xtalEta[25] = {1.4425, 1.3567, 1.2711, 1.1855, 1.10, 1.01, 0.92, 0.83, 0.7468,
416  0.6612, 0.5756, 0.4897, 0.3985, 0.3117, 0.2250, 0.1384, 0.0487, 1.546,
417  1.651, 1.771, 1.908, 2.071, 2.267, 2.516, 2.8};
418 
419  Double_t zero[25] = {0.026}; //interval/sqrt(12)
420 
421  for (int j = 0; j < 25; j++)
422  h2_coeffVsEtaGrouped_->Fill(xtalEta[j], mean[j]);
423 
424  // for(int sho = 0; sho <25; sho++)
425  //cout<<"xtalEta[j] "<< xtalEta[sho]<<" mean[j] "<<mean[sho]<<" err[j] "<<meanErr[sho]<<std::endl;
426 
427  TProfile* px = h2_coeffVsEta_->ProfileX("coeffVsEtaProfile");
428  px->SetXTitle("Eta channel");
429  px->SetYTitle("recalibCoeff");
430  px->Write();
431 
432  h2_coeffVsEta_->Write();
433  h2_coeffVsEtaGrouped_->Write();
434  h2_zMassVsLoop_->Write();
435  h2_zMassDiffVsLoop_->Write();
436  h2_zWidthVsLoop_->Write();
437  h2_coeffVsLoop_->Write();
438  h2_miscalRecal_->Write();
439  h1_mc_->Write();
440  h2_miscalRecalEB_->Write();
441  h1_mcEB_->Write();
442  h2_miscalRecalEE_->Write();
443  h1_mcEE_->Write();
444 
445  h2_residualSigma_->Write();
446 
448 
449  double weightSumMeanBarrel = 0.;
450  double weightSumMeanEndcap = 0.;
451 
452  for (int iIteration = 0; iIteration < theAlgorithm_->getNumberOfIterations(); iIteration++)
453  for (int iChannel = 0; iChannel < theAlgorithm_->getNumberOfChannels(); iChannel++) {
454  if (iIteration == (theAlgorithm_->getNumberOfIterations() - 1)) {
455  if (iChannel < 170)
456  weightSumMeanBarrel += algoHistos->weightedRescaleFactor[iIteration][iChannel]->Integral() / 170.;
457 
458  if (iChannel >= 170)
459  weightSumMeanEndcap += algoHistos->weightedRescaleFactor[iIteration][iChannel]->Integral() / 78.;
460 
461  h1_occupancyVsEta_->Fill((Double_t)ringNumberCorrector(iChannel),
462  algoHistos->weightedRescaleFactor[iIteration][iChannel]->Integral());
463 
464  h1_occupancy_->Fill(algoHistos->weightedRescaleFactor[iIteration][iChannel]->Integral());
465 
466  if (iChannel < 170)
467  h1_occupancyBarrel_->Fill(algoHistos->weightedRescaleFactor[iIteration][iChannel]->Integral());
468 
469  if (iChannel >= 170)
470  h1_occupancyEndcap_->Fill(algoHistos->weightedRescaleFactor[iIteration][iChannel]->Integral());
471 
472 #ifdef DEBUG
473  std::cout << "Writing weighted integral for channel " << ringNumberCorrector(iChannel) << " ,value "
474  << algoHistos->weightedRescaleFactor[iIteration][iChannel]->Integral() << std::endl;
475 #endif
476  }
477  }
478 
479  // std::cout<<"Done! Closing output file... "<<std::endl;
480 
481  h1_weightSumMeanBarrel_->Fill(weightSumMeanBarrel);
482  h1_weightSumMeanEndcap_->Fill(weightSumMeanEndcap);
483 
484  std::cout << "Weight sum mean on channels in Barrel is :" << weightSumMeanBarrel << std::endl;
485  std::cout << "Weight sum mean on channels in Endcap is :" << weightSumMeanEndcap << std::endl;
486 
487  h1_weightSumMeanBarrel_->Write();
488  h1_weightSumMeanEndcap_->Write();
489 
490  h1_occupancyVsEta_->Write();
491  h1_occupancy_->Write();
492  h1_occupancyBarrel_->Write();
493  h1_occupancyEndcap_->Write();
494 
495  myTree->Write();
496 
497  TGraphErrors* graph = new TGraphErrors(25, xtalEta, mean, zero, meanErr);
498  graph->Draw("APL");
499  graph->Write();
500 
501  double zero50[50] = {0.};
502 
503  TGraphErrors* residualSigmaGraph = new TGraphErrors(50, loopArray, sigmaArray, zero50, sigmaErrorArray);
504  residualSigmaGraph->SetName("residualSigmaGraph");
505  residualSigmaGraph->Draw("APL");
506  residualSigmaGraph->Write();
507 
508  TGraphErrors* coefficientDistanceAtIterationGraph =
509  new TGraphErrors(50, loopArray, coefficientDistanceAtIteration, zero50, zero50);
510  coefficientDistanceAtIterationGraph->SetName("coefficientDistanceAtIterationGraph");
511  coefficientDistanceAtIterationGraph->Draw("APL");
512  coefficientDistanceAtIterationGraph->Write();
513 
514  Float_t noError[250] = {0.};
515 
516  Float_t ringInd[250];
517  for (int i = 0; i < 250; i++)
518  ringInd[i] = ringNumberCorrector(i);
519 
520  TGraphErrors* graphCoeff =
521  new TGraphErrors(theAlgorithm_->getNumberOfChannels(), ringInd, calibCoeff, noError, calibCoeffError);
522  graphCoeff->SetName("graphCoeff");
523  graphCoeff->Draw("APL");
524  graphCoeff->Write();
525 
526  // outputFile_->Write();//this automatically writes all histos on file
527 
528  h1_ZCandMult_->Write();
529  h1_reco_ZMass_->Write();
530 
531  h1_reco_ZMassCorr_->Write();
532  h1_reco_ZMassCorrBB_->Write();
533  h1_reco_ZMassCorrEE_->Write();
534 
535  outputFile_->Close();
536 
541 
542  // myZeeRescaleFactorPlots_ = new ZeeRescaleFactorPlots("zeeRescaleFactorPlots.root");
543  //myZeeRescaleFactorPlots_->writeHistograms( theAlgorithm_ );
544 
545  // delete myZeeRescaleFactorPlots_;
546 }
547 
548 //_____________________________________________________________________________
549 // Called at each event
550 //________________________________________
551 
553  using namespace edm;
554 
555 #ifdef DEBUG
556  std::cout << "[ZeeCalibration] Entering duringLoop" << std::endl;
557 #endif
558 
559  // code that used to be in beginJob
560  if (isfirstcall_) {
561  //inizializzare la geometria di ecal
562  const auto& geometry = iSetup.getData(geometryToken_);
564 
565  myZeePlots_ = new ZeePlots("zeePlots.root");
566  // myZeeRescaleFactorPlots_ = new ZeeRescaleFactorPlots("zeeRescaleFactorPlots.root");
567 
568  // go to *OUR* rootfile and book histograms
569  outputFile_->cd();
570  bookHistograms();
571 
572  std::cout << "[ZeeCalibration::beginOfJob] Histograms booked " << std::endl;
573 
574  loopFlag_ = 0;
575 
576  //Read miscalibration map if requested
577  CaloMiscalibMapEcal* miscalibMap = nullptr;
578  if (!barrelfile_.empty() || !barrelfile_.empty()) {
579  miscalibMap = new CaloMiscalibMapEcal();
580  miscalibMap->prefillMap();
581  }
582 
583  if (!barrelfile_.empty()) {
584  MiscalibReaderFromXMLEcalBarrel barrelreader_(*miscalibMap);
585  barrelreader_.parseXMLMiscalibFile(barrelfile_);
586 #ifdef DEBUG
587  std::cout << "[ZeeCalibration::beginOfJob] Parsed EB miscal file" << std::endl;
588 #endif
589  }
590 
591  if (!endcapfile_.empty()) {
592  MiscalibReaderFromXMLEcalEndcap endcapreader_(*miscalibMap);
593  endcapreader_.parseXMLMiscalibFile(endcapfile_);
594 #ifdef DEBUG
595  std::cout << "[ZeeCalibration::beginOfJob] Parsed EE miscal file" << std::endl;
596 #endif
597  }
598 
599  std::cout << " theAlgorithm_->getNumberOfChannels() " << theAlgorithm_->getNumberOfChannels() << std::endl;
600 
602  for (int k = 0; k < theAlgorithm_->getNumberOfChannels(); k++) {
603  calibCoeff[k] = 1.;
604  calibCoeffError[k] = 0.;
605 
606  std::vector<DetId> ringIds;
607 
608  if (calibMode_ == "RING")
610 
611  if (calibMode_ == "MODULE")
613 
614  if (calibMode_ == "ABS_SCALE" || calibMode_ == "ETA_ET_MODE")
616 
617  if (miscalibMap) {
618  initCalibCoeff[k] = 0.;
619  for (unsigned int iid = 0; iid < ringIds.size(); ++iid) {
620  float miscalib = *(miscalibMap->get().getMap().find(ringIds[iid]));
621  // float miscalib=miscalibMap->get().getMap().find(ringIds[iid])->second; ////////AP
622  initCalibCoeff[k] += miscalib;
623  }
624  initCalibCoeff[k] /= (float)ringIds.size();
625  std::cout << k << " " << initCalibCoeff[k] << " " << ringIds.size() << std::endl;
626  } else {
627  initCalibCoeff[k] = 1.;
628  }
629  }
630 
631  ical = std::make_shared<EcalIntercalibConstants>();
632 
633  for (int k = 0; k < theAlgorithm_->getNumberOfChannels(); k++) {
634  // std::vector<DetId> ringIds = EcalRingCalibrationTools::getDetIdsInRing(k);
635 
636  std::vector<DetId> ringIds;
637 
638  if (calibMode_ == "RING")
640 
641  if (calibMode_ == "MODULE")
643 
644  if (calibMode_ == "ABS_SCALE" || calibMode_ == "ETA_ET_MODE")
646 
647  for (unsigned int iid = 0; iid < ringIds.size(); ++iid) {
648  // ical->setValue( ringIds[iid], 1. * initCalibCoeff[k] );
649 
650  if (ringIds[iid].subdetId() == EcalBarrel) {
651  EBDetId myEBDetId(ringIds[iid]);
652  h2_xtalMiscalibCoeffBarrel_->SetBinContent(
653  myEBDetId.ieta() + 86,
654  myEBDetId.iphi(),
655  *(miscalibMap->get().getMap().find(ringIds[iid]))); //fill TH2 with miscalibCoeff
656  }
657 
658  if (ringIds[iid].subdetId() == EcalEndcap) {
659  EEDetId myEEDetId(ringIds[iid]);
660  if (myEEDetId.zside() < 0)
661  h2_xtalMiscalibCoeffEndcapMinus_->SetBinContent(
662  myEEDetId.ix(),
663  myEEDetId.iy(),
664  *(miscalibMap->get().getMap().find(ringIds[iid]))); //fill TH2 with miscalibCoeff
665 
666  if (myEEDetId.zside() > 0)
667  h2_xtalMiscalibCoeffEndcapPlus_->SetBinContent(
668  myEEDetId.ix(),
669  myEEDetId.iy(),
670  *(miscalibMap->get().getMap().find(ringIds[iid]))); //fill TH2 with miscalibCoeff
671  }
672 
673  ical->setValue(ringIds[iid], *(miscalibMap->get().getMap().find(ringIds[iid])));
674  }
675 
676  read_events = 0;
677  init_ = false;
678  }
679  isfirstcall_ = false;
680  } // if isfirstcall
681 
683 
684  for (unsigned int iHLT = 0; iHLT < 200; ++iHLT) {
685  aHLTResults[iHLT] = false;
686  }
687 
688 #ifdef DEBUG
689  std::cout << "[ZeeCalibration::duringLoop] Done with initializing aHLTresults[] " << std::endl;
690 #endif
691 
692  edm::Handle<edm::TriggerResults> hltTriggerResultHandle;
693  iEvent.getByToken(trigResultsToken_, hltTriggerResultHandle);
694 
695  if (!hltTriggerResultHandle.isValid()) {
696  //std::cout << "invalid handle for HLT TriggerResults" << std::endl;
697  } else {
698  hltCount = hltTriggerResultHandle->size();
699 
700  if (loopFlag_ == 0)
701  myZeePlots_->fillHLTInfo(hltTriggerResultHandle);
702 
703 #ifdef DEBUG
704  std::cout << "[ZeeCalibration::duringLoop] Done with myZeePlots_->fillHLTInfo(hltTriggerResultHandle); "
705  << std::endl;
706 #endif
707 
708  for (int i = 0; i < hltCount; i++) {
709  aHLTResults[i] = hltTriggerResultHandle->accept(i);
710 
711  //HLT bit 32 = HLT1Electron
712  //HLT bit 34 = HLT2Electron
713  //HLT bit 35 = HLT2ElectronRelaxed
714  }
715 
716  if (!aHLTResults[32] && !aHLTResults[34] && !aHLTResults[35])
717  return kContinue;
718  }
719 
720 #ifdef DEBUG
721  std::cout << "[ZeeCalibration::duringLoop] End HLT section" << std::endl;
722 #endif
723 
725 
726  std::vector<HepMC::GenParticle*> mcEle;
727 
728  float myGenZMass(-1);
729 
730  if (!mcProducer_.empty()) {
731  //DUMP GENERATED Z MASS - BEGIN
732  Handle<HepMCProduct> hepProd;
733  iEvent.getByToken(hepMCToken_, hepProd);
734 
735  const HepMC::GenEvent* myGenEvent = hepProd->GetEvent();
736 
737  if (loopFlag_ == 0)
738  myZeePlots_->fillZMCInfo(&(*myGenEvent));
739 
740 #ifdef DEBUG
741  std::cout << "[ZeeCalibration::duringLoop] Done with myZeePlots_->fillZMCInfo( & (*myGenEvent) ); " << std::endl;
742 #endif
743 
744  for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
745  ++p) {
746  //return a pointer to MC Z in the event
747  if ((*p)->pdg_id() == 23 && (*p)->status() == 2) {
748  myGenZMass = (*p)->momentum().m();
749  }
750  }
751  //DUMP GENERATED Z MASS - END
752 
753  if (loopFlag_ == 0)
754  myZeePlots_->fillEleMCInfo(&(*myGenEvent));
755 
756  //loop over MC positrons and find the closest MC positron in (eta,phi) phace space - begin
757  HepMC::GenParticle MCele;
758 
759  for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
760  ++p) {
761  if (abs((*p)->pdg_id()) == 11) {
762  mcEle.push_back((*p));
763  MCele = *(*p);
764  }
765  }
766 
767  if (mcEle.size() == 2 && fabs(mcEle[0]->momentum().eta()) < 2.4 && fabs(mcEle[1]->momentum().eta()) < 2.4) {
768  NEVT++;
769 
770  if (fabs(mcEle[0]->momentum().eta()) < 1.479 && fabs(mcEle[1]->momentum().eta()) < 1.479)
771  MCZBB++;
772 
773  if ((fabs(mcEle[0]->momentum().eta()) > 1.479 && fabs(mcEle[1]->momentum().eta()) < 1.479) ||
774  (fabs(mcEle[0]->momentum().eta()) < 1.479 && fabs(mcEle[1]->momentum().eta()) > 1.479))
775  MCZEB++;
776 
777  if (fabs(mcEle[0]->momentum().eta()) > 1.479 && fabs(mcEle[1]->momentum().eta()) > 1.479)
778  MCZEE++;
779  }
780  }
781 
782  // Get EBRecHits
784  iEvent.getByToken(ebRecHitToken_, phits);
785  if (!phits.isValid()) {
786  std::cerr << "Error! can't get the product EBRecHitCollection " << std::endl;
787  }
788  const EBRecHitCollection* hits = phits.product(); // get a ptr to the product
789 
790  // Get EERecHits
792  iEvent.getByToken(eeRecHitToken_, ephits);
793  if (!ephits.isValid()) {
794  std::cerr << "Error! can't get the product EERecHitCollection " << std::endl;
795  }
796  const EERecHitCollection* ehits = ephits.product(); // get a ptr to the product
797 
798  //Get Hybrid SuperClusters
800  iEvent.getByToken(scToken_, pSuperClusters);
801  if (!pSuperClusters.isValid()) {
802  std::cerr << "Error! can't get the product SuperClusterCollection " << std::endl;
803  }
804  const reco::SuperClusterCollection* scCollection = pSuperClusters.product();
805 
806 #ifdef DEBUG
807  std::cout << "scCollection->size()" << scCollection->size() << std::endl;
808  for (reco::SuperClusterCollection::const_iterator scIt = scCollection->begin(); scIt != scCollection->end(); scIt++) {
809  std::cout << scIt->energy() << std::endl;
810  }
811 #endif
812 
813  //Get Island SuperClusters
814  Handle<reco::SuperClusterCollection> pIslandSuperClusters;
815  iEvent.getByToken(islandSCToken_, pIslandSuperClusters);
816  if (!pIslandSuperClusters.isValid()) {
817  std::cerr << "Error! can't get the product IslandSuperClusterCollection " << std::endl;
818  }
819  const reco::SuperClusterCollection* scIslandCollection = pIslandSuperClusters.product();
820 
821 #ifdef DEBUG
822  std::cout << "scCollection->size()" << scIslandCollection->size() << std::endl;
823 #endif
824 
825  if ((scCollection->size() + scIslandCollection->size()) < 2)
826  return kContinue;
827 
828  // Get Electrons
830  iEvent.getByToken(gsfElectronToken_, pElectrons);
831  if (!pElectrons.isValid()) {
832  std::cerr << "Error! can't get the product ElectronCollection " << std::endl;
833  }
835 
836  /*
837  //reco-mc association map
838  std::map<HepMC::GenParticle*,const reco::PixelMatchGsfElectron*> myMCmap;
839 
840  fillMCmap(&(*electronCollection),mcEle,myMCmap);
841 
842  fillEleInfo(mcEle,myMCmap);
843  */
844 
845  if (electronCollection->size() < 2)
846  return kContinue;
847 
848  if (!hits && !ehits) {
849  std::cout << "!hits" << std::endl;
850  return kContinue;
851  }
852 
853  if (hits->empty() && ehits->empty()) {
854  std::cout << "hits->size() == 0" << std::endl;
855  return kContinue;
856  }
857 
858  if (!electronCollection) {
859  std::cout << "!electronCollection" << std::endl;
860  return kContinue;
861  }
862 
863  if (electronCollection->empty()) {
864  std::cout << "electronCollection->size() == 0" << std::endl;
865  return kContinue;
866  }
867 
871 
872  read_events++;
873 
874  // std::cout << "read_events = " << read_events << std::endl;
875 
877 
878 #ifdef DEBUG
879  std::cout << " Starting with myZeePlots_->fillEleInfo(electronCollection); " << std::endl;
880 #endif
881 
882  if (loopFlag_ == 0)
884 
885 #ifdef DEBUG
886  std::cout << " Done with myZeePlots_->fillEleInfo(electronCollection); " << std::endl;
887 #endif
888 
889  //FILL an electron vector - end
890  //###################################Electron-SC association: begin#####################################################
891  //Filling new ElectronCollection with new SC ref and calibElectron container
892  std::vector<calib::CalibElectron> calibElectrons;
893  //std::map< const calib::CalibElectron* , const reco::SuperCluster* > eleScMap;
894 
895  //#####################################Electron-SC association map: end#####################################################
896  for (unsigned int e_it = 0; e_it != electronCollection->size(); e_it++) {
897  calibElectrons.push_back(calib::CalibElectron(&((*electronCollection)[e_it]), hits, ehits));
898 #ifdef DEBUG
899  std::cout << calibElectrons.back().getRecoElectron()->superCluster()->energy() << " "
900  << calibElectrons.back().getRecoElectron()->energy() << std::endl;
901 #endif
902  // h1_recoEleEnergy_->Fill(calibElectrons.back().getRecoElectron()->superCluster()->energy());
903  }
904  // if (iLoop == 0)
905  //fillCalibElectrons(calibElectrons);
906 
907 #ifdef DEBUG
908  std::cout << "Filled histos" << std::endl;
909 #endif
910 
911  //COMBINATORY FOR Z MASS - begin
912  std::vector<std::pair<calib::CalibElectron*, calib::CalibElectron*> > zeeCandidates;
913  int myBestZ = -1;
914 
915  mass = -1.;
916  double DeltaMinvMin(5000.);
917 
918  if (calibElectrons.size() < 2)
919  return kContinue;
920 
921  for (unsigned int e_it = 0; e_it != calibElectrons.size() - 1; e_it++) {
922  for (unsigned int p_it = e_it + 1; p_it != calibElectrons.size(); p_it++) {
923 #ifdef DEBUG
924  std::cout << e_it << " " << calibElectrons[e_it].getRecoElectron()->charge() << " " << p_it << " "
925  << calibElectrons[p_it].getRecoElectron()->charge() << std::endl;
926 #endif
927  if (calibElectrons[e_it].getRecoElectron()->charge() * calibElectrons[p_it].getRecoElectron()->charge() != -1)
928  continue;
929 
931  std::pair<calib::CalibElectron*, calib::CalibElectron*>(&(calibElectrons[e_it]), &(calibElectrons[p_it])));
932 
933  if (mass < 0)
934  continue;
935 
936 #ifdef DEBUG
937  std::cout << "#######################mass " << mass << std::endl;
938 #endif
939 
940  zeeCandidates.push_back(
941  std::pair<calib::CalibElectron*, calib::CalibElectron*>(&(calibElectrons[e_it]), &(calibElectrons[p_it])));
942  double DeltaMinv = fabs(mass - MZ);
943 
944  if (DeltaMinv < DeltaMinvMin) {
945  DeltaMinvMin = DeltaMinv;
946  myBestZ = zeeCandidates.size() - 1;
947  }
948  }
949  }
950 
951  // h_DeltaZMassDistr_[loopFlag_]->Fill( (mass-MZ) / MZ );
952 
953  // zeeCa->Fill(zeeCandidates);
954  //
955  h1_ZCandMult_->Fill(zeeCandidates.size());
956 
957  if (zeeCandidates.empty() || myBestZ == -1)
958  return kContinue;
959 
960  if (loopFlag_ == 0)
961  myZeePlots_->fillZInfo(zeeCandidates[myBestZ]);
962 
963 #ifdef DEBUG
964  std::cout << "Found ZCandidates " << myBestZ << std::endl;
965 #endif
966 
967  // h1_zMassResol_ ->Fill(mass-myGenZMass);
968 
970 
971  int class1 = zeeCandidates[myBestZ].first->getRecoElectron()->classification();
972  int class2 = zeeCandidates[myBestZ].second->getRecoElectron()->classification();
973 
974  h1_eleClasses_->Fill(class1);
975  h1_eleClasses_->Fill(class2);
976 
977  std::cout << "BEFORE " << std::endl;
978 
979  // myZeePlots_->fillEleClassesPlots( zeeCandidates[myBestZ].first );
980  //myZeePlots_->fillEleClassesPlots( zeeCandidates[myBestZ].second );
981 
982  std::cout << "AFTER " << std::endl;
983 
985 
986  if (class1 < 100)
987  // h1_Elec_->Fill(1);
989  if (class1 >= 100)
991 
992  if (class2 < 100)
994  if (class2 >= 100)
996 
997  if (class1 == 0)
999  if (class1 == 100)
1001  if (class1 == 10 || class1 == 20)
1003  if (class1 == 110 || class1 == 120)
1005  if (class1 >= 30 && class1 <= 34)
1007  if (class1 >= 130 && class1 <= 134)
1009  if (class1 == 40)
1011  if (class1 == 140)
1013 
1014  if (class2 == 0)
1016  if (class2 == 100)
1018  if (class2 == 10 || class2 == 20)
1020  if (class2 == 110 || class2 == 120)
1022  if (class2 >= 30 && class2 <= 34)
1024  if (class2 >= 130 && class2 <= 134)
1026  if (class2 == 40)
1028  if (class2 == 140)
1030 
1032 
1034 
1035  DetId firstElehottestDetId =
1037  zeeCandidates[myBestZ].first->getRecoElectron()->superCluster()->seed()->hitsAndFractions(), hits, ehits)
1038  .first;
1039  DetId secondElehottestDetId =
1041  zeeCandidates[myBestZ].second->getRecoElectron()->superCluster()->seed()->hitsAndFractions(), hits, ehits)
1042  .first;
1043 
1044  bool firstElectronIsOnModuleBorder(false);
1045  bool secondElectronIsOnModuleBorder(false);
1046 
1048 
1049  if (class1 < 100) {
1050  if (firstElehottestDetId.subdetId() == EcalBarrel)
1051  firstElectronIsOnModuleBorder = xtalIsOnModuleBorder(firstElehottestDetId);
1052 
1054 
1055  if (firstElehottestDetId.subdetId() == EcalBarrel && !firstElectronIsOnModuleBorder)
1057  }
1058 
1059  if (class2 < 100) {
1060  if (secondElehottestDetId.subdetId() == EcalBarrel)
1061  secondElectronIsOnModuleBorder = xtalIsOnModuleBorder(secondElehottestDetId);
1062 
1064 
1065  if (secondElehottestDetId.subdetId() == EcalBarrel && !secondElectronIsOnModuleBorder)
1067  }
1068 
1069  if (class1 < 100) {
1070  if (firstElehottestDetId.subdetId() == EcalBarrel && firstElectronIsOnModuleBorder) {
1071  h1_borderElectronClassification_->Fill(class1);
1072  return kContinue;
1073  }
1074  }
1075 
1076  if (class2 < 100) {
1077  if (secondElehottestDetId.subdetId() == EcalBarrel && secondElectronIsOnModuleBorder) {
1078  h1_borderElectronClassification_->Fill(class2);
1079  return kContinue;
1080  }
1081  }
1082 
1085 
1086  if (class1 < 100 && class2 < 100) {
1087  BBZN++;
1088  if (class1 == 0 && class2 == 0)
1089  BBZN_gg++;
1090  if (class1 < 21 && class2 < 21)
1091  BBZN_tt++;
1092  if (class1 < 21 || class2 < 21)
1093  BBZN_t0++;
1094  }
1095 
1096  if (class1 >= 100 && class2 >= 100) {
1097  EEZN++;
1098  if (class1 == 100 && class2 == 100)
1099  EEZN_gg++;
1100  if (class1 < 121 && class2 < 121)
1101  EEZN_tt++;
1102  if (class1 < 121 || class2 < 121)
1103  EEZN_t0++;
1104  }
1105 
1106  if ((class1 < 100 && class2 >= 100) || (class2 < 100 && class1 >= 100)) {
1107  EBZN++;
1108  if ((class1 == 0 && class2 == 100) || (class2 == 0 && class1 == 100))
1109  EBZN_gg++;
1110  if ((class1 < 21 && class2 < 121) || (class2 < 21 && class1 < 121))
1111  EBZN_tt++;
1112  if (class2 < 21 || class1 < 21 || class2 < 121 || class1 < 121)
1113  EBZN_t0++;
1114  }
1115 
1117 
1118  if (myBestZ == -1)
1119  return kContinue;
1120 
1121  bool invMassBool = ((mass > minInvMassCut_) && (mass < maxInvMassCut_));
1122 
1123  bool selectionBool = false;
1124  //0 = all electrons (but no crack)
1125 
1127 
1128  float theta1 = 2. * atan(exp(-zeeCandidates[myBestZ].first->getRecoElectron()->superCluster()->eta()));
1129  bool ET_1 = ((zeeCandidates[myBestZ].first->getRecoElectron()->superCluster()->energy() * sin(theta1)) > 20.);
1130 
1131  float theta2 = 2. * atan(exp(-zeeCandidates[myBestZ].second->getRecoElectron()->superCluster()->eta()));
1132  bool ET_2 = ((zeeCandidates[myBestZ].second->getRecoElectron()->superCluster()->energy() * sin(theta2)) > 20.);
1133 
1134  bool HoE_1 = (zeeCandidates[myBestZ].first->getRecoElectron()->hadronicOverEm() < 0.115);
1135  bool HoE_2 = (zeeCandidates[myBestZ].second->getRecoElectron()->hadronicOverEm() < 0.115);
1136 
1137  bool DeltaPhiIn_1 = (zeeCandidates[myBestZ].first->getRecoElectron()->deltaPhiSuperClusterTrackAtVtx() < 0.090);
1138  bool DeltaPhiIn_2 = (zeeCandidates[myBestZ].second->getRecoElectron()->deltaPhiSuperClusterTrackAtVtx() < 0.090);
1139 
1140  bool DeltaEtaIn_1 = (zeeCandidates[myBestZ].first->getRecoElectron()->deltaEtaSuperClusterTrackAtVtx() < 0.0090);
1141  bool DeltaEtaIn_2 = (zeeCandidates[myBestZ].second->getRecoElectron()->deltaEtaSuperClusterTrackAtVtx() < 0.0090);
1142 
1144 
1145  if (!(invMassBool && ET_1 && ET_2 && HoE_1 && HoE_2 && DeltaPhiIn_1 && DeltaPhiIn_2 && DeltaEtaIn_1 && DeltaEtaIn_2))
1146  return kContinue;
1148 
1149  h1_eventsAfterEWKSelection_->Fill(1);
1150 
1152 
1153  int c1st = zeeCandidates[myBestZ].first->getRecoElectron()->classification();
1154  int c2nd = zeeCandidates[myBestZ].second->getRecoElectron()->classification();
1155  if (electronSelection_ == 0)
1156  selectionBool = (myBestZ != -1 && c1st != 40 && c1st != 40 && c2nd != 40 && c2nd != 140);
1157 
1158  //1 = all electrons are Golden, BB or Narrow
1159 
1160  if (electronSelection_ == 1)
1161  selectionBool =
1162  (myBestZ != -1 && (c1st == 0 || c1st == 10 || c1st == 20 || c1st == 100 || c1st == 110 || c1st == 120) &&
1163  (c2nd == 0 || c2nd == 10 || c2nd == 20 || c2nd == 100 || c2nd == 110 || c2nd == 120));
1164 
1165  //2 = all electrons are Golden
1166  if (electronSelection_ == 2)
1167  selectionBool = (myBestZ != -1 && (c1st == 0 || c1st == 100) && (c2nd == 0 || c2nd == 100));
1168  //3 = all electrons are showering
1169  if (electronSelection_ == 3)
1170  selectionBool = (myBestZ != -1 && ((c1st >= 30 && c1st <= 34) || ((c1st >= 130 && c1st <= 134))) &&
1171  ((c2nd >= 30 && c2nd <= 34) || ((c2nd >= 130 && c2nd <= 134)))
1172 
1173  );
1174 
1175  //4 = all Barrel electrons are Golden, BB or Narrow; take all Endcap electrons
1176 
1177  if (electronSelection_ == 4)
1178  selectionBool = (myBestZ != -1 && (
1179 
1180  ((c1st == 0 || c1st == 10 || c1st == 20) && c2nd >= 100 && c2nd != 140)
1181 
1182  ||
1183 
1184  ((c2nd == 0 || c2nd == 10 || c2nd == 20) && c1st >= 100 && c1st != 140)
1185 
1186  ));
1187 
1188  //5 = all Endcap electrons (but no crack)
1189 
1190  if (electronSelection_ == 5)
1191  selectionBool = (myBestZ != -1 && c1st >= 100 && c2nd >= 100 && c1st != 140 && c2nd != 140);
1192 
1193  //6 = all Barrel electrons (but no crack)
1194 
1195  if (electronSelection_ == 6)
1196  selectionBool = (myBestZ != -1 && c1st < 100 && c2nd < 100 && c1st != 40 && c2nd != 40);
1197 
1198  //7 = this eliminates the events which have 1 ele in the Barrel and 1 in the Endcap
1199 
1200  if (electronSelection_ == 7)
1201  selectionBool = (myBestZ != -1 && !(c1st < 100 && c2nd >= 100) && !(c1st >= 100 && c2nd < 100));
1202 
1203  float ele1EnergyCorrection(1.);
1204  float ele2EnergyCorrection(1.);
1205 
1206  if (invMassBool && selectionBool && wantEtaCorrection_) {
1207  ele1EnergyCorrection = getEtaCorrection(zeeCandidates[myBestZ].first->getRecoElectron());
1208  ele2EnergyCorrection = getEtaCorrection(zeeCandidates[myBestZ].second->getRecoElectron());
1209  }
1210 
1211  if (invMassBool && selectionBool) {
1213  ZeeKinematicTools::cosThetaElectrons_SC(zeeCandidates[myBestZ], ele1EnergyCorrection, ele2EnergyCorrection));
1215  ZeeKinematicTools::cosThetaElectrons_TK(zeeCandidates[myBestZ], ele1EnergyCorrection, ele2EnergyCorrection));
1217  ZeeKinematicTools::cosThetaElectrons_SC(zeeCandidates[myBestZ], ele1EnergyCorrection, ele2EnergyCorrection) /
1219  zeeCandidates[myBestZ], ele1EnergyCorrection, ele2EnergyCorrection) -
1220  1.);
1221 
1222  if (!mcProducer_.empty()) {
1223  h1_zMassResol_->Fill(mass - myGenZMass);
1224 
1225  //reco-mc association map - begin
1226 
1227  std::map<HepMC::GenParticle*, const reco::GsfElectron*> myMCmap;
1228 
1229  std::vector<const reco::GsfElectron*> dauElectronCollection;
1230 
1231  dauElectronCollection.push_back(zeeCandidates[myBestZ].first->getRecoElectron());
1232  dauElectronCollection.push_back(zeeCandidates[myBestZ].second->getRecoElectron());
1233 
1234  fillMCmap(&dauElectronCollection, mcEle, myMCmap);
1235  fillEleInfo(mcEle, myMCmap);
1236  //h_DiffZMassDistr_[loopFlag_]->Fill( (mass-myGenZMass) );
1237  }
1238 
1239  //PUT f(eta) IN OUR Zee ALGORITHM
1240  theAlgorithm_->addEvent(zeeCandidates[myBestZ].first,
1241  zeeCandidates[myBestZ].second,
1242  MZ * sqrt(ele1EnergyCorrection * ele2EnergyCorrection));
1243 
1244  h1_reco_ZMass_->Fill(ZeeKinematicTools::calculateZMass_withTK(zeeCandidates[myBestZ]));
1245 
1247  zeeCandidates[myBestZ], ele1EnergyCorrection, ele2EnergyCorrection));
1248 
1249  int c1st = zeeCandidates[myBestZ].first->getRecoElectron()->classification();
1250  int c2nd = zeeCandidates[myBestZ].second->getRecoElectron()->classification();
1251  if (c1st < 100 && c2nd < 100)
1253  zeeCandidates[myBestZ], ele1EnergyCorrection, ele2EnergyCorrection));
1254 
1255  if (c1st >= 100 && c2nd >= 100)
1257  zeeCandidates[myBestZ], ele1EnergyCorrection, ele2EnergyCorrection));
1258 
1260  zeeCandidates[myBestZ], ele1EnergyCorrection, ele2EnergyCorrection);
1261 
1263  zeeCandidates[myBestZ], ele1EnergyCorrection, ele2EnergyCorrection) -
1264  myGenZMass;
1265 
1266  // h_ZMassDistr_[loopFlag_]->Fill(ZeeKinematicTools::calculateZMassWithCorrectedElectrons_withTK(zeeCandidates[myBestZ],ele1EnergyCorrection,ele2EnergyCorrection));
1267 
1268  myTree->Fill();
1269  }
1270 
1271 #ifdef DEBUG
1272  std::cout << "Added event to algorithm" << std::endl;
1273 #endif
1274 
1275  return kContinue;
1276 }
1277 //end of ZeeCalibration::duringLoop
1278 
1279 // Called at beginning of loop
1280 void ZeeCalibration::startingNewLoop(unsigned int iLoop) {
1281  std::cout << "[ZeeCalibration] Starting loop number " << iLoop << std::endl;
1282 
1284 
1285  resetVariables();
1286 
1287  resetHistograms();
1288 
1289 #ifdef DEBUG
1290  std::cout << "[ZeeCalibration] exiting from startingNewLoop" << std::endl;
1291 #endif
1292 }
1293 
1294 // Called at end of loop
1296  double par[3];
1297  double errpar[3];
1298  double zChi2;
1299  int zIters;
1300 
1301  ZIterativeAlgorithmWithFit::gausfit(h1_reco_ZMass_, par, errpar, 2., 2., &zChi2, &zIters);
1302 
1303  h2_zMassVsLoop_->Fill(loopFlag_, par[1]);
1304 
1305  h2_zMassDiffVsLoop_->Fill(loopFlag_, (par[1] - MZ) / MZ);
1306 
1307  h2_zWidthVsLoop_->Fill(loopFlag_, par[2]);
1308 
1310 
1311  std::cout << "[ZeeCalibration] Ending loop " << iLoop << std::endl;
1312  //RUN the algorithm
1314 
1315  const std::vector<float>& optimizedCoefficients = theAlgorithm_->getOptimizedCoefficients();
1316  const std::vector<float>& optimizedCoefficientsError = theAlgorithm_->getOptimizedCoefficientsError();
1317  //const std::vector<float>& weightSum = theAlgorithm_->getWeightSum();
1318  const std::vector<float>& optimizedChi2 = theAlgorithm_->getOptimizedChiSquare();
1319  const std::vector<int>& optimizedIterations = theAlgorithm_->getOptimizedIterations();
1320 
1321  //#ifdef DEBUG
1322  std::cout << "Optimized coefficients " << optimizedCoefficients.size() << std::endl;
1323  //#endif
1324 
1325  // h2_coeffVsLoop_->Fill(loopFlag_, optimizedCoefficients[75]); //show the evolution of just 1 ring coefficient (well chosen...)
1326 
1328  for (unsigned int ieta = 0; ieta < optimizedCoefficients.size(); ieta++) {
1329  NewCalibCoeff[ieta] = calibCoeff[ieta] * optimizedCoefficients[ieta];
1330 
1331  h2_chi2_[loopFlag_]->Fill(ringNumberCorrector(ieta), optimizedChi2[ieta]);
1332  h2_iterations_[loopFlag_]->Fill(ringNumberCorrector(ieta), optimizedIterations[ieta]);
1333  }
1335 
1337  computeCoefficientDistanceAtIteration(calibCoeff, NewCalibCoeff, optimizedCoefficients.size());
1338 
1339  std::cout << "Iteration # : " << loopFlag_ << " CoefficientDistanceAtIteration "
1340  << coefficientDistanceAtIteration[loopFlag_] << std::endl;
1341  std::cout << "size " << optimizedCoefficients.size() << std::endl;
1342 
1343  for (unsigned int ieta = 0; ieta < optimizedCoefficients.size(); ieta++) {
1344  calibCoeff[ieta] *= optimizedCoefficients[ieta];
1346  calibCoeff[ieta] * sqrt(pow(optimizedCoefficientsError[ieta] / optimizedCoefficients[ieta], 2) +
1348  //calibCoeffError[ieta] = optimizedCoefficientsError[ieta];
1349 
1350 #ifdef DEBUG
1351  std::cout << ieta << " " << optimizedCoefficients[ieta] << std::endl;
1352 #endif
1353 
1354  std::vector<DetId> ringIds;
1355 
1356  if (calibMode_ == "RING")
1358 
1359  if (calibMode_ == "MODULE")
1361 
1362  if (calibMode_ == "ABS_SCALE" || calibMode_ == "ETA_ET_MODE")
1364 
1365  for (unsigned int iid = 0; iid < ringIds.size(); ++iid) {
1366  if (ringIds[iid].subdetId() == EcalBarrel) {
1367  EBDetId myEBDetId(ringIds[iid]);
1368  h2_xtalRecalibCoeffBarrel_[loopFlag_]->SetBinContent(
1369  myEBDetId.ieta() + 86,
1370  myEBDetId.iphi(),
1371  100 * (calibCoeff[ieta] * initCalibCoeff[ieta] - 1.)); //fill TH2 with recalibCoeff
1372  }
1373 
1374  if (ringIds[iid].subdetId() == EcalEndcap) {
1375  EEDetId myEEDetId(ringIds[iid]);
1376  if (myEEDetId.zside() < 0)
1378  myEEDetId.ix(),
1379  myEEDetId.iy(),
1380  100 * (calibCoeff[ieta] * initCalibCoeff[ieta] - 1.)); //fill TH2 with recalibCoeff
1381 
1382  if (myEEDetId.zside() > 0)
1383  h2_xtalRecalibCoeffEndcapPlus_[loopFlag_]->SetBinContent(
1384  myEEDetId.ix(),
1385  myEEDetId.iy(),
1386  100 * (calibCoeff[ieta] * initCalibCoeff[ieta] - 1.)); //fill TH2 with recalibCoeff
1387  }
1388 
1389  ical->setValue(ringIds[iid], *(ical->getMap().find(ringIds[iid])) * optimizedCoefficients[ieta]);
1390  }
1391  }
1392 
1394 
1395  for (int k = 0; k < theAlgorithm_->getNumberOfChannels(); k++) {
1396  bool isNearCrack =
1397  (abs(ringNumberCorrector(k)) == 1 || abs(ringNumberCorrector(k)) == 25 || abs(ringNumberCorrector(k)) == 26 ||
1398  abs(ringNumberCorrector(k)) == 45 || abs(ringNumberCorrector(k)) == 46 || abs(ringNumberCorrector(k)) == 65 ||
1399  abs(ringNumberCorrector(k)) == 66 || abs(ringNumberCorrector(k)) == 85 || abs(ringNumberCorrector(k)) == 86 ||
1400  abs(ringNumberCorrector(k)) == 124);
1401 
1402  if (!isNearCrack) {
1403  // h2_miscalRecalParz_[iLoop]->Fill( initCalibCoeff[k], 1./calibCoeff[k] );
1404  h1_mcParz_[iLoop]->Fill(initCalibCoeff[k] * calibCoeff[k] - 1.);
1405 
1406  if (k < 170) {
1407  //h2_miscalRecalEBParz_[iLoop]->Fill( initCalibCoeff[k], 1./calibCoeff[k] );
1408  h1_mcEBParz_[iLoop]->Fill(initCalibCoeff[k] * calibCoeff[k] - 1.);
1409  }
1410 
1411  if (k >= 170) {
1412  //h2_miscalRecalEEParz_[iLoop]->Fill( initCalibCoeff[k], 1./calibCoeff[k] );
1413  h1_mcEEParz_[iLoop]->Fill(initCalibCoeff[k] * calibCoeff[k] - 1.);
1414  }
1415  }
1416  }
1417 
1419  double parResidual[3];
1420  double errparResidual[3];
1421  double zResChi2;
1422  int zResIters;
1423 
1424  ZIterativeAlgorithmWithFit::gausfit(h1_mcParz_[iLoop], parResidual, errparResidual, 3., 3., &zResChi2, &zResIters);
1425  //h1_mcParz_[iLoop]->Fit("gaus");
1426 
1427  h2_residualSigma_->Fill(loopFlag_ + 1, parResidual[2]);
1428  loopArray[loopFlag_] = loopFlag_ + 1;
1429  sigmaArray[loopFlag_] = parResidual[2];
1430  sigmaErrorArray[loopFlag_] = errparResidual[2];
1431 
1432  std::cout << "Fit on residuals, sigma is " << parResidual[2] << " +/- " << errparResidual[2] << std::endl;
1433 
1435  outputFile_->cd();
1436 
1437  // h2_miscalRecalParz_[iLoop]->Write();
1438  h1_mcParz_[iLoop]->Write();
1439 
1440  //h2_miscalRecalEBParz_[iLoop]->Write();
1441  h1_mcEBParz_[iLoop]->Write();
1442 
1443  //h2_miscalRecalEEParz_[iLoop]->Write();
1444  h1_mcEEParz_[iLoop]->Write();
1448 
1450 
1451  loopFlag_++;
1452 
1453 #ifdef DEBUG
1454  std::cout << " loopFlag_ is " << loopFlag_ << std::endl;
1455 #endif
1456 
1457  if (iLoop == theMaxLoops - 1 || iLoop >= theMaxLoops)
1458  return kStop;
1459  else
1460  return kContinue;
1461 }
1462 
1464  h1_eventsBeforeEWKSelection_ = new TH1F("h1_eventsBeforeEWKSelection", "h1_eventsBeforeEWKSelection", 5, 0, 5);
1465  h1_eventsAfterEWKSelection_ = new TH1F("h1_eventsAfterEWKSelection", "h1_eventsAfterEWKSelection", 5, 0, 5);
1466 
1468  new TH1F("h1_eventsBeforeBorderSelection", "h1_eventsBeforeBorderSelection", 5, 0, 5);
1469  h1_eventsAfterBorderSelection_ = new TH1F("h1_eventsAfterBorderSelection", "h1_eventsAfterBorderSelection", 5, 0, 5);
1470 
1471  h1_seedOverSC_ = new TH1F("h1_seedOverSC", "h1_seedOverSC", 400, 0., 2.);
1472 
1474 
1476  new TH1F("h1_borderElectronClassification", "h1_borderElectronClassification", 55, -5, 50);
1477  h1_preshowerOverSC_ = new TH1F("h1_preshowerOverSC", "h1_preshowerOverSC", 400, 0., 1.);
1478 
1479  h2_fEtaBarrelGood_ = new TH2F("fEtaBarrelGood", "fEtaBarrelGood", 800, -4., 4., 800, 0.8, 1.2);
1480  h2_fEtaBarrelGood_->SetXTitle("Eta");
1481  h2_fEtaBarrelGood_->SetYTitle("1/fEtaBarrelGood");
1482 
1483  h2_fEtaBarrelBad_ = new TH2F("fEtaBarrelBad", "fEtaBarrelBad", 800, -4., 4., 800, 0.8, 1.2);
1484  h2_fEtaBarrelBad_->SetXTitle("Eta");
1485  h2_fEtaBarrelBad_->SetYTitle("1/fEtaBarrelBad");
1486 
1487  h2_fEtaEndcapGood_ = new TH2F("fEtaEndcapGood", "fEtaEndcapGood", 800, -4., 4., 800, 0.8, 1.2);
1488  h2_fEtaEndcapGood_->SetXTitle("Eta");
1489  h2_fEtaEndcapGood_->SetYTitle("1/fEtaEndcapGood");
1490 
1491  h2_fEtaEndcapBad_ = new TH2F("fEtaEndcapBad", "fEtaEndcapBad", 800, -4., 4., 800, 0.8, 1.2);
1492  h2_fEtaEndcapBad_->SetXTitle("Eta");
1493  h2_fEtaEndcapBad_->SetYTitle("1/fEtaEndcapBad");
1494 
1495  for (int i = 0; i < 2; i++) {
1496  char histoName[50];
1497 
1498  sprintf(histoName, "h_eleEffEta_%d", i);
1499  h_eleEffEta_[i] = new TH1F(histoName, histoName, 150, 0., 2.7);
1500  h_eleEffEta_[i]->SetXTitle("|#eta|");
1501 
1502  sprintf(histoName, "h_eleEffPhi_%d", i);
1503  h_eleEffPhi_[i] = new TH1F(histoName, histoName, 400, -4., 4.);
1504  h_eleEffPhi_[i]->SetXTitle("Phi");
1505 
1506  sprintf(histoName, "h_eleEffPt_%d", i);
1507  h_eleEffPt_[i] = new TH1F(histoName, histoName, 200, 0., 200.);
1508  h_eleEffPt_[i]->SetXTitle("p_{T}(GeV/c)");
1509  }
1510 
1512  new TH2F("h2_xtalMiscalibCoeffBarrel", "h2_xtalMiscalibCoeffBarrel", 171, -85, 85, 360, 0, 360);
1514  new TH2F("h2_xtalMiscalibCoeffEndcapMinus", "h2_xtalMiscalibCoeffEndcapMinus", 100, 0, 100, 100, 0, 100);
1516  new TH2F("h2_xtalMiscalibCoeffEndcapPlus", "h2_xtalMiscalibCoeffEndcapPlus", 100, 0, 100, 100, 0, 100);
1517 
1518  h2_xtalMiscalibCoeffBarrel_->SetXTitle("ieta");
1519  h2_xtalMiscalibCoeffBarrel_->SetYTitle("iphi");
1520 
1521  h2_xtalMiscalibCoeffEndcapMinus_->SetXTitle("ix");
1522  h2_xtalMiscalibCoeffEndcapMinus_->SetYTitle("iy");
1523 
1524  for (int i = 0; i < 25; i++) {
1525  char histoName[50];
1526  sprintf(histoName, "h_ESCEtrueVsEta_%d", i);
1527 
1528  h_ESCEtrueVsEta_[i] = new TH2F(histoName, histoName, 150, 0., 2.7, 300, 0., 1.5);
1529  h_ESCEtrueVsEta_[i]->SetXTitle("|#eta|");
1530  h_ESCEtrueVsEta_[i]->SetYTitle("E_{SC,raw}/E_{MC}");
1531 
1532  sprintf(histoName, "h_ESCEtrue_%d", i);
1533 
1534  h_ESCEtrue_[i] = new TH1F(histoName, histoName, 300, 0., 1.5);
1535 
1536  sprintf(histoName, "h2_chi2_%d", i);
1537  h2_chi2_[i] = new TH2F(histoName, histoName, 1000, -150, 150, 1000, -1, 5);
1538 
1539  sprintf(histoName, "h2_iterations_%d", i);
1540  h2_iterations_[i] = new TH2F(histoName, histoName, 1000, -150, 150, 1000, -1, 15);
1541 
1542  sprintf(histoName, "h_ESCcorrEtrueVsEta_%d", i);
1543 
1544  h_ESCcorrEtrueVsEta_[i] = new TH2F(histoName, histoName, 150, 0., 2.7, 300, 0., 1.5);
1545  h_ESCcorrEtrueVsEta_[i]->SetXTitle("|#eta|");
1546  h_ESCcorrEtrueVsEta_[i]->SetYTitle("E_{SC,#eta-corr}/E_{MC}");
1547 
1548  sprintf(histoName, "h_ESCcorrEtrue_%d", i);
1549 
1550  h_ESCcorrEtrue_[i] = new TH1F(histoName, histoName, 300, 0., 1.5);
1551 
1552  sprintf(histoName, "h2_xtalRecalibCoeffBarrel_%d", i);
1553  h2_xtalRecalibCoeffBarrel_[i] = new TH2F(histoName, histoName, 171, -85, 85, 360, 0, 360);
1554 
1555  h2_xtalRecalibCoeffBarrel_[i]->SetXTitle("ieta");
1556  h2_xtalRecalibCoeffBarrel_[i]->SetYTitle("iphi");
1557 
1558  sprintf(histoName, "h2_xtalRecalibCoeffEndcapMinus_%d", i);
1559  h2_xtalRecalibCoeffEndcapMinus_[i] = new TH2F(histoName, histoName, 100, 0, 100, 100, 0, 100);
1560  h2_xtalRecalibCoeffEndcapMinus_[i]->SetXTitle("ix");
1561  h2_xtalRecalibCoeffEndcapMinus_[i]->SetYTitle("iy");
1562 
1563  sprintf(histoName, "h2_xtalRecalibCoeffEndcapPlus_%d", i);
1564  h2_xtalRecalibCoeffEndcapPlus_[i] = new TH2F(histoName, histoName, 100, 0, 100, 100, 0, 100);
1565  h2_xtalRecalibCoeffEndcapPlus_[i]->SetXTitle("ix");
1566  h2_xtalRecalibCoeffEndcapPlus_[i]->SetYTitle("iy");
1567  }
1568 
1569  /*
1570  for (int i=0;i<15;i++)
1571  {
1572 
1573  char histoName[50];
1574 
1575  sprintf(histoName,"h_DiffZMassDistr_%d",i);
1576  h_DiffZMassDistr_[i] = new TH1F(histoName,histoName, 400, -20., 20.);
1577  h_DiffZMassDistr_[i]->SetXTitle("M_{Z, reco} - M_{Z, MC}");
1578  h_DiffZMassDistr_[i]->SetYTitle("events");
1579 
1580  sprintf(histoName,"h_ZMassDistr_%d",i);
1581  h_ZMassDistr_[i] = new TH1F(histoName,histoName, 200, 0., 150.);
1582  h_ZMassDistr_[i]->SetXTitle("RecoZmass (GeV)");
1583  h_ZMassDistr_[i]->SetYTitle("events");
1584 
1585  }
1586  */
1587 
1588  h1_zMassResol_ = new TH1F("zMassResol", "zMassResol", 200, -50., 50.);
1589  h1_zMassResol_->SetXTitle("M_{Z, reco} - M_{Z, MC}");
1590  h1_zMassResol_->SetYTitle("events");
1591 
1592  h1_eleEtaResol_ = new TH1F("eleEtaResol", "eleEtaResol", 100, -0.01, 0.01);
1593  h1_eleEtaResol_->SetXTitle("#eta_{reco} - #eta_{MC}");
1594  h1_eleEtaResol_->SetYTitle("events");
1595 
1596  h1_electronCosTheta_TK_ = new TH1F("electronCosTheta_TK", "electronCosTheta_TK", 100, -1, 1);
1597  h1_electronCosTheta_TK_->SetXTitle("cos #theta_{12}");
1598  h1_electronCosTheta_TK_->SetYTitle("events");
1599 
1600  h1_electronCosTheta_SC_ = new TH1F("electronCosTheta_SC", "electronCosTheta_SC", 100, -1, 1);
1601  h1_electronCosTheta_SC_->SetXTitle("cos #theta_{12}");
1602  h1_electronCosTheta_SC_->SetYTitle("events");
1603 
1604  h1_electronCosTheta_SC_TK_ = new TH1F("electronCosTheta_SC_TK", "electronCosTheta_SC_TK", 200, -0.1, 0.1);
1605  h1_electronCosTheta_SC_TK_->SetXTitle("cos #theta_{12}^{SC}/ cos #theta_{12}^{TK} - 1");
1606  h1_electronCosTheta_SC_TK_->SetYTitle("events");
1607 
1608  h1_elePhiResol_ = new TH1F("elePhiResol", "elePhiResol", 100, -0.01, 0.01);
1609  h1_elePhiResol_->SetXTitle("#phi_{reco} - #phi_{MC}");
1610  h1_elePhiResol_->SetYTitle("events");
1611 
1612  h1_zEtaResol_ = new TH1F("zEtaResol", "zEtaResol", 200, -1., 1.);
1613  h1_zEtaResol_->SetXTitle("#eta_{Z, reco} - #eta_{Z, MC}");
1614  h1_zEtaResol_->SetYTitle("events");
1615 
1616  h1_zPhiResol_ = new TH1F("zPhiResol", "zPhiResol", 200, -1., 1.);
1617  h1_zPhiResol_->SetXTitle("#phi_{Z, reco} - #phi_{Z, MC}");
1618  h1_zPhiResol_->SetYTitle("events");
1619 
1620  h1_nEleReco_ = new TH1F("nEleReco", "Number of reco electrons", 10, -0.5, 10.5);
1621  h1_nEleReco_->SetXTitle("nEleReco");
1622  h1_nEleReco_->SetYTitle("events");
1623 
1624  // h1_occupancyVsEta_ = new TH1F("occupancyVsEta","occupancyVsEta",EcalRingCalibrationTools::N_RING_TOTAL,0,(float)EcalRingCalibrationTools::N_RING_TOTAL);
1625 
1626  h1_occupancyVsEta_ = new TH1F("occupancyVsEta", "occupancyVsEta", 249, -124, 124);
1627  h1_occupancyVsEta_->SetYTitle("Weighted electron statistics");
1628  h1_occupancyVsEta_->SetXTitle("Eta channel");
1629 
1630  h1_weightSumMeanBarrel_ = new TH1F("weightSumMeanBarrel", "weightSumMeanBarrel", 10000, 0, 10000);
1631  h1_weightSumMeanEndcap_ = new TH1F("weightSumMeanEndcap", "weightSumMeanEndcap", 10000, 0, 10000);
1632 
1633  h1_occupancy_ = new TH1F("occupancy", "occupancy", 1000, 0, 10000);
1634  h1_occupancy_->SetXTitle("Weighted electron statistics");
1635 
1636  h1_occupancyBarrel_ = new TH1F("occupancyBarrel", "occupancyBarrel", 1000, 0, 10000);
1637  h1_occupancyBarrel_->SetXTitle("Weighted electron statistics");
1638 
1639  h1_occupancyEndcap_ = new TH1F("occupancyEndcap", "occupancyEndcap", 1000, 0, 10000);
1640  h1_occupancyEndcap_->SetXTitle("Weighted electron statistics");
1641 
1642  h1_eleClasses_ = new TH1F("eleClasses", "eleClasses", 301, -1, 300);
1643  h1_eleClasses_->SetXTitle("classCode");
1644  h1_eleClasses_->SetYTitle("#");
1645 
1647 
1649 
1651 
1653 
1654  h1_ZCandMult_ = new TH1F("ZCandMult", "Multiplicity of Z candidates in one event", 10, -0.5, 10.5);
1655  h1_ZCandMult_->SetXTitle("ZCandMult");
1656 
1657  h1_reco_ZMass_ = new TH1F("reco_ZMass", "Inv. mass of 2 reco Electrons", 200, 0., 150.);
1658  h1_reco_ZMass_->SetXTitle("reco_ZMass (GeV)");
1659  h1_reco_ZMass_->SetYTitle("events");
1660 
1661  h1_reco_ZMassCorr_ = new TH1F("reco_ZMassCorr", "Inv. mass of 2 corrected reco Electrons", 200, 0., 150.);
1662  h1_reco_ZMassCorr_->SetXTitle("reco_ZMass (GeV)");
1663  h1_reco_ZMassCorr_->SetYTitle("events");
1664 
1665  h1_reco_ZMassCorrBB_ = new TH1F("reco_ZMassCorrBB", "Inv. mass of 2 corrected reco Electrons", 200, 0., 150.);
1666  h1_reco_ZMassCorrBB_->SetXTitle("reco_ZMass (GeV)");
1667  h1_reco_ZMassCorrBB_->SetYTitle("events");
1668 
1669  h1_reco_ZMassCorrEE_ = new TH1F("reco_ZMassCorrEE", "Inv. mass of 2 corrected reco Electrons", 200, 0., 150.);
1670  h1_reco_ZMassCorrEE_->SetXTitle("reco_ZMass (GeV)");
1671  h1_reco_ZMassCorrEE_->SetYTitle("events");
1672 
1673  // h2_coeffVsEta_= new TH2F("h2_calibCoeffVsEta","h2_calibCoeffVsEta",EcalRingCalibrationTools::N_RING_TOTAL,0, (double)EcalRingCalibrationTools::N_RING_TOTAL, 200, 0., 2.);
1674 
1675  h2_coeffVsEta_ = new TH2F("h2_calibCoeffVsEta", "h2_calibCoeffVsEta", 249, -124, 125, 200, 0., 2.);
1676  h2_coeffVsEta_->SetXTitle("Eta channel");
1677  h2_coeffVsEta_->SetYTitle("recalibCoeff");
1678 
1680  new TH2F("h2_calibCoeffVsEtaGrouped", "h2_calibCoeffVsEtaGrouped", 200, 0., 3., 200, 0.6, 1.4);
1681  h2_coeffVsEtaGrouped_->SetXTitle("|#eta|");
1682  h2_coeffVsEtaGrouped_->SetYTitle("recalibCoeff");
1683 
1684  h2_zMassVsLoop_ = new TH2F("h2_zMassVsLoop", "h2_zMassVsLoop", 1000, 0, 40, 90, 80., 95.);
1685 
1686  h2_zMassDiffVsLoop_ = new TH2F("h2_zMassDiffVsLoop", "h2_zMassDiffVsLoop", 1000, 0, 40, 100, -1., 1.);
1687  h2_zMassDiffVsLoop_->SetXTitle("Iteration");
1688  h2_zMassDiffVsLoop_->SetYTitle("M_{Z, reco peak} - M_{Z, true}");
1689 
1690  h2_zWidthVsLoop_ = new TH2F("h2_zWidthVsLoop", "h2_zWidthVsLoop", 1000, 0, 40, 100, 0., 10.);
1691 
1692  h2_coeffVsLoop_ = new TH2F("h2_coeffVsLoop", "h2_coeffVsLoop", 1000, 0, 40, 100, 0., 2.);
1693 
1694  h2_residualSigma_ = new TH2F("h2_residualSigma", "h2_residualSigma", 1000, 0, 40, 100, 0., .5);
1695 
1696  h2_miscalRecal_ = new TH2F("h2_miscalRecal", "h2_miscalRecal", 500, 0., 2., 500, 0., 2.);
1697  h2_miscalRecal_->SetXTitle("initCalibCoeff");
1698  h2_miscalRecal_->SetYTitle("1/RecalibCoeff");
1699 
1700  h2_miscalRecalEB_ = new TH2F("h2_miscalRecalEB", "h2_miscalRecalEB", 500, 0., 2., 500, 0., 2.);
1701  h2_miscalRecalEB_->SetXTitle("initCalibCoeff");
1702  h2_miscalRecalEB_->SetYTitle("1/RecalibCoeff");
1703 
1704  h2_miscalRecalEE_ = new TH2F("h2_miscalRecalEE", "h2_miscalRecalEE", 500, 0., 2., 500, 0., 2.);
1705  h2_miscalRecalEE_->SetXTitle("initCalibCoeff");
1706  h2_miscalRecalEE_->SetYTitle("1/RecalibCoeff");
1707 
1708  h1_mc_ = new TH1F("h1_residualMiscalib", "h1_residualMiscalib", 200, -0.2, 0.2);
1709  h1_mcEB_ = new TH1F("h1_residualMiscalibEB", "h1_residualMiscalibEB", 200, -0.2, 0.2);
1710  h1_mcEE_ = new TH1F("h1_residualMiscalibEE", "h1_residualMiscalibEE", 200, -0.2, 0.2);
1711 
1712  for (int i = 0; i < 25; i++) {
1713  char histoName[50];
1714  /*
1715  sprintf(histoName,"h2_miscalRecalParz_%d",i);
1716  h2_miscalRecalParz_[i] = new TH2F(histoName,histoName,500, 0., 2., 500, 0., 2.);
1717  h2_miscalRecalParz_[i]->SetXTitle("initCalibCoeff");
1718  h2_miscalRecalParz_[i]->SetYTitle("1/recalibCoeff");
1719 
1720  sprintf(histoName,"h2_miscalRecalEBParz_%d",i);
1721  h2_miscalRecalEBParz_[i] = new TH2F(histoName,histoName,500, 0., 2., 500, 0., 2.);
1722  h2_miscalRecalEBParz_[i]->SetXTitle("initCalibCoeff");
1723  h2_miscalRecalEBParz_[i]->SetYTitle("1/recalibCoeff");
1724 
1725  sprintf(histoName,"h2_miscalRecalEEParz_%d",i);
1726  h2_miscalRecalEEParz_[i] = new TH2F(histoName,histoName,500, 0., 2., 500, 0., 2.);
1727  h2_miscalRecalEEParz_[i]->SetXTitle("initCalibCoeff");
1728  h2_miscalRecalEEParz_[i]->SetYTitle("1/recalibCoeff");
1729  */
1730 
1731  sprintf(histoName, "h1_residualMiscalibParz_%d", i);
1732  h1_mcParz_[i] = new TH1F(histoName, histoName, 200, -0.2, 0.2);
1733  sprintf(histoName, "h1_residualMiscalibEBParz_%d", i);
1734  h1_mcEBParz_[i] = new TH1F(histoName, histoName, 200, -0.2, 0.2);
1735  sprintf(histoName, "h1_residualMiscalibEEParz_%d", i);
1736  h1_mcEEParz_[i] = new TH1F(histoName, histoName, 200, -0.2, 0.2);
1737  }
1738 }
1739 
1740 double ZeeCalibration::fEtaBarrelBad(double scEta) const {
1741  float p0 = 1.00153e+00;
1742  float p1 = 3.29331e-02;
1743  float p2 = 1.21187e-03;
1744 
1745  double x = (double)fabs(scEta);
1746 
1747  return 1. / (p0 + p1 * x * x + p2 * x * x * x * x);
1748 }
1749 
1750 double ZeeCalibration::fEtaEndcapGood(double scEta) const {
1751  // f(eta) for the first 3 classes (100, 110 and 120)
1752  // Ivica's new corrections 01/06
1753  float p0 = 1.06819e+00;
1754  float p1 = -1.53189e-02;
1755  float p2 = 4.01707e-04;
1756 
1757  double x = (double)fabs(scEta);
1758 
1759  return 1. / (p0 + p1 * x * x + p2 * x * x * x * x);
1760 }
1761 
1762 double ZeeCalibration::fEtaEndcapBad(double scEta) const {
1763  float p0 = 1.17382e+00;
1764  float p1 = -6.52319e-02;
1765  float p2 = 6.26108e-03;
1766 
1767  double x = (double)fabs(scEta);
1768 
1769  return 1. / (p0 + p1 * x * x + p2 * x * x * x * x);
1770 }
1771 
1772 double ZeeCalibration::fEtaBarrelGood(double scEta) const {
1773  float p0 = 9.99782e-01;
1774  float p1 = 1.26983e-02;
1775  float p2 = 2.16344e-03;
1776 
1777  double x = (double)fabs(scEta);
1778 
1779  return 1. / (p0 + p1 * x * x + p2 * x * x * x * x);
1780 }
1781 
1783 
1784 void ZeeCalibration::fillMCmap(const std::vector<const reco::GsfElectron*>* electronCollection,
1785  const std::vector<HepMC::GenParticle*>& mcEle,
1786  std::map<HepMC::GenParticle*, const reco::GsfElectron*>& myMCmap) {
1787  for (unsigned int i = 0; i < mcEle.size(); i++) {
1788  float minDR = 0.1;
1789  const reco::GsfElectron* myMatchEle = nullptr;
1790  for (unsigned int j = 0; j < electronCollection->size(); j++) {
1791  float dr = EvalDR(mcEle[i]->momentum().pseudoRapidity(),
1792  (*(*electronCollection)[j]).eta(),
1793  mcEle[i]->momentum().phi(),
1794  (*(*electronCollection)[j]).phi());
1795  if (dr < minDR) {
1796  myMatchEle = (*electronCollection)[j];
1797  minDR = dr;
1798  }
1799  }
1800  myMCmap.insert(std::pair<HepMC::GenParticle*, const reco::GsfElectron*>(mcEle[i], myMatchEle));
1801  }
1802 }
1803 
1804 float ZeeCalibration::EvalDR(float Eta, float Eta_ref, float Phi, float Phi_ref) {
1805  if (Phi < 0)
1806  Phi = 2 * TMath::Pi() + Phi;
1807  if (Phi_ref < 0)
1808  Phi_ref = 2 * TMath::Pi() + Phi_ref;
1809  float DPhi = Phi - Phi_ref;
1810  if (fabs(DPhi) > TMath::Pi())
1811  DPhi = 2 * TMath::Pi() - fabs(DPhi);
1812 
1813  float DEta = Eta - Eta_ref;
1814 
1815  float DR = sqrt(DEta * DEta + DPhi * DPhi);
1816  return DR;
1817 }
1818 
1819 float ZeeCalibration::EvalDPhi(float Phi, float Phi_ref) {
1820  if (Phi < 0)
1821  Phi = 2 * TMath::Pi() + Phi;
1822  if (Phi_ref < 0)
1823  Phi_ref = 2 * TMath::Pi() + Phi_ref;
1824  return (Phi - Phi_ref);
1825 }
1826 
1827 void ZeeCalibration::fillEleInfo(std::vector<HepMC::GenParticle*>& mcEle,
1828  std::map<HepMC::GenParticle*, const reco::GsfElectron*>& associationMap) {
1829  for (unsigned int i = 0; i < mcEle.size(); i++) {
1830  h_eleEffEta_[0]->Fill(fabs(mcEle[i]->momentum().pseudoRapidity()));
1831  h_eleEffPhi_[0]->Fill(mcEle[i]->momentum().phi());
1832  h_eleEffPt_[0]->Fill(mcEle[i]->momentum().perp());
1833 
1834  std::map<HepMC::GenParticle*, const reco::GsfElectron*>::const_iterator mIter = associationMap.find(mcEle[i]);
1835  if (mIter == associationMap.end())
1836  continue;
1837 
1838  if ((*mIter).second) {
1839  const reco::GsfElectron* myEle = (*mIter).second;
1840 
1841  h_eleEffEta_[1]->Fill(fabs(mcEle[i]->momentum().pseudoRapidity()));
1842  h_eleEffPhi_[1]->Fill(mcEle[i]->momentum().phi());
1843  h_eleEffPt_[1]->Fill(mcEle[i]->momentum().perp());
1844  h1_eleEtaResol_->Fill(myEle->eta() - mcEle[i]->momentum().eta());
1845  h1_elePhiResol_->Fill(myEle->phi() - mcEle[i]->momentum().phi());
1846 
1847  const reco::SuperCluster* mySC = &(*(myEle->superCluster()));
1848  if (/*fabs(mySC->position().eta()) < 2.4*/ true) {
1849  // if(myEle->classification()>=100)std::cout<<"mySC->preshowerEnergy()"<<mySC->preshowerEnergy()<<std::endl;
1850 
1851  h_ESCEtrue_[loopFlag_]->Fill(mySC->energy() / mcEle[i]->momentum().e());
1852  h_ESCEtrueVsEta_[loopFlag_]->Fill(fabs(mySC->position().eta()), mySC->energy() / mcEle[i]->momentum().e());
1853 
1854  double corrSCenergy = (mySC->energy()) / getEtaCorrection(myEle);
1855  h_ESCcorrEtrue_[loopFlag_]->Fill(corrSCenergy / mcEle[i]->momentum().e());
1856  h_ESCcorrEtrueVsEta_[loopFlag_]->Fill(fabs(mySC->position().eta()), corrSCenergy / mcEle[i]->momentum().e());
1857 
1858  // std::vector<DetId> mySCRecHits = mySC->seed()->getHitsByDetId();
1859 
1860  h1_seedOverSC_->Fill(mySC->seed()->energy() / mySC->energy());
1861  h1_preshowerOverSC_->Fill(mySC->preshowerEnergy() / mySC->energy());
1862  }
1863  }
1864  }
1865 }
1866 
1868  int index = -999;
1869 
1870  if (calibMode_ == "RING") {
1871  if (k >= 0 && k <= 84)
1872  index = k - 85;
1873 
1874  if (k >= 85 && k <= 169)
1875  index = k - 84;
1876 
1877  if (k >= 170 && k <= 208)
1878  index = -k + 84;
1879 
1880  if (k >= 209 && k <= 247)
1881  index = k - 123;
1882 
1883  }
1884 
1885  else if (calibMode_ == "MODULE") {
1886  if (k >= 0 && k <= 71)
1887  index = k - 72;
1888 
1889  if (k >= 72 && k <= 143)
1890  index = k - 71;
1891  }
1892  return index;
1893 }
1894 
1896  double correction(1.);
1897 
1898  int c = ele->classification();
1899  if (c == 0 || c == 10 || c == 20)
1900  correction = fEtaBarrelGood(ele->superCluster()->eta());
1901 
1902  if (c == 100 || c == 110 || c == 120)
1903  correction = fEtaEndcapGood(ele->superCluster()->eta());
1904 
1905  if (c == 30 || c == 31 || c == 32 || c == 33 || c == 34)
1906  correction = fEtaBarrelBad(ele->superCluster()->eta());
1907 
1908  if (c == 130 || c == 131 || c == 132 || c == 133 || c == 134)
1909  correction = fEtaEndcapBad(ele->superCluster()->eta());
1910 
1911  return correction;
1912 }
1913 
1914 std::pair<DetId, double> ZeeCalibration::getHottestDetId(const std::vector<std::pair<DetId, float> >& mySCRecHits,
1915  const EBRecHitCollection* ebhits,
1916  const EERecHitCollection* eehits) {
1917  double maxEnergy = -9999.;
1918  const EcalRecHit* hottestRecHit = nullptr;
1919 
1920  std::pair<DetId, double> myPair(DetId(0), -9999.);
1921 
1922  for (std::vector<std::pair<DetId, float> >::const_iterator idIt = mySCRecHits.begin(); idIt != mySCRecHits.end();
1923  idIt++) {
1924  if (idIt->first.subdetId() == EcalBarrel) {
1925  hottestRecHit = &(*(ebhits->find((*idIt).first)));
1926 
1927  if (hottestRecHit == &(*(ebhits->end()))) {
1928  std::cout << "@@@@@@@@@@@@@@@@@@@@@@@@@@@ NO RECHIT FOUND SHOULD NEVER HAPPEN" << std::endl;
1929  continue;
1930  }
1931  } else if (idIt->first.subdetId() == EcalEndcap) {
1932  hottestRecHit = &(*(eehits->find((*idIt).first)));
1933  if (hottestRecHit == &(*(eehits->end()))) {
1934  std::cout << "@@@@@@@@@@@@@@@@@@@@@@@@@@@ NO RECHIT FOUND SHOULD NEVER HAPPEN" << std::endl;
1935  continue;
1936  }
1937  }
1938 
1939  //std::cout<<"[getHottestDetId] hottestRecHit->energy() "<<hottestRecHit->energy()<<std::endl;
1940 
1941  if (hottestRecHit && hottestRecHit->energy() > maxEnergy) {
1942  maxEnergy = hottestRecHit->energy();
1943 
1944  myPair.first = hottestRecHit->id();
1945  myPair.second = maxEnergy;
1946  }
1947 
1948  } //end loop to find hottest RecHit
1949 
1950  //std::cout<<"[ZeeCalibration::getHottestDetId] going to return..."<<std::endl;
1951 
1952  return myPair;
1953 }
1954 
1956  bool myBool(false);
1957 
1958  short ieta = myEBDetId.ieta();
1959  short iphi = myEBDetId.iphi();
1960 
1961  // std::cout<<"[xtalIsOnModuleBorder] ieta: "<<ieta<<" iphi "<<iphi<<std::endl;
1962 
1963  myBool = (abs(ieta) == 1 || abs(ieta) == 25 || abs(ieta) == 26 || abs(ieta) == 45 || abs(ieta) == 46 ||
1964  abs(ieta) == 65 || abs(ieta) == 66 || abs(ieta) == 85);
1965 
1966  for (int i = 0; i < 19; i++) {
1967  if (iphi == (20 * i + 1) || iphi == 20 * i)
1968  myBool = true;
1969  }
1970 
1971  return myBool;
1972 }
1973 
1974 float ZeeCalibration::computeCoefficientDistanceAtIteration(float v1[250], float v2[250], int size) {
1975  float dist(0.);
1976 
1977  for (int i = 0; i < size; i++) {
1978  // std::cout<< "[ZeeCalibration::computeCoefficientDistanceAtIteration] Adding term "<<pow( v1[i]-v2[i], 2 )<<" from v1 "<<v1[i]<<" and v2 "<<v2[i]<<std::endl;
1979 
1980  bool isNearCrack = false;
1981 
1982  if (calibMode_ == "RING") { //exclude non-calibrated rings from computation
1983 
1984  isNearCrack = (abs(ringNumberCorrector(i)) == 1 || abs(ringNumberCorrector(i)) == 25 ||
1985  abs(ringNumberCorrector(i)) == 26 || abs(ringNumberCorrector(i)) == 45 ||
1986  abs(ringNumberCorrector(i)) == 46 || abs(ringNumberCorrector(i)) == 65 ||
1987  abs(ringNumberCorrector(i)) == 66 || abs(ringNumberCorrector(i)) == 85 ||
1988  abs(ringNumberCorrector(i)) == 86 || abs(ringNumberCorrector(i)) == 124);
1989  }
1990 
1991  if (!isNearCrack)
1992  dist += pow(v1[i] - v2[i], 2);
1993  }
1994 
1995  dist = sqrt(dist) / size;
1996 
1997  return dist;
1998 }
1999 
2001  BBZN = 0;
2002  EBZN = 0;
2003  EEZN = 0;
2004  BBZN_gg = 0;
2005  EBZN_gg = 0;
2006  EEZN_gg = 0;
2007 
2008  BBZN_tt = 0;
2009  EBZN_tt = 0;
2010  EEZN_tt = 0;
2011 
2012  BBZN_t0 = 0;
2013  EBZN_t0 = 0;
2014  EEZN_t0 = 0;
2015 
2018 
2027 
2030 
2031  return;
2032 }
2033 
2036  h1_eventsAfterEWKSelection_->Reset();
2037 
2040 
2041  for (int i = 0; i < 2; i++) {
2042  h_eleEffEta_[i]->Reset();
2043  h_eleEffPhi_[i]->Reset();
2044  h_eleEffPt_[i]->Reset();
2045  }
2046 
2047  h1_seedOverSC_->Reset();
2048  h1_preshowerOverSC_->Reset();
2049 
2050  h1_eleEtaResol_->Reset();
2051  h1_elePhiResol_->Reset();
2052 
2053  h1_zMassResol_->Reset();
2054 
2055  h1_electronCosTheta_TK_->Reset();
2056  h1_electronCosTheta_SC_->Reset();
2057  h1_electronCosTheta_SC_TK_->Reset();
2058 
2059  h2_fEtaBarrelGood_->Reset();
2060  h2_fEtaBarrelBad_->Reset();
2061  h2_fEtaEndcapGood_->Reset();
2062  h2_fEtaEndcapBad_->Reset();
2063  h1_eleClasses_->Reset();
2064 
2065  h1_ZCandMult_->Reset();
2066  h1_reco_ZMass_->Reset();
2067  h1_reco_ZMassCorr_->Reset();
2068  h1_reco_ZMassCorrBB_->Reset();
2069  h1_reco_ZMassCorrEE_->Reset();
2070  h1_occupancyVsEta_->Reset();
2071  h1_occupancy_->Reset();
2072  h1_occupancyBarrel_->Reset();
2073  h1_occupancyEndcap_->Reset();
2074 
2075  return;
2076 }
2077 
2079  std::cout << "[ CHECK ON BARREL ELECTRON NUMBER ]"
2080  << " first " << BARREL_ELECTRONS_BEFORE_BORDER_CUT << " second " << TOTAL_ELECTRONS_IN_BARREL << std::endl;
2081 
2082  std::cout << "[ EFFICIENCY OF THE BORDER SELECTION ]"
2084 
2085  std::cout << "[ EFFICIENCY OF THE GOLDEN SELECTION ] BARREL: "
2087  << " ENDCAP: " << (float)GOLDEN_ELECTRONS_IN_ENDCAP / (float)TOTAL_ELECTRONS_IN_ENDCAP << std::endl;
2088 
2089  std::cout << "[ EFFICIENCY OF THE SILVER SELECTION ] BARREL: "
2091  << " ENDCAP: " << (float)SILVER_ELECTRONS_IN_ENDCAP / (float)TOTAL_ELECTRONS_IN_ENDCAP << std::endl;
2092 
2093  std::cout << "[ EFFICIENCY OF THE SHOWER SELECTION ] BARREL: "
2095  << " ENDCAP: " << (float)SHOWER_ELECTRONS_IN_ENDCAP / (float)TOTAL_ELECTRONS_IN_ENDCAP << std::endl;
2096 
2097  std::cout << "[ EFFICIENCY OF THE CRACK SELECTION ] BARREL: "
2099  << " ENDCAP: " << (float)CRACK_ELECTRONS_IN_ENDCAP / (float)TOTAL_ELECTRONS_IN_ENDCAP << std::endl;
2100 
2101  std::ofstream fout("ZeeStatistics.txt");
2102 
2103  if (!fout) {
2104  std::cout << "Cannot open output file.\n";
2105  }
2106 
2107  fout << "ZeeStatistics" << std::endl;
2108 
2109  fout << "##########################RECO#########################" << std::endl;
2110  fout << "##################Zee with Barrel-Barrel electrons: " << BBZN << std::endl;
2111 
2112  fout << "Golden-Golden fraction: " << (float)BBZN_gg / BBZN << " 3-3 fraction is " << (float)BBZN_tt / BBZN
2113  << " 3-whatever fraction is " << (float)BBZN_t0 / BBZN << std::endl;
2114  fout << "##################Zee with Barrel-Endcap electrons: " << EBZN << std::endl;
2115  fout << "Golden-Golden fraction: " << (float)EBZN_gg / EBZN << " 3-3 fraction is " << (float)EBZN_tt / EBZN
2116  << " 3-whatever fraction is " << (float)EBZN_t0 / EBZN << std::endl;
2117  fout << "##################Zee with Endcap-Endcap electrons: " << EEZN << std::endl;
2118  fout << "Golden-Golden fraction: " << (float)EEZN_gg / EEZN << " 3-3 fraction is " << (float)EEZN_tt / EEZN
2119  << " 3-whatever fraction is " << (float)EEZN_t0 / EEZN << std::endl;
2120 
2121  fout << "\n" << std::endl;
2122 
2123  fout << "##########################GEN#########################" << std::endl;
2124  fout << "##################Zee with Barrel-Barrel electrons: " << (float)MCZBB / NEVT << std::endl;
2125  fout << "##################Zee with Barrel-Endcap electrons: " << (float)MCZEB / NEVT << std::endl;
2126  fout << "##################Zee with Endcap-Endcap electrons: " << (float)MCZEE / NEVT << std::endl;
2127 
2128  fout.close();
2129 }
unsigned int etBins_
size
Write out results.
const math::XYZPoint & position() const
cluster centroid position
Definition: CaloCluster.h:154
bool accept() const
Has at least one path accepted the event?
const double Pi
auto setWhatProduced(T *iThis, const es::Label &iLabel={})
Definition: ESProducer.h:163
TH2F * h2_fEtaBarrelBad_
TH1F * h1_mcEEParz_[25]
ESGetTokenH3DDVariant esConsumes(std::string const &Record, edm::ConsumesCollector &)
Definition: DeDxTools.cc:283
const edm::EDGetTokenT< edm::HepMCProduct > hepMCToken_
std::string calibMode_
std::string outputFileName_
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
static const int MIN_IPHI
Definition: EBDetId.h:135
void setBinRange(int, double, double, int, double, double)
static void setCaloGeometry(const CaloGeometry *geometry)
const std::vector< float > & getOptimizedCoefficientsError() const
T const & getData(const ESGetToken< T, R > &iToken) const noexcept(false)
Definition: EventSetup.h:119
TH1F * h_eleEffPhi_[2]
const edm::ESGetToken< CaloGeometry, CaloGeometryRecord > geometryToken_
int BARREL_ELECTRONS_AFTER_BORDER_CUT
void fillEleInfo(std::vector< HepMC::GenParticle *> &a, std::map< HepMC::GenParticle *, const reco::GsfElectron *> &b)
TH2F * h2_miscalRecal_
TH1F * h1_eventsBeforeEWKSelection_
float computeCoefficientDistanceAtIteration(float v1[250], float v2[250], int size)
float calibCoeff[250]
void fillEleInfo(const reco::GsfElectronCollection *)
Definition: ZeePlots.cc:220
void fillEleMCInfo(const HepMC::GenEvent *)
Definition: ZeePlots.cc:178
const edm::EDGetTokenT< EERecHitCollection > eeRecHitToken_
unsigned int theMaxLoops
int GOLDEN_ELECTRONS_IN_ENDCAP
void fillZMCInfo(const HepMC::GenEvent *)
Definition: ZeePlots.cc:132
int iphi() const
get the crystal iphi
Definition: EBDetId.h:51
TH1F * h_ESCEtrue_[25]
TH1F * h1_eventsAfterBorderSelection_
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
T const * product() const
Definition: Handle.h:70
int ix() const
Definition: EEDetId.h:77
bool parseXMLMiscalibFile(std::string configFile)
void fillZInfo(std::pair< calib::CalibElectron *, calib::CalibElectron *> myZeeCandidate)
Definition: ZeePlots.cc:100
TH1F * h1_occupancyEndcap_
#define MZ
bool addEvent(calib::CalibElectron *, calib::CalibElectron *, float)
void writeMCEleHistograms()
Definition: ZeePlots.cc:255
Status duringLoop(const edm::Event &, const edm::EventSetup &) override
Called at each event.
TH1F * h1_eventsBeforeBorderSelection_
TH1F * h1_electronCosTheta_SC_
void writeMCZHistograms()
Definition: ZeePlots.cc:122
TH1F * h1_reco_ZMassCorr_
static void gausfit(TH1F *histoou, double *par, double *errpar, float nsigmalow, float nsigmaup, double *mychi2, int *iterations)
int GOLDEN_ELECTRONS_IN_BARREL
TH2F * h2_coeffVsEtaGrouped_
TH1F * h1_occupancyVsEta_
std::string endcapfile_
std::vector< GsfElectron > GsfElectronCollection
collection of GsfElectron objects
static std::vector< DetId > getDetIdsInRing(short aRingIndex)
Retrieve the DetIds in a phi-ring.
static const int IX_MIN
Definition: EEDetId.h:290
TH1F * h1_occupancyBarrel_
static const int IY_MIN
Definition: EEDetId.h:294
~ZeeCalibration() override
Destructor.
TH2F * h2_miscalRecalEE_
const edm::EDGetTokenT< edm::TriggerResults > trigResultsToken_
Classification classification() const
Definition: GsfElectron.h:805
static bool validDetId(int i, int j)
check if a valid index combination
Definition: EBDetId.h:118
float EvalDR(float Eta, float Eta_ref, float Phi, float Phi_ref)
int ieta() const
get the crystal ieta
Definition: EBDetId.h:49
T getUntrackedParameter(std::string const &, T const &) const
void bookEleHistograms()
Definition: ZeePlots.cc:195
U second(std::pair< T, U > const &p)
TH1F * h1_electronCosTheta_TK_
TH2F * h_ESCcorrEtrueVsEta_[25]
unsigned int size() const
Get number of paths stored.
int iEvent
Definition: GenABIO.cc:224
TH2F * h2_fEtaEndcapBad_
const EcalIntercalibConstants & get()
int ringNumberCorrector(int k)
TupleMultiplicity< TrackerTraits > const *__restrict__ TrackingRecHitSoAConstView< TrackerTraits > TrackerTraits::tindex_type *__restrict__ double *__restrict__ phits
int SILVER_ELECTRONS_IN_BARREL
std::vector< SuperCluster > SuperClusterCollection
collection of SuperCluser objectr
void beginOfJob() override
Called at beginning of job.
double sigmaArray[50]
int TOTAL_ELECTRONS_IN_ENDCAP
int BARREL_ELECTRONS_BEFORE_BORDER_CUT
T sqrt(T t)
Definition: SSEVec.h:19
void fillHLTInfo(edm::Handle< edm::TriggerResults >)
Definition: ZeePlots.cc:299
int SHOWER_ELECTRONS_IN_ENDCAP
static EcalIndexingTools * getInstance()
static float calculateZMassWithCorrectedElectrons_withTK(const std::pair< calib::CalibElectron *, calib::CalibElectron *> &aZCandidate, float ele1EnergyCorrection, float ele2EnergyCorrection)
void writeZHistograms()
Definition: ZeePlots.cc:110
static float cosThetaElectrons_TK(const std::pair< calib::CalibElectron *, calib::CalibElectron *> &aZCandidate, float ele1EnergyCorrection, float ele2EnergyCorrection)
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
TH1F * h1_mcEBParz_[25]
static float calculateZMass_withTK(const std::pair< calib::CalibElectron *, calib::CalibElectron *> &aZCandidate)
ZIterativeAlgorithmWithFit * theAlgorithm_
const edm::EDGetTokenT< reco::SuperClusterCollection > islandSCToken_
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:48
TH1F * h1_elePhiResol_
unsigned int electronSelection_
bool xtalIsOnModuleBorder(EBDetId myEBDetId)
TH1F * h1_eleEtaResol_
TH2F * h2_xtalRecalibCoeffBarrel_[25]
float EvalDPhi(float Phi, float Phi_ref)
TH1F * h1_eventsAfterEWKSelection_
T perp() const
Magnitude of transverse component.
static std::vector< DetId > getDetIdsInModule(short int)
void bookZMCHistograms()
Definition: ZeePlots.cc:50
TH1F * h1_preshowerOverSC_
const HepMC::GenEvent * GetEvent() const
Definition: HepMCProduct.h:37
const_iterator find(uint32_t rawId) const
TH1F * h1_electronCosTheta_SC_TK_
double coefficientDistanceAtIteration[50]
double fEtaBarrelBad(double scEta) const
void fillMCmap(const std::vector< const reco::GsfElectron *> *electronCollection, const std::vector< HepMC::GenParticle *> &mcEle, std::map< HepMC::GenParticle *, const reco::GsfElectron *> &myMCmap)
TH1F * h_eleEffPt_[2]
const_iterator end() const
static const int IX_MAX
Definition: EEDetId.h:298
TH2F * h2_zMassDiffVsLoop_
double energy() const
cluster energy
Definition: CaloCluster.h:149
void startingNewLoop(unsigned int iLoop) override
Called at beginning of loop.
Definition: DetId.h:17
TH1F * h1_mcParz_[25]
const edm::EDGetTokenT< reco::SuperClusterCollection > scToken_
static const int MAX_IPHI
Definition: EBDetId.h:137
ZeePlots * myZeePlots_
int zside() const
Definition: EEDetId.h:71
double fEtaEndcapGood(double scEta) const
void writeEleHistograms()
Definition: ZeePlots.cc:240
unsigned int etaBins_
edm::ParameterSet theParameterSet
const ZIterativeAlgorithmWithFitPlots * getHistos() const
double fEtaBarrelGood(double scEta) const
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:57
static bool validDetId(int crystal_ix, int crystal_iy, int iz)
Definition: EEDetId.h:248
const std::vector< int > & getOptimizedIterations() const
TH2F * h2_miscalRecalEB_
TH1F * h_ESCcorrEtrue_[25]
std::pair< DetId, double > getHottestDetId(const std::vector< std::pair< DetId, float > > &mySCRecHits, const EBRecHitCollection *ebhits, const EERecHitCollection *eehits)
double loopArray[50]
TH2F * h2_xtalMiscalibCoeffEndcapPlus_
static const int MAX_IETA
Definition: EBDetId.h:136
const edm::EDGetTokenT< EBRecHitCollection > ebRecHitToken_
Status endOfLoop(const edm::EventSetup &, unsigned int iLoop) override
Called at end of loop.
const std::string mcProducer_
static std::vector< DetId > getDetIdsInECAL()
DetId id() const
get the id
Definition: EcalRecHit.h:77
bool isValid() const
Definition: HandleBase.h:70
float initCalibCoeff[250]
const std::vector< float > & getOptimizedCoefficients() const
int SILVER_ELECTRONS_IN_ENDCAP
TH1F * h1_reco_ZMassCorrBB_
const edm::EDGetTokenT< reco::GsfElectronCollection > gsfElectronToken_
const CaloClusterPtr & seed() const
seed BasicCluster
Definition: SuperCluster.h:79
void writeLine(EBDetId const &, float)
iterator find(key_type k)
virtual std::shared_ptr< EcalIntercalibConstants > produceEcalIntercalibConstants(const EcalIntercalibConstantsRcd &iRecord)
Produce Ecal interCalibrations.
fixed size matrix
HLT enums.
TH1F * h1_reco_ZMassCorrEE_
TH2F * h2_xtalMiscalibCoeffBarrel_
TH2F * h2_chi2_[25]
TH2F * h2_coeffVsLoop_
const std::vector< float > & getOptimizedChiSquare() const
std::string barrelfile_
TH1F * h1_weightSumMeanBarrel_
TH2F * h_ESCEtrueVsEta_[25]
TH2F * h2_zMassVsLoop_
void bookZHistograms()
Definition: ZeePlots.cc:74
float DPhi(float phi1, float phi2)
Definition: L1Clustering.h:127
static const int IY_MAX
Definition: EEDetId.h:302
TH2F * h2_iterations_[25]
void bookHLTHistograms()
Definition: ZeePlots.cc:270
ZeeCalibration(const edm::ParameterSet &iConfig)
Constructor.
float NewCalibCoeff[250]
float calibCoeffError[250]
TH2F * h2_xtalRecalibCoeffEndcapPlus_[25]
float energy() const
Definition: EcalRecHit.h:68
int CRACK_ELECTRONS_IN_ENDCAP
int CRACK_ELECTRONS_IN_BARREL
double fEtaEndcapBad(double scEta) const
TH2F * h2_fEtaEndcapGood_
TH2F * h2_xtalRecalibCoeffEndcapMinus_[25]
void bookEleMCHistograms()
Definition: ZeePlots.cc:156
double phi() const final
momentum azimuthal angle
double getEtaCorrection(const reco::GsfElectron *)
double preshowerEnergy() const
energy deposited in preshower
Definition: SuperCluster.h:63
static float cosThetaElectrons_SC(const std::pair< calib::CalibElectron *, calib::CalibElectron *> &aZCandidate, float ele1EnergyCorrection, float ele2EnergyCorrection)
bool aHLTResults[200]
TH2F * h2_fEtaBarrelGood_
TH1F * h_eleEffEta_[2]
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
TH1F * h1_borderElectronClassification_
SuperClusterRef superCluster() const override
reference to a SuperCluster
Definition: GsfElectron.h:155
TH2F * h2_xtalMiscalibCoeffEndcapMinus_
TH2F * h2_zWidthVsLoop_
void endOfJob() override
Called at end of job.
double sigmaErrorArray[50]
TH2F * h2_residualSigma_
int SHOWER_ELECTRONS_IN_BARREL
ib
Definition: cuy.py:661
int iy() const
Definition: EEDetId.h:83
int TOTAL_ELECTRONS_IN_BARREL
TH1F * h1_weightSumMeanEndcap_
double eta() const final
momentum pseudorapidity
std::shared_ptr< EcalIntercalibConstants > ical