Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-09-03 22:27:57

0001 
0002 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0003 
0004 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0005 #include "DataFormats/DetId/interface/DetId.h"
0006 #include "Calibration/Tools/interface/calibXMLwriter.h"
0007 #include "Calibration/Tools/interface/CalibrationCluster.h"
0008 #include "Calibration/Tools/interface/HouseholderDecomposition.h"
0009 #include "Calibration/Tools/interface/MinL3Algorithm.h"
0010 #include "Calibration/EcalCalibAlgos/interface/ElectronCalibrationUniv.h"
0011 #include "FWCore/Utilities/interface/isFinite.h"
0012 #include "TF1.h"
0013 #include "TRandom.h"
0014 
0015 #include <iostream>
0016 #include <string>
0017 #include <stdexcept>
0018 #include <vector>
0019 
0020 ElectronCalibrationUniv::ElectronCalibrationUniv(const edm::ParameterSet& iConfig)
0021     : ebRecHitLabel_(iConfig.getParameter<edm::InputTag>("ebRecHitsLabel")),
0022       eeRecHitLabel_(iConfig.getParameter<edm::InputTag>("eeRecHitsLabel")),
0023       electronLabel_(iConfig.getParameter<edm::InputTag>("electronLabel")),
0024       rootfile_(iConfig.getParameter<std::string>("rootfile")),
0025       calibAlgo_(iConfig.getParameter<std::string>("CALIBRATION_ALGO")),
0026       miscalibfile_(iConfig.getParameter<std::string>("miscalibfile")),
0027       miscalibfileEndCap_(iConfig.getParameter<std::string>("miscalibfileEndCap")),
0028       keventweight_(iConfig.getParameter<int>("keventweight")),
0029       elePt_(iConfig.getParameter<double>("ElePt")),
0030       maxeta_(iConfig.getParameter<int>("maxeta")),
0031       mineta_(iConfig.getParameter<int>("mineta")),
0032       maxphi_(iConfig.getParameter<int>("maxphi")),
0033       minphi_(iConfig.getParameter<int>("minphi")),
0034       cut1_(iConfig.getParameter<double>("cut1")),
0035       cut2_(iConfig.getParameter<double>("cut2")),
0036       cut3_(iConfig.getParameter<double>("cut3")),
0037       numevent_(iConfig.getParameter<int>("numevent")),
0038       cutEPCalo1_(iConfig.getParameter<double>("cutEPCaloMin")),
0039       cutEPCalo2_(iConfig.getParameter<double>("cutEPCaloMax")),
0040       cutEPin1_(iConfig.getParameter<double>("cutEPinMin")),
0041       cutEPin2_(iConfig.getParameter<double>("cutEPinMax")),
0042       cutCalo1_(iConfig.getParameter<double>("cutCaloMin")),
0043       cutCalo2_(iConfig.getParameter<double>("cutCaloMax")),
0044       cutESeed_(iConfig.getParameter<double>("cutESeed")),
0045       clusterSize_(iConfig.getParameter<int>("Clustersize")),
0046       elecclass_(iConfig.getParameter<int>("elecclass")),
0047       ebRecHitToken_(consumes<EBRecHitCollection>(ebRecHitLabel_)),
0048       eeRecHitToken_(consumes<EERecHitCollection>(eeRecHitLabel_)),
0049       gsfElectronToken_(consumes<reco::GsfElectronCollection>(electronLabel_)),
0050       topologyToken_(esConsumes<edm::Transition::BeginRun>()) {}
0051 
0052 ElectronCalibrationUniv::~ElectronCalibrationUniv() {}
0053 
0054 //========================================================================
0055 void ElectronCalibrationUniv::beginJob() {
0056   //========================================================================
0057   f = new TFile(rootfile_.c_str(), "RECREATE");
0058   f->cd();
0059   EventsAfterCuts = new TH1F("EventsAfterCuts", "Events After Cuts", 30, 0, 30);
0060 
0061   // Book histograms
0062   e9 = new TH1F("e9", "E9 energy", 300, 0., 150.);
0063   e25 = new TH1F("e25", "E25 energy", 300, 0., 150.);
0064   scE = new TH1F("scE", "SC energy", 300, 0., 150.);
0065   trP = new TH1F("trP", "Trk momentum", 300, 0., 150.);
0066   EoP = new TH1F("EoP", "EoP", 600, 0., 3.);
0067   EoP_all = new TH1F("EoP_all", "EoP_all", 600, 0., 3.);
0068 
0069   calibs = new TH1F("calib", "Calibration constants", 800, 0.5, 2.);
0070   calibsEndCapMinus = new TH1F("calibEndCapMinus", "Calibration constants EE-", 800, 0.5, 2.);
0071   calibsEndCapPlus = new TH1F("calibEndCapPlus", "Calibration constants EE+", 800, 0.5, 2.);
0072 
0073   e25OverScE = new TH1F("e25OverscE", "E25 / SC energy", 400, 0., 2.);
0074   E25oP = new TH1F("E25oP", "E25 / P", 750, 0., 1.5);
0075 
0076   Map = new TH2F("Map", "Nb Events in Crystal", 173, -86, 86, 362, 0, 361);
0077   e9Overe25 = new TH1F("e9Overe25", "E9 / E25", 400, 0., 2.);
0078   Map3Dcalib = new TH2F("3Dcalib", "3Dcalib", 173, -86, 86, 362, 0, 361);
0079   Map3DcalibEndCapMinus = new TH2F("3DcalibEndCapMinus", "3Dcalib EE-", 100, 0, 100, 100, 0, 100);
0080   Map3DcalibEndCapPlus = new TH2F("3DcalibEndCapPlus", "3Dcalib EE+", 100, 0, 100, 100, 0, 100);
0081 
0082   MapCor1 = new TH2F("MapCor1", "Correlation E25/Pcalo versus E25/Pin", 100, 0., 5., 100, 0., 5.);
0083   MapCor2 = new TH2F("MapCor2", "Correlation E25/Pcalo versus E/P", 100, 0., 5., 100, 0., 5.);
0084   MapCor3 = new TH2F("MapCor3", "Correlation E25/Pcalo versus Pout/Pin", 100, 0., 5., 100, 0., 5.);
0085   MapCor4 = new TH2F("MapCor4", "Correlation E25/Pcalo versus E25/highestP", 100, 0., 5., 100, 0., 5.);
0086   MapCor5 = new TH2F("MapCor5", "Correlation E25/Pcalo versus Pcalo/Pout", 100, 0., 5., 100, 0., 5.);
0087   MapCor6 = new TH2F("MapCor6", "Correlation Pout/Pin versus E25/Pin", 100, 0., 5., 100, 0., 5.);
0088   MapCor7 = new TH2F("MapCor7", "Correlation Pout/Pin versus Pcalo/Pout", 100, 0., 5., 100, 0., 5.);
0089   MapCor8 = new TH2F("MapCor8", "Correlation E25/Pin versus Pcalo/Pout", 100, 0., 5., 100, 0., 5.);
0090   MapCor9 = new TH2F("MapCor9", "Correlation  E25/Pcalo versus Eseed/Pout", 100, 0., 5., 100, 0., 5.);
0091   MapCor10 = new TH2F("MapCor10", "Correlation Eseed/Pout versus Pout/Pin", 100, 0., 5., 100, 0., 5.);
0092   MapCor11 = new TH2F("MapCor11", "Correlation Eseed/Pout versus E25/Pin", 100, 0., 5., 100, 0., 5.);
0093   //   MapCorCalib = new TH2F ("MapCorCalib", "Correlation Miscalibration versus Calibration constants", 500, 0.5,1.5, 500, 0.5, 1.5);
0094 
0095   E25oPvsEta = new TH2F("E25oPvsEta", "E/P vs Eta", 173, -86, 86, 600, 0.7, 1.3);
0096   E25oPvsEtaEndCapMinus = new TH2F("E25oPvsEtaEndCapMinus", "E/P vs R EE-", 100, 0, 100, 600, 0.7, 1.3);
0097   E25oPvsEtaEndCapPlus = new TH2F("E25oPvsEtaEndCapPlus", "E/P vs R EE+", 100, 0, 100, 600, 0.7, 1.3);
0098 
0099   PinMinPout = new TH1F("PinMinPout", "(Pin - Pout)/Pin", 600, -2.0, 2.0);
0100 
0101   calibinter = new TH1F("calibinter", "internal calibration constants", 800, 0.5, 2.);
0102   PinOverPout = new TH1F("PinOverPout", "pinOverpout", 600, 0., 3.);
0103   eSeedOverPout = new TH1F("eSeedOverPout", "eSeedOverpout ", 600, 0., 3.);
0104   //   MisCalibs = new TH1F("MisCalibs","Miscalibration constants",800,0.5,2.);
0105   //   RatioCalibs = new TH1F("RatioCalibs","Ratio in Calibration Constants", 800, 0.5, 2.0);
0106   //   DiffCalibs = new TH1F("DiffCalibs", "Difference in Calibration constants", 800, -1.0,1.0);
0107   calibinterEndCapMinus = new TH1F("calibinterEndCapMinus", "internal calibration constants", 800, 0.5, 2.);
0108   calibinterEndCapPlus = new TH1F("calibinterEndCapPlus", "internal calibration constants", 800, 0.5, 2.);
0109   //   MisCalibsEndCapMinus = new TH1F("MisCalibsEndCapMinus","Miscalibration constants",800,0.5,2.);
0110   //   MisCalibsEndCapPlus = new TH1F("MisCalibsEndCapPlus","Miscalibration constants",800,0.5,2.);
0111   //   RatioCalibsEndCapMinus = new TH1F("RatioCalibsEndCapMinus","Ratio in Calibration Constants", 800, 0.5, 2.0);
0112   //   RatioCalibsEndCapPlus = new TH1F("RatioCalibsEndCapPlus","Ratio in Calibration Constants", 800, 0.5, 2.0);
0113   //   DiffCalibsEndCapMinus = new TH1F("DiffCalibsEndCapMinus", "Difference in Calibration constants", 800, -1.0,1.0);
0114   //   DiffCalibsEndCapPlus = new TH1F("DiffCalibsEndCapPlus", "Difference in Calibration constants", 800, -1.0,1.0);
0115   Error1 = new TH1F("Error1", "DeltaP/Pin", 800, -1.0, 1.0);
0116   Error2 = new TH1F("Error2", "DeltaP/Pout", 800, -1.0, 1.0);
0117   Error3 = new TH1F("Error3", "DeltaP/Pcalo", 800, -1.0, 1.0);
0118   eSeedOverPout2 = new TH1F("eSeedOverPout2", "eSeedOverpout (No Supercluster)", 600, 0., 4.);
0119   hadOverEm = new TH1F("hadOverEm", "Had/EM distribution", 600, -2., 2.);
0120 
0121   // Book histograms
0122   Map3DcalibNoCuts = new TH2F("3DcalibNoCuts", "3Dcalib (Before Cuts)", 173, -86, 86, 362, 0, 361);
0123   e9NoCuts = new TH1F("e9NoCuts", "E9 energy (Before Cuts)", 300, 0., 150.);
0124   e25NoCuts = new TH1F("e25NoCuts", "E25 energy (Before Cuts)", 300, 0., 150.);
0125   scENoCuts = new TH1F("scENoCuts", "SC energy (Before Cuts)", 300, 0., 150.);
0126   trPNoCuts = new TH1F("trPNoCuts", "Trk momentum (Before Cuts)", 300, 0., 150.);
0127   EoPNoCuts = new TH1F("EoPNoCuts", "EoP (Before Cuts)", 600, 0., 3.);
0128   calibsNoCuts = new TH1F("calibNoCuts", "Calibration constants (Before Cuts)", 800, 0., 2.);
0129   e25OverScENoCuts = new TH1F("e25OverscENoCuts", "E25 / SC energy (Before Cuts)", 400, 0., 2.);
0130   E25oPNoCuts = new TH1F("E25oPNoCuts", "E25 / P (Before Cuts)", 750, 0., 1.5);
0131   MapEndCapMinus = new TH2F("MapEndCapMinus", "Nb Events in Crystal (EndCap)", 100, 0, 100, 100, 0, 100);
0132   MapEndCapPlus = new TH2F("MapEndCapPlus", "Nb Events in Crystal (EndCap)", 100, 0, 100, 100, 0, 100);
0133   e9Overe25NoCuts = new TH1F("e9Overe25NoCuts", "E9 / E25 (Before Cuts)", 400, 0., 2.);
0134   PinOverPoutNoCuts = new TH1F("PinOverPoutNoCuts", "pinOverpout (Before Cuts)", 600, 0., 3.);
0135   eSeedOverPoutNoCuts = new TH1F(" eSeedOverPoutNoCuts", "eSeedOverpout (Before Cuts) ", 600, 0., 4.);
0136   PinMinPoutNoCuts = new TH1F("PinMinPoutNoCuts", "(Pin - Pout)/Pin (Before Cuts)", 600, -2.0, 2.0);
0137 
0138   //   RatioCalibsNoCuts = new TH1F("RatioCalibsNoCuts","Ratio in Calibration Constants (Before Cuts)", 800, 0.5, 2.0);
0139   //   DiffCalibsNoCuts = new TH1F("DiffCalibsNoCuts", "Difference in Calibration constants (Before Cuts)", 800, -1.0,1.0);
0140   calibinterNoCuts = new TH1F("calibinterNoCuts", "internal calibration constants", 2000, 0.5, 2.);
0141 
0142   MapCor1NoCuts =
0143       new TH2F("MapCor1NoCuts", "Correlation E25/PatCalo versus E25/Pin (Before Cuts)", 100, 0., 5., 100, 0., 5.);
0144   MapCor2NoCuts =
0145       new TH2F("MapCor2NoCuts", "Correlation E25/PatCalo versus E/P (Before Cuts)", 100, 0., 5., 100, 0., 5.);
0146   MapCor3NoCuts =
0147       new TH2F("MapCor3NoCuts", "Correlation E25/PatCalo versus Pout/Pin (Before Cuts)", 100, 0., 5., 100, 0., 5.);
0148   MapCor4NoCuts =
0149       new TH2F("MapCor4NoCuts", "Correlation E25/PatCalo versus E25/highestP (Before Cuts)", 100, 0., 5., 100, 0., 5.);
0150   MapCor5NoCuts =
0151       new TH2F("MapCor5NoCuts", "Correlation E25/Pcalo versus Pcalo/Pout (Before Cuts)", 100, 0., 5., 100, 0., 5.);
0152   MapCor6NoCuts =
0153       new TH2F("MapCor6NoCuts", "Correlation Pout/Pin versus E25/Pin (Before Cuts)", 100, 0., 5., 100, 0., 5.);
0154   MapCor7NoCuts =
0155       new TH2F("MapCor7NoCuts", "Correlation Pout/Pin versus Pcalo/Pout (Before Cuts)", 100, 0., 5., 100, 0., 5.);
0156   MapCor8NoCuts =
0157       new TH2F("MapCor8NoCuts", "Correlation E25/Pin versus Pcalo/Pout (Before Cuts)", 100, 0., 5., 100, 0., 5.);
0158   MapCor9NoCuts =
0159       new TH2F("MapCor9NoCuts", "Correlation  E25/Pcalo versus Eseed/Pout (Before Cuts)", 100, 0., 5., 100, 0., 5.);
0160   MapCor10NoCuts =
0161       new TH2F("MapCor10NoCuts", "Correlation Eseed/Pout versus Pout/Pin (Before Cuts)", 100, 0., 5., 100, 0., 5.);
0162   MapCor11NoCuts =
0163       new TH2F("MapCor11NoCuts", "Correlation Eseed/Pout versus E25/Pin (Before Cuts)", 100, 0., 5., 100, 0., 5.);
0164   //   MapCorCalibEndCapMinus = new TH2F ("MapCorCalibEndCapMinus", "Correlation Miscalibration versus Calibration constants (EndCap)",  500, 0.5,1.5, 500, 0.5, 1.5);
0165   //   MapCorCalibEndCapPlus = new TH2F ("MapCorCalibEndCapPlus", "Correlation Miscalibration versus Calibration constants (EndCap)",  500, 0.5,1.5, 500, 0.5, 1.5);
0166 
0167   Error1NoCuts = new TH1F("Eror1NoCuts", "DeltaP/Pin (Before Cuts)", 800, -1.0, 1.0);
0168   Error2NoCuts = new TH1F("Error2NoCuts", "DeltaP/Pout (Before Cuts)", 800, -1.0, 1.0);
0169   Error3NoCuts = new TH1F("Error3NoCuts", "DeltaP/Pcalo (Before Cuts)", 800, -1.0, 1.0);
0170   eSeedOverPout2NoCuts = new TH1F("eSeedOverPout2NoCuts", "eSeedOverpout (No Supercluster, Before Cuts)", 600, 0., 4.);
0171   hadOverEmNoCuts = new TH1F("hadOverEmNoCuts", "Had/EM distribution (Before Cuts)", 600, -2., 2.);
0172 
0173   //Book histograms after ESeed cut
0174   MapCor1ESeed =
0175       new TH2F("MapCor1ESeed", "Correlation E25/Pcalo versus E25/Pin (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
0176   MapCor2ESeed =
0177       new TH2F("MapCor2ESeed", "Correlation E25/Pcalo versus E/P (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
0178   MapCor3ESeed = new TH2F(
0179       "MapCor3ESeed", "Correlation E25/Pcalo versus Pout/Pin (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
0180   MapCor4ESeed = new TH2F(
0181       "MapCor4ESeed", "Correlation E25/Pcalo versus E25/highestP (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
0182   MapCor5ESeed = new TH2F(
0183       "MapCor5ESeed", "Correlation E25/Pcalo versus Pcalo/Pout (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
0184   MapCor6ESeed =
0185       new TH2F("MapCor6ESeed", "Correlation Pout/Pin versus E25/Pin (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
0186   MapCor7ESeed = new TH2F(
0187       "MapCor7ESeed", "Correlation Pout/Pin versus Pcalo/Pout (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
0188   MapCor8ESeed = new TH2F(
0189       "MapCor8ESeed", "Correlation E25/Pin versus Pcalo/Pout (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
0190   MapCor9ESeed = new TH2F(
0191       "MapCor9ESeed", "Correlation  E25/Pcalo versus Eseed/Pout (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
0192   MapCor10ESeed = new TH2F(
0193       "MapCor10ESeed", "Correlation Eseed/Pout versus Pout/Pin (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
0194   MapCor11ESeed = new TH2F(
0195       "MapCor11ESeed", "Correlation Eseed/Pout versus E25/Pin (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
0196 
0197   eSeedOverPout2ESeed =
0198       new TH1F("eSeedOverPout2ESeed", "eSeedOverpout (No Supercluster, after Eseed/Pout cut)", 600, 0., 4.);
0199 
0200   hadOverEmESeed = new TH1F("hadOverEmESeed", "Had/EM distribution (after Eseed/Pout cut)", 600, -2., 2.);
0201 
0202   //Book histograms without any cut
0203   GeneralMap = new TH2F("GeneralMap", "Map without any cuts", 173, -86, 86, 362, 0, 361);
0204   GeneralMapEndCapMinus = new TH2F("GeneralMapEndCapMinus", "Map without any cuts", 100, 0, 100, 100, 0, 100);
0205   GeneralMapEndCapPlus = new TH2F("GeneralMapEndCapPlus", "Map without any cuts", 100, 0, 100, 100, 0, 100);
0206   GeneralMapBeforePt = new TH2F("GeneralMapBeforePt", "Map without any cuts", 173, -86, 86, 362, 0, 361);
0207   GeneralMapEndCapMinusBeforePt =
0208       new TH2F("GeneralMapEndCapMinusBeforePt", "Map without any cuts", 100, 0, 100, 100, 0, 100);
0209   GeneralMapEndCapPlusBeforePt =
0210       new TH2F("GeneralMapEndCapPlusBeforePt", "Map without any cuts", 100, 0, 100, 100, 0, 100);
0211 
0212   calibClusterSize = clusterSize_;
0213   etaMin = int(mineta_);
0214   etaMax = int(maxeta_);
0215   phiMin = int(minphi_);
0216   phiMax = int(maxphi_);
0217   if (calibAlgo_ == "L3") {
0218     MyL3Algo1 = new MinL3Algorithm(keventweight_, calibClusterSize, etaMin, etaMax, phiMin, phiMax);
0219   } else {
0220     if (calibAlgo_ == "L3Univ") {
0221       UnivL3 = new MinL3AlgoUniv<DetId>(keventweight_);
0222     } else {
0223       if (calibAlgo_ == "HH" || calibAlgo_ == "HHReg") {
0224         MyHH = new HouseholderDecomposition(calibClusterSize, etaMin, etaMax, phiMin, phiMax);
0225       } else {
0226         std::cout << " Name of Algorithm is not recognize " << calibAlgo_
0227                   << " Should be either L3, HH or HHReg. Abort! " << std::endl;
0228       }
0229     }
0230   }
0231   read_events = 0;
0232 }
0233 
0234 //========================================================================
0235 void ElectronCalibrationUniv::beginRun(edm::Run const&, edm::EventSetup const& iSetup) {
0236   //========================================================================
0237 
0238   //To Deal with Geometry:
0239   theCaloTopology_ = &iSetup.getData(topologyToken_);
0240 }
0241 
0242 //========================================================================
0243 void ElectronCalibrationUniv::endRun(edm::Run const&, edm::EventSetup const& iSetup) {}
0244 //========================================================================
0245 
0246 //========================================================================
0247 
0248 void ElectronCalibrationUniv::endJob() {
0249   //========================================================================
0250 
0251   f->cd();
0252   time_t start, end;
0253   time_t cpu_time_used;
0254   start = time(nullptr);
0255 
0256   //In order to do only one loop to use properly looper properties, ask only for 1 iterations!
0257   int nIterations = 10;
0258   if (calibAlgo_ == "L3") {
0259     solution = MyL3Algo1->iterate(EventMatrix, MaxCCeta, MaxCCphi, EnergyVector, nIterations);
0260   } else {
0261     if (calibAlgo_ == "L3Univ") {
0262       //Univsolution= UnivL3->getSolution();
0263       //     std::cout<<" Should derive solution "<<EnergyVector.size()<<std::endl;
0264       Univsolution = UnivL3->iterate(EventMatrix, UnivEventIds, EnergyVector, nIterations);
0265       //std::cout<<" solution size "<<Univsolution.size()<<std::endl;
0266     } else {
0267       if (calibAlgo_ == "HH") {
0268         solution = MyHH->iterate(EventMatrix, MaxCCeta, MaxCCphi, EnergyVector, 1, false);
0269       } else {
0270         if (calibAlgo_ == "HHReg") {
0271           solution = MyHH->runRegional(EventMatrix, MaxCCeta, MaxCCphi, EnergyVector, 2);
0272         } else {
0273           std::cout << " Calibration not run due to problem in Algo Choice..." << std::endl;
0274           return;
0275         }
0276       }
0277     }
0278   }
0279   end = time(nullptr);
0280   cpu_time_used = end - start;
0281   //     std::cout<<"222 solution size "<<Univsolution.size()<<std::endl;
0282 
0283   calibXMLwriter write_calibrations;
0284 
0285   //   FILE* MisCalib;
0286   //   //char* calibfile="miscalibfile";
0287   //   MisCalib = fopen(miscalibfile_.c_str(),"r");
0288 
0289   //   int fileStatus=0;
0290   //   int eta=-1;
0291   //   int phi=-1;
0292   //   float coeff=-1;
0293 
0294   std::map<EBDetId, float> OldCoeff;
0295 
0296   //  while(fileStatus != EOF) {
0297   //    fileStatus = fscanf(MisCalib,"%d %d %f\n",  &eta,&phi,&coeff);
0298   //    if(eta!=-1&&phi!=-1&& coeff!=-1){
0299   //      //     std::cout<<" We have read correctly the coefficient " << coeff << " corresponding to eta "<<eta<<" and  phi "<<phi<<std::endl;
0300   //      OldCoeff.insert(std::make_pair(EBDetId(eta,phi,EBDetId::ETAPHIMODE),coeff ));
0301   //    }
0302   //  }
0303 
0304   //  fclose(MisCalib);
0305   //   FILE* MisCalibEndCap;
0306   //   //char* calibfile="miscalibfile";
0307   //   MisCalibEndCap = fopen(miscalibfileEndCap_.c_str(),"r");
0308 
0309   //   int fileStatus2=0;
0310   //   int X=-1;
0311   //   int Y=-1;
0312   //   float coeff2=-1;
0313   std::map<EEDetId, float> OldCoeffEndCap;
0314 
0315   //  while(fileStatus2 != EOF) {
0316   //    fileStatus2 = fscanf(MisCalibEndCap,"%d %d %f\n",  &X,&Y,&coeff2);
0317   //    if(X!=-1&&Y!=-1&& coeff2!=-1){
0318   //      //     std::cout<<" We have read correctly the coefficient " << coeff << " corresponding to eta "<<eta<<" and  phi "<<phi<<std::endl;
0319   //      if(TestEEvalidDetId(X,Y,1)){
0320   //        OldCoeffEndCap.insert(std::make_pair(EEDetId(X,Y,1,EEDetId::XYMODE),coeff2 ));
0321   //      }
0322   //    }
0323   //  }
0324 
0325   // fclose(MisCalibEndCap);
0326   std::map<DetId, float>::const_iterator itmap;
0327   for (itmap = Univsolution.begin(); itmap != Univsolution.end(); itmap++) {
0328     const DetId Id(itmap->first);
0329     if (Id.subdetId() == 1) {
0330       const EBDetId IChannelDetId(itmap->first);
0331       if (IChannelDetId.ieta() < mineta_) {
0332         continue;
0333       }
0334       if (IChannelDetId.ieta() > maxeta_) {
0335         continue;
0336       }
0337       if (IChannelDetId.iphi() < minphi_) {
0338         continue;
0339       }
0340       if (IChannelDetId.iphi() > maxphi_) {
0341         continue;
0342       }
0343       //      float Compare=1;
0344       //      std::map<EBDetId,float>::iterator iter = OldCoeff.find(itmap->first);
0345       //      if( iter != OldCoeff.end() )Compare = iter->second;
0346       Map3Dcalib->Fill(IChannelDetId.ieta(), IChannelDetId.iphi(), itmap->second);
0347       calibs->Fill(itmap->second);
0348       //DiffCalibs->Fill(newCalibs[icry]-miscalib[IChannelDetId.ieta()-1][IChannelDetId.iphi()-21]);
0349       //RatioCalibs->Fill(newCalibs[icry]/miscalib[IChannelDetId.ieta()-1][IChannelDetId.iphi()-21]);
0350       if (IChannelDetId.ieta() < mineta_ + 2) {
0351         continue;
0352       }
0353       if (IChannelDetId.ieta() > maxeta_ - 2) {
0354         continue;
0355       }
0356       if (IChannelDetId.iphi() < minphi_ + 2) {
0357         continue;
0358       }
0359       if (IChannelDetId.iphi() > maxphi_ - 2) {
0360         continue;
0361       }
0362       write_calibrations.writeLine(IChannelDetId, itmap->second);
0363       calibinter->Fill(itmap->second);
0364       //       MapCorCalib->Fill(itmap->second,Compare);
0365       //       DiffCalibs->Fill(itmap->second-Compare);
0366       //       RatioCalibs->Fill(itmap->second*Compare);
0367     } else {
0368       const EEDetId IChannelDetId(itmap->first);
0369       //       if (IChannelDetId.ix()<0 ){continue;}
0370       //       if (IChannelDetId.ix()>100 ){continue;}
0371       //       if (IChannelDetId.iy()<0 ){continue;}
0372       //       if (IChannelDetId.iy()>100 ){continue;}
0373       //     std::map<EEDetId,float>::iterator iter = OldCoeffEndCap.find(itmap->first);
0374       //      float Compare=1;
0375       //      if( iter != OldCoeffEndCap.end() )Compare = iter->second;
0376       if (IChannelDetId.zside() < 0) {
0377         Map3DcalibEndCapMinus->Fill(IChannelDetId.ix(), IChannelDetId.iy(), itmap->second);
0378         calibsEndCapMinus->Fill(itmap->second);
0379         calibinterEndCapMinus->Fill(itmap->second);
0380         //  DiffCalibsEndCapMinus->Fill(itmap->second-Compare);
0381         //  RatioCalibsEndCapMinus->Fill(itmap->second*Compare);
0382         //  MapCorCalibEndCapMinus->Fill(itmap->second,Compare);
0383       } else {
0384         Map3DcalibEndCapPlus->Fill(IChannelDetId.ix(), IChannelDetId.iy(), itmap->second);
0385         calibsEndCapPlus->Fill(itmap->second);
0386         calibinterEndCapPlus->Fill(itmap->second);
0387         //  DiffCalibsEndCapPlus->Fill(itmap->second-Compare);
0388         //  RatioCalibsEndCapPlus->Fill(itmap->second*Compare);
0389         //  MapCorCalibEndCapPlus->Fill(itmap->second,Compare);
0390       }
0391       write_calibrations.writeLine(IChannelDetId, itmap->second);
0392     }
0393   }
0394   EventsAfterCuts->Write();
0395 
0396   // Book histograms
0397   e25->Write();
0398   e9->Write();
0399   scE->Write();
0400   trP->Write();
0401   EoP->Write();
0402   EoP_all->Write();
0403   calibs->Write();
0404   calibsEndCapMinus->Write();
0405   calibsEndCapPlus->Write();
0406   e9Overe25->Write();
0407   e25OverScE->Write();
0408   Map->Write();
0409   E25oP->Write();
0410 
0411   PinOverPout->Write();
0412   eSeedOverPout->Write();
0413   //   MisCalibs->Write();
0414   //   RatioCalibs->Write();
0415   //   DiffCalibs->Write();
0416   //   RatioCalibsNoCuts->Write();
0417   //   DiffCalibsNoCuts->Write();
0418   //   MisCalibsEndCapMinus->Write();
0419   //   MisCalibsEndCapPlus->Write();
0420   //   RatioCalibsEndCapMinus->Write();
0421   //   RatioCalibsEndCapPlus->Write();
0422   //   DiffCalibsEndCapMinus->Write();
0423   //   DiffCalibsEndCapPlus->Write();
0424 
0425   e25NoCuts->Write();
0426   e9NoCuts->Write();
0427   scENoCuts->Write();
0428   trPNoCuts->Write();
0429   EoPNoCuts->Write();
0430   calibsNoCuts->Write();
0431   e9Overe25NoCuts->Write();
0432   e25OverScENoCuts->Write();
0433   MapEndCapMinus->Write();
0434   MapEndCapPlus->Write();
0435   E25oPNoCuts->Write();
0436   Map3Dcalib->Write();
0437   Map3DcalibEndCapMinus->Write();
0438   Map3DcalibEndCapPlus->Write();
0439   Map3DcalibNoCuts->Write();
0440   calibinter->Write();
0441   calibinterEndCapMinus->Write();
0442   calibinterEndCapPlus->Write();
0443   calibinterNoCuts->Write();
0444   PinOverPoutNoCuts->Write();
0445   eSeedOverPoutNoCuts->Write();
0446 
0447   GeneralMap->Write();
0448   GeneralMapEndCapMinus->Write();
0449   GeneralMapEndCapPlus->Write();
0450   GeneralMapBeforePt->Write();
0451   GeneralMapEndCapMinusBeforePt->Write();
0452   GeneralMapEndCapPlusBeforePt->Write();
0453 
0454   MapCor1->Write();
0455   MapCor2->Write();
0456   MapCor3->Write();
0457   MapCor4->Write();
0458   MapCor5->Write();
0459   MapCor6->Write();
0460   MapCor7->Write();
0461   MapCor8->Write();
0462   MapCor9->Write();
0463   MapCor10->Write();
0464   MapCor11->Write();
0465   //  MapCorCalib->Write();
0466 
0467   MapCor1NoCuts->Write();
0468   MapCor2NoCuts->Write();
0469   MapCor3NoCuts->Write();
0470   MapCor4NoCuts->Write();
0471   MapCor5NoCuts->Write();
0472   MapCor6NoCuts->Write();
0473   MapCor7NoCuts->Write();
0474   MapCor8NoCuts->Write();
0475   MapCor9NoCuts->Write();
0476   MapCor10NoCuts->Write();
0477   MapCor11NoCuts->Write();
0478   //   MapCorCalibEndCapMinus->Write();
0479   //   MapCorCalibEndCapPlus->Write();
0480 
0481   MapCor1ESeed->Write();
0482   MapCor2ESeed->Write();
0483   MapCor3ESeed->Write();
0484   MapCor4ESeed->Write();
0485   MapCor5ESeed->Write();
0486   MapCor6ESeed->Write();
0487   MapCor7ESeed->Write();
0488   MapCor8ESeed->Write();
0489   MapCor9ESeed->Write();
0490   MapCor10ESeed->Write();
0491   MapCor11ESeed->Write();
0492 
0493   E25oPvsEta->Write();
0494   E25oPvsEtaEndCapMinus->Write();
0495   E25oPvsEtaEndCapPlus->Write();
0496 
0497   PinMinPout->Write();
0498   PinMinPoutNoCuts->Write();
0499 
0500   Error1->Write();
0501   Error2->Write();
0502   Error3->Write();
0503   Error1NoCuts->Write();
0504   Error2NoCuts->Write();
0505   Error3NoCuts->Write();
0506 
0507   eSeedOverPout2->Write();
0508   eSeedOverPout2NoCuts->Write();
0509   eSeedOverPout2ESeed->Write();
0510 
0511   hadOverEm->Write();
0512   hadOverEmNoCuts->Write();
0513   hadOverEmESeed->Write();
0514 
0515   f->Write();
0516 
0517   f->Close();
0518   //   if(MyL3Algo1)delete MyL3Algo1;
0519   //   if(UnivL3)delete UnivL3;
0520   //   if(MyHH)delete MyHH;
0521   //  delete f;
0522   ////////////////////////       FINAL STATISTICS           ////////////////////
0523 
0524   std::cout << " " << std::endl;
0525   std::cout << "************* STATISTICS **************" << std::endl;
0526   std::cout << " Events Studied " << read_events << std::endl;
0527   std::cout << "Timing info:" << std::endl;
0528   std::cout << "CPU time usage  -- calibrating: " << cpu_time_used << " sec." << std::endl;
0529 
0530   /////////////////////////////////////////////////////////////////////////////
0531 }
0532 
0533 DetId ElectronCalibrationUniv::findMaxHit(const std::vector<DetId>& v1,
0534                                           const EBRecHitCollection* EBhits,
0535                                           const EERecHitCollection* EEhits) {
0536   //=================================================================================
0537 
0538   double currEnergy = 0.;
0539   DetId maxHit;
0540 
0541   for (std::vector<DetId>::const_iterator idsIt = v1.begin(); idsIt != v1.end(); ++idsIt) {
0542     if (idsIt->subdetId() == 1) {
0543       EBRecHitCollection::const_iterator itrechit;
0544       itrechit = EBhits->find(*idsIt);
0545       if (itrechit == EBhits->end()) {
0546         std::cout << "ElectronCalibration::findMaxHit: rechit not found! " << (EBDetId)(*idsIt) << std::endl;
0547         continue;
0548       }
0549       if (itrechit->energy() > currEnergy) {
0550         currEnergy = itrechit->energy();
0551         maxHit = *idsIt;
0552       }
0553     } else {
0554       EERecHitCollection::const_iterator itrechit;
0555       itrechit = EEhits->find(*idsIt);
0556       if (itrechit == EEhits->end()) {
0557         std::cout << "ElectronCalibration::findMaxHit: rechit not found! idsIt = " << (EEDetId)(*idsIt) << std::endl;
0558         continue;
0559       }
0560 
0561       if (itrechit->energy() > currEnergy) {
0562         currEnergy = itrechit->energy();
0563         maxHit = *idsIt;
0564       }
0565     }
0566   }
0567 
0568   return maxHit;
0569 }
0570 
0571 bool ElectronCalibrationUniv::TestEEvalidDetId(int crystal_ix, int crystal_iy, int iz) {
0572   bool valid = false;
0573   if (crystal_ix < 1 || crystal_ix > 100 || crystal_iy < 1 || crystal_iy > 100 || abs(iz) != 1) {
0574     return valid;
0575   }
0576   if ((crystal_ix >= 1 && crystal_ix <= 3 && (crystal_iy <= 40 || crystal_iy > 60)) ||
0577       (crystal_ix >= 4 && crystal_ix <= 5 && (crystal_iy <= 35 || crystal_iy > 65)) ||
0578       (crystal_ix >= 6 && crystal_ix <= 8 && (crystal_iy <= 25 || crystal_iy > 75)) ||
0579       (crystal_ix >= 9 && crystal_ix <= 13 && (crystal_iy <= 20 || crystal_iy > 80)) ||
0580       (crystal_ix >= 14 && crystal_ix <= 15 && (crystal_iy <= 15 || crystal_iy > 85)) ||
0581       (crystal_ix >= 16 && crystal_ix <= 20 && (crystal_iy <= 13 || crystal_iy > 87)) ||
0582       (crystal_ix >= 21 && crystal_ix <= 25 && (crystal_iy <= 8 || crystal_iy > 92)) ||
0583       (crystal_ix >= 26 && crystal_ix <= 35 && (crystal_iy <= 5 || crystal_iy > 95)) ||
0584       (crystal_ix >= 36 && crystal_ix <= 39 && (crystal_iy <= 3 || crystal_iy > 97)) ||
0585       (crystal_ix >= 98 && crystal_ix <= 100 && (crystal_iy <= 40 || crystal_iy > 60)) ||
0586       (crystal_ix >= 96 && crystal_ix <= 97 && (crystal_iy <= 35 || crystal_iy > 65)) ||
0587       (crystal_ix >= 93 && crystal_ix <= 95 && (crystal_iy <= 25 || crystal_iy > 75)) ||
0588       (crystal_ix >= 88 && crystal_ix <= 92 && (crystal_iy <= 20 || crystal_iy > 80)) ||
0589       (crystal_ix >= 86 && crystal_ix <= 87 && (crystal_iy <= 15 || crystal_iy > 85)) ||
0590       (crystal_ix >= 81 && crystal_ix <= 85 && (crystal_iy <= 13 || crystal_iy > 87)) ||
0591       (crystal_ix >= 76 && crystal_ix <= 80 && (crystal_iy <= 8 || crystal_iy > 92)) ||
0592       (crystal_ix >= 66 && crystal_ix <= 75 && (crystal_iy <= 5 || crystal_iy > 95)) ||
0593       (crystal_ix >= 62 && crystal_ix <= 65 && (crystal_iy <= 3 || crystal_iy > 97)) ||
0594       ((crystal_ix == 40 || crystal_ix == 61) &&
0595        ((crystal_iy >= 46 && crystal_iy <= 55) || crystal_iy <= 3 || crystal_iy > 97)) ||
0596       ((crystal_ix == 41 || crystal_ix == 60) && crystal_iy >= 44 && crystal_iy <= 57) ||
0597       ((crystal_ix == 42 || crystal_ix == 59) && crystal_iy >= 43 && crystal_iy <= 58) ||
0598       ((crystal_ix == 43 || crystal_ix == 58) && crystal_iy >= 42 && crystal_iy <= 59) ||
0599       ((crystal_ix == 44 || crystal_ix == 45 || crystal_ix == 57 || crystal_ix == 56) && crystal_iy >= 41 &&
0600        crystal_iy <= 60) ||
0601       (crystal_ix >= 46 && crystal_ix <= 55 && crystal_iy >= 40 && crystal_iy <= 61)) {
0602     return valid;
0603   }
0604   valid = true;
0605   return valid;
0606 }
0607 
0608 //=================================================================================
0609 void ElectronCalibrationUniv::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0610   //=================================================================================
0611   using namespace edm;
0612 
0613   // Get EBRecHits
0614   edm::Handle<EBRecHitCollection> EBphits;
0615   iEvent.getByToken(ebRecHitToken_, EBphits);
0616   if (!EBphits.isValid()) {
0617     std::cerr << "Error! can't get the product EBRecHitCollection: " << std::endl;
0618   }
0619   const EBRecHitCollection* EBhits = EBphits.product();  // get a ptr to the product
0620 
0621   // Get EERecHits
0622   edm::Handle<EERecHitCollection> EEphits;
0623 
0624   iEvent.getByToken(eeRecHitToken_, EEphits);
0625   if (!EEphits.isValid()) {
0626     std::cerr << "Error! can't get the product EERecHitCollection: " << std::endl;
0627   }
0628   const EERecHitCollection* EEhits = EEphits.product();  // get a ptr to the product
0629 
0630   // Get pixelElectrons
0631   edm::Handle<reco::GsfElectronCollection> pElectrons;
0632   iEvent.getByToken(gsfElectronToken_, pElectrons);
0633   if (!pElectrons.isValid()) {
0634     std::cerr << "Error! can't get the product ElectronCollection: " << std::endl;
0635   }
0636   const reco::GsfElectronCollection* electronCollection = pElectrons.product();
0637   read_events++;
0638   if (read_events % 1000 == 0)
0639     std::cout << "read_events = " << read_events << std::endl;
0640 
0641   EventsAfterCuts->Fill(1);
0642   if (!EBhits || !EEhits)
0643     return;
0644   EventsAfterCuts->Fill(2);
0645   if (EBhits->empty() && EEhits->empty())
0646     return;
0647   EventsAfterCuts->Fill(3);
0648   if (!electronCollection)
0649     return;
0650   EventsAfterCuts->Fill(4);
0651   if (electronCollection->empty())
0652     return;
0653 
0654   //    ////////////////Need to recalibrate the events (copy code from EcalRecHitRecalib):
0655 
0656   ////////////////////////////////////////////////////////////////////////////////////////
0657   ///                          START HERE....
0658   ///////////////////////////////////////////////////////////////////////////////////////
0659 
0660   reco::GsfElectronCollection::const_iterator eleIt = electronCollection->begin();
0661 
0662   reco::GsfElectron highPtElectron;
0663 
0664   float highestElePt = 0.;
0665   bool found = false;
0666   for (eleIt = electronCollection->begin(); eleIt != electronCollection->end(); eleIt++) {
0667     if (fabs(eleIt->eta()) > 2.4)
0668       continue;
0669     //     if(eleIt->eta()<0.0) continue;
0670 
0671     if (eleIt->pt() > highestElePt) {
0672       highestElePt = eleIt->pt();
0673       highPtElectron = *eleIt;
0674       found = true;
0675       //       std::cout<<" eleIt->pt( "<<eleIt->pt()<<" eleIt->eta() "<<eleIt->eta()<<std::endl;
0676     }
0677   }
0678   EventsAfterCuts->Fill(5);
0679   if (!found)
0680     return;
0681 
0682   const reco::SuperCluster& sc = *(highPtElectron.superCluster());
0683   //  if(fabs(sc.eta())>1.479){std::cout<<" SC not in Barrel "<<sc.eta()<<std::endl;;}
0684   //  const std::vector<DetId> & v1 = sc.getHitsByDetId();
0685 
0686   std::vector<DetId> v1;
0687   //Loop to fill the vector of DetIds
0688   for (std::vector<std::pair<DetId, float> >::const_iterator idsIt = sc.hitsAndFractions().begin();
0689        idsIt != sc.hitsAndFractions().end();
0690        ++idsIt) {
0691     v1.push_back(idsIt->first);
0692   }
0693 
0694   DetId maxHitId;
0695 
0696   maxHitId = findMaxHit(v1, (EBhits), (EEhits));
0697   //maxHitId = findMaxHit(v1,EBhits,EEhits);
0698 
0699   EventsAfterCuts->Fill(6);
0700   if (maxHitId.null()) {
0701     std::cout << " Null " << std::endl;
0702     return;
0703   }
0704 
0705   int maxCC_Eta = 0;
0706   int maxCC_Phi = 0;
0707   int Zside = 0;
0708   if (maxHitId.subdetId() != 1) {
0709     maxCC_Eta = ((EEDetId)maxHitId).ix();
0710     maxCC_Phi = ((EEDetId)maxHitId).iy();
0711     Zside = ((EEDetId)maxHitId).zside();
0712     //    std::cout<<" ++++++++ Zside "<<Zside<<std::endl;
0713   } else {
0714     maxCC_Eta = ((EBDetId)maxHitId).ieta();
0715     maxCC_Phi = ((EBDetId)maxHitId).iphi();
0716   }
0717 
0718   //   if(maxCC_Eta>maxeta_ ) ;
0719   //   if(maxCC_Eta<mineta_ )  ;
0720 
0721   // number of events per crystal is set
0722   //   eventcrystal[maxCC_Eta][maxCC_Phi]+=1;
0723   //   if(eventcrystal[maxCC_Eta][maxCC_Phi] > numevent_) ;
0724 
0725   // fill cluster energy
0726   std::vector<float> energy;
0727   float energy3x3 = 0.;
0728   float energy5x5 = 0.;
0729   //Should be moved to cfg file!
0730   int ClusterSize = clusterSize_;
0731 
0732   const CaloSubdetectorTopology* topology = theCaloTopology_->getSubdetectorTopology(DetId::Ecal, maxHitId.subdetId());
0733   std::vector<DetId> NxNaroundMax = topology->getWindow(maxHitId, ClusterSize, ClusterSize);
0734   //ToCompute 3x3
0735   std::vector<DetId> S9aroundMax = topology->getWindow(maxHitId, 3, 3);
0736 
0737   EventsAfterCuts->Fill(7);
0738   if ((int)NxNaroundMax.size() != ClusterSize * ClusterSize)
0739     return;
0740   EventsAfterCuts->Fill(8);
0741   if (S9aroundMax.size() != 9)
0742     return;
0743 
0744   //   std::cout<<" ******** New Event "<<std::endl;
0745 
0746   EventsAfterCuts->Fill(9);
0747   for (int icry = 0; icry < ClusterSize * ClusterSize; icry++) {
0748     if (NxNaroundMax[icry].subdetId() == EcalBarrel) {
0749       EBRecHitCollection::const_iterator itrechit;
0750       itrechit = EBhits->find(NxNaroundMax[icry]);
0751       if (itrechit == EBhits->end()) {
0752         //  std::cout << "EB DetId not in e25" << std::endl;
0753         energy.push_back(0.);
0754         energy5x5 += 0.;
0755         continue;
0756       }
0757 
0758       if (edm::isNotFinite(itrechit->energy())) {
0759         std::cout << " nan energy " << std::endl;
0760         return;
0761       }
0762       energy.push_back(itrechit->energy());
0763       energy5x5 += itrechit->energy();
0764 
0765       //Summing in 3x3 to cut later on:
0766       for (int tt = 0; tt < 9; tt++) {
0767         if (NxNaroundMax[icry] == S9aroundMax[tt])
0768           energy3x3 += itrechit->energy();
0769       }
0770     } else {
0771       EERecHitCollection::const_iterator itrechit;
0772 
0773       itrechit = EEhits->find(NxNaroundMax[icry]);
0774 
0775       if (itrechit == EEhits->end()) {
0776         //  std::cout << "EE DetId not in e25" << std::endl;
0777         //  std::cout<<" ******** putting 0 "<<std::endl;
0778         energy.push_back(0.);
0779         energy5x5 += 0.;
0780         continue;
0781       }
0782 
0783       if (edm::isNotFinite(itrechit->energy())) {
0784         std::cout << " nan energy " << std::endl;
0785         return;
0786       }
0787       energy.push_back(itrechit->energy());
0788       energy5x5 += itrechit->energy();
0789 
0790       //Summing in 3x3 to cut later on:
0791       for (int tt = 0; tt < 9; tt++) {
0792         if (NxNaroundMax[icry] == S9aroundMax[tt])
0793           energy3x3 += itrechit->energy();
0794       }
0795     }
0796   }
0797   //  if((read_events-50)%10000 ==0)cout << "++++++++++++ENERGY 5x5 " <<  energy5x5 << std::endl;
0798   EventsAfterCuts->Fill(10);
0799   //  std::cout<<" ******** NxNaroundMax.size() "<<NxNaroundMax.size()<<std::endl;
0800   //  std::cout<<" ******** energy.size() "<<energy.size()<<std::endl;
0801   if ((int)energy.size() != ClusterSize * ClusterSize)
0802     return;
0803 
0804   if (maxHitId.subdetId() == EcalBarrel) {
0805     GeneralMapBeforePt->Fill(maxCC_Eta, maxCC_Phi);
0806   } else {
0807     if (Zside < 0) {
0808       GeneralMapEndCapMinusBeforePt->Fill(maxCC_Eta, maxCC_Phi);
0809     } else {
0810       GeneralMapEndCapPlusBeforePt->Fill(maxCC_Eta, maxCC_Phi);
0811     }
0812   }
0813 
0814   EventsAfterCuts->Fill(11);
0815   if (highestElePt < elePt_)
0816     return;
0817 
0818   if (maxHitId.subdetId() == EcalBarrel) {
0819     GeneralMap->Fill(maxCC_Eta, maxCC_Phi);
0820   } else {
0821     if (Zside < 0) {
0822       GeneralMapEndCapMinus->Fill(maxCC_Eta, maxCC_Phi);
0823     } else {
0824       GeneralMapEndCapPlus->Fill(maxCC_Eta, maxCC_Phi);
0825     }
0826   }
0827 
0828   EventsAfterCuts->Fill(12);
0829   if (highPtElectron.classification() != elecclass_ && elecclass_ != -1)
0830     return;
0831 
0832   float Ptrack_in =
0833       sqrt(pow(highPtElectron.trackMomentumAtVtx().X(), 2) + pow(highPtElectron.trackMomentumAtVtx().Y(), 2) +
0834            pow(highPtElectron.trackMomentumAtVtx().Z(), 2));
0835 
0836   float UncorrectedPatCalo =
0837       sqrt(pow(highPtElectron.trackMomentumAtCalo().X(), 2) + pow(highPtElectron.trackMomentumAtCalo().Y(), 2) +
0838            pow(highPtElectron.trackMomentumAtCalo().Z(), 2));
0839 
0840   float Ptrack_out =
0841       sqrt(pow(highPtElectron.trackMomentumOut().X(), 2) + pow(highPtElectron.trackMomentumOut().Y(), 2) +
0842            pow(highPtElectron.trackMomentumOut().Z(), 2));
0843 
0844   e9NoCuts->Fill(energy3x3);
0845   e25NoCuts->Fill(energy5x5);
0846   e9Overe25NoCuts->Fill(energy3x3 / energy5x5);
0847   scENoCuts->Fill(sc.energy());
0848 
0849   trPNoCuts->Fill(UncorrectedPatCalo);
0850 
0851   EoPNoCuts->Fill(highPtElectron.eSuperClusterOverP());
0852   e25OverScENoCuts->Fill(energy5x5 / sc.energy());
0853 
0854   E25oPNoCuts->Fill(energy5x5 / UncorrectedPatCalo);
0855 
0856   PinOverPoutNoCuts->Fill(
0857       sqrt(pow(highPtElectron.trackMomentumAtVtx().X(), 2) + pow(highPtElectron.trackMomentumAtVtx().Y(), 2) +
0858            pow(highPtElectron.trackMomentumAtVtx().Z(), 2)) /
0859       sqrt(pow(highPtElectron.trackMomentumOut().X(), 2) + pow(highPtElectron.trackMomentumOut().Y(), 2) +
0860            pow(highPtElectron.trackMomentumOut().Z(), 2)));
0861   eSeedOverPoutNoCuts->Fill(highPtElectron.eSuperClusterOverP());
0862 
0863   MapCor1NoCuts->Fill(energy5x5 / UncorrectedPatCalo, energy5x5 / Ptrack_in);
0864   MapCor2NoCuts->Fill(energy5x5 / UncorrectedPatCalo, highPtElectron.eSuperClusterOverP());
0865   MapCor3NoCuts->Fill(energy5x5 / UncorrectedPatCalo, Ptrack_out / Ptrack_in);
0866   MapCor4NoCuts->Fill(energy5x5 / UncorrectedPatCalo, energy5x5 / highPtElectron.p());
0867   MapCor5NoCuts->Fill(energy5x5 / UncorrectedPatCalo, UncorrectedPatCalo / Ptrack_out);
0868   MapCor6NoCuts->Fill(Ptrack_out / Ptrack_in, energy5x5 / Ptrack_in);
0869   MapCor7NoCuts->Fill(Ptrack_out / Ptrack_in, UncorrectedPatCalo / Ptrack_out);
0870   MapCor8NoCuts->Fill(energy5x5 / Ptrack_in, UncorrectedPatCalo / Ptrack_out);
0871   MapCor9NoCuts->Fill(energy5x5 / UncorrectedPatCalo, highPtElectron.eSeedClusterOverPout());
0872   MapCor10NoCuts->Fill(highPtElectron.eSeedClusterOverPout(), Ptrack_out / Ptrack_in);
0873   MapCor11NoCuts->Fill(highPtElectron.eSeedClusterOverPout(), energy5x5 / Ptrack_in);
0874 
0875   PinMinPoutNoCuts->Fill((Ptrack_in - Ptrack_out) / Ptrack_in);
0876 
0877   Error1NoCuts->Fill(highPtElectron.trackMomentumError() / Ptrack_in);
0878   Error2NoCuts->Fill(highPtElectron.trackMomentumError() / Ptrack_out);
0879   Error3NoCuts->Fill(highPtElectron.trackMomentumError() / UncorrectedPatCalo);
0880   eSeedOverPout2NoCuts->Fill(highPtElectron.eSeedClusterOverPout());
0881 
0882   hadOverEmNoCuts->Fill(highPtElectron.hadronicOverEm());
0883 
0884   //Cuts!
0885   if ((energy3x3 / energy5x5) < cut1_)
0886     return;
0887   if ((Ptrack_out / Ptrack_in) < cut2_ || (Ptrack_out / Ptrack_in) > cut3_)
0888     return;
0889   if ((energy5x5 / Ptrack_in) < cutEPin1_ || (energy5x5 / Ptrack_in) > cutEPin2_)
0890     return;
0891   //    if(!highPtElectron.ecalDriven())return;
0892   //    if(!highPtElectron.passingCutBasedPreselection())return;
0893 
0894   // //  Apply Pietro cuts:
0895   //    EventsAfterCuts->Fill(13);
0896   //    //Module 1
0897   //    if(maxHitId.subdetId() == EcalBarrel){
0898   //      //Module 1
0899   //      if(maxCC_Eta <= 25){
0900   //        if(highPtElectron.eSuperClusterOverP()>1.05 || highPtElectron.eSuperClusterOverP()<0.95)return ;
0901   //        if(highPtElectron.eSeedClusterOverPout()>1.4 || highPtElectron.eSeedClusterOverPout()<0.90)return ;
0902   //        if((Ptrack_in- Ptrack_out) / Ptrack_in <-0.05 || (Ptrack_in- Ptrack_out) / Ptrack_in >0.2)return ;
0903   //      }else{
0904   //        //Module 2
0905   //        if( maxCC_Eta > 25&& maxCC_Eta <= 45){
0906   //          if(highPtElectron.eSuperClusterOverP()>1.05 || highPtElectron.eSuperClusterOverP()<0.95)return ;
0907   //          if(highPtElectron.eSeedClusterOverPout()>1.25 || highPtElectron.eSeedClusterOverPout()<0.90)return ;
0908   //          if((Ptrack_in- Ptrack_out) / Ptrack_in <-0.05 || (Ptrack_in- Ptrack_out) / Ptrack_in >0.2)return ;
0909   //        }else{
0910   //        //Module 3
0911   //          if( maxCC_Eta > 45&& maxCC_Eta <= 65){
0912   //        if(highPtElectron.eSuperClusterOverP()>1.05 || highPtElectron.eSuperClusterOverP()<0.95)return ;
0913   //        if(highPtElectron.eSeedClusterOverPout()>1.15 || highPtElectron.eSeedClusterOverPout()<0.90)return ;
0914   //        if((Ptrack_in- Ptrack_out) / Ptrack_in <-0.05 || (Ptrack_in- Ptrack_out) / Ptrack_in >0.15)return ;
0915   //          }else{
0916   //          if( maxCC_Eta > 65&& maxCC_Eta <= 85){
0917   //        if(highPtElectron.eSuperClusterOverP()>1.05 || highPtElectron.eSuperClusterOverP()<0.95)return ;
0918   //        if(highPtElectron.eSeedClusterOverPout()>1.15 || highPtElectron.eSeedClusterOverPout()<0.90)return ;
0919   //        if((Ptrack_in- Ptrack_out) / Ptrack_in <-0.05 || (Ptrack_in- Ptrack_out) / Ptrack_in >0.15)return ;
0920   //          }else{
0921   //        return;
0922   //          }
0923   //          }
0924   //        }
0925   //      }
0926   //    }else{
0927   //      //EndCapMinus Side:
0928   //      //EndCapPlus Side:
0929   //      int iR = sqrt((maxCC_Eta-50)*(maxCC_Eta-50) + (maxCC_Phi-50)*(maxCC_Phi-50));
0930   //      if( iR >= 22&& iR < 27){
0931   //        if(highPtElectron.eSuperClusterOverP()>1.05 || highPtElectron.eSuperClusterOverP()<0.95)return ;
0932   //        if(highPtElectron.eSeedClusterOverPout()>1.15 || highPtElectron.eSeedClusterOverPout()<0.90)return ;
0933   //        if((Ptrack_in- Ptrack_out) / Ptrack_in <-0.05 || (Ptrack_in- Ptrack_out) / Ptrack_in >0.2)return ;
0934   //      }else{
0935   //        if( iR >= 27&& iR < 32){
0936   //          if(highPtElectron.eSuperClusterOverP()>1.1 || highPtElectron.eSuperClusterOverP()<0.95)return ;
0937   //          if(highPtElectron.eSeedClusterOverPout()>1.25 || highPtElectron.eSeedClusterOverPout()<0.90)return ;
0938   //          if((Ptrack_in- Ptrack_out) / Ptrack_in <-0.05 || (Ptrack_in- Ptrack_out) / Ptrack_in >0.2)return ;
0939   //        }else{
0940   //          if( iR >= 32&& iR < 37){
0941   //        if(highPtElectron.eSuperClusterOverP()>1.05 || highPtElectron.eSuperClusterOverP()<0.95)return ;
0942   //        if(highPtElectron.eSeedClusterOverPout()>1.15 || highPtElectron.eSeedClusterOverPout()<0.90)return ;
0943   //        if((Ptrack_in- Ptrack_out) / Ptrack_in <-0.05 || (Ptrack_in- Ptrack_out) / Ptrack_in >0.2)return ;
0944   //          }else{
0945   //        if( iR >= 37&& iR < 42){
0946   //          if(highPtElectron.eSuperClusterOverP()>1.1 || highPtElectron.eSuperClusterOverP()<0.95)return ;
0947   //          if(highPtElectron.eSeedClusterOverPout()>1.15 || highPtElectron.eSeedClusterOverPout()<0.90)return ;
0948   //          if((Ptrack_in- Ptrack_out) / Ptrack_in <-0.05 || (Ptrack_in- Ptrack_out) / Ptrack_in >0.15)return ;
0949   //        }else{
0950   //          if( iR >= 42){
0951   //            if(highPtElectron.eSuperClusterOverP()>1.05 || highPtElectron.eSuperClusterOverP()<0.95)return ;
0952   //            if(highPtElectron.eSeedClusterOverPout()>1.15 || highPtElectron.eSeedClusterOverPout()<0.90)return ;
0953   //          if((Ptrack_in- Ptrack_out) / Ptrack_in <-0.05 || (Ptrack_in- Ptrack_out) / Ptrack_in >0.15)return ;
0954   //          }
0955   //        }
0956   //          }
0957   //        }
0958   //      }
0959   //    }
0960 
0961   if (maxHitId.subdetId() == EcalBarrel) {
0962     E25oPvsEta->Fill(maxCC_Eta, energy5x5 / UncorrectedPatCalo);
0963   } else {
0964     float Radius = sqrt((maxCC_Eta) * (maxCC_Eta) + (maxCC_Phi) * (maxCC_Phi));
0965     if (Zside < 0) {
0966       E25oPvsEtaEndCapMinus->Fill(Radius, energy5x5 / UncorrectedPatCalo);
0967     } else {
0968       E25oPvsEtaEndCapPlus->Fill(Radius, energy5x5 / UncorrectedPatCalo);
0969     }
0970   }
0971   e9->Fill(energy3x3);
0972   e25->Fill(energy5x5);
0973   e9Overe25->Fill(energy3x3 / energy5x5);
0974   scE->Fill(sc.energy());
0975   trP->Fill(UncorrectedPatCalo);
0976 
0977   EoP->Fill(highPtElectron.eSuperClusterOverP());
0978   e25OverScE->Fill(energy5x5 / sc.energy());
0979 
0980   E25oP->Fill(energy5x5 / UncorrectedPatCalo);
0981 
0982   if (maxHitId.subdetId() == EcalBarrel) {
0983     Map->Fill(maxCC_Eta, maxCC_Phi);
0984   } else {
0985     if (Zside < 0) {
0986       MapEndCapMinus->Fill(maxCC_Eta, maxCC_Phi);
0987     } else {
0988       MapEndCapPlus->Fill(maxCC_Eta, maxCC_Phi);
0989     }
0990   }
0991 
0992   PinOverPout->Fill(sqrt(pow(highPtElectron.trackMomentumAtVtx().X(), 2) +
0993                          pow(highPtElectron.trackMomentumAtVtx().Y(), 2) +
0994                          pow(highPtElectron.trackMomentumAtVtx().Z(), 2)) /
0995                     sqrt(pow(highPtElectron.trackMomentumOut().X(), 2) + pow(highPtElectron.trackMomentumOut().Y(), 2) +
0996                          pow(highPtElectron.trackMomentumOut().Z(), 2)));
0997   eSeedOverPout->Fill(highPtElectron.eSuperClusterOverP());
0998 
0999   MapCor1->Fill(energy5x5 / UncorrectedPatCalo, energy5x5 / Ptrack_in);
1000   MapCor2->Fill(energy5x5 / UncorrectedPatCalo, highPtElectron.eSuperClusterOverP());
1001   MapCor3->Fill(energy5x5 / UncorrectedPatCalo, Ptrack_out / Ptrack_in);
1002   MapCor4->Fill(energy5x5 / UncorrectedPatCalo, energy5x5 / highPtElectron.p());
1003   MapCor5->Fill(energy5x5 / UncorrectedPatCalo, UncorrectedPatCalo / Ptrack_out);
1004   MapCor6->Fill(Ptrack_out / Ptrack_in, energy5x5 / Ptrack_in);
1005   MapCor7->Fill(Ptrack_out / Ptrack_in, UncorrectedPatCalo / Ptrack_out);
1006   MapCor8->Fill(energy5x5 / Ptrack_in, UncorrectedPatCalo / Ptrack_out);
1007   MapCor9->Fill(energy5x5 / UncorrectedPatCalo, highPtElectron.eSeedClusterOverPout());
1008   MapCor10->Fill(highPtElectron.eSeedClusterOverPout(), Ptrack_out / Ptrack_in);
1009   MapCor11->Fill(highPtElectron.eSeedClusterOverPout(), energy5x5 / Ptrack_in);
1010 
1011   PinMinPout->Fill((Ptrack_in - Ptrack_out) / Ptrack_in);
1012 
1013   Error1->Fill(highPtElectron.trackMomentumError() / Ptrack_in);
1014   Error2->Fill(highPtElectron.trackMomentumError() / Ptrack_out);
1015   Error3->Fill(highPtElectron.trackMomentumError() / UncorrectedPatCalo);
1016 
1017   eSeedOverPout2->Fill(highPtElectron.eSeedClusterOverPout());
1018   hadOverEm->Fill(highPtElectron.hadronicOverEm());
1019 
1020   UnivEventIds.push_back(NxNaroundMax);
1021   EventMatrix.push_back(energy);
1022   EnergyVector.push_back(UncorrectedPatCalo);
1023 
1024   EventsAfterCuts->Fill(14);
1025 
1026   if (!highPtElectron.ecalDrivenSeed())
1027     EventsAfterCuts->Fill(15);
1028 
1029   return;
1030 }
1031 
1032 DEFINE_FWK_MODULE(ElectronCalibrationUniv);