File indexing completed on 2023-03-17 11:18:31
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
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
0060 produces<HTTTopJetTagInfoCollection>();
0061
0062
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
0105 vector<fastjet::PseudoJet> inclusiveJets = fjClusterSeq_->inclusive_jets(minFatjetPt_);
0106
0107 if (verbose_)
0108 cout << "Getting central jets" << endl;
0109
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
0140 auto tagInfos = std::make_unique<HTTTopJetTagInfoCollection>();
0141
0142
0143 for (size_t ij = 0; ij != fjJets_.size(); ij++) {
0144 HTTTopJetProperties properties;
0145 HTTTopJetTagInfo tagInfo;
0146
0147
0148
0149
0150
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);
0213
0214
0215 FastjetJetProducer::fillDescriptionsFromFastJetProducer(desc);
0216
0217 desc.add<std::string>("jetCollInstanceName", "SubJets");
0218 VirtualJetProducer::fillDescriptionsFromVirtualJetProducer(desc);
0219
0220 descriptions.add("HTTTopJetProducer", desc);
0221 }
0222
0223
0224 DEFINE_FWK_MODULE(HTTTopJetProducer);