File indexing completed on 2024-04-06 11:59:20
0001 #include "Calibration/IsolatedParticles/interface/MatrixECALDetIds.h"
0002 #include "Calibration/IsolatedParticles/interface/MatrixHCALDetIds.h"
0003 #include "Calibration/IsolatedParticles/interface/ChargeIsolation.h"
0004 #include "Calibration/IsolatedParticles/interface/GenSimInfo.h"
0005 #include "Calibration/IsolatedParticles/interface/DebugInfo.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007
0008 namespace spr {
0009
0010 void eGenSimInfo(const DetId& coreDet,
0011 HepMC::GenEvent::particle_const_iterator trkItr,
0012 std::vector<spr::propagatedGenTrackID>& trackIds,
0013 const CaloGeometry* geo,
0014 const CaloTopology* caloTopology,
0015 int ieta,
0016 int iphi,
0017 genSimInfo& info,
0018 bool debug) {
0019 if (debug)
0020 edm::LogVerbatim("IsoTrack") << "eGenSimInfo:: For track " << (*trkItr)->momentum().rho() << "/"
0021 << (*trkItr)->momentum().eta() << "/" << (*trkItr)->momentum().phi()
0022 << " with ieta:iphi " << ieta << ":" << iphi;
0023
0024 std::vector<DetId> vdets = spr::matrixECALIds(coreDet, ieta, iphi, geo, caloTopology, false);
0025 if (debug)
0026 spr::debugEcalDets(0, vdets);
0027 spr::cGenSimInfo(vdets, trkItr, trackIds, true, info, debug);
0028 }
0029
0030 void eGenSimInfo(const DetId& coreDet,
0031 HepMC::GenEvent::particle_const_iterator trkItr,
0032 std::vector<spr::propagatedGenTrackID>& trackIds,
0033 const CaloGeometry* geo,
0034 const CaloTopology* caloTopology,
0035 double dR,
0036 const GlobalVector& trackMom,
0037 spr::genSimInfo& info,
0038 bool debug) {
0039 if (debug)
0040 edm::LogVerbatim("IsoTrack") << "eGenSimInfo:: For track " << (*trkItr)->momentum().rho() << "/"
0041 << (*trkItr)->momentum().eta() << "/" << (*trkItr)->momentum().phi()
0042 << " with dR,tMom " << dR << " " << trackMom;
0043
0044 std::vector<DetId> vdets = spr::matrixECALIds(coreDet, dR, trackMom, geo, caloTopology, false);
0045 if (debug)
0046 spr::debugEcalDets(0, vdets);
0047 spr::cGenSimInfo(vdets, trkItr, trackIds, true, info, debug);
0048 }
0049
0050 void eGenSimInfo(const DetId& coreDet,
0051 reco::GenParticleCollection::const_iterator trkItr,
0052 std::vector<spr::propagatedGenParticleID>& trackIds,
0053 const CaloGeometry* geo,
0054 const CaloTopology* caloTopology,
0055 int ieta,
0056 int iphi,
0057 genSimInfo& info,
0058 bool debug) {
0059 if (debug)
0060 edm::LogVerbatim("IsoTrack") << "eGenSimInfo:: For track " << trkItr->momentum().R() << "/"
0061 << trkItr->momentum().eta() << "/" << trkItr->momentum().phi() << " with ieta:iphi "
0062 << ieta << ":" << iphi;
0063
0064 std::vector<DetId> vdets = spr::matrixECALIds(coreDet, ieta, iphi, geo, caloTopology, false);
0065 if (debug)
0066 spr::debugEcalDets(0, vdets);
0067 spr::cGenSimInfo(vdets, trkItr, trackIds, true, info, debug);
0068 }
0069
0070 void eGenSimInfo(const DetId& coreDet,
0071 reco::GenParticleCollection::const_iterator trkItr,
0072 std::vector<spr::propagatedGenParticleID>& trackIds,
0073 const CaloGeometry* geo,
0074 const CaloTopology* caloTopology,
0075 double dR,
0076 const GlobalVector& trackMom,
0077 spr::genSimInfo& info,
0078 bool debug) {
0079 if (debug)
0080 edm::LogVerbatim("IsoTrack") << "eGenSimInfo:: For track " << trkItr->momentum().R() << "/"
0081 << trkItr->momentum().eta() << "/" << trkItr->momentum().phi() << " with dR,tMom "
0082 << dR << " " << trackMom;
0083
0084 std::vector<DetId> vdets = spr::matrixECALIds(coreDet, dR, trackMom, geo, caloTopology, false);
0085 if (debug)
0086 spr::debugEcalDets(0, vdets);
0087 spr::cGenSimInfo(vdets, trkItr, trackIds, true, info, debug);
0088 }
0089
0090 void hGenSimInfo(const DetId& coreDet,
0091 HepMC::GenEvent::particle_const_iterator trkItr,
0092 std::vector<spr::propagatedGenTrackID>& trackIds,
0093 const HcalTopology* topology,
0094 int ieta,
0095 int iphi,
0096 genSimInfo& info,
0097 bool includeHO,
0098 bool debug) {
0099 if (debug)
0100 edm::LogVerbatim("IsoTrack") << "hGenSimInfo:: For track " << (*trkItr)->momentum().rho() << "/"
0101 << (*trkItr)->momentum().eta() << "/" << (*trkItr)->momentum().phi()
0102 << " with ieta:iphi " << ieta << ":" << iphi;
0103
0104 std::vector<DetId> dets;
0105 dets.push_back(coreDet);
0106 std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ieta, iphi, includeHO, false);
0107 if (debug)
0108 spr::debugHcalDets(0, vdets);
0109 spr::cGenSimInfo(vdets, trkItr, trackIds, false, info, debug);
0110 }
0111
0112 void hGenSimInfo(const DetId& coreDet,
0113 HepMC::GenEvent::particle_const_iterator trkItr,
0114 std::vector<spr::propagatedGenTrackID>& trackIds,
0115 const CaloGeometry* geo,
0116 const HcalTopology* topology,
0117 double dR,
0118 const GlobalVector& trackMom,
0119 spr::genSimInfo& info,
0120 bool includeHO,
0121 bool debug) {
0122 if (debug)
0123 edm::LogVerbatim("IsoTrack") << "hGenSimInfo:: For track " << (*trkItr)->momentum().rho() << "/"
0124 << (*trkItr)->momentum().eta() << "/" << (*trkItr)->momentum().phi()
0125 << " with dR,tMom " << dR << " " << trackMom;
0126
0127 std::vector<DetId> vdets = spr::matrixHCALIds(coreDet, geo, topology, dR, trackMom, includeHO, false);
0128 if (debug)
0129 spr::debugHcalDets(0, vdets);
0130 spr::cGenSimInfo(vdets, trkItr, trackIds, false, info, debug);
0131 }
0132
0133 void hGenSimInfo(const DetId& coreDet,
0134 reco::GenParticleCollection::const_iterator trkItr,
0135 std::vector<spr::propagatedGenParticleID>& trackIds,
0136 const HcalTopology* topology,
0137 int ieta,
0138 int iphi,
0139 genSimInfo& info,
0140 bool includeHO,
0141 bool debug) {
0142 if (debug)
0143 edm::LogVerbatim("IsoTrack") << "hGenSimInfo:: For track " << trkItr->momentum().R() << "/"
0144 << trkItr->momentum().eta() << "/" << trkItr->momentum().phi() << " with ieta:iphi "
0145 << ieta << ":" << iphi;
0146
0147 std::vector<DetId> dets;
0148 dets.push_back(coreDet);
0149 std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ieta, iphi, includeHO, false);
0150 if (debug)
0151 spr::debugHcalDets(0, vdets);
0152 spr::cGenSimInfo(vdets, trkItr, trackIds, false, info, debug);
0153 }
0154
0155 void hGenSimInfo(const DetId& coreDet,
0156 reco::GenParticleCollection::const_iterator trkItr,
0157 std::vector<spr::propagatedGenParticleID>& trackIds,
0158 const CaloGeometry* geo,
0159 const HcalTopology* topology,
0160 double dR,
0161 const GlobalVector& trackMom,
0162 spr::genSimInfo& info,
0163 bool includeHO,
0164 bool debug) {
0165 if (debug)
0166 edm::LogVerbatim("IsoTrack") << "hGenSimInfo:: For track " << trkItr->momentum().R() << "/"
0167 << trkItr->momentum().eta() << "/" << trkItr->momentum().phi() << " with dR,tMom "
0168 << dR << " " << trackMom;
0169
0170 std::vector<DetId> vdets = spr::matrixHCALIds(coreDet, geo, topology, dR, trackMom, includeHO, false);
0171 if (debug)
0172 spr::debugHcalDets(0, vdets);
0173 spr::cGenSimInfo(vdets, trkItr, trackIds, false, info, debug);
0174 }
0175
0176 void cGenSimInfo(std::vector<DetId>& vdets,
0177 HepMC::GenEvent::particle_const_iterator trkItr,
0178 std::vector<spr::propagatedGenTrackID>& trackIds,
0179 bool ifECAL,
0180 spr::genSimInfo& info,
0181 bool debug) {
0182 info.maxNearP = -1.0;
0183 info.cHadronEne = info.nHadronEne = info.eleEne = info.muEne = info.photonEne = 0.0;
0184 info.isChargedIso = true;
0185 for (int i = 0; i < 3; ++i)
0186 info.cHadronEne_[i] = 0.0;
0187 for (unsigned int i = 0; i < trackIds.size(); ++i) {
0188 HepMC::GenEvent::particle_const_iterator trkItr2 = trackIds[i].trkItr;
0189
0190 if ((trkItr2 != trkItr) && trackIds[i].ok) {
0191 int charge = trackIds[i].charge;
0192 int pdgid = trackIds[i].pdgId;
0193 double p = (*trkItr2)->momentum().rho();
0194 bool isolat = false;
0195 if (ifECAL) {
0196 const DetId anyCell = trackIds[i].detIdECAL;
0197 isolat = spr::chargeIsolation(anyCell, vdets);
0198 } else {
0199 const DetId anyCell = trackIds[i].detIdHCAL;
0200 isolat = spr::chargeIsolation(anyCell, vdets);
0201 }
0202 if (!isolat)
0203 spr::cGenSimInfo(charge, pdgid, p, info, debug);
0204 }
0205 }
0206
0207 if (debug) {
0208 edm::LogVerbatim("IsoTrack") << "Isolation variables: isChargedIso :" << info.isChargedIso << " maxNearP "
0209 << info.maxNearP << " Energy e/mu/g/ch/nh " << info.eleEne << "," << info.muEne
0210 << "," << info.photonEne << "," << info.cHadronEne << "," << info.nHadronEne
0211 << " charge " << info.cHadronEne_[0] << "," << info.cHadronEne_[1] << ","
0212 << info.cHadronEne_[2];
0213 }
0214 }
0215
0216 void cGenSimInfo(std::vector<DetId>& vdets,
0217 reco::GenParticleCollection::const_iterator trkItr,
0218 std::vector<spr::propagatedGenParticleID>& trackIds,
0219 bool ifECAL,
0220 spr::genSimInfo& info,
0221 bool debug) {
0222 info.maxNearP = -1.0;
0223 info.cHadronEne = info.nHadronEne = info.eleEne = info.muEne = info.photonEne = 0.0;
0224 info.isChargedIso = true;
0225 for (int i = 0; i < 3; ++i)
0226 info.cHadronEne_[i] = 0.0;
0227 for (unsigned int i = 0; i < trackIds.size(); ++i) {
0228 reco::GenParticleCollection::const_iterator trkItr2 = trackIds[i].trkItr;
0229
0230 if ((trkItr2 != trkItr) && trackIds[i].ok) {
0231 int charge = trackIds[i].charge;
0232 int pdgid = trackIds[i].pdgId;
0233 double p = trkItr2->momentum().R();
0234 bool isolat = false;
0235 if (ifECAL) {
0236 const DetId anyCell = trackIds[i].detIdECAL;
0237 isolat = spr::chargeIsolation(anyCell, vdets);
0238 } else {
0239 const DetId anyCell = trackIds[i].detIdHCAL;
0240 isolat = spr::chargeIsolation(anyCell, vdets);
0241 }
0242 if (!isolat)
0243 spr::cGenSimInfo(charge, pdgid, p, info, debug);
0244 }
0245 }
0246
0247 if (debug) {
0248 edm::LogVerbatim("IsoTrack") << "Isolation variables: isChargedIso :" << info.isChargedIso << " maxNearP "
0249 << info.maxNearP << " Energy e/mu/g/ch/nh " << info.eleEne << "," << info.muEne
0250 << "," << info.photonEne << "," << info.cHadronEne << "," << info.nHadronEne
0251 << " charge " << info.cHadronEne_[0] << "," << info.cHadronEne_[1] << ","
0252 << info.cHadronEne_[2];
0253 }
0254 }
0255
0256 void cGenSimInfo(int charge, int pdgid, double p, spr::genSimInfo& info, bool) {
0257 if (pdgid == 22)
0258 info.photonEne += p;
0259 else if (pdgid == 11)
0260 info.eleEne += p;
0261 else if (pdgid == 13)
0262 info.muEne += p;
0263 else if (std::abs(charge) > 0) {
0264 info.isChargedIso = false;
0265 info.cHadronEne += p;
0266 if (p > 1.0)
0267 info.cHadronEne_[0] += p;
0268 if (p > 2.0)
0269 info.cHadronEne_[1] += p;
0270 if (p > 3.0)
0271 info.cHadronEne_[2] += p;
0272 if (info.maxNearP < p)
0273 info.maxNearP = p;
0274 } else if (std::abs(charge) == 0) {
0275 info.nHadronEne += p;
0276 }
0277 }
0278 }