25 throw(std::runtime_error(
"\n\nEgammaObjects: Only 11 = Photon and 22 = Electron are acceptable parictleIDs! Exiting...\n\n"));
93 TH1::SetDefaultSumw2(
true);
102 =
new TH1D(
"hist_Et_",(
"Et Distribution of "+
particleString).c_str(),
105 =
new TH1D(
"hist_EtOverTruth_",(
"Reco Et over True Et of "+
particleString).c_str(),
115 =
new TH1D(
"hist_E_",(
"E Distribution of "+
particleString).c_str(),
118 =
new TH1D(
"hist_EOverTruth_",(
"Reco E over True E of "+
particleString).c_str(),
128 =
new TH1D(
"hist_Eta_",(
"Eta Distribution of "+
particleString).c_str(),
131 =
new TH1D(
"hist_EtaOverTruth_",(
"Reco Eta over True Eta of "+
particleString).c_str(),
141 =
new TH1D(
"hist_Phi_",(
"Phi Distribution of "+
particleString).c_str(),
144 =
new TH1D(
"hist_PhiOverTruth_",(
"Reco Phi over True Phi of "+
particleString).c_str(),
153 std::string recoParticleName;
155 if(
particleID == 22 ) recoParticleName =
"Higgs";
156 else if(
particleID == 11 ) recoParticleName =
"Z";
159 =
new TH1D(
"hist_All_recoMass_",(recoParticleName+
" Mass from "+
particleString+
" in All Regions").c_str(),
162 =
new TH1D(
"hist_BarrelOnly_recoMass_",(recoParticleName+
" Mass from "+
particleString+
" in Barrel").c_str(),
165 =
new TH1D(
"hist_EndcapOnly_recoMass_",(recoParticleName+
" Mass from "+
particleString+
" in EndCap").c_str(),
168 =
new TH1D(
"hist_Mixed_recoMass_",(recoParticleName+
" Mass from "+
particleString+
" in Split Detectors").c_str(),
172 =
new TH1D(
"hist_recoMass_withBackgroud_NoEtCut_",(recoParticleName+
" Mass from "+
particleString+
" with Background, No Et Cut").c_str(),
175 =
new TH1D(
"hist_recoMass_withBackgroud_5EtCut_",(recoParticleName+
" Mass from "+
particleString+
" with Background, 5 Et Cut").c_str(),
178 =
new TH1D(
"hist_recoMass_withBackgroud_10EtCut_",(recoParticleName+
" Mass from "+
particleString+
" with Background, 10 Et Cut").c_str(),
181 =
new TH1D(
"hist_recoMass_withBackgroud_20EtCut_",(recoParticleName+
" Mass from "+
particleString+
" with Background, 20 Et Cut").c_str(),
188 =
new TH2D(
"_TEMP_scatterPlot_EtOverTruthVsEt_",
"_TEMP_scatterPlot_EtOverTruthVsEt_",
191 =
new TH2D(
"_TEMP_scatterPlot_EtOverTruthVsE_",
"_TEMP_scatterPlot_EtOverTruthVsE_",
194 =
new TH2D(
"_TEMP_scatterPlot_EtOverTruthVsEta_",
"_TEMP_scatterPlot_EtOverTruthVsEta_",
197 =
new TH2D(
"_TEMP_scatterPlot_EtOverTruthVsPhi_",
"_TEMP_scatterPlot_EtOverTruthVsPhi_",
201 =
new TH2D(
"_TEMP_scatterPlot_EOverTruthVsEt_",
"_TEMP_scatterPlot_EOverTruthVsEt_",
204 =
new TH2D(
"_TEMP_scatterPlot_EOverTruthVsE_",
"_TEMP_scatterPlot_EOverTruthVsE_",
207 =
new TH2D(
"_TEMP_scatterPlot_EOverTruthVsEta_",
"_TEMP_scatterPlot_EOverTruthVsEta_",
210 =
new TH2D(
"_TEMP_scatterPlot_EOverTruthVsPhi_",
"_TEMP_scatterPlot_EOverTruthVsPhi_",
214 =
new TH2D(
"_TEMP_scatterPlot_deltaEtaVsEt_",
"_TEMP_scatterPlot_deltaEtaVsEt_",
217 =
new TH2D(
"_TEMP_scatterPlot_deltaEtaVsE_",
"_TEMP_scatterPlot_deltaEtaVsE_",
220 =
new TH2D(
"_TEMP_scatterPlot_deltaEtaVsEta_",
"_TEMP_scatterPlot_deltaEtaVsEta_",
223 =
new TH2D(
"_TEMP_scatterPlot_deltaEtaVsPhi_",
"_TEMP_scatterPlot_deltaEtaVsPhi_",
227 =
new TH2D(
"_TEMP_scatterPlot_deltaPhiVsEt_",
"_TEMP_scatterPlot_deltaPhiVsEt_",
230 =
new TH2D(
"_TEMP_scatterPlot_deltaPhiVsE_",
"_TEMP_scatterPlot_deltaPhiVsE_",
233 =
new TH2D(
"_TEMP_scatterPlot_deltaPhiVsEta_",
"_TEMP_scatterPlot_deltaPhiVsEta_",
236 =
new TH2D(
"_TEMP_scatterPlot_deltaPhiVsPhi_",
"_TEMP_scatterPlot_deltaPhiVsPhi_",
255 std::vector<reco::Photon> photonsMCMatched;
257 for(reco::PhotonCollection::const_iterator aClus = photons->begin(); aClus != photons->end(); aClus++)
259 if(aClus->et() >=
EtCut)
262 hist_E_->Fill(aClus->energy());
268 for(
int firstPhoton = 0, numPhotons = photons->size(); firstPhoton < numPhotons - 1; firstPhoton++)
269 for(
int secondPhoton = firstPhoton + 1; secondPhoton < numPhotons; secondPhoton++)
278 if(pOne.
et() >= 5 && pTwo.
et() >= 5)
281 if(pOne.
et() >= 10 && pTwo.
et() >= 10)
284 if(pOne.
et() >= 20 && pTwo.
et() >= 20)
294 const HepMC::GenEvent*
genEvent = pMCTruth->GetEvent();
296 for(HepMC::GenEvent::particle_const_iterator currentParticle = genEvent->particles_begin();
297 currentParticle != genEvent->particles_end(); currentParticle++ )
299 if(
abs((*currentParticle)->pdg_id())==22 && (*currentParticle)->status()==1
300 && (*currentParticle)->momentum().e()/cosh(
ecalEta((*currentParticle)->momentum().eta(), (*currentParticle)->production_vertex()->position().z()/10.,
301 (*currentParticle)->production_vertex()->position().perp()/10.)) >=
EtCut)
303 HepMC::FourVector vtx = (*currentParticle)->production_vertex()->position();
304 double phiTrue = (*currentParticle)->momentum().phi();
305 double etaTrue =
ecalEta((*currentParticle)->momentum().eta(), vtx.z()/10., vtx.perp()/10.);
306 double eTrue = (*currentParticle)->momentum().e();
307 double etTrue = (*currentParticle)->momentum().e()/cosh(etaTrue);
309 double etaCurrent, etaFound = -999;
310 double phiCurrent, phiFound = -999;
311 double etCurrent, etFound = -999;
312 double eCurrent, eFound = -999;
316 double closestParticleDistance = 999;
318 for(reco::PhotonCollection::const_iterator aClus = photons->begin(); aClus != photons->end(); aClus++)
320 if(aClus->et() >
EtCut)
322 etaCurrent = aClus->
eta();
323 phiCurrent = aClus->phi();
324 etCurrent = aClus->et();
325 eCurrent = aClus->energy();
327 double deltaPhi = phiCurrent-phiTrue;
332 if(deltaR < closestParticleDistance)
336 etaFound = etaCurrent;
337 phiFound = phiCurrent;
338 closestParticleDistance =
deltaR;
339 bestMatchPhoton = *aClus;
344 if(closestParticleDistance < 0.05 && etFound/etTrue > .5 && etFound/etTrue < 1.5)
380 photonsMCMatched.push_back(bestMatchPhoton);
390 if(photonsMCMatched.size() == 2)
417 std::vector<reco::GsfElectron> electronsMCMatched;
419 for(reco::GsfElectronCollection::const_iterator aClus = electrons->begin(); aClus != electrons->end(); aClus++)
421 if(aClus->et() >=
EtCut)
424 hist_E_->Fill(aClus->energy());
430 for(
int firstElectron = 0, numElectrons = electrons->size(); firstElectron < numElectrons - 1; firstElectron++)
431 for(
int secondElectron = firstElectron + 1; secondElectron < numElectrons; secondElectron++)
440 if(eOne.
et() >= 5 && eTwo.
et() >= 5)
443 if(eOne.
et() >= 10 && eTwo.
et() >= 10)
446 if(eOne.
et() >= 20 && eTwo.
et() >= 20)
456 const HepMC::GenEvent*
genEvent = pMCTruth->GetEvent();
457 for(HepMC::GenEvent::particle_const_iterator currentParticle = genEvent->particles_begin();
458 currentParticle != genEvent->particles_end(); currentParticle++ )
460 if(
abs((*currentParticle)->pdg_id())==11 && (*currentParticle)->status()==1
461 && (*currentParticle)->momentum().e()/cosh((*currentParticle)->momentum().eta()) >=
EtCut)
463 HepMC::FourVector vtx = (*currentParticle)->production_vertex()->position();
464 double phiTrue = (*currentParticle)->momentum().phi();
465 double etaTrue = (*currentParticle)->momentum().eta();
466 double eTrue = (*currentParticle)->momentum().e();
467 double etTrue = (*currentParticle)->momentum().e()/cosh(etaTrue);
469 double etaCurrent, etaFound = -999;
470 double phiCurrent, phiFound = -999;
471 double etCurrent, etFound = -999;
472 double eCurrent, eFound = -999;
476 double closestParticleDistance = 999;
478 for(reco::GsfElectronCollection::const_iterator aClus = electrons->begin(); aClus != electrons->end(); aClus++)
480 if(aClus->et() >
EtCut)
482 etaCurrent = aClus->
eta();
483 phiCurrent = aClus->phi();
484 etCurrent = aClus->et();
485 eCurrent = aClus->energy();
487 double deltaPhi = phiCurrent-phiTrue;
492 if(deltaR < closestParticleDistance)
496 etaFound = etaCurrent;
497 phiFound = phiCurrent;
498 closestParticleDistance =
deltaR;
499 bestMatchElectron = *aClus;
504 if(closestParticleDistance < 0.05 && etFound/etTrue > .5 && etFound/etTrue < 1.5)
540 electronsMCMatched.push_back(bestMatchElectron);
550 if(electronsMCMatched.size() == 2)
592 const float R_ECAL = 136.5;
596 if(EtaParticle != 0.)
599 float ZEcal = (R_ECAL-plane_Radius)*sinh(EtaParticle)+Zvertex;
601 if(ZEcal != 0.0) Theta = atan(R_ECAL/ZEcal);
602 if(Theta<0.0) Theta = Theta+
Geom::pi() ;
606 if(
std::abs(ETA) > etaBarrelEndcap )
609 if(EtaParticle<0.0 ) Zend = -Zend ;
610 float Zlen = Zend - Zvertex ;
611 float RR = Zlen/sinh(EtaParticle);
612 Theta = atan((RR+plane_Radius)/Zend);
613 if(Theta<0.0) Theta = Theta+
Geom::pi() ;
614 ETA = -
log(
tan(0.5*Theta));
621 edm::LogWarning(
"") <<
"[EgammaObjects::ecalEta] Warning: Eta equals to zero, not correcting" ;
805 hist_Et_->GetXaxis()->SetTitle(
"Et (GeV)");
806 hist_Et_->GetYaxis()->SetTitle(
"# per Et Bin");
830 hist_E_->GetXaxis()->SetTitle(
"E (GeV)");
831 hist_E_->GetYaxis()->SetTitle(
"# per E Bin");
835 hist_EEfficiency_->GetYaxis()->SetTitle(
"# True Reconstructed/# True per E Bin");
855 hist_Eta_->GetXaxis()->SetTitle(
"#eta (Radians)");
856 hist_Eta_->GetYaxis()->SetTitle(
"# per Eta Bin");
880 hist_Phi_->GetXaxis()->SetTitle(
"#phi (Radians)");
881 hist_Phi_->GetYaxis()->SetTitle(
"# per Phi Bin");
905 std::string recoParticleName;
907 if(
particleID == 22 ) recoParticleName =
"Reco Higgs";
908 else if(
particleID == 11 ) recoParticleName =
"Reco Z";
1062 std::string recoParticleName;
1064 if(
particleID == 22 ) recoParticleName =
"HiggsRecoMass";
1065 else if(
particleID == 11 ) recoParticleName =
"ZRecoMass";
TH1D * hist_resolutionEVsPhi_
TH2D * _TEMP_scatterPlot_EtOverTruthVsPhi_
double hist_min_deltaEta_
T getParameter(std::string const &) const
double hist_min_EtaOverTruth_
TH2D * _TEMP_scatterPlot_deltaEtaVsEta_
TH1D * hist_resolutionEVsE_
TH2D * _TEMP_scatterPlot_deltaPhiVsPhi_
TH1D * hist_PhiEfficiency_
TH1D * hist_deltaEtaVsEta_
float ecalEta(float EtaParticle, float Zvertex, float plane_Radius)
void analyzePhotons(const edm::Event &, const edm::EventSetup &)
TH1D * hist_EtOverTruthVsEt_
TH1D * hist_resolutionEVsEta_
reco::SuperClusterRef superCluster() const
Retrieve photonCore attributes.
virtual double et() const
transverse energy
TH1D * hist_deltaEtaVsPhi_
TH1D * hist_EOverTruthVsEt_
double deltaPhi(float phi1, float phi2)
TH1D * hist_recoMass_withBackgroud_20EtCut_
TH2D * _TEMP_scatterPlot_EtOverTruthVsE_
TH2D * _TEMP_scatterPlot_deltaPhiVsE_
void getEfficiencyHistosViaDividing()
TH1D * hist_Mixed_recoMass_
double hist_min_PhiOverTruth_
TH1D * hist_resolutionPhiVsPhi_
TH1D * hist_resolutionEtVsEta_
edm::InputTag MCTruthCollection_
TH2D * _TEMP_scatterPlot_EOverTruthVsEta_
TH2D * _TEMP_scatterPlot_deltaPhiVsEt_
SuperClusterRef superCluster() const
reference to a SuperCluster
TH1D * hist_ENumRecoOverNumTrue_
TH1D * hist_deltaPhiVsPhi_
double hist_max_EtaOverTruth_
void loadHistoParameters(const edm::ParameterSet &ps)
TH1D * hist_EOverTruthVsE_
TH1D * hist_EtNumRecoOverNumTrue_
std::vector< GsfElectron > GsfElectronCollection
collection of GsfElectron objects
TH1D * hist_resolutionEVsEt_
void createTempHistoObjects()
TH1D * hist_resolutionEtVsPhi_
std::string particleString
virtual double eta() const
momentum pseudorapidity
int hist_bins_EOverTruth_
TH1D * hist_EOverTruthVsPhi_
TH2D * _TEMP_scatterPlot_EtOverTruthVsEta_
TH1D * hist_EtOverTruthVsEta_
double hist_min_EtOverTruth_
double findRecoMass(reco::Photon pOne, reco::Photon pTwo)
TH1D * hist_EOverTruthVsEta_
void getDeltaResHistosViaSlicing()
TH1D * hist_resolutionEtVsE_
TH2D * _TEMP_scatterPlot_deltaEtaVsE_
double hist_max_recoMass_
TH1D * hist_PhiNumRecoOverNumTrue_
TH1D * hist_PhiOverTruth_
double hist_min_deltaPhi_
TH1D * hist_recoMass_withBackgroud_5EtCut_
TH2D * _TEMP_scatterPlot_deltaEtaVsPhi_
TH1D * hist_resolutionEtVsEt_
TH2D * _TEMP_scatterPlot_EOverTruthVsEt_
double hist_max_deltaEta_
Cos< T >::type cos(const T &t)
TH1D * hist_EtOverTruthVsE_
TH1D * hist_resolutionEtaVsE_
Tan< T >::type tan(const T &t)
TH1D * hist_BarrelOnly_recoMass_
TH1D * hist_EndcapOnly_recoMass_
TH1D * hist_deltaPhiVsEt_
TH1D * hist_EtaOverTruth_
TH1D * hist_deltaPhiVsEta_
TH1D * hist_deltaEtaVsEt_
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
double hist_max_deltaPhi_
int hist_bins_EtOverTruth_
TH1D * hist_recoMass_withBackgroud_10EtCut_
double deltaR(double eta1, double eta2, double phi1, double phi2)
edm::InputTag RecoCollection_
TH1D * hist_resolutionPhiVsEt_
double hist_max_PhiOverTruth_
void analyzeElectrons(const edm::Event &, const edm::EventSetup &)
TH2D * _TEMP_scatterPlot_deltaPhiVsEta_
TH2D * _TEMP_scatterPlot_EOverTruthVsE_
double hist_max_EtOverTruth_
Log< T >::type log(const T &t)
EgammaObjects(const edm::ParameterSet &)
TH1D * hist_All_recoMass_
TH2D * _TEMP_scatterPlot_deltaEtaVsEt_
math::XYZPoint caloPosition() const
TH2D * _TEMP_scatterPlot_EOverTruthVsPhi_
std::vector< Photon > PhotonCollection
collectin of Photon objects
TH1D * hist_resolutionEtaVsPhi_
static const float etaBarrelEndcap
T const * product() const
TH1D * hist_EtOverTruthVsPhi_
static const float Z_Endcap
double hist_min_recoMass_
TH1D * hist_EtaEfficiency_
TH1D * hist_resolutionPhiVsE_
double hist_min_EOverTruth_
void loadCMSSWObjects(const edm::ParameterSet &ps)
TH1D * hist_recoMass_withBackgroud_NoEtCut_
TH2D * _TEMP_scatterPlot_EtOverTruthVsEt_
int hist_bins_PhiOverTruth_
static const float R_ECAL
double hist_max_EOverTruth_
int hist_bins_EtaOverTruth_
void createBookedHistoObjects()
TH1D * hist_resolutionEtaVsEta_
TH1D * hist_resolutionPhiVsEta_
TH1D * hist_EtEfficiency_
virtual void analyze(const edm::Event &, const edm::EventSetup &)
TH1D * hist_EtaNumRecoOverNumTrue_
Power< A, B >::type pow(const A &a, const B &b)
TH1D * hist_resolutionEtaVsEt_