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