/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/GeneratorInterface/PartonShowerVeto/src/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