JetProducers/src/JetMatchingTools.cc

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 }  // namespace
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 /// get towers contributing to CaloJet
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 /// get CaloRecHits contributing to the tower
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 /// get cells contributing to the tower
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 /// get PCaloHits contributing to the detId
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 /// GEANT track ID
0230 int JetMatchingTools::getTrackId(const PCaloHit& fHit) { return fHit.geantTrackId(); }
0231 /// convert trackId to SimTrack
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 /// Generator ID
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 /// GenParticle
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]);  // knowhow: index is shifted by 1
0268 }
0269
0270 /// GenParticles for CaloJet
0271 std::vector<const reco::GenParticle*> JetMatchingTools::getGenParticles(const reco::CaloJet& fJet, bool fVerbose) {
0272   std::set<const reco::GenParticle*> result;
0273   // follow the chain
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 /// GenParticles for GenJet
0306 std::vector<const reco::GenParticle*> JetMatchingTools::getGenParticles(const reco::GenJet& fJet) {
0307   return fJet.getGenConstituents();
0308 }
0309
0310 /// energy in broken links
0311 double JetMatchingTools::lostEnergyFraction(const reco::CaloJet& fJet) {
0312   double totalEnergy = 0;
0313   double lostEnergy = 0;
0314   // follow the chain
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 /// energy overlap
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 }