Page MenuHomeHEPForge
Files F8725547

No OneTemporary
Actions

View Options

diff --git a/Shower/Dipole/Base/DipoleEventRecord.cc b/Shower/Dipole/Base/DipoleEventRecord.cc
--- a/Shower/Dipole/Base/DipoleEventRecord.cc
+++ b/Shower/Dipole/Base/DipoleEventRecord.cc
@@ -1,1493 +1,1493 @@
// -*- C++ -*-
//
// DipoleEventRecord.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig is licenced under version 2 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 DipoleEventRecord class.
//
#include "DipoleEventRecord.h"
#include "Herwig/Shower/Dipole/Utility/DipolePartonSplitter.h"
#include "Herwig/Shower/ShowerHandler.h"
#include "ThePEG/PDT/DecayMode.h"
#include "Herwig/Decay/HwDecayerBase.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/PDF/PartonExtractor.h"
#include "Herwig/Shower/RealEmissionProcess.h"
#include <boost/utility.hpp>
#include <algorithm>
#include <iterator>
using namespace Herwig;
PList DipoleEventRecord::colourOrdered(PPair & in,
PList & out) {
PList colour_ordered;
size_t done_size = out.size();
if (in.first->coloured())
++done_size;
if (in.second && in.second->coloured())
++done_size;
while (colour_ordered.size() != done_size) {
PPtr current;
// start with singlets, as long as we have some
if (find(colour_ordered.begin(),colour_ordered.end(),in.first) ==
colour_ordered.end() && in.first->coloured()) {
if (!in.first->hasColour() || !in.first->hasAntiColour())
current = in.first;
}
if (!current) {
for (PList::iterator p = out.begin();
p != out.end(); ++p) {
if (find(colour_ordered.begin(),colour_ordered.end(),*p) ==
colour_ordered.end() && (**p).coloured()) {
if (!(**p).hasColour() || !(**p).hasAntiColour()) {
current = *p;
break;
}
}
}
}
if (!current) {
if (in.second && find(colour_ordered.begin(),colour_ordered.end(),in.second) ==
colour_ordered.end() && in.second->coloured()) {
if (!in.second->hasColour() || !in.second->hasAntiColour())
current = in.second;
}
}
// then go on with anything else
if (!current) {
if (find(colour_ordered.begin(),colour_ordered.end(),in.first) ==
colour_ordered.end() && in.first->coloured()) {
current = in.first;
}
}
if (!current) {
for (PList::iterator p = out.begin();
p != out.end(); ++p) {
if (find(colour_ordered.begin(),colour_ordered.end(),*p) ==
colour_ordered.end() && (**p).coloured()) {
current = *p;
break;
}
}
}
if (!current) {
if (in.second && find(colour_ordered.begin(),colour_ordered.end(),in.second) ==
colour_ordered.end() && in.second->coloured()) {
current = in.second;
}
}
assert(current);
PPtr next;
Ptr<ColourLine>::ptr walk_the_line;
while (true) {
if (!walk_the_line) {
if (current->hasColour()) {
walk_the_line = current->colourLine();
}
else if (current->hasAntiColour()) {
walk_the_line = current->antiColourLine();
}
}
if (!next)
for (tPVector::const_iterator p = walk_the_line->coloured().begin();
p != walk_the_line->coloured().end(); ++p) {
if (*p == current)
continue;
if (find(out.begin(),out.end(),*p) != out.end() ||
*p == in.first ||
(in.second && *p == in.second)) {
next = *p;
if (next->hasColour() && next->hasAntiColour()) {
walk_the_line = walk_the_line == next->colourLine() ? next->antiColourLine() : next->colourLine();
}
break;
}
}
if (!next)
for (tPVector::const_iterator p = walk_the_line->antiColoured().begin();
p != walk_the_line->antiColoured().end(); ++p) {
if (*p == current)
continue;
if (find(out.begin(),out.end(),*p) != out.end() ||
*p == in.first ||
(in.second && *p == in.second)) {
next = *p;
if (next->hasColour() && next->hasAntiColour()) {
walk_the_line = walk_the_line == next->colourLine() ? next->antiColourLine() : next->colourLine();
}
break;
}
}
assert(next);
colour_ordered.push_back(current);
current = next;
// done if next is not a gluon or next is already in colour_ordered
if ((current->hasColour() && !current->hasAntiColour()) ||
(!current->hasColour() && current->hasAntiColour())) {
colour_ordered.push_back(current);
break;
}
if (next->hasColour() && next->hasAntiColour()) {
if (find(colour_ordered.begin(),colour_ordered.end(),next) != colour_ordered.end())
break;
}
next = PPtr();
}
}
return colour_ordered;
}
void DipoleEventRecord::popChain() {
assert(!theChains.empty());
theDoneChains.push_back(DipoleChain());
theDoneChains.back().dipoles().splice(theDoneChains.back().dipoles().begin(),theChains.front().dipoles());
theChains.pop_front();
}
void DipoleEventRecord::popChain(list<DipoleChain>::iterator ch) {
assert(!theChains.empty());
theDoneChains.push_back(DipoleChain());
theDoneChains.back().dipoles().splice(theDoneChains.back().dipoles().begin(),ch->dipoles());
theChains.erase(ch);
}
void DipoleEventRecord::popChains(const list<list<DipoleChain>::iterator>& chs) {
assert(!theChains.empty());
for ( list<list<DipoleChain>::iterator>::const_iterator ch =
chs.begin(); ch != chs.end(); ++ch ) {
theDoneChains.push_back(DipoleChain());
theDoneChains.back().dipoles().splice(theDoneChains.back().dipoles().begin(),(*ch)->dipoles());
}
for ( list<list<DipoleChain>::iterator>::const_iterator ch =
chs.begin(); ch != chs.end(); ++ch )
theChains.erase(*ch);
}
DipoleIndex
DipoleEventRecord::mergeIndex(list<Dipole>::iterator firstDipole, const pair<bool,bool>& whichFirst,
list<Dipole>::iterator secondDipole, const pair<bool,bool>& whichSecond) const {
tcPDPtr emitterData =
whichFirst.first ? firstDipole->leftParticle()->dataPtr() : firstDipole->rightParticle()->dataPtr();
tcPDPtr spectatorData =
whichSecond.first ? secondDipole->leftParticle()->dataPtr() : secondDipole->rightParticle()->dataPtr();
const PDF& emitterPDF =
whichFirst.first ? firstDipole->leftPDF() : firstDipole->rightPDF();
const PDF& spectatorPDF =
whichSecond.first ? secondDipole->leftPDF() : secondDipole->rightPDF();
return DipoleIndex(emitterData,spectatorData,emitterPDF,spectatorPDF);
}
SubleadingSplittingInfo
DipoleEventRecord::mergeSplittingInfo(list<DipoleChain>::iterator firstChain, list<Dipole>::iterator firstDipole,
const pair<bool,bool>& whichFirst,
list<DipoleChain>::iterator secondChain, list<Dipole>::iterator secondDipole,
const pair<bool,bool>& whichSecond) const {
SubleadingSplittingInfo res;
res.index(mergeIndex(firstDipole,whichFirst,secondDipole,whichSecond));
res.emitter(whichFirst.first ? firstDipole->leftParticle() : firstDipole->rightParticle());
res.spectator(whichSecond.first ? secondDipole->leftParticle() : secondDipole->rightParticle());
res.emitterX(whichFirst.first ? firstDipole->leftFraction() : firstDipole->rightFraction());
res.spectatorX(whichSecond.first ? secondDipole->leftFraction() : secondDipole->rightFraction());
res.configuration(whichFirst);
res.spectatorConfiguration(whichSecond);
res.emitterChain(firstChain);
res.emitterDipole(firstDipole);
res.spectatorChain(secondChain);
res.spectatorDipole(secondDipole);
return res;
}
void DipoleEventRecord::getSubleadingSplittings(list<SubleadingSplittingInfo>& res) {
static pair<bool,bool> left(true,false);
static pair<bool,bool> right(false,true);
res.clear();
for ( list<DipoleChain>::iterator cit = theChains.begin();
cit != theChains.end(); ++cit ) {
for ( list<Dipole>::iterator dit = cit->dipoles().begin();
dit != cit->dipoles().end(); ++dit ) {
for ( list<Dipole>::iterator djt = dit;
djt != cit->dipoles().end(); ++djt ) {
res.push_back(mergeSplittingInfo(cit,dit,left,cit,djt,left));
res.push_back(mergeSplittingInfo(cit,dit,right,cit,djt,right));
if ( dit != djt ) {
res.push_back(mergeSplittingInfo(cit,dit,left,cit,djt,right));
res.push_back(mergeSplittingInfo(cit,dit,right,cit,djt,left));
}
}
}
list<DipoleChain>::iterator cjt = cit; ++cjt;
for ( ; cjt != theChains.end(); ++cjt ) {
for ( list<Dipole>::iterator dit = cit->dipoles().begin();
dit != cit->dipoles().end(); ++dit ) {
for ( list<Dipole>::iterator djt = cjt->dipoles().begin();
djt != cjt->dipoles().end(); ++djt ) {
res.push_back(mergeSplittingInfo(cit,dit,left,cjt,djt,left));
res.push_back(mergeSplittingInfo(cit,dit,right,cjt,djt,right));
res.push_back(mergeSplittingInfo(cit,dit,left,cjt,djt,right));
res.push_back(mergeSplittingInfo(cit,dit,right,cjt,djt,left));
}
}
}
}
}
void DipoleEventRecord::splitSubleading(SubleadingSplittingInfo& dsplit,
pair<list<Dipole>::iterator,list<Dipole>::iterator>& childIterators,
DipoleChain*& firstChain, DipoleChain*& secondChain) {
if ( dsplit.emitterDipole() == dsplit.spectatorDipole() ) {
assert(dsplit.emitterChain() == dsplit.spectatorChain());
split(dsplit.emitterDipole(),dsplit.emitterChain(),dsplit,
childIterators,firstChain,secondChain,false);
} else {
// first need to recoil, then split
recoil(dsplit.spectatorDipole(),dsplit.spectatorChain(),dsplit);
split(dsplit.emitterDipole(),dsplit.emitterChain(),dsplit,
childIterators,firstChain,secondChain,true);
}
}
void DipoleEventRecord::findChains(const PList& ordered, const bool decay) {
theChains.clear();
theDoneChains.clear();
DipoleChain current_chain;
// this whole thing needs to have a more elegant implementation at some point
bool startIsTriplet =
(ordered.front()->hasColour() && !ordered.front()->hasAntiColour()) ||
(!ordered.front()->hasColour() && ordered.front()->hasAntiColour());
bool endIsTriplet =
(ordered.back()->hasColour() && !ordered.back()->hasAntiColour()) ||
(!ordered.back()->hasColour() && ordered.back()->hasAntiColour());
if (!( ordered.size() == 2 && startIsTriplet && endIsTriplet)) {
PList::const_iterator theStart = ordered.begin();
bool onceMore = false;
for (PList::const_iterator p = ordered.begin();
p != ordered.end(); ++p) {
PList::const_iterator next_it =
p != --ordered.end() ? std::next(p) : ordered.begin();
if (!DipolePartonSplitter::colourConnected(*p,*next_it)) {
// it may have happened that we need to close the chain due to another
// chain starting right now; see the above global comment for this fix
bool startIsOctet =
(**theStart).hasColour() && (**theStart).hasAntiColour();
bool endIsOctet =
(**p).hasColour() && (**p).hasAntiColour();
if ( DipolePartonSplitter::colourConnected(*p,*theStart) &&
startIsOctet && endIsOctet ) {
swap(next_it,theStart);
onceMore = true;
} else {
theStart = next_it;
current_chain.check();
theChains.push_back(current_chain);
current_chain.dipoles().clear();
continue;
}
}
pair<bool,bool> initial_state (false,false);
initial_state.first = (*p == incoming().first || *p == incoming().second);
initial_state.second = (*next_it == incoming().first || *next_it == incoming().second);
pair<int,int> which_in (-1,-1);
if (initial_state.first)
which_in.first = *p == incoming().first ? 0 : 1;
if (initial_state.second)
which_in.second = *next_it == incoming().first ? 0 : 1;
pair<double,double> xs (1.,1.);
if (initial_state.first)
xs.first = *p == incoming().first ? fractions().first : fractions().second;
if (initial_state.second)
xs.second = *next_it == incoming().first ? fractions().first : fractions().second;
pair<PDF,PDF> pdf;
if ( which_in.first == 0 )
pdf.first = pdfs().first;
else if ( which_in.first == 1 )
pdf.first = pdfs().second;
if ( which_in.second == 0 )
pdf.second = pdfs().first;
else if ( which_in.second == 1 )
pdf.second = pdfs().second;
// In the case of a decay process register which
// parton is incoming to the decay
pair<bool,bool> decayed_parton (false,false);
if (decay) {
decayed_parton.first = (*p == currentDecay()->incoming()[0].first);
decayed_parton.second = (*next_it == currentDecay()->incoming()[0].first);
}
current_chain.dipoles().push_back(Dipole({*p,*next_it},pdf,xs,decayed_parton));
if ( onceMore ) {
next_it = theStart;
current_chain.check();
theChains.push_back(current_chain);
current_chain.dipoles().clear();
onceMore = false;
}
}
}
else {
// treat 2 -> singlet, singlet -> 2 and 1 + singlet -> 1 + singlet special
// to prevent duplicate dipole
assert(DipolePartonSplitter::colourConnected(ordered.front(),ordered.back()));
pair<bool,bool> initial_state (false,false);
initial_state.first = (ordered.front() == incoming().first || ordered.front() == incoming().second);
initial_state.second = (ordered.back() == incoming().first || ordered.back() == incoming().second);
pair<int,int> which_in (-1,-1);
if (initial_state.first)
which_in.first = ordered.front() == incoming().first ? 0 : 1;
if (initial_state.second)
which_in.second = ordered.back() == incoming().first ? 0 : 1;
pair<double,double> xs (1.,1.);
if (initial_state.first)
xs.first = ordered.front() == incoming().first ? fractions().first : fractions().second;
if (initial_state.second)
xs.second = ordered.back() == incoming().first ? fractions().first : fractions().second;
pair<PDF,PDF> pdf;
if ( which_in.first == 0 )
pdf.first = pdfs().first;
else if ( which_in.first == 1 )
pdf.first = pdfs().second;
if ( which_in.second == 0 )
pdf.second = pdfs().first;
else if ( which_in.second == 1 )
pdf.second = pdfs().second;
// In the case of a decay process register which
// parton is incoming to the decay
pair<bool,bool> decayed_parton (false,false);
if (decay) {
decayed_parton.first = (ordered.front() == currentDecay()->incoming()[0].first);
decayed_parton.second = (ordered.back() == currentDecay()->incoming()[0].first);
}
current_chain.dipoles().push_back(Dipole({ordered.front(),ordered.back()},pdf,xs,decayed_parton));
}
if (!current_chain.dipoles().empty()) {
current_chain.check();
theChains.push_back(current_chain);
}
}
const map<PPtr,PPtr>&
DipoleEventRecord::prepare(tSubProPtr subpro,
tStdXCombPtr xc,
const pair<PDF,PDF>& pdf,tPPair beam,
bool dipoles) {
// set the subprocess
subProcess(subpro);
// clear the event record
outgoing().clear();
theHard.clear();
theOriginals.clear();
theDecays.clear();
theCurrentDecay = PerturbativeProcessPtr();
// extract incoming particles
PPair in = subpro->incoming();
// get the incoming momentum fractions
// don't take these from the XComb as it may be null
pair<double,double> xs;
ThePEG::Direction<0> dir(true);
xs.first = in.first->momentum().dirPlus()/beam.first->momentum().dirPlus();
dir.reverse();
xs.second = in.second->momentum().dirPlus()/beam.second->momentum().dirPlus();
xcombPtr(xc);
pdfs() = pdf;
fractions() = xs;
// use ShowerHandler to split up the hard process
PerturbativeProcessPtr hard;
DecayProcessMap decay;
ShowerHandler::currentHandler()->splitHardProcess(tPVector(subpro->outgoing().begin(),
subpro->outgoing().end()),
hard,decay);
// vectors for originals and copies of the particles
vector<PPtr> original;
vector<PPtr> copies;
// fill originals
for(unsigned int ix=0;ix<2;++ix)
original.push_back(hard->incoming()[ix].first);
for(unsigned int ix=0;ix<hard->outgoing().size();++ix)
original.push_back(hard->outgoing()[ix].first);
for(DecayProcessMap::const_iterator it=decay.begin();it!=decay.end();++it) {
fillFromDecays(it->second, original);
}
// and make copies
for ( vector<PPtr>::const_iterator p = original.begin();
p != original.end(); ++p ) {
PPtr copy = new_ptr(Particle(**p));
copies.push_back(copy);
theOriginals[*p] = copy;
}
// isolate the colour of the copies from the originals
colourIsolate(original,copies);
// set the incoming particles
incoming().first = copies[0];
ParticleVector children = incoming().first->children();
for ( ParticleVector::const_iterator c = children.begin();
c != children.end(); ++c )
incoming().first->abandonChild(*c);
incoming().second = copies[1];
children = incoming().second->children();
for ( ParticleVector::const_iterator c = children.begin();
c != children.end(); ++c )
incoming().second->abandonChild(*c);
// set the outgoing particles for the hard process
for(unsigned int ix=0;ix<hard->outgoing().size();++ix) {
if(hard->outgoing()[ix].first->coloured())
outgoing().push_back(theOriginals[hard->outgoing()[ix].first]);
else
theHard.push_back(theOriginals[hard->outgoing()[ix].first]);
}
if ( dipoles ) {
PList cordered = colourOrdered(incoming(),outgoing());
findChains(cordered,false);
}
// sort out the decays
for(auto const & dec : decay) {
// If the decay particle is in original it needs
// to be added to the decays and the decay needs to be
// changed to the copied particles.
if ( theOriginals.find(dec.second->incoming()[0].first) != theOriginals.end() ) {
theDecays[theOriginals[dec.second->incoming()[0].first]] = dec.second;
PerturbativeProcessPtr decayProc = theDecays[theOriginals[dec.second->incoming()[0].first]];
separateDecay(decayProc);
}
else {
assert( find( copies.begin(), copies.end(), dec.second->incoming()[0].first ) != copies.end() );
theDecays[dec.second->incoming()[0].first] = dec.second;
}
}
PList::const_iterator XFirst, XLast;
if ( !theHard.empty() ) {
XFirst = theHard.begin();
XLast = theHard.end();
} else {
XFirst = outgoing().begin();
XLast = outgoing().end();
}
thePX = (**XFirst).momentum();
++XFirst;
for ( ; XFirst != XLast; ++XFirst )
thePX += (**XFirst).momentum();
identifyEventType();
return theOriginals;
}
void DipoleEventRecord::slimprepare(tSubProPtr subpro,
tStdXCombPtr xc,
const pair<PDF,PDF>& pdf,tPPair beam,
bool dipoles) {
// set the subprocess
subProcess(subpro);
// clear the event record
outgoing().clear();
theHard.clear();
theOriginals.clear();
theDecays.clear();
theCurrentDecay = PerturbativeProcessPtr();
// extract incoming particles
PPair in = subpro->incoming();
// get the beam
// get the incoming momentum fractions
// don't take these from the XComb as it may be null
pair<double,double> xs;
ThePEG::Direction<0> dir(true);
xs.first = in.first->momentum().dirPlus()/beam.first->momentum().dirPlus();
dir.reverse();
xs.second = in.second->momentum().dirPlus()/beam.second->momentum().dirPlus();
xcombPtr(xc);
pdfs() = pdf;
fractions() = xs;
incoming() = in;
for(unsigned int ix=0;ix<subpro->outgoing().size();++ix) {
if(subpro->outgoing()[ix]->coloured())
outgoing().push_back(subpro->outgoing()[ix]);
}
if ( dipoles ) {
PList cordered = colourOrdered(incoming(),outgoing());
findChains(cordered,false);
}
}
void DipoleEventRecord::fillFromDecays(PerturbativeProcessPtr decayProc, vector<PPtr>& original) {
// Loop over the outgoing of the given perturbative process
for ( auto const & outIt : decayProc->outgoing() ) {
// Add the outgoing particle to the vector of original particles
original.push_back(outIt.first);
// Iterate through the outgoing
if ( outIt.second )
fillFromDecays( outIt.second, original);
}
}
void DipoleEventRecord::separateDecay(PerturbativeProcessPtr decayProc) {
// Iteratively replace all entries in the incoming
// with their copies.
for ( auto & inIt : decayProc->incoming() ) {
if ( theOriginals.find( inIt.first ) != theOriginals.end() )
inIt.first = theOriginals[inIt.first];
}
// Iteratively replace all entries in the outgoing
// with their copies.
for ( auto & outIt : decayProc->outgoing()) {
if ( theOriginals.count( outIt.first ) )
outIt.first = theOriginals[outIt.first];
if ( outIt.second )
separateDecay(outIt.second);
}
}
void DipoleEventRecord::clear() {
ShowerEventRecord::clear();
theDecays.clear();
theHard.clear();
theChains.clear();
theDoneChains.clear();
theOriginals.clear();
}
pair<PVector,PVector> DipoleEventRecord::tmpupdate(DipoleSplittingInfo& dsplit) {
PVector inc;
PVector out;
tcPPtr IF = incoming().first;
tcPPtr IS = incoming().second;
tcPPtr DE = dsplit.emitter();
tcPPtr DS = dsplit.spectator();
if ( IF != DE && IF != DS ) {
PPtr p = IF->data().produceParticle(IF->momentum());
inc.push_back(p);
}
else if ( IF == DE ) inc.push_back( dsplit.splitEmitter() );
else if ( IF == DS ) inc.push_back( dsplit.splitSpectator() );
if ( IS != DE && IS != DS ) {
PPtr p = IS->data().produceParticle(IS->momentum());
inc.push_back(p);
}
else if ( IS == DE ) inc.push_back( dsplit.splitEmitter() );
else if ( IS == DS ) inc.push_back( dsplit.splitSpectator() );
if ( IF != DE && IS != DE)
out.push_back( dsplit.splitEmitter());
if ( IF != DS && IS != DS)
out.push_back( dsplit.splitSpectator());
out.push_back( dsplit.emission());
for ( tcPPtr h : theHard ){
PPtr p = h->data().produceParticle(h->momentum());
if ( dsplit.splittingKinematics()->doesTransform() ) {
p->set5Momentum( dsplit.splittingKinematics()->transform(p->momentum()) );
}
out.push_back(p);
}
for ( tcPPtr p : outgoing() )
if ( p != DE &&
p != DS &&
p != dsplit.emission() ){
PPtr ou = p->data().produceParticle(p->momentum());;
if ( dsplit.splittingKinematics()->doesTransform() ){
ou->set5Momentum( dsplit.splittingKinematics()->transform(ou->momentum()) );
}
out.push_back(ou);
}
return {inc,out};
}
void DipoleEventRecord::update(DipoleSplittingInfo& dsplit) {
if ( incoming().first == dsplit.emitter() ) {
intermediates().push_back(dsplit.emitter());
incoming().first = dsplit.splitEmitter();
fractions().first /= dsplit.lastEmitterZ();
} else if ( incoming().first == dsplit.spectator() ) {
intermediates().push_back(dsplit.spectator());
incoming().first = dsplit.splitSpectator();
fractions().first /= dsplit.lastSpectatorZ();
}
if ( incoming().second == dsplit.emitter() ) {
intermediates().push_back(dsplit.emitter());
incoming().second = dsplit.splitEmitter();
fractions().second /= dsplit.lastEmitterZ();
} else if ( incoming().second == dsplit.spectator() ) {
intermediates().push_back(dsplit.spectator());
incoming().second = dsplit.splitSpectator();
fractions().second /= dsplit.lastSpectatorZ();
}
PList::iterator pos;
pos = find(outgoing().begin(), outgoing().end(), dsplit.emitter());
if (pos != outgoing().end()) {
intermediates().push_back(*pos);
*pos = dsplit.splitEmitter();
}
pos = find(outgoing().begin(), outgoing().end(), dsplit.spectator());
if (pos != outgoing().end()) {
intermediates().push_back(*pos);
*pos = dsplit.splitSpectator();
}
outgoing().push_back(dsplit.emission());
if (dsplit.splittingKinematics()->doesTransform()) {
for (PList::iterator p = intermediates().begin();
p != intermediates().end(); ++p) {
(**p).set5Momentum(dsplit.splittingKinematics()->transform((**p).momentum()));
}
for (PList::iterator h = theHard.begin();
h != theHard.end(); ++h) {
(**h).set5Momentum(dsplit.splittingKinematics()->transform((**h).momentum()));
}
for (PList::iterator p = outgoing().begin();
p != outgoing().end(); ++p) {
if ((*p) != dsplit.splitEmitter() &&
(*p) != dsplit.splitSpectator() &&
(*p) != dsplit.emission())
(**p).set5Momentum(dsplit.splittingKinematics()->transform((**p).momentum()));
}
}
// Handle updates related to decays
// Showering of decay processes
// Treat the evolution of the incoming
// decayed particle as in backward evolution
if ( dsplit.isDecayProc() ) {
// Create a pointer to the decay process
PerturbativeProcessPtr decayProc = currentDecay();
// Add the emission to the outgoing of the decay process
decayProc->outgoing().push_back( {dsplit.emission(), PerturbativeProcessPtr() });
// Bools to be used throughout
const bool decayedEmtr = dsplit.index().incomingDecayEmitter();
const bool decayedSpec = dsplit.index().incomingDecaySpectator();
/*
In the current implementation, **following the hard process**
all particles in theDecays evolve independently
e.g. if we have W -> XYZ where all X, Y and Z need to be
showered and decayed, we only identify them as needing decaying
(and hence put them in theDecays) AFTER showering the decay of W.
Hence, XYZ are not even in theDecays until W has been fully
showered and then they are decayed and showered completely independently
KEY POINT - Never need to update other entries of theDecays
Note: The PPtr in theDecays should remain unchanged and all changes
should be made to the relative PerturbativeProcess.
*/
// Splittings from dipoles in the decay process which
// do not have the decayed parton as emitter or spectator.
// Update the decay process in theDecays
if ( !decayedEmtr && !decayedSpec ) {
// Find and replace the old spectator and
// emitter in the outgoing of the decay process
bool decayProcEm = false;
bool decayProcSp = false;
for ( auto & outIt : decayProc->outgoing() ) {
if ( !decayProcEm && outIt.first == dsplit.emitter() ) {
outIt = {dsplit.splitEmitter(), PerturbativeProcessPtr()};
decayProcEm = true;
}
if ( !decayProcSp && outIt.first == dsplit.spectator() ) {
outIt = {dsplit.splitSpectator(), PerturbativeProcessPtr() };
decayProcSp = true;
}
if ( decayProcEm && decayProcSp )
break;
}
// Test that nothing strange is happening
assert( (decayProcEm && decayProcSp) );
return;
}
// The spectator is the decayed particle
else if ( decayedSpec ) {
// Update the dipole event record intermediates
intermediates().push_back(dsplit.splitSpectator());
// Update the the decayProcess incoming
decayProc->incoming().clear();
decayProc->incoming().push_back({dsplit.splitSpectator(),decayProc});
// Update the decay process outgoing
// Replace the old emitter with the new emitter
for ( auto & outEmtrIt : decayProc->outgoing() ) {
if ( outEmtrIt.first == dsplit.emitter() ){
outEmtrIt = {dsplit.splitEmitter(), PerturbativeProcessPtr() };
break;
}
}
// Perform the recoil transformation
// Find all particles in the recoil system
PList recoilSystem;
for ( auto const & outIt : decayProc->outgoing() ) {
if ( outIt.first != dsplit.splitEmitter() && outIt.first != dsplit.emission() ) {
recoilSystem.push_back(outIt.first);
}
}
dsplit.splittingKinematics()->decayRecoil( recoilSystem );
return;
}
// The emitter is the decayed particle
else {
throw Exception()
<< "DipoleEventRecord: The emitter as a decayed particle is currently not implemented."
<< Exception::runerror;
assert( currentDecay()->incoming()[0].first == dsplit.emitter() && decayedEmtr && !decayedSpec );
// Update the dipole event record intermediates
intermediates().push_back(dsplit.splitEmitter());
// Update the the decayProcess incoming
decayProc->incoming().clear();
decayProc->incoming().push_back({dsplit.splitEmitter(),decayProc});
// Update the decay process outgoing
// Replace the old spectator with the new spectator
for (auto & outSpecIt : decayProc->outgoing() ) {
if ( outSpecIt.first == dsplit.spectator() ){
outSpecIt = { dsplit.splitSpectator(), PerturbativeProcessPtr() };
break;
}
}
// Perform the recoil transformation
assert(dsplit.splittingKinematics()->isDecay());
// Find all particles in the recoil system
PList recoilSystem;
for ( auto const & outIt : decayProc->outgoing() ) {
if ( outIt.first != dsplit.splitSpectator() && outIt.first != dsplit.emission() ) {
recoilSystem.push_back(outIt.first);
}
}
dsplit.splittingKinematics()->decayRecoil( recoilSystem );
return;
}
}
}
void
DipoleEventRecord::split(list<Dipole>::iterator dip,
list<DipoleChain>::iterator ch,
DipoleSplittingInfo& dsplit,
pair<list<Dipole>::iterator,list<Dipole>::iterator>& childIterators,
DipoleChain*& firstChain, DipoleChain*& secondChain,
bool colourSpectator) {
static DipoleChain empty;
pair<Dipole,Dipole> children = dip->split(dsplit,colourSpectator);
list<Dipole>::iterator breakup =
ch->insertSplitting(dip,children,childIterators);
if ( breakup == ch->dipoles().end() ) {
firstChain = &(*ch);
secondChain = &empty;
}
else {
DipoleChain other;
other.dipoles().splice(other.dipoles().end(),ch->dipoles(),breakup,ch->dipoles().end());
chains().push_back(other);
firstChain = &(*ch);
secondChain = &(chains().back());
// explicitly fix iterators in case the splice implementation
// at hand does invalidate iterators (the SGI docu says, it doesn't,
// but it seems that this behaviour is not part of the standard)
childIterators.first = --firstChain->dipoles().end();
childIterators.second = secondChain->dipoles().begin(); }
if ( !colourSpectator ) {
update(dsplit); // otherwise done by recoil(...)
}
}
pair<PVector,PVector> DipoleEventRecord::tmpsplit(list<Dipole>::iterator dip,
list<DipoleChain>::iterator ,
DipoleSplittingInfo& dsplit,
pair<list<Dipole>::iterator,list<Dipole>::iterator>& ,
DipoleChain*& , DipoleChain*& ,
bool colourSpectator) {
dip->tmpsplit(dsplit,colourSpectator);
return tmpupdate(dsplit); // otherwise done by recoil(...)
}
void DipoleEventRecord::recoil(list<Dipole>::iterator dip,
list<DipoleChain>::iterator ch,
DipoleSplittingInfo& dsplit) {
dip->recoil(dsplit);
ch->updateDipole(dip);
update(dsplit);
}
list<pair<list<Dipole>::iterator,list<DipoleChain>::iterator> >
DipoleEventRecord::inDipoles() {
list<pair<list<Dipole>::iterator,list<DipoleChain>::iterator> > res;
for ( list<DipoleChain>::iterator chit = theDoneChains.begin();
chit != theDoneChains.end(); ++chit ) {
bool haveOne = false;
for ( list<Dipole>::iterator dit = chit->dipoles().begin();
dit != chit->dipoles().end(); ++dit ) {
if ( dit->leftPDF().pdf() || dit->rightPDF().pdf() ) {
haveOne = true;
break;
}
}
if ( haveOne ) {
theChains.splice(theChains.begin(),theDoneChains,chit);
for ( list<Dipole>::iterator dit = theChains.front().dipoles().begin();
dit != theChains.front().dipoles().end(); ++dit ) {
if ( dit->leftPDF().pdf() || dit->rightPDF().pdf() ) {
res.push_back({dit,theChains.begin()});
}
}
}
}
return res;
}
void DipoleEventRecord::transform(const SpinOneLorentzRotation& rot) {
Lorentz5Momentum tmp;
for (PList::iterator p = intermediates().begin();
p != intermediates().end(); ++p) {
tmp = (**p).momentum(); tmp = rot * tmp;
(**p).set5Momentum(tmp);
}
for (PList::iterator h = theHard.begin();
h != theHard.end(); ++h) {
tmp = (**h).momentum(); tmp = rot * tmp;
(**h).set5Momentum(tmp);
}
for (PList::iterator p = outgoing().begin();
p != outgoing().end(); ++p) {
tmp = (**p).momentum(); tmp = rot * tmp;
(**p).set5Momentum(tmp);
}
}
tPPair DipoleEventRecord::fillEventRecord(StepPtr step, bool firstInteraction, bool) {
PPtr inSubPro = subProcess()->incoming().first;
PPtr inParticle;
if ( !(inSubPro->parents().empty()) )
inParticle = inSubPro->parents()[0];
else
inParticle = inSubPro;
PPtr inParton = theOriginals[inSubPro];
theOriginals.erase(inSubPro);
updateColour(incoming().first,true);
if ( inParticle != inSubPro )
inParticle->abandonChild(inSubPro);
inParton->addChild(inSubPro);
if ( inParticle != inSubPro )
inParticle->addChild(incoming().first);
intermediates().push_back(inSubPro);
intermediates().push_back(inParton);
// Repeat all the above for the second incoming particle
inSubPro = subProcess()->incoming().second;
if ( !(inSubPro->parents().empty()) )
inParticle = inSubPro->parents()[0];
else
inParticle = inSubPro;
inParton = theOriginals[inSubPro];
theOriginals.erase(inSubPro);
updateColour(incoming().second,true);
if ( inParticle != inSubPro )
inParticle->abandonChild(inSubPro);
inParton->addChild(inSubPro);
if ( inParticle != inSubPro )
inParticle->addChild(incoming().second);
intermediates().push_back(inSubPro);
intermediates().push_back(inParton);
// theOriginals is populated in ::prepare and contains all of the incoming and outgoing particles of the original hard process
// Here outgoing particles from theOriginals are added into the intermediates()
while ( !theOriginals.empty() ) {
PPtr outSubPro = theOriginals.begin()->first;
PPtr outParton = theOriginals.begin()->second;
// workaround for OS X Mavericks LLVM libc++
#ifdef _LIBCPP_VERSION
map<PPtr,PPtr>::const_iterator beg = theOriginals.begin();
#else
map<PPtr,PPtr>::iterator beg = theOriginals.begin();
#endif
theOriginals.erase(beg);
updateColour(outParton,true);
outSubPro->addChild(outParton);
intermediates().push_back(outSubPro);
}
// Update the intermediates of the step
step->addIntermediates(intermediates().begin(),intermediates().end());
for (auto const & p : outgoing())
step->addDecayProduct( p );
for (auto const & p : theHard)
step->addDecayProduct( p );
if ( firstInteraction &&
(incoming().first->coloured() ||
incoming().second->coloured() ) ) {
ShowerHandler::currentHandler()->lastExtractor()
->newRemnants(subProcess()->incoming(),incoming(),step);
}
step->addIntermediate(incoming().first);
step->addIntermediate(incoming().second);
return incoming();
}
bool DipoleEventRecord::prepareDecay( PerturbativeProcessPtr decayProc ) {
// Create objects containing the incoming and outgoing partons,
// required as inputs for colourOrdered.
PList out;
for( auto const & dec : decayProc->outgoing()) {
if(dec.first->coloured()) {
out.push_back(dec.first);
}
}
// Only need to shower if we have coloured outgoing particles
if ( out.empty() )
return false;
else {
// For the incoming, use a PPair containing the incoming and a null pointer
PPair in;
in.first = decayProc->incoming()[0].first;
// Create an ordered list of particles
PList cordered;
cordered = colourOrdered(in,out);
// Find the dipole chains for this decay
findChains(cordered,true);
return true;
}
}
Energy DipoleEventRecord::decay(PPtr incoming, bool& powhegEmission) {
// get the process
PerturbativeProcessPtr process = theDecays[incoming];
assert(process);
//tDMPtr decayMode = new_ptr(DecayMode());
tDMPtr decayMode = DMPtr();
// Do not decay particles that have already been decayed
// Note the herwig decayer deals with colour connections
if ( process->outgoing().empty() ) {
process->incoming()[0].first = incoming;
DecayProcessMap decay;
// Decay the particle, returning a pointer to the decay mode
- decayMode = ShowerHandler::currentHandler()->decay(process,decay);
+ decayMode = ShowerHandler::currentHandler()->decay(process,decay,true);
}
// Sort out the colour connections of particles already decayed
else {
// sort out the colour of the incoming
map<tColinePtr,tColinePtr> cmap;
if(incoming->colourLine())
cmap[process->incoming()[0].first->colourLine()] = incoming->colourLine();
if(incoming->antiColourLine())
cmap[process->incoming()[0].first->antiColourLine()] = incoming->antiColourLine();
// fix colours of outgoing
for(auto const & outg : process->outgoing()) {
map<tColinePtr,tColinePtr>::iterator it =
cmap.find(outg.first->colourLine());
if(it!=cmap.end()) {
ColinePtr c1=outg.first->colourLine();
c1->removeColoured(outg.first);
it->second->addColoured(outg.first);
}
it = cmap.find(outg.first->antiColourLine());
if(it!=cmap.end()) {
ColinePtr c1=outg.first->antiColourLine();
c1->removeAntiColoured(outg.first);
it->second->addAntiColoured(outg.first);
}
}
// swap the incoming
process->incoming()[0].first = incoming;
}
// Set the scale of all particles involved in the decay process to the
// mass of the decaying particle
// Initialise the scale for the evolution of
// the parton shower following the decay
Energy showerScale = ZERO;
// Set the scale for the evolution of the shower
showerScale = process->incoming()[0].first->momentum().m();
Energy2 decayScaleSqr = sqr( showerScale );
process->incoming()[0].first->scale( decayScaleSqr );
for(auto & outg : process->outgoing()) {
outg.first->scale( decayScaleSqr );
}
// Update the decaying particle in the process and the event
PList::iterator posOut = find(outgoing().begin(), outgoing().end(), incoming);
PList::iterator posHard = find(hard().begin(), hard().end(), incoming);
assert((posOut!=outgoing().end() && posHard==hard().end()) ||
(posOut==outgoing().end() && posHard!=hard().end()) );
if ( posOut!=outgoing().end() ) {
outgoing().erase(posOut);
}
else {
hard().erase(posHard);
}
intermediates().push_back(process->incoming()[0].first);
// Populate the children of the incoming
for(auto const & outg : process->outgoing()) {
PPtr outgoing = outg.first;
process->incoming()[0].first->addChild(outgoing);
}
// If a decayed particle is not decayed above,
// e.g. a W in a 3-body top decay, find its decaymode.
if ( powhegEmission && !decayMode ) {
string tag = incoming->dataPtr()->name() + "->";
// Must use OrderedParticles for a tag search
ShowerHandler::OrderedParticles decayOut;
for(auto const & outg : process->outgoing()) {
decayOut.insert(outg.first->dataPtr());
}
// Construct the tag
for(auto const & dec : decayOut) {
if( dec!=*decayOut.begin() ) tag += ",";
tag +=dec->name();
}
tag += ";";
// Find the decay mode
decayMode = ShowerHandler::currentHandler()->findDecayMode(tag);
}
// Perform the powheg emission
if ( powhegEmission ) {
if ( decayMode ) {
HwDecayerBasePtr decayer;
decayer = dynamic_ptr_cast<HwDecayerBasePtr>(decayMode->decayer());
if ( decayer->hasPOWHEGCorrection() ) {
// Construct a real emission process and populate its
// incoming and outcoming prior to any powheg emission
RealEmissionProcessPtr born = new_ptr( RealEmissionProcess() );
born->bornIncoming().push_back( incoming );
for(auto const & outg : process->outgoing()) {
born->bornOutgoing().push_back(outg.first);
}
// Generate any powheg emission, returning 'real'
RealEmissionProcessPtr real = decayer->generateHardest( born );
// If no emission is generated the new incoming and
// outgoing containers will be empty.
// Only do something if an emission is generated.
if ( real && real->incoming().size() != 0 ) {
assert ( real->incoming().size() == 1 );
// Update the decay process
// Note: Do not use the new incoming particle
PPtr oldEmitter;
PPtr newEmitter;
// Use the name recoiler to avoid confusion with
// the spectator in the POWHEGDecayer
// i.e. the recoiler can be coloured or non-coloured
PPtr oldRecoiler;
PPtr newRecoiler;
if ( real->emitter() == 1 ) {
oldEmitter = real->bornOutgoing()[0];
oldRecoiler = real->bornOutgoing()[1];
newEmitter = real->outgoing()[0];
newRecoiler = real->outgoing()[1];
}
else if ( real->emitter() == 2) {
oldEmitter = real->bornOutgoing()[1];
oldRecoiler = real->bornOutgoing()[0];
newEmitter = real->outgoing()[1];
newRecoiler = real->outgoing()[0];
}
PPtr emitted = real->outgoing()[ real->emitted()-1];
// Update the scales
newRecoiler->scale(oldRecoiler->scale());
showerScale = real->pT()[ShowerInteraction::QCD];
newEmitter->scale(sqr(showerScale));
emitted->scale(sqr(showerScale));
// Update the colour flow of the new outgoing particles
// Note the emitted and newEmitter are already colour
// connected by the powheg emission function
emitted->incomingColour(oldEmitter, oldEmitter->id()<0);
if ( newRecoiler->coloured() )
newRecoiler->incomingColour(oldRecoiler, oldRecoiler->id()<0);
// Update the children of the outgoing
oldRecoiler->addChild( newRecoiler );
oldEmitter->addChild( newEmitter );
oldEmitter->addChild( emitted );
// Note: The particles in the pert proc outgoing and both outgoing
// vectors of the real emission proc are in the same order
for(unsigned int ix=0;ix<real->bornOutgoing().size();++ix) {
// Update the decay process
assert(process->outgoing()[ix].first == real->bornOutgoing()[ix]);
process->outgoing()[ix].first = real->outgoing()[ix];
// Add the outgoing from the born
// decay to the event intermediates
intermediates().push_back(real->bornOutgoing()[ix]);
}
// Add the emitted to the outgoing of the decay process
process->outgoing().push_back( { emitted, PerturbativeProcessPtr() } );
}
// No powheg emission occurred:
else
powhegEmission = false;
}
// No powheg emission occurred:
else
powhegEmission = false;
}
// No powheg emission occurred:
else
powhegEmission = false;
}
// Copy the outgoing from the decay
// process to the event record
for(auto const & outg : process->outgoing()) {
if ( outg.first->coloured() )
outgoing().push_back(outg.first);
else
hard().push_back(outg.first);
}
return showerScale;
}
void DipoleEventRecord::updateDecayMom( PPtr decayParent, PerturbativeProcessPtr decayProc ) {
// Only particles that have already been decayed
// should be passed to this function
assert( !(decayProc->outgoing().empty()) );
// Create a list of the children to update their momenta
PList children;
for ( auto const & outg : decayProc->outgoing() ) {
children.push_back( outg.first );
}
// Boost the children
PList::iterator beginChildren = children.begin();
PList::iterator endChildren = children.end();
ThePEG::UtilityBase::setMomentum(beginChildren, endChildren, decayParent->momentum().vect() );
}
void DipoleEventRecord::updateDecayChainMom( PPtr decayParent, PerturbativeProcessPtr decayProc ) {
// Note - this updates the momenta of the
// outgoing of the given decay process
// Update the momenta of the outgoing from this decay
updateDecayMom( decayParent, decayProc );
// Iteratively update the momenta of the rest of the decay chain
for ( auto & outg : decayProc->outgoing() ) {
// If a child has a corresponding pert proc
// then it has decay products
if ( outg.second ) {
for ( auto & dec : theDecays ) {
if(dec.second==outg.second) {
dec.first->setMomentum(outg.first->momentum());
break;
}
}
// Iteratively update any decay products
if ( !outg.second->outgoing().empty() )
updateDecayChainMom( outg.first, outg.second );
}
}
}
void DipoleEventRecord::updateDecays(PerturbativeProcessPtr decayProc, bool iterate) {
// Note - This does not update the momenta of the outgoing
// of decayProc.
// i.e. it is for use following the (non-)showering
// of a decay when the daughter momentum are correct.
// With iterate = true, this updates the rest of the decay chain.
// Loop over the outgoing from this decay
for ( auto & outg : decayProc->outgoing() ) {
if ( outg.second && !outg.second->outgoing().empty() ) {
// Outgoing particles which have already been decayed
PPtr newDecayed = outg.first;
PerturbativeProcessPtr newDecayProc = outg.second;
// Update the outgoing momenta from this decay
updateDecayMom( newDecayed, newDecayProc);
// If this decay is already in theDecays then erase it
for ( auto const & dec : theDecays ) {
if(dec.second==newDecayProc) {
theDecays.erase(dec.first);
break;
}
}
// Add to theDecays
theDecays[newDecayed] = newDecayProc;
//assert(theDecays[newDecayed]->incoming()[0].second==decayProc);
// Iteratively update theDecays from the decay chain
if ( iterate )
updateDecays( newDecayProc );
}
// Deal with any outgoing which need to be decayed
else if ( ShowerHandler::currentHandler()->decaysInShower(outg.first->id()) ) {
PerturbativeProcessPtr newDecay=new_ptr(PerturbativeProcess());
newDecay->incoming().push_back({ outg.first , decayProc } );
theDecays[outg.first] = newDecay;
}
}
}
void DipoleEventRecord::debugLastEvent(ostream& os) const {
bool first = ShowerHandler::currentHandler()->firstInteraction();
os << "--- DipoleEventRecord ----------------------------------------------------------\n";
os << " the " << (first ? "hard" : "secondary") << " subprocess is:\n"
<< (*subProcess());
os << " using PDF's " << pdfs().first.pdf() << " and "
<< pdfs().second.pdf() << "\n";
os << " chains showering currently:\n";
for ( list<DipoleChain>::const_iterator chit = theChains.begin();
chit != theChains.end(); ++chit )
os << (*chit);
os << " chains which finished showering:\n";
for ( list<DipoleChain>::const_iterator chit = theDoneChains.begin();
chit != theDoneChains.end(); ++chit )
os << (*chit);
os << "--------------------------------------------------------------------------------\n";
os << flush;
}
diff --git a/Shower/ShowerHandler.cc b/Shower/ShowerHandler.cc
--- a/Shower/ShowerHandler.cc
+++ b/Shower/ShowerHandler.cc
@@ -1,1081 +1,1092 @@
// -*- C++ -*-
//
// ShowerHandler.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 The Herwig Collaboration
//
// Herwig is licenced under version 2 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/Shower/QTilde/Base/ShowerTree.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() )
particlesDecayInShower_.insert(inputparticlesDecayInShower_.begin(),
inputparticlesDecayInShower_.end());
ShowerTree::_vmin2 = vMin_;
ShowerTree::_spaceTime = includeSpaceTime_;
}
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),
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());
}
ShowerTree::_vmin2 = vMin_;
ShowerTree::_spaceTime = includeSpaceTime_;
}
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_;
}
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_;
}
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 interfaceRestrictPhasespaceOn
(interfaceRestrictPhasespace,
"On",
"Perform phasespace restrictions",
true);
static SwitchOption interfaceRestrictPhasespaceOff
(interfaceRestrictPhasespace,
"Off",
"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 interfaceDoFSROn
(interfaceDoFSR,
"Yes",
"Switch on final state radiation.",
true);
static SwitchOption interfaceDoFSROff
(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 interfaceDoISROn
(interfaceDoISR,
"Yes",
"Switch on initial state radiation.",
true);
static SwitchOption interfaceDoISROff
(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);
}
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();
}
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(CurrentGenerator::current().currentEventHandler()->currentStep()->
find(orig)&&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(in,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) const {
+ 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) {
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,803 +1,807 @@
// -*- C++ -*-
//
// ShowerHandler.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 The Herwig Collaboration
//
// Herwig is licenced under version 2 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_;
}
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 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_; }
/**
* 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;
/**
- * Decay a particle
+ * 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) const;
+ 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);
};
/**
* 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;
//@}
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:
/**
* 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:
/** @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 &);
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_;
//@}
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 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 */

File Metadata

Mime Type
text/x-diff
Expires
Tue, Jan 21, 1:48 AM (1 d, 17 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4243561
Default Alt Text
(106 KB)

Event Timeline