JetProducers/src/JetSpecific.cc

0001 ////////////////////////////////////////////////////////////////////////////////
0002 //
0003 // JetSpecific
0004 // -----------
0005 //
0006 ////////////////////////////////////////////////////////////////////////////////
0007
0008 #include "RecoJets/JetProducers/interface/JetSpecific.h"
0009
0010 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0011 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0012 #include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
0013 #include "DataFormats/RecoCandidate/interface/RecoCaloTowerCandidate.h"
0014 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0015 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0016 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0017 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0018 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0019 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0020 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0021
0022 #include <cmath>
0023
0024 ////////////////////////////////////////////////////////////////////////////////
0025 // implementation of global functions
0026 ////////////////////////////////////////////////////////////////////////////////
0027
0028 //______________________________________________________________________________
0029 // Overloaded methods to write out specific types
0030
0031 // CaloJet
0032 void reco::writeSpecific(reco::CaloJet& jet,
0033                          reco::Particle::LorentzVector const& p4,
0034                          reco::Particle::Point const& point,
0035                          std::vector<reco::CandidatePtr> const& constituents,
0036                          CaloGeometry const& geometry,
0037                          HcalTopology const& topology) {
0038   const CaloSubdetectorGeometry* towerGeometry = geometry.getSubdetectorGeometry(DetId::Calo, CaloTowerDetId::SubdetId);
0039
0040   // Make the specific
0041   reco::CaloJet::Specific specific;
0042   makeSpecific(constituents, towerGeometry, &specific, topology);
0043   // Set the calo jet
0044   jet = reco::CaloJet(p4, point, specific, constituents);
0045 }
0046
0047 // BasicJet
0048 void reco::writeSpecific(reco::BasicJet& jet,
0049                          reco::Particle::LorentzVector const& p4,
0050                          reco::Particle::Point const& point,
0051                          std::vector<reco::CandidatePtr> const& constituents) {
0052   jet = reco::BasicJet(p4, point, constituents);
0053 }
0054
0055 // GenJet
0056 void reco::writeSpecific(reco::GenJet& jet,
0057                          reco::Particle::LorentzVector const& p4,
0058                          reco::Particle::Point const& point,
0059                          std::vector<reco::CandidatePtr> const& constituents) {
0060   // Make the specific
0061   reco::GenJet::Specific specific;
0062   makeSpecific(constituents, &specific);
0063   // Set to the jet
0064   jet = reco::GenJet(p4, point, specific, constituents);
0065 }
0066
0067 // PFJet
0068 void reco::writeSpecific(reco::PFJet& jet,
0069                          reco::Particle::LorentzVector const& p4,
0070                          reco::Particle::Point const& point,
0071                          std::vector<reco::CandidatePtr> const& constituents,
0072                          edm::ValueMap<float> const* weights) {
0073   // Make the specific
0074   reco::PFJet::Specific specific;
0075   makeSpecific(constituents, &specific, weights);
0076   // now make jet charge
0077   int charge = 0.;
0078   for (std::vector<reco::CandidatePtr>::const_iterator ic = constituents.begin(), icend = constituents.end();
0079        ic != icend;
0080        ++ic) {
0081     float weight = (weights != nullptr) ? (*weights)[*ic] : 1.0;
0082     charge += (*ic)->charge() * weight;
0083   }
0084   jet = reco::PFJet(p4, point, specific, constituents);
0085   jet.setCharge(charge);
0086 }
0087
0088 // TrackJet
0089 void reco::writeSpecific(reco::TrackJet& jet,
0090                          reco::Particle::LorentzVector const& p4,
0091                          reco::Particle::Point const& point,
0092                          std::vector<reco::CandidatePtr> const& constituents) {
0093   jet = reco::TrackJet(p4, point, constituents);
0094 }
0095
0096 // PFClusterJet
0097 void reco::writeSpecific(reco::PFClusterJet& jet,
0098                          reco::Particle::LorentzVector const& p4,
0099                          reco::Particle::Point const& point,
0100                          std::vector<reco::CandidatePtr> const& constituents) {
0101   jet = reco::PFClusterJet(p4, point, constituents);
0102 }
0103
0104 //______________________________________________________________________________
0105 bool reco::makeSpecific(std::vector<reco::CandidatePtr> const& towers,
0106                         const CaloSubdetectorGeometry* towerGeometry,
0107                         CaloJet::Specific* caloJetSpecific,
0108                         const HcalTopology& topology) {
0109   if (nullptr == caloJetSpecific)
0110     return false;
0111
0112   // 1.- Loop over the tower Ids,
0113   // 2.- Get the corresponding CaloTower
0114   // 3.- Calculate the different CaloJet specific quantities
0115   std::vector<double> eECal_i;
0116   std::vector<double> eHCal_i;
0117   double eInHad = 0.;
0118   double eInEm = 0.;
0119   double eInHO = 0.;
0120   double eInHB = 0.;
0121   double eInHE = 0.;
0122   double eHadInHF = 0.;
0123   double eEmInHF = 0.;
0124   double eInEB = 0.;
0125   double eInEE = 0.;
0126   double jetArea = 0.;
0127
0128   std::vector<reco::CandidatePtr>::const_iterator itTower;
0129   for (itTower = towers.begin(); itTower != towers.end(); ++itTower) {
0130     if (itTower->isNull() || !itTower->isAvailable()) {
0131       edm::LogWarning("DataNotFound") << " JetSpecific: Tower is invalid\n";
0132       continue;
0133     }
0134     const CaloTower* tower = dynamic_cast<const CaloTower*>(itTower->get());
0135     if (tower) {
0136       //Array of energy in EM Towers:
0137       eECal_i.push_back(tower->emEnergy());
0138       eInEm += tower->emEnergy();
0139       //Array of energy in HCAL Towers:
0140       eHCal_i.push_back(tower->hadEnergy());
0141       eInHad += tower->hadEnergy();
0142
0143       //  figure out contributions
0144       switch (reco::hcalSubdetector(tower->id().ieta(), topology)) {
0145         case HcalBarrel:
0146           eInHB += tower->hadEnergy();
0147           eInHO += tower->outerEnergy();
0148           eInEB += tower->emEnergy();
0149           break;
0150         case HcalEndcap:
0151           eInHE += tower->hadEnergy();
0152           eInEE += tower->emEnergy();
0153           break;
0154         case HcalForward:
0155           eHadInHF += tower->hadEnergy();
0156           eEmInHF += tower->emEnergy();
0157           break;
0158         default:
0159           break;
0160       }
0161       // get area of the tower (++ minus --)
0162       auto geometry = towerGeometry->getGeometry(tower->id());
0163       if (geometry) {
0164         jetArea += geometry->etaSpan() * geometry->phiSpan();
0165       } else {
0166         edm::LogWarning("DataNotFound") << "reco::makeCaloJetSpecific: Geometry for cell " << tower->id()
0167                                         << " can not be found. Ignoring cell\n";
0168       }
0169     } else {
0170       edm::LogWarning("DataNotFound") << "reco::makeCaloJetSpecific: Constituent is not of "
0171                                       << "CaloTower type\n";
0172     }
0173   }
0174
0175   double towerEnergy = eInHad + eInEm;
0176   caloJetSpecific->mHadEnergyInHO = eInHO;
0177   caloJetSpecific->mHadEnergyInHB = eInHB;
0178   caloJetSpecific->mHadEnergyInHE = eInHE;
0179   caloJetSpecific->mHadEnergyInHF = eHadInHF;
0180   caloJetSpecific->mEmEnergyInHF = eEmInHF;
0181   caloJetSpecific->mEmEnergyInEB = eInEB;
0182   caloJetSpecific->mEmEnergyInEE = eInEE;
0183   if (towerEnergy > 0) {
0184     caloJetSpecific->mEnergyFractionHadronic = eInHad / towerEnergy;
0185     caloJetSpecific->mEnergyFractionEm = eInEm / towerEnergy;
0186   } else {  // HO only jet
0187     caloJetSpecific->mEnergyFractionHadronic = 1.;
0188     caloJetSpecific->mEnergyFractionEm = 0.;
0189   }
0190   caloJetSpecific->mTowersArea = jetArea;
0191   caloJetSpecific->mMaxEInEmTowers = 0;
0192   caloJetSpecific->mMaxEInHadTowers = 0;
0193
0194   //Sort the arrays
0195   sort(eECal_i.begin(), eECal_i.end(), std::greater<double>());
0196   sort(eHCal_i.begin(), eHCal_i.end(), std::greater<double>());
0197
0198   if (!towers.empty()) {
0199     //Highest value in the array is the first element of the array
0200     caloJetSpecific->mMaxEInEmTowers = eECal_i.front();
0201     caloJetSpecific->mMaxEInHadTowers = eHCal_i.front();
0202   }
0203
0204   return true;
0205 }
0206
0207 //______________________________________________________________________________
0208 bool reco::makeSpecific(std::vector<reco::CandidatePtr> const& particles,
0209                         PFJet::Specific* pfJetSpecific,
0210                         edm::ValueMap<float> const* weights) {
0211   if (nullptr == pfJetSpecific)
0212     return false;
0213
0214   // 1.- Loop over PFCandidates,
0215   // 2.- Get the corresponding PFCandidate
0216   // 3.- Calculate the different PFJet specific quantities
0217
0218   float chargedHadronEnergy = 0.;
0219   float neutralHadronEnergy = 0.;
0220   float photonEnergy = 0.;
0221   float electronEnergy = 0.;
0222   float muonEnergy = 0.;
0223   float HFHadronEnergy = 0.;
0224   float HFEMEnergy = 0.;
0225   float chargedHadronMultiplicity = 0;
0226   float neutralHadronMultiplicity = 0;
0227   float photonMultiplicity = 0;
0228   float electronMultiplicity = 0;
0229   float muonMultiplicity = 0;
0230   float HFHadronMultiplicity = 0;
0231   float HFEMMultiplicity = 0;
0232
0233   float chargedEmEnergy = 0.;
0234   float neutralEmEnergy = 0.;
0235   float chargedMuEnergy = 0.;
0236   float chargedMultiplicity = 0;
0237   float neutralMultiplicity = 0;
0238
0239   float HOEnergy = 0.;
0240
0241   std::vector<reco::CandidatePtr>::const_iterator itParticle;
0242   for (itParticle = particles.begin(); itParticle != particles.end(); ++itParticle) {
0243     if (itParticle->isNull() || !itParticle->isAvailable()) {
0244       edm::LogWarning("DataNotFound") << " JetSpecific: PF Particle is invalid\n";
0245       continue;
0246     }
0247     const Candidate* pfCand = itParticle->get();
0248     if (pfCand) {
0249       const PFCandidate* pfCandCast = dynamic_cast<const PFCandidate*>(pfCand);
0250       float weight = (weights != nullptr) ? (*weights)[*itParticle] : 1.0;
0251       if (pfCandCast)
0252         HOEnergy += pfCandCast->hoEnergy() * weight;
0253
0254       switch (std::abs(pfCand->pdgId())) {
0255         case 211:  //PFCandidate::h:       // charged hadron
0256           chargedHadronEnergy += pfCand->energy() * weight;
0257           chargedHadronMultiplicity += weight;
0258           chargedMultiplicity += weight;
0259           break;
0260
0261         case 130:  //PFCandidate::h0 :    // neutral hadron
0262           neutralHadronEnergy += pfCand->energy() * weight;
0263           neutralHadronMultiplicity += weight;
0264           neutralMultiplicity += weight;
0265           break;
0266
0267         case 22:  //PFCandidate::gamma:   // photon
0268           photonEnergy += pfCand->energy() * weight;
0269           photonMultiplicity += weight;
0270           neutralEmEnergy += pfCand->energy() * weight;
0271           neutralMultiplicity += weight;
0272           break;
0273
0274         case 11:  // PFCandidate::e:       // electron
0275           electronEnergy += pfCand->energy() * weight;
0276           electronMultiplicity += weight;
0277           chargedEmEnergy += pfCand->energy() * weight;
0278           chargedMultiplicity += weight;
0279           break;
0280
0281         case 13:  //PFCandidate::mu:      // muon
0282           muonEnergy += pfCand->energy() * weight;
0283           muonMultiplicity += weight;
0284           chargedMuEnergy += pfCand->energy() * weight;
0285           chargedMultiplicity += weight;
0286           break;
0287
0288         case 1:  // PFCandidate::h_HF :    // hadron in HF
0289           HFHadronEnergy += pfCand->energy() * weight;
0290           HFHadronMultiplicity += weight;
0291           neutralHadronEnergy += pfCand->energy() * weight;
0292           neutralMultiplicity += weight;
0293           break;
0294
0295         case 2:  //PFCandidate::egamma_HF :    // electromagnetic in HF
0296           HFEMEnergy += pfCand->energy() * weight;
0297           HFEMMultiplicity += weight;
0298           neutralEmEnergy += pfCand->energy() * weight;
0299           neutralMultiplicity += weight;
0300           break;
0301
0302         default:
0303           edm::LogWarning("DataNotFound")
0304               << "reco::makePFJetSpecific: Unknown PFCandidate::ParticleType: " << pfCand->pdgId() << " is ignored\n";
0305           break;
0306       }
0307     } else {
0308       edm::LogWarning("DataNotFound") << "reco::makePFJetSpecific: Referred constituent is not "
0309                                       << "a PFCandidate\n";
0310     }
0311   }
0312
0313   pfJetSpecific->mChargedHadronEnergy = chargedHadronEnergy;
0314   pfJetSpecific->mNeutralHadronEnergy = neutralHadronEnergy;
0315   pfJetSpecific->mPhotonEnergy = photonEnergy;
0316   pfJetSpecific->mElectronEnergy = electronEnergy;
0317   pfJetSpecific->mMuonEnergy = muonEnergy;
0318   pfJetSpecific->mHFHadronEnergy = HFHadronEnergy;
0319   pfJetSpecific->mHFEMEnergy = HFEMEnergy;
0320
0321   pfJetSpecific->mChargedHadronMultiplicity = std::round(chargedHadronMultiplicity);
0322   pfJetSpecific->mNeutralHadronMultiplicity = std::round(neutralHadronMultiplicity);
0323   pfJetSpecific->mPhotonMultiplicity = std::round(photonMultiplicity);
0324   pfJetSpecific->mElectronMultiplicity = std::round(electronMultiplicity);
0325   pfJetSpecific->mMuonMultiplicity = std::round(muonMultiplicity);
0326   pfJetSpecific->mHFHadronMultiplicity = std::round(HFHadronMultiplicity);
0327   pfJetSpecific->mHFEMMultiplicity = std::round(HFEMMultiplicity);
0328
0329   pfJetSpecific->mChargedEmEnergy = chargedEmEnergy;
0330   pfJetSpecific->mChargedMuEnergy = chargedMuEnergy;
0331   pfJetSpecific->mNeutralEmEnergy = neutralEmEnergy;
0332   pfJetSpecific->mChargedMultiplicity = std::round(chargedMultiplicity);
0333   pfJetSpecific->mNeutralMultiplicity = std::round(neutralMultiplicity);
0334
0335   pfJetSpecific->mHOEnergy = HOEnergy;
0336
0337   return true;
0338 }
0339
0340 //______________________________________________________________________________
0341 bool reco::makeSpecific(std::vector<reco::CandidatePtr> const& mcparticles, GenJet::Specific* genJetSpecific) {
0342   if (nullptr == genJetSpecific)
0343     return false;
0344
0345   std::vector<reco::CandidatePtr>::const_iterator itMcParticle = mcparticles.begin();
0346   for (; itMcParticle != mcparticles.end(); ++itMcParticle) {
0347     if (itMcParticle->isNull() || !itMcParticle->isAvailable()) {
0348       edm::LogWarning("DataNotFound") << " JetSpecific: MC Particle is invalid\n";
0349       continue;
0350     }
0351
0352     const Candidate* candidate = itMcParticle->get();
0353     if (candidate->hasMasterClone())
0354       candidate = candidate->masterClone().get();
0355     //const GenParticle* genParticle = GenJet::genParticle(candidate);
0356
0357     if (candidate) {
0358       double e = candidate->energy();
0359
0360       // Legacy calo-like definitions
0361       switch (std::abs(candidate->pdgId())) {
0362         case 22:  // photon
0363         case 11:  // e
0364           genJetSpecific->m_EmEnergy += e;
0365           break;
0366         case 211:   // pi
0367         case 321:   // K
0368         case 130:   // KL
0369         case 2212:  // p
0370         case 2112:  // n
0371           genJetSpecific->m_HadEnergy += e;
0372           break;
0373         case 13:  // muon
0374         case 12:  // nu_e
0375         case 14:  // nu_mu
0376         case 16:  // nu_tau
0377           genJetSpecific->m_InvisibleEnergy += e;
0378           break;
0379         default:
0380           genJetSpecific->m_AuxiliaryEnergy += e;
0381       }
0382
0383       // PF-like definitions
0384       switch (std::abs(candidate->pdgId())) {
0385         case 11:  //electron
0386           genJetSpecific->m_ChargedEmEnergy += e;
0387           ++(genJetSpecific->m_ChargedEmMultiplicity);
0388           break;
0389         case 13:  // muon
0390           genJetSpecific->m_MuonEnergy += e;
0391           ++(genJetSpecific->m_MuonMultiplicity);
0392           break;
0393         case 211:   //pi+-
0394         case 321:   //K
0395         case 2212:  //p
0396         case 3222:  //Sigma+
0397         case 3112:  //Sigma-
0398         case 3312:  //Xi-
0399         case 3334:  //Omega-
0400           genJetSpecific->m_ChargedHadronEnergy += e;
0401           ++(genJetSpecific->m_ChargedHadronMultiplicity);
0402           break;
0403         case 310:   //KS0
0404         case 130:   //KL0
0405         case 3122:  //Lambda0
0406         case 3212:  //Sigma0
0407         case 3322:  //Xi0
0408         case 2112:  //n0
0409           genJetSpecific->m_NeutralHadronEnergy += e;
0410           ++(genJetSpecific->m_NeutralHadronMultiplicity);
0411           break;
0412         case 22:  //photon
0413           genJetSpecific->m_NeutralEmEnergy += e;
0414           ++(genJetSpecific->m_NeutralEmMultiplicity);
0415           break;
0416       }
0417     }  // end if found a candidate
0418     else {
0419       edm::LogWarning("DataNotFound") << "reco::makeGenJetSpecific: Referred  GenParticleCandidate "
0420                                       << "is not available in the event\n";
0421     }
0422   }  // end for loop over MC particles
0423
0424   return true;
0425 }
0426
0427 //______________________________________________________________________________
0428 HcalSubdetector reco::hcalSubdetector(int iEta, const HcalTopology& topology) {
0429   int eta = std::abs(iEta);
0430   if (eta <= topology.lastHBRing())
0431     return HcalBarrel;
0432   else if (eta <= topology.lastHERing())
0433     return HcalEndcap;
0434   else if (eta <= topology.lastHFRing())
0435     return HcalForward;
0436   return HcalEmpty;
0437 }