Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-09-07 04:37:34

0001 #include <memory>
0002 
0003 #include "RecoJets/JetProducers/interface/NjettinessAdder.h"
0004 
0005 #include "FWCore/Framework/interface/MakerMacros.h"
0006 
0007 NjettinessAdder::NjettinessAdder(const edm::ParameterSet& iConfig)
0008     : src_(iConfig.getParameter<edm::InputTag>("src")),
0009       src_token_(consumes<edm::View<reco::Jet>>(src_)),
0010       Njets_(iConfig.getParameter<std::vector<unsigned>>("Njets")),
0011       measureDefinition_(iConfig.getParameter<unsigned>("measureDefinition")),
0012       beta_(iConfig.getParameter<double>("beta")),
0013       R0_(iConfig.getParameter<double>("R0")),
0014       Rcutoff_(iConfig.getParameter<double>("Rcutoff")),
0015       axesDefinition_(iConfig.getParameter<unsigned>("axesDefinition")),
0016       nPass_(iConfig.getParameter<int>("nPass")),
0017       akAxesR0_(iConfig.getParameter<double>("akAxesR0")) {
0018   edm::InputTag srcWeights = iConfig.getParameter<edm::InputTag>("srcWeights");
0019   if (!srcWeights.label().empty())
0020     input_weights_token_ = consumes<edm::ValueMap<float>>(srcWeights);
0021 
0022   for (std::vector<unsigned>::const_iterator n = Njets_.begin(); n != Njets_.end(); ++n) {
0023     std::ostringstream tauN_str;
0024     tauN_str << "tau" << *n;
0025 
0026     produces<edm::ValueMap<float>>(tauN_str.str());
0027   }
0028 
0029   // Get the measure definition
0030   fastjet::contrib::NormalizedMeasure normalizedMeasure(beta_, R0_);
0031   fastjet::contrib::UnnormalizedMeasure unnormalizedMeasure(beta_);
0032   fastjet::contrib::OriginalGeometricMeasure geometricMeasure(beta_);  // changed in 1.020
0033   fastjet::contrib::NormalizedCutoffMeasure normalizedCutoffMeasure(beta_, R0_, Rcutoff_);
0034   fastjet::contrib::UnnormalizedCutoffMeasure unnormalizedCutoffMeasure(beta_, Rcutoff_);
0035   //fastjet::contrib::GeometricCutoffMeasure     geometricCutoffMeasure   (beta_,Rcutoff_); // removed in 1.020
0036 
0037   fastjet::contrib::MeasureDefinition const* measureDef = nullptr;
0038   switch (measureDefinition_) {
0039     case UnnormalizedMeasure:
0040       measureDef = &unnormalizedMeasure;
0041       break;
0042     case OriginalGeometricMeasure:
0043       measureDef = &geometricMeasure;
0044       break;  // changed in 1.020
0045     case NormalizedCutoffMeasure:
0046       measureDef = &normalizedCutoffMeasure;
0047       break;
0048     case UnnormalizedCutoffMeasure:
0049       measureDef = &unnormalizedCutoffMeasure;
0050       break;
0051     //case GeometricCutoffMeasure : measureDef = &geometricCutoffMeasure; break; // removed in 1.020
0052     case NormalizedMeasure:
0053     default:
0054       measureDef = &normalizedMeasure;
0055       break;
0056   }
0057 
0058   // Get the axes definition
0059   fastjet::contrib::KT_Axes kt_axes;
0060   fastjet::contrib::CA_Axes ca_axes;
0061   fastjet::contrib::AntiKT_Axes antikt_axes(akAxesR0_);
0062   fastjet::contrib::WTA_KT_Axes wta_kt_axes;
0063   fastjet::contrib::WTA_CA_Axes wta_ca_axes;
0064   fastjet::contrib::OnePass_KT_Axes onepass_kt_axes;
0065   fastjet::contrib::OnePass_CA_Axes onepass_ca_axes;
0066   fastjet::contrib::OnePass_AntiKT_Axes onepass_antikt_axes(akAxesR0_);
0067   fastjet::contrib::OnePass_WTA_KT_Axes onepass_wta_kt_axes;
0068   fastjet::contrib::OnePass_WTA_CA_Axes onepass_wta_ca_axes;
0069   fastjet::contrib::MultiPass_Axes multipass_axes(nPass_);
0070 
0071   fastjet::contrib::AxesDefinition const* axesDef = nullptr;
0072   switch (axesDefinition_) {
0073     case KT_Axes:
0074     default:
0075       axesDef = &kt_axes;
0076       break;
0077     case CA_Axes:
0078       axesDef = &ca_axes;
0079       break;
0080     case AntiKT_Axes:
0081       axesDef = &antikt_axes;
0082       break;
0083     case WTA_KT_Axes:
0084       axesDef = &wta_kt_axes;
0085       break;
0086     case WTA_CA_Axes:
0087       axesDef = &wta_ca_axes;
0088       break;
0089     case OnePass_KT_Axes:
0090       axesDef = &onepass_kt_axes;
0091       break;
0092     case OnePass_CA_Axes:
0093       axesDef = &onepass_ca_axes;
0094       break;
0095     case OnePass_AntiKT_Axes:
0096       axesDef = &onepass_antikt_axes;
0097       break;
0098     case OnePass_WTA_KT_Axes:
0099       axesDef = &onepass_wta_kt_axes;
0100       break;
0101     case OnePass_WTA_CA_Axes:
0102       axesDef = &onepass_wta_ca_axes;
0103       break;
0104     case MultiPass_Axes:
0105       axesDef = &multipass_axes;
0106       break;
0107   };
0108 
0109   routine_ = std::make_unique<fastjet::contrib::Njettiness>(*axesDef, *measureDef);
0110 }
0111 
0112 void NjettinessAdder::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0113   // read input collection
0114   edm::Handle<edm::View<reco::Jet>> jets;
0115   iEvent.getByToken(src_token_, jets);
0116 
0117   // Get Weights Collection
0118   if (!input_weights_token_.isUninitialized())
0119     weightsHandle_ = &iEvent.get(input_weights_token_);
0120 
0121   for (std::vector<unsigned>::const_iterator n = Njets_.begin(); n != Njets_.end(); ++n) {
0122     std::ostringstream tauN_str;
0123     tauN_str << "tau" << *n;
0124 
0125     // prepare room for output
0126     std::vector<float> tauN;
0127     tauN.reserve(jets->size());
0128 
0129     for (typename edm::View<reco::Jet>::const_iterator jetIt = jets->begin(); jetIt != jets->end(); ++jetIt) {
0130       edm::Ptr<reco::Jet> jetPtr = jets->ptrAt(jetIt - jets->begin());
0131 
0132       float t = getTau(*n, jetPtr);
0133 
0134       tauN.push_back(t);
0135     }
0136 
0137     auto outT = std::make_unique<edm::ValueMap<float>>();
0138     edm::ValueMap<float>::Filler fillerT(*outT);
0139     fillerT.insert(jets, tauN.begin(), tauN.end());
0140     fillerT.fill();
0141 
0142     iEvent.put(std::move(outT), tauN_str.str());
0143   }
0144 }
0145 
0146 float NjettinessAdder::getTau(unsigned num, const edm::Ptr<reco::Jet>& object) const {
0147   std::vector<fastjet::PseudoJet> fjParticles;
0148   for (unsigned k = 0; k < object->numberOfDaughters(); ++k) {
0149     const reco::CandidatePtr& dp = object->daughterPtr(k);
0150     if (dp.isNonnull() && dp.isAvailable()) {
0151       // Here, the daughters are the "end" node, so this is a PFJet
0152       if (dp->numberOfDaughters() == 0) {
0153         if (object->isWeighted()) {
0154           if (input_weights_token_.isUninitialized())
0155             throw cms::Exception("MissingConstituentWeight")
0156                 << "ECFAdder: No weights (e.g. PUPPI) given for weighted jet collection" << std::endl;
0157           float w = (*weightsHandle_)[dp];
0158           fjParticles.push_back(fastjet::PseudoJet(dp->px() * w, dp->py() * w, dp->pz() * w, dp->energy() * w));
0159         } else
0160           fjParticles.push_back(fastjet::PseudoJet(dp->px(), dp->py(), dp->pz(), dp->energy()));
0161       } else {  // Otherwise, this is a BasicJet, so you need to descend further.
0162         auto subjet = dynamic_cast<reco::Jet const*>(dp.get());
0163         for (unsigned l = 0; l < subjet->numberOfDaughters(); ++l) {
0164           if (subjet != nullptr) {
0165             const reco::CandidatePtr& ddp = subjet->daughterPtr(l);
0166             if (subjet->isWeighted()) {
0167               if (input_weights_token_.isUninitialized())
0168                 throw cms::Exception("MissingConstituentWeight")
0169                     << "ECFAdder: No weights (e.g. PUPPI) given for weighted jet collection" << std::endl;
0170               float w = (*weightsHandle_)[ddp];
0171               fjParticles.push_back(fastjet::PseudoJet(ddp->px() * w, ddp->py() * w, ddp->pz() * w, ddp->energy() * w));
0172             } else
0173               fjParticles.push_back(fastjet::PseudoJet(ddp->px(), ddp->py(), ddp->pz(), ddp->energy()));
0174           } else {
0175             edm::LogWarning("MissingJetConstituent")
0176                 << "BasicJet constituent required for N-jettiness computation is missing!";
0177           }
0178         }
0179       }  // end if basic jet
0180     }  // end if daughter pointer is nonnull and available
0181     else
0182       edm::LogWarning("MissingJetConstituent") << "Jet constituent required for N-jettiness computation is missing!";
0183   }
0184 
0185   return routine_->getTau(num, fjParticles);
0186 }
0187 
0188 DEFINE_FWK_MODULE(NjettinessAdder);