CMS 3D CMS Logo

/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/RecoHI/HiEgammaAlgos/src/HiPhotonType.cc

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    // Check if a given particle is isolated.
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       //      cout <<    "No mom for this particle.." << endl;
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   // Check if a given particle is isolated.
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; //get rid of muons and neutrinos
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   // isolation cuts
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   // Check if a given particle is isolated. 
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       //int apid= abs(p.pdgId());
00156       //  if(apid>11 &&  apid<20)
00157       //  continue; //get rid of muons and neutrinos 
00158 
00159         if(getDeltaR(p,pp)>cone)
00160           continue;
00161         double et=p.et();
00162         //  double pt = p.pt();
00163         EtCone+=et;
00164 
00165     }
00166 
00167     EtPhoton = pp.et();
00168 
00169     // isolation cuts 
00170 
00171       if(EtCone-EtPhoton > 5)
00172         return kFALSE;
00173       //if(PtMaxHadron > 4.5+EtPhoton/40)
00174       //  return kFALSE;
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   // Check if a given particle is isolated. 
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       //  double pt = p.pt();           
00206       EtCone+=et;
00207 
00208     }
00209     EtPhoton = pp.et();
00210 
00211     // isolation cuts   
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