Skip to content

Commit

Permalink
Merge pull request #3328 from mkirsano/rapidityuserhooks_3
Browse files Browse the repository at this point in the history 
Generators -- Add rapidityuserhooks and some cleanup
  • Loading branch information
@ktf
ktf committed Apr 16, 2014
2 parents 6cffc70 + ef1cca2 commit d618342
Show file tree
Hide file tree
Showing 15 changed files with 217 additions and 90 deletions.
98 changes: 50 additions & 48 deletions GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc
View file
Expand Up @@ -89,11 +89,13 @@ class Pythia8Hadronizer : public BaseHadronizer, public Py8InterfaceBase {
double fBeam1PZ;
double fBeam2PZ;

// Reweight user hook
// Reweight user hooks
//
UserHooks* fReweightUserHook;
UserHooks* fReweightRapUserHook;
UserHooks* fReweightPtHatRapUserHook;

// PS matching protot6ype
// PS matching prototype
//
JetMatchingHook* fJetMatchingHook;

Expand Down Expand Up @@ -123,7 +125,7 @@ Pythia8Hadronizer::Pythia8Hadronizer(const edm::ParameterSet &params) :
comEnergy(params.getParameter<double>("comEnergy")),
LHEInputFileName(params.getUntrackedParameter<string>("LHEInputFileName","")),
fInitialState(PP),
fReweightUserHook(0),
fReweightUserHook(0),fReweightRapUserHook(0),fReweightPtHatRapUserHook(0),
fJetMatchingHook(0),
fEmissionVetoHook(0),fEmissionVetoHook1(0)
{
Expand All @@ -135,11 +137,9 @@ Pythia8Hadronizer::Pythia8Hadronizer(const edm::ParameterSet &params) :
if ( fInitialState == PP )
{
fInitialState = PPbar;
edm::LogInfo("GeneratorInterface|Pythia6Interface")
<< "Pythia6 will be initialized for PROTON-ANTIPROTON INITIAL STATE. "
<< "This is a user-request change from the DEFAULT PROTON-PROTON initial state." << std::endl;
std::cout << "Pythia6 will be initialized for PROTON-ANTIPROTON INITIAL STATE." << std::endl;
std::cout << "This is a user-request change from the DEFAULT PROTON-PROTON initial state." << std::endl;
edm::LogImportant("GeneratorInterface|Pythia8Interface")
<< "Pythia8 will be initialized for PROTON-ANTIPROTON INITIAL STATE. "
<< "This is a user-request change from the DEFAULT PROTON-PROTON initial state.";
}
else
{
Expand All @@ -151,11 +151,9 @@ Pythia8Hadronizer::Pythia8Hadronizer(const edm::ParameterSet &params) :
if ( fInitialState == PP )
{
fInitialState = ElectronPositron;
edm::LogInfo("GeneratorInterface|Pythia6Interface")
<< "Pythia6 will be initialized for ELECTRON-POSITRON INITIAL STATE. "
<< "This is a user-request change from the DEFAULT PROTON-PROTON initial state." << std::endl;
std::cout << "Pythia6 will be initialized for ELECTRON-POSITRON INITIAL STATE." << std::endl;
std::cout << "This is a user-request change from the DEFAULT PROTON-PROTON initial state." << std::endl;
edm::LogInfo("GeneratorInterface|Pythia8Interface")
<< "Pythia8 will be initialized for ELECTRON-POSITRON INITIAL STATE. "
<< "This is a user-request change from the DEFAULT PROTON-PROTON initial state.";
}
else
{
Expand Down Expand Up @@ -188,6 +186,34 @@ Pythia8Hadronizer::Pythia8Hadronizer(const edm::ParameterSet &params) :
//
if( params.exists( "reweightGen" ) )
fReweightUserHook = new PtHatReweightUserHook();
if( params.exists( "reweightGenRap" ) )
{
edm::LogInfo("Pythia8Interface") << "Start setup for reweightGenRap";
edm::ParameterSet rgrParams =
params.getParameter<edm::ParameterSet>("reweightGenRap");
fReweightRapUserHook =
new RapReweightUserHook(rgrParams.getParameter<std::string>("yLabSigmaFunc"),
rgrParams.getParameter<double>("yLabPower"),
rgrParams.getParameter<std::string>("yCMSigmaFunc"),
rgrParams.getParameter<double>("yCMPower"),
rgrParams.getParameter<double>("pTHatMin"),
rgrParams.getParameter<double>("pTHatMax"));
edm::LogInfo("Pythia8Interface") << "End setup for reweightGenRap";
}
if( params.exists( "reweightGenPtHatRap" ) )
{
edm::LogInfo("Pythia8Interface") << "Start setup for reweightGenPtHatRap";
edm::ParameterSet rgrParams =
params.getParameter<edm::ParameterSet>("reweightGenPtHatRap");
fReweightPtHatRapUserHook =
new PtHatRapReweightUserHook(rgrParams.getParameter<std::string>("yLabSigmaFunc"),
rgrParams.getParameter<double>("yLabPower"),
rgrParams.getParameter<std::string>("yCMSigmaFunc"),
rgrParams.getParameter<double>("yCMPower"),
rgrParams.getParameter<double>("pTHatMin"),
rgrParams.getParameter<double>("pTHatMax"));
edm::LogInfo("Pythia8Interface") << "End setup for reweightGenPtHatRap";
}

if( params.exists( "useUserHook" ) )
throw edm::Exception(edm::errors::Configuration,"Pythia8Interface")
Expand Down Expand Up @@ -241,28 +267,22 @@ Pythia8Hadronizer::Pythia8Hadronizer(const edm::ParameterSet &params) :

int NHooks=0;
if(fReweightUserHook) NHooks++;
if(fReweightRapUserHook) NHooks++;
if(fReweightPtHatRapUserHook) NHooks++;
if(fJetMatchingHook) NHooks++;
if(fEmissionVetoHook) NHooks++;
if(fEmissionVetoHook1) NHooks++;
if(NHooks > 1)
throw edm::Exception(edm::errors::Configuration,"Pythia8Interface")
<<" Too many User Hooks. \n Please choose one from: reweightGen, jetMatching, emissionVeto \n";

<<" Too many User Hooks. \n Please choose one from: reweightGen, reweightGenRap, reweightGenPtHatRap, jetMatching, emissionVeto, emissionVeto1 \n";
if(fReweightUserHook) fMasterGen->setUserHooksPtr(fReweightUserHook);
if(fReweightRapUserHook) fMasterGen->setUserHooksPtr(fReweightRapUserHook);
if(fReweightPtHatRapUserHook) fMasterGen->setUserHooksPtr(fReweightPtHatRapUserHook);
if(fJetMatchingHook) fMasterGen->setUserHooksPtr(fJetMatchingHook);
if(fEmissionVetoHook || fEmissionVetoHook1) {
cout << "Turning on Emission Veto Hook";
if(fEmissionVetoHook1) cout << " 1";
cout << endl;
int nversion = (int)(1000.*(fMasterGen->settings.parm("Pythia:versionNumber") - 8.));
if(nversion < 157) {
cout << "obsolete pythia8 version for this Emission Veto code" << endl;
cout << "Please update pythia8 version using the instructions here:" << endl;
cout << "https://twiki.cern.ch/twiki/bin/view/CMS/Pythia8Interface" << endl;
cout << "or try to use tag V00-01-28 of this interface" << endl;
throw edm::Exception(edm::errors::Configuration,"Pythia8Interface")
<<" Obsolete pythia8 version for this Emission Veto code\n";
}
std::cout << "Turning on Emission Veto Hook";
if(fEmissionVetoHook1) std::cout << " 1";
std::cout << std::endl;
if(fEmissionVetoHook) fMasterGen->setUserHooksPtr(fEmissionVetoHook);
if(fEmissionVetoHook1) fMasterGen->setUserHooksPtr(fEmissionVetoHook1);
}
Expand Down Expand Up @@ -325,13 +345,11 @@ bool Pythia8Hadronizer::initializeForExternalPartons()

std::cout << "Initializing for external partons" << std::endl;

// pythiaEvent = &pythia->event;

if(LHEInputFileName != string()) {

cout << endl;
cout << "LHE Input File Name = " << LHEInputFileName << endl;
cout << endl;
std::cout << std::endl;
std::cout << "LHE Input File Name = " << LHEInputFileName << std::endl;
std::cout << std::endl;
fMasterGen->init(LHEInputFileName);

} else {
Expand Down Expand Up @@ -392,20 +410,6 @@ bool Pythia8Hadronizer::initializeForExternalPartons()
return true;
}

#if 0
// naive Pythia8 HepMC status fixup
static int getStatus(const HepMC::GenParticle *p)
{
int status = p->status();
if (status > 0)
return status;
else if (status > -30 && status < 0)
return 3;
else
return 2;
}
#endif


void Pythia8Hadronizer::statistics()
{
Expand Down Expand Up @@ -489,7 +493,6 @@ bool Pythia8Hadronizer::residualDecay()
if ( !fDecayer->event[nentries-1].mayDecay() ) continue;
fDecayer->next();
int nentries1 = fDecayer->event.size();
// fDecayer->event.list(std::cout);
if ( nentries1 <= nentries ) continue; //same number of particles, no decays...

part->set_status(2);
Expand Down Expand Up @@ -524,7 +527,6 @@ bool Pythia8Hadronizer::residualDecay()
}

event().get()->add_vertex( DecVtx );
// fCurrentEventState->add_vertex( DecVtx );

}
}
Expand Down
94 changes: 94 additions & 0 deletions GeneratorInterface/Pythia8Interface/plugins/ReweightUserHooks.h
View file
@@ -1,4 +1,5 @@
#include "Pythia8/Pythia.h"
#include "TF1.h"

class PtHatReweightUserHook : public Pythia8::UserHooks
{
Expand All @@ -24,3 +25,96 @@ class PtHatReweightUserHook : public Pythia8::UserHooks
private:
double pt, power;
};


class RapReweightUserHook : public Pythia8::UserHooks
{
public:
RapReweightUserHook(const std::string& _yLabsigma_func, double _yLab_power,
const std::string& _yCMsigma_func, double _yCM_power,
double _pTHatMin, double _pTHatMax) :
yLabsigma_func(_yLabsigma_func), yCMsigma_func(_yCMsigma_func),
yLab_power(_yLab_power), yCM_power(_yCM_power),
pTHatMin(_pTHatMin), pTHatMax(_pTHatMax)
{
// empirical parametrizations defined in configuration file
yLabsigma = TF1("yLabsigma", yLabsigma_func.c_str(), pTHatMin, pTHatMax);
yCMsigma = TF1("yCMsigma", yLabsigma_func.c_str(), pTHatMin, pTHatMax);
}
virtual ~RapReweightUserHook() {}

virtual bool canBiasSelection() { return true; }

virtual double biasSelectionBy(const Pythia8::SigmaProcess* sigmaProcessPtr,
const Pythia8::PhaseSpace* phaseSpacePtr, bool inEvent)
{
//the variable selBias of the base class should be used;
if ((sigmaProcessPtr->nFinal() == 2)) {
double x1 = phaseSpacePtr->x1();
double x2 = phaseSpacePtr->x2();
double yLab = 0.5*log(x1/x2);
double yCM = 0.5*log( phaseSpacePtr->tHat() / phaseSpacePtr->uHat() );
double pTHat = phaseSpacePtr->pTHat();
double sigmaLab = yLabsigma.Eval(pTHat);
double sigmaCM = yCMsigma.Eval(pTHat);
// empirical reweighting function
selBias = exp( pow(fabs(yLab),yLab_power)/(2*sigmaLab*sigmaLab) +
pow(fabs(yCM),yCM_power)/(2*sigmaCM*sigmaCM) );
return selBias;
}
selBias = 1.;
return selBias;
}

private:
std::string yLabsigma_func, yCMsigma_func;
double yLab_power, yCM_power, pTHatMin, pTHatMax;
TF1 yLabsigma, yCMsigma;
};


class PtHatRapReweightUserHook : public Pythia8::UserHooks
{
public:
PtHatRapReweightUserHook(const std::string& _yLabsigma_func, double _yLab_power,
const std::string& _yCMsigma_func, double _yCM_power,
double _pTHatMin, double _pTHatMax,
double _pt = 15, double _power = 4.5) :
yLabsigma_func(_yLabsigma_func), yCMsigma_func(_yCMsigma_func),
yLab_power(_yLab_power), yCM_power(_yCM_power),
pTHatMin(_pTHatMin), pTHatMax(_pTHatMax), pt(_pt), power(_power)
{
// empirical parametrizations defined in configuration file
yLabsigma = TF1("yLabsigma", yLabsigma_func.c_str(), pTHatMin, pTHatMax);
yCMsigma = TF1("yCMsigma", yLabsigma_func.c_str(), pTHatMin, pTHatMax);
}
virtual ~PtHatRapReweightUserHook() {}

virtual bool canBiasSelection() { return true; }

virtual double biasSelectionBy(const Pythia8::SigmaProcess* sigmaProcessPtr,
const Pythia8::PhaseSpace* phaseSpacePtr, bool inEvent)
{
//the variable selBias of the base class should be used;
if ((sigmaProcessPtr->nFinal() == 2)) {
double x1 = phaseSpacePtr->x1();
double x2 = phaseSpacePtr->x2();
double yLab = 0.5*log(x1/x2);
double yCM = 0.5*log( phaseSpacePtr->tHat() / phaseSpacePtr->uHat() );
double pTHat = phaseSpacePtr->pTHat();
double sigmaLab = yLabsigma.Eval(pTHat);
double sigmaCM = yCMsigma.Eval(pTHat);
// empirical reweighting function
selBias = pow(pTHat / pt, power) * exp( pow(fabs(yLab),yLab_power)/(2*sigmaLab*sigmaLab) +
pow(fabs(yCM),yCM_power)/(2*sigmaCM*sigmaCM) );
return selBias;
}
selBias = 1.;
return selBias;
}

private:
std::string yLabsigma_func, yCMsigma_func;
double yLab_power, yCM_power, pTHatMin, pTHatMax, pt, power;
TF1 yLabsigma, yCMsigma;
};
5 changes: 0 additions & 5 deletions GeneratorInterface/Pythia8Interface/test/Py8EGun_Z2tau_cfg.py
View file
Expand Up @@ -47,12 +47,7 @@
)
)

# The following three lines reduce the clutter of repeated printouts
# of the same exception message.
process.load("FWCore.MessageLogger.MessageLogger_cfi")
process.MessageLogger.destinations = ['cerr']
process.MessageLogger.statistics = []
process.MessageLogger.fwkJobReports = []

process.RandomNumberGeneratorService = cms.Service("RandomNumberGeneratorService",
generator = cms.PSet(
Expand Down
17 changes: 9 additions & 8 deletions GeneratorInterface/Pythia8Interface/test/Py8Had_MGFastJet_cfg.py
View file
Expand Up @@ -57,14 +57,15 @@
)
)

#process.MessageLogger = cms.Service("MessageLogger",
# cout = cms.untracked.PSet(
# default = cms.untracked.PSet(
# limit = cms.untracked.int32(0)
# )
# ),
# destinations = cms.untracked.vstring('cout')
#)
process.load("FWCore.MessageLogger.MessageLogger_cfi")
process.MessageLogger = cms.Service("MessageLogger",
cout = cms.untracked.PSet(
default = cms.untracked.PSet(
limit = cms.untracked.int32(2)
)
),
destinations = cms.untracked.vstring('cout')
)

process.RandomNumberGeneratorService = cms.Service("RandomNumberGeneratorService",
generator = cms.PSet(
Expand Down
17 changes: 9 additions & 8 deletions GeneratorInterface/Pythia8Interface/test/Py8Had_mgmatching_cfg.py
View file
Expand Up @@ -58,14 +58,15 @@
)
)

#process.MessageLogger = cms.Service("MessageLogger",
# cout = cms.untracked.PSet(
# default = cms.untracked.PSet(
# limit = cms.untracked.int32(0)
# )
# ),
# destinations = cms.untracked.vstring('cout')
#)
process.load("FWCore.MessageLogger.MessageLogger_cfi")
process.MessageLogger = cms.Service("MessageLogger",
cout = cms.untracked.PSet(
default = cms.untracked.PSet(
limit = cms.untracked.int32(2)
)
),
destinations = cms.untracked.vstring('cout')
)

process.RandomNumberGeneratorService = cms.Service("RandomNumberGeneratorService",
generator = cms.PSet(
Expand Down
13 changes: 8 additions & 5 deletions GeneratorInterface/Pythia8Interface/test/Py8JetGun_cfg.py
View file
Expand Up @@ -33,12 +33,15 @@
)
)

# The following three lines reduce the clutter of repeated printouts
# of the same exception message.
process.load("FWCore.MessageLogger.MessageLogger_cfi")
process.MessageLogger.destinations = ['cerr']
process.MessageLogger.statistics = []
process.MessageLogger.fwkJobReports = []
process.MessageLogger = cms.Service("MessageLogger",
cout = cms.untracked.PSet(
default = cms.untracked.PSet(
limit = cms.untracked.int32(2)
)
),
destinations = cms.untracked.vstring('cout')
)

process.RandomNumberGeneratorService = cms.Service("RandomNumberGeneratorService",
generator = cms.PSet(
Expand Down
11 changes: 8 additions & 3 deletions GeneratorInterface/Pythia8Interface/test/Py8_Z2tautau_cfg.py
View file
Expand Up @@ -92,9 +92,14 @@
# The following three lines reduce the clutter of repeated printouts
# of the same exception message.
process.load("FWCore.MessageLogger.MessageLogger_cfi")
process.MessageLogger.destinations = ['cerr']
process.MessageLogger.statistics = []
process.MessageLogger.fwkJobReports = []
process.MessageLogger = cms.Service("MessageLogger",
cout = cms.untracked.PSet(
default = cms.untracked.PSet(
limit = cms.untracked.int32(2)
)
),
destinations = cms.untracked.vstring('cout')
)

process.RandomNumberGeneratorService = cms.Service("RandomNumberGeneratorService",
generator = cms.PSet(
Expand Down

0 comments on commit d618342

Please sign in to comment.