diff --git a/Decay/Radiation/QEDRadiationHandler.cc b/Decay/Radiation/QEDRadiationHandler.cc
--- a/Decay/Radiation/QEDRadiationHandler.cc
+++ b/Decay/Radiation/QEDRadiationHandler.cc
@@ -1,174 +1,174 @@
// -*- C++ -*-
//
// QEDRadiationHandler.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the QEDRadiationHandler class.
//
#include "QEDRadiationHandler.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "DecayRadiationGenerator.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Handlers/EventHandler.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "Herwig/Decay/DecayIntegrator.h"
#include "Herwig/Decay/DecayPhaseSpaceMode.h"
#include "ThePEG/PDT/DecayMode.h"
using namespace Herwig;
namespace {
/**
* A struct to order the particles in the same way as in the DecayMode's
*/
struct ParticleOrdering {
/**
* Operator for the ordering
* @param p1 The first ParticleData object
* @param p2 The second ParticleData object
*/
- bool operator()(cPDPtr p1, cPDPtr p2) {
+ bool operator()(cPDPtr p1, cPDPtr p2) const {
return abs(p1->id()) > abs(p2->id()) ||
( abs(p1->id()) == abs(p2->id()) && p1->id() > p2->id() ) ||
( p1->id() == p2->id() && p1->fullName() > p2->fullName() );
}
};
/**
* A set of ParticleData objects ordered as for the DecayMode's
*/
typedef multiset<cPDPtr,ParticleOrdering> OrderedParticles;
}
QEDRadiationHandler::QEDRadiationHandler() {
// only include electroweak gauge bosons
_decayingParticles.push_back( 22);
_decayingParticles.push_back( 23);
_decayingParticles.push_back( 24);
_decayingParticles.push_back(-24);
// only include the charged leptons
_decayProducts.push_back( 11);
_decayProducts.push_back( 13);
_decayProducts.push_back( 15);
_decayProducts.push_back(-11);
_decayProducts.push_back(-13);
_decayProducts.push_back(-15);
}
void QEDRadiationHandler::
handle(EventHandler & eh, const tPVector & tagged,
const Hint &) {
// find the potential decaying particles to be considered
set<tPPtr> parents;
for(unsigned int ix=0;ix<tagged.size();++ix) {
long id=tagged[ix]->id();
if(find(_decayProducts.begin(),_decayProducts.end(),id)!=_decayProducts.end()) {
PPtr par=tagged[ix]->parents()[0];
id=par->id();
if(tagged[ix]->parents()[0]->mass()>ZERO&&
find(_decayingParticles.begin(),_decayingParticles.end(),id)!=
_decayingParticles.end()) parents.insert(par);
}
}
set<tPPtr>::const_iterator sit;
StepPtr step=eh.currentStep();
// loop over parents
for(sit=parents.begin();sit!=parents.end();++sit) {
// extract children
ParticleVector children=(**sit).children();
// store number of children
unsigned int initsize = children.size();
// boost the decay to the parent rest frame
Boost boost = - (**sit).momentum().boostVector();
(**sit).deepBoost(boost);
// construct the tag for the decay mode to find the decayer
OrderedParticles products;
for(unsigned int ix=0;ix<children.size();++ix)
products.insert(children[ix]->dataPtr());
string tag = (**sit).dataPtr()->name() + "->";
unsigned int iprod=0;
for(OrderedParticles::const_iterator it = products.begin();
it != products.end(); ++it) {
++iprod;
tag += (**it).name();
if(iprod != initsize) tag += ",";
}
tag += ";";
tDMPtr dm = generator()->findDecayMode(tag);
tDecayIntegratorPtr decayer;
if(dm) {
tDecayerPtr dtemp = dm->decayer();
decayer = dynamic_ptr_cast<tDecayIntegratorPtr>(dtemp);
bool cc;
tPDVector ctemp;
for(unsigned int ix=0;ix<children.size();++ix)
ctemp.push_back(const_ptr_cast<tPDPtr>(children[ix]->dataPtr()));
unsigned int imode = decayer->modeNumber(cc,(**sit).dataPtr(),ctemp);
decayer->me2(imode,**sit,children,DecayIntegrator::Initialize);
}
// generate photons
ParticleVector newchildren =
_generator->generatePhotons(**sit,children,decayer);
// if photons produced add as children and to step
for(unsigned int ix=initsize;ix<newchildren.size();++ix) {
(**sit).addChild(newchildren[ix]);
step->addDecayProduct(newchildren[ix]);
}
(**sit).deepBoost(-boost);
}
}
void QEDRadiationHandler::persistentOutput(PersistentOStream & os) const {
os << _generator << _decayingParticles << _decayProducts;
}
void QEDRadiationHandler::persistentInput(PersistentIStream & is, int) {
is >> _generator >> _decayingParticles >> _decayProducts;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<QEDRadiationHandler,StepHandler>
describeHerwigQEDRadiationHandler("Herwig::QEDRadiationHandler", "Herwig.so");
void QEDRadiationHandler::Init() {
static ClassDocumentation<QEDRadiationHandler> documentation
("The QEDRadiationHandler class is designed to be used as a PostSubProcessHandler"
"so that the same approach as for radiation in decays can be used for resonances"
"produced as part of the hard process");
static Reference<QEDRadiationHandler,DecayRadiationGenerator>
interfaceRadiationGenerator
("RadiationGenerator",
"Reference to the DecayRadiationGenerator",
&QEDRadiationHandler::_generator, false, false, true, false, false);
static ParVector<QEDRadiationHandler,long> interfaceDecayingParticles
("DecayingParticles",
"List of PDF codes of the particles which should have radiation"
" generated for them.",
&QEDRadiationHandler::_decayingParticles, -1, long(24), 0, 0,
false, false, Interface::nolimits);
static ParVector<QEDRadiationHandler,long> interfaceDecayProducts
("DecayProducts",
"List of PDG codes of the particles which should be present"
" as decay products for the radiation to be generated.",
&QEDRadiationHandler::_decayProducts, -1, long(11), 0, 0,
false, false, Interface::nolimits);
}
diff --git a/Models/General/ModelGenerator.cc b/Models/General/ModelGenerator.cc
--- a/Models/General/ModelGenerator.cc
+++ b/Models/General/ModelGenerator.cc
@@ -1,556 +1,556 @@
// -*- C++ -*-
//
// ModelGenerator.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the ModelGenerator class.
//
#include "ModelGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/DecayMode.h"
#include "ThePEG/Repository/CurrentGenerator.h"
#include "BSMWidthGenerator.h"
#include "Herwig/PDT/GenericMassGenerator.h"
#include "Herwig/Decay/DecayIntegrator.h"
#include "ThePEG/Repository/BaseRepository.h"
using namespace Herwig;
IBPtr ModelGenerator::clone() const {
return new_ptr(*this);
}
IBPtr ModelGenerator::fullclone() const {
return new_ptr(*this);
}
void ModelGenerator::persistentOutput(PersistentOStream & os) const {
os << hardProcessConstructors_ << _theDecayConstructor << particles_
<< offshell_ << Offsel_ << BRnorm_ << twoBodyOnly_ << howOffShell_
<< Npoints_ << Iorder_ << BWshape_ << brMin_ << decayOutput_;
}
void ModelGenerator::persistentInput(PersistentIStream & is, int) {
is >> hardProcessConstructors_ >> _theDecayConstructor >> particles_
>> offshell_ >> Offsel_ >> BRnorm_ >> twoBodyOnly_ >> howOffShell_
>> Npoints_ >> Iorder_ >> BWshape_ >> brMin_ >> decayOutput_;
}
// Static variable needed for the type description system in ThePEG.
DescribeClass<ModelGenerator,Interfaced>
describeThePEGModelGenerator("Herwig::ModelGenerator", "Herwig.so");
void ModelGenerator::Init() {
static ClassDocumentation<ModelGenerator> documentation
("This class controls the the use of BSM physics.",
"BSM physics was produced using the algorithm of "
"\\cite{Gigg:2007cr,Gigg:2008yc}",
"\\bibitem{Gigg:2007cr} M.~Gigg and P.~Richardson, \n"
"Eur.\\ Phys.\\ J.\\ C {\\bf 51} (2007) 989.\n"
"%%CITATION = EPHJA,C51,989;%%\n"
" %\\cite{Gigg:2008yc}\n"
"\\bibitem{Gigg:2008yc}\n"
" M.~A.~Gigg and P.~Richardson,\n"
" %``Simulation of Finite Width Effects in Physics Beyond the Standard Model,''\n"
" arXiv:0805.3037 [hep-ph].\n"
" %%CITATION = ARXIV:0805.3037;%%\n"
);
static RefVector<ModelGenerator,HardProcessConstructor>
interfaceHardProcessConstructors
("HardProcessConstructors",
"The objects to construct hard processes",
&ModelGenerator::hardProcessConstructors_, -1,
false, false, true, false, false);
static Reference<ModelGenerator,Herwig::DecayConstructor>
interfaceDecayConstructor
("DecayConstructor",
"Pointer to DecayConstructor helper class",
&ModelGenerator::_theDecayConstructor, false, false, true, false);
static RefVector<ModelGenerator,ThePEG::ParticleData> interfaceModelParticles
("DecayParticles",
"ParticleData pointers to the particles requiring spin correlation "
"decayers. If decay modes do not exist they will also be created.",
&ModelGenerator::particles_, -1, false, false, true, false);
static RefVector<ModelGenerator,ParticleData> interfaceOffshell
("Offshell",
"The particles to treat as off-shell",
&ModelGenerator::offshell_, -1, false, false, true, false);
static Switch<ModelGenerator,int> interfaceWhichOffshell
("WhichOffshell",
"A switch to determine which particles to create mass and width "
"generators for.",
&ModelGenerator::Offsel_, 0, false, false);
static SwitchOption interfaceWhichOffshellSelected
(interfaceWhichOffshell,
"Selected",
"Only create mass and width generators for the particles specified",
0);
static SwitchOption interfaceWhichOffshellAll
(interfaceWhichOffshell,
"All",
"Treat all particles specified in the DecayParticles "
"list as off-shell",
1);
static Switch<ModelGenerator,bool> interfaceBRNormalize
("BRNormalize",
"Whether to normalize the partial widths to BR*total width for an "
"on-shell particle",
&ModelGenerator::BRnorm_, true, false, false);
static SwitchOption interfaceBRNormalizeNormalize
(interfaceBRNormalize,
"Yes",
"Normalize the partial widths",
true);
static SwitchOption interfaceBRNormalizeNoNormalize
(interfaceBRNormalize,
"No",
"Do not normalize the partial widths",
false);
static Parameter<ModelGenerator,int> interfacePoints
("InterpolationPoints",
"Number of points to use for interpolation tables when needed",
&ModelGenerator::Npoints_, 10, 5, 1000,
false, false, true);
static Parameter<ModelGenerator,unsigned int>
interfaceInterpolationOrder
("InterpolationOrder", "The interpolation order for the tables",
&ModelGenerator::Iorder_, 1, 1, 5,
false, false, Interface::limited);
static Switch<ModelGenerator,int> interfaceBreitWignerShape
("BreitWignerShape",
"Controls the shape of the mass distribution generated",
&ModelGenerator::BWshape_, 0, false, false);
static SwitchOption interfaceBreitWignerShapeDefault
(interfaceBreitWignerShape,
"Default",
"Running width with q in numerator and denominator width factor",
0);
static SwitchOption interfaceBreitWignerShapeFixedWidth
(interfaceBreitWignerShape,
"FixedWidth",
"Use a fixed width",
1);
static SwitchOption interfaceBreitWignerShapeNoq
(interfaceBreitWignerShape,
"Noq",
"Use M rather than q in the numerator and denominator width factor",
2);
static SwitchOption interfaceBreitWignerShapeNoNumerator
(interfaceBreitWignerShape,
"NoNumerator",
"Neglect the numerator factors",
3);
static Switch<ModelGenerator,bool> interfaceTwoBodyOnly
("TwoBodyOnly",
"Whether to use only two-body or all modes in the running width calculation",
&ModelGenerator::twoBodyOnly_, false, false, false);
static SwitchOption interfaceTwoBodyOnlyYes
(interfaceTwoBodyOnly,
"Yes",
"Only use two-body modes",
true);
static SwitchOption interfaceTwoBodyOnlyNo
(interfaceTwoBodyOnly,
"No",
"Use all modes",
false);
static Parameter<ModelGenerator,double> interfaceMinimumBR
("MinimumBR",
"The minimum branching fraction to include",
&ModelGenerator::brMin_, 1e-6, 0.0, 1.0,
false, false, Interface::limited);
static Switch<ModelGenerator,unsigned int> interfaceDecayOutput
("DecayOutput",
"Option to control the output of the decay mode information",
&ModelGenerator::decayOutput_, 1, false, false);
static SwitchOption interfaceDecayOutputNone
(interfaceDecayOutput,
"None",
"No output",
0);
static SwitchOption interfaceDecayOutputPlain
(interfaceDecayOutput,
"Plain",
"Default plain text output",
1);
static SwitchOption interfaceDecayOutputSLHA
(interfaceDecayOutput,
"SLHA",
"Output in the Susy Les Houches Accord format",
2);
static Parameter<ModelGenerator,double> interfaceMinimumWidthFraction
("MinimumWidthFraction",
"Minimum fraction of the particle's mass the width can be"
" for the off-shell treatment.",
&ModelGenerator::minWidth_, 1e-6, 1e-15, 1.,
false, false, Interface::limited);
static Parameter<ModelGenerator,double> interfaceHowMuchOffShell
("HowMuchOffShell",
"The multiple of the particle's width by which it is allowed to be off-shell",
&ModelGenerator::howOffShell_, 5., 0.0, 100.,
false, false, Interface::limited);
}
namespace {
/// Helper function for sorting by mass
inline bool massIsLess(tcPDPtr a, tcPDPtr b) {
return a->mass() < b->mass();
}
// Helper function to find minimum possible mass of a particle
inline Energy minimumMass(tcPDPtr parent) {
Energy output(Constants::MaxEnergy);
for(set<tDMPtr>::const_iterator dit = parent->decayModes().begin();
dit != parent->decayModes().end(); ++dit) {
Energy outMass(ZERO);
for(unsigned int ix=0;ix<(**dit).orderedProducts().size();++ix) {
outMass += (**dit).orderedProducts()[ix]->massMin();
}
output = min(output,outMass);
}
return output;
}
}
void ModelGenerator::doinit() {
useMe();
Interfaced::doinit();
// make sure the model is initialized
Ptr<Herwig::StandardModel>::pointer model
= dynamic_ptr_cast<Ptr<Herwig::StandardModel>::pointer>(generator()->standardModel());
model->init();
// and the vertices
for(size_t iv = 0; iv < model->numberOfVertices(); ++iv)
model->vertex(iv)->init();
// uniq and sort DecayParticles list by mass
set<PDPtr> tmp(particles_.begin(),particles_.end());
particles_.assign(tmp.begin(),tmp.end());
sort(particles_.begin(),particles_.end(),massIsLess);
//create decayers and decaymodes (if necessary)
if( _theDecayConstructor ) {
_theDecayConstructor->init();
_theDecayConstructor->createDecayers(particles_,brMin_);
}
// write out decays with spin correlations
ostream & os = CurrentGenerator::current().misc();
ofstream ofs;
if ( decayOutput_ > 1 ) {
string filename
= CurrentGenerator::current().filename() + "-BR.spc";
ofs.open(filename.c_str());
}
if(decayOutput_!=0) {
if(decayOutput_==1) {
os << "# The decay modes listed below will have spin\n"
<< "# correlations included when they are generated.\n#\n#";
}
else {
ofs << "# Herwig decay tables in SUSY Les Houches accord format\n";
ofs << "Block DCINFO # Program information\n";
ofs << "1 Herwig # Decay Calculator\n";
ofs << "2 " << generator()->strategy()->versionstring()
<< " # Version number\n";
}
}
//create mass and width generators for the requested particles
set<PDPtr> offShell;
if( Offsel_ == 0 ) offShell = set<PDPtr>(offshell_.begin() ,offshell_.end() );
else offShell = set<PDPtr>(particles_.begin(),particles_.end());
for(PDVector::iterator pit = particles_.begin();
pit != particles_.end(); ++pit) {
tPDPtr parent = *pit;
// Check decays for ones where quarks cannot be put on constituent
// mass-shell
checkDecays(parent);
parent->reset();
// Now switch off the modes if needed
for(DecaySet::const_iterator it=parent->decayModes().begin();
it!=parent->decayModes().end();++it) {
if( _theDecayConstructor->disableDecayMode((**it).tag()) ) {
DMPtr mode = *it;
generator()->preinitInterface(mode, "Active", "set", "No");
if(mode->CC()) {
DMPtr CCmode = mode->CC();
generator()->preinitInterface(CCmode, "Active", "set", "No");
}
}
}
parent->update();
if( parent->CC() ) parent->CC()->synchronize();
if( parent->decaySelector().empty() ) {
parent->stable(true);
parent->width(ZERO);
parent->widthCut(ZERO);
parent->massGenerator(tGenericMassGeneratorPtr());
parent->widthGenerator(tGenericWidthGeneratorPtr());
}
else {
if(parent->mass()*minWidth_>parent->width()) {
parent->massGenerator(tGenericMassGeneratorPtr());
parent->widthGenerator(tGenericWidthGeneratorPtr());
}
else {
if( offShell.find(*pit) != offShell.end() ) {
createWidthGenerator(*pit);
}
else {
parent->massGenerator(tGenericMassGeneratorPtr());
parent->widthGenerator(tGenericWidthGeneratorPtr());
}
}
}
// set the offshellness
Energy minMass = minimumMass(parent);
Energy offShellNess = howOffShell_*parent->width();
if(minMass>parent->mass()-offShellNess) {
offShellNess = parent->mass()-minMass;
}
parent->widthCut(offShellNess);
if( parent->massGenerator() ) {
parent->massGenerator()->reset();
if(decayOutput_==1)
os << "# " <<parent->PDGName() << " will be considered off-shell.\n#\n";
}
if( parent->widthGenerator() ) parent->widthGenerator()->reset();
}
// loop again to initialise mass and width generators
// switch off modes and write output
for(PDVector::iterator pit = particles_.begin();
pit != particles_.end(); ++pit) {
tPDPtr parent = *pit;
if(parent->widthGenerator())
parent->widthGenerator()->init();
if(parent->massGenerator())
parent->massGenerator()->init();
// output the modes if needed
if( !parent->decaySelector().empty() ) {
if ( decayOutput_ == 2 )
writeDecayModes(ofs, parent);
else
writeDecayModes(os, parent);
}
}
//Now construct hard processes given that we know which
//objects have running widths
for(unsigned int ix=0;ix<hardProcessConstructors_.size();++ix) {
hardProcessConstructors_[ix]->init();
hardProcessConstructors_[ix]->constructDiagrams();
}
}
void ModelGenerator::checkDecays(PDPtr parent) {
if( parent->stable() ) {
if(parent->coloured())
cerr << "Warning: No decays for coloured particle " << parent->PDGName() << "\n\n"
<< "have been calcluated in BSM model.\n"
<< "This may cause problems in the hadronization phase.\n"
<< "You may have forgotten to switch on the decay mode calculation using\n"
<< " set TwoBodyDC:CreateDecayModes Yes\n"
<< " set ThreeBodyDC:CreateDecayModes Yes\n"
<< " set WeakDecayConstructor:CreateDecayModes Yes\n"
<< "or the decays of this particle are missing from your\n"
<< "input spectrum and decay file in the SLHA format.\n\n";
return;
}
DecaySet::iterator dit = parent->decayModes().begin();
DecaySet::iterator dend = parent->decayModes().end();
Energy oldwidth(parent->width()), newwidth(ZERO);
bool rescalebrat(false);
double brsum(0.);
for(; dit != dend; ++dit ) {
if( !(**dit).on() ) continue;
Energy release((**dit).parent()->mass());
tPDVector::const_iterator pit = (**dit).orderedProducts().begin();
tPDVector::const_iterator pend =(**dit).orderedProducts().end();
for( ; pit != pend; ++pit ) {
release -= (**pit).constituentMass();
}
if( (**dit).brat() < brMin_ || release < ZERO ) {
if( release < ZERO )
cerr << "Warning: The shower cannot be generated using this decay "
<< (**dit).tag() << " because it is too close to threshold.\nIt "
<< "will be switched off and the branching fractions of the "
<< "remaining modes rescaled.\n";
rescalebrat = true;
generator()->preinitInterface(*dit, "Active", "set", "No");
generator()->preinitInterface(*dit, "BranchingRatio",
"set", "0.0");
DecayIntegratorPtr decayer = dynamic_ptr_cast<DecayIntegratorPtr>((**dit).decayer());
if(decayer) {
generator()->preinitInterface(decayer->fullName(), "Initialize", "set","0");
}
}
else {
brsum += (**dit).brat();
newwidth += (**dit).brat()*oldwidth;
}
}
// if no modes left set stable
if(newwidth==ZERO) {
parent->stable(true);
parent->width(ZERO);
parent->widthCut(ZERO);
parent->massGenerator(tGenericMassGeneratorPtr());
parent->widthGenerator(tGenericWidthGeneratorPtr());
}
// otherwise rescale if needed
else if( ( rescalebrat || abs(brsum - 1.) > 1e-12 ) && !parent->decayModes().empty()) {
dit = parent->decayModes().begin();
dend = parent->decayModes().end();
double factor = oldwidth/newwidth;
brsum = 0.;
for( ; dit != dend; ++dit ) {
if( !(**dit).on() ) continue;
double newbrat = ((**dit).brat())*factor;
brsum += newbrat;
ostringstream brf;
brf << setprecision(13) << newbrat;
generator()->preinitInterface(*dit, "BranchingRatio",
"set", brf.str());
}
parent->width(newwidth);
if( newwidth > ZERO ) parent->cTau(hbarc/newwidth);
}
}
namespace {
struct DecayModeOrdering {
- bool operator()(tcDMPtr m1, tcDMPtr m2) {
+ bool operator()(tcDMPtr m1, tcDMPtr m2) const {
if(m1->brat()!=m2->brat()) {
return m1->brat()>m2->brat();
}
else {
if(m1->products().size()==m2->products().size()) {
ParticleMSet::const_iterator it1=m1->products().begin();
ParticleMSet::const_iterator it2=m2->products().begin();
do {
if((**it1).id()!=(**it2).id()) {
return (**it1).id()>(**it2).id();
}
++it1;
++it2;
}
while(it1!=m1->products().end()&&
it2!=m2->products().end());
assert(false);
}
else
return m1->products().size()<m2->products().size();
}
return false;
}
};
}
void ModelGenerator::writeDecayModes(ostream & os, tcPDPtr parent) const {
if(decayOutput_==0) return;
set<tcDMPtr,DecayModeOrdering> modes(parent->decayModes().begin(),
parent->decayModes().end());
if(decayOutput_==1) {
os << " Parent: " << parent->PDGName() << " Mass (GeV): "
<< parent->mass()/GeV << " Total Width (GeV): "
<< parent->width()/GeV << endl;
os << std::left << std::setw(40) << '#'
<< std::left << std::setw(20) << "Partial Width/GeV"
<< std::left << std::setw(20) << "BR" << "Yes/No\n";
for(set<tcDMPtr,DecayModeOrdering>::iterator dit=modes.begin();
dit!=modes.end();++dit)
os << std::left << std::setw(40) << (**dit).tag()
<< std::left << std::setw(20) << (**dit).brat()*parent->width()/GeV
<< std::left << std::setw(20) << (**dit).brat()
<< ((**dit).on() ? "Yes" : "No" ) << '\n';
os << "#\n#";
}
else if(decayOutput_==2) {
os << "# \t PDG \t Width\n";
os << "DECAY\t" << parent->id() << "\t" << parent->width()/GeV << "\t # " << parent->PDGName() << "\n";
for(set<tcDMPtr,DecayModeOrdering>::iterator dit=modes.begin();
dit!=modes.end();++dit) {
os << "\t" << std::left << std::setw(10)
<< (**dit).brat() << "\t" << (**dit).orderedProducts().size()
<< "\t";
for(unsigned int ix=0;ix<(**dit).orderedProducts().size();++ix)
os << std::right << std::setw(10)
<< (**dit).orderedProducts()[ix]->id() ;
for(unsigned int ix=(**dit).orderedProducts().size();ix<4;++ix)
os << "\t";
os << "# " << (**dit).tag() << "\n";
}
}
}
void ModelGenerator::createWidthGenerator(tPDPtr p) {
string wn = p->fullName() + string("-WGen");
string mn = p->fullName() + string("-MGen");
GenericMassGeneratorPtr mgen = dynamic_ptr_cast<GenericMassGeneratorPtr>
(generator()->preinitCreate("Herwig::GenericMassGenerator", mn));
BSMWidthGeneratorPtr wgen = dynamic_ptr_cast<BSMWidthGeneratorPtr>
(generator()->preinitCreate("Herwig::BSMWidthGenerator", wn));
//set the particle interface
mgen->particle(p);
wgen->particle(p);
//set the generator interfaces in the ParticleData object
generator()->preinitInterface(p, "Mass_generator","set", mn);
generator()->preinitInterface(p, "Width_generator","set", wn);
//allow the branching fraction of this particle type to vary
p->variableRatio(true);
if( p->CC() ) p->CC()->variableRatio(true);
//initialize the generators
generator()->preinitInterface(mgen, "Initialize", "set", "Yes");
generator()->preinitInterface(wgen, "Initialize", "set", "Yes");
string norm = BRnorm_ ? "Yes" : "No";
generator()->preinitInterface(wgen, "BRNormalize", "set", norm);
string twob = twoBodyOnly_ ? "Yes" : "No";
generator()->preinitInterface(wgen, "TwoBodyOnly", "set", twob);
ostringstream os;
os << Npoints_;
generator()->preinitInterface(wgen, "Points", "set", os.str());
os.str("");
os << Iorder_;
generator()->preinitInterface(wgen, "InterpolationOrder", "set",
os.str());
os.str("");
os << BWshape_;
generator()->preinitInterface(mgen, "BreitWignerShape", "set",
os.str());
}
diff --git a/Models/General/PrototypeVertex.h b/Models/General/PrototypeVertex.h
--- a/Models/General/PrototypeVertex.h
+++ b/Models/General/PrototypeVertex.h
@@ -1,477 +1,477 @@
// -*- C++ -*-
//
// PrototypeVertex.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_PrototypeVertex_H
#define HERWIG_PrototypeVertex_H
#include <stack>
#include "ThePEG/Helicity/Vertex/VertexBase.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Utilities/EnumIO.h"
#include "NBodyDecayConstructorBase.fh"
//
// This is the declaration of the PrototypeVertex class.
//
namespace Herwig {
using namespace ThePEG;
using Helicity::VertexBasePtr;
class PrototypeVertex;
ThePEG_DECLARE_POINTERS(Herwig::PrototypeVertex,PrototypeVertexPtr);
/** Pair of int,double */
typedef pair<unsigned int, double> CFPair;
/**
* A struct to order the particles in the same way as in the DecayModes
*/
struct ParticleOrdering {
/**
* Operator for the ordering
* @param p1 The first ParticleData object
* @param p2 The second ParticleData object
*/
- bool operator() (tcPDPtr p1, tcPDPtr p2) {
+ bool operator() (tcPDPtr p1, tcPDPtr p2) const {
return abs(p1->id()) > abs(p2->id()) ||
( abs(p1->id()) == abs(p2->id()) && p1->id() > p2->id() ) ||
( p1->id() == p2->id() && p1->fullName() > p2->fullName() );
}
};
/**
* A struct to order the particles in the same way as in the DecayMode's
*/
struct VertexOrdering {
/**
* Operator for the ordering
* @param p1 The first ParticleData object
* @param p2 The second ParticleData object
*/
bool operator() (const pair< tPDPtr, PrototypeVertexPtr > & p1,
const pair< tPDPtr, PrototypeVertexPtr > & p2) const {
return abs(p1.first->id()) > abs(p2.first->id()) ||
( abs(p1.first->id()) == abs(p2.first->id()) && p1.first->id() > p2.first->id() ) ||
( p1.first->id() == p2.first->id() && p1.first->fullName() > p2.first->fullName() );
}
};
typedef multiset<pair< tPDPtr, PrototypeVertexPtr >,VertexOrdering > OrderedVertices;
/**
* A set of ParticleData objects ordered as for the DecayMode's
*/
typedef multiset<PDPtr,ParticleOrdering> OrderedParticles;
/**
* Storage of a potenital n-body decay
*/
class PrototypeVertex : public Base {
public:
/**
* Default Constructor
*/
PrototypeVertex() : npart(0), possibleOnShell(false) {}
/**
* Constructor
*/
PrototypeVertex(tPDPtr in, OrderedVertices out,
VertexBasePtr v, int n) :
incoming(in), outgoing(out), vertex(v), npart(n),
possibleOnShell(false) {}
/**
* Incoming particle
*/
tPDPtr incoming;
/**
* Outgoing particles
*/
OrderedVertices outgoing;
/**
* The vertex for the interaction
*/
VertexBasePtr vertex;
/**
* The parent of the vertex
*/
tPrototypeVertexPtr parent;
/**
* Number of particles
*/
unsigned int npart;
/**
* Outgoing particles
*/
mutable OrderedParticles outPart;
/**
* Can have on-shell intermediates
*/
bool possibleOnShell;
/**
* Increment the number of particles
*/
void incrementN(int in) {
npart += in;
if(parent) parent->incrementN(in);
}
/**
* Mass of the incoming particle
*/
Energy incomingMass() {
return incoming->mass();
}
/**
* Total mass of all the outgoing particles
*/
Energy outgoingMass() {
Energy mass(ZERO);
for(OrderedVertices::const_iterator it = outgoing.begin();
it!=outgoing.end();++it) {
mass += it->second ?
it->second->outgoingMass() : it->first->mass();
}
return mass;
}
/**
* Total constituent mass of all the outgoing particles
*/
Energy outgoingConstituentMass() {
Energy mass(ZERO);
for(OrderedVertices::const_iterator it = outgoing.begin();
it!=outgoing.end();++it) {
mass += it->second ?
it->second->outgoingConstituentMass() : it->first->constituentMass();
}
return mass;
}
/**
* Check the external particles
*/
bool checkExternal(bool first=true) {
if(outPart.empty()) setOutgoing();
if(first&&outPart.find(incoming)!=outPart.end()) return false;
bool output = true;
for(OrderedVertices::const_iterator it = outgoing.begin();
it!=outgoing.end();++it) {
if(it->second&& !it->second->checkExternal(false)) output = false;
}
return output;
}
/**
* Set the outgoing particles
*/
void setOutgoing() const {
assert(outPart.empty());
for(OrderedVertices::const_iterator it = outgoing.begin();
it!=outgoing.end();++it) {
if(it->second) {
it->second->setOutgoing();
outPart.insert(it->second->outPart.begin(),
it->second->outPart.end());
}
else
outPart.insert(it->first);
}
}
/**
* Are there potential on-shell intermediates?
*/
bool canBeOnShell(unsigned int opt,Energy maxMass,bool first);
/**
* Check if same external particles
*/
bool sameDecay(const PrototypeVertex & x) const;
/**
* Create a \f$1\to2\f$ prototype
*/
static void createPrototypes(tPDPtr inpart, VertexBasePtr vertex,
std::stack<PrototypeVertexPtr> & prototypes,
NBodyDecayConstructorBasePtr decayCon);
/**
* Expand the prototypes by adding more legs
*/
static void expandPrototypes(PrototypeVertexPtr proto, VertexBasePtr vertex,
std::stack<PrototypeVertexPtr> & prototypes,
const set<PDPtr> & excluded,
NBodyDecayConstructorBasePtr decayCon);
/**
* Copy the whole structure with a new branching
*/
static PrototypeVertexPtr replicateTree(PrototypeVertexPtr parent,
PrototypeVertexPtr oldChild,
PrototypeVertexPtr & newChild);
};
/**
* Output to a stream
*/
inline ostream & operator<<(ostream & os, const PrototypeVertex & diag) {
os << diag.incoming->PDGName() << " -> ";
bool seq=false;
for(OrderedVertices::const_iterator it = diag.outgoing.begin();
it!=diag.outgoing.end();++it) {
os << it->first->PDGName() << " ";
if(it->second) seq = true;
}
os << " decays via "
<< diag.vertex->fullName() << " in a "
<< diag.npart << "-body decay\n";
if(!seq) return os;
os << "Followed by\n";
for(OrderedVertices::const_iterator it = diag.outgoing.begin();
it!=diag.outgoing.end();++it) {
if(it->second) os << *it->second;
}
return os;
}
/**
* Test whether two diagrams are identical.
*/
inline bool operator==(const PrototypeVertex & x, const PrototypeVertex & y) {
if(x.incoming != y.incoming) return false;
if(x.vertex != y.vertex) return false;
if(x.npart != y.npart) return false;
if(x.outgoing.empty()&&y.outgoing.empty()) return true;
if(x.outgoing.size() != y.outgoing.size()) return false;
OrderedVertices::const_iterator xt = x.outgoing.begin();
OrderedVertices::const_iterator yt = y.outgoing.begin();
while(xt!=x.outgoing.end()) {
if(xt->first != yt->first) return false;
// special for identical particles
OrderedVertices::const_iterator lxt = x.outgoing.lower_bound(*xt);
OrderedVertices::const_iterator uxt = x.outgoing.upper_bound(*xt);
--uxt;
// just one particle
if(lxt==uxt) {
if(xt->second && yt->second) {
if(*(xt->second)==*(yt->second)) {
++xt;
++yt;
continue;
}
else return false;
}
else if(xt->second || yt->second)
return false;
++xt;
++yt;
}
// identical particles
else {
++uxt;
OrderedVertices::const_iterator lyt = y.outgoing.lower_bound(*xt);
OrderedVertices::const_iterator uyt = y.outgoing.upper_bound(*xt);
unsigned int nx=0;
for(OrderedVertices::const_iterator ixt=lxt;ixt!=uxt;++ixt) {++nx;}
unsigned int ny=0;
for(OrderedVertices::const_iterator iyt=lyt;iyt!=uyt;++iyt) {++ny;}
if(nx!=ny) return false;
vector<bool> matched(ny,false);
for(OrderedVertices::const_iterator ixt=lxt;ixt!=uxt;++ixt) {
bool found = false;
unsigned int iy=0;
for(OrderedVertices::const_iterator iyt=lyt;iyt!=uyt;++iyt) {
if(matched[iy]) {
++iy;
continue;
}
if( (!ixt->second &&!iyt->second) ||
( ixt->second&&iyt->second &&
*(ixt->second)==*(iyt->second)) ) {
matched[iy] = true;
found = true;
break;
}
++iy;
}
if(!found) return false;
}
xt=uxt;
yt=uyt;
}
}
return true;
}
/**
* A simple vertex for the N-body diagram
*/
struct NBVertex {
/**
* Constructor taking a prototype vertex as the arguments
*/
NBVertex(PrototypeVertexPtr proto = PrototypeVertexPtr() );
/**
* Incoming particle
*/
tPDPtr incoming;
/**
* Outgoing particles
*/
mutable OrderedParticles outgoing;
/**
* The vertices
*/
list<pair<PDPtr,NBVertex> > vertices;
/**
* The vertex
*/
VertexBasePtr vertex;
};
/**
* The NBDiagram struct contains information about a \f$1\to n\f$ decay
* that has been automatically generated.
*/
struct NBDiagram : public NBVertex {
/**
* Constructor taking a prototype vertex as the arguments*/
NBDiagram(PrototypeVertexPtr proto=PrototypeVertexPtr());
/**
* The type of channel
*/
vector<unsigned int> channelType;
/** Store colour flow at \f$N_c=3\f$ information */
mutable vector<CFPair> colourFlow;
/** Store colour flow at \f$N_c=\infty\f$ information */
mutable vector<CFPair> largeNcColourFlow;
};
/**
* Output operator to allow the structure to be persistently written
* @param os The output stream
* @param x The NBVertex
*/
inline PersistentOStream & operator<<(PersistentOStream & os,
const NBVertex & x) {
os << x.incoming << x.outgoing << x.vertices << x.vertex;
return os;
}
/**
* Input operator to allow persistently written data to be read in
* @param is The input stream
* @param x The NBVertex
*/
inline PersistentIStream & operator>>(PersistentIStream & is,
NBVertex & x) {
is >> x.incoming >> x.outgoing >> x.vertices >> x.vertex;
return is;
}
/**
* Output operator to allow the structure to be persistently written
* @param os The output stream
* @param x The NBDiagram
*/
inline PersistentOStream & operator<<(PersistentOStream & os,
const NBDiagram & x) {
os << x.incoming << x.channelType << x.outgoing << x.vertices << x.vertex
<< x.colourFlow << x.largeNcColourFlow;
return os;
}
/**
* Input operator to allow persistently written data to be read in
* @param is The input stream
* @param x The NBDiagram
*/
inline PersistentIStream & operator>>(PersistentIStream & is,
NBDiagram & x) {
is >> x.incoming >> x.channelType >> x.outgoing >> x.vertices >> x.vertex
>> x.colourFlow >> x.largeNcColourFlow;
return is;
}
/**
* Output a NBVertex to a stream
*/
inline ostream & operator<<(ostream & os, const NBVertex & vertex) {
os << vertex.incoming->PDGName() << " -> ";
bool seq=false;
for(list<pair<PDPtr,NBVertex> >::const_iterator it=vertex.vertices.begin();
it!=vertex.vertices.end();++it) {
os << it->first->PDGName() << " ";
if(it->second.incoming) seq = true;
}
os << "via vertex " << vertex.vertex->fullName() << "\n";
if(!seq) return os;
os << "Followed by\n";
for(list<pair<PDPtr,NBVertex> >::const_iterator it=vertex.vertices.begin();
it!=vertex.vertices.end();++it) {
if(it->second.incoming) os << it->second;
}
return os;
}
/**
* Output a NBDiagram to a stream
*/
inline ostream & operator<<(ostream & os, const NBDiagram & diag) {
os << diag.incoming->PDGName() << " -> ";
for(OrderedParticles::const_iterator it=diag.outgoing.begin();
it!=diag.outgoing.end();++it) {
os << (**it).PDGName() << " ";
}
os << " has order ";
for(unsigned int ix=0;ix<diag.channelType.size();++ix)
os << diag.channelType[ix] << " ";
os << "\n";
os << "First decay " << diag.incoming->PDGName() << " -> ";
bool seq=false;
for(list<pair<PDPtr,NBVertex> >::const_iterator it=diag.vertices.begin();
it!=diag.vertices.end();++it) {
os << it->first->PDGName() << " ";
if(it->second.incoming) seq = true;
}
os << "via vertex " << diag.vertex->fullName() << "\n";
if(!seq) return os;
os << "Followed by\n";
for(list<pair<PDPtr,NBVertex> >::const_iterator it=diag.vertices.begin();
it!=diag.vertices.end();++it) {
if(it->second.incoming) os << it->second;
}
return os;
}
}
#endif /* HERWIG_PrototypeVertex_H */
diff --git a/Shower/ShowerHandler.cc b/Shower/ShowerHandler.cc
--- a/Shower/ShowerHandler.cc
+++ b/Shower/ShowerHandler.cc
@@ -1,1130 +1,1130 @@
// -*- C++ -*-
//
// ShowerHandler.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the ShowerHandler class.
//
#include "ShowerHandler.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Command.h"
#include "ThePEG/PDF/PartonExtractor.h"
#include "ThePEG/PDF/PartonBinInstance.h"
#include "Herwig/PDT/StandardMatchers.h"
#include "ThePEG/Cuts/Cuts.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include "ThePEG/Utilities/Throw.h"
#include "ThePEG/Utilities/StringUtils.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "Herwig/Utilities/EnumParticles.h"
#include "Herwig/PDF/MPIPDF.h"
#include "Herwig/PDF/MinBiasPDF.h"
#include "ThePEG/Handlers/EventHandler.h"
#include "Herwig/PDF/HwRemDecayer.h"
#include <cassert>
#include "ThePEG/Utilities/DescribeClass.h"
#include "Herwig/Decay/DecayIntegrator.h"
#include "Herwig/Decay/DecayPhaseSpaceMode.h"
using namespace Herwig;
DescribeClass<ShowerHandler,CascadeHandler>
describeShowerHandler ("Herwig::ShowerHandler","HwShower.so");
ShowerHandler::~ShowerHandler() {}
tShowerHandlerPtr ShowerHandler::currentHandler_ = tShowerHandlerPtr();
void ShowerHandler::doinit() {
CascadeHandler::doinit();
// copy particles to decay before showering from input vector to the
// set used in the simulation
if ( particlesDecayInShower_.empty() ) {
for(unsigned int ix=0;ix<inputparticlesDecayInShower_.size();++ix)
particlesDecayInShower_.insert(abs(inputparticlesDecayInShower_[ix]));
}
if ( profileScales() ) {
if ( profileScales()->unrestrictedPhasespace() &&
restrictPhasespace() ) {
generator()->log()
<< "ShowerApproximation warning: The scale profile chosen requires an unrestricted phase space,\n"
<< "however, the phase space was set to be restricted. Will switch to unrestricted phase space.\n"
<< flush;
restrictPhasespace_ = false;
}
}
}
IBPtr ShowerHandler::clone() const {
return new_ptr(*this);
}
IBPtr ShowerHandler::fullclone() const {
return new_ptr(*this);
}
ShowerHandler::ShowerHandler() :
maxtry_(10),maxtryMPI_(10),maxtryDP_(10),maxtryDecay_(100),
factorizationScaleFactor_(1.0),
renormalizationScaleFactor_(1.0),
hardScaleFactor_(1.0),
restrictPhasespace_(true), maxPtIsMuF_(false),
spinOpt_(1), pdfFreezingScale_(2.5*GeV),
doFSR_(true), doISR_(true),
splitHardProcess_(true),
includeSpaceTime_(false), vMin_(0.1*GeV2),
reweight_(1.0) {
inputparticlesDecayInShower_.push_back( 6 ); // top
inputparticlesDecayInShower_.push_back( 23 ); // Z0
inputparticlesDecayInShower_.push_back( 24 ); // W+/-
inputparticlesDecayInShower_.push_back( 25 ); // h0
}
void ShowerHandler::doinitrun(){
CascadeHandler::doinitrun();
//can't use isMPIOn here, because the EventHandler is not set at that stage
if(MPIHandler_) {
MPIHandler_->initialize();
if(MPIHandler_->softInt())
remDec_->initSoftInteractions(MPIHandler_->Ptmin(), MPIHandler_->beta());
}
}
void ShowerHandler::dofinish() {
CascadeHandler::dofinish();
if(MPIHandler_) MPIHandler_->finalize();
}
void ShowerHandler::persistentOutput(PersistentOStream & os) const {
os << remDec_ << ounit(pdfFreezingScale_,GeV) << maxtry_
<< maxtryMPI_ << maxtryDP_ << maxtryDecay_
<< inputparticlesDecayInShower_
<< particlesDecayInShower_ << MPIHandler_ << PDFA_ << PDFB_
<< PDFARemnant_ << PDFBRemnant_
<< includeSpaceTime_ << ounit(vMin_,GeV2)
<< factorizationScaleFactor_ << renormalizationScaleFactor_
<< hardScaleFactor_
<< restrictPhasespace_ << maxPtIsMuF_ << hardScaleProfile_
<< showerVariations_ << doFSR_ << doISR_ << splitHardProcess_
<< spinOpt_ << useConstituentMasses_;
}
void ShowerHandler::persistentInput(PersistentIStream & is, int) {
is >> remDec_ >> iunit(pdfFreezingScale_,GeV) >> maxtry_
>> maxtryMPI_ >> maxtryDP_ >> maxtryDecay_
>> inputparticlesDecayInShower_
>> particlesDecayInShower_ >> MPIHandler_ >> PDFA_ >> PDFB_
>> PDFARemnant_ >> PDFBRemnant_
>> includeSpaceTime_ >> iunit(vMin_,GeV2)
>> factorizationScaleFactor_ >> renormalizationScaleFactor_
>> hardScaleFactor_
>> restrictPhasespace_ >> maxPtIsMuF_ >> hardScaleProfile_
>> showerVariations_ >> doFSR_ >> doISR_ >> splitHardProcess_
>> spinOpt_ >> useConstituentMasses_;
}
void ShowerHandler::Init() {
static ClassDocumentation<ShowerHandler> documentation
("Main driver class for the showering.");
static Reference<ShowerHandler,HwRemDecayer>
interfaceRemDecayer("RemDecayer",
"A reference to the Remnant Decayer object",
&Herwig::ShowerHandler::remDec_,
false, false, true, false);
static Parameter<ShowerHandler,Energy> interfacePDFFreezingScale
("PDFFreezingScale",
"The PDF freezing scale",
&ShowerHandler::pdfFreezingScale_, GeV, 2.5*GeV, 2.0*GeV, 10.0*GeV,
false, false, Interface::limited);
static Parameter<ShowerHandler,unsigned int> interfaceMaxTry
("MaxTry",
"The maximum number of attempts for the main showering loop",
&ShowerHandler::maxtry_, 10, 1, 100,
false, false, Interface::limited);
static Parameter<ShowerHandler,unsigned int> interfaceMaxTryMPI
("MaxTryMPI",
"The maximum number of regeneration attempts for an additional scattering",
&ShowerHandler::maxtryMPI_, 10, 0, 100,
false, false, Interface::limited);
static Parameter<ShowerHandler,unsigned int> interfaceMaxTryDP
("MaxTryDP",
"The maximum number of regeneration attempts for an additional hard scattering",
&ShowerHandler::maxtryDP_, 10, 0, 100,
false, false, Interface::limited);
static ParVector<ShowerHandler,long> interfaceDecayInShower
("DecayInShower",
"PDG codes of the particles to be decayed in the shower",
&ShowerHandler::inputparticlesDecayInShower_, -1, 0l, -10000000l, 10000000l,
false, false, Interface::limited);
static Reference<ShowerHandler,UEBase> interfaceMPIHandler
("MPIHandler",
"The object that administers all additional scatterings.",
&ShowerHandler::MPIHandler_, false, false, true, true);
static Reference<ShowerHandler,PDFBase> interfacePDFA
("PDFA",
"The PDF for beam particle A. Overrides the particle's own PDF setting."
"By default used for both the shower and forced splitting in the remnant",
&ShowerHandler::PDFA_, false, false, true, true, false);
static Reference<ShowerHandler,PDFBase> interfacePDFB
("PDFB",
"The PDF for beam particle B. Overrides the particle's own PDF setting."
"By default used for both the shower and forced splitting in the remnant",
&ShowerHandler::PDFB_, false, false, true, true, false);
static Reference<ShowerHandler,PDFBase> interfacePDFARemnant
("PDFARemnant",
"The PDF for beam particle A used to generate forced splittings of the remnant."
" This overrides both the particle's own PDF setting and the value set by PDFA if used.",
&ShowerHandler::PDFARemnant_, false, false, true, true, false);
static Reference<ShowerHandler,PDFBase> interfacePDFBRemnant
("PDFBRemnant",
"The PDF for beam particle B used to generate forced splittings of the remnant."
" This overrides both the particle's own PDF setting and the value set by PDFB if used.",
&ShowerHandler::PDFBRemnant_, false, false, true, true, false);
static Switch<ShowerHandler,bool> interfaceIncludeSpaceTime
("IncludeSpaceTime",
"Whether to include the model for the calculation of space-time distances",
&ShowerHandler::includeSpaceTime_, false, false, false);
static SwitchOption interfaceIncludeSpaceTimeYes
(interfaceIncludeSpaceTime,
"Yes",
"Include the model",
true);
static SwitchOption interfaceIncludeSpaceTimeNo
(interfaceIncludeSpaceTime,
"No",
"Only include the displacement from the particle-s lifetime for decaying particles",
false);
static Parameter<ShowerHandler,Energy2> interfaceMinimumVirtuality
("MinimumVirtuality",
"The minimum virtuality for the space-time model",
&ShowerHandler::vMin_, GeV2, 0.1*GeV2, 0.0*GeV2, 1000.0*GeV2,
false, false, Interface::limited);
static Parameter<ShowerHandler,double> interfaceFactorizationScaleFactor
("FactorizationScaleFactor",
"The factorization scale factor.",
&ShowerHandler::factorizationScaleFactor_, 1.0, 0.0, 0,
false, false, Interface::lowerlim);
static Parameter<ShowerHandler,double> interfaceRenormalizationScaleFactor
("RenormalizationScaleFactor",
"The renormalization scale factor.",
&ShowerHandler::renormalizationScaleFactor_, 1.0, 0.0, 0,
false, false, Interface::lowerlim);
static Parameter<ShowerHandler,double> interfaceHardScaleFactor
("HardScaleFactor",
"The hard scale factor.",
&ShowerHandler::hardScaleFactor_, 1.0, 0.0, 0,
false, false, Interface::lowerlim);
static Parameter<ShowerHandler,unsigned int> interfaceMaxTryDecay
("MaxTryDecay",
"The maximum number of attempts to generate a decay",
&ShowerHandler::maxtryDecay_, 200, 10, 0,
false, false, Interface::lowerlim);
static Reference<ShowerHandler,HardScaleProfile> interfaceHardScaleProfile
("HardScaleProfile",
"The hard scale profile to use.",
&ShowerHandler::hardScaleProfile_, false, false, true, true, false);
static Switch<ShowerHandler,bool> interfaceMaxPtIsMuF
("MaxPtIsMuF",
"",
&ShowerHandler::maxPtIsMuF_, false, false, false);
static SwitchOption interfaceMaxPtIsMuFYes
(interfaceMaxPtIsMuF,
"Yes",
"",
true);
static SwitchOption interfaceMaxPtIsMuFNo
(interfaceMaxPtIsMuF,
"No",
"",
false);
static Switch<ShowerHandler,bool> interfaceRestrictPhasespace
("RestrictPhasespace",
"Switch on or off phasespace restrictions",
&ShowerHandler::restrictPhasespace_, true, false, false);
static SwitchOption interfaceRestrictPhasespaceYes
(interfaceRestrictPhasespace,
"Yes",
"Perform phasespace restrictions",
true);
static SwitchOption interfaceRestrictPhasespaceNo
(interfaceRestrictPhasespace,
"No",
"Do not perform phasespace restrictions",
false);
static Command<ShowerHandler> interfaceAddVariation
("AddVariation",
"Add a shower variation.",
&ShowerHandler::doAddVariation, false);
static Switch<ShowerHandler,bool> interfaceDoFSR
("DoFSR",
"Switch on or off final state radiation.",
&ShowerHandler::doFSR_, true, false, false);
static SwitchOption interfaceDoFSRYes
(interfaceDoFSR,
"Yes",
"Switch on final state radiation.",
true);
static SwitchOption interfaceDoFSRNo
(interfaceDoFSR,
"No",
"Switch off final state radiation.",
false);
static Switch<ShowerHandler,bool> interfaceDoISR
("DoISR",
"Switch on or off initial state radiation.",
&ShowerHandler::doISR_, true, false, false);
static SwitchOption interfaceDoISRYes
(interfaceDoISR,
"Yes",
"Switch on initial state radiation.",
true);
static SwitchOption interfaceDoISRNo
(interfaceDoISR,
"No",
"Switch off initial state radiation.",
false);
static Switch<ShowerHandler,bool> interfaceSplitHardProcess
("SplitHardProcess",
"Whether or not to try and split the hard process into production and decay processes",
&ShowerHandler::splitHardProcess_, true, false, false);
static SwitchOption interfaceSplitHardProcessYes
(interfaceSplitHardProcess,
"Yes",
"Split the hard process",
true);
static SwitchOption interfaceSplitHardProcessNo
(interfaceSplitHardProcess,
"No",
"Don't split the hard process",
false);
static Switch<ShowerHandler,unsigned int> interfaceSpinCorrelations
("SpinCorrelations",
"Treatment of spin correlations in the parton shower",
&ShowerHandler::spinOpt_, 1, false, false);
static SwitchOption interfaceSpinCorrelationsNo
(interfaceSpinCorrelations,
"No",
"No spin correlations",
0);
static SwitchOption interfaceSpinCorrelationsSpin
(interfaceSpinCorrelations,
"Yes",
"Include the azimuthal spin correlations",
1);
static Switch<ShowerHandler,bool> interfaceUseConstituentMasses
("UseConstituentMasses",
"Whether or not to use constituent masses for the reconstruction of the particle after showering.",
&ShowerHandler::useConstituentMasses_, true, false, false);
static SwitchOption interfaceUseConstituentMassesYes
(interfaceUseConstituentMasses,
"Yes",
"Use constituent masses.",
true);
static SwitchOption interfaceUseConstituentMassesNo
(interfaceUseConstituentMasses,
"No",
"Don't use constituent masses.",
false);
}
Energy ShowerHandler::hardScale() const {
assert(false);
}
void ShowerHandler::cascade() {
useMe();
// Initialise the weights in the event object
// so that any variations are output regardless of
// whether showering occurs for the given event
initializeWeights();
// get the PDF's from ThePEG (if locally overridden use the local versions)
tcPDFPtr first = PDFA_ ? tcPDFPtr(PDFA_) : firstPDF().pdf();
tcPDFPtr second = PDFB_ ? tcPDFPtr(PDFB_) : secondPDF().pdf();
resetPDFs(make_pair(first,second));
// set the PDFs for the remnant
if( ! rempdfs_.first)
rempdfs_.first = PDFARemnant_ ? PDFPtr(PDFARemnant_) : const_ptr_cast<PDFPtr>(first);
if( ! rempdfs_.second)
rempdfs_.second = PDFBRemnant_ ? PDFPtr(PDFBRemnant_) : const_ptr_cast<PDFPtr>(second);
// get the incoming partons
tPPair incomingPartons =
eventHandler()->currentCollision()->primarySubProcess()->incoming();
// and the parton bins
PBIPair incomingBins =
make_pair(lastExtractor()->partonBinInstance(incomingPartons.first),
lastExtractor()->partonBinInstance(incomingPartons.second));
// and the incoming hadrons
tPPair incomingHadrons =
eventHandler()->currentCollision()->incoming();
remnantDecayer()->setHadronContent(incomingHadrons);
// check if incoming hadron == incoming parton
// and get the incoming hadron if exists or parton otherwise
incoming_ = make_pair(incomingBins.first ?
incomingBins.first ->particle() : incomingPartons.first,
incomingBins.second ?
incomingBins.second->particle() : incomingPartons.second);
// check the collision is of the beam particles
// and if not boost collision to the right frame
// i.e. the hadron-hadron CMF of the collision
bool btotal(false);
LorentzRotation rtotal;
if(incoming_.first != incomingHadrons.first ||
incoming_.second != incomingHadrons.second ) {
btotal = true;
boostCollision(false);
}
// set the current ShowerHandler
setCurrentHandler();
// first shower the hard process
try {
SubProPtr sub = eventHandler()->currentCollision()->primarySubProcess();
incomingPartons = cascade(sub,lastXCombPtr());
}
catch(ShowerTriesVeto &veto){
throw Exception() << "Failed to generate the shower after "
<< veto.tries
<< " attempts in ShowerHandler::cascade()"
<< Exception::eventerror;
}
if(showerHardProcessVeto()) throw Veto();
// if a non-hadron collision return (both incoming non-hadronic)
if( ( !incomingBins.first||
!isResolvedHadron(incomingBins.first ->particle()))&&
( !incomingBins.second||
!isResolvedHadron(incomingBins.second->particle()))) {
// boost back to lab if needed
if(btotal) boostCollision(true);
// perform the reweighting for the hard process shower
combineWeights();
// unset the current ShowerHandler
unSetCurrentHandler();
return;
}
// get the remnants for hadronic collision
pair<tRemPPtr,tRemPPtr> remnants(getRemnants(incomingBins));
// set the starting scale of the forced splitting to the PDF freezing scale
remnantDecayer()->initialize(remnants, incoming_, *currentStep(), pdfFreezingScale());
// do the first forcedSplitting
try {
remnantDecayer()->doSplit(incomingPartons, make_pair(rempdfs_.first,rempdfs_.second), true);
}
catch (ExtraScatterVeto) {
throw Exception() << "Remnant extraction failed in "
<< "ShowerHandler::cascade() from primary interaction"
<< Exception::eventerror;
}
// perform the reweighting for the hard process shower
combineWeights();
// if no MPI return
if( !isMPIOn() ) {
remnantDecayer()->finalize();
// boost back to lab if needed
if(btotal) boostCollision(true);
// unset the current ShowerHandler
unSetCurrentHandler();
return;
}
// generate the multiple scatters use modified pdf's now:
setMPIPDFs();
// additional "hard" processes
unsigned int tries(0);
// This is the loop over additional hard scatters (most of the time
// only one, but who knows...)
for(unsigned int i=1; i <= getMPIHandler()->additionalHardProcs(); i++){
//counter for regeneration
unsigned int multSecond = 0;
// generate the additional scatters
while( multSecond < getMPIHandler()->multiplicity(i) ) {
// generate the hard scatter
tStdXCombPtr lastXC = getMPIHandler()->generate(i);
SubProPtr sub = lastXC->construct();
// add to the Step
newStep()->addSubProcess(sub);
// increment the counters
tries++;
multSecond++;
if(tries == maxtryDP_)
throw Exception() << "Failed to establish the requested number "
<< "of additional hard processes. If this error "
<< "occurs often, your selection of additional "
<< "scatter is probably unphysical"
<< Exception::eventerror;
// Generate the shower. If not possible veto the event
try {
incomingPartons = cascade(sub,lastXC);
}
catch(ShowerTriesVeto &veto){
throw Exception() << "Failed to generate the shower of "
<< "a secondary hard process after "
<< veto.tries
<< " attempts in Evolver::showerHardProcess()"
<< Exception::eventerror;
}
try {
// do the forcedSplitting
remnantDecayer()->doSplit(incomingPartons, make_pair(remmpipdfs_.first,remmpipdfs_.second), false);
}
catch(ExtraScatterVeto){
//remove all particles associated with the subprocess
newStep()->removeParticle(incomingPartons.first);
newStep()->removeParticle(incomingPartons.second);
//remove the subprocess from the list
newStep()->removeSubProcess(sub);
//regenerate the scattering
multSecond--;
continue;
}
// connect with the remnants but don't set Remnant colour,
// because that causes problems due to the multiple colour lines.
if ( !remnants.first ->extract(incomingPartons.first , false) ||
!remnants.second->extract(incomingPartons.second, false) )
throw Exception() << "Remnant extraction failed in "
<< "ShowerHandler::cascade() for additional scatter"
<< Exception::runerror;
}
// perform the reweighting for the additional hard scatter shower
combineWeights();
}
// the underlying event processes
unsigned int ptveto(1), veto(0);
unsigned int max(getMPIHandler()->multiplicity());
for(unsigned int i=0; i<max; i++) {
// check how often this scattering has been regenerated
if(veto > maxtryMPI_) break;
//generate PSpoint
tStdXCombPtr lastXC = getMPIHandler()->generate();
SubProPtr sub = lastXC->construct();
//If Algorithm=1 additional scatters of the signal type
// with pt > ptmin have to be vetoed
//with probability 1/(m+1), where m is the number of occurances in this event
if( getMPIHandler()->Algorithm() == 1 ){
//get the pT
Energy pt = sub->outgoing().front()->momentum().perp();
if(pt > getMPIHandler()->PtForVeto() && UseRandom::rnd() < 1./(ptveto+1) ){
ptveto++;
i--;
continue;
}
}
// add to the SubProcess to the step
newStep()->addSubProcess(sub);
// Run the Shower. If not possible veto the scattering
try {
incomingPartons = cascade(sub,lastXC);
}
// discard this extra scattering, but try the next one
catch(ShowerTriesVeto) {
newStep()->removeSubProcess(sub);
//regenerate the scattering
veto++;
i--;
continue;
}
try{
//do the forcedSplitting
remnantDecayer()->doSplit(incomingPartons, make_pair(remmpipdfs_.first,remmpipdfs_.second), false);
}
catch (ExtraScatterVeto) {
//remove all particles associated with the subprocess
newStep()->removeParticle(incomingPartons.first);
newStep()->removeParticle(incomingPartons.second);
//remove the subprocess from the list
newStep()->removeSubProcess(sub);
//regenerate the scattering
veto++;
i--;
continue;
}
//connect with the remnants but don't set Remnant colour,
//because that causes problems due to the multiple colour lines.
if ( !remnants.first ->extract(incomingPartons.first , false) ||
!remnants.second->extract(incomingPartons.second, false) )
throw Exception() << "Remnant extraction failed in "
<< "ShowerHandler::cascade() for MPI hard scattering"
<< Exception::runerror;
//reset veto counter
veto = 0;
// perform the reweighting for the MPI process shower
combineWeights();
}
// finalize the remnants
remnantDecayer()->finalize(getMPIHandler()->colourDisrupt(),
getMPIHandler()->softMultiplicity());
// boost back to lab if needed
if(btotal) boostCollision(true);
// unset the current ShowerHandler
unSetCurrentHandler();
getMPIHandler()->clean();
resetPDFs(make_pair(first,second));
}
void ShowerHandler::initializeWeights() {
if ( !showerVariations().empty() ) {
tEventPtr event = eventHandler()->currentEvent();
for ( map<string,ShowerVariation>::const_iterator var =
showerVariations().begin();
var != showerVariations().end(); ++var ) {
// Check that this is behaving as intended
//map<string,double>::iterator wi = event->optionalWeights().find(var->first);
//assert(wi == event->optionalWeights().end() );
event->optionalWeights()[var->first] = 1.0;
currentWeights_[var->first] = 1.0;
}
}
reweight_ = 1.0;
}
void ShowerHandler::resetWeights() {
for ( map<string,double>::iterator w = currentWeights_.begin();
w != currentWeights_.end(); ++w ) {
w->second = 1.0;
}
reweight_ = 1.0;
}
void ShowerHandler::combineWeights() {
tEventPtr event = eventHandler()->currentEvent();
for ( map<string,double>::const_iterator w =
currentWeights_.begin(); w != currentWeights_.end(); ++w ) {
map<string,double>::iterator ew = event->optionalWeights().find(w->first);
if ( ew != event->optionalWeights().end() )
ew->second *= w->second;
else {
assert(false && "Weight name unknown.");
//event->optionalWeights()[w->first] = w->second;
}
}
if ( reweight_ != 1.0 ) {
Ptr<StandardEventHandler>::tptr eh =
dynamic_ptr_cast<Ptr<StandardEventHandler>::tptr>(eventHandler());
if ( !eh ) {
throw Exception() << "ShowerHandler::combineWeights() : Cross section reweighting "
<< "through the shower is currently only available with standard "
<< "event generators" << Exception::runerror;
}
eh->reweight(reweight_);
}
}
string ShowerHandler::doAddVariation(string in) {
if ( in.empty() )
return "expecting a name and a variation specification";
string name = StringUtils::car(in);
ShowerVariation var;
string res = var.fromInFile(StringUtils::cdr(in));
if ( res.empty() ) {
if ( !var.firstInteraction && !var.secondaryInteractions ) {
// TODO what about decay showers?
return "variation does not apply to any shower";
}
if ( var.renormalizationScaleFactor == 1.0 &&
var.factorizationScaleFactor == 1.0 ) {
return "variation does not vary anything";
}
/*
Repository::clog() << "adding a variation with tag '" << name << "' using\nxir = "
<< var.renormalizationScaleFactor
<< " xif = "
<< var.factorizationScaleFactor
<< "\napplying to:\n"
<< "first interaction = " << var.firstInteraction << " "
<< "secondary interactions = " << var.secondaryInteractions << "\n"
<< flush;
*/
showerVariations()[name] = var;
}
return res;
}
tPPair ShowerHandler::cascade(tSubProPtr, XCPtr) {
assert(false);
}
ShowerHandler::RemPair
ShowerHandler::getRemnants(PBIPair incomingBins) {
RemPair remnants;
// first beam particle
if(incomingBins.first&&!incomingBins.first->remnants().empty()) {
remnants.first =
dynamic_ptr_cast<tRemPPtr>(incomingBins.first->remnants()[0] );
if(remnants.first) {
ParticleVector children=remnants.first->children();
for(unsigned int ix=0;ix<children.size();++ix) {
if(children[ix]->dataPtr()==remnants.first->dataPtr())
remnants.first = dynamic_ptr_cast<RemPPtr>(children[ix]);
}
//remove existing colour lines from the remnants
if(remnants.first->colourLine())
remnants.first->colourLine()->removeColoured(remnants.first);
if(remnants.first->antiColourLine())
remnants.first->antiColourLine()->removeAntiColoured(remnants.first);
}
}
// seconnd beam particle
if(incomingBins.second&&!incomingBins. second->remnants().empty()) {
remnants.second =
dynamic_ptr_cast<tRemPPtr>(incomingBins.second->remnants()[0] );
if(remnants.second) {
ParticleVector children=remnants.second->children();
for(unsigned int ix=0;ix<children.size();++ix) {
if(children[ix]->dataPtr()==remnants.second->dataPtr())
remnants.second = dynamic_ptr_cast<RemPPtr>(children[ix]);
}
//remove existing colour lines from the remnants
if(remnants.second->colourLine())
remnants.second->colourLine()->removeColoured(remnants.second);
if(remnants.second->antiColourLine())
remnants.second->antiColourLine()->removeAntiColoured(remnants.second);
}
}
assert(remnants.first || remnants.second);
return remnants;
}
namespace {
void addChildren(tPPtr in,set<tPPtr> & particles) {
particles.insert(in);
for(unsigned int ix=0;ix<in->children().size();++ix)
addChildren(in->children()[ix],particles);
}
}
void ShowerHandler::boostCollision(bool boost) {
// calculate boost from lab to rest
if(!boost) {
Lorentz5Momentum ptotal=incoming_.first ->momentum()+incoming_.second->momentum();
boost_ = LorentzRotation(-ptotal.boostVector());
Axis axis((boost_*incoming_.first ->momentum()).vect().unit());
if(axis.perp2()>0.) {
double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
boost_.rotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
}
}
// first call performs the boost and second inverse
// get the particles to be boosted
set<tPPtr> particles;
addChildren(incoming_.first,particles);
addChildren(incoming_.second,particles);
// apply the boost
for(set<tPPtr>::const_iterator cit=particles.begin();
cit!=particles.end();++cit) {
(*cit)->transform(boost_);
}
if(!boost) boost_.invert();
}
void ShowerHandler::setMPIPDFs() {
if ( !mpipdfs_.first ) {
// first have to check for MinBiasPDF
tcMinBiasPDFPtr first = dynamic_ptr_cast<tcMinBiasPDFPtr>(firstPDF().pdf());
if(first)
mpipdfs_.first = new_ptr(MPIPDF(first->originalPDF()));
else
mpipdfs_.first = new_ptr(MPIPDF(firstPDF().pdf()));
}
if ( !mpipdfs_.second ) {
tcMinBiasPDFPtr second = dynamic_ptr_cast<tcMinBiasPDFPtr>(secondPDF().pdf());
if(second)
mpipdfs_.second = new_ptr(MPIPDF(second->originalPDF()));
else
mpipdfs_.second = new_ptr(MPIPDF(secondPDF().pdf()));
}
if( !remmpipdfs_.first ) {
tcMinBiasPDFPtr first = dynamic_ptr_cast<tcMinBiasPDFPtr>(rempdfs_.first);
if(first)
remmpipdfs_.first = new_ptr(MPIPDF(first->originalPDF()));
else
remmpipdfs_.first = new_ptr(MPIPDF(rempdfs_.first));
}
if( !remmpipdfs_.second ) {
tcMinBiasPDFPtr second = dynamic_ptr_cast<tcMinBiasPDFPtr>(rempdfs_.second);
if(second)
remmpipdfs_.second = new_ptr(MPIPDF(second->originalPDF()));
else
remmpipdfs_.second = new_ptr(MPIPDF(rempdfs_.second));
}
// reset the PDFs stored in the base class
resetPDFs(mpipdfs_);
}
bool ShowerHandler::isResolvedHadron(tPPtr particle) {
if(!HadronMatcher::Check(particle->data())) return false;
for(unsigned int ix=0;ix<particle->children().size();++ix) {
if(particle->children()[ix]->id()==ParticleID::Remnant) return true;
}
return false;
}
namespace {
bool decayProduct(tSubProPtr subProcess,
tPPtr particle) {
// must be time-like and not incoming
if(particle->momentum().m2()<=ZERO||
particle == subProcess->incoming().first||
particle == subProcess->incoming().second) return false;
// if only 1 outgoing and this is it
if(subProcess->outgoing().size()==1 &&
subProcess->outgoing()[0]==particle) return true;
// must not be the s-channel intermediate otherwise
if(find(subProcess->incoming().first->children().begin(),
subProcess->incoming().first->children().end(),particle)!=
subProcess->incoming().first->children().end()&&
find(subProcess->incoming().second->children().begin(),
subProcess->incoming().second->children().end(),particle)!=
subProcess->incoming().second->children().end()&&
subProcess->incoming().first ->children().size()==1&&
subProcess->incoming().second->children().size()==1)
return false;
// if non-coloured this is enough
if(!particle->dataPtr()->coloured()) return true;
// if coloured must be unstable
if(particle->dataPtr()->stable()) return false;
// must not have same particle type as a child
int id = particle->id();
for(unsigned int ix=0;ix<particle->children().size();++ix)
if(particle->children()[ix]->id()==id) return false;
// otherwise its a decaying particle
return true;
}
PPtr findParent(PPtr original, bool & isHard,
set<PPtr> outgoingset,
tSubProPtr subProcess) {
PPtr parent=original;
isHard |=(outgoingset.find(original) != outgoingset.end());
if(!original->parents().empty()) {
PPtr orig=original->parents()[0];
if(decayProduct(subProcess,orig))
parent=findParent(orig,isHard,outgoingset,subProcess);
}
return parent;
}
}
void ShowerHandler::findDecayProducts(PPtr in,PerturbativeProcessPtr hard,
DecayProcessMap & decay) const {
ParticleVector children=in->children();
for(ParticleVector::const_iterator it=children.begin(); it!=children.end();++it) {
// if decayed or should be decayed in shower make the PerturbaitveProcess
bool radiates = false;
if(!(**it).children().empty()) {
// remove d,u,s,c,b quarks and leptons other than on-shell taus
if( StandardQCDPartonMatcher::Check((**it).id()) ||
( LeptonMatcher::Check((**it).id()) && !(abs((**it).id())==ParticleID::tauminus &&
abs((**it).mass()-(**it).dataPtr()->mass())<MeV))) {
radiates = true;
}
else {
bool foundParticle(false),foundGauge(false);
for(unsigned int iy=0;iy<(**it).children().size();++iy) {
if((**it).children()[iy]->id()==(**it).id()) {
foundParticle = true;
}
else if((**it).children()[iy]->id()==ParticleID::g ||
(**it).children()[iy]->id()==ParticleID::gamma) {
foundGauge = true;
}
}
radiates = foundParticle && foundGauge;
}
}
if(radiates) {
findDecayProducts(*it,hard,decay);
}
else if(!(**it).children().empty()||
(decaysInShower((**it).id())&&!(**it).dataPtr()->stable())) {
createDecayProcess(*it,hard,decay);
}
else {
hard->outgoing().push_back(make_pair(*it,PerturbativeProcessPtr()));
}
}
}
void ShowerHandler::splitHardProcess(tPVector tagged, PerturbativeProcessPtr & hard,
DecayProcessMap & decay) const {
// temporary storage of the particles
set<PPtr> hardParticles;
// tagged particles in a set
set<PPtr> outgoingset(tagged.begin(),tagged.end());
bool isHard=false;
// loop over the tagged particles
for (tParticleVector::const_iterator taggedP = tagged.begin();
taggedP != tagged.end(); ++taggedP) {
// skip remnants
if (eventHandler()->currentCollision()&&
eventHandler()->currentCollision()->isRemnant(*taggedP)) continue;
// find the parent and whether its a decaying particle
bool isDecayProd=false;
// check if hard
isHard |=(outgoingset.find(*taggedP) != outgoingset.end());
if(splitHardProcess_) {
tPPtr parent = *taggedP;
// check if from s channel decaying colourless particle
while(parent&&!parent->parents().empty()&&!isDecayProd) {
parent = parent->parents()[0];
if(parent == subProcess_->incoming().first ||
parent == subProcess_->incoming().second ) break;
isDecayProd = decayProduct(subProcess_,parent);
}
if (isDecayProd)
hardParticles.insert(findParent(parent,isHard,outgoingset,subProcess_));
}
if (!isDecayProd)
hardParticles.insert(*taggedP);
}
// there must be something to shower
if(hardParticles.empty())
throw Exception() << "No particles to shower in "
<< "ShowerHandler::splitHardProcess()"
<< Exception::eventerror;
// must be a hard process
if(!isHard)
throw Exception() << "Starting on decay not yet implemented in "
<< "ShowerHandler::splitHardProcess()"
<< Exception::runerror;
// create the hard process
hard = new_ptr(PerturbativeProcess());
// incoming particles
hard->incoming().push_back(make_pair(subProcess_->incoming().first ,PerturbativeProcessPtr()));
hard->incoming().push_back(make_pair(subProcess_->incoming().second,PerturbativeProcessPtr()));
// outgoing particles
for(set<PPtr>::const_iterator it=hardParticles.begin();it!=hardParticles.end();++it) {
// if decayed or should be decayed in shower make the tree
PPtr orig = *it;
bool radiates = false;
if(!orig->children().empty()) {
// remove d,u,s,c,b quarks and leptons other than on-shell taus
if( StandardQCDPartonMatcher::Check(orig->id()) ||
( LeptonMatcher::Check(orig->id()) &&
!(abs(orig->id())==ParticleID::tauminus && abs(orig->mass()-orig->dataPtr()->mass())<MeV))) {
radiates = true;
}
else {
bool foundParticle(false),foundGauge(false);
for(unsigned int iy=0;iy<orig->children().size();++iy) {
if(orig->children()[iy]->id()==orig->id()) {
foundParticle = true;
}
else if(orig->children()[iy]->id()==ParticleID::g ||
orig->children()[iy]->id()==ParticleID::gamma) {
foundGauge = true;
}
}
radiates = foundParticle && foundGauge;
}
}
if(radiates) {
findDecayProducts(orig,hard,decay);
}
else if(!(**it).children().empty()||
(decaysInShower((**it).id())&&!(**it).dataPtr()->stable())) {
createDecayProcess(*it,hard,decay);
}
else {
hard->outgoing().push_back(make_pair(*it,PerturbativeProcessPtr()));
}
}
}
void ShowerHandler::createDecayProcess(PPtr in,PerturbativeProcessPtr hard, DecayProcessMap & decay) const {
// there must be an incoming particle
assert(in);
// create the new process and connect with the parent
PerturbativeProcessPtr newDecay=new_ptr(PerturbativeProcess());
newDecay->incoming().push_back(make_pair(in,hard));
Energy width=in->dataPtr()->generateWidth(in->mass());
decay.insert(make_pair(width,newDecay));
hard->outgoing().push_back(make_pair(in,newDecay));
// we need to deal with the decay products if decayed
ParticleVector children = in->children();
if(!children.empty()) {
for(ParticleVector::const_iterator it = children.begin();
it!= children.end(); ++it) {
// if decayed or should be decayed in shower make the tree
in->abandonChild(*it);
bool radiates = false;
if(!(**it).children().empty()) {
if(StandardQCDPartonMatcher::Check((**it).id())||
(LeptonMatcher::Check((**it).id())&& !(abs((**it).id())==ParticleID::tauminus &&
abs((**it).mass()-(**it).dataPtr()->mass())<MeV))) {
radiates = true;
}
else {
bool foundParticle(false),foundGauge(false);
for(unsigned int iy=0;iy<(**it).children().size();++iy) {
if((**it).children()[iy]->id()==(**it).id()) {
foundParticle = true;
}
else if((**it).children()[iy]->id()==ParticleID::g ||
(**it).children()[iy]->id()==ParticleID::gamma) {
foundGauge = true;
}
}
radiates = foundParticle && foundGauge;
}
// finally assume all non-decaying particles are in this class
// pr 27/11/15 not sure about this bit
// if(!radiates) {
// radiates = !decaysInShower((**it).id());
// }
}
if(radiates) {
findDecayProducts(*it,newDecay,decay);
}
else if(!(**it).children().empty()||
(decaysInShower((**it).id())&&!(**it).dataPtr()->stable())) {
createDecayProcess(*it,newDecay,decay);
}
else {
newDecay->outgoing().push_back(make_pair(*it,PerturbativeProcessPtr()));
}
}
}
}
tDMPtr ShowerHandler::decay(PerturbativeProcessPtr process,
DecayProcessMap & decayMap,
bool radPhotons ) const {
PPtr parent = process->incoming()[0].first;
assert(parent);
if(parent->spinInfo()) parent->spinInfo()->decay(true);
unsigned int ntry = 0;
ParticleVector children;
tDMPtr dm = DMPtr();
while (true) {
// exit if fails
if (++ntry>=maxtryDecay_)
throw Exception() << "Failed to perform decay in ShowerHandler::decay()"
<< " after " << maxtryDecay_
<< " attempts for " << parent->PDGName()
<< Exception::eventerror;
// select decay mode
dm = parent->data().selectMode(*parent);
if(!dm)
throw Exception() << "Failed to select decay mode in ShowerHandler::decay()"
<< "for " << parent->PDGName()
<< Exception::eventerror;
if(!dm->decayer())
throw Exception() << "No Decayer for selected decay mode "
<< " in ShowerHandler::decay()"
<< Exception::runerror;
// start of try block
try {
children = dm->decayer()->decay(*dm, *parent);
// if no children have another go
if(children.empty()) continue;
if(radPhotons){
// generate radiation in the decay
tDecayIntegratorPtr hwdec=dynamic_ptr_cast<tDecayIntegratorPtr>(dm->decayer());
if (hwdec && hwdec->canGeneratePhotons())
children = hwdec->generatePhotons(*parent,children);
}
// set up parent
parent->decayMode(dm);
// add children
for (unsigned int i = 0, N = children.size(); i < N; ++i ) {
children[i]->setLabVertex(parent->labDecayVertex());
//parent->addChild(children[i]);
}
// if succeeded break out of loop
break;
}
catch(Veto) {
}
}
assert(!children.empty());
for(ParticleVector::const_iterator it = children.begin();
it!= children.end(); ++it) {
if(!(**it).children().empty()||
(decaysInShower((**it).id())&&!(**it).dataPtr()->stable())) {
createDecayProcess(*it,process,decayMap);
}
else {
process->outgoing().push_back(make_pair(*it,PerturbativeProcessPtr()));
}
}
return dm;
}
// Note: The tag must be constructed from an ordered particle container.
tDMPtr ShowerHandler::findDecayMode(const string & tag) const {
static map<string,DMPtr> cache;
map<string,DMPtr>::const_iterator pos = cache.find(tag);
if ( pos != cache.end() )
return pos->second;
tDMPtr dm = CurrentGenerator::current().findDecayMode(tag);
cache[tag] = dm;
return dm;
}
/**
* Operator for the particle ordering
* @param p1 The first ParticleData object
* @param p2 The second ParticleData object
*/
-bool ShowerHandler::ParticleOrdering::operator() (tcPDPtr p1, tcPDPtr p2) {
+bool ShowerHandler::ParticleOrdering::operator() (tcPDPtr p1, tcPDPtr p2) const {
return abs(p1->id()) > abs(p2->id()) ||
( abs(p1->id()) == abs(p2->id()) && p1->id() > p2->id() ) ||
( p1->id() == p2->id() && p1->fullName() > p2->fullName() );
}
diff --git a/Shower/ShowerHandler.h b/Shower/ShowerHandler.h
--- a/Shower/ShowerHandler.h
+++ b/Shower/ShowerHandler.h
@@ -1,898 +1,898 @@
// -*- C++ -*-
//
// ShowerHandler.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_ShowerHandler_H
#define HERWIG_ShowerHandler_H
//
// This is the declaration of the ShowerHandler class.
//
#include "ThePEG/Handlers/EventHandler.h"
#include "ThePEG/Handlers/CascadeHandler.h"
#include "ShowerVariation.h"
#include "Herwig/PDF/HwRemDecayer.fh"
#include "ThePEG/EventRecord/RemnantParticle.fh"
#include "UEBase.h"
#include "PerturbativeProcess.h"
#include "Herwig/MatrixElement/Matchbox/Matching/HardScaleProfile.h"
#include "ShowerHandler.fh"
namespace Herwig {
using namespace ThePEG;
/** \ingroup Shower
*
* This class is the main driver of the shower: it is responsible for
* the proper handling of all other specific collaborating classes
* and for the storing of the produced particles in the event record.
*
* @see \ref ShowerHandlerInterfaces "The interfaces"
*
* @see ThePEG::CascadeHandler
* @see MPIHandler
* @see HwRemDecayer
*/
class ShowerHandler: public CascadeHandler {
public:
/**
* Typedef for a pair of ThePEG::RemnantParticle pointers.
*/
typedef pair<tRemPPtr, tRemPPtr> RemPair;
public:
/**
* Default constructor
*/
ShowerHandler();
/**
* Destructor
*/
virtual ~ShowerHandler();
public:
/**
* The main method which manages the multiple interactions and starts
* the shower by calling cascade(sub, lastXC).
*/
virtual void cascade();
/**
* pointer to "this", the current ShowerHandler.
*/
static const tShowerHandlerPtr currentHandler() {
assert(currentHandler_);
return currentHandler_;
}
/**
* pointer to "this", the current ShowerHandler.
*/
static bool currentHandlerIsSet() {
return currentHandler_;
}
public:
/**
* Hook to allow vetoing of event after showering hard sub-process
* as in e.g. MLM merging.
*/
virtual bool showerHardProcessVeto() const { return false; }
/**
* Return true, if this cascade handler will perform reshuffling from hard
* process masses.
*/
virtual bool isReshuffling() const { return true; }
/**
* Return true, if this cascade handler will put the final state
* particles to their constituent mass. If false the nominal mass is used.
*/
virtual bool retConstituentMasses() const { return useConstituentMasses_; }
/**
* Return true, if the shower handler can generate a truncated
* shower for POWHEG style events generated using Matchbox
*/
virtual bool canHandleMatchboxTrunc() const { return false; }
/**
* Get the PDF freezing scale
*/
Energy pdfFreezingScale() const { return pdfFreezingScale_; }
/**
* Get the local PDFs.
*/
PDFPtr getPDFA() const {return PDFA_;}
/**
* Get the local PDFs.
*/
PDFPtr getPDFB() const {return PDFB_;}
/**
* Return true if currently the primary subprocess is showered.
*/
bool firstInteraction() const {
if (!eventHandler()->currentCollision())return true;
return ( subProcess_ ==
eventHandler()->currentCollision()->primarySubProcess() );
}
/**
* Return the remnant decayer.
*/
tHwRemDecPtr remnantDecayer() const { return remDec_; }
/**
* Split the hard process into production and decays
* @param tagged The tagged particles from the StepHandler
* @param hard The hard perturbative process
* @param decay The decay particles
*/
void splitHardProcess(tPVector tagged, PerturbativeProcessPtr & hard,
DecayProcessMap & decay) const;
/**
* Information if the Showerhandler splits the hard process.
*/
bool doesSplitHardProcess()const {return splitHardProcess_;}
/**
* Decay a particle.
* radPhotons switches the generation of photon
* radiation on/off.
* Required for Dipole Shower but not QTilde Shower.
*/
tDMPtr decay(PerturbativeProcessPtr,
DecayProcessMap & decay,
bool radPhotons = false) const;
/**
* Cached lookup of decay modes.
* Generator::findDecayMode() is not efficient.
*/
tDMPtr findDecayMode(const string & tag) const;
/**
* A struct to order the particles in the same way as in the DecayMode's
*/
struct ParticleOrdering {
- bool operator() (tcPDPtr p1, tcPDPtr p2);
+ bool operator() (tcPDPtr p1, tcPDPtr p2) const;
};
/**
* A container for ordered particles required
* for constructing tags for decay mode lookup.
*/
typedef multiset<tcPDPtr,ParticleOrdering> OrderedParticles;
public:
/**
* @name Switches for initial- and final-state radiation
*/
//@{
/**
* Switch for any radiation
*/
bool doRadiation() const {return doFSR_ || doISR_;}
/**
* Switch on or off final state radiation.
*/
bool doFSR() const { return doFSR_;}
/**
* Switch on or off initial state radiation.
*/
bool doISR() const { return doISR_;}
//@}
public:
/**
* @name Switches for scales
*/
//@{
/**
* Return true if maximum pt should be deduced from the factorization scale
*/
bool hardScaleIsMuF() const { return maxPtIsMuF_; }
/**
* The factorization scale factor.
*/
double factorizationScaleFactor() const {
return factorizationScaleFactor_;
}
/**
* The renormalization scale factor.
*/
double renFac() const {
return renormalizationScaleFactor_;
}
/**
* The factorization scale factor.
*/
double facFac() const {
return factorizationScaleFactor_;
}
/**
* The renormalization scale factor.
*/
double renormalizationScaleFactor() const {
return renormalizationScaleFactor_;
}
/**
* The scale factor for the hard scale
*/
double hardScaleFactor() const {
return hardScaleFactor_;
}
/**
* Return true, if the phase space restrictions of the dipole shower should
* be applied.
*/
bool restrictPhasespace() const { return restrictPhasespace_; }
/**
* Return profile scales
*/
Ptr<HardScaleProfile>::tptr profileScales() const { return hardScaleProfile_; }
/**
* Return the relevant hard scale to be used in the profile scales
*/
virtual Energy hardScale() const;
/**
* Return information about shower phase space choices
*/
virtual int showerPhaseSpaceOption() const {
assert(false && "not implemented in general");
return -1;
}
//@}
public:
/**
* Access the shower variations
*/
map<string,ShowerVariation>& showerVariations() {
return showerVariations_;
}
/**
* Return the shower variations
*/
const map<string,ShowerVariation>& showerVariations() const {
return showerVariations_;
}
/**
* Access the current Weights
*/
map<string,double>& currentWeights() {
return currentWeights_;
}
/**
* Return the current Weights
*/
const map<string,double>& currentWeights() const {
return currentWeights_;
}
/**
* Change the current reweighting factor
*/
void reweight(double w) {
reweight_ = w;
}
/**
* Return the current reweighting factor
*/
double reweight() const {
return reweight_;
}
public :
/**
* Access to switches for spin correlations
*/
//@{
/**
* Spin Correlations
*/
unsigned int spinCorrelations() const {
return spinOpt_;
}
/**
* Any correlations
*/
virtual bool correlations() const {
return spinOpt_!=0;
}
//@}
public:
/**
* struct that is used to catch exceptions which are thrown
* due to energy conservation issues of additional scatters
*/
struct ExtraScatterVeto {};
/**
* struct that is used to catch exceptions which are thrown
* due to fact that the Shower has been invoked more than
* a defined threshold on a certain configuration
*/
struct ShowerTriesVeto {
/** variable to store the number of attempts */
const int tries;
/** constructor */
ShowerTriesVeto(int t) : tries(t) {}
};
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Functions to perform the cascade
*/
//@{
/**
* The main method which manages the showering of a subprocess.
*/
virtual tPPair cascade(tSubProPtr sub, XCPtr xcomb);
/**
* Set up for the cascade
*/
void prepareCascade(tSubProPtr sub) {
current_ = currentStep();
subProcess_ = sub;
}
/**
* Boost all the particles in the collision so that the collision always occurs
* in the rest frame with the incoming particles along the z axis
*/
void boostCollision(bool boost);
//@}
protected:
/**
* Set/unset the current shower handler
*/
//@{
/**
* Set the current handler
*/
void setCurrentHandler() {
currentHandler_ = tShowerHandlerPtr(this);
}
/**
* Unset the current handler
*/
void unSetCurrentHandler() {
currentHandler_ = tShowerHandlerPtr();
}
//@}
protected:
/**
* @name Members relating to the underlying event and MPI
*/
//@{
/**
* Return true if multiple parton interactions are switched on
* and can be used for this beam setup.
*/
bool isMPIOn() const {
return MPIHandler_ && MPIHandler_->beamOK();
}
/**
* Access function for the MPIHandler, it should only be called after
* checking with isMPIOn.
*/
tUEBasePtr getMPIHandler() const {
assert(MPIHandler_);
return MPIHandler_;
}
/**
* Is a beam particle where hadronic structure is resolved
*/
bool isResolvedHadron(tPPtr);
/**
* Get the remnants from the ThePEG::PartonBinInstance es and
* do some checks.
*/
RemPair getRemnants(PBIPair incbins);
/**
* Reset the PDF's after the hard collision has been showered
*/
void setMPIPDFs();
//@}
public:
/**
* Check if a particle decays in the shower
* @param id The PDG code for the particle
*/
bool decaysInShower(long id) const {
return ( particlesDecayInShower_.find( abs(id) ) !=
particlesDecayInShower_.end() );
}
protected:
/**
* Members to handle splitting up of hard process and decays
*/
//@{
/**
* Find decay products from the hard process and create decay processes
* @param parent The parent particle
* @param hard The hard process
* @param decay The decay processes
*/
void findDecayProducts(PPtr parent, PerturbativeProcessPtr hard, DecayProcessMap & decay) const;
/**
* Find decay products from the hard process and create decay processes
* @param parent The parent particle
* @param hard The parent hard process
* @param decay The decay processes
*/
void createDecayProcess(PPtr parent,PerturbativeProcessPtr hard, DecayProcessMap & decay) const;
//@}
/**
* @name Functions to return information relevant to the process being showered
*/
//@{
/**
* Return the currently used SubProcess.
*/
tSubProPtr currentSubProcess() const {
assert(subProcess_);
return subProcess_;
}
/**
* Access to the incoming beam particles
*/
tPPair incomingBeams() const {
return incoming_;
}
//@}
protected:
/**
* Weight handling for shower variations
*/
//@
/**
* Combine the variation weights which have been encountered
*/
void combineWeights();
/**
* Initialise the weights in currentEvent()
*/
void initializeWeights();
/**
* Reset the current weights
*/
void resetWeights();
//@}
protected:
/**
* Return the maximum number of attempts for showering
* a given subprocess.
*/
unsigned int maxtry() const { return maxtry_; }
protected:
/**
* Parameters for the space-time model
*/
//@{
/**
* Whether or not to include spa-cetime distances in the shower
*/
bool includeSpaceTime() const {return includeSpaceTime_;}
/**
* The minimum virtuality for the space-time model
*/
Energy2 vMin() const {return vMin_;}
//@}
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
/**
* Finalize this object. Called in the run phase just after a
* run has ended. Used eg. to write out statistics.
*/
virtual void dofinish();
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
ShowerHandler & operator=(const ShowerHandler &) = delete;
private:
/**
* pointer to "this", the current ShowerHandler.
*/
static tShowerHandlerPtr currentHandler_;
/**
* a MPIHandler to administer the creation of several (semihard)
* partonic interactions.
*/
UEBasePtr MPIHandler_;
/**
* Pointer to the HwRemDecayer
*/
HwRemDecPtr remDec_;
private:
/**
* Maximum tries for various stages of the showering process
*/
//@{
/**
* Maximum number of attempts for the
* main showering loop
*/
unsigned int maxtry_;
/**
* Maximum number of attempts for the regeneration of an additional
* scattering, before the number of scatters is reduced.
*/
unsigned int maxtryMPI_;
/**
* Maximum number of attempts for the regeneration of an additional
* hard scattering, before this event is vetoed.
*/
unsigned int maxtryDP_;
/**
* Maximum number of attempts to generate a decay
*/
unsigned int maxtryDecay_;
//@}
private:
/**
* Factors for the various scales
*/
//@{
/**
* The factorization scale factor.
*/
double factorizationScaleFactor_;
/**
* The renormalization scale factor.
*/
double renormalizationScaleFactor_;
/**
* The scale factor for the hard scale
*/
double hardScaleFactor_;
/**
* True, if the phase space restrictions of the dipole shower should
* be applied.
*/
bool restrictPhasespace_;
/**
* True if maximum pt should be deduced from the factorization scale
*/
bool maxPtIsMuF_;
/**
* The profile scales
*/
Ptr<HardScaleProfile>::ptr hardScaleProfile_;
//@}
/**
* Option to include spin correlations
*/
unsigned int spinOpt_;
private:
/**
* Storage of information about the current event
*/
//@{
/**
* The incoming beam particles for the current collision
*/
tPPair incoming_;
/**
* Boost to get back to the lab
*/
LorentzRotation boost_;
/**
* Const pointer to the currently handeled ThePEG::SubProcess
*/
tSubProPtr subProcess_;
/**
* Const pointer to the current step
*/
tcStepPtr current_;
//@}
private:
/**
* PDFs to be used for the various stages and related parameters
*/
//@{
/**
* The PDF freezing scale
*/
Energy pdfFreezingScale_;
/**
* PDFs to be used for the various stages and related parameters
*/
//@{
/**
* The PDF for beam particle A. Overrides the particle's own PDF setting.
*/
PDFPtr PDFA_;
/**
* The PDF for beam particle B. Overrides the particle's own PDF setting.
*/
PDFPtr PDFB_;
/**
* The PDF for beam particle A for remnant splitting. Overrides the particle's own PDF setting.
*/
PDFPtr PDFARemnant_;
/**
* The PDF for beam particle B for remnant splitting. Overrides the particle's own PDF setting.
*/
PDFPtr PDFBRemnant_;
/**
* The MPI PDF's to be used for secondary scatters.
*/
pair <PDFPtr, PDFPtr> mpipdfs_;
/**
* The MPI PDF's to be used for secondary scatters.
*/
pair <PDFPtr, PDFPtr> rempdfs_;
/**
* The MPI PDF's to be used for secondary scatters.
*/
pair <PDFPtr, PDFPtr> remmpipdfs_;
//@}
private:
/**
* @name Parameters for initial- and final-state radiation
*/
//@{
/**
* Switch on or off final state radiation.
*/
bool doFSR_;
/**
* Switch on or off initial state radiation.
*/
bool doISR_;
//@}
private:
/**
* @name Parameters for particle decays
*/
//@{
/**
* Whether or not to split into hard and decay trees
*/
bool splitHardProcess_;
/**
* PDG codes of the particles which decay during showering
* this is fast storage for use during running
*/
set<long> particlesDecayInShower_;
/**
* PDG codes of the particles which decay during showering
* this is a vector that is interfaced so they can be changed
*/
vector<long> inputparticlesDecayInShower_;
//@}
private:
/**
* Parameters for the space-time model
*/
//@{
/**
* Whether or not to include spa-cetime distances in the shower
*/
bool includeSpaceTime_;
/**
* The minimum virtuality for the space-time model
*/
Energy2 vMin_;
//@}
private:
/**
* Parameters for the constituent mass treatment.
*/
//@{
// True if shower should return constituent masses.
bool useConstituentMasses_=true;
//@}
private:
/**
* Parameters relevant for reweight and variations
*/
//@{
/**
* The shower variations
*/
map<string,ShowerVariation> showerVariations_;
/**
* Command to add a shower variation
*/
string doAddVariation(string);
/**
* A reweighting factor applied by the showering
*/
double reweight_;
/**
* The shower variation weights
*/
map<string,double> currentWeights_;
//@}
};
}
#endif /* HERWIG_ShowerHandler_H */