JetProducers/plugins/CATopJetProducer.cc

0001 #include "FWCore/Framework/interface/MakerMacros.h"
0002 #include "RecoJets/JetProducers/plugins/CATopJetProducer.h"
0003
0004 #include <memory>
0005
0006 #include "RecoJets/JetProducers/plugins/FastjetJetProducer.h"
0007
0008 using namespace edm;
0009 using namespace cms;
0010 using namespace reco;
0011 using namespace std;
0012
0013 CATopJetProducer::CATopJetProducer(edm::ParameterSet const& conf)
0014     : FastjetJetProducer(conf),
0015       tagAlgo_(conf.getParameter<int>("tagAlgo")),
0016       ptMin_(conf.getParameter<double>("jetPtMin")),
0017       centralEtaCut_(conf.getParameter<double>("centralEtaCut")),
0018       verbose_(conf.getParameter<bool>("verbose")) {
0019   if (tagAlgo_ == CA_TOPTAGGER) {
0020     legacyCMSTopTagger_ = std::make_unique<CATopJetAlgorithm>(
0021         src_,
0022         conf.getParameter<bool>("verbose"),
0023         conf.getParameter<int>("algorithm"),         // 0 = KT, 1 = CA, 2 = anti-KT
0024         conf.getParameter<int>("useAdjacency"),      // choose adjacency requirement:
0025                                                      //  0 = no adjacency
0026                                                      //  1 = deltar adjacency
0027                                                      //  2 = modified adjacency
0028                                                      //  3 = calotower neirest neigbor based adjacency (untested)
0029         conf.getParameter<double>("centralEtaCut"),  // eta for defining "central" jets
0030         conf.getParameter<double>("jetPtMin"),       // min jet pt
0031         conf.getParameter<std::vector<double>>(
0032             "sumEtBins"),  // sumEt bins over which cuts may vary. vector={bin 0 lower bound, bin 1 lower bound, ...}
0033         conf.getParameter<std::vector<double>>("rBins"),  // Jet distance paramter R. R values depend on sumEt bins.
0034         conf.getParameter<std::vector<double>>("ptFracBins"),  // fraction of hard jet pt that subjet must have (deltap)
0035         conf.getParameter<std::vector<double>>(
0036             "deltarBins"),  // Applicable only if useAdjacency=1. deltar adjacency values for each sumEtBin
0037         conf.getParameter<std::vector<double>>(
0038             "nCellBins"),  // Applicable only if useAdjacency=3. number of cells to consider two subjets adjacent
0039         conf.getParameter<double>("inputEtMin"),  // seed threshold - NOT USED
0040         conf.getParameter<bool>(
0041             "useMaxTower"),  // use max tower as adjacency criterion, otherwise use centroid - NOT USED
0042         conf.getParameter<double>("sumEtEtaCut"),  // eta for event SumEt - NOT USED
0043         conf.getParameter<double>("etFrac")  // fraction of event sumEt / 2 for a jet to be considered "hard" -NOT USED
0044     );
0045   } else if (tagAlgo_ == FJ_CMS_TOPTAG) {
0046     fjCMSTopTagger_ = std::make_unique<fastjet::CMSTopTagger>(conf.getParameter<double>("ptFrac"),
0047                                                               conf.getParameter<double>("rFrac"),
0048                                                               conf.getParameter<double>("adjacencyParam"));
0049   } else if (tagAlgo_ == FJ_JHU_TOPTAG) {
0050     fjJHUTopTagger_ = std::make_unique<fastjet::JHTopTagger>(conf.getParameter<double>("ptFrac"),
0051                                                              conf.getParameter<double>("deltaRCut"),
0052                                                              conf.getParameter<double>("cosThetaWMax"));
0053   } else if (tagAlgo_ == FJ_NSUB_TAG) {
0054     fastjet::JetDefinition::Plugin* plugin = new fastjet::SISConePlugin(0.6, 0.75);
0055     fastjet::JetDefinition NsubJetDef(plugin);
0056     fjNSUBTagger_ = std::make_unique<fastjet::RestFrameNSubjettinessTagger>(NsubJetDef,
0057                                                                             conf.getParameter<double>("tau2Cut"),
0058                                                                             conf.getParameter<double>("cosThetaSCut"),
0059                                                                             conf.getParameter<bool>("useExclusive"));
0060   }
0061 }
0062
0063 void CATopJetProducer::produce(edm::Event& e, const edm::EventSetup& c) { FastjetJetProducer::produce(e, c); }
0064
0065 void CATopJetProducer::runAlgorithm(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0066   if (!doAreaFastjet_ && !doRhoFastjet_) {
0067     fjClusterSeq_ = std::make_shared<fastjet::ClusterSequence>(fjInputs_, *fjJetDefinition_);
0068   } else if (voronoiRfact_ <= 0) {
0069     fjClusterSeq_ =
0070         ClusterSequencePtr(new fastjet::ClusterSequenceArea(fjInputs_, *fjJetDefinition_, *fjAreaDefinition_));
0071   } else {
0072     fjClusterSeq_ = ClusterSequencePtr(
0073         new fastjet::ClusterSequenceVoronoiArea(fjInputs_, *fjJetDefinition_, fastjet::VoronoiAreaSpec(voronoiRfact_)));
0074   }
0075
0076   if (tagAlgo_ == CA_TOPTAGGER) {
0077     (*legacyCMSTopTagger_).run(fjInputs_, fjJets_, fjClusterSeq_);
0078
0079   } else {
0080     //Run the jet clustering
0081     vector<fastjet::PseudoJet> inclusiveJets = fjClusterSeq_->inclusive_jets(ptMin_);
0082
0083     if (verbose_)
0084       cout << "Getting central jets" << endl;
0085     // Find the transient central jets
0086     vector<fastjet::PseudoJet> centralJets;
0087     for (unsigned int i = 0; i < inclusiveJets.size(); i++) {
0088       if (inclusiveJets[i].perp() > ptMin_ && fabs(inclusiveJets[i].rapidity()) < centralEtaCut_) {
0089         centralJets.push_back(inclusiveJets[i]);
0090       }
0091     }
0092
0093     fastjet::CMSTopTagger& CMSTagger = *fjCMSTopTagger_;
0094     fastjet::JHTopTagger& JHUTagger = *fjJHUTopTagger_;
0095     fastjet::RestFrameNSubjettinessTagger& NSUBTagger = *fjNSUBTagger_;
0096
0097     vector<fastjet::PseudoJet>::iterator jetIt = centralJets.begin(), centralJetsEnd = centralJets.end();
0098     if (verbose_)
0099       cout << "Loop over jets" << endl;
0100     for (; jetIt != centralJetsEnd; ++jetIt) {
0101       if (verbose_)
0102         cout << "CMS FJ jet pt: " << (*jetIt).perp() << endl;
0103
0104       fastjet::PseudoJet taggedJet;
0105       if (tagAlgo_ == FJ_CMS_TOPTAG)
0106         taggedJet = CMSTagger.result(*jetIt);
0107       else if (tagAlgo_ == FJ_JHU_TOPTAG)
0108         taggedJet = JHUTagger.result(*jetIt);
0109       else if (tagAlgo_ == FJ_NSUB_TAG)
0110         taggedJet = NSUBTagger.result(*jetIt);
0111       else
0112         cout << "NOT A VALID TAGGING ALGORITHM CHOICE!" << endl;
0113
0114       if (taggedJet != 0)
0115         fjJets_.push_back(taggedJet);
0116     }
0117   }
0118 }
0119
0120 void CATopJetProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0121   edm::ParameterSetDescription desc;
0122   /// Cambridge-Aachen top jet producer parameters
0123   desc.add<int>("tagAlgo", 0);             // choice of top tagging algorithm
0124   desc.add<double>("centralEtaCut", 2.5);  // eta for defining "central" jets
0125   desc.add<bool>("verbose", false);
0126   desc.add<string>("jetCollInstanceName", "caTopSubJets");  // subjet collection
0127   desc.add<int>("algorithm", 1);                            // 0 = KT, 1 = CA, 2 = anti-KT
0128   desc.add<int>("useAdjacency", 2);                         // veto adjacent subjets:
0129                                                             // 0,   no adjacency
0130                                                             //  1,  deltar adjacency
0131                                                             //  2,  modified adjacency
0132                                                             //  3,  calotower neirest neigbor based adjacency (untested)
0133   vector<double> sumEtBinsDefault = {0., 1600., 2600.};
0134   desc.add<vector<double>>(
0135       "sumEtBins",
0136       sumEtBinsDefault);  // sumEt bins over which cuts vary. vector={bin 0 lower bound, bin 1 lower bound, ...}
0137   vector<double> rBinsDefault(3, 0.8);
0138   desc.add<vector<double>>("rBins", rBinsDefault);  // Jet distance paramter R. R values depend on sumEt bins.
0139   vector<double> ptFracBinsDefault(3, 0.05);
0140   desc.add<vector<double>>("ptFracBins", ptFracBinsDefault);  // minimum fraction of central jet pt for subjets (deltap)
0141   vector<double> deltarBinsDefault(3, 0.019);
0142   desc.add<vector<double>>(
0143       "deltarBins", deltarBinsDefault);  // Applicable only if useAdjacency=1. deltar adjacency values for each sumEtBin
0144   vector<double> nCellBinsDefault(3, 1.9);
0145   desc.add<vector<double>>(
0146       "nCellBins",
0147       nCellBinsDefault);  // Applicable only if useAdjacency=3. number of cells apart for two subjets to be considered "independent"
0148   desc.add<bool>("useMaxTower", false);  // use max tower in adjacency criterion, otherwise use centroid - NOT USED
0149   desc.add<double>("sumEtEtaCut", 3.0);  // eta for event SumEt - NOT USED
0150   desc.add<double>("etFrac", 0.7);       // fraction of event sumEt / 2 for a jet to be considered "hard" - NOT USED
0151   desc.add<double>("ptFrac", 0.05);
0152   desc.add<double>("rFrac", 0.);
0153   desc.add<double>("adjacencyParam", 0.);
0154   desc.add<double>("deltaRCut", 0.19);
0155   desc.add<double>("cosThetaWMax", 0.7);
0156   desc.add<double>("tau2Cut", 0.);
0157   desc.add<double>("cosThetaSCut", 0.);
0158   desc.add<bool>("useExclusive", false);
0159   ////// From FastjetJetProducer
0160   FastjetJetProducer::fillDescriptionsFromFastJetProducer(desc);
0161   ///// From VirtualJetProducer
0162   VirtualJetProducer::fillDescriptionsFromVirtualJetProducer(desc);
0163   /////////////////////
0164   descriptions.add("CATopJetProducer", desc);
0165 }
0166 //define this as a plug-in
0167 DEFINE_FWK_MODULE(CATopJetProducer);