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
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
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 }