PartonShowerVeto/src/JetMatchingMGFastJet.cc

0001 #include "GeneratorInterface/PartonShowerVeto/interface/JetMatchingMGFastJet.h"
0002 #include "FWCore/Utilities/interface/Exception.h"
0003 #include <iomanip>
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0006 #include "SimDataFormats/GeneratorProducts/interface/LHECommonBlocks.h"
0007 #include "GeneratorInterface/LHEInterface/interface/LHERunInfo.h"
0008 #include "GeneratorInterface/LHEInterface/interface/LHEEvent.h"
0009
0010 #include <boost/algorithm/string/trim.hpp>
0011
0012 namespace gen {
0013
0014   extern "C" {
0015
0016 #define PARAMLEN 20
0017   namespace {
0018     struct Param {
0019       Param(const std::string &str) {
0020         int len = std::min(PARAMLEN, (int)str.length());
0021         std::memcpy(value, str.c_str(), len);
0022         std::memset(value + len, ' ', PARAMLEN - len);
0023       }
0024
0025       char value[PARAMLEN];
0026     };
0027   }  // namespace
0028
0029   extern struct UPPRIV {
0030     int lnhin, lnhout;
0031     int mscal, ievnt;
0032     int ickkw, iscale;
0033   } uppriv_;
0034
0035   extern struct MEMAIN {
0036     double etcjet, rclmax, etaclmax, qcut, showerkt, clfact;
0037     int maxjets, minjets, iexcfile, ktsche;
0038     int mektsc, nexcres, excres[30];
0039     int nqmatch, excproc, iexcproc[1000], iexcval[1000];
0040     bool nosingrad, jetprocs;
0041   } memain_;
0042
0043   extern struct OUTTREE { int flag; } outtree_;
0044   }
0045
0046   template <typename T>
0047   T JetMatchingMGFastJet::parseParameter(const std::string &value) {
0048     std::istringstream ss(value);
0049     T result;
0050     ss >> result;
0051     return result;
0052   }
0053
0054   template <>
0055   std::string JetMatchingMGFastJet::parseParameter(const std::string &value) {
0056     std::string result;
0057     if (!result.empty() && result[0] == '\'')
0058       result = result.substr(1);
0059     if (!result.empty() && result[result.length() - 1] == '\'')
0060       result.resize(result.length() - 1);
0061     return result;
0062   }
0063
0064   template <>
0065   bool JetMatchingMGFastJet::parseParameter(const std::string &value_) {
0066     std::string value(value_);
0067     std::transform(value.begin(), value.end(), value.begin(), (int (*)(int))std::toupper);
0068     return value == "T" || value == "Y" || value == "True" || value == "1" || value == ".TRUE.";
0069   }
0070
0071   template <typename T>
0072   T JetMatchingMGFastJet::getParameter(const std::map<std::string, std::string> &params,
0073                                        const std::string &var,
0074                                        const T &defValue) {
0075     std::map<std::string, std::string>::const_iterator pos = params.find(var);
0076     if (pos == params.end())
0077       return defValue;
0078     return parseParameter<T>(pos->second);
0079   }
0080
0081   template <typename T>
0082   T JetMatchingMGFastJet::getParameter(const std::string &var, const T &defValue) const {
0083     return getParameter(mgParams, var, defValue);
0084   }
0085
0086   static std::map<std::string, std::string> parseHeader(const std::vector<std::string> &header) {
0087     std::map<std::string, std::string> params;
0088
0089     for (std::vector<std::string>::const_iterator iter = header.begin(); iter != header.end(); ++iter) {
0090       std::string line = *iter;
0091       if (line.empty() || line[0] == '#')
0092         continue;
0093
0094       std::string::size_type pos = line.find('!');
0095       if (pos != std::string::npos)
0096         line.resize(pos);
0097
0098       pos = line.find('=');
0099       if (pos == std::string::npos)
0100         continue;
0101
0102       std::string var = boost::algorithm::trim_copy(line.substr(pos + 1));
0103       std::string value = boost::algorithm::trim_copy(line.substr(0, pos));
0104
0105       params[var] = value;
0106     }
0107
0108     return params;
0109   }
0110
0111   template <typename T>
0112   void JetMatchingMGFastJet::updateOrDie(const std::map<std::string, std::string> &params,
0113                                          T &param,
0114                                          const std::string &name) {
0115     if (param < 0) {
0116       param = getParameter(params, name, param);
0117     }
0118     if (param < 0)
0119       throw cms::Exception("Generator|PartonShowerVeto") << "The MGParamCMS header does not specify the jet "
0120                                                             "matching parameter \""
0121                                                          << name
0122                                                          << "\", but it "
0123                                                             "is requested by the CMSSW configuration."
0124                                                          << std::endl;
0125   }
0126
0127   JetMatchingMGFastJet::JetMatchingMGFastJet(const edm::ParameterSet &params)
0128       : JetMatching(params), runInitialized(false), fJetFinder(nullptr), fIsInit(false) {
0129     std::string mode = params.getParameter<std::string>("mode");
0130     if (mode == "inclusive") {
0131       soup = false;
0132       exclusive = false;
0133     } else if (mode == "exclusive") {
0134       soup = false;
0135       exclusive = true;
0136     } else if (mode == "auto")
0137       soup = true;
0138     else
0139       throw cms::Exception("Generator|LHEInterface") << "MGFastJet jet matching scheme requires \"mode\" "
0140                                                         "parameter to be set to either \"inclusive\", "
0141                                                         "\"exclusive\" or \"auto\"."
0142                                                      << std::endl;
0143
0144     memain_.etcjet = 0.;
0145     memain_.rclmax = 0.0;
0146     memain_.clfact = 0.0;
0147     memain_.ktsche = 0.0;
0148     memain_.etaclmax = params.getParameter<double>("MEMAIN_etaclmax");
0149     memain_.qcut = params.getParameter<double>("MEMAIN_qcut");
0150     memain_.minjets = params.getParameter<int>("MEMAIN_minjets");
0151     memain_.maxjets = params.getParameter<int>("MEMAIN_maxjets");
0152     memain_.showerkt = params.getParameter<double>("MEMAIN_showerkt");
0153     memain_.nqmatch = params.getParameter<int>("MEMAIN_nqmatch");
0154     outtree_.flag = params.getParameter<int>("outTree_flag");
0155     std::string list_excres = params.getParameter<std::string>("MEMAIN_excres");
0156     std::vector<std::string> elems;
0157     std::stringstream ss(list_excres);
0158     std::string item;
0159     int index = 0;
0160     while (std::getline(ss, item, ',')) {
0161       elems.push_back(item);
0162       memain_.excres[index] = std::atoi(item.c_str());
0163       index++;
0164     }
0165     memain_.nexcres = index;
0166
0167     fDJROutFlag = params.getParameter<int>("outTree_flag");
0168   }
0169
0170   double JetMatchingMGFastJet::getJetEtaMax() const { return memain_.etaclmax; }
0171
0172   bool JetMatchingMGFastJet::initAfterBeams() {
0173     if (fIsInit)
0174       return true;
0175
0176     //
0177     doMerge = uppriv_.ickkw;
0178     // doMerge = true;
0179     qCut = memain_.qcut;  //
0180     nQmatch = memain_.nqmatch;
0181     clFact = 1.;  // Default value
0182                   // NOTE: ME2pythia seems to default to 1.5 - need to check !!!
0183                   // in general, needs to read key ALPSFACT from LHE file - fix CMSSW code !!!
0184     nJetMin = memain_.minjets;
0185     nJetMax = memain_.maxjets;
0186
0187     etaJetMax = memain_.etaclmax;
0188
0189     coneRadius = 1.0;
0190     jetAlgoPower = 1;  //   this is the kT algorithm !!!
0191
0192     // Matching procedure
0193     //
0194     qCutSq = pow(qCut, 2);
0195     // this should be something like memaev_.iexc
0196     fExcLocal = true;
0197
0198     // If not merging, then done (?????)
0199     //
0200     // if (!doMerge) return true;
0201
0202     // Initialise chosen jet algorithm.
0203     //
0204     fJetFinder = new fastjet::JetDefinition(fastjet::kt_algorithm, coneRadius);
0205     fClusJets.clear();
0206     fPtSortedJets.clear();
0207
0208     fIsInit = true;
0209
0210     return true;
0211   }
0212
0213   void JetMatchingMGFastJet::beforeHadronisation(const lhef::LHEEvent *lhee) {
0214     if (!runInitialized)
0215       throw cms::Exception("Generator|PartonShowerVeto") << "Run not initialized in JetMatchingMGFastJet" << std::endl;
0216
0217     for (int i = 0; i < 3; i++) {
0218       typeIdx[i].clear();
0219     }
0220
0221     // Sort original process final state into light/heavy jets and 'other'.
0222     // Criteria:
0223     //   1 <= ID <= 5 and massless, or ID == 21 --> light jet (typeIdx[0])
0224     //   4 <= ID <= 6 and massive               --> heavy jet (typeIdx[1])
0225     //   All else                               --> other     (typeIdx[2])
0226     //
0227     const lhef::HEPEUP &hepeup = *lhee->getHEPEUP();
0228     int idx = 2;
0229     for (int i = 0; i < hepeup.NUP; i++) {
0230       if (hepeup.ISTUP[i] < 0)
0231         continue;
0232       if (hepeup.MOTHUP[i].first != 1 && hepeup.MOTHUP[i].second != 2)
0233         continue;  // this way we skip entries that come
0234                    // from resonance decays;
0235                    // we only take those that descent
0236                    // directly from "incoming partons"
0237       idx = 2;
0238       if (hepeup.IDUP[i] == ID_GLUON || (std::abs(hepeup.IDUP[i]) <= nQmatch))  // light jet
0239                                                                                 // light jet
0240         idx = 0;
0241       else if (std::abs(hepeup.IDUP[i]) > nQmatch && std::abs(hepeup.IDUP[i]) <= ID_TOP)  // heavy jet
0242         idx = 1;
0243       // Store
0244       typeIdx[idx].push_back(i);
0245     }
0246
0247     // NOTE: In principle, I should use exclusive, inclusive, or soup !!!
0248     // should be like this:
0249     if (soup) {
0250       int NPartons = typeIdx[0].size();
0251       fExcLocal = (NPartons < nJetMax);
0252     } else
0253       fExcLocal = exclusive;
0254
0255     return;
0256   }
0257
0258   void JetMatchingMGFastJet::init(const lhef::LHERunInfo *runInfo) {
0259     if (fIsInit)
0260       return;
0261
0262     // read MGFastJet run card
0263
0264     std::map<std::string, std::string> parameters;
0265
0266     std::vector<std::string> header = runInfo->findHeader("MGRunCard");
0267     if (header.empty())
0268       throw cms::Exception("Generator|PartonShowerVeto") << "In order to use MGFastJet jet matching, "
0269                                                             "the input file has to contain the corresponding "
0270                                                             "MadGraph headers."
0271                                                          << std::endl;
0272
0273     mgParams = parseHeader(header);
0274
0275     // set variables in common block
0276
0277     std::vector<Param> params;
0278     std::vector<Param> values;
0279     for (std::map<std::string, std::string>::const_iterator iter = mgParams.begin(); iter != mgParams.end(); ++iter) {
0280       params.push_back(" " + iter->first);
0281       values.push_back(iter->second);
0282     }
0283
0284     // set MG matching parameters
0285
0286     uppriv_.ickkw = getParameter<int>("ickkw", 0);
0287     memain_.mektsc = getParameter<int>("ktscheme", 0);
0288
0289     header = runInfo->findHeader("MGParamCMS");
0290
0291     std::map<std::string, std::string> mgInfoCMS = parseHeader(header);
0292
0293     for (std::map<std::string, std::string>::const_iterator iter = mgInfoCMS.begin(); iter != mgInfoCMS.end(); ++iter) {
0294       std::cout << "mgInfoCMS: " << iter->first << " " << iter->second << std::endl;
0295     }
0296
0297     updateOrDie(mgInfoCMS, memain_.etaclmax, "etaclmax");
0298     updateOrDie(mgInfoCMS, memain_.qcut, "qcut");
0299     updateOrDie(mgInfoCMS, memain_.minjets, "minjets");
0300     updateOrDie(mgInfoCMS, memain_.maxjets, "maxjets");
0301     updateOrDie(mgInfoCMS, memain_.showerkt, "showerkt");
0302     updateOrDie(mgInfoCMS, memain_.nqmatch, "nqmatch");
0303
0304     runInitialized = true;
0305
0306     initAfterBeams();
0307
0308     fIsInit = true;
0309
0310     return;
0311   }
0312
0313   int JetMatchingMGFastJet::match(const lhef::LHEEvent *partonLevel, const std::vector<fastjet::PseudoJet> *jetInput) {
0314     // Number of hard partons
0315     //
0316     int NPartons = typeIdx[0].size();
0317
0318     fClusJets.clear();
0319     fPtSortedJets.clear();
0320
0321     int ClusSeqNJets = 0;
0322
0323     fastjet::ClusterSequence ClusSequence(*jetInput, *fJetFinder);
0324
0325     if (fExcLocal) {
0326       fClusJets = ClusSequence.exclusive_jets(qCutSq);
0327     } else {
0328       fClusJets = ClusSequence.inclusive_jets(qCut);
0329     }
0330
0331     ClusSeqNJets = fClusJets.size();
0332
0333     if (ClusSeqNJets < NPartons)
0334       return LESS_JETS;
0335
0336     double localQcutSq = qCutSq;
0337
0338     if (fExcLocal)  // exclusive
0339     {
0340       if (ClusSeqNJets > NPartons)
0341         return MORE_JETS;
0342     } else  // inclusive
0343     {
0344       fPtSortedJets = fastjet::sorted_by_pt(*jetInput);
0345       localQcutSq = std::max(qCutSq, fPtSortedJets[0].pt2());
0346       fClusJets = ClusSequence.exclusive_jets(NPartons);  // override
0347       ClusSeqNJets = NPartons;
0348     }
0349
0350     if (clFact != 0)
0351       localQcutSq *= pow(clFact, 2);
0352
0353     std::vector<fastjet::PseudoJet> MatchingInput;
0354
0355     std::vector<bool> jetAssigned;
0356     jetAssigned.assign(fClusJets.size(), false);
0357
0358     int counter = 0;
0359
0360     const lhef::HEPEUP &hepeup = *partonLevel->getHEPEUP();
0361
0362     while (counter < NPartons) {
0363       MatchingInput.clear();
0364
0365       for (int i = 0; i < ClusSeqNJets; i++) {
0366         if (jetAssigned[i])
0367           continue;
0368         if (i == NPartons)
0369           break;
0370         //
0371         // this looks "awkward" but this way we do NOT pass in cluster_hist_index
0372         // which would actually indicate position in "history" of the "master" ClusSeq
0373         //
0374         fastjet::PseudoJet exjet = fClusJets[i];
0375         MatchingInput.push_back(fastjet::PseudoJet(exjet.px(), exjet.py(), exjet.pz(), exjet.e()));
0376         MatchingInput.back().set_user_index(i);
0377       }
0378
0379       int idx = typeIdx[0][counter];
0380       MatchingInput.push_back(
0381           fastjet::PseudoJet(hepeup.PUP[idx][0], hepeup.PUP[idx][1], hepeup.PUP[idx][2], hepeup.PUP[idx][3]));
0382
0383       //
0384       // in principle, one can use ClusterSequence::n_particles()
0385       //
0386       int NBefore = MatchingInput.size();
0387
0388       // Create new clustering object - which includes the 1st clustering run !!!
0389       // NOTE-1: it better be a new object for each try, or history will be "too crowded".
0390       // NOTE-2: when created, the it ALWAYS makes at least 1 clustering step, so in most
0391       //         cases at least 1 new jet object will be added; although in some cases
0392       //         no jet is added, if the system doesn't cluster to anything (for example,
0393       //         input jet(s) and parton(s) are going opposite directions)
0394       //
0395       fastjet::ClusterSequence ClusSeq(MatchingInput, *fJetFinder);
0396
0397       const std::vector<fastjet::PseudoJet> &output = ClusSeq.jets();
0398       int NClusJets = output.size() - NBefore;
0399
0400       //
0401       // JVY - I think this is the right idea:
0402       // at least 1 (one) new jet needs to be added
0403       // however, need to double check details and refine, especially for inclusive mode
0404       //
0405       //
0406       if (NClusJets < 1) {
0407         return UNMATCHED_PARTON;
0408       }
0409       //
0410       // very unlikely case but let's do it just to be safe
0411       //
0412       if (NClusJets >= NBefore) {
0413         return MORE_JETS;
0414       }
0415
0416       // Now browse history and see how close the clustering distance
0417       //
0418       // NOTE: Remember, there maybe more than one new jet in the list (for example,
0419       //       for process=2,3,...);
0420       //       in this case we take the ** first ** one that satisfies the distance/cut,
0421       //       which is ** typically ** also the best one
0422       //
0423       bool matched = false;
0424       const std::vector<fastjet::ClusterSequence::history_element> &history = ClusSeq.history();
0425
0426       // for ( unsigned int i=nBefore; i<history.size(); i++ )
0427       for (unsigned int i = NBefore; i < output.size(); i++) {
0428         int hidx = output[i].cluster_sequence_history_index();
0429         double dNext = history[hidx].dij;
0430         if (dNext < localQcutSq) {
0431           //
0432           // the way we form input, parent1 is always jet, and parent2 can be any,
0433           // but if it's a good match/cluster, it should point at the parton
0434           //
0435           int parent1 = history[hidx].parent1;
0436           int parent2 = history[hidx].parent2;
0437           if (parent1 < 0 || parent2 < 0)
0438             break;  // bad cluster, no match
0439           //
0440           // pull up jet's "global" index
0441           //
0442           int pidx = MatchingInput[parent1].user_index();
0443           jetAssigned[pidx] = true;
0444           matched = true;
0445           break;
0446         }
0447       }
0448       if (!matched) {
0449         return UNMATCHED_PARTON;
0450       }
0451
0452       counter++;
0453     }
0454
0455     // Now save some info for DJR analysis (if requested).
0456     // This mimics what is done in ME2pythia.f
0457     // Basically, NJets and matching scale for these 4 cases:
0458     // 1->0, 2->1, 3->2, and 4->3
0459     //
0460     if (fDJROutFlag > 0) {
0461       std::vector<double> MergingScale;
0462       MergingScale.clear();
0463       for (int nj = 0; nj < 4; nj++) {
0464         double dmscale2 = ClusSequence.exclusive_dmerge(nj);
0465         double dmscale = sqrt(dmscale2);
0466         MergingScale.push_back(dmscale);
0467       }
0468       fDJROutput.open("events.tree", std::ios_base::app);
0469       double dNJets = (double)NPartons;
0470       fDJROutput << " " << dNJets << " " << MergingScale[0] << " " << MergingScale[1] << " " << MergingScale[2] << " "
0471                  << MergingScale[3] << std::endl;
0472       fDJROutput.close();
0473     }
0474
0475     return NONE;
0476   }
0477
0478 }  // namespace gen