/data/refman/pasoursint/CMSSW_5_3_1/src/GeneratorInterface/LHEInterface/plugins/JetMatchingMadgraph.cc

Go to the documentation of this file.

#include <functional>
#include <algorithm>
#include <iostream>
#include <sstream>
#include <vector>
#include <memory>
#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/LHEInterface/interface/JetMatching.h"

namespace lhef {

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, clfact;
                int     maxjets, minjets, iexcfile, ktsche;
                int     nexcres, excres[30];
        } memain_;

        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"

class JetMatchingMadgraph : public JetMatching {
    public:
        JetMatchingMadgraph(const edm::ParameterSet &params);
        ~JetMatchingMadgraph();

    private:
        void init(const boost::shared_ptr<LHERunInfo> &runInfo);
        void beforeHadronisation(const boost::shared_ptr<LHEEvent> &event);

        double match(const HepMC::GenEvent *partonLevel,
                     const HepMC::GenEvent *finalState,
                     bool showeredFinalState);

        std::set<std::string> capabilities() const;

        template<typename T>
        static T parseParameter(const std::string &value);
        template<typename T>
        T getParameter(const std::string &var, const T &defValue = T()) const;

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

        bool                                    runInitialized;
        bool                                    eventInitialized;
        bool                                    soup;
        bool                                    exclusive;
};

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 == "1" || value == ".TRUE.";
}

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

} // namespace lhef

using namespace lhef;

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_.iexcfile = 0;
        memain_.ktsche = 0;
        memain_.etaclmax = params.getParameter<double>("etaclmax");
        memain_.qcut = params.getParameter<double>("qcut");
        memain_.minjets = params.getParameter<int>("minjets");
        memain_.maxjets = params.getParameter<int>("maxjets");
}

JetMatchingMadgraph::~JetMatchingMadgraph()
{
}

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

// implements the Madgraph method - use ME2pythia.f

void JetMatchingMadgraph::init(const boost::shared_ptr<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|LHEInterface")
                        << "In order to use MadGraph jet matching, "
                           "the input file has to contain the corresponding "
                           "MadGraph headers." << std::endl;

        mgParams.clear();
        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));

                mgParams[var] = value;
        }

        // 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);

        // run Fortran initialization code

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

void JetMatchingMadgraph::beforeHadronisation(
                                const boost::shared_ptr<LHEEvent> &event)
{
        if (!runInitialized)
                throw cms::Exception("Generator|LHEInterface")
                        << "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 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;
}

double JetMatchingMadgraph::match(const HepMC::GenEvent *partonLevel,
                                  const HepMC::GenEvent *finalState,
                                  bool showeredFinalState)
{
        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);
        eventInitialized = false;

        return veto ? 0.0 : 1.0;
}

DEFINE_LHE_JETMATCHING_PLUGIN(JetMatchingMadgraph);