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
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
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
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
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
0241 ReducedMap = calibCluster.getMap(etaMin, etaMax, phiMin, phiMax);
0242
0243 oldCalibs.resize(ReducedMap.size(), 0.);
0244
0245
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
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
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();
0465
0466
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
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
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
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
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
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
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
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
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
0658
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
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 }
