JetProducers/plugins/CompoundJetProducer.cc

0001
0002 #include "FWCore/Framework/interface/MakerMacros.h"
0003 #include "RecoJets/JetProducers/plugins/CompoundJetProducer.h"
0004 #include "FWCore/Utilities/interface/Exception.h"
0005 #include "RecoJets/JetProducers/interface/JetSpecific.h"
0006 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0007 #include "Geometry/Records/interface/HcalRecNumberingRecord.h"
0008
0009 using namespace std;
0010 using namespace reco;
0011 using namespace edm;
0012 using namespace cms;
0013
0014 CompoundJetProducer::CompoundJetProducer(edm::ParameterSet const& conf) : VirtualJetProducer(conf) {
0015   produces<reco::BasicJetCollection>();
0016   // the subjet collections are set through the config file in the "jetCollInstanceName" field.
0017 }
0018
0019 void CompoundJetProducer::inputTowers() {
0020   fjCompoundJets_.clear();
0021   VirtualJetProducer::inputTowers();
0022 }
0023
0024 void CompoundJetProducer::output(edm::Event& iEvent, edm::EventSetup const& iSetup) {
0025   // Write jets and constitutents. Will use fjJets_.
0026   switch (jetTypeE) {
0027     case JetType::CaloJet:
0028       writeCompoundJets<reco::CaloJet>(iEvent, iSetup);
0029       break;
0030     case JetType::PFJet:
0031       writeCompoundJets<reco::PFJet>(iEvent, iSetup);
0032       break;
0033     case JetType::GenJet:
0034       writeCompoundJets<reco::GenJet>(iEvent, iSetup);
0035       break;
0036     case JetType::BasicJet:
0037       writeCompoundJets<reco::BasicJet>(iEvent, iSetup);
0038       break;
0039     default:
0040       throw cms::Exception("InvalidInput") << "invalid jet type in CompoundJetProducer\n";
0041       break;
0042   };
0043 }
0044
0045 /// function template to write out the outputs
0046 template <class T>
0047 void CompoundJetProducer::writeCompoundJets(edm::Event& iEvent, edm::EventSetup const& iSetup) {
0048   // get a list of output jets
0049   auto jetCollection = std::make_unique<reco::BasicJetCollection>();
0050   // get a list of output subjets
0051   auto subjetCollection = std::make_unique<std::vector<T>>();
0052
0053   // This will store the handle for the subjets after we write them
0054   edm::OrphanHandle<std::vector<T>> subjetHandleAfterPut;
0055   // this is the mapping of subjet to hard jet
0056   std::vector<std::vector<int>> indices;
0057   // this is the list of hardjet 4-momenta
0058   std::vector<math::XYZTLorentzVector> p4_hardJets;
0059   // this is the hardjet areas
0060   std::vector<double> area_hardJets;
0061
0062   [[maybe_unused]] const CaloGeometry* pGeometry = nullptr;
0063   [[maybe_unused]] const HcalTopology* pTopology = nullptr;
0064   if constexpr (std::is_same_v<T, reco::CaloJet>) {
0065     pGeometry = &getGeometry(iSetup);
0066     pTopology = &getTopology(iSetup);
0067   }
0068
0069   // Loop over the hard jets
0070   std::vector<CompoundPseudoJet>::const_iterator it = fjCompoundJets_.begin(), iEnd = fjCompoundJets_.end(),
0071                                                  iBegin = fjCompoundJets_.begin();
0072   indices.resize(fjCompoundJets_.size());
0073   for (; it != iEnd; ++it) {
0074     int jetIndex = it - iBegin;
0075     fastjet::PseudoJet localJet = it->hardJet();
0076     // Get the 4-vector for the hard jet
0077     p4_hardJets.push_back(math::XYZTLorentzVector(localJet.px(), localJet.py(), localJet.pz(), localJet.e()));
0078     area_hardJets.push_back(it->hardJetArea());
0079
0080     // create the subjet list
0081     std::vector<CompoundPseudoSubJet>::const_iterator itSubJetBegin = it->subjets().begin(), itSubJet = itSubJetBegin,
0082                                                       itSubJetEnd = it->subjets().end();
0083     for (; itSubJet != itSubJetEnd; ++itSubJet) {
0084       fastjet::PseudoJet subjet = itSubJet->subjet();
0085       math::XYZTLorentzVector p4Subjet(subjet.px(), subjet.py(), subjet.pz(), subjet.e());
0086       reco::Particle::Point point(0, 0, 0);
0087
0088       // This will hold ptr's to the subjets
0089       std::vector<reco::CandidatePtr> subjetConstituents;
0090
0091       // Get the transient subjet constituents from fastjet
0092       std::vector<int> const& subjetFastjetConstituentIndices = itSubJet->constituents();
0093       std::vector<int>::const_iterator fastSubIt = subjetFastjetConstituentIndices.begin(),
0094                                        transConstEnd = subjetFastjetConstituentIndices.end();
0095       for (; fastSubIt != transConstEnd; ++fastSubIt) {
0096         // Add a ptr to this constituent
0097         if (*fastSubIt < static_cast<int>(inputs_.size()))
0098           subjetConstituents.push_back(inputs_[*fastSubIt]);
0099       }
0100
0101       // This holds the subjet-to-hardjet mapping
0102       indices[jetIndex].push_back(subjetCollection->size());
0103
0104       // Add the concrete subjet type to the subjet list to write to event record
0105       T jet;
0106       if constexpr (std::is_same_v<T, reco::CaloJet>) {
0107         reco::writeSpecific(jet, p4Subjet, point, subjetConstituents, *pGeometry, *pTopology);
0108       } else {
0109         reco::writeSpecific(jet, p4Subjet, point, subjetConstituents);
0110       }
0111       jet.setJetArea(itSubJet->subjetArea());
0112       subjetCollection->push_back(jet);
0113     }
0114   }
0115   // put subjets into event record
0116   subjetHandleAfterPut = iEvent.put(std::move(subjetCollection), jetCollInstanceName_);
0117
0118   // Now create the hard jets with ptr's to the subjets as constituents
0119   std::vector<math::XYZTLorentzVector>::const_iterator ip4 = p4_hardJets.begin(), ip4Begin = p4_hardJets.begin(),
0120                                                        ip4End = p4_hardJets.end();
0121
0122   for (; ip4 != ip4End; ++ip4) {
0123     int p4_index = ip4 - ip4Begin;
0124     std::vector<int>& ind = indices[p4_index];
0125     std::vector<reco::CandidatePtr> i_hardJetConstituents;
0126     // Add the subjets to the hard jet
0127     for (std::vector<int>::const_iterator isub = ind.begin(); isub != ind.end(); ++isub) {
0128       reco::CandidatePtr candPtr(subjetHandleAfterPut, *isub, false);
0129       i_hardJetConstituents.push_back(candPtr);
0130     }
0131     reco::Particle::Point point(0, 0, 0);
0132     reco::BasicJet toput(*ip4, point, i_hardJetConstituents);
0133     toput.setJetArea(area_hardJets[ip4 - ip4Begin]);
0134     jetCollection->push_back(toput);
0135   }
0136
0137   // put hard jets into event record
0138   iEvent.put(std::move(jetCollection));
0139 }