HiEgammaAlgos/src/HiPhotonType.cc

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   // Check if a given particle is isolated.
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     //      cout <<    "No mom for this particle.." << endl;
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   // Check if a given particle is isolated.
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;  //get rid of muons and neutrinos
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   // isolation cuts
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   // Check if a given particle is isolated.
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     //int apid= abs(p.pdgId());
0132     //  if(apid>11 &&  apid<20)
0133     //  continue; //get rid of muons and neutrinos
0134
0135     if (getDeltaR(p, pp) > cone)
0136       continue;
0137     double et = p.et();
0138     //  double pt = p.pt();
0139     EtCone += et;
0140   }
0141
0142   EtPhoton = pp.et();
0143
0144   // isolation cuts
0145
0146   if (EtCone - EtPhoton > 5)
0147     return kFALSE;
0148   //if(PtMaxHadron > 4.5+EtPhoton/40)
0149   //  return kFALSE;
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   // Check if a given particle is isolated.
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     //  double pt = p.pt();
0178     EtCone += et;
0179   }
0180   EtPhoton = pp.et();
0181
0182   // isolation cuts
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 }