CMS 3D CMS Logo

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