diff --git a/EventRecord/Particle.cc b/EventRecord/Particle.cc
--- a/EventRecord/Particle.cc
+++ b/EventRecord/Particle.cc
@@ -1,521 +1,524 @@
// -*- C++ -*-
//
// Particle.cc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2017 Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined member functions of
// the Particle class.
//
#include "Particle.h"
#include "ThePEG/EventRecord/Step.h"
#include "ThePEG/EventRecord/Event.h"
#include "ThePEG/EventRecord/ColourLine.h"
#include "ThePEG/Utilities/Rebinder.h"
#include "ThePEG/Config/algorithm.h"
#include "ThePEG/EventRecord/ParticleTraits.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/DecayMode.h"
#include <iostream>
#include <iomanip>
#include <cctype>
#ifdef ThePEG_TEMPLATES_IN_CC_FILE
#include "Particle.tcc"
#endif
using namespace ThePEG;
Particle::ParticleRep::ParticleRep(const ParticleRep & p)
: theParents(p.theParents), theChildren(p.theChildren),
thePrevious(p.thePrevious), theNext(p.theNext),
theBirthStep(p.theBirthStep), theVertex(p.theVertex),
theLifeLength(p.theLifeLength), theScale(p.theScale),
theVetoScale(p.theVetoScale), theNumber(p.theNumber),
theExtraInfo(p.theExtraInfo.size()) {
if ( p.theColourInfo )
theColourInfo = dynamic_ptr_cast<CBPtr>(p.theColourInfo->clone());
if ( p.theSpinInfo )
theSpinInfo = dynamic_ptr_cast<SpinPtr>(p.theSpinInfo->clone());
for ( int i = 0, N = p.theExtraInfo.size(); i < N; ++i )
theExtraInfo[i] = p.theExtraInfo[i]->clone();
}
Particle::Particle(const Particle & p)
: Base(p), theData(p.theData), theMomentum(p.theMomentum), theRep(p.theRep) {
if ( p.theRep ) {
theRep = new ParticleRep(*p.theRep);
theRep->theParents.clear();
}
}
Particle::~Particle() {
if ( theRep ) {
if ( colourLine() ) colourLine()->removeColoured(this);
if ( antiColourLine() ) antiColourLine()->removeAntiColoured(this);
delete theRep;
}
theRep = 0;
theData = cEventPDPtr();
}
void Particle::initFull() {
if ( theRep ) return;
theRep = new ParticleRep;
Energy width = data().generateWidth(mass());
if ( width > ZERO ) {
Time lifetime = data().generateLifeTime(mass(), width);
theRep->theLifeLength.setTau(lifetime);
theRep->theLifeLength.
setVect((momentum().vect()*(lifetime /
max(mass(), Constants::epsilon*GeV))));
theRep->theLifeLength.rescaleEnergy();
}
}
PPtr Particle::clone() const {
return ptr_new<PPtr>(*this);
}
void Particle::rebind(const EventTranslationMap & trans) {
for ( ParticleVector::iterator pit = rep().theChildren.begin();
pit != rep().theChildren.end(); ++pit ) *pit = trans.translate(*pit);
for ( tParticleVector::iterator pit = rep().theParents.begin();
pit != rep().theParents.end(); ++pit ) *pit = trans.translate(*pit);
rep().thePrevious = trans.translate(rep().thePrevious);
rep().theNext = trans.translate(rep().theNext);
if ( hasColourInfo() ) colourInfo()->rebind(trans);
if ( spinInfo() ) spinInfo()->rebind(trans);
rep().theBirthStep = trans.translate(rep().theBirthStep);
for ( EIVector::const_iterator ie = rep().theExtraInfo.begin();
ie != rep().theExtraInfo.end(); ++ie ) (**ie).rebind(trans);
}
tParticleSet Particle::siblings() const {
tParticleSet theSiblings;
for ( tParticleVector::const_iterator pit = parents().begin();
pit != parents().end(); ++pit )
theSiblings.insert((*pit)->children().begin(), (*pit)->children().end());
theSiblings.erase(const_cast<Particle *>(this));
return theSiblings;
}
void Particle::colourNeighbour(tPPtr p, bool anti) {
tColinePtr line = colourLine(!anti);
if ( !line ) line = ColourLine::create(this, !anti);
line->addColoured(p, anti);
}
void Particle::outgoingColour(tPPtr p, bool anti) {
tColinePtr line = colourLine(anti);
if ( !line ) line = ColourLine::create(this, anti);
line->addColoured(p, anti);
}
tPPtr Particle::incomingColour(bool anti) const {
if ( !hasColourInfo() ) return tPPtr();
tColinePtr line = colourLine(anti);
if ( !line ) return tPPtr();
for ( int i = 0, N = parents().size(); i < N; ++i )
if ( parents()[i]->hasColourLine(line, anti) ) return parents()[i];
return tPPtr();
}
tPPtr Particle::outgoingColour(bool anti) const {
if ( !hasColourInfo() ) return tPPtr();
tColinePtr line = colourLine(anti);
if ( !line ) return tPPtr();
for ( int i = 0, N = children().size(); i < N; ++i )
if ( children()[i]->hasColourLine(line, anti) ) return children()[i];
return tPPtr();
}
LorentzPoint Particle::labVertex() const {
LorentzPoint r(rep().theBirthStep && rep().theBirthStep->collision()?
vertex() + rep().theBirthStep->collision()->vertex():
vertex());
return r;
}
void Particle::setLabVertex(const LorentzPoint & p) {
rep().theVertex = ( rep().theBirthStep && rep().theBirthStep->collision()?
p - rep().theBirthStep->collision()->vertex() : p );
}
void Particle::transform(const LorentzRotation & r) {
if ( hasRep() && spinInfo() ) spinInfo()->transform(momentum(), r);
theMomentum.transform(r);
if ( !hasRep() ) return;
rep().theVertex.transform(r);
rep().theLifeLength.transform(r);
}
void Particle::deepTransform(const LorentzRotation & r) {
transform(r);
if ( !theRep ) return;
for ( int i = 0, N = children().size(); i < N; ++i )
rep().theChildren[i]->deepTransform(r);
if ( rep().theNext ) rep().theNext->deepTransform(r);
}
void Particle::rotateX(double a) {
LorentzRotation r;
r.rotateX(a);
transform(r);
}
void Particle::deepRotateX(double a) {
LorentzRotation r;
r.rotateX(a);
deepTransform(r);
}
void Particle::rotateY(double a) {
LorentzRotation r;
r.rotateY(a);
transform(r);
}
void Particle::deepRotateY(double a) {
LorentzRotation r;
r.rotateY(a);
deepTransform(r);
}
void Particle::rotateZ(double a) {
LorentzRotation r;
r.rotateZ(a);
transform(r);
}
void Particle::deepRotateZ(double a) {
LorentzRotation r;
r.rotateZ(a);
deepTransform(r);
}
void Particle::rotate(double a, const Axis & axis) {
LorentzRotation r;
r.rotate(a, axis);
transform(r);
}
void Particle::deepRotate(double a, const Axis & axis) {
LorentzRotation r;
r.rotate(a, axis);
deepTransform(r);
}
string Particle::outputFormat =
"%n3%s10 %i7 %p[,]0 %c(,) %^^0%vv0 %>>0%<>0 %l{,}0\n"
" %x10.3%y10.3%z10.3%e10.3%m10.3\n";
int getNumber(string::const_iterator & pos, int def) {
if ( !isdigit(*pos) ) return def;
def = *pos++ - '0';
while ( isdigit(*pos) ) def = 10*def + *pos++ - '0';
return def;
}
void writePrecision(ostream & os, string::const_iterator & pos,
int defw, int defp, double x) {
defw = getNumber(pos, defw);
if ( *pos == '.' ) defp = getNumber(++pos, defp);
int oldp = os.precision();
os << setprecision(defp) << setw(defw) << x << setprecision(oldp);
}
void writeStringAdjusted(ostream & os, bool left, int w, string str) {
while ( !left && w-- > int(str.size()) ) os << ' ';
os << str;
while ( left && w-- > int(str.size()) ) os << ' ';
}
template <typename Container>
void writeParticleRanges(ostream & os, const Container & co, char sep, int w) {
set<int> cnum;
for ( typename Container::const_iterator it = co.begin();
it != co.end(); ++it) cnum.insert((**it).number());
bool elipsis = false;
int last = -10;
for ( set<int>::iterator it = cnum.begin(); it != cnum.end(); ++it) {
int n = *it;
int next = 0;
set<int>::iterator itn = it;
if ( ++itn != cnum.end() ) next = *itn;
bool writeit = true;
bool writesep = false;
if ( elipsis && ( n != last + 1 || n != next - 1 ) )
elipsis = false;
else if ( !elipsis && n == last + 1 && n == next -1 ) {
os << "..";
elipsis = true;
writeit = false;
}
else if ( elipsis && n == last + 1 && n == next -1 )
writeit = false;
else if ( it != cnum.begin() )
writesep = true;
if ( writeit ) {
if ( writesep ) os << sep;
os << setw(w) << n;
}
last = n;
}
}
ostream & ThePEG::operator<<(ostream & os, const Particle & p) {
return p.print(os, p.birthStep());
}
ostream & Particle::print(ostream & os, tcStepPtr step) const {
if ( !step ) step = birthStep();
tCollPtr coll = step? step->collision(): tCollPtr();
tEventPtr event = coll? coll->event(): tEventPtr();
string::const_iterator pos = Particle::outputFormat.begin();
ios::fmtflags saveflags = os.setf(ios::fixed, ios::floatfield);
while ( pos != Particle::outputFormat.end() ) {
if ( *pos == '%' && ++pos != Particle::outputFormat.end() ) {
bool left = false;
if ( *pos == '-' ) {
left = true;
os.setf(ios::left, ios::adjustfield);
++pos;
} else {
os.setf(ios::right, ios::adjustfield);
}
char mark;
char open;
char close;
char sep;
int w;
string str;
string fill;
if ( pos == Particle::outputFormat.end() ) break;
bool fullColour = false;
switch ( *pos ) {
case 'n':
os << setw(getNumber(++pos, 3)) << number();
break;
case 'i':
os << setw(getNumber(++pos, 8)) << id();
break;
case 's':
writeStringAdjusted(os, left, getNumber(++pos, 8), PDGName());
break;
case 'x':
writePrecision(os, ++pos, 10, 3, momentum().x()/GeV);
break;
case 'y':
writePrecision(os, ++pos, 10, 3, momentum().y()/GeV);
break;
case 'z':
writePrecision(os, ++pos, 10, 3, momentum().z()/GeV);
break;
case 'e':
writePrecision(os, ++pos, 10, 3, momentum().e()/GeV);
break;
case 'm':
writePrecision(os, ++pos, 10, 3, momentum().mass()/GeV);
break;
case 'P':
fullColour = true;
+ [[fallthrough]];
case 'p':
open = *++pos;
sep = *++pos;
close = *++pos;
w = getNumber(++pos, 0);
if ( parents().empty() ) break;
if ( open ) os << open;
writeParticleRanges(os, parents(), sep, w);
if ( fullColour && hasColourInfo() &&
( incomingColour() || incomingAntiColour() ) ) {
if ( close ) os << open;
if ( incomingColour() )
os << "+" << incomingColour()->number();
if ( incomingAntiColour() )
os << "-" << incomingAntiColour()->number();
if ( close ) os << close;
}
if ( close ) os << close;
break;
case 'l':
open = *++pos;
sep = *++pos;
close = *++pos;
w = getNumber(++pos, 0);
if ( hasColourInfo() &&
( colourLine() || antiColourLine() ) && event) {
if ( open ) os << open;
vector<tcColinePtr> clines = colourInfo()->colourLines();
for ( int i = 0, N = clines.size(); i < N; ++i ) {
if ( i > 0 && sep ) os << sep;
clines[i]->write(os, event, false);
}
vector<tcColinePtr> aclines = colourInfo()->antiColourLines();
for ( int i = 0, N = aclines.size(); i < N; ++i ) {
if ( ( i > 0 || clines.size() ) && sep ) os << sep;
aclines[i]->write(os, event, true);
}
if ( close ) os << close;
}
break;
case 'C':
fullColour = true;
+ [[fallthrough]];
case 'c':
open = *++pos;
sep = *++pos;
close = *++pos;
w = getNumber(++pos, 0);
if ( children().empty() ) break;
if ( open ) os << open;
writeParticleRanges(os, children(), sep, w);
if ( fullColour && hasColourInfo() &&
( outgoingColour() || outgoingAntiColour() ) ) {
if ( close ) os << open;
if ( outgoingColour() )
os << "+" << outgoingColour()->number();
if ( outgoingAntiColour() )
os << "-" << outgoingAntiColour()->number();
if ( close ) os << close;
}
if ( close ) os << close;
break;
case '>':
mark = *++pos;
w = getNumber(++pos, 0);
if ( hasColourInfo() && step && step->colourNeighbour(this) ) {
os << setw(w-1) << step->colourNeighbour(this)->number() << mark;
}
break;
case '<':
mark = *++pos;
w = getNumber(++pos, 0);
if ( hasColourInfo() && step && step->antiColourNeighbour(this) ) {
int n = step->antiColourNeighbour(this)->number();
ostringstream oss;
oss << mark << n;
writeStringAdjusted(os, left, w, oss.str());
}
break;
case 'v':
mark = *++pos;
w = getNumber(++pos, 0);
if ( next() ) {
if ( left && mark ) os << mark;
os << setw(w) << next()->number();
if ( !left && mark ) os << mark;
}
break;
case '^':
mark = *++pos;
w = getNumber(++pos, 0);
if ( previous() ) {
if ( left && mark ) os << mark;
os << setw(w) << previous()->number();
if ( !left && mark ) os << mark;
}
break;
case 'd':
switch ( *++pos ) {
case 'x':
writePrecision(os, ++pos, 10, 3, lifeLength().x()/mm);
break;
case 'y':
writePrecision(os, ++pos, 10, 3, lifeLength().y()/mm);
break;
case 'z':
writePrecision(os, ++pos, 10, 3, lifeLength().z()/mm);
break;
case 't':
writePrecision(os, ++pos, 10, 3, lifeLength().e()/mm);
break;
case 'T':
writePrecision(os, ++pos, 10, 3, lifeLength().tau()/mm);
break;
}
break;
case 'V':
switch ( *++pos ) {
case 'x':
writePrecision(os, ++pos, 10, 3, vertex().x()/mm);
break;
case 'y':
writePrecision(os, ++pos, 10, 3, vertex().y()/mm);
break;
case 'z':
writePrecision(os, ++pos, 10, 3, vertex().z()/mm);
break;
case 't':
writePrecision(os, ++pos, 10, 3, vertex().e()/mm);
break;
}
+ break;
case 'L':
switch ( *++pos ) {
case 'x':
writePrecision(os, ++pos, 10, 3, labVertex().x()/mm);
break;
case 'y':
writePrecision(os, ++pos, 10, 3, labVertex().y()/mm);
break;
case 'z':
writePrecision(os, ++pos, 10, 3, labVertex().z()/mm);
break;
case 't':
writePrecision(os, ++pos, 10, 3, labVertex().e()/mm);
break;
}
break;
default:
os << *pos++;
}
} else {
if ( pos != Particle::outputFormat.end() ) os << *pos++;
}
}
os.flags(saveflags);
return os;
}
void Particle::debugme() const {
cerr << *this;
EventRecordBase::debugme();
}
void Particle::persistentOutput(PersistentOStream & os) const {
EventConfig::putParticleData(os, theData);
os << ounit(theMomentum, GeV) << bool( theRep != 0 );
if ( !theRep ) return;
os << rep().theParents << rep().theChildren
<< rep().thePrevious << rep().theNext << rep().theBirthStep
<< ounit(rep().theVertex, mm) << ounit(rep().theLifeLength, mm)
<< ounit(rep().theScale, GeV2) << ounit(rep().theVetoScale, GeV2)
<< rep().theNumber << rep().theDecayMode
<< rep().theColourInfo << rep().theSpinInfo << rep().theExtraInfo;
}
void Particle::persistentInput(PersistentIStream & is, int) {
bool readRep;
EventConfig::getParticleData(is, theData);
is >> iunit(theMomentum, GeV) >> readRep;
if ( !readRep ) return;
if ( !hasRep() ) theRep = new ParticleRep;
is >> rep().theParents >> rep().theChildren
>> rep().thePrevious >> rep().theNext >> rep().theBirthStep
>> iunit(rep().theVertex, mm) >> iunit(rep().theLifeLength, mm)
>> iunit(rep().theScale, GeV2) >> iunit(rep().theVetoScale, GeV2)
>> rep().theNumber >> rep().theDecayMode
>> rep().theColourInfo >> rep().theSpinInfo >> rep().theExtraInfo;
}
ClassDescription<Particle> Particle::initParticle;
void Particle::Init() {}
diff --git a/Helicity/Vertex/Tensor/SSTVertex.cc b/Helicity/Vertex/Tensor/SSTVertex.cc
--- a/Helicity/Vertex/Tensor/SSTVertex.cc
+++ b/Helicity/Vertex/Tensor/SSTVertex.cc
@@ -1,161 +1,161 @@
// -*- C++ -*-
//
// SSTVertex.cc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2003-2017 Peter Richardson, Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the SSTVertex class.
//
#include "SSTVertex.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
using namespace ThePEG;
using namespace Helicity;
// The following static variable is needed for the type
// description system in ThePEG.
DescribeAbstractNoPIOClass<SSTVertex,AbstractSSTVertex>
describeThePEGSSTVertex("ThePEG::SSTVertex", "libThePEG.so");
void SSTVertex::Init() {
static ClassDocumentation<SSTVertex> documentation
("The SSTVertex class is the implementation of the "
"helicity amplitude calculation for the scalar-scalar-tensor vertex"
", all such vertices should inherit from it");
}
// evaluate the vertex
Complex SSTVertex::evaluate(Energy2 q2, const ScalarWaveFunction & sca1,
const ScalarWaveFunction & sca2,
const TensorWaveFunction & ten) {
// obtain the coupling
setCoupling(q2,sca1.particle(),sca2.particle(),ten.particle());
Complex ii(0.,1.);
// evaluate the trace of the tensor
Complex trace = ten.tt()-ten.xx()-ten.yy()-ten.zz();
// dot product of the two momenta
Energy2 dot = sca1.momentum()*sca2.momentum();
Energy mass = sca1.particle()->mass();
// second term
complex<Energy2> second =
+2.*ten.tt()*sca1.e()*sca2.e() +2.*ten.xx()*sca1.px()*sca2.px()
+2.*ten.yy()*sca1.py()*sca2.py()+2.*ten.zz()*sca1.pz()*sca2.pz()
-(ten.tx()+ten.xt())*(sca1.e()*sca2.px() +sca1.px()*sca2.e())
-(ten.ty()+ten.yt())*(sca1.e()*sca2.py() +sca1.py()*sca2.e())
-(ten.tz()+ten.zt())*(sca1.e()*sca2.pz() +sca1.pz()*sca2.e())
+(ten.xy()+ten.yx())*(sca1.py()*sca2.px()+sca1.px()*sca2.py())
+(ten.xz()+ten.zx())*(sca1.pz()*sca2.px()+sca1.px()*sca2.pz())
+(ten.yz()+ten.zy())*(sca1.pz()*sca2.py()+sca1.py()*sca2.pz());
// return the answer
return -0.5*ii*norm()*UnitRemoval::InvE2*
(trace*(mass*mass-dot)+second)*sca1.wave()*sca2.wave();
}
// off-shell tensor
TensorWaveFunction SSTVertex::evaluate(Energy2 q2, int iopt, tcPDPtr out,
const ScalarWaveFunction & sca1,
const ScalarWaveFunction & sca2,
complex<Energy> mass, complex<Energy> width) {
// obtain the coupling
setCoupling(q2,sca1.particle(),sca2.particle(),out);
// array for the tensor
Complex ten[4][4];
// calculate the outgoing momentum
Lorentz5Momentum pout = sca1.momentum()+sca2.momentum();
// prefactor
if(mass.real() < ZERO) mass = out->mass();
complex<Energy2> mass2 = sqr(mass);
Energy2 p2 = pout.m2();
Complex fact = 0.25*norm()*sca1.wave()*sca2.wave()*
propagator(iopt,p2,out,mass,width);
// dot products we need
Energy2 dot12 = sca1.momentum()*sca2.momentum();
Energy2 dot1 = sca1.momentum()*pout;
Energy2 dot2 = pout*sca2.momentum();
// the vectors that we need for the tensor
LorentzPolarizationVectorE vec1,vec2;
Complex a,b;
Energy2 mphi2 = sqr(sca1.particle()->mass());
// massive case
if(mass.real()!=ZERO) {
Complex norm1=dot1/mass2;
Complex norm2=dot2/mass2;
- a = UnitRemoval::InvE2 * ((mphi2+dot12)*(Complex(2.*p2/mass2)-5.)
- +4.*(dot12-dot1*dot2/mass2))/3.;
+ a = Complex(UnitRemoval::InvE2 * ((mphi2+dot12)*(Complex(2.*p2/mass2)-5.)
+ +4.*(dot12-dot1*dot2/mass2)))/3.;
b = -(-(mphi2+dot12)*(2.+p2/mass2)+4.*(dot12-dot1*(dot2/mass2)))/3./mass2;
vec1 = LorentzPolarizationVectorE(sca1.momentum()) -
LorentzPolarizationVectorE(norm1 * pout);
vec2 = LorentzPolarizationVectorE(sca2.momentum()) -
LorentzPolarizationVectorE(norm2 * pout);
}
// massless case
else {
a = UnitRemoval::InvE2 * (-5.*(mphi2+dot12)+4.*dot12)/3.;
b = 0.;
vec1 = sca1.momentum();
vec2 = sca2.momentum();
}
// calculate the wavefunction
complex<Energy> vec1_tmp[4] = {vec1.x(), vec1.y(), vec1.z(), vec1.t()};
complex<Energy> vec2_tmp[4] = {vec2.x(), vec2.y(), vec2.z(), vec2.t()};
complex<Energy> pout_tmp[4] = {pout.x(), pout.y(), pout.z(), pout.t()};
for(int ix=0;ix<4;++ix) {
for(int iy=0;iy<4;++iy) {
complex<Energy2> temp = -2.*( vec1_tmp[ix]*vec2_tmp[iy] +
vec1_tmp[ix]*vec2_tmp[iy])
-b*pout_tmp[ix]*pout_tmp[iy];
ten[ix][iy]= UnitRemoval::InvE2 * temp;
}
}
ten[3][3]=ten[3][3]-a;
for(int ix=0;ix<3;++ix){ten[ix][ix]=ten[ix][ix]+a;}
// prefactor
for(int ix=0;ix<4;++ix){for(int iy=0;iy<4;++iy){ten[ix][iy]=fact*ten[ix][iy];}}
// return the wavefunction
return TensorWaveFunction(pout,out,
ten[0][0],ten[0][1],ten[0][2],ten[0][3],
ten[1][0],ten[1][1],ten[1][2],ten[1][3],
ten[2][0],ten[2][1],ten[2][2],ten[2][3],
ten[3][0],ten[3][1],ten[3][2],ten[3][3]);
}
// off-shell scalar
ScalarWaveFunction SSTVertex::evaluate(Energy2 q2,int iopt, tcPDPtr out,
const ScalarWaveFunction & sca,
const TensorWaveFunction & ten,
complex<Energy> mass, complex<Energy> width) {
// obtain the coupling
setCoupling(q2,sca.particle(),out,ten.particle());
// calculate the outgoing momentum
Lorentz5Momentum pout = sca.momentum()+ten.momentum();
// prefactors
if(mass.real() < ZERO) mass = out->mass();
complex<Energy2> mass2 = sqr(mass);
Energy2 p2 = pout.m2();
Complex fact = 0.5*norm()*sca.wave()*propagator(iopt,p2,out,mass,width);
// trace of the tensor
Complex trace1 = ten.tt()-ten.xx()-ten.yy()-ten.zz();
// dot product of the two momenta
Energy2 dot = sca.momentum()*pout;
// first term
complex<Energy2> trace = trace1*(mass2-dot);
// second term
complex<Energy2> second =
+2.*ten.tt()*sca.e()*pout.e() +2.*ten.xx()*sca.px()*pout.x()
+2.*ten.yy()*sca.py()*pout.y()+2.*ten.zz()*sca.pz()*pout.z()
-(ten.tx()+ten.xt())*( sca.e()*pout.x()+sca.px()*pout.e())
-(ten.ty()+ten.yt())*( sca.e()*pout.y()+sca.py()*pout.e())
-(ten.tz()+ten.zt())*( sca.e()*pout.z()+sca.pz()*pout.e())
+(ten.xy()+ten.yx())*(sca.py()*pout.x()+sca.px()*pout.y())
+(ten.xz()+ten.zx())*(sca.pz()*pout.x()+sca.px()*pout.z())
+(ten.yz()+ten.zy())*(sca.py()*pout.z()+sca.pz()*pout.y());
// put it all together
second = fact*(trace+second);
Complex result = second * UnitRemoval::InvE2;
return ScalarWaveFunction(pout,out,result);
}
diff --git a/Helicity/Vertex/Tensor/VVTVertex.cc b/Helicity/Vertex/Tensor/VVTVertex.cc
--- a/Helicity/Vertex/Tensor/VVTVertex.cc
+++ b/Helicity/Vertex/Tensor/VVTVertex.cc
@@ -1,276 +1,277 @@
// -*- C++ -*-
//
// VVTVertex.cc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 2003-2017 Peter Richardson, Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the VVTVertex class.
//
#include "VVTVertex.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
using namespace ThePEG;
using namespace Helicity;
// The following static variable is needed for the type
// description system in ThePEG.
DescribeAbstractNoPIOClass<VVTVertex,AbstractVVTVertex>
describeThePEGVVTVertex("ThePEG::VVTVertex", "libThePEG.so");
void VVTVertex::Init() {
static ClassDocumentation<VVTVertex> documentation
("The VVTVertex class is the implementation of"
"the vector-vector tensor vertices for helicity "
"amplitude calculations. All such vertices should inherit"
"from it.");
}
// function to evaluate the vertex
Complex VVTVertex::evaluate(Energy2 q2, const VectorWaveFunction & vec1,
const VectorWaveFunction & vec2,
const TensorWaveFunction & ten,
Energy vmass) {
// set the couplings
setCoupling(q2,vec1.particle(),vec2.particle(),ten.particle());
// mass of the vector
if(vmass<ZERO) vmass = vec1.particle()->mass();
// mass+k1.k2
Energy2 mdot = vec1.momentum()*vec2.momentum();
if(vmass!=ZERO) mdot += sqr(vmass);
// dot product of wavefunctions and momenta
Complex dotv1v2 = vec1.wave().dot(vec2.wave());
complex<Energy> dotk1v2 = vec1.momentum()*vec2.wave();
complex<Energy> dotk2v1 = vec1.wave()*vec2.momentum();
// dot product of wavefunctions and momenta with the tensor
LorentzPolarizationVector tv1 = ten.wave().postDot(vec1.wave());
LorentzPolarizationVector tv2 = ten.wave().postDot(vec2.wave());
LorentzPolarizationVectorE tk1 = ten.wave().postDot(vec1.momentum());
LorentzPolarizationVectorE tk2 = ten.wave().postDot(vec2.momentum());
Complex tenv1v2 =
vec1.wave ().dot(tv2) + vec2.wave ().dot(tv1);
complex<Energy> tenk1v2 =
vec1.momentum().dot(tv2) + vec2.wave ().dot(tk1);
complex<Energy> tenk2v1 =
vec2.momentum().dot(tv1) + vec1.wave ().dot(tk2);
complex<Energy2> tenk1k2 =
vec2.momentum().dot(tk1) + vec1.momentum().dot(tk2);
// trace of the tensor
Complex trace = ten.tt()-ten.xx()-ten.yy()-ten.zz();
// evaluate the vertex
return -0.5*Complex(0.,1.)*norm()*UnitRemoval::InvE2 *
(trace*(dotk1v2*dotk2v1-dotv1v2*mdot)
+mdot*tenv1v2-dotk2v1*tenk1v2
-dotk1v2*tenk2v1+dotv1v2*tenk1k2);
}
// evaluate an off-shell vector
VectorWaveFunction VVTVertex::evaluate(Energy2 q2, int iopt, tcPDPtr out,
const VectorWaveFunction & vec,
const TensorWaveFunction & ten,
complex<Energy> mass,
complex<Energy> width) {
// evaluate the couplings
setCoupling(q2,vec.particle(),out,ten.particle());
// outgoing momentum
Lorentz5Momentum pout = ten.momentum()+vec.momentum();
// normalisation factor
if(mass.real() < ZERO) mass = out->mass();
complex<Energy2> mass2 = sqr(mass);
Energy2 p2 = pout.m2();
Complex fact = -0.5*norm()*propagator(iopt,p2,out,mass,width);
// dot product of wavefunctions and momenta
complex<Energy2> dotk1k2 = vec.momentum()*pout;
complex<Energy> dotk2v1 = vec.wave() *pout;
// mass-k1.k2
complex<Energy2> mdot = -dotk1k2;
if(mass.real()!=ZERO) mdot += mass2;
// components of the tensor
Complex tentx = ten.tx()+ten.xt();
Complex tenty = ten.ty()+ten.yt();
Complex tentz = ten.tz()+ten.zt();
Complex tenxy = ten.xy()+ten.yx();
Complex tenxz = ten.xz()+ten.zx();
Complex tenyz = ten.yz()+ten.zy();
// dot product of momenta and polarization vectors with the tensor
complex<Energy> tenk2v1 =
2.*(+ten.tt()*vec.t()*pout.e() +ten.xx()*vec.x()*pout.x()
+ten.yy()*vec.y()*pout.y()+ten.zz()*vec.z()*pout.z())
-tentx*(vec.t()*pout.x()+vec.x()*pout.e())
-tenty*(vec.t()*pout.y()+vec.y()*pout.e())
-tentz*(vec.t()*pout.z()+vec.z()*pout.e())
+tenxy*(vec.x()*pout.y()+vec.y()*pout.x())
+tenxz*(vec.x()*pout.z()+vec.z()*pout.x())
+tenyz*(vec.y()*pout.z()+vec.z()*pout.y());
complex<Energy2> tenk1k2 =
2.*(+ten.tt()*vec.e()*pout.e() +ten.xx()*vec.px()*pout.x()
+ten.yy()*vec.py()*pout.y()+ten.zz()*vec.pz()*pout.z())
-tentx*(vec.e()*pout.x() +vec.px()*pout.e())
-tenty*(vec.e()*pout.y() +vec.py()*pout.e())
-tentz*(vec.e()*pout.z() +vec.pz()*pout.e())
+tenxy*(vec.px()*pout.y()+vec.py()*pout.x())
+tenxz*(vec.px()*pout.z()+vec.pz()*pout.x())
+tenyz*(vec.py()*pout.z()+vec.pz()*pout.y());
// trace of the tensor
Complex trace = ten.tt()-ten.xx()-ten.yy()-ten.zz();
// compute the vector
Complex vec1[4];
vec1[0] = UnitRemoval::InvE2*
(mdot*(+ tentx* vec.t()-2.*ten.xx()*vec.x()
- tenxy* vec.y()- tenxz* vec.z()-trace*vec.x())
+(tenk2v1-trace*dotk2v1)*vec.px()-tenk1k2*vec.x()
+dotk2v1*(+ tentx* vec.e() -2.*ten.xx()*vec.px()
- tenxy* vec.py()- tenxz* vec.pz()));
vec1[1] = UnitRemoval::InvE2*(
mdot *
(+ tenty* vec.t()
- tenxy* vec.x()
-2.*ten.yy()*vec.y()
- tenyz* vec.z()-trace*vec.y()
)
+(tenk2v1-trace*dotk2v1)*vec.py()
-tenk1k2*vec.y()
+dotk2v1*(+ tenty* vec.e()
- tenxy* vec.px()
-2.*ten.yy()*vec.py()
- tenyz* vec.pz()));
vec1[2] = UnitRemoval::InvE2*
(mdot*
(+ tentz* vec.t()- tenxz* vec.x()
- tenyz* vec.y()-2.*ten.zz()*vec.z()-trace*vec.z())
+(tenk2v1-trace*dotk2v1)*vec.pz()-tenk1k2*vec.z()
+dotk2v1*(+ tentz* vec.e() - tenxz *vec.px()
- tenyz* vec.py()-2.*ten.zz()*vec.pz()));
vec1[3] = UnitRemoval::InvE2*
(mdot*(+2.*ten.tt()*vec.t()- tentx* vec.x()
- tenty* vec.y()- tentz* vec.z()-trace*vec.t())
+(tenk2v1-trace*dotk2v1)*vec.e() -tenk1k2*vec.t()
+dotk2v1*(+2.*ten.tt()*vec.e() - tentx* vec.px()
- tenty* vec.py()- tentz* vec.pz()));
// now add the piece for massive bosons
if(mass.real()!=ZERO) {
// DGRELL unit problem?
Complex dot = tenk2v1 * UnitRemoval::InvE
- dotk1k2 * trace * UnitRemoval::InvE2;
vec1[0] -= dot*pout.x() * UnitRemoval::InvE;
vec1[1] -= dot*pout.y() * UnitRemoval::InvE;
vec1[2] -= dot*pout.z() * UnitRemoval::InvE;
vec1[3] -= dot*pout.e() * UnitRemoval::InvE;
}
// return the VectorWaveFunction
for(int ix=0;ix<4;++ix){vec1[ix]=vec1[ix]*fact;}
return VectorWaveFunction(pout,out,vec1[0],vec1[1],vec1[2],vec1[3]);
}
// off-shell tensor
TensorWaveFunction VVTVertex::evaluate(Energy2 q2, int iopt,tcPDPtr out,
const VectorWaveFunction & vec1,
const VectorWaveFunction & vec2,
Energy vmass,
complex<Energy> tmass,
complex<Energy> width) {
// coupling
setCoupling(q2,vec1.particle(),vec2.particle(),out);
// momenta of the outgoing tensor
// outgoing momentum
Lorentz5Momentum pout= vec1.momentum()+vec2.momentum();
// overall normalisation
if(tmass.real() < ZERO) tmass = out->mass();
complex<Energy2> tmass2 = sqr(tmass);
if(vmass<ZERO) vmass = vec1.particle()->mass();
Energy2 vmass2 = sqr(vmass);
Energy2 p2 = pout.m2();
Complex fact = 0.25*norm()*propagator(iopt,p2,out,tmass,width);
// dot products we need to construct the tensor
complex<Energy2> dotk1k2 = vec1.momentum()*vec2.momentum();
complex<Energy> dotv1k2 = vec1.wave()*vec2.momentum();
Complex dotv1v2 = vec1.wave().dot(vec2.wave());
complex<Energy> dotk1v2 = vec1.momentum()*vec2.wave();
complex<Energy2> dotkk1 = vec1.momentum()*pout;
complex<Energy2> dotkk2 = vec2.momentum()*pout;
complex<Energy> dotkv1 = vec1.wave()*pout;
complex<Energy> dotkv2 = vec2.wave()*pout;
// dot product ma^2+k1.k2
complex<Energy2> mdot = vmass2+dotk1k2;
// vectors to help construct the tensor
Complex vecv1[4],vecv2[4];
complex<Energy> veck1[4],veck2[4];
complex<InvEnergy2> tmass2inv(ZERO);
if(tmass.real()> ZERO) tmass2inv = 1./tmass2;
vecv1[0]=vec1.x() -pout.x()*dotkv1*tmass2inv;
vecv2[0]=vec2.x() -pout.x()*dotkv2*tmass2inv;
veck1[0]=vec1.px()-pout.x()*dotkk1*tmass2inv;
veck2[0]=vec2.px()-pout.x()*dotkk2*tmass2inv;
vecv1[1]=vec1.y() -pout.y()*dotkv1*tmass2inv;
vecv2[1]=vec2.y() -pout.y()*dotkv2*tmass2inv;
veck1[1]=vec1.py()-pout.y()*dotkk1*tmass2inv;
veck2[1]=vec2.py()-pout.y()*dotkk2*tmass2inv;
vecv1[2]=vec1.z() -pout.z()*dotkv1*tmass2inv;
vecv2[2]=vec2.z() -pout.z()*dotkv2*tmass2inv;
veck1[2]=vec1.pz()-pout.z()*dotkk1*tmass2inv;
veck2[2]=vec2.pz()-pout.z()*dotkk2*tmass2inv;
vecv1[3]=vec1.t() -pout.e()*dotkv1*tmass2inv;
vecv2[3]=vec2.t() -pout.e()*dotkv2*tmass2inv;
veck1[3]=vec1.e() -pout.e()*dotkk1*tmass2inv;
veck2[3]=vec2.e() -pout.e()*dotkk2*tmass2inv;
// coefficient of g(nu,mu)-k^muk^nu/m^2
+ Complex omp2 = 1.-Complex(p2*tmass2inv);
Complex coeff1 = UnitRemoval::InvE2 *
(
+4./3.*mdot*(-2.*dotv1v2+Complex((dotkv1*dotkv2+p2*dotv1v2)*tmass2inv))
+4./3.*(dotv1k2*(dotk1v2-dotkk1*dotkv2*tmass2inv)
+dotk1v2*(dotv1k2-dotkk2*dotkv1*tmass2inv)
-dotv1v2*(dotk1k2-dotkk1*dotkk2*tmass2inv)
- +(1.-p2*tmass2inv)*dotk1v2*dotv1k2)
+ +omp2*dotk1v2*dotv1k2)
);
// coefficient of g(nu,mu)
Complex coeff2 = UnitRemoval::InvE2 *
(
- 2.*mdot*(1.-p2*tmass2inv)*dotv1v2
- -2.*(1.-p2*tmass2inv)*dotk1v2*dotv1k2
+ 2.*mdot*omp2*dotv1v2
+ -2.*omp2*dotk1v2*dotv1k2
);
// construct the tensor
Complex ten[4][4];
const Energy pout_tmp[4] = {pout.x(), pout.y(), pout.z(), pout.e()};
for(int ix=0;ix<4;++ix) {
for(int iy=0;iy<4;++iy) {
complex<Energy2> temp;
temp = 2.*mdot* (vecv1[ix]*vecv2[iy]+vecv1[iy]*vecv2[ix]);
temp -= 2.*dotv1k2*(veck1[ix]*vecv2[iy]+veck1[iy]*vecv2[ix]);
temp -= 2.*dotk1v2*(veck2[ix]*vecv1[iy]+veck2[iy]*vecv1[ix]);
temp += 2.*dotv1v2*(veck1[ix]*veck2[iy]+veck1[iy]*veck2[ix]);
ten[ix][iy] = UnitRemoval::InvE2 * temp
-coeff1*tmass2inv*pout_tmp[ix]*pout_tmp[iy];
}
}
// add the g(mu,nu) term
ten[0][0] = ten[0][0]-(coeff1+coeff2);
ten[1][1] = ten[1][1]-(coeff1+coeff2);
ten[2][2] = ten[2][2]-(coeff1+coeff2);
ten[3][3] = ten[3][3]+(coeff1+coeff2);
// overall coefficent
for(int ix=0;ix<4;++ix) {
for(int iy=0;iy<4;++iy) {
ten[ix][iy] = fact*ten[ix][iy];
}
}
// return the wavefunction
return TensorWaveFunction(pout,out,
ten[0][0],ten[0][1],ten[0][2],ten[0][3],
ten[1][0],ten[1][1],ten[1][2],ten[1][3],
ten[2][0],ten[2][1],ten[2][2],ten[2][3],
ten[3][0],ten[3][1],ten[3][2],ten[3][3]);
}
diff --git a/MatrixElement/ME2to2Base.cc b/MatrixElement/ME2to2Base.cc
--- a/MatrixElement/ME2to2Base.cc
+++ b/MatrixElement/ME2to2Base.cc
@@ -1,202 +1,202 @@
// -*- C++ -*-
//
// ME2to2Base.cc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2017 Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the ME2to2Base class.
//
#include "ME2to2Base.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Utilities/SimplePhaseSpace.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include "ThePEG/Cuts/Cuts.h"
using namespace ThePEG;
ME2to2Base::~ME2to2Base() {}
Energy2 ME2to2Base::scale() const {
switch ( scaleChoice() ) {
case 1:
return -tHat()*uHat()/(tHat() + uHat());
default:
return tHat()*uHat()/sHat();
}
}
void ME2to2Base::setKinematics() {
MEBase::setKinematics();
theLastTHat = (meMomenta()[0] - meMomenta()[2]).m2();
theLastUHat = (meMomenta()[1] - meMomenta()[2]).m2();
theLastPhi = meMomenta()[2].phi();
}
bool ME2to2Base::generateKinematics(const double * r) {
// generate the masses of the particles
for ( int i = 2, N = meMomenta().size(); i < N; ++i ) {
meMomenta()[i] = Lorentz5Momentum(mePartonData()[i]->generateMass());
}
double ctmin = -1.0;
double ctmax = 1.0;
Energy q = ZERO;
try {
q = SimplePhaseSpace::
getMagnitude(sHat(), meMomenta()[2].mass(), meMomenta()[3].mass());
- } catch ( ImpossibleKinematics ) {
+ } catch ( ImpossibleKinematics & e ) {
return false;
}
Energy e = sqrt(sHat())/2.0;
Energy2 m22 = meMomenta()[2].mass2();
Energy2 m32 = meMomenta()[3].mass2();
Energy2 e0e2 = 2.0*e*sqrt(sqr(q) + m22);
Energy2 e1e2 = 2.0*e*sqrt(sqr(q) + m22);
Energy2 e0e3 = 2.0*e*sqrt(sqr(q) + m32);
Energy2 e1e3 = 2.0*e*sqrt(sqr(q) + m32);
Energy2 pq = 2.0*e*q;
Energy2 thmin = lastCuts().minTij(mePartonData()[0], mePartonData()[2]);
if ( thmin > ZERO ) ctmax = min(ctmax, (e0e2 - m22 - thmin)/pq);
thmin = lastCuts().minTij(mePartonData()[1], mePartonData()[2]);
if ( thmin > ZERO ) ctmin = max(ctmin, (thmin + m22 - e1e2)/pq);
thmin = lastCuts().minTij(mePartonData()[1], mePartonData()[3]);
if ( thmin > ZERO ) ctmax = min(ctmax, (e1e3 - m32 - thmin)/pq);
thmin = lastCuts().minTij(mePartonData()[0], mePartonData()[3]);
if ( thmin > ZERO ) ctmin = max(ctmin, (thmin + m32 - e0e3)/pq);
Energy ptmin = max(lastCuts().minKT(mePartonData()[2]),
lastCuts().minKT(mePartonData()[3]));
if ( ptmin > ZERO ) {
double ctm = 1.0 - sqr(ptmin/q);
if ( ctm <= 0.0 ) return false;
ctmin = max(ctmin, -sqrt(ctm));
ctmax = min(ctmax, sqrt(ctm));
}
double ymin2 = lastCuts().minYStar(mePartonData()[2]);
double ymax2 = lastCuts().maxYStar(mePartonData()[2]);
double ymin3 = lastCuts().minYStar(mePartonData()[3]);
double ymax3 = lastCuts().maxYStar(mePartonData()[3]);
double ytot = lastCuts().Y() + lastCuts().currentYHat();
if ( ymin2 + ytot > -0.9*Constants::MaxRapidity )
ctmin = max(ctmin, sqrt(sqr(q) + m22)*tanh(ymin2)/q);
if ( ymax2 + ytot < 0.9*Constants::MaxRapidity )
ctmax = min(ctmax, sqrt(sqr(q) + m22)*tanh(ymax2)/q);
if ( ymin3 + ytot > -0.9*Constants::MaxRapidity )
ctmax = min(ctmax, sqrt(sqr(q) + m32)*tanh(-ymin3)/q);
if ( ymax3 + ytot < 0.9*Constants::MaxRapidity )
ctmin = max(ctmin, sqrt(sqr(q) + m32)*tanh(-ymax3)/q);
if ( ctmin >= ctmax ) return false;
double cth = getCosTheta(ctmin, ctmax, r);
Energy pt = q*sqrt(1.0-sqr(cth));
theLastPhi = rnd(2.0*Constants::pi);
meMomenta()[2].setVect(Momentum3( pt*sin(phi()), pt*cos(phi()), q*cth));
meMomenta()[3].setVect(Momentum3(-pt*sin(phi()), -pt*cos(phi()), -q*cth));
meMomenta()[2].rescaleEnergy();
meMomenta()[3].rescaleEnergy();
vector<LorentzMomentum> out(2);
out[0] = meMomenta()[2];
out[1] = meMomenta()[3];
tcPDVector tout(2);
tout[0] = mePartonData()[2];
tout[1] = mePartonData()[3];
if ( !lastCuts().passCuts(tout, out, mePartonData()[0], mePartonData()[1]) )
return false;
theLastTHat = pq*cth + m22 - e0e2;
theLastUHat = m22 + m32 - sHat() - theLastTHat;
jacobian((pq/sHat())*Constants::pi*jacobian());
return true;
}
double ME2to2Base::getCosTheta(double ctmin, double ctmax, const double * r) {
double cth = 0.0;
static const double eps = 1.0e-6;
if ( 1.0 + ctmin <= eps && 1.0 - ctmax <= eps ) {
jacobian(ctmax - ctmin);
cth = ctmin + (*r)*jacobian();
} else if ( 1.0 + ctmin <= eps ) {
cth = 1.0 - (1.0 - ctmax)*pow((1.0 - ctmin)/(1.0 - ctmax), *r);
jacobian(log((1.0 - ctmin)/(1.0 - ctmax))*(1.0 - cth));
} else if ( 1.0 - ctmax <= eps ) {
cth = -1.0 + (1.0 + ctmin)*pow((1.0 + ctmax)/(1.0 + ctmin), *r);
jacobian(log((1.0 + ctmax)/(1.0 + ctmin))*(1.0 + cth));
} else {
double zmin = 0.5*(1.0 - ctmax);
double zmax = 0.5*(1.0 - ctmin);
double A1 = -ctmin/(zmax*(1.0-zmax));
double A0 = -ctmax/(zmin*(1.0-zmin));
double A = *r*(A1 - A0) + A0;
double z = A < 2.0? 2.0/(sqrt(sqr(A) + 4.0) + 2 - A):
0.5*(A - 2.0 + sqrt(sqr(A) + 4.0))/A;
cth = 1.0 - 2.0*z;
jacobian(2.0*(A1 - A0)*sqr(z)*sqr(1.0 - z)/(sqr(z) + sqr(1.0 - z)));
}
return cth;
}
CrossSection ME2to2Base::dSigHatDR() const {
return me2()*jacobian()/(16.0*sqr(Constants::pi)*sHat())*sqr(hbarc);
}
void ME2to2Base::persistentOutput(PersistentOStream & os) const {
os << theScaleChoice << ounit(theLastTHat, GeV2) << ounit(theLastUHat, GeV2)
<< theLastPhi;
}
void ME2to2Base::persistentInput(PersistentIStream & is, int) {
is >> theScaleChoice >> iunit(theLastTHat, GeV2) >> iunit(theLastUHat, GeV2)
>> theLastPhi;
}
AbstractClassDescription<ME2to2Base> ME2to2Base::initME2to2Base;
// Definition of the static class description member.
Switch<ME2to2Base,int> & ME2to2Base::interfaceScaleChoice() {
static Switch<ME2to2Base,int> dummy
("ScaleChoice",
"Different options for calculating the scale of the generated "
"hard sub-process.",
&ME2to2Base::theScaleChoice, 0, false, false);
return dummy;
}
void ME2to2Base::Init() {
static ClassDocumentation<ME2to2Base> documentation
("The ThePEG::ME2to2Base class may be used as a base class "
"for all \\f$2\\rightarrow 2\\f$ matrix elements.");
static SwitchOption interfaceScaleChoice0
(interfaceScaleChoice(),
"that.uhat/shat", "\\f$\\hat{t}\\hat{u}/\\hat{s}\\f$", 0);
static SwitchOption interfaceScaleChoice1
(interfaceScaleChoice(),
"that.uhat/(that+uhat)",
"\\f$-\\hat{t}\\hat{u}/(\\hat{t}+\\hat{u})\\f$", 1);
}
diff --git a/PDF/PartonExtractor.cc b/PDF/PartonExtractor.cc
--- a/PDF/PartonExtractor.cc
+++ b/PDF/PartonExtractor.cc
@@ -1,653 +1,653 @@
// -*- C++ -*-
//
// PartonExtractor.cc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2017 Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the PartonExtractor class.
//
#include "PartonExtractor.h"
#include "ThePEG/Handlers/XComb.h"
#include "ThePEG/Cuts/Cuts.h"
#include "ThePEG/PDF/NoPDF.h"
#include "ThePEG/PDF/RemnantHandler.h"
#include "ThePEG/PDF/BeamParticleData.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/EventRecord/Step.h"
#include "ThePEG/EventRecord/SubProcess.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/SimplePhaseSpace.h"
#include "ThePEG/Utilities/UtilityBase.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/PDT/EnumParticles.h"
using namespace ThePEG;
PartonExtractor::PartonExtractor()
: theMaxTries(100), flatSHatY(false) {}
PartonExtractor::~PartonExtractor() {}
IBPtr PartonExtractor::clone() const {
return new_ptr(*this);
}
IBPtr PartonExtractor::fullclone() const {
return new_ptr(*this);
}
PartonPairVec PartonExtractor::
getPartons(Energy maxEnergy, const cPDPair & incoming,
const Cuts & kc) const {
PartonPairVec result;
PartonVector first;
PDFCuts cuts1(kc, true, maxEnergy);
PBPtr p1 =
new_ptr(PartonBin(PDPtr(), PBPtr(), incoming.first, PDFPtr(), cuts1));
addPartons(p1, cuts1, theFirstPDF, first);
PartonVector second;
PDFCuts cuts2(kc, false, maxEnergy);
PBPtr p2 =
new_ptr(PartonBin(PDPtr(), PBPtr(), incoming.second, PDFPtr(), cuts2));
addPartons(p2, cuts2, theSecondPDF, second);
for ( PartonVector::iterator it1 = first.begin();
it1 != first.end(); ++it1 )
for ( PartonVector::iterator it2 = second.begin();
it2 != second.end(); ++it2 )
result.push_back(PBPair(*it1, *it2));
// We add the original parton bins as well to avoid them being
// deleted.
result.push_back(PBPair(p1, p2));
return result;
}
void PartonExtractor::
addPartons(tPBPtr incoming, const PDFCuts & cuts, tcPDFPtr pdf,
PartonVector & pbins) const {
if(!pdf) pdf = getPDF(incoming->parton());
if ( dynamic_ptr_cast<Ptr<NoPDF>::tcp>(pdf) ||
incoming->parton() == incoming->particle() ) {
pbins.push_back(incoming);
return;
}
cPDVector partons = pdf->partons(incoming->parton());
for ( int i = 0, N = partons.size(); i < N; ++i ) {
PBPtr pb =
new_ptr(PartonBin(incoming->parton(), incoming, partons[i], pdf, cuts));
incoming->addOutgoing(pb);
addPartons(pb, cuts, PDFPtr(), pbins);
}
}
tcPDFPtr PartonExtractor::getPDF(tcPDPtr particle) const {
for ( vector<PDFPtr>::const_iterator it = theSpecialDensities.begin();
it != theSpecialDensities.end(); ++it )
if ( (**it).canHandle(particle) ) return *it;
Ptr<BeamParticleData>::tcp p =
dynamic_ptr_cast<Ptr<BeamParticleData>::tcp>(particle);
if ( !p || !p->pdf() ) return noPDF();
return p->pdf();
}
void PartonExtractor::select(tXCombPtr newXComb) {
theLastXComb = newXComb;
}
tPBIPtr PartonExtractor::partonBinInstance(tcPPtr p) const {
PartonBinInstanceMap::const_iterator it = partonBinInstances().find(p);
return it == partonBinInstances().end()? PBIPtr(): it->second;
}
void PartonExtractor::
colourConnect(tPPtr particle, tPPtr parton, const tPVector & remnants) const {
// Sorry cannot handle coloured resolved particles.
if ( particle->coloured() ) throw RemColException(*this);
// First connect the loose colour line from the extacted parton.
if ( parton->hasColour() )
findConnect(parton->colourLine(), parton, true,
remnants.rbegin(), remnants.rend());
// First connect the loose anti-colour line from the extacted parton.
if ( parton->hasAntiColour() )
findConnect(parton->antiColourLine(), parton, false,
remnants.begin(), remnants.end());
// Go through the rest of the remnants and create new colour lines
// if needed. Go through it forwards and backwards to catch possible
// inconsistencies.
for ( tPVector::const_iterator it = remnants.begin();
it != remnants.end(); ++it ) {
if ( (**it).hasAntiColour() && !(**it).antiColourLine() )
findConnect(ColourLine::create(*it, true), *it, false,
it + 1, remnants.end());
}
for ( tPVector::const_reverse_iterator it = remnants.rbegin();
it != remnants.rend(); ++it ) {
if ( (**it).hasColour() && !(**it).colourLine() )
findConnect(ColourLine::create(*it), *it, true, it + 1, remnants.rend());
}
}
Energy2 PartonExtractor::newScale() {
return lastScale();
}
pair<int,int> PartonExtractor::nDims(const PBPair & pbins) {
// if photon from a lepton or proton generate scale
bool genscale[2]={false,false};
for(unsigned int ix=0;ix<2;++ix) {
PBPtr bin = ix==0 ? pbins.first : pbins.second;
if (!bin || !bin->particle() || !bin->parton()) continue;
if(bin->pdf()->partons(bin->particle()).size()==1 &&
bin->particle()->id()!=bin->parton()->id())
genscale[ix]=true;
}
return make_pair(pbins.first ->nDim(genscale[0]),
pbins.second->nDim(genscale[1]));
}
void PartonExtractor::prepare(const PBIPair & pbins) {
partonBinInstances().clear();
pbins.first->prepare();
pbins.second->prepare();
}
void PartonExtractor::updatePartonBinInstances(const PBIPair & pbins) {
partonBinInstances().clear();
tPBIPtr current = pbins.first;
while ( current->incoming() ) {
partonBinInstances()[current->parton()] = current;
current = current->incoming();
}
current = pbins.second;
while ( current->incoming() ) {
partonBinInstances()[current->parton()] = current;
current = current->incoming();
}
}
bool PartonExtractor::
generateL(const PBIPair & pbins, const double * r1, const double * r2) {
Direction<0> dir(true);
generateL(*pbins.first, r1);
dir.reverse();
generateL(*pbins.second, r2);
if ( !flatSHatY || pbins.first->hasPoleIn1() || pbins.second->hasPoleIn1() )
return true;
Energy2 shmax = lastCuts().sHatMax();
Energy2 shmin = lastCuts().sHatMin();
Energy2 sh = shmin*pow(shmax/shmin, *r1);
double ymax = lastCuts().yHatMax();
double ymin = lastCuts().yHatMin();
double km = log(shmax/shmin);
ymax = min(ymax, log(lastCuts().x1Max()*sqrt(lastS()/sh)));
ymin = max(ymin, -log(lastCuts().x2Max()*sqrt(lastS()/sh)));
double y = ymin + (*r2)*(ymax - ymin);
double l1 = 0.5*log(lastS()/sh) - y;
double l2 = 0.5*log(lastS()/sh) + y;
pbins.first->li(l1 - pbins.first->l() + pbins.first->li());
pbins.first->l(l1);
pbins.first->jacobian(km*(ymax - ymin));
pbins.second->li(l2 - pbins.second->l() + pbins.second->li());
pbins.second->l(l2);
pbins.second->jacobian(1.0);
return ( pbins.first->li() >= 0.0 && pbins.second->li() >= 0.0 );
}
Energy2 PartonExtractor::
generateSHat(Energy2, const PBIPair & pbins,
const double * r1, const double * r2,
bool haveMEPartons) {
Direction<0> dir(true);
if(pbins.first->bin()->pdfDim()<=1) pbins.first->scale(-lastScale());
if ( !generate(*pbins.first, r1, lastSHat(),
pbins.first->getFirst()->parton()->momentum(),
haveMEPartons) )
return -1.0*GeV2;
dir.reverse();
if(pbins.second->bin()->pdfDim()<=1) pbins.second->scale(-lastScale());
if ( !generate(*pbins.second, r2, lastSHat(),
pbins.second->getFirst()->parton()->momentum(),
haveMEPartons) )
return -1.0*GeV2;
return (pbins.first->parton()->momentum() +
pbins.second->parton()->momentum()).m2();
}
void PartonExtractor::
generateL(PartonBinInstance & pb, const double * r) {
if ( !pb.incoming() ) return;
pb.parton(pb.partonData()->produceParticle(Lorentz5Momentum()));
generateL(*pb.incoming(), r + pb.bin()->pdfDim() + pb.bin()->remDim());
pb.particle(pb.incoming()->parton());
if ( pb.li() >= 0 ) return;
double jac = 1.0;
if ( pb.bin()->pdfDim() )
pb.li(pb.pdf()->flattenL(pb.particleData(), pb.partonData(),
pb.bin()->cuts(), *r++, jac));
pb.scale(-1.0*GeV2);
if ( pb.bin()->pdfDim() > 1 )
pb.scale(pb.pdf()->flattenScale(pb.particleData(), pb.partonData(),
pb.bin()->cuts(), pb.li(), *r++, jac)
*pb.bin()->cuts().scaleMaxL(pb.li()));
pb.jacobian(jac);
pb.l(pb.incoming()->l() + pb.li());
}
bool PartonExtractor::
generate(PartonBinInstance & pb, const double * r,
Energy2 shat, const Lorentz5Momentum & first,
bool haveMEPartons) {
if ( !pb.incoming() ) return true;
if ( !generate(*pb.incoming(), r + pb.bin()->pdfDim() + pb.bin()->remDim(),
shat/pb.xi(), first) )
return false;
pb.remnantWeight(1.0);
pb.parton()->setMomentum
(pb.remnantHandler()->generate(pb, r + pb.bin()->pdfDim(), pb.scale(), shat,
pb.particle()->momentum(),haveMEPartons));
if ( pb.remnantWeight() <= 0.0 ) return false;
partonBinInstances()[pb.parton()] = &pb;
return true;
}
void PartonExtractor::
constructRemnants(const PBIPair & pbins, tSubProPtr sub, tStepPtr step) const {
partonBinInstances().clear();
LorentzMomentum k1 = pbins.first->parton()->momentum();
LorentzMomentum k2 = pbins.second->parton()->momentum();
LorentzMomentum Ph = k1 + k2;
LorentzMomentum Phold = Ph;
LorentzRotation Rh = Utilities::getBoostToCM(make_pair(k1, k2));
bool pickside = rndbool();
if ( pickside && pbins.first->incoming() ) {
Direction<0> dir(true);
constructRemnants(*pbins.first, Ph, k2);
construct(*pbins.first, step, false);
}
if ( pbins.second->incoming() ) {
Direction<0> dir(false);
constructRemnants(*pbins.second, Ph, pbins.first->parton()->momentum());
construct(*pbins.second, step, false);
}
if ( (!pickside) && pbins.first->incoming() ) {
Direction<0> dir(true);
constructRemnants(*pbins.first, Ph, pbins.second->parton()->momentum());
construct(*pbins.first, step, false);
}
// LorentzRotation rot = Utilities::transformToMomentum(Phold, Ph);
k1 = pbins.first->parton()->momentum();
k2 = pbins.second->parton()->momentum();
LorentzRotation rot = Utilities::getBoostFromCM(make_pair(k1, k2))*Rh;
Utilities::transform(sub->outgoing(), rot);
Utilities::transform(sub->intermediates(), rot);
Ph = k1 + k2;
if ( abs(Ph.m2() - Phold.m2())/Phold.m2() > 0.000001 )
cerr << Ph.m2()/GeV2 << " was (" << Phold.m2()/GeV2 << ")" << endl;
}
void PartonExtractor::
constructRemnants(PartonBinInstance & pb, LorentzMomentum & Ph,
const LorentzMomentum & k) const {
LorentzMomentum P = pb.particle()->momentum();
DVector r = UseRandom::rndvec(pb.bin()->remDim());
if ( r.empty() ) r.push_back(0.0);
pb.parton()->setMomentum(pb.remnantHandler()->
generate(pb, &r[0], pb.scale(), Ph.m2(), P));
if ( pb.remnantWeight() <= 0.0 ) throw Veto();
pb.remnantHandler()->boostRemnants(pb);
LorentzMomentum Pr = Utilities::sumMomentum(pb.remnants());
transformRemnants(Ph, Pr, k, pb.particle()->momentum());
pb.parton()->setMomentum(pb.particle()->momentum() - Pr);
try {
Utilities::setMomentum(pb.remnants().begin(),
pb.remnants().end(),
static_cast<const LorentzMomentum &>(Pr));
}
catch ( ThePEG::Exception & e) {
throw e;
}
catch ( ThePEG::Veto ) {
throw;
}
catch ( std::exception & e ) {
throw Exception() << "Caught non-ThePEG exception " << e.what() << "in "
<< "PartonExtractor::constructRemnants"
<< Exception::eventerror;
}
partonBinInstances()[pb.parton()] = &pb;
if ( !pb.incoming()->incoming() ) return;
// We get here if we need to construct remnants recursively.
LorentzMomentum Phnew = Ph + Pr;
constructRemnants(*pb.incoming(), Phnew, k);
LorentzRotation rot = Utilities::transformToMomentum(Ph + Pr, Phnew);
Utilities::transform(pb.remnants(), rot);
Ph.transform(rot);
}
LorentzRotation PartonExtractor::
boostRemnants(PBIPair & bins, LorentzMomentum k1, LorentzMomentum k2,
bool side1, bool side2) const {
if ( !side1 && !side2 ) return LorentzRotation();
LorentzMomentum P1 =
bins.first? LorentzMomentum(bins.first->parton()->momentum()): k1;
LorentzMomentum Pr1;
if ( side1 ) {
P1 = bins.first->particle()->momentum();
Pr1 = Utilities::sumMomentum(bins.first->remnants());
}
LorentzMomentum P2 =
bins.second? LorentzMomentum(bins.second->parton()->momentum()): k2;
LorentzMomentum Pr2;
if ( side2 ) {
P2 = bins.second->particle()->momentum();
Pr2 = Utilities::sumMomentum(bins.second->remnants());
}
LorentzRotation Rh = Utilities::getBoostToCM(make_pair(k1, k2));
LorentzMomentum Ph = k1 + k2;
// LorentzMomentum Phold = Ph;
bool otherside = rndbool();
if ( otherside && side2 ){
Direction<0> dir(false);
transformRemnants(Ph, Pr2, k1, P2);
k2 = P2 - Pr2;
}
if ( side1 ){
Direction<0> dir(true);
transformRemnants(Ph, Pr1, k2, P1);
k1 = P1 - Pr1;
}
if ( side2 && !otherside ) {
Direction<0> dir(false);
transformRemnants(Ph, Pr2, k1, P2);
k2 = P2 - Pr2;
}
if ( bins.first ) {
if ( bins.first->remnants().size() == 1 )
bins.first->remnants()[0]->setMomentum(Pr1);
else
Utilities::setMomentum(bins.first->remnants().begin(),
bins.first->remnants().end(),
static_cast<const LorentzMomentum &>(Pr1));
bins.first->parton()->setMomentum(k1);
}
if ( bins.second ) {
if ( bins.second->remnants().size() == 1 )
bins.second->remnants()[0]->setMomentum(Pr2);
else
Utilities::setMomentum(bins.second->remnants().begin(),
bins.second->remnants().end(),
static_cast<const LorentzMomentum &>(Pr2));
bins.second->parton()->setMomentum(k2);
}
Rh.transform(Utilities::getBoostFromCM(make_pair(k1, k2)));
// LorentzMomentum phh = Rh*Phold;
return Rh;
// return Utilities::transformToMomentum(Phold, Ph);
}
void PartonExtractor::
transformRemnants(LorentzMomentum & Ph, LorentzMomentum & Pr,
const LorentzMomentum & k, const LorentzMomentum & P) const {
// don't do this for very soft remnants, as
// we may run into numerical troubles; threshold
// needs to become a parameter at some point
if ( Pr.vect().mag2()/k.vect().mag2() < 1e-10 &&
sqr(Pr.e()/k.e()) < 1e-10 )
return;
TransverseMomentum pt = Pr;
try {
if ( Direction<0>::pos() )
SimplePhaseSpace::CMS(Pr, Ph, (P + k).m2(), 1.0, 0.0);
else
SimplePhaseSpace::CMS(Ph, Pr, (k + P).m2(), 1.0, 0.0);
LorentzRotation rpt;
if ( sqr(Pr.z()) > ZERO ) rpt.rotateY(asin(pt.pt()/Pr.z()));
rpt.rotateZ(pt.phi());
rpt = Direction<0>::pos()?
Utilities::getBoostFromCM(make_pair(P, k))*rpt:
Utilities::getBoostFromCM(make_pair(k, P))*rpt;
Ph.transform(rpt);
Pr.transform(rpt);
- } catch ( ImpossibleKinematics ) {}
+ } catch ( ImpossibleKinematics & e ) {}
}
double PartonExtractor::fullFn(const PBIPair & pbins, Energy2 scale,
pair<bool,bool> noLastPDF) {
if(pbins.first->bin()->pdfDim()<=1) pbins.first->scale(scale);
if(pbins.second->bin()->pdfDim()<=1) pbins.second->scale(scale);
return fullFn(*pbins.first,noLastPDF.first)*fullFn(*pbins.second,noLastPDF.second);
}
double PartonExtractor::fullFn(const PartonBinInstance & pb,
bool noLastPDF) {
if ( !pb.incoming() ) return 1.0;
if (noLastPDF)
return
fullFn(*pb.incoming(),false) * pb.jacobian() *
pb.remnantWeight() * exp(-pb.li());
return fullFn(*pb.incoming(),false) * pb.jacobian() * pb.remnantWeight() *
pb.pdf()->xfl(pb.particleData(), pb.partonData(), pb.scale(),
pb.li(), pb.incoming()->scale());
}
void PartonExtractor::
construct(const PBIPair & pbins, tStepPtr step) const {
// if a long chain we need to break some mother/child relationships
if(pbins.first->incoming()) {
if(pbins.first->incoming()->incoming()) {
if(!pbins.first->parton()->parents().empty()) {
tParticleVector parents=pbins.first->parton()->parents();
tPPtr parton = pbins.first->parton();
for(unsigned int ix=0;ix<parents.size();++ix) parents[ix]->abandonChild(parton);
}
}
}
if(pbins.second->incoming()) {
if(pbins.second->incoming()->incoming()) {
if(!pbins.second->parton()->parents().empty()) {
tParticleVector parents=pbins.second->parton()->parents();
tPPtr parton = pbins.second->parton();
for(unsigned int ix=0;ix<parents.size();++ix) parents[ix]->abandonChild(parton);
}
}
}
Direction<0> dir(true);
construct(*pbins.first, step);
dir.reverse();
construct(*pbins.second, step);
}
void PartonExtractor::
construct(PartonBinInstance & pb, tStepPtr step, bool boost) const {
if ( !pb.incoming() ) return;
if ( boost ) pb.remnantHandler()->boostRemnants(pb);
if ( pb.incoming()->incoming() ) {
step->insertIntermediate(pb.particle(),pb.incoming()->particle(),pb.parton());
}
tPVector rem(pb.remnants().begin(), pb.remnants().end());
if ( !step->addDecayProduct(pb.particle(), rem.begin(), rem.end(), false) )
{}
colourConnect(pb.particle(), pb.parton(), rem);
construct(*pb.incoming(), step);
}
PBIPair PartonExtractor::newRemnants(tPPair oldp, tPPair newp, tStepPtr step) {
PBIPair pb;
Direction<0> dir(true);
pb.first = newRemnants(partonBinInstance(oldp.first),
newp.first, newp.second->momentum());
dir.reverse();
pb.second = newRemnants(partonBinInstance(oldp.second),
newp.second, newp.first->momentum());
addNewRemnants(partonBinInstance(oldp.first), pb.first, step);
addNewRemnants(partonBinInstance(oldp.second), pb.second, step);
return pb;
}
PBIPtr PartonExtractor::
newRemnants(tPBIPtr oldpb, tPPtr newp, const LorentzMomentum & k) {
if ( ! oldpb || !oldpb->incoming() ) return oldpb;
Energy2 shat = (k + newp->momentum()).m2();
// Loop over all possible PartonBin sisters to find the one
// corresponding to the newly extracted parton.
const PartonBin::PBVector & sisters = oldpb->incoming()->bin()->outgoing();
for ( int i = 0, N = sisters.size(); i < N; ++i )
if ( sisters[i]->parton() == newp->dataPtr() ) {
// Setup necessary info in new PartonBinInstance object.
PBIPtr newpb = new_ptr(PartonBinInstance(sisters[i], oldpb->incoming()));
newpb->particle(oldpb->particle());
newpb->parton(newp);
newpb->li(log(oldpb->particle()->momentum().dirPlus()/
newp->momentum().dirPlus()));
newpb->l(oldpb->l() - oldpb->li() + newpb->li());
Energy2 sc = -newp->scale();
newpb->scale(newp->scale());
if ( oldpb->incoming()->incoming() )
sc = -newpb->particle()->momentum().m2();
// Now we can construct the new remnants.
newpb->remnantWeight(1.0);
if ( !newpb->remnantHandler()->
recreateRemnants(*newpb, oldpb->parton(), newp, newpb->li(),
sc, shat, newpb->particle()->momentum()) )
throw Veto();
if ( newpb->remnantWeight() <= 0.0 ) throw Veto();
return newpb;
}
throw Veto();
}
void PartonExtractor::
addNewRemnants(tPBIPtr oldpb, tPBIPtr newpb, tStepPtr step) {
if ( oldpb == newpb ) return;
if ( oldpb->parton() != newpb->parton() ) {
step->removeDecayProduct(newpb->particle(), oldpb->parton());
if ( !step->addDecayProduct(newpb->particle(), newpb->parton()) )
throw Veto();
}
tPVector rem(newpb->remnants().begin(), newpb->remnants().end());
colourConnect(newpb->particle(), newpb->parton(), rem);
partonBinInstances()[newpb->parton()] = newpb;
if ( !step->addDecayProduct(oldpb->remnants().begin(),
oldpb->remnants().end(),
rem.begin(), rem.end()) )
throw Veto();
}
void PartonExtractor::persistentOutput(PersistentOStream & os) const {
os << theLastXComb << theSpecialDensities << theNoPDF << theMaxTries
<< flatSHatY << theFirstPDF << theSecondPDF;
}
void PartonExtractor::persistentInput(PersistentIStream & is, int) {
is >> theLastXComb >> theSpecialDensities >> theNoPDF >> theMaxTries
>> flatSHatY >> theFirstPDF >> theSecondPDF;
}
ClassDescription<PartonExtractor> PartonExtractor::initPartonExtractor;
void PartonExtractor::Init() {
static ClassDocumentation<PartonExtractor> documentation
("There is no documentation for the ThePEG::PartonExtractor class");
static RefVector<PartonExtractor,PDFBase> interfaceSpecialDensities
("SpecialDensities",
"A list of parton density objects to be used for incoming particles "
"overriding possible densities given for particles of the "
"BeamParticleData class.",
&PartonExtractor::theSpecialDensities, 0, false, false, true, false);
static Reference<PartonExtractor,PDFBase> interfaceNoPDF
("NoPDF",
"A fixed reference to a NoPDF object to be used for particles without "
"substructure.",
&PartonExtractor::theNoPDF, true, true, true, false);
static Parameter<PartonExtractor,int> interfaceMaxTries
("MaxTries",
"The maximum number of attempts allowed when trying to generate "
"remnants.",
&PartonExtractor::theMaxTries, 100, 1, 1000, false, false, true);
static Switch<PartonExtractor,bool> interfaceFlatSHatY
("FlatSHatY",
"The possibility to override the l-generation in the PDFs and generate "
"a flat distribution in \\f$\\log(\\hat{s})\\f$ and \\f$y\\f$. This only "
"applies if the parton densities do not have poles in \\f$x=1\\f$.",
&PartonExtractor::flatSHatY, false, false, false);
static SwitchOption interfaceFlatSHatY0
(interfaceFlatSHatY,
"Off", "Use the l-generation defined by the PDFs", false);
static SwitchOption interfaceFlatSHatY1
(interfaceFlatSHatY,
"On", "Generate flat rapidity and \\f$\\log(\\hat{s})\\f$", true);
static SwitchOption interfaceFlatSHatNo
(interfaceFlatSHatY,
"No", "Use the l-generation defined by the PDFs", false);
static SwitchOption interfaceFlatSHatYes
(interfaceFlatSHatY,
"Yes", "Generate flat rapidity and \\f$\\log(\\hat{s})\\f$", true);
static Reference<PartonExtractor,PDFBase> interfaceFirstPDF
("FirstPDF",
"PDF to override the default PDF for the first beam particle",
&PartonExtractor::theFirstPDF, false, false, true, true, false);
static Reference<PartonExtractor,PDFBase> interfaceSecondPDF
("SecondPDF",
"PDF to override the default PDF for the second beam particle",
&PartonExtractor::theSecondPDF, false, false, true, true, false);
}
RemColException::RemColException(const PartonExtractor & pe) {
theMessage << "Parton extractor '" << pe.name() << "' failed to connect "
<< "the colours of the outgoing partons and the remnants.";
severity(maybeabort);
}
void PartonExtractor::dofinish() {
partonBinInstances().clear();
HandlerBase::dofinish();
}
diff --git a/PDT/DecayMode.cc b/PDT/DecayMode.cc
--- a/PDT/DecayMode.cc
+++ b/PDT/DecayMode.cc
@@ -1,691 +1,692 @@
// -*- C++ -*-
//
// DecayMode.cc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2017 Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the DecayMode class.
//
#include "DecayMode.h"
#include "DecayMode.xh"
#include "ThePEG/Repository/Repository.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Utilities/Rebinder.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/EnumIO.h"
using namespace ThePEG;
DecayMode::DecayMode()
: theBrat(0.0), isOn(false) {}
DecayMode::DecayMode(tPDPtr newParticle, double newBrat, bool newOn)
: theBrat(newBrat), isOn(newOn), theParent(newParticle) {}
DecayMode::DecayMode(const DecayMode & dm)
: Interfaced(dm), theTag(dm.theTag), theBrat(dm.theBrat), isOn(dm.isOn),
theParent(dm.theParent), theProducts(dm.theProducts),
theOrderedProducts(dm.theOrderedProducts),
theCascadeProducts(dm.theCascadeProducts), theMatchers(dm.theMatchers),
theWildMatcher(dm.theWildMatcher), theExcluded(dm.theExcluded),
theOverlap(dm.theOverlap), theDecayer(dm.theDecayer),
theAntiPartner(dm.theAntiPartner), theLinks(dm.theLinks) {}
DecayMode::~DecayMode() {}
DMPtr DecayMode::Create(tPDPtr newParent, double newBrat, bool newOn) {
DMPtr dm = new_ptr(DecayMode(newParent, newBrat, newOn));
Repository::Register(dm, newParent->fullName() + "/NEWMODE");
if ( !newParent->CC() ) return dm;
DMPtr adm = new_ptr(DecayMode(newParent->CC(), newBrat, newOn));
Repository::Register(adm, newParent->CC()->fullName() + "/NEWMODE");
dm->theAntiPartner = adm;
adm->theAntiPartner = dm;
return dm;
}
void DecayMode::readSetup(istream & is) {
string decnam;
is >> theBrat >> ienum(isOn) >> decnam;
if ( decnam.empty() ) return;
BaseRepository::DirectoryAppend(decnam);
setDecayer(BaseRepository::GetObject<tDecayerPtr>(decnam));
}
IBPtr DecayMode::clone() const {
return dmclone();
}
DMPtr DecayMode::dmclone() const {
return new_ptr(*this);
}
IBPtr DecayMode::fullclone() const {
DMPtr dm = dmclone();
Repository::Register(dm);
if ( !CC() ) return dm;
DMPtr adm = CC()->dmclone();
Repository::Register(adm);
dm->theAntiPartner = adm;
adm->theAntiPartner = dm;
return dm;
}
void DecayMode::doupdate() {
Interfaced::doupdate();
bool redo = touched();
UpdateChecker::check(decayer(), redo);
if ( !redo ) return;
if ( theBrat > 0.0 && isOn && !decayer()->accept(*this) )
throw DecModNoAccept(tag(), decayer()->name());
}
void DecayMode::rebind(const TranslationMap & trans) {
try {
theParent = trans.alwaysTranslate(theParent);
ParticleMSet newProducts;
trans.alwaysTranslate(inserter(newProducts),
products().begin(), products().end());
products().swap(newProducts);
tPDVector newOrdered;
trans.alwaysTranslate(inserter(newOrdered),
orderedProducts().begin(), orderedProducts().end());
theOrderedProducts.swap(newOrdered);
ModeMSet newCasc;
trans.alwaysTranslate(inserter(newCasc),
cascadeProducts().begin(), cascadeProducts().end());
cascadeProducts().swap(newCasc);
MatcherMSet newMatchers;
trans.alwaysTranslate(inserter(newMatchers),
productMatchers().begin(), productMatchers().end());
productMatchers().swap(newMatchers);
wildProductMatcher() = trans.alwaysTranslate(wildProductMatcher());
ParticleMSet newExclude;
trans.alwaysTranslate(inserter(newExclude),
excluded().begin(), excluded().end());
excluded().swap(newExclude);
theAntiPartner = trans.alwaysTranslate(CC());
for ( int i = 0, N = theLinks.size(); i < N; ++i ) {
theLinks[i].first = trans.alwaysTranslate(theLinks[i].first);
theLinks[i].second = trans.alwaysTranslate(theLinks[i].second);
}
}
catch (IBPtr ip) {
throw DecModRebind(name(), ip->name());
}
catch (...) {
throw DecModRebind(name(), "<unknown>");
}
}
IVector DecayMode::getReferences() {
IVector ret;
ret.push_back(theParent);
ret.insert(ret.end(), products().begin(), products().end());
ret.insert(ret.end(), cascadeProducts().begin(), cascadeProducts().end());
ret.insert(ret.end(), productMatchers().begin(), productMatchers().end());
if ( wildProductMatcher() ) ret.push_back(wildProductMatcher());
ret.insert(ret.end(), excluded().begin(), excluded().end());
if ( CC() ) ret.push_back(CC());
return ret;
}
bool DecayMode::addOverlap(tcDMPtr d) {
bool inc = includes(*d);
if ( !inc ) return false;
if ( find(theOverlap.begin(), theOverlap.end(), d) != theOverlap.end() )
return true;
theOverlap.push_back(d);
return true;
}
void DecayMode::resetOverlap() {
theOverlap.clear();
}
bool DecayMode::
compareId(const ParticleMSet & s1, const ParticleMSet & si2) const {
if ( generator() ) return s1 == si2;
ParticleMSet s2 = si2;
for ( ParticleMSet::const_iterator p1 = s1.begin(); p1 != s1.end(); ++p1 ) {
ParticleMSet::const_iterator p2 = s2.begin();
while ( p2 != s2.end() && (**p2).id() != (**p1).id() ) ++p2;
if ( p2 == s2.end() ) return false;
s2.erase(p2);
}
return s2.empty();
}
ParticleMSet::const_iterator
DecayMode::findId(const ParticleMSet & s, const ParticleData & p) const {
for ( ParticleMSet::const_iterator pit = s.begin(); pit != s.end(); ++pit )
if ( (**pit).id() == p.id() ) return pit;
return s.end();
}
struct IdCmp {
bool operator()(tcPDPtr p1, tcPDPtr p2) { return p1->id() == p2->id(); }
};
bool DecayMode::includes(const DecayMode & d) const {
// Fast check for ordinary decay modes.
if ( cascadeProducts().empty() && productMatchers().empty() &&
excluded().empty() && !wildProductMatcher() &&
d.cascadeProducts().empty() && d.productMatchers().empty() &&
d.excluded().empty() && !d.wildProductMatcher() ) {
if ( links().size() != d.links().size() ) return false;
if ( !compareId(products(), d.products()) ) return false;
LinkVector dlinks = d.links();
for ( int i = 0, N = links().size(); i < N; ++i ) {
for ( int j = 0, M = dlinks.size(); j < M; ++j ) {
if ( ( links()[i].first->id() == dlinks[j].first->id() &&
links()[i].second->id() == dlinks[j].second->id() ) ||
( links()[i].first->id() == dlinks[j].second->id() &&
links()[i].second->id() == dlinks[j].first->id() ) ) {
dlinks.erase(dlinks.begin() + j);
break;
}
}
return false;
}
return dlinks.empty();
}
// First check that none of the excluded products in this are
// present in the other.
ParticleMSet::const_iterator pit;
for ( pit = excluded().begin(); pit != excluded().end(); ++pit ) {
if ( findId(d.products(), **pit ) != d.products().end() ) return false;
for ( ModeMSet::const_iterator mit = d.cascadeProducts().begin();
mit != d.cascadeProducts().end(); ++mit )
if ( (**pit).id() == (**mit).parent()->id() ) return false;
}
// Check that all cascade decays in this overlaps with one in the
// other. Save the ones that are left
ModeMSet cascleft = d.cascadeProducts();
for ( ModeMSet::iterator mit = cascadeProducts().begin();
mit != cascadeProducts().end(); ++mit ) {
ModeMSet::iterator mit2 = cascleft.begin();
while ( mit2 != cascleft.end() && !(**mit).includes(**mit2) ) ++mit2;
if ( mit2 == cascleft.end() ) return false;
}
// Check that all cascade product parents in the other matches
// something in this. Otherwise expand the cascade product.
ParticleMSet partleft = d.products();
MatcherMSet matchleft = d.productMatchers();
ParticleMSet excludeleft = d.excluded();
MatcherMSet wildleft;
if ( d.wildProductMatcher() ) wildleft.insert(wildProductMatcher());
ParticleMSet part = products();
MatcherMSet match = productMatchers();
while ( cascleft.size() ) {
ModeMSet::iterator cdmit = cascleft.begin();
cDMPtr cdm = *cdmit;
ParticleMSet::iterator pit = findId(part, *(cdm->parent()));
if ( pit != part.end() ) {
cascleft.erase(cdmit);
part.erase(pit);
} else {
MatcherMSet::iterator mit = match.begin();
while ( mit != match.end() && !(**mit).matches(*(cdm->parent())) ) ++mit;
if ( mit != match.end() ) {
cascleft.erase(cdmit);
match.erase(mit);
} else {
if ( wildProductMatcher() &&
wildProductMatcher()->matches(*(cdm->parent())) ) {
cascleft.erase(cdmit);
} else {
cascleft.erase(cdmit);
partleft.insert(cdm->products().begin(), cdm->products().end());
matchleft.insert(cdm->productMatchers().begin(),
cdm->productMatchers().end());
if ( cdm->wildProductMatcher() )
wildleft.insert(cdm->wildProductMatcher());
excludeleft.insert(cdm->excluded().begin(), cdm->excluded().end());
cascleft.insert(cdm->cascadeProducts().begin(),
cdm->cascadeProducts().end());
}
}
}
}
// Check that all excluded left in the other are absent in this.
if ( find_first_of(excludeleft.begin(), excludeleft.end(), part.begin(),
part.end(), IdCmp()) != excludeleft.end() ) return false;
// Now all particles and matches left in this must match something
// in the other.
pit = part.begin();
while ( pit != part.end() ) {
ParticleMSet::iterator pit2 = findId(partleft, **pit++);
if ( pit2 == partleft.end() ) return false;
partleft.erase(pit2);
}
MatcherMSet::const_iterator pmit = match.begin();
while ( pmit != match.end() ) {
ParticleMSet::iterator pit2 = partleft.begin();
while ( pit2 != partleft.end() && ! (**pmit).matches(**pit2) ) ++pit2;
if ( pit2 != partleft.end() ) {
partleft.erase(pit2);
} else {
MatcherMSet::iterator pmit2 = matchleft.begin();
while ( pmit2 != matchleft.end() && ! (**pmit).matches(**pmit2) ) ++pmit2;
if ( pmit2 != matchleft.end() ) {
matchleft.erase(pmit2);
} else
return false;
}
}
// Now all particles and matchers left in the other must be matched
// by the wild match in this.
if ( wildProductMatcher() ) {
pit = partleft.begin();
while ( pit != partleft.end() )
if ( !(wildProductMatcher()->matches(**pit++)) ) return false;
pmit = matchleft.begin();
while (pmit != matchleft.end() )
if ( !(wildProductMatcher()->matches(**pmit++)) ) return false;
pmit = wildleft.begin();
while (pmit != wildleft.end() )
if ( !(wildProductMatcher()->matches(**pmit++)) ) return false;
} else
return partleft.empty() && matchleft.empty() && wildleft.empty();
return true;
}
DMPtr DecayMode::clone(tPDPtr pd) const {
DMPtr dm = dmclone();
dm->theParent = pd;
Repository::Register(dm, pd->fullName() + "/" + dm->name());
if ( !theDecayer || !theDecayer->accept(*dm) ) dm->isOn = false;
if ( pd->CC() ) {
DMPtr adm = CC()? CC()->dmclone(): dmclone();
adm->theParent = pd->CC();
Repository::Register(adm, pd->CC()->fullName() + "/" + adm->name());
dm->theAntiPartner = adm;
adm->theAntiPartner = dm;
if ( !adm->theDecayer->accept(*adm) ) adm->isOn = false;
} else
dm->theAntiPartner = DMPtr();
return dm;
}
void DecayMode::synchronize() {
if ( !CC() ) return;
theBrat = CC()->theBrat;
isOn = CC()->isOn;
theDecayer = CC()->theDecayer;
}
struct ParticleOrdering {
bool 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() );
}
};
struct ModeOrdering {
bool operator()(tcDMPtr d1, tcDMPtr d2) {
ParticleOrdering ord;
return ord(d1->parent(), d2->parent()) ||
( !ord(d2->parent(), d1->parent()) &&
( d1->tag() < d2->tag() ||
( d1->tag() == d2->tag() && d1->fullName() < d2->fullName() ) ) );
}
};
struct MatcherOrdering {
bool operator()(tcPMPtr m1, tcPMPtr m2) {
return m1->name() < m2->name() ||
( m1->name() == m2->name() && m1->fullName() < m2->fullName() );
}
};
PVector DecayMode::produceProducts() const {
PVector ret;
for ( int i = 0, N = orderedProducts().size(); i < N; ++i )
ret.push_back(orderedProducts()[i]->produceParticle());
return ret;
}
string DecayMode::makeTag() const {
string ret;
typedef multiset<tcPDPtr,ParticleOrdering> OrderedParticles;
typedef multiset<tcPMPtr,MatcherOrdering> OrderedMatchers;
typedef multiset<tcDMPtr,ModeOrdering> OrderedModes;
LinkVector dlinks = links();
ret = theParent->PDGName() + "->";
if ( dlinks.empty() ) {
OrderedParticles prod(products().begin(), products().end());
for (OrderedParticles::iterator pit = prod.begin();
pit != prod.end(); ++pit )
ret += (**pit).PDGName() + ",";
} else {
unsigned int dl = 0;
for ( int i = 0, N = orderedProducts().size(); i < N; ++i ) {
if ( dl < dlinks.size() && orderedProducts()[i] == dlinks[dl].first ) {
ret += orderedProducts()[i]->PDGName() + "=";
++dl;
} else
ret += orderedProducts()[i]->PDGName() + ",";
}
}
OrderedModes casc(cascadeProducts().begin(), cascadeProducts().end());
for ( OrderedModes::iterator dmit = casc.begin();dmit != casc.end(); ++dmit )
ret += "[" + (**dmit).tag() + "],";
OrderedMatchers match(productMatchers().begin(), productMatchers().end());
for ( OrderedMatchers::iterator mit = match.begin();
mit != match.end(); ++mit ) ret += "?" +(**mit).name() + ",";
if ( theWildMatcher )
ret += "*" + theWildMatcher->name() + ",";
OrderedParticles ex(excluded().begin(), excluded().end());
for ( OrderedParticles::iterator pit = ex.begin(); pit != ex.end(); ++pit )
ret += "!" + (**pit).PDGName() + ",";
ret[ret.size()-1] = ';';
return ret;
}
void DecayMode::brat(double newBrat) {
theBrat = newBrat;
if ( theBrat <= 0.0 ) switchOff();
if ( CC() && parent()->synchronized() ) CC()->theBrat = newBrat;
}
double DecayMode::brat() const {
return isOn? theDecayer->brat(*this, *theParent, theBrat): 0.0;
}
double DecayMode::brat(const Particle & p) const {
return isOn && p.dataPtr() == parent()?
theDecayer->brat(*this, p, theBrat): 0.0;
}
void DecayMode::switchOn() {
isOn = true;
if ( CC() && parent()->synchronized() ) CC()->isOn = true;
}
void DecayMode::switchOff() {
isOn = false;
if ( CC() && parent()->synchronized() ) CC()->isOn = false;
}
void DecayMode::addProduct(tPDPtr pd) {
products().insert(pd);
theOrderedProducts.push_back(pd);
if ( CC() ) {
CC()->products().insert(pd->CC()? pd->CC(): pd);
CC()->theOrderedProducts.push_back(pd->CC()? pd->CC(): pd);
}
resetTag();
}
void DecayMode::addLink(tPDPtr a, tPDPtr b) {
theLinks.push_back(make_pair(a, b));
if ( CC() ) CC()->theLinks.push_back(make_pair(a->CC()? a->CC(): a,
b->CC()? b->CC(): b));
resetTag();
}
void DecayMode::addCascadeProduct(tDMPtr dm) {
cascadeProducts().insert(dm);
if ( CC() ) CC()->cascadeProducts().insert(dm->CC()? dm->CC(): dm);
resetTag();
}
void DecayMode::addProductMatcher( tPMPtr pm) {
productMatchers().insert(pm);
if ( CC() ) CC()->productMatchers().insert(pm->CC()? pm->CC(): pm);
resetTag();
}
void DecayMode::setWildMatcher(tPMPtr pm) {
wildProductMatcher() = pm;
if ( CC() ) CC()->wildProductMatcher() = pm->CC()? pm->CC(): pm;
resetTag();
}
void DecayMode::addExcluded(tPDPtr pd) {
excluded().insert(pd);
if ( CC() ) CC()->excluded().insert(pd->CC()? pd->CC(): pd);
resetTag();
}
void DecayMode::decayer(tDecayerPtr dec) {
if ( !dec || !dec->accept(*this) )
throw DecModSetupNoAccept(tag(), dec->name());
if ( CC() && parent()->synchronized() ) {
if ( !dec->accept(*CC()) )
throw DecModSetupNoAccept(CC()->tag(), dec->name());
CC()->theDecayer = dec;
}
theDecayer = dec;
}
DMPtr DecayMode::constructDecayMode(string & tag, vector<DMPtr> * save) {
DMPtr rdm;
DMPtr adm;
int level = 0;
string::size_type end = 0;
while ( end < tag.size() && ( tag[end] != ']' || level ) ) {
switch ( tag[end++] ) {
case '[':
++level;
break;
case ']':
--level;
break;
}
}
string::size_type next = tag.find("->");
if ( next == string::npos ) return rdm;
if ( tag.find(';') == string::npos ) return rdm;
tPDPtr pd = Repository::findParticle(tag.substr(0,next));
if ( !pd ) return rdm;
rdm = ptr_new<DMPtr>();
rdm->parent(pd);
if ( pd->CC() ) {
adm = ptr_new<DMPtr>();
adm->parent(pd->CC());
rdm->theAntiPartner = adm;
adm->theAntiPartner = rdm;
}
bool error = false;
tag = tag.substr(next+2);
tPDPtr lastprod;
bool dolink = false;
do {
switch ( tag[0] ) {
case '[':
{
tag = tag.substr(1);
DMPtr cdm = constructDecayMode(tag, save);
if ( save ) save->push_back(cdm);
if ( cdm ) rdm->addCascadeProduct(cdm);
else error = true;
} break;
case '=':
dolink = true;
+ [[fallthrough]];
case ',':
case ']':
tag = tag.substr(1);
break;
case '?':
{
next = min(tag.find(','), tag.find(';'));
tPMPtr pm = Repository::findMatcher(tag.substr(1,next-1));
if ( pm ) rdm->addProductMatcher(pm);
else error = true;
tag = tag.substr(next);
} break;
case '!':
{
next = min(tag.find(','), tag.find(';'));
tPDPtr pd = Repository::findParticle(tag.substr(1,next-1));
if ( pd ) rdm->addExcluded(pd);
else error = true;
tag = tag.substr(next);
} break;
case '*':
{
next = min(tag.find(','), tag.find(';'));
tPMPtr pm = Repository::findMatcher(tag.substr(1,next-1));
if ( pm ) rdm->setWildMatcher(pm);
else error = true;
tag = tag.substr(next);
} break;
default:
{
next = min(tag.find('='), min(tag.find(','), tag.find(';')));
tPDPtr pdp = Repository::findParticle(tag.substr(0,next));
if ( pdp ) rdm->addProduct(pdp);
else error = true;
tag = tag.substr(next);
if ( dolink && lastprod ) {
rdm->addLink(lastprod, pdp);
dolink = false;
}
lastprod = pdp;
} break;
}
} while ( tag[0] != ';' && tag.size() );
if ( tag[0] != ';' || error ) {
return DMPtr();
}
tag = tag.substr(1);
for ( DecaySet::const_iterator dit = pd->decayModes().begin();
dit != pd->decayModes().end(); ++dit )
if ( (**dit).tag() == rdm->tag() ) return *dit;
if ( save ) {
save->push_back(rdm);
save->push_back(adm);
} else {
pd->addDecayMode(rdm);
Repository::Register(rdm, pd->fullName() + "/" + rdm->tag());
if ( adm )
Repository::Register(adm, pd->CC()->fullName() + "/" + adm->tag());
}
return rdm;
}
void DecayMode::persistentOutput(PersistentOStream & os) const {
multiset<tcPDPtr,ParticleOrdering> prod(products().begin(), products().end());
multiset<tcDMPtr,ModeOrdering>
casc(cascadeProducts().begin(), cascadeProducts().end());
multiset<tcPMPtr,MatcherOrdering>
match(productMatchers().begin(), productMatchers().end());
multiset<tcPDPtr,ParticleOrdering> ex(excluded().begin(), excluded().end());
multiset<tcDMPtr,ModeOrdering> ovlap(overlap().begin(), overlap().end());
os << theTag << theBrat << isOn << theParent << prod << theOrderedProducts
<< casc << match << theWildMatcher << ex << ovlap << theDecayer
<< theAntiPartner << theLinks;
}
void DecayMode::persistentInput(PersistentIStream & is, int) {
is >> theTag >> theBrat >> isOn >> theParent >> theProducts
>> theOrderedProducts >> theCascadeProducts >> theMatchers
>> theWildMatcher >> theExcluded >> theOverlap >> theDecayer
>> theAntiPartner >> theLinks;
}
ClassDescription<DecayMode> DecayMode::initDecayMode;
void DecayMode::setOn(long i) {
isOn = i;
}
long DecayMode::getOn() const {
return isOn;
}
void DecayMode::setDecayer(DecayerPtr dp) {
decayer(dp);
}
void DecayMode::Init() {
static ClassDocumentation<DecayMode> documentation
("Represents a specific decay channel of a particle.");
static Parameter<DecayMode,double> interfaceBrat
("BranchingRatio",
"The branching fraction for this decay mode. Note that if the sum of "
"branching ratios for one particle is always renormalized to 1. Also, "
"the decaying particle may change this branching ratio if it has a "
"ThePEG::WidthGenerator object assigned to it. ",
&DecayMode::theBrat, 0.0, 0.0, 1.0, false, false, true);
interfaceBrat.setHasDefault(false);
static Switch<DecayMode> interfaceOn
("OnOff",
"Indicates if the decay mode is switched on or off.",
0, 0, false, false, &DecayMode::setOn, &DecayMode::getOn);
static SwitchOption interfaceOnYes
(interfaceOn, "On", "The decay channel is switched on.", 1);
static SwitchOption interfaceOnNo
(interfaceOn, "Off", "The decay channel is switched off.", 0);
interfaceOn.setHasDefault(false);
static Switch<DecayMode> interfaceActive
("Active",
"Indicates if the decay mode is switched on or off.",
0, 0, false, false, &DecayMode::setOn, &DecayMode::getOn);
static SwitchOption interfaceActiveYes
(interfaceActive, "Yes", "The decay channel is switched on.", 1);
static SwitchOption interfaceActiveNo
(interfaceActive, "No", "The decay channel is switched off.", 0);
interfaceActive.setHasDefault(false);
static Reference<DecayMode,Decayer> interfaceDecayer
("Decayer",
"The ThePEG::Decayer object responsible for performing this decay.",
&DecayMode::theDecayer, false, false, true, false,
&DecayMode::setDecayer);
interfaceBrat.rank(10);
interfaceDecayer.rank(9);
interfaceOn.rank(8);
interfaceActive.rank(8);
}
DecModNoAccept::DecModNoAccept(string tag, string dec) {
theMessage << "The Decayer '" << dec << "' is not capable to "
<< "perform the decay in the DecayMode '" << tag << "'.";
severity(warning);
}
DecModSetupNoAccept::DecModSetupNoAccept(string tag, string dec) {
theMessage << "The Decayer '" << dec << "' is not capable to "
<< "perform the decay in the DecayMode '" << tag << "'.";
severity(warning);
}
DecModRebind::DecModRebind(string tag, string obj) {
theMessage << "'Rebind' of DecayMode '" << tag << "' failed because "
<< "the object '" << obj << "' refered to lacked a translation.";
severity(abortnow);
}
diff --git a/PDT/ParticleData.cc b/PDT/ParticleData.cc
--- a/PDT/ParticleData.cc
+++ b/PDT/ParticleData.cc
@@ -1,1038 +1,1038 @@
// -*- C++ -*-
//
// ParticleData.cc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2017 Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the ParticleData class.
//
#include "ParticleData.h"
#include "ParticleData.xh"
#include "ThePEG/PDT/DecayMode.h"
#include "ThePEG/Utilities/HoldFlag.h"
#include "ThePEG/Utilities/Rebinder.h"
#include "ThePEG/Utilities/StringUtils.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Command.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Repository/Repository.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Config/algorithm.h"
#include "ThePEG/Utilities/Exception.h"
#include "ThePEG/Utilities/EnumIO.h"
#include "ThePEG/Repository/UseRandom.h"
namespace ThePEG {
ParticleData::ParticleData()
: theId(0), thePDGName(""), theMass(-1.0*GeV), theWidth(-1.0*GeV),
theHardProcessMass(-1.0*GeV), hardProcessMassSet(false),
theHardProcessWidth(-1.0*GeV), hardProcessWidthSet(false),
theWidthUpCut(-1.0*GeV), theWidthLoCut(-1.0*GeV), theCTau(-1.0*mm),
theCharge(PDT::ChargeUnknown),
theSpin(PDT::SpinUnknown), theColour(PDT::ColourUnknown), isStable(true),
theVariableRatio(false), syncAnti(false), theDefMass(-1.0*GeV),
theDefWidth(-1.0*GeV), theDefCut(-1.0*GeV), theDefCTau(-1.0*mm),
theDefCharge(PDT::ChargeUnknown), theDefSpin(PDT::SpinUnknown),
theDefColour(PDT::ColourUnknown) {}
ParticleData::
ParticleData(PID newId, const string & newPDGName)
: theId(newId), thePDGName(newPDGName), theMass(-1.0*GeV), theWidth(-1.0*GeV),
theHardProcessMass(-1.0*GeV), hardProcessMassSet(false),
theHardProcessWidth(-1.0*GeV), hardProcessWidthSet(false),
theWidthUpCut(-1.0*GeV), theWidthLoCut(-1.0*GeV), theCTau(-1.0*mm),
theCharge(PDT::ChargeUnknown),
theSpin(PDT::SpinUnknown), theColour(PDT::ColourUnknown), isStable(true),
theVariableRatio(false), syncAnti(false), theDefMass(-1.0*GeV),
theDefWidth(-1.0*GeV), theDefCut(-1.0*GeV), theDefCTau(-1.0*mm),
theDefCharge(PDT::ChargeUnknown), theDefSpin(PDT::SpinUnknown),
theDefColour(PDT::ColourUnknown) {}
ParticleData::~ParticleData() {}
PDPtr ParticleData::Create(PID newId, const string & newPDGName) {
return new_ptr(ParticleData(newId, newPDGName));
}
PDPair ParticleData::
Create(PID newId, const string & newPDGName, const string & newAntiPDGName) {
PDPair pap;
pap.first = new_ptr(ParticleData(newId, newPDGName));
pap.second = new_ptr(ParticleData(-newId, newAntiPDGName));
antiSetup(pap);
return pap;
}
void ParticleData::readSetup(istream & is) {
long id;
is >> id >> thePDGName >> iunit(theDefMass, GeV) >> iunit(theDefWidth, GeV)
>> iunit(theDefCut, GeV) >> iunit(theDefCTau, mm) >> ienum(theDefCharge)
>> ienum(theDefColour) >> ienum(theDefSpin) >> ienum(isStable);
theId = id;
theMass = theDefMass;
theWidth = theDefWidth;
theWidthUpCut = theDefCut;
theWidthLoCut = theDefCut;
theCTau = theDefCTau;
theCharge = theDefCharge;
theColour = theDefColour;
theSpin = theDefSpin;
if ( PDGName() == "-" ) thePDGName = name();
return;
}
void ParticleData::antiSetup(const PDPair & pap) {
pap.first->theAntiPartner = pap.second;
pap.second->theAntiPartner = pap.first;
pap.first->syncAnti = pap.second->syncAnti = true;
}
PDPtr ParticleData::pdclone() const {
return new_ptr(*this);
}
IBPtr ParticleData::clone() const {
return pdclone();
}
IBPtr ParticleData::fullclone() const {
PDPtr pd = pdclone();
Repository::Register(pd);
pd->theDecaySelector.clear();
pd->theDecayModes.clear();
pd->isStable = true;
PDPtr apd;
if ( CC() ) {
apd = CC()->pdclone();
Repository::Register(apd);
apd->theDecaySelector.clear();
apd->theDecayModes.clear();
apd->isStable = true;
pd->theAntiPartner = apd;
apd->theAntiPartner = pd;
pd->syncAnti = syncAnti;
apd->syncAnti = CC()->syncAnti;
}
HoldFlag<> dosync(pd->syncAnti, true);
for ( DecaySet::const_iterator it = theDecayModes.begin();
it != theDecayModes.end(); ++it )
pd->addDecayMode(*it);
return pd;
}
Energy ParticleData::mass(Energy mi) {
theMass = mi;
if ( synchronized() && CC() ) CC()->theMass = theMass;
return theMass;
}
Energy ParticleData::width(Energy wi) {
theWidth = wi;
if ( synchronized() && CC() ) CC()->theWidth = theWidth;
return theWidth;
}
Energy ParticleData::widthUpCut(Energy wci) {
theWidthUpCut = wci;
if ( synchronized() && CC() ) CC()->theWidthUpCut = theWidthUpCut;
return theWidthUpCut;
}
Energy ParticleData::widthLoCut(Energy wci) {
theWidthLoCut = wci;
if ( synchronized() && CC() ) CC()->theWidthLoCut = theWidthLoCut;
return theWidthLoCut;
}
Length ParticleData::cTau(Length ti) {
theCTau = ti;
if ( synchronized() && CC() ) CC()->theCTau = theCTau;
return theCTau;
}
PDT::Charge ParticleData::iCharge(PDT::Charge ci) {
theCharge = ci;
if ( synchronized() && CC() ) CC()->theCharge = PDT::Charge(-ci);
return theCharge;
}
PDT::Spin ParticleData::iSpin(PDT::Spin si) {
theSpin = si;
if ( synchronized() && CC() ) CC()->theSpin = si;
return si;
}
PDT::Colour ParticleData::iColour(PDT::Colour ci) {
theColour = ci;
if ( synchronized() && CC() ) CC()->theColour = PDT::Colour(-ci);
return theColour;
}
void ParticleData::stable(bool s) {
isStable = s;
if ( synchronized() && CC() ) CC()->isStable = s;
}
void ParticleData::synchronized(bool h) {
syncAnti = h;
if ( CC() ) CC()->syncAnti = h;
}
void ParticleData::variableRatio(bool varRatio) {
theVariableRatio=varRatio;
}
void ParticleData::addDecayMode(tDMPtr dm) {
if ( member(theDecayModes, dm) ) return;
cPDPtr parent = dm->parent();
if ( !parent ) parent = this;
if ( parent != this ) {
dm = dm->clone(this);
}
theDecayModes.insert(dm);
theDecaySelector.insert(dm->brat(), dm);
if ( CC() ) {
if ( !synchronized() ) dm->CC()->switchOff();
CC()->theDecayModes.insert(dm->CC());
CC()->theDecaySelector.insert(dm->CC()->brat(), dm->CC());
}
}
void ParticleData::removeDecayMode(tDMPtr dm) {
theDecayModes.erase(theDecayModes.find(dm));
if(theDecayModes.empty()) isStable = true;
theDecaySelector.erase(dm);
if ( !CC() ) return;
CC()->theDecayModes.erase(dm->CC());
if(CC()->theDecayModes.empty()) CC()->isStable = true;
CC()->theDecaySelector.erase(dm->CC());
}
void ParticleData::synchronize() {
if ( !CC() ) return;
isStable = CC()->isStable;
theMass = CC()->theMass;
theHardProcessMass = CC()->theHardProcessMass;
hardProcessMassSet = CC()->hardProcessMassSet;
theWidth = CC()->theWidth;
theHardProcessWidth = CC()->theHardProcessWidth;
hardProcessWidthSet = CC()->hardProcessWidthSet;
theWidthUpCut = CC()->theWidthUpCut;
theWidthLoCut = CC()->theWidthLoCut;
theCTau = CC()->theCTau;
theCharge = PDT::Charge(-CC()->theCharge);
theSpin = CC()->theSpin;
theColour = PDT::antiColour(CC()->theColour);
theMassGenerator = CC()->theMassGenerator;
theWidthGenerator = CC()->theWidthGenerator;
syncAnti = CC()->syncAnti;
theDecaySelector.clear();
for ( DecaySet::iterator it = theDecayModes.begin();
it != theDecayModes.end(); ++it ) {
(*it)->synchronize();
theDecaySelector.insert((*it)->brat(), *it);
}
}
void ParticleData::doupdate() {
Interfaced::doupdate();
bool redo = touched();
for_each(theDecayModes, UpdateChecker(redo));
UpdateChecker::check(theMassGenerator, redo);
UpdateChecker::check(theWidthGenerator, redo);
if ( !redo ) return;
theDecaySelector.clear();
for ( DecaySet::const_iterator dit = theDecayModes.begin();
dit != theDecayModes.end(); ++dit ) {
tDMPtr dm = *dit;
dm->resetOverlap();
for ( DecaySet::const_iterator dit2 = theDecayModes.begin();
dit2 != theDecayModes.end(); ++dit2 )
if ( dit2 != dit ) dm->addOverlap(dm);
if ( dm->brat() > 0.0 ) theDecaySelector.insert(dm->brat(), dm);
}
if ( theMassGenerator && !theMassGenerator->accept(*this) )
throw UpdateException();
if ( theWidthGenerator &&
!theWidthGenerator->accept(*this) )
throw UpdateException();
if ( theWidthGenerator ) theDecaySelector = theWidthGenerator->rate(*this);
touch();
}
tDMPtr ParticleData::selectMode(Particle & p) const {
if ( &(p.data()) != this ) return tDMPtr();
try {
if ( !theWidthGenerator || !theVariableRatio )
return theDecaySelector.select(UseRandom::current());
DecaySelector local;
if ( theWidthGenerator )
local = theWidthGenerator->rate(p);
else
for ( DecaySet::const_iterator mit = theDecayModes.begin();
mit != theDecayModes.end(); ++mit )
local.insert((*mit)->brat(p), *mit);
return local.select(UseRandom::current());
}
- catch (range_error) {
+ catch (range_error & e) {
return tDMPtr();
}
}
void ParticleData::rebind(const TranslationMap & trans) {
if ( CC() ) theAntiPartner = trans.translate(theAntiPartner);
DecaySet newModes;
DecaySelector newSelector;
for ( DecaySet::iterator it = theDecayModes.begin();
it != theDecayModes.end(); ++it ) {
DMPtr dm;
dm = trans.translate(*it);
if ( !dm ) throw RebindException();
newModes.insert(dm);
newSelector.insert(dm->brat(), dm);
}
theDecayModes.swap(newModes);
theDecaySelector.swap(newSelector);
}
IVector ParticleData::getReferences() {
IVector refs = Interfaced::getReferences();
if ( CC() ) refs.push_back(CC());
refs.insert(refs.end(), theDecayModes.begin(), theDecayModes.end());
return refs;
}
void ParticleData::massGenerator(tMassGenPtr mg) {
if ( mg && !mg->accept(*this) ) return;
if ( mg && synchronized() && CC() && !mg->accept(*CC()) ) return;
theMassGenerator = mg;
if ( synchronized() && CC() ) CC()->theMassGenerator = mg;
}
void ParticleData::widthGenerator(tWidthGeneratorPtr newGen) {
if ( newGen && !newGen->accept(*this) ) return;
if ( newGen && synchronized() && CC() && !newGen->accept(*CC()) ) return;
theWidthGenerator = newGen;
if ( synchronized() && CC() ) CC()->theWidthGenerator = newGen;
}
Energy ParticleData::generateMass() const {
return massGenerator()? massGenerator()->mass(*this): mass();
}
Energy ParticleData::generateWidth(Energy m) const {
return widthGenerator()? widthGenerator()->width(*this, m): width();
}
Length ParticleData::generateLifeTime(Energy m, Energy w) const {
return widthGenerator() ?
widthGenerator()->lifeTime(*this, m, w) :
UseRandom::rndExp(cTau());
}
PPtr ParticleData::produceParticle(const Lorentz5Momentum & pp) const {
PPtr p = new_ptr(Particle(this));
p->set5Momentum(pp);
return p;
}
PPtr ParticleData::produceParticle(const LorentzMomentum & pp) const {
PPtr p(produceParticle(Lorentz5Momentum(pp)));
return p;
}
PPtr ParticleData::produceParticle(const LorentzMomentum & pp, Energy m) const {
PPtr p(produceParticle(Lorentz5Momentum(pp, m)));
return p;
}
PPtr ParticleData::produceParticle(Energy m, const Momentum3 & pp) const {
PPtr p(produceParticle(Lorentz5Momentum(m, pp)));
return p;
}
PPtr ParticleData::produceParticle(const Momentum3 & pp) const {
PPtr p(produceParticle(Lorentz5Momentum(generateMass(), pp)));
return p;
}
PPtr ParticleData::
produceParticle(Energy plus, Energy minus, Energy px, Energy py) const {
PPtr p(produceParticle(LorentzMomentum(px, py, 0.5*(plus-minus),
0.5*(plus+minus))));
return p;
}
void ParticleData::setMass(Energy mi) {
mass(mi);
}
void ParticleData::setHardProcessMass(Energy mi) {
theHardProcessMass = mi;
hardProcessMassSet = true;
ParticleData * apd = CC().operator->();
if ( synchronized() && apd ) {
apd->theHardProcessMass = theHardProcessMass;
apd->hardProcessMassSet = true;
}
}
void ParticleData::setHardProcessWidth(Energy mi) {
theHardProcessWidth = mi;
hardProcessWidthSet = true;
ParticleData * apd = CC().operator->();
if ( synchronized() && apd ) {
apd->theHardProcessWidth = theHardProcessWidth;
apd->hardProcessWidthSet = true;
}
}
string ParticleData::doUnsetHardProcessMass(string) {
hardProcessMassSet = false;
theHardProcessMass = -1.*GeV;
return "";
}
string ParticleData::doAdjustNominalMass(string) {
if ( hardProcessMassSet )
setMass(theHardProcessMass);
return "";
}
string ParticleData::doUnsetHardProcessWidth(string) {
hardProcessWidthSet = false;
theHardProcessWidth = -1.*GeV;
return "";
}
Energy ParticleData::defMass() const {
return theDefMass;
}
void ParticleData::setWidth(Energy wi) {
width(wi);
}
Energy ParticleData::getWidth() const {
return width();
}
Energy ParticleData::defWidth() const {
return theDefWidth;
}
void ParticleData::setCut(Energy ci) {
widthCut(ci);
}
Energy ParticleData::getCut() const {
return (theWidthUpCut >= ZERO && theWidthLoCut >= ZERO)?
max(theWidthUpCut, theWidthLoCut): min(theWidthUpCut, theWidthLoCut);
}
Energy ParticleData::defCut() const {
return theDefCut;
}
void ParticleData::setUpCut(Energy ci) {
widthUpCut(ci);
}
Energy ParticleData::getUpCut() const {
return theWidthUpCut;
}
void ParticleData::setLoCut(Energy ci) {
widthLoCut(ci);
}
Energy ParticleData::getLoCut() const {
return theWidthLoCut;
}
void ParticleData::setCTau(Length ti) {
cTau(ti);
}
Length ParticleData::getCTau() const {
return cTau();
}
Length ParticleData::defCTau() const {
return theDefCTau;
}
void ParticleData::setStable(long is) {
stable(is);
}
long ParticleData::getStable() const {
return stable();
}
void ParticleData::setSync(long is) {
synchronized(is);
}
long ParticleData::getSync() const {
return synchronized();
}
void ParticleData::setVariableRatio(long is) {
variableRatio(is);
}
long ParticleData::getVariableRatio() const {
return variableRatio();
}
string ParticleData::doSync(string) {
synchronize();
return "";
}
void ParticleData::setMassGenerator(MassGenPtr gi) {
massGenerator(gi);
}
void ParticleData::setWidthGenerator(WidthGeneratorPtr wg) {
widthGenerator(wg);
}
void ParticleData::setColour(long c) {
theColour = PDT::Colour(c);
}
long ParticleData::getColour() const {
return theColour;
}
long ParticleData::defColour() const {
return theDefColour;
}
void ParticleData::setCharge(int c) {
theCharge = PDT::Charge(c);
}
string ParticleData::ssetCharge(string arg) {
istringstream is(arg);
long i;
if ( is >> i ) {
theCharge = PDT::Charge(i);
return "New charge is " + arg;
}
if ( arg == "unknown" )
theCharge = PDT::ChargeUnknown;
else if ( arg == "charged" )
theCharge = PDT::Charged;
else if ( arg == "positive" )
theCharge = PDT::Positive;
else if ( arg == "negative" )
theCharge = PDT::Negative;
else throw ParticleChargeCommand(*this, arg);
return "New charge is " + arg;
}
int ParticleData::getCharge() const {
return theCharge;
}
int ParticleData::defCharge() const {
return theDefCharge;
}
void ParticleData::setSpin(int s) {
theSpin = PDT::Spin(s);
}
int ParticleData::getSpin() const {
return theSpin;
}
int ParticleData::defSpin() const {
return theDefSpin;
}
ClassDescription<ParticleData> ParticleData::initParticleData;
struct ParticleOrdering {
bool 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() );
}
};
struct ModeOrdering {
bool operator()(const tcDMPtr & d1, const tcDMPtr & d2) {
ParticleOrdering ord;
return ord(d1->parent(), d2->parent()) ||
( !ord(d2->parent(), d1->parent()) &&
( d1->tag() < d2->tag() ||
( d1->tag() == d2->tag() && d1->fullName() < d2->fullName() ) ) );
}
};
void ParticleData::persistentOutput(PersistentOStream & os) const {
multiset<tcDMPtr,ModeOrdering>
modes(theDecayModes.begin(), theDecayModes.end());
os << long(theId) << thePDGName << ounit(theMass, GeV) << ounit(theWidth, GeV)
<< ounit(theHardProcessMass,GeV) << hardProcessMassSet
<< ounit(theHardProcessWidth,GeV) << hardProcessWidthSet
<< ounit(theWidthUpCut, GeV) << ounit(theWidthLoCut, GeV)
<< ounit(theCTau, mm) << oenum(theCharge) << oenum(theSpin)
<< oenum(theColour);
os << theMassGenerator << isStable << modes << theDecaySelector
<< theWidthGenerator << theVariableRatio << theAntiPartner << syncAnti
<< ounit(theDefMass, GeV) << ounit(theDefWidth, GeV)
<< ounit(theDefCut, GeV) << ounit(theDefCTau, mm) << oenum(theDefColour)
<< oenum(theDefCharge) << oenum(theDefSpin);
}
void ParticleData::persistentInput(PersistentIStream & is, int) {
long id;
is >> id >> thePDGName >> iunit(theMass, GeV) >> iunit(theWidth, GeV)
>> iunit(theHardProcessMass,GeV) >> hardProcessMassSet
>> iunit(theHardProcessWidth,GeV) >> hardProcessWidthSet
>> iunit(theWidthUpCut, GeV) >> iunit(theWidthLoCut, GeV)
>> iunit(theCTau, mm) >> ienum(theCharge) >> ienum(theSpin)
>> ienum(theColour) >> theMassGenerator >> isStable
>> theDecayModes >> theDecaySelector >> theWidthGenerator >> theVariableRatio
>> theAntiPartner >> syncAnti >> iunit(theDefMass, GeV)
>> iunit(theDefWidth, GeV) >> iunit(theDefCut, GeV)
>> iunit(theDefCTau, mm) >> ienum(theDefColour) >> ienum(theDefCharge)
>> ienum(theDefSpin);
theId = id;
}
void ParticleData::Init() {
static ClassDocumentation<ParticleData> documentation
("There is no documentation for the ThePEG::ParticleData class");
static Parameter<ParticleData,Energy> interfaceMass
("NominalMass",
"The nominal mass in GeV of the particle. The actual mass "
"of a particle instance is generated depending on the "
"nominal mass and the width and is generated by the "
"<interface>Mass_generator</interface> object associated with the "
"particle.",
&ParticleData::theMass, GeV, ZERO, ZERO, Constants::MaxEnergy,
false, false, Interface::lowerlim,
&ParticleData::setMass, 0, 0, 0, &ParticleData::defMass);
static Parameter<ParticleData,Energy> interfaceHardProcessMass
("HardProcessMass",
"The mass in GeV of the particle to be used in calculating hard process cross sections.",
&ParticleData::theHardProcessMass, GeV, ZERO, ZERO, Constants::MaxEnergy,
false, false, Interface::lowerlim,
&ParticleData::setHardProcessMass, 0, 0, 0, 0);
static Parameter<ParticleData,Energy> interfaceDefMass
("DefaultMass",
"The default nominal mass in GeV of the particle. The actual mass "
"of a particle instance is generated depending on the "
"nominal mass and the width and is generated by the "
"<interface>Mass_generator</interface> object associated with the "
"particle.",
&ParticleData::theDefMass, GeV, ZERO, ZERO, Constants::MaxEnergy,
false, true, Interface::lowerlim);
interfaceDefMass.setHasDefault(false);
static Parameter<ParticleData,Energy> interfaceWidth
("Width",
"The width of the particle in GeV.",
0, GeV, ZERO, ZERO, Constants::MaxEnergy,
false, false, Interface::lowerlim,
&ParticleData::setWidth, &ParticleData::getWidth,
0, 0, &ParticleData::defWidth);
static Parameter<ParticleData,Energy> interfaceHardProcessWidth
("HardProcessWidth",
"The width in GeV of the particle to be used in calculating hard process cross sections.",
&ParticleData::theHardProcessWidth, GeV, ZERO, ZERO, Constants::MaxEnergy,
false, false, Interface::lowerlim,
&ParticleData::setHardProcessWidth, 0, 0, 0, 0);
static Parameter<ParticleData,Energy> interfaceDefWidth
("DefaultWidth",
"The default width of the particle in GeV.",
&ParticleData::theDefWidth, GeV, ZERO, ZERO, Constants::MaxEnergy,
false, true, Interface::lowerlim);
interfaceDefWidth.setHasDefault(false);
static Parameter<ParticleData,Energy> interfaceWidthUpCut
("WidthUpCut",
"The upper hard cutoff in GeV in generated mass, which is the maximum "
"allowed upwards deviation from the nominal mass. A negative value "
"corresponds to no limit.",
0, GeV, ZERO, -1.0*GeV, Constants::MaxEnergy,
false, false, Interface::lowerlim,
&ParticleData::setUpCut, &ParticleData::getUpCut,
0, 0, &ParticleData::defCut);
static Parameter<ParticleData,Energy> interfaceWidthLoCut
("WidthLoCut",
"The lower hard cutoff in GeV in generated mass, which is the maximum "
"allowed downwards deviation from the nominal mass. A negative value "
"corresponds to no limit.",
0, GeV, ZERO, -1.0*GeV, Constants::MaxEnergy,
false, false, Interface::lowerlim,
&ParticleData::setLoCut, &ParticleData::getLoCut,
0, 0, &ParticleData::defCut);
static Parameter<ParticleData,Energy> interfaceWidthCut
("WidthCut",
"The hard cutoff in GeV in generated mass, which is the maximum "
"allowed deviation from the nominal mass. Sets both the upper and lower "
"cut. (The displayed value is the maximium of the upper and lower cut.) "
"A negative value corresponds to no limit.",
0, GeV, ZERO, -1.0*GeV, Constants::MaxEnergy,
false, false, Interface::lowerlim,
&ParticleData::setCut, &ParticleData::getCut,
0, 0, &ParticleData::defCut);
interfaceWidthCut.setHasDefault(false);
static Parameter<ParticleData,Energy> interfaceDefWidthCut
("DefaultWidthCut",
"The default hard cutoff in GeV in generated mass, which is the maximum "
"allowed deviation from the nominal mass. For the actual cutoff, the "
"upper and lower cut can be set separately.",
&ParticleData::theDefCut, GeV, ZERO, ZERO, Constants::MaxEnergy,
false, true, Interface::lowerlim);
interfaceDefWidthCut.setHasDefault(false);
static Parameter<ParticleData,Length> interfaceCTau
("LifeTime",
"c times the average lifetime of the particle measuerd in mm."
"The actual lifetime of a particle instance is generated "
"from this number by the <interface>Mass_generator</interface> "
"object associated with the particle.",
0, mm, ZERO, ZERO, Constants::MaxLength,
false, false, Interface::lowerlim,
&ParticleData::setCTau, &ParticleData::getCTau,
0, 0, &ParticleData::defCTau);
interfaceCTau.setHasDefault(false);
static Parameter<ParticleData,Length> interfaceDefCTau
("DefaultLifeTime",
"c times the default average lifetime of the particle measuerd in mm."
"The actual lifetime of a particle instance is generated "
"from this number by the <interface>Mass_generator</interface> "
"object associated with the particle.",
&ParticleData::theDefCTau, mm, ZERO, ZERO, Constants::MaxLength,
false, true, Interface::lowerlim);
interfaceDefCTau.setHasDefault(false);
static Switch<ParticleData> interfaceColour
("Colour",
"The colour quantum number of this particle type.",
0, -1, false, false, &ParticleData::setColour, &ParticleData::getColour,
&ParticleData::defColour);
static SwitchOption interfaceColourUndefined
(interfaceColour, "Undefined", "The coulur is undefined.", -1);
static SwitchOption interfaceColourNeutral
(interfaceColour, "Neutral", "This particle is colour neutral.", 0);
static SwitchOption interfaceColour3
(interfaceColour, "Triplet", "This particle is a colour triplet.", 3);
static SwitchOption interfaceColour3bar
(interfaceColour, "AntiTriplet",
"This particle is a colour anti-triplet.", -3);
static SwitchOption interfaceColour6
(interfaceColour, "Sextet", "This particle is a colour sextet.", 6);
static SwitchOption interfaceColour6bar
(interfaceColour, "AntiSextet",
"This particle is a colour anti-sextet.", -6);
static SwitchOption interfaceColour8
(interfaceColour, "Octet", "This particle is a colour octet.", 8);
static Switch<ParticleData,PDT::Colour> interfaceDefColour
("DefaultColour",
"The default colour quantum number of this particle type.",
&ParticleData::theDefColour, PDT::Colour(-1), false, true);
static SwitchOption interfaceDefColourUndefined
(interfaceDefColour, "Undefined", "The coulur is undefined.", -1);
static SwitchOption interfaceDefColourNeutral
(interfaceDefColour, "Neutral", "This particle is colour neutral.", 0);
static SwitchOption interfaceDefColour3
(interfaceDefColour, "Triplet", "This particle is a colour triplet.", 3);
static SwitchOption interfaceDefColour3bar
(interfaceDefColour, "AntiTriplet",
"This particle is a colour anti-triplet.", -3);
static SwitchOption interfaceDefColour6
(interfaceDefColour, "Sextet", "This particle is a colour sextet.", 6);
static SwitchOption interfaceDefColour6bar
(interfaceDefColour, "AntiSextet",
"This particle is a colour anti-sextet.", -6);
static SwitchOption interfaceDefColour8
(interfaceDefColour, "Octet", "This particle is a colour octet.", 8);
interfaceDefColour.setHasDefault(false);
static Parameter<ParticleData, int> interfaceCharge
("Charge",
"The charge of this particle in units of e/3. "
"See also the command interface <interface>SetCharge</interface>.",
0, 0, -24, 24, false, false, true,
&ParticleData::setCharge, &ParticleData::getCharge, 0, 0,
&ParticleData::defCharge);
static Parameter<ParticleData, PDT::Charge> interfaceDefCharge
("DefaultCharge",
"The default charge of this particle in units of e/3. "
"See also the command interface <interface>SetCharge</interface>.",
&ParticleData::theDefCharge, PDT::Charge(0), PDT::Charge(-24),
PDT::Charge(24), false, true, true);
interfaceDefCharge.setHasDefault(false);
static Command<ParticleData> interfaceSetCharge
("SetCharge",
"Set the charge of this particle. The argument should be given as an "
"interger giving three times the unit charge, or 'unknown', "
"'charged', 'positive' or 'negative'", &ParticleData::ssetCharge);
static Parameter<ParticleData, int> interfaceSpin
("Spin",
"The spin quantim number of this particle on the form 2j+1.",
0, 0, 0, 9, false, false, true,
&ParticleData::setSpin, &ParticleData::getSpin, 0, 0,
&ParticleData::defSpin);
static Parameter<ParticleData, PDT::Spin> interfaceDefSpin
("DefaultSpin",
"The default spin quantim number of this particle on the form 2j+1.",
&ParticleData::theDefSpin, PDT::Spin(0), PDT::Spin(0), PDT::Spin(9),
false, true, true);
interfaceDefSpin.setHasDefault(false);
static Switch<ParticleData> interfaceStable
("Stable",
"Indicates if the particle is stable or not.",
0, 0, false, false,
&ParticleData::setStable, &ParticleData::getStable, 0);
static SwitchOption interfaceStableYes
(interfaceStable,
"Stable",
"This particle is stable",
1);
static SwitchOption interfaceStableNo
(interfaceStable,
"Unstable",
"This particle is not stable",
0);
interfaceStable.setHasDefault(false);
static Switch<ParticleData> interfaceVariableRatio
("VariableRatio",
"Indicates if the branching ratios of the particle are allowed"
" to vary for given Particle instances depending on the mass of the instance.",
0, 0, false, false,
&ParticleData::setVariableRatio, &ParticleData::getVariableRatio, 0);
static SwitchOption interfaceVariableRatioYes
(interfaceVariableRatio,
"Yes",
"The branching ratio varies.",
1);
static SwitchOption interfaceVariableRatioNo
(interfaceVariableRatio,
"No",
"The branching ratio does not vary.",
0);
static Switch<ParticleData> interfaceSync
("Synchronized",
"Indicates if the changes to this particle is propagated to "
"its anti-partner or not. Note that setting this switch does not "
"actually synchronize the properties with the anti-partner, "
"it only assures that following changes are propagated. "
"To sync the particle with its anti-particle, use the "
"<interface>Synchronize</interface> command.",
0, 1, false, false,
&ParticleData::setSync, &ParticleData::getSync, 0);
static SwitchOption interfaceSyncYes
(interfaceSync,
"Synchronized",
"Changes to this particle will propagate to its "
"anti-partner",
1);
static SwitchOption interfaceSyncNo
(interfaceSync,
"Not_synchronized",
"Changes to this particle will propagate to its "
"anti-partner",
0);
interfaceSync.setHasDefault(false);
static Command<ParticleData> interfaceSynchronize
("Synchronize",
"Synchronizes this particle so that all its properties "
"correspond to those of its anti-partner",
&ParticleData::doSync, false);
static Reference<ParticleData,MassGenerator> interfaceMassGenerator
("Mass_generator",
"An object derived from the ThePEG::MassGenerator"
"class, which is able to generate a mass for a given "
"particle instance",
&ParticleData::theMassGenerator, false, false, true, true,
&ParticleData::setMassGenerator, 0, 0);
static Reference<ParticleData,WidthGenerator> interfaceWidthGenerator
("Width_generator",
"An object derived from the ThePEG::WidthGenerator class, "
"which is able to calculate the full and partial widths for"
"this particle type and for a given instance of this "
"particle type.",
&ParticleData::theWidthGenerator, false, false, true, true,
&ParticleData::setWidthGenerator, 0, 0);
static RefVector<ParticleData,DecayMode> interfaceDecayModes
("DecayModes",
"The list of decay modes defined for this particle type.",
0, -1, false, false, false, false, 0,
&ParticleData::insDecayModes, &ParticleData::delDecayModes,
&ParticleData::getDecayModes);
static Command<ParticleData> interfaceSelectDecayModes
("SelectDecayModes",
"Only the decay modes which are given as (white-space separated) "
"decay tags will be switched on, all others will be switched off. "
"If no argument or 'none' is given, all decay modes are switched off. "
"If the argument is 'all', all decay modes are switched on.",
&ParticleData::doSelectDecayModes, false);
static Command<ParticleData> interfacePrintDecayModes
("PrintDecayModes",
"Print all decay modes of this particle.",
&ParticleData::doPrintDecayModes, true);
static Command<ParticleData> interfaceUnsetHardProcessMass
("UnsetHardProcessMass",
"Unset a previously set hard process mass.",
&ParticleData::doUnsetHardProcessMass, false);
static Command<ParticleData> interfaceAdjustNominalMass
("AdjustNominalMass",
"Unset a previously set hard process mass.",
&ParticleData::doAdjustNominalMass, false);
static Command<ParticleData> interfaceUnsetHardProcessWidth
("UnsetHardProcessWidth",
"Unset a previously set hard process width.",
&ParticleData::doUnsetHardProcessWidth, false);
interfaceStable.rank(14);
interfaceDecayModes.rank(13);
interfaceMass.rank(12);
interfaceWidth.rank(11);
interfaceWidthCut.rank(10);
interfaceCTau.rank(9);
interfaceMassGenerator.rank(8);
interfaceWidthGenerator.rank(7);
interfaceWidthUpCut.rank(-0.1);
interfaceWidthLoCut.rank(-0.1);
}
string ParticleData::doPrintDecayModes(string) {
multimap<double,tDMPtr, std::greater<double> > sorted;
for ( DecaySet::iterator it = decayModes().begin();
it != decayModes().end(); ++it )
sorted.insert(make_pair((**it).brat(), *it));
ostringstream os;
for ( multimap<double,tDMPtr,
std::greater<double> >::iterator it = sorted.begin();
it != sorted.end(); ++it )
os << it->second->tag()
<< (it->second->on()? " ": " (off) ")
<< it->first << endl;
return os.str();
}
string ParticleData::doSelectDecayModes(string args) {
DecaySet on;
while ( !args.empty() ) {
string arg = StringUtils::car(args);
if ( arg == "all" ) {
on = decayModes();
break;
}
if ( arg == "none" ) {
on.clear();
break;
}
string name = arg;
args = StringUtils::cdr(args);
if ( arg.empty() ) continue;
if ( arg[0] != '/' ) arg = fullName() + "/" + arg;
DMPtr dm = Repository::GetPtr<DMPtr>(arg);
if ( !dm ) return "Error: No decay mode with tag '" + name + "' exists.";
on.insert(dm);
}
for ( DecaySet::iterator it = decayModes().begin();
it != decayModes().end(); ++it ) {
if ( on.find(*it) != on.end() ) {
(**it).switchOn();
on.erase(*it);
} else {
(**it).switchOff();
}
}
if ( !on.empty() )
return "Error: decay mode '" + (**on.begin()).tag() + "'was not available.";
return "";
}
void ParticleData::insDecayModes(DMPtr dm, int) {
addDecayMode(dm);
}
void ParticleData::delDecayModes(int i) {
vector<DMPtr> mv = getDecayModes();
if ( i >= 0 && static_cast<unsigned int>(i) < mv.size() ) removeDecayMode(mv[i]);
}
vector<DMPtr> ParticleData::getDecayModes() const {
return vector<DMPtr>(theDecayModes.begin(), theDecayModes.end());
}
ParticleChargeCommand::
ParticleChargeCommand(const ParticleData & pd, string arg) {
theMessage << "Cannot set the charge of particle '" << pd.name()
<< "' to '" << arg << "'.";
severity(warning);
}
void ParticleData::doinit() {
Interfaced::doinit();
if( theMassGenerator ) theMassGenerator->init();
if( theWidthGenerator ) theWidthGenerator->init();
}
void ParticleData::doinitrun() {
Interfaced::doinitrun();
if( theMassGenerator ) theMassGenerator->initrun();
if( theWidthGenerator ) theWidthGenerator->initrun();
}
}
diff --git a/PDT/QuarksToHadronsDecayer.cc b/PDT/QuarksToHadronsDecayer.cc
--- a/PDT/QuarksToHadronsDecayer.cc
+++ b/PDT/QuarksToHadronsDecayer.cc
@@ -1,235 +1,235 @@
// -*- C++ -*-
//
// QuarksToHadronsDecayer.cc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2017 Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the QuarksToHadronsDecayer class.
//
#include "QuarksToHadronsDecayer.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/PDT/DecayMode.h"
#include "ThePEG/PDT/StandardMatchers.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/SimplePhaseSpace.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace ThePEG;
QuarksToHadronsDecayer::~QuarksToHadronsDecayer() {}
IBPtr QuarksToHadronsDecayer::clone() const {
return new_ptr(*this);
}
IBPtr QuarksToHadronsDecayer::fullclone() const {
return new_ptr(*this);
}
bool QuarksToHadronsDecayer::accept(const DecayMode & dm) const {
int col = 0;
int acol = 0;
if ( !dm.productMatchers().empty() ) {
for ( MatcherMSet::const_iterator it = dm.productMatchers().begin();
it != dm.productMatchers().end(); ++it ) {
if ( typeid(**it) == typeid(MatchLightQuark) ) ++col;
else if ( typeid(**it) == typeid(MatchLightAntiQuark) ) ++acol;
else return false;
}
if ( col != 1 || col != acol ) return false;
}
if ( dm.orderedProducts().size() + col + acol < 2 ||
!dm.cascadeProducts().empty() || dm.wildProductMatcher() ) return false;
for ( int i = 0, N = dm.orderedProducts().size(); i < N; ++i ) {
if ( DiquarkMatcher::Check(*dm.orderedProducts()[i]) ) {
if ( i + 1 != N ) return false;
if ( dm.orderedProducts()[i]->id() < 0 ) ++col;
else ++acol;
}
if ( QuarkMatcher::Check(*dm.orderedProducts()[i]) ) {
if ( dm.orderedProducts()[i]->id() > 0 ) ++col;
else ++acol;
}
}
if ( acol != col || col < 1 || col > 2 ) return false;
return true;
}
PVector QuarksToHadronsDecayer::decay(const DecayMode & dm,
const Particle & parent) const {
PVector children;
tcPDVector quarks;
if ( !dm.productMatchers().empty() ) {
tcPDPtr pd = getParticleData(flavourGenerator()->selectQuark());
quarks.push_back(pd);
quarks.push_back(pd->CC());
}
Energy summq = ZERO;
Energy summp = ZERO;
tPDVector prods = dm.orderedProducts();
for ( int i = 0, N = prods.size(); i < N; ++i )
if ( QuarkMatcher::Check(*prods[i]) || DiquarkMatcher::Check(*prods[i])) {
quarks.push_back(prods[i]);
summq += quarks.back()->mass();
} else {
children.push_back(prods[i]->produceParticle());
summp += children.back()->mass();
}
Energy summh = ZERO;
PVector hadrons;
if ( !quarks.empty() ) do {
hadrons = getHadrons(getN(parent.mass(), summq, quarks.size()), quarks);
summh = ZERO;
for ( int i = 0, N = hadrons.size(); i < N; ++i )
summh += hadrons[i]->mass();
} while ( hadrons.empty() || summp + summh >= parent.mass() );
children.insert(children.end(), hadrons.begin(), hadrons.end());
distribute(parent, children);
finalBoost(parent, children);
setScales(parent, children);
return children;
}
int QuarksToHadronsDecayer::getN(Energy m0, Energy summq, int Nq) const {
int Nh = fixedN();
if ( Nh >= 2 ) return Nh;
double c = c1()*log((m0 - summq)/c2()) + c3();
if ( c < 0.0 ) return minN();
while ( true ) {
using namespace Constants;
Nh = int(0.5 + double(Nq)/4.0 + c +
sqrt(-2.0*c*log(max(1.0e-10, rnd())))*sin(2.0*pi*rnd()));
if ( Nh >= minN() ) return Nh;
}
}
PVector QuarksToHadronsDecayer::
getHadrons(int Nh, tcPDVector quarks) const {
PVector hadrons;
Nh -= quarks.size()/2;
while ( Nh-- > 0 ) {
int i = irnd(quarks.size() - 1);
tcPDPair hq = flavourGenerator()->alwaysGenerateHadron(quarks[i]);
hadrons.push_back(hq.first->produceParticle());
quarks[i] = hq.second;
}
if ( DiquarkMatcher::Check(*quarks[0]) && DiquarkMatcher::Check(*quarks[1]) )
return PVector();
tcPDPtr h = flavourGenerator()->alwaysGetHadron(quarks[0], quarks[1]);
hadrons.push_back(h->produceParticle());
if ( quarks.size() <= 2 ) return hadrons;
if ( DiquarkMatcher::Check(*quarks[2]) && DiquarkMatcher::Check(*quarks[3]) )
return PVector();
h = flavourGenerator()->alwaysGetHadron(quarks[2], quarks[3]);
hadrons.push_back(h->produceParticle());
return hadrons;
}
void QuarksToHadronsDecayer::
distribute(const Particle & parent, PVector & children) const {
do {
try {
SimplePhaseSpace::CMSn(children, parent.mass());
}
- catch ( ImpossibleKinematics ) {
+ catch ( ImpossibleKinematics & e) {
children.clear();
return;
}
} while ( reweight(parent, children) < rnd() );
}
double QuarksToHadronsDecayer::
reweight(const Particle &, const PVector &) const {
return 1.0;
}
void QuarksToHadronsDecayer::persistentOutput(PersistentOStream & os) const {
os << theFixedN << theMinN << theC1 << ounit(theC2,GeV) << theC3 << theFlavourGenerator;
}
void QuarksToHadronsDecayer::persistentInput(PersistentIStream & is, int) {
is >> theFixedN >> theMinN >> theC1 >> iunit(theC2,GeV) >> theC3 >> theFlavourGenerator;
}
ClassDescription<QuarksToHadronsDecayer> QuarksToHadronsDecayer::initQuarksToHadronsDecayer;
// Definition of the static class description member.
void QuarksToHadronsDecayer::Init() {
static ClassDocumentation<QuarksToHadronsDecayer> documentation
("This class decays particles to nq (2 or 4) quarks which then are "
"decayes to hadrons according to phase space. The number of final "
"hadrons can either be given by a fixed number or as a Gaussian "
"multiplicity distribution centered around c+nq/4+c3 and a width "
"sqrt(c), where c = c1 log((m - summ)/c2), m is the mass of the "
"decaying particle, summ the sum of the quark masses and ci real "
"parameters.");
static Parameter<QuarksToHadronsDecayer,int> interfaceFixedN
("FixedN",
"The fixed number of hadrons to be produced. If less than 2, the "
"number is instead given by a gaussian multiplicity distribution.",
&QuarksToHadronsDecayer::theFixedN, 0, 0, 10,
true, false, true);
static Parameter<QuarksToHadronsDecayer,int> interfaceMinN
("MinN",
"The minimum hadrons to be produced.",
&QuarksToHadronsDecayer::theMinN, 2, 2, 10,
true, false, true);
static Parameter<QuarksToHadronsDecayer,double> interfaceC1
("C1",
"The c1 parameter of the gaussian multiplicity distribution centered "
"around c1 log((m - summ)/c2) +c3.",
&QuarksToHadronsDecayer::theC1, 4.5, 0.0, 10.0,
true, false, true);
static Parameter<QuarksToHadronsDecayer,Energy> interfaceC2
("C2",
"The c2 parameter of the gaussian multiplicity distribution centered "
"around c1 log((m - summ)/c2) +c3.",
&QuarksToHadronsDecayer::theC2, GeV, 0.7*GeV, ZERO, 10.0*GeV,
true, false, true);
static Parameter<QuarksToHadronsDecayer,double> interfaceC3
("C3",
"The c3 parameter of the gaussian multiplicity distribution centered "
"around c1 log((m - summ)/c2) +c3.",
&QuarksToHadronsDecayer::theC3, 0.0, 0.0, 10.0,
true, false, true);
static Reference<QuarksToHadronsDecayer,FlavourGenerator>
interfaceFlavourGenerator
("FlavourGenerator",
"The object in charge of generating hadrons spieces from given quark "
"flavours.",
&QuarksToHadronsDecayer::theFlavourGenerator,
true, false, true, false, true);
interfaceFixedN.rank(10);
interfaceMinN.rank(9);
interfaceFlavourGenerator.rank(8);
interfaceMinN.setHasDefault(false);;
}
diff --git a/Repository/BaseRepository.cc b/Repository/BaseRepository.cc
--- a/Repository/BaseRepository.cc
+++ b/Repository/BaseRepository.cc
@@ -1,978 +1,977 @@
// -*- C++ -*-
//
// BaseRepository.cc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2017 Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the BaseRepository class.
//
// macro is passed in from -D compile flag
#ifndef THEPEG_PKGDATADIR
#error Makefile.am needs to define THEPEG_PKGDATADIR
#endif
#include "BaseRepository.h"
#include "ThePEG/Config/algorithm.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/InterfaceBase.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Command.h"
#include "ThePEG/Utilities/ClassDescription.h"
#include "ThePEG/Utilities/DescriptionList.h"
#include "ThePEG/Utilities/HoldFlag.h"
#include "ThePEG/Utilities/TypeInfo.h"
#include "ThePEG/Utilities/DynamicLoader.h"
#include "ThePEG/Utilities/StringUtils.h"
#include "ThePEG/Utilities/Throw.h"
#include "ThePEG/PDT/DecayMode.h"
#ifdef ThePEG_TEMPLATES_IN_CC_FILE
#include "BaseRepository.tcc"
#endif
using namespace ThePEG;
ostream *& BaseRepository::coutp() {
static ostream * theCout = &std::cout;
return theCout;
}
ostream *& BaseRepository::cerrp() {
static ostream * theCerr = &std::cerr;
return theCerr;
}
ostream *& BaseRepository::clogp() {
static ostream * theClog = &std::clog;
return theClog;
}
bool & BaseRepository::updating() {
static bool theBool = false;
return theBool;
}
ObjectMap & BaseRepository::objects() {
static ObjectMap theObjectMap;
return theObjectMap;
}
ObjectSet & BaseRepository::allObjects() {
static ObjectSet theObjectSet;
return theObjectSet;
}
BaseRepository::TypeInterfaceMap & BaseRepository::interfaces() {
static TypeInterfaceMap theInterfaceMap;
return theInterfaceMap;
}
BaseRepository::TypeDocumentationMap & BaseRepository::documentations() {
static TypeDocumentationMap theDocumentationMap;
return theDocumentationMap;
}
BaseRepository::DirectorySet & BaseRepository::directories() {
- static char root[1][2] = {"/"};
- static DirectorySet theDirectories(root, root+1);
+ static DirectorySet theDirectories = {"/"};
return theDirectories;
}
vector<string> & BaseRepository::globalLibraries() {
static vector<string> theGlobalLibraries;
return theGlobalLibraries;
}
stack<string> & BaseRepository::currentReadDirStack() {
static stack<string> theCurrentReadDirStack;
if ( theCurrentReadDirStack.empty() ) theCurrentReadDirStack.push("");
return theCurrentReadDirStack;
}
vector<string> & BaseRepository::readDirs() {
// macro is passed in from -D compile flag
static vector<string>
theReadDirs(1, THEPEG_PKGDATADIR);
return theReadDirs;
}
const vector<string> & BaseRepository::listReadDirs() {
return BaseRepository::readDirs();
}
void BaseRepository::prependReadDir(string dir) {
readDirs().insert(readDirs().begin(), dir);
}
void BaseRepository::prependReadDir(const std::vector<std::string>& dirs) {
readDirs().insert(readDirs().begin(), dirs.begin(), dirs.end());
}
void BaseRepository::appendReadDir(string dir) {
readDirs().push_back(dir);
}
void BaseRepository::appendReadDir(const std::vector<std::string>& dirs) {
readDirs().insert(readDirs().end(), dirs.begin(), dirs.end());
}
BaseRepository::StringVector & BaseRepository::directoryStack() {
static StringVector theDirectoryStack(1, "/");
return theDirectoryStack;
}
void BaseRepository::Register(const InterfaceBase & ib, const type_info & i) {
const ClassDescriptionBase * db = DescriptionList::find(i);
if ( db ) interfaces()[db].insert(&ib);
}
void BaseRepository::
Register(const ClassDocumentationBase & cd, const type_info & i) {
const ClassDescriptionBase * db = DescriptionList::find(i);
if ( db ) documentations()[db] = &cd;
}
void BaseRepository::Register(IBPtr ip, string newName) {
DirectoryAppend(newName);
ip->name(newName);
Register(ip);
}
void BaseRepository::Register(IBPtr ip) {
if ( !ip || member(allObjects(), ip) ) return;
if ( member(objects(), ip->fullName()) )
throw RepoNameExistsException(ip->fullName());
objects()[ip->fullName()] = ip;
allObjects().insert(ip);
ip->clear();
ip->update();
ip->touch();
}
void BaseRepository::DirectoryAppend(string & name) {
if ( name == "." ) name = directoryStack().back();
if ( name[0] != '/' ) name = directoryStack().back() + name;
}
void BaseRepository::CreateDirectory(string name) {
DirectoryAppend(name);
if ( name[name.size()-1] != '/' ) name += "/";
if ( member(directories(), name) ) return;
directories().insert(name);
name = name.substr(0, name.size() - 1);
name = name.substr(0, name.rfind('/'));
if ( name.size() ) CreateDirectory(name);
}
void BaseRepository::CheckObjectDirectory(string name) {
if ( name[name.size() - 1] != '/' )
name = name.substr(0, name.rfind('/') + 1);
CheckDirectory(name);
}
void BaseRepository::CheckDirectory(string name) {
DirectoryAppend(name);
if ( name[name.size()-1] != '/' ) name += "/";
if ( member(directories(), name) ) return;
throw RepositoryNoDirectory(name);
}
void BaseRepository::ChangeDirectory(string name) {
DirectoryAppend(name);
if ( name[name.size()-1] != '/' ) name += "/";
if ( member(directories(), name) ) {
directoryStack().back() = name;
return;
}
throw RepositoryNoDirectory(name);
}
void BaseRepository::PushDirectory(string name) {
DirectoryAppend(name);
if ( name[name.size()-1] != '/' ) name += "/";
if ( member(directories(), name) ) {
directoryStack().push_back(name);
return;
}
throw RepositoryNoDirectory(name);
}
void BaseRepository::PopDirectory() {
if ( directoryStack().size() > 1 ) directoryStack().pop_back();
}
IBPtr BaseRepository::GetPointer(string name) {
ObjectMap::iterator it = objects().find(name);
return it == objects().end()? IBPtr(): it->second;
}
IVector BaseRepository::SearchDirectory(string name, string className) {
IVector ret;
DirectoryAppend(name);
const ClassDescriptionBase * cdb = 0;
if ( className.size() ) {
cdb = DescriptionList::find(className);
if ( !cdb ) return ret;
}
if ( name[name.size()-1] != '/' ) name += "/";
string::size_type size = name.size();
for ( ObjectMap::const_iterator i = objects().begin();
i != objects().end(); ++i ) {
if ( cdb && !DescriptionList::find(typeid(*(i->second)))->isA(*cdb) )
continue;
if ( i->first.substr(0, size) == name ) ret.push_back(i->second);
}
return ret;
}
IVector BaseRepository::GetObjectsReferringTo(IBPtr obj) {
IVector ret;
for ( ObjectMap::const_iterator i = objects().begin();
i != objects().end(); ++i ) {
if ( obj == i->second ) continue;
IVector ov = DirectReferences(i->second);
if ( member(ov, obj) ) ret.push_back(i->second);
}
return ret;
}
IVector BaseRepository::DirectReferences(IBPtr obj) {
IVector ov = obj->getReferences();
InterfaceMap interfaceMap = getInterfaces(typeid(*obj));
for ( InterfaceMap::iterator iit = interfaceMap.begin();
iit != interfaceMap.end(); ++iit ) {
IVector ovi = iit->second->getReferences(*obj);
ov.insert(ov.end(), ovi.begin(), ovi.end());
}
return ov;
}
void BaseRepository::
addReferences(tIBPtr obj, ObjectSet & refs) {
if ( !obj ) return;
refs.insert(obj);
IVector ov = obj->getReferences();
for ( IVector::const_iterator it = ov.begin(); it != ov.end(); ++it )
if ( !member(refs, *it) ) addReferences(*it, refs);
InterfaceMap interfaceMap = getInterfaces(typeid(*obj));
for ( InterfaceMap::iterator iit = interfaceMap.begin();
iit != interfaceMap.end(); ++iit ) {
IVector ov = iit->second->getReferences(*obj);
for ( IVector::const_iterator it = ov.begin(); it != ov.end(); ++it )
if ( !member(refs, *it) ) addReferences(*it, refs);
}
}
void BaseRepository::
addInterfaces(const ClassDescriptionBase & db,
InterfaceMap & interfaceMap, bool all) {
for ( ClassDescriptionBase::DescriptionVector::const_iterator it =
db.descriptions().begin(); it != db.descriptions().end(); ++it )
if ( *it ) addInterfaces(**it, interfaceMap, all);
TypeInterfaceMap::const_iterator cit = interfaces().find(&db);
if ( cit == interfaces().end() ) return;
for ( InterfaceSet::const_iterator iit = (cit->second).begin();
iit != (cit->second).end(); ++iit ) {
string n = (**iit).name();
while ( all && member(interfaceMap, n) ) n = "+" + n;
interfaceMap[n] = *iit;
}
}
InterfaceMap BaseRepository::getInterfaces(const type_info & ti, bool all) {
InterfaceMap interfaceMap;
const ClassDescriptionBase * db = DescriptionList::find(ti);
if ( !db ) return interfaceMap;
addInterfaces(*db, interfaceMap, all);
return interfaceMap;
}
void BaseRepository::
rebind(InterfacedBase & i, const TranslationMap & trans,
const IVector & defaults) {
InterfaceMap interfaceMap = getInterfaces(typeid(i), true);
for ( InterfaceMap::iterator iit = interfaceMap.begin();
iit != interfaceMap.end(); ++iit )
iit->second->rebind(i, trans, defaults);
i.rebind(trans);
}
void BaseRepository::update() {
for_each(allObjects(), mem_fun(&InterfacedBase::update));
clearAll(allObjects());
}
template <typename Set1, typename Set2>
bool overlap(const Set1 & s1, const Set2 & s2) {
typename Set1::const_iterator i1 = s1.begin();
typename Set2::const_iterator i2 = s2.begin();
while ( i1 != s1.end() && i2 != s2.end() ) {
if ( *i1 == *i2 ) return true;
if ( *i1 < *i2 ) {
i1 = s1.lower_bound(*i2);
if ( *i1 == *i2 ) return true;
++i1;
} else {
i2 = s2.lower_bound(*i1);
if ( *i1 == *i2 ) return true;
++i2;
}
}
return false;
}
void BaseRepository::remove(tIBPtr ip) {
ObjectMap::iterator it = objects().find(ip->fullName());
if ( it == objects().end() || ip != it->second ) return;
objects().erase(it);
allObjects().erase(ip);
}
string BaseRepository::remove(const ObjectSet & rmset) {
ObjectSet refset;
for ( ObjectSet::const_iterator oi = rmset.begin();
oi != rmset.end(); ++oi ) {
IVector ov = GetObjectsReferringTo(*oi);
refset.insert(ov.begin(), ov.end());
}
for ( ObjectSet::iterator oi = rmset.begin(); oi != rmset.end(); ++oi )
refset.erase(*oi);
if ( refset.empty() ) {
for ( ObjectSet::iterator oi = rmset.begin(); oi != rmset.end(); ++oi )
remove(*oi);
return "";
}
string ret = "Error: cannot remove the objects because the following "
"objects refers to some of them:\n";
for ( ObjectSet::iterator oi = refset.begin(); oi != refset.end(); ++oi )
ret += (**oi).fullName() + "\n";
return ret;
}
void BaseRepository::rename(tIBPtr ip, string newName) {
ObjectSet::iterator it = allObjects().find(ip);
if ( it == allObjects().end() ) {
Register(ip, newName);
return;
}
ObjectMap::iterator mit = objects().find(ip->fullName());
if ( mit == objects().end() || mit->second != ip )
throw RepoNameException(ip->fullName());
objects().erase(mit);
ip->name(newName);
while ( member(objects(), ip->fullName()) ) ip->name(ip->fullName() + "#");
objects()[ip->fullName()] = ip;
}
const InterfaceBase * BaseRepository::FindInterface(IBPtr ip, string name) {
InterfaceMap imap = getInterfaces(typeid(*ip), false);
InterfaceMap::iterator it = imap.find(name);
return it == imap.end()? 0: it->second;
}
const ClassDocumentationBase * BaseRepository::getDocumentation(tcIBPtr ip) {
TypeDocumentationMap::const_iterator cdoc =
documentations().find(DescriptionList::find(typeid(*ip)));
return cdoc != documentations().end()? cdoc->second: 0;
}
string BaseRepository::getModelDescription(tcIBPtr ip) {
const ClassDocumentationBase * cd = getDocumentation(ip);
return cd? cd->modelDescription(): string("");
}
string BaseRepository::getModelReferences(tcIBPtr ip) {
const ClassDocumentationBase * cd = getDocumentation(ip);
return cd? cd->modelReferences(): string("");
}
IBPtr BaseRepository::TraceObject(string path) {
DirectoryAppend(path);
string::size_type colon = path.find(':');
IBPtr ip = GetPointer(path.substr(0, colon));
if ( !ip ) {
// Do special check if this is a decay mode.
string name = path.substr(0, colon);
string::size_type slash = name.rfind('/');
if ( slash != string::npos ) name = name.substr(slash + 1);
if ( name.find("->") != string::npos && name[name.length() - 1] == ';' ) {
vector<DMPtr> save;
DMPtr dm = DecayMode::constructDecayMode(name, &save);
if ( dm ) ip = dynamic_ptr_cast<DMPtr>(GetPointer(path.substr(0, slash + 1)
+ dm->tag()));
if ( ip ) Throw<Exception>()
<< "Warning: rewriting DecayMode name '"
<< path.substr(0, colon).substr(slash + 1) << "' to '"
<< ip->name() << Exception::warning;
}
}
while ( colon != string::npos ) {
if ( !ip ) throw RepositoryNotFound(path);
path = path.substr(colon+1);
colon = path.find(':');
string::size_type bra = path.find('[');
const InterfaceBase * ifb =
FindInterface(ip, path.substr(0, min(colon, bra)));
const ReferenceBase * rb = dynamic_cast<const ReferenceBase *>(ifb);
if ( rb ) {
ip = rb->get(*ip);
continue;
}
const RefVectorBase * rvb = dynamic_cast<const RefVectorBase *>(ifb);
if ( rvb ) {
unsigned int place = 0;
if ( bra < colon ) {
string::size_type ket = path.find(']');
place = atoi(path.substr(bra + 1,ket - bra - 1).c_str());
}
IVector iv = rvb->get(*ip);
if ( place >= iv.size() ) throw RepositoryNotFound(path);
ip = iv[place];
continue;
}
const CommandBase * cb = dynamic_cast<const CommandBase *>(ifb);
if ( cb ) {
string::size_type ket = path.find(']');
string newobj = cb->cmd(*ip, path.substr(bra + 1,ket - bra - 1));
ip = GetPointer(newobj);
continue;
}
throw RepositoryNotFound(path);
}
if ( !ip ) throw RepositoryNotFound(path);
return ip;
}
IBPtr BaseRepository::getObjectFromNoun(string noun) {
string::size_type colon = noun.rfind(':');
return TraceObject(noun.substr(0, colon));
}
string BaseRepository::getInterfaceFromNoun(string noun) {
string::size_type colon = noun.rfind(':');
string interface = noun.substr(colon+1);
string::size_type bra = interface.find('[');
if ( bra != string::npos ) return interface.substr(0, bra);
else return interface;
}
string BaseRepository::getPosArgFromNoun(string noun) {
string::size_type colon = noun.rfind(':');
string interface = noun.substr(colon+1);
string::size_type bra = interface.find('[');
if ( bra != string::npos ) {
string::size_type ket = interface.find(']');
return interface.substr(bra + 1,ket - bra - 1);
}
return "";
}
string BaseRepository::
GetInterfacedBaseClasses(const ClassDescriptionBase * cdb) {
if ( !cdb || cdb->name() == "ThePEG::Interfaced" ||
cdb->name() == "ThePEG::InterfacedBase" ) return "";
string ret = cdb->name() + "\n";
for ( int i = 0, N = cdb->descriptions().size(); i < N; ++i )
ret += GetInterfacedBaseClasses(cdb->descriptions()[i]);
return ret;
}
struct InterfaceOrder {
bool operator()(const InterfaceBase * x, const InterfaceBase * y) const {
return x->rank() > y->rank() ||
( x->rank() == y->rank() && x->name() < y->name() );
}
};
void BaseRepository::readSetup(tIBPtr ip, istream & is) {
ip->setup(is);
}
string BaseRepository::exec(string command, ostream &) {
string verb = StringUtils::car(command);
command = StringUtils::cdr(command);
if ( verb.empty() || verb[0] == '#' ) return "";
try {
if ( verb == "DISABLEREADONLY" ) {
InterfaceBase::NoReadOnly = true;
return "";
}
if ( verb == "ENABLEREADONLY" ) {
InterfaceBase::NoReadOnly = false;
return "";
}
if ( verb == "cd" || verb == "pushd" || verb == "mkdir") {
string dir = StringUtils::car(command);
if ( verb == "cd" )
ChangeDirectory(dir);
else if ( verb == "pushd" )
PushDirectory(dir);
else
CreateDirectory(dir);
return "";
}
if ( verb == "popd" ) {
PopDirectory();
return "";
}
if ( verb == "pwd" ) return directoryStack().back();
if ( verb == "dirs" ) {
string ret;
for ( StringVector::reverse_iterator it = directoryStack().rbegin();
it != directoryStack().rend(); ++it ) ret += *it;
return ret;
}
if ( verb == "cp" || verb == "mv" ) {
string oldname = StringUtils::car(command);
DirectoryAppend(oldname);
IBPtr obj = GetPointer(oldname);
if ( !obj ) return "Error: No object named '" + oldname + "' available.";
command = StringUtils::cdr(command);
string newname = StringUtils::car(command);
DirectoryAppend(newname);
if ( newname[newname.size() - 1] == '/' ) newname += obj->name();
if ( verb == "cp" ) obj = obj->fullclone();
rename(obj, newname);
return "";
}
if ( verb == "check" ) {
string name = StringUtils::car(command);
if ( directories().find(name) != directories().end() ) return name;
if ( objects().find(name) != objects().end() ) return name;
return "Not found";
}
if ( verb == "ls" ) {
string className;
string dir = StringUtils::car(command);
if ( dir.size() ) {
PushDirectory(dir);
command = StringUtils::cdr(command);
className = StringUtils::car(command);
}
string ret;
string thisdir = directoryStack().back();
for ( DirectorySet::iterator it = directories().begin();
it != directories().end(); ++it ) {
string d = *it;
if ( d.size() <= thisdir.size() ) continue;
string d0 = d.substr(0, thisdir.size());
string d1 = d.substr(thisdir.size());
if ( d0 == thisdir && d1.find('/') == d1.size() - 1 ) {
if ( className.size() && SearchDirectory(d, className).empty() )
continue;
ret += (dir.size()? d: d1) + "\n";
}
}
for ( ObjectMap::iterator it = objects().begin();
it != objects().end(); ++it ) {
if ( className.size() ) {
const ClassDescriptionBase * cdb = DescriptionList::find(className);
if ( cdb &&
!DescriptionList::find(typeid(*(it->second)))->isA(*cdb) )
continue;
}
if ( thisdir + it->second->name() == it->first )
ret += (dir.size()? it->first: it->second->name()) + '\n';
}
if ( dir.size() ) PopDirectory();
return ret;
}
if ( verb == "library" ) {
string library = StringUtils::car(command);
if ( library.empty() ) return "Error: No library specified.";
if ( !DynamicLoader::load(library) )
return "Error: Could not load library " + library +
"\n - " + DynamicLoader::lastErrorMessage;
return "";
}
if ( verb == "globallibrary" ) {
string library = StringUtils::car(command);
if ( library.empty() ) return "Error: No library specified.";
if ( !DynamicLoader::load(library) )
return "Error: Could not load library " + library +
"\n - " + DynamicLoader::lastErrorMessage;
globalLibraries().push_back(library);
return "";
}
if ( verb == "rmgloballibrary" ) {
string library = StringUtils::car(command);
if ( library.empty() ) return "Error: No library specified.";
vector<string>::iterator it;
while ( (it = find(globalLibraries(), library)) != globalLibraries().end() )
globalLibraries().erase(it);
return "";
}
if ( verb == "appendpath" ) {
string path = StringUtils::car(command);
if ( !path.empty() ) DynamicLoader::appendPath(path);
return "";
}
if ( verb == "lspaths" ) {
string paths;
for ( int i = 0, N = DynamicLoader::allPaths().size(); i < N; ++i )
paths += DynamicLoader::allPaths()[i] + "\n";
return paths;
}
if ( verb == "prependpath" ) {
string path = StringUtils::car(command);
if ( !path.empty() ) DynamicLoader::prependPath(path);
return "";
}
if ( verb == "create" ) {
string className = StringUtils::car(command);
command = StringUtils::cdr(command);
string name = StringUtils::car(command);
const ClassDescriptionBase * db = DescriptionList::find(className);
command = StringUtils::cdr(command);
while ( !db && command.length() ) {
string library = StringUtils::car(command);
command = StringUtils::cdr(command);
DynamicLoader::load(library);
db = DescriptionList::find(className);
}
if ( !db ) {
string msg = "Error: " + className + ": No such class found.";
if ( !DynamicLoader::lastErrorMessage.empty() )
msg += "\nerror message from dynamic loader:\n" +
DynamicLoader::lastErrorMessage;
return msg;
}
IBPtr obj = dynamic_ptr_cast<IBPtr>(db->create());
if ( !obj ) return "Error: Could not create object of class "+className;
if ( name.empty() ) return "Error: No name specified.";
Register(obj, name);
return "";
}
if ( verb == "setup" ) {
string name = StringUtils::car(command);
DirectoryAppend(name);
IBPtr obj = GetPointer(name);
if ( !obj ) return "Error: Could not find object named " + name;
istringstream is(StringUtils::cdr(command));
readSetup(obj, is);
return "";
}
if ( verb == "rm" ) {
ObjectSet rmset;
while ( !command.empty() ) {
string name = StringUtils::car(command);
DirectoryAppend(name);
IBPtr obj = GetPointer(name);
if ( !obj ) return "Error: Could not find object named " + name;
rmset.insert(obj);
command = StringUtils::cdr(command);
}
return remove(rmset);
}
if ( verb == "rmdir" || verb == "rrmdir" ) {
string dir = StringUtils::car(command);
DirectoryAppend(dir);
if ( dir[dir.size() - 1] != '/' ) dir += '/';
if ( !member(directories(), dir) )
return verb == "rmdir"? "Error: No such directory.": "";
IVector ov = SearchDirectory(dir);
if ( ov.size() && verb == "rmdir" )
return "Error: Cannot remove a non-empty directory. "
"(Use rrmdir do remove all object and subdirectories.)";
ObjectSet rmset(ov.begin(), ov.end());
string ret = remove(rmset);
if ( !ret.empty() ) return ret;
StringVector dirs(directories().begin(), directories().end());
for ( int i = 0, N = dirs.size(); i < N; ++ i )
if ( dirs[i].substr(0, dir.size()) == dir )
directories().erase(dirs[i]);
for ( int i = 0, N = directoryStack().size(); i < N; ++i )
if ( directoryStack()[i].substr(0, dir.size()) == dir )
directoryStack()[i] = '/';
return "";
}
if ( verb == "rcp" ) {
string name = StringUtils::car(command);
DirectoryAppend(name);
string newName = StringUtils::car(StringUtils::cdr(command));
if ( newName.empty() )
return "Error: No destination directory specified.";
DirectoryAppend(newName);
CreateDirectory(newName);
if ( newName[newName.size() - 1] != '/' ) newName += '/';
IBPtr obj = GetPointer(name);
if ( name[name.size() - 1] != '/' ) name += '/';
IVector ov = SearchDirectory(name);
ov.push_back(obj);
if ( ov.empty() ) return "Error: No such object or directory.";
ObjectSet toclone;
for ( IVector::iterator i = ov.begin(); i != ov.end(); ++i ) {
toclone.insert(*i);
addReferences(*i, toclone);
}
for ( ObjectSet::iterator i = toclone.begin(); i != toclone.end(); ++i )
Register((**i).clone(), newName + (**i).name());
return "";
}
if ( verb == "rebind" ) {
// For all objects in the repository, replace any references to
// the first object given with references to the second
// object. The two objects will change names
IBPtr ip1 = TraceObject(StringUtils::car(command));
string newname = StringUtils::car(StringUtils::cdr(command));
DirectoryAppend(newname);
IBPtr ip2 = GetPointer(newname);
if ( !ip2 ) {
ip2 = ip1->fullclone();
rename(ip2, newname);
}
TranslationMap trans;
trans[ip1] = ip2;
IVector objs = GetObjectsReferringTo(ip1);
for ( int i = 0, N = objs.size(); i < N; ++i ) rebind(*objs[i], trans, IVector());
}
if ( verb == "doxygendump" ) {
string spacename = StringUtils::car(command);
command = StringUtils::cdr(command);
string filename = StringUtils::car(command);
ofstream os(filename.c_str());
for ( TypeDocumentationMap::const_iterator it = documentations().begin();
it != documentations().end(); ++it ) {
const ClassDescriptionBase & db = *(it->first);
string classname = db.name();
if ( classname.substr(0, spacename.length()) != spacename ) continue;
string briefname = classname.substr(spacename.length());
os << "/** \\page " << briefname << "Interfaces "
<< "Interfaces defined for the " << classname << " class.\n\n"
<< "\\par Brief class description:\n";
string doc = it->second->documentation();
if ( doc.substr(0,25) == "There is no documentation" )
os << "See " << classname << "\n\n";
else
os << doc << "<br>See also " << classname << "\n\n";
TypeInterfaceMap::const_iterator isit = interfaces().find(it->first);
if ( isit == interfaces().end() || isit->second.empty() ) {
os << "There are no interfaces declared for this class.\n\n";
} else {
const InterfaceSet & ints = isit->second;
for ( InterfaceSet::const_iterator iit = ints.begin();
iit != ints.end(); ++iit )
(**iit).doxygenDescription(os);
}
string baserefs = "";
int nbases = 0;
for ( int ib = 0, N = db.descriptions().size(); ib < N; ++ib ) {
if ( documentations().find(db.descriptions()[ib]) ==
documentations().end() ) continue;
const ClassDescriptionBase & bdb = *db.descriptions()[ib];
if ( nbases ) baserefs += " and ";
string briefname = bdb.name().substr(bdb.name().rfind("::") + 2);
baserefs += "\\ref " + briefname +
"Interfaces \"" + bdb.name() + "\"";
++nbases;
}
if ( nbases == 1 )
os << "<hr>There may be interfaces inherited from the "
<< baserefs << " class.";
else if ( nbases > 1 )
os << "<hr>There may be interfaces inherited from the "
<< "following classes: " << baserefs << ".";
os << "\n\n*/\n\n";
}
return "";
}
if ( verb == "mset" || verb == "msetdef" || verb == "minsert" ||
verb == "mdo" || verb == "mget" || verb == "mdef" || verb == "mmin" ||
verb == "mmax" || verb == "merase" || verb == "msend" ) {
if ( verb == "msend" ) verb = "mdo";
string dir = StringUtils::car(command);
command = StringUtils::cdr(command);
string className = StringUtils::car(command);
command = StringUtils::cdr(command);
string interface = StringUtils::car(command);
string arguments = StringUtils::cdr(command);
string::size_type bra = interface.find('[');
if ( bra != string::npos ) {
string::size_type ket = interface.find(']');
arguments = interface.substr(bra + 1,ket - bra - 1) + " " + arguments;
interface = interface.substr(0, bra);
}
IVector ov = SearchDirectory(dir, className);
if ( ov.empty() ) return "Error: no matching objects found.";
string ret;
verb = verb.substr(1);
for ( IVector::size_type i = 0; i < ov.size(); ++i ) {
const InterfaceBase * ifb = FindInterface(ov[i], interface);
if ( !ifb ) continue;
string mess = ifb->exec(*ov[i], verb, arguments);
if ( !mess.empty() ) ret += ov[i]->fullName() + ": " + mess + "\n";
}
return ret.substr(0, ret.size() - 1);
}
if ( verb == "set" || verb == "setdef" || verb == "insert" ||
verb == "do" || verb == "get" || verb == "def" || verb == "min" ||
verb == "max" || verb == "describe" || verb == "fulldescribe" ||
verb == "erase" || verb == "clear" || verb == "send" || verb == "newdef" ) {
if ( verb == "send" ) verb = "do";
if ( verb == "newdef" && !InterfaceBase::NoReadOnly )
return "Error: The default value of an interface is a read-only "
"entity. Use the command 'DISABLEREADONLY' to override.";
string noun = StringUtils::car(command);
string arguments = getPosArgFromNoun(noun) + " "
+ StringUtils::cdr(command);
IBPtr ip = getObjectFromNoun(noun);
const InterfaceBase * ifb = FindInterface(ip, getInterfaceFromNoun(noun));
if ( !ifb && verb != "describe" && verb != "fulldescribe" ) {
string ret = "Error: The interface '" + noun + "' was not found.\n";
ret += "Valid interfaces:\n";
InterfaceMap imap = getInterfaces(typeid(*ip));
for ( InterfaceMap::iterator it = imap.begin(); it != imap.end(); ++it )
ret += "* " + it->second->name() + "\n";
return ret;
}
if ( verb == "describe" ) {
if ( ifb ) return ifb->description();
const ClassDescriptionBase * cd = DescriptionList::find(typeid(*ip));
string ret = "Object '" + ip->name() + "' of class '" +
cd->name() + "':\n";
TypeDocumentationMap::const_iterator cdoc = documentations().find(cd);
if ( cdoc != documentations().end() )
ret += cdoc->second->documentation() + "\n";
ret +="Interfaces:\n";
InterfaceMap imap = getInterfaces(typeid(*ip));
for ( InterfaceMap::iterator it = imap.begin(); it != imap.end(); ++it )
ret += "* " + it->second->name() + "\n";
return ret;
} else if ( verb == "fulldescribe" ) {
if ( ifb ) return ifb->fullDescription(*ip);
ostringstream ret;
const ClassDescriptionBase * cd = DescriptionList::find(typeid(*ip));
TypeDocumentationMap::const_iterator cdoc = documentations().find(cd);
ret << ip->fullName() << endl << cd->name() << endl;
if ( cdoc != documentations().end() )
ret << cdoc->second->documentation() << endl;
ret << "Interfaces:" << endl;
InterfaceMap imap = getInterfaces(typeid(*ip));
typedef set<const InterfaceBase *, InterfaceOrder> InterfaceSet;
InterfaceSet iset;
for ( InterfaceMap::iterator it = imap.begin(); it != imap.end(); ++it )
iset.insert(it->second);
double rank = 1.0;
for ( InterfaceSet::iterator it = iset.begin();
it != iset.end(); ++it ) {
if ( rank >= 0.0 && (**it).rank() < 0.0 ) ret << "0" << endl;
rank = (**it).rank();
ret << (**it).type() << " " << (**it).name() << endl;
}
return ret.str();
}
else
return ifb->exec(*ip, verb, arguments);
}
if ( verb == "baseclasses" ) {
string className = StringUtils::car(command);
const ClassDescriptionBase * cdb = 0;
if ( className.size() ) {
cdb = DescriptionList::find(className);
if ( !cdb ) return "Error: no class '" + className + "' found.";
}
return GetInterfacedBaseClasses(cdb);
}
if ( verb == "describeclass" ) {
string className = StringUtils::car(command);
const ClassDescriptionBase * cdb = 0;
if ( className.size() ) {
cdb = DescriptionList::find(className);
if ( !cdb ) return "Error: no class '" + className + "' found.";
}
TypeDocumentationMap::const_iterator cdoc = documentations().find(cdb);
if ( cdoc != documentations().end() )
return cdoc->second->documentation() + "\n";
else
return "";
}
if ( verb == "lsclass" ) {
string className = StringUtils::car(command);
const ClassDescriptionBase * cdb = 0;
if ( className.size() ) {
cdb = DescriptionList::find(className);
if ( !cdb ) return "Error: no class '" + className + "' found.";
}
vector<const ClassDescriptionBase *> classes;
if ( cdb && !cdb->abstract() ) classes.push_back(cdb);
for ( DescriptionList::DescriptionMap::const_iterator
it = DescriptionList::all().begin();
it != DescriptionList::all().end(); ++it ) {
if ( it->second == cdb || it->second->abstract() ) continue;
if ( cdb && !it->second->isA(*cdb) ) continue;
classes.push_back(it->second);
}
if ( classes.empty() ) return "Error: no classes found.";
string ret;
for ( int i = 0, N = classes.size(); i < N; ++i )
ret += classes[i]->name() + "\n";
return ret;
}
}
catch (const Exception & e) {
e.handle();
return "Error: " + e.message();
}
return "Error: Unrecognized command '" + verb + "'.";
}
BadClassClone::BadClassClone(const InterfacedBase & o) {
theMessage << "Could not clone the object '" << o.name()
<< "' of class '" << TypeInfo::name(o)
<< "' because the class does not"
<< " implement a working 'clone' method.";
severity(abortnow);
}
BadClone::BadClone(const InterfacedBase & o) {
theMessage << "Could not clone the object '" << o.name()
<< "' of class '" << TypeInfo::name(o)
<< "' because the clone method threw an unknown exception.";
severity(abortnow);
}
RepoNameException::RepoNameException(string name) {
theMessage << "The object '" << name << "' is present in the Repository but "
<< "under a different name. This means that the name of the "
<< "object has been illegally changed outside of the Repository.";
severity(abortnow);
}
RepoNameExistsException::RepoNameExistsException(string name) {
theMessage
<< "The object '" << name
<< "' was not created as another object with that name already exists.";
severity(warning);
}
RepositoryNoDirectory::RepositoryNoDirectory(string name) {
theMessage << "The directory '" << name << "' does not exist.";
severity(warning);
}
RepositoryNotFound::RepositoryNotFound(string name) {
theMessage << "There was no object named '" << name << "' in the repository.";
severity(warning);
}
RepositoryClassMisMatch::
RepositoryClassMisMatch(const InterfacedBase & o, string name) {
theMessage << "The requested object '" << o.fullName() << "' was not of the "
<< "specified type (" << name << ").";
severity(warning);
}
diff --git a/Repository/EventGenerator.cc b/Repository/EventGenerator.cc
--- a/Repository/EventGenerator.cc
+++ b/Repository/EventGenerator.cc
@@ -1,1414 +1,1415 @@
// -*- C++ -*-
//
// EventGenerator.cc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2017 Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the EventGenerator class.
//
#include "EventGenerator.h"
#include "EventGenerator.xh"
#include "ThePEG/Handlers/EventHandler.h"
#include "Repository.h"
#include "ThePEG/Utilities/HoldFlag.h"
#include "ThePEG/Utilities/Debug.h"
#include "ThePEG/Utilities/DebugItem.h"
#include "ThePEG/Interface/Interfaced.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Command.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/PDT/ParticleData.h"
#include "ThePEG/PDT/MatcherBase.h"
#include "ThePEG/PDT/DecayMode.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Repository/Strategy.h"
#include "ThePEG/Repository/CurrentGenerator.h"
#include "ThePEG/Handlers/AnalysisHandler.h"
#include "ThePEG/Analysis/FactoryBase.h"
#include "ThePEG/Handlers/EventManipulator.h"
#include "ThePEG/Handlers/LuminosityFunction.h"
#include "ThePEG/MatrixElement/MEBase.h"
#include "ThePEG/EventRecord/Event.h"
#include "ThePEG/Handlers/SubProcessHandler.h"
#include "ThePEG/Handlers/CascadeHandler.h"
#include "ThePEG/Handlers/HadronizationHandler.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Config/algorithm.h"
#include "ThePEG/Utilities/DynamicLoader.h"
#include <cstdlib>
#include "ThePEG/Repository/Main.h"
#include <csignal>
#ifdef ThePEG_TEMPLATES_IN_CC_FILE
#include "EventGenerator.tcc"
#endif
using namespace ThePEG;
namespace {
volatile sig_atomic_t THEPEG_SIGNAL_STATE = 0;
}
// signal handler function
// very restricted in what it is allowed do
// without causing undefined behaviour
extern "C" {
void thepegSignalHandler(int id) {
THEPEG_SIGNAL_STATE=id;
signal(id,SIG_DFL);
}
}
void EventGenerator::checkSignalState() {
if (THEPEG_SIGNAL_STATE) {
log() << "Caught signal " << THEPEG_SIGNAL_STATE << ". Exiting ..." << std::endl;
finalize();
exit(0);
}
}
EventGenerator::EventGenerator()
: thePath("."), theNumberOfEvents(1000), theQuickSize(7000),
preinitializing(false), ieve(0), weightSum(0.0),
theDebugLevel(0), logNonDefault(-1), printEvent(0), dumpPeriod(0),
keepAllDumps(false),
debugEvent(0), maxWarnings(10), maxErrors(10), theCurrentRandom(0),
theCurrentGenerator(0), useStdout(false), theIntermediateOutput(false) {}
EventGenerator::EventGenerator(const EventGenerator & eg)
: Interfaced(eg), theDefaultObjects(eg.theDefaultObjects),
theLocalParticles(eg.theLocalParticles),
theStandardModel(eg.theStandardModel),
theStrategy(eg.theStrategy), theRandom(eg.theRandom),
theEventHandler(eg.theEventHandler),
theAnalysisHandlers(eg.theAnalysisHandlers),
theHistogramFactory(eg.theHistogramFactory),
theEventManipulator(eg.theEventManipulator),
thePath(eg.thePath), theRunName(eg.theRunName),
theNumberOfEvents(eg.theNumberOfEvents), theObjects(eg.theObjects),
theObjectMap(eg.theObjectMap),
theParticles(eg.theParticles), theQuickParticles(eg.theQuickParticles),
theQuickSize(eg.theQuickSize), preinitializing(false),
theMatchers(eg.theMatchers),
usedObjects(eg.usedObjects), ieve(eg.ieve), weightSum(eg.weightSum),
theDebugLevel(eg.theDebugLevel), logNonDefault(eg.logNonDefault),
printEvent(eg.printEvent), dumpPeriod(eg.dumpPeriod),
keepAllDumps(eg.keepAllDumps),
debugEvent(eg.debugEvent),
maxWarnings(eg.maxWarnings), maxErrors(eg.maxErrors), theCurrentRandom(0),
theCurrentGenerator(0),
theCurrentEventHandler(eg.theCurrentEventHandler),
theCurrentStepHandler(eg.theCurrentStepHandler),
useStdout(eg.useStdout),
theIntermediateOutput(eg.theIntermediateOutput) {}
EventGenerator::~EventGenerator() {
if ( theCurrentRandom ) delete theCurrentRandom;
if ( theCurrentGenerator ) delete theCurrentGenerator;
}
IBPtr EventGenerator::clone() const {
return new_ptr(*this);
}
IBPtr EventGenerator::fullclone() const {
return new_ptr(*this);
}
tcEventPtr EventGenerator::currentEvent() const {
return eventHandler()->currentEvent();
}
CrossSection EventGenerator::histogramScale() const {
return eventHandler()->histogramScale();
}
CrossSection EventGenerator::integratedXSec() const {
return eventHandler()->integratedXSec();
}
CrossSection EventGenerator::integratedXSecErr() const {
return eventHandler()->integratedXSecErr();
}
void EventGenerator::setSeed(long seed) {
random().setSeed(seed);
ostringstream s;
s << seed;
const InterfaceBase * ifb = BaseRepository::FindInterface(theRandom, "Seed");
ifb->exec(*theRandom, "set", s.str());
}
void
EventGenerator::setup(string newRunName,
ObjectSet & newObjects,
ParticleMap & newParticles,
MatcherSet & newMatchers) {
HoldFlag<int> debug(Debug::level, Debug::isset? Debug::level: theDebugLevel);
theRunName = newRunName;
theObjects.swap(newObjects);
theParticles.swap(newParticles);
theMatchers.swap(newMatchers);
theObjectMap.clear();
for ( ObjectSet::const_iterator it = objects().begin();
it != objects().end(); ++it ) theObjectMap[(**it).fullName()] = *it;
UseRandom currentRandom(theRandom);
CurrentGenerator currentGenerator(this);
// Force update of all objects and then reset.
touch();
for_each(theObjects, mem_fun(&InterfacedBase::touch));
update();
for_each(theObjects, mem_fun(&InterfacedBase::update));
clear();
BaseRepository::clearAll(theObjects);
init();
}
IBPtr EventGenerator::getPointer(string name) const {
ObjectMap::const_iterator it = objectMap().find(name);
if ( it == objectMap().end() ) return IBPtr();
else return it->second;
}
void EventGenerator::openOutputFiles() {
if ( !useStdout ) {
logfile().open((filename() + ".log").c_str());
theOutFileName = filename() + ".out";
outfile().open(theOutFileName.c_str());
outfile().close();
theOutStream.str("");
}
out() << Repository::banner() << endl;
log() << Repository::banner() << endl;
}
void EventGenerator::closeOutputFiles() {
flushOutputFile();
if ( !useStdout ) logfile().close();
}
void EventGenerator::flushOutputFile() {
if ( !useStdout ) {
outfile().open(theOutFileName.c_str(), ios::out|ios::app);
outfile() << theOutStream.str();
outfile().close();
} else
BaseRepository::cout() << theOutStream.str();
theOutStream.str("");
}
void EventGenerator::doinit() {
HoldFlag<int> debug(Debug::level, Debug::isset? Debug::level: theDebugLevel);
// First initialize base class and random number generator.
Interfaced::doinit();
random().init();
// Make random generator and this available in standard static
// classes.
UseRandom useRandom(theRandom);
CurrentGenerator currentGenerator(this);
// First initialize all objects which have requested this by
// implementing a InterfacedBase::preInitialize() function which
// returns true.
while ( true ) {
HoldFlag<bool> hold(preinitializing, true);
ObjectSet preinits;
for ( ObjectSet::iterator it = objects().begin();
it != objects().end(); ++it )
if ( (**it).preInitialize() &&
(**it).state() == InterfacedBase::uninitialized )
preinits.insert(*it);
if ( preinits.empty() ) break;
for_each(preinits, mem_fun(&InterfacedBase::init));
}
// Initialize the quick access to particles.
theQuickParticles.clear();
theQuickParticles.resize(2*theQuickSize);
for ( ParticleMap::const_iterator pit = theParticles.begin();
pit != theParticles.end(); ++pit )
if ( abs(pit->second->id()) < theQuickSize )
theQuickParticles[pit->second->id()+theQuickSize] = pit->second;
// Then call the init method for all objects. Start with the
// standard model and the strategy.
standardModel()->init();
if ( strategy() ) strategy()->init();
eventHandler()->init();
// initialize particles first
for(ParticleMap::const_iterator pit = particles().begin();
pit != particles().end(); ++pit) pit->second->init();
for_each(objects(), mem_fun(&InterfacedBase::init));
// Then initialize the Event Handler calculating initial cross
// sections and stuff.
eventHandler()->initialize();
}
void EventGenerator::doinitrun() {
HoldFlag<int> debug(Debug::level, Debug::isset? Debug::level: theDebugLevel);
signal(SIGHUP, thepegSignalHandler);
signal(SIGINT, thepegSignalHandler);
signal(SIGTERM,thepegSignalHandler);
currentEventHandler(eventHandler());
Interfaced::doinitrun();
random().initrun();
// Then call the init method for all objects. Start with the
// standard model and the strategy.
standardModel()->initrun();
if ( strategy() ) {
strategy()->initrun();
if ( ! strategy()->versionstring().empty() ) {
out() << ">> " << strategy()->versionstring() << '\n' << endl;
log() << ">> " << strategy()->versionstring() << '\n' << endl;
}
}
// initialize particles first
for(ParticleMap::const_iterator pit = particles().begin();
pit != particles().end(); ++pit) {
pit->second->initrun();
}
eventHandler()->initrun();
for_each(objects(), mem_fun(&InterfacedBase::initrun));
if ( logNonDefault > 0 || ( ThePEG_DEBUG_LEVEL && logNonDefault == 0 ) ) {
vector< pair<IBPtr, const InterfaceBase *> > changed =
Repository::getNonDefaultInterfaces(objects());
if ( changed.size() ) {
log() << string(78, '=') << endl
<< "The following interfaces have non-default values (default):"
<< endl << string(78, '-') << endl;
for ( int i = 0, N = changed.size(); i < N; ++i ) {
log() << changed[i].first->fullName() << ":"
<< changed[i].second->name() << " = "
<< changed[i].second->exec(*changed[i].first, "notdef", "")
<< endl;
}
log() << string(78,'=') << endl;
}
}
weightSum = 0.0;
}
PDPtr EventGenerator::getParticleData(PID id) const {
long newId = id;
if ( abs(newId) < theQuickSize && theQuickParticles.size() )
return theQuickParticles[newId+theQuickSize];
ParticleMap::const_iterator it = theParticles.find(newId);
if ( it == theParticles.end() ) return PDPtr();
return it->second;
}
PPtr EventGenerator::getParticle(PID newId) const {
tcPDPtr pd = getParticleData(newId);
if ( !pd ) return PPtr();
return pd->produceParticle();
}
void EventGenerator::finalize() {
UseRandom currentRandom(theRandom);
CurrentGenerator currentGenerator(this);
finish();
finally();
}
void EventGenerator::dofinish() {
HoldFlag<int> debug(Debug::level, Debug::isset? Debug::level: theDebugLevel);
// first write out statistics from the event handler.
eventHandler()->statistics(out());
// Call the finish method for all other objects.
for_each(objects(), mem_fun(&InterfacedBase::finish));
if ( theExceptions.empty() ) {
log() << "No exceptions reported in this run.\n";
} else {
log() << "\nThe following exception classes were reported in this run:\n";
for ( ExceptionMap::iterator it = theExceptions.begin();
it != theExceptions.end(); ++it ) {
string severity;
switch ( it->first.second ) {
case Exception::info : severity="info"; break;
case Exception::warning : severity="warning"; break;
case Exception::setuperror : severity="setuperror"; break;
case Exception::eventerror : severity="eventerror"; break;
case Exception::runerror : severity="runerror"; break;
case Exception::maybeabort : severity="maybeabort"; break;
case Exception::abortnow : severity="abortnow"; break;
default : severity="unknown";
}
log() << it->first.first << ' ' << severity
<< " (" << it->second << " times)\n";
}
}
theExceptions.clear();
const string & msg = theMiscStream.str();
if ( ! msg.empty() ) {
log() << endl
<< "Miscellaneous output from modules to the standard output:\n\n"
<< msg;
theMiscStream.str("");
}
flushOutputFile();
}
void EventGenerator::finally() {
generateReferences();
closeOutputFiles();
if ( theCurrentRandom ) delete theCurrentRandom;
if ( theCurrentGenerator ) delete theCurrentGenerator;
theCurrentRandom = 0;
theCurrentGenerator = 0;
}
void EventGenerator::initialize(bool initOnly) {
UseRandom currentRandom(theRandom);
CurrentGenerator currentGenerator(this);
doInitialize(initOnly);
}
bool EventGenerator::loadMain(string file) {
initialize();
UseRandom currentRandom(theRandom);
CurrentGenerator currentGenerator(this);
Main::eventGenerator(this);
bool ok = DynamicLoader::load(file);
finish();
finally();
return ok;
}
void EventGenerator::go(long next, long maxevent, bool tics) {
UseRandom currentRandom(theRandom);
CurrentGenerator currentGenerator(this);
doGo(next, maxevent, tics);
}
EventPtr EventGenerator::shoot() {
static DebugItem debugfpu("ThePEG::FPU", 1);
if ( debugfpu ) Debug::unmaskFpuErrors();
UseRandom currentRandom(theRandom);
CurrentGenerator currentGenerator(this);
checkSignalState();
EventPtr event = doShoot();
if ( event ) weightSum += event->weight();
DebugItem::tic();
return event;
}
EventPtr EventGenerator::doShoot() {
EventPtr event;
if ( N() >= 0 && ++ieve > N() ) return event;
HoldFlag<int> debug(Debug::level, Debug::isset? Debug::level: theDebugLevel);
do {
int state = 0;
int loop = 1;
eventHandler()->clearEvent();
try {
do {
// Generate a full event or part of an event
if ( eventHandler()->empty() ) event = eventHandler()->generateEvent();
else event = eventHandler()->continueEvent();
if ( eventHandler()->empty() ) loop = -loop;
// Analyze the possibly uncomplete event
for ( AnalysisVector::iterator it = analysisHandlers().begin();
it != analysisHandlers().end(); ++it )
(**it).analyze(event, ieve, loop, state);
// Manipulate the current event, possibly deleting some steps
// and telling the event handler to redo them.
if ( manipulator() )
state = manipulator()->manipulate(eventHandler(), event);
// If the event was not completed, continue generation and continue.
loop = abs(loop) + 1;
} while ( !eventHandler()->empty() );
}
catch (Exception & ex) {
if ( logException(ex, eventHandler()->currentEvent()) ) throw;
}
catch (...) {
event = eventHandler()->currentEvent();
if ( event )
log() << *event;
else
log() << "An exception occurred before any event object was created!";
log() << endl;
dump();
throw;
}
if ( ThePEG_DEBUG_LEVEL ) {
if ( ( ThePEG_DEBUG_LEVEL == Debug::printEveryEvent ||
ieve < printEvent ) && event ) log() << *event;
if ( debugEvent > 0 && ieve + 1 >= debugEvent )
Debug::level = Debug::full;
}
} while ( !event );
// If scheduled, dump a clean state between events
if ( ThePEG_DEBUG_LEVEL && dumpPeriod > 0 && ieve%dumpPeriod == 0 ) {
eventHandler()->clearEvent();
eventHandler()->clean();
dump();
}
return event;
}
EventPtr EventGenerator::doGenerateEvent(tEventPtr e) {
if ( N() >= 0 && ++ieve > N() ) return EventPtr();
EventPtr event = e;
try {
event = eventHandler()->generateEvent(e);
}
catch (Exception & ex) {
if ( logException(ex, eventHandler()->currentEvent()) ) throw;
}
catch (...) {
event = eventHandler()->currentEvent();
if ( !event ) event = e;
log() << *event << endl;
dump();
throw;
}
return event;
}
EventPtr EventGenerator::doGenerateEvent(tStepPtr s) {
if ( N() >= 0 && ++ieve > N() ) return EventPtr();
EventPtr event;
try {
event = eventHandler()->generateEvent(s);
}
catch (Exception & ex) {
if ( logException(ex, eventHandler()->currentEvent()) ) throw;
}
catch (...) {
event = eventHandler()->currentEvent();
if ( event ) log() << *event << endl;
dump();
throw;
}
return event;
}
EventPtr EventGenerator::generateEvent(Event & e) {
UseRandom currentRandom(theRandom);
CurrentGenerator currentGenerator(this);
EventPtr event = doGenerateEvent(tEventPtr(&e));
if ( event ) weightSum += event->weight();
return event;
}
EventPtr EventGenerator::generateEvent(Step & s) {
UseRandom currentRandom(theRandom);
CurrentGenerator currentGenerator(this);
EventPtr event = doGenerateEvent(tStepPtr(&s));
if ( event ) weightSum += event->weight();
return event;
}
Energy EventGenerator::maximumCMEnergy() const {
tcEHPtr eh = eventHandler();
return eh->lumiFnPtr()? eh->lumiFn().maximumCMEnergy(): ZERO;
}
void EventGenerator::doInitialize(bool initOnly) {
if ( !initOnly )
openOutputFiles();
init();
if ( !initOnly )
initrun();
if ( !ThePEG_DEBUG_LEVEL ) Exception::noabort = true;
}
void EventGenerator::doGo(long next, long maxevent, bool tics) {
if ( maxevent >= 0 ) N(maxevent);
if ( next >= 0 ) {
if ( tics )
cerr << "event> " << setw(9) << "init\r" << flush;
initialize();
ieve = next-1;
} else {
openOutputFiles();
}
if ( tics ) tic();
try {
while ( shoot() ) {
if ( tics ) tic();
}
}
catch ( ... ) {
finish();
throw;
}
finish();
finally();
}
void EventGenerator::tic(long currev, long totev) const {
if ( !currev ) currev = ieve;
if ( !totev ) totev = N();
long i = currev;
long n = totev;
bool skip = currev%(max(totev/100, 1L));
if ( i > n/2 ) i = n-i;
while ( skip && i >= 10 && !(i%10) ) i /= 10;
if ( i == 1 || i == 2 || i == 5 ) skip = false;
if (!theIntermediateOutput) { //default
if ( skip ) return;
cerr << "event> " << setw(8) << currev << " " << setw(8) << totev << "\r";
}
else if (theIntermediateOutput) {
if ( skip && currev%10000!=0) return;
cerr << "event> " << setw(9) << right << currev << "/" << totev
<< "; xs = " << integratedXSec()/picobarn << " pb +- "
<< integratedXSecErr()/picobarn << " pb" << endl;
}
cerr.flush();
if ( currev == totev ) cerr << endl;
}
void EventGenerator::dump() const {
if ( dumpPeriod > -1 ) {
string dumpfile;
if ( keepAllDumps ) {
ostringstream number;
number << ieve;
dumpfile = filename() + "-" + number.str() + ".dump";
}
else
dumpfile = filename() + ".dump";
PersistentOStream file(dumpfile, globalLibraries());
file << tcEGPtr(this);
}
}
void EventGenerator::use(const Interfaced & i) {
IBPtr ip = getPtr(i);
if ( ip ) usedObjects.insert(ip);
}
void EventGenerator::generateReferences() {
typedef map<string,string> StringMap;
StringMap references;
// First get all model descriptions and model references from the
// used objects. Put them in a map indexed by the description to
// avoid duplicates.
for ( ObjectSet::iterator it = usedObjects.begin();
it != usedObjects.end(); ++it ) {
if ( *it == strategy() ) continue;
string desc = Repository::getModelDescription(*it);
if ( desc.empty() ) continue;
if ( dynamic_ptr_cast<cEHPtr>(*it) ) desc = "A " + desc;
else if ( dynamic_ptr_cast<cSMPtr>(*it) ) desc = "B " + desc;
else if ( dynamic_ptr_cast<cMEPtr>(*it) ) desc = "C " + desc;
else if ( dynamic_ptr_cast<cCascHdlPtr>(*it) ) desc = "D " + desc;
else if ( dynamic_ptr_cast<cHadrHdlPtr>(*it) ) desc = "E " + desc;
else if ( dynamic_ptr_cast<cStepHdlPtr>(*it) ) desc = "F " + desc;
else if ( dynamic_ptr_cast<cDecayerPtr>(*it) ) desc = "Y " + desc;
else if ( dynamic_ptr_cast<cAnalysisHdlPtr>(*it) ) desc = "Z " + desc;
else if ( dynamic_ptr_cast<Ptr<HandlerBase>::const_pointer>(*it) )
desc = "G " + desc;
else desc = "H " + desc;
references[desc] = Repository::getModelReferences(*it);
}
// Now get the main strategy description which should put first and
// remove it from the map.
string stratdesc;
string stratref;
if ( strategy() ) {
stratdesc = Repository::getModelDescription(strategy());
stratref = Repository::getModelReferences(strategy());
references.erase(stratdesc);
}
// Open the file and write out an appendix header
if ( !useStdout )
reffile().open((filename() + ".tex").c_str());
ref() << "\\documentclass{article}\n"
<< "\\usepackage{graphics}\n"
<< "\\begin{document}\n"
<< "\\appendix\n"
<< "\\section[xxx]{\\textsc{ThePEG} version " << Repository::version()
<< " \\cite{ThePEG} Run Information}\n"
<< "Run name: \\textbf{" << runName()
<< "}:\\\\\n";
if ( !stratdesc.empty() )
ref() << "This run was generated using " << stratdesc
<< " and the following models:\n";
else
ref() << "The following models were used:\n";
ref() << "\\begin{itemize}\n";
// Write out all descriptions.
for ( StringMap::iterator it = references.begin();
it != references.end(); ++it )
ref() << "\\item " << it->first.substr(2) << endl;
// Write out thebibliography header and all references.
ref() << "\\end{itemize}\n\n"
<< "\\begin{thebibliography}{99}\n"
<< "\\bibitem{ThePEG} L.~L\\\"onnblad, "
<< "Comput.~Phys.~Commun.\\ {\\bf 118} (1999) 213.\n";
if ( !stratref.empty() ) ref() << stratref << '\n';
for ( StringMap::iterator it = references.begin();
it != references.end(); ++it )
ref() << it->second << '\n';
ref() << "\\end{thebibliography}\n"
<< "\\end{document}" << endl;
if ( !useStdout )
reffile().close();
}
void EventGenerator::strategy(StrategyPtr s) {
theStrategy = s;
}
int EventGenerator::count(const Exception & ex) {
return ++theExceptions[make_pair(StringUtils::typeName(typeid(ex)),
ex.severity())];
}
void EventGenerator::printException(const Exception & ex) {
switch ( ex.severity() ) {
case Exception::info:
log() << "* An information";
break;
case Exception::warning:
log() << "* A warning";
break;
case Exception::setuperror:
log() << "** A setup";
break;
case Exception::eventerror:
log() << "** An event";
break;
case Exception::runerror:
log() << "*** An run";
break;
case Exception::maybeabort:
case Exception::abortnow:
log() << "**** A serious";
break;
default:
log() << "**** An unknown";
break;
}
if ( ieve > 0 )
log() << " exception of type " << StringUtils::typeName(typeid(ex))
<< " occurred while generating event number "
<< ieve << ": \n" << ex.message() << endl;
else
log() << " exception occurred in the initialization of "
<< name() << ": \n" << ex.message() << endl;
if ( ex.severity() == Exception::eventerror )
log() << "The event will be discarded." << endl;
}
void EventGenerator::logWarning(const Exception & ex) {
if ( ex.severity() != Exception::info &&
ex.severity() != Exception::warning ) throw ex;
ex.handle();
int c = count(ex);
if ( c > maxWarnings ) return;
printException(ex);
if ( c == maxWarnings )
log() << "No more warnings of this kind will be reported." << endl;
}
bool EventGenerator::
logException(const Exception & ex, tcEventPtr event) {
bool noEvent = !event;
ex.handle();
int c = count(ex);
if ( c <= maxWarnings ) {
printException(ex);
if ( c == maxWarnings )
log() << "No more warnings of this kind will be reported." << endl;
}
if ( ex.severity() == Exception::info ||
ex.severity() == Exception::warning ) {
ex.handle();
return false;
}
if ( ex.severity() == Exception::eventerror ) {
if ( c < maxErrors || maxErrors <= 0 ) {
ex.handle();
if ( ThePEG_DEBUG_LEVEL > 0 && !noEvent ) log() << *event;
return false;
}
if ( c > maxErrors ) printException(ex);
log() << "Too many (" << c << ") exceptions of this kind has occurred. "
"Execution will be stopped.\n";
} else {
log() << "This exception is too serious. Execution will be stopped.\n";
}
if ( !noEvent ) log() << *event;
else log()
<< "An exception occurred before any event object was created!\n";
dump();
return true;
}
struct MatcherOrdering {
bool operator()(tcPMPtr m1, tcPMPtr m2) {
return m1->name() < m2->name() ||
( m1->name() == m2->name() && m1->fullName() < m2->fullName() );
}
};
struct ObjectOrdering {
bool operator()(tcIBPtr i1, tcIBPtr i2) {
return i1->fullName() < i2->fullName();
}
};
void EventGenerator::persistentOutput(PersistentOStream & os) const {
set<tcPMPtr,MatcherOrdering> match(theMatchers.begin(), theMatchers.end());
set<tcIBPtr,ObjectOrdering> usedset(usedObjects.begin(), usedObjects.end());
os << theDefaultObjects << theLocalParticles << theStandardModel
<< theStrategy << theRandom << theEventHandler << theAnalysisHandlers
<< theHistogramFactory << theEventManipulator << thePath << theRunName
<< theNumberOfEvents << theObjectMap << theParticles
<< theQuickParticles << theQuickSize << match << usedset
<< ieve << weightSum << theDebugLevel << logNonDefault << printEvent
<< dumpPeriod << keepAllDumps << debugEvent
<< maxWarnings << maxErrors << theCurrentEventHandler
<< theCurrentStepHandler << useStdout << theIntermediateOutput << theMiscStream.str()
<< Repository::listReadDirs();
}
void EventGenerator::persistentInput(PersistentIStream & is, int) {
string dummy;
vector<string> readdirs;
theGlobalLibraries = is.globalLibraries();
is >> theDefaultObjects >> theLocalParticles >> theStandardModel
>> theStrategy >> theRandom >> theEventHandler >> theAnalysisHandlers
>> theHistogramFactory >> theEventManipulator >> thePath >> theRunName
>> theNumberOfEvents >> theObjectMap >> theParticles
>> theQuickParticles >> theQuickSize >> theMatchers >> usedObjects
>> ieve >> weightSum >> theDebugLevel >> logNonDefault >> printEvent
>> dumpPeriod >> keepAllDumps >> debugEvent
>> maxWarnings >> maxErrors >> theCurrentEventHandler
>> theCurrentStepHandler >> useStdout >> theIntermediateOutput >> dummy
>> readdirs;
theMiscStream.str(dummy);
theMiscStream.seekp(0, std::ios::end);
theObjects.clear();
for ( ObjectMap::iterator it = theObjectMap.begin();
it != theObjectMap.end(); ++it ) theObjects.insert(it->second);
Repository::appendReadDir(readdirs);
}
void EventGenerator::setLocalParticles(PDPtr pd, int) {
localParticles()[pd->id()] = pd;
}
void EventGenerator::insLocalParticles(PDPtr pd, int) {
localParticles()[pd->id()] = pd;
}
void EventGenerator::delLocalParticles(int place) {
ParticleMap::iterator it = localParticles().begin();
while ( place-- && it != localParticles().end() ) ++it;
if ( it != localParticles().end() ) localParticles().erase(it);
}
vector<PDPtr> EventGenerator::getLocalParticles() const {
vector<PDPtr> ret;
for ( ParticleMap::const_iterator it = localParticles().begin();
it != localParticles().end(); ++it ) ret.push_back(it->second);
return ret;
}
void EventGenerator::setPath(string newPath) {
if ( std::system(("mkdir -p " + newPath).c_str()) ) throw EGNoPath(newPath);
if ( std::system(("touch " + newPath + "/.ThePEG").c_str()) )
throw EGNoPath(newPath);
if ( std::system(("rm -f " + newPath + "/.ThePEG").c_str()) )
throw EGNoPath(newPath);
thePath = newPath;
}
string EventGenerator::defPath() const {
char * env = std::getenv("ThePEG_RUN_DIR");
if ( env ) return string(env);
return string(".");
}
ostream & EventGenerator::out() {
return theOutStream;
}
ostream & EventGenerator::log() {
return logfile().is_open()? logfile(): BaseRepository::cout();
}
ostream & EventGenerator::ref() {
return reffile().is_open()? reffile(): BaseRepository::cout();
}
string EventGenerator::doSaveRun(string runname) {
runname = StringUtils::car(runname);
if ( runname.empty() ) runname = theRunName;
if ( runname.empty() ) runname = name();
EGPtr eg = Repository::makeRun(this, runname);
string file = eg->filename() + ".run";
PersistentOStream os(file);
os << eg;
if ( !os ) return "Error: Save failed! (I/O error)";
return "";
}
string EventGenerator::doMakeRun(string runname) {
runname = StringUtils::car(runname);
if ( runname.empty() ) runname = theRunName;
if ( runname.empty() ) runname = name();
Repository::makeRun(this, runname);
return "";
}
bool EventGenerator::preinitRegister(IPtr obj, string fullname) {
if ( !preinitializing ) throw InitException()
<< "Tried to register a new object in the initialization of an "
<< "EventGenerator outside of the pre-initialization face. "
<< "The preinitRegister() can only be called from a doinit() function "
<< "in an object for which preInitialize() returns true.";
if ( objectMap().find(fullname) != objectMap().end() ) return false;
obj->name(fullname);
objectMap()[fullname] = obj;
objects().insert(obj);
obj->theGenerator = this;
PDPtr pd = dynamic_ptr_cast<PDPtr>(obj);
if ( pd ) theParticles[pd->id()] = pd;
PMPtr pm = dynamic_ptr_cast<PMPtr>(obj);
if ( pm ) theMatchers.insert(pm);
return true;
}
bool EventGenerator::preinitRemove(IPtr obj) {
bool deleted=true;
if(theObjects.find(obj)!=theObjects.end()) {
theObjects.erase(obj);
}
else
deleted = false;
if(theObjectMap.find(obj->fullName())!=theObjectMap.end())
theObjectMap.erase(obj->fullName());
else
deleted = false;
return deleted;
}
IPtr EventGenerator::
preinitCreate(string classname, string fullname, string libraries) {
if ( !preinitializing ) throw InitException()
<< "Tried to create a new object in the initialization of an "
<< "EventGenerator outside of the pre-initialization face. "
<< "The preinitCreate() can only be called from a doinit() function "
<< "in an object for which preInitialize() returns true.";
if ( objectMap().find(fullname) != objectMap().end() ) return IPtr();
const ClassDescriptionBase * db = DescriptionList::find(classname);
while ( !db && libraries.length() ) {
string library = StringUtils::car(libraries);
libraries = StringUtils::cdr(libraries);
DynamicLoader::load(library);
db = DescriptionList::find(classname);
}
if ( !db ) return IPtr();
IPtr obj = dynamic_ptr_cast<IPtr>(db->create());
if ( !obj ) return IPtr();
if ( !preinitRegister(obj, fullname) ) return IPtr();
return obj;
}
string EventGenerator::
preinitInterface(IPtr obj, string ifcname, string cmd, string value) {
if ( !preinitializing ) throw InitException()
<< "Tried to manipulate an external object in the initialization of an "
<< "EventGenerator outside of the pre-initialization face. "
<< "The preinitSet() can only be called from a doinit() function "
<< "in an object for which preInitialize() returns true.";
if ( !obj ) return "Error: No object found.";
const InterfaceBase * ifc = Repository::FindInterface(obj, ifcname);
if ( !ifc ) return "Error: No such interface found.";
try {
return ifc->exec(*obj, cmd, value);
}
catch ( const InterfaceException & ex) {
ex.handle();
return "Error: " + ex.message();
}
}
string EventGenerator::
preinitInterface(IPtr obj, string ifcname, int index,
string cmd, string value) {
ostringstream os;
os << index;
return preinitInterface(obj, ifcname, cmd, os.str() + " " + value);
}
string EventGenerator::
preinitInterface(string fullname, string ifcname, string cmd, string value) {
return preinitInterface(getObject<Interfaced>(fullname), ifcname, cmd, value);
}
string EventGenerator::
preinitInterface(string fullname, string ifcname, int index,
string cmd, string value) {
return preinitInterface(getObject<Interfaced>(fullname), ifcname, index,
cmd, value);
}
tDMPtr EventGenerator::findDecayMode(string tag) const {
for ( ObjectSet::const_iterator it = objects().begin();
it != objects().end(); ++it ) {
tDMPtr dm = dynamic_ptr_cast<tDMPtr>(*it);
if ( dm && dm->tag() == tag ) return dm;
}
return tDMPtr();
}
tDMPtr EventGenerator::preinitCreateDecayMode(string tag) {
return constructDecayMode(tag);
}
DMPtr EventGenerator::constructDecayMode(string & tag) {
DMPtr rdm;
DMPtr adm;
int level = 0;
string::size_type end = 0;
while ( end < tag.size() && ( tag[end] != ']' || level ) ) {
switch ( tag[end++] ) {
case '[':
++level;
break;
case ']':
--level;
break;
}
}
rdm = findDecayMode(tag.substr(0,end));
if ( rdm ) return rdm;
string::size_type next = tag.find("->");
if ( next == string::npos ) return rdm;
if ( tag.find(';') == string::npos ) return rdm;
tPDPtr pd = getObject<ParticleData>(tag.substr(0,next));
if ( !pd ) pd = findParticle(tag.substr(0,next));
if ( !pd ) return rdm;
rdm = ptr_new<DMPtr>();
rdm->parent(pd);
if ( pd->CC() ) {
adm = ptr_new<DMPtr>();
adm->parent(pd->CC());
rdm->theAntiPartner = adm;
adm->theAntiPartner = rdm;
}
bool error = false;
tag = tag.substr(next+2);
tPDPtr lastprod;
bool dolink = false;
do {
switch ( tag[0] ) {
case '[':
{
tag = tag.substr(1);
tDMPtr cdm = constructDecayMode(tag);
if ( cdm ) rdm->addCascadeProduct(cdm);
else error = true;
} break;
case '=':
dolink = true;
+ [[fallthrough]];
case ',':
case ']':
tag = tag.substr(1);
break;
case '?':
{
next = min(tag.find(','), tag.find(';'));
tPMPtr pm = findMatcher(tag.substr(1,next-1));
if ( pm ) rdm->addProductMatcher(pm);
else error = true;
tag = tag.substr(next);
} break;
case '!':
{
next = min(tag.find(','), tag.find(';'));
tPDPtr pd = findParticle(tag.substr(1,next-1));
if ( pd ) rdm->addExcluded(pd);
else error = true;
tag = tag.substr(next);
} break;
case '*':
{
next = min(tag.find(','), tag.find(';'));
tPMPtr pm = findMatcher(tag.substr(1,next-1));
if ( pm ) rdm->setWildMatcher(pm);
else error = true;
tag = tag.substr(next);
} break;
default:
{
next = min(tag.find('='), min(tag.find(','), tag.find(';')));
tPDPtr pdp = findParticle(tag.substr(0,next));
if ( pdp ) rdm->addProduct(pdp);
else error = true;
tag = tag.substr(next);
if ( dolink && lastprod ) {
rdm->addLink(lastprod, pdp);
dolink = false;
}
lastprod = pdp;
} break;
}
} while ( tag[0] != ';' && tag.size() );
if ( tag[0] != ';' || error ) {
return DMPtr();
}
tag = tag.substr(1);
DMPtr ndm = findDecayMode(rdm->tag());
if ( ndm ) return ndm;
pd->addDecayMode(rdm);
if ( !preinitRegister(rdm, pd->fullName() + "/" + rdm->tag()) )
return DMPtr();
if ( adm ) {
preinitRegister(adm, pd->CC()->fullName() + "/" + adm->tag());
rdm->CC(adm);
adm->CC(rdm);
}
return rdm;
}
tPDPtr EventGenerator::findParticle(string pdgname) const {
for ( ParticleMap::const_iterator it = particles().begin();
it != particles().end(); ++it )
if ( it->second->PDGName() == pdgname ) return it->second;
return tPDPtr();
}
tPMPtr EventGenerator::findMatcher(string name) const {
for ( MatcherSet::const_iterator it = matchers().begin();
it != matchers().end(); ++it )
if ( (**it).name() == name ) return *it;
return tPMPtr();
}
ClassDescription<EventGenerator> EventGenerator::initEventGenerator;
void EventGenerator::Init() {
static ClassDocumentation<EventGenerator> documentation
("This is the main class used to administer an event generation run. "
"The actual generation of each event is handled by the assigned "
"<interface>EventHandler</interface> object. When the event generator"
"is properly set up it can be initialized with the command "
"<interface>MakeRun</interface> and/or saved to a file with the command "
"<interface>SaveRun</interface>. If saved to a file, the event generator "
"can be read into another program to produce events. The file can also "
"be read into the <tt>runThePEG</tt> program where a number of events "
"determined by the parameter <interface>NumberOfEvents</interface> is "
"generated with each event analysed by the list of assigned "
"<interface>AnalysisHandlers</interface>.");
static Reference<EventGenerator,StandardModelBase> interfaceStandardModel
("StandardModelParameters",
"The ThePEG::StandardModelBase object to be used to access standard "
"model parameters in this run.",
&EventGenerator::theStandardModel, false, false, true, false);
static Reference<EventGenerator,EventHandler> interfaceEventHandler
("EventHandler",
"The ThePEG::EventHandler object to be used to generate the "
"individual events in this run.",
&EventGenerator::theEventHandler, false, false, true, false);
static RefVector<EventGenerator,AnalysisHandler> interfaceAnalysisHandlers
("AnalysisHandlers",
"ThePEG::AnalysisHandler objects to be used to analyze the produced "
"events in this run.",
&EventGenerator::theAnalysisHandlers, 0, true, false, true, false);
static Reference<EventGenerator,FactoryBase> interfaceHistogramFactory
("HistogramFactory",
"An associated factory object for handling histograms to be used by "
"<interface>AnalysisHandlers</interface>.",
&EventGenerator::theHistogramFactory, true, false, true, true, true);
static Reference<EventGenerator,EventManipulator> interfaceEventManip
("EventManipulator",
"An ThePEG::EventManipulator called each time the generation of an "
"event is stopped. The ThePEG::EventManipulator object is able to "
"manipulate the generated event, as opposed to an "
"ThePEG::AnalysisHandler which may only look at the event.",
&EventGenerator::theEventManipulator, true, false, true, true);
static RefVector<EventGenerator,ParticleData> interfaceLocalParticles
("LocalParticles",
"Special versions of ThePEG::ParticleData objects to be used "
"in this run. Note that to delete an object, its number in the list "
"should be given, rather than its id number.",
0, 0, false, false, true, false,
&EventGenerator::setLocalParticles, &EventGenerator::insLocalParticles,
&EventGenerator::delLocalParticles, &EventGenerator::getLocalParticles);
static RefVector<EventGenerator,Interfaced> interfaceDefaultObjects
("DefaultObjects",
"A vector of pointers to default objects. In a ThePEG::Reference or "
"ThePEG::RefVector interface with the defaultIfNull() flag set, if a "
"null pointer is encountered this vector is gone through until an "
"acceptable object is found in which case the null pointer is replaced "
"by a pointer to this object.",
&EventGenerator::theDefaultObjects, 0, true, false, true, false, false);
static Reference<EventGenerator,Strategy> interfaceStrategy
("Strategy",
"An ThePEG::Strategy with additional ThePEG::ParticleData objects to "
"be used in this run.",
&EventGenerator::theStrategy, false, false, true, true);
static Reference<EventGenerator,RandomGenerator> interfaceRandomGenerator
("RandomNumberGenerator",
"An ThePEG::RandomGenerator object which should typically interaface to "
"a CLHEP Random object. This will be the default random number generator "
"for the run, but individual objects may use their own random generator "
"if they wish.",
&EventGenerator::theRandom, true, false, true, false);
static Parameter<EventGenerator,string> interfacePath
("Path",
"The directory where the output files are put.",
&EventGenerator::thePath, ".", true, false,
&EventGenerator::setPath, 0, &EventGenerator::defPath);
interfacePath.directoryType();
static Parameter<EventGenerator,string> interfaceRunName
("RunName",
"The name of this run. This name will be used in the output filenames. "
"The files wil be placed in the directory specified by the "
"<interface>Path</interface> parameter"
"If empty the name of the event generator will be used instead.",
&EventGenerator::theRunName, "", true, false,
0, 0, &EventGenerator::name);
static Parameter<EventGenerator,long> interfaceNumberOfEvents
("NumberOfEvents",
"The number of events to be generated in this run. If less than zero, "
"the number of events is unlimited",
&EventGenerator::theNumberOfEvents, 1000, -1, Constants::MaxInt,
true, false, Interface::lowerlim);
static Parameter<EventGenerator,int> interfaceDebugLevel
("DebugLevel",
"The level of debug information sent out to the log file in the run. "
"Level 0 only gives a limited ammount of warnings and error messages. "
"Level 1 will print the first few events. "
"Level 5 will print every event. "
"Level 9 will print every step in every event.",
&EventGenerator::theDebugLevel, 0, 0, 9, true, false, true);
static Parameter<EventGenerator,int> interfacePrintEvent
("PrintEvent",
"If the debug level is above zero, print the first 'PrintEvent' events.",
&EventGenerator::printEvent, 0, 0, 1000, true, false, Interface::lowerlim);
static Parameter<EventGenerator,long> interfaceDumpPeriod
("DumpPeriod",
"If the debug level is above zero, dump the full state of the run every "
"'DumpPeriod' events. Set it to -1 to disable dumping even in the case of errors.",
&EventGenerator::dumpPeriod, 0, -1, Constants::MaxInt,
true, false, Interface::lowerlim);
static Switch<EventGenerator,bool> interfaceKeepAllDumps
("KeepAllDumps",
"Whether all dump files should be kept, labelled by event number.",
&EventGenerator::keepAllDumps, false, true, false);
static SwitchOption interfaceKeepAllDumpsYes
(interfaceKeepAllDumps,
"Yes",
"Keep all dump files, labelled by event number.",
true);
static SwitchOption interfaceKeepAllDumpsNo
(interfaceKeepAllDumps,
"No",
"Keep only the latest dump file.",
false);
static Parameter<EventGenerator,long> interfaceDebugEvent
("DebugEvent",
"If the debug level is above zero, step up to the highest debug level "
"befor event number 'DebugEvent'.",
&EventGenerator::debugEvent, 0, 0, Constants::MaxInt,
true, false, Interface::lowerlim);
static Parameter<EventGenerator,int> interfaceMaxWarnings
("MaxWarnings",
"The maximum number of warnings of each type which will be printed.",
&EventGenerator::maxWarnings,
10, 1, 100, true, false, Interface::lowerlim);
static Parameter<EventGenerator,int> interfaceMaxErrors
("MaxErrors",
"The maximum number of errors of each type which will be tolerated. "
"If more errors are reported, the run will be aborted.",
&EventGenerator::maxErrors,
10, -1, 100000, true, false, Interface::lowerlim);
static Parameter<EventGenerator,long> interfaceQuickSize
("QuickSize",
"The max absolute id number of particle data objects which are accessed "
"quickly through a vector indexed by the id number.",
&EventGenerator::theQuickSize,
7000, 0, 50000, true, false, Interface::lowerlim);
static Command<EventGenerator> interfaceSaveRun
("SaveRun",
"Isolate, initialize and save this event generator to a file, from which "
"it can be read in and run in another program. If an agument is given "
"this is used as the run name, otherwise the run name is taken from the "
"<interface>RunName</interface> parameter.",
&EventGenerator::doSaveRun, true);
static Command<EventGenerator> interfaceMakeRun
("MakeRun",
"Isolate and initialize this event generator and give it a run name. "
"If no argument is given, the run name is taken from the "
"<interface>RunName</interface> parameter.",
&EventGenerator::doMakeRun, true);
interfaceEventHandler.rank(11.0);
interfaceSaveRun.rank(10.0);
interfaceMakeRun.rank(9.0);
interfaceRunName.rank(8.0);
interfaceNumberOfEvents.rank(7.0);
interfaceAnalysisHandlers.rank(6.0);
static Switch<EventGenerator,bool> interfaceUseStdout
("UseStdout",
"Redirect the logging and output to stdout instead of files.",
&EventGenerator::useStdout, false, true, false);
static SwitchOption interfaceUseStdoutYes
(interfaceUseStdout,
"Yes",
"Use stdout instead of log files.",
true);
static SwitchOption interfaceUseStdoutNo
(interfaceUseStdout,
"No",
"Use log files.",
false);
static Switch<EventGenerator,int> interfaceLogNonDefault
("LogNonDefault",
"Controls the printout of important interfaces which has been changed from their default values.",
&EventGenerator::logNonDefault, -1, true, false);
static SwitchOption interfaceLogNonDefaultYes
(interfaceLogNonDefault,
"Yes",
"Always print changed interfaces.",
1);
static SwitchOption interfaceLogNonDefaultOnDebug
(interfaceLogNonDefault,
"OnDebug",
"Only print changed interfaces if debugging is turned on.",
0);
static SwitchOption interfaceLogNonDefaultNo
(interfaceLogNonDefault,
"No",
"Don't print changed interfaces.",
-1);
interfaceLogNonDefault.setHasDefault(false);
static Switch<EventGenerator,bool> interfaceIntermediateOutput
("IntermediateOutput",
"Modified event number count with the number of events processed so far, "
"which updates at least every 10000 events, together with the corresponding "
"intermediate estimate for the cross section plus the integration error.",
&EventGenerator::theIntermediateOutput, false, true, false);
static SwitchOption interfaceIntermediateOutputYes
(interfaceIntermediateOutput,
"Yes",
"Show the modified event number count with the number of events processed so far, "
"plus further information on the intermediate cross section estimate.",
true);
static SwitchOption interfaceIntermediateOutputNo
(interfaceIntermediateOutput,
"No",
"Show the usual event number count with the number of events processed so far, "
"but no further information on the intermediate cross section estimate.",
false);
}
EGNoPath::EGNoPath(string path) {
theMessage << "Cannot set the directory path for output files to '" << path
<< "' because the directory did not exist and could not be "
<< "created.";
severity(warning);
}