Page MenuHomeHEPForge
Files F7879073

No OneTemporary
Actions

View Options

diff --git a/Handlers/StdDependentXComb.cc b/Handlers/StdDependentXComb.cc
--- a/Handlers/StdDependentXComb.cc
+++ b/Handlers/StdDependentXComb.cc
@@ -1,164 +1,165 @@
// -*- C++ -*-
//
// StdDependentXComb.cc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2007 Leif Lonnblad
// Copyright (C) 2009-2010 Simon Platzer
//
// ThePEG 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 StdDependentXComb class.
//
#include "StdDependentXComb.h"
#include "StdXCombGroup.h"
#include "ThePEG/MatrixElement/MEGroup.h"
#include "ThePEG/Cuts/Cuts.h"
#include "ThePEG/Handlers/CascadeHandler.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace ThePEG;
StdDependentXComb::StdDependentXComb(tStdXCombPtr newHead,
const PBPair & newPartonBins, tMEPtr newME,
const DiagramVector & newDiagrams)
: StandardXComb(newHead->maxEnergy(),newHead->particles(),
newHead->eventHandlerPtr(),
const_ptr_cast<tSubHdlPtr>(newHead->subProcessHandler()),
newHead->pExtractor(),newHead->CKKWHandler(),
newPartonBins,newHead->cuts(),newME,newDiagrams,newHead->mirror(),
newHead), resetIncoming(true) {}
StdDependentXComb::StdDependentXComb()
: StandardXComb(), resetIncoming(true) {}
StdDependentXComb::~StdDependentXComb() { }
tSubProPtr StdDependentXComb::construct() {
// first get the meMomenta in their CMS, as this may
// not be the case
Boost cm = (meMomenta()[0] + meMomenta()[1]).findBoostToCM();
if ( cm.mag2() > Constants::epsilon ) {
for ( vector<Lorentz5Momentum>::iterator m = meMomenta().begin();
m != meMomenta().end(); ++m ) {
*m = m->boost(cm);
}
}
tSubProPtr sub = StandardXComb::construct();
sub->head(head()->subProcess());
sub->groupWeight(lastCrossSection()/head()->lastCrossSection());
return sub;
}
void StdDependentXComb::setProcess() {
meMomenta().resize(mePartonData().size());
}
void StdDependentXComb::setPartonBinInstances(Energy2 scale) {
PBIPair newBins;
Direction<0> dir(true);
newBins.first =
new_ptr(PartonBinInstance(lastPartons().first,partonBins().first,scale));
dir.reverse();
newBins.second =
new_ptr(PartonBinInstance(lastPartons().second,partonBins().second,scale));
resetPartonBinInstances(newBins);
setPartonBinInfo();
lastPartons().first->scale(partonBinInstances().first->scale());
lastPartons().second->scale(partonBinInstances().second->scale());
}
void StdDependentXComb::setIncomingPartons() {
if ( !resetIncoming )
return;
resetIncoming = false;
clean();
createPartonBinInstances();
setPartonBinInfo();
lastParticles(head()->lastParticles());
lastPartons(make_pair(mePartonData()[0]->produceParticle(Lorentz5Momentum()),
mePartonData()[1]->produceParticle(Lorentz5Momentum())));
Lorentz5Momentum pFirst = meMomenta()[0];
Lorentz5Momentum pSecond = meMomenta()[1];
if ( head()->matrixElement()->wantCMS() ) {
Boost toLab = (head()->lastPartons().first->momentum() +
head()->lastPartons().second->momentum()).boostVector();
if ( toLab.mag2() > Constants::epsilon ) {
pFirst.boost(toLab);
pSecond.boost(toLab);
}
}
lastPartons().first->set5Momentum(pFirst);
lastPartons().second->set5Momentum(pSecond);
lastS((lastParticles().first->momentum() +
lastParticles().second->momentum()).m2());
lastSHat((lastPartons().first->momentum() +
lastPartons().second->momentum()).m2());
lastP1P2(make_pair(0.0, 0.0));
double x1 =
lastPartons().first->momentum().plus()/
lastParticles().first->momentum().plus();
double x2 =
lastPartons().second->momentum().minus()/
lastParticles().second->momentum().minus();
lastX1X2(make_pair(x1,x2));
lastY(log(lastX1()/lastX2())*0.5);
}
CrossSection StdDependentXComb::dSigDR() {
setProcess();
setIncomingPartons();
lastScale(head()->lastScale());
lastAlphaS(head()->lastAlphaS());
lastAlphaEM(head()->lastAlphaEM());
- if ( !willPassCuts() ) {
+ if ( (!willPassCuts() && !matrixElement()->headCuts()) ||
+ !matrixElement()->apply() ) {
subProcess(SubProPtr());
lastCrossSection(ZERO);
return ZERO;
}
lastPDFWeight(head()->lastPDFWeight());
CrossSection xsec = matrixElement()->dSigHatDR() * lastPDFWeight();
subProcess(SubProPtr());
lastCrossSection(xsec);
return xsec;
}
void StdDependentXComb::Init() {}
void StdDependentXComb::persistentOutput(PersistentOStream &) const {
}
void StdDependentXComb::persistentInput(PersistentIStream &, int) {
}
ClassDescription<StdDependentXComb> StdDependentXComb::initStdDependentXComb;
diff --git a/Handlers/StdXCombGroup.cc b/Handlers/StdXCombGroup.cc
--- a/Handlers/StdXCombGroup.cc
+++ b/Handlers/StdXCombGroup.cc
@@ -1,366 +1,362 @@
// -*- C++ -*-
//
// StdXCombGroup.cc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2007 Leif Lonnblad
// Copyright (C) 2009-2010 Simon Platzer
//
// ThePEG 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 StdXCombGroup class.
//
#include "StdXCombGroup.h"
#include "ThePEG/MatrixElement/MEGroup.h"
#include "ThePEG/Handlers/StandardEventHandler.h"
#include "ThePEG/Handlers/SubProcessHandler.h"
#include "ThePEG/Cuts/Cuts.h"
#include "ThePEG/PDF/PartonExtractor.h"
#include "ThePEG/Utilities/Debug.h"
#include "ThePEG/Utilities/Maths.h"
#include "ThePEG/PDT/ParticleData.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Utilities/SimplePhaseSpace.h"
#include "ThePEG/Utilities/UtilityBase.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/EventRecord/SubProcessGroup.h"
#include "ThePEG/Vectors/LorentzRotation.h"
#include "ThePEG/MatrixElement/MEBase.h"
#include "ThePEG/MatrixElement/ColourLines.h"
#include "ThePEG/Handlers/LuminosityFunction.h"
#include "ThePEG/Handlers/CascadeHandler.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/EventRecord/TmpTransform.h"
using namespace ThePEG;
StdXCombGroup::StdXCombGroup(Energy newMaxEnergy, const cPDPair & inc,
tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler,
tPExtrPtr newExtractor, tCascHdlPtr newCKKW,
const PBPair & newPartonBins, tCutsPtr newCuts, tMEGroupPtr newME,
const DiagramVector & newDiagrams, bool mir)
: StandardXComb(newMaxEnergy,inc,newEventHandler,newSubProcessHandler,
newExtractor, newCKKW, newPartonBins, newCuts,
newME, newDiagrams, mir),
theMEGroup(newME), theDependent(), theLastHeadCrossSection(ZERO) {}
StdXCombGroup::StdXCombGroup()
: StandardXComb(), theDependent() {}
void StdXCombGroup::build(const PartonPairVec& allPBins) {
for ( MEVector::const_iterator me = theMEGroup->dependent().begin();
me != theMEGroup->dependent().end(); ++me ) {
StdDependentXCombPtr dep =
theMEGroup->makeDependentXComb(this,diagrams().front()->partons(),*me,allPBins);
theDependent.push_back(dep);
}
}
StdXCombGroup::~StdXCombGroup() { }
CrossSection StdXCombGroup::dSigDR(const pair<double,double> ll, int nr, const double * r) {
matrixElement()->flushCaches();
if ( matrixElement()->keepRandomNumbers() ) {
lastRandomNumbers().resize(nDim());
copy(r,r+nDim(),lastRandomNumbers().begin());
}
pExtractor()->select(this);
setPartonBinInfo();
lastP1P2(ll);
lastS(sqr(maxEnergy())/exp(lastP1() + lastP2()));
meMomenta().resize(mePartonData().size());
matrixElement()->setXComb(this);
PPair partons;
if ( !matrixElement()->haveX1X2() ) {
if ( !pExtractor()->generateL(partonBinInstances(),
r, r + nr - partonDims.second) ) {
lastCrossSection(ZERO);
return ZERO;
}
partons = make_pair(partonBinInstances().first->parton(),
partonBinInstances().second->parton());
lastSHat(lastS()/exp(partonBinInstances().first->l() +
partonBinInstances().second->l()));
meMomenta()[0] = partons.first->momentum();
meMomenta()[1] = partons.second->momentum();
} else {
if ( !matrixElement()->generateKinematics(r + partonDims.first) ) {
lastCrossSection(ZERO);
return ZERO;
}
lastSHat((meMomenta()[0]+meMomenta()[1]).m2());
matrixElement()->setKinematics();
lastScale(matrixElement()->scale());
partons.first = mePartonData()[0]->produceParticle(meMomenta()[0]);
partons.second = mePartonData()[1]->produceParticle(meMomenta()[1]);
Direction<0> dir(true);
partonBinInstances().first =
new_ptr(PartonBinInstance(lastParticles().first,partons.first,
partonBins().first,lastScale()));
dir.reverse();
partonBinInstances().second =
new_ptr(PartonBinInstance(lastParticles().second,partons.second,
partonBins().second,lastScale()));
}
lastPartons(partons);
if ( lastSHat() < cuts()->sHatMin() ) {
lastCrossSection(ZERO);
return ZERO;
}
lastY(0.5*(partonBinInstances().second->l() -
partonBinInstances().first->l()));
if ( !cuts()->initSubProcess(lastSHat(), lastY(), mirror()) ) {
lastCrossSection(ZERO);
return ZERO;
}
if ( mirror() ) swap(meMomenta()[0], meMomenta()[1]);
if ( matrixElement()->wantCMS() &&
!matrixElement()->haveX1X2() )
SimplePhaseSpace::CMS(meMomenta()[0], meMomenta()[1], lastSHat());
Energy summ = ZERO;
if ( meMomenta().size() == 3 ) {
if ( !matrixElement()->haveX1X2() )
meMomenta()[2] = Lorentz5Momentum(sqrt(lastSHat()));
} else {
for ( int i = 2, N = meMomenta().size(); i < N; ++i ) {
if ( !matrixElement()->haveX1X2() )
meMomenta()[i] = Lorentz5Momentum(mePartonData()[i]->mass());
summ += mePartonData()[i]->massMin();
}
if ( sqr(summ) >= lastSHat() ) {
lastCrossSection(ZERO);
return ZERO;
}
}
matrixElement()->setXComb(this);
if ( !matrixElement()->haveX1X2() )
lastScale(max(lastSHat()/4.0, cuts()->scaleMin()));
lastSHat(pExtractor()->generateSHat(lastS(), partonBinInstances(),
r, r + nr - partonDims.second,
matrixElement()->haveX1X2()));
if ( !cuts()->sHat(lastSHat()) ) {
lastCrossSection(ZERO);
return ZERO;
}
r += partonDims.first;
lastX1X2(make_pair(lastPartons().first->momentum().plus()/
lastParticles().first->momentum().plus(),
lastPartons().second->momentum().minus()/
lastParticles().second->momentum().minus()));
if ( !cuts()->x1(lastX1()) || !cuts()->x2(lastX2()) ) {
lastCrossSection(ZERO);
return ZERO;
}
lastY((lastPartons().first->momentum() +
lastPartons().second->momentum()).rapidity());
if ( !cuts()->yHat(lastY()) ) {
lastCrossSection(ZERO);
return ZERO;
}
if ( !cuts()->initSubProcess(lastSHat(), lastY(), mirror()) ) {
lastCrossSection(ZERO);
return ZERO;
}
meMomenta()[0] = lastPartons().first->momentum();
meMomenta()[1] = lastPartons().second->momentum();
if ( mirror() ) swap(meMomenta()[0], meMomenta()[1]);
if ( matrixElement()->wantCMS() &&
!matrixElement()->haveX1X2() )
SimplePhaseSpace::CMS(meMomenta()[0], meMomenta()[1], lastSHat());
if ( meMomenta().size() == 3 ) {
if ( !matrixElement()->haveX1X2() )
meMomenta()[2] = Lorentz5Momentum(sqrt(lastSHat()));
} else {
if ( sqr(summ) >= lastSHat() ) {
lastCrossSection(ZERO);
return ZERO;
}
}
matrixElement()->setXComb(this);
if ( !matrixElement()->haveX1X2() ) {
if ( !matrixElement()->generateKinematics(r) ) {
lastCrossSection(ZERO);
return ZERO;
}
}
lastScale(matrixElement()->scale());
if ( !cuts()->scale(lastScale()) ) {
lastCrossSection(ZERO);
return ZERO;
}
pair<bool,bool> evalPDFS =
make_pair(matrixElement()->havePDFWeight1(),
matrixElement()->havePDFWeight2());
if ( mirror() )
swap(evalPDFS.first,evalPDFS.second);
lastPDFWeight(pExtractor()->fullFn(partonBinInstances(), lastScale(),
evalPDFS));
if ( lastPDFWeight() == 0.0 ) {
lastCrossSection(ZERO);
return ZERO;
}
matrixElement()->setKinematics();
CrossSection xsec = matrixElement()->dSigHatDR() * lastPDFWeight();
bool noHeadPass = !willPassCuts() || xsec == ZERO;
if ( noHeadPass )
lastCrossSection(ZERO);
lastAlphaS(matrixElement()->alphaS());
lastAlphaEM(matrixElement()->alphaEM());
subProcess(SubProPtr());
if ( CKKWHandler() && matrixElement()->maxMultCKKW() > 0 &&
matrixElement()->maxMultCKKW() > matrixElement()->minMultCKKW() ) {
newSubProcess(theMEGroup->subProcessGroups());
CKKWHandler()->setXComb(this);
xsec *= CKKWHandler()->reweightCKKW(matrixElement()->minMultCKKW(),
matrixElement()->maxMultCKKW());
}
if ( matrixElement()->reweighted() ) {
newSubProcess(theMEGroup->subProcessGroups());
xsec *= matrixElement()->reWeight() * matrixElement()->preWeight();
}
lastHeadCrossSection(xsec);
CrossSection depxsec = ZERO;
if ( !theMEGroup->mcSumDependent() ) {
for ( StdDepXCVector::const_iterator dep = theDependent.begin();
dep != theDependent.end(); ++dep ) {
if ( !(*dep) )
continue;
if ( !(**dep).matrixElement()->apply() )
continue;
if ( noHeadPass && (**dep).matrixElement()->headCuts() )
continue;
depxsec += (**dep).dSigDR();
}
} else {
if ( !noHeadPass || !theMEGroup->lastDependentXComb()->matrixElement()->headCuts() )
depxsec = theDependent.size()*theMEGroup->lastDependentXComb()->dSigDR();
}
if ( xsec != ZERO ) {
double rw = 1.0 + depxsec/xsec;
xsec *= rw;
} else {
xsec = depxsec;
}
lastCrossSection(xsec);
if ( xsec != ZERO )
theMEGroup->lastEventStatistics();
- if ( !theMEGroup->subProcessGroups() ) {
- lastHeadCrossSection(xsec);
- }
-
return xsec;
}
void StdXCombGroup::newSubProcess(bool) {
// subprocess selection goes here
// if me group returns an associated
// selector
StandardXComb::newSubProcess(theMEGroup->subProcessGroups());
if ( !theMEGroup->subProcessGroups() )
return;
subProcess()->groupWeight(lastHeadCrossSection()/lastCrossSection());
Ptr<SubProcessGroup>::tptr group =
dynamic_ptr_cast<Ptr<SubProcessGroup>::tptr>(subProcess());
assert(group);
for ( StdDepXCVector::iterator dep = theDependent.begin();
dep != theDependent.end(); ++dep ) {
if ( !(*dep) )
continue;
if ( (**dep).lastCrossSection() == ZERO ||
!(**dep).matrixElement()->apply() )
continue;
tSubProPtr ds;
try {
ds = (**dep).construct();
} catch(Veto&) {
continue;
}
if ( ds )
group->add(ds);
}
}
tSubProPtr StdXCombGroup::construct() {
if ( !theMEGroup->subProcessGroups() )
return StandardXComb::construct();
matrixElement()->setXComb(this);
setPartonBinInfo();
matrixElement()->setKinematics();
newSubProcess(true);
TmpTransform<tSubProPtr>
tmp(subProcess(), Utilities::getBoostToCM(subProcess()->incoming()));
if ( !cuts()->passCuts(*subProcess()) ) {
subProcess()->groupWeight(0.0);
}
return subProcess();
}
void StdXCombGroup::Init() {}
void StdXCombGroup::persistentOutput(PersistentOStream & os) const {
os << theDependent << theMEGroup << ounit(theLastHeadCrossSection,nanobarn);
}
void StdXCombGroup::persistentInput(PersistentIStream & is, int) {
is >> theDependent >> theMEGroup >> iunit(theLastHeadCrossSection,nanobarn);
}
ClassDescription<StdXCombGroup> StdXCombGroup::initStdXCombGroup;
diff --git a/Handlers/StdXCombGroup.h b/Handlers/StdXCombGroup.h
--- a/Handlers/StdXCombGroup.h
+++ b/Handlers/StdXCombGroup.h
@@ -1,187 +1,187 @@
// -*- C++ -*-
//
// StdXCombGroup.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2007 Leif Lonnblad
// Copyright (C) 2009-2010 Simon Platzer
//
// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_StdXCombGroup_H
#define ThePEG_StdXCombGroup_H
// This is the declaration of the StdXCombGroup class.
#include "StandardXComb.h"
#include "StdDependentXComb.h"
#include "StdXCombGroup.fh"
#include "ThePEG/MatrixElement/MEGroup.fh"
namespace ThePEG {
/**
* The StdXCombGroup class represents a 'head' XComb object
* in association with a group of dependent XComb objects.
*
* @see StdDependendXComb
* @see MEGroup
*/
class StdXCombGroup: public StandardXComb {
/** MEBase needs to be a friend. */
friend class MEBase;
public:
/** @name Standard constructors and destructors. */
//@{
/**
* Standard constructor.
*/
StdXCombGroup(Energy newMaxEnergy, const cPDPair & inc,
tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler,
tPExtrPtr newExtractor, tCascHdlPtr newCKKW,
const PBPair & newPartonBins, tCutsPtr newCuts, tMEGroupPtr newME,
const DiagramVector & newDiagrams, bool mir);
/**
* Default constructor.
*/
StdXCombGroup();
/**
* Destructor.
*/
virtual ~StdXCombGroup();
public:
/**
* Generate a phase space point from a vector \a r of \a nr numbers
* in the interval ]0,1[ and return the corresponding differential
* cross section.
*/
virtual CrossSection dSigDR(const pair<double,double> ll, int nr, const double * r);
/**
+ * Return the cross section calculated from the head matrix element
+ */
+ CrossSection lastHeadCrossSection() const { return theLastHeadCrossSection; }
+
+ /**
* Visit the dependent XComb objects
*/
const StdDepXCVector& dependent() const { return theDependent; }
/**
* Return the matrix element group steered by this
* XComb group.
*/
tcMEGroupPtr meGroup() const { return theMEGroup; }
/**
* Initialize this XComb group
*/
void build(const PartonPairVec& allPBins);
/**
* Construct a sub-process object from the information available.
*/
virtual tSubProPtr construct();
protected:
/**
* Construct the corresponding SubProcess object if it hasn't been
* done before.
*/
virtual void newSubProcess(bool);
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);
//@}
/**
* Standard Init function used to initialize the interface.
*/
static void Init();
protected:
/**
* Set the cross section calculated from the head matrix element
*/
void lastHeadCrossSection(CrossSection xs) { theLastHeadCrossSection = xs; }
- /**
- * Return the cross section calculated from the head matrix element
- */
- CrossSection lastHeadCrossSection() const { return theLastHeadCrossSection; }
-
private:
/**
* The MEGroup object
*/
MEGroupPtr theMEGroup;
/**
* The dependent XComb objects
*/
StdDepXCVector theDependent;
/**
* The cross section calculated from the head matrix element
*/
CrossSection theLastHeadCrossSection;
private:
/**
* Describe a concrete class with persistent data.
*/
static ClassDescription<StdXCombGroup> initStdXCombGroup;
/**
* Private and non-existent assignment operator.
*/
StdXCombGroup & operator=(const StdXCombGroup &);
};
/** @cond TRAITSPECIALIZATIONS */
/**
* This template specialization informs ThePEG about the base class of
* StdXCombGroup.
*/
template <>
struct BaseClassTrait<StdXCombGroup,1> {
/** Typedef of the base class of StdXCombGroup. */
typedef StandardXComb NthBase;
};
/**
* This template specialization informs ThePEG about the name of the
* StdXCombGroup class.
*/
template <>
struct ClassTraits<StdXCombGroup>:
public ClassTraitsBase<StdXCombGroup> {
/** Return the class name. */
static string className() { return "ThePEG::StdXCombGroup"; }
};
/** @endcond */
}
#endif /* ThePEG_StdXCombGroup_H */
diff --git a/MatrixElement/ColourLines.cc b/MatrixElement/ColourLines.cc
--- a/MatrixElement/ColourLines.cc
+++ b/MatrixElement/ColourLines.cc
@@ -1,138 +1,143 @@
// -*- C++ -*-
//
// ColourLines.cc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2011 Leif Lonnblad
//
// ThePEG 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 ColourLines class.
//
#include "ColourLines.h"
#include "ColourLines.xh"
#include "ThePEG/EventRecord/ColourLine.h"
#include "ThePEG/EventRecord/MultiColour.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Utilities/StringUtils.h"
using namespace ThePEG;
ColourLines::ColourLines(string s) {
+ reset(s);
+}
+
+void ColourLines::reset(string s) {
+ theLines.clear();
while ( true ) {
string line = StringUtils::car(s, ",");
line = StringUtils::stripws(line);
Line l;
while (line!="") {
string loc = StringUtils::car(line);
string first = StringUtils::car(loc,":");
string second = StringUtils::cdr(loc,":");
if(second!="") {
int i;
istringstream is(first);
is >> i;
int j;
istringstream is2(second);
is2 >> j;
l.push_back(make_pair(i,j));
}
else {
int i;
istringstream is(first);
is >> i;
l.push_back(make_pair(i,0));
}
line = StringUtils::cdr(line);
};
if ( l.empty() ) return;
theLines.push_back(l);
s = StringUtils::cdr(s, ",");
}
}
void ColourLines::connect(const tPVector & partons) const {
VertexVector sinks;
VertexVector sources;
long np = partons.size();
// Create each line and connect the specified partons to them. Save
// all lines coming from a source or ending in a sink.
for ( LineVector::size_type il = 0; il < theLines.size(); ++il ) {
const Line & line = theLines[il];
ColinePtr cline = new_ptr(ColourLine());
for ( Line::size_type i = 0; i < line.size(); ++i ) {
if ( line[i].first > np ) {
// this is a colour source.
int is = line[i].first - np;
sources.resize(is);
sources[is - 1].push_back(cline);
} else if ( -line[i].first > np ) {
// this is a colour sink.
int is = -line[i].first - np;
sources.resize(is);
sources[is - 1].push_back(cline);
} else if ( line[i].first > 0 ) {
// This is a coloured particle.
if ( !partons[line[i].first - 1]->hasColour() )
throw ColourGeometryException(partons, line);
if(line[i].second==0) {
cline->addColoured(partons[line[i].first - 1]);
}
else {
Ptr<MultiColour>::pointer colour =
dynamic_ptr_cast<Ptr<MultiColour>::pointer>(partons[line[i].first - 1]->colourInfo());
assert(colour);
colour->colourLine(cline,line[i].second);
}
} else {
if ( !partons[-line[i].first - 1]->hasAntiColour() )
throw ColourGeometryException(partons, line);
if(line[i].second==0) {
cline->addAntiColoured(partons[-line[i].first - 1]);
}
else {
Ptr<MultiColour>::pointer colour =
dynamic_ptr_cast<Ptr<MultiColour>::pointer>(partons[-line[i].first - 1]->colourInfo());
assert(colour);
colour->antiColourLine(cline,line[i].second);
}
}
}
}
// Now connect up all lines steming from sources.
for ( VertexVector::size_type i = 0; i < sources.size(); ++i ) {
if ( sources[i].empty() ) continue;
if ( sources[i].size() != 3 ) throw ColourGeometryException(partons,
vector<pair<int,int> >() );
sources[i][0]->setSourceNeighbours(sources[i][1], sources[i][2]);
}
// Now connect up all lines ending in sinks.
for ( VertexVector::size_type i = 0; i < sinks.size(); ++i ) {
if ( sinks[i].empty() ) continue;
if ( sinks[i].size() != 3 ) throw ColourGeometryException(partons,
vector<pair<int,int> >());
sinks[i][0]->setSinkNeighbours(sinks[i][1], sinks[i][2]);
}
}
ColourGeometryException::
ColourGeometryException(const tPVector & p, const vector<pair<int,int> > & c) {
if ( c.empty() )
theMessage << "The number of colour lines steming from one colour source "
<< "or ending in one colour sink was not equal to three.\n";
else {
theMessage << "Cannot connect the following partons:\n";
for ( unsigned i = 0; i < p.size(); ++i )
theMessage << " " << p[i]->PDGName();
theMessage << "\n to the following colour line:\n";
for ( unsigned i = 0; i < c.size(); ++i ) theMessage << " (" << c[i].first << ","
<< c[i].second << ") ";
theMessage << endl;
}
severity(maybeabort);
}
diff --git a/MatrixElement/ColourLines.h b/MatrixElement/ColourLines.h
--- a/MatrixElement/ColourLines.h
+++ b/MatrixElement/ColourLines.h
@@ -1,78 +1,89 @@
// -*- C++ -*-
//
// ColourLines.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2011 Leif Lonnblad
//
// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_ColourLines_H
#define ThePEG_ColourLines_H
// This is the declaration of the ColourLines class.
#include "ThePEG/Config/ThePEG.h"
namespace ThePEG {
/**
* The ColourLines class defines the colour flow in a SubProcess. It
* defines a number of colour lines and specifies which particles are
* connected to them.
*
*/
class ColourLines: public Base {
public:
/** A single colour line */
typedef vector<pair<int,int> > Line;
/** A vector of colour lines. */
typedef vector<Line> LineVector;
/** A vector of <code>ColourLine</code>. */
typedef vector<ColinePtr> Vertex;
/** A vector of vertices. */
typedef vector<Vertex> VertexVector;
public:
/** @name Standard constructors and destructors. */
//@{
/**
* Default constructor.
*/
ColourLines() {}
/**
* The standard constructor. The string \a s should contain a
* comma-separated sequence of integers. Each sequence of numbers
* indicates a colour line and the integer represents a parton
* connected to it. If the integer is negative, it means that the
* line is the corresponding partons anti-colour. Note that the
* partons are numbered from 1: The first spacelike particle is 1, the second
* is 2 and the internal time-like and outgoing are numbered after all
* the spacelike particles.
*/
ColourLines(string s);
//@}
+ /**
+ * Reset this ColourLines object. The string \a s should contain a
+ * comma-separated sequence of integers. Each sequence of numbers
+ * indicates a colour line and the integer represents a parton
+ * connected to it. If the integer is negative, it means that the
+ * line is the corresponding partons anti-colour. Note that the
+ * partons are numbered from 1: The first incoming is 1, the second
+ * is 2 and the internal and outgoing are numbered 3 and upwards.
+ */
+ void reset(string s);
+
public:
/**
* Create the corresponding <code>ColourLine</code>s and connect the
* given \a partons. The partons are assumed to be in the same order
* as the numbers specified in the constructor.
*/
void connect(const tPVector & partons) const;
private:
/**
* The vector of colour lines.
*/
LineVector theLines;
};
}
#endif /* ThePEG_ColourLines_H */

File Metadata

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

Event Timeline