Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:07:21

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