Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Redesign jet matching mg fast jet #1708

Merged
merged 2 commits into from Jan 15, 2014
Merged

Redesign jet matching mg fast jet #1708

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
4c9d0b5
Redesign JetMatchingFastJet
mkirsano Dec 6, 2013
263f1ae
Redesign JetMatchingFastJet
mkirsano Dec 6, 2013
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
62 changes: 38 additions & 24 deletions GeneratorInterface/PartonShowerVeto/interface/JetMatchingMGFastJet.h
View file
Open in desktop
@@ -1,7 +1,6 @@
#ifndef gen_JetMatchingMGFastJet_h
#define gen_JetMatchingMGFastJet_h


//
// Julia V. Yarba, Feb.10, 2013
//
Expand All @@ -10,11 +9,8 @@
// somewhat differently, and is also using FastJet package
// instead of Pythia8's native SlowJet
//
// At this point, we inherit from JetMatchingMadgraph,
// mainly to use parameters input machinery
//

#include "GeneratorInterface/PartonShowerVeto/interface/JetMatchingMadgraph.h"
#include "GeneratorInterface/PartonShowerVeto/interface/JetMatching.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "SimDataFormats/GeneratorProducts/interface/LHECommonBlocks.h"
Expand All @@ -33,43 +29,57 @@

namespace gen
{
class JetMatchingMGFastJet : public JetMatchingMadgraph
class JetMatchingMGFastJet : public JetMatching
{

public:

JetMatchingMGFastJet(const edm::ParameterSet& params) : JetMatchingMadgraph(params),
fJetFinder(0),
fIsInit(false)
{ fDJROutFlag = params.getParameter<int>("outTree_flag"); }
JetMatchingMGFastJet(const edm::ParameterSet& params);

~JetMatchingMGFastJet() { if (fJetFinder) delete fJetFinder; }

const std::vector<int>* getPartonList() { return typeIdx; }

protected:
virtual void init( const lhef::LHERunInfo* runInfo ) { if (fIsInit) return; JetMatchingMadgraph::init(runInfo);
initAfterBeams(); fIsInit=true; return; }
bool initAfterBeams();
void beforeHadronisation( const lhef::LHEEvent* );
void beforeHadronisationExec() { return; }

virtual void init( const lhef::LHERunInfo* runInfo );

virtual bool initAfterBeams();
virtual void beforeHadronisation( const lhef::LHEEvent* );
virtual void beforeHadronisationExec() { return; }

int match( const lhef::LHEEvent* partonLevel, const std::vector<fastjet::PseudoJet>* jetInput );

virtual int match( const lhef::LHEEvent* partonLevel, const std::vector<fastjet::PseudoJet>* jetInput );

virtual double getJetEtaMax() const;

private:

//parameters staff from Madgraph

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

template<typename T>
static void updateOrDie(
const std::map<std::string, std::string> &params,
T &param, const std::string &name);

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

// ----------------------------

enum vetoStatus { NONE, LESS_JETS, MORE_JETS, HARD_JET, UNMATCHED_PARTON };
enum partonTypes { ID_TOP=6, ID_GLUON=21, ID_PHOTON=22 };

double qCut, qCutSq;
double clFact;
int nQmatch;

//
// these 2 below are legacy but keep here for now
//
//int ktScheme;
//int showerKt;

// Master switch for merging
bool doMerge;

Expand All @@ -81,13 +91,17 @@ class JetMatchingMGFastJet : public JetMatchingMadgraph
double coneRadius, etaJetMax ;

// Merging procedure control flag(s)
// (there're also inclusive, exclusive, and soup/auto in JetMatchingMadgraph)
// (there're also inclusive, exclusive, and soup/auto in JetMatchingMGFastJet)
//
bool fExcLocal; // this is similar to memaev_.iexc

// Sort final-state of incoming process into light/heavy jets and 'other'
std::vector < int > typeIdx[3];

bool runInitialized;
bool soup;
bool exclusive;

//
// FastJets tool(s)
//
Expand Down
240 changes: 238 additions & 2 deletions GeneratorInterface/PartonShowerVeto/src/JetMatchingMGFastJet.cc
View file
Open in desktop
@@ -1,4 +1,3 @@

#include "GeneratorInterface/PartonShowerVeto/interface/JetMatchingMGFastJet.h"
#include "FWCore/Utilities/interface/Exception.h"
#include <iomanip>
Expand All @@ -8,11 +7,30 @@
#include "GeneratorInterface/LHEInterface/interface/LHERunInfo.h"
#include "GeneratorInterface/LHEInterface/interface/LHEEvent.h"

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


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 struct UPPRIV {
int lnhin, lnhout;
int mscal, ievnt;
Expand All @@ -27,8 +45,165 @@ extern "C" {
bool nosingrad,jetprocs;
} memain_;

extern struct OUTTREE {
int flag;
} outtree_;

}


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

template<>
std::string JetMatchingMGFastJet::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 JetMatchingMGFastJet::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 JetMatchingMGFastJet::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 JetMatchingMGFastJet::getParameter(const std::string &var,
const T &defValue) const
{
return getParameter(mgParams, var, defValue);
}


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 JetMatchingMGFastJet::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;
}


JetMatchingMGFastJet::JetMatchingMGFastJet(const edm::ParameterSet &params) :
JetMatching(params),
runInitialized(false),
fJetFinder(0),
fIsInit(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")
<< "MGFastJet 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;

fDJROutFlag = params.getParameter<int>("outTree_flag");
}


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


bool JetMatchingMGFastJet::initAfterBeams()
{

Expand Down Expand Up @@ -78,7 +253,7 @@ void JetMatchingMGFastJet::beforeHadronisation( const lhef::LHEEvent* lhee )

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

for (int i = 0; i < 3; i++)
Expand Down Expand Up @@ -125,6 +300,67 @@ void JetMatchingMGFastJet::beforeHadronisation( const lhef::LHEEvent* lhee )

}


void JetMatchingMGFastJet::init(const lhef::LHERunInfo* runInfo)
{
if (fIsInit) return;

// read MGFastJet 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 MGFastJet 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");

runInitialized = true;

initAfterBeams();

fIsInit=true;

return;
}


int JetMatchingMGFastJet::match( const lhef::LHEEvent* partonLevel,
const std::vector<fastjet::PseudoJet>* jetInput )
{
Expand Down