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