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