Page MenuHomeHEPForge
Files F7878686

No OneTemporary
Actions

View Options

diff --git a/DipoleShower/Kinematics/FFMassiveKinematics.cc b/DipoleShower/Kinematics/FFMassiveKinematics.cc
--- a/DipoleShower/Kinematics/FFMassiveKinematics.cc
+++ b/DipoleShower/Kinematics/FFMassiveKinematics.cc
@@ -1,350 +1,343 @@
// -*- C++ -*-
//
// FFMassiveKinematics.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the FFMassiveKinematics class.
//
#include "FFMassiveKinematics.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "Herwig++/DipoleShower/Base/DipoleSplittingInfo.h"
#include "Herwig++/DipoleShower/Kernels/DipoleSplittingKernel.h"
// TODO: remove after verification
// only for checking for NaN or inf
#include <gsl/gsl_math.h>
using namespace Herwig;
FFMassiveKinematics::FFMassiveKinematics()
: DipoleSplittingKinematics() {}
FFMassiveKinematics::~FFMassiveKinematics() {}
IBPtr FFMassiveKinematics::clone() const {
return new_ptr(*this);
}
IBPtr FFMassiveKinematics::fullclone() const {
return new_ptr(*this);
}
pair<double,double> FFMassiveKinematics::kappaSupport(const DipoleSplittingInfo&) const {
return make_pair(0.0,1.0);
}
pair<double,double> FFMassiveKinematics::xiSupport(const DipoleSplittingInfo& split) const {
double c = sqrt(1.-4.*sqr(IRCutoff()/generator()->maximumCMEnergy()));
if ( split.index().emitterData()->id() == ParticleID::g ) {
if ( split.emissionData()->id() != ParticleID::g )
return make_pair(0.5*(1.-c),0.5*(1.+c));
double b = log((1.+c)/(1.-c));
return make_pair(-b,b);
}
return make_pair(-log(0.5*(1.+c)),-log(0.5*(1.-c)));
}
Energy FFMassiveKinematics::dipoleScale(const Lorentz5Momentum& pEmitter,
const Lorentz5Momentum& pSpectator) const {
return (pEmitter+pSpectator).m();
}
Energy FFMassiveKinematics::ptMax(Energy dScale,
double, double,
const DipoleIndex& ind,
const DipoleSplittingKernel& split) const {
double mui = split.emitter(ind)->mass() / dScale;
double mu = split.emission(ind)->mass() / dScale;
double muj = split.spectator(ind)->mass() / dScale;
double mui2 = sqr( mui ), mu2 = sqr( mu ), muj2 = sqr( muj );
-
- if(mui>(1.+mu-muj)&&
- ( mui2*mui2 + mu2*mu2 + sqr(1.-muj)*sqr(1.-muj) - 2.*( mui2*mu2+mui2*sqr(1.-muj)+mu2*sqr(1.-muj) ))>0
- ) {
- return rootOfKallen( mui2, mu2, sqr(1.-muj) ) / ( 2.-2.*sqrt(muj2) ) * dScale;
- }
- else
- return ZERO;
+ return rootOfKallen( mui2, mu2, sqr(1.-muj) ) / ( 2.-2.*sqrt(muj2) ) * dScale;
}
Energy FFMassiveKinematics::QMax(Energy dScale,
double, double,
const DipoleIndex& ind,
const DipoleSplittingKernel&) const {
assert(false && "implementation missing");
double Muj = ind.spectatorData()->mass() / dScale;
return dScale * ( 1.-2.*Muj+sqr(Muj) );
}
Energy FFMassiveKinematics::PtFromQ(Energy scale, const DipoleSplittingInfo& split) const {
// from Martin's thesis
double z = split.lastZ();
Energy mi = split.emitterData()->mass();
Energy m = split.emissionData()->mass();
Energy2 pt2 = z*(1.-z)*sqr(scale) - (1-z)*sqr(mi) - z*sqr(m);
assert(pt2 >= ZERO);
return sqrt(pt2);
}
Energy FFMassiveKinematics::QFromPt(Energy scale, const DipoleSplittingInfo& split) const {
// from Martin's thesis
double z = split.lastZ();
Energy mi = split.emitterData()->mass();
Energy m = split.emissionData()->mass();
Energy2 Q2 = (sqr(scale) + (1-z)*sqr(mi) + z*sqr(m))/(z*(1.-z));
return sqrt(Q2);
}
double FFMassiveKinematics::ptToRandom(Energy pt, Energy,
double,double,
const DipoleIndex&,
const DipoleSplittingKernel&) const {
return log(pt/IRCutoff()) / log(0.5 * generator()->maximumCMEnergy()/IRCutoff());
}
// own kinematics close to Dinsdale,Ternick,Weinzierl
bool FFMassiveKinematics::generateSplitting(double kappa, double xi, double rphi,
DipoleSplittingInfo& info,
const DipoleSplittingKernel&) {
Energy pt = IRCutoff() * pow(0.5 * generator()->maximumCMEnergy()/IRCutoff(),kappa);
if ( pt > info.hardPt() || pt < IRCutoff() ) {
jacobian(0.0);
return false;
}
double z;
double mapZJacobian;
if ( info.index().emitterData()->id() == ParticleID::g ) {
if ( info.emissionData()->id() != ParticleID::g ) {
z = xi;
mapZJacobian = 1.;
} else {
z = exp(xi)/(1.+exp(xi));
mapZJacobian = z*(1.-z);
}
} else {
z = 1.-exp(-xi);
mapZJacobian = 1.-z;
}
// masses
double mui2 = sqr( info.emitterData()->mass() / info.scale() );
double mu2 = sqr( info.emissionData()->mass() / info.scale() );
double muj2 = sqr( info.spectatorData()->mass() / info.scale() );
double Mui2 = 0.;
if ( info.emitterData()->id() + info.emissionData()->id() == 0 ) Mui2 = 0.; // gluon
else Mui2 = mui2; // (anti)quark
double Muj2 = muj2;
// new: 2011-08-31
// 2011-11-08: this does happen
if( sqrt(mui2)+sqrt(mu2)+sqrt(muj2) > 1. ){
jacobian(0.0);
return false;
}
double bar = 1.-mui2-mu2-muj2;
double y = ( sqr( pt / info.scale() ) + sqr(1.-z)*mui2 + z*z*mu2 ) /
(z*(1.-z)*bar);
// do this here for simplicity
Energy ptmax1 = rootOfKallen( mui2, mu2, sqr(1.-sqrt(muj2)) ) /
( 2.-2.*sqrt(muj2) ) * info.scale();
Energy auxHardPt = ptmax1 > info.hardPt() ? info.hardPt() : ptmax1;
// 2011-11-09
assert(ptmax1>=info.hardPt());
// phasespace constraint to incorporate ptMax
double zp1 = ( 1.+mui2-mu2+muj2-2.*sqrt(muj2) +
rootOfKallen(mui2,mu2,sqr(1-sqrt(muj2))) *
sqrt( 1.-sqr(pt/auxHardPt) ) ) /
( 2.*sqr(1.-sqrt(muj2)) );
double zm1 = ( 1.+mui2-mu2+muj2-2.*sqrt(muj2) -
rootOfKallen(mui2,mu2,sqr(1-sqrt(muj2))) *
sqrt( 1.-sqr(pt/auxHardPt) ) ) /
( 2.*sqr(1.-sqrt(muj2)) );
if ( z > zp1 || z < zm1 ||
pt > auxHardPt) {
jacobian(0.0);
return false;
}
// kinematic phasespace boundaries for (y,z)
// same as in Dittmaier hep-ph/9904440v2 (equivalent to CS)
double ym = 2.*sqrt(mui2)*sqrt(mu2)/bar;
double yp = 1. - 2.*sqrt(muj2)*(1.-sqrt(muj2))/bar;
if ( y < ym || y > yp ) {
jacobian(0.0);
return false;
}
double zm = ( (2.*mui2+bar*y)*(1.-y) - sqrt(y*y-ym*ym)*sqrt(sqr(2.*muj2+bar-bar*y)-4.*muj2) ) /
( 2.*(1.-y)*(mui2+mu2+bar*y) );
double zp = ( (2.*mui2+bar*y)*(1.-y) + sqrt(y*y-ym*ym)*sqrt(sqr(2.*muj2+bar-bar*y)-4.*muj2) ) /
( 2.*(1.-y)*(mui2+mu2+bar*y) );
if ( z < zm || z > zp ) {
jacobian(0.0);
return false;
}
double phi = 2.*Constants::pi*rphi;
jacobian( 2. * mapZJacobian * (1.-y) *
log(0.5 * generator()->maximumCMEnergy()/IRCutoff()) *
bar / rootOfKallen(1.,Mui2,Muj2) );
lastPt(pt);
lastZ(z);
lastPhi(phi);
if ( theMCCheck )
theMCCheck->book(1.,1.,info.scale(),info.hardPt(),pt,z,jacobian());
return true;
}
// kinematics close to Dinsdale,Ternick,Weinzierl
// revised 2011-08-22
// revised 2011-11-06
void FFMassiveKinematics::generateKinematics(const Lorentz5Momentum& pEmitter,
const Lorentz5Momentum& pSpectator,
const DipoleSplittingInfo& dInfo) {
double z = dInfo.lastZ();
Energy pt = dInfo.lastPt();
// masses
double mui2 = sqr( dInfo.emitterData()->mass() / dInfo.scale() );
double mu2 = sqr( dInfo.emissionData()->mass() / dInfo.scale() );
double muj2 = sqr( dInfo.spectatorData()->mass() / dInfo.scale() );
double y = ( sqr( pt / dInfo.scale() ) + sqr(1.-z)*mui2 + z*z*mu2 ) /
(z*(1.-z)*(1.-mui2-mu2-muj2));
Energy2 sbar = sqr(dInfo.scale()) *(1.-mui2-mu2-muj2);
// CMF: particle energies
Energy Ei = ( sbar*(1.-(1.-z)*(1.-y)) + 2.*sqr(dInfo.scale())*mui2 ) / (2.*dInfo.scale());
Energy E = ( sbar*(1.- z *(1.-y)) + 2.*sqr(dInfo.scale())*mu2 ) / (2.*dInfo.scale());
Energy Ej = ( sbar*(1.- y ) + 2.*sqr(dInfo.scale())*muj2 ) / (2.*dInfo.scale());
// CMF: momenta in z-direction (axis of pEmitter & pSpectator)
Energy qi3 = (2.*Ei*Ej-z*(1.-y)*sbar ) / 2./sqrt(Ej*Ej-sqr(dInfo.scale())*muj2);
Energy q3 = (2.*E *Ej-(1.-z)*(1.-y)*sbar) / 2./sqrt(Ej*Ej-sqr(dInfo.scale())*muj2);
Energy qj3 = sqrt( sqr(Ej) - sqr(dInfo.scale())*muj2 );
// get z axis in the dipole's CMF which is parallel to pSpectator
Boost toCMF = (pEmitter+pSpectator).findBoostToCM();
Lorentz5Momentum pjAux = pSpectator; pjAux.boost(toCMF);
ThreeVector<double> pjAxis = pjAux.vect().unit();
// set the momenta in this special reference frame
// note that pt might in some cases differ from the physical pt!
// phi is defined exactly as in getKt
Energy ptResc = sqrt( sqr(Ei)-sqr(dInfo.scale())*mui2-sqr(qi3) );
Lorentz5Momentum em ( ptResc*cos(dInfo.lastPhi()), -ptResc*sin(dInfo.lastPhi()), qi3, Ei );
Lorentz5Momentum emm ( -ptResc*cos(dInfo.lastPhi()), ptResc*sin(dInfo.lastPhi()), q3, E );
Lorentz5Momentum spe ( 0.*GeV, 0.*GeV, qj3, Ej );
// output the mismatch between pt and physical pt
// ofstream output1("ptDiffOnPtAxis-uub-m.dat",ofstream::app);
// ofstream output2("ptDiffOnCosAxis-uub-m.dat",ofstream::app);
// if( abs(dInfo.spectatorData()->id())==5 && dInfo.emitterData()->id()+dInfo.emissionData()->id()==0 &&
// abs(dInfo.emitterData()->id())==1 ) {
// output1 << pt/dInfo.scale() << " " << abs(ptResc-pt)/(ptResc+pt) << " " << endl;
// output2 << em.vect().unit()*emm.vect().unit() << " " << abs(ptResc-pt)/(ptResc+pt) << " " << endl;
// }
// output1.close(); output2.close();
// rotate back
em.rotateUz (pjAxis);
emm.rotateUz(pjAxis);
spe.rotateUz(pjAxis);
// boost back
em.boost (-toCMF);
emm.boost(-toCMF);
spe.boost(-toCMF);
// mass shells, rescale energy
em.setMass(dInfo.scale()*sqrt(mui2));
em.rescaleEnergy();
emm.setMass(dInfo.scale()*sqrt(mu2));
emm.rescaleEnergy();
spe.setMass(dInfo.scale()*sqrt(muj2));
spe.rescaleEnergy();
// book
emitterMomentum(em);
emissionMomentum(emm);
spectatorMomentum(spe);
// TODO: remove
// 2011-11-09: never occurred
if(em.t()/GeV>=0. && emm.t()/GeV>=0. && spe.t()/GeV>=0.);
else cout << "FFMassiveKinematics::generateKinematics momenta corrupt" << endl;
// 2011-11-03 LEP run with full masses:
// x,y,t no problem
// z order > 5.e-7 happend 41 times in 10000 runs
// maximum was 2.5e-6
// double order=2.4e-6;
// Lorentz5Momentum pDif=em+emm+spe-pEmitter-pSpectator, pSum=em+emm+spe+pEmitter+pSpectator;
// if(abs(pDif.x()/(pSum.x()==ZERO ? 1.*GeV : pSum.x())) <order || (pDif-pSum).x()==ZERO);
// else cout << "FFMassiveKinematics::generateKinematics momenta corrupt: x " <<
// abs(pDif.x()/(pSum.x()==ZERO ? 1.*GeV : pSum.x())) << endl <<
// " " << em/GeV << " " << emm/GeV << " " << spe/GeV << endl <<
// " " << pEmitter/GeV << " " << pSpectator/GeV << endl;
// if(abs(pDif.y()/(pSum.y()==ZERO ? 1.*GeV : pSum.y())) <order || (pDif-pSum).y()==ZERO);
// else cout << "FFMassiveKinematics::generateKinematics momenta corrupt: y " <<
// abs(pDif.y()/(pSum.y()==ZERO ? 1.*GeV : pSum.y())) << endl;
// if(abs(pDif.z()/(pSum.z()==ZERO ? 1.*GeV : pSum.z())) <order | (pDif-pSum).z()==ZERO);
// else cout << "FFMassiveKinematics::generateKinematics momenta corrupt: z " <<
// abs(pDif.z()/(pSum.z()==ZERO ? 1.*GeV : pSum.z())) << endl <<
// " " << em/GeV << " " << emm/GeV << " " << spe/GeV << endl <<
// " " << pEmitter/GeV << " " << pSpectator/GeV << endl;
// if(abs(pDif.t()/(pSum.t()==ZERO ? 1.*GeV : pSum.t())) <order || (pDif-pSum).t()==ZERO);
// else cout << "FFMassiveKinematics::generateKinematics momenta corrupt: t " <<
// // abs(pDif.t()/(pSum.t()==ZERO ? 1.*GeV : pSum.t())) << endl;
// cout << endl;
}
// If needed, insert default implementations of function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void FFMassiveKinematics::persistentOutput(PersistentOStream & ) const {
}
void FFMassiveKinematics::persistentInput(PersistentIStream & , int) {
}
ClassDescription<FFMassiveKinematics> FFMassiveKinematics::initFFMassiveKinematics;
// Definition of the static class description member.
void FFMassiveKinematics::Init() {
static ClassDocumentation<FFMassiveKinematics> documentation
("FFMassiveKinematics implements massive splittings "
"off a final-final dipole.");
}
diff --git a/DipoleShower/Kinematics/FFMassiveKinematics.h b/DipoleShower/Kinematics/FFMassiveKinematics.h
--- a/DipoleShower/Kinematics/FFMassiveKinematics.h
+++ b/DipoleShower/Kinematics/FFMassiveKinematics.h
@@ -1,261 +1,262 @@
// -*- C++ -*-
//
// FFMassiveKinematics.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_FFMassiveKinematics_H
#define HERWIG_FFMassiveKinematics_H
//
// This is the declaration of the FFMassiveKinematics class.
//
#include "DipoleSplittingKinematics.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup DipoleShower
* \author Simon Platzer
*
* \brief FFMassiveKinematics implements massive splittings
* off a final-final dipole.
*
*/
class FFMassiveKinematics: public DipoleSplittingKinematics {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
FFMassiveKinematics();
/**
* The destructor.
*/
virtual ~FFMassiveKinematics();
//@}
public:
/**
* Return the boundaries in between the evolution
* variable random number is to be sampled; the lower
* cuoff is assumed to correspond to the infrared cutoff.
*/
virtual pair<double,double> kappaSupport(const DipoleSplittingInfo& dIndex) const;
/**
* Return the boundaries in between the momentum
* fraction random number is to be sampled.
*/
virtual pair<double,double> xiSupport(const DipoleSplittingInfo& dIndex) const;
/**
* Return the boundaries on the momentum fraction
*/
virtual pair<double,double> zBoundaries(Energy,
const DipoleSplittingInfo&,
const DipoleSplittingKernel&) const {
return make_pair(0.0,1.0);
}
/**
* Return the dipole scale associated to the
* given pair of emitter and spectator. This
* should be the invariant mass or absolute value
* final/final or initial/initial and the absolute
* value of the momentum transfer for intial/final or
* final/initial dipoles.
*/
virtual Energy dipoleScale(const Lorentz5Momentum& pEmitter,
const Lorentz5Momentum& pSpectator) const;
/**
* Return the maximum pt for the given dipole scale.
*/
virtual Energy ptMax(Energy dScale,
double emX, double specX,
const DipoleIndex& ind,
const DipoleSplittingKernel& split) const;
/**
* Return the maximum virtuality for the given dipole scale.
*/
virtual Energy QMax(Energy dScale,
double emX, double specX,
const DipoleIndex& dIndex,
const DipoleSplittingKernel& split) const;
/**
* Return the pt given a virtuality.
*/
virtual Energy PtFromQ(Energy scale, const DipoleSplittingInfo&) const;
/**
* Return the virtuality given a pt.
*/
virtual Energy QFromPt(Energy scale, const DipoleSplittingInfo&) const;
/**
* Return the random number associated to
* the given pt.
*/
virtual double ptToRandom(Energy pt, Energy dScale,
double emX, double specX,
const DipoleIndex& dIndex,
const DipoleSplittingKernel& split) const;
/**
* Generate splitting variables given three random numbers
* and the momentum fractions of the emitter and spectator.
* Return true on success.
*/
virtual bool generateSplitting(double kappa, double xi, double phi,
DipoleSplittingInfo& dIndex,
const DipoleSplittingKernel& split);
/**
* Generate the full kinematics given emitter and
* spectator momentum and a previously completeted
* DipoleSplittingInfo object.
*/
virtual void generateKinematics(const Lorentz5Momentum& pEmitter,
const Lorentz5Momentum& pSpectator,
const DipoleSplittingInfo& dInfo);
public:
/**
* Triangular / Kallen function
*/
template <class T>
inline double rootOfKallen (T a, T b, T c) const {
- return sqrt( a*a + b*b + c*c - 2.*( a*b+a*c+b*c ) ); }
+ double sres=a*a + b*b + c*c - 2.*( a*b+a*c+b*c );
+ return sres>0.?sqrt( sres ):0.; }
/**
* Perform a rotation on both momenta such that the first one will
* point along the (positive) z axis. Rotate back to the original
* reference frame by applying rotateUz(returnedVector) to each momentum.
*/
ThreeVector<double> rotateToZ (Lorentz5Momentum& pTarget, Lorentz5Momentum& p1){
ThreeVector<double> oldAxis = pTarget.vect().unit();
double ct = oldAxis.z(); double st = sqrt( 1.-sqr(ct) ); // cos,sin(theta)
double cp = oldAxis.x()/st; double sp = oldAxis.y()/st; // cos,sin(phi)
pTarget.setZ( pTarget.vect().mag() ); pTarget.setX( 0.*GeV ); pTarget.setY( 0.*GeV );
Lorentz5Momentum p1old = p1;
p1.setX( sp*p1old.x() - cp*p1old.y() );
p1.setY( ct*cp*p1old.x() + ct*sp*p1old.y() - st*p1old.z() );
p1.setZ( st*cp*p1old.x() + st*sp*p1old.y() + ct*p1old.z() );
return oldAxis;
}
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);
//@}
/**
* 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();
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;
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<FFMassiveKinematics> initFFMassiveKinematics;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
FFMassiveKinematics & operator=(const FFMassiveKinematics &);
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of FFMassiveKinematics. */
template <>
struct BaseClassTrait<Herwig::FFMassiveKinematics,1> {
/** Typedef of the first base class of FFMassiveKinematics. */
typedef Herwig::DipoleSplittingKinematics NthBase;
};
/** This template specialization informs ThePEG about the name of
* the FFMassiveKinematics class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::FFMassiveKinematics>
: public ClassTraitsBase<Herwig::FFMassiveKinematics> {
/** Return a platform-independent class name */
static string className() { return "Herwig::FFMassiveKinematics"; }
/**
* The name of a file containing the dynamic library where the class
* FFMassiveKinematics is implemented. It may also include several, space-separated,
* libraries if the class FFMassiveKinematics depends on other classes (base classes
* excepted). In this case the listed libraries will be dynamically
* linked in the order they are specified.
*/
static string library() { return "HwDipoleShower.so"; }
};
/** @endcond */
}
#endif /* HERWIG_FFMassiveKinematics_H */
diff --git a/DipoleShower/Kinematics/FIMassiveKinematics.cc b/DipoleShower/Kinematics/FIMassiveKinematics.cc
--- a/DipoleShower/Kinematics/FIMassiveKinematics.cc
+++ b/DipoleShower/Kinematics/FIMassiveKinematics.cc
@@ -1,282 +1,279 @@
// -*- C++ -*-
//
// FIMassiveKinematics.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the FIMassiveKinematics class.
//
#include "FIMassiveKinematics.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "Herwig++/DipoleShower/Base/DipoleSplittingInfo.h"
#include "Herwig++/DipoleShower/Kernels/DipoleSplittingKernel.h"
using namespace Herwig;
FIMassiveKinematics::FIMassiveKinematics()
: DipoleSplittingKinematics() {}
FIMassiveKinematics::~FIMassiveKinematics() {}
IBPtr FIMassiveKinematics::clone() const {
return new_ptr(*this);
}
IBPtr FIMassiveKinematics::fullclone() const {
return new_ptr(*this);
}
pair<double,double> FIMassiveKinematics::kappaSupport(const DipoleSplittingInfo&) const {
return make_pair(0.0,1.0);
}
pair<double,double> FIMassiveKinematics::xiSupport(const DipoleSplittingInfo& split) const {
double c = sqrt(1.-4.*sqr(IRCutoff()/generator()->maximumCMEnergy()));
if ( split.index().emitterData()->id() == ParticleID::g ) {
if ( split.emissionData()->id() != ParticleID::g )
return make_pair(0.5*(1.-c),0.5*(1.+c));
double b = log((1.+c)/(1.-c));
return make_pair(-b,b);
}
return make_pair(-log(0.5*(1.+c)),-log(0.5*(1.-c)));
}
// sbar
Energy FIMassiveKinematics::dipoleScale(const Lorentz5Momentum& pEmitter,
const Lorentz5Momentum& pSpectator) const {
return sqrt(2.*(pEmitter*pSpectator));
}
Energy FIMassiveKinematics::ptMax(Energy dScale,
double, double specX,
const DipoleIndex& ind,
const DipoleSplittingKernel& split) const {
Energy mi = split.emitter(ind)->mass(), m = split.emission(ind)->mass();
Energy2 mi2 = sqr(mi), m2 = sqr(m);
// Energy2 Mi2 = split.emitter(int)->id() + split.emission(int)->id() == 0 ?
// 0.*GeV2 : mi2;
Energy2 Mi2 = mi2 == m2 ? 0.*GeV2 : mi2;
// s^star/x
Energy2 s = sqr(dScale) * (1.-specX)/specX + Mi2;
Energy roots = sqrt(s);
- if(roots>mi+m)
- return .5 * sqrt(s) * rootOfKallen( s/s, mi2/s, m2/s );
- else
- return ZERO;
+ return .5 * sqrt(s) * rootOfKallen( s/s, mi2/s, m2/s );
}
// what is this? in FF it is given by y+*dScale = sqrt( 2qi*q / bar )->max
Energy FIMassiveKinematics::QMax(Energy dScale,
double, double specX,
const DipoleIndex&,
const DipoleSplittingKernel&) const {
assert(false && "implementation missing");
cout << "FIMassiveKinematics::QMax called.\n" << flush;
// this is sqrt( 2qi*q ) -> max;
return dScale * sqrt((1.-specX)/specX);
}
Energy FIMassiveKinematics::PtFromQ(Energy scale, const DipoleSplittingInfo& split) const {
// from Martin's thesis
double z = split.lastZ();
Energy mi = split.emitterData()->mass();
Energy m = split.emissionData()->mass();
Energy2 pt2 = z*(1.-z)*sqr(scale) - (1-z)*sqr(mi) - z*sqr(m);
assert(pt2 >= ZERO);
return sqrt(pt2);
}
Energy FIMassiveKinematics::QFromPt(Energy scale, const DipoleSplittingInfo& split) const {
// from Martin's thesis
double z = split.lastZ();
Energy mi = split.emitterData()->mass();
Energy m = split.emissionData()->mass();
Energy2 Q2 = (sqr(scale) + (1-z)*sqr(mi) + z*sqr(m))/(z*(1.-z));
return sqrt(Q2);
}
double FIMassiveKinematics::ptToRandom(Energy pt, Energy,
double,double,
const DipoleIndex&,
const DipoleSplittingKernel&) const {
return log(pt/IRCutoff()) / log(0.5 * generator()->maximumCMEnergy()/IRCutoff());
}
bool FIMassiveKinematics::generateSplitting(double kappa, double xi, double rphi,
DipoleSplittingInfo& info,
const DipoleSplittingKernel&) {
if ( info.spectatorX() < xMin() ) {
jacobian(0.0);
return false;
}
Energy pt = IRCutoff() * pow(0.5 * generator()->maximumCMEnergy()/IRCutoff(),kappa);
if ( pt > info.hardPt() || pt < IRCutoff() ) {
jacobian(0.0);
return false;
}
double z;
double mapZJacobian;
if ( info.index().emitterData()->id() == ParticleID::g ) {
if ( info.emissionData()->id() != ParticleID::g ) {
z = xi;
mapZJacobian = 1.;
} else {
z = exp(xi)/(1.+exp(xi));
mapZJacobian = z*(1.-z);
}
} else {
z = 1.-exp(-xi);
mapZJacobian = 1.-z;
}
// double s = z*(1.-z);
// double xs = info.spectatorX();
// double x = 1. / ( 1. + sqr(pt/info.scale()) / s );
// double zp = 0.5*(1.+sqrt(1.-sqr(pt/info.hardPt())));
// double zm = 0.5*(1.-sqrt(1.-sqr(pt/info.hardPt())));
Energy2 mi2 = sqr(info.emitterData()->mass());
Energy2 m2 = sqr(info.emissionData()->mass());
Energy2 Mi2 = info.emitterData()->id()+info.emissionData()->id() == 0 ?
0.*GeV2 : mi2;
// s^star/x
Energy2 s = sqr(info.scale()) * (1.-info.spectatorX())/info.spectatorX() + Mi2;
double xs = info.spectatorX();
double x = 1. / ( 1. +
( sqr(pt) + (1.-z)*mi2 + z*m2 - z*(1.-z)*Mi2 ) /
( z*(1.-z)*s ) );
double zm1 = .5*( 1.+(mi2-m2)/s - rootOfKallen(s/s,mi2/s,m2/s) *
sqrt( 1.-sqr(pt/info.hardPt()) ) );
double zp1 = .5*( 1.+(mi2-m2)/s + rootOfKallen(s/s,mi2/s,m2/s) *
sqrt( 1.-sqr(pt/info.hardPt()) ) );
if ( // pt < IRCutoff() ||
// pt > info.hardPt() ||
z > zp1 || z < zm1 ||
x < xs ) {
jacobian(0.0);
return false;
}
// additional purely kinematic constraints
double mui2 = x*mi2/sqr(info.scale());
double mu2 = x*m2/sqr(info.scale());
double Mui2 = x*Mi2/sqr(info.scale());
double xp = 1. + Mui2 - sqr(sqrt(mui2)+sqrt(mu2));
double root = sqr(1.-x+Mui2-mui2-mu2)-4.*mui2*mu2;
if( root < 0. && root>-1e-10 )
root = 0.;
else if (root <0. ) {
jacobian(0.0);
return false;
}
root = sqrt(root);
double zm = .5*( 1.-x+Mui2+mui2-mui2 - root ) / (1.-x+Mui2);
double zp = .5*( 1.-x+Mui2+mui2-mui2 + root ) / (1.-x+Mui2);
if (x > xp ||
z > zp || z < zm ) {
jacobian(0.0);
return false;
}
double phi = 2.*Constants::pi*rphi;
jacobian( 2. * mapZJacobian * log(0.5 * generator()->maximumCMEnergy()/IRCutoff()));
lastPt(pt);
lastZ(z);
lastPhi(phi);
lastSpectatorZ(x);
if ( theMCCheck )
theMCCheck->book(1.,info.spectatorX(),info.scale(),info.hardPt(),pt,z,jacobian());
return true;
}
void FIMassiveKinematics::generateKinematics(const Lorentz5Momentum& pEmitter,
const Lorentz5Momentum& pSpectator,
const DipoleSplittingInfo& dInfo) {
Energy pt = dInfo.lastPt();
double z = dInfo.lastZ();
Lorentz5Momentum kt =
getKt (pSpectator, pEmitter, pt, dInfo.lastPhi(),true);
Energy2 mi2 = sqr(dInfo.emitterData()->mass());
Energy2 m2 = sqr(dInfo.emissionData()->mass());
Energy2 Mi2 = dInfo.emitterData()->id() + dInfo.emissionData()->id() == 0 ?
0.*GeV2 : mi2;
double xInv = ( 1. +
(pt*pt+(1.-z)*mi2+z*m2-z*(1.-z)*Mi2) /
(z*(1.-z)*sqr(dInfo.scale())) );
Lorentz5Momentum em = z*pEmitter +
(sqr(pt)+mi2-z*z*Mi2)/(z*sqr(dInfo.scale()))*pSpectator + kt;
em.setMass(sqrt(mi2));
em.rescaleEnergy();
Lorentz5Momentum emm = (1.-z)*pEmitter +
(pt*pt+m2-sqr(1.-z)*Mi2)/((1.-z)*sqr(dInfo.scale()))*pSpectator - kt;
emm.setMass(sqrt(m2));
emm.rescaleEnergy();
Lorentz5Momentum spe = xInv*pSpectator;
spe.setMass(ZERO);
spe.rescaleEnergy();
emitterMomentum(em);
emissionMomentum(emm);
spectatorMomentum(spe);
}
// If needed, insert default implementations of function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void FIMassiveKinematics::persistentOutput(PersistentOStream & ) const {
}
void FIMassiveKinematics::persistentInput(PersistentIStream & , int) {
}
ClassDescription<FIMassiveKinematics> FIMassiveKinematics::initFIMassiveKinematics;
// Definition of the static class description member.
void FIMassiveKinematics::Init() {
static ClassDocumentation<FIMassiveKinematics> documentation
("FIMassiveKinematics implements massless splittings "
"off a final-initial dipole.");
}
diff --git a/DipoleShower/Kinematics/FIMassiveKinematics.h b/DipoleShower/Kinematics/FIMassiveKinematics.h
--- a/DipoleShower/Kinematics/FIMassiveKinematics.h
+++ b/DipoleShower/Kinematics/FIMassiveKinematics.h
@@ -1,244 +1,246 @@
// -*- C++ -*-
//
// FILightKinematics.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_FILightKinematics_H
#define HERWIG_FILightKinematics_H
//
// This is the declaration of the FILightKinematics class.
//
#include "DipoleSplittingKinematics.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup DipoleShower
* \author Simon Platzer, Martin Stoll
*
* \brief FIMassiveKinematics implements massless splittings
* off a final-initial dipole.
*
*/
class FIMassiveKinematics: public DipoleSplittingKinematics {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
FIMassiveKinematics();
/**
* The destructor.
*/
virtual ~FIMassiveKinematics();
//@}
public:
/**
* Return the boundaries in between the evolution
* variable random number is to be sampled; the lower
* cuoff is assumed to correspond to the infrared cutoff.
*/
virtual pair<double,double> kappaSupport(const DipoleSplittingInfo& dIndex) const;
/**
* Return the boundaries in between the momentum
* fraction random number is to be sampled.
*/
virtual pair<double,double> xiSupport(const DipoleSplittingInfo& dIndex) const;
/**
* Return the boundaries on the momentum fraction
*/
virtual pair<double,double> zBoundaries(Energy,
const DipoleSplittingInfo&,
const DipoleSplittingKernel&) const {
return make_pair(0.0,1.0);
}
/**
* Return the dipole scale associated to the
* given pair of emitter and spectator. This
* should be the invariant mass or absolute value
* final/final or initial/initial and the absolute
* value of the momentum transfer for intial/final or
* final/initial dipoles.
*/
virtual Energy dipoleScale(const Lorentz5Momentum& pEmitter,
const Lorentz5Momentum& pSpectator) const;
/**
* Return the maximum pt for the given dipole scale.
*/
virtual Energy ptMax(Energy dScale,
double emX, double specX,
const DipoleIndex& dIndex,
const DipoleSplittingKernel&) const;
/**
* Return the maximum virtuality for the given dipole scale.
*/
virtual Energy QMax(Energy dScale,
double emX, double specX,
const DipoleIndex& dIndex,
const DipoleSplittingKernel&) const;
/**
* Return the pt given a virtuality.
*/
virtual Energy PtFromQ(Energy scale, const DipoleSplittingInfo&) const;
/**
* Return the virtuality given a pt.
*/
virtual Energy QFromPt(Energy scale, const DipoleSplittingInfo&) const;
/**
* Return the random number associated to
* the given pt.
*/
virtual double ptToRandom(Energy pt, Energy dScale,
double emX, double specX,
const DipoleIndex& dIndex,
const DipoleSplittingKernel&) const;
/**
* Generate splitting variables given three random numbers
* and the momentum fractions of the emitter and spectator.
* Return true on success.
*/
virtual bool generateSplitting(double kappa, double xi, double phi,
DipoleSplittingInfo& dIndex,
const DipoleSplittingKernel&);
/**
* Generate the full kinematics given emitter and
* spectator momentum and a previously completeted
* DipoleSplittingInfo object.
*/
virtual void generateKinematics(const Lorentz5Momentum& pEmitter,
const Lorentz5Momentum& pSpectator,
const DipoleSplittingInfo& dInfo);
public:
/**
* Triangular / Kallen function
*/
template <class T>
inline double rootOfKallen (T a, T b, T c) const {
- return sqrt( a*a + b*b + c*c - 2.*( a*b+a*c+b*c ) ); }
+ double sres=a*a + b*b + c*c - 2.*( a*b+a*c+b*c );
+ return sres>0.?sqrt( sres ):0.; }
+
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);
//@}
/**
* 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();
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;
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<FIMassiveKinematics> initFIMassiveKinematics;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
FIMassiveKinematics & operator=(const FIMassiveKinematics &);
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of FIMassiveKinematics. */
template <>
struct BaseClassTrait<Herwig::FIMassiveKinematics,1> {
/** Typedef of the first base class of FIMassiveKinematics. */
typedef Herwig::DipoleSplittingKinematics NthBase;
};
/** This template specialization informs ThePEG about the name of
* the FIMassiveKinematics class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::FIMassiveKinematics>
: public ClassTraitsBase<Herwig::FIMassiveKinematics> {
/** Return a platform-independent class name */
static string className() { return "Herwig::FIMassiveKinematics"; }
/**
* The name of a file containing the dynamic library where the class
* FIMassiveKinematics is implemented. It may also include several, space-separated,
* libraries if the class FIMassiveKinematics depends on other classes (base classes
* excepted). In this case the listed libraries will be dynamically
* linked in the order they are specified.
*/
static string library() { return "HwDipoleShower.so"; }
};
/** @endcond */
}
#endif /* HERWIG_FIMassiveKinematics_H */

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 6:39 PM (1 d, 16 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805641
Default Alt Text
(36 KB)

Event Timeline