00001
00002 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00003
00004 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
00005 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00006 #include "DataFormats/DetId/interface/DetId.h"
00007 #include "CondFormats/EcalObjects/interface/EcalIntercalibConstants.h"
00008 #include "CondFormats/DataRecord/interface/EcalIntercalibConstantsRcd.h"
00009 #include "Calibration/Tools/interface/calibXMLwriter.h"
00010 #include "Calibration/Tools/interface/CalibrationCluster.h"
00011 #include "Calibration/Tools/interface/HouseholderDecomposition.h"
00012 #include "Calibration/Tools/interface/MinL3Algorithm.h"
00013 #include "Calibration/EcalCalibAlgos/interface/ElectronCalibration.h"
00014 #include "DataFormats/EgammaCandidates/interface/Electron.h"
00015 #include "DataFormats/EgammaCandidates/interface/ElectronFwd.h"
00016 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
00017 #include "DataFormats/TrackReco/interface/Track.h"
00018 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00019 #include "DataFormats/TrackReco/interface/TrackExtraFwd.h"
00020 #include "TFile.h"
00021 #include "TH1.h"
00022 #include "TH2.h"
00023 #include "TF1.h"
00024 #include "TRandom.h"
00025
00026 #include <iostream>
00027 #include <string>
00028 #include <stdexcept>
00029 #include <vector>
00030
00031 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
00032 #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
00033
00034 ElectronCalibration::ElectronCalibration(const edm::ParameterSet& iConfig)
00035 {
00036
00037 rootfile_ = iConfig.getParameter<std::string>("rootfile");
00038 recHitLabel_ = iConfig.getParameter< edm::InputTag > ("ebRecHitsLabel");
00039 electronLabel_ = iConfig.getParameter< edm::InputTag > ("electronLabel");
00040 trackLabel_ = iConfig.getParameter< edm::InputTag > ("trackLabel");
00041 calibAlgo_ = iConfig.getParameter<std::string>("CALIBRATION_ALGO");
00042 std::cout << " The used Algorithm is " << calibAlgo_ << std::endl;
00043 keventweight_ = iConfig.getParameter<int>("keventweight");
00044 ClusterSize_ = iConfig.getParameter<int>("Clustersize");
00045 ElePt_ = iConfig.getParameter<double>("ElePt");
00046 maxeta_ = iConfig.getParameter<int>("maxeta");
00047 mineta_ = iConfig.getParameter<int>("mineta");
00048 maxphi_ = iConfig.getParameter<int>("maxphi");
00049 minphi_ = iConfig.getParameter<int>("minphi");
00050 cut1_ = iConfig.getParameter<double>("cut1");
00051 cut2_ = iConfig.getParameter<double>("cut2");
00052 cut3_ = iConfig.getParameter<double>("cut3");
00053 elecclass_ = iConfig.getParameter<int>("elecclass");
00054 std::cout << " The electronclass is " << elecclass_ <<std::endl;
00055 numevent_ = iConfig.getParameter<int>("numevent");
00056 miscalibfile_ = iConfig.getParameter<std::string>("miscalibfile");
00057
00058 cutEPCalo1_ = iConfig.getParameter<double>("cutEPCaloMin");
00059 cutEPCalo2_ = iConfig.getParameter<double>("cutEPCaloMax");
00060 cutEPin1_ = iConfig.getParameter<double>("cutEPinMin");
00061 cutEPin2_ = iConfig.getParameter<double>("cutEPinMax");
00062 cutCalo1_ = iConfig.getParameter<double>("cutCaloMin");
00063 cutCalo2_ = iConfig.getParameter<double>("cutCaloMax");
00064
00065 cutESeed_ = iConfig.getParameter<double>("cutESeed");
00066 }
00067
00068
00069 ElectronCalibration::~ElectronCalibration()
00070 {
00071
00072
00073 }
00074
00075
00076 void ElectronCalibration::beginJob() {
00077
00078 f = new TFile(rootfile_.c_str(),"RECREATE");
00079
00080
00081 e9 = new TH1F("e9","E9 energy", 300, 0., 150.);
00082 e25 = new TH1F("e25","E25 energy", 300, 0., 150.);
00083 scE = new TH1F("scE","SC energy", 300, 0., 150.);
00084 trP = new TH1F("trP","Trk momentum", 300, 0., 150.);
00085 EoP = new TH1F("EoP","EoP", 600, 0., 3.);
00086 EoP_all = new TH1F("EoP_all","EoP_all",600, 0., 3.);
00087
00088 if (elecclass_ ==0 || elecclass_ == -1) {
00089 calibs = new TH1F("calib","Calibration constants", 4000, 0.5, 2.);
00090 }else{
00091 calibs = new TH1F("calib","Calibration constants", 800, 0.5, 2.);
00092 }
00093
00094 e25OverScE = new TH1F("e25OverscE","E25 / SC energy", 400, 0., 2.);
00095 E25oP = new TH1F("E25oP","E25 / P", 1200, 0., 1.5);
00096
00097 Map = new TH2F("Map","Nb Events in Crystal",85,1, 85,70 ,5, 75);
00098 e9Overe25 = new TH1F("e9Overe25","E9 / E25", 400, 0., 2.);
00099 Map3Dcalib = new TH2F("3Dcalib", "3Dcalib",85 ,1 ,85,70, 5, 75 );
00100
00101 MapCor1 = new TH2F ("MapCor1", "Correlation E25/Pcalo versus E25/Pin",100 ,0. ,5. ,100,0.,5. );
00102 MapCor2 = new TH2F ("MapCor2", "Correlation E25/Pcalo versus E/P",100 ,0. ,5. ,100,0.,5. );
00103 MapCor3 = new TH2F ("MapCor3", "Correlation E25/Pcalo versus Pout/Pin",100 ,0. ,5. ,100,0.,5. );
00104 MapCor4 = new TH2F ("MapCor4", "Correlation E25/Pcalo versus E25/highestP",100 ,0. ,5. ,100,0.,5. );
00105 MapCor5 = new TH2F ("MapCor5", "Correlation E25/Pcalo versus Pcalo/Pout",100 ,0. ,5. ,100,0.,5. );
00106 MapCor6 = new TH2F ("MapCor6", "Correlation Pout/Pin versus E25/Pin",100 ,0. ,5. ,100,0.,5. );
00107 MapCor7 = new TH2F ("MapCor7", "Correlation Pout/Pin versus Pcalo/Pout",100 ,0. ,5. ,100,0.,5. );
00108 MapCor8 = new TH2F ("MapCor8", "Correlation E25/Pin versus Pcalo/Pout",100 ,0. ,5. ,100,0.,5. );
00109 MapCor9 = new TH2F ("MapCor9", "Correlation E25/Pcalo versus Eseed/Pout",100 ,0. ,5. ,100,0.,5. );
00110 MapCor10 = new TH2F ("MapCor10", "Correlation Eseed/Pout versus Pout/Pin",100 ,0. ,5. ,100,0.,5. );
00111 MapCor11 = new TH2F ("MapCor11", "Correlation Eseed/Pout versus E25/Pin",100 ,0. ,5. ,100,0.,5. );
00112 MapCorCalib = new TH2F ("MapCorCalib", "Correlation Miscalibration versus Calibration constants", 100, 0.5,1.5, 100, 0.5, 1.5);
00113
00114 PinMinPout = new TH1F("PinMinPout","(Pin - Pout)/Pin",600,-2.0,2.0);
00115
00116 if(elecclass_ == 0 || elecclass_ == -1) {
00117 calibinter = new TH1F("calibinter", "internal calibration constants", 2000 , 0.5,2.);
00118 PinOverPout= new TH1F("PinOverPout", "pinOverpout", 600,0., 3.);
00119 eSeedOverPout= new TH1F("eSeedOverPout", "eSeedOverpout ", 600, 0., 3.);
00120 MisCalibs = new TH1F("MisCalibs","Miscalibration constants",4000,0.5,2.);
00121 RatioCalibs = new TH1F("RatioCalibs","Ratio in Calibration Constants", 4000, 0.5, 2.0);
00122 DiffCalibs = new TH1F("DiffCalibs", "Difference in Calibration constants", 4000, -1.0,1.0);
00123 }else {
00124 calibinter = new TH1F("calibinter", "internal calibration constants",400 , 0.5,2.);
00125 PinOverPout= new TH1F("PinOverPout", "pinOverpout", 600,0., 3.);
00126 eSeedOverPout= new TH1F("eSeedOverPout", "eSeedOverpout ", 600, 0., 3.);
00127 MisCalibs = new TH1F("MisCalibs","Miscalibration constants",800,0.5,2.);
00128 RatioCalibs = new TH1F("RatioCalibs","Ratio in Calibration Constants", 800, 0.5, 2.0);
00129 DiffCalibs = new TH1F("DiffCalibs", "Difference in Calibration constants", 800, -1.0,1.0);
00130 }
00131 Error1 = new TH1F ("Error1","DeltaP/Pin",800 ,-1.0,1.0 );
00132 Error2 = new TH1F ("Error2","DeltaP/Pout",800 ,-1.0,1.0 );
00133 Error3 = new TH1F ("Error3","DeltaP/Pcalo",800 ,-1.0,1.0 );
00134 eSeedOverPout2= new TH1F("eSeedOverPout2", "eSeedOverpout (No Supercluster)", 600, 0., 4.);
00135 hadOverEm= new TH1F("hadOverEm", "Had/EM distribution", 600, -2., 2.);
00136
00137
00138 Map3DcalibNoCuts = new TH2F("3DcalibNoCuts", "3Dcalib (Before Cuts)",85 ,1 ,85,70, 5, 75 );
00139 e9NoCuts = new TH1F("e9NoCuts","E9 energy (Before Cuts)",300, 0., 150.);
00140 e25NoCuts = new TH1F("e25NoCuts","E25 energy (Before Cuts)", 300, 0., 150.);
00141 scENoCuts = new TH1F("scENoCuts","SC energy (Before Cuts)", 300, 0., 150.);
00142 trPNoCuts = new TH1F("trPNoCuts","Trk momentum (Before Cuts)", 300, 0., 150.);
00143 EoPNoCuts = new TH1F("EoPNoCuts","EoP (Before Cuts)", 600, 0., 3.);
00144 if (elecclass_ ==0 || elecclass_ == -1){
00145 calibsNoCuts = new TH1F("calibNoCuts","Calibration constants (Before Cuts)", 4000, 0., 2.);
00146 }else{
00147 calibsNoCuts = new TH1F("calibNoCuts","Calibration constants (Before Cuts)", 800, 0., 2.);
00148 }
00149 e25OverScENoCuts = new TH1F("e25OverscENoCuts","E25 / SC energy (Before Cuts)", 400, 0., 2.);
00150 E25oPNoCuts = new TH1F("E25oPNoCuts","E25 / P (Before Cuts)", 1200, 0., 1.5);
00151 MapNoCuts = new TH2F("MapNoCuts","Nb Events in Crystal (Before Cuts)",85,1, 85,70 ,5, 75);
00152 e9Overe25NoCuts = new TH1F("e9Overe25NoCuts","E9 / E25 (Before Cuts)", 400, 0., 2.);
00153 PinOverPoutNoCuts = new TH1F("PinOverPoutNoCuts", "pinOverpout (Before Cuts)", 600,0., 3.);
00154 eSeedOverPoutNoCuts = new TH1F(" eSeedOverPoutNoCuts", "eSeedOverpout (Before Cuts) ", 600, 0., 4.);
00155 PinMinPoutNoCuts = new TH1F("PinMinPoutNoCuts","(Pin - Pout)/Pin (Before Cuts)",600,-2.0,2.0);
00156
00157 RatioCalibsNoCuts = new TH1F("RatioCalibsNoCuts","Ratio in Calibration Constants (Before Cuts)", 4000, 0.5, 2.0);
00158 DiffCalibsNoCuts = new TH1F("DiffCalibsNoCuts", "Difference in Calibration constants (Before Cuts)", 4000, -1.0,1.0);
00159 calibinterNoCuts = new TH1F("calibinterNoCuts", "internal calibration constants", 2000 , 0.5,2.);
00160
00161 MapCor1NoCuts = new TH2F ("MapCor1NoCuts", "Correlation E25/PatCalo versus E25/Pin (Before Cuts)",100 ,0. ,5. ,100,0.,5. );
00162 MapCor2NoCuts = new TH2F ("MapCor2NoCuts", "Correlation E25/PatCalo versus E/P (Before Cuts)",100 ,0. ,5. ,100,0.,5. );
00163 MapCor3NoCuts = new TH2F ("MapCor3NoCuts", "Correlation E25/PatCalo versus Pout/Pin (Before Cuts)",100 ,0. ,5. ,100,0.,5. );
00164 MapCor4NoCuts = new TH2F ("MapCor4NoCuts", "Correlation E25/PatCalo versus E25/highestP (Before Cuts)",100 ,0. ,5. ,100,0.,5. );
00165 MapCor5NoCuts = new TH2F ("MapCor5NoCuts", "Correlation E25/Pcalo versus Pcalo/Pout (Before Cuts)",100 ,0. ,5. ,100,0.,5. );
00166 MapCor6NoCuts = new TH2F ("MapCor6NoCuts", "Correlation Pout/Pin versus E25/Pin (Before Cuts)",100 ,0. ,5. ,100,0.,5. );
00167 MapCor7NoCuts = new TH2F ("MapCor7NoCuts", "Correlation Pout/Pin versus Pcalo/Pout (Before Cuts)",100 ,0. ,5. ,100,0.,5. );
00168 MapCor8NoCuts = new TH2F ("MapCor8NoCuts", "Correlation E25/Pin versus Pcalo/Pout (Before Cuts)",100 ,0. ,5. ,100,0.,5. );
00169 MapCor9NoCuts = new TH2F ("MapCor9NoCuts", "Correlation E25/Pcalo versus Eseed/Pout (Before Cuts)",100 ,0. ,5. ,100,0.,5. );
00170 MapCor10NoCuts = new TH2F ("MapCor10NoCuts", "Correlation Eseed/Pout versus Pout/Pin (Before Cuts)",100 ,0. ,5. ,100,0.,5. );
00171 MapCor11NoCuts = new TH2F ("MapCor11NoCuts", "Correlation Eseed/Pout versus E25/Pin (Before Cuts)",100 ,0. ,5. ,100,0.,5. );
00172 MapCorCalibNoCuts = new TH2F ("MapCorCalibNoCuts", "Correlation Miscalibration versus Calibration constants (Before Cuts)", 100, 0., 3., 100, 0., 3.);
00173
00174 Error1NoCuts = new TH1F ("Eror1NoCuts","DeltaP/Pin (Before Cuts)",800 ,-1.0,1.0 );
00175 Error2NoCuts = new TH1F ("Error2NoCuts","DeltaP/Pout (Before Cuts)",800 ,-1.0,1.0 );
00176 Error3NoCuts = new TH1F ("Error3NoCuts","DeltaP/Pcalo (Before Cuts)",800 ,-1.0, 1.0);
00177 eSeedOverPout2NoCuts= new TH1F("eSeedOverPout2NoCuts", "eSeedOverpout (No Supercluster, Before Cuts)", 600, 0., 4.);
00178 hadOverEmNoCuts= new TH1F("hadOverEmNoCuts", "Had/EM distribution (Before Cuts)", 600, -2., 2.);
00179
00180
00181 MapCor1ESeed = new TH2F ("MapCor1ESeed", "Correlation E25/Pcalo versus E25/Pin (after Eseed/Pout cut)",100 ,0. ,5. ,100,0.,5. );
00182 MapCor2ESeed = new TH2F ("MapCor2ESeed", "Correlation E25/Pcalo versus E/P (after Eseed/Pout cut)",100 ,0. ,5. ,100,0.,5. );
00183 MapCor3ESeed = new TH2F ("MapCor3ESeed", "Correlation E25/Pcalo versus Pout/Pin (after Eseed/Pout cut)",100 ,0. ,5. ,100,0.,5. );
00184 MapCor4ESeed = new TH2F ("MapCor4ESeed", "Correlation E25/Pcalo versus E25/highestP (after Eseed/Pout cut)",100 ,0. ,5. ,100,0.,5. );
00185 MapCor5ESeed = new TH2F ("MapCor5ESeed", "Correlation E25/Pcalo versus Pcalo/Pout (after Eseed/Pout cut)",100 ,0. ,5. ,100,0.,5. );
00186 MapCor6ESeed = new TH2F ("MapCor6ESeed", "Correlation Pout/Pin versus E25/Pin (after Eseed/Pout cut)",100 ,0. ,5. ,100,0.,5. );
00187 MapCor7ESeed = new TH2F ("MapCor7ESeed", "Correlation Pout/Pin versus Pcalo/Pout (after Eseed/Pout cut)",100 ,0. ,5. ,100,0.,5. );
00188 MapCor8ESeed = new TH2F ("MapCor8ESeed", "Correlation E25/Pin versus Pcalo/Pout (after Eseed/Pout cut)",100 ,0. ,5. ,100,0.,5. );
00189 MapCor9ESeed = new TH2F ("MapCor9ESeed", "Correlation E25/Pcalo versus Eseed/Pout (after Eseed/Pout cut)",100 ,0. ,5. ,100,0.,5. );
00190 MapCor10ESeed = new TH2F ("MapCor10ESeed", "Correlation Eseed/Pout versus Pout/Pin (after Eseed/Pout cut)",100 ,0. ,5. ,100,0.,5. );
00191 MapCor11ESeed = new TH2F ("MapCor11ESeed", "Correlation Eseed/Pout versus E25/Pin (after Eseed/Pout cut)",100 ,0. ,5. ,100,0.,5. );
00192
00193 eSeedOverPout2ESeed= new TH1F("eSeedOverPout2ESeed", "eSeedOverpout (No Supercluster, after Eseed/Pout cut)", 600, 0., 4.);
00194
00195 hadOverEmESeed= new TH1F("hadOverEmESeed", "Had/EM distribution (after Eseed/Pout cut)", 600, -2., 2.);
00196
00197
00198 GeneralMap = new TH2F("GeneralMap","Map without any cuts",85,1,85,70,5,75);
00199
00200 calibClusterSize=ClusterSize_;
00201 etaMin = mineta_;
00202 etaMax = maxeta_;
00203 phiMin = minphi_;
00204 phiMax = maxphi_;
00205 if(calibAlgo_=="L3") {
00206 MyL3Algo1 = new MinL3Algorithm(keventweight_,calibClusterSize, etaMin, etaMax, phiMin, phiMax);
00207 }else{
00208 if(calibAlgo_=="HH" || calibAlgo_=="HHReg"){
00209 MyHH = new HouseholderDecomposition(calibClusterSize, etaMin,etaMax, phiMin, phiMax);
00210 }else{
00211 std::cout<<" Name of Algorithm is not recognize "<<calibAlgo_<<" Should be either L3, HH or HHReg. Abort! "<<std::endl;
00212 }
00213 }
00214 read_events=0;
00215
00216
00217 ReducedMap = calibCluster.getMap(etaMin, etaMax, phiMin, phiMax);
00218
00219 oldCalibs.resize(ReducedMap.size(),0.);
00220
00221
00222 for (int phi=0; phi<360; phi++){for (int eta=0; eta<171; eta++){eventcrystal[eta][phi]=0;}}
00223
00224
00225 std::cout<<" Begin JOB "<<std::endl;
00226 }
00227
00228
00229
00230
00231 void
00232 ElectronCalibration::endJob() {
00233
00234
00235 int nIterations =10;
00236 if(calibAlgo_=="L3"){
00237 solution = MyL3Algo1->iterate(EventMatrix, MaxCCeta, MaxCCphi, EnergyVector, nIterations);
00238 }else{
00239 if(calibAlgo_=="HH"){
00240 solution = MyHH->iterate(EventMatrix, MaxCCeta, MaxCCphi,EnergyVector, 1,false);
00241 }else{
00242 if(calibAlgo_=="HHReg"){solution = MyHH->runRegional(EventMatrix, MaxCCeta, MaxCCphi,EnergyVector, 2);
00243 }else{
00244 std::cout<<" Calibration not run due to problem in Algo Choice..."<<std::endl;return;
00245 }
00246 }
00247 }
00248 for (int ii=0;ii<(int)solution.size();ii++)
00249 {
00250 std::cout << "solution[" << ii << "] = " << solution[ii] << std::endl;
00251 calibs->Fill(solution[ii]);
00252 }
00253
00254 newCalibs.resize(ReducedMap.size(),0.);
00255
00256 calibXMLwriter write_calibrations;
00257
00258 FILE* MisCalib;
00259 MisCalib = fopen(miscalibfile_.c_str(),"r");
00260 int fileStatus=1;
00261 int eta=-1;
00262 int phi=-1;
00263 float coeff=-1;
00264
00265 std::map<EBDetId,float> OldCoeff;
00266
00267 while(fileStatus != EOF) {
00268 fileStatus = fscanf(MisCalib,"%d %d %f\n", &eta,&phi,&coeff);
00269 if(eta!=-1&&phi!=-1&& coeff!=-1){
00270
00271 OldCoeff.insert(std::make_pair(EBDetId(eta,phi,EBDetId::ETAPHIMODE),coeff ));
00272 }
00273 }
00274
00275 fclose(MisCalib);
00276
00277 int icry=0;
00278 CalibrationCluster::CalibMap::iterator itmap;
00279 for (itmap=ReducedMap.begin(); itmap != ReducedMap.end();itmap++){
00280
00281 newCalibs[icry] = solution[icry];
00282
00283 write_calibrations.writeLine(itmap->first,newCalibs[icry]);
00284 float Compare =1.;
00285 std::map<EBDetId,float>::iterator iter = OldCoeff.find(itmap->first);
00286 if( iter != OldCoeff.end() )Compare = iter->second;
00287
00288 if((itmap->first).ieta()>mineta_ && (itmap->first).ieta()<maxeta_ && (itmap->first).iphi()>minphi_ && (itmap->first).iphi()<maxphi_){
00289 Map3Dcalib->Fill((itmap->first).ieta(),(itmap->first).iphi(),newCalibs[icry]*Compare ) ;
00290 MisCalibs->Fill(Compare);
00291
00292 }
00293 if((itmap->first).ieta()< mineta_+2){icry++; continue;}
00294 if((itmap->first).ieta()> maxeta_-2){icry++; continue;}
00295 if((itmap->first).iphi()< minphi_+2){icry++; continue;}
00296 if((itmap->first).iphi()> maxphi_-2){icry++; continue;}
00297
00298 calibinter->Fill(newCalibs[icry]);
00299 DiffCalibs->Fill(newCalibs[icry]-1./Compare);
00300 RatioCalibs->Fill(newCalibs[icry]*Compare);
00301 MapCorCalib->Fill(1./Compare, newCalibs[icry]);
00302 icry++;
00303 }
00304
00305 if(calibAlgo_=="L3"){
00306 solutionNoCuts = MyL3Algo1->iterate(EventMatrixNoCuts, MaxCCetaNoCuts, MaxCCphiNoCuts,EnergyVectorNoCuts,nIterations);
00307 }else{
00308 if(calibAlgo_=="HH"){
00309 solutionNoCuts = MyHH->iterate(EventMatrixNoCuts, MaxCCetaNoCuts, MaxCCphiNoCuts, EnergyVectorNoCuts, 1,false);
00310 }else{
00311 if(calibAlgo_=="HHReg"){
00312 solutionNoCuts = MyHH->runRegional(EventMatrixNoCuts, MaxCCetaNoCuts, MaxCCphiNoCuts,EnergyVectorNoCuts, 2);
00313 }else{
00314 std::cout<<" Calibration not run due to problem in AlgoChoice..."<<std::endl;return;
00315 }
00316 }
00317 }
00318 for (int ii=0;ii<(int)solutionNoCuts.size();ii++){
00319 calibsNoCuts->Fill(solutionNoCuts[ii]);
00320 }
00321 int icryp=0;
00322 CalibrationCluster::CalibMap::iterator itmapp;
00323 for (itmapp=ReducedMap.begin(); itmapp != ReducedMap.end();itmapp++){
00324
00325 newCalibs[icryp] = solutionNoCuts[icryp];
00326 float Compare2 =1.;
00327 std::map<EBDetId,float>::iterator iter2 = OldCoeff.find(itmapp->first);
00328 if( iter2 != OldCoeff.end() )Compare2 = iter2->second;
00329
00330 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) ;
00331 if ((itmapp->first).ieta()< mineta_+2){icryp++; continue;}
00332 if ((itmapp->first).ieta()> maxeta_-2){icryp++; continue;}
00333 if ((itmapp->first).iphi()< minphi_+2){icryp++; continue;}
00334 if ((itmapp->first).iphi()> maxphi_-2){icryp++; continue;}
00335 calibinterNoCuts->Fill(newCalibs[icryp]);
00336 DiffCalibsNoCuts->Fill(newCalibs[icryp]-1./(Compare2));
00337 RatioCalibsNoCuts->Fill(newCalibs[icryp]*Compare2);
00338 MapCorCalibNoCuts->Fill(1./Compare2 ,newCalibs[icryp]);
00339 icryp++;
00340 }
00341
00342
00343
00345
00346 std::cout << " " << std::endl;
00347 std::cout << "************* STATISTICS **************" << std::endl;
00348 std::cout << " Events Studied "<<read_events<< std::endl;
00349
00351
00352 f->Write();
00353
00354 f->Close();
00355 }
00356
00357
00358
00359 EBDetId ElectronCalibration::findMaxHit(edm::Handle<EBRecHitCollection> & phits) {
00360
00361
00362 EcalRecHitCollection ecrh = *phits;
00363 EcalRecHitCollection::iterator it;
00364 int count=0;
00365 EBDetId save;
00366 float en_save=0;
00367 for (it = ecrh.begin(); it != ecrh.end(); it++)
00368 {
00369 EBDetId p = EBDetId(it->id().rawId());
00370 if(it->energy()> en_save){
00371 en_save=it->energy();
00372 save=p;
00373 }
00374 count++;
00375 }
00376 return save;
00377
00378 }
00379
00380
00381 EBDetId ElectronCalibration::findMaxHit2(const std::vector<DetId> & v1,const EBRecHitCollection* hits) {
00382
00383
00384 double currEnergy = 0.;
00385 EBDetId maxHit;
00386
00387 for( std::vector<DetId>::const_iterator idsIt = v1.begin(); idsIt != v1.end(); ++idsIt) {
00388 if(idsIt->subdetId()!=1) continue;
00389 EBRecHitCollection::const_iterator itrechit;
00390 itrechit = hits->find(*idsIt);
00391
00392 if(itrechit == hits->end()){
00393 std::cout << "ElectronCalibration::findMaxHit2: rechit not found! " << std::endl;
00394 continue;
00395 }
00396 if(itrechit->energy() > currEnergy) {
00397 currEnergy=itrechit->energy();
00398 maxHit= *idsIt;
00399 }
00400 }
00401
00402 return maxHit;
00403 }
00404
00405
00406
00407 void ElectronCalibration::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup){
00408
00409 using namespace edm;
00410
00411
00412 Handle<EBRecHitCollection> phits;
00413 iEvent.getByLabel( recHitLabel_, phits);
00414 if (!phits.isValid()) {
00415 std::cerr << "Error! can't get the product EBRecHitCollection: " << std::endl;
00416 }
00417
00418 const EBRecHitCollection* hits = phits.product();
00419
00420
00421 Handle<reco::GsfElectronCollection> pElectrons;
00422
00423 iEvent.getByLabel(electronLabel_, pElectrons);
00424 if (!pElectrons.isValid()) {
00425 std::cerr << "Error! can't get the product ElectronCollection: " << std::endl;
00426 }
00427
00428 const reco::GsfElectronCollection* electronCollection = pElectrons.product();
00429 read_events++;
00430 if(read_events%1000 ==0)std::cout << "read_events = " << read_events << std::endl;
00431
00432 if(!hits)return;
00433 if(hits->size() == 0)return;
00434 if(!electronCollection)return;
00435 if(electronCollection->size() == 0)return;
00436
00437
00439
00441 reco::GsfElectronCollection::const_iterator eleIt = electronCollection->begin();
00442
00443 reco::GsfElectron highPtElectron;
00444
00445 float highestElePt=0.;
00446 bool found=false;
00447 for (eleIt=electronCollection->begin(); eleIt!=electronCollection->end(); eleIt++) {
00448
00449 if(fabs(eleIt->eta())>(maxeta_+3) * 0.0175) continue;
00450 if(eleIt->eta()<(mineta_-3) * 0.0175) continue;
00451
00452 if(eleIt->pt()>highestElePt) {
00453 highestElePt=eleIt->pt();
00454 highPtElectron = *eleIt;
00455 found =true;
00456 }
00457
00458 }
00459 if(highestElePt<ElePt_)return;
00460 if(!found) return;
00461 const reco::SuperCluster & sc = *(highPtElectron.superCluster()) ;
00462 if(fabs(sc.eta())>(maxeta_+3) * 0.0175){
00463 std::cout<<"++++ Problem with electron, electron eta is "<< highPtElectron.eta()<<" while SC is "<<sc.eta()<<std::endl;return;
00464 }
00465
00466
00467
00468 std::vector<DetId> v1;
00469
00470 for (std::vector<std::pair<DetId,float> >::const_iterator idsIt = sc.hitsAndFractions().begin();
00471 idsIt != sc.hitsAndFractions().end ();++idsIt)
00472 {v1.push_back(idsIt->first);
00473 }
00474
00475
00476 EBDetId maxHitId;
00477
00478 maxHitId = findMaxHit2(v1,hits);
00479
00480 if(maxHitId.null()){std::cout<<" Null "<<std::endl; return;}
00481
00482 int maxCC_Eta = maxHitId.ieta();
00483 int maxCC_Phi = maxHitId.iphi();
00484
00485 if(maxCC_Eta>maxeta_ )return;
00486 if(maxCC_Eta<mineta_ )return;
00487 if(maxCC_Phi>maxphi_ ) return;
00488 if(maxCC_Phi<minphi_ ) return;
00489
00490
00491 if(numevent_>0){
00492 eventcrystal[maxCC_Eta+85][maxCC_Phi-1]+=1;
00493 if (eventcrystal[maxCC_Eta+85][maxCC_Phi-1] > numevent_) return;
00494 }
00495
00496 std::vector<EBDetId> Xtals5x5 = calibCluster.get5x5Id(maxHitId);
00497
00498 if((int)Xtals5x5.size()!=ClusterSize_*ClusterSize_)return;
00499
00500
00501 std::vector<float> energy;
00502 float energy3x3=0.;
00503 float energy5x5=0.;
00504
00505 for (int icry=0;icry<ClusterSize_*ClusterSize_;icry++){
00506
00507 EBRecHitCollection::const_iterator itrechit;
00508 if(Xtals5x5[icry].subdetId()!=1) continue;
00509 itrechit = hits->find(Xtals5x5[icry]);
00510 if(itrechit==hits->end())
00511 { std::cout << "DetId not is e25" << std::endl;
00512 continue;
00513 }
00514
00515 if (std::isnan(itrechit->energy())) return;
00516 energy.push_back(itrechit->energy());
00517 energy5x5 += energy[icry];
00518
00519 if ( icry == 6 || icry == 7 || icry == 8 ||
00520 icry == 11 || icry == 12 || icry ==13 ||
00521 icry == 16 || icry == 17 || icry ==18 )
00522 {
00523 energy3x3+=energy[icry];
00524 }
00525
00526 }
00527 if((int)energy.size()!=ClusterSize_*ClusterSize_) return;
00528
00529
00530 GeneralMap->Fill(maxCC_Eta,maxCC_Phi);
00531
00532 EoP_all->Fill(highPtElectron.eSuperClusterOverP());
00533
00534 if(highPtElectron.classification()==elecclass_ || elecclass_== -1 ){
00535
00536 float Ptrack_in=sqrt( pow(highPtElectron.trackMomentumAtVtx().X(),2) +pow(highPtElectron.trackMomentumAtVtx().Y(),2) + pow(highPtElectron.trackMomentumAtVtx().Z(),2) );
00537
00538 float UncorrectedPatCalo = sqrt(pow(highPtElectron.trackMomentumAtCalo().X(),2)+pow(highPtElectron.trackMomentumAtCalo().Y(),2)+pow(highPtElectron.trackMomentumAtCalo().Z(),2));
00539
00540 float Ptrack_out = sqrt( pow(highPtElectron.trackMomentumOut().X(),2)+ pow(highPtElectron.trackMomentumOut().Y(),2)+ pow(highPtElectron.trackMomentumOut().Z(),2) );
00541
00542
00543 EventMatrixNoCuts.push_back(energy);
00544 EnergyVectorNoCuts.push_back(UncorrectedPatCalo);
00545
00546 MaxCCetaNoCuts.push_back(maxCC_Eta);
00547 MaxCCphiNoCuts.push_back(maxCC_Phi);
00548
00549 WeightVectorNoCuts.push_back(energy5x5/UncorrectedPatCalo);
00550
00551
00552 e9NoCuts->Fill(energy3x3);
00553 e25NoCuts->Fill(energy5x5);
00554 e9Overe25NoCuts->Fill(energy3x3/energy5x5);
00555 scENoCuts->Fill(sc.energy());
00556
00557 trPNoCuts->Fill(UncorrectedPatCalo);
00558
00559 EoPNoCuts->Fill(highPtElectron.eSuperClusterOverP());
00560 e25OverScENoCuts->Fill(energy5x5/sc.energy());
00561
00562 E25oPNoCuts->Fill(energy5x5/UncorrectedPatCalo);
00563
00564 MapNoCuts->Fill(maxCC_Eta,maxCC_Phi);
00565 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) ) );
00566 eSeedOverPoutNoCuts->Fill(highPtElectron.eSuperClusterOverP());
00567
00568 MapCor1NoCuts->Fill(energy5x5/UncorrectedPatCalo,energy5x5/Ptrack_in);
00569 MapCor2NoCuts->Fill(energy5x5/UncorrectedPatCalo,highPtElectron.eSuperClusterOverP());
00570 MapCor3NoCuts->Fill(energy5x5/UncorrectedPatCalo,Ptrack_out/Ptrack_in);
00571 MapCor4NoCuts->Fill(energy5x5/UncorrectedPatCalo,energy5x5/highPtElectron.p());
00572 MapCor5NoCuts->Fill(energy5x5/UncorrectedPatCalo,UncorrectedPatCalo/Ptrack_out);
00573 MapCor6NoCuts->Fill(Ptrack_out/Ptrack_in,energy5x5/Ptrack_in);
00574 MapCor7NoCuts->Fill(Ptrack_out/Ptrack_in,UncorrectedPatCalo/Ptrack_out);
00575 MapCor8NoCuts->Fill(energy5x5/Ptrack_in,UncorrectedPatCalo/Ptrack_out);
00576 MapCor9NoCuts->Fill(energy5x5/UncorrectedPatCalo,highPtElectron.eSeedClusterOverPout());
00577 MapCor10NoCuts->Fill(highPtElectron.eSeedClusterOverPout(),Ptrack_out/Ptrack_in);
00578 MapCor11NoCuts->Fill(highPtElectron.eSeedClusterOverPout(),energy5x5/Ptrack_in);
00579
00580 PinMinPoutNoCuts->Fill((Ptrack_in-Ptrack_out)/Ptrack_in);
00581
00582 Error1NoCuts-> Fill(highPtElectron.trackMomentumError()/Ptrack_in);
00583 Error2NoCuts->Fill(highPtElectron.trackMomentumError()/Ptrack_out);
00584 Error3NoCuts->Fill(highPtElectron.trackMomentumError()/UncorrectedPatCalo);
00585 eSeedOverPout2NoCuts->Fill(highPtElectron.eSeedClusterOverPout());
00586
00587 hadOverEmNoCuts->Fill(highPtElectron.hadronicOverEm());
00588
00589
00590
00591 if((energy3x3/energy5x5)<cut1_)return;
00592
00593 if((Ptrack_out/Ptrack_in)< cut2_ || (Ptrack_out/Ptrack_in)> cut3_ )return;
00594 if((energy5x5/Ptrack_in)< cutEPin1_ || (energy5x5/Ptrack_in)> cutEPin2_ )return;
00595
00596 e9->Fill(energy3x3);
00597 e25->Fill(energy5x5);
00598 e9Overe25->Fill(energy3x3/energy5x5);
00599 scE->Fill(sc.energy());
00600 trP->Fill(UncorrectedPatCalo);
00601
00602 EoP->Fill(highPtElectron.eSuperClusterOverP());
00603 e25OverScE->Fill(energy5x5/sc.energy());
00604
00605 E25oP->Fill(energy5x5/UncorrectedPatCalo);
00606
00607 Map->Fill(maxCC_Eta,maxCC_Phi);
00608 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) ) );
00609 eSeedOverPout->Fill(highPtElectron.eSuperClusterOverP());
00610
00611 MapCor1->Fill(energy5x5/UncorrectedPatCalo,energy5x5/Ptrack_in);
00612 MapCor2->Fill(energy5x5/UncorrectedPatCalo,highPtElectron.eSuperClusterOverP());
00613 MapCor3->Fill(energy5x5/UncorrectedPatCalo,Ptrack_out/Ptrack_in);
00614 MapCor4->Fill(energy5x5/UncorrectedPatCalo,energy5x5/highPtElectron.p());
00615 MapCor5->Fill(energy5x5/UncorrectedPatCalo,UncorrectedPatCalo/Ptrack_out);
00616 MapCor6->Fill(Ptrack_out/Ptrack_in,energy5x5/Ptrack_in);
00617 MapCor7->Fill(Ptrack_out/Ptrack_in,UncorrectedPatCalo/Ptrack_out);
00618 MapCor8->Fill(energy5x5/Ptrack_in,UncorrectedPatCalo/Ptrack_out);
00619 MapCor9->Fill(energy5x5/UncorrectedPatCalo,highPtElectron.eSeedClusterOverPout());
00620 MapCor10->Fill(highPtElectron.eSeedClusterOverPout(),Ptrack_out/Ptrack_in);
00621 MapCor11->Fill(highPtElectron.eSeedClusterOverPout(),energy5x5/Ptrack_in);
00622
00623 PinMinPout->Fill((Ptrack_in-Ptrack_out)/Ptrack_in);
00624
00625 Error1-> Fill(highPtElectron.trackMomentumError()/Ptrack_in);
00626 Error2->Fill(highPtElectron.trackMomentumError()/Ptrack_out);
00627 Error3->Fill(highPtElectron.trackMomentumError()/UncorrectedPatCalo);
00628
00629 eSeedOverPout2->Fill(highPtElectron.eSeedClusterOverPout());
00630 hadOverEm->Fill(highPtElectron.hadronicOverEm());
00631
00632 EventMatrix.push_back(energy);
00633 EnergyVector.push_back(UncorrectedPatCalo);
00634 MaxCCeta.push_back(maxCC_Eta);
00635 MaxCCphi.push_back(maxCC_Phi);
00636
00637 WeightVector.push_back(energy5x5/UncorrectedPatCalo);
00638
00639
00640 if(highPtElectron.eSeedClusterOverPout()< cutESeed_ ) return;
00641
00642 MapCor1ESeed->Fill(energy5x5/UncorrectedPatCalo,energy5x5/Ptrack_in);
00643 MapCor2ESeed->Fill(energy5x5/UncorrectedPatCalo,highPtElectron.eSuperClusterOverP());
00644 MapCor3ESeed->Fill(energy5x5/UncorrectedPatCalo,Ptrack_out/Ptrack_in);
00645 MapCor4ESeed->Fill(energy5x5/UncorrectedPatCalo,energy5x5/highPtElectron.p());
00646 MapCor5ESeed->Fill(energy5x5/UncorrectedPatCalo,UncorrectedPatCalo/Ptrack_out);
00647 MapCor6ESeed->Fill(Ptrack_out/Ptrack_in,energy5x5/Ptrack_in);
00648 MapCor7ESeed->Fill(Ptrack_out/Ptrack_in,UncorrectedPatCalo/Ptrack_out);
00649 MapCor8ESeed->Fill(energy5x5/Ptrack_in,UncorrectedPatCalo/Ptrack_out);
00650 MapCor9ESeed->Fill(energy5x5/UncorrectedPatCalo,highPtElectron.eSeedClusterOverPout());
00651 MapCor10ESeed->Fill(highPtElectron.eSeedClusterOverPout(),Ptrack_out/Ptrack_in);
00652 MapCor11ESeed->Fill(highPtElectron.eSeedClusterOverPout(),energy5x5/Ptrack_in);
00653
00654 eSeedOverPout2ESeed->Fill(highPtElectron.eSeedClusterOverPout());
00655
00656 hadOverEmESeed->Fill(highPtElectron.hadronicOverEm());
00657
00658 }else{return;}
00659 }
00660
1.7.3