CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/Calibration/IsolatedParticles/src/GenSimInfo.cc

Go to the documentation of this file.
00001 #include "Calibration/IsolatedParticles/interface/MatrixECALDetIds.h"
00002 #include "Calibration/IsolatedParticles/interface/MatrixHCALDetIds.h"
00003 #include "Calibration/IsolatedParticles/interface/ChargeIsolation.h"
00004 #include "Calibration/IsolatedParticles/interface/GenSimInfo.h"
00005 #include "Calibration/IsolatedParticles/interface/DebugInfo.h"
00006 
00007 #include<iostream>
00008 
00009 namespace spr{
00010 
00011   void eGenSimInfo(const DetId& coreDet, HepMC::GenEvent::particle_const_iterator trkItr, std::vector<spr::propagatedGenTrackID>& trackIds, const CaloGeometry* geo, const CaloTopology* caloTopology, int ieta, int iphi, genSimInfo & info, bool debug) {
00012     
00013     if (debug) std::cout << "eGenSimInfo:: For track " << (*trkItr)->momentum().rho() << "/" << (*trkItr)->momentum().eta() << "/" << (*trkItr)->momentum().phi() << " with ieta:iphi " << ieta << ":" << iphi << std::endl;
00014 
00015     std::vector<DetId> vdets = spr::matrixECALIds(coreDet, ieta, iphi, geo, caloTopology, false);
00016     if (debug) spr::debugEcalDets(0, vdets);
00017     spr::cGenSimInfo(vdets, trkItr, trackIds, true, info, debug);
00018   }
00019 
00020   void eGenSimInfo(const DetId& coreDet, HepMC::GenEvent::particle_const_iterator trkItr, std::vector<spr::propagatedGenTrackID>& trackIds, const CaloGeometry* geo, const CaloTopology* caloTopology, double dR, const GlobalVector& trackMom, spr::genSimInfo & info, bool debug) {
00021 
00022     if (debug) std::cout << "eGenSimInfo:: For track " << (*trkItr)->momentum().rho() << "/" << (*trkItr)->momentum().eta() << "/" << (*trkItr)->momentum().phi() << " with dR,tMom " << dR << " " << trackMom << std::endl;
00023 
00024     std::vector<DetId> vdets = spr::matrixECALIds(coreDet, dR, trackMom, geo, caloTopology, false);
00025     if (debug) spr::debugEcalDets(0, vdets);
00026     spr::cGenSimInfo(vdets, trkItr, trackIds, true, info, debug);
00027   }
00028 
00029   void eGenSimInfo(const DetId& coreDet, reco::GenParticleCollection::const_iterator trkItr, std::vector<spr::propagatedGenParticleID>& trackIds, const CaloGeometry* geo, const CaloTopology* caloTopology, int ieta, int iphi, genSimInfo & info, bool debug) {
00030     
00031     if (debug) std::cout << "eGenSimInfo:: For track " << trkItr->momentum().R() << "/" << trkItr->momentum().eta() << "/" << trkItr->momentum().phi() << " with ieta:iphi " << ieta << ":" << iphi << std::endl;
00032 
00033     std::vector<DetId> vdets = spr::matrixECALIds(coreDet, ieta, iphi, geo, caloTopology, false);
00034     if (debug) spr::debugEcalDets(0, vdets);
00035     spr::cGenSimInfo(vdets, trkItr, trackIds, true, info, debug);
00036   }
00037 
00038   void eGenSimInfo(const DetId& coreDet, reco::GenParticleCollection::const_iterator trkItr, std::vector<spr::propagatedGenParticleID>& trackIds, const CaloGeometry* geo, const CaloTopology* caloTopology, double dR, const GlobalVector& trackMom, spr::genSimInfo & info, bool debug) {
00039 
00040     if (debug) std::cout << "eGenSimInfo:: For track " << trkItr->momentum().R() << "/" << trkItr->momentum().eta() << "/" << trkItr->momentum().phi() << " with dR,tMom " << dR << " " << trackMom << std::endl;
00041 
00042     std::vector<DetId> vdets = spr::matrixECALIds(coreDet, dR, trackMom, geo, caloTopology, false);
00043     if (debug) spr::debugEcalDets(0, vdets);
00044     spr::cGenSimInfo(vdets, trkItr, trackIds, true, info, debug);
00045   }
00046 
00047   void hGenSimInfo(const DetId& coreDet, HepMC::GenEvent::particle_const_iterator trkItr, std::vector<spr::propagatedGenTrackID>& trackIds, const HcalTopology* topology, int ieta, int iphi, genSimInfo & info, bool includeHO, bool debug) {
00048     
00049     if (debug) std::cout << "hGenSimInfo:: For track " << (*trkItr)->momentum().rho() << "/" << (*trkItr)->momentum().eta() << "/" << (*trkItr)->momentum().phi() << " with ieta:iphi " << ieta << ":" << iphi << std::endl;
00050 
00051     std::vector<DetId> dets;
00052     dets.push_back(coreDet);
00053     std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ieta, iphi, includeHO, false);
00054     if (debug) spr::debugHcalDets(0, vdets);
00055     spr::cGenSimInfo(vdets, trkItr, trackIds, false, info, debug);
00056   }
00057 
00058   void hGenSimInfo(const DetId& coreDet, HepMC::GenEvent::particle_const_iterator trkItr, std::vector<spr::propagatedGenTrackID>& trackIds, const CaloGeometry* geo, const HcalTopology* topology, double dR, const GlobalVector& trackMom, spr::genSimInfo & info, bool includeHO, bool debug) {
00059     
00060     if (debug) std::cout << "hGenSimInfo:: For track " << (*trkItr)->momentum().rho() << "/" << (*trkItr)->momentum().eta() << "/" << (*trkItr)->momentum().phi() << " with dR,tMom " << dR << " " << trackMom << std::endl;
00061 
00062     std::vector<DetId> vdets = spr::matrixHCALIds(coreDet, geo, topology, dR, trackMom, includeHO, false);
00063     if (debug) spr::debugHcalDets(0, vdets);
00064     spr::cGenSimInfo(vdets, trkItr, trackIds, false, info, debug);
00065   }
00066 
00067   void hGenSimInfo(const DetId& coreDet, reco::GenParticleCollection::const_iterator trkItr, std::vector<spr::propagatedGenParticleID>& trackIds, const HcalTopology* topology, int ieta, int iphi, genSimInfo & info, bool includeHO, bool debug) {
00068     
00069     if (debug) std::cout << "hGenSimInfo:: For track " << trkItr->momentum().R() << "/" << trkItr->momentum().eta() << "/" << trkItr->momentum().phi() << " with ieta:iphi " << ieta << ":" << iphi << std::endl;
00070 
00071     std::vector<DetId> dets;
00072     dets.push_back(coreDet);
00073     std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ieta, iphi, includeHO, false);
00074     if (debug) spr::debugHcalDets(0, vdets);
00075     spr::cGenSimInfo(vdets, trkItr, trackIds, false, info, debug);
00076   }
00077 
00078   void hGenSimInfo(const DetId& coreDet, reco::GenParticleCollection::const_iterator trkItr, std::vector<spr::propagatedGenParticleID>& trackIds, const CaloGeometry* geo, const HcalTopology* topology, double dR, const GlobalVector& trackMom, spr::genSimInfo & info, bool includeHO, bool debug) {
00079     
00080     if (debug) std::cout << "hGenSimInfo:: For track " << trkItr->momentum().R() << "/" << trkItr->momentum().eta() << "/" << trkItr->momentum().phi() << " with dR,tMom " << dR << " " << trackMom << std::endl;
00081 
00082     std::vector<DetId> vdets = spr::matrixHCALIds(coreDet, geo, topology, dR, trackMom, includeHO, false);
00083     if (debug) spr::debugHcalDets(0, vdets);
00084     spr::cGenSimInfo(vdets, trkItr, trackIds, false, info, debug);
00085   }
00086 
00087   void cGenSimInfo(std::vector<DetId>& vdets, HepMC::GenEvent::particle_const_iterator trkItr, std::vector<spr::propagatedGenTrackID>& trackIds, bool ifECAL, spr::genSimInfo & info, bool debug) {
00088 
00089     info.maxNearP=-1.0;
00090     info.cHadronEne=info.nHadronEne=info.eleEne=info.muEne=info.photonEne=0.0;
00091     info.isChargedIso=true;
00092     for (int i=0; i<3; ++i) info.cHadronEne_[i]=0.0;
00093     for (unsigned int i=0; i<trackIds.size(); ++i) {
00094       HepMC::GenEvent::particle_const_iterator trkItr2 = trackIds[i].trkItr;
00095       // avoid the track under consideration
00096       if ( (trkItr2 != trkItr) && trackIds[i].ok) {
00097         int charge = trackIds[i].charge;
00098         int pdgid  = trackIds[i].pdgId;
00099         double p   = (*trkItr2)->momentum().rho();
00100         bool isolat= false;
00101         if (ifECAL) {
00102           const DetId anyCell = trackIds[i].detIdECAL;
00103           isolat              = spr::chargeIsolation(anyCell,vdets);
00104         } else {
00105           const DetId anyCell = trackIds[i].detIdHCAL;
00106           isolat              = spr::chargeIsolation(anyCell,vdets);
00107         }
00108         if (!isolat) spr::cGenSimInfo(charge, pdgid, p, info, debug);
00109       }
00110     }
00111     if (debug) {
00112       std::cout << "Isolation variables: isChargedIso :" << info.isChargedIso 
00113                 << " maxNearP " << info.maxNearP << " Energy e/mu/g/ch/nh "
00114                 << info.eleEne << "," << info.muEne << "," << info.photonEne
00115                 << "," << info.cHadronEne << "," << info.nHadronEne 
00116                 << " charge " << info.cHadronEne_[0] << "," 
00117                 << info.cHadronEne_[1] << "," << info.cHadronEne_[2] 
00118                 << std::endl;
00119     }
00120   }
00121 
00122   void cGenSimInfo(std::vector<DetId>& vdets, reco::GenParticleCollection::const_iterator trkItr, std::vector<spr::propagatedGenParticleID>& trackIds, bool ifECAL, spr::genSimInfo & info, bool debug) {
00123 
00124     info.maxNearP=-1.0;
00125     info.cHadronEne=info.nHadronEne=info.eleEne=info.muEne=info.photonEne=0.0;
00126     info.isChargedIso=true;
00127     for (int i=0; i<3; ++i) info.cHadronEne_[i]=0.0;
00128     for (unsigned int i=0; i<trackIds.size(); ++i) {
00129       reco::GenParticleCollection::const_iterator trkItr2 = trackIds[i].trkItr;
00130       // avoid the track under consideration
00131       if ( (trkItr2 != trkItr) && trackIds[i].ok) {
00132         int charge = trackIds[i].charge;
00133         int pdgid  = trackIds[i].pdgId;
00134         double p   = trkItr2->momentum().R();
00135         bool isolat= false;
00136         if (ifECAL) {
00137           const DetId anyCell = trackIds[i].detIdECAL;
00138           isolat              = spr::chargeIsolation(anyCell,vdets);
00139         } else {
00140           const DetId anyCell = trackIds[i].detIdHCAL;
00141           isolat              = spr::chargeIsolation(anyCell,vdets);
00142         }
00143         if (!isolat) spr::cGenSimInfo(charge, pdgid, p, info, debug);
00144       }
00145     }
00146 
00147     if (debug) {
00148       std::cout << "Isolation variables: isChargedIso :" << info.isChargedIso 
00149                 << " maxNearP " << info.maxNearP << " Energy e/mu/g/ch/nh "
00150                 << info.eleEne << "," << info.muEne << "," << info.photonEne
00151                 << "," << info.cHadronEne << "," << info.nHadronEne 
00152                 << " charge " << info.cHadronEne_[0] << "," 
00153                 << info.cHadronEne_[1] << "," << info.cHadronEne_[2] 
00154                 << std::endl;
00155     }
00156   }
00157 
00158   void cGenSimInfo(int charge, int pdgid, double p, spr::genSimInfo & info, bool debug) {
00159 
00160     if      (pdgid==22 ) info.photonEne  += p;
00161     else if (pdgid==11)  info.eleEne     += p;
00162     else if (pdgid==13)  info.muEne      += p;
00163     else if (std::abs(charge)>0)   { 
00164       info.isChargedIso = false; 
00165       info.cHadronEne  += p;
00166       if (p>1.0) info.cHadronEne_[0] += p;
00167       if (p>2.0) info.cHadronEne_[1] += p;
00168       if (p>3.0) info.cHadronEne_[2] += p;
00169       if (info.maxNearP<p) info.maxNearP=p;
00170     } else if (std::abs(charge)==0)  { 
00171       info.nHadronEne += p;
00172     }
00173   }
00174 }