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