Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 12:47:58

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