JetProducers/plugins/HTTTopJetProducer.cc

0001 #include "FWCore/Framework/interface/MakerMacros.h"
0002 #include "DataFormats/BTauReco/interface/HTTTopJetTagInfo.h"
0003 #include "DataFormats/JetReco/interface/BasicJetCollection.h"
0004 #include "HTTTopJetProducer.h"
0005
0006 #include <memory>
0007
0008 using namespace edm;
0009 using namespace cms;
0010 using namespace reco;
0011 using namespace std;
0012
0013 HTTTopJetProducer::HTTTopJetProducer(edm::ParameterSet const& conf)
0014     : FastjetJetProducer(conf),
0015       optimalR_(false),
0016       qJets_(false),
0017       minFatjetPt_(200.),
0018       minSubjetPt_(20.),
0019       minCandPt_(200.),
0020       maxFatjetAbsEta_(2.5),
0021       subjetMass_(30.),
0022       muCut_(0.8),
0023       filtR_(0.3),
0024       filtN_(5),
0025       mode_(0),
0026       minCandMass_(150.),
0027       maxCandMass_(200.),
0028       massRatioWidth_(15),
0029       minM23Cut_(0.35),
0030       minM13Cut_(0.2),
0031       maxM13Cut_(1.3),
0032       maxR_(1.5),
0033       minR_(0.5),
0034       rejectMinR_(false),
0035       verbose_(false) {
0036   // Read in all the options from the configuration
0037   optimalR_ = conf.getParameter<bool>("optimalR");
0038   qJets_ = conf.getParameter<bool>("qJets");
0039   minFatjetPt_ = conf.getParameter<double>("minFatjetPt");
0040   minSubjetPt_ = conf.getParameter<double>("minSubjetPt");
0041   minCandPt_ = conf.getParameter<double>("minCandPt");
0042   maxFatjetAbsEta_ = conf.getParameter<double>("maxFatjetAbsEta");
0043   subjetMass_ = conf.getParameter<double>("subjetMass");
0044   muCut_ = conf.getParameter<double>("muCut");
0045   filtR_ = conf.getParameter<double>("filtR");
0046   filtN_ = conf.getParameter<int>("filtN");
0047   mode_ = conf.getParameter<int>("mode");
0048   minCandMass_ = conf.getParameter<double>("minCandMass");
0049   maxCandMass_ = conf.getParameter<double>("maxCandMass");
0050   massRatioWidth_ = conf.getParameter<double>("massRatioWidth");
0051   minM23Cut_ = conf.getParameter<double>("minM23Cut");
0052   minM13Cut_ = conf.getParameter<double>("minM13Cut");
0053   maxM13Cut_ = conf.getParameter<double>("maxM13Cut");
0054   maxR_ = conf.getParameter<double>("maxR");
0055   minR_ = conf.getParameter<double>("minR");
0056   rejectMinR_ = conf.getParameter<bool>("rejectMinR");
0057   verbose_ = conf.getParameter<bool>("verbose");
0058
0059   // Create the tagger-wrapper
0060   produces<HTTTopJetTagInfoCollection>();
0061
0062   // Signal to the VirtualJetProducer that we have to add HTT information
0063   fromHTTTopJetProducer_ = true;
0064
0065   fjHEPTopTagger_ = std::make_unique<fastjet::HEPTopTaggerV2>(optimalR_,
0066                                                               qJets_,
0067                                                               minSubjetPt_,
0068                                                               minCandPt_,
0069                                                               subjetMass_,
0070                                                               muCut_,
0071                                                               filtR_,
0072                                                               filtN_,
0073                                                               mode_,
0074                                                               minCandMass_,
0075                                                               maxCandMass_,
0076                                                               massRatioWidth_,
0077                                                               minM23Cut_,
0078                                                               minM13Cut_,
0079                                                               maxM13Cut_,
0080                                                               rejectMinR_);
0081 }
0082
0083 void HTTTopJetProducer::produce(edm::Event& e, const edm::EventSetup& c) {
0084   if (qJets_) {
0085     edm::Service<edm::RandomNumberGenerator> rng;
0086     CLHEP::HepRandomEngine* engine = &rng->getEngine(e.streamID());
0087     fjHEPTopTagger_->set_rng(engine);
0088   }
0089
0090   FastjetJetProducer::produce(e, c);
0091 }
0092
0093 void HTTTopJetProducer::runAlgorithm(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0094   if (!doAreaFastjet_ && !doRhoFastjet_) {
0095     fjClusterSeq_ = std::make_shared<fastjet::ClusterSequence>(fjInputs_, *fjJetDefinition_);
0096   } else if (voronoiRfact_ <= 0) {
0097     fjClusterSeq_ =
0098         ClusterSequencePtr(new fastjet::ClusterSequenceArea(fjInputs_, *fjJetDefinition_, *fjAreaDefinition_));
0099   } else {
0100     fjClusterSeq_ = ClusterSequencePtr(
0101         new fastjet::ClusterSequenceVoronoiArea(fjInputs_, *fjJetDefinition_, fastjet::VoronoiAreaSpec(voronoiRfact_)));
0102   }
0103
0104   //Run the jet clustering
0105   vector<fastjet::PseudoJet> inclusiveJets = fjClusterSeq_->inclusive_jets(minFatjetPt_);
0106
0107   if (verbose_)
0108     cout << "Getting central jets" << endl;
0109   // Find the transient central jets
0110   vector<fastjet::PseudoJet> centralJets;
0111   for (unsigned int i = 0; i < inclusiveJets.size(); i++) {
0112     if (inclusiveJets[i].perp() > minFatjetPt_ && fabs(inclusiveJets[i].rapidity()) < maxFatjetAbsEta_) {
0113       centralJets.push_back(inclusiveJets[i]);
0114     }
0115   }
0116
0117   fastjet::HEPTopTaggerV2& HEPTagger = *fjHEPTopTagger_;
0118
0119   vector<fastjet::PseudoJet>::iterator jetIt = centralJets.begin(), centralJetsEnd = centralJets.end();
0120   if (verbose_)
0121     cout << "Loop over jets" << endl;
0122   for (; jetIt != centralJetsEnd; ++jetIt) {
0123     if (verbose_)
0124       cout << "CMS FJ jet pt: " << (*jetIt).perp() << endl;
0125
0126     fastjet::PseudoJet taggedJet;
0127
0128     taggedJet = HEPTagger.result(*jetIt);
0129
0130     if (taggedJet != 0) {
0131       fjJets_.push_back(taggedJet);
0132     }
0133   }
0134 }
0135
0136 void HTTTopJetProducer::addHTTTopJetTagInfoCollection(edm::Event& iEvent,
0137                                                       const edm::EventSetup& iSetup,
0138                                                       edm::OrphanHandle<reco::BasicJetCollection>& oh) {
0139   // Set up output list
0140   auto tagInfos = std::make_unique<HTTTopJetTagInfoCollection>();
0141
0142   // Loop over jets
0143   for (size_t ij = 0; ij != fjJets_.size(); ij++) {
0144     HTTTopJetProperties properties;
0145     HTTTopJetTagInfo tagInfo;
0146
0147     // Black magic:
0148     // In the standard CA treatment the RefToBase is made from the handle directly
0149     // Since we only have a OrphanHandle (the JetCollection is created by this process)
0150     // we have to take the detour via the Ref
0151     edm::Ref<reco::BasicJetCollection> ref(oh, ij);
0152     edm::RefToBase<reco::Jet> rtb(ref);
0153
0154     fastjet::HEPTopTaggerV2Structure* s = (fastjet::HEPTopTaggerV2Structure*)fjJets_[ij].structure_non_const_ptr();
0155
0156     properties.fjMass = s->fj_mass();
0157     properties.fjPt = s->fj_pt();
0158     properties.fjEta = s->fj_eta();
0159     properties.fjPhi = s->fj_phi();
0160
0161     properties.topMass = s->top_mass();
0162     properties.unfilteredMass = s->unfiltered_mass();
0163     properties.prunedMass = s->pruned_mass();
0164     properties.fRec = s->fRec();
0165     properties.massRatioPassed = s->mass_ratio_passed();
0166
0167     properties.ropt = s->ropt();
0168     properties.roptCalc = s->roptCalc();
0169     properties.ptForRoptCalc = s->ptForRoptCalc();
0170
0171     properties.tau1Unfiltered = s->Tau1Unfiltered();
0172     properties.tau2Unfiltered = s->Tau2Unfiltered();
0173     properties.tau3Unfiltered = s->Tau3Unfiltered();
0174     properties.tau1Filtered = s->Tau1Filtered();
0175     properties.tau2Filtered = s->Tau2Filtered();
0176     properties.tau3Filtered = s->Tau3Filtered();
0177     properties.qWeight = s->qweight();
0178     properties.qEpsilon = s->qepsilon();
0179     properties.qSigmaM = s->qsigmaM();
0180
0181     tagInfo.insert(rtb, properties);
0182     tagInfos->push_back(tagInfo);
0183   }
0184
0185   iEvent.put(std::move(tagInfos));
0186 };
0187
0188 void HTTTopJetProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0189   edm::ParameterSetDescription desc;
0190
0191   desc.add<bool>("optimalR", true);
0192   desc.add<bool>("qJets", false);
0193   desc.add<double>("minFatjetPt", 200.);
0194   desc.add<double>("minSubjetPt", 0.);
0195   desc.add<double>("minCandPt", 0.);
0196   desc.add<double>("maxFatjetAbsEta", 99.);
0197   desc.add<double>("subjetMass", 30.);
0198   desc.add<double>("filtR", 0.3);
0199   desc.add<int>("filtN", 5);
0200   desc.add<int>("mode", 4);
0201   desc.add<double>("minCandMass", 0.);
0202   desc.add<double>("maxCandMass", 999999.);
0203   desc.add<double>("massRatioWidth", 999999.);
0204   desc.add<double>("minM23Cut", 0.);
0205   desc.add<double>("minM13Cut", 0.);
0206   desc.add<double>("maxM13Cut", 999999.);
0207   desc.add<double>("maxR", 1.5);
0208   desc.add<double>("minR", 0.5);
0209   desc.add<bool>("rejectMinR", false);
0210   desc.add<bool>("verbose", false);
0211
0212   desc.add<int>("algorithm", 1);  // where is it needed?
0213
0214   ////// From FastjetJetProducer
0215   FastjetJetProducer::fillDescriptionsFromFastJetProducer(desc);
0216   ///// From VirtualJetProducer
0217   desc.add<std::string>("jetCollInstanceName", "SubJets");
0218   VirtualJetProducer::fillDescriptionsFromVirtualJetProducer(desc);
0219
0220   descriptions.add("HTTTopJetProducer", desc);
0221 }
0222
0223 //define this as a plug-in
0224 DEFINE_FWK_MODULE(HTTTopJetProducer);