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