PartonShowerVeto/src/JetMatchingMadgraph.cc

0001 #include <functional>
0002 #include <algorithm>
0003 #include <iostream>
0004 #include <sstream>
0005 #include <vector>
0006 #include <memory>
0007 #include <cstring>
0008 #include <string>
0009 #include <cctype>
0010 #include <map>
0011 #include <set>
0012
0013 #include <boost/algorithm/string/trim.hpp>
0014
0015 // #include <HepMC/GenEvent.h>
0016
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0019
0020 #include "GeneratorInterface/LHEInterface/interface/LHERunInfo.h"
0021 #include "GeneratorInterface/LHEInterface/interface/LHEEvent.h"
0022 #include "GeneratorInterface/PartonShowerVeto/interface/JetMatchingMadgraph.h"
0023
0024 namespace gen {
0025
0026   extern "C" {
0027 #define PARAMLEN 20
0028   namespace {
0029     struct Param {
0030       Param(const std::string &str) {
0031         int len = std::min(PARAMLEN, (int)str.length());
0032         std::memcpy(value, str.c_str(), len);
0033         std::memset(value + len, ' ', PARAMLEN - len);
0034       }
0035
0036       char value[PARAMLEN];
0037     };
0038   }  // namespace
0039   extern void mginit_(int *npara, Param *params, Param *values);
0040   extern void mgevnt_(void);
0041   extern void mgveto_(int *veto);
0042
0043   extern struct UPPRIV {
0044     int lnhin, lnhout;
0045     int mscal, ievnt;
0046     int ickkw, iscale;
0047   } uppriv_;
0048
0049   extern struct MEMAIN {
0050     double etcjet, rclmax, etaclmax, qcut, showerkt, clfact;
0051     int maxjets, minjets, iexcfile, ktsche;
0052     int mektsc, nexcres, excres[30];
0053     int nqmatch, excproc, iexcproc[1000], iexcval[1000];
0054     bool nosingrad, jetprocs;
0055   } memain_;
0056
0057   extern struct OUTTREE {
0058     int flag;
0059   } outtree_;
0060
0061   extern struct MEMAEV {
0062     double ptclus[20];
0063     int nljets, iexc, ifile;
0064   } memaev_;
0065
0066   extern struct PYPART {
0067     int npart, npartd, ipart[1000];
0068     double ptpart[1000];
0069   } pypart_;
0070   }  // extern "C"
0071
0072   template <typename T>
0073   T JetMatchingMadgraph::parseParameter(const std::string &value) {
0074     std::istringstream ss(value);
0075     T result;
0076     ss >> result;
0077     return result;
0078   }
0079
0080   template <>
0081   std::string JetMatchingMadgraph::parseParameter(const std::string &value) {
0082     std::string result;
0083     if (!result.empty() && result[0] == '\'')
0084       result = result.substr(1);
0085     if (!result.empty() && result[result.length() - 1] == '\'')
0086       result.resize(result.length() - 1);
0087     return result;
0088   }
0089
0090   template <>
0091   bool JetMatchingMadgraph::parseParameter(const std::string &value_) {
0092     std::string value(value_);
0093     std::transform(value.begin(), value.end(), value.begin(), (int (*)(int))std::toupper);
0094     return value == "T" || value == "Y" || value == "True" || value == "1" || value == ".TRUE.";
0095   }
0096
0097   template <typename T>
0098   T JetMatchingMadgraph::getParameter(const std::map<std::string, std::string> &params,
0099                                       const std::string &var,
0100                                       const T &defValue) {
0101     std::map<std::string, std::string>::const_iterator pos = params.find(var);
0102     if (pos == params.end())
0103       return defValue;
0104     return parseParameter<T>(pos->second);
0105   }
0106
0107   template <typename T>
0108   T JetMatchingMadgraph::getParameter(const std::string &var, const T &defValue) const {
0109     return getParameter(mgParams, var, defValue);
0110   }
0111
0112   JetMatchingMadgraph::JetMatchingMadgraph(const edm::ParameterSet &params)
0113       : JetMatching(params), runInitialized(false) {
0114     std::string mode = params.getParameter<std::string>("mode");
0115     if (mode == "inclusive") {
0116       soup = false;
0117       exclusive = false;
0118     } else if (mode == "exclusive") {
0119       soup = false;
0120       exclusive = true;
0121     } else if (mode == "auto")
0122       soup = true;
0123     else
0124       throw cms::Exception("Generator|LHEInterface") << "Madgraph jet matching scheme requires \"mode\" "
0125                                                         "parameter to be set to either \"inclusive\", "
0126                                                         "\"exclusive\" or \"auto\"."
0127                                                      << std::endl;
0128
0129     memain_.etcjet = 0.;
0130     memain_.rclmax = 0.0;
0131     memain_.clfact = 0.0;
0132     memain_.ktsche = 0.0;
0133     memain_.etaclmax = params.getParameter<double>("MEMAIN_etaclmax");
0134     memain_.qcut = params.getParameter<double>("MEMAIN_qcut");
0135     memain_.minjets = params.getParameter<int>("MEMAIN_minjets");
0136     memain_.maxjets = params.getParameter<int>("MEMAIN_maxjets");
0137     memain_.showerkt = params.getParameter<double>("MEMAIN_showerkt");
0138     memain_.nqmatch = params.getParameter<int>("MEMAIN_nqmatch");
0139     outtree_.flag = params.getParameter<int>("outTree_flag");
0140     std::string list_excres = params.getParameter<std::string>("MEMAIN_excres");
0141     std::vector<std::string> elems;
0142     std::stringstream ss(list_excres);
0143     std::string item;
0144     int index = 0;
0145     while (std::getline(ss, item, ',')) {
0146       elems.push_back(item);
0147       memain_.excres[index] = std::atoi(item.c_str());
0148       index++;
0149     }
0150     memain_.nexcres = index;
0151   }
0152
0153   JetMatchingMadgraph::~JetMatchingMadgraph() {}
0154
0155   double JetMatchingMadgraph::getJetEtaMax() const { return memain_.etaclmax; }
0156
0157   std::set<std::string> JetMatchingMadgraph::capabilities() const {
0158     std::set<std::string> result;
0159     result.insert("psFinalState");
0160     result.insert("hepevt");
0161     result.insert("pythia6");
0162     return result;
0163   }
0164
0165   static std::map<std::string, std::string> parseHeader(const std::vector<std::string> &header) {
0166     std::map<std::string, std::string> params;
0167
0168     for (std::vector<std::string>::const_iterator iter = header.begin(); iter != header.end(); ++iter) {
0169       std::string line = *iter;
0170       if (line.empty() || line[0] == '#')
0171         continue;
0172
0173       std::string::size_type pos = line.find('!');
0174       if (pos != std::string::npos)
0175         line.resize(pos);
0176
0177       pos = line.find('=');
0178       if (pos == std::string::npos)
0179         continue;
0180
0181       std::string var = boost::algorithm::trim_copy(line.substr(pos + 1));
0182       std::string value = boost::algorithm::trim_copy(line.substr(0, pos));
0183
0184       params[var] = value;
0185     }
0186
0187     return params;
0188   }
0189
0190   template <typename T>
0191   void JetMatchingMadgraph::updateOrDie(const std::map<std::string, std::string> &params,
0192                                         T &param,
0193                                         const std::string &name) {
0194     if (param < 0) {
0195       param = getParameter(params, name, param);
0196     }
0197     if (param < 0)
0198       throw cms::Exception("Generator|PartonShowerVeto") << "The MGParamCMS header does not specify the jet "
0199                                                             "matching parameter \""
0200                                                          << name
0201                                                          << "\", but it "
0202                                                             "is requested by the CMSSW configuration."
0203                                                          << std::endl;
0204   }
0205
0206   // implements the Madgraph method - use ME2pythia.f
0207
0208   void JetMatchingMadgraph::init(const lhef::LHERunInfo *runInfo) {
0209     // read MadGraph run card
0210
0211     std::map<std::string, std::string> parameters;
0212
0213     std::vector<std::string> header = runInfo->findHeader("MGRunCard");
0214     if (header.empty())
0215       throw cms::Exception("Generator|PartonShowerVeto") << "In order to use MadGraph jet matching, "
0216                                                             "the input file has to contain the corresponding "
0217                                                             "MadGraph headers."
0218                                                          << std::endl;
0219
0220     mgParams = parseHeader(header);
0221
0222     // set variables in common block
0223
0224     std::vector<Param> params;
0225     std::vector<Param> values;
0226     for (std::map<std::string, std::string>::const_iterator iter = mgParams.begin(); iter != mgParams.end(); ++iter) {
0227       params.push_back(" " + iter->first);
0228       values.push_back(iter->second);
0229     }
0230
0231     // set MG matching parameters
0232
0233     uppriv_.ickkw = getParameter<int>("ickkw", 0);
0234     memain_.mektsc = getParameter<int>("ktscheme", 0);
0235
0236     header = runInfo->findHeader("MGParamCMS");
0237
0238     std::map<std::string, std::string> mgInfoCMS = parseHeader(header);
0239
0240     for (std::map<std::string, std::string>::const_iterator iter = mgInfoCMS.begin(); iter != mgInfoCMS.end(); ++iter) {
0241       std::cout << "mgInfoCMS: " << iter->first << " " << iter->second << std::endl;
0242     }
0243
0244     updateOrDie(mgInfoCMS, memain_.etaclmax, "etaclmax");
0245     updateOrDie(mgInfoCMS, memain_.qcut, "qcut");
0246     updateOrDie(mgInfoCMS, memain_.minjets, "minjets");
0247     updateOrDie(mgInfoCMS, memain_.maxjets, "maxjets");
0248     updateOrDie(mgInfoCMS, memain_.showerkt, "showerkt");
0249     updateOrDie(mgInfoCMS, memain_.nqmatch, "nqmatch");
0250
0251     // run Fortran initialization code
0252
0253     int nparam = params.size();
0254     mginit_(&nparam, &params.front(), &values.front());
0255     runInitialized = true;
0256   }
0257
0258   //void JetMatchingMadgraph::beforeHadronisation(
0259   //                const std::shared_ptr<lhef::LHEEvent> &event)
0260
0261   void JetMatchingMadgraph::beforeHadronisation(const lhef::LHEEvent *event) {
0262     if (!runInitialized)
0263       throw cms::Exception("Generator|PartonShowerVeto") << "Run not initialized in JetMatchingMadgraph" << std::endl;
0264
0265     if (uppriv_.ickkw) {
0266       std::vector<std::string> comments = event->getComments();
0267       if (comments.size() == 1) {
0268         std::istringstream ss(comments[0].substr(1));
0269         for (int i = 0; i < 1000; i++) {
0270           double pt;
0271           ss >> pt;
0272           if (!ss.good())
0273             break;
0274           pypart_.ptpart[i] = pt;
0275         }
0276       } else {
0277         edm::LogWarning("Generator|LHEInterface") << "Expected exactly one comment line per "
0278                                                      "event containing MadGraph parton scale "
0279                                                      "information."
0280                                                   << std::endl;
0281
0282         const lhef::HEPEUP *hepeup = event->getHEPEUP();
0283         for (int i = 2; i < hepeup->NUP; i++) {
0284           double mt2 = hepeup->PUP[i][0] * hepeup->PUP[i][0] + hepeup->PUP[i][1] * hepeup->PUP[i][1] +
0285                        hepeup->PUP[i][4] * hepeup->PUP[i][4];
0286           pypart_.ptpart[i - 2] = std::sqrt(mt2);
0287         }
0288       }
0289     }
0290
0291     // mgevnt_();
0292     eventInitialized = true;
0293   }
0294
0295   void JetMatchingMadgraph::beforeHadronisationExec() {
0296     mgevnt_();
0297     eventInitialized = true;
0298     return;
0299   }
0300
0301   /*
0302 int JetMatchingMadgraph::match(const HepMC::GenEvent *partonLevel,
0303                                   const HepMC::GenEvent *finalState,
0304                                   bool showeredFinalState)
0305 */
0306   int JetMatchingMadgraph::match(const lhef::LHEEvent *partonLevel, const std::vector<fastjet::PseudoJet> *jetInput) {
0307     /*
0308     if (!showeredFinalState)
0309         throw cms::Exception("Generator|LHEInterface")
0310             << "MadGraph matching expected parton shower "
0311                "final state." << std::endl;
0312 */
0313
0314     if (!runInitialized)
0315       throw cms::Exception("Generator|LHEInterface") << "Run not initialized in JetMatchingMadgraph" << std::endl;
0316
0317     if (!eventInitialized)
0318       throw cms::Exception("Generator|LHEInterface") << "Event not initialized in JetMatchingMadgraph" << std::endl;
0319
0320     if (soup)
0321       memaev_.iexc = (memaev_.nljets < memain_.maxjets);
0322     else
0323       memaev_.iexc = exclusive;
0324
0325     int veto = 0;
0326     mgveto_(&veto);
0327     fMatchingStatus = true;
0328     eventInitialized = false;
0329
0330     return veto;
0331   }
0332
0333 }  // end namespace gen