File indexing completed on 2024-04-06 12:25:17
0001 #include "RecoHI/HiEgammaAlgos/interface/HiPhotonType.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0004 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0005 #include "DataFormats/Common/interface/Handle.h"
0006
0007 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0008 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
0009 #include "DataFormats/Candidate/interface/Candidate.h"
0010
0011 #include <vector>
0012
0013 using namespace edm;
0014 using namespace reco;
0015
0016 HiPhotonType::HiPhotonType(edm::Handle<GenParticleCollection> inputHandle) {
0017 using namespace std;
0018
0019 const GenParticleCollection *collection1 = inputHandle.product();
0020 mcisocut = HiGammaJetSignalDef(collection1);
0021 PI = 3.141592653589793238462643383279502884197169399375105820974945;
0022 }
0023
0024 bool HiPhotonType::IsIsolated(const reco::GenParticle &pp) {
0025 using namespace std;
0026 using namespace edm;
0027 using namespace reco;
0028
0029
0030 return mcisocut.IsIsolatedPP(pp);
0031 }
0032
0033 bool HiPhotonType::IsPrompt(const reco::GenParticle &pp) {
0034 using namespace std;
0035 if (pp.pdgId() != 22)
0036 return false;
0037
0038 if (pp.mother() == nullptr) {
0039
0040 return false;
0041 } else {
0042 if (pp.mother()->pdgId() == 22) {
0043 cout << " found a prompt photon" << endl;
0044 return true;
0045 } else
0046 return false;
0047 }
0048 }
0049
0050 HiGammaJetSignalDef::HiGammaJetSignalDef() : fSigParticles(nullptr) {
0051 PI = 3.141592653589793238462643383279502884197169399375105820974945;
0052 }
0053 HiGammaJetSignalDef::HiGammaJetSignalDef(const reco::GenParticleCollection *sigParticles) : fSigParticles(nullptr) {
0054 using namespace std;
0055
0056 fSigParticles = sigParticles;
0057 PI = 3.141592653589793238462643383279502884197169399375105820974945;
0058 }
0059
0060 bool HiGammaJetSignalDef::IsIsolated(const reco::GenParticle &pp) {
0061 using namespace std;
0062 using namespace edm;
0063 using namespace reco;
0064
0065
0066 double EtCone = 0;
0067 double EtPhoton = 0;
0068 double PtMaxHadron = 0;
0069 double cone = 0.5;
0070 const int maxindex = (int)fSigParticles->size();
0071 for (int i = 0; i < maxindex; ++i) {
0072 const GenParticle &p = (*fSigParticles)[i];
0073
0074 if (p.status() != 1)
0075 continue;
0076 if (p.collisionId() != pp.collisionId())
0077 continue;
0078
0079 int apid = abs(p.pdgId());
0080 if (apid > 11 && apid < 20)
0081 continue;
0082
0083 if (getDeltaR(p, pp) > cone)
0084 continue;
0085
0086 double et = p.et();
0087 double pt = p.pt();
0088 EtCone += et;
0089 bool isHadron = false;
0090 if (fabs(pp.pdgId()) > 100 && fabs(pp.pdgId()) != 310)
0091 isHadron = true;
0092
0093 if (apid > 100 && apid != 310 && pt > PtMaxHadron) {
0094 if ((isHadron == true) && (fabs(pp.pt() - pt) < 0.001) && (pp.pdgId() == p.pdgId()))
0095 continue;
0096 else
0097 PtMaxHadron = pt;
0098 }
0099 }
0100 EtPhoton = pp.et();
0101
0102
0103 if (EtCone - EtPhoton > 5 + EtPhoton / 20)
0104 return kFALSE;
0105 if (PtMaxHadron > 4.5 + EtPhoton / 40)
0106 return kFALSE;
0107
0108 return kTRUE;
0109 }
0110
0111 bool HiGammaJetSignalDef::IsIsolatedPP(const reco::GenParticle &pp) {
0112 using namespace std;
0113 using namespace edm;
0114 using namespace reco;
0115
0116
0117
0118 double EtCone = 0;
0119 double EtPhoton = 0;
0120 double cone = 0.4;
0121
0122 const int maxindex = (int)fSigParticles->size();
0123 for (int i = 0; i < maxindex; ++i) {
0124 const GenParticle &p = (*fSigParticles)[i];
0125
0126 if (p.status() != 1)
0127 continue;
0128 if (p.collisionId() != pp.collisionId())
0129 continue;
0130
0131
0132
0133
0134
0135 if (getDeltaR(p, pp) > cone)
0136 continue;
0137 double et = p.et();
0138
0139 EtCone += et;
0140 }
0141
0142 EtPhoton = pp.et();
0143
0144
0145
0146 if (EtCone - EtPhoton > 5)
0147 return kFALSE;
0148
0149
0150
0151 return kTRUE;
0152 }
0153 bool HiGammaJetSignalDef::IsIsolatedJP(const reco::GenParticle &pp) {
0154 using namespace std;
0155 using namespace edm;
0156 using namespace reco;
0157
0158
0159
0160 double EtCone = 0;
0161 double EtPhoton = 0;
0162 double cone = 0.4;
0163
0164 const int maxindex = (int)fSigParticles->size();
0165 for (int i = 0; i < maxindex; ++i) {
0166 const GenParticle &p = (*fSigParticles)[i];
0167
0168 if (p.status() != 1)
0169 continue;
0170 if (p.collisionId() != pp.collisionId())
0171 continue;
0172
0173 if (getDeltaR(p, pp) > cone)
0174 continue;
0175
0176 double et = p.et();
0177
0178 EtCone += et;
0179 }
0180 EtPhoton = pp.et();
0181
0182
0183
0184 if (EtCone - EtPhoton > 5)
0185 return kFALSE;
0186 return kTRUE;
0187 }
0188
0189 double HiGammaJetSignalDef::getDeltaR(const reco::Candidate &track1, const reco::Candidate &track2) {
0190 double dEta = track1.eta() - track2.eta();
0191 double dPhi = track1.phi() - track2.phi();
0192
0193 while (dPhi >= PI)
0194 dPhi -= (2.0 * PI);
0195 while (dPhi < (-1.0 * PI))
0196 dPhi += (2.0 * PI);
0197
0198 return sqrt(dEta * dEta + dPhi * dPhi);
0199 }
0200
0201 double HiGammaJetSignalDef::getDeltaPhi(const reco::Candidate &track1, const reco::Candidate &track2) {
0202 double dPhi = track1.phi() - track2.phi();
0203
0204 while (dPhi >= PI)
0205 dPhi -= (2.0 * PI);
0206 while (dPhi < (-1.0 * PI))
0207 dPhi += (2.0 * PI);
0208
0209 return dPhi;
0210 }