Page MenuHomeHEPForge
Files F7879343

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,1498 +1,1496 @@
- // -*- 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.
- //
+// -*- 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;
+ ++done_size;
if (in.second && in.second->coloured())
- ++done_size;
+ ++done_size;
while (colour_ordered.size() != done_size) {
PPtr current;
- // start with singlets, as long as we have some
+ // 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;
+ 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;
+ current = in.second;
}
}
- // then go on with anything else
+ // 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;
- }
- }
+ 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;
- }
- }
+ 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
+ // 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;
+ 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 ) {
+ 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);
+ 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();
+ whichFirst.first ? firstDipole->leftParticle()->dataPtr() : firstDipole->rightParticle()->dataPtr();
tcPDPtr spectatorData =
- whichSecond.first ? secondDipole->leftParticle()->dataPtr() : secondDipole->rightParticle()->dataPtr();
+ whichSecond.first ? secondDipole->leftParticle()->dataPtr() : secondDipole->rightParticle()->dataPtr();
const PDF& emitterPDF =
- whichFirst.first ? firstDipole->leftPDF() : firstDipole->rightPDF();
+ whichFirst.first ? firstDipole->leftPDF() : firstDipole->rightPDF();
const PDF& spectatorPDF =
- whichSecond.first ? secondDipole->leftPDF() : secondDipole->rightPDF();
+ 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 ) {
+ cit != theChains.end(); ++cit ) {
for ( list<Dipole>::iterator dit = cit->dipoles().begin();
- dit != cit->dipoles().end(); ++dit ) {
+ dit != cit->dipoles().end(); ++dit ) {
for ( list<Dipole>::iterator djt = dit;
- djt != cit->dipoles().end(); ++djt ) {
+ 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 ) {
+ dit != cit->dipoles().end(); ++dit ) {
for ( list<Dipole>::iterator djt = cjt->dipoles().begin();
- djt != cjt->dipoles().end(); ++djt ) {
+ 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
+ // 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
+ // 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());
+ (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());
+ (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
+ // 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();
+ (**theStart).hasColour() && (**theStart).hasAntiColour();
bool endIsOctet =
- (**p).hasColour() && (**p).hasAntiColour();
+ (**p).hasColour() && (**p).hasAntiColour();
if ( DipolePartonSplitter::colourConnected(*p,*theStart) &&
- startIsOctet && endIsOctet ) {
+ 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;
+ which_in.first = *p == incoming().first ? 0 : 1;
if (initial_state.second)
- which_in.second = *next_it == incoming().first ? 0 : 1;
+ 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;
+ xs.first = *p == incoming().first ? fractions().first : fractions().second;
if (initial_state.second)
- xs.second = *next_it == incoming().first ? fractions().first : fractions().second;
+ xs.second = *next_it == incoming().first ? fractions().first : fractions().second;
pair<PDF,PDF> pdf;
if ( which_in.first == 0 )
- pdf.first = pdfs().first;
+ pdf.first = pdfs().first;
else if ( which_in.first == 1 )
- pdf.first = pdfs().second;
+ pdf.first = pdfs().second;
if ( which_in.second == 0 )
- pdf.second = pdfs().first;
+ pdf.second = pdfs().first;
else if ( which_in.second == 1 )
- pdf.second = pdfs().second;
+ pdf.second = pdfs().second;
- // In the case of a decay process register which
- // parton is incoming to the decay
+ // 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
+ // 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;
+ which_in.first = ordered.front() == incoming().first ? 0 : 1;
if (initial_state.second)
- which_in.second = ordered.back() == incoming().first ? 0 : 1;
+ 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;
+ xs.first = ordered.front() == incoming().first ? fractions().first : fractions().second;
if (initial_state.second)
- xs.second = ordered.back() == incoming().first ? fractions().first : fractions().second;
+ xs.second = ordered.back() == incoming().first ? fractions().first : fractions().second;
pair<PDF,PDF> pdf;
if ( which_in.first == 0 )
- pdf.first = pdfs().first;
+ pdf.first = pdfs().first;
else if ( which_in.first == 1 )
- pdf.first = pdfs().second;
+ pdf.first = pdfs().second;
if ( which_in.second == 0 )
- pdf.second = pdfs().first;
+ pdf.second = pdfs().first;
else if ( which_in.second == 1 )
- pdf.second = pdfs().second;
+ pdf.second = pdfs().second;
- // In the case of a decay process register which
- // parton is incoming to the decay
+ // 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
+ // set the subprocess
subProcess(subpro);
- // clear the event record
+ // clear the event record
outgoing().clear();
theHard.clear();
theOriginals.clear();
theDecays.clear();
theCurrentDecay = PerturbativeProcessPtr();
- // extract incoming particles
+ // extract incoming particles
PPair in = subpro->incoming();
- // get the incoming momentum fractions
- // don't take these from the XComb as it may be null
+ // 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
+ // 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
+ // vectors for originals and copies of the particles
vector<PPtr> original;
vector<PPtr> copies;
- // fill originals
+ // fill originals
for(unsigned int ix=0;ix<2;++ix)
- original.push_back(hard->incoming()[ix].first);
+ original.push_back(hard->incoming()[ix].first);
for(unsigned int ix=0;ix<hard->outgoing().size();++ix)
- original.push_back(hard->outgoing()[ix].first);
+ 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
+ // and make copies
for ( vector<PPtr>::const_iterator p = original.begin();
- p != original.end(); ++p ) {
+ 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
+ // isolate the colour of the copies from the originals
colourIsolate(original,copies);
- // set the incoming particles
+ // 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);
+ 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);
+ c != children.end(); ++c )
+ incoming().second->abandonChild(*c);
- // set the outgoing particles for the hard process
+ // 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]);
+ outgoing().push_back(theOriginals[hard->outgoing()[ix].first]);
else
- theHard.push_back(theOriginals[hard->outgoing()[ix].first]);
+ theHard.push_back(theOriginals[hard->outgoing()[ix].first]);
}
if ( dipoles ) {
PList cordered = colourOrdered(incoming(),outgoing());
findChains(cordered,false);
}
- // sort out the decays
+ // 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 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();
+ 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
+ // set the subprocess
subProcess(subpro);
- // clear the event record
+ // clear the event record
outgoing().clear();
theHard.clear();
theOriginals.clear();
theDecays.clear();
theCurrentDecay = PerturbativeProcessPtr();
- // extract incoming particles
+ // 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
+ // 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]);
+ 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
+ // 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
+ // Add the outgoing particle to the vector of original particles
original.push_back(outIt.first);
- // Iterate through the outgoing
+ // Iterate through the outgoing
if ( outIt.second )
- fillFromDecays( outIt.second, original);
+ fillFromDecays( outIt.second, original);
}
}
void DipoleEventRecord::separateDecay(PerturbativeProcessPtr decayProc) {
- // Iteratively replace all entries in the incoming
- // with their copies.
+ // 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];
+ inIt.first = theOriginals[inIt.first];
}
- // Iteratively replace all entries in the outgoing
- // with their copies.
+ // 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];
+ outIt.first = theOriginals[outIt.first];
if ( outIt.second )
- separateDecay(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());
+ out.push_back( dsplit.splitEmitter());
if ( IF != DS && IS != DS)
- out.push_back( dsplit.splitSpectator());
+ 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() ){
+ 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()) );
+ PPtr ou = p->data().produceParticle(p->momentum());;
+ if ( dsplit.splittingKinematics()->doesTransform() ){
+ ou->set5Momentum( dsplit.splittingKinematics()->transform(ou->momentum()) );
+ }
+ out.push_back(ou);
}
- 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()));
+ (**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
+ // 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
+ // Create a pointer to the decay process
PerturbativeProcessPtr decayProc = currentDecay();
- // Add the emission to the outgoing of the decay process
+ // Add the emission to the outgoing of the decay process
decayProc->outgoing().push_back( {dsplit.emission(), PerturbativeProcessPtr() });
- // Bools to be used throughout
+ // 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
+ /*
+ 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.
- */
+ 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
+ // 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
+ // 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;
+ break;
}
- // Test that nothing strange is happening
- assert( (decayProcEm && decayProcSp) );
+ // Test that nothing strange is happening
+ assert( (decayProcEm && decayProcSp) );
return;
}
- // The spectator is the decayed particle
+ // The spectator is the decayed particle
else if ( decayedSpec ) {
- // Update the dipole event record intermediates
+ // Update the dipole event record intermediates
intermediates().push_back(dsplit.splitSpectator());
- // Update the the decayProcess incoming
+ // 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
+ // 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
+ // 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
+ // The emitter is the decayed particle
else {
throw Exception()
- << "DipoleEventRecord: The emitter as a decayed particle is currently not implemented."
- << Exception::runerror;
+ << "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
+ // Update the dipole event record intermediates
intermediates().push_back(dsplit.splitEmitter());
- // Update the the decayProcess incoming
+ // 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
+ // 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
+ // Perform the recoil transformation
assert(dsplit.splittingKinematics()->isDecay());
- // Find all particles in the recoil system
+ // 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);
+ 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)
+ // 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 ) {
+ chit != theDoneChains.end(); ++chit ) {
bool haveOne = false;
for ( list<Dipole>::iterator dit = chit->dipoles().begin();
- dit != chit->dipoles().end(); ++dit ) {
+ 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 ) {
+ 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];
+ inParticle = inSubPro->parents()[0];
else
- inParticle = inSubPro;
+ inParticle = inSubPro;
PPtr inParton = theOriginals[inSubPro];
theOriginals.erase(inSubPro);
updateColour(incoming().first,true);
if ( inParticle != inSubPro )
- inParticle->abandonChild(inSubPro);
+ inParticle->abandonChild(inSubPro);
inParton->addChild(inSubPro);
if ( inParticle != inSubPro )
- inParticle->addChild(incoming().first);
+ inParticle->addChild(incoming().first);
intermediates().push_back(inSubPro);
intermediates().push_back(inParton);
- // Repeat all the above for the second incoming particle
+ // Repeat all the above for the second incoming particle
inSubPro = subProcess()->incoming().second;
if ( !(inSubPro->parents().empty()) )
- inParticle = inSubPro->parents()[0];
+ inParticle = inSubPro->parents()[0];
else
- inParticle = inSubPro;
+ inParticle = inSubPro;
inParton = theOriginals[inSubPro];
theOriginals.erase(inSubPro);
updateColour(incoming().second,true);
if ( inParticle != inSubPro )
- inParticle->abandonChild(inSubPro);
+ inParticle->abandonChild(inSubPro);
inParton->addChild(inSubPro);
if ( inParticle != inSubPro )
- inParticle->addChild(incoming().second);
+ 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()
+ // 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++
+ // 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
+ // Update the intermediates of the step
step->addIntermediates(intermediates().begin(),intermediates().end());
for (auto const & p : outgoing())
- step->addDecayProduct( p );
+ 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);
- }
+ (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.
+ // 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
+ // Only need to shower if we have coloured outgoing particles
if ( out.empty() )
- return false;
+ return false;
else {
- // For the incoming, use a PPair containing the incoming and a null pointer
+ // 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
+ // Create an ordered list of particles
PList cordered;
cordered = colourOrdered(in,out);
- // Find the dipole chains for this decay
+ // Find the dipole chains for this decay
findChains(cordered,true);
return true;
}
}
Energy DipoleEventRecord::decay(PPtr incoming, bool& powhegEmission) {
- // get the process
+ // get the process
PerturbativeProcessPtr process = theDecays[incoming];
assert(process);
- //tDMPtr decayMode = new_ptr(DecayMode());
+ //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
+ // 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
+ // Decay the particle, returning a pointer to the decay mode
decayMode = ShowerHandler::currentHandler()->decay(process,decay);
}
- // Sort out the colour connections of particles already decayed
+ // Sort out the colour connections of particles already decayed
else {
- // sort out the colour of the incoming
+ // sort out the colour of the incoming
map<tColinePtr,tColinePtr> cmap;
if(incoming->colourLine())
- cmap[process->incoming()[0].first->colourLine()] = 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
+ 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());
+ 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
+ // 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
+ // 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
+ // 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
+ // 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
+ // 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
+ // 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 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
+ // 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(ShowerHandler::OrderedParticles::const_iterator it=decayOut.begin(); it!=decayOut.end();++it) {
- if(it!=decayOut.begin()) tag += ",";
- tag +=(**it).name();
+ // Construct the tag
+ for(auto const & dec : decayOut) {
+ if(dec=*decayOut.begin()) tag += ",";
+ tag +=dec->name();
}
tag += ";";
- // Find the decay mode
+ // Find the decay mode
decayMode = ShowerHandler::currentHandler()->findDecayMode(tag);
}
- // Perform the powheg emission
+ // Perform the powheg emission
if ( powhegEmission ) {
if ( decayMode ) {
- HwDecayerBasePtr decayer;
- decayer = dynamic_ptr_cast<HwDecayerBasePtr>(decayMode->decayer());
+ HwDecayerBasePtr decayer;
+ decayer = dynamic_ptr_cast<HwDecayerBasePtr>(decayMode->decayer());
- if ( decayer->hasPOWHEGCorrection() ) {
+ 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 );
+ // 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);
- }
+ for(auto const & outg : process->outgoing()) {
+ born->bornOutgoing().push_back(outg.first);
+ }
- // Generate any powheg emission, returning 'real'
+ // 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 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
+ // 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
+ // 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
+ // 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
+ // 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);
+ newRecoiler->incomingColour(oldRecoiler, oldRecoiler->id()<0);
- // Update the children of the outgoing
+ // 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
+ // 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
+ // 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
+ // 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
+ // Add the emitted to the outgoing of the decay process
process->outgoing().push_back( { emitted, PerturbativeProcessPtr() } );
}
- // No powheg emission occurred:
+ // No powheg emission occurred:
else
- powhegEmission = false;
+ powhegEmission = false;
}
- // No powheg emission occurred:
+ // No powheg emission occurred:
else
- powhegEmission = false;
+ powhegEmission = false;
}
- // No powheg emission occurred:
+ // No powheg emission occurred:
else
- powhegEmission = false;
+ powhegEmission = false;
}
- // Copy the outgoing from the decay
- // process to the event record
+ // 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);
+ outgoing().push_back(outg.first);
else
- hard().push_back(outg.first);
+ 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
+ // 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
+ // Create a list of the children to update their momenta
PList children;
- for ( vector<pair<PPtr,PerturbativeProcessPtr> >::iterator outIt = decayProc->outgoing().begin();
- outIt != decayProc->outgoing().end(); ++outIt ) {
- children.push_back( outIt->first );
+ for ( auto const & out : decayProc->outgoing() ) {
+ children.push_back( out.first );
}
- // Boost the children
+ // 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
+ // Note - this updates the momenta of the
+ // outgoing of the given decay process
- // Update the momenta of the outgoing from this decay
+ // Update the momenta of the outgoing from this decay
updateDecayMom( decayParent, decayProc );
- // Iteratively update the momenta of the rest of the decay chain
- for ( vector<pair<PPtr,PerturbativeProcessPtr> >::iterator outIt = decayProc->outgoing().begin();
- outIt != decayProc->outgoing().end(); ++outIt ) {
-
- // If a child has a corresponding pert proc
- // then it has decay products
- if ( outIt->second ) {
+ // Iteratively update the momenta of the rest of the decay chain
+ for ( auto const & out : decayProc->outgoing() ) {
+
+ // If a child has a corresponding pert proc
+ // then it has decay products
+ if ( out.second ) {
for(map<PPtr,PerturbativeProcessPtr>::iterator it=theDecays.begin();
it!=theDecays.end();++it) {
- if(it->second==outIt->second) {
- it->first->setMomentum(outIt->first->momentum());
+ if(it->second==out.second) {
+ it->first->setMomentum(out.first->momentum());
break;
}
}
- // Iteratively update any decay products
- if ( !outIt->second->outgoing().empty() )
- updateDecayChainMom( outIt->first, outIt->second );
+ // Iteratively update any decay products
+ if ( !out.second->outgoing().empty() )
+ updateDecayChainMom( out.first, out.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.
+ // 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 ( vector<pair<PPtr,PerturbativeProcessPtr> >::iterator outIt = decayProc->outgoing().begin(); outIt!=decayProc->outgoing().end(); ++outIt ) {
+ // Loop over the outgoing from this decay
+ for ( auto const & out : decayProc->outgoing() ) {
- if ( outIt->second && !outIt->second->outgoing().empty() ) {
- // Outgoing particles which have already been decayed
- PPtr newDecayed = outIt->first;
- PerturbativeProcessPtr newDecayProc = outIt->second;
+ if ( out.second && !out.second->outgoing().empty() ) {
+ // Outgoing particles which have already been decayed
+ PPtr newDecayed = out.first;
+ PerturbativeProcessPtr newDecayProc = out.second;
- // Update the outgoing momenta from this decay
+ // Update the outgoing momenta from this decay
updateDecayMom( newDecayed, newDecayProc);
- // If this decay is already in theDecays then erase it
+ // If this decay is already in theDecays then erase it
for(map<PPtr,PerturbativeProcessPtr>::const_iterator it=theDecays.begin();
it!=theDecays.end();++it) {
- if(it->second==outIt->second) {
+ if(it->second==out.second) {
theDecays.erase(it->first);
break;
}
}
- // Add to theDecays
+ // Add to theDecays
theDecays[newDecayed] = newDecayProc;
- //assert(theDecays[newDecayed]->incoming()[0].second==decayProc);
+ //assert(theDecays[newDecayed]->incoming()[0].second==decayProc);
- // Iteratively update theDecays from the decay chain
+ // Iteratively update theDecays from the decay chain
if ( iterate )
- updateDecays( newDecayProc );
+ updateDecays( newDecayProc );
}
- // Deal with any outgoing which need to be decayed
- else if ( ShowerHandler::currentHandler()->decaysInShower(outIt->first->id()) ) {
+ // Deal with any outgoing which need to be decayed
+ else if ( ShowerHandler::currentHandler()->decaysInShower(out.first->id()) ) {
PerturbativeProcessPtr newDecay=new_ptr(PerturbativeProcess());
- newDecay->incoming().push_back({ outIt->first , decayProc } );
- theDecays[outIt->first] = newDecay;
+ newDecay->incoming().push_back({ out.first , decayProc } );
+ theDecays[out.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());
+ << (*subProcess());
os << " using PDF's " << pdfs().first.pdf() << " and "
- << pdfs().second.pdf() << "\n";
+ << pdfs().second.pdf() << "\n";
os << " chains showering currently:\n";
for ( list<DipoleChain>::const_iterator chit = theChains.begin();
- chit != theChains.end(); ++chit )
- os << (*chit);
+ 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);
+ chit != theDoneChains.end(); ++chit )
+ os << (*chit);
os << "--------------------------------------------------------------------------------\n";
os << flush;
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 8:01 PM (1 d, 7 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3799762
Default Alt Text
(66 KB)

Event Timeline