EcalCalibAlgos/src/ZeeCalibration.cc

0001 #include <fstream>
0002 #include <iostream>
0003 #include <utility>
0004
0005 #include <TBranch.h>
0006 #include <TCanvas.h>
0007 #include <TF1.h>
0008 #include <TProfile.h>
0009 #include <TRandom.h>
0010
0011 #include <CLHEP/Vector/LorentzVector.h>
0012
0013 #include "CalibCalorimetry/CaloMiscalibTools/interface/CaloMiscalibMapEcal.h"
0014 #include "CalibCalorimetry/CaloMiscalibTools/interface/MiscalibReaderFromXMLEcalBarrel.h"
0015 #include "CalibCalorimetry/CaloMiscalibTools/interface/MiscalibReaderFromXMLEcalEndcap.h"
0016 #include "Calibration/EcalCalibAlgos/interface/ZeeCalibration.h"
0017 #include "Calibration/EcalCalibAlgos/interface/ZeeKinematicTools.h"
0018 #include "Calibration/Tools/interface/CalibrationCluster.h"
0019 #include "Calibration/Tools/interface/EcalIndexingTools.h"
0020 #include "Calibration/Tools/interface/EcalRingCalibrationTools.h"
0021 #include "Calibration/Tools/interface/HouseholderDecomposition.h"
0022 #include "Calibration/Tools/interface/MinL3Algorithm.h"
0023 #include "Calibration/Tools/interface/calibXMLwriter.h"
0024 #include "CondFormats/DataRecord/interface/EcalIntercalibConstantsRcd.h"
0025 #include "CondFormats/EcalObjects/interface/EcalIntercalibConstants.h"
0026 #include "DataFormats/CaloRecHit/interface/CaloRecHit.h"
0027 #include "DataFormats/Common/interface/Handle.h"
0028 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0029 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0030 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0031 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0032 #include "FWCore/Framework/interface/TriggerNamesService.h"
0033 #include "FWCore/ServiceRegistry/interface/Service.h"
0034 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0035
0036 #define MZ 91.1876
0037
0038 #define DEBUG 1
0039
0040 ZeeCalibration::ZeeCalibration(const edm::ParameterSet& iConfig)
0041     : hlTriggerResults_(iConfig.getParameter<edm::InputTag>("HLTriggerResults")),
0042       mcProducer_(iConfig.getUntrackedParameter<std::string>("mcProducer", "")),
0043       rechitProducer_(iConfig.getParameter<std::string>("rechitProducer")),
0044       rechitCollection_(iConfig.getParameter<std::string>("rechitCollection")),
0045       erechitProducer_(iConfig.getParameter<std::string>("erechitProducer")),
0046       erechitCollection_(iConfig.getParameter<std::string>("erechitCollection")),
0047       scProducer_(iConfig.getParameter<std::string>("scProducer")),
0048       scCollection_(iConfig.getParameter<std::string>("scCollection")),
0049       scIslandProducer_(iConfig.getParameter<std::string>("scIslandProducer")),
0050       scIslandCollection_(iConfig.getParameter<std::string>("scIslandCollection")),
0051       electronProducer_(iConfig.getParameter<std::string>("electronProducer")),
0052       electronCollection_(iConfig.getParameter<std::string>("electronCollection")),
0053       trigResultsToken_(consumes<edm::TriggerResults>(hlTriggerResults_)),
0054       hepMCToken_(consumes<edm::HepMCProduct>(edm::InputTag(mcProducer_))),
0055       ebRecHitToken_(consumes<EBRecHitCollection>(edm::InputTag(rechitProducer_, rechitCollection_))),
0056       eeRecHitToken_(consumes<EERecHitCollection>(edm::InputTag(erechitProducer_, erechitCollection_))),
0057       scToken_(consumes<reco::SuperClusterCollection>(edm::InputTag(scProducer_, scCollection_))),
0058       islandSCToken_(consumes<reco::SuperClusterCollection>(edm::InputTag(scIslandProducer_, scIslandCollection_))),
0059       gsfElectronToken_(consumes<reco::GsfElectronCollection>(edm::InputTag(electronProducer_, electronCollection_))),
0060       geometryToken_(esConsumes()) {
0061 #ifdef DEBUG
0062   std::cout << "[ZeeCalibration] Starting the ctor" << std::endl;
0063 #endif
0064
0065   theMaxLoops = iConfig.getUntrackedParameter<unsigned int>("maxLoops", 0);
0066
0067   wantEtaCorrection_ = iConfig.getUntrackedParameter<bool>("wantEtaCorrection", true);
0068
0069   outputFileName_ = iConfig.getParameter<std::string>("outputFile");
0070
0071   minInvMassCut_ = iConfig.getUntrackedParameter<double>("minInvMassCut", 70.);
0072   maxInvMassCut_ = iConfig.getUntrackedParameter<double>("maxInvMassCut", 110.);
0073
0074   calibMode_ = iConfig.getUntrackedParameter<std::string>("ZCalib_CalibType");
0075
0076   outputFile_ = TFile::Open(outputFileName_.c_str(), "RECREATE");  // open output file to store histograms
0077
0078   myTree = new TTree("myTree", "myTree");
0079   //  myTree->Branch("zMass","zMass", &mass);
0080   myTree->Branch("zMass", &mass4tree, "mass/F");
0081   myTree->Branch("zMassDiff", &massDiff4tree, "massDiff/F");
0082
0083   barrelfile_ = iConfig.getUntrackedParameter<std::string>("initialMiscalibrationBarrel", "");
0084   endcapfile_ = iConfig.getUntrackedParameter<std::string>("initialMiscalibrationEndcap", "");
0085
0086   electronSelection_ =
0087       iConfig.getUntrackedParameter<unsigned int>("electronSelection", 0);  //option for electron selection
0088
0089   etaBins_ = iConfig.getUntrackedParameter<unsigned int>("etaBins", 10);
0090   etBins_ = iConfig.getUntrackedParameter<unsigned int>("etBins", 10);
0091
0092   etaMin_ = iConfig.getUntrackedParameter<double>("etaMin", 0.);
0093   etMin_ = iConfig.getUntrackedParameter<double>("etMin", 0.);
0094   etaMax_ = iConfig.getUntrackedParameter<double>("etaMax", 3.);
0095   etMax_ = iConfig.getUntrackedParameter<double>("etMax", 100.);
0096
0097   //  new ZeePlots("zeePlots.root");
0098   //  ZeePlots->bookHistos();
0099
0100   //ZeeCalibrationPLots("zeeCalibPlots");
0101   //ZeecaPlots->bookHistos(maxsIter);
0102
0103   theParameterSet = iConfig;
0104   EcalIndexingTools* myIndexTool = nullptr;
0105
0106   myIndexTool = EcalIndexingTools::getInstance();
0107
0108   myIndexTool->setBinRange(etaBins_, etaMin_, etaMax_, etBins_, etMin_, etMax_);
0109
0110   //creating the algorithm
0111   theAlgorithm_ = new ZIterativeAlgorithmWithFit(iConfig);
0112
0113   // Tell the framework what data is being produced
0114   //setWhatProduced(this);
0115   setWhatProduced(this, &ZeeCalibration::produceEcalIntercalibConstants);
0116   findingRecord<EcalIntercalibConstantsRcd>();
0117
0118   for (int i = 0; i < 50; i++) {
0119     coefficientDistanceAtIteration[i] = -1.;
0120     loopArray[i] = -1.;
0121     sigmaArray[i] = -1.;
0122     sigmaErrorArray[i] = -1.;
0123   }
0124
0125 #ifdef DEBUG
0126   std::cout << "[ZeeCalibration] Done with the ctor" << std::endl;
0127 #endif
0128 }
0129
0130 ZeeCalibration::~ZeeCalibration() {
0131   //   if (theAlgorithm_)
0132   //     delete theAlgorithm_;
0133 }
0134
0135 //_____________________________________________________________________________
0136 // Produce EcalIntercalibConstants
0137 std::shared_ptr<EcalIntercalibConstants> ZeeCalibration::produceEcalIntercalibConstants(
0138     const EcalIntercalibConstantsRcd& iRecord) {
0139   std::cout << "@SUB=ZeeCalibration::produceEcalIntercalibConstants" << std::endl;
0140   return ical;
0141 }
0142
0143 void ZeeCalibration::beginOfJob() { isfirstcall_ = true; }
0144
0145 //========================================================================
0146 void ZeeCalibration::endOfJob() {
0147   printStatistics();
0148
0149   if (calibMode_ != "ETA_ET_MODE") {
0150     ///if not ETA_ET MODE - begin
0151
0152     //Writing out calibration coefficients
0153     calibXMLwriter* barrelWriter = new calibXMLwriter(EcalBarrel);
0154     for (int ieta = -EBDetId::MAX_IETA; ieta <= EBDetId::MAX_IETA; ++ieta) {
0155       if (ieta == 0)
0156         continue;
0157       for (int iphi = EBDetId::MIN_IPHI; iphi <= EBDetId::MAX_IPHI; ++iphi) {
0158         // make an EBDetId since we need EBDetId::rawId() to be used as the key for the pedestals
0159         if (EBDetId::validDetId(ieta, iphi)) {
0160           EBDetId ebid(ieta, iphi);
0161           barrelWriter->writeLine(ebid, *(ical->getMap().find(ebid.rawId())));
0162         }
0163       }
0164     }
0165
0166     calibXMLwriter* endcapWriter = new calibXMLwriter(EcalEndcap);
0167     for (int iX = EEDetId::IX_MIN; iX <= EEDetId::IX_MAX; ++iX) {
0168       for (int iY = EEDetId::IY_MIN; iY <= EEDetId::IY_MAX; ++iY) {
0169         // make an EEDetId since we need EEDetId::rawId() to be used as the key for the pedestals
0170         if (EEDetId::validDetId(iX, iY, 1)) {
0171           EEDetId eeid(iX, iY, 1);
0172           endcapWriter->writeLine(eeid, *(ical->getMap().find(eeid.rawId())));
0173         }
0174         if (EEDetId::validDetId(iX, iY, -1)) {
0175           EEDetId eeid(iX, iY, -1);
0176           endcapWriter->writeLine(eeid, *(ical->getMap().find(eeid.rawId())));
0177         }
0178       }
0179     }
0180
0181   }  ///if not ETA_ET MODE - end
0182
0183   std::cout << "Writing  histos..." << std::endl;
0184   outputFile_->cd();
0185
0186   //  zeeplts->Write();
0187
0188   h1_eventsBeforeEWKSelection_->Write();
0189   h1_eventsAfterEWKSelection_->Write();
0190
0191   h1_eventsBeforeBorderSelection_->Write();
0192   h1_eventsAfterBorderSelection_->Write();
0193
0194   h1_borderElectronClassification_->Write();
0195
0196   h2_xtalMiscalibCoeffBarrel_->Write();
0197   h2_xtalMiscalibCoeffEndcapMinus_->Write();
0198   h2_xtalMiscalibCoeffEndcapPlus_->Write();
0199
0200   h1_electronCosTheta_SC_->Write();
0201   h1_electronCosTheta_TK_->Write();
0202   h1_electronCosTheta_SC_TK_->Write();
0203
0204   h1_zMassResol_->Write();
0205   h1_zEtaResol_->Write();
0206   h1_zPhiResol_->Write();
0207   h1_eleEtaResol_->Write();
0208   h1_elePhiResol_->Write();
0209   h1_seedOverSC_->Write();
0210   h1_preshowerOverSC_->Write();
0211
0212   for (unsigned int i = 0; i < 25; i++) {
0213     if (i < theMaxLoops) {
0214       h_ESCEtrueVsEta_[i]->Write();
0215       h_ESCEtrue_[i]->Write();
0216
0217       h_ESCcorrEtrueVsEta_[i]->Write();
0218       h_ESCcorrEtrue_[i]->Write();
0219
0220       h2_chi2_[i]->Write();
0221       h2_iterations_[i]->Write();
0222
0223       //      h_DiffZMassDistr_[i]->Write();
0224
0225       //h_ZMassDistr_[i]->Write();
0226     }
0227   }
0228
0229   h2_fEtaBarrelGood_->Write();
0230   h2_fEtaBarrelBad_->Write();
0231   h2_fEtaEndcapGood_->Write();
0232   h2_fEtaEndcapBad_->Write();
0233   h1_eleClasses_->Write();
0234
0235   h_eleEffEta_[0]->Write();
0236   h_eleEffPhi_[0]->Write();
0237   h_eleEffPt_[0]->Write();
0238
0239   h_eleEffEta_[1]->Write();
0240   h_eleEffPhi_[1]->Write();
0241   h_eleEffPt_[1]->Write();
0242
0243   int j = 0;
0244
0245   int flag = 0;
0246
0247   Double_t mean[25] = {0.};
0248   Double_t num[25] = {0.};
0249   Double_t meanErr[25] = {0.};
0250   Float_t rms[25] = {0.};
0251   Float_t tempRms[10][25];
0252
0253   for (int ia = 0; ia < 10; ia++) {
0254     for (int ib = 0; ib < 25; ib++) {
0255       tempRms[ia][ib] = 0.;
0256     }
0257   }
0258
0259   int aa = 0;
0260
0261   for (int k = 0; k < theAlgorithm_->getNumberOfChannels(); k++) {
0262     //////grouped
0263     bool isNearCrack = false;
0264
0265     if (calibMode_ == "RING") {
0266       isNearCrack = (abs(ringNumberCorrector(k)) == 1 || abs(ringNumberCorrector(k)) == 25 ||
0267                      abs(ringNumberCorrector(k)) == 26 || abs(ringNumberCorrector(k)) == 45 ||
0268                      abs(ringNumberCorrector(k)) == 46 || abs(ringNumberCorrector(k)) == 65 ||
0269                      abs(ringNumberCorrector(k)) == 66 || abs(ringNumberCorrector(k)) == 85 ||
0270                      abs(ringNumberCorrector(k)) == 86 || abs(ringNumberCorrector(k)) == 124);
0271     }
0272
0273     if (k < 85) {
0274       if ((k + 1) % 5 != 0) {
0275         if (!isNearCrack) {
0276           mean[j] += calibCoeff[k];
0277           mean[j] += calibCoeff[169 - k];
0278
0279           num[j] += 2.;
0280
0281           //meanErr[j]+= calibCoeffError[k];
0282           //meanErr[j]+= calibCoeffError[169 - k];
0283
0284           meanErr[j] += 1. / pow(calibCoeffError[k], 2);
0285           meanErr[j] += 1. / pow(calibCoeffError[169 - k], 2);
0286
0287           tempRms[aa][j] += calibCoeff[k];
0288           aa++;
0289           tempRms[aa][j] += calibCoeff[169 - k];
0290           aa++;
0291         }
0292       }
0293
0294       else {
0295         if (!isNearCrack) {
0296           mean[j] += calibCoeff[k];
0297           mean[j] += calibCoeff[169 - k];
0298
0299           num[j] += 2.;
0300
0301           //meanErr[j]+= calibCoeffError[k];
0302           //meanErr[j]+= calibCoeffError[169 - k];
0303
0304           meanErr[j] += 1. / pow(calibCoeffError[k], 2);
0305           meanErr[j] += 1. / pow(calibCoeffError[169 - k], 2);
0306
0307           tempRms[aa][j] += calibCoeff[k];
0308           aa++;
0309           tempRms[aa][j] += calibCoeff[169 - k];
0310           aa++;
0311         }
0312         j++;
0313         aa = 0;
0314       }
0315     }
0316     //EE begin
0317
0318     if (k >= 170 && k <= 204) {
0319       if (flag < 4) {
0320         //make groups of 5 Xtals in #eta
0321         mean[j] += calibCoeff[k] / 10.;
0322         mean[j] += calibCoeff[k + 39] / 10.;
0323
0324         meanErr[j] += calibCoeffError[k] / 30.;
0325         meanErr[j] += calibCoeffError[k + 39] / 30.;
0326
0327         tempRms[aa][j] += calibCoeff[k];
0328         aa++;
0329         tempRms[aa][j] += calibCoeff[k + 39];
0330         aa++;
0331
0332         flag++;
0333       }
0334
0335       else if (flag == 4) {
0336         //make groups of 5 Xtals in #eta
0337         mean[j] += calibCoeff[k] / 10.;
0338         mean[j] += calibCoeff[k + 39] / 10.;
0339
0340         meanErr[j] += calibCoeffError[k] / 30.;
0341         meanErr[j] += calibCoeffError[k + 39] / 30.;
0342
0343         tempRms[aa][j] += calibCoeff[k];
0344         aa++;
0345         tempRms[aa][j] += calibCoeff[k + 39];
0346         aa++;
0347
0348         flag = 0;
0349         //    std::cout<<" index(>85) "<<k<<" j is "<<j<<" mean[j] is "<<mean[j]<<std::endl;
0350         j++;
0351         aa = 0;
0352       }
0353     }
0354     if (k >= 205 && k <= 208) {
0355       mean[j] += calibCoeff[k] / 8.;
0356       mean[j] += calibCoeff[k + 39] / 8.;
0357
0358       meanErr[j] += calibCoeffError[k] / 30.;
0359       meanErr[j] += calibCoeffError[k + 39] / 30.;
0360
0361       tempRms[aa][j] += calibCoeff[k];
0362       aa++;
0363       tempRms[aa][j] += calibCoeff[k + 39];
0364       aa++;
0365     }
0366     //EE end
0367
0368     /*
0369       for(int jj =0; jj< 25; jj++){
0370       if(meanErr[jj] > 0.)
0371     std::cout<<" meanErr[jj] before sqrt: "<<meanErr[jj]<<std::endl;
0372
0373     meanErr[jj] = 1./sqrt( meanErr[jj] );
0374
0375     std::cout<<" meanErr[jj] after sqrt: "<<meanErr[jj]<<std::endl;
0376       }
0377       */
0378
0379     if (!isNearCrack) {
0380       h2_coeffVsEta_->Fill(ringNumberCorrector(k), calibCoeff[k]);
0381       h2_miscalRecal_->Fill(initCalibCoeff[k], 1. / calibCoeff[k]);
0382       h1_mc_->Fill(initCalibCoeff[k] * calibCoeff[k] - 1.);
0383
0384       if (k < 170) {
0385         h2_miscalRecalEB_->Fill(initCalibCoeff[k], 1. / calibCoeff[k]);
0386         h1_mcEB_->Fill(initCalibCoeff[k] * calibCoeff[k] - 1.);
0387       }
0388
0389       if (k >= 170) {
0390         h2_miscalRecalEE_->Fill(initCalibCoeff[k], 1. / calibCoeff[k]);
0391         h1_mcEE_->Fill(initCalibCoeff[k] * calibCoeff[k] - 1.);
0392       }
0393     }
0394   }
0395
0396   for (int ic = 0; ic < 17; ic++) {
0397     mean[ic] = mean[ic] / num[ic];  //find mean of recalib coeff on group of rings
0398     //meanErr[ic] = meanErr[ic] / ( sqrt( num[ic] ) * num[ic] ); //find mean of recalib coeff on group of rings
0399     meanErr[ic] = 1. / sqrt(meanErr[ic]);  //find mean of recalib coeff on group of rings
0400   }
0401
0402   //build array of RMS
0403   for (int ic = 0; ic < 25; ic++) {
0404     for (int id = 0; id < 10; id++) {
0405       if (tempRms[id][ic] > 0.) {
0406         rms[ic] += (tempRms[id][ic] - mean[j]) * (tempRms[id][ic] - mean[j]);
0407       }
0408     }
0409     rms[ic] /= 10.;  //this is approximate
0410     rms[ic] = sqrt(rms[ic]);
0411   }
0412
0413   //build array of RMS
0414
0415   Double_t xtalEta[25] = {1.4425, 1.3567, 1.2711, 1.1855, 1.10,   1.01,   0.92,   0.83,   0.7468,
0416                           0.6612, 0.5756, 0.4897, 0.3985, 0.3117, 0.2250, 0.1384, 0.0487, 1.546,
0417                           1.651,  1.771,  1.908,  2.071,  2.267,  2.516,  2.8};
0418
0419   Double_t zero[25] = {0.026};  //interval/sqrt(12)
0420
0421   for (int j = 0; j < 25; j++)
0422     h2_coeffVsEtaGrouped_->Fill(xtalEta[j], mean[j]);
0423
0424   //  for(int sho = 0; sho <25; sho++)
0425   //cout<<"xtalEta[j] "<< xtalEta[sho]<<" mean[j]  "<<mean[sho]<<"  err[j] "<<meanErr[sho]<<std::endl;
0426
0427   TProfile* px = h2_coeffVsEta_->ProfileX("coeffVsEtaProfile");
0428   px->SetXTitle("Eta channel");
0429   px->SetYTitle("recalibCoeff");
0430   px->Write();
0431
0432   h2_coeffVsEta_->Write();
0433   h2_coeffVsEtaGrouped_->Write();
0434   h2_zMassVsLoop_->Write();
0435   h2_zMassDiffVsLoop_->Write();
0436   h2_zWidthVsLoop_->Write();
0437   h2_coeffVsLoop_->Write();
0438   h2_miscalRecal_->Write();
0439   h1_mc_->Write();
0440   h2_miscalRecalEB_->Write();
0441   h1_mcEB_->Write();
0442   h2_miscalRecalEE_->Write();
0443   h1_mcEE_->Write();
0444
0445   h2_residualSigma_->Write();
0446
0447   const ZIterativeAlgorithmWithFit::ZIterativeAlgorithmWithFitPlots* algoHistos = theAlgorithm_->getHistos();
0448
0449   double weightSumMeanBarrel = 0.;
0450   double weightSumMeanEndcap = 0.;
0451
0452   for (int iIteration = 0; iIteration < theAlgorithm_->getNumberOfIterations(); iIteration++)
0453     for (int iChannel = 0; iChannel < theAlgorithm_->getNumberOfChannels(); iChannel++) {
0454       if (iIteration == (theAlgorithm_->getNumberOfIterations() - 1)) {
0455         if (iChannel < 170)
0456           weightSumMeanBarrel += algoHistos->weightedRescaleFactor[iIteration][iChannel]->Integral() / 170.;
0457
0458         if (iChannel >= 170)
0459           weightSumMeanEndcap += algoHistos->weightedRescaleFactor[iIteration][iChannel]->Integral() / 78.;
0460
0461         h1_occupancyVsEta_->Fill((Double_t)ringNumberCorrector(iChannel),
0462                                  algoHistos->weightedRescaleFactor[iIteration][iChannel]->Integral());
0463
0464         h1_occupancy_->Fill(algoHistos->weightedRescaleFactor[iIteration][iChannel]->Integral());
0465
0466         if (iChannel < 170)
0467           h1_occupancyBarrel_->Fill(algoHistos->weightedRescaleFactor[iIteration][iChannel]->Integral());
0468
0469         if (iChannel >= 170)
0470           h1_occupancyEndcap_->Fill(algoHistos->weightedRescaleFactor[iIteration][iChannel]->Integral());
0471
0472 #ifdef DEBUG
0473         std::cout << "Writing weighted integral for channel " << ringNumberCorrector(iChannel) << " ,value "
0474                   << algoHistos->weightedRescaleFactor[iIteration][iChannel]->Integral() << std::endl;
0475 #endif
0476       }
0477     }
0478
0479   //  std::cout<<"Done! Closing output file... "<<std::endl;
0480
0481   h1_weightSumMeanBarrel_->Fill(weightSumMeanBarrel);
0482   h1_weightSumMeanEndcap_->Fill(weightSumMeanEndcap);
0483
0484   std::cout << "Weight sum mean on channels in Barrel is :" << weightSumMeanBarrel << std::endl;
0485   std::cout << "Weight sum mean on channels in Endcap is :" << weightSumMeanEndcap << std::endl;
0486
0487   h1_weightSumMeanBarrel_->Write();
0488   h1_weightSumMeanEndcap_->Write();
0489
0490   h1_occupancyVsEta_->Write();
0491   h1_occupancy_->Write();
0492   h1_occupancyBarrel_->Write();
0493   h1_occupancyEndcap_->Write();
0494
0495   myTree->Write();
0496
0497   TGraphErrors* graph = new TGraphErrors(25, xtalEta, mean, zero, meanErr);
0498   graph->Draw("APL");
0499   graph->Write();
0500
0501   double zero50[50] = {0.};
0502
0503   TGraphErrors* residualSigmaGraph = new TGraphErrors(50, loopArray, sigmaArray, zero50, sigmaErrorArray);
0504   residualSigmaGraph->SetName("residualSigmaGraph");
0505   residualSigmaGraph->Draw("APL");
0506   residualSigmaGraph->Write();
0507
0508   TGraphErrors* coefficientDistanceAtIterationGraph =
0509       new TGraphErrors(50, loopArray, coefficientDistanceAtIteration, zero50, zero50);
0510   coefficientDistanceAtIterationGraph->SetName("coefficientDistanceAtIterationGraph");
0511   coefficientDistanceAtIterationGraph->Draw("APL");
0512   coefficientDistanceAtIterationGraph->Write();
0513
0514   Float_t noError[250] = {0.};
0515
0516   Float_t ringInd[250];
0517   for (int i = 0; i < 250; i++)
0518     ringInd[i] = ringNumberCorrector(i);
0519
0520   TGraphErrors* graphCoeff =
0521       new TGraphErrors(theAlgorithm_->getNumberOfChannels(), ringInd, calibCoeff, noError, calibCoeffError);
0522   graphCoeff->SetName("graphCoeff");
0523   graphCoeff->Draw("APL");
0524   graphCoeff->Write();
0525
0526   //   outputFile_->Write();//this automatically writes all histos on file
0527
0528   h1_ZCandMult_->Write();
0529   h1_reco_ZMass_->Write();
0530
0531   h1_reco_ZMassCorr_->Write();
0532   h1_reco_ZMassCorrBB_->Write();
0533   h1_reco_ZMassCorrEE_->Write();
0534
0535   outputFile_->Close();
0536
0537   myZeePlots_->writeEleHistograms();
0538   myZeePlots_->writeMCEleHistograms();
0539   myZeePlots_->writeZHistograms();
0540   myZeePlots_->writeMCZHistograms();
0541
0542   // myZeeRescaleFactorPlots_ = new ZeeRescaleFactorPlots("zeeRescaleFactorPlots.root");
0543   //myZeeRescaleFactorPlots_->writeHistograms( theAlgorithm_ );
0544
0545   //  delete myZeeRescaleFactorPlots_;
0546 }
0547
0548 //_____________________________________________________________________________
0549 // Called at each event
0550 //________________________________________
0551
0552 edm::EDLooper::Status ZeeCalibration::duringLoop(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0553   using namespace edm;
0554
0555 #ifdef DEBUG
0556   std::cout << "[ZeeCalibration] Entering duringLoop" << std::endl;
0557 #endif
0558
0559   // code that used to be in beginJob
0560   if (isfirstcall_) {
0561     //inizializzare la geometria di ecal
0562     const auto& geometry = iSetup.getData(geometryToken_);
0563     EcalRingCalibrationTools::setCaloGeometry(&geometry);
0564
0565     myZeePlots_ = new ZeePlots("zeePlots.root");
0566     //  myZeeRescaleFactorPlots_ = new ZeeRescaleFactorPlots("zeeRescaleFactorPlots.root");
0567
0568     // go to *OUR* rootfile and book histograms
0569     outputFile_->cd();
0570     bookHistograms();
0571
0572     std::cout << "[ZeeCalibration::beginOfJob] Histograms booked " << std::endl;
0573
0574     loopFlag_ = 0;
0575
0576     //Read miscalibration map if requested
0577     CaloMiscalibMapEcal* miscalibMap = nullptr;
0578     if (!barrelfile_.empty() || !barrelfile_.empty()) {
0579       miscalibMap = new CaloMiscalibMapEcal();
0580       miscalibMap->prefillMap();
0581     }
0582
0583     if (!barrelfile_.empty()) {
0584       MiscalibReaderFromXMLEcalBarrel barrelreader_(*miscalibMap);
0585       barrelreader_.parseXMLMiscalibFile(barrelfile_);
0586 #ifdef DEBUG
0587       std::cout << "[ZeeCalibration::beginOfJob] Parsed EB miscal file" << std::endl;
0588 #endif
0589     }
0590
0591     if (!endcapfile_.empty()) {
0592       MiscalibReaderFromXMLEcalEndcap endcapreader_(*miscalibMap);
0593       endcapreader_.parseXMLMiscalibFile(endcapfile_);
0594 #ifdef DEBUG
0595       std::cout << "[ZeeCalibration::beginOfJob] Parsed EE miscal file" << std::endl;
0596 #endif
0597     }
0598
0599     std::cout << "  theAlgorithm_->getNumberOfChannels() " << theAlgorithm_->getNumberOfChannels() << std::endl;
0600
0601     ////////////////////set miscalibration
0602     for (int k = 0; k < theAlgorithm_->getNumberOfChannels(); k++) {
0603       calibCoeff[k] = 1.;
0604       calibCoeffError[k] = 0.;
0605
0606       std::vector<DetId> ringIds;
0607
0608       if (calibMode_ == "RING")
0609         ringIds = EcalRingCalibrationTools::getDetIdsInRing(k);
0610
0611       if (calibMode_ == "MODULE")
0612         ringIds = EcalRingCalibrationTools::getDetIdsInModule(k);
0613
0614       if (calibMode_ == "ABS_SCALE" || calibMode_ == "ETA_ET_MODE")
0615         ringIds = EcalRingCalibrationTools::getDetIdsInECAL();
0616
0617       if (miscalibMap) {
0618         initCalibCoeff[k] = 0.;
0619         for (unsigned int iid = 0; iid < ringIds.size(); ++iid) {
0620           float miscalib = *(miscalibMap->get().getMap().find(ringIds[iid]));
0621           //          float miscalib=miscalibMap->get().getMap().find(ringIds[iid])->second; ////////AP
0622           initCalibCoeff[k] += miscalib;
0623         }
0624         initCalibCoeff[k] /= (float)ringIds.size();
0625         std::cout << k << " " << initCalibCoeff[k] << " " << ringIds.size() << std::endl;
0626       } else {
0627         initCalibCoeff[k] = 1.;
0628       }
0629     }
0630
0631     ical = std::make_shared<EcalIntercalibConstants>();
0632
0633     for (int k = 0; k < theAlgorithm_->getNumberOfChannels(); k++) {
0634       //      std::vector<DetId> ringIds = EcalRingCalibrationTools::getDetIdsInRing(k);
0635
0636       std::vector<DetId> ringIds;
0637
0638       if (calibMode_ == "RING")
0639         ringIds = EcalRingCalibrationTools::getDetIdsInRing(k);
0640
0641       if (calibMode_ == "MODULE")
0642         ringIds = EcalRingCalibrationTools::getDetIdsInModule(k);
0643
0644       if (calibMode_ == "ABS_SCALE" || calibMode_ == "ETA_ET_MODE")
0645         ringIds = EcalRingCalibrationTools::getDetIdsInECAL();
0646
0647       for (unsigned int iid = 0; iid < ringIds.size(); ++iid) {
0648         //  ical->setValue( ringIds[iid], 1. * initCalibCoeff[k] );
0649
0650         if (ringIds[iid].subdetId() == EcalBarrel) {
0651           EBDetId myEBDetId(ringIds[iid]);
0652           h2_xtalMiscalibCoeffBarrel_->SetBinContent(
0653               myEBDetId.ieta() + 86,
0654               myEBDetId.iphi(),
0655               *(miscalibMap->get().getMap().find(ringIds[iid])));  //fill TH2 with miscalibCoeff
0656         }
0657
0658         if (ringIds[iid].subdetId() == EcalEndcap) {
0659           EEDetId myEEDetId(ringIds[iid]);
0660           if (myEEDetId.zside() < 0)
0661             h2_xtalMiscalibCoeffEndcapMinus_->SetBinContent(
0662                 myEEDetId.ix(),
0663                 myEEDetId.iy(),
0664                 *(miscalibMap->get().getMap().find(ringIds[iid])));  //fill TH2 with miscalibCoeff
0665
0666           if (myEEDetId.zside() > 0)
0667             h2_xtalMiscalibCoeffEndcapPlus_->SetBinContent(
0668                 myEEDetId.ix(),
0669                 myEEDetId.iy(),
0670                 *(miscalibMap->get().getMap().find(ringIds[iid])));  //fill TH2 with miscalibCoeff
0671         }
0672
0673         ical->setValue(ringIds[iid], *(miscalibMap->get().getMap().find(ringIds[iid])));
0674       }
0675
0676       read_events = 0;
0677       init_ = false;
0678     }
0679     isfirstcall_ = false;
0680   }  // if isfirstcall
0681
0682   ////////////////////////////////////////////////////////////////////////////HLT begin
0683
0684   for (unsigned int iHLT = 0; iHLT < 200; ++iHLT) {
0685     aHLTResults[iHLT] = false;
0686   }
0687
0688 #ifdef DEBUG
0689   std::cout << "[ZeeCalibration::duringLoop] Done with initializing aHLTresults[] " << std::endl;
0690 #endif
0691
0692   edm::Handle<edm::TriggerResults> hltTriggerResultHandle;
0693   iEvent.getByToken(trigResultsToken_, hltTriggerResultHandle);
0694
0695   if (!hltTriggerResultHandle.isValid()) {
0696     //std::cout << "invalid handle for HLT TriggerResults" << std::endl;
0697   } else {
0698     hltCount = hltTriggerResultHandle->size();
0699
0700     if (loopFlag_ == 0)
0701       myZeePlots_->fillHLTInfo(hltTriggerResultHandle);
0702
0703 #ifdef DEBUG
0704     std::cout << "[ZeeCalibration::duringLoop] Done with myZeePlots_->fillHLTInfo(hltTriggerResultHandle); "
0705               << std::endl;
0706 #endif
0707
0708     for (int i = 0; i < hltCount; i++) {
0709       aHLTResults[i] = hltTriggerResultHandle->accept(i);
0710
0711       //HLT bit 32 = HLT1Electron
0712       //HLT bit 34 = HLT2Electron
0713       //HLT bit 35 = HLT2ElectronRelaxed
0714     }
0715
0716     if (!aHLTResults[32] && !aHLTResults[34] && !aHLTResults[35])
0717       return kContinue;
0718   }
0719
0720 #ifdef DEBUG
0721   std::cout << "[ZeeCalibration::duringLoop] End HLT section" << std::endl;
0722 #endif
0723
0724   //////////////////////////////////////////////////////////////////////HLT end
0725
0726   std::vector<HepMC::GenParticle*> mcEle;
0727
0728   float myGenZMass(-1);
0729
0730   if (!mcProducer_.empty()) {
0731     //DUMP GENERATED Z MASS - BEGIN
0732     Handle<HepMCProduct> hepProd;
0733     iEvent.getByToken(hepMCToken_, hepProd);
0734
0735     const HepMC::GenEvent* myGenEvent = hepProd->GetEvent();
0736
0737     if (loopFlag_ == 0)
0738       myZeePlots_->fillZMCInfo(&(*myGenEvent));
0739
0740 #ifdef DEBUG
0741     std::cout << "[ZeeCalibration::duringLoop] Done with myZeePlots_->fillZMCInfo( & (*myGenEvent) ); " << std::endl;
0742 #endif
0743
0744     for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0745          ++p) {
0746       //return a pointer to MC Z in the event
0747       if ((*p)->pdg_id() == 23 && (*p)->status() == 2) {
0748         myGenZMass = (*p)->momentum().m();
0749       }
0750     }
0751     //DUMP GENERATED Z MASS - END
0752
0753     if (loopFlag_ == 0)
0754       myZeePlots_->fillEleMCInfo(&(*myGenEvent));
0755
0756     //loop over MC positrons and find the closest MC positron in (eta,phi) phace space - begin
0757     HepMC::GenParticle MCele;
0758
0759     for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0760          ++p) {
0761       if (abs((*p)->pdg_id()) == 11) {
0762         mcEle.push_back((*p));
0763         MCele = *(*p);
0764       }
0765     }
0766
0767     if (mcEle.size() == 2 && fabs(mcEle[0]->momentum().eta()) < 2.4 && fabs(mcEle[1]->momentum().eta()) < 2.4) {
0768       NEVT++;
0769
0770       if (fabs(mcEle[0]->momentum().eta()) < 1.479 && fabs(mcEle[1]->momentum().eta()) < 1.479)
0771         MCZBB++;
0772
0773       if ((fabs(mcEle[0]->momentum().eta()) > 1.479 && fabs(mcEle[1]->momentum().eta()) < 1.479) ||
0774           (fabs(mcEle[0]->momentum().eta()) < 1.479 && fabs(mcEle[1]->momentum().eta()) > 1.479))
0775         MCZEB++;
0776
0777       if (fabs(mcEle[0]->momentum().eta()) > 1.479 && fabs(mcEle[1]->momentum().eta()) > 1.479)
0778         MCZEE++;
0779     }
0780   }
0781
0782   // Get EBRecHits
0783   Handle<EBRecHitCollection> phits;
0784   iEvent.getByToken(ebRecHitToken_, phits);
0785   if (!phits.isValid()) {
0786     std::cerr << "Error! can't get the product EBRecHitCollection " << std::endl;
0787   }
0788   const EBRecHitCollection* hits = phits.product();  // get a ptr to the product
0789
0790   // Get EERecHits
0791   Handle<EERecHitCollection> ephits;
0792   iEvent.getByToken(eeRecHitToken_, ephits);
0793   if (!ephits.isValid()) {
0794     std::cerr << "Error! can't get the product EERecHitCollection " << std::endl;
0795   }
0796   const EERecHitCollection* ehits = ephits.product();  // get a ptr to the product
0797
0798   //Get Hybrid SuperClusters
0799   Handle<reco::SuperClusterCollection> pSuperClusters;
0800   iEvent.getByToken(scToken_, pSuperClusters);
0801   if (!pSuperClusters.isValid()) {
0802     std::cerr << "Error! can't get the product SuperClusterCollection " << std::endl;
0803   }
0804   const reco::SuperClusterCollection* scCollection = pSuperClusters.product();
0805
0806 #ifdef DEBUG
0807   std::cout << "scCollection->size()" << scCollection->size() << std::endl;
0808   for (reco::SuperClusterCollection::const_iterator scIt = scCollection->begin(); scIt != scCollection->end(); scIt++) {
0809     std::cout << scIt->energy() << std::endl;
0810   }
0811 #endif
0812
0813   //Get Island SuperClusters
0814   Handle<reco::SuperClusterCollection> pIslandSuperClusters;
0815   iEvent.getByToken(islandSCToken_, pIslandSuperClusters);
0816   if (!pIslandSuperClusters.isValid()) {
0817     std::cerr << "Error! can't get the product IslandSuperClusterCollection " << std::endl;
0818   }
0819   const reco::SuperClusterCollection* scIslandCollection = pIslandSuperClusters.product();
0820
0821 #ifdef DEBUG
0822   std::cout << "scCollection->size()" << scIslandCollection->size() << std::endl;
0823 #endif
0824
0825   if ((scCollection->size() + scIslandCollection->size()) < 2)
0826     return kContinue;
0827
0828   // Get Electrons
0829   Handle<reco::GsfElectronCollection> pElectrons;
0830   iEvent.getByToken(gsfElectronToken_, pElectrons);
0831   if (!pElectrons.isValid()) {
0832     std::cerr << "Error! can't get the product ElectronCollection " << std::endl;
0833   }
0834   const reco::GsfElectronCollection* electronCollection = pElectrons.product();
0835
0836   /*
0837   //reco-mc association map
0838   std::map<HepMC::GenParticle*,const reco::PixelMatchGsfElectron*> myMCmap;
0839
0840     fillMCmap(&(*electronCollection),mcEle,myMCmap);
0841
0842     fillEleInfo(mcEle,myMCmap);
0843   */
0844
0845   if (electronCollection->size() < 2)
0846     return kContinue;
0847
0848   if (!hits && !ehits) {
0849     std::cout << "!hits" << std::endl;
0850     return kContinue;
0851   }
0852
0853   if (hits->empty() && ehits->empty()) {
0854     std::cout << "hits->size() == 0" << std::endl;
0855     return kContinue;
0856   }
0857
0858   if (!electronCollection) {
0859     std::cout << "!electronCollection" << std::endl;
0860     return kContinue;
0861   }
0862
0863   if (electronCollection->empty()) {
0864     std::cout << "electronCollection->size() == 0" << std::endl;
0865     return kContinue;
0866   }
0867
0868   ///////////////////////////////////////////////////////////////////////////////////////
0869   ///                          START HERE....
0870   ///////////////////////////////////////////////////////////////////////////////////////
0871
0872   read_events++;
0873
0874   //  std::cout << "read_events = " << read_events << std::endl;
0875
0876   ////////////////////////CHOOSE EVENTS WITH 2 OR MORE ELECTRONS ////////////////////
0877
0878 #ifdef DEBUG
0879   std::cout << " Starting with myZeePlots_->fillEleInfo(electronCollection); " << std::endl;
0880 #endif
0881
0882   if (loopFlag_ == 0)
0883     myZeePlots_->fillEleInfo(electronCollection);
0884
0885 #ifdef DEBUG
0886   std::cout << " Done with myZeePlots_->fillEleInfo(electronCollection); " << std::endl;
0887 #endif
0888
0889   //FILL an electron vector - end
0890   //###################################Electron-SC association: begin#####################################################
0891   //Filling new ElectronCollection with new SC ref and calibElectron container
0892   std::vector<calib::CalibElectron> calibElectrons;
0893   //std::map< const calib::CalibElectron* , const reco::SuperCluster* > eleScMap;
0894
0895   //#####################################Electron-SC association map: end#####################################################
0896   for (unsigned int e_it = 0; e_it != electronCollection->size(); e_it++) {
0897     calibElectrons.push_back(calib::CalibElectron(&((*electronCollection)[e_it]), hits, ehits));
0898 #ifdef DEBUG
0899     std::cout << calibElectrons.back().getRecoElectron()->superCluster()->energy() << " "
0900               << calibElectrons.back().getRecoElectron()->energy() << std::endl;
0901 #endif
0902     //      h1_recoEleEnergy_->Fill(calibElectrons.back().getRecoElectron()->superCluster()->energy());
0903   }
0904   //  if (iLoop == 0)
0905   //fillCalibElectrons(calibElectrons);
0906
0907 #ifdef DEBUG
0908   std::cout << "Filled histos" << std::endl;
0909 #endif
0910
0911   //COMBINATORY FOR Z MASS - begin
0912   std::vector<std::pair<calib::CalibElectron*, calib::CalibElectron*> > zeeCandidates;
0913   int myBestZ = -1;
0914
0915   mass = -1.;
0916   double DeltaMinvMin(5000.);
0917
0918   if (calibElectrons.size() < 2)
0919     return kContinue;
0920
0921   for (unsigned int e_it = 0; e_it != calibElectrons.size() - 1; e_it++) {
0922     for (unsigned int p_it = e_it + 1; p_it != calibElectrons.size(); p_it++) {
0923 #ifdef DEBUG
0924       std::cout << e_it << " " << calibElectrons[e_it].getRecoElectron()->charge() << " " << p_it << " "
0925                 << calibElectrons[p_it].getRecoElectron()->charge() << std::endl;
0926 #endif
0927       if (calibElectrons[e_it].getRecoElectron()->charge() * calibElectrons[p_it].getRecoElectron()->charge() != -1)
0928         continue;
0929
0930       mass = ZeeKinematicTools::calculateZMass_withTK(
0931           std::pair<calib::CalibElectron*, calib::CalibElectron*>(&(calibElectrons[e_it]), &(calibElectrons[p_it])));
0932
0933       if (mass < 0)
0934         continue;
0935
0936 #ifdef DEBUG
0937       std::cout << "#######################mass " << mass << std::endl;
0938 #endif
0939
0940       zeeCandidates.push_back(
0941           std::pair<calib::CalibElectron*, calib::CalibElectron*>(&(calibElectrons[e_it]), &(calibElectrons[p_it])));
0942       double DeltaMinv = fabs(mass - MZ);
0943
0944       if (DeltaMinv < DeltaMinvMin) {
0945         DeltaMinvMin = DeltaMinv;
0946         myBestZ = zeeCandidates.size() - 1;
0947       }
0948     }
0949   }
0950
0951   //  h_DeltaZMassDistr_[loopFlag_]->Fill( (mass-MZ) / MZ );
0952
0953   //  zeeCa->Fill(zeeCandidates);
0954   //
0955   h1_ZCandMult_->Fill(zeeCandidates.size());
0956
0957   if (zeeCandidates.empty() || myBestZ == -1)
0958     return kContinue;
0959
0960   if (loopFlag_ == 0)
0961     myZeePlots_->fillZInfo(zeeCandidates[myBestZ]);
0962
0963 #ifdef DEBUG
0964   std::cout << "Found ZCandidates " << myBestZ << std::endl;
0965 #endif
0966
0967   //  h1_zMassResol_ ->Fill(mass-myGenZMass);
0968
0969   /////////////////////////////DUMP ELECTRON CLASS
0970
0971   int class1 = zeeCandidates[myBestZ].first->getRecoElectron()->classification();
0972   int class2 = zeeCandidates[myBestZ].second->getRecoElectron()->classification();
0973
0974   h1_eleClasses_->Fill(class1);
0975   h1_eleClasses_->Fill(class2);
0976
0977   std::cout << "BEFORE " << std::endl;
0978
0979   //  myZeePlots_->fillEleClassesPlots( zeeCandidates[myBestZ].first );
0980   //myZeePlots_->fillEleClassesPlots( zeeCandidates[myBestZ].second );
0981
0982   std::cout << "AFTER " << std::endl;
0983
0984   ///////////////////////ELECTRON CLASS STATISTICS
0985
0986   if (class1 < 100)
0987     //    h1_Elec_->Fill(1);
0988     TOTAL_ELECTRONS_IN_BARREL++;
0989   if (class1 >= 100)
0990     TOTAL_ELECTRONS_IN_ENDCAP++;
0991
0992   if (class2 < 100)
0993     TOTAL_ELECTRONS_IN_BARREL++;
0994   if (class2 >= 100)
0995     TOTAL_ELECTRONS_IN_ENDCAP++;
0996
0997   if (class1 == 0)
0998     GOLDEN_ELECTRONS_IN_BARREL++;
0999   if (class1 == 100)
1000     GOLDEN_ELECTRONS_IN_ENDCAP++;
1001   if (class1 == 10 || class1 == 20)
1002     SILVER_ELECTRONS_IN_BARREL++;
1003   if (class1 == 110 || class1 == 120)
1004     SILVER_ELECTRONS_IN_ENDCAP++;
1005   if (class1 >= 30 && class1 <= 34)
1006     SHOWER_ELECTRONS_IN_BARREL++;
1007   if (class1 >= 130 && class1 <= 134)
1008     SHOWER_ELECTRONS_IN_ENDCAP++;
1009   if (class1 == 40)
1010     CRACK_ELECTRONS_IN_BARREL++;
1011   if (class1 == 140)
1012     CRACK_ELECTRONS_IN_ENDCAP++;
1013
1014   if (class2 == 0)
1015     GOLDEN_ELECTRONS_IN_BARREL++;
1016   if (class2 == 100)
1017     GOLDEN_ELECTRONS_IN_ENDCAP++;
1018   if (class2 == 10 || class2 == 20)
1019     SILVER_ELECTRONS_IN_BARREL++;
1020   if (class2 == 110 || class2 == 120)
1021     SILVER_ELECTRONS_IN_ENDCAP++;
1022   if (class2 >= 30 && class2 <= 34)
1023     SHOWER_ELECTRONS_IN_BARREL++;
1024   if (class2 >= 130 && class2 <= 134)
1025     SHOWER_ELECTRONS_IN_ENDCAP++;
1026   if (class2 == 40)
1027     CRACK_ELECTRONS_IN_BARREL++;
1028   if (class2 == 140)
1029     CRACK_ELECTRONS_IN_ENDCAP++;
1030
1031   /////////////////////////////////
1032
1033   ///////////////////////////////////////EXCLUDE ELECTRONS HAVING HOTTEST XTAL WHICH IS A BORDER XTAL - begin
1034
1035   DetId firstElehottestDetId =
1036       getHottestDetId(
1037           zeeCandidates[myBestZ].first->getRecoElectron()->superCluster()->seed()->hitsAndFractions(), hits, ehits)
1038           .first;
1039   DetId secondElehottestDetId =
1040       getHottestDetId(
1041           zeeCandidates[myBestZ].second->getRecoElectron()->superCluster()->seed()->hitsAndFractions(), hits, ehits)
1042           .first;
1043
1044   bool firstElectronIsOnModuleBorder(false);
1045   bool secondElectronIsOnModuleBorder(false);
1046
1047   h1_eventsBeforeBorderSelection_->Fill(1);
1048
1049   if (class1 < 100) {
1050     if (firstElehottestDetId.subdetId() == EcalBarrel)
1051       firstElectronIsOnModuleBorder = xtalIsOnModuleBorder(firstElehottestDetId);
1052
1053     BARREL_ELECTRONS_BEFORE_BORDER_CUT++;
1054
1055     if (firstElehottestDetId.subdetId() == EcalBarrel && !firstElectronIsOnModuleBorder)
1056       BARREL_ELECTRONS_AFTER_BORDER_CUT++;
1057   }
1058
1059   if (class2 < 100) {
1060     if (secondElehottestDetId.subdetId() == EcalBarrel)
1061       secondElectronIsOnModuleBorder = xtalIsOnModuleBorder(secondElehottestDetId);
1062
1063     BARREL_ELECTRONS_BEFORE_BORDER_CUT++;
1064
1065     if (secondElehottestDetId.subdetId() == EcalBarrel && !secondElectronIsOnModuleBorder)
1066       BARREL_ELECTRONS_AFTER_BORDER_CUT++;
1067   }
1068
1069   if (class1 < 100) {
1070     if (firstElehottestDetId.subdetId() == EcalBarrel && firstElectronIsOnModuleBorder) {
1071       h1_borderElectronClassification_->Fill(class1);
1072       return kContinue;
1073     }
1074   }
1075
1076   if (class2 < 100) {
1077     if (secondElehottestDetId.subdetId() == EcalBarrel && secondElectronIsOnModuleBorder) {
1078       h1_borderElectronClassification_->Fill(class2);
1079       return kContinue;
1080     }
1081   }
1082
1083   h1_eventsAfterBorderSelection_->Fill(1);
1084   ///////////////////////////////////////EXCLUDE ELECTRONS HAVING HOTTEST XTAL WHICH IS A BORDER XTAL - end
1085
1086   if (class1 < 100 && class2 < 100) {
1087     BBZN++;
1088     if (class1 == 0 && class2 == 0)
1089       BBZN_gg++;
1090     if (class1 < 21 && class2 < 21)
1091       BBZN_tt++;
1092     if (class1 < 21 || class2 < 21)
1093       BBZN_t0++;
1094   }
1095
1096   if (class1 >= 100 && class2 >= 100) {
1097     EEZN++;
1098     if (class1 == 100 && class2 == 100)
1099       EEZN_gg++;
1100     if (class1 < 121 && class2 < 121)
1101       EEZN_tt++;
1102     if (class1 < 121 || class2 < 121)
1103       EEZN_t0++;
1104   }
1105
1106   if ((class1 < 100 && class2 >= 100) || (class2 < 100 && class1 >= 100)) {
1107     EBZN++;
1108     if ((class1 == 0 && class2 == 100) || (class2 == 0 && class1 == 100))
1109       EBZN_gg++;
1110     if ((class1 < 21 && class2 < 121) || (class2 < 21 && class1 < 121))
1111       EBZN_tt++;
1112     if (class2 < 21 || class1 < 21 || class2 < 121 || class1 < 121)
1113       EBZN_t0++;
1114   }
1115
1116   ///////////////////////////ELECTRON SELECTION///////////////////////////////
1117
1118   if (myBestZ == -1)
1119     return kContinue;
1120
1121   bool invMassBool = ((mass > minInvMassCut_) && (mass < maxInvMassCut_));
1122
1123   bool selectionBool = false;
1124   //0 = all electrons (but no crack)
1125
1126   ////////EWK selection - begin
1127
1128   float theta1 = 2. * atan(exp(-zeeCandidates[myBestZ].first->getRecoElectron()->superCluster()->eta()));
1129   bool ET_1 = ((zeeCandidates[myBestZ].first->getRecoElectron()->superCluster()->energy() * sin(theta1)) > 20.);
1130
1131   float theta2 = 2. * atan(exp(-zeeCandidates[myBestZ].second->getRecoElectron()->superCluster()->eta()));
1132   bool ET_2 = ((zeeCandidates[myBestZ].second->getRecoElectron()->superCluster()->energy() * sin(theta2)) > 20.);
1133
1134   bool HoE_1 = (zeeCandidates[myBestZ].first->getRecoElectron()->hadronicOverEm() < 0.115);
1135   bool HoE_2 = (zeeCandidates[myBestZ].second->getRecoElectron()->hadronicOverEm() < 0.115);
1136
1137   bool DeltaPhiIn_1 = (zeeCandidates[myBestZ].first->getRecoElectron()->deltaPhiSuperClusterTrackAtVtx() < 0.090);
1138   bool DeltaPhiIn_2 = (zeeCandidates[myBestZ].second->getRecoElectron()->deltaPhiSuperClusterTrackAtVtx() < 0.090);
1139
1140   bool DeltaEtaIn_1 = (zeeCandidates[myBestZ].first->getRecoElectron()->deltaEtaSuperClusterTrackAtVtx() < 0.0090);
1141   bool DeltaEtaIn_2 = (zeeCandidates[myBestZ].second->getRecoElectron()->deltaEtaSuperClusterTrackAtVtx() < 0.0090);
1142
1143   h1_eventsBeforeEWKSelection_->Fill(1);
1144
1145   if (!(invMassBool && ET_1 && ET_2 && HoE_1 && HoE_2 && DeltaPhiIn_1 && DeltaPhiIn_2 && DeltaEtaIn_1 && DeltaEtaIn_2))
1146     return kContinue;
1147   ////////EWK selection - end
1148
1149   h1_eventsAfterEWKSelection_->Fill(1);
1150
1151   ///////////////////////////////////////EXCLUDE ELECTRONS HAVING HOTTEST XTAL WHICH IS A BORDER XTAL
1152
1153   int c1st = zeeCandidates[myBestZ].first->getRecoElectron()->classification();
1154   int c2nd = zeeCandidates[myBestZ].second->getRecoElectron()->classification();
1155   if (electronSelection_ == 0)
1156     selectionBool = (myBestZ != -1 && c1st != 40 && c1st != 40 && c2nd != 40 && c2nd != 140);
1157
1158   //1 = all electrons are Golden, BB or Narrow
1159
1160   if (electronSelection_ == 1)
1161     selectionBool =
1162         (myBestZ != -1 && (c1st == 0 || c1st == 10 || c1st == 20 || c1st == 100 || c1st == 110 || c1st == 120) &&
1163          (c2nd == 0 || c2nd == 10 || c2nd == 20 || c2nd == 100 || c2nd == 110 || c2nd == 120));
1164
1165   //2 = all electrons are Golden
1166   if (electronSelection_ == 2)
1167     selectionBool = (myBestZ != -1 && (c1st == 0 || c1st == 100) && (c2nd == 0 || c2nd == 100));
1168   //3 = all electrons are showering
1169   if (electronSelection_ == 3)
1170     selectionBool = (myBestZ != -1 && ((c1st >= 30 && c1st <= 34) || ((c1st >= 130 && c1st <= 134))) &&
1171                      ((c2nd >= 30 && c2nd <= 34) || ((c2nd >= 130 && c2nd <= 134)))
1172
1173     );
1174
1175   //4 = all Barrel electrons are Golden, BB or Narrow; take all Endcap electrons
1176
1177   if (electronSelection_ == 4)
1178     selectionBool = (myBestZ != -1 && (
1179
1180                                           ((c1st == 0 || c1st == 10 || c1st == 20) && c2nd >= 100 && c2nd != 140)
1181
1182                                           ||
1183
1184                                           ((c2nd == 0 || c2nd == 10 || c2nd == 20) && c1st >= 100 && c1st != 140)
1185
1186                                               ));
1187
1188   //5 = all Endcap electrons (but no crack)
1189
1190   if (electronSelection_ == 5)
1191     selectionBool = (myBestZ != -1 && c1st >= 100 && c2nd >= 100 && c1st != 140 && c2nd != 140);
1192
1193   //6 = all Barrel electrons (but no crack)
1194
1195   if (electronSelection_ == 6)
1196     selectionBool = (myBestZ != -1 && c1st < 100 && c2nd < 100 && c1st != 40 && c2nd != 40);
1197
1198   //7 = this eliminates the events which have 1 ele in the Barrel and 1 in the Endcap
1199
1200   if (electronSelection_ == 7)
1201     selectionBool = (myBestZ != -1 && !(c1st < 100 && c2nd >= 100) && !(c1st >= 100 && c2nd < 100));
1202
1203   float ele1EnergyCorrection(1.);
1204   float ele2EnergyCorrection(1.);
1205
1206   if (invMassBool && selectionBool && wantEtaCorrection_) {
1207     ele1EnergyCorrection = getEtaCorrection(zeeCandidates[myBestZ].first->getRecoElectron());
1208     ele2EnergyCorrection = getEtaCorrection(zeeCandidates[myBestZ].second->getRecoElectron());
1209   }
1210
1211   if (invMassBool && selectionBool) {
1212     h1_electronCosTheta_SC_->Fill(
1213         ZeeKinematicTools::cosThetaElectrons_SC(zeeCandidates[myBestZ], ele1EnergyCorrection, ele2EnergyCorrection));
1214     h1_electronCosTheta_TK_->Fill(
1215         ZeeKinematicTools::cosThetaElectrons_TK(zeeCandidates[myBestZ], ele1EnergyCorrection, ele2EnergyCorrection));
1216     h1_electronCosTheta_SC_TK_->Fill(
1217         ZeeKinematicTools::cosThetaElectrons_SC(zeeCandidates[myBestZ], ele1EnergyCorrection, ele2EnergyCorrection) /
1218             ZeeKinematicTools::cosThetaElectrons_TK(
1219                 zeeCandidates[myBestZ], ele1EnergyCorrection, ele2EnergyCorrection) -
1220         1.);
1221
1222     if (!mcProducer_.empty()) {
1223       h1_zMassResol_->Fill(mass - myGenZMass);
1224
1225       //reco-mc association map - begin
1226
1227       std::map<HepMC::GenParticle*, const reco::GsfElectron*> myMCmap;
1228
1229       std::vector<const reco::GsfElectron*> dauElectronCollection;
1230
1231       dauElectronCollection.push_back(zeeCandidates[myBestZ].first->getRecoElectron());
1232       dauElectronCollection.push_back(zeeCandidates[myBestZ].second->getRecoElectron());
1233
1234       fillMCmap(&dauElectronCollection, mcEle, myMCmap);
1235       fillEleInfo(mcEle, myMCmap);
1236       //h_DiffZMassDistr_[loopFlag_]->Fill( (mass-myGenZMass) );
1237     }
1238
1239     //PUT f(eta) IN OUR Zee ALGORITHM
1240     theAlgorithm_->addEvent(zeeCandidates[myBestZ].first,
1241                             zeeCandidates[myBestZ].second,
1242                             MZ * sqrt(ele1EnergyCorrection * ele2EnergyCorrection));
1243
1244     h1_reco_ZMass_->Fill(ZeeKinematicTools::calculateZMass_withTK(zeeCandidates[myBestZ]));
1245
1246     h1_reco_ZMassCorr_->Fill(ZeeKinematicTools::calculateZMassWithCorrectedElectrons_withTK(
1247         zeeCandidates[myBestZ], ele1EnergyCorrection, ele2EnergyCorrection));
1248
1249     int c1st = zeeCandidates[myBestZ].first->getRecoElectron()->classification();
1250     int c2nd = zeeCandidates[myBestZ].second->getRecoElectron()->classification();
1251     if (c1st < 100 && c2nd < 100)
1252       h1_reco_ZMassCorrBB_->Fill(ZeeKinematicTools::calculateZMassWithCorrectedElectrons_withTK(
1253           zeeCandidates[myBestZ], ele1EnergyCorrection, ele2EnergyCorrection));
1254
1255     if (c1st >= 100 && c2nd >= 100)
1256       h1_reco_ZMassCorrEE_->Fill(ZeeKinematicTools::calculateZMassWithCorrectedElectrons_withTK(
1257           zeeCandidates[myBestZ], ele1EnergyCorrection, ele2EnergyCorrection));
1258
1259     mass4tree = ZeeKinematicTools::calculateZMassWithCorrectedElectrons_withTK(
1260         zeeCandidates[myBestZ], ele1EnergyCorrection, ele2EnergyCorrection);
1261
1262     massDiff4tree = ZeeKinematicTools::calculateZMassWithCorrectedElectrons_withTK(
1263                         zeeCandidates[myBestZ], ele1EnergyCorrection, ele2EnergyCorrection) -
1264                     myGenZMass;
1265
1266     //      h_ZMassDistr_[loopFlag_]->Fill(ZeeKinematicTools::calculateZMassWithCorrectedElectrons_withTK(zeeCandidates[myBestZ],ele1EnergyCorrection,ele2EnergyCorrection));
1267
1268     myTree->Fill();
1269   }
1270
1271 #ifdef DEBUG
1272   std::cout << "Added event to algorithm" << std::endl;
1273 #endif
1274
1275   return kContinue;
1276 }
1277 //end of ZeeCalibration::duringLoop
1278
1279 // Called at beginning of loop
1280 void ZeeCalibration::startingNewLoop(unsigned int iLoop) {
1281   std::cout << "[ZeeCalibration] Starting loop number " << iLoop << std::endl;
1282
1283   theAlgorithm_->resetIteration();
1284
1285   resetVariables();
1286
1287   resetHistograms();
1288
1289 #ifdef DEBUG
1290   std::cout << "[ZeeCalibration] exiting from startingNewLoop" << std::endl;
1291 #endif
1292 }
1293
1294 // Called at end of loop
1295 edm::EDLooper::Status ZeeCalibration::endOfLoop(const edm::EventSetup& iSetup, unsigned int iLoop) {
1296   double par[3];
1297   double errpar[3];
1298   double zChi2;
1299   int zIters;
1300
1301   ZIterativeAlgorithmWithFit::gausfit(h1_reco_ZMass_, par, errpar, 2., 2., &zChi2, &zIters);
1302
1303   h2_zMassVsLoop_->Fill(loopFlag_, par[1]);
1304
1305   h2_zMassDiffVsLoop_->Fill(loopFlag_, (par[1] - MZ) / MZ);
1306
1307   h2_zWidthVsLoop_->Fill(loopFlag_, par[2]);
1308
1309   //////////////////FIT Z PEAK
1310
1311   std::cout << "[ZeeCalibration] Ending loop " << iLoop << std::endl;
1312   //RUN the algorithm
1313   theAlgorithm_->iterate();
1314
1315   const std::vector<float>& optimizedCoefficients = theAlgorithm_->getOptimizedCoefficients();
1316   const std::vector<float>& optimizedCoefficientsError = theAlgorithm_->getOptimizedCoefficientsError();
1317   //const std::vector<float>& weightSum = theAlgorithm_->getWeightSum();
1318   const std::vector<float>& optimizedChi2 = theAlgorithm_->getOptimizedChiSquare();
1319   const std::vector<int>& optimizedIterations = theAlgorithm_->getOptimizedIterations();
1320
1321   //#ifdef DEBUG
1322   std::cout << "Optimized coefficients " << optimizedCoefficients.size() << std::endl;
1323   //#endif
1324
1325   //  h2_coeffVsLoop_->Fill(loopFlag_, optimizedCoefficients[75]); //show the evolution of just 1 ring coefficient (well chosen...)
1326
1327   //////////////define NewCalibCoeff - begin
1328   for (unsigned int ieta = 0; ieta < optimizedCoefficients.size(); ieta++) {
1329     NewCalibCoeff[ieta] = calibCoeff[ieta] * optimizedCoefficients[ieta];
1330
1331     h2_chi2_[loopFlag_]->Fill(ringNumberCorrector(ieta), optimizedChi2[ieta]);
1332     h2_iterations_[loopFlag_]->Fill(ringNumberCorrector(ieta), optimizedIterations[ieta]);
1333   }
1334   //////////////define NewCalibCoeff - end
1335
1336   coefficientDistanceAtIteration[loopFlag_] =
1337       computeCoefficientDistanceAtIteration(calibCoeff, NewCalibCoeff, optimizedCoefficients.size());
1338
1339   std::cout << "Iteration # : " << loopFlag_ << " CoefficientDistanceAtIteration "
1340             << coefficientDistanceAtIteration[loopFlag_] << std::endl;
1341   std::cout << "size " << optimizedCoefficients.size() << std::endl;
1342
1343   for (unsigned int ieta = 0; ieta < optimizedCoefficients.size(); ieta++) {
1344     calibCoeff[ieta] *= optimizedCoefficients[ieta];
1345     calibCoeffError[ieta] =
1346         calibCoeff[ieta] * sqrt(pow(optimizedCoefficientsError[ieta] / optimizedCoefficients[ieta], 2) +
1347                                 pow(calibCoeffError[ieta] / calibCoeff[ieta], 2));
1348     //calibCoeffError[ieta] = optimizedCoefficientsError[ieta];
1349
1350 #ifdef DEBUG
1351     std::cout << ieta << " " << optimizedCoefficients[ieta] << std::endl;
1352 #endif
1353
1354     std::vector<DetId> ringIds;
1355
1356     if (calibMode_ == "RING")
1357       ringIds = EcalRingCalibrationTools::getDetIdsInRing(ieta);
1358
1359     if (calibMode_ == "MODULE")
1360       ringIds = EcalRingCalibrationTools::getDetIdsInModule(ieta);
1361
1362     if (calibMode_ == "ABS_SCALE" || calibMode_ == "ETA_ET_MODE")
1363       ringIds = EcalRingCalibrationTools::getDetIdsInECAL();
1364
1365     for (unsigned int iid = 0; iid < ringIds.size(); ++iid) {
1366       if (ringIds[iid].subdetId() == EcalBarrel) {
1367         EBDetId myEBDetId(ringIds[iid]);
1368         h2_xtalRecalibCoeffBarrel_[loopFlag_]->SetBinContent(
1369             myEBDetId.ieta() + 86,
1370             myEBDetId.iphi(),
1371             100 * (calibCoeff[ieta] * initCalibCoeff[ieta] - 1.));  //fill TH2 with recalibCoeff
1372       }
1373
1374       if (ringIds[iid].subdetId() == EcalEndcap) {
1375         EEDetId myEEDetId(ringIds[iid]);
1376         if (myEEDetId.zside() < 0)
1377           h2_xtalRecalibCoeffEndcapMinus_[loopFlag_]->SetBinContent(
1378               myEEDetId.ix(),
1379               myEEDetId.iy(),
1380               100 * (calibCoeff[ieta] * initCalibCoeff[ieta] - 1.));  //fill TH2 with recalibCoeff
1381
1382         if (myEEDetId.zside() > 0)
1383           h2_xtalRecalibCoeffEndcapPlus_[loopFlag_]->SetBinContent(
1384               myEEDetId.ix(),
1385               myEEDetId.iy(),
1386               100 * (calibCoeff[ieta] * initCalibCoeff[ieta] - 1.));  //fill TH2 with recalibCoeff
1387       }
1388
1389       ical->setValue(ringIds[iid], *(ical->getMap().find(ringIds[iid])) * optimizedCoefficients[ieta]);
1390     }
1391   }
1392
1393   /////////dump residual miscalibration at each loop
1394
1395   for (int k = 0; k < theAlgorithm_->getNumberOfChannels(); k++) {
1396     bool isNearCrack =
1397         (abs(ringNumberCorrector(k)) == 1 || abs(ringNumberCorrector(k)) == 25 || abs(ringNumberCorrector(k)) == 26 ||
1398          abs(ringNumberCorrector(k)) == 45 || abs(ringNumberCorrector(k)) == 46 || abs(ringNumberCorrector(k)) == 65 ||
1399          abs(ringNumberCorrector(k)) == 66 || abs(ringNumberCorrector(k)) == 85 || abs(ringNumberCorrector(k)) == 86 ||
1400          abs(ringNumberCorrector(k)) == 124);
1401
1402     if (!isNearCrack) {
1403       //    h2_miscalRecalParz_[iLoop]->Fill( initCalibCoeff[k], 1./calibCoeff[k] );
1404       h1_mcParz_[iLoop]->Fill(initCalibCoeff[k] * calibCoeff[k] - 1.);
1405
1406       if (k < 170) {
1407         //h2_miscalRecalEBParz_[iLoop]->Fill( initCalibCoeff[k], 1./calibCoeff[k] );
1408         h1_mcEBParz_[iLoop]->Fill(initCalibCoeff[k] * calibCoeff[k] - 1.);
1409       }
1410
1411       if (k >= 170) {
1412         //h2_miscalRecalEEParz_[iLoop]->Fill( initCalibCoeff[k], 1./calibCoeff[k] );
1413         h1_mcEEParz_[iLoop]->Fill(initCalibCoeff[k] * calibCoeff[k] - 1.);
1414       }
1415     }
1416   }
1417
1418   /////////////////////////
1419   double parResidual[3];
1420   double errparResidual[3];
1421   double zResChi2;
1422   int zResIters;
1423
1424   ZIterativeAlgorithmWithFit::gausfit(h1_mcParz_[iLoop], parResidual, errparResidual, 3., 3., &zResChi2, &zResIters);
1425   //h1_mcParz_[iLoop]->Fit("gaus");
1426
1427   h2_residualSigma_->Fill(loopFlag_ + 1, parResidual[2]);
1428   loopArray[loopFlag_] = loopFlag_ + 1;
1429   sigmaArray[loopFlag_] = parResidual[2];
1430   sigmaErrorArray[loopFlag_] = errparResidual[2];
1431
1432   std::cout << "Fit on residuals, sigma is " << parResidual[2] << " +/- " << errparResidual[2] << std::endl;
1433
1434   /////////////////////
1435   outputFile_->cd();
1436
1437   //  h2_miscalRecalParz_[iLoop]->Write();
1438   h1_mcParz_[iLoop]->Write();
1439
1440   //h2_miscalRecalEBParz_[iLoop]->Write();
1441   h1_mcEBParz_[iLoop]->Write();
1442
1443   //h2_miscalRecalEEParz_[iLoop]->Write();
1444   h1_mcEEParz_[iLoop]->Write();
1445   h2_xtalRecalibCoeffBarrel_[loopFlag_]->Write();
1446   h2_xtalRecalibCoeffEndcapPlus_[loopFlag_]->Write();
1447   h2_xtalRecalibCoeffEndcapMinus_[loopFlag_]->Write();
1448
1449   /////////dump residual miscalibration at each loop
1450
1451   loopFlag_++;
1452
1453 #ifdef DEBUG
1454   std::cout << " loopFlag_ is " << loopFlag_ << std::endl;
1455 #endif
1456
1457   if (iLoop == theMaxLoops - 1 || iLoop >= theMaxLoops)
1458     return kStop;
1459   else
1460     return kContinue;
1461 }
1462
1463 void ZeeCalibration::bookHistograms() {
1464   h1_eventsBeforeEWKSelection_ = new TH1F("h1_eventsBeforeEWKSelection", "h1_eventsBeforeEWKSelection", 5, 0, 5);
1465   h1_eventsAfterEWKSelection_ = new TH1F("h1_eventsAfterEWKSelection", "h1_eventsAfterEWKSelection", 5, 0, 5);
1466
1467   h1_eventsBeforeBorderSelection_ =
1468       new TH1F("h1_eventsBeforeBorderSelection", "h1_eventsBeforeBorderSelection", 5, 0, 5);
1469   h1_eventsAfterBorderSelection_ = new TH1F("h1_eventsAfterBorderSelection", "h1_eventsAfterBorderSelection", 5, 0, 5);
1470
1471   h1_seedOverSC_ = new TH1F("h1_seedOverSC", "h1_seedOverSC", 400, 0., 2.);
1472
1473   myZeePlots_->bookHLTHistograms();
1474
1475   h1_borderElectronClassification_ =
1476       new TH1F("h1_borderElectronClassification", "h1_borderElectronClassification", 55, -5, 50);
1477   h1_preshowerOverSC_ = new TH1F("h1_preshowerOverSC", "h1_preshowerOverSC", 400, 0., 1.);
1478
1479   h2_fEtaBarrelGood_ = new TH2F("fEtaBarrelGood", "fEtaBarrelGood", 800, -4., 4., 800, 0.8, 1.2);
1480   h2_fEtaBarrelGood_->SetXTitle("Eta");
1481   h2_fEtaBarrelGood_->SetYTitle("1/fEtaBarrelGood");
1482
1483   h2_fEtaBarrelBad_ = new TH2F("fEtaBarrelBad", "fEtaBarrelBad", 800, -4., 4., 800, 0.8, 1.2);
1484   h2_fEtaBarrelBad_->SetXTitle("Eta");
1485   h2_fEtaBarrelBad_->SetYTitle("1/fEtaBarrelBad");
1486
1487   h2_fEtaEndcapGood_ = new TH2F("fEtaEndcapGood", "fEtaEndcapGood", 800, -4., 4., 800, 0.8, 1.2);
1488   h2_fEtaEndcapGood_->SetXTitle("Eta");
1489   h2_fEtaEndcapGood_->SetYTitle("1/fEtaEndcapGood");
1490
1491   h2_fEtaEndcapBad_ = new TH2F("fEtaEndcapBad", "fEtaEndcapBad", 800, -4., 4., 800, 0.8, 1.2);
1492   h2_fEtaEndcapBad_->SetXTitle("Eta");
1493   h2_fEtaEndcapBad_->SetYTitle("1/fEtaEndcapBad");
1494
1495   for (int i = 0; i < 2; i++) {
1496     char histoName[50];
1497
1498     sprintf(histoName, "h_eleEffEta_%d", i);
1499     h_eleEffEta_[i] = new TH1F(histoName, histoName, 150, 0., 2.7);
1500     h_eleEffEta_[i]->SetXTitle("|#eta|");
1501
1502     sprintf(histoName, "h_eleEffPhi_%d", i);
1503     h_eleEffPhi_[i] = new TH1F(histoName, histoName, 400, -4., 4.);
1504     h_eleEffPhi_[i]->SetXTitle("Phi");
1505
1506     sprintf(histoName, "h_eleEffPt_%d", i);
1507     h_eleEffPt_[i] = new TH1F(histoName, histoName, 200, 0., 200.);
1508     h_eleEffPt_[i]->SetXTitle("p_{T}(GeV/c)");
1509   }
1510
1511   h2_xtalMiscalibCoeffBarrel_ =
1512       new TH2F("h2_xtalMiscalibCoeffBarrel", "h2_xtalMiscalibCoeffBarrel", 171, -85, 85, 360, 0, 360);
1513   h2_xtalMiscalibCoeffEndcapMinus_ =
1514       new TH2F("h2_xtalMiscalibCoeffEndcapMinus", "h2_xtalMiscalibCoeffEndcapMinus", 100, 0, 100, 100, 0, 100);
1515   h2_xtalMiscalibCoeffEndcapPlus_ =
1516       new TH2F("h2_xtalMiscalibCoeffEndcapPlus", "h2_xtalMiscalibCoeffEndcapPlus", 100, 0, 100, 100, 0, 100);
1517
1518   h2_xtalMiscalibCoeffBarrel_->SetXTitle("ieta");
1519   h2_xtalMiscalibCoeffBarrel_->SetYTitle("iphi");
1520
1521   h2_xtalMiscalibCoeffEndcapMinus_->SetXTitle("ix");
1522   h2_xtalMiscalibCoeffEndcapMinus_->SetYTitle("iy");
1523
1524   for (int i = 0; i < 25; i++) {
1525     char histoName[50];
1526     sprintf(histoName, "h_ESCEtrueVsEta_%d", i);
1527
1528     h_ESCEtrueVsEta_[i] = new TH2F(histoName, histoName, 150, 0., 2.7, 300, 0., 1.5);
1529     h_ESCEtrueVsEta_[i]->SetXTitle("|#eta|");
1530     h_ESCEtrueVsEta_[i]->SetYTitle("E_{SC,raw}/E_{MC}");
1531
1532     sprintf(histoName, "h_ESCEtrue_%d", i);
1533
1534     h_ESCEtrue_[i] = new TH1F(histoName, histoName, 300, 0., 1.5);
1535
1536     sprintf(histoName, "h2_chi2_%d", i);
1537     h2_chi2_[i] = new TH2F(histoName, histoName, 1000, -150, 150, 1000, -1, 5);
1538
1539     sprintf(histoName, "h2_iterations_%d", i);
1540     h2_iterations_[i] = new TH2F(histoName, histoName, 1000, -150, 150, 1000, -1, 15);
1541
1542     sprintf(histoName, "h_ESCcorrEtrueVsEta_%d", i);
1543
1544     h_ESCcorrEtrueVsEta_[i] = new TH2F(histoName, histoName, 150, 0., 2.7, 300, 0., 1.5);
1545     h_ESCcorrEtrueVsEta_[i]->SetXTitle("|#eta|");
1546     h_ESCcorrEtrueVsEta_[i]->SetYTitle("E_{SC,#eta-corr}/E_{MC}");
1547
1548     sprintf(histoName, "h_ESCcorrEtrue_%d", i);
1549
1550     h_ESCcorrEtrue_[i] = new TH1F(histoName, histoName, 300, 0., 1.5);
1551
1552     sprintf(histoName, "h2_xtalRecalibCoeffBarrel_%d", i);
1553     h2_xtalRecalibCoeffBarrel_[i] = new TH2F(histoName, histoName, 171, -85, 85, 360, 0, 360);
1554
1555     h2_xtalRecalibCoeffBarrel_[i]->SetXTitle("ieta");
1556     h2_xtalRecalibCoeffBarrel_[i]->SetYTitle("iphi");
1557
1558     sprintf(histoName, "h2_xtalRecalibCoeffEndcapMinus_%d", i);
1559     h2_xtalRecalibCoeffEndcapMinus_[i] = new TH2F(histoName, histoName, 100, 0, 100, 100, 0, 100);
1560     h2_xtalRecalibCoeffEndcapMinus_[i]->SetXTitle("ix");
1561     h2_xtalRecalibCoeffEndcapMinus_[i]->SetYTitle("iy");
1562
1563     sprintf(histoName, "h2_xtalRecalibCoeffEndcapPlus_%d", i);
1564     h2_xtalRecalibCoeffEndcapPlus_[i] = new TH2F(histoName, histoName, 100, 0, 100, 100, 0, 100);
1565     h2_xtalRecalibCoeffEndcapPlus_[i]->SetXTitle("ix");
1566     h2_xtalRecalibCoeffEndcapPlus_[i]->SetYTitle("iy");
1567   }
1568
1569   /*
1570   for (int i=0;i<15;i++)
1571     {
1572
1573       char histoName[50];
1574
1575       sprintf(histoName,"h_DiffZMassDistr_%d",i);
1576       h_DiffZMassDistr_[i] = new TH1F(histoName,histoName, 400, -20., 20.);
1577       h_DiffZMassDistr_[i]->SetXTitle("M_{Z, reco} - M_{Z, MC}");
1578       h_DiffZMassDistr_[i]->SetYTitle("events");
1579
1580       sprintf(histoName,"h_ZMassDistr_%d",i);
1581       h_ZMassDistr_[i] = new TH1F(histoName,histoName, 200, 0., 150.);
1582       h_ZMassDistr_[i]->SetXTitle("RecoZmass (GeV)");
1583       h_ZMassDistr_[i]->SetYTitle("events");
1584
1585     }
1586   */
1587
1588   h1_zMassResol_ = new TH1F("zMassResol", "zMassResol", 200, -50., 50.);
1589   h1_zMassResol_->SetXTitle("M_{Z, reco} - M_{Z, MC}");
1590   h1_zMassResol_->SetYTitle("events");
1591
1592   h1_eleEtaResol_ = new TH1F("eleEtaResol", "eleEtaResol", 100, -0.01, 0.01);
1593   h1_eleEtaResol_->SetXTitle("#eta_{reco} - #eta_{MC}");
1594   h1_eleEtaResol_->SetYTitle("events");
1595
1596   h1_electronCosTheta_TK_ = new TH1F("electronCosTheta_TK", "electronCosTheta_TK", 100, -1, 1);
1597   h1_electronCosTheta_TK_->SetXTitle("cos #theta_{12}");
1598   h1_electronCosTheta_TK_->SetYTitle("events");
1599
1600   h1_electronCosTheta_SC_ = new TH1F("electronCosTheta_SC", "electronCosTheta_SC", 100, -1, 1);
1601   h1_electronCosTheta_SC_->SetXTitle("cos #theta_{12}");
1602   h1_electronCosTheta_SC_->SetYTitle("events");
1603
1604   h1_electronCosTheta_SC_TK_ = new TH1F("electronCosTheta_SC_TK", "electronCosTheta_SC_TK", 200, -0.1, 0.1);
1605   h1_electronCosTheta_SC_TK_->SetXTitle("cos #theta_{12}^{SC}/ cos #theta_{12}^{TK} - 1");
1606   h1_electronCosTheta_SC_TK_->SetYTitle("events");
1607
1608   h1_elePhiResol_ = new TH1F("elePhiResol", "elePhiResol", 100, -0.01, 0.01);
1609   h1_elePhiResol_->SetXTitle("#phi_{reco} - #phi_{MC}");
1610   h1_elePhiResol_->SetYTitle("events");
1611
1612   h1_zEtaResol_ = new TH1F("zEtaResol", "zEtaResol", 200, -1., 1.);
1613   h1_zEtaResol_->SetXTitle("#eta_{Z, reco} - #eta_{Z, MC}");
1614   h1_zEtaResol_->SetYTitle("events");
1615
1616   h1_zPhiResol_ = new TH1F("zPhiResol", "zPhiResol", 200, -1., 1.);
1617   h1_zPhiResol_->SetXTitle("#phi_{Z, reco} - #phi_{Z, MC}");
1618   h1_zPhiResol_->SetYTitle("events");
1619
1620   h1_nEleReco_ = new TH1F("nEleReco", "Number of reco electrons", 10, -0.5, 10.5);
1621   h1_nEleReco_->SetXTitle("nEleReco");
1622   h1_nEleReco_->SetYTitle("events");
1623
1624   //  h1_occupancyVsEta_ = new TH1F("occupancyVsEta","occupancyVsEta",EcalRingCalibrationTools::N_RING_TOTAL,0,(float)EcalRingCalibrationTools::N_RING_TOTAL);
1625
1626   h1_occupancyVsEta_ = new TH1F("occupancyVsEta", "occupancyVsEta", 249, -124, 124);
1627   h1_occupancyVsEta_->SetYTitle("Weighted electron statistics");
1628   h1_occupancyVsEta_->SetXTitle("Eta channel");
1629
1630   h1_weightSumMeanBarrel_ = new TH1F("weightSumMeanBarrel", "weightSumMeanBarrel", 10000, 0, 10000);
1631   h1_weightSumMeanEndcap_ = new TH1F("weightSumMeanEndcap", "weightSumMeanEndcap", 10000, 0, 10000);
1632
1633   h1_occupancy_ = new TH1F("occupancy", "occupancy", 1000, 0, 10000);
1634   h1_occupancy_->SetXTitle("Weighted electron statistics");
1635
1636   h1_occupancyBarrel_ = new TH1F("occupancyBarrel", "occupancyBarrel", 1000, 0, 10000);
1637   h1_occupancyBarrel_->SetXTitle("Weighted electron statistics");
1638
1639   h1_occupancyEndcap_ = new TH1F("occupancyEndcap", "occupancyEndcap", 1000, 0, 10000);
1640   h1_occupancyEndcap_->SetXTitle("Weighted electron statistics");
1641
1642   h1_eleClasses_ = new TH1F("eleClasses", "eleClasses", 301, -1, 300);
1643   h1_eleClasses_->SetXTitle("classCode");
1644   h1_eleClasses_->SetYTitle("#");
1645
1646   myZeePlots_->bookZMCHistograms();
1647
1648   myZeePlots_->bookZHistograms();
1649
1650   myZeePlots_->bookEleMCHistograms();
1651
1652   myZeePlots_->bookEleHistograms();
1653
1654   h1_ZCandMult_ = new TH1F("ZCandMult", "Multiplicity of Z candidates in one event", 10, -0.5, 10.5);
1655   h1_ZCandMult_->SetXTitle("ZCandMult");
1656
1657   h1_reco_ZMass_ = new TH1F("reco_ZMass", "Inv. mass of 2 reco Electrons", 200, 0., 150.);
1658   h1_reco_ZMass_->SetXTitle("reco_ZMass (GeV)");
1659   h1_reco_ZMass_->SetYTitle("events");
1660
1661   h1_reco_ZMassCorr_ = new TH1F("reco_ZMassCorr", "Inv. mass of 2 corrected reco Electrons", 200, 0., 150.);
1662   h1_reco_ZMassCorr_->SetXTitle("reco_ZMass (GeV)");
1663   h1_reco_ZMassCorr_->SetYTitle("events");
1664
1665   h1_reco_ZMassCorrBB_ = new TH1F("reco_ZMassCorrBB", "Inv. mass of 2 corrected reco Electrons", 200, 0., 150.);
1666   h1_reco_ZMassCorrBB_->SetXTitle("reco_ZMass (GeV)");
1667   h1_reco_ZMassCorrBB_->SetYTitle("events");
1668
1669   h1_reco_ZMassCorrEE_ = new TH1F("reco_ZMassCorrEE", "Inv. mass of 2 corrected reco Electrons", 200, 0., 150.);
1670   h1_reco_ZMassCorrEE_->SetXTitle("reco_ZMass (GeV)");
1671   h1_reco_ZMassCorrEE_->SetYTitle("events");
1672
1673   //  h2_coeffVsEta_= new TH2F("h2_calibCoeffVsEta","h2_calibCoeffVsEta",EcalRingCalibrationTools::N_RING_TOTAL,0, (double)EcalRingCalibrationTools::N_RING_TOTAL, 200, 0., 2.);
1674
1675   h2_coeffVsEta_ = new TH2F("h2_calibCoeffVsEta", "h2_calibCoeffVsEta", 249, -124, 125, 200, 0., 2.);
1676   h2_coeffVsEta_->SetXTitle("Eta channel");
1677   h2_coeffVsEta_->SetYTitle("recalibCoeff");
1678
1679   h2_coeffVsEtaGrouped_ =
1680       new TH2F("h2_calibCoeffVsEtaGrouped", "h2_calibCoeffVsEtaGrouped", 200, 0., 3., 200, 0.6, 1.4);
1681   h2_coeffVsEtaGrouped_->SetXTitle("|#eta|");
1682   h2_coeffVsEtaGrouped_->SetYTitle("recalibCoeff");
1683
1684   h2_zMassVsLoop_ = new TH2F("h2_zMassVsLoop", "h2_zMassVsLoop", 1000, 0, 40, 90, 80., 95.);
1685
1686   h2_zMassDiffVsLoop_ = new TH2F("h2_zMassDiffVsLoop", "h2_zMassDiffVsLoop", 1000, 0, 40, 100, -1., 1.);
1687   h2_zMassDiffVsLoop_->SetXTitle("Iteration");
1688   h2_zMassDiffVsLoop_->SetYTitle("M_{Z, reco peak} - M_{Z, true}");
1689
1690   h2_zWidthVsLoop_ = new TH2F("h2_zWidthVsLoop", "h2_zWidthVsLoop", 1000, 0, 40, 100, 0., 10.);
1691
1692   h2_coeffVsLoop_ = new TH2F("h2_coeffVsLoop", "h2_coeffVsLoop", 1000, 0, 40, 100, 0., 2.);
1693
1694   h2_residualSigma_ = new TH2F("h2_residualSigma", "h2_residualSigma", 1000, 0, 40, 100, 0., .5);
1695
1696   h2_miscalRecal_ = new TH2F("h2_miscalRecal", "h2_miscalRecal", 500, 0., 2., 500, 0., 2.);
1697   h2_miscalRecal_->SetXTitle("initCalibCoeff");
1698   h2_miscalRecal_->SetYTitle("1/RecalibCoeff");
1699
1700   h2_miscalRecalEB_ = new TH2F("h2_miscalRecalEB", "h2_miscalRecalEB", 500, 0., 2., 500, 0., 2.);
1701   h2_miscalRecalEB_->SetXTitle("initCalibCoeff");
1702   h2_miscalRecalEB_->SetYTitle("1/RecalibCoeff");
1703
1704   h2_miscalRecalEE_ = new TH2F("h2_miscalRecalEE", "h2_miscalRecalEE", 500, 0., 2., 500, 0., 2.);
1705   h2_miscalRecalEE_->SetXTitle("initCalibCoeff");
1706   h2_miscalRecalEE_->SetYTitle("1/RecalibCoeff");
1707
1708   h1_mc_ = new TH1F("h1_residualMiscalib", "h1_residualMiscalib", 200, -0.2, 0.2);
1709   h1_mcEB_ = new TH1F("h1_residualMiscalibEB", "h1_residualMiscalibEB", 200, -0.2, 0.2);
1710   h1_mcEE_ = new TH1F("h1_residualMiscalibEE", "h1_residualMiscalibEE", 200, -0.2, 0.2);
1711
1712   for (int i = 0; i < 25; i++) {
1713     char histoName[50];
1714     /*
1715          sprintf(histoName,"h2_miscalRecalParz_%d",i);
1716          h2_miscalRecalParz_[i] = new TH2F(histoName,histoName,500, 0., 2., 500, 0., 2.);
1717          h2_miscalRecalParz_[i]->SetXTitle("initCalibCoeff");
1718          h2_miscalRecalParz_[i]->SetYTitle("1/recalibCoeff");
1719
1720          sprintf(histoName,"h2_miscalRecalEBParz_%d",i);
1721          h2_miscalRecalEBParz_[i] = new TH2F(histoName,histoName,500, 0., 2., 500, 0., 2.);
1722          h2_miscalRecalEBParz_[i]->SetXTitle("initCalibCoeff");
1723          h2_miscalRecalEBParz_[i]->SetYTitle("1/recalibCoeff");
1724
1725       sprintf(histoName,"h2_miscalRecalEEParz_%d",i);
1726       h2_miscalRecalEEParz_[i] = new TH2F(histoName,histoName,500, 0., 2., 500, 0., 2.);
1727       h2_miscalRecalEEParz_[i]->SetXTitle("initCalibCoeff");
1728       h2_miscalRecalEEParz_[i]->SetYTitle("1/recalibCoeff");
1729       */
1730
1731     sprintf(histoName, "h1_residualMiscalibParz_%d", i);
1732     h1_mcParz_[i] = new TH1F(histoName, histoName, 200, -0.2, 0.2);
1733     sprintf(histoName, "h1_residualMiscalibEBParz_%d", i);
1734     h1_mcEBParz_[i] = new TH1F(histoName, histoName, 200, -0.2, 0.2);
1735     sprintf(histoName, "h1_residualMiscalibEEParz_%d", i);
1736     h1_mcEEParz_[i] = new TH1F(histoName, histoName, 200, -0.2, 0.2);
1737   }
1738 }
1739
1740 double ZeeCalibration::fEtaBarrelBad(double scEta) const {
1741   float p0 = 1.00153e+00;
1742   float p1 = 3.29331e-02;
1743   float p2 = 1.21187e-03;
1744
1745   double x = (double)fabs(scEta);
1746
1747   return 1. / (p0 + p1 * x * x + p2 * x * x * x * x);
1748 }
1749
1750 double ZeeCalibration::fEtaEndcapGood(double scEta) const {
1751   // f(eta) for the first 3 classes (100, 110 and 120)
1752   // Ivica's new corrections 01/06
1753   float p0 = 1.06819e+00;
1754   float p1 = -1.53189e-02;
1755   float p2 = 4.01707e-04;
1756
1757   double x = (double)fabs(scEta);
1758
1759   return 1. / (p0 + p1 * x * x + p2 * x * x * x * x);
1760 }
1761
1762 double ZeeCalibration::fEtaEndcapBad(double scEta) const {
1763   float p0 = 1.17382e+00;
1764   float p1 = -6.52319e-02;
1765   float p2 = 6.26108e-03;
1766
1767   double x = (double)fabs(scEta);
1768
1769   return 1. / (p0 + p1 * x * x + p2 * x * x * x * x);
1770 }
1771
1772 double ZeeCalibration::fEtaBarrelGood(double scEta) const {
1773   float p0 = 9.99782e-01;
1774   float p1 = 1.26983e-02;
1775   float p2 = 2.16344e-03;
1776
1777   double x = (double)fabs(scEta);
1778
1779   return 1. / (p0 + p1 * x * x + p2 * x * x * x * x);
1780 }
1781
1782 //////////////////////////////////new part
1783
1784 void ZeeCalibration::fillMCmap(const std::vector<const reco::GsfElectron*>* electronCollection,
1785                                const std::vector<HepMC::GenParticle*>& mcEle,
1786                                std::map<HepMC::GenParticle*, const reco::GsfElectron*>& myMCmap) {
1787   for (unsigned int i = 0; i < mcEle.size(); i++) {
1788     float minDR = 0.1;
1789     const reco::GsfElectron* myMatchEle = nullptr;
1790     for (unsigned int j = 0; j < electronCollection->size(); j++) {
1791       float dr = EvalDR(mcEle[i]->momentum().pseudoRapidity(),
1792                         (*(*electronCollection)[j]).eta(),
1793                         mcEle[i]->momentum().phi(),
1794                         (*(*electronCollection)[j]).phi());
1795       if (dr < minDR) {
1796         myMatchEle = (*electronCollection)[j];
1797         minDR = dr;
1798       }
1799     }
1800     myMCmap.insert(std::pair<HepMC::GenParticle*, const reco::GsfElectron*>(mcEle[i], myMatchEle));
1801   }
1802 }
1803
1804 float ZeeCalibration::EvalDR(float Eta, float Eta_ref, float Phi, float Phi_ref) {
1805   if (Phi < 0)
1806     Phi = 2 * TMath::Pi() + Phi;
1807   if (Phi_ref < 0)
1808     Phi_ref = 2 * TMath::Pi() + Phi_ref;
1809   float DPhi = Phi - Phi_ref;
1810   if (fabs(DPhi) > TMath::Pi())
1811     DPhi = 2 * TMath::Pi() - fabs(DPhi);
1812
1813   float DEta = Eta - Eta_ref;
1814
1815   float DR = sqrt(DEta * DEta + DPhi * DPhi);
1816   return DR;
1817 }
1818
1819 float ZeeCalibration::EvalDPhi(float Phi, float Phi_ref) {
1820   if (Phi < 0)
1821     Phi = 2 * TMath::Pi() + Phi;
1822   if (Phi_ref < 0)
1823     Phi_ref = 2 * TMath::Pi() + Phi_ref;
1824   return (Phi - Phi_ref);
1825 }
1826
1827 void ZeeCalibration::fillEleInfo(std::vector<HepMC::GenParticle*>& mcEle,
1828                                  std::map<HepMC::GenParticle*, const reco::GsfElectron*>& associationMap) {
1829   for (unsigned int i = 0; i < mcEle.size(); i++) {
1830     h_eleEffEta_[0]->Fill(fabs(mcEle[i]->momentum().pseudoRapidity()));
1831     h_eleEffPhi_[0]->Fill(mcEle[i]->momentum().phi());
1832     h_eleEffPt_[0]->Fill(mcEle[i]->momentum().perp());
1833
1834     std::map<HepMC::GenParticle*, const reco::GsfElectron*>::const_iterator mIter = associationMap.find(mcEle[i]);
1835     if (mIter == associationMap.end())
1836       continue;
1837
1838     if ((*mIter).second) {
1839       const reco::GsfElectron* myEle = (*mIter).second;
1840
1841       h_eleEffEta_[1]->Fill(fabs(mcEle[i]->momentum().pseudoRapidity()));
1842       h_eleEffPhi_[1]->Fill(mcEle[i]->momentum().phi());
1843       h_eleEffPt_[1]->Fill(mcEle[i]->momentum().perp());
1844       h1_eleEtaResol_->Fill(myEle->eta() - mcEle[i]->momentum().eta());
1845       h1_elePhiResol_->Fill(myEle->phi() - mcEle[i]->momentum().phi());
1846
1847       const reco::SuperCluster* mySC = &(*(myEle->superCluster()));
1848       if (/*fabs(mySC->position().eta()) < 2.4*/ true) {
1849         //      if(myEle->classification()>=100)std::cout<<"mySC->preshowerEnergy()"<<mySC->preshowerEnergy()<<std::endl;
1850
1851         h_ESCEtrue_[loopFlag_]->Fill(mySC->energy() / mcEle[i]->momentum().e());
1852         h_ESCEtrueVsEta_[loopFlag_]->Fill(fabs(mySC->position().eta()), mySC->energy() / mcEle[i]->momentum().e());
1853
1854         double corrSCenergy = (mySC->energy()) / getEtaCorrection(myEle);
1855         h_ESCcorrEtrue_[loopFlag_]->Fill(corrSCenergy / mcEle[i]->momentum().e());
1856         h_ESCcorrEtrueVsEta_[loopFlag_]->Fill(fabs(mySC->position().eta()), corrSCenergy / mcEle[i]->momentum().e());
1857
1858         //        std::vector<DetId> mySCRecHits = mySC->seed()->getHitsByDetId();
1859
1860         h1_seedOverSC_->Fill(mySC->seed()->energy() / mySC->energy());
1861         h1_preshowerOverSC_->Fill(mySC->preshowerEnergy() / mySC->energy());
1862       }
1863     }
1864   }
1865 }
1866
1867 int ZeeCalibration::ringNumberCorrector(int k) {
1868   int index = -999;
1869
1870   if (calibMode_ == "RING") {
1871     if (k >= 0 && k <= 84)
1872       index = k - 85;
1873
1874     if (k >= 85 && k <= 169)
1875       index = k - 84;
1876
1877     if (k >= 170 && k <= 208)
1878       index = -k + 84;
1879
1880     if (k >= 209 && k <= 247)
1881       index = k - 123;
1882
1883   }
1884
1885   else if (calibMode_ == "MODULE") {
1886     if (k >= 0 && k <= 71)
1887       index = k - 72;
1888
1889     if (k >= 72 && k <= 143)
1890       index = k - 71;
1891   }
1892   return index;
1893 }
1894
1895 double ZeeCalibration::getEtaCorrection(const reco::GsfElectron* ele) {
1896   double correction(1.);
1897
1898   int c = ele->classification();
1899   if (c == 0 || c == 10 || c == 20)
1900     correction = fEtaBarrelGood(ele->superCluster()->eta());
1901
1902   if (c == 100 || c == 110 || c == 120)
1903     correction = fEtaEndcapGood(ele->superCluster()->eta());
1904
1905   if (c == 30 || c == 31 || c == 32 || c == 33 || c == 34)
1906     correction = fEtaBarrelBad(ele->superCluster()->eta());
1907
1908   if (c == 130 || c == 131 || c == 132 || c == 133 || c == 134)
1909     correction = fEtaEndcapBad(ele->superCluster()->eta());
1910
1911   return correction;
1912 }
1913
1914 std::pair<DetId, double> ZeeCalibration::getHottestDetId(const std::vector<std::pair<DetId, float> >& mySCRecHits,
1915                                                          const EBRecHitCollection* ebhits,
1916                                                          const EERecHitCollection* eehits) {
1917   double maxEnergy = -9999.;
1918   const EcalRecHit* hottestRecHit = nullptr;
1919
1920   std::pair<DetId, double> myPair(DetId(0), -9999.);
1921
1922   for (std::vector<std::pair<DetId, float> >::const_iterator idIt = mySCRecHits.begin(); idIt != mySCRecHits.end();
1923        idIt++) {
1924     if (idIt->first.subdetId() == EcalBarrel) {
1925       hottestRecHit = &(*(ebhits->find((*idIt).first)));
1926
1927       if (hottestRecHit == &(*(ebhits->end()))) {
1928         std::cout << "@@@@@@@@@@@@@@@@@@@@@@@@@@@ NO RECHIT FOUND SHOULD NEVER HAPPEN" << std::endl;
1929         continue;
1930       }
1931     } else if (idIt->first.subdetId() == EcalEndcap) {
1932       hottestRecHit = &(*(eehits->find((*idIt).first)));
1933       if (hottestRecHit == &(*(eehits->end()))) {
1934         std::cout << "@@@@@@@@@@@@@@@@@@@@@@@@@@@ NO RECHIT FOUND SHOULD NEVER HAPPEN" << std::endl;
1935         continue;
1936       }
1937     }
1938
1939     //std::cout<<"[getHottestDetId] hottestRecHit->energy() "<<hottestRecHit->energy()<<std::endl;
1940
1941     if (hottestRecHit && hottestRecHit->energy() > maxEnergy) {
1942       maxEnergy = hottestRecHit->energy();
1943
1944       myPair.first = hottestRecHit->id();
1945       myPair.second = maxEnergy;
1946     }
1947
1948   }  //end loop to find hottest RecHit
1949
1950   //std::cout<<"[ZeeCalibration::getHottestDetId] going to return..."<<std::endl;
1951
1952   return myPair;
1953 }
1954
1955 bool ZeeCalibration::xtalIsOnModuleBorder(EBDetId myEBDetId) {
1956   bool myBool(false);
1957
1958   short ieta = myEBDetId.ieta();
1959   short iphi = myEBDetId.iphi();
1960
1961   //  std::cout<<"[xtalIsOnModuleBorder] ieta: "<<ieta<<" iphi "<<iphi<<std::endl;
1962
1963   myBool = (abs(ieta) == 1 || abs(ieta) == 25 || abs(ieta) == 26 || abs(ieta) == 45 || abs(ieta) == 46 ||
1964             abs(ieta) == 65 || abs(ieta) == 66 || abs(ieta) == 85);
1965
1966   for (int i = 0; i < 19; i++) {
1967     if (iphi == (20 * i + 1) || iphi == 20 * i)
1968       myBool = true;
1969   }
1970
1971   return myBool;
1972 }
1973
1974 float ZeeCalibration::computeCoefficientDistanceAtIteration(float v1[250], float v2[250], int size) {
1975   float dist(0.);
1976
1977   for (int i = 0; i < size; i++) {
1978     //    std::cout<< "[ZeeCalibration::computeCoefficientDistanceAtIteration] Adding term "<<pow( v1[i]-v2[i], 2 )<<" from v1 "<<v1[i]<<" and v2 "<<v2[i]<<std::endl;
1979
1980     bool isNearCrack = false;
1981
1982     if (calibMode_ == "RING") {  //exclude non-calibrated rings from computation
1983
1984       isNearCrack = (abs(ringNumberCorrector(i)) == 1 || abs(ringNumberCorrector(i)) == 25 ||
1985                      abs(ringNumberCorrector(i)) == 26 || abs(ringNumberCorrector(i)) == 45 ||
1986                      abs(ringNumberCorrector(i)) == 46 || abs(ringNumberCorrector(i)) == 65 ||
1987                      abs(ringNumberCorrector(i)) == 66 || abs(ringNumberCorrector(i)) == 85 ||
1988                      abs(ringNumberCorrector(i)) == 86 || abs(ringNumberCorrector(i)) == 124);
1989     }
1990
1991     if (!isNearCrack)
1992       dist += pow(v1[i] - v2[i], 2);
1993   }
1994
1995   dist = sqrt(dist) / size;
1996
1997   return dist;
1998 }
1999
2000 void ZeeCalibration::resetVariables() {
2001   BBZN = 0;
2002   EBZN = 0;
2003   EEZN = 0;
2004   BBZN_gg = 0;
2005   EBZN_gg = 0;
2006   EEZN_gg = 0;
2007
2008   BBZN_tt = 0;
2009   EBZN_tt = 0;
2010   EEZN_tt = 0;
2011
2012   BBZN_t0 = 0;
2013   EBZN_t0 = 0;
2014   EEZN_t0 = 0;
2015
2016   TOTAL_ELECTRONS_IN_BARREL = 0;
2017   TOTAL_ELECTRONS_IN_ENDCAP = 0;
2018
2019   GOLDEN_ELECTRONS_IN_BARREL = 0;
2020   GOLDEN_ELECTRONS_IN_ENDCAP = 0;
2021   SILVER_ELECTRONS_IN_BARREL = 0;
2022   SILVER_ELECTRONS_IN_ENDCAP = 0;
2023   SHOWER_ELECTRONS_IN_BARREL = 0;
2024   SHOWER_ELECTRONS_IN_ENDCAP = 0;
2025   CRACK_ELECTRONS_IN_BARREL = 0;
2026   CRACK_ELECTRONS_IN_ENDCAP = 0;
2027
2028   BARREL_ELECTRONS_BEFORE_BORDER_CUT = 0;
2029   BARREL_ELECTRONS_AFTER_BORDER_CUT = 0;
2030
2031   return;
2032 }
2033
2034 void ZeeCalibration::resetHistograms() {
2035   h1_eventsBeforeEWKSelection_->Reset();
2036   h1_eventsAfterEWKSelection_->Reset();
2037
2038   h1_eventsBeforeBorderSelection_->Reset();
2039   h1_eventsAfterBorderSelection_->Reset();
2040
2041   for (int i = 0; i < 2; i++) {
2042     h_eleEffEta_[i]->Reset();
2043     h_eleEffPhi_[i]->Reset();
2044     h_eleEffPt_[i]->Reset();
2045   }
2046
2047   h1_seedOverSC_->Reset();
2048   h1_preshowerOverSC_->Reset();
2049
2050   h1_eleEtaResol_->Reset();
2051   h1_elePhiResol_->Reset();
2052
2053   h1_zMassResol_->Reset();
2054
2055   h1_electronCosTheta_TK_->Reset();
2056   h1_electronCosTheta_SC_->Reset();
2057   h1_electronCosTheta_SC_TK_->Reset();
2058
2059   h2_fEtaBarrelGood_->Reset();
2060   h2_fEtaBarrelBad_->Reset();
2061   h2_fEtaEndcapGood_->Reset();
2062   h2_fEtaEndcapBad_->Reset();
2063   h1_eleClasses_->Reset();
2064
2065   h1_ZCandMult_->Reset();
2066   h1_reco_ZMass_->Reset();
2067   h1_reco_ZMassCorr_->Reset();
2068   h1_reco_ZMassCorrBB_->Reset();
2069   h1_reco_ZMassCorrEE_->Reset();
2070   h1_occupancyVsEta_->Reset();
2071   h1_occupancy_->Reset();
2072   h1_occupancyBarrel_->Reset();
2073   h1_occupancyEndcap_->Reset();
2074
2075   return;
2076 }
2077
2078 void ZeeCalibration::printStatistics() {
2079   std::cout << "[ CHECK ON BARREL ELECTRON NUMBER ]"
2080             << " first " << BARREL_ELECTRONS_BEFORE_BORDER_CUT << " second " << TOTAL_ELECTRONS_IN_BARREL << std::endl;
2081
2082   std::cout << "[ EFFICIENCY OF THE BORDER SELECTION ]"
2083             << (float)BARREL_ELECTRONS_AFTER_BORDER_CUT / (float)BARREL_ELECTRONS_BEFORE_BORDER_CUT << std::endl;
2084
2085   std::cout << "[ EFFICIENCY OF THE GOLDEN SELECTION ] BARREL: "
2086             << (float)GOLDEN_ELECTRONS_IN_BARREL / (float)TOTAL_ELECTRONS_IN_BARREL
2087             << " ENDCAP: " << (float)GOLDEN_ELECTRONS_IN_ENDCAP / (float)TOTAL_ELECTRONS_IN_ENDCAP << std::endl;
2088
2089   std::cout << "[ EFFICIENCY OF THE SILVER SELECTION ] BARREL: "
2090             << (float)SILVER_ELECTRONS_IN_BARREL / (float)TOTAL_ELECTRONS_IN_BARREL
2091             << " ENDCAP: " << (float)SILVER_ELECTRONS_IN_ENDCAP / (float)TOTAL_ELECTRONS_IN_ENDCAP << std::endl;
2092
2093   std::cout << "[ EFFICIENCY OF THE SHOWER SELECTION ] BARREL: "
2094             << (float)SHOWER_ELECTRONS_IN_BARREL / (float)TOTAL_ELECTRONS_IN_BARREL
2095             << " ENDCAP: " << (float)SHOWER_ELECTRONS_IN_ENDCAP / (float)TOTAL_ELECTRONS_IN_ENDCAP << std::endl;
2096
2097   std::cout << "[ EFFICIENCY OF THE CRACK SELECTION ] BARREL: "
2098             << (float)CRACK_ELECTRONS_IN_BARREL / (float)TOTAL_ELECTRONS_IN_BARREL
2099             << " ENDCAP: " << (float)CRACK_ELECTRONS_IN_ENDCAP / (float)TOTAL_ELECTRONS_IN_ENDCAP << std::endl;
2100
2101   std::ofstream fout("ZeeStatistics.txt");
2102
2103   if (!fout) {
2104     std::cout << "Cannot open output file.\n";
2105   }
2106
2107   fout << "ZeeStatistics" << std::endl;
2108
2109   fout << "##########################RECO#########################" << std::endl;
2110   fout << "##################Zee with Barrel-Barrel electrons: " << BBZN << std::endl;
2111
2112   fout << "Golden-Golden fraction: " << (float)BBZN_gg / BBZN << " 3-3 fraction is " << (float)BBZN_tt / BBZN
2113        << " 3-whatever fraction is " << (float)BBZN_t0 / BBZN << std::endl;
2114   fout << "##################Zee with Barrel-Endcap electrons: " << EBZN << std::endl;
2115   fout << "Golden-Golden fraction: " << (float)EBZN_gg / EBZN << " 3-3 fraction is " << (float)EBZN_tt / EBZN
2116        << " 3-whatever fraction is " << (float)EBZN_t0 / EBZN << std::endl;
2117   fout << "##################Zee with Endcap-Endcap electrons: " << EEZN << std::endl;
2118   fout << "Golden-Golden fraction: " << (float)EEZN_gg / EEZN << " 3-3 fraction is " << (float)EEZN_tt / EEZN
2119        << " 3-whatever fraction is " << (float)EEZN_t0 / EEZN << std::endl;
2120
2121   fout << "\n" << std::endl;
2122
2123   fout << "##########################GEN#########################" << std::endl;
2124   fout << "##################Zee with Barrel-Barrel electrons: " << (float)MCZBB / NEVT << std::endl;
2125   fout << "##################Zee with Barrel-Endcap electrons: " << (float)MCZEB / NEVT << std::endl;
2126   fout << "##################Zee with Endcap-Endcap electrons: " << (float)MCZEE / NEVT << std::endl;
2127
2128   fout.close();
2129 }