CMS 3D CMS Logo
ElectronCalibrationUniv.cc
Go to the documentation of this file.
1 
2 #include "FWCore/ParameterSet/interface/ParameterSet.h"
3 
4 #include "DataFormats/EcalDetId/interface/EBDetId.h"
5 #include "DataFormats/DetId/interface/DetId.h"
6 #include "Calibration/Tools/interface/calibXMLwriter.h"
7 #include "Calibration/Tools/interface/CalibrationCluster.h"
8 #include "Calibration/Tools/interface/HouseholderDecomposition.h"
9 #include "Calibration/Tools/interface/MinL3Algorithm.h"
10 #include "Calibration/EcalCalibAlgos/interface/ElectronCalibrationUniv.h"
11 #include "FWCore/Utilities/interface/isFinite.h"
12 #include "TF1.h"
13 #include "TRandom.h"
14 
15 #include <iostream>
16 #include <string>
17 #include <stdexcept>
18 #include <vector>
19 
20 ElectronCalibrationUniv::ElectronCalibrationUniv(const edm::ParameterSet& iConfig)
21  : ebRecHitLabel_(iConfig.getParameter<edm::InputTag>("ebRecHitsLabel")),
22  eeRecHitLabel_(iConfig.getParameter<edm::InputTag>("eeRecHitsLabel")),
23  electronLabel_(iConfig.getParameter<edm::InputTag>("electronLabel")),
24  rootfile_(iConfig.getParameter<std::string>("rootfile")),
25  calibAlgo_(iConfig.getParameter<std::string>("CALIBRATION_ALGO")),
26  miscalibfile_(iConfig.getParameter<std::string>("miscalibfile")),
27  miscalibfileEndCap_(iConfig.getParameter<std::string>("miscalibfileEndCap")),
28  keventweight_(iConfig.getParameter<int>("keventweight")),
29  elePt_(iConfig.getParameter<double>("ElePt")),
30  maxeta_(iConfig.getParameter<int>("maxeta")),
31  mineta_(iConfig.getParameter<int>("mineta")),
32  maxphi_(iConfig.getParameter<int>("maxphi")),
33  minphi_(iConfig.getParameter<int>("minphi")),
34  cut1_(iConfig.getParameter<double>("cut1")),
35  cut2_(iConfig.getParameter<double>("cut2")),
36  cut3_(iConfig.getParameter<double>("cut3")),
37  numevent_(iConfig.getParameter<int>("numevent")),
38  cutEPCalo1_(iConfig.getParameter<double>("cutEPCaloMin")),
39  cutEPCalo2_(iConfig.getParameter<double>("cutEPCaloMax")),
40  cutEPin1_(iConfig.getParameter<double>("cutEPinMin")),
41  cutEPin2_(iConfig.getParameter<double>("cutEPinMax")),
42  cutCalo1_(iConfig.getParameter<double>("cutCaloMin")),
43  cutCalo2_(iConfig.getParameter<double>("cutCaloMax")),
44  cutESeed_(iConfig.getParameter<double>("cutESeed")),
45  clusterSize_(iConfig.getParameter<int>("Clustersize")),
46  elecclass_(iConfig.getParameter<int>("elecclass")),
47  ebRecHitToken_(consumes<EBRecHitCollection>(ebRecHitLabel_)),
48  eeRecHitToken_(consumes<EERecHitCollection>(eeRecHitLabel_)),
49  gsfElectronToken_(consumes<reco::GsfElectronCollection>(electronLabel_)),
50  topologyToken_(esConsumes<edm::Transition::BeginRun>()) {}
51 
52 ElectronCalibrationUniv::~ElectronCalibrationUniv() {}
53 
54 //========================================================================
55 void ElectronCalibrationUniv::beginJob() {
56  //========================================================================
57  f = new TFile(rootfile_.c_str(), "RECREATE");
58  f->cd();
59  EventsAfterCuts = new TH1F("EventsAfterCuts", "Events After Cuts", 30, 0, 30);
60 
61  // Book histograms
62  e9 = new TH1F("e9", "E9 energy", 300, 0., 150.);
63  e25 = new TH1F("e25", "E25 energy", 300, 0., 150.);
64  scE = new TH1F("scE", "SC energy", 300, 0., 150.);
65  trP = new TH1F("trP", "Trk momentum", 300, 0., 150.);
66  EoP = new TH1F("EoP", "EoP", 600, 0., 3.);
67  EoP_all = new TH1F("EoP_all", "EoP_all", 600, 0., 3.);
68 
69  calibs = new TH1F("calib", "Calibration constants", 800, 0.5, 2.);
70  calibsEndCapMinus = new TH1F("calibEndCapMinus", "Calibration constants EE-", 800, 0.5, 2.);
71  calibsEndCapPlus = new TH1F("calibEndCapPlus", "Calibration constants EE+", 800, 0.5, 2.);
72 
73  e25OverScE = new TH1F("e25OverscE", "E25 / SC energy", 400, 0., 2.);
74  E25oP = new TH1F("E25oP", "E25 / P", 750, 0., 1.5);
75 
76  Map = new TH2F("Map", "Nb Events in Crystal", 173, -86, 86, 362, 0, 361);
77  e9Overe25 = new TH1F("e9Overe25", "E9 / E25", 400, 0., 2.);
78  Map3Dcalib = new TH2F("3Dcalib", "3Dcalib", 173, -86, 86, 362, 0, 361);
79  Map3DcalibEndCapMinus = new TH2F("3DcalibEndCapMinus", "3Dcalib EE-", 100, 0, 100, 100, 0, 100);
80  Map3DcalibEndCapPlus = new TH2F("3DcalibEndCapPlus", "3Dcalib EE+", 100, 0, 100, 100, 0, 100);
81 
82  MapCor1 = new TH2F("MapCor1", "Correlation E25/Pcalo versus E25/Pin", 100, 0., 5., 100, 0., 5.);
83  MapCor2 = new TH2F("MapCor2", "Correlation E25/Pcalo versus E/P", 100, 0., 5., 100, 0., 5.);
84  MapCor3 = new TH2F("MapCor3", "Correlation E25/Pcalo versus Pout/Pin", 100, 0., 5., 100, 0., 5.);
85  MapCor4 = new TH2F("MapCor4", "Correlation E25/Pcalo versus E25/highestP", 100, 0., 5., 100, 0., 5.);
86  MapCor5 = new TH2F("MapCor5", "Correlation E25/Pcalo versus Pcalo/Pout", 100, 0., 5., 100, 0., 5.);
87  MapCor6 = new TH2F("MapCor6", "Correlation Pout/Pin versus E25/Pin", 100, 0., 5., 100, 0., 5.);
88  MapCor7 = new TH2F("MapCor7", "Correlation Pout/Pin versus Pcalo/Pout", 100, 0., 5., 100, 0., 5.);
89  MapCor8 = new TH2F("MapCor8", "Correlation E25/Pin versus Pcalo/Pout", 100, 0., 5., 100, 0., 5.);
90  MapCor9 = new TH2F("MapCor9", "Correlation E25/Pcalo versus Eseed/Pout", 100, 0., 5., 100, 0., 5.);
91  MapCor10 = new TH2F("MapCor10", "Correlation Eseed/Pout versus Pout/Pin", 100, 0., 5., 100, 0., 5.);
92  MapCor11 = new TH2F("MapCor11", "Correlation Eseed/Pout versus E25/Pin", 100, 0., 5., 100, 0., 5.);
93  // MapCorCalib = new TH2F ("MapCorCalib", "Correlation Miscalibration versus Calibration constants", 500, 0.5,1.5, 500, 0.5, 1.5);
94 
95  E25oPvsEta = new TH2F("E25oPvsEta", "E/P vs Eta", 173, -86, 86, 600, 0.7, 1.3);
96  E25oPvsEtaEndCapMinus = new TH2F("E25oPvsEtaEndCapMinus", "E/P vs R EE-", 100, 0, 100, 600, 0.7, 1.3);
97  E25oPvsEtaEndCapPlus = new TH2F("E25oPvsEtaEndCapPlus", "E/P vs R EE+", 100, 0, 100, 600, 0.7, 1.3);
98 
99  PinMinPout = new TH1F("PinMinPout", "(Pin - Pout)/Pin", 600, -2.0, 2.0);
100 
101  calibinter = new TH1F("calibinter", "internal calibration constants", 800, 0.5, 2.);
102  PinOverPout = new TH1F("PinOverPout", "pinOverpout", 600, 0., 3.);
103  eSeedOverPout = new TH1F("eSeedOverPout", "eSeedOverpout ", 600, 0., 3.);
104  // MisCalibs = new TH1F("MisCalibs","Miscalibration constants",800,0.5,2.);
105  // RatioCalibs = new TH1F("RatioCalibs","Ratio in Calibration Constants", 800, 0.5, 2.0);
106  // DiffCalibs = new TH1F("DiffCalibs", "Difference in Calibration constants", 800, -1.0,1.0);
107  calibinterEndCapMinus = new TH1F("calibinterEndCapMinus", "internal calibration constants", 800, 0.5, 2.);
108  calibinterEndCapPlus = new TH1F("calibinterEndCapPlus", "internal calibration constants", 800, 0.5, 2.);
109  // MisCalibsEndCapMinus = new TH1F("MisCalibsEndCapMinus","Miscalibration constants",800,0.5,2.);
110  // MisCalibsEndCapPlus = new TH1F("MisCalibsEndCapPlus","Miscalibration constants",800,0.5,2.);
111  // RatioCalibsEndCapMinus = new TH1F("RatioCalibsEndCapMinus","Ratio in Calibration Constants", 800, 0.5, 2.0);
112  // RatioCalibsEndCapPlus = new TH1F("RatioCalibsEndCapPlus","Ratio in Calibration Constants", 800, 0.5, 2.0);
113  // DiffCalibsEndCapMinus = new TH1F("DiffCalibsEndCapMinus", "Difference in Calibration constants", 800, -1.0,1.0);
114  // DiffCalibsEndCapPlus = new TH1F("DiffCalibsEndCapPlus", "Difference in Calibration constants", 800, -1.0,1.0);
115  Error1 = new TH1F("Error1", "DeltaP/Pin", 800, -1.0, 1.0);
116  Error2 = new TH1F("Error2", "DeltaP/Pout", 800, -1.0, 1.0);
117  Error3 = new TH1F("Error3", "DeltaP/Pcalo", 800, -1.0, 1.0);
118  eSeedOverPout2 = new TH1F("eSeedOverPout2", "eSeedOverpout (No Supercluster)", 600, 0., 4.);
119  hadOverEm = new TH1F("hadOverEm", "Had/EM distribution", 600, -2., 2.);
120 
121  // Book histograms
122  Map3DcalibNoCuts = new TH2F("3DcalibNoCuts", "3Dcalib (Before Cuts)", 173, -86, 86, 362, 0, 361);
123  e9NoCuts = new TH1F("e9NoCuts", "E9 energy (Before Cuts)", 300, 0., 150.);
124  e25NoCuts = new TH1F("e25NoCuts", "E25 energy (Before Cuts)", 300, 0., 150.);
125  scENoCuts = new TH1F("scENoCuts", "SC energy (Before Cuts)", 300, 0., 150.);
126  trPNoCuts = new TH1F("trPNoCuts", "Trk momentum (Before Cuts)", 300, 0., 150.);
127  EoPNoCuts = new TH1F("EoPNoCuts", "EoP (Before Cuts)", 600, 0., 3.);
128  calibsNoCuts = new TH1F("calibNoCuts", "Calibration constants (Before Cuts)", 800, 0., 2.);
129  e25OverScENoCuts = new TH1F("e25OverscENoCuts", "E25 / SC energy (Before Cuts)", 400, 0., 2.);
130  E25oPNoCuts = new TH1F("E25oPNoCuts", "E25 / P (Before Cuts)", 750, 0., 1.5);
131  MapEndCapMinus = new TH2F("MapEndCapMinus", "Nb Events in Crystal (EndCap)", 100, 0, 100, 100, 0, 100);
132  MapEndCapPlus = new TH2F("MapEndCapPlus", "Nb Events in Crystal (EndCap)", 100, 0, 100, 100, 0, 100);
133  e9Overe25NoCuts = new TH1F("e9Overe25NoCuts", "E9 / E25 (Before Cuts)", 400, 0., 2.);
134  PinOverPoutNoCuts = new TH1F("PinOverPoutNoCuts", "pinOverpout (Before Cuts)", 600, 0., 3.);
135  eSeedOverPoutNoCuts = new TH1F(" eSeedOverPoutNoCuts", "eSeedOverpout (Before Cuts) ", 600, 0., 4.);
136  PinMinPoutNoCuts = new TH1F("PinMinPoutNoCuts", "(Pin - Pout)/Pin (Before Cuts)", 600, -2.0, 2.0);
137 
138  // RatioCalibsNoCuts = new TH1F("RatioCalibsNoCuts","Ratio in Calibration Constants (Before Cuts)", 800, 0.5, 2.0);
139  // DiffCalibsNoCuts = new TH1F("DiffCalibsNoCuts", "Difference in Calibration constants (Before Cuts)", 800, -1.0,1.0);
140  calibinterNoCuts = new TH1F("calibinterNoCuts", "internal calibration constants", 2000, 0.5, 2.);
141 
142  MapCor1NoCuts =
143  new TH2F("MapCor1NoCuts", "Correlation E25/PatCalo versus E25/Pin (Before Cuts)", 100, 0., 5., 100, 0., 5.);
144  MapCor2NoCuts =
145  new TH2F("MapCor2NoCuts", "Correlation E25/PatCalo versus E/P (Before Cuts)", 100, 0., 5., 100, 0., 5.);
146  MapCor3NoCuts =
147  new TH2F("MapCor3NoCuts", "Correlation E25/PatCalo versus Pout/Pin (Before Cuts)", 100, 0., 5., 100, 0., 5.);
148  MapCor4NoCuts =
149  new TH2F("MapCor4NoCuts", "Correlation E25/PatCalo versus E25/highestP (Before Cuts)", 100, 0., 5., 100, 0., 5.);
150  MapCor5NoCuts =
151  new TH2F("MapCor5NoCuts", "Correlation E25/Pcalo versus Pcalo/Pout (Before Cuts)", 100, 0., 5., 100, 0., 5.);
152  MapCor6NoCuts =
153  new TH2F("MapCor6NoCuts", "Correlation Pout/Pin versus E25/Pin (Before Cuts)", 100, 0., 5., 100, 0., 5.);
154  MapCor7NoCuts =
155  new TH2F("MapCor7NoCuts", "Correlation Pout/Pin versus Pcalo/Pout (Before Cuts)", 100, 0., 5., 100, 0., 5.);
156  MapCor8NoCuts =
157  new TH2F("MapCor8NoCuts", "Correlation E25/Pin versus Pcalo/Pout (Before Cuts)", 100, 0., 5., 100, 0., 5.);
158  MapCor9NoCuts =
159  new TH2F("MapCor9NoCuts", "Correlation E25/Pcalo versus Eseed/Pout (Before Cuts)", 100, 0., 5., 100, 0., 5.);
160  MapCor10NoCuts =
161  new TH2F("MapCor10NoCuts", "Correlation Eseed/Pout versus Pout/Pin (Before Cuts)", 100, 0., 5., 100, 0., 5.);
162  MapCor11NoCuts =
163  new TH2F("MapCor11NoCuts", "Correlation Eseed/Pout versus E25/Pin (Before Cuts)", 100, 0., 5., 100, 0., 5.);
164  // MapCorCalibEndCapMinus = new TH2F ("MapCorCalibEndCapMinus", "Correlation Miscalibration versus Calibration constants (EndCap)", 500, 0.5,1.5, 500, 0.5, 1.5);
165  // MapCorCalibEndCapPlus = new TH2F ("MapCorCalibEndCapPlus", "Correlation Miscalibration versus Calibration constants (EndCap)", 500, 0.5,1.5, 500, 0.5, 1.5);
166 
167  Error1NoCuts = new TH1F("Eror1NoCuts", "DeltaP/Pin (Before Cuts)", 800, -1.0, 1.0);
168  Error2NoCuts = new TH1F("Error2NoCuts", "DeltaP/Pout (Before Cuts)", 800, -1.0, 1.0);
169  Error3NoCuts = new TH1F("Error3NoCuts", "DeltaP/Pcalo (Before Cuts)", 800, -1.0, 1.0);
170  eSeedOverPout2NoCuts = new TH1F("eSeedOverPout2NoCuts", "eSeedOverpout (No Supercluster, Before Cuts)", 600, 0., 4.);
171  hadOverEmNoCuts = new TH1F("hadOverEmNoCuts", "Had/EM distribution (Before Cuts)", 600, -2., 2.);
172 
173  //Book histograms after ESeed cut
174  MapCor1ESeed =
175  new TH2F("MapCor1ESeed", "Correlation E25/Pcalo versus E25/Pin (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
176  MapCor2ESeed =
177  new TH2F("MapCor2ESeed", "Correlation E25/Pcalo versus E/P (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
178  MapCor3ESeed = new TH2F(
179  "MapCor3ESeed", "Correlation E25/Pcalo versus Pout/Pin (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
180  MapCor4ESeed = new TH2F(
181  "MapCor4ESeed", "Correlation E25/Pcalo versus E25/highestP (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
182  MapCor5ESeed = new TH2F(
183  "MapCor5ESeed", "Correlation E25/Pcalo versus Pcalo/Pout (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
184  MapCor6ESeed =
185  new TH2F("MapCor6ESeed", "Correlation Pout/Pin versus E25/Pin (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
186  MapCor7ESeed = new TH2F(
187  "MapCor7ESeed", "Correlation Pout/Pin versus Pcalo/Pout (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
188  MapCor8ESeed = new TH2F(
189  "MapCor8ESeed", "Correlation E25/Pin versus Pcalo/Pout (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
190  MapCor9ESeed = new TH2F(
191  "MapCor9ESeed", "Correlation E25/Pcalo versus Eseed/Pout (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
192  MapCor10ESeed = new TH2F(
193  "MapCor10ESeed", "Correlation Eseed/Pout versus Pout/Pin (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
194  MapCor11ESeed = new TH2F(
195  "MapCor11ESeed", "Correlation Eseed/Pout versus E25/Pin (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
196 
197  eSeedOverPout2ESeed =
198  new TH1F("eSeedOverPout2ESeed", "eSeedOverpout (No Supercluster, after Eseed/Pout cut)", 600, 0., 4.);
199 
200  hadOverEmESeed = new TH1F("hadOverEmESeed", "Had/EM distribution (after Eseed/Pout cut)", 600, -2., 2.);
201 
202  //Book histograms without any cut
203  GeneralMap = new TH2F("GeneralMap", "Map without any cuts", 173, -86, 86, 362, 0, 361);
204  GeneralMapEndCapMinus = new TH2F("GeneralMapEndCapMinus", "Map without any cuts", 100, 0, 100, 100, 0, 100);
205  GeneralMapEndCapPlus = new TH2F("GeneralMapEndCapPlus", "Map without any cuts", 100, 0, 100, 100, 0, 100);
206  GeneralMapBeforePt = new TH2F("GeneralMapBeforePt", "Map without any cuts", 173, -86, 86, 362, 0, 361);
207  GeneralMapEndCapMinusBeforePt =
208  new TH2F("GeneralMapEndCapMinusBeforePt", "Map without any cuts", 100, 0, 100, 100, 0, 100);
209  GeneralMapEndCapPlusBeforePt =
210  new TH2F("GeneralMapEndCapPlusBeforePt", "Map without any cuts", 100, 0, 100, 100, 0, 100);
211 
212  calibClusterSize = clusterSize_;
213  etaMin = int(mineta_);
214  etaMax = int(maxeta_);
215  phiMin = int(minphi_);
216  phiMax = int(maxphi_);
217  if (calibAlgo_ == "L3") {
218  MyL3Algo1 = new MinL3Algorithm(keventweight_, calibClusterSize, etaMin, etaMax, phiMin, phiMax);
219  } else {
220  if (calibAlgo_ == "L3Univ") {
221  UnivL3 = new MinL3AlgoUniv<DetId>(keventweight_);
222  } else {
223  if (calibAlgo_ == "HH" || calibAlgo_ == "HHReg") {
224  MyHH = new HouseholderDecomposition(calibClusterSize, etaMin, etaMax, phiMin, phiMax);
225  } else {
226  std::cout << " Name of Algorithm is not recognize " << calibAlgo_
227  << " Should be either L3, HH or HHReg. Abort! " << std::endl;
228  }
229  }
230  }
231  read_events = 0;
232 }
233 
234 //========================================================================
235 void ElectronCalibrationUniv::beginRun(edm::Run const&, edm::EventSetup const& iSetup) {
236  //========================================================================
237 
238  //To Deal with Geometry:
239  theCaloTopology_ = &iSetup.getData(topologyToken_);
240 }
241 
242 //========================================================================
243 void ElectronCalibrationUniv::endRun(edm::Run const&, edm::EventSetup const& iSetup) {}
244 //========================================================================
245 
246 //========================================================================
247 
248 void ElectronCalibrationUniv::endJob() {
249  //========================================================================
250 
251  f->cd();
252  time_t start, end;
253  time_t cpu_time_used;
254  start = time(nullptr);
255 
256  //In order to do only one loop to use properly looper properties, ask only for 1 iterations!
257  int nIterations = 10;
258  if (calibAlgo_ == "L3") {
259  solution = MyL3Algo1->iterate(EventMatrix, MaxCCeta, MaxCCphi, EnergyVector, nIterations);
260  } else {
261  if (calibAlgo_ == "L3Univ") {
262  //Univsolution= UnivL3->getSolution();
263  // std::cout<<" Should derive solution "<<EnergyVector.size()<<std::endl;
264  Univsolution = UnivL3->iterate(EventMatrix, UnivEventIds, EnergyVector, nIterations);
265  //std::cout<<" solution size "<<Univsolution.size()<<std::endl;
266  } else {
267  if (calibAlgo_ == "HH") {
268  solution = MyHH->iterate(EventMatrix, MaxCCeta, MaxCCphi, EnergyVector, 1, false);
269  } else {
270  if (calibAlgo_ == "HHReg") {
271  solution = MyHH->runRegional(EventMatrix, MaxCCeta, MaxCCphi, EnergyVector, 2);
272  } else {
273  std::cout << " Calibration not run due to problem in Algo Choice..." << std::endl;
274  return;
275  }
276  }
277  }
278  }
279  end = time(nullptr);
280  cpu_time_used = end - start;
281  // std::cout<<"222 solution size "<<Univsolution.size()<<std::endl;
282 
283  calibXMLwriter write_calibrations;
284 
285  // FILE* MisCalib;
286  // //char* calibfile="miscalibfile";
287  // MisCalib = fopen(miscalibfile_.c_str(),"r");
288 
289  // int fileStatus=0;
290  // int eta=-1;
291  // int phi=-1;
292  // float coeff=-1;
293 
294  std::map<EBDetId, float> OldCoeff;
295 
296  // while(fileStatus != EOF) {
297  // fileStatus = fscanf(MisCalib,"%d %d %f\n", &eta,&phi,&coeff);
298  // if(eta!=-1&&phi!=-1&& coeff!=-1){
299  // // std::cout<<" We have read correctly the coefficient " << coeff << " corresponding to eta "<<eta<<" and phi "<<phi<<std::endl;
300  // OldCoeff.insert(std::make_pair(EBDetId(eta,phi,EBDetId::ETAPHIMODE),coeff ));
301  // }
302  // }
303 
304  // fclose(MisCalib);
305  // FILE* MisCalibEndCap;
306  // //char* calibfile="miscalibfile";
307  // MisCalibEndCap = fopen(miscalibfileEndCap_.c_str(),"r");
308 
309  // int fileStatus2=0;
310  // int X=-1;
311  // int Y=-1;
312  // float coeff2=-1;
313  std::map<EEDetId, float> OldCoeffEndCap;
314 
315  // while(fileStatus2 != EOF) {
316  // fileStatus2 = fscanf(MisCalibEndCap,"%d %d %f\n", &X,&Y,&coeff2);
317  // if(X!=-1&&Y!=-1&& coeff2!=-1){
318  // // std::cout<<" We have read correctly the coefficient " << coeff << " corresponding to eta "<<eta<<" and phi "<<phi<<std::endl;
319  // if(TestEEvalidDetId(X,Y,1)){
320  // OldCoeffEndCap.insert(std::make_pair(EEDetId(X,Y,1,EEDetId::XYMODE),coeff2 ));
321  // }
322  // }
323  // }
324 
325  // fclose(MisCalibEndCap);
326  std::map<DetId, float>::const_iterator itmap;
327  for (itmap = Univsolution.begin(); itmap != Univsolution.end(); itmap++) {
328  const DetId Id(itmap->first);
329  if (Id.subdetId() == 1) {
330  const EBDetId IChannelDetId(itmap->first);
331  if (IChannelDetId.ieta() < mineta_) {
332  continue;
333  }
334  if (IChannelDetId.ieta() > maxeta_) {
335  continue;
336  }
337  if (IChannelDetId.iphi() < minphi_) {
338  continue;
339  }
340  if (IChannelDetId.iphi() > maxphi_) {
341  continue;
342  }
343  // float Compare=1;
344  // std::map<EBDetId,float>::iterator iter = OldCoeff.find(itmap->first);
345  // if( iter != OldCoeff.end() )Compare = iter->second;
346  Map3Dcalib->Fill(IChannelDetId.ieta(), IChannelDetId.iphi(), itmap->second);
347  calibs->Fill(itmap->second);
348  //DiffCalibs->Fill(newCalibs[icry]-miscalib[IChannelDetId.ieta()-1][IChannelDetId.iphi()-21]);
349  //RatioCalibs->Fill(newCalibs[icry]/miscalib[IChannelDetId.ieta()-1][IChannelDetId.iphi()-21]);
350  if (IChannelDetId.ieta() < mineta_ + 2) {
351  continue;
352  }
353  if (IChannelDetId.ieta() > maxeta_ - 2) {
354  continue;
355  }
356  if (IChannelDetId.iphi() < minphi_ + 2) {
357  continue;
358  }
359  if (IChannelDetId.iphi() > maxphi_ - 2) {
360  continue;
361  }
362  write_calibrations.writeLine(IChannelDetId, itmap->second);
363  calibinter->Fill(itmap->second);
364  // MapCorCalib->Fill(itmap->second,Compare);
365  // DiffCalibs->Fill(itmap->second-Compare);
366  // RatioCalibs->Fill(itmap->second*Compare);
367  } else {
368  const EEDetId IChannelDetId(itmap->first);
369  // if (IChannelDetId.ix()<0 ){continue;}
370  // if (IChannelDetId.ix()>100 ){continue;}
371  // if (IChannelDetId.iy()<0 ){continue;}
372  // if (IChannelDetId.iy()>100 ){continue;}
373  // std::map<EEDetId,float>::iterator iter = OldCoeffEndCap.find(itmap->first);
374  // float Compare=1;
375  // if( iter != OldCoeffEndCap.end() )Compare = iter->second;
376  if (IChannelDetId.zside() < 0) {
377  Map3DcalibEndCapMinus->Fill(IChannelDetId.ix(), IChannelDetId.iy(), itmap->second);
378  calibsEndCapMinus->Fill(itmap->second);
379  calibinterEndCapMinus->Fill(itmap->second);
380  // DiffCalibsEndCapMinus->Fill(itmap->second-Compare);
381  // RatioCalibsEndCapMinus->Fill(itmap->second*Compare);
382  // MapCorCalibEndCapMinus->Fill(itmap->second,Compare);
383  } else {
384  Map3DcalibEndCapPlus->Fill(IChannelDetId.ix(), IChannelDetId.iy(), itmap->second);
385  calibsEndCapPlus->Fill(itmap->second);
386  calibinterEndCapPlus->Fill(itmap->second);
387  // DiffCalibsEndCapPlus->Fill(itmap->second-Compare);
388  // RatioCalibsEndCapPlus->Fill(itmap->second*Compare);
389  // MapCorCalibEndCapPlus->Fill(itmap->second,Compare);
390  }
391  write_calibrations.writeLine(IChannelDetId, itmap->second);
392  }
393  }
394  EventsAfterCuts->Write();
395 
396  // Book histograms
397  e25->Write();
398  e9->Write();
399  scE->Write();
400  trP->Write();
401  EoP->Write();
402  EoP_all->Write();
403  calibs->Write();
404  calibsEndCapMinus->Write();
405  calibsEndCapPlus->Write();
406  e9Overe25->Write();
407  e25OverScE->Write();
408  Map->Write();
409  E25oP->Write();
410 
411  PinOverPout->Write();
412  eSeedOverPout->Write();
413  // MisCalibs->Write();
414  // RatioCalibs->Write();
415  // DiffCalibs->Write();
416  // RatioCalibsNoCuts->Write();
417  // DiffCalibsNoCuts->Write();
418  // MisCalibsEndCapMinus->Write();
419  // MisCalibsEndCapPlus->Write();
420  // RatioCalibsEndCapMinus->Write();
421  // RatioCalibsEndCapPlus->Write();
422  // DiffCalibsEndCapMinus->Write();
423  // DiffCalibsEndCapPlus->Write();
424 
425  e25NoCuts->Write();
426  e9NoCuts->Write();
427  scENoCuts->Write();
428  trPNoCuts->Write();
429  EoPNoCuts->Write();
430  calibsNoCuts->Write();
431  e9Overe25NoCuts->Write();
432  e25OverScENoCuts->Write();
433  MapEndCapMinus->Write();
434  MapEndCapPlus->Write();
435  E25oPNoCuts->Write();
436  Map3Dcalib->Write();
437  Map3DcalibEndCapMinus->Write();
438  Map3DcalibEndCapPlus->Write();
439  Map3DcalibNoCuts->Write();
440  calibinter->Write();
441  calibinterEndCapMinus->Write();
442  calibinterEndCapPlus->Write();
443  calibinterNoCuts->Write();
444  PinOverPoutNoCuts->Write();
445  eSeedOverPoutNoCuts->Write();
446 
447  GeneralMap->Write();
448  GeneralMapEndCapMinus->Write();
449  GeneralMapEndCapPlus->Write();
450  GeneralMapBeforePt->Write();
451  GeneralMapEndCapMinusBeforePt->Write();
452  GeneralMapEndCapPlusBeforePt->Write();
453 
454  MapCor1->Write();
455  MapCor2->Write();
456  MapCor3->Write();
457  MapCor4->Write();
458  MapCor5->Write();
459  MapCor6->Write();
460  MapCor7->Write();
461  MapCor8->Write();
462  MapCor9->Write();
463  MapCor10->Write();
464  MapCor11->Write();
465  // MapCorCalib->Write();
466 
467  MapCor1NoCuts->Write();
468  MapCor2NoCuts->Write();
469  MapCor3NoCuts->Write();
470  MapCor4NoCuts->Write();
471  MapCor5NoCuts->Write();
472  MapCor6NoCuts->Write();
473  MapCor7NoCuts->Write();
474  MapCor8NoCuts->Write();
475  MapCor9NoCuts->Write();
476  MapCor10NoCuts->Write();
477  MapCor11NoCuts->Write();
478  // MapCorCalibEndCapMinus->Write();
479  // MapCorCalibEndCapPlus->Write();
480 
481  MapCor1ESeed->Write();
482  MapCor2ESeed->Write();
483  MapCor3ESeed->Write();
484  MapCor4ESeed->Write();
485  MapCor5ESeed->Write();
486  MapCor6ESeed->Write();
487  MapCor7ESeed->Write();
488  MapCor8ESeed->Write();
489  MapCor9ESeed->Write();
490  MapCor10ESeed->Write();
491  MapCor11ESeed->Write();
492 
493  E25oPvsEta->Write();
494  E25oPvsEtaEndCapMinus->Write();
495  E25oPvsEtaEndCapPlus->Write();
496 
497  PinMinPout->Write();
498  PinMinPoutNoCuts->Write();
499 
500  Error1->Write();
501  Error2->Write();
502  Error3->Write();
503  Error1NoCuts->Write();
504  Error2NoCuts->Write();
505  Error3NoCuts->Write();
506 
507  eSeedOverPout2->Write();
508  eSeedOverPout2NoCuts->Write();
509  eSeedOverPout2ESeed->Write();
510 
511  hadOverEm->Write();
512  hadOverEmNoCuts->Write();
513  hadOverEmESeed->Write();
514 
515  f->Write();
516 
517  f->Close();
518  // if(MyL3Algo1)delete MyL3Algo1;
519  // if(UnivL3)delete UnivL3;
520  // if(MyHH)delete MyHH;
521  // delete f;
523 
524  std::cout << " " << std::endl;
525  std::cout << "************* STATISTICS **************" << std::endl;
526  std::cout << " Events Studied " << read_events << std::endl;
527  std::cout << "Timing info:" << std::endl;
528  std::cout << "CPU time usage -- calibrating: " << cpu_time_used << " sec." << std::endl;
529 
531 }
532 
533 DetId ElectronCalibrationUniv::findMaxHit(const std::vector<DetId>& v1,
534  const EBRecHitCollection* EBhits,
535  const EERecHitCollection* EEhits) {
536  //=================================================================================
537 
538  double currEnergy = 0.;
539  DetId maxHit;
540 
541  for (std::vector<DetId>::const_iterator idsIt = v1.begin(); idsIt != v1.end(); ++idsIt) {
542  if (idsIt->subdetId() == 1) {
543  EBRecHitCollection::const_iterator itrechit;
544  itrechit = EBhits->find(*idsIt);
545  if (itrechit == EBhits->end()) {
546  std::cout << "ElectronCalibration::findMaxHit: rechit not found! " << (EBDetId)(*idsIt) << std::endl;
547  continue;
548  }
549  if (itrechit->energy() > currEnergy) {
550  currEnergy = itrechit->energy();
551  maxHit = *idsIt;
552  }
553  } else {
554  EERecHitCollection::const_iterator itrechit;
555  itrechit = EEhits->find(*idsIt);
556  if (itrechit == EEhits->end()) {
557  std::cout << "ElectronCalibration::findMaxHit: rechit not found! idsIt = " << (EEDetId)(*idsIt) << std::endl;
558  continue;
559  }
560 
561  if (itrechit->energy() > currEnergy) {
562  currEnergy = itrechit->energy();
563  maxHit = *idsIt;
564  }
565  }
566  }
567 
568  return maxHit;
569 }
570 
571 bool ElectronCalibrationUniv::TestEEvalidDetId(int crystal_ix, int crystal_iy, int iz) {
572  bool valid = false;
573  if (crystal_ix < 1 || crystal_ix > 100 || crystal_iy < 1 || crystal_iy > 100 || abs(iz) != 1) {
574  return valid;
575  }
576  if ((crystal_ix >= 1 && crystal_ix <= 3 && (crystal_iy <= 40 || crystal_iy > 60)) ||
577  (crystal_ix >= 4 && crystal_ix <= 5 && (crystal_iy <= 35 || crystal_iy > 65)) ||
578  (crystal_ix >= 6 && crystal_ix <= 8 && (crystal_iy <= 25 || crystal_iy > 75)) ||
579  (crystal_ix >= 9 && crystal_ix <= 13 && (crystal_iy <= 20 || crystal_iy > 80)) ||
580  (crystal_ix >= 14 && crystal_ix <= 15 && (crystal_iy <= 15 || crystal_iy > 85)) ||
581  (crystal_ix >= 16 && crystal_ix <= 20 && (crystal_iy <= 13 || crystal_iy > 87)) ||
582  (crystal_ix >= 21 && crystal_ix <= 25 && (crystal_iy <= 8 || crystal_iy > 92)) ||
583  (crystal_ix >= 26 && crystal_ix <= 35 && (crystal_iy <= 5 || crystal_iy > 95)) ||
584  (crystal_ix >= 36 && crystal_ix <= 39 && (crystal_iy <= 3 || crystal_iy > 97)) ||
585  (crystal_ix >= 98 && crystal_ix <= 100 && (crystal_iy <= 40 || crystal_iy > 60)) ||
586  (crystal_ix >= 96 && crystal_ix <= 97 && (crystal_iy <= 35 || crystal_iy > 65)) ||
587  (crystal_ix >= 93 && crystal_ix <= 95 && (crystal_iy <= 25 || crystal_iy > 75)) ||
588  (crystal_ix >= 88 && crystal_ix <= 92 && (crystal_iy <= 20 || crystal_iy > 80)) ||
589  (crystal_ix >= 86 && crystal_ix <= 87 && (crystal_iy <= 15 || crystal_iy > 85)) ||
590  (crystal_ix >= 81 && crystal_ix <= 85 && (crystal_iy <= 13 || crystal_iy > 87)) ||
591  (crystal_ix >= 76 && crystal_ix <= 80 && (crystal_iy <= 8 || crystal_iy > 92)) ||
592  (crystal_ix >= 66 && crystal_ix <= 75 && (crystal_iy <= 5 || crystal_iy > 95)) ||
593  (crystal_ix >= 62 && crystal_ix <= 65 && (crystal_iy <= 3 || crystal_iy > 97)) ||
594  ((crystal_ix == 40 || crystal_ix == 61) &&
595  ((crystal_iy >= 46 && crystal_iy <= 55) || crystal_iy <= 3 || crystal_iy > 97)) ||
596  ((crystal_ix == 41 || crystal_ix == 60) && crystal_iy >= 44 && crystal_iy <= 57) ||
597  ((crystal_ix == 42 || crystal_ix == 59) && crystal_iy >= 43 && crystal_iy <= 58) ||
598  ((crystal_ix == 43 || crystal_ix == 58) && crystal_iy >= 42 && crystal_iy <= 59) ||
599  ((crystal_ix == 44 || crystal_ix == 45 || crystal_ix == 57 || crystal_ix == 56) && crystal_iy >= 41 &&
600  crystal_iy <= 60) ||
601  (crystal_ix >= 46 && crystal_ix <= 55 && crystal_iy >= 40 && crystal_iy <= 61)) {
602  return valid;
603  }
604  valid = true;
605  return valid;
606 }
607 
608 //=================================================================================
609 void ElectronCalibrationUniv::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
610  //=================================================================================
611  using namespace edm;
612 
613  // Get EBRecHits
614  edm::Handle<EBRecHitCollection> EBphits;
615  iEvent.getByToken(ebRecHitToken_, EBphits);
616  if (!EBphits.isValid()) {
617  std::cerr << "Error! can't get the product EBRecHitCollection: " << std::endl;
618  }
619  const EBRecHitCollection* EBhits = EBphits.product(); // get a ptr to the product
620 
621  // Get EERecHits
622  edm::Handle<EERecHitCollection> EEphits;
623 
624  iEvent.getByToken(eeRecHitToken_, EEphits);
625  if (!EEphits.isValid()) {
626  std::cerr << "Error! can't get the product EERecHitCollection: " << std::endl;
627  }
628  const EERecHitCollection* EEhits = EEphits.product(); // get a ptr to the product
629 
630  // Get pixelElectrons
631  edm::Handle<reco::GsfElectronCollection> pElectrons;
632  iEvent.getByToken(gsfElectronToken_, pElectrons);
633  if (!pElectrons.isValid()) {
634  std::cerr << "Error! can't get the product ElectronCollection: " << std::endl;
635  }
636  const reco::GsfElectronCollection* electronCollection = pElectrons.product();
637  read_events++;
638  if (read_events % 1000 == 0)
639  std::cout << "read_events = " << read_events << std::endl;
640 
641  EventsAfterCuts->Fill(1);
642  if (!EBhits || !EEhits)
643  return;
644  EventsAfterCuts->Fill(2);
645  if (EBhits->empty() && EEhits->empty())
646  return;
647  EventsAfterCuts->Fill(3);
648  if (!electronCollection)
649  return;
650  EventsAfterCuts->Fill(4);
651  if (electronCollection->empty())
652  return;
653 
654  // ////////////////Need to recalibrate the events (copy code from EcalRecHitRecalib):
655 
659 
660  reco::GsfElectronCollection::const_iterator eleIt = electronCollection->begin();
661 
662  reco::GsfElectron highPtElectron;
663 
664  float highestElePt = 0.;
665  bool found = false;
666  for (eleIt = electronCollection->begin(); eleIt != electronCollection->end(); eleIt++) {
667  if (fabs(eleIt->eta()) > 2.4)
668  continue;
669  // if(eleIt->eta()<0.0) continue;
670 
671  if (eleIt->pt() > highestElePt) {
672  highestElePt = eleIt->pt();
673  highPtElectron = *eleIt;
674  found = true;
675  // std::cout<<" eleIt->pt( "<<eleIt->pt()<<" eleIt->eta() "<<eleIt->eta()<<std::endl;
676  }
677  }
678  EventsAfterCuts->Fill(5);
679  if (!found)
680  return;
681 
682  const reco::SuperCluster& sc = *(highPtElectron.superCluster());
683  // if(fabs(sc.eta())>1.479){std::cout<<" SC not in Barrel "<<sc.eta()<<std::endl;;}
684  // const std::vector<DetId> & v1 = sc.getHitsByDetId();
685 
686  std::vector<DetId> v1;
687  //Loop to fill the vector of DetIds
688  for (std::vector<std::pair<DetId, float> >::const_iterator idsIt = sc.hitsAndFractions().begin();
689  idsIt != sc.hitsAndFractions().end();
690  ++idsIt) {
691  v1.push_back(idsIt->first);
692  }
693 
694  DetId maxHitId;
695 
696  maxHitId = findMaxHit(v1, (EBhits), (EEhits));
697  //maxHitId = findMaxHit(v1,EBhits,EEhits);
698 
699  EventsAfterCuts->Fill(6);
700  if (maxHitId.null()) {
701  std::cout << " Null " << std::endl;
702  return;
703  }
704 
705  int maxCC_Eta = 0;
706  int maxCC_Phi = 0;
707  int Zside = 0;
708  if (maxHitId.subdetId() != 1) {
709  maxCC_Eta = ((EEDetId)maxHitId).ix();
710  maxCC_Phi = ((EEDetId)maxHitId).iy();
711  Zside = ((EEDetId)maxHitId).zside();
712  // std::cout<<" ++++++++ Zside "<<Zside<<std::endl;
713  } else {
714  maxCC_Eta = ((EBDetId)maxHitId).ieta();
715  maxCC_Phi = ((EBDetId)maxHitId).iphi();
716  }
717 
718  // if(maxCC_Eta>maxeta_ ) ;
719  // if(maxCC_Eta<mineta_ ) ;
720 
721  // number of events per crystal is set
722  // eventcrystal[maxCC_Eta][maxCC_Phi]+=1;
723  // if(eventcrystal[maxCC_Eta][maxCC_Phi] > numevent_) ;
724 
725  // fill cluster energy
726  std::vector<float> energy;
727  float energy3x3 = 0.;
728  float energy5x5 = 0.;
729  //Should be moved to cfg file!
730  int ClusterSize = clusterSize_;
731 
732  const CaloSubdetectorTopology* topology = theCaloTopology_->getSubdetectorTopology(DetId::Ecal, maxHitId.subdetId());
733  std::vector<DetId> NxNaroundMax = topology->getWindow(maxHitId, ClusterSize, ClusterSize);
734  //ToCompute 3x3
735  std::vector<DetId> S9aroundMax = topology->getWindow(maxHitId, 3, 3);
736 
737  EventsAfterCuts->Fill(7);
738  if ((int)NxNaroundMax.size() != ClusterSize * ClusterSize)
739  return;
740  EventsAfterCuts->Fill(8);
741  if (S9aroundMax.size() != 9)
742  return;
743 
744  // std::cout<<" ******** New Event "<<std::endl;
745 
746  EventsAfterCuts->Fill(9);
747  for (int icry = 0; icry < ClusterSize * ClusterSize; icry++) {
748  if (NxNaroundMax[icry].subdetId() == EcalBarrel) {
749  EBRecHitCollection::const_iterator itrechit;
750  itrechit = EBhits->find(NxNaroundMax[icry]);
751  if (itrechit == EBhits->end()) {
752  // std::cout << "EB DetId not in e25" << std::endl;
753  energy.push_back(0.);
754  energy5x5 += 0.;
755  continue;
756  }
757 
758  if (edm::isNotFinite(itrechit->energy())) {
759  std::cout << " nan energy " << std::endl;
760  return;
761  }
762  energy.push_back(itrechit->energy());
763  energy5x5 += itrechit->energy();
764 
765  //Summing in 3x3 to cut later on:
766  for (int tt = 0; tt < 9; tt++) {
767  if (NxNaroundMax[icry] == S9aroundMax[tt])
768  energy3x3 += itrechit->energy();
769  }
770  } else {
771  EERecHitCollection::const_iterator itrechit;
772 
773  itrechit = EEhits->find(NxNaroundMax[icry]);
774 
775  if (itrechit == EEhits->end()) {
776  // std::cout << "EE DetId not in e25" << std::endl;
777  // std::cout<<" ******** putting 0 "<<std::endl;
778  energy.push_back(0.);
779  energy5x5 += 0.;
780  continue;
781  }
782 
783  if (edm::isNotFinite(itrechit->energy())) {
784  std::cout << " nan energy " << std::endl;
785  return;
786  }
787  energy.push_back(itrechit->energy());
788  energy5x5 += itrechit->energy();
789 
790  //Summing in 3x3 to cut later on:
791  for (int tt = 0; tt < 9; tt++) {
792  if (NxNaroundMax[icry] == S9aroundMax[tt])
793  energy3x3 += itrechit->energy();
794  }
795  }
796  }
797  // if((read_events-50)%10000 ==0)cout << "++++++++++++ENERGY 5x5 " << energy5x5 << std::endl;
798  EventsAfterCuts->Fill(10);
799  // std::cout<<" ******** NxNaroundMax.size() "<<NxNaroundMax.size()<<std::endl;
800  // std::cout<<" ******** energy.size() "<<energy.size()<<std::endl;
801  if ((int)energy.size() != ClusterSize * ClusterSize)
802  return;
803 
804  if (maxHitId.subdetId() == EcalBarrel) {
805  GeneralMapBeforePt->Fill(maxCC_Eta, maxCC_Phi);
806  } else {
807  if (Zside < 0) {
808  GeneralMapEndCapMinusBeforePt->Fill(maxCC_Eta, maxCC_Phi);
809  } else {
810  GeneralMapEndCapPlusBeforePt->Fill(maxCC_Eta, maxCC_Phi);
811  }
812  }
813 
814  EventsAfterCuts->Fill(11);
815  if (highestElePt < elePt_)
816  return;
817 
818  if (maxHitId.subdetId() == EcalBarrel) {
819  GeneralMap->Fill(maxCC_Eta, maxCC_Phi);
820  } else {
821  if (Zside < 0) {
822  GeneralMapEndCapMinus->Fill(maxCC_Eta, maxCC_Phi);
823  } else {
824  GeneralMapEndCapPlus->Fill(maxCC_Eta, maxCC_Phi);
825  }
826  }
827 
828  EventsAfterCuts->Fill(12);
829  if (highPtElectron.classification() != elecclass_ && elecclass_ != -1)
830  return;
831 
832  float Ptrack_in =
833  sqrt(pow(highPtElectron.trackMomentumAtVtx().X(), 2) + pow(highPtElectron.trackMomentumAtVtx().Y(), 2) +
834  pow(highPtElectron.trackMomentumAtVtx().Z(), 2));
835 
836  float UncorrectedPatCalo =
837  sqrt(pow(highPtElectron.trackMomentumAtCalo().X(), 2) + pow(highPtElectron.trackMomentumAtCalo().Y(), 2) +
838  pow(highPtElectron.trackMomentumAtCalo().Z(), 2));
839 
840  float Ptrack_out =
841  sqrt(pow(highPtElectron.trackMomentumOut().X(), 2) + pow(highPtElectron.trackMomentumOut().Y(), 2) +
842  pow(highPtElectron.trackMomentumOut().Z(), 2));
843 
844  e9NoCuts->Fill(energy3x3);
845  e25NoCuts->Fill(energy5x5);
846  e9Overe25NoCuts->Fill(energy3x3 / energy5x5);
847  scENoCuts->Fill(sc.energy());
848 
849  trPNoCuts->Fill(UncorrectedPatCalo);
850 
851  EoPNoCuts->Fill(highPtElectron.eSuperClusterOverP());
852  e25OverScENoCuts->Fill(energy5x5 / sc.energy());
853 
854  E25oPNoCuts->Fill(energy5x5 / UncorrectedPatCalo);
855 
856  PinOverPoutNoCuts->Fill(
857  sqrt(pow(highPtElectron.trackMomentumAtVtx().X(), 2) + pow(highPtElectron.trackMomentumAtVtx().Y(), 2) +
858  pow(highPtElectron.trackMomentumAtVtx().Z(), 2)) /
859  sqrt(pow(highPtElectron.trackMomentumOut().X(), 2) + pow(highPtElectron.trackMomentumOut().Y(), 2) +
860  pow(highPtElectron.trackMomentumOut().Z(), 2)));
861  eSeedOverPoutNoCuts->Fill(highPtElectron.eSuperClusterOverP());
862 
863  MapCor1NoCuts->Fill(energy5x5 / UncorrectedPatCalo, energy5x5 / Ptrack_in);
864  MapCor2NoCuts->Fill(energy5x5 / UncorrectedPatCalo, highPtElectron.eSuperClusterOverP());
865  MapCor3NoCuts->Fill(energy5x5 / UncorrectedPatCalo, Ptrack_out / Ptrack_in);
866  MapCor4NoCuts->Fill(energy5x5 / UncorrectedPatCalo, energy5x5 / highPtElectron.p());
867  MapCor5NoCuts->Fill(energy5x5 / UncorrectedPatCalo, UncorrectedPatCalo / Ptrack_out);
868  MapCor6NoCuts->Fill(Ptrack_out / Ptrack_in, energy5x5 / Ptrack_in);
869  MapCor7NoCuts->Fill(Ptrack_out / Ptrack_in, UncorrectedPatCalo / Ptrack_out);
870  MapCor8NoCuts->Fill(energy5x5 / Ptrack_in, UncorrectedPatCalo / Ptrack_out);
871  MapCor9NoCuts->Fill(energy5x5 / UncorrectedPatCalo, highPtElectron.eSeedClusterOverPout());
872  MapCor10NoCuts->Fill(highPtElectron.eSeedClusterOverPout(), Ptrack_out / Ptrack_in);
873  MapCor11NoCuts->Fill(highPtElectron.eSeedClusterOverPout(), energy5x5 / Ptrack_in);
874 
875  PinMinPoutNoCuts->Fill((Ptrack_in - Ptrack_out) / Ptrack_in);
876 
877  Error1NoCuts->Fill(highPtElectron.trackMomentumError() / Ptrack_in);
878  Error2NoCuts->Fill(highPtElectron.trackMomentumError() / Ptrack_out);
879  Error3NoCuts->Fill(highPtElectron.trackMomentumError() / UncorrectedPatCalo);
880  eSeedOverPout2NoCuts->Fill(highPtElectron.eSeedClusterOverPout());
881 
882  hadOverEmNoCuts->Fill(highPtElectron.hadronicOverEm());
883 
884  //Cuts!
885  if ((energy3x3 / energy5x5) < cut1_)
886  return;
887  if ((Ptrack_out / Ptrack_in) < cut2_ || (Ptrack_out / Ptrack_in) > cut3_)
888  return;
889  if ((energy5x5 / Ptrack_in) < cutEPin1_ || (energy5x5 / Ptrack_in) > cutEPin2_)
890  return;
891  // if(!highPtElectron.ecalDriven())return;
892  // if(!highPtElectron.passingCutBasedPreselection())return;
893 
894  // // Apply Pietro cuts:
895  // EventsAfterCuts->Fill(13);
896  // //Module 1
897  // if(maxHitId.subdetId() == EcalBarrel){
898  // //Module 1
899  // if(maxCC_Eta <= 25){
900  // if(highPtElectron.eSuperClusterOverP()>1.05 || highPtElectron.eSuperClusterOverP()<0.95)return ;
901  // if(highPtElectron.eSeedClusterOverPout()>1.4 || highPtElectron.eSeedClusterOverPout()<0.90)return ;
902  // if((Ptrack_in- Ptrack_out) / Ptrack_in <-0.05 || (Ptrack_in- Ptrack_out) / Ptrack_in >0.2)return ;
903  // }else{
904  // //Module 2
905  // if( maxCC_Eta > 25&& maxCC_Eta <= 45){
906  // if(highPtElectron.eSuperClusterOverP()>1.05 || highPtElectron.eSuperClusterOverP()<0.95)return ;
907  // if(highPtElectron.eSeedClusterOverPout()>1.25 || highPtElectron.eSeedClusterOverPout()<0.90)return ;
908  // if((Ptrack_in- Ptrack_out) / Ptrack_in <-0.05 || (Ptrack_in- Ptrack_out) / Ptrack_in >0.2)return ;
909  // }else{
910  // //Module 3
911  // if( maxCC_Eta > 45&& maxCC_Eta <= 65){
912  // if(highPtElectron.eSuperClusterOverP()>1.05 || highPtElectron.eSuperClusterOverP()<0.95)return ;
913  // if(highPtElectron.eSeedClusterOverPout()>1.15 || highPtElectron.eSeedClusterOverPout()<0.90)return ;
914  // if((Ptrack_in- Ptrack_out) / Ptrack_in <-0.05 || (Ptrack_in- Ptrack_out) / Ptrack_in >0.15)return ;
915  // }else{
916  // if( maxCC_Eta > 65&& maxCC_Eta <= 85){
917  // if(highPtElectron.eSuperClusterOverP()>1.05 || highPtElectron.eSuperClusterOverP()<0.95)return ;
918  // if(highPtElectron.eSeedClusterOverPout()>1.15 || highPtElectron.eSeedClusterOverPout()<0.90)return ;
919  // if((Ptrack_in- Ptrack_out) / Ptrack_in <-0.05 || (Ptrack_in- Ptrack_out) / Ptrack_in >0.15)return ;
920  // }else{
921  // return;
922  // }
923  // }
924  // }
925  // }
926  // }else{
927  // //EndCapMinus Side:
928  // //EndCapPlus Side:
929  // int iR = sqrt((maxCC_Eta-50)*(maxCC_Eta-50) + (maxCC_Phi-50)*(maxCC_Phi-50));
930  // if( iR >= 22&& iR < 27){
931  // if(highPtElectron.eSuperClusterOverP()>1.05 || highPtElectron.eSuperClusterOverP()<0.95)return ;
932  // if(highPtElectron.eSeedClusterOverPout()>1.15 || highPtElectron.eSeedClusterOverPout()<0.90)return ;
933  // if((Ptrack_in- Ptrack_out) / Ptrack_in <-0.05 || (Ptrack_in- Ptrack_out) / Ptrack_in >0.2)return ;
934  // }else{
935  // if( iR >= 27&& iR < 32){
936  // if(highPtElectron.eSuperClusterOverP()>1.1 || highPtElectron.eSuperClusterOverP()<0.95)return ;
937  // if(highPtElectron.eSeedClusterOverPout()>1.25 || highPtElectron.eSeedClusterOverPout()<0.90)return ;
938  // if((Ptrack_in- Ptrack_out) / Ptrack_in <-0.05 || (Ptrack_in- Ptrack_out) / Ptrack_in >0.2)return ;
939  // }else{
940  // if( iR >= 32&& iR < 37){
941  // if(highPtElectron.eSuperClusterOverP()>1.05 || highPtElectron.eSuperClusterOverP()<0.95)return ;
942  // if(highPtElectron.eSeedClusterOverPout()>1.15 || highPtElectron.eSeedClusterOverPout()<0.90)return ;
943  // if((Ptrack_in- Ptrack_out) / Ptrack_in <-0.05 || (Ptrack_in- Ptrack_out) / Ptrack_in >0.2)return ;
944  // }else{
945  // if( iR >= 37&& iR < 42){
946  // if(highPtElectron.eSuperClusterOverP()>1.1 || highPtElectron.eSuperClusterOverP()<0.95)return ;
947  // if(highPtElectron.eSeedClusterOverPout()>1.15 || highPtElectron.eSeedClusterOverPout()<0.90)return ;
948  // if((Ptrack_in- Ptrack_out) / Ptrack_in <-0.05 || (Ptrack_in- Ptrack_out) / Ptrack_in >0.15)return ;
949  // }else{
950  // if( iR >= 42){
951  // if(highPtElectron.eSuperClusterOverP()>1.05 || highPtElectron.eSuperClusterOverP()<0.95)return ;
952  // if(highPtElectron.eSeedClusterOverPout()>1.15 || highPtElectron.eSeedClusterOverPout()<0.90)return ;
953  // if((Ptrack_in- Ptrack_out) / Ptrack_in <-0.05 || (Ptrack_in- Ptrack_out) / Ptrack_in >0.15)return ;
954  // }
955  // }
956  // }
957  // }
958  // }
959  // }
960 
961  if (maxHitId.subdetId() == EcalBarrel) {
962  E25oPvsEta->Fill(maxCC_Eta, energy5x5 / UncorrectedPatCalo);
963  } else {
964  float Radius = sqrt((maxCC_Eta) * (maxCC_Eta) + (maxCC_Phi) * (maxCC_Phi));
965  if (Zside < 0) {
966  E25oPvsEtaEndCapMinus->Fill(Radius, energy5x5 / UncorrectedPatCalo);
967  } else {
968  E25oPvsEtaEndCapPlus->Fill(Radius, energy5x5 / UncorrectedPatCalo);
969  }
970  }
971  e9->Fill(energy3x3);
972  e25->Fill(energy5x5);
973  e9Overe25->Fill(energy3x3 / energy5x5);
974  scE->Fill(sc.energy());
975  trP->Fill(UncorrectedPatCalo);
976 
977  EoP->Fill(highPtElectron.eSuperClusterOverP());
978  e25OverScE->Fill(energy5x5 / sc.energy());
979 
980  E25oP->Fill(energy5x5 / UncorrectedPatCalo);
981 
982  if (maxHitId.subdetId() == EcalBarrel) {
983  Map->Fill(maxCC_Eta, maxCC_Phi);
984  } else {
985  if (Zside < 0) {
986  MapEndCapMinus->Fill(maxCC_Eta, maxCC_Phi);
987  } else {
988  MapEndCapPlus->Fill(maxCC_Eta, maxCC_Phi);
989  }
990  }
991 
992  PinOverPout->Fill(sqrt(pow(highPtElectron.trackMomentumAtVtx().X(), 2) +
993  pow(highPtElectron.trackMomentumAtVtx().Y(), 2) +
994  pow(highPtElectron.trackMomentumAtVtx().Z(), 2)) /
995  sqrt(pow(highPtElectron.trackMomentumOut().X(), 2) + pow(highPtElectron.trackMomentumOut().Y(), 2) +
996  pow(highPtElectron.trackMomentumOut().Z(), 2)));
997  eSeedOverPout->Fill(highPtElectron.eSuperClusterOverP());
998 
999  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);
start
Definition: start.py:1
ElectronCalibrationUniv::ElectronCalibrationUniv
ElectronCalibrationUniv(const edm::ParameterSet &)
Definition: ElectronCalibrationUniv.cc:20
ProducerED_cfi.InputTag
InputTag
Definition: ProducerED_cfi.py:5
deDxTools::esConsumes
ESGetTokenH3DDVariant esConsumes(std::string const &Record, edm::ConsumesCollector &)
Definition: DeDxTools.cc:283
reco::CaloCluster::hitsAndFractions
const std::vector< std::pair< DetId, float > > & hitsAndFractions() const
Definition: CaloCluster.h:209
edm::EventSetup::getData
T const & getData(const ESGetToken< T, R > &iToken) const noexcept(false)
Definition: EventSetup.h:119
reco::GsfElectron::ecalDrivenSeed
bool ecalDrivenSeed() const
Definition: GsfElectron.h:158
ElectronCalibrationUniv::topologyToken_
const edm::ESGetToken< CaloTopology, CaloTopologyRecord > topologyToken_
Definition: ElectronCalibrationUniv.h:122
PixelTestBeamValidation_cfi.ClusterSize
ClusterSize
Definition: PixelTestBeamValidation_cfi.py:93
ElectronCalibrationUniv::EventMatrix
std::vector< std::vector< float > > EventMatrix
Definition: ElectronCalibrationUniv.h:108
edm::isNotFinite
constexpr bool isNotFinite(T x)
Definition: isFinite.h:9
EBDetId::iphi
int iphi() const
get the crystal iphi
Definition: EBDetId.h:51
reco::GsfElectron::trackMomentumError
float trackMomentumError() const
Definition: GsfElectron.h:884
edm::Handle::product
T const * product() const
Definition: Handle.h:70
EEDetId::ix
int ix() const
Definition: EEDetId.h:77
edm::SortedCollection< EcalRecHit >::const_iterator
std::vector< EcalRecHit >::const_iterator const_iterator
Definition: SortedCollection.h:80
ElectronCalibrationUniv::analyze
void analyze(const edm::Event &, const edm::EventSetup &) override
Definition: ElectronCalibrationUniv.cc:609
ElectronCalibrationUniv::UnivL3
MinL3AlgoUniv< DetId > * UnivL3
Definition: ElectronCalibrationUniv.h:117
conifer::pow
constexpr int pow(int x)
Definition: conifer.h:24
ElectronCalibrationUniv::theCaloTopology_
const CaloTopology * theCaloTopology_
Definition: ElectronCalibrationUniv.h:123
reco::GsfElectronCollection
std::vector< GsfElectron > GsfElectronCollection
collection of GsfElectron objects
Definition: GsfElectronFwd.h:14
hcalRecHitTable_cff.time
time
Definition: hcalRecHitTable_cff.py:14
hcalRecHitTable_cff.energy
energy
Definition: hcalRecHitTable_cff.py:13
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
ElectronCalibrationUniv::UnivEventIds
std::vector< std::vector< DetId > > UnivEventIds
Definition: ElectronCalibrationUniv.h:129
reco::GsfElectron::eSuperClusterOverP
float eSuperClusterOverP() const
Definition: GsfElectron.h:221
reco::GsfElectron::classification
Classification classification() const
Definition: GsfElectron.h:805
EBDetId::ieta
int ieta() const
get the crystal ieta
Definition: EBDetId.h:49
ElectronCalibrationUniv::ebRecHitToken_
const edm::EDGetTokenT< EBRecHitCollection > ebRecHitToken_
Definition: ElectronCalibrationUniv.h:119
reco::GsfElectron::trackMomentumAtCalo
math::XYZVectorF trackMomentumAtCalo() const
Definition: GsfElectron.h:269
reco::GsfElectron::eSeedClusterOverPout
float eSeedClusterOverPout() const
Definition: GsfElectron.h:223
DetId::null
constexpr bool null() const
is this a null id ?
Definition: DetId.h:59
pdwgLeptonRecoSkim_cfi.electronCollection
electronCollection
Definition: pdwgLeptonRecoSkim_cfi.py:5
command_line.start
start
Definition: command_line.py:165
reco::GsfElectron::trackMomentumOut
math::XYZVectorF trackMomentumOut() const
Definition: GsfElectron.h:270
ElectronCalibrationUniv::TestEEvalidDetId
bool TestEEvalidDetId(int crystal_ix, int crystal_iy, int iz)
Definition: ElectronCalibrationUniv.cc:571
iEvent
int iEvent
Definition: GenABIO.cc:224
tt
Definition: TTTypes.h:54
reco::LeafCandidate::p
double p() const final
magnitude of momentum vector
Definition: LeafCandidate.h:123
validateGeometry_cfg.valid
valid
Definition: validateGeometry_cfg.py:21
HouseholderDecomposition::iterate
std::vector< float > iterate(const std::vector< std::vector< float > > &eventMatrix, const std::vector< int > &VmaxCeta, const std::vector< int > &VmaxCphi, const std::vector< float > &energyVector, const int &nIter, const bool &normalizeFlag=false)
Definition: HouseholderDecomposition.cc:123
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
edm::Transition
Transition
Definition: Transition.h:12
reco::GsfElectron::trackMomentumAtVtx
math::XYZVectorF trackMomentumAtVtx() const
Definition: GsfElectron.h:268
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
DetId::subdetId
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:48
CaloSubdetectorTopology::getWindow
virtual std::vector< DetId > getWindow(const DetId &id, const int &northSouthSize, const int &eastWestSize) const
Definition: CaloSubdetectorTopology.cc:4
ElectronCalibrationUniv::Univsolution
std::map< DetId, float > Univsolution
Definition: ElectronCalibrationUniv.h:130
ElectronCalibrationUniv::MyHH
HouseholderDecomposition * MyHH
Definition: ElectronCalibrationUniv.h:185
edm::SortedCollection::end
const_iterator end() const
Definition: SortedCollection.h:267
reco::CaloCluster::energy
double energy() const
cluster energy
Definition: CaloCluster.h:148
DetId
Definition: DetId.h:17
mps_fire.end
end
Definition: mps_fire.py:242
HouseholderDecomposition::runRegional
std::vector< float > runRegional(const std::vector< std::vector< float > > &eventMatrix, const std::vector< int > &VmaxCeta, const std::vector< int > &VmaxCphi, const std::vector< float > &energyVector, const int &nIter, const int &regLength=5)
Definition: HouseholderDecomposition.cc:30
EEDetId::zside
int zside() const
Definition: EEDetId.h:71
MTVHistoProducerAlgoForTrackerBlock_cfi.maxHit
maxHit
Definition: MTVHistoProducerAlgoForTrackerBlock_cfi.py:37
CaloTopology::getSubdetectorTopology
const CaloSubdetectorTopology * getSubdetectorTopology(const DetId &id) const
access the subdetector Topology for the given subdetector directly
Definition: CaloTopology.cc:17
ElectronCalibrationUniv::gsfElectronToken_
const edm::EDGetTokenT< reco::GsfElectronCollection > gsfElectronToken_
Definition: ElectronCalibrationUniv.h:121
DMR_cfg.cerr
cerr
Definition: DMR_cfg.py:78
edm::HandleBase::isValid
bool isValid() const
Definition: HandleBase.h:70
ElectronCalibrationUniv::endRun
void endRun(edm::Run const &, edm::EventSetup const &) override
Definition: ElectronCalibrationUniv.cc:243
calibXMLwriter::writeLine
void writeLine(EBDetId const &, float)
Definition: calibXMLwriter.cc:32
edm::SortedCollection::find
iterator find(key_type k)
Definition: SortedCollection.h:240
reco
fixed size matrix
Definition: AlignmentAlgorithmBase.h:46
edm
HLT enums.
Definition: AlignableModifier.h:19
gather_cfg.cout
cout
Definition: gather_cfg.py:144
MinL3AlgoUniv::iterate
IDmap iterate(const std::vector< std::vector< float > > &eventMatrix, const std::vector< std::vector< IDdet > > &idMatrix, const std::vector< float > &energyVector, const int &nIter, const bool &normalizeFlag=false)
Definition: MinL3AlgoUniv.h:76
ElectronCalibrationUniv::EnergyVector
std::vector< float > EnergyVector
Definition: ElectronCalibrationUniv.h:106
reco::GsfElectron::hadronicOverEm
float hadronicOverEm() const
Definition: GsfElectron.h:500
reco::GsfElectron::superCluster
SuperClusterRef superCluster() const override
reference to a SuperCluster
Definition: GsfElectron.h:155
MinL3Algorithm::iterate
std::vector< float > iterate(const std::vector< std::vector< float > > &eventMatrix, const std::vector< int > &VmaxCeta, const std::vector< int > &VmaxCphi, const std::vector< float > &energyVector, const int &nIter, const bool &normalizeFlag=false)
Definition: MinL3Algorithm.cc:35
edm::Run
Definition: Run.h:45
ElectronCalibrationUniv::eeRecHitToken_
const edm::EDGetTokenT< EERecHitCollection > eeRecHitToken_
Definition: ElectronCalibrationUniv.h:120
ElectronCalibrationUniv::beginRun
void beginRun(edm::Run const &, edm::EventSetup const &) override
Definition: ElectronCalibrationUniv.cc:235
EEDetId::iy
int iy() const
Definition: EEDetId.h:83
ElectronCalibrationUniv::findMaxHit
DetId findMaxHit(const std::vector< DetId > &v1, const EBRecHitCollection *EBhits, const EERecHitCollection *EEhits)
Definition: ElectronCalibrationUniv.cc:533
Generated for CMSSW Reference Manual by doxygen 1.8.14