Go to the documentation of this file.00001 #include "RecoHI/HiEgammaAlgos/interface/HiPhotonType.h"
00002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00003 #include "Geometry/Records/interface/IdealGeometryRecord.h"
00004 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00005 #include "DataFormats/Common/interface/Handle.h"
00006
00007 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
00008 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
00009 #include "DataFormats/Candidate/interface/Candidate.h"
00010
00011 #include <vector>
00012
00013 using namespace edm;
00014 using namespace reco;
00015
00016 HiPhotonType::HiPhotonType(edm::Handle<GenParticleCollection> inputHandle)
00017 {
00018 using namespace std;
00019
00020 const GenParticleCollection *collection1 = inputHandle.product();
00021 mcisocut = HiGammaJetSignalDef(collection1);
00022 PI = 3.141592653589793238462643383279502884197169399375105820974945;
00023
00024 }
00025
00026
00027 bool HiPhotonType::IsIsolated(const reco::GenParticle &pp)
00028 {
00029 using namespace std;
00030 using namespace edm;
00031 using namespace reco;
00032
00033
00034 return mcisocut.IsIsolatedPP(pp);
00035 }
00036
00037
00038
00039 bool HiPhotonType::IsPrompt(const reco::GenParticle &pp)
00040 {
00041 using namespace std;
00042 if ( pp.pdgId() != 22)
00043 return false;
00044
00045 if ( pp.mother() ==0 )
00046 {
00047
00048 return false;
00049 }
00050 else
00051 {
00052 if (pp.mother()->pdgId() == 22)
00053 {
00054 cout << " found a prompt photon" << endl;
00055 return true;
00056 }
00057 else
00058 return false;
00059 }
00060 }
00061
00062
00063 HiGammaJetSignalDef::HiGammaJetSignalDef () :
00064 fSigParticles(0)
00065 { PI = 3.141592653589793238462643383279502884197169399375105820974945;
00066 }
00067 HiGammaJetSignalDef::HiGammaJetSignalDef (const reco::GenParticleCollection *sigParticles) :
00068 fSigParticles(0)
00069 {
00070 using namespace std;
00071
00072 fSigParticles = sigParticles;
00073 PI = 3.141592653589793238462643383279502884197169399375105820974945;
00074
00075 }
00076
00077 bool HiGammaJetSignalDef::IsIsolated(const reco::GenParticle &pp)
00078 {
00079 using namespace std;
00080 using namespace edm;
00081 using namespace reco;
00082
00083
00084 double EtCone = 0;
00085 double EtPhoton = 0;
00086 double PtMaxHadron = 0;
00087 double cone = 0.5;
00088 const int maxindex = (int)fSigParticles->size();
00089 for(int i=0; i < maxindex ; ++i) {
00090
00091 const GenParticle &p=(*fSigParticles)[i];
00092
00093 if(p.status()!=1)
00094 continue;
00095 if (p.collisionId() != pp.collisionId())
00096 continue;
00097
00098 int apid= abs(p.pdgId());
00099 if(apid>11 && apid<20)
00100 continue;
00101
00102 if(getDeltaR(p,pp)>cone)
00103 continue;
00104
00105 double et=p.et();
00106 double pt = p.pt();
00107 EtCone+=et;
00108 bool isHadron = false;
00109 if ( fabs(pp.pdgId())>100 && fabs(pp.pdgId()) != 310)
00110 isHadron = true;
00111
00112 if(apid>100 && apid!=310 && pt>PtMaxHadron)
00113 {
00114 if ( (isHadron == true) && (fabs(pp.pt()-pt)<0.001) && (pp.pdgId()==p.pdgId()) )
00115 continue;
00116 else
00117 PtMaxHadron=pt;
00118 }
00119
00120 }
00121 EtPhoton = pp.et();
00122
00123
00124 if(EtCone-EtPhoton > 5+EtPhoton/20)
00125 return kFALSE;
00126 if(PtMaxHadron > 4.5+EtPhoton/40)
00127 return kFALSE;
00128
00129 return kTRUE;
00130
00131 }
00132
00133 bool HiGammaJetSignalDef::IsIsolatedPP(const reco::GenParticle &pp)
00134 {
00135 using namespace std;
00136 using namespace edm;
00137 using namespace reco;
00138
00139
00140
00141 double EtCone = 0;
00142 double EtPhoton = 0;
00143 double cone = 0.4;
00144
00145 const int maxindex = (int)fSigParticles->size();
00146 for(int i=0; i < maxindex ; ++i) {
00147
00148 const GenParticle &p=(*fSigParticles)[i];
00149
00150 if(p.status()!=1)
00151 continue;
00152 if (p.collisionId() != pp.collisionId())
00153 continue;
00154
00155
00156
00157
00158
00159 if(getDeltaR(p,pp)>cone)
00160 continue;
00161 double et=p.et();
00162
00163 EtCone+=et;
00164
00165 }
00166
00167 EtPhoton = pp.et();
00168
00169
00170
00171 if(EtCone-EtPhoton > 5)
00172 return kFALSE;
00173
00174
00175
00176 return kTRUE;
00177
00178 }
00179 bool HiGammaJetSignalDef::IsIsolatedJP(const reco::GenParticle &pp)
00180 {
00181 using namespace std;
00182 using namespace edm;
00183 using namespace reco;
00184
00185
00186
00187 double EtCone = 0;
00188 double EtPhoton = 0;
00189 double cone = 0.4;
00190
00191 const int maxindex = (int)fSigParticles->size();
00192 for(int i=0; i < maxindex ; ++i) {
00193
00194 const GenParticle &p=(*fSigParticles)[i];
00195
00196 if(p.status()!=1)
00197 continue;
00198 if (p.collisionId() != pp.collisionId())
00199 continue;
00200
00201 if(getDeltaR(p,pp)>cone)
00202 continue;
00203
00204 double et=p.et();
00205
00206 EtCone+=et;
00207
00208 }
00209 EtPhoton = pp.et();
00210
00211
00212
00213 if(EtCone-EtPhoton > 5)
00214 return kFALSE;
00215 return kTRUE;
00216
00217 }
00218
00219 double HiGammaJetSignalDef::getDeltaR(const reco::Candidate &track1, const reco::Candidate &track2)
00220 {
00221 double dEta = track1.eta() - track2.eta();
00222 double dPhi = track1.phi() - track2.phi();
00223
00224 while(dPhi >= PI) dPhi -= (2.0*PI);
00225 while(dPhi < (-1.0*PI)) dPhi += (2.0*PI);
00226
00227 return sqrt(dEta*dEta+dPhi*dPhi);
00228 }
00229
00230 double HiGammaJetSignalDef::getDeltaPhi(const reco::Candidate &track1, const reco::Candidate &track2)
00231 {
00232 double dPhi = track1.phi() - track2.phi();
00233
00234 while(dPhi >= PI) dPhi -= (2.0*PI);
00235 while(dPhi < (-1.0*PI)) dPhi += (2.0*PI);
00236
00237 return dPhi;
00238 }
00239
00240
00241