CMS 3D CMS Logo

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