JetMatchingMadgraph.cc

Go to the documentation of this file.

#include <functional>
#include <algorithm>
#include <iostream>
#include <sstream>
#include <vector>
#include <memory>
#include <cstring>
#include <string>
#include <cctype>
#include <map>
#include <set>

#include <boost/shared_ptr.hpp>
#include <boost/algorithm/string/trim.hpp>

// #include <HepMC/GenEvent.h>

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "GeneratorInterface/LHEInterface/interface/LHERunInfo.h"
#include "GeneratorInterface/LHEInterface/interface/LHEEvent.h"
#include "GeneratorInterface/PartonShowerVeto/interface/JetMatchingMadgraph.h"

namespace gen {

extern "C" {
    #define PARAMLEN 20
    namespace {
        struct Param {
            Param(const std::string &str)
            {
                int len = std::min(PARAMLEN,
                                   (int)str.length());
                std::memcpy(value, str.c_str(), len);
                std::memset(value + len, ' ', PARAMLEN - len);
            }

            char value[PARAMLEN];
        };
    }
    extern void mginit_(int *npara, Param *params, Param *values);
    extern void mgevnt_(void);
    extern void mgveto_(int *veto);

    extern struct UPPRIV {
        int lnhin, lnhout;
        int mscal, ievnt;
        int ickkw, iscale;
    } uppriv_;

    extern struct MEMAIN {
        double  etcjet, rclmax, etaclmax, qcut, showerkt, clfact;
        int maxjets, minjets, iexcfile, ktsche;
        int mektsc,nexcres, excres[30];
            int nqmatch,excproc,iexcproc[1000],iexcval[1000];
                bool    nosingrad,jetprocs;
    } memain_;
    
    extern struct OUTTREE {
        int flag;
    } outtree_;
    
    extern struct MEMAEV {
        double  ptclus[20];
        int nljets, iexc, ifile;
    } memaev_;

    extern struct PYPART {
        int npart, npartd, ipart[1000];
        double  ptpart[1000];
    } pypart_;
} // extern "C"


template<typename T>
T JetMatchingMadgraph::parseParameter(const std::string &value)
{
    std::istringstream ss(value);
    T result;
    ss >> result;
    return result;
}

template<>
std::string JetMatchingMadgraph::parseParameter(const std::string &value)
{   
    std::string result;
    if (!result.empty() && result[0] == '\'')
        result = result.substr(1);
    if (!result.empty() && result[result.length() - 1] == '\'')
        result.resize(result.length() - 1);
    return result;
}

template<>
bool JetMatchingMadgraph::parseParameter(const std::string &value_)
{
    std::string value(value_);
    std::transform(value.begin(), value.end(),
                   value.begin(), (int(*)(int))std::toupper);
    return value == "T" || value == "Y" || value=="True" ||
           value == "1" || value == ".TRUE.";
}

template<typename T>
T JetMatchingMadgraph::getParameter(
            const std::map<std::string, std::string> &params,
            const std::string &var, const T &defValue)
{
    std::map<std::string, std::string>::const_iterator pos =
                            params.find(var);
    if (pos == params.end())
        return defValue;
    return parseParameter<T>(pos->second);
}

template<typename T>
T JetMatchingMadgraph::getParameter(const std::string &var,
                                    const T &defValue) const
{ 
    return getParameter(mgParams, var, defValue); 
}

JetMatchingMadgraph::JetMatchingMadgraph(const edm::ParameterSet &params) :
    JetMatching(params),
    runInitialized(false)
{
    std::string mode = params.getParameter<std::string>("mode");
    if (mode == "inclusive") {
        soup = false;
        exclusive = false;
    } else if (mode == "exclusive") {
        soup = false;
        exclusive = true;
    } else if (mode == "auto")
        soup = true;
    else
        throw cms::Exception("Generator|LHEInterface")
            << "Madgraph jet matching scheme requires \"mode\" "
               "parameter to be set to either \"inclusive\", "
               "\"exclusive\" or \"auto\"." << std::endl;

    memain_.etcjet = 0.;
    memain_.rclmax = 0.0;
    memain_.clfact = 0.0;
    memain_.ktsche = 0.0;
    memain_.etaclmax = params.getParameter<double>("MEMAIN_etaclmax");
    memain_.qcut = params.getParameter<double>("MEMAIN_qcut");
    memain_.minjets = params.getParameter<int>("MEMAIN_minjets");
    memain_.maxjets = params.getParameter<int>("MEMAIN_maxjets");
    memain_.showerkt = params.getParameter<double>("MEMAIN_showerkt");
    memain_.nqmatch = params.getParameter<int>("MEMAIN_nqmatch");
    outtree_.flag = params.getParameter<int>("outTree_flag");
    std::string list_excres = params.getParameter<std::string>("MEMAIN_excres");
    std::vector<std::string> elems;
    std::stringstream ss(list_excres);
    std::string item;
    int index=0;
    while(std::getline(ss, item, ',')) {
            elems.push_back(item);
        memain_.excres[index]=std::atoi(item.c_str());
        index++;
    }
    memain_.nexcres=index;

    

}

JetMatchingMadgraph::~JetMatchingMadgraph()
{
}

double JetMatchingMadgraph::getJetEtaMax() const
{   
   return memain_.etaclmax;
}

std::set<std::string> JetMatchingMadgraph::capabilities() const
{
    std::set<std::string> result;
    result.insert("psFinalState");
    result.insert("hepevt");
    result.insert("pythia6");
    return result;
}

static std::map<std::string, std::string>
parseHeader(const std::vector<std::string> &header)
{
    std::map<std::string, std::string> params;

    for(std::vector<std::string>::const_iterator iter = header.begin();
        iter != header.end(); ++iter) {
        std::string line = *iter;
        if (line.empty() || line[0] == '#')
            continue;

        std::string::size_type pos = line.find('!');
        if (pos != std::string::npos)
            line.resize(pos);

        pos = line.find('=');
        if (pos == std::string::npos)
            continue;

        std::string var =
            boost::algorithm::trim_copy(line.substr(pos + 1));
        std::string value = 
            boost::algorithm::trim_copy(line.substr(0, pos));

        params[var] = value;
    }

    return params;
}

template<typename T>
void JetMatchingMadgraph::updateOrDie(
            const std::map<std::string, std::string> &params,
            T &param, const std::string &name)
{
    if (param < 0){
        param = getParameter(params, name, param);
    }
    if (param < 0)
        throw cms::Exception("Generator|PartonShowerVeto")
            << "The MGParamCMS header does not specify the jet "
               "matching parameter \"" << name << "\", but it "
               "is requested by the CMSSW configuration."
            << std::endl;
}

// implements the Madgraph method - use ME2pythia.f

void JetMatchingMadgraph::init(const lhef::LHERunInfo* runInfo)
{
    // read MadGraph run card

    std::map<std::string, std::string> parameters;

    std::vector<std::string> header = runInfo->findHeader("MGRunCard");
    if (header.empty())
        throw cms::Exception("Generator|PartonShowerVeto")
            << "In order to use MadGraph jet matching, "
               "the input file has to contain the corresponding "
               "MadGraph headers." << std::endl;

    mgParams = parseHeader(header);

    // set variables in common block

    std::vector<Param> params;
    std::vector<Param> values;
    for(std::map<std::string, std::string>::const_iterator iter =
            mgParams.begin(); iter != mgParams.end(); ++iter) {
        params.push_back(" " + iter->first);
        values.push_back(iter->second);

    }

    // set MG matching parameters

    uppriv_.ickkw = getParameter<int>("ickkw", 0);
    memain_.mektsc = getParameter<int>("ktscheme", 0);

    header = runInfo->findHeader("MGParamCMS");
    
    std::map<std::string, std::string> mgInfoCMS = parseHeader(header);

    for(std::map<std::string, std::string>::const_iterator iter =
                        mgInfoCMS.begin(); iter != mgInfoCMS.end(); ++iter) {
                        std::cout<<"mgInfoCMS: "<<iter->first<<" "<<iter->second<<std::endl;

        }


    updateOrDie(mgInfoCMS, memain_.etaclmax, "etaclmax");
    updateOrDie(mgInfoCMS, memain_.qcut, "qcut");
    updateOrDie(mgInfoCMS, memain_.minjets, "minjets");
    updateOrDie(mgInfoCMS, memain_.maxjets, "maxjets");
    updateOrDie(mgInfoCMS, memain_.showerkt, "showerkt");
    updateOrDie(mgInfoCMS, memain_.nqmatch, "nqmatch");

    // run Fortran initialization code

    int nparam = params.size();
    mginit_(&nparam, &params.front(), &values.front());
    runInitialized = true;
}

//void JetMatchingMadgraph::beforeHadronisation(
//              const boost::shared_ptr<lhef::LHEEvent> &event)

void JetMatchingMadgraph::beforeHadronisation(const lhef::LHEEvent* event)
{
    if (!runInitialized)
        throw cms::Exception("Generator|PartonShowerVeto")
            << "Run not initialized in JetMatchingMadgraph"
            << std::endl;


    if (uppriv_.ickkw) 
    {
        std::vector<std::string> comments = event->getComments();
        if (comments.size() == 1) 
        {
            std::istringstream ss(comments[0].substr(1));
            for(int i = 0; i < 1000; i++) 
            {
                double pt;
                ss >> pt;
                if (!ss.good())
                    break;
                pypart_.ptpart[i] = pt;
            }
        } 
        else 
        {
            edm::LogWarning("Generator|LHEInterface")
                << "Expected exactly one comment line per "
                   "event containing MadGraph parton scale "
                   "information."
                << std::endl;

            const lhef::HEPEUP *hepeup = event->getHEPEUP();
            for(int i = 2; i < hepeup->NUP; i++) 
            {
                double mt2 =
                    hepeup->PUP[i][0] * hepeup->PUP[i][0] +
                    hepeup->PUP[i][1] * hepeup->PUP[i][1] +
                    hepeup->PUP[i][4] * hepeup->PUP[i][4];
                pypart_.ptpart[i - 2] = std::sqrt(mt2);
            }
        }
    }

    // mgevnt_();
    eventInitialized = true;
}

void JetMatchingMadgraph::beforeHadronisationExec()
{
    mgevnt_();
    eventInitialized = true;
        return;
}

/*
int JetMatchingMadgraph::match(const HepMC::GenEvent *partonLevel,
                                  const HepMC::GenEvent *finalState,
                                  bool showeredFinalState)
*/
int JetMatchingMadgraph::match( const lhef::LHEEvent* partonLevel, const std::vector<fastjet::PseudoJet>* jetInput )
{
/*
    if (!showeredFinalState)
        throw cms::Exception("Generator|LHEInterface")
            << "MadGraph matching expected parton shower "
               "final state." << std::endl;
*/

    if (!runInitialized)
        throw cms::Exception("Generator|LHEInterface")
            << "Run not initialized in JetMatchingMadgraph"
            << std::endl;

    if (!eventInitialized)
        throw cms::Exception("Generator|LHEInterface")
            << "Event not initialized in JetMatchingMadgraph"
            << std::endl;

    if (soup)
        memaev_.iexc = (memaev_.nljets < memain_.maxjets);
    else
        memaev_.iexc = exclusive;

    int veto = 0;
    mgveto_(&veto);
    fMatchingStatus=true;
    eventInitialized = false;

    return veto;
}

} // end namespace gen