File indexing completed on 2024-04-06 12:25:35
0001 #include <set>
0002
0003 #include "RecoJets/JetProducers/interface/JetMatchingTools.h"
0004
0005 #include "FWCore/Framework/interface/Event.h"
0006 #include "DataFormats/Common/interface/Handle.h"
0007
0008 #include "DataFormats/JetReco/interface/CaloJet.h"
0009 #include "DataFormats/JetReco/interface/GenJet.h"
0010 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0011 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0012 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0013 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0014
0015 using namespace edm;
0016
0017 namespace {
0018 template <class T>
0019 const typename T::value_type* getHit(const T& fCollection, DetId fId) {
0020 typename T::const_iterator hit = fCollection.find(fId);
0021 if (hit != fCollection.end())
0022 return &*hit;
0023 return nullptr;
0024 }
0025
0026 std::vector<const PCaloHit*> getSimHits(const PCaloHitContainer& fCollection, DetId fId) {
0027 std::vector<const PCaloHit*> result;
0028 for (unsigned i = 0; i < fCollection.size(); ++i) {
0029 if (fCollection[i].id() == fId.rawId()) {
0030 result.push_back(&(fCollection[i]));
0031 }
0032 }
0033 return result;
0034 }
0035 }
0036
0037 JetMatchingTools::JetMatchingTools(const edm::Event& fEvent, edm::ConsumesCollector&& iC)
0038 : mEvent(&fEvent),
0039 mEBRecHitCollection(nullptr),
0040 mEERecHitCollection(nullptr),
0041 mHBHERecHitCollection(nullptr),
0042 mHORecHitCollection(nullptr),
0043 mHFRecHitCollection(nullptr),
0044 mEBSimHitCollection(nullptr),
0045 mEESimHitCollection(nullptr),
0046 mHcalSimHitCollection(nullptr),
0047 mSimTrackCollection(nullptr),
0048 mSimVertexCollection(nullptr),
0049 mGenParticleCollection(nullptr) {
0050 input_ebrechits_token_ = iC.mayConsume<EBRecHitCollection>(edm::InputTag("ecalRecHit:EcalRecHitsEB"));
0051 input_eerechits_token_ = iC.mayConsume<EERecHitCollection>(edm::InputTag("ecalRecHit:EcalRecHitsEE"));
0052 input_hbherechits_token_ = iC.mayConsume<HBHERecHitCollection>(edm::InputTag("hbhereco"));
0053 input_horechits_token_ = iC.mayConsume<HORecHitCollection>(edm::InputTag("horeco"));
0054 input_hfrechits_token_ = iC.mayConsume<HFRecHitCollection>(edm::InputTag("hfreco"));
0055 input_pcalohits_ebcal_token_ = iC.mayConsume<edm::PCaloHitContainer>(edm::InputTag("g4SimHits:EcalHitsEB"));
0056 input_pcalohits_eecal_token_ = iC.mayConsume<edm::PCaloHitContainer>(edm::InputTag("g4SimHits:EcalHitsEE"));
0057 input_pcalohits_hcal_token_ = iC.mayConsume<edm::PCaloHitContainer>(edm::InputTag("g4SimHits:HcalHits"));
0058 input_simtrack_token_ = iC.mayConsume<edm::SimTrackContainer>(edm::InputTag("g4SimHits"));
0059 input_simvertex_token_ = iC.mayConsume<edm::SimVertexContainer>(edm::InputTag("g4SimHits"));
0060 input_cands_token_ = iC.mayConsume<reco::CandidateCollection>(edm::InputTag("genParticleCandidates"));
0061 }
0062
0063 JetMatchingTools::~JetMatchingTools() {}
0064
0065 const EBRecHitCollection* JetMatchingTools::getEBRecHitCollection() {
0066 if (!mEBRecHitCollection) {
0067 edm::Handle<EBRecHitCollection> recHits;
0068 mEvent->getByToken(input_ebrechits_token_, recHits);
0069 mEBRecHitCollection = &*recHits;
0070 }
0071 return mEBRecHitCollection;
0072 }
0073 const EERecHitCollection* JetMatchingTools::getEERecHitCollection() {
0074 if (!mEERecHitCollection) {
0075 edm::Handle<EERecHitCollection> recHits;
0076 mEvent->getByToken(input_eerechits_token_, recHits);
0077 mEERecHitCollection = &*recHits;
0078 }
0079 return mEERecHitCollection;
0080 }
0081 const HBHERecHitCollection* JetMatchingTools::getHBHERecHitCollection() {
0082 if (!mHBHERecHitCollection) {
0083 edm::Handle<HBHERecHitCollection> recHits;
0084 mEvent->getByToken(input_hbherechits_token_, recHits);
0085 mHBHERecHitCollection = &*recHits;
0086 }
0087 return mHBHERecHitCollection;
0088 }
0089 const HORecHitCollection* JetMatchingTools::getHORecHitCollection() {
0090 if (!mHORecHitCollection) {
0091 edm::Handle<HORecHitCollection> recHits;
0092 mEvent->getByToken(input_horechits_token_, recHits);
0093 mHORecHitCollection = &*recHits;
0094 }
0095 return mHORecHitCollection;
0096 }
0097 const HFRecHitCollection* JetMatchingTools::getHFRecHitCollection() {
0098 if (!mHFRecHitCollection) {
0099 edm::Handle<HFRecHitCollection> recHits;
0100 mEvent->getByToken(input_hfrechits_token_, recHits);
0101 mHFRecHitCollection = &*recHits;
0102 }
0103 return mHFRecHitCollection;
0104 }
0105 const PCaloHitContainer* JetMatchingTools::getEBSimHitCollection() {
0106 if (!mEBSimHitCollection) {
0107 edm::Handle<PCaloHitContainer> simHits;
0108 mEvent->getByToken(input_pcalohits_ebcal_token_, simHits);
0109 mEBSimHitCollection = &*simHits;
0110 }
0111 return mEBSimHitCollection;
0112 }
0113 const PCaloHitContainer* JetMatchingTools::getEESimHitCollection() {
0114 if (!mEESimHitCollection) {
0115 edm::Handle<PCaloHitContainer> simHits;
0116 mEvent->getByToken(input_pcalohits_eecal_token_, simHits);
0117 mEESimHitCollection = &*simHits;
0118 }
0119 return mEESimHitCollection;
0120 }
0121 const PCaloHitContainer* JetMatchingTools::getHcalSimHitCollection() {
0122 if (!mHcalSimHitCollection) {
0123 edm::Handle<PCaloHitContainer> simHits;
0124 mEvent->getByToken(input_pcalohits_hcal_token_, simHits);
0125 mHcalSimHitCollection = &*simHits;
0126 }
0127 return mHcalSimHitCollection;
0128 }
0129 const SimTrackContainer* JetMatchingTools::getSimTrackCollection() {
0130 if (!mSimTrackCollection) {
0131 edm::Handle<SimTrackContainer> simHits;
0132 mEvent->getByToken(input_simtrack_token_, simHits);
0133 mSimTrackCollection = &*simHits;
0134 }
0135 return mSimTrackCollection;
0136 }
0137 const SimVertexContainer* JetMatchingTools::getSimVertexCollection() {
0138 if (!mSimVertexCollection) {
0139 edm::Handle<SimVertexContainer> simHits;
0140 mEvent->getByToken(input_simvertex_token_, simHits);
0141 mSimVertexCollection = &*simHits;
0142 }
0143 return mSimVertexCollection;
0144 }
0145 const reco::CandidateCollection* JetMatchingTools::getGenParticlesCollection() {
0146 if (!mGenParticleCollection) {
0147 edm::Handle<reco::CandidateCollection> handle;
0148 mEvent->getByToken(input_cands_token_, handle);
0149 mGenParticleCollection = &*handle;
0150 }
0151 return mGenParticleCollection;
0152 }
0153
0154
0155 std::vector<const CaloTower*> JetMatchingTools::getConstituents(const reco::CaloJet& fJet) {
0156 std::vector<const CaloTower*> result;
0157 std::vector<CaloTowerPtr> constituents = fJet.getCaloConstituents();
0158 for (unsigned i = 0; i < constituents.size(); ++i)
0159 result.push_back(&*(constituents[i]));
0160 return result;
0161 }
0162
0163
0164 std::vector<JetMatchingTools::JetConstituent> JetMatchingTools::getConstituentHits(const CaloTower& fTower) {
0165 std::vector<JetConstituent> result;
0166
0167 for (unsigned i = 0; i < fTower.constituentsSize(); ++i) {
0168 DetId id = fTower.constituent(i);
0169
0170 if (id.det() == DetId::Ecal) {
0171 const EcalRecHit* hit = nullptr;
0172
0173 if ((EcalSubdetector)id.subdetId() == EcalBarrel) {
0174 hit = getHit(*getEBRecHitCollection(), id);
0175 } else if ((EcalSubdetector)id.subdetId() == EcalEndcap) {
0176 hit = getHit(*getEERecHitCollection(), id);
0177 }
0178
0179 assert(hit != nullptr);
0180 if (hit)
0181 result.push_back(JetConstituent(*hit));
0182 else
0183 std::cerr << "Can not find rechit for id " << id.rawId() << std::endl;
0184 } else if (id.det() == DetId::Hcal) {
0185 const CaloRecHit* hit = nullptr;
0186
0187 if ((HcalSubdetector)id.subdetId() == HcalBarrel || (HcalSubdetector)id.subdetId() == HcalEndcap) {
0188 hit = getHit(*getHBHERecHitCollection(), id);
0189 } else if ((HcalSubdetector)id.subdetId() == HcalOuter) {
0190 hit = getHit(*getHORecHitCollection(), id);
0191 }
0192 if ((HcalSubdetector)id.subdetId() == HcalForward) {
0193 hit = getHit(*getHFRecHitCollection(), id);
0194 }
0195
0196 if (hit)
0197 result.push_back(JetConstituent(*hit));
0198 else
0199 std::cerr << "Can not find rechit for id " << id.rawId() << std::endl;
0200 }
0201 }
0202
0203 return result;
0204 }
0205
0206
0207 std::vector<DetId> JetMatchingTools::getConstituentIds(const CaloTower& fTower) {
0208 std::vector<DetId> result;
0209 for (unsigned i = 0; i < fTower.constituentsSize(); ++i) {
0210 DetId id = fTower.constituent(i);
0211 result.push_back(id);
0212 }
0213 return result;
0214 }
0215
0216 std::vector<const PCaloHit*> JetMatchingTools::getPCaloHits(DetId fId) {
0217 std::vector<const PCaloHit*> result;
0218 if (fId.det() == DetId::Ecal) {
0219 if ((EcalSubdetector)fId.subdetId() == EcalBarrel) {
0220 result = getSimHits(*getEBSimHitCollection(), fId);
0221 } else if ((EcalSubdetector)fId.subdetId() == EcalEndcap) {
0222 result = getSimHits(*getEESimHitCollection(), fId);
0223 }
0224 } else if (fId.det() == DetId::Hcal) {
0225 result = getSimHits(*getHcalSimHitCollection(), fId);
0226 }
0227 return result;
0228 }
0229
0230 int JetMatchingTools::getTrackId(const PCaloHit& fHit) { return fHit.geantTrackId(); }
0231
0232 const SimTrack* JetMatchingTools::getTrack(unsigned fSimTrackId) {
0233 for (unsigned i = 0; i < getSimTrackCollection()->size(); ++i) {
0234 if ((*getSimTrackCollection())[i].trackId() == fSimTrackId)
0235 return &(*getSimTrackCollection())[i];
0236 }
0237 return nullptr;
0238 }
0239
0240 int JetMatchingTools::generatorId(unsigned fSimTrackId) {
0241 const SimTrack* track = getTrack(fSimTrackId);
0242 if (!track)
0243 return -1;
0244 while (track->noGenpart()) {
0245 if (track->noVertex()) {
0246 std::cerr << "JetMatchingTools::generatorId-> No vertex for track " << *track << std::endl;
0247 return -1;
0248 }
0249 const SimVertex* vertex = &((*getSimVertexCollection())[track->vertIndex()]);
0250 if (vertex->noParent()) {
0251 std::cerr << "JetMatchingTools::generatorId-> No track for vertex " << *vertex << std::endl;
0252 return -1;
0253 }
0254 track = getTrack(vertex->parentIndex());
0255 }
0256 return track->genpartIndex();
0257 }
0258
0259
0260 const reco::GenParticle* JetMatchingTools::getGenParticle(int fGeneratorId) {
0261 if (fGeneratorId > int(getGenParticlesCollection()->size())) {
0262 std::cerr << "JetMatchingTools::getGenParticle-> requested index " << fGeneratorId
0263 << " is grater then container size " << getGenParticlesCollection()->size() << std::endl;
0264 return nullptr;
0265 }
0266 return reco::GenJet::genParticle(
0267 &(*getGenParticlesCollection())[fGeneratorId - 1]);
0268 }
0269
0270
0271 std::vector<const reco::GenParticle*> JetMatchingTools::getGenParticles(const reco::CaloJet& fJet, bool fVerbose) {
0272 std::set<const reco::GenParticle*> result;
0273
0274 std::vector<const CaloTower*> towers = getConstituents(fJet);
0275 for (unsigned itower = 0; itower < towers.size(); ++itower) {
0276 std::vector<DetId> detids = getConstituentIds(*(towers[itower]));
0277 for (unsigned iid = 0; iid < detids.size(); ++iid) {
0278 std::vector<const PCaloHit*> phits = getPCaloHits(detids[iid]);
0279 for (unsigned iphit = 0; iphit < phits.size(); ++iphit) {
0280 int trackId = getTrackId(*(phits[iphit]));
0281 if (trackId >= 0) {
0282 int genId = generatorId(trackId);
0283 if (genId >= 0) {
0284 const reco::GenParticle* genPart = getGenParticle(genId);
0285 if (genPart) {
0286 result.insert(genPart);
0287 } else if (fVerbose) {
0288 std::cerr << "JetMatchingTools::getGenParticles-> Can not convert genId " << genId << " to GenParticle"
0289 << std::endl;
0290 }
0291 } else if (fVerbose) {
0292 std::cerr << "JetMatchingTools::getGenParticles-> Can not convert trackId " << trackId << " to genId"
0293 << std::endl;
0294 }
0295 } else if (fVerbose) {
0296 std::cerr << "JetMatchingTools::getGenParticles-> Unknown trackId for PCaloHit " << *(phits[iphit])
0297 << std::endl;
0298 }
0299 }
0300 }
0301 }
0302 return std::vector<const reco::GenParticle*>(result.begin(), result.end());
0303 }
0304
0305
0306 std::vector<const reco::GenParticle*> JetMatchingTools::getGenParticles(const reco::GenJet& fJet) {
0307 return fJet.getGenConstituents();
0308 }
0309
0310
0311 double JetMatchingTools::lostEnergyFraction(const reco::CaloJet& fJet) {
0312 double totalEnergy = 0;
0313 double lostEnergy = 0;
0314
0315 std::vector<const CaloTower*> towers = getConstituents(fJet);
0316 for (unsigned itower = 0; itower < towers.size(); ++itower) {
0317 std::vector<JetConstituent> recHits = getConstituentHits(*(towers[itower]));
0318 for (unsigned ihit = 0; ihit < recHits.size(); ++ihit) {
0319 double foundSimEnergy = 0;
0320 double lostSimEnergy = 0;
0321 std::vector<const PCaloHit*> phits = getPCaloHits(recHits[ihit].id);
0322 for (unsigned iphit = 0; iphit < phits.size(); ++iphit) {
0323 double simEnergy = phits[iphit]->energy();
0324 int trackId = getTrackId(*(phits[iphit]));
0325 if (trackId < 0 || generatorId(trackId) < 0)
0326 lostSimEnergy += simEnergy;
0327 else
0328 foundSimEnergy += simEnergy;
0329 }
0330 if (foundSimEnergy > 0 || lostSimEnergy > 0) {
0331 totalEnergy += recHits[ihit].energy;
0332 lostEnergy += recHits[ihit].energy * lostSimEnergy / (foundSimEnergy + lostSimEnergy);
0333 }
0334 }
0335 }
0336 return lostEnergy / totalEnergy;
0337 }
0338
0339
0340 double JetMatchingTools::overlapEnergyFraction(const std::vector<const reco::GenParticle*>& fObject,
0341 const std::vector<const reco::GenParticle*>& fReference) const {
0342 if (fObject.empty())
0343 return 0;
0344 double totalEnergy = 0;
0345 double overlapEnergy = 0;
0346 for (unsigned i = 0; i < fObject.size(); ++i) {
0347 totalEnergy += fObject[i]->energy();
0348 if (find(fReference.begin(), fReference.end(), fObject[i]) != fReference.end())
0349 overlapEnergy += fObject[i]->energy();
0350 }
0351 return overlapEnergy / totalEnergy;
0352 }