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