Page MenuHomeHEPForge
Files F8725611

No OneTemporary
Actions

View Options

diff --git a/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.cc b/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.cc
--- a/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.cc
+++ b/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.cc
@@ -1,145 +1,157 @@
// -*- C++ -*-
//
// MatchboxPhasespace.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MatchboxPhasespace class.
//
#include "MatchboxPhasespace.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
using namespace Herwig;
MatchboxPhasespace::MatchboxPhasespace() {}
MatchboxPhasespace::~MatchboxPhasespace() {}
void MatchboxPhasespace::cloneDependencies(const std::string&) {
}
void MatchboxPhasespace::dumpInfo(const string& prefix) const {
generator()->log() << prefix << fullName()
<< " [" << this << "]\n";
generator()->log() << prefix << " | XComb " << lastXCombPtr()
<< " for ";
if ( lastXCombPtr() ) {
for ( cPDVector::const_iterator p = lastXComb().mePartonData().begin();
p != lastXComb().mePartonData().end(); ++p ) {
generator()->log() << (**p).PDGName() << " ";
}
}
generator()->log() << "\n";
}
pair<double,Lorentz5Momentum>
MatchboxPhasespace::timeLikeWeight(const Tree2toNDiagram& diag,
- int branch) const {
+ int branch, double flatCut) const {
pair<int,int> children = diag.children(branch);
if ( children.first == -1 ) {
return make_pair(1.,meMomenta()[diag.externalId(branch)]);
}
pair<double,Lorentz5Momentum> res
- = timeLikeWeight(diag,children.first);
+ = timeLikeWeight(diag,children.first,flatCut);
pair<double,Lorentz5Momentum> other
- = timeLikeWeight(diag,children.second);
+ = timeLikeWeight(diag,children.second,flatCut);
res.first *= other.first;
res.second += other.second;
Energy2 mass2 = sqr(diag.allPartons()[branch]->mass());
Energy2 width2 = sqr(diag.allPartons()[branch]->width());
- res.first /=
- sqr((res.second.m2()-mass2)/lastSHat()) +
- mass2*width2/sqr(lastSHat());
+ if ( width2 == ZERO ) {
+ if ( abs((res.second.m2()-mass2)/lastSHat()) > flatCut )
+ res.first /=
+ abs((res.second.m2()-mass2)/lastSHat());
+ } else {
+ res.first /=
+ sqr((res.second.m2()-mass2)/lastSHat()) +
+ mass2*width2/sqr(lastSHat());
+ }
return res;
}
double MatchboxPhasespace::spaceLikeWeight(const Tree2toNDiagram& diag,
const Lorentz5Momentum& incoming,
- int branch) const {
+ int branch, double flatCut) const {
if ( branch == -1 )
return 1.;
pair<int,int> children = diag.children(branch);
pair<double,Lorentz5Momentum> res =
- timeLikeWeight(diag,children.second);
+ timeLikeWeight(diag,children.second,flatCut);
if ( children.first == diag.nSpace() - 1 ) {
return res.first;
}
res.second = incoming - res.second;
Energy2 mass2 = sqr(diag.allPartons()[branch]->mass());
Energy2 width2 = sqr(diag.allPartons()[branch]->width());
- res.first /=
- sqr((res.second.m2()-mass2)/lastSHat()) +
- mass2*width2/sqr(lastSHat());
+ if ( width2 == ZERO ) {
+ if ( abs((res.second.m2()-mass2)/lastSHat()) > flatCut )
+ res.first /=
+ abs((res.second.m2()-mass2)/lastSHat());
+ } else {
+ res.first /=
+ sqr((res.second.m2()-mass2)/lastSHat()) +
+ mass2*width2/sqr(lastSHat());
+ }
return
- res.first * spaceLikeWeight(diag,res.second,children.first);
+ res.first * spaceLikeWeight(diag,res.second,children.first,flatCut);
}
-void MatchboxPhasespace::fillDiagramWeights() {
+void MatchboxPhasespace::fillDiagramWeights(double flatCut) {
diagramWeights.clear();
for ( StandardXComb::DiagramVector::const_iterator d =
lastXComb().diagrams().begin(); d != lastXComb().diagrams().end(); ++d ) {
diagramWeights[(**d).id()] =
- spaceLikeWeight(dynamic_cast<const Tree2toNDiagram&>(**d),meMomenta()[0],0);
+ spaceLikeWeight(dynamic_cast<const Tree2toNDiagram&>(**d),meMomenta()[0],0,flatCut);
}
}
Selector<MEBase::DiagramIndex>
MatchboxPhasespace::selectDiagrams(const MEBase::DiagramVector& diags) const {
Selector<MEBase::DiagramIndex> ret;
for ( MEBase::DiagramIndex d = 0; d < diags.size(); ++d ) {
ret.insert(diagramWeight(dynamic_cast<const Tree2toNDiagram&>(*diags[d])),d);
}
return ret;
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeAbstractNoPIOClass<MatchboxPhasespace,HandlerBase>
describeMatchboxPhasespace("Herwig::MatchboxPhasespace", "HwMatchbox.so");
void MatchboxPhasespace::Init() {
static ClassDocumentation<MatchboxPhasespace> documentation
("MatchboxPhasespace defines an abstract interface to a phase "
"space generator.");
}
diff --git a/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.h b/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.h
--- a/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.h
+++ b/MatrixElement/Matchbox/Phasespace/MatchboxPhasespace.h
@@ -1,206 +1,206 @@
// -*- C++ -*-
//
// MatchboxPhasespace.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_MatchboxPhasespace_H
#define HERWIG_MatchboxPhasespace_H
//
// This is the declaration of the MatchboxPhasespace class.
//
#include "ThePEG/Handlers/StandardXComb.h"
#include "ThePEG/Handlers/HandlerBase.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief Wrap around a vector of random numbers to behave as a stream
* of those.
*/
struct StreamingRnd {
/**
* The random numbers
*/
const double* numbers;
/**
* The number of random numbers available.
*/
size_t nRnd;
/**
* Default constructor.
*/
StreamingRnd()
: numbers(0), nRnd(0) {}
/**
* Construct from random numbers.
*/
explicit StreamingRnd(const double* newNumbers,
size_t n)
: numbers(newNumbers), nRnd(n) {}
/**
* Return next random number
*/
inline double operator()() {
assert(numbers && nRnd > 0);
const double ret = numbers[0];
++numbers; --nRnd;
return ret;
}
};
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief MatchboxPhasespace defines an abstract interface to a phase
* space generator.
*
*/
class MatchboxPhasespace:
public HandlerBase, public LastXCombInfo<StandardXComb> {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
MatchboxPhasespace();
/**
* The destructor.
*/
virtual ~MatchboxPhasespace();
//@}
public:
/**
* Prepare a phase space generator for the given xcomb object.
*/
virtual void prepare(tStdXCombPtr, bool verbose = false) = 0;
/**
* Generate a phase space point and return its weight.
*/
virtual double generateKinematics(const double*,
vector<Lorentz5Momentum>& momenta) = 0;
/**
* Return the number of random numbers required to produce a given
* multiplicity final state.
*/
virtual int nDim(int nFinal) const = 0;
/**
* Return true, if this phasespace generator will generate incoming
* partons itself.
*/
virtual bool haveX1X2() const { return false; }
/**
* Return true, if this phase space generator expects
* the incoming partons in their center-of-mass system
*/
virtual bool wantCMS() const { return true; }
/**
* Fill a diagram selector for the last phase space point.
*/
virtual Selector<MEBase::DiagramIndex> selectDiagrams(const MEBase::DiagramVector&) const;
/**
* Return the momentum and weight appropriate to the given timelike
* branch of the diagram.
*/
pair<double,Lorentz5Momentum> timeLikeWeight(const Tree2toNDiagram& diag,
- int branch) const;
+ int branch, double flatCut) const;
/**
* Return the weight appropriate to the given spacelike branch of
* the diagram.
*/
double spaceLikeWeight(const Tree2toNDiagram& diag,
const Lorentz5Momentum& incoming,
- int branch) const;
+ int branch, double flatCut) const;
/**
* Return the weight appropriate to the given diagram.
*/
double diagramWeight(const Tree2toNDiagram& diag) const {
assert( !diagramWeights.empty() );
return diagramWeights.find(diag.id())->second;
}
/**
* Fill the diagram weights.
*/
- void fillDiagramWeights();
+ void fillDiagramWeights(double flatCut = 0.0);
/**
* Clone this phase space generator.
*/
Ptr<MatchboxPhasespace>::ptr cloneMe() const {
return dynamic_ptr_cast<Ptr<MatchboxPhasespace>::ptr>(clone());
}
/**
* Clone the dependencies, using a given prefix.
*/
virtual void cloneDependencies(const std::string& prefix = "");
/**
* Dump xcomb hierarchies.
*/
void dumpInfo(const string& prefix = "") const;
public:
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
private:
/**
* The diagram weights indexed by diagram id.
*/
map<int,double> diagramWeights;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MatchboxPhasespace & operator=(const MatchboxPhasespace &);
};
}
#endif /* HERWIG_MatchboxPhasespace_H */
diff --git a/MatrixElement/Matchbox/Phasespace/PhasespaceHelpers.cc b/MatrixElement/Matchbox/Phasespace/PhasespaceHelpers.cc
--- a/MatrixElement/Matchbox/Phasespace/PhasespaceHelpers.cc
+++ b/MatrixElement/Matchbox/Phasespace/PhasespaceHelpers.cc
@@ -1,514 +1,516 @@
// -*- C++ -*-
//
// PhasespaceHelpers.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#include "PhasespaceHelpers.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/RandomHelpers.h"
using namespace Herwig;
using namespace Herwig::PhasespaceHelpers;
using namespace Herwig::RandomHelpers;
Energy PhasespaceInfo::generateMass(tcPDPtr data,
const pair<Energy,Energy>& range) {
double xlow = sqr(range.first)/sHat;
if ( range.first < ZERO )
xlow = -xlow;
double xup = sqr(range.second)/sHat;
if ( range.second < ZERO )
xup = -xup;
double mu = sqr(data->mass())/sHat;
double gamma = sqr(data->width())/sHat;
double r = rnd();
pair<double,double> event;
if ( gamma == 0. ) {
if ( mu < xlow || mu > xup ) {
if ( abs(xlow-mu) < xc )
xlow = mu;
if ( abs(xup-mu) < xc )
xup = mu;
}
if ( mu < xlow || mu > xup ) {
event =
generate(inverse(mu,xlow,xup),r);
} else {
pair<double,double> pLeft(xlow,xlow < mu-x0 ? mu-x0 : xlow);
pair<double,double> pRight(xup > mu+x0 ? mu+x0 : xup,xup);
pair<double,double> fLeft(pLeft.second < mu-x0 ? mu-x0 : pLeft.second,mu-xc);
if ( fLeft.first >= fLeft.second )
fLeft.first = fLeft.second;
pair<double,double> fRight(mu+xc,pRight.first > mu+x0 ? mu+x0 : pRight.first);
if ( fRight.first >= fRight.second )
fRight.second = fRight.first;
if ( pLeft.first != pLeft.second &&
fLeft.first != fLeft.second &&
fRight.first != fRight.second &&
pRight.first != pRight.second ) {
event =
generate((piecewise(),
inverse(mu,pLeft.first,pLeft.second),
match(flat(fLeft.first,fLeft.second))) +
match((piecewise(),
flat(fRight.first,fRight.second),
match(inverse(mu,pRight.first,pRight.second)))),
r);
} else if ( pLeft.first == pLeft.second &&
fLeft.first != fLeft.second &&
fRight.first != fRight.second &&
pRight.first != pRight.second ) {
event =
generate(flat(fLeft.first,fLeft.second) +
match((piecewise(),
flat(fRight.first,fRight.second),
match(inverse(mu,pRight.first,pRight.second)))),
r);
} else if ( pLeft.first != pLeft.second &&
fLeft.first != fLeft.second &&
fRight.first != fRight.second &&
pRight.first == pRight.second ) {
event =
generate((piecewise(),
inverse(mu,pLeft.first,pLeft.second),
match(flat(fLeft.first,fLeft.second))) +
match(flat(fRight.first,fRight.second)),
r);
} else if ( pLeft.first == pLeft.second &&
fLeft.first == fLeft.second &&
fRight.first != fRight.second &&
pRight.first != pRight.second ) {
event =
generate((piecewise(),
flat(fRight.first,fRight.second),
match(inverse(mu,pRight.first,pRight.second))),
r);
} else if ( pLeft.first != pLeft.second &&
fLeft.first != fLeft.second &&
fRight.first == fRight.second &&
pRight.first == pRight.second ) {
event =
generate((piecewise(),
inverse(mu,pLeft.first,pLeft.second),
match(flat(fLeft.first,fLeft.second))),
r);
} else if ( pLeft.first == pLeft.second &&
fLeft.first != fLeft.second &&
fRight.first != fRight.second &&
pRight.first == pRight.second ) {
event =
generate(flat(fLeft.first,fLeft.second) +
match(flat(fRight.first,fRight.second)),
r);
} else if ( pLeft.first == pLeft.second &&
fLeft.first == fLeft.second &&
fRight.first != fRight.second &&
pRight.first == pRight.second ) {
event =
generate(flat(fRight.first,fRight.second),
r);
} else if ( pLeft.first == pLeft.second &&
fLeft.first != fLeft.second &&
fRight.first == fRight.second &&
pRight.first == pRight.second ) {
event =
generate(flat(fLeft.first,fLeft.second),
r);
} else if ( pLeft.first == pLeft.second &&
fLeft.first == fLeft.second &&
fRight.first == fRight.second &&
pRight.first == pRight.second ) {
throw Veto();
} else assert(false);
}
} else {
event = generate(breitWigner(mu,gamma,xlow,xup),r);
}
weight *= event.second;
Energy res = sqrt(abs(event.first)*sHat);
if ( event.first < 0. )
res = -res;
return res;
}
Lorentz5Momentum PhasespaceInfo::generateKt(const Lorentz5Momentum& p1,
const Lorentz5Momentum& p2,
Energy pt) {
- if ( (p1+p2).m2() <= ZERO ) {
- cerr << "cannot boost ... " << ((p1+p2).m2()/GeV2) << "\n";
- throw Veto();
- }
- Boost beta = (p1+p2).findBoostToCM();
-
- Lorentz5Momentum p1c = p1;
+ double phi = 2.*Constants::pi*rnd();
- if (beta.mag2() > Constants::epsilon) {
- p1c.boost(beta);
- }
-
- Lorentz5Momentum k (0.*GeV,0.*GeV,0.*GeV,0.*GeV);
-
- double ct = p1c.vect().unit().z();
- double st = sqrt(1.-ct*ct);
-
- double phi = 2.*Constants::pi*rnd();
- weight *= 2.*Constants::pi;
- double cphi = cos(phi);
- double sphi = sqrt(1.-cphi*cphi);
- if (phi > Constants::pi) sphi = -sphi;
-
- // adding an output somewhere around here results in the correct behaviour
- // when compiling with g++-4.6.1 -O3
- // cerr << "bla\n" << flush;
- if (st > Constants::epsilon) {
- double cchi = p1c.vect().unit().x()/st;
- double schi = p1c.vect().unit().y()/st;
- k.setX((cphi*cchi*ct-sphi*schi)*pt);
- k.setY((cphi*schi*ct+sphi*cchi)*pt);
- k.setZ(-cphi*st*pt);
- } else {
- k.setX(pt*cphi);
- k.setY(pt*sphi);
- k.setZ(0.*GeV);
- }
-
- if (beta.mag2() > Constants::epsilon)
- k.boost(-beta);
-
- return k;
+ Lorentz5Momentum P = p1 + p2;
+
+ Energy2 Q2 = abs(P.m2());
+
+ Lorentz5Momentum Q =
+ Lorentz5Momentum(ZERO,ZERO,ZERO,sqrt(Q2),sqrt(Q2));
+
+ bool boost =
+ abs((P-Q).vect().mag2()/GeV2) > 1e-10 ||
+ abs((P-Q).t()/GeV) > 1e-5;
+
+ Lorentz5Momentum inFrame1;
+ if ( boost )
+ inFrame1 = p1 + ((P*p1-Q*p1)/(P*Q-Q.mass2()))*(P-Q);
+ else
+ inFrame1 = p1;
+
+ Energy ptx = inFrame1.x();
+ Energy pty = inFrame1.y();
+ Energy q = 2.*inFrame1.z();
+
+ Energy Qp = sqrt(4.*(sqr(ptx)+sqr(pty))+sqr(q));
+ Energy Qy = sqrt(4.*sqr(pty)+sqr(q));
+
+ double cPhi = cos(phi);
+ double sPhi = sqrt(1.-sqr(cPhi));
+ if ( phi > Constants::pi )
+ sPhi = -sPhi;
+
+ Lorentz5Momentum kt;
+
+ kt.setT(ZERO);
+ kt.setX(pt*Qy*cPhi/Qp);
+ kt.setY(-pt*(4*ptx*pty*cPhi/Qp+q*sPhi)/Qy);
+ kt.setZ(2.*pt*(-ptx*q*cPhi/Qp + pty*sPhi)/Qy);
+
+ if ( boost )
+ kt = kt + ((P*kt-Q*kt)/(P*Q-Q.mass2()))*(P-Q);
+ kt.setMass(-pt);
+ kt.rescaleRho();
+
+ return kt;
}
void PhasespaceTree::setup(const Tree2toNDiagram& diag,
int pos) {
pair<int,int> dchildren =
diag.children(pos);
data = diag.allPartons()[pos];
spacelike = pos < diag.nSpace();
if ( pos == 0 )
externalId = 0;
if ( dchildren.first == -1 ) {
externalId = diag.externalId(pos);
leafs.insert(externalId);
return;
}
children.push_back(PhasespaceTree());
children.back().setup(diag,dchildren.first);
children.push_back(PhasespaceTree());
children.back().setup(diag,dchildren.second);
if ( !children[0].children.empty() &&
children[1].children.empty() &&
!spacelike )
swap(children[0],children[1]);
if ( spacelike &&
!children[0].spacelike )
swap(children[0],children[1]);
copy(children[0].leafs.begin(),children[0].leafs.end(),
inserter(leafs,leafs.begin()));
copy(children[1].leafs.begin(),children[1].leafs.end(),
inserter(leafs,leafs.begin()));
}
void PhasespaceTree::init(const vector<Lorentz5Momentum>& meMomenta) {
if ( children.empty() ) {
massRange.first = meMomenta[externalId].mass();
massRange.second = massRange.first;
momentum.setMass(meMomenta[externalId].mass());
if ( externalId == 1 )
momentum = meMomenta[1];
return;
}
children[0].init(meMomenta);
children[1].init(meMomenta);
if ( !children[0].spacelike &&
!children[1].spacelike ) {
massRange.first =
children[0].massRange.first +
children[1].massRange.first;
}
}
void PhasespaceTree::generateKinematics(PhasespaceInfo& info,
vector<Lorentz5Momentum>& meMomenta) {
if ( externalId == 0 ) {
init(meMomenta);
momentum = meMomenta[0];
}
if ( children.empty() ) {
if ( externalId != 1 )
meMomenta[externalId] = momentum;
return;
}
// s-channel
if ( externalId == 0 &&
children[0].externalId == 1 ) {
children[1].momentum = meMomenta[0] + meMomenta[1];
children[1].momentum.setMass(info.sqrtSHat);
children[1].momentum.rescaleEnergy();
children[1].generateKinematics(info,meMomenta);
return;
}
if ( !spacelike ) {
Energy mij = momentum.mass();
Energy mi,mj;
// work out the mass for the first child
if ( !children[0].children.empty() ) {
Energy sumOthers = ZERO;
for ( size_t k = 2; k < meMomenta.size(); ++k )
if ( children[1].leafs.find(k) != children[1].leafs.end() )
sumOthers += meMomenta[k].mass();
children[0].massRange.second = momentum.mass() - sumOthers;
if ( children[0].massRange.second < children[0].massRange.first )
throw Veto();
if ( children[0].massRange.second > momentum.mass() )
throw Veto();
mi = info.generateMass(children[0].data,children[0].massRange);
children[0].momentum.setMass(mi);
} else {
mi = children[0].momentum.mass();
}
// work out the mass for the second child
if ( !children[1].children.empty() ) {
children[1].massRange.second = momentum.mass()-children[0].momentum.mass();
if ( children[1].massRange.second < children[1].massRange.first )
throw Veto();
mj = info.generateMass(children[1].data,children[1].massRange);
children[1].momentum.setMass(mj);
} else {
mj = children[1].momentum.mass();
}
Energy2 mij2 = sqr(mij);
Energy2 mi2 = sqr(mi);
Energy2 mj2 = sqr(mj);
// perform the decay
Energy4 lambda2 = sqr(mij2-mi2-mj2)-4.*mi2*mj2;
if ( lambda2 < ZERO )
throw Veto();
Energy2 lambda = sqrt(lambda2);
double phi = 2.*Constants::pi*info.rnd();
double cosPhi = cos(phi);
double sinPhi = sqrt(1.-sqr(cosPhi));
if ( phi > Constants::pi )
sinPhi = -sinPhi;
info.weight *= Constants::pi*lambda/(2.*mij2);
double cosTheta = 2.*info.rnd() - 1.;
double sinTheta = sqrt(1.-sqr(cosTheta));
Energy p = lambda/(2.*mij);
children[0].momentum.setX(p*cosPhi*sinTheta);
children[0].momentum.setY(p*sinPhi*sinTheta);
children[0].momentum.setZ(p*cosTheta);
children[0].momentum.rescaleEnergy();
if ( momentum.m2() <= ZERO ) {
cerr << "cannot boost in decay ... " << (momentum.m2()/GeV2) << "\n";
throw Veto();
}
Boost out = momentum.boostVector();
if ( out.mag2() > Constants::epsilon ) {
children[0].momentum.boost(out);
}
children[1].momentum = momentum - children[0].momentum;
children[1].momentum.setMass(mj);
children[1].momentum.rescaleEnergy();
// go on with next branchings
children[0].generateKinematics(info,meMomenta);
children[1].generateKinematics(info,meMomenta);
return;
}
// get the minimum mass of the `W' system
Energy Wmin = ZERO;
PhasespaceTree* current = &children[0];
while ( !(current->children.empty()) ) {
Wmin += current->children[1].massRange.first;
current = &(current->children[0]);
}
// get the CM energy avaialble
Energy2 s = (momentum+meMomenta[1]).m2();
if ( s <= ZERO )
throw Veto();
// generate a mass for the timelike child
Energy mi;
if ( !children[1].children.empty() ) {
children[1].massRange.second = sqrt(s)-Wmin;
if ( children[1].massRange.second < children[1].massRange.first )
throw Veto();
mi = info.generateMass(children[1].data,children[1].massRange);
children[1].momentum.setMass(mi);
} else {
mi = children[1].momentum.mass();
}
Energy2 mi2 = sqr(mi);
// wether or not this is the last 2->2 scatter
bool lastScatter
= children[0].children[0].children.empty();
// `W' mass relevant for the other boundaries
Energy MW = Wmin;
// generate a mass for second outgoing leg, if needed
if ( lastScatter )
if ( !children[0].children[1].children.empty() ) {
// get the maximum `W' mass
Energy Wmax = sqrt(s)-children[1].momentum.mass();
children[0].children[1].massRange.second = Wmax;
if ( children[0].children[1].massRange.second <
children[0].children[1].massRange.first )
throw Veto();
MW = info.generateMass(children[0].children[1].data,
children[0].children[1].massRange);
children[0].children[1].momentum.setMass(MW);
}
Energy2 MW2 = sqr(MW);
Energy ma = momentum.mass();
Energy2 ma2 = sqr(ma);
if ( ma < ZERO )
ma2 = -ma2;
Energy mb = meMomenta[1].mass();
Energy2 mb2 = sqr(mb);
// pick the ys variable
Energy2 ys = ZERO;
if ( !lastScatter ) {
ys = info.rnd()*(sqr(sqrt(s)-mi)-MW2);
info.weight *= (sqr(sqrt(s)-mi)-MW2)/info.sHat;
}
Energy4 lambda2 = sqr(s-ma2-mb2)-4.*ma2*mb2;
if ( lambda2 <= ZERO ) {
throw Veto();
}
Energy2 lambda = sqrt(lambda2);
info.weight *= info.sHat/(4.*lambda);
// get the boundaries on the momentum transfer
Energy4 rho2 = sqr(s-ys-MW2-mi2)-4.*mi2*(ys+MW2);
if ( rho2 < ZERO )
throw Veto();
Energy2 rho = sqrt(rho2);
Energy4 tau2 =
ys*(ma2-mb2+s)
- sqr(s)+s*(ma2+mb2+mi2+MW2)-(mi2-MW2)*(ma2-mb2);
pair<Energy2,Energy2> tBounds
((tau2-rho*lambda)/(2.*s),(tau2+rho*lambda)/(2.*s));
children[0].massRange.first = sqrt(abs(tBounds.first));
if ( tBounds.first < ZERO )
children[0].massRange.first = -children[0].massRange.first;
children[0].massRange.second = sqrt(abs(tBounds.second));
if ( tBounds.second < ZERO )
children[0].massRange.second = -children[0].massRange.second;
// generate a momentum transfer
Energy mai = info.generateMass(children[0].data,children[0].massRange);
children[0].momentum.setMass(mai);
Energy2 t = sqr(mai);
if ( mai < ZERO )
t = -t;
Energy2 u = -s -t + ys + ma2 + mb2 + mi2 + MW2;
Energy2 st = s - ma2 - mb2;
Energy2 tt = t - mi2 - ma2;
Energy2 ut = u - mi2 - mb2;
// get the timelike momentum
double xa = (-st*ut+2.*mb2*tt)/lambda2;
double xb = (-st*tt+2.*ma2*ut)/lambda2;
Energy2 pt2 = (st*tt*ut-ma2*sqr(ut)-mb2*sqr(tt)-mi2*sqr(st)+4.*ma2*mb2*mi2)/lambda2;
if ( pt2 < ZERO )
throw Veto();
Energy pt = sqrt(pt2);
children[1].momentum =
xa*momentum + xb*meMomenta[1] + info.generateKt(momentum,meMomenta[1],pt);
children[1].momentum.setMass(mi);
children[1].momentum.rescaleEnergy();
children[0].momentum =
momentum - children[1].momentum;
children[0].momentum.setMass(mai);
bool changeSign = false;
if ( children[0].momentum.t() < ZERO && mai > ZERO ) changeSign = true;
if ( mai < ZERO )
children[0].momentum.rescaleRho();
else
children[0].momentum.rescaleEnergy();
if ( changeSign ) children[0].momentum.setT(-children[0].momentum.t());
children[1].generateKinematics(info,meMomenta);
if ( !lastScatter ) {
children[0].generateKinematics(info,meMomenta);
} else {
children[0].children[1].momentum =
meMomenta[1] + children[0].momentum;
children[0].children[1].momentum.setMass(MW);
children[0].children[1].momentum.rescaleEnergy();
children[0].children[1].generateKinematics(info,meMomenta);
}
}
void PhasespaceTree::put(PersistentOStream& os) const {
os << children.size();
if ( !children.empty() ) {
children[0].put(os);
children[1].put(os);
}
os << data << externalId << leafs << spacelike;
}
void PhasespaceTree::get(PersistentIStream& is) {
size_t nc; is >> nc;
assert(nc == 0 || nc == 2);
if ( nc == 2 ) {
children.resize(2,PhasespaceTree());
children[0].get(is);
children[1].get(is);
}
is >> data >> externalId >> leafs >> spacelike;
}
diff --git a/MatrixElement/Matchbox/Phasespace/TreePhasespace.cc b/MatrixElement/Matchbox/Phasespace/TreePhasespace.cc
--- a/MatrixElement/Matchbox/Phasespace/TreePhasespace.cc
+++ b/MatrixElement/Matchbox/Phasespace/TreePhasespace.cc
@@ -1,170 +1,170 @@
// -*- C++ -*-
//
// TreePhasespace.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the TreePhasespace class.
//
#include "TreePhasespace.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
using namespace Herwig::PhasespaceHelpers;
TreePhasespace::TreePhasespace()
: x0(0.01), xc(1e-4) {
lastPhasespaceInfo.x0 = x0;
lastPhasespaceInfo.xc = xc;
}
TreePhasespace::~TreePhasespace() {}
IBPtr TreePhasespace::clone() const {
return new_ptr(*this);
}
IBPtr TreePhasespace::fullclone() const {
return new_ptr(*this);
}
void TreePhasespace::prepare(tStdXCombPtr xco, bool) {
theLastXComb = xco;
lastChannelsIterator = channelMap().find(lastXCombPtr());
if ( lastChannelsIterator == channelMap().end() ) {
map<Ptr<Tree2toNDiagram>::ptr,PhasespaceTree> channels;
for ( StandardXComb::DiagramVector::const_iterator d =
lastXComb().diagrams().begin(); d != lastXComb().diagrams().end(); ++d ) {
PhasespaceTree tree;
Ptr<Tree2toNDiagram>::ptr diag =
dynamic_ptr_cast<Ptr<Tree2toNDiagram>::ptr>(*d);
tree.setup(*diag);
channels[diag] = tree;
}
channelMap()[lastXCombPtr()] = channels;
lastChannelsIterator = channelMap().find(lastXCombPtr());
}
lastPhasespaceInfo.sHat = lastXComb().lastSHat();
lastPhasespaceInfo.sqrtSHat = sqrt(lastXComb().lastSHat());
lastPhasespaceInfo.weight = 1.;
}
double TreePhasespace::generateKinematics(const double* random,
vector<Lorentz5Momentum>& momenta) {
size_t nchannels = lastXComb().diagrams().size();
map<Ptr<Tree2toNDiagram>::ptr,PhasespaceHelpers::PhasespaceTree>::iterator ds =
lastChannels().begin();
advance(ds,(size_t)(random[0]*nchannels));
Ptr<Tree2toNDiagram>::ptr channel = ds->first;
++random;
lastPhasespaceInfo.rnd.numbers = random;
lastPhasespaceInfo.rnd.nRnd = 3*momenta.size() - 10;
cPDVector::const_iterator pd = mePartonData().begin();
vector<Lorentz5Momentum>::iterator p = momenta.begin();
for ( ; pd != mePartonData().end(); ++pd, ++p )
p->setMass((**pd).mass());
try {
lastChannels()[channel].generateKinematics(lastPhasespaceInfo,momenta);
} catch (Veto) {
return 0.;
}
- fillDiagramWeights();
+ fillDiagramWeights(x0);
double sum = 0.;
for ( map<Ptr<Tree2toNDiagram>::ptr,PhasespaceHelpers::PhasespaceTree>::const_iterator d
= lastChannels().begin(); d != lastChannels().end(); ++d )
sum += diagramWeight(*(d->first));
double piWeight = pow(2.*Constants::pi,(double)(3*(momenta.size()-2)-4));
return nchannels*lastPhasespaceInfo.weight*diagramWeight(*channel)/(sum*piWeight);
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void TreePhasespace::doinit() {
MatchboxPhasespace::doinit();
lastPhasespaceInfo.x0 = x0;
lastPhasespaceInfo.xc = xc;
}
void TreePhasespace::doinitrun() {
MatchboxPhasespace::doinitrun();
lastPhasespaceInfo.x0 = x0;
lastPhasespaceInfo.xc = xc;
}
void TreePhasespace::persistentOutput(PersistentOStream & os) const {
os << theChannelMap << x0 << xc;
}
void TreePhasespace::persistentInput(PersistentIStream & is, int) {
is >> theChannelMap >> x0 >> xc;
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeClass<TreePhasespace,MatchboxPhasespace>
describeHerwigTreePhasespace("Herwig::TreePhasespace", "HwMatchbox.so");
void TreePhasespace::Init() {
static ClassDocumentation<TreePhasespace> documentation
("TreePhasespace is a multichannel phasespace generator "
"adapting to singularity structures as determined from the matrix "
"elements diagrams.");
static Reference<TreePhasespace,TreePhasespaceChannels> interfaceChannelMap
("ChannelMap",
"Set the object storing the channels.",
&TreePhasespace::theChannelMap, false, false, true, false, false);
static Parameter<TreePhasespace,double> interfaceX0
("X0",
"Set the cut below which flat virtuality sampling is imposed.",
&TreePhasespace::x0, 0.01, 0.0, 0,
false, false, Interface::lowerlim);
static Parameter<TreePhasespace,double> interfaceXC
("XC",
"Set the cut below which no virtualities are generated.",
&TreePhasespace::xc, 1e-4, 0.0, 0,
false, false, Interface::lowerlim);
}

File Metadata

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

Event Timeline