Page MenuHomeHEPForge
Files F8309033

No OneTemporary
Actions

View Options

diff --git a/Handlers/StandardXComb.cc b/Handlers/StandardXComb.cc
--- a/Handlers/StandardXComb.cc
+++ b/Handlers/StandardXComb.cc
@@ -1,701 +1,798 @@
// -*- C++ -*-
//
// StandardXComb.cc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2011 Leif Lonnblad
// Copyright (C) 2009-2011 Simon Platzer
//
// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the StandardXComb class.
//
#include "StandardXComb.h"
#include "StdXCombGroup.h"
#include "ThePEG/Handlers/StandardEventHandler.h"
#include "ThePEG/Handlers/SubProcessHandler.h"
#include "ThePEG/Cuts/Cuts.h"
#include "ThePEG/PDF/PartonExtractor.h"
#include "ThePEG/Utilities/Debug.h"
#include "ThePEG/Utilities/Maths.h"
#include "ThePEG/PDT/ParticleData.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Utilities/SimplePhaseSpace.h"
#include "ThePEG/Utilities/UtilityBase.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/EventRecord/SubProcessGroup.h"
#include "ThePEG/Vectors/LorentzRotation.h"
#include "ThePEG/MatrixElement/MEBase.h"
#include "ThePEG/MatrixElement/ColourLines.h"
#include "ThePEG/Handlers/LuminosityFunction.h"
#include "ThePEG/Handlers/CascadeHandler.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/EventRecord/TmpTransform.h"
#ifdef ThePEG_TEMPLATES_IN_CC_FILE
#include "StandardXComb.tcc"
#endif
using namespace ThePEG;
StandardXComb::StandardXComb()
: XComb(), isMirror(false), theNDim(0),
partonDims(make_pair(0, 0)), theKinematicsGenerated(false),
theLastDiagramIndex(0), theLastPDFWeight(0.0),
theLastCrossSection(ZERO), theLastJacobian(1.0), theLastME2(-1.0), theLastPreweight(1.0),
theLastMECrossSection(ZERO), theLastMEPDFWeight(1.0), theLastMECouplings(1.0),
- checkedCuts(false), passedCuts(false), theCutWeight(1.0) {}
+ checkedCuts(false), passedCuts(false), theCutWeight(1.0), theNeedsReshuffling(false) {}
StandardXComb::
StandardXComb(Energy newMaxEnergy, const cPDPair & inc,
tEHPtr newEventHandler, tSubHdlPtr newSubProcessHandler,
tPExtrPtr newExtractor, tCascHdlPtr newCKKW,
const PBPair & newPartonBins, tCutsPtr newCuts,
tMEPtr newME, const DiagramVector & newDiagrams, bool mir,
tStdXCombPtr newHead)
: XComb(newMaxEnergy, inc, newEventHandler,
newExtractor, newCKKW, newPartonBins, newCuts),
theSubProcessHandler(newSubProcessHandler), theME(newME),
- theDiagrams(newDiagrams), isMirror(mir), theKinematicsGenerated(false),
- theLastDiagramIndex(0), theLastPDFWeight(0.0),
- theLastCrossSection(ZERO), theLastME2(-1.0), theLastPreweight(1.0), theLastMECrossSection(ZERO),
+ theDiagrams(newDiagrams), isMirror(mir), theNDim(0), partonDims(0,0),
+ theKinematicsGenerated(false), theLastDiagramIndex(0), theLastPDFWeight(0.0),
+ theLastCrossSection(ZERO), theLastJacobian(0.0), theLastME2(-1.0),
+ theLastPreweight(1.0), theLastMECrossSection(ZERO),
theLastMEPDFWeight(1.0), theLastMECouplings(1.0), theHead(newHead),
- checkedCuts(false), passedCuts(false) {
+ checkedCuts(false), passedCuts(false), theCutWeight(1.0), theNeedsReshuffling(false) {
partonDims = pExtractor()->nDims(partonBins());
if ( matrixElement()->haveX1X2() ) {
partonDims.first = 0;
partonDims.second = 0;
}
theNDim = matrixElement()->nDim() + partonDims.first + partonDims.second;
mePartonData() = lastDiagram()->partons();
+ checkReshufflingNeeds();
}
StandardXComb::StandardXComb(tMEPtr me, const tPVector & parts,
DiagramIndex indx)
: theME(me), isMirror(false), theNDim(0), partonDims(make_pair(0, 0)),
theKinematicsGenerated(false),
theLastDiagramIndex(0), theLastPDFWeight(0.0), theLastCrossSection(ZERO),
theLastME2(-1.0), theLastPreweight(1.0), theLastMECrossSection(ZERO), theLastMEPDFWeight(1.0),
theLastMECouplings(1.0),
- checkedCuts(false), passedCuts(false) {
+ checkedCuts(false), passedCuts(false), theCutWeight(1.0), theNeedsReshuffling(false) {
subProcess(new_ptr(SubProcess(make_pair(parts[0], parts[1]),
tCollPtr(), me)));
for ( int i = 0, N = parts.size(); i < N; ++i ) {
subProcess()->addOutgoing(parts[i], false);
theMEPartonData.push_back(parts[i]->dataPtr());
theMEMomenta.push_back(parts[i]->momentum());
}
lastSHat((meMomenta()[0] + meMomenta()[1]).m2());
string tag = me->diagrams()[indx]->getTag();
for ( int i = 0, N = me->diagrams().size(); i < N; ++i )
if ( me->diagrams()[i]->getTag() == tag )
theDiagrams.push_back(me->diagrams()[i]);
+ checkReshufflingNeeds();
}
StandardXComb::StandardXComb(tStdXCombPtr newHead,
const PBPair & newPartonBins, tMEPtr newME,
const DiagramVector & newDiagrams)
: XComb(newHead->maxEnergy(), newHead->particles(),
newHead->eventHandlerPtr(), newHead->pExtractor(),
newHead->CKKWHandler(), newPartonBins, newHead->cuts()),
theSubProcessHandler(const_ptr_cast<tSubHdlPtr>(newHead->subProcessHandler())),
- theME(newME), theDiagrams(newDiagrams), isMirror(newHead->mirror()),
- theKinematicsGenerated(false), theLastDiagramIndex(0), theLastPDFWeight(0.0),
- theLastCrossSection(ZERO), theLastME2(-1.0), theLastPreweight(1.0), theLastMECrossSection(ZERO),
+ theME(newME), theDiagrams(newDiagrams), isMirror(newHead->mirror()), theNDim(0),
+ partonDims(0,0), theKinematicsGenerated(false), theLastDiagramIndex(0), theLastPDFWeight(0.0),
+ theLastCrossSection(ZERO), theLastJacobian(0.0), theLastME2(-1.0),
+ theLastPreweight(1.0), theLastMECrossSection(ZERO),
theLastMEPDFWeight(1.0), theLastMECouplings(1.0), theHead(newHead),
- checkedCuts(false), passedCuts(false) {
+ checkedCuts(false), passedCuts(false), theCutWeight(1.0), theNeedsReshuffling(false) {
partonDims = pExtractor()->nDims(partonBins());
if ( matrixElement()->haveX1X2() ) {
partonDims.first = 0;
partonDims.second = 0;
}
theNDim = matrixElement()->nDim() + partonDims.first + partonDims.second;
mePartonData() = lastDiagram()->partons();
+ checkReshufflingNeeds();
/* // wait for C++11
StandardXComb(newHead->maxEnergy(),newHead->particles(),
newHead->eventHandlerPtr(),
const_ptr_cast<tSubHdlPtr>(newHead->subProcessHandler()),
newHead->pExtractor(),newHead->CKKWHandler(),
newPartonBins,newHead->cuts(),newME,newDiagrams,newHead->mirror(),
newHead);
*/
}
StandardXComb::~StandardXComb() {}
void StandardXComb::recreatePartonBinInstances(Energy2 scale) {
PBIPair newBins;
Direction<0> dir(true);
newBins.first =
new_ptr(PartonBinInstance(lastPartons().first,partonBins().first,scale));
dir.reverse();
newBins.second =
new_ptr(PartonBinInstance(lastPartons().second,partonBins().second,scale));
resetPartonBinInstances(newBins);
setPartonBinInfo();
lastPartons().first->scale(partonBinInstances().first->scale());
lastPartons().second->scale(partonBinInstances().second->scale());
}
void StandardXComb::refillPartonBinInstances(const double* r) {
pExtractor()->select(this);
pExtractor()->updatePartonBinInstances(partonBinInstances());
pExtractor()->generateSHat(lastS(), partonBinInstances(),
r, r + nDim() - partonDims.second,true);
}
void StandardXComb::setIncomingPartons(tStdXCombPtr labHead) {
if ( lastPartons().first )
return;
if ( !labHead )
labHead = head();
createPartonBinInstances();
lastParticles(labHead->lastParticles());
setPartonBinInfo();
lastPartons(make_pair(mePartonData()[0]->produceParticle(Lorentz5Momentum()),
mePartonData()[1]->produceParticle(Lorentz5Momentum())));
Lorentz5Momentum pFirst = meMomenta()[0];
Lorentz5Momentum pSecond = meMomenta()[1];
if ( labHead->matrixElement()->wantCMS() ) {
Boost toLab = (labHead->lastPartons().first->momentum() +
labHead->lastPartons().second->momentum()).boostVector();
if ( toLab.mag2() > Constants::epsilon ) {
pFirst.boost(toLab);
pSecond.boost(toLab);
}
}
lastPartons().first->set5Momentum(pFirst);
lastPartons().second->set5Momentum(pSecond);
partonBinInstances().first->parton(lastPartons().first);
partonBinInstances().second->parton(lastPartons().second);
lastS((lastParticles().first->momentum() +
lastParticles().second->momentum()).m2());
lastSHat((lastPartons().first->momentum() +
lastPartons().second->momentum()).m2());
lastP1P2(make_pair(labHead->lastP1(),labHead->lastP2()));
double x1 =
lastPartons().first->momentum().plus()/
lastParticles().first->momentum().plus();
double x2 =
lastPartons().second->momentum().minus()/
lastParticles().second->momentum().minus();
lastX1X2(make_pair(x1,x2));
lastY((lastPartons().first->momentum()+
lastPartons().second->momentum()).rapidity());
}
void StandardXComb::fill(const PPair& newParticles,
const PPair& newPartons,
const vector<Lorentz5Momentum>& newMEMomenta,
const DVector& newLastRandomNumbers) {
lastParticles(newParticles);
lastP1P2(make_pair(0.0, 0.0));
lastS((lastParticles().first->momentum() +
lastParticles().second->momentum()).m2());
lastPartons(newPartons);
lastSHat((lastPartons().first->momentum() +
lastPartons().second->momentum()).m2());
lastX1X2(make_pair(lastPartons().first->momentum().plus()/
lastParticles().first->momentum().plus(),
lastPartons().second->momentum().minus()/
lastParticles().second->momentum().minus()));
lastY((lastPartons().first->momentum() +
lastPartons().second->momentum()).rapidity());
meMomenta().resize(newMEMomenta.size());
copy(newMEMomenta.begin(),newMEMomenta.end(),
meMomenta().begin());
lastRandomNumbers().resize(newLastRandomNumbers.size());
copy(newLastRandomNumbers.begin(),newLastRandomNumbers.end(),
lastRandomNumbers().begin());
}
bool StandardXComb::checkInit() {
Energy summin = ZERO;
for ( int i = 2, N = mePartonData().size(); i < N; ++i ) {
summin += mePartonData()[i]->massMin();
}
return ( summin < min(maxEnergy(), cuts()->mHatMax()) );
}
bool StandardXComb::willPassCuts() {
if ( checkedCuts )
return passedCuts;
checkedCuts = true;
if ( !head() ) {
if ( !cuts()->initSubProcess(lastSHat(), lastY(), mirror()) ) {
passedCuts = false;
theCutWeight = 0.0;
return false;
}
} else {
cuts()->initSubProcess(lastSHat(), lastY(), mirror());
}
// check for exceptional configurations which may happen in NLO
// dijets subtraction with an extremely soft incoming parton giving
// rise to about lightlike CM momentum
if ( (meMomenta()[0]+meMomenta()[1]).m2() <= ZERO ) {
passedCuts = false;
return false;
}
tcPDVector outdata(mePartonData().begin()+2,mePartonData().end());
vector<LorentzMomentum> outmomenta(meMomenta().begin()+2,meMomenta().end());
Boost tocm = (meMomenta()[0]+meMomenta()[1]).findBoostToCM();
if ( tocm.mag2() > Constants::epsilon ) {
for ( vector<LorentzMomentum>::iterator p = outmomenta.begin();
p != outmomenta.end(); ++p ) {
p->boost(tocm);
}
}
if ( !cuts()->passCuts(outdata,outmomenta,mePartonData()[0],mePartonData()[1]) ) {
theCutWeight = cuts()->cutWeight();
passedCuts = false;
return false;
}
theCutWeight = cuts()->cutWeight();
passedCuts = true;
return true;
}
void StandardXComb::clean() {
XComb::clean();
theLastPDFWeight = 0.0;
theLastCrossSection = ZERO;
theLastJacobian = 0.0;
theLastME2 = 0.0;
theLastPreweight = 1.0;
theLastMECrossSection = ZERO;
theLastMEPDFWeight = 0.0;
theLastMECouplings = 0.0;
theProjectors.clear();
theProjector = StdXCombPtr();
theKinematicsGenerated = false;
checkedCuts = false;
passedCuts = false;
theCutWeight = 1.0;
matrixElement()->flushCaches();
}
CrossSection StandardXComb::dSigDR(const double * r) {
matrixElement()->setXComb(this);
if ( !matrixElement()->apply() ) {
subProcess(SubProPtr());
lastCrossSection(ZERO);
return ZERO;
}
meMomenta().resize(mePartonData().size());
if ( !matrixElement()->generateKinematics(r) ) {
subProcess(SubProPtr());
lastCrossSection(ZERO);
return ZERO;
}
setIncomingPartons();
lastScale(matrixElement()->scale());
lastAlphaS(matrixElement()->alphaS());
lastAlphaEM(matrixElement()->alphaEM());
if ( (!willPassCuts() &&
!matrixElement()->headCuts() &&
!matrixElement()->ignoreCuts()) ||
!matrixElement()->apply() ) {
subProcess(SubProPtr());
lastCrossSection(ZERO);
return ZERO;
}
lastPDFWeight(head()->lastPDFWeight());
matrixElement()->setKinematics();
CrossSection xsec = matrixElement()->dSigHatDR() * lastPDFWeight();
xsec *= cutWeight();
subProcess(SubProPtr());
lastCrossSection(xsec);
return xsec;
}
CrossSection StandardXComb::
dSigDR(const pair<double,double> ll, int nr, const double * r) {
if ( matrixElement()->keepRandomNumbers() ) {
lastRandomNumbers().resize(nDim());
copy(r,r+nDim(),lastRandomNumbers().begin());
}
pExtractor()->select(this);
setPartonBinInfo();
lastP1P2(ll);
lastS(sqr(maxEnergy())/exp(lastP1() + lastP2()));
meMomenta().resize(mePartonData().size());
matrixElement()->setXComb(this);
PPair partons;
if ( !matrixElement()->haveX1X2() ) {
if ( !pExtractor()->generateL(partonBinInstances(),
r, r + nr - partonDims.second) ) {
lastCrossSection(ZERO);
return ZERO;
}
partons = make_pair(partonBinInstances().first->parton(),
partonBinInstances().second->parton());
lastSHat(lastS()/exp(partonBinInstances().first->l() +
partonBinInstances().second->l()));
meMomenta()[0] = partons.first->momentum();
meMomenta()[1] = partons.second->momentum();
} else {
if ( !matrixElement()->generateKinematics(r + partonDims.first) ) {
lastCrossSection(ZERO);
return ZERO;
}
lastSHat((meMomenta()[0]+meMomenta()[1]).m2());
matrixElement()->setKinematics();
lastScale(matrixElement()->scale());
partons.first = mePartonData()[0]->produceParticle(meMomenta()[0]);
partons.second = mePartonData()[1]->produceParticle(meMomenta()[1]);
Direction<0> dir(true);
partonBinInstances().first =
new_ptr(PartonBinInstance(lastParticles().first,partons.first,
partonBins().first,lastScale()));
dir.reverse();
partonBinInstances().second =
new_ptr(PartonBinInstance(lastParticles().second,partons.second,
partonBins().second,lastScale()));
}
lastPartons(partons);
if ( lastSHat() < cuts()->sHatMin() ) {
lastCrossSection(ZERO);
return ZERO;
}
lastY(0.5*(partonBinInstances().second->l() -
partonBinInstances().first->l()));
if ( !cuts()->initSubProcess(lastSHat(), lastY(), mirror()) ) {
lastCrossSection(ZERO);
return ZERO;
}
if ( mirror() ) swap(meMomenta()[0], meMomenta()[1]);
if ( matrixElement()->wantCMS() &&
!matrixElement()->haveX1X2() )
SimplePhaseSpace::CMS(meMomenta()[0], meMomenta()[1], lastSHat());
Energy summ = ZERO;
if ( meMomenta().size() == 3 ) {
if ( !matrixElement()->haveX1X2() )
meMomenta()[2] = Lorentz5Momentum(sqrt(lastSHat()));
} else {
for ( int i = 2, N = meMomenta().size(); i < N; ++i ) {
if ( !matrixElement()->haveX1X2() )
meMomenta()[i] = Lorentz5Momentum(mePartonData()[i]->mass());
summ += mePartonData()[i]->massMin();
}
if ( sqr(summ) >= lastSHat() ) {
lastCrossSection(ZERO);
return ZERO;
}
}
if ( !matrixElement()->haveX1X2() )
lastScale(max(lastSHat()/4.0, cuts()->scaleMin()));
lastSHat(pExtractor()->generateSHat(lastS(), partonBinInstances(),
r, r + nr - partonDims.second,
matrixElement()->haveX1X2()));
if ( !cuts()->sHat(lastSHat()) ) {
lastCrossSection(ZERO);
return ZERO;
}
r += partonDims.first;
lastX1X2(make_pair(lastPartons().first->momentum().plus()/
lastParticles().first->momentum().plus(),
lastPartons().second->momentum().minus()/
lastParticles().second->momentum().minus()));
if ( !cuts()->x1(lastX1()) || !cuts()->x2(lastX2()) ) {
lastCrossSection(ZERO);
return ZERO;
}
lastY((lastPartons().first->momentum() +
lastPartons().second->momentum()).rapidity());
if ( !cuts()->yHat(lastY()) ) {
lastCrossSection(ZERO);
return ZERO;
}
if ( !cuts()->initSubProcess(lastSHat(), lastY(), mirror()) ) {
lastCrossSection(ZERO);
return ZERO;
}
meMomenta()[0] = lastPartons().first->momentum();
meMomenta()[1] = lastPartons().second->momentum();
if ( mirror() ) swap(meMomenta()[0], meMomenta()[1]);
if ( matrixElement()->wantCMS() &&
!matrixElement()->haveX1X2() )
SimplePhaseSpace::CMS(meMomenta()[0], meMomenta()[1], lastSHat());
if ( meMomenta().size() == 3 ) {
if ( !matrixElement()->haveX1X2() )
meMomenta()[2] = Lorentz5Momentum(sqrt(lastSHat()));
} else {
if ( sqr(summ) >= lastSHat() ) {
lastCrossSection(ZERO);
return ZERO;
}
}
if ( !matrixElement()->haveX1X2() ) {
if ( !matrixElement()->generateKinematics(r) ) {
lastCrossSection(ZERO);
return ZERO;
}
}
lastScale(matrixElement()->scale());
if ( !cuts()->scale(lastScale()) ) {
lastCrossSection(ZERO);
return ZERO;
}
// get information on cuts; we don't take this into account here for
// reasons of backward compatibility but this will change eventually
willPassCuts();
pair<bool,bool> evalPDFS =
make_pair(matrixElement()->havePDFWeight1(),
matrixElement()->havePDFWeight2());
if ( mirror() )
swap(evalPDFS.first,evalPDFS.second);
lastPDFWeight(pExtractor()->fullFn(partonBinInstances(), lastScale(),
evalPDFS));
if ( lastPDFWeight() == 0.0 ) {
lastCrossSection(ZERO);
return ZERO;
}
matrixElement()->setKinematics();
CrossSection xsec = matrixElement()->dSigHatDR() * lastPDFWeight();
if ( xsec == ZERO ) {
lastCrossSection(ZERO);
return ZERO;
}
xsec *= cutWeight();
lastAlphaS (matrixElement()->orderInAlphaS () >0 ?
matrixElement()->alphaS() : -1.);
lastAlphaEM(matrixElement()->orderInAlphaEW() >0 ?
matrixElement()->alphaEM() : -1.);
matrixElement()->fillProjectors();
if ( !projectors().empty() ) {
lastProjector(projectors().select(UseRandom::rnd()));
}
subProcess(SubProPtr());
if ( CKKWHandler() && matrixElement()->maxMultCKKW() > 0 &&
matrixElement()->maxMultCKKW() > matrixElement()->minMultCKKW() ) {
newSubProcess();
CKKWHandler()->setXComb(this);
xsec *= CKKWHandler()->reweightCKKW(matrixElement()->minMultCKKW(),
matrixElement()->maxMultCKKW());
}
if ( matrixElement()->reweighted() ) {
newSubProcess();
xsec *= matrixElement()->reWeight() * matrixElement()->preWeight();
}
lastCrossSection(xsec);
return xsec;
}
void StandardXComb::newSubProcess(bool group) {
if ( subProcess() ) return;
if ( head() && matrixElement()->wantCMS() ) {
// first get the meMomenta in their CMS, as this may
// not be the case
Boost cm = (meMomenta()[0] + meMomenta()[1]).findBoostToCM();
if ( cm.mag2() > Constants::epsilon ) {
for ( vector<Lorentz5Momentum>::iterator m = meMomenta().begin();
m != meMomenta().end(); ++m ) {
*m = m->boost(cm);
}
}
}
if ( !lastProjector() ) {
if ( !group )
subProcess(new_ptr(SubProcess(lastPartons(), tCollPtr(), matrixElement())));
else
subProcess(new_ptr(SubProcessGroup(lastPartons(), tCollPtr(), matrixElement())));
lastDiagramIndex(matrixElement()->diagram(diagrams()));
const ColourLines & cl = matrixElement()->selectColourGeometry(lastDiagram());
Lorentz5Momentum p1 = lastPartons().first->momentum();
Lorentz5Momentum p2 = lastPartons().second->momentum();
tPPair inc = lastPartons();
if ( mirror() ) swap(inc.first, inc.second);
if ( matrixElement()->wantCMS() &&
!matrixElement()->haveX1X2() ) {
LorentzRotation r = Utilities::boostToCM(inc);
lastDiagram()->construct(subProcess(), *this, cl);
subProcess()->transform(r.inverse());
lastPartons().first->set5Momentum(p1);
lastPartons().second->set5Momentum(p2);
} else {
lastDiagram()->construct(subProcess(), *this, cl);
}
lastPartons().first ->scale(partonBinInstances().first ->scale());
lastPartons().second->scale(partonBinInstances().second->scale());
for ( int i = 0, N = subProcess()->outgoing().size(); i < N; ++i )
subProcess()->outgoing()[i]->scale(lastScale());
// construct the spin information for the interaction
matrixElement()->constructVertex(subProcess(),&cl);
// set veto scales
matrixElement()->setVetoScales(subProcess());
} else {
lastProjector()->newSubProcess();
subProcess(lastProjector()->subProcess());
lastPartons(lastProjector()->lastPartons());
lastSHat((lastPartons().first->momentum() +
lastPartons().second->momentum()).m2());
lastX1X2(make_pair(lastPartons().first->momentum().plus()/
lastParticles().first->momentum().plus(),
lastPartons().second->momentum().minus()/
lastParticles().second->momentum().minus()));
lastY(log(lastX1()/lastX2())*0.5);
partonBinInstances().first->parton(lastPartons().first);
partonBinInstances().second->parton(lastPartons().second);
if ( !matrixElement()->keepRandomNumbers() )
throw Exception() << "Matrix element needs to request random number storage "
<< "for creating projected subprocesses."
<< Exception::runerror;
const double * r = &lastRandomNumbers()[0];
pExtractor()->generateSHat(lastS(), partonBinInstances(),
r, r + nDim() - partonDims.second,true);
pExtractor()->updatePartonBinInstances(partonBinInstances());
}
}
+void StandardXComb::checkReshufflingNeeds() {
+
+ theNeedsReshuffling = false;
+
+ for ( cPDVector::const_iterator p = mePartonData().begin() + 2;
+ p != mePartonData().end(); ++p ) {
+ theNeedsReshuffling |= ( (**p).mass() != (**p).hardProcessMass() );
+ }
+
+}
+
+double StandardXComb::reshuffleEquation(double k,
+ const vector<pair<Energy2,Energy2> >& coeffs,
+ Energy2 Q2) const {
+ double res = 0.;
+ for ( vector<pair<Energy2,Energy2> >::const_iterator c = coeffs.begin();
+ c != coeffs.end(); ++c ) {
+ double x = sqr(k)*c->first/Q2+c->second/Q2;
+ if ( x < 0.0 )
+ throw Veto();
+ res += sqrt(x);
+ }
+ return res - 1.;
+}
+
+double StandardXComb::solveReshuffleEquation(const vector<pair<Energy2,Energy2> >& coeffs,
+ Energy2 Q2) const {
+ // shamelessly adapted from Herwig++'s QTildeReconstructor
+ // in principle we could have included the exact solution for two
+ // outgoing particles, but for this to be simple we would require
+ // boosting to the CM so the gain is questionable
+ double k1 = 0.,k2 = 1.,k = 0.;
+ if ( reshuffleEquation(k1,coeffs,Q2) < 0.0 ) {
+ while ( reshuffleEquation(k2, coeffs, Q2) < 0.0 ) {
+ k1 = k2;
+ k2 *= 2;
+ }
+ while ( abs( (k1 - k2)/(k1 + k2) ) > 1.e-10 ) {
+ if( reshuffleEquation(k2,coeffs,Q2) == 0.0 ) {
+ return k2;
+ } else {
+ k = (k1+k2)/2.;
+ if ( reshuffleEquation(k,coeffs,Q2) > 0.0 ) {
+ k2 = k;
+ } else {
+ k1 = k;
+ }
+ }
+ }
+ return k1;
+ }
+ throw Veto();
+ return 1.;
+}
+
+void StandardXComb::reshuffle(vector<Lorentz5Momentum>& pout) const {
+
+ Lorentz5Momentum Q(ZERO,ZERO,ZERO,ZERO);
+ for ( vector<Lorentz5Momentum>::const_iterator p = pout.begin();
+ p != pout.end(); ++p )
+ Q += *p;
+
+ if ( Q.m2() <= ZERO )
+ throw Veto();
+
+ LorentzVector<double> Qhat = Q/sqrt(Q.m2());
+
+ vector<Energy> pQhat;
+ vector<pair<Energy2,Energy2> > coeffs;
+
+ vector<Lorentz5Momentum>::const_iterator p = pout.begin();
+ cPDVector::const_iterator pd = mePartonData().begin() + 2;
+
+ for ( ; p != pout.end(); ++p, ++pd ) {
+ pQhat.push_back((*p)*Qhat);
+ coeffs.push_back(make_pair(sqr(pQhat.back())-sqr((**pd).hardProcessMass()),
+ sqr((**pd).mass())));
+ }
+
+ double k = solveReshuffleEquation(coeffs,Q.m2());
+
+ vector<Lorentz5Momentum>::iterator px = pout.begin();
+ vector<Energy>::const_iterator pQ = pQhat.begin();
+ vector<pair<Energy2,Energy2> >::const_iterator coeff = coeffs.begin();
+
+ for ( ; px != pout.end(); ++px, ++pQ, ++coeff ) {
+ *px = k*((*px) - (*pQ) * Qhat) + sqrt(sqr(k)*coeff->first + coeff->second)*Qhat;
+ px->rescaleMass();
+ }
+
+}
+
tSubProPtr StandardXComb::construct() {
matrixElement()->setXComb(this);
if ( !head() ) {
if ( !cuts()->initSubProcess(lastSHat(), lastY()) ) throw Veto();
} else {
cuts()->initSubProcess(lastSHat(), lastY());
}
if ( head() )
if ( !matrixElement()->apply() )
return tSubProPtr();
setPartonBinInfo();
matrixElement()->setKinematics();
newSubProcess();
TmpTransform<tSubProPtr>
tmp(subProcess(), Utilities::getBoostToCM(subProcess()->incoming()));
if ( !cuts()->passCuts(*subProcess()) ) throw Veto();
if ( head() ) {
subProcess()->head(head()->subProcess());
if(lastCrossSection()!=ZERO&&head()->lastCrossSection()!=ZERO)
subProcess()->groupWeight(lastCrossSection()/head()->lastCrossSection());
else
subProcess()->groupWeight(0.);
}
return subProcess();
}
void StandardXComb::Init() {}
void StandardXComb::persistentOutput(PersistentOStream & os) const {
os << theSubProcessHandler << theME << theStats
<< theDiagrams << isMirror << theNDim << partonDims
<< theLastDiagramIndex << theMEInfo << theLastRandomNumbers << theMEPartonData
<< theLastPDFWeight << ounit(theLastCrossSection,nanobarn) << theLastJacobian
<< theLastME2 << theLastPreweight
<< ounit(theLastMECrossSection,nanobarn) << theLastMEPDFWeight << theLastMECouplings
<< theHead << theProjectors << theProjector << theKinematicsGenerated
- << checkedCuts << passedCuts << theCutWeight;
+ << checkedCuts << passedCuts << theCutWeight << theNeedsReshuffling;
}
void StandardXComb::persistentInput(PersistentIStream & is, int) {
is >> theSubProcessHandler >> theME >> theStats
>> theDiagrams >> isMirror >> theNDim >> partonDims
>> theLastDiagramIndex >> theMEInfo >> theLastRandomNumbers >> theMEPartonData
>> theLastPDFWeight >> iunit(theLastCrossSection,nanobarn) >> theLastJacobian
>> theLastME2 >> theLastPreweight
>> iunit(theLastMECrossSection,nanobarn) >> theLastMEPDFWeight >> theLastMECouplings
>> theHead >> theProjectors >> theProjector >> theKinematicsGenerated
- >> checkedCuts >> passedCuts >> theCutWeight;
+ >> checkedCuts >> passedCuts >> theCutWeight >> theNeedsReshuffling;
}
ClassDescription<StandardXComb> StandardXComb::initStandardXComb;
diff --git a/Handlers/StandardXComb.h b/Handlers/StandardXComb.h
--- a/Handlers/StandardXComb.h
+++ b/Handlers/StandardXComb.h
@@ -1,710 +1,744 @@
// -*- C++ -*-
//
// StandardXComb.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2011 Leif Lonnblad
// Copyright (C) 2009-2011 Simon Platzer
//
// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_StandardXComb_H
#define ThePEG_StandardXComb_H
// This is the declaration of the StandardXComb class.
#include "ThePEG/Config/ThePEG.h"
#include "SubProcessHandler.fh"
#include "ThePEG/PDF/PartonExtractor.fh"
#include "ThePEG/PDF/PartonBin.h"
#include "ThePEG/PDF/PartonBinInstance.h"
#include "ThePEG/Utilities/VSelector.h"
#include "ThePEG/Utilities/ClassDescription.h"
#include "ThePEG/Utilities/Maths.h"
#include "ThePEG/Utilities/XSecStat.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/MatrixElement/MEBase.h"
#include "ThePEG/Handlers/XComb.h"
#include "ThePEG/Handlers/StandardEventHandler.h"
#include "ThePEG/Handlers/SubProcessHandler.fh"
#include "StandardXComb.fh"
namespace ThePEG {
/**
* The StandardXComb class inherits from the more general XComb class
* which stores all information about the generation of a hard
* sub-proces for a given pair of incoming particles, a pair of
* extracted partons, etc. This class stores more information related
* to thestandard process generation scheme in ThePEG, such as the
* PartonExtractor and MEBase object used. It also does some of the
* administration of the process generation.
*
* The main function is dSigDR() which returns the differential cross
* section w.r.t. a given vector of random numbers in the interval
* ]0,1[. In the initialization this is used to pre-sample the phase
* space. In the generation phase it is used to give the cross section
* for a phase space point, and if this StandardXComb is chosen the
* construct() function is called to generate the actual sub-process.
*
* @see ParonExtractor
* @see MEBase
* @see Cuts
* @see StdXCombGroup
*/
class StandardXComb: public XComb {
public:
/** A vector of DiagramBase objects. */
typedef MEBase::DiagramVector DiagramVector;
/** A vector of indices. */
typedef MEBase::DiagramIndex DiagramIndex;
/** MEBase needs to be a friend. */
friend class MEBase;
public:
/** @name Standard constructors and destructors. */
//@{
/**
* Standard constructor.
*/
StandardXComb(Energy newMaxEnergy, const cPDPair & inc,
tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler,
tPExtrPtr newExtractor, tCascHdlPtr newCKKW,
const PBPair & newPartonBins, tCutsPtr newCuts, tMEPtr newME,
const DiagramVector & newDiagrams, bool mir,
tStdXCombPtr newHead = tStdXCombPtr());
/**
* Constructor given a head xcomb.
*/
StandardXComb(tStdXCombPtr newHead,
const PBPair & newPartonBins, tMEPtr newME,
const DiagramVector & newDiagrams);
/**
* Default constructor.
*/
StandardXComb();
/**
* Destructor.
*/
virtual ~StandardXComb();
/**
* Constructor used by MEBase to create a temporary object to store info.
*/
StandardXComb(tMEPtr me, const tPVector & parts, DiagramIndex i);
//@}
/** @name Utilities for incoming partons. */
//@{
/**
* Properly setup the PartonBinInstance objects provided a sub
* process has been constructed using this XComb.
*/
void recreatePartonBinInstances(Energy2 scale);
/**
* Fill the variables needed to generate remnants; momenta will be
* used from the partons set in this xcomb, but random numbers need
* to be provided to (re)generate variables not fixed by the
* incoming partons.
*/
void refillPartonBinInstances(const double* r);
/**
* Setup information on incoming partons depending
* on the information previously supplied through the
* choice of diagram and incoming momenta in the first
* two entries of meMomenta(). Partons are not actually
* extracted from the incoming particles, though a subprocess
* detached from the current Event may be created.
*/
void setIncomingPartons(tStdXCombPtr labHead = tStdXCombPtr());
/**
* Fill phase space information as far as possible
*/
void fill(const PPair& newParticles,
const PPair& newPartons,
const vector<Lorentz5Momentum>& newMEMomenta,
const DVector& newLastRandomNumbers = DVector());
//@}
/** @name Access the assigned objects used in the generation. */
//@{
/**
* Return a pointer to the corresponding sub-process handler. May be
* null if the standard process generation in ThePEG was not used.
*/
tcSubHdlPtr subProcessHandler() const { return theSubProcessHandler; }
/**
* The matrix element to be used.
*/
tMEPtr matrixElement() const { return theME; }
/**
* Return a pointer to the head XComb this XComb
* depends on. May return NULL, if this is not a
* member of a XComb group.
*/
tStdXCombPtr head() const { return theHead; }
/**
* Set the head XComb pointer.
*/
void head(tStdXCombPtr headXC) { theHead = headXC; }
/**
* Return a selector object of xcombs to choose subprocesses
* different than the one currently integrated.
*/
Selector<tStdXCombPtr>& projectors() { return theProjectors; }
/**
* Return a selector object of xcombs to choose subprocesses
* different than the one currently integrated.
*/
const Selector<tStdXCombPtr>& projectors() const { return theProjectors; }
/**
* Return a pointer to a projector xcomb which will generate a subprocess
* different from the one just integrated.
*/
tStdXCombPtr lastProjector() const { return theProjector; }
/**
* Set a pointer to a projector xcomb which will generate a subprocess
* different from the one just integrated.
*/
void lastProjector(tStdXCombPtr pxc) { theProjector = pxc; }
//@}
/** @name Main functions used for the generation. */
//@{
/**
* Try to determine if this subprocess is at all possible.
*/
virtual bool checkInit();
/**
* The number of dimensions of the phase space used to generate this
* process.
*/
virtual int nDim() const { return theNDim; }
/**
* Return the parton extraction dimensions
*/
const pair<int,int>& partonDimensions() const { return partonDims; }
/**
* Return true, if the current configuration will pass the cuts
*/
bool willPassCuts();
/**
* Return the cut weight encountered from fuzzy cuts
*/
double cutWeight() const { return theCutWeight; }
/**
* Reset all saved data about last generated phasespace point;
*/
virtual void clean();
/**
* Return true, if kinematics have already been generated
*/
bool kinematicsGenerated() const { return theKinematicsGenerated; }
/**
* Indicate that kinematics have been generated
*/
void didGenerateKinematics() { theKinematicsGenerated = true; }
/**
* Generate a phase space point from a vector \a r of \a nr numbers
* in the interval ]0,1[ and return the corresponding differential
* cross section.
*/
virtual CrossSection dSigDR(const pair<double,double> ll, int nr, const double * r);
/**
* If this XComb has a head XComb, return the cross section
* differential in the variables previously supplied. The PDF weight
* is taken from the lastPDFWeight supplied by the head XComb
* object.
*/
CrossSection dSigDR(const double * r);
/**
* Return the PDF weight used in the last call to dSigDR
*/
double lastPDFWeight() const { return theLastPDFWeight; }
/**
* Return the cross section calculated in the last call to dSigDR
*/
CrossSection lastCrossSection() const { return theLastCrossSection; }
/**
+ * Check if a reshuffling is required when constructing the hard
+ * subprocess.
+ */
+ void checkReshufflingNeeds();
+
+ /**
+ * Return true if a reshuffling is required when constructing the hard
+ * subprocess.
+ */
+ bool needsReshuffling() const { return theNeedsReshuffling; }
+
+ /**
+ * Perform the reshuffling from hardProcessMass to mass values,
+ * given outgoing momenta
+ */
+ void reshuffle(vector<Lorentz5Momentum>&) const;
+
+ /**
* Construct a sub-process object from the information available.
*/
virtual tSubProPtr construct();
//@}
/** @name Functions used for collecting statistics. */
//@{
/**
* The statistics object for this XComb.
*/
virtual const XSecStat & stats() const { return theStats; }
/**
* Select the current event. It will later be rejected with a
* probability given by \a weight.
*/
virtual void select(double weight) { theStats.select(weight); }
/**
* Accept the current event assuming it was previously selcted.
*/
virtual void accept() { theStats.accept(); }
/**
* Reweight a selected and accepted event.
*/
void reweight(double oldWeight, double newWeight) {
theStats.reweight(oldWeight,newWeight);
}
/**
* Reject the current event assuming it was previously accepted. If
* weighted events are produced, the \a weight should be the same as
* the previous call to select(double).
*/
virtual void reject(double weight = 1.0) { theStats.reject(weight); }
/**
* Reset statistics.
*/
virtual void reset() { theStats.reset(); }
//@}
/** @name Access information used by the MEBase object. */
//@{
/**
* The diagrams used by the matrix element.
*/
const DiagramVector & diagrams() const { return theDiagrams; }
/**
* True if the TreeDiagram's for this matrix element should in fact
* be mirrored before used to create an actual sub-rocess.
*/
bool mirror() const { return isMirror; }
/**
* Return the momenta of the partons to be used by the matrix
* element object, in the order specified by the TreeDiagram objects
* given by the matrix element.
*/
const vector<Lorentz5Momentum> & meMomenta() const { return theMEMomenta; }
/**
* Return the last selected diagram.
*/
tcDiagPtr lastDiagram() const { return diagrams()[lastDiagramIndex()]; }
/**
* Return the parton types to be used by the matrix element object,
* in the order specified by the TreeDiagram objects given by the
* matrix element.
*/
const cPDVector & mePartonData() const { return theMEPartonData; }
/**
* Return the index of the last selected diagram.
*/
DiagramIndex lastDiagramIndex() const { return theLastDiagramIndex; }
/**
* Get information saved by the matrix element in the calculation of
* the cross section to be used later when selecting diagrams and
* colour flow.
*/
const DVector & meInfo() const { return theMEInfo; }
/**
* Set information saved by the matrix element in the calculation of
* the cross section to be used later when selecting diagrams and
* colour flow.
*/
void meInfo(const DVector & info) { theMEInfo = info; }
/**
* Return the random numbers used to generate the
* last phase space point, if the matrix element
* requested so.
*/
const DVector& lastRandomNumbers() const { return theLastRandomNumbers; }
/**
* Get the last jacobian obtained when generating the kinematics
* for the call to dSigHatDR.
*/
double jacobian() const { return theLastJacobian; }
/**
* Return the matrix element squared as calculated
* for the last phase space point. This may optionally
* be used by a matrix element for caching.
*/
double lastME2() const { return theLastME2; }
/**
* Return the last preweight factor
*/
double lastPreweight() const { return theLastPreweight; }
/**
* Return the partonic cross section as calculated
* for the last phase space point. This may optionally
* be used by a matrix element for caching.
*/
CrossSection lastMECrossSection() const { return theLastMECrossSection; }
/**
* Return the PDF weight as calculated
* for the last phase space point, if the matrix
* element does supply PDF weights. This may optionally
* be used by a matrix element for caching.
*/
double lastMEPDFWeight() const { return theLastMEPDFWeight; }
/**
* Return the coupling factor as calculated for the lats phase space
* point.
*/
double lastMECouplings() const { return theLastMECouplings; }
//@}
/**
* Construct the corresponding SubProcess object if it hasn't been
* done before.
*/
virtual void newSubProcess(bool group = false);
/**
* Return the momenta of the partons to be used by the matrix
* element object, in the order specified by the TreeDiagram objects
* given by the matrix element.
*/
vector<Lorentz5Momentum> & meMomenta() { return theMEMomenta; }
/**
* Access the random numbers used to generate the
* last phase space point, if the matrix element
* requested so.
*/
DVector& lastRandomNumbers() { return theLastRandomNumbers; }
/**
* Return the parton types to be used by the matrix element object,
* in the order specified by the TreeDiagram objects given by the
* matrix element.
*/
cPDVector & mePartonData() { return theMEPartonData; }
/**
* Set the last selected diagram.
*/
void lastDiagramIndex(DiagramIndex i) { theLastDiagramIndex = i; }
/**
* Set the PDF weight used in the last call to dSigDR
*/
void lastPDFWeight(double w) { theLastPDFWeight = w; }
/**
* Set the cross section calculated in the last call to dSigDR
*/
void lastCrossSection(CrossSection s) { theLastCrossSection = s; }
/**
* Set the last jacobian obtained when generating the kinematics for
* the call to dSigHatDR.
*/
void jacobian(double j) { theLastJacobian = j; }
/**
* Set the matrix element squared as calculated
* for the last phase space point. This may optionally
* be used by a matrix element for caching.
*/
void lastME2(double v) { theLastME2 = v; }
/**
* Set the last preweight factor
*/
void lastPreweight(double w) { theLastPreweight = w; }
/**
* Set the partonic cross section as calculated
* for the last phase space point. This may optionally
* be used by a matrix element for caching.
*/
void lastMECrossSection(CrossSection v) { theLastMECrossSection = v; }
/**
* Set the PDF weight as calculated
* for the last phase space point, if the matrix
* element does supply PDF weights. This may optionally
* be used by a matrix element for caching.
*/
void lastMEPDFWeight(double v) { theLastMEPDFWeight = v; }
/**
* Set the coupling factor
*/
void lastMECouplings(double v) { theLastMECouplings = v; }
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* Standard Init function used to initialize the interface.
*/
static void Init();
private:
/**
* The corresponding sub-process handler
*/
tSubHdlPtr theSubProcessHandler;
/**
* The matrix element to be used.
*/
tMEPtr theME;
/**
* Statistics gathering for this XComb.
*/
XSecStat theStats;
/**
* The diagrams used by the matrix element.
*/
DiagramVector theDiagrams;
/**
* True if the TreeDiagram's for this matrix element should in fact
* be mirrored before used to create an actual sub-rocess.
*/
bool isMirror;
/**
* The number of dimensions of the phase space used to generate this
* process.
*/
int theNDim;
protected:
/**
* The number of dimensions of the phase space used for each of the
* incoming partons.
*/
pair<int,int> partonDims;
private:
/**
* True, if kinematics have already been generated
*/
bool theKinematicsGenerated;
/**
* The momenta of the partons to be used by the matrix element
* object, in the order specified by the TreeDiagram objects given
* by the matrix element.
*/
vector<Lorentz5Momentum> theMEMomenta;
/**
* The parton types to be used by the matrix element object, in the
* order specified by the TreeDiagram objects given by the matrix
* element.
*/
cPDVector theMEPartonData;
/**
* The last selected tree diagram.
*/
DiagramIndex theLastDiagramIndex;
/**
* Information saved by the matrix element in the calculation of the
* cross section to be used later when selecting diagrams and colour
* flow.
*/
DVector theMEInfo;
/**
* The random numbers used to generate the
* last phase space point, if the matrix element
* requested so.
*/
DVector theLastRandomNumbers;
/**
* The PDF weight used in the last call to dSigDR
*/
double theLastPDFWeight;
/**
* The cross section calculated in the last call to dSigDR
*/
CrossSection theLastCrossSection;
/**
* Save the last jacobian obtained when generating the kinematics for
* the call to dSigHatDR.
*/
double theLastJacobian;
/**
* The matrix element squared as calculated
* for the last phase space point. This may optionally
* be used by a matrix element for caching.
*/
double theLastME2;
/**
* The last preweight factor
*/
double theLastPreweight;
/**
* The partonic cross section as calculated
* for the last phase space point. This may optionally
* be used by a matrix element for caching.
*/
CrossSection theLastMECrossSection;
/**
* The PDF weight as calculated
* for the last phase space point, if the matrix
* element does supply PDF weights. This may optionally
* be used by a matrix element for caching.
*/
double theLastMEPDFWeight;
/**
* The coupling factor
*/
double theLastMECouplings;
/**
* A pointer to the head XComb this XComb
* depends on. May return NULL, if this is not a
* member of a XComb group.
*/
tStdXCombPtr theHead;
/**
* A selector object of xcombs to choose subprocesses
* different than the one currently integrated.
*/
Selector<tStdXCombPtr> theProjectors;
/**
* A pointer to a projector xcomb which will generate a subprocess
* different from the one just integrated.
*/
tStdXCombPtr theProjector;
/**
* True, if cuts have already been checked
*/
bool checkedCuts;
/**
* The result of the last call to willPassCuts
*/
bool passedCuts;
/**
* The cut weight encountered from fuzzy cuts
*/
double theCutWeight;
+ /**
+ * True if a reshuffling is required when constructing the hard
+ * subprocess.
+ */
+ bool theNeedsReshuffling;
+
+ /**
+ * Calculate the reshuffling equation given the coefficients
+ */
+ double reshuffleEquation(double, const vector<pair<Energy2,Energy2> >&, Energy2) const;
+
+ /**
+ * Solve the reshuffling equation given the coefficients
+ */
+ double solveReshuffleEquation(const vector<pair<Energy2,Energy2> >&, Energy2) const;
+
private:
/**
* Describe a concrete class with persistent data.
*/
static ClassDescription<StandardXComb> initStandardXComb;
/**
* Private and non-existent assignment operator.
*/
StandardXComb & operator=(const StandardXComb &);
};
/** @cond TRAITSPECIALIZATIONS */
/**
* This template specialization informs ThePEG about the base class of
* StandardXComb.
*/
template <>
struct BaseClassTrait<StandardXComb,1>: public ClassTraitsType {
/** Typedef of the base class of StandardXComb. */
typedef XComb NthBase;
};
/**
* This template specialization informs ThePEG about the name of the
* StandardXComb class.
*/
template <>
struct ClassTraits<StandardXComb>:
public ClassTraitsBase<StandardXComb> {
/** Return the class name. */
static string className() { return "ThePEG::StandardXComb"; }
};
/** @endcond */
}
#endif /* ThePEG_StandardXComb_H */
diff --git a/MatrixElement/Tree2toNDiagram.cc b/MatrixElement/Tree2toNDiagram.cc
--- a/MatrixElement/Tree2toNDiagram.cc
+++ b/MatrixElement/Tree2toNDiagram.cc
@@ -1,440 +1,440 @@
// -*- C++ -*-
//
// Tree2toNDiagram.cc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2011 Leif Lonnblad
//
// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the Tree2toNDiagram class.
//
#include "Tree2toNDiagram.h"
#include "ThePEG/EventRecord/SubProcess.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Utilities/Rebinder.h"
#include "ThePEG/Utilities/UtilityBase.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include "ThePEG/PDT/EnumParticles.h"
using namespace ThePEG;
Tree2toNDiagram::~Tree2toNDiagram() {}
Tree2toNDiagram & Tree2toNDiagram::add(tcPDPtr pd) {
if ( thePartons.size() < theNSpace ) addSpacelike(pd);
else addTimelike(pd, nextOrig);
return *this;
}
void Tree2toNDiagram::addTimelike(tcPDPtr pd, size_type orig) {
if ( allPartons().size() < theNSpace ||
orig >= allPartons().size())
throw Tree2toNDiagramError();
thePartons.push_back(pd);
theParents.push_back(orig);
}
tPVector Tree2toNDiagram::
construct(SubProPtr sp, const StandardXComb & xc, const ColourLines & cl) const {
tPVector out;
vector<Lorentz5Momentum> pout(xc.meMomenta().begin() + 2,
xc.meMomenta().end());
-// Utilities::transform(pout.begin(), pout.end(),
-// Utilities::getBoostFromCM(xc.lastPartons()));
+ if ( xc.needsReshuffling() )
+ xc.reshuffle(pout);
tPPair in = xc.lastPartons();
if ( xc.mirror() ) swap(in.first, in.second);
tPVector ret;
if ( in.first->dataPtr() != allPartons()[0] ||
in.second->dataPtr() != allPartons()[nSpace() - 1] )
throw Tree2toNDiagramError();
PVector slike;
slike.push_back(in.first);
for ( int i = 1; i < nSpace() - 1; ++i )
slike.push_back(allPartons()[i]->produceParticle());
slike.push_back(in.second);
ret = tPVector(slike.begin(), slike.end());
for ( size_type i = 1; i < slike.size() - 1; ++i ) {
slike[i-1]->addChild(slike[i]);
sp->addIntermediate(slike[xc.mirror()? i: slike.size() - 1 - i], false);
}
int io = pout.size();
PVector tlike(allPartons().size() - nSpace());
ParticleSet done;
for ( int i = allPartons().size() - 1; i >= nSpace(); --i ) {
int it = i - nSpace();
pair<int,int> ch = children(i);
bool iso = ch.first < 0;
if ( iso ) {
tlike[it] = allPartons()[i]->produceParticle(pout[--io]);
done.insert(tlike[it]);
} else {
Lorentz5Momentum p = tlike[ch.first - nSpace()]->momentum() +
tlike[ch.second - nSpace()]->momentum();
tlike[it] = allPartons()[i]->produceParticle(p);
}
if ( parent(i) < nSpace() ) {
slike[parent(i)]->addChild(tlike[it]);
if ( parent(i) == nSpace() - 2 )
slike[parent(i) + 1]->addChild(tlike[it]);
}
if ( !iso ) {
tlike[it]->addChild(tlike[ch.first - nSpace()]);
tlike[it]->addChild(tlike[ch.second - nSpace()]);
}
if ( iso ) out.push_back(tlike[it]);
else sp->addIntermediate(tlike[it], false);
}
ret.insert(ret.end(), tlike.begin(), tlike.end());
for ( int i = 0, N = out.size(); i < N; ++i )
sp->addOutgoing(out[xc.mirror()? i: out.size() - i - 1], false);
for ( PVector::size_type i = 0; i < slike.size() - 2; ++i ) {
pair<int,int> ch = children(i);
slike[ch.first]->set5Momentum(slike[i]->momentum() -
tlike[ch.second - nSpace()]->momentum());
}
cl.connect(ret);
return out;
}
tcPDVector Tree2toNDiagram::outgoing() const {
tcPDVector pdv;
for ( size_type i = nSpace(); i < allPartons().size(); ++i )
if ( children(i).first < 0 ) pdv.push_back(allPartons()[i]);
return pdv;
}
tcPDVector Tree2toNDiagram::external() const {
tcPDVector pdv;
pdv.push_back(allPartons()[0]);
pdv.push_back(allPartons()[nSpace() - 1]);
for ( size_type i = nSpace(); i < allPartons().size(); ++i )
if ( children(i).first < 0 ) pdv.push_back(allPartons()[i]);
return pdv;
}
tcPDPair Tree2toNDiagram::incoming() const {
return tcPDPair(allPartons()[0], allPartons()[nSpace() - 1]);
}
pair<int,int> Tree2toNDiagram::children(int ii) const {
pair<int,int> ret = make_pair(-1, -1);
for ( size_type i = 0; i < theParents.size(); ++i ) {
if ( parent(i) == ii ) {
if ( ret.first < 0 ) ret.first = i;
else if ( ret.second < 0 ) ret.second = i;
else throw Tree2toNDiagramError();
}
}
return ret;
}
void Tree2toNDiagram::check() {
vector< pair<int,int> > children(allPartons().size(), make_pair(-1, -1));
theNOutgoing = 0;
for ( size_type i = nSpace(); i < allPartons().size(); ++i ) {
if ( children[parent(i)].first < 0 ) children[parent(i)].first = i;
else if ( children[parent(i)].second < 0 ) children[parent(i)].second = i;
else throw Tree2toNDiagramError();
}
for ( size_type i = nSpace(); i < allPartons().size(); ++i ) {
if ( children[i].first < 0 && children[i].second < 0 ) ++theNOutgoing;
else if ( children[i].first < 0 || children[i].second < 0 )
throw Tree2toNDiagramError();
}
cPDVector parts(2);
parts[0] = incoming().first;
parts[1] = incoming().second;
tcPDVector out(outgoing());
parts.insert(parts.end(), out.begin(), out.end());
partons(2, parts, nextOrig + 1);
}
bool Tree2toNDiagram::isSame (tcDiagPtr diag) const {
Ptr<Tree2toNDiagram>::tcptr cmp =
dynamic_ptr_cast<Ptr<Tree2toNDiagram>::tcptr>( diag );
if ( !cmp )
return false;
return equals(cmp) && external() == cmp->external();
}
bool Tree2toNDiagram::isSame (tcDiagPtr diag, map<int,int>& remap) const {
Ptr<Tree2toNDiagram>::tcptr cmp =
dynamic_ptr_cast<Ptr<Tree2toNDiagram>::tcptr>( diag );
if ( !cmp )
return false;
remap.clear();
remap[0] = 0;
return equals(cmp,remap);
}
bool Tree2toNDiagram::equals(Ptr<Tree2toNDiagram>::tcptr diag,
int start, int startCmp) const {
if ( start < 0 && startCmp < 0 )
return true;
if ( allPartons()[start] != diag->allPartons()[startCmp] )
return false;
pair<int,int> ch = children(start);
pair<int,int> chCmp = diag->children(startCmp);
bool match =
equals(diag,ch.first,chCmp.first) &&
equals(diag,ch.second,chCmp.second);
// also try swapped outgoing legs on same vertex
if ( !match && start > nSpace() - 1 &&
children(ch.first).first < 0 && children(ch.second).first < 0 &&
diag->children(chCmp.first).first < 0 && diag->children(chCmp.second).first < 0 )
match =
equals(diag,ch.first,chCmp.second) &&
equals(diag,ch.second,chCmp.first);
return match;
}
bool Tree2toNDiagram::equals(Ptr<Tree2toNDiagram>::tcptr diag,
map<int,int>& remap,
int start, int startCmp) const {
if ( start < 0 && startCmp < 0 )
return true;
if ( allPartons()[start] != diag->allPartons()[startCmp] )
return false;
pair<int,int> ch = children(start);
pair<int,int> chCmp = diag->children(startCmp);
if ( ch.first < 0 && chCmp.first < 0 ) {
remap[externalId(start)] = diag->externalId(startCmp);
}
bool match =
equals(diag,remap,ch.first,chCmp.first) &&
equals(diag,remap,ch.second,chCmp.second);
// also try swapped outgoing legs on same vertex
if ( !match && start > nSpace() - 1 &&
children(ch.first).first < 0 && children(ch.second).first < 0 &&
diag->children(chCmp.first).first < 0 && diag->children(chCmp.second).first < 0 )
match =
equals(diag,remap,ch.first,chCmp.second) &&
equals(diag,remap,ch.second,chCmp.first);
return match;
}
int Tree2toNDiagram::externalId(int id) const {
if ( id < 0 )
return -1;
if ( id == 0 )
return 0;
if ( id == nSpace() - 1 )
return 1;
int k = 1;
for ( size_type i = nSpace(); i < allPartons().size(); ++i ) {
if ( children(i).first < 0 ) ++k;
if ( i == size_type(id) )
break;
}
return k;
}
int Tree2toNDiagram::diagramId(int id) const {
if ( id < 0 )
return -1;
if ( id == 0 ) return 0;
if ( id == 1 ) return nSpace() - 1;
int k = 1;
size_type i = nSpace();
for ( ; i < allPartons().size(); ++i ) {
if ( children(i).first < 0 ) ++k;
if ( k == id )
break;
}
return i;
}
int Tree2toNDiagram::mergeEmission(int emitter, int id, map<int,int>& remap) {
if ( id < 2 )
return -1;
if ( remap.find(emitter) != remap.end() ) {
remap.erase(emitter);
}
if ( remap.find(id) != remap.end() ) {
remap.erase(id);
}
for ( map<int,int>::iterator rm = remap.begin();
rm != remap.end(); ++rm ) {
if ( rm->first == 0 || rm->first == 1 ) {
rm->second = rm->first;
} else {
rm->second = diagramId(rm->first);
}
}
// translate to diagram id
int did = diagramId(id);
int demitter = diagramId(emitter);
if ( children(did) != make_pair(-1,-1) )
return -1;
// now get the parent
int p = parent(did);
int npos = -1;
if ( p == 0 || p == nSpace() - 2 ) {
npos = ( p == 0 ? 0 : 1 );
} else if ( p >= nSpace() ) {
if ( id > emitter )
npos = emitter;
else
npos = emitter - 1;
}
pair<int,int> remove;
size_type theNSpaceBackup = theNSpace;
int theNOutgoingBackup = theNOutgoing;
int nextOrigBackup = nextOrig;
cPDVector thePartonsBackup = thePartons;
vector<int> theParentsBackup = theParents;
int deltaFlow = 0;
if ( npos == 1 ) {
if ( thePartons[did]->CC() )
deltaFlow -= ( thePartons[did]->id() < 0 ? -1 : 1 );
if ( thePartons[nSpace()-1]->CC() )
deltaFlow += ( thePartons[nSpace()-1]->id() < 0 ? -1 : 1 );
}
// emitted from spacelike
if ( p == 0 || p == nSpace() - 2 ) {
if ( p == 0 && p != demitter )
return -1;
if ( p == nSpace() - 2 && demitter != nSpace()-1 )
return -1;
if ( p == 0 )
remove = make_pair(p,did);
else
remove = make_pair(nSpace()-1,did);
--theNSpace;
--theNOutgoing;
} else if ( p >= nSpace() ) {
remove = children(p);
if ( remove.first != demitter )
swap(remove.first,remove.second);
if ( remove != make_pair(demitter,did) )
return -1;
--theNOutgoing;
} else {
return -1;
}
if ( remove.first > remove.second )
swap(remove.first,remove.second);
for ( map<int,int>::iterator rm = remap.begin();
rm != remap.end(); ++rm ) {
if ( rm->first > 1 ) {
if ( rm->second > remove.first &&
rm->second < remove.second )
rm->second -= 1;
else if ( rm->second > remove.second )
rm->second -= 2;
}
}
for ( unsigned int k = remove.first + 1; k < theParents.size(); ++k ) {
if ( theParents[k] >= remove.first &&
theParents[k] < remove.second &&
theParents[k] >= 0 )
theParents[k] -= 1;
else if ( theParents[k] > remove.second && theParents[k] > 0 )
theParents[k] -= 2;
}
thePartons.erase(thePartons.begin() + remove.second);
theParents.erase(theParents.begin() + remove.second);
thePartons.erase(thePartons.begin() + remove.first);
theParents.erase(theParents.begin() + remove.first);
if ( npos > 1 )
if ( npos != externalId(p) ) {
pair<int,int> swapDiagIds(p,diagramId(npos));
swap(thePartons[swapDiagIds.first],thePartons[swapDiagIds.second]);
swap(theParents[swapDiagIds.first],theParents[swapDiagIds.second]);
for ( map<int,int>::iterator rm = remap.begin();
rm != remap.end(); ++rm ) {
if ( rm->first > 1 ) {
if ( rm->second == swapDiagIds.first ) {
rm->second = swapDiagIds.second;
} else if ( rm->second == swapDiagIds.second ) {
rm->second = swapDiagIds.first;
}
}
}
}
for ( map<int,int>::iterator rm = remap.begin();
rm != remap.end(); ++rm ) {
if ( rm->first > 1 ) {
rm->second = externalId(rm->second);
}
}
if ( npos == 1 ) {
if ( thePartons[nSpace()-1]->CC() )
deltaFlow -= ( thePartons[nSpace()-1]->id() < 0 ? -1 : 1 );
if ( deltaFlow != 0 )
thePartons[nSpace()-1] = thePartons[nSpace()-1]->CC();
}
try {
check();
} catch (Tree2toNDiagramError&) {
theNSpace = theNSpaceBackup;
theNOutgoing = theNOutgoingBackup;
nextOrig = nextOrigBackup;
thePartons = thePartonsBackup;
theParents = theParentsBackup;
return -1;
}
return npos;
}
ClassDescription<Tree2toNDiagram> Tree2toNDiagram::initTree2toNDiagram;
void Tree2toNDiagram::persistentInput(PersistentIStream & is, int) {
is >> theNSpace >> theNOutgoing >> thePartons >> theParents >> nextOrig;
}
void Tree2toNDiagram::persistentOutput(PersistentOStream & os) const {
os << theNSpace << theNOutgoing << thePartons << theParents << nextOrig;
}
Tree2toNDiagramError::Tree2toNDiagramError() {
theMessage << "An error occurred while setting up a diagram of class "
<< "'Tree2toNDiagram'.";
severity(abortnow);
}
diff --git a/PDT/Matcher.cc b/PDT/Matcher.cc
--- a/PDT/Matcher.cc
+++ b/PDT/Matcher.cc
@@ -1,67 +1,69 @@
// -*- C++ -*-
//
// Matcher.cc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2011 Leif Lonnblad
//
// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the Matcher class.
//
#include "Matcher.h"
#include "StandardMatchers.h"
#ifdef ThePEG_TEMPLATES_IN_CC_FILE
#include "Matcher.tcc"
#endif
#define THEPEG_MATCH_DESC(T) \
/** \
* This template specialization registers the Matcher \
*/ \
template <> \
NoPIOClassDescription<T> T::initMatcher \
= NoPIOClassDescription<T>(); \
namespace ThePEG {
THEPEG_MATCH_DESC(MatchAny)
THEPEG_MATCH_DESC(MatchStandardQCDParton)
THEPEG_MATCH_DESC(MatchLightAntiQuark)
THEPEG_MATCH_DESC(MatchLightQuark)
THEPEG_MATCH_DESC(MatchLepton)
THEPEG_MATCH_DESC(MatchDiquark)
THEPEG_MATCH_DESC(MatchMeson)
THEPEG_MATCH_DESC(MatchBaryon)
THEPEG_MATCH_DESC(MatchNegative)
THEPEG_MATCH_DESC(MatchNeutral)
THEPEG_MATCH_DESC(MatchPositive)
THEPEG_MATCH_DESC(MatchCharged)
+ THEPEG_MATCH_DESC(MatchNeutrino)
}
using namespace ThePEG;
namespace {
void dummy() {
static MatchAny m00;
static MatchStandardQCDParton m01;
static MatchLightAntiQuark m02;
static MatchLightQuark m03;
static MatchLepton m04;
static MatchDiquark m05;
static MatchMeson m06;
static MatchBaryon m07;
static MatchNegative m08;
static MatchNeutral m09;
static MatchPositive m11;
static MatchCharged m12;
+ static MatchNeutrino m13;
}
}
diff --git a/PDT/ParticleData.cc b/PDT/ParticleData.cc
--- a/PDT/ParticleData.cc
+++ b/PDT/ParticleData.cc
@@ -1,954 +1,1034 @@
// -*- C++ -*-
//
// ParticleData.cc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2011 Leif Lonnblad
//
// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the 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::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) {
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) {
theMass = mi;
ParticleData * apd = CC().operator->();
if ( synchronized() && apd ) apd->theMass = theMass;
}
+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/ParticleData.h b/PDT/ParticleData.h
--- a/PDT/ParticleData.h
+++ b/PDT/ParticleData.h
@@ -1,901 +1,963 @@
// -*- C++ -*-
//
// ParticleData.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2011 Leif Lonnblad
//
// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_ParticleData_H
#define ThePEG_ParticleData_H
// This is the declaration of the ParticleData class.
#include "ThePEG/Config/ThePEG.h"
#include "ThePEG/PDT/PDT.h"
#include "ThePEG/PDT/PID.h"
#include "ThePEG/Vectors/LorentzVector.h"
#include "ThePEG/Vectors/ThreeVector.h"
#include "ThePEG/Interface/Interfaced.h"
#include "ThePEG/Utilities/Selector.h"
#include "ThePEG/PDT/WidthGenerator.h"
#include "ThePEG/PDT/MassGenerator.h"
#include "ThePEG/PDT/DecayMode.fh"
#include "ThePEG/Utilities/ClassTraits.h"
#include "ThePEG/Utilities/ClassDescription.h"
namespace ThePEG {
/**
* ParticleData inherits from InterfacedBase and represents the
* properties of a particle type. It is also able to produce instances
* of this Particle type and, among other things, to decay them.
*
* @see \ref ParticleDataInterfaces "The interfaces"
* defined for ParticleData.
*/
class ParticleData: public Interfaced {
public:
/** The Repository is a friend. */
friend class Repository;
/** The EventGenerator is a friend. */
friend class EventGenerator;
/** DecayMode is a friend. */
friend class DecayMode;
/** A selector of DecayMode objects. */
typedef Selector<tDMPtr> DecaySelector;
public:
/** @name Standard constructors and destructors. */
//@{
/**
* Default constructor.
*/
ParticleData();
/**
* Destructor.
*/
virtual ~ParticleData();
//@}
/** @name The Create methods are special interfaces for ParticleData
classes. */
//@{
/**
* Create a Particle which is its own anti-particle.
*/
static PDPtr Create(PID newId, const string & newPDGName);
/**
* Create a particle - anti particle pair.
*/
static PDPair Create(PID newId, const string & newPDGName, const string & newAntiPDGName);
//@}
public:
/** @name Acces to number and name. */
//@{
/**
* Return the PDG id number.
*/
long id() const { return theId; }
/**
* Return the generic PDG name. Note that this is not really
* standardised.
*/
const string & PDGName() const { return thePDGName; }
/**
* Return the generic PDG name. Note that this is not really
* standardised.
*/
const string & genericName() const { return thePDGName; }
//@}
/** @name Functions used for producing Particle instances. */
//@{
/**
* Produce a particle specifying momentum.
*/
PPtr produceParticle(const Lorentz5Momentum &) const;
/**
* Produce a particle specifying momentum.
*/
PPtr produceParticle(const LorentzMomentum &) const;
/**
* Produce a particle specifying 4-momentum and a mass.
*/
PPtr produceParticle(const LorentzMomentum &, Energy m) const;
/**
* Produce a particle specifying 3-momentum.
*/
PPtr produceParticle(const Momentum3 & pp = Momentum3()) const;
/**
* Produce a particle specifying mass and 3-momentum.
*/
PPtr produceParticle(Energy m, const Momentum3 & pp = Momentum3()) const;
/**
* Produce a particle specifying light-cone momentum components and
* transverse momentum components.
*/
PPtr produceParticle(Energy plus, Energy minus, Energy px, Energy py) const;
/**
* Generate a mass for an instance of this particle type.
*/
Energy generateMass() const;
/**
* Generate a width for an instance of this particle type. Given a
* \a mass of an instance of this particle type, calculate its width.
*/
Energy generateWidth(Energy mass) const;
/**
* Generate a mass for an instance of this particle type. Given a \a
* mass and a \a width of an instance of this particle type,
* generate a life time.
*/
Length generateLifeTime(Energy mass, Energy width) const;
// Given a mass and a width of an instance of this particle type,
// generate a life time.
//@}
/** @name Access the decay modes. */
//@{
/**
* Return the nominal decay selector for this particle. Ie. the
* decay modes weighted by their nominal branching ratios.
*/
const DecaySelector & decaySelector() const { return theDecaySelector; }
/**
* Selects a decay mode randomly according to the branching
* ratios. The nominal branching ratios may be changed for the
* particular Particle instance \a p, iether by an assigned
* WidthGenerator or the respective Decayers.
*/
tDMPtr selectMode(Particle & p) const;
/**
* Access all the decay modes, including those which are
* switched off, or have zero branching ratio
*/
const DecaySet & decayModes() const { return theDecayModes; }
//@}
/**
* Return the nominal mass.
*/
Energy mass() const { return theMass; }
/**
+ * Return the mass to be used when evaluating hard process cross sections.
+ */
+ Energy hardProcessMass() const {
+ return
+ hardProcessMassSet ? theHardProcessMass : mass();
+ }
+
+ /**
* Return the maximum possible mass of this particle type.
*/
Energy massMax() const { return mass() + widthUpCut(); }
/**
* Return the minimum possible mass of this particle type.
*/
Energy massMin() const { return max(mass() - widthLoCut(), ZERO); }
/**
* Return the constituent mass of this particle if relevant. This
* version simply returns the nominal mass.
*/
virtual Energy constituentMass() const { return mass(); }
/**
* Set the width.
*/
Energy width(Energy);
/**
* Get the width. If no width is specified, it is calculated from
* the lifetime.
*/
Energy width() const {
return theWidth >= ZERO ? theWidth :
( theCTau > Length() ? hbarc/theCTau :
( theCTau == Length() ? Constants::MaxEnergy : ZERO ) );
}
/**
* Set the width cut. Both upper and lower cut is set.
*/
Energy widthCut(Energy wci) {
widthUpCut(wci);
return widthLoCut(wci);
}
/**
* Get the width cut.
*/
Energy widthCut() const { return max(widthUpCut(), widthLoCut()); }
/**
* Set the upper width cut.
*/
Energy widthUpCut(Energy);
/**
* Get the upper width cut.
*/
Energy widthUpCut() const {
return theWidthUpCut >= ZERO? theWidthUpCut: Constants::MaxEnergy;
}
/**
* Set the lower width cut.
*/
Energy widthLoCut(Energy);
/**
* Get the lower width cut.
*/
Energy widthLoCut() const {
return theWidthLoCut >= ZERO? theWidthLoCut: Constants::MaxEnergy;
}
/**
* Set the life time cTau.
*/
Length cTau(Length);
/**
* Get the life time cTau cTau. If no life time is specified, it is
* calculated from the width. If the width is also not specified,
* the lifetime is assumed to be zero for ustable particles and
* infinite for stable ones.
*/
Length cTau() const {
return theCTau > Length() ? theCTau :
( theWidth > ZERO ? hbarc/theWidth :
( stable() ? Constants::MaxLength : Length() ) );
}
/**
+ * Return the width to be used when evaluating hard process cross sections.
+ */
+ Energy hardProcessWidth() const {
+ return
+ hardProcessWidthSet ? theHardProcessWidth : width();
+ }
+
+ /**
* Set the charge. The charge should be given
* in units of e/3 using the PDT::Charge enum.
*/
PDT::Charge iCharge(PDT::Charge);
/**
* Get the charge. The charge is returned in standard units and in
* iCharge the charge is returned in units of e/3.
*/
Charge charge() const { return eplus*double(theCharge)/3.0; }
/**
* Get the charge. The charge is returned in units of e/3.
*/
PDT::Charge iCharge() const { return theCharge; }
/**
* Return true if charged.
*/
bool charged() const { return PDT::charged(theCharge); }
/**
* Return true if positively charged.
*/
bool positive() const { return PDT::positive(theCharge); }
/**
* Return true if negatively charged.
*/
bool negative() const { return PDT::negative(theCharge); }
/**
* Set the spin. The spin should be given as 2J+1 (in units of
* hbar/2) using the PDT::Spin enum.
*/
PDT::Spin iSpin(PDT::Spin);
/**
* Get the spin.The spin is returned in standard units.
*/
AngularMomentum spin() const { return hbar_Planck*double(theSpin-1)*0.5; }
/**
* Get the spin. The spin is returned as 2J+1 in units of hbar/2.
*/
PDT::Spin iSpin() const { return theSpin; }
/**
* Set the colour of the particle in units of PDT::Colour.
*/
PDT::Colour iColour(PDT::Colour);
/**
* Get the colour of the particle in units of PDT::Colour.
*/
PDT::Colour iColour() const { return theColour; }
/**
* Return true if coloured.
*/
bool coloured() const { return PDT::coloured(iColour()); }
/**
* Return true if (\a anti) coloured or colour-octet.
*/
bool hasColour(bool anti = false) const {
return anti? hasAntiColour():
( iColour() == PDT::Colour3 || iColour() == PDT::Colour6 ||
iColour() == PDT::Colour8 );
}
/**
* Return true if anti coloured or colour-octet.
*/
bool hasAntiColour() const {
return iColour() == PDT::Colour3bar || iColour() == PDT::Colour6bar ||
iColour() == PDT::Colour8;
}
/**
* Specify if particle is to be considered stable according to \a
* stab.
*/
void stable(bool stab);
/**
* Return true if particle is to be considered stable. If the decay
* table is empty the function always returns true, even if the
* member variable is false.
*/
bool stable() const { return isStable; }
/**
* Get the pointer to the corresponding anti partner.
*/
tPDPtr CC() const { return theAntiPartner; }
/**
* Specify if the anti partner chould be changed automatically when
* this object is changed according to \a sync.
*/
void synchronized(bool sync);
/**
* Return true if the anti partner chould be changed automatically
* when this object is changed.
*/
bool synchronized() const { return syncAnti; }
/**
* If there is an anti-partner, update this object to have correct
* anti-properties.
*/
void synchronize();
/**
* Set the mass generator object.
*/
void massGenerator(tMassGenPtr);
/**
* Get the mass generator object.
*/
tMassGenPtr massGenerator() const { return theMassGenerator; }
/**
* Set the width generator object.
*/
void widthGenerator(tWidthGeneratorPtr);
/**
* Get the width generator object.
*/
tWidthGeneratorPtr widthGenerator() const { return theWidthGenerator; }
/**
* Specify if the branching ratio of the Particle instances should vary with their
* masses.
*/
void variableRatio(bool varRatio);
/**
* Return true if the branching ratio should vary with the mass of the Particle
* instance.
*/
bool variableRatio() const { return theVariableRatio; }
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
/**
* Special clone function used by the Repository. Also makes copies
* the decay modes and the anti-partner if it exists and if
* synchronized() is true.
*/
virtual PDPtr pdclone() const;
/**
* Protected constructor only to be used by subclasses or by the
* Create method.
*/
ParticleData(PID newId, const string & newPDGName);
/**
* Read setup info from a standard stream. The following information
* must be supplied in a white-space separated list: PDG number,
* generic name, default mass (GeV), default width (GeV), width cut
* (GeV), the lifetime ctau (mm), the charge, the colour, the spin,
* stable (true) or not (false). Note that if a minus sign is given
* instead of a generic name, the name of the object will be used
* instead.
*/
virtual void readSetup(istream & is);
/**
* Used by subclasses or by the Create method to setup
* anti-relationship.
*/
static void antiSetup(const PDPair & pap);
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Check sanity of the object during the setup phase.
*/
virtual void doupdate();
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
/**
* Rebind pointer to other Interfaced objects. Called in the setup phase
* after all objects used in an EventGenerator has been cloned so that
* the pointers will refer to the cloned objects afterwards.
* @param trans a TranslationMap relating the original objects to
* their respective clones.
* @throws RebindException if no cloned object was found for a given
* pointer.
*/
virtual void rebind(const TranslationMap & trans)
;
/**
* Return a vector of all pointers to Interfaced objects used in this
* object.
* @return a vector of pointers.
*/
virtual IVector getReferences();
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
//@}
protected:
/**
* Add a decay mode for this particle.
*/
void addDecayMode(tDMPtr);
/**
* Remove a decay mode for this particle.
*/
void removeDecayMode(tDMPtr);
private:
/**
* Id number according to the STDHEP/PDG standard.
*/
PID theId;
/**
* Name and Id number according to the STDHEP/PDG standard.
*/
string thePDGName;
/**
* Nominal mass.
*/
Energy theMass;
/**
* Width.
*/
Energy theWidth;
/**
+ * The mass to be used when evaluating hard process cross sections.
+ */
+ Energy theHardProcessMass;
+
+ /**
+ * True, if a hard process mass has been set.
+ */
+ bool hardProcessMassSet;
+
+ /**
+ * The width to be used when evaluating hard process cross sections.
+ */
+ Energy theHardProcessWidth;
+
+ /**
+ * True, if a hard process width has been set.
+ */
+ bool hardProcessWidthSet;
+
+ /**
* Upper width cut.
*/
Energy theWidthUpCut;
/**
* Lower width cut.
*/
Energy theWidthLoCut;
/**
* Lifetime.
*/
Length theCTau;
/**
* Three times the charge.
*/
PDT::Charge theCharge;
/**
* 2 times the spin plus one.
*/
PDT::Spin theSpin;
/**
* The colour for this particle.
*/
PDT::Colour theColour;
/**
* A pointer to an object capable to generate a mass for a particle
* of this type.
*/
MassGenPtr theMassGenerator;
/**
* True if the particle is considered stable.
*/
bool isStable;
/**
* A selector of decay modes weighted by the nominal branching
* ratios.
*/
DecaySelector theDecaySelector;
/**
* The set of all decay modes.
*/
DecaySet theDecayModes;
/**
* A pointer to an object capable to generate the branching
* fractions for different decay modes for this particle type. The
* object will be asked to generate branching fractions every time
* the ParticleData object it updated and will modify the branching
* fractions for every particle instance if variableRatio is true.
*/
WidthGeneratorPtr theWidthGenerator;
/**
* Determine whether the branching fractions are allowed to change
* on a particle-by-particle basis.
*/
bool theVariableRatio;
/**
* Pointer to the object corresponding to the antiparticle. Set to
* null if it is its own antiparticle.
*/
tPDPtr theAntiPartner;
/**
* If syncAnti is true all changes to this object will be transfered
* to the antiParticle.
*/
bool syncAnti;
/**
* Helper variable to keep track of the default mass.
*/
Energy theDefMass;
/**
* Helper variable to keep track of the default width.
*/
Energy theDefWidth;
/**
* Helper variable to keep track of the default width cut.
*/
Energy theDefCut;
/**
* Helper variable to keep track of the default lifetime.
*/
Length theDefCTau;
/**
* Helper variable to keep track of the default charge.
*/
PDT::Charge theDefCharge;
/**
* Helper variable to keep track of the default spin.
*/
PDT::Spin theDefSpin;
/**
* Helper variable to keep track of the default colour.
*/
PDT::Colour theDefColour;
/**
* Utility function for the interface.
*/
void setMass(Energy);
/**
* Utility function for the interface.
*/
+ void setHardProcessMass(Energy);
+
+ /**
+ * Reset the hard process mass
+ */
+ string doUnsetHardProcessMass(string);
+
+ /**
+ * Adjust the nominal mass to the hard process mass if a reshuffling
+ * is not desirable.
+ */
+ string doAdjustNominalMass(string);
+
+ /**
+ * Utility function for the interface.
+ */
Energy defMass() const;
/**
* Utility function for the interface.
*/
void setWidth(Energy);
/**
* Utility function for the interface.
*/
+ void setHardProcessWidth(Energy);
+
+ /**
+ * Reset the hard process mass
+ */
+ string doUnsetHardProcessWidth(string);
+
+ /**
+ * Utility function for the interface.
+ */
Energy getWidth() const;
/**
* Utility function for the interface.
*/
Energy defWidth() const;
/**
* Utility function for the interface.
*/
void setCut(Energy);
/**
* Utility function for the interface.
*/
Energy getCut() const;
/**
* Utility function for the interface.
*/
Energy defCut() const;
/**
* Utility function for the interface.
*/
void setUpCut(Energy);
/**
* Utility function for the interface.
*/
Energy getUpCut() const;
/**
* Utility function for the interface.
*/
void setLoCut(Energy);
/**
* Utility function for the interface.
*/
Energy getLoCut() const;
/**
* Utility function for the interface.
*/
void setCTau(Length);
/**
* Utility function for the interface.
*/
Length getCTau() const;
/**
* Utility function for the interface.
*/
Length defCTau() const;
/**
* Utility function for the interface.
*/
void setStable(long);
/**
* Utility function for the interface.
*/
long getStable() const;
/**
* Utility function for the interface.
*/
void setSync(long);
/**
* Utility function for the interface.
*/
long getSync() const;
/**
* Utility function for the interface.
*/
void setVariableRatio(long);
/**
* Utility function for the interface.
*/
long getVariableRatio() const;
/**
* Utility function for the interface.
*/
string doSync(string);
/**
* Utility function for the interface.
*/
void setMassGenerator(MassGenPtr);
/**
* Utility function for the interface.
*/
void setWidthGenerator(WidthGeneratorPtr);
/**
* Utility function for the interface.
*/
void setCharge(int);
/**
* Utility function for the interface.
*/
string ssetCharge(string);
/**
* Utility function for the interface.
*/
int getCharge() const;
/**
* Utility function for the interface.
*/
int defCharge() const;
/**
* Utility function for the interface.
*/
void setSpin(int);
/**
* Utility function for the interface.
*/
int getSpin() const;
/**
* Utility function for the interface.
*/
int defSpin() const;
/**
* Utility function for the interface.
*/
void setColour(long);
/**
* Utility function for the interface.
*/
long getColour() const;
/**
* Utility function for the interface.
*/
long defColour() const;
/**
* Utility function for the interface.
*/
void insDecayModes(DMPtr dm, int);
/**
* Utility function for the interface.
*/
void delDecayModes(int i);
/**
* Utility function for the interface.
*/
vector<DMPtr> getDecayModes() const;
/**
* Utility function for the interface.
*/
string doSelectDecayModes(string);
/**
* Utility function for the interface.
*/
string doPrintDecayModes(string);
/**
* Describe a concrete class with persistent data.
*/
static ClassDescription<ParticleData> initParticleData;
};
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the base classes
* of ParticleData. */
template <>
struct BaseClassTrait<ParticleData,1>: public ClassTraitsType {
/** Typedef of the first base class of ParticleData. */
typedef Interfaced NthBase;
};
/** This template specialization informs ThePEG about the name of the
* ParticleData class. */
template <>
struct ClassTraits<ParticleData>: public ClassTraitsBase<ParticleData> {
/** Return a platform-independent class name */
static string className() { return "ThePEG::ParticleData"; }
};
/** @endcond */
}
#endif /* ThePEG_ParticleData_H */
diff --git a/PDT/StandardMatchers.h b/PDT/StandardMatchers.h
--- a/PDT/StandardMatchers.h
+++ b/PDT/StandardMatchers.h
@@ -1,338 +1,359 @@
// -*- C++ -*-
//
// StandardMatchers.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2011 Leif Lonnblad
//
// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_StandardMatchers_H
#define ThePEG_StandardMatchers_H
// This is the declaration of the AnyMatcher,
#include "Matcher.h"
#include "ThePEG/PDT/EnumParticles.h"
namespace ThePEG {
/** \file StandardMatchers.h
*
* This file declare a set of standard matcher classes. The
* ChargedMatcher, NegativeMatcher, PositiveMatcher, NeutralMatcher,
* BaryonMatcher, MesonMatcher, DiquarkMatcher, LeptonMatcher,
* LightAntiQuarkMatcher, LightQuarkMatcher and
* StandardQCDPartonMatcher classes can be used by themselves (with
* their static functions) or together with the Matcher class to
* define Interfaced objects of the MatcherBase type to be used in the
* Repository. Suitable typedefs are declared for the latter.
*
* @see Matcher
* @see MatcherBase
*/
/**
* A Matcher class which matches any particle.
*/
struct AnyMatcher: public MatcherType {
/** Typedef the class matching the complex conjugate particles. */
typedef AnyMatcher CC;
/** The main static function to check if a given particle type \a pd
matches. */
static bool Check(const ParticleData & pd) {
return bool(pd.id());
}
/** A simplified but unique class name. */
static string className() { return "Any"; }
};
/** Gives a MatcherBase class based on AnyMatcher. */
typedef Matcher<AnyMatcher> MatchAny;
/**
* A Matcher class which matches any charged particle.
*/
struct ChargedMatcher: public MatcherType {
/** Typedef the class matching the complex conjugate particles. */
typedef ChargedMatcher CC;
/** The main static function to check if a given particle type \a pd
matches. */
static bool Check(const ParticleData & pd) {
return PDT::charged(pd.iCharge());
}
/** A simplified but unique class name. */
static string className() { return "Charged"; }
};
/** Gives a MatcherBase class based on ChargedMatcher. */
typedef Matcher<ChargedMatcher> MatchCharged;
struct NegativeMatcher;
/**
* A Matcher class which matches any positively charged particle.
*/
struct PositiveMatcher: public MatcherType {
/** Typedef the class matching the complex conjugate particles. */
typedef NegativeMatcher CC;
/** The main static function to check if a given particle type \a pd
matches. */
static bool Check(const ParticleData & pd) {
return PDT::positive(pd.iCharge());
}
/** A simplified but unique class name. */
static string className() { return "Positive"; }
};
/** Gives a MatcherBase class based on PositiveMatcher. */
typedef Matcher<PositiveMatcher> MatchPositive;
/**
* A Matcher class which matches any uncharged particle.
*/
struct NeutralMatcher: public MatcherType {
/** Typedef the class matching the complex conjugate particles. */
typedef NeutralMatcher CC;
/** The main static function to check if a given particle type \a pd
matches. */
static bool Check(const ParticleData & pd) {
return pd.iCharge() == PDT::Charge0;
}
/** A simplified but unique class name. */
static string className() { return "Neutral"; }
};
/** Gives a MatcherBase class based on NeutralMatcher. */
typedef Matcher<NeutralMatcher> MatchNeutral;
/**
* A Matcher class which matches any negatively charged particle.
*/
struct NegativeMatcher: public MatcherType {
/** Typedef the class matching the complex conjugate particles. */
typedef PositiveMatcher CC;
/** The main static function to check if a given particle type \a pd
matches. */
static bool Check(const ParticleData & pd) {
return PDT::negative(pd.iCharge());
}
/** A simplified but unique class name. */
static string className() { return "Negative"; }
};
/** Gives a MatcherBase class based on NegativeMatcher. */
typedef Matcher<NegativeMatcher> MatchNegative;
/**
* A Matcher class which matches any baryon.
*/
struct BaryonMatcher: public MatcherType {
/** Typedef the class matching the complex conjugate particles. */
typedef BaryonMatcher CC;
/** The main static function to check if a given particle type \a pd
matches. */
static bool Check(const ParticleData & pd) { return Check(pd.id()); }
/** The main static function to check if a given particle with type
\a id matches. */
static bool Check(long id) {
return (id/10)%10 && (id/100)%10 && (id/1000)%10;
}
/** A simplified but unique class name. */
static string className() { return "Baryon"; }
};
/** Gives a MatcherBase class based on BaryonMatcher. */
typedef Matcher<BaryonMatcher> MatchBaryon;
/**
* A Matcher class which matches any meson.
*/
struct MesonMatcher: public MatcherType {
/** Typedef the class matching the complex conjugate particles. */
typedef MesonMatcher CC;
/** The main static function to check if a given particle type \a pd
matches. */
static bool Check(const ParticleData & pd) { return Check(pd.id()); }
/** The main static function to check if a given particle with type
\a id matches. */
static bool Check(long id) {
return (id/10)%10 && (id/100)%10 && (id/1000)%10 == 0;
}
/** A simplified but unique class name. */
static string className() { return "Meson"; }
};
/** Gives a MatcherBase class based on MesonMatcher. */
typedef Matcher<MesonMatcher> MatchMeson;
/**
* A Matcher class which matches any (anti-)diquark.
*/
struct DiquarkMatcher: public MatcherType {
/** Typedef the class matching the complex conjugate particles. */
typedef DiquarkMatcher CC;
/** The main static function to check if a given particle type \a pd
matches. */
static bool Check(const ParticleData & pd) { return Check(pd.id()); }
/** The main static function to check if a given particle with type
\a id matches. */
static bool Check(long id) {
return id/10 && (id/10)%10 == 0 && (id/100)%10 && (id/1000)%10;
}
/** A simplified but unique class name. */
static string className() { return "Diquark"; }
};
/** Gives a MatcherBase class based on DiquarkMatcher. */
typedef Matcher<DiquarkMatcher> MatchDiquark;
/**
* A Matcher class which matches any (anti-)quark.
*/
struct QuarkMatcher: public MatcherType {
/** Typedef the class matching the complex conjugate particles. */
typedef QuarkMatcher CC;
/** The main static function to check if a given particle type \a pd
matches. */
static bool Check(const ParticleData & pd) { return Check(pd.id()); }
/** The main static function to check if a given particle with type
\a id matches. */
static bool Check(long id) {
return id && abs(id) < 10;
}
/** A simplified but unique class name. */
static string className() { return "Quark"; }
};
/** Gives a MatcherBase class based on QuarkMatcher. */
typedef Matcher<QuarkMatcher> MatchQuark;
/**
* A Matcher class which matches any lepton.
*/
struct LeptonMatcher: public MatcherType {
/** Typedef the class matching the complex conjugate particles. */
typedef LeptonMatcher CC;
/** The main static function to check if a given particle type \a pd
matches. */
static bool Check(const ParticleData & pd) { return Check(pd.id()); }
/** The main static function to check if a given particle with type
\a id matches. */
static bool Check(long id) {
return abs(id) > 10 && abs(id) <= 20;
}
/** A simplified but unique class name. */
static string className() { return "Lepton"; }
};
/** Gives a MatcherBase class based on LeptonMatcher. */
typedef Matcher<LeptonMatcher> MatchLepton;
+/**
+ * A Matcher class which matches any neutrino.
+ */
+struct NeutrinoMatcher: public MatcherType {
+ /** Typedef the class matching the complex conjugate particles. */
+ typedef NeutrinoMatcher CC;
+ /** The main static function to check if a given particle type \a pd
+ matches. */
+ static bool Check(const ParticleData & pd) { return Check(pd.id()); }
+ /** The main static function to check if a given particle with type
+ \a id matches. */
+ static bool Check(long id) {
+ return abs(id)-10 >= 2 && abs(id)-10 <= 8 && abs(id)%2 == 0;
+ }
+ /** A simplified but unique class name. */
+ static string className() { return "Neutrino"; }
+};
+
+/** Gives a MatcherBase class based on VectorMesonMatcher. */
+typedef Matcher<NeutrinoMatcher> MatchNeutrino;
+
struct LightAntiQuarkMatcher;
/**
* A Matcher class which matches any light quark (d,u or s).
*/
struct LightQuarkMatcher: public MatcherType {
/** Typedef the class matching the complex conjugate particles. */
typedef LightAntiQuarkMatcher CC;
/** The main static function to check if a given particle type \a pd
matches. */
static bool Check(const ParticleData & pd) { return Check(pd.id()); }
/** The main static function to check if a given particle with type
\a id matches. */
static bool Check(long id) {
return id > 0 && id < 4 ;
}
/** A simplified but unique class name. */
static string className() { return "LightQuark"; }
};
/** Gives a MatcherBase class based on LightQuarkMatcher. */
typedef Matcher<LightQuarkMatcher> MatchLightQuark;
/**
* A Matcher class which matches any light anti-quark
* (\f$\bar{\mbox{d}}\f$,\f$\bar{\mbox{u}}\f$ or
* \f$\bar{\mbox{s}}\f$).
*/
struct LightAntiQuarkMatcher: public MatcherType {
/** Typedef the class matching the complex conjugate particles. */
typedef LightQuarkMatcher CC;
/** The main static function to check if a given particle type \a pd
matches. */
static bool Check(const ParticleData & pd) { return Check(pd.id()); }
/** The main static function to check if a given particle with type
\a id matches. */
static bool Check(long id) {
return id < 0 && id > -4 ;
}
/** A simplified but unique class name. */
static string className() { return "LightAntiQuark"; }
};
/** Gives a MatcherBase class based on LightAntiQuarkMatcher. */
typedef Matcher<LightAntiQuarkMatcher> MatchLightAntiQuark;
/**
* A Matcher class which matches any standard QCD parton, ie. gluons
* and quarks up to bottom.
*/
struct StandardQCDPartonMatcher: public MatcherType {
/** Typedef the class matching the complex conjugate particles. */
typedef StandardQCDPartonMatcher CC;
/** The main static function to check if a given particle type \a pd
matches. */
static bool Check(const ParticleData & pd) { return Check(pd.id()); }
/** The main static function to check if a given particle with type
\a id matches. */
static bool Check(long id) {
return id && ( abs(id) <= 5 || id == ParticleID::g );
}
/** A simplified but unique class name. */
static string className() { return "StandardQCDParton"; }
};
/** Gives a MatcherBase class based on StandardQCDPartonMatcher. */
typedef Matcher<StandardQCDPartonMatcher> MatchStandardQCDParton;
/**
* A Matcher class which matches any pseudo scalar meson.
*/
struct PseudoScalarMesonMatcher: public MatcherType {
/** Typedef the class matching the complex conjugate particles. */
typedef PseudoScalarMesonMatcher CC;
/** The main static function to check if a given particle type \a pd
matches. */
static bool Check(const ParticleData & pd) { return Check(pd.id()); }
/** The main static function to check if a given particle with type
\a id matches. */
static bool Check(long id) {
return ( (abs(id)/1000)%1 == 0 && abs(id) > 100 && abs(id)%10 == 1 ) ||
( id == ParticleID::K_L0 || id == ParticleID::K_S0 );
}
/** A simplified but unique class name. */
static string className() { return "PseudoScalarMeson"; }
};
/** Gives a MatcherBase class based on PseudoScalarMesonMatcher. */
typedef Matcher<PseudoScalarMesonMatcher> MatchPseudoScalarMeson;
/**
* A Matcher class which matches any vector meson.
*/
struct VectorMesonMatcher: public MatcherType {
/** Typedef the class matching the complex conjugate particles. */
typedef VectorMesonMatcher CC;
/** The main static function to check if a given particle type \a pd
matches. */
static bool Check(const ParticleData & pd) { return Check(pd.id()); }
/** The main static function to check if a given particle with type
\a id matches. */
static bool Check(long id) {
return (abs(id)/1000)%1 == 0 && abs(id) > 100 && abs(id)%10 == 3;
}
/** A simplified but unique class name. */
static string className() { return "VectorMeson"; }
};
/** Gives a MatcherBase class based on VectorMesonMatcher. */
typedef Matcher<VectorMesonMatcher> MatchVectorMeson;
}
#endif /* ThePEG_StandardMatchers_H */
diff --git a/Utilities/DynamicLoader.cc b/Utilities/DynamicLoader.cc
--- a/Utilities/DynamicLoader.cc
+++ b/Utilities/DynamicLoader.cc
@@ -1,156 +1,159 @@
// -*- C++ -*-
//
// DynamicLoader.cc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2011 Leif Lonnblad
//
// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the DynamicLoader class.
//
// macro is passed in from -D compile flag
#ifndef THEPEG_PKGLIBDIR
#error Makefile.am needs to define THEPEG_PKGLIBDIR
#endif
#include "DynamicLoader.h"
#include "ThePEG/Utilities/StringUtils.h"
#include "ThePEG/Config/algorithm.h"
#include "config.h"
#ifdef ThePEG_HAS_DLOPEN
#include <dlfcn.h>
#endif
#include <cstdlib>
#ifdef ThePEG_HAS_FENV
#include <fenv.h>
#endif
using namespace ThePEG;
void DynamicLoader::dlname(string sofile) {
if ( StringUtils::suffix(sofile) == "so" ) {
string lafile = StringUtils::remsuf(sofile) + ".la";
ifstream is(lafile.c_str());
string line;
while ( getline(is, line) ) {
if ( line.find("dlname='") != string::npos ) {
int pos = line.find('\'') + 1;
int l = line.rfind('\'') - pos;
sofile = StringUtils::basename(sofile);
versionMap[sofile] = line.substr(pos, l);
}
}
} else if ( StringUtils::suffix(StringUtils::remsuf(sofile)) == "so" ) {
versionMap[StringUtils::basename(StringUtils::remsuf(sofile))] =
StringUtils::basename(sofile);
}
}
string DynamicLoader::dlnameversion(string libs) {
string ret;
do {
string soname = StringUtils::car(libs);
string dir = StringUtils::dirname(soname);
if ( dir.length() ) dir += '/';
libs = StringUtils::cdr(libs);
if ( versionMap.find(StringUtils::basename(soname)) != versionMap.end() )
ret += dir + versionMap[StringUtils::basename(soname)] + " ";
else
ret += soname + " ";
} while ( libs.length() );
return StringUtils::stripws(ret);
}
bool DynamicLoader::loadcmd(string file) {
#ifdef ThePEG_HAS_DLOPEN
dlname(file);
#ifdef ThePEG_HAS_FENV
int oldfpe = fegetexcept();
#endif
bool ret = dlopen(file.c_str(), RTLD_LAZY|RTLD_GLOBAL) != 0;
#ifdef ThePEG_HAS_FENV
feenableexcept(oldfpe);
#endif
if ( !ret ) lastErrorMessage += string(dlerror()) + string("\n");
return ret;
#else
#error ThePEG can only be run on platforms which support
#error dynamic loading.
return false;
#endif
}
void DynamicLoader::appendPath(string path) {
if ( path[path.size()-1] != '/' ) path += '/';
paths.push_back(path);
apppaths.push_back(path);
}
void DynamicLoader::prependPath(string path) {
if ( path[path.size()-1] != '/' ) path += '/';
paths.insert(paths.begin(), path);
prepaths.push_back(path);
}
bool DynamicLoader::load(string name) {
lastErrorMessage = "";
static set<string> loaded;
if ( loaded.find(name) != loaded.end() ) return true;
loaded.insert(name);
bool success = false;
- if ( name[0] == '/' ) success = loadcmd(name);
+ const string name_dylib = StringUtils::remsuf(name) + ".dylib";
+ if ( name[0] == '/' ) {
+ success = loadcmd(name) || loadcmd(name_dylib);
+ }
else {
for ( unsigned i = 0; i < paths.size(); ++i ) {
string path = paths[i];
if ( path[path.size() - 1] != '/' ) path += '/';
- if ( loadcmd(path + name) ) {
+ if ( loadcmd(path + name) || loadcmd(path + name_dylib) ) {
success = true;
break;
}
}
}
- if ( success || loadcmd(name) ) {
+ if ( success || loadcmd(name) || loadcmd(name_dylib) ) {
lastErrorMessage = "";
return true;
}
loaded.erase(name);
return false;
}
const vector<string> & DynamicLoader::appendedPaths() {
return apppaths;
}
const vector<string> & DynamicLoader::prependedPaths() {
return prepaths;
}
const vector<string> & DynamicLoader::allPaths() {
return paths;
}
vector<string> DynamicLoader::paths = DynamicLoader::defaultPaths();
vector<string> DynamicLoader::prepaths = vector<string>();
vector<string> DynamicLoader::apppaths = vector<string>();
string DynamicLoader::lastErrorMessage;
map<string,string> DynamicLoader::versionMap;
vector<string> DynamicLoader::defaultPaths() {
vector<string> vpaths;
// macro is passed in from -D compile flag
string instpath = THEPEG_PKGLIBDIR;
vpaths.push_back(instpath);
vpaths.push_back(".");
return vpaths;
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 2:01 PM (11 h, 40 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023053
Default Alt Text
(136 KB)

Event Timeline