Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-08-30 22:38:32

0001 
0002 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0003 
0004 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0005 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0006 #include "DataFormats/DetId/interface/DetId.h"
0007 #include "CondFormats/EcalObjects/interface/EcalIntercalibConstants.h"
0008 #include "CondFormats/DataRecord/interface/EcalIntercalibConstantsRcd.h"
0009 #include "Calibration/Tools/interface/calibXMLwriter.h"
0010 #include "Calibration/Tools/interface/CalibrationCluster.h"
0011 #include "Calibration/Tools/interface/HouseholderDecomposition.h"
0012 #include "Calibration/Tools/interface/MinL3Algorithm.h"
0013 #include "Calibration/EcalCalibAlgos/interface/ElectronCalibration.h"
0014 #include "DataFormats/EgammaCandidates/interface/Electron.h"
0015 #include "DataFormats/EgammaCandidates/interface/ElectronFwd.h"
0016 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0017 #include "DataFormats/TrackReco/interface/Track.h"
0018 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0019 #include "DataFormats/TrackReco/interface/TrackExtraFwd.h"
0020 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0021 #include "FWCore/Utilities/interface/isFinite.h"
0022 #include "TFile.h"
0023 #include "TH1.h"
0024 #include "TH2.h"
0025 #include "TF1.h"
0026 #include "TRandom.h"
0027 
0028 #include <iostream>
0029 #include <string>
0030 #include <stdexcept>
0031 #include <vector>
0032 
0033 ElectronCalibration::ElectronCalibration(const edm::ParameterSet& iConfig) {
0034   rootfile_ = iConfig.getParameter<std::string>("rootfile");
0035   recHitToken_ = consumes<EBRecHitCollection>(iConfig.getParameter<edm::InputTag>("ebRecHitsLabel"));
0036   electronToken_ = consumes<reco::GsfElectronCollection>(iConfig.getParameter<edm::InputTag>("electronLabel"));
0037   trackLabel_ = iConfig.getParameter<edm::InputTag>("trackLabel");
0038   calibAlgo_ = iConfig.getParameter<std::string>("CALIBRATION_ALGO");
0039   edm::LogVerbatim("EcalCalib") << " The used Algorithm is  " << calibAlgo_;
0040   keventweight_ = iConfig.getParameter<int>("keventweight");
0041   ClusterSize_ = iConfig.getParameter<int>("Clustersize");
0042   ElePt_ = iConfig.getParameter<double>("ElePt");
0043   maxeta_ = iConfig.getParameter<int>("maxeta");
0044   mineta_ = iConfig.getParameter<int>("mineta");
0045   maxphi_ = iConfig.getParameter<int>("maxphi");
0046   minphi_ = iConfig.getParameter<int>("minphi");
0047   cut1_ = iConfig.getParameter<double>("cut1");
0048   cut2_ = iConfig.getParameter<double>("cut2");
0049   cut3_ = iConfig.getParameter<double>("cut3");
0050   elecclass_ = iConfig.getParameter<int>("elecclass");
0051   edm::LogVerbatim("EcalCalib") << " The electronclass is " << elecclass_;
0052   numevent_ = iConfig.getParameter<int>("numevent");
0053   miscalibfile_ = iConfig.getParameter<std::string>("miscalibfile");
0054 
0055   cutEPCalo1_ = iConfig.getParameter<double>("cutEPCaloMin");
0056   cutEPCalo2_ = iConfig.getParameter<double>("cutEPCaloMax");
0057   cutEPin1_ = iConfig.getParameter<double>("cutEPinMin");
0058   cutEPin2_ = iConfig.getParameter<double>("cutEPinMax");
0059   cutCalo1_ = iConfig.getParameter<double>("cutCaloMin");
0060   cutCalo2_ = iConfig.getParameter<double>("cutCaloMax");
0061 
0062   cutESeed_ = iConfig.getParameter<double>("cutESeed");
0063 }
0064 
0065 ElectronCalibration::~ElectronCalibration() {}
0066 
0067 //========================================================================
0068 void ElectronCalibration::beginJob() {
0069   //========================================================================
0070   f = new TFile(rootfile_.c_str(), "RECREATE");
0071 
0072   // Book histograms
0073   e9 = new TH1F("e9", "E9 energy", 300, 0., 150.);
0074   e25 = new TH1F("e25", "E25 energy", 300, 0., 150.);
0075   scE = new TH1F("scE", "SC energy", 300, 0., 150.);
0076   trP = new TH1F("trP", "Trk momentum", 300, 0., 150.);
0077   EoP = new TH1F("EoP", "EoP", 600, 0., 3.);
0078   EoP_all = new TH1F("EoP_all", "EoP_all", 600, 0., 3.);
0079 
0080   if (elecclass_ == 0 || elecclass_ == -1) {
0081     calibs = new TH1F("calib", "Calibration constants", 4000, 0.5, 2.);
0082   } else {
0083     calibs = new TH1F("calib", "Calibration constants", 800, 0.5, 2.);
0084   }
0085 
0086   e25OverScE = new TH1F("e25OverscE", "E25 / SC energy", 400, 0., 2.);
0087   E25oP = new TH1F("E25oP", "E25 / P", 1200, 0., 1.5);
0088 
0089   Map = new TH2F("Map", "Nb Events in Crystal", 85, 1, 85, 70, 5, 75);
0090   e9Overe25 = new TH1F("e9Overe25", "E9 / E25", 400, 0., 2.);
0091   Map3Dcalib = new TH2F("3Dcalib", "3Dcalib", 85, 1, 85, 70, 5, 75);
0092 
0093   MapCor1 = new TH2F("MapCor1", "Correlation E25/Pcalo versus E25/Pin", 100, 0., 5., 100, 0., 5.);
0094   MapCor2 = new TH2F("MapCor2", "Correlation E25/Pcalo versus E/P", 100, 0., 5., 100, 0., 5.);
0095   MapCor3 = new TH2F("MapCor3", "Correlation E25/Pcalo versus Pout/Pin", 100, 0., 5., 100, 0., 5.);
0096   MapCor4 = new TH2F("MapCor4", "Correlation E25/Pcalo versus E25/highestP", 100, 0., 5., 100, 0., 5.);
0097   MapCor5 = new TH2F("MapCor5", "Correlation E25/Pcalo versus Pcalo/Pout", 100, 0., 5., 100, 0., 5.);
0098   MapCor6 = new TH2F("MapCor6", "Correlation Pout/Pin versus E25/Pin", 100, 0., 5., 100, 0., 5.);
0099   MapCor7 = new TH2F("MapCor7", "Correlation Pout/Pin versus Pcalo/Pout", 100, 0., 5., 100, 0., 5.);
0100   MapCor8 = new TH2F("MapCor8", "Correlation E25/Pin versus Pcalo/Pout", 100, 0., 5., 100, 0., 5.);
0101   MapCor9 = new TH2F("MapCor9", "Correlation  E25/Pcalo versus Eseed/Pout", 100, 0., 5., 100, 0., 5.);
0102   MapCor10 = new TH2F("MapCor10", "Correlation Eseed/Pout versus Pout/Pin", 100, 0., 5., 100, 0., 5.);
0103   MapCor11 = new TH2F("MapCor11", "Correlation Eseed/Pout versus E25/Pin", 100, 0., 5., 100, 0., 5.);
0104   MapCorCalib =
0105       new TH2F("MapCorCalib", "Correlation Miscalibration versus Calibration constants", 100, 0.5, 1.5, 100, 0.5, 1.5);
0106 
0107   PinMinPout = new TH1F("PinMinPout", "(Pin - Pout)/Pin", 600, -2.0, 2.0);
0108 
0109   if (elecclass_ == 0 || elecclass_ == -1) {
0110     calibinter = new TH1F("calibinter", "internal calibration constants", 2000, 0.5, 2.);
0111     PinOverPout = new TH1F("PinOverPout", "pinOverpout", 600, 0., 3.);
0112     eSeedOverPout = new TH1F("eSeedOverPout", "eSeedOverpout ", 600, 0., 3.);
0113     MisCalibs = new TH1F("MisCalibs", "Miscalibration constants", 4000, 0.5, 2.);
0114     RatioCalibs = new TH1F("RatioCalibs", "Ratio in Calibration Constants", 4000, 0.5, 2.0);
0115     DiffCalibs = new TH1F("DiffCalibs", "Difference in Calibration constants", 4000, -1.0, 1.0);
0116   } else {
0117     calibinter = new TH1F("calibinter", "internal calibration constants", 400, 0.5, 2.);
0118     PinOverPout = new TH1F("PinOverPout", "pinOverpout", 600, 0., 3.);
0119     eSeedOverPout = new TH1F("eSeedOverPout", "eSeedOverpout ", 600, 0., 3.);
0120     MisCalibs = new TH1F("MisCalibs", "Miscalibration constants", 800, 0.5, 2.);
0121     RatioCalibs = new TH1F("RatioCalibs", "Ratio in Calibration Constants", 800, 0.5, 2.0);
0122     DiffCalibs = new TH1F("DiffCalibs", "Difference in Calibration constants", 800, -1.0, 1.0);
0123   }
0124   Error1 = new TH1F("Error1", "DeltaP/Pin", 800, -1.0, 1.0);
0125   Error2 = new TH1F("Error2", "DeltaP/Pout", 800, -1.0, 1.0);
0126   Error3 = new TH1F("Error3", "DeltaP/Pcalo", 800, -1.0, 1.0);
0127   eSeedOverPout2 = new TH1F("eSeedOverPout2", "eSeedOverpout (No Supercluster)", 600, 0., 4.);
0128   hadOverEm = new TH1F("hadOverEm", "Had/EM distribution", 600, -2., 2.);
0129 
0130   // Book histograms
0131   Map3DcalibNoCuts = new TH2F("3DcalibNoCuts", "3Dcalib (Before Cuts)", 85, 1, 85, 70, 5, 75);
0132   e9NoCuts = new TH1F("e9NoCuts", "E9 energy (Before Cuts)", 300, 0., 150.);
0133   e25NoCuts = new TH1F("e25NoCuts", "E25 energy (Before Cuts)", 300, 0., 150.);
0134   scENoCuts = new TH1F("scENoCuts", "SC energy (Before Cuts)", 300, 0., 150.);
0135   trPNoCuts = new TH1F("trPNoCuts", "Trk momentum (Before Cuts)", 300, 0., 150.);
0136   EoPNoCuts = new TH1F("EoPNoCuts", "EoP (Before Cuts)", 600, 0., 3.);
0137   if (elecclass_ == 0 || elecclass_ == -1) {
0138     calibsNoCuts = new TH1F("calibNoCuts", "Calibration constants (Before Cuts)", 4000, 0., 2.);
0139   } else {
0140     calibsNoCuts = new TH1F("calibNoCuts", "Calibration constants (Before Cuts)", 800, 0., 2.);
0141   }
0142   e25OverScENoCuts = new TH1F("e25OverscENoCuts", "E25 / SC energy (Before Cuts)", 400, 0., 2.);
0143   E25oPNoCuts = new TH1F("E25oPNoCuts", "E25 / P (Before Cuts)", 1200, 0., 1.5);
0144   MapNoCuts = new TH2F("MapNoCuts", "Nb Events in Crystal (Before Cuts)", 85, 1, 85, 70, 5, 75);
0145   e9Overe25NoCuts = new TH1F("e9Overe25NoCuts", "E9 / E25 (Before Cuts)", 400, 0., 2.);
0146   PinOverPoutNoCuts = new TH1F("PinOverPoutNoCuts", "pinOverpout (Before Cuts)", 600, 0., 3.);
0147   eSeedOverPoutNoCuts = new TH1F(" eSeedOverPoutNoCuts", "eSeedOverpout (Before Cuts) ", 600, 0., 4.);
0148   PinMinPoutNoCuts = new TH1F("PinMinPoutNoCuts", "(Pin - Pout)/Pin (Before Cuts)", 600, -2.0, 2.0);
0149 
0150   RatioCalibsNoCuts = new TH1F("RatioCalibsNoCuts", "Ratio in Calibration Constants (Before Cuts)", 4000, 0.5, 2.0);
0151   DiffCalibsNoCuts = new TH1F("DiffCalibsNoCuts", "Difference in Calibration constants (Before Cuts)", 4000, -1.0, 1.0);
0152   calibinterNoCuts = new TH1F("calibinterNoCuts", "internal calibration constants", 2000, 0.5, 2.);
0153 
0154   MapCor1NoCuts =
0155       new TH2F("MapCor1NoCuts", "Correlation E25/PatCalo versus E25/Pin (Before Cuts)", 100, 0., 5., 100, 0., 5.);
0156   MapCor2NoCuts =
0157       new TH2F("MapCor2NoCuts", "Correlation E25/PatCalo versus E/P (Before Cuts)", 100, 0., 5., 100, 0., 5.);
0158   MapCor3NoCuts =
0159       new TH2F("MapCor3NoCuts", "Correlation E25/PatCalo versus Pout/Pin (Before Cuts)", 100, 0., 5., 100, 0., 5.);
0160   MapCor4NoCuts =
0161       new TH2F("MapCor4NoCuts", "Correlation E25/PatCalo versus E25/highestP (Before Cuts)", 100, 0., 5., 100, 0., 5.);
0162   MapCor5NoCuts =
0163       new TH2F("MapCor5NoCuts", "Correlation E25/Pcalo versus Pcalo/Pout (Before Cuts)", 100, 0., 5., 100, 0., 5.);
0164   MapCor6NoCuts =
0165       new TH2F("MapCor6NoCuts", "Correlation Pout/Pin versus E25/Pin (Before Cuts)", 100, 0., 5., 100, 0., 5.);
0166   MapCor7NoCuts =
0167       new TH2F("MapCor7NoCuts", "Correlation Pout/Pin versus Pcalo/Pout (Before Cuts)", 100, 0., 5., 100, 0., 5.);
0168   MapCor8NoCuts =
0169       new TH2F("MapCor8NoCuts", "Correlation E25/Pin versus Pcalo/Pout (Before Cuts)", 100, 0., 5., 100, 0., 5.);
0170   MapCor9NoCuts =
0171       new TH2F("MapCor9NoCuts", "Correlation  E25/Pcalo versus Eseed/Pout (Before Cuts)", 100, 0., 5., 100, 0., 5.);
0172   MapCor10NoCuts =
0173       new TH2F("MapCor10NoCuts", "Correlation Eseed/Pout versus Pout/Pin (Before Cuts)", 100, 0., 5., 100, 0., 5.);
0174   MapCor11NoCuts =
0175       new TH2F("MapCor11NoCuts", "Correlation Eseed/Pout versus E25/Pin (Before Cuts)", 100, 0., 5., 100, 0., 5.);
0176   MapCorCalibNoCuts = new TH2F("MapCorCalibNoCuts",
0177                                "Correlation Miscalibration versus Calibration constants (Before Cuts)",
0178                                100,
0179                                0.,
0180                                3.,
0181                                100,
0182                                0.,
0183                                3.);
0184 
0185   Error1NoCuts = new TH1F("Eror1NoCuts", "DeltaP/Pin (Before Cuts)", 800, -1.0, 1.0);
0186   Error2NoCuts = new TH1F("Error2NoCuts", "DeltaP/Pout (Before Cuts)", 800, -1.0, 1.0);
0187   Error3NoCuts = new TH1F("Error3NoCuts", "DeltaP/Pcalo (Before Cuts)", 800, -1.0, 1.0);
0188   eSeedOverPout2NoCuts = new TH1F("eSeedOverPout2NoCuts", "eSeedOverpout (No Supercluster, Before Cuts)", 600, 0., 4.);
0189   hadOverEmNoCuts = new TH1F("hadOverEmNoCuts", "Had/EM distribution (Before Cuts)", 600, -2., 2.);
0190 
0191   //Book histograms after ESeed cut
0192   MapCor1ESeed =
0193       new TH2F("MapCor1ESeed", "Correlation E25/Pcalo versus E25/Pin (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
0194   MapCor2ESeed =
0195       new TH2F("MapCor2ESeed", "Correlation E25/Pcalo versus E/P (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
0196   MapCor3ESeed = new TH2F(
0197       "MapCor3ESeed", "Correlation E25/Pcalo versus Pout/Pin (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
0198   MapCor4ESeed = new TH2F(
0199       "MapCor4ESeed", "Correlation E25/Pcalo versus E25/highestP (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
0200   MapCor5ESeed = new TH2F(
0201       "MapCor5ESeed", "Correlation E25/Pcalo versus Pcalo/Pout (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
0202   MapCor6ESeed =
0203       new TH2F("MapCor6ESeed", "Correlation Pout/Pin versus E25/Pin (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
0204   MapCor7ESeed = new TH2F(
0205       "MapCor7ESeed", "Correlation Pout/Pin versus Pcalo/Pout (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
0206   MapCor8ESeed = new TH2F(
0207       "MapCor8ESeed", "Correlation E25/Pin versus Pcalo/Pout (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
0208   MapCor9ESeed = new TH2F(
0209       "MapCor9ESeed", "Correlation  E25/Pcalo versus Eseed/Pout (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
0210   MapCor10ESeed = new TH2F(
0211       "MapCor10ESeed", "Correlation Eseed/Pout versus Pout/Pin (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
0212   MapCor11ESeed = new TH2F(
0213       "MapCor11ESeed", "Correlation Eseed/Pout versus E25/Pin (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
0214 
0215   eSeedOverPout2ESeed =
0216       new TH1F("eSeedOverPout2ESeed", "eSeedOverpout (No Supercluster, after Eseed/Pout cut)", 600, 0., 4.);
0217 
0218   hadOverEmESeed = new TH1F("hadOverEmESeed", "Had/EM distribution (after Eseed/Pout cut)", 600, -2., 2.);
0219 
0220   //Book histograms without any cut
0221   GeneralMap = new TH2F("GeneralMap", "Map without any cuts", 85, 1, 85, 70, 5, 75);
0222 
0223   calibClusterSize = ClusterSize_;
0224   etaMin = mineta_;
0225   etaMax = maxeta_;
0226   phiMin = minphi_;
0227   phiMax = maxphi_;
0228   if (calibAlgo_ == "L3") {
0229     MyL3Algo1 = new MinL3Algorithm(keventweight_, calibClusterSize, etaMin, etaMax, phiMin, phiMax);
0230   } else {
0231     if (calibAlgo_ == "HH" || calibAlgo_ == "HHReg") {
0232       MyHH = new HouseholderDecomposition(calibClusterSize, etaMin, etaMax, phiMin, phiMax);
0233     } else {
0234       edm::LogError("EcalCalib") << " Name of Algorithm is not recognize " << calibAlgo_
0235                                  << " Should be either L3, HH or HHReg. Abort! ";
0236     }
0237   }
0238   read_events = 0;
0239 
0240   // get Region to be calibrated
0241   ReducedMap = calibCluster.getMap(etaMin, etaMax, phiMin, phiMax);
0242 
0243   oldCalibs.resize(ReducedMap.size(), 0.);
0244 
0245   // table is set to zero
0246   for (int phi = 0; phi < 360; phi++) {
0247     for (int eta = 0; eta < 171; eta++) {
0248       eventcrystal[eta][phi] = 0;
0249     }
0250   }
0251 
0252   edm::LogVerbatim("EcalCalib") << " Begin JOB ";
0253 }
0254 
0255 //========================================================================
0256 
0257 void ElectronCalibration::endJob() {
0258   //========================================================================
0259 
0260   int nIterations = 10;
0261   if (calibAlgo_ == "L3") {
0262     solution = MyL3Algo1->iterate(EventMatrix, MaxCCeta, MaxCCphi, EnergyVector, nIterations);
0263   } else {
0264     if (calibAlgo_ == "HH") {
0265       solution = MyHH->iterate(EventMatrix, MaxCCeta, MaxCCphi, EnergyVector, 1, false);
0266     } else {
0267       if (calibAlgo_ == "HHReg") {
0268         solution = MyHH->runRegional(EventMatrix, MaxCCeta, MaxCCphi, EnergyVector, 2);
0269       } else {
0270         edm::LogError("EcalCalib") << " Calibration not run due to problem in Algo Choice...";
0271         return;
0272       }
0273     }
0274   }
0275   for (int ii = 0; ii < (int)solution.size(); ii++) {
0276     edm::LogVerbatim("EcalCalib") << "solution[" << ii << "] = " << solution[ii];
0277     calibs->Fill(solution[ii]);
0278   }
0279 
0280   newCalibs.resize(ReducedMap.size(), 0.);
0281 
0282   calibXMLwriter write_calibrations;
0283 
0284   FILE* MisCalib;
0285   MisCalib = fopen(miscalibfile_.c_str(), "r");
0286   int fileStatus = 1;
0287   int eta = -1;
0288   int phi = -1;
0289   float coeff = -1;
0290 
0291   std::map<EBDetId, float> OldCoeff;
0292 
0293   while (fileStatus != EOF) {
0294     fileStatus = fscanf(MisCalib, "%d %d %f\n", &eta, &phi, &coeff);
0295     if (eta != -1 && phi != -1 && coeff != -1) {
0296       OldCoeff.insert(std::make_pair(EBDetId(eta, phi, EBDetId::ETAPHIMODE), coeff));
0297     }
0298   }
0299 
0300   fclose(MisCalib);
0301 
0302   int icry = 0;
0303   CalibrationCluster::CalibMap::iterator itmap;
0304   for (itmap = ReducedMap.begin(); itmap != ReducedMap.end(); itmap++) {
0305     newCalibs[icry] = solution[icry];
0306 
0307     write_calibrations.writeLine(itmap->first, newCalibs[icry]);
0308     float Compare = 1.;
0309     std::map<EBDetId, float>::iterator iter = OldCoeff.find(itmap->first);
0310     if (iter != OldCoeff.end())
0311       Compare = iter->second;
0312 
0313     if ((itmap->first).ieta() > mineta_ && (itmap->first).ieta() < maxeta_ && (itmap->first).iphi() > minphi_ &&
0314         (itmap->first).iphi() < maxphi_) {
0315       Map3Dcalib->Fill((itmap->first).ieta(), (itmap->first).iphi(), newCalibs[icry] * Compare);
0316       MisCalibs->Fill(Compare);
0317     }
0318     if ((itmap->first).ieta() < mineta_ + 2) {
0319       icry++;
0320       continue;
0321     }
0322     if ((itmap->first).ieta() > maxeta_ - 2) {
0323       icry++;
0324       continue;
0325     }
0326     if ((itmap->first).iphi() < minphi_ + 2) {
0327       icry++;
0328       continue;
0329     }
0330     if ((itmap->first).iphi() > maxphi_ - 2) {
0331       icry++;
0332       continue;
0333     }
0334 
0335     calibinter->Fill(newCalibs[icry]);
0336     DiffCalibs->Fill(newCalibs[icry] - 1. / Compare);
0337     RatioCalibs->Fill(newCalibs[icry] * Compare);
0338     MapCorCalib->Fill(1. / Compare, newCalibs[icry]);
0339     icry++;
0340   }
0341 
0342   if (calibAlgo_ == "L3") {
0343     solutionNoCuts =
0344         MyL3Algo1->iterate(EventMatrixNoCuts, MaxCCetaNoCuts, MaxCCphiNoCuts, EnergyVectorNoCuts, nIterations);
0345   } else {
0346     if (calibAlgo_ == "HH") {
0347       solutionNoCuts = MyHH->iterate(EventMatrixNoCuts, MaxCCetaNoCuts, MaxCCphiNoCuts, EnergyVectorNoCuts, 1, false);
0348     } else {
0349       if (calibAlgo_ == "HHReg") {
0350         solutionNoCuts = MyHH->runRegional(EventMatrixNoCuts, MaxCCetaNoCuts, MaxCCphiNoCuts, EnergyVectorNoCuts, 2);
0351       } else {
0352         edm::LogError("EcalCalib") << " Calibration not run due to problem in AlgoChoice...";
0353         return;
0354       }
0355     }
0356   }
0357   for (int ii = 0; ii < (int)solutionNoCuts.size(); ii++) {
0358     calibsNoCuts->Fill(solutionNoCuts[ii]);
0359   }
0360   int icryp = 0;
0361   CalibrationCluster::CalibMap::iterator itmapp;
0362   for (itmapp = ReducedMap.begin(); itmapp != ReducedMap.end(); itmapp++) {
0363     newCalibs[icryp] = solutionNoCuts[icryp];
0364     float Compare2 = 1.;
0365     std::map<EBDetId, float>::iterator iter2 = OldCoeff.find(itmapp->first);
0366     if (iter2 != OldCoeff.end())
0367       Compare2 = iter2->second;
0368 
0369     if ((itmapp->first).ieta() > mineta_ && (itmapp->first).ieta() < maxeta_ && (itmapp->first).iphi() > minphi_ &&
0370         (itmapp->first).iphi() < maxphi_)
0371       Map3DcalibNoCuts->Fill((itmapp->first).ieta(), (itmapp->first).iphi(), newCalibs[icryp] * Compare2);
0372     if ((itmapp->first).ieta() < mineta_ + 2) {
0373       icryp++;
0374       continue;
0375     }
0376     if ((itmapp->first).ieta() > maxeta_ - 2) {
0377       icryp++;
0378       continue;
0379     }
0380     if ((itmapp->first).iphi() < minphi_ + 2) {
0381       icryp++;
0382       continue;
0383     }
0384     if ((itmapp->first).iphi() > maxphi_ - 2) {
0385       icryp++;
0386       continue;
0387     }
0388     calibinterNoCuts->Fill(newCalibs[icryp]);
0389     DiffCalibsNoCuts->Fill(newCalibs[icryp] - 1. / (Compare2));
0390     RatioCalibsNoCuts->Fill(newCalibs[icryp] * Compare2);
0391     MapCorCalibNoCuts->Fill(1. / Compare2, newCalibs[icryp]);
0392     icryp++;
0393   }
0394 
0395   ////////////////////////       FINAL STATISTICS           ////////////////////
0396 
0397   edm::LogVerbatim("EcalCalib") << "\n************* STATISTICS **************";
0398   edm::LogVerbatim("EcalCalib") << " Events Studied " << read_events;
0399 
0400   /////////////////////////////////////////////////////////////////////////////
0401 
0402   f->Write();
0403 
0404   f->Close();
0405 }
0406 
0407 //=================================================================================
0408 EBDetId ElectronCalibration::findMaxHit(edm::Handle<EBRecHitCollection>& phits) {
0409   //=================================================================================
0410 
0411   EcalRecHitCollection ecrh = *phits;
0412   EcalRecHitCollection::iterator it;
0413   int count = 0;
0414   EBDetId save;
0415   float en_save = 0;
0416   for (it = ecrh.begin(); it != ecrh.end(); it++) {
0417     EBDetId p = EBDetId(it->id().rawId());
0418     if (it->energy() > en_save) {
0419       en_save = it->energy();
0420       save = p;
0421     }
0422     count++;
0423   }
0424   return save;
0425 }
0426 
0427 //=================================================================================
0428 EBDetId ElectronCalibration::findMaxHit2(const std::vector<DetId>& v1, const EBRecHitCollection* hits) {
0429   //=================================================================================
0430 
0431   double currEnergy = 0.;
0432   EBDetId maxHit;
0433 
0434   for (std::vector<DetId>::const_iterator idsIt = v1.begin(); idsIt != v1.end(); ++idsIt) {
0435     if (idsIt->subdetId() != 1)
0436       continue;
0437     EBRecHitCollection::const_iterator itrechit;
0438     itrechit = hits->find(*idsIt);
0439 
0440     if (itrechit == hits->end()) {
0441       edm::LogVerbatim("EcalCalib") << "ElectronCalibration::findMaxHit2: rechit not found! ";
0442       continue;
0443     }
0444     if (itrechit->energy() > currEnergy) {
0445       currEnergy = itrechit->energy();
0446       maxHit = *idsIt;
0447     }
0448   }
0449 
0450   return maxHit;
0451 }
0452 
0453 //=================================================================================
0454 void ElectronCalibration::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0455   //=================================================================================
0456   using namespace edm;
0457 
0458   // Get EBRecHits
0459   const Handle<EBRecHitCollection>& phits = iEvent.getHandle(recHitToken_);
0460   if (!phits.isValid()) {
0461     edm::LogError("EcalCalib") << "Error! can't get the product EBRecHitCollection: ";
0462   }
0463 
0464   const EBRecHitCollection* hits = phits.product();  // get a ptr to the product
0465 
0466   // Get pixelElectrons
0467   const Handle<reco::GsfElectronCollection>& pElectrons = iEvent.getHandle(electronToken_);
0468   if (!pElectrons.isValid()) {
0469     edm::LogError("EcalCalib") << "Error! can't get the product ElectronCollection: ";
0470   }
0471 
0472   const reco::GsfElectronCollection* electronCollection = pElectrons.product();
0473   read_events++;
0474   if (read_events % 1000 == 0)
0475     edm::LogVerbatim("EcalCalib") << "read_events = " << read_events << std::endl;
0476 
0477   if (!hits)
0478     return;
0479   if (hits->empty())
0480     return;
0481   if (!electronCollection)
0482     return;
0483   if (electronCollection->empty())
0484     return;
0485 
0486   ////////////////////////////////////////////////////////////////////////////////////////
0487   //                          START HERE....
0488   ///////////////////////////////////////////////////////////////////////////////////////
0489   reco::GsfElectronCollection::const_iterator eleIt = electronCollection->begin();
0490 
0491   reco::GsfElectron highPtElectron;
0492 
0493   float highestElePt = 0.;
0494   bool found = false;
0495   for (eleIt = electronCollection->begin(); eleIt != electronCollection->end(); eleIt++) {
0496     //Comments
0497     if (fabs(eleIt->eta()) > (maxeta_ + 3) * 0.0175)
0498       continue;
0499     if (eleIt->eta() < (mineta_ - 3) * 0.0175)
0500       continue;
0501 
0502     if (eleIt->pt() > highestElePt) {
0503       highestElePt = eleIt->pt();
0504       highPtElectron = *eleIt;
0505       found = true;
0506     }
0507   }
0508   if (highestElePt < ElePt_)
0509     return;
0510   if (!found)
0511     return;
0512   const reco::SuperCluster& sc = *(highPtElectron.superCluster());
0513   if (fabs(sc.eta()) > (maxeta_ + 3) * 0.0175) {
0514     edm::LogVerbatim("EcalCalib") << "++++ Problem with electron, electron eta is " << highPtElectron.eta()
0515                                   << " while SC is " << sc.eta() << std::endl;
0516     return;
0517   }
0518 
0519   std::vector<DetId> v1;
0520   //Loop to fill the vector of DetIds
0521   for (std::vector<std::pair<DetId, float> >::const_iterator idsIt = sc.hitsAndFractions().begin();
0522        idsIt != sc.hitsAndFractions().end();
0523        ++idsIt) {
0524     v1.push_back(idsIt->first);
0525   }
0526 
0527   //Change function name
0528   EBDetId maxHitId;
0529 
0530   maxHitId = findMaxHit2(v1, hits);
0531 
0532   if (maxHitId.null()) {
0533     edm::LogVerbatim("EcalCalib") << " Null " << std::endl;
0534     return;
0535   }
0536 
0537   int maxCC_Eta = maxHitId.ieta();
0538   int maxCC_Phi = maxHitId.iphi();
0539 
0540   if (maxCC_Eta > maxeta_)
0541     return;
0542   if (maxCC_Eta < mineta_)
0543     return;
0544   if (maxCC_Phi > maxphi_)
0545     return;
0546   if (maxCC_Phi < minphi_)
0547     return;
0548 
0549   // number of events per crystal is set
0550   if (numevent_ > 0) {
0551     eventcrystal[maxCC_Eta + 85][maxCC_Phi - 1] += 1;
0552     if (eventcrystal[maxCC_Eta + 85][maxCC_Phi - 1] > numevent_)
0553       return;
0554   }
0555 
0556   std::vector<EBDetId> Xtals5x5 = calibCluster.get5x5Id(maxHitId);
0557 
0558   if ((int)Xtals5x5.size() != ClusterSize_ * ClusterSize_)
0559     return;
0560 
0561   // fill cluster energy
0562   std::vector<float> energy;
0563   float energy3x3 = 0.;
0564   float energy5x5 = 0.;
0565 
0566   for (int icry = 0; icry < ClusterSize_ * ClusterSize_; icry++) {
0567     EBRecHitCollection::const_iterator itrechit;
0568     if (Xtals5x5[icry].subdetId() != 1)
0569       continue;
0570     itrechit = hits->find(Xtals5x5[icry]);
0571     if (itrechit == hits->end()) {
0572       edm::LogVerbatim("EcalCalib") << "DetId not is e25";
0573       continue;
0574     }
0575 
0576     if (edm::isNotFinite(itrechit->energy()))
0577       return;
0578     energy.push_back(itrechit->energy());
0579     energy5x5 += energy[icry];
0580 
0581     if (icry == 6 || icry == 7 || icry == 8 || icry == 11 || icry == 12 || icry == 13 || icry == 16 || icry == 17 ||
0582         icry == 18) {
0583       energy3x3 += energy[icry];
0584     }
0585   }
0586   if ((int)energy.size() != ClusterSize_ * ClusterSize_)
0587     return;
0588   // Once we have the matrix 5x5, we have to correct for gaps/cracks/umbrella and maincontainement
0589 
0590   GeneralMap->Fill(maxCC_Eta, maxCC_Phi);
0591 
0592   EoP_all->Fill(highPtElectron.eSuperClusterOverP());
0593 
0594   if (highPtElectron.classification() == elecclass_ || elecclass_ == -1) {
0595     float Ptrack_in =
0596         sqrt(pow(highPtElectron.trackMomentumAtVtx().X(), 2) + pow(highPtElectron.trackMomentumAtVtx().Y(), 2) +
0597              pow(highPtElectron.trackMomentumAtVtx().Z(), 2));
0598 
0599     float UncorrectedPatCalo =
0600         sqrt(pow(highPtElectron.trackMomentumAtCalo().X(), 2) + pow(highPtElectron.trackMomentumAtCalo().Y(), 2) +
0601              pow(highPtElectron.trackMomentumAtCalo().Z(), 2));
0602 
0603     float Ptrack_out =
0604         sqrt(pow(highPtElectron.trackMomentumOut().X(), 2) + pow(highPtElectron.trackMomentumOut().Y(), 2) +
0605              pow(highPtElectron.trackMomentumOut().Z(), 2));
0606 
0607     EventMatrixNoCuts.push_back(energy);
0608     EnergyVectorNoCuts.push_back(UncorrectedPatCalo);
0609 
0610     MaxCCetaNoCuts.push_back(maxCC_Eta);
0611     MaxCCphiNoCuts.push_back(maxCC_Phi);
0612 
0613     WeightVectorNoCuts.push_back(energy5x5 / UncorrectedPatCalo);
0614 
0615     //---------------------------------------------------No Cuts-------------------------------------------------------
0616     e9NoCuts->Fill(energy3x3);
0617     e25NoCuts->Fill(energy5x5);
0618     e9Overe25NoCuts->Fill(energy3x3 / energy5x5);
0619     scENoCuts->Fill(sc.energy());
0620 
0621     trPNoCuts->Fill(UncorrectedPatCalo);
0622 
0623     EoPNoCuts->Fill(highPtElectron.eSuperClusterOverP());
0624     e25OverScENoCuts->Fill(energy5x5 / sc.energy());
0625 
0626     E25oPNoCuts->Fill(energy5x5 / UncorrectedPatCalo);
0627 
0628     MapNoCuts->Fill(maxCC_Eta, maxCC_Phi);
0629     PinOverPoutNoCuts->Fill(
0630         sqrt(pow(highPtElectron.trackMomentumAtVtx().X(), 2) + pow(highPtElectron.trackMomentumAtVtx().Y(), 2) +
0631              pow(highPtElectron.trackMomentumAtVtx().Z(), 2)) /
0632         sqrt(pow(highPtElectron.trackMomentumOut().X(), 2) + pow(highPtElectron.trackMomentumOut().Y(), 2) +
0633              pow(highPtElectron.trackMomentumOut().Z(), 2)));
0634     eSeedOverPoutNoCuts->Fill(highPtElectron.eSuperClusterOverP());
0635 
0636     MapCor1NoCuts->Fill(energy5x5 / UncorrectedPatCalo, energy5x5 / Ptrack_in);
0637     MapCor2NoCuts->Fill(energy5x5 / UncorrectedPatCalo, highPtElectron.eSuperClusterOverP());
0638     MapCor3NoCuts->Fill(energy5x5 / UncorrectedPatCalo, Ptrack_out / Ptrack_in);
0639     MapCor4NoCuts->Fill(energy5x5 / UncorrectedPatCalo, energy5x5 / highPtElectron.p());
0640     MapCor5NoCuts->Fill(energy5x5 / UncorrectedPatCalo, UncorrectedPatCalo / Ptrack_out);
0641     MapCor6NoCuts->Fill(Ptrack_out / Ptrack_in, energy5x5 / Ptrack_in);
0642     MapCor7NoCuts->Fill(Ptrack_out / Ptrack_in, UncorrectedPatCalo / Ptrack_out);
0643     MapCor8NoCuts->Fill(energy5x5 / Ptrack_in, UncorrectedPatCalo / Ptrack_out);
0644     MapCor9NoCuts->Fill(energy5x5 / UncorrectedPatCalo, highPtElectron.eSeedClusterOverPout());
0645     MapCor10NoCuts->Fill(highPtElectron.eSeedClusterOverPout(), Ptrack_out / Ptrack_in);
0646     MapCor11NoCuts->Fill(highPtElectron.eSeedClusterOverPout(), energy5x5 / Ptrack_in);
0647 
0648     PinMinPoutNoCuts->Fill((Ptrack_in - Ptrack_out) / Ptrack_in);
0649 
0650     Error1NoCuts->Fill(highPtElectron.trackMomentumError() / Ptrack_in);
0651     Error2NoCuts->Fill(highPtElectron.trackMomentumError() / Ptrack_out);
0652     Error3NoCuts->Fill(highPtElectron.trackMomentumError() / UncorrectedPatCalo);
0653     eSeedOverPout2NoCuts->Fill(highPtElectron.eSeedClusterOverPout());
0654 
0655     hadOverEmNoCuts->Fill(highPtElectron.hadronicOverEm());
0656 
0657     //------------------------------------------------Cuts-----------------------------------------------------
0658     //Cuts!
0659     if ((energy3x3 / energy5x5) < cut1_)
0660       return;
0661 
0662     if ((Ptrack_out / Ptrack_in) < cut2_ || (Ptrack_out / Ptrack_in) > cut3_)
0663       return;
0664     if ((energy5x5 / Ptrack_in) < cutEPin1_ || (energy5x5 / Ptrack_in) > cutEPin2_)
0665       return;
0666 
0667     e9->Fill(energy3x3);
0668     e25->Fill(energy5x5);
0669     e9Overe25->Fill(energy3x3 / energy5x5);
0670     scE->Fill(sc.energy());
0671     trP->Fill(UncorrectedPatCalo);
0672 
0673     EoP->Fill(highPtElectron.eSuperClusterOverP());
0674     e25OverScE->Fill(energy5x5 / sc.energy());
0675 
0676     E25oP->Fill(energy5x5 / UncorrectedPatCalo);
0677 
0678     Map->Fill(maxCC_Eta, maxCC_Phi);
0679     PinOverPout->Fill(
0680         sqrt(pow(highPtElectron.trackMomentumAtVtx().X(), 2) + pow(highPtElectron.trackMomentumAtVtx().Y(), 2) +
0681              pow(highPtElectron.trackMomentumAtVtx().Z(), 2)) /
0682         sqrt(pow(highPtElectron.trackMomentumOut().X(), 2) + pow(highPtElectron.trackMomentumOut().Y(), 2) +
0683              pow(highPtElectron.trackMomentumOut().Z(), 2)));
0684     eSeedOverPout->Fill(highPtElectron.eSuperClusterOverP());
0685 
0686     MapCor1->Fill(energy5x5 / UncorrectedPatCalo, energy5x5 / Ptrack_in);
0687     MapCor2->Fill(energy5x5 / UncorrectedPatCalo, highPtElectron.eSuperClusterOverP());
0688     MapCor3->Fill(energy5x5 / UncorrectedPatCalo, Ptrack_out / Ptrack_in);
0689     MapCor4->Fill(energy5x5 / UncorrectedPatCalo, energy5x5 / highPtElectron.p());
0690     MapCor5->Fill(energy5x5 / UncorrectedPatCalo, UncorrectedPatCalo / Ptrack_out);
0691     MapCor6->Fill(Ptrack_out / Ptrack_in, energy5x5 / Ptrack_in);
0692     MapCor7->Fill(Ptrack_out / Ptrack_in, UncorrectedPatCalo / Ptrack_out);
0693     MapCor8->Fill(energy5x5 / Ptrack_in, UncorrectedPatCalo / Ptrack_out);
0694     MapCor9->Fill(energy5x5 / UncorrectedPatCalo, highPtElectron.eSeedClusterOverPout());
0695     MapCor10->Fill(highPtElectron.eSeedClusterOverPout(), Ptrack_out / Ptrack_in);
0696     MapCor11->Fill(highPtElectron.eSeedClusterOverPout(), energy5x5 / Ptrack_in);
0697 
0698     PinMinPout->Fill((Ptrack_in - Ptrack_out) / Ptrack_in);
0699 
0700     Error1->Fill(highPtElectron.trackMomentumError() / Ptrack_in);
0701     Error2->Fill(highPtElectron.trackMomentumError() / Ptrack_out);
0702     Error3->Fill(highPtElectron.trackMomentumError() / UncorrectedPatCalo);
0703 
0704     eSeedOverPout2->Fill(highPtElectron.eSeedClusterOverPout());
0705     hadOverEm->Fill(highPtElectron.hadronicOverEm());
0706 
0707     EventMatrix.push_back(energy);
0708     EnergyVector.push_back(UncorrectedPatCalo);
0709     MaxCCeta.push_back(maxCC_Eta);
0710     MaxCCphi.push_back(maxCC_Phi);
0711 
0712     WeightVector.push_back(energy5x5 / UncorrectedPatCalo);
0713 
0714     //-------------------------------------------------------Extra Cut-----------------------------------------------------
0715     if (highPtElectron.eSeedClusterOverPout() < cutESeed_)
0716       return;
0717 
0718     MapCor1ESeed->Fill(energy5x5 / UncorrectedPatCalo, energy5x5 / Ptrack_in);
0719     MapCor2ESeed->Fill(energy5x5 / UncorrectedPatCalo, highPtElectron.eSuperClusterOverP());
0720     MapCor3ESeed->Fill(energy5x5 / UncorrectedPatCalo, Ptrack_out / Ptrack_in);
0721     MapCor4ESeed->Fill(energy5x5 / UncorrectedPatCalo, energy5x5 / highPtElectron.p());
0722     MapCor5ESeed->Fill(energy5x5 / UncorrectedPatCalo, UncorrectedPatCalo / Ptrack_out);
0723     MapCor6ESeed->Fill(Ptrack_out / Ptrack_in, energy5x5 / Ptrack_in);
0724     MapCor7ESeed->Fill(Ptrack_out / Ptrack_in, UncorrectedPatCalo / Ptrack_out);
0725     MapCor8ESeed->Fill(energy5x5 / Ptrack_in, UncorrectedPatCalo / Ptrack_out);
0726     MapCor9ESeed->Fill(energy5x5 / UncorrectedPatCalo, highPtElectron.eSeedClusterOverPout());
0727     MapCor10ESeed->Fill(highPtElectron.eSeedClusterOverPout(), Ptrack_out / Ptrack_in);
0728     MapCor11ESeed->Fill(highPtElectron.eSeedClusterOverPout(), energy5x5 / Ptrack_in);
0729 
0730     eSeedOverPout2ESeed->Fill(highPtElectron.eSeedClusterOverPout());
0731 
0732     hadOverEmESeed->Fill(highPtElectron.hadronicOverEm());
0733 
0734   } else {
0735     return;
0736   }
0737 }