CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
ElectronCalibration.cc
Go to the documentation of this file.
1 
3 
21 #include "TFile.h"
22 #include "TH1.h"
23 #include "TH2.h"
24 #include "TF1.h"
25 #include "TRandom.h"
26 
27 #include <iostream>
28 #include <string>
29 #include <stdexcept>
30 #include <vector>
31 
34 
36  rootfile_ = iConfig.getParameter<std::string>("rootfile");
37  recHitLabel_ = iConfig.getParameter<edm::InputTag>("ebRecHitsLabel");
38  electronLabel_ = iConfig.getParameter<edm::InputTag>("electronLabel");
39  trackLabel_ = iConfig.getParameter<edm::InputTag>("trackLabel");
40  calibAlgo_ = iConfig.getParameter<std::string>("CALIBRATION_ALGO");
41  std::cout << " The used Algorithm is " << calibAlgo_ << std::endl;
42  keventweight_ = iConfig.getParameter<int>("keventweight");
43  ClusterSize_ = iConfig.getParameter<int>("Clustersize");
44  ElePt_ = iConfig.getParameter<double>("ElePt");
45  maxeta_ = iConfig.getParameter<int>("maxeta");
46  mineta_ = iConfig.getParameter<int>("mineta");
47  maxphi_ = iConfig.getParameter<int>("maxphi");
48  minphi_ = iConfig.getParameter<int>("minphi");
49  cut1_ = iConfig.getParameter<double>("cut1");
50  cut2_ = iConfig.getParameter<double>("cut2");
51  cut3_ = iConfig.getParameter<double>("cut3");
52  elecclass_ = iConfig.getParameter<int>("elecclass");
53  std::cout << " The electronclass is " << elecclass_ << std::endl;
54  numevent_ = iConfig.getParameter<int>("numevent");
55  miscalibfile_ = iConfig.getParameter<std::string>("miscalibfile");
56 
57  cutEPCalo1_ = iConfig.getParameter<double>("cutEPCaloMin");
58  cutEPCalo2_ = iConfig.getParameter<double>("cutEPCaloMax");
59  cutEPin1_ = iConfig.getParameter<double>("cutEPinMin");
60  cutEPin2_ = iConfig.getParameter<double>("cutEPinMax");
61  cutCalo1_ = iConfig.getParameter<double>("cutCaloMin");
62  cutCalo2_ = iConfig.getParameter<double>("cutCaloMax");
63 
64  cutESeed_ = iConfig.getParameter<double>("cutESeed");
65 }
66 
68 
69 //========================================================================
71  //========================================================================
72  f = new TFile(rootfile_.c_str(), "RECREATE");
73 
74  // Book histograms
75  e9 = new TH1F("e9", "E9 energy", 300, 0., 150.);
76  e25 = new TH1F("e25", "E25 energy", 300, 0., 150.);
77  scE = new TH1F("scE", "SC energy", 300, 0., 150.);
78  trP = new TH1F("trP", "Trk momentum", 300, 0., 150.);
79  EoP = new TH1F("EoP", "EoP", 600, 0., 3.);
80  EoP_all = new TH1F("EoP_all", "EoP_all", 600, 0., 3.);
81 
82  if (elecclass_ == 0 || elecclass_ == -1) {
83  calibs = new TH1F("calib", "Calibration constants", 4000, 0.5, 2.);
84  } else {
85  calibs = new TH1F("calib", "Calibration constants", 800, 0.5, 2.);
86  }
87 
88  e25OverScE = new TH1F("e25OverscE", "E25 / SC energy", 400, 0., 2.);
89  E25oP = new TH1F("E25oP", "E25 / P", 1200, 0., 1.5);
90 
91  Map = new TH2F("Map", "Nb Events in Crystal", 85, 1, 85, 70, 5, 75);
92  e9Overe25 = new TH1F("e9Overe25", "E9 / E25", 400, 0., 2.);
93  Map3Dcalib = new TH2F("3Dcalib", "3Dcalib", 85, 1, 85, 70, 5, 75);
94 
95  MapCor1 = new TH2F("MapCor1", "Correlation E25/Pcalo versus E25/Pin", 100, 0., 5., 100, 0., 5.);
96  MapCor2 = new TH2F("MapCor2", "Correlation E25/Pcalo versus E/P", 100, 0., 5., 100, 0., 5.);
97  MapCor3 = new TH2F("MapCor3", "Correlation E25/Pcalo versus Pout/Pin", 100, 0., 5., 100, 0., 5.);
98  MapCor4 = new TH2F("MapCor4", "Correlation E25/Pcalo versus E25/highestP", 100, 0., 5., 100, 0., 5.);
99  MapCor5 = new TH2F("MapCor5", "Correlation E25/Pcalo versus Pcalo/Pout", 100, 0., 5., 100, 0., 5.);
100  MapCor6 = new TH2F("MapCor6", "Correlation Pout/Pin versus E25/Pin", 100, 0., 5., 100, 0., 5.);
101  MapCor7 = new TH2F("MapCor7", "Correlation Pout/Pin versus Pcalo/Pout", 100, 0., 5., 100, 0., 5.);
102  MapCor8 = new TH2F("MapCor8", "Correlation E25/Pin versus Pcalo/Pout", 100, 0., 5., 100, 0., 5.);
103  MapCor9 = new TH2F("MapCor9", "Correlation E25/Pcalo versus Eseed/Pout", 100, 0., 5., 100, 0., 5.);
104  MapCor10 = new TH2F("MapCor10", "Correlation Eseed/Pout versus Pout/Pin", 100, 0., 5., 100, 0., 5.);
105  MapCor11 = new TH2F("MapCor11", "Correlation Eseed/Pout versus E25/Pin", 100, 0., 5., 100, 0., 5.);
106  MapCorCalib =
107  new TH2F("MapCorCalib", "Correlation Miscalibration versus Calibration constants", 100, 0.5, 1.5, 100, 0.5, 1.5);
108 
109  PinMinPout = new TH1F("PinMinPout", "(Pin - Pout)/Pin", 600, -2.0, 2.0);
110 
111  if (elecclass_ == 0 || elecclass_ == -1) {
112  calibinter = new TH1F("calibinter", "internal calibration constants", 2000, 0.5, 2.);
113  PinOverPout = new TH1F("PinOverPout", "pinOverpout", 600, 0., 3.);
114  eSeedOverPout = new TH1F("eSeedOverPout", "eSeedOverpout ", 600, 0., 3.);
115  MisCalibs = new TH1F("MisCalibs", "Miscalibration constants", 4000, 0.5, 2.);
116  RatioCalibs = new TH1F("RatioCalibs", "Ratio in Calibration Constants", 4000, 0.5, 2.0);
117  DiffCalibs = new TH1F("DiffCalibs", "Difference in Calibration constants", 4000, -1.0, 1.0);
118  } else {
119  calibinter = new TH1F("calibinter", "internal calibration constants", 400, 0.5, 2.);
120  PinOverPout = new TH1F("PinOverPout", "pinOverpout", 600, 0., 3.);
121  eSeedOverPout = new TH1F("eSeedOverPout", "eSeedOverpout ", 600, 0., 3.);
122  MisCalibs = new TH1F("MisCalibs", "Miscalibration constants", 800, 0.5, 2.);
123  RatioCalibs = new TH1F("RatioCalibs", "Ratio in Calibration Constants", 800, 0.5, 2.0);
124  DiffCalibs = new TH1F("DiffCalibs", "Difference in Calibration constants", 800, -1.0, 1.0);
125  }
126  Error1 = new TH1F("Error1", "DeltaP/Pin", 800, -1.0, 1.0);
127  Error2 = new TH1F("Error2", "DeltaP/Pout", 800, -1.0, 1.0);
128  Error3 = new TH1F("Error3", "DeltaP/Pcalo", 800, -1.0, 1.0);
129  eSeedOverPout2 = new TH1F("eSeedOverPout2", "eSeedOverpout (No Supercluster)", 600, 0., 4.);
130  hadOverEm = new TH1F("hadOverEm", "Had/EM distribution", 600, -2., 2.);
131 
132  // Book histograms
133  Map3DcalibNoCuts = new TH2F("3DcalibNoCuts", "3Dcalib (Before Cuts)", 85, 1, 85, 70, 5, 75);
134  e9NoCuts = new TH1F("e9NoCuts", "E9 energy (Before Cuts)", 300, 0., 150.);
135  e25NoCuts = new TH1F("e25NoCuts", "E25 energy (Before Cuts)", 300, 0., 150.);
136  scENoCuts = new TH1F("scENoCuts", "SC energy (Before Cuts)", 300, 0., 150.);
137  trPNoCuts = new TH1F("trPNoCuts", "Trk momentum (Before Cuts)", 300, 0., 150.);
138  EoPNoCuts = new TH1F("EoPNoCuts", "EoP (Before Cuts)", 600, 0., 3.);
139  if (elecclass_ == 0 || elecclass_ == -1) {
140  calibsNoCuts = new TH1F("calibNoCuts", "Calibration constants (Before Cuts)", 4000, 0., 2.);
141  } else {
142  calibsNoCuts = new TH1F("calibNoCuts", "Calibration constants (Before Cuts)", 800, 0., 2.);
143  }
144  e25OverScENoCuts = new TH1F("e25OverscENoCuts", "E25 / SC energy (Before Cuts)", 400, 0., 2.);
145  E25oPNoCuts = new TH1F("E25oPNoCuts", "E25 / P (Before Cuts)", 1200, 0., 1.5);
146  MapNoCuts = new TH2F("MapNoCuts", "Nb Events in Crystal (Before Cuts)", 85, 1, 85, 70, 5, 75);
147  e9Overe25NoCuts = new TH1F("e9Overe25NoCuts", "E9 / E25 (Before Cuts)", 400, 0., 2.);
148  PinOverPoutNoCuts = new TH1F("PinOverPoutNoCuts", "pinOverpout (Before Cuts)", 600, 0., 3.);
149  eSeedOverPoutNoCuts = new TH1F(" eSeedOverPoutNoCuts", "eSeedOverpout (Before Cuts) ", 600, 0., 4.);
150  PinMinPoutNoCuts = new TH1F("PinMinPoutNoCuts", "(Pin - Pout)/Pin (Before Cuts)", 600, -2.0, 2.0);
151 
152  RatioCalibsNoCuts = new TH1F("RatioCalibsNoCuts", "Ratio in Calibration Constants (Before Cuts)", 4000, 0.5, 2.0);
153  DiffCalibsNoCuts = new TH1F("DiffCalibsNoCuts", "Difference in Calibration constants (Before Cuts)", 4000, -1.0, 1.0);
154  calibinterNoCuts = new TH1F("calibinterNoCuts", "internal calibration constants", 2000, 0.5, 2.);
155 
156  MapCor1NoCuts =
157  new TH2F("MapCor1NoCuts", "Correlation E25/PatCalo versus E25/Pin (Before Cuts)", 100, 0., 5., 100, 0., 5.);
158  MapCor2NoCuts =
159  new TH2F("MapCor2NoCuts", "Correlation E25/PatCalo versus E/P (Before Cuts)", 100, 0., 5., 100, 0., 5.);
160  MapCor3NoCuts =
161  new TH2F("MapCor3NoCuts", "Correlation E25/PatCalo versus Pout/Pin (Before Cuts)", 100, 0., 5., 100, 0., 5.);
162  MapCor4NoCuts =
163  new TH2F("MapCor4NoCuts", "Correlation E25/PatCalo versus E25/highestP (Before Cuts)", 100, 0., 5., 100, 0., 5.);
164  MapCor5NoCuts =
165  new TH2F("MapCor5NoCuts", "Correlation E25/Pcalo versus Pcalo/Pout (Before Cuts)", 100, 0., 5., 100, 0., 5.);
166  MapCor6NoCuts =
167  new TH2F("MapCor6NoCuts", "Correlation Pout/Pin versus E25/Pin (Before Cuts)", 100, 0., 5., 100, 0., 5.);
168  MapCor7NoCuts =
169  new TH2F("MapCor7NoCuts", "Correlation Pout/Pin versus Pcalo/Pout (Before Cuts)", 100, 0., 5., 100, 0., 5.);
170  MapCor8NoCuts =
171  new TH2F("MapCor8NoCuts", "Correlation E25/Pin versus Pcalo/Pout (Before Cuts)", 100, 0., 5., 100, 0., 5.);
172  MapCor9NoCuts =
173  new TH2F("MapCor9NoCuts", "Correlation E25/Pcalo versus Eseed/Pout (Before Cuts)", 100, 0., 5., 100, 0., 5.);
175  new TH2F("MapCor10NoCuts", "Correlation Eseed/Pout versus Pout/Pin (Before Cuts)", 100, 0., 5., 100, 0., 5.);
177  new TH2F("MapCor11NoCuts", "Correlation Eseed/Pout versus E25/Pin (Before Cuts)", 100, 0., 5., 100, 0., 5.);
178  MapCorCalibNoCuts = new TH2F("MapCorCalibNoCuts",
179  "Correlation Miscalibration versus Calibration constants (Before Cuts)",
180  100,
181  0.,
182  3.,
183  100,
184  0.,
185  3.);
186 
187  Error1NoCuts = new TH1F("Eror1NoCuts", "DeltaP/Pin (Before Cuts)", 800, -1.0, 1.0);
188  Error2NoCuts = new TH1F("Error2NoCuts", "DeltaP/Pout (Before Cuts)", 800, -1.0, 1.0);
189  Error3NoCuts = new TH1F("Error3NoCuts", "DeltaP/Pcalo (Before Cuts)", 800, -1.0, 1.0);
190  eSeedOverPout2NoCuts = new TH1F("eSeedOverPout2NoCuts", "eSeedOverpout (No Supercluster, Before Cuts)", 600, 0., 4.);
191  hadOverEmNoCuts = new TH1F("hadOverEmNoCuts", "Had/EM distribution (Before Cuts)", 600, -2., 2.);
192 
193  //Book histograms after ESeed cut
194  MapCor1ESeed =
195  new TH2F("MapCor1ESeed", "Correlation E25/Pcalo versus E25/Pin (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
196  MapCor2ESeed =
197  new TH2F("MapCor2ESeed", "Correlation E25/Pcalo versus E/P (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
198  MapCor3ESeed = new TH2F(
199  "MapCor3ESeed", "Correlation E25/Pcalo versus Pout/Pin (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
200  MapCor4ESeed = new TH2F(
201  "MapCor4ESeed", "Correlation E25/Pcalo versus E25/highestP (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
202  MapCor5ESeed = new TH2F(
203  "MapCor5ESeed", "Correlation E25/Pcalo versus Pcalo/Pout (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
204  MapCor6ESeed =
205  new TH2F("MapCor6ESeed", "Correlation Pout/Pin versus E25/Pin (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
206  MapCor7ESeed = new TH2F(
207  "MapCor7ESeed", "Correlation Pout/Pin versus Pcalo/Pout (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
208  MapCor8ESeed = new TH2F(
209  "MapCor8ESeed", "Correlation E25/Pin versus Pcalo/Pout (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
210  MapCor9ESeed = new TH2F(
211  "MapCor9ESeed", "Correlation E25/Pcalo versus Eseed/Pout (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
212  MapCor10ESeed = new TH2F(
213  "MapCor10ESeed", "Correlation Eseed/Pout versus Pout/Pin (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
214  MapCor11ESeed = new TH2F(
215  "MapCor11ESeed", "Correlation Eseed/Pout versus E25/Pin (after Eseed/Pout cut)", 100, 0., 5., 100, 0., 5.);
216 
218  new TH1F("eSeedOverPout2ESeed", "eSeedOverpout (No Supercluster, after Eseed/Pout cut)", 600, 0., 4.);
219 
220  hadOverEmESeed = new TH1F("hadOverEmESeed", "Had/EM distribution (after Eseed/Pout cut)", 600, -2., 2.);
221 
222  //Book histograms without any cut
223  GeneralMap = new TH2F("GeneralMap", "Map without any cuts", 85, 1, 85, 70, 5, 75);
224 
226  etaMin = mineta_;
227  etaMax = maxeta_;
228  phiMin = minphi_;
229  phiMax = maxphi_;
230  if (calibAlgo_ == "L3") {
232  } else {
233  if (calibAlgo_ == "HH" || calibAlgo_ == "HHReg") {
235  } else {
236  std::cout << " Name of Algorithm is not recognize " << calibAlgo_ << " Should be either L3, HH or HHReg. Abort! "
237  << std::endl;
238  }
239  }
240  read_events = 0;
241 
242  // get Region to be calibrated
244 
245  oldCalibs.resize(ReducedMap.size(), 0.);
246 
247  // table is set to zero
248  for (int phi = 0; phi < 360; phi++) {
249  for (int eta = 0; eta < 171; eta++) {
250  eventcrystal[eta][phi] = 0;
251  }
252  }
253 
254  std::cout << " Begin JOB " << std::endl;
255 }
256 
257 //========================================================================
258 
260  //========================================================================
261 
262  int nIterations = 10;
263  if (calibAlgo_ == "L3") {
265  } else {
266  if (calibAlgo_ == "HH") {
268  } else {
269  if (calibAlgo_ == "HHReg") {
271  } else {
272  std::cout << " Calibration not run due to problem in Algo Choice..." << std::endl;
273  return;
274  }
275  }
276  }
277  for (int ii = 0; ii < (int)solution.size(); ii++) {
278  std::cout << "solution[" << ii << "] = " << solution[ii] << std::endl;
279  calibs->Fill(solution[ii]);
280  }
281 
282  newCalibs.resize(ReducedMap.size(), 0.);
283 
284  calibXMLwriter write_calibrations;
285 
286  FILE* MisCalib;
287  MisCalib = fopen(miscalibfile_.c_str(), "r");
288  int fileStatus = 1;
289  int eta = -1;
290  int phi = -1;
291  float coeff = -1;
292 
293  std::map<EBDetId, float> OldCoeff;
294 
295  while (fileStatus != EOF) {
296  fileStatus = fscanf(MisCalib, "%d %d %f\n", &eta, &phi, &coeff);
297  if (eta != -1 && phi != -1 && coeff != -1) {
298  // std::cout<<" We have read correctly the coefficient " << coeff << " corresponding to eta "<<eta<<" and phi "<<phi<<std::endl;
299  OldCoeff.insert(std::make_pair(EBDetId(eta, phi, EBDetId::ETAPHIMODE), coeff));
300  }
301  }
302 
303  fclose(MisCalib);
304 
305  int icry = 0;
306  CalibrationCluster::CalibMap::iterator itmap;
307  for (itmap = ReducedMap.begin(); itmap != ReducedMap.end(); itmap++) {
308  newCalibs[icry] = solution[icry];
309 
310  write_calibrations.writeLine(itmap->first, newCalibs[icry]);
311  float Compare = 1.;
312  std::map<EBDetId, float>::iterator iter = OldCoeff.find(itmap->first);
313  if (iter != OldCoeff.end())
314  Compare = iter->second;
315 
316  if ((itmap->first).ieta() > mineta_ && (itmap->first).ieta() < maxeta_ && (itmap->first).iphi() > minphi_ &&
317  (itmap->first).iphi() < maxphi_) {
318  Map3Dcalib->Fill((itmap->first).ieta(), (itmap->first).iphi(), newCalibs[icry] * Compare);
319  MisCalibs->Fill(Compare);
320  }
321  if ((itmap->first).ieta() < mineta_ + 2) {
322  icry++;
323  continue;
324  }
325  if ((itmap->first).ieta() > maxeta_ - 2) {
326  icry++;
327  continue;
328  }
329  if ((itmap->first).iphi() < minphi_ + 2) {
330  icry++;
331  continue;
332  }
333  if ((itmap->first).iphi() > maxphi_ - 2) {
334  icry++;
335  continue;
336  }
337 
338  calibinter->Fill(newCalibs[icry]);
339  DiffCalibs->Fill(newCalibs[icry] - 1. / Compare);
340  RatioCalibs->Fill(newCalibs[icry] * Compare);
341  MapCorCalib->Fill(1. / Compare, newCalibs[icry]);
342  icry++;
343  }
344 
345  if (calibAlgo_ == "L3") {
348  } else {
349  if (calibAlgo_ == "HH") {
351  } else {
352  if (calibAlgo_ == "HHReg") {
354  } else {
355  std::cout << " Calibration not run due to problem in AlgoChoice..." << std::endl;
356  return;
357  }
358  }
359  }
360  for (int ii = 0; ii < (int)solutionNoCuts.size(); ii++) {
362  }
363  int icryp = 0;
364  CalibrationCluster::CalibMap::iterator itmapp;
365  for (itmapp = ReducedMap.begin(); itmapp != ReducedMap.end(); itmapp++) {
366  newCalibs[icryp] = solutionNoCuts[icryp];
367  float Compare2 = 1.;
368  std::map<EBDetId, float>::iterator iter2 = OldCoeff.find(itmapp->first);
369  if (iter2 != OldCoeff.end())
370  Compare2 = iter2->second;
371 
372  if ((itmapp->first).ieta() > mineta_ && (itmapp->first).ieta() < maxeta_ && (itmapp->first).iphi() > minphi_ &&
373  (itmapp->first).iphi() < maxphi_)
374  Map3DcalibNoCuts->Fill((itmapp->first).ieta(), (itmapp->first).iphi(), newCalibs[icryp] * Compare2);
375  if ((itmapp->first).ieta() < mineta_ + 2) {
376  icryp++;
377  continue;
378  }
379  if ((itmapp->first).ieta() > maxeta_ - 2) {
380  icryp++;
381  continue;
382  }
383  if ((itmapp->first).iphi() < minphi_ + 2) {
384  icryp++;
385  continue;
386  }
387  if ((itmapp->first).iphi() > maxphi_ - 2) {
388  icryp++;
389  continue;
390  }
391  calibinterNoCuts->Fill(newCalibs[icryp]);
392  DiffCalibsNoCuts->Fill(newCalibs[icryp] - 1. / (Compare2));
393  RatioCalibsNoCuts->Fill(newCalibs[icryp] * Compare2);
394  MapCorCalibNoCuts->Fill(1. / Compare2, newCalibs[icryp]);
395  icryp++;
396  }
397 
399 
400  std::cout << " " << std::endl;
401  std::cout << "************* STATISTICS **************" << std::endl;
402  std::cout << " Events Studied " << read_events << std::endl;
403 
405 
406  f->Write();
407 
408  f->Close();
409 }
410 
411 //=================================================================================
413  //=================================================================================
414 
415  EcalRecHitCollection ecrh = *phits;
417  int count = 0;
418  EBDetId save;
419  float en_save = 0;
420  for (it = ecrh.begin(); it != ecrh.end(); it++) {
421  EBDetId p = EBDetId(it->id().rawId());
422  if (it->energy() > en_save) {
423  en_save = it->energy();
424  save = p;
425  }
426  count++;
427  }
428  return save;
429 }
430 
431 //=================================================================================
432 EBDetId ElectronCalibration::findMaxHit2(const std::vector<DetId>& v1, const EBRecHitCollection* hits) {
433  //=================================================================================
434 
435  double currEnergy = 0.;
436  EBDetId maxHit;
437 
438  for (std::vector<DetId>::const_iterator idsIt = v1.begin(); idsIt != v1.end(); ++idsIt) {
439  if (idsIt->subdetId() != 1)
440  continue;
442  itrechit = hits->find(*idsIt);
443 
444  if (itrechit == hits->end()) {
445  std::cout << "ElectronCalibration::findMaxHit2: rechit not found! " << std::endl;
446  continue;
447  }
448  if (itrechit->energy() > currEnergy) {
449  currEnergy = itrechit->energy();
450  maxHit = *idsIt;
451  }
452  }
453 
454  return maxHit;
455 }
456 
457 //=================================================================================
459  //=================================================================================
460  using namespace edm;
461 
462  // Get EBRecHits
464  iEvent.getByLabel(recHitLabel_, phits);
465  if (!phits.isValid()) {
466  std::cerr << "Error! can't get the product EBRecHitCollection: " << std::endl;
467  }
468 
469  const EBRecHitCollection* hits = phits.product(); // get a ptr to the product
470 
471  // Get pixelElectrons
473 
474  iEvent.getByLabel(electronLabel_, pElectrons);
475  if (!pElectrons.isValid()) {
476  std::cerr << "Error! can't get the product ElectronCollection: " << std::endl;
477  }
478 
479  const reco::GsfElectronCollection* electronCollection = pElectrons.product();
480  read_events++;
481  if (read_events % 1000 == 0)
482  std::cout << "read_events = " << read_events << std::endl;
483 
484  if (!hits)
485  return;
486  if (hits->empty())
487  return;
488  if (!electronCollection)
489  return;
490  if (electronCollection->empty())
491  return;
492 
494  // START HERE....
496  reco::GsfElectronCollection::const_iterator eleIt = electronCollection->begin();
497 
498  reco::GsfElectron highPtElectron;
499 
500  float highestElePt = 0.;
501  bool found = false;
502  for (eleIt = electronCollection->begin(); eleIt != electronCollection->end(); eleIt++) {
503  //Comments
504  if (fabs(eleIt->eta()) > (maxeta_ + 3) * 0.0175)
505  continue;
506  if (eleIt->eta() < (mineta_ - 3) * 0.0175)
507  continue;
508 
509  if (eleIt->pt() > highestElePt) {
510  highestElePt = eleIt->pt();
511  highPtElectron = *eleIt;
512  found = true;
513  }
514  }
515  if (highestElePt < ElePt_)
516  return;
517  if (!found)
518  return;
519  const reco::SuperCluster& sc = *(highPtElectron.superCluster());
520  if (fabs(sc.eta()) > (maxeta_ + 3) * 0.0175) {
521  std::cout << "++++ Problem with electron, electron eta is " << highPtElectron.eta() << " while SC is " << sc.eta()
522  << std::endl;
523  return;
524  }
525  // std::cout << "track eta = " << highPtElectron.eta() << std::endl;
526  // std::cout << "track phi = " << highPtElectron.phi() << std::endl;
527 
528  std::vector<DetId> v1;
529  //Loop to fill the vector of DetIds
530  for (std::vector<std::pair<DetId, float> >::const_iterator idsIt = sc.hitsAndFractions().begin();
531  idsIt != sc.hitsAndFractions().end();
532  ++idsIt) {
533  v1.push_back(idsIt->first);
534  }
535 
536  //getHitsByDetId(); //Change function name
537  EBDetId maxHitId;
538 
539  maxHitId = findMaxHit2(v1, hits);
540 
541  if (maxHitId.null()) {
542  std::cout << " Null " << std::endl;
543  return;
544  }
545 
546  int maxCC_Eta = maxHitId.ieta();
547  int maxCC_Phi = maxHitId.iphi();
548 
549  if (maxCC_Eta > maxeta_)
550  return;
551  if (maxCC_Eta < mineta_)
552  return;
553  if (maxCC_Phi > maxphi_)
554  return;
555  if (maxCC_Phi < minphi_)
556  return;
557 
558  // number of events per crystal is set
559  if (numevent_ > 0) {
560  eventcrystal[maxCC_Eta + 85][maxCC_Phi - 1] += 1;
561  if (eventcrystal[maxCC_Eta + 85][maxCC_Phi - 1] > numevent_)
562  return;
563  }
564 
565  std::vector<EBDetId> Xtals5x5 = calibCluster.get5x5Id(maxHitId);
566 
567  if ((int)Xtals5x5.size() != ClusterSize_ * ClusterSize_)
568  return;
569 
570  // fill cluster energy
571  std::vector<float> energy;
572  float energy3x3 = 0.;
573  float energy5x5 = 0.;
574 
575  for (int icry = 0; icry < ClusterSize_ * ClusterSize_; icry++) {
577  if (Xtals5x5[icry].subdetId() != 1)
578  continue;
579  itrechit = hits->find(Xtals5x5[icry]);
580  if (itrechit == hits->end()) {
581  std::cout << "DetId not is e25" << std::endl;
582  continue;
583  }
584 
585  if (edm::isNotFinite(itrechit->energy()))
586  return;
587  energy.push_back(itrechit->energy());
588  energy5x5 += energy[icry];
589 
590  if (icry == 6 || icry == 7 || icry == 8 || icry == 11 || icry == 12 || icry == 13 || icry == 16 || icry == 17 ||
591  icry == 18) {
592  energy3x3 += energy[icry];
593  }
594  }
595  if ((int)energy.size() != ClusterSize_ * ClusterSize_)
596  return;
597  //Once we have the matrix 5x5, we have to correct for gaps/cracks/umbrella and maincontainement
598 
599  GeneralMap->Fill(maxCC_Eta, maxCC_Phi);
600 
601  EoP_all->Fill(highPtElectron.eSuperClusterOverP());
602 
603  if (highPtElectron.classification() == elecclass_ || elecclass_ == -1) {
604  float Ptrack_in =
605  sqrt(pow(highPtElectron.trackMomentumAtVtx().X(), 2) + pow(highPtElectron.trackMomentumAtVtx().Y(), 2) +
606  pow(highPtElectron.trackMomentumAtVtx().Z(), 2));
607 
608  float UncorrectedPatCalo =
609  sqrt(pow(highPtElectron.trackMomentumAtCalo().X(), 2) + pow(highPtElectron.trackMomentumAtCalo().Y(), 2) +
610  pow(highPtElectron.trackMomentumAtCalo().Z(), 2));
611 
612  float Ptrack_out =
613  sqrt(pow(highPtElectron.trackMomentumOut().X(), 2) + pow(highPtElectron.trackMomentumOut().Y(), 2) +
614  pow(highPtElectron.trackMomentumOut().Z(), 2));
615 
616  EventMatrixNoCuts.push_back(energy);
617  EnergyVectorNoCuts.push_back(UncorrectedPatCalo);
618 
619  MaxCCetaNoCuts.push_back(maxCC_Eta);
620  MaxCCphiNoCuts.push_back(maxCC_Phi);
621 
622  WeightVectorNoCuts.push_back(energy5x5 / UncorrectedPatCalo);
623 
624  //---------------------------------------------------No Cuts-------------------------------------------------------
625  e9NoCuts->Fill(energy3x3);
626  e25NoCuts->Fill(energy5x5);
627  e9Overe25NoCuts->Fill(energy3x3 / energy5x5);
628  scENoCuts->Fill(sc.energy());
629 
630  trPNoCuts->Fill(UncorrectedPatCalo);
631 
632  EoPNoCuts->Fill(highPtElectron.eSuperClusterOverP());
633  e25OverScENoCuts->Fill(energy5x5 / sc.energy());
634 
635  E25oPNoCuts->Fill(energy5x5 / UncorrectedPatCalo);
636 
637  MapNoCuts->Fill(maxCC_Eta, maxCC_Phi);
638  PinOverPoutNoCuts->Fill(
639  sqrt(pow(highPtElectron.trackMomentumAtVtx().X(), 2) + pow(highPtElectron.trackMomentumAtVtx().Y(), 2) +
640  pow(highPtElectron.trackMomentumAtVtx().Z(), 2)) /
641  sqrt(pow(highPtElectron.trackMomentumOut().X(), 2) + pow(highPtElectron.trackMomentumOut().Y(), 2) +
642  pow(highPtElectron.trackMomentumOut().Z(), 2)));
643  eSeedOverPoutNoCuts->Fill(highPtElectron.eSuperClusterOverP());
644 
645  MapCor1NoCuts->Fill(energy5x5 / UncorrectedPatCalo, energy5x5 / Ptrack_in);
646  MapCor2NoCuts->Fill(energy5x5 / UncorrectedPatCalo, highPtElectron.eSuperClusterOverP());
647  MapCor3NoCuts->Fill(energy5x5 / UncorrectedPatCalo, Ptrack_out / Ptrack_in);
648  MapCor4NoCuts->Fill(energy5x5 / UncorrectedPatCalo, energy5x5 / highPtElectron.p());
649  MapCor5NoCuts->Fill(energy5x5 / UncorrectedPatCalo, UncorrectedPatCalo / Ptrack_out);
650  MapCor6NoCuts->Fill(Ptrack_out / Ptrack_in, energy5x5 / Ptrack_in);
651  MapCor7NoCuts->Fill(Ptrack_out / Ptrack_in, UncorrectedPatCalo / Ptrack_out);
652  MapCor8NoCuts->Fill(energy5x5 / Ptrack_in, UncorrectedPatCalo / Ptrack_out);
653  MapCor9NoCuts->Fill(energy5x5 / UncorrectedPatCalo, highPtElectron.eSeedClusterOverPout());
654  MapCor10NoCuts->Fill(highPtElectron.eSeedClusterOverPout(), Ptrack_out / Ptrack_in);
655  MapCor11NoCuts->Fill(highPtElectron.eSeedClusterOverPout(), energy5x5 / Ptrack_in);
656 
657  PinMinPoutNoCuts->Fill((Ptrack_in - Ptrack_out) / Ptrack_in);
658 
659  Error1NoCuts->Fill(highPtElectron.trackMomentumError() / Ptrack_in);
660  Error2NoCuts->Fill(highPtElectron.trackMomentumError() / Ptrack_out);
661  Error3NoCuts->Fill(highPtElectron.trackMomentumError() / UncorrectedPatCalo);
662  eSeedOverPout2NoCuts->Fill(highPtElectron.eSeedClusterOverPout());
663 
664  hadOverEmNoCuts->Fill(highPtElectron.hadronicOverEm());
665 
666  //------------------------------------------------Cuts-----------------------------------------------------
667  //Cuts!
668  if ((energy3x3 / energy5x5) < cut1_)
669  return;
670 
671  if ((Ptrack_out / Ptrack_in) < cut2_ || (Ptrack_out / Ptrack_in) > cut3_)
672  return;
673  if ((energy5x5 / Ptrack_in) < cutEPin1_ || (energy5x5 / Ptrack_in) > cutEPin2_)
674  return;
675 
676  e9->Fill(energy3x3);
677  e25->Fill(energy5x5);
678  e9Overe25->Fill(energy3x3 / energy5x5);
679  scE->Fill(sc.energy());
680  trP->Fill(UncorrectedPatCalo);
681 
682  EoP->Fill(highPtElectron.eSuperClusterOverP());
683  e25OverScE->Fill(energy5x5 / sc.energy());
684 
685  E25oP->Fill(energy5x5 / UncorrectedPatCalo);
686 
687  Map->Fill(maxCC_Eta, maxCC_Phi);
688  PinOverPout->Fill(
689  sqrt(pow(highPtElectron.trackMomentumAtVtx().X(), 2) + pow(highPtElectron.trackMomentumAtVtx().Y(), 2) +
690  pow(highPtElectron.trackMomentumAtVtx().Z(), 2)) /
691  sqrt(pow(highPtElectron.trackMomentumOut().X(), 2) + pow(highPtElectron.trackMomentumOut().Y(), 2) +
692  pow(highPtElectron.trackMomentumOut().Z(), 2)));
693  eSeedOverPout->Fill(highPtElectron.eSuperClusterOverP());
694 
695  MapCor1->Fill(energy5x5 / UncorrectedPatCalo, energy5x5 / Ptrack_in);
696  MapCor2->Fill(energy5x5 / UncorrectedPatCalo, highPtElectron.eSuperClusterOverP());
697  MapCor3->Fill(energy5x5 / UncorrectedPatCalo, Ptrack_out / Ptrack_in);
698  MapCor4->Fill(energy5x5 / UncorrectedPatCalo, energy5x5 / highPtElectron.p());
699  MapCor5->Fill(energy5x5 / UncorrectedPatCalo, UncorrectedPatCalo / Ptrack_out);
700  MapCor6->Fill(Ptrack_out / Ptrack_in, energy5x5 / Ptrack_in);
701  MapCor7->Fill(Ptrack_out / Ptrack_in, UncorrectedPatCalo / Ptrack_out);
702  MapCor8->Fill(energy5x5 / Ptrack_in, UncorrectedPatCalo / Ptrack_out);
703  MapCor9->Fill(energy5x5 / UncorrectedPatCalo, highPtElectron.eSeedClusterOverPout());
704  MapCor10->Fill(highPtElectron.eSeedClusterOverPout(), Ptrack_out / Ptrack_in);
705  MapCor11->Fill(highPtElectron.eSeedClusterOverPout(), energy5x5 / Ptrack_in);
706 
707  PinMinPout->Fill((Ptrack_in - Ptrack_out) / Ptrack_in);
708 
709  Error1->Fill(highPtElectron.trackMomentumError() / Ptrack_in);
710  Error2->Fill(highPtElectron.trackMomentumError() / Ptrack_out);
711  Error3->Fill(highPtElectron.trackMomentumError() / UncorrectedPatCalo);
712 
713  eSeedOverPout2->Fill(highPtElectron.eSeedClusterOverPout());
714  hadOverEm->Fill(highPtElectron.hadronicOverEm());
715 
716  EventMatrix.push_back(energy);
717  EnergyVector.push_back(UncorrectedPatCalo);
718  MaxCCeta.push_back(maxCC_Eta);
719  MaxCCphi.push_back(maxCC_Phi);
720 
721  WeightVector.push_back(energy5x5 / UncorrectedPatCalo);
722 
723  //-------------------------------------------------------Extra Cut-----------------------------------------------------
724  if (highPtElectron.eSeedClusterOverPout() < cutESeed_)
725  return;
726 
727  MapCor1ESeed->Fill(energy5x5 / UncorrectedPatCalo, energy5x5 / Ptrack_in);
728  MapCor2ESeed->Fill(energy5x5 / UncorrectedPatCalo, highPtElectron.eSuperClusterOverP());
729  MapCor3ESeed->Fill(energy5x5 / UncorrectedPatCalo, Ptrack_out / Ptrack_in);
730  MapCor4ESeed->Fill(energy5x5 / UncorrectedPatCalo, energy5x5 / highPtElectron.p());
731  MapCor5ESeed->Fill(energy5x5 / UncorrectedPatCalo, UncorrectedPatCalo / Ptrack_out);
732  MapCor6ESeed->Fill(Ptrack_out / Ptrack_in, energy5x5 / Ptrack_in);
733  MapCor7ESeed->Fill(Ptrack_out / Ptrack_in, UncorrectedPatCalo / Ptrack_out);
734  MapCor8ESeed->Fill(energy5x5 / Ptrack_in, UncorrectedPatCalo / Ptrack_out);
735  MapCor9ESeed->Fill(energy5x5 / UncorrectedPatCalo, highPtElectron.eSeedClusterOverPout());
736  MapCor10ESeed->Fill(highPtElectron.eSeedClusterOverPout(), Ptrack_out / Ptrack_in);
737  MapCor11ESeed->Fill(highPtElectron.eSeedClusterOverPout(), energy5x5 / Ptrack_in);
738 
739  eSeedOverPout2ESeed->Fill(highPtElectron.eSeedClusterOverPout());
740 
741  hadOverEmESeed->Fill(highPtElectron.hadronicOverEm());
742 
743  } else {
744  return;
745  }
746 }
CalibrationCluster::CalibMap ReducedMap
constexpr bool isNotFinite(T x)
Definition: isFinite.h:9
std::vector< float > WeightVector
edm::InputTag recHitLabel_
void analyze(const edm::Event &, const edm::EventSetup &) override
HouseholderDecomposition * MyHH
std::vector< EcalRecHit >::const_iterator const_iterator
CalibMap getMap(int, int, int, int)
std::vector< float > solutionNoCuts
MinL3Algorithm * MyL3Algo1
void beginJob() override
std::vector< int > MaxCCphiNoCuts
std::vector< int > MaxCCetaNoCuts
int ii
Definition: cuy.py:589
std::vector< GsfElectron > GsfElectronCollection
collection of GsfElectron objects
std::vector< float > solution
edm::InputTag electronLabel_
std::vector< float > newCalibs
std::vector< float > EnergyVector
int iEvent
Definition: GenABIO.cc:224
ElectronCalibration(const edm::ParameterSet &)
EBDetId findMaxHit(edm::Handle< EBRecHitCollection > &)
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)
std::vector< int > MaxCCphi
T sqrt(T t)
Definition: SSEVec.h:19
if(conf_.getParameter< bool >("UseStripCablingDB"))
std::vector< float > WeightVectorNoCuts
std::vector< std::vector< float > > EventMatrixNoCuts
CalibrationCluster calibCluster
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:500
static const int ETAPHIMODE
Definition: EBDetId.h:158
std::vector< EcalRecHit >::iterator iterator
const_iterator end() const
std::vector< std::vector< float > > EventMatrix
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)
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
EBDetId findMaxHit2(const std::vector< DetId > &, const EBRecHitCollection *)
void writeLine(EBDetId const &, float)
iterator find(key_type k)
list save
Definition: cuy.py:1164
std::vector< float > EnergyVectorNoCuts
tuple cout
Definition: gather_cfg.py:144
std::vector< EBDetId > get5x5Id(EBDetId const &)
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
std::vector< float > oldCalibs
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)
const_iterator begin() const
caConstants::TupleMultiplicity const *__restrict__ HitsOnGPU const *__restrict__ tindex_type *__restrict__ double *__restrict__ phits
std::vector< int > MaxCCeta