Page MenuHomeHEPForge
Files F8310253

No OneTemporary
Actions

View Options

diff --git a/Shower/Base/ShowerParticle.cc b/Shower/Base/ShowerParticle.cc
--- a/Shower/Base/ShowerParticle.cc
+++ b/Shower/Base/ShowerParticle.cc
@@ -1,45 +1,431 @@
// -*- C++ -*-
//
// ShowerParticle.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 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 ShowerParticle class.
//
#include "ShowerParticle.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
+#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
+#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
+#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
+#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
+#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
#include <ostream>
using namespace Herwig;
PPtr ShowerParticle::clone() const {
return new_ptr(*this);
}
PPtr ShowerParticle::fullclone() const {
return new_ptr(*this);
}
ClassDescription<ShowerParticle> ShowerParticle::initShowerParticle;
// Definition of the static class description member.
void ShowerParticle::vetoEmission(ShowerPartnerType::Type, Energy scale) {
scales_.QED = min(scale,scales_.QED );
scales_.QED_noAO = min(scale,scales_.QED_noAO );
scales_.QCD_c = min(scale,scales_.QCD_c );
scales_.QCD_c_noAO = min(scale,scales_.QCD_c_noAO );
scales_.QCD_ac = min(scale,scales_.QCD_ac );
scales_.QCD_ac_noAO = min(scale,scales_.QCD_ac_noAO);
scales_.EW = min(scale,scales_.EW );
}
void ShowerParticle::addPartner(EvolutionPartner in ) {
partners_.push_back(in);
}
+
+namespace {
+
+LorentzRotation boostToShower(const vector<Lorentz5Momentum> & basis,
+ ShowerKinematics::Frame frame,
+ Lorentz5Momentum & porig) {
+ LorentzRotation output;
+ if(frame==ShowerKinematics::BackToBack) {
+ // we are doing the evolution in the back-to-back frame
+ // work out the boostvector
+ Boost boostv(-(basis[0]+basis[1]).boostVector());
+ // momentum of the parton
+ Lorentz5Momentum ptest(basis[0]);
+ // construct the Lorentz boost
+ output = LorentzRotation(boostv);
+ ptest *= output;
+ Axis axis(ptest.vect().unit());
+ // now rotate so along the z axis as needed for the splitting functions
+ if(axis.perp2()>1e-10) {
+ double sinth(sqrt(1.-sqr(axis.z())));
+ output.rotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
+ }
+ else if(axis.z()<0.) {
+ output.rotate(Constants::pi,Axis(1.,0.,0.));
+ }
+ porig = output*basis[0];
+ porig.setX(ZERO);
+ porig.setY(ZERO);
+ }
+ else {
+ output = LorentzRotation(-basis[0].boostVector());
+ porig = output*basis[0];
+ porig.setX(ZERO);
+ porig.setY(ZERO);
+ porig.setZ(ZERO);
+ }
+ return output;
+}
+
+RhoDMatrix bosonMapping(ShowerParticle & particle,
+ const Lorentz5Momentum & porig,
+ VectorSpinPtr vspin,
+ const LorentzRotation & rot) {
+ // rotate the original basis
+ vector<LorentzPolarizationVector> sbasis;
+ for(unsigned int ix=0;ix<3;++ix) {
+ sbasis.push_back(vspin->getProductionBasisState(ix));
+ sbasis.back().transform(rot);
+ }
+ // splitting basis
+ vector<LorentzPolarizationVector> fbasis;
+ bool massless(particle.id()==ParticleID::g||particle.id()==ParticleID::gamma);
+ VectorWaveFunction wave(porig,particle.dataPtr(),outgoing);
+ for(unsigned int ix=0;ix<3;++ix) {
+ if(massless&&ix==1) {
+ fbasis.push_back(LorentzPolarizationVector());
+ }
+ else {
+ wave.reset(ix);
+ fbasis.push_back(wave.wave());
+ }
+ }
+ // work out the mapping
+ RhoDMatrix mapping=RhoDMatrix(PDT::Spin1,false);
+ for(unsigned int ix=0;ix<3;++ix) {
+ for(unsigned int iy=0;iy<3;++iy) {
+ mapping(ix,iy)= sbasis[iy].dot(fbasis[ix].conjugate());
+ if(particle.id()<0)
+ mapping(ix,iy)=conj(mapping(ix,iy));
+ }
+ }
+ // \todo need to fix this
+ mapping = RhoDMatrix(PDT::Spin1,false);
+ if(massless) {
+ mapping(0,0) = 1.;
+ mapping(2,2) = 1.;
+ }
+ else {
+ mapping(0,0) = 1.;
+ mapping(1,1) = 1.;
+ mapping(2,2) = 1.;
+ }
+ return mapping;
+}
+
+RhoDMatrix fermionMapping(ShowerParticle & particle,
+ const Lorentz5Momentum & porig,
+ FermionSpinPtr fspin,
+ const LorentzRotation & rot) {
+ // extract the original basis states
+ vector<LorentzSpinor<SqrtEnergy> > sbasis;
+ for(unsigned int ix=0;ix<2;++ix) {
+ sbasis.push_back(fspin->getProductionBasisState(ix));
+ sbasis.back().transform(rot);
+ }
+ // calculate the states in the splitting basis
+ vector<LorentzSpinor<SqrtEnergy> > fbasis;
+ SpinorWaveFunction wave(porig,particle.dataPtr(),
+ particle.id()>0 ? incoming : outgoing);
+ for(unsigned int ix=0;ix<2;++ix) {
+ wave.reset(ix);
+ fbasis.push_back(wave.dimensionedWave());
+ }
+ RhoDMatrix mapping=RhoDMatrix(PDT::Spin1Half,false);
+ for(unsigned int ix=0;ix<2;++ix) {
+ if(fbasis[0].s2()==complex<SqrtEnergy>()) {
+ mapping(ix,0) = sbasis[ix].s3()/fbasis[0].s3();
+ mapping(ix,1) = sbasis[ix].s2()/fbasis[1].s2();
+ }
+ else {
+ mapping(ix,0) = sbasis[ix].s2()/fbasis[0].s2();
+ mapping(ix,1) = sbasis[ix].s3()/fbasis[1].s3();
+ }
+ }
+ return mapping;
+}
+
+FermionSpinPtr createFermionSpinInfo(ShowerParticle & particle,
+ const Lorentz5Momentum & porig,
+ const LorentzRotation & rot,
+ Helicity::Direction dir) {
+ // calculate the splitting basis for the branching
+ // and rotate back to construct the basis states
+ LorentzRotation rinv = rot.inverse();
+ SpinorWaveFunction wave;
+ if(particle.id()>0)
+ wave=SpinorWaveFunction(porig,particle.dataPtr(),incoming);
+ else
+ wave=SpinorWaveFunction(porig,particle.dataPtr(),outgoing);
+ FermionSpinPtr fspin = new_ptr(FermionSpinInfo(particle.momentum(),dir==outgoing));
+ for(unsigned int ix=0;ix<2;++ix) {
+ wave.reset(ix);
+ LorentzSpinor<SqrtEnergy> basis = wave.dimensionedWave();
+ basis.transform(rinv);
+ fspin->setBasisState(ix,basis);
+ fspin->setDecayState(ix,basis);
+ }
+ particle.spinInfo(fspin);
+ return fspin;
+}
+
+VectorSpinPtr createVectorSpinInfo(ShowerParticle & particle,
+ const Lorentz5Momentum & porig,
+ const LorentzRotation & rot,
+ Helicity::Direction dir) {
+ // calculate the splitting basis for the branching
+ // and rotate back to construct the basis states
+ LorentzRotation rinv = rot.inverse();
+ bool massless(particle.id()==ParticleID::g||particle.id()==ParticleID::gamma);
+ VectorWaveFunction wave(porig,particle.dataPtr(),dir);
+ VectorSpinPtr vspin = new_ptr(VectorSpinInfo(particle.momentum(),dir==outgoing));
+ for(unsigned int ix=0;ix<3;++ix) {
+ LorentzPolarizationVector basis;
+ if(massless&&ix==1) {
+ basis = LorentzPolarizationVector();
+ }
+ else {
+ wave.reset(ix);
+ basis = wave.wave();
+ }
+ basis *= rinv;
+ vspin->setBasisState(ix,basis);
+ vspin->setDecayState(ix,basis);
+ }
+ particle.spinInfo(vspin);
+ vspin-> DMatrix() = RhoDMatrix(PDT::Spin1);
+ vspin->rhoMatrix() = RhoDMatrix(PDT::Spin1);
+ if(massless) {
+ vspin-> DMatrix()(0,0) = 0.5;
+ vspin->rhoMatrix()(0,0) = 0.5;
+ vspin-> DMatrix()(2,2) = 0.5;
+ vspin->rhoMatrix()(2,2) = 0.5;
+ }
+ return vspin;
+}
+}
+
+RhoDMatrix ShowerParticle::extractRhoMatrix(ShoKinPtr kinematics,bool forward) {
+ // get the spin density matrix and the mapping
+ RhoDMatrix mapping;
+ SpinPtr inspin;
+ bool needMapping = getMapping(inspin,mapping,kinematics);
+ // set the decayed flag
+ inspin->decay();
+ // get the spin density matrix
+ RhoDMatrix rho = forward ? inspin->rhoMatrix() : inspin->DMatrix();
+ // map to the shower basis if needed
+ if(needMapping) {
+ RhoDMatrix rhop(rho.iSpin(),false);
+ for(int ixa=0;ixa<rho.iSpin();++ixa) {
+ for(int ixb=0;ixb<rho.iSpin();++ixb) {
+ for(int iya=0;iya<rho.iSpin();++iya) {
+ for(int iyb=0;iyb<rho.iSpin();++iyb) {
+ rhop(ixa,ixb) += rho(iya,iyb)*mapping(iya,ixa)*conj(mapping(iyb,ixb));
+ }
+ }
+ }
+ }
+ rhop.normalize();
+ rho = rhop;
+ }
+ return rho;
+}
+
+bool ShowerParticle::getMapping(SpinPtr & output, RhoDMatrix & mapping,
+ ShoKinPtr showerkin) {
+ // if the particle is not from the hard process
+ if(!this->perturbative()) {
+ // mapping is the identity
+ output=this->spinInfo();
+ mapping=RhoDMatrix(this->dataPtr()->iSpin());
+ if(output) {
+ return false;
+ }
+ else {
+ Lorentz5Momentum porig;
+ LorentzRotation rot = boostToShower(showerkin->getBasis(),showerkin->frame(),porig);
+ Helicity::Direction dir = this->isFinalState() ? outgoing : incoming;
+ if(this->dataPtr()->iSpin()==PDT::Spin0) {
+ assert(false);
+ }
+ else if(this->dataPtr()->iSpin()==PDT::Spin1Half) {
+ output = createFermionSpinInfo(*this,porig,rot,dir);
+ }
+ else if(this->dataPtr()->iSpin()==PDT::Spin1) {
+ output = createVectorSpinInfo(*this,porig,rot,dir);
+ }
+ else {
+ assert(false);
+ }
+ return false;
+ }
+ }
+ // if particle is final-state and is from the hard process
+ else if(this->isFinalState()) {
+ assert(this->perturbative()==1 || this->perturbative()==2);
+ // get transform to shower frame
+ Lorentz5Momentum porig;
+ LorentzRotation rot = boostToShower(showerkin->getBasis(),showerkin->frame(),porig);
+ // the rest depends on the spin of the particle
+ PDT::Spin spin(this->dataPtr()->iSpin());
+ mapping=RhoDMatrix(spin,false);
+ // do the spin dependent bit
+ if(spin==PDT::Spin0) {
+ ScalarSpinPtr sspin=dynamic_ptr_cast<ScalarSpinPtr>(this->spinInfo());
+ if(!sspin) {
+ ScalarWaveFunction::constructSpinInfo(this,outgoing,true);
+ }
+ output=this->spinInfo();
+ return false;
+ }
+ else if(spin==PDT::Spin1Half) {
+ FermionSpinPtr fspin=dynamic_ptr_cast<FermionSpinPtr>(this->spinInfo());
+ // spin info exists get information from it
+ if(fspin) {
+ output=fspin;
+ mapping = fermionMapping(*this,porig,fspin,rot);
+ return true;
+ }
+ // spin info does not exist create it
+ else {
+ output = createFermionSpinInfo(*this,porig,rot,outgoing);
+ return false;
+ }
+ }
+ else if(spin==PDT::Spin1) {
+ VectorSpinPtr vspin=dynamic_ptr_cast<VectorSpinPtr>(this->spinInfo());
+ // spin info exists get information from it
+ if(vspin) {
+ output=vspin;
+ mapping = bosonMapping(*this,porig,vspin,rot);
+ return true;
+ }
+ else {
+ output = createVectorSpinInfo(*this,porig,rot,outgoing);
+ return false;
+ }
+ }
+ // not scalar/fermion/vector
+ else
+ assert(false);
+ }
+ // incoming to hard process
+ else if(this->perturbative()==1 && !this->isFinalState()) {
+ // get the basis vectors
+ // get transform to shower frame
+ Lorentz5Momentum porig;
+ LorentzRotation rot = boostToShower(showerkin->getBasis(),showerkin->frame(),porig);
+ porig *= this->x();
+ // the rest depends on the spin of the particle
+ PDT::Spin spin(this->dataPtr()->iSpin());
+ mapping=RhoDMatrix(spin);
+ // do the spin dependent bit
+ if(spin==PDT::Spin0) {
+ cerr << "testing spin 0 not yet implemented " << endl;
+ assert(false);
+ }
+ // spin-1/2
+ else if(spin==PDT::Spin1Half) {
+ FermionSpinPtr fspin=dynamic_ptr_cast<FermionSpinPtr>(this->spinInfo());
+ // spin info exists get information from it
+ if(fspin) {
+ output=fspin;
+ mapping = fermionMapping(*this,porig,fspin,rot);
+ return true;
+ }
+ // spin info does not exist create it
+ else {
+ output = createFermionSpinInfo(*this,porig,rot,incoming);
+ return false;
+ }
+ }
+ // spin-1
+ else if(spin==PDT::Spin1) {
+ VectorSpinPtr vspin=dynamic_ptr_cast<VectorSpinPtr>(this->spinInfo());
+ // spinInfo exists map it
+ if(vspin) {
+ output=vspin;
+ mapping = bosonMapping(*this,porig,vspin,rot);
+ return true;
+ }
+ // create the spininfo
+ else {
+ output = createVectorSpinInfo(*this,porig,rot,incoming);
+ return false;
+ }
+ }
+ assert(false);
+ }
+ // incoming to decay
+ else if(this->perturbative() == 2 && !this->isFinalState()) {
+ // get the basis vectors
+ Lorentz5Momentum porig;
+ LorentzRotation rot=boostToShower(showerkin->getBasis(),
+ showerkin->frame(),porig);
+ // the rest depends on the spin of the particle
+ PDT::Spin spin(this->dataPtr()->iSpin());
+ mapping=RhoDMatrix(spin);
+ // do the spin dependent bit
+ if(spin==PDT::Spin0) {
+ cerr << "testing spin 0 not yet implemented " << endl;
+ assert(false);
+ }
+ // spin-1/2
+ else if(spin==PDT::Spin1Half) {
+ // FermionSpinPtr fspin=dynamic_ptr_cast<FermionSpinPtr>(this->spinInfo());
+ // // spin info exists get information from it
+ // if(fspin) {
+ // output=fspin;
+ // mapping = fermionMapping(*this,porig,fspin,rot);
+ // return true;
+ // // spin info does not exist create it
+ // else {
+ // output = createFermionSpinInfo(*this,porig,rot,incoming);
+ // return false;
+ // }
+ // }
+ assert(false);
+ }
+ // // spin-1
+ // else if(spin==PDT::Spin1) {
+ // VectorSpinPtr vspin=dynamic_ptr_cast<VectorSpinPtr>(this->spinInfo());
+ // // spinInfo exists map it
+ // if(vspin) {
+ // output=vspin;
+ // mapping = bosonMapping(*this,porig,vspin,rot);
+ // return true;
+ // }
+ // // create the spininfo
+ // else {
+ // output = createVectorSpinInfo(*this,porig,rot,incoming);
+ // return false;
+ // }
+ // }
+ // assert(false);
+ assert(false);
+ }
+ else
+ assert(false);
+ return true;
+}
diff --git a/Shower/Base/ShowerParticle.h b/Shower/Base/ShowerParticle.h
--- a/Shower/Base/ShowerParticle.h
+++ b/Shower/Base/ShowerParticle.h
@@ -1,462 +1,481 @@
// -*- C++ -*-
//
// ShowerParticle.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 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_ShowerParticle_H
#define HERWIG_ShowerParticle_H
//
// This is the declaration of the ShowerParticle class.
//
#include "ThePEG/EventRecord/Particle.h"
#include "Herwig/Shower/SplittingFunctions/SplittingFunction.fh"
#include "Herwig/Shower/ShowerConfig.h"
#include "ShowerKinematics.h"
#include "ShowerParticle.fh"
#include <iosfwd>
namespace Herwig {
using namespace ThePEG;
/** \ingroup Shower
* This class represents a particle in the showering process.
* It inherits from the Particle class of ThePEG and has some
* specifics information useful only during the showering process.
*
* Notice that:
* - for forward evolution, it is clear what is meant by parent/child;
* for backward evolution, however, it depends whether we want
* to keep a physical picture or a Monte-Carlo effective one.
* In the former case, an incoming particle (emitting particle)
* splits into an emitted particle and the emitting particle after
* the emission: the latter two are then children of the
* emitting particle, the parent. In the Monte-Carlo effective
* picture, we have that the particle close to the hard subprocess,
* with higher (space-like) virtuality, splits into an emitted particle
* and the emitting particle at lower virtuality: the latter two are,
* in this case, the children of the first one, the parent. However we
* choose a more physical picture where the new emitting particle is the
* parented of the emitted final-state particle and the original emitting
* particle.
* - the pointer to a SplitFun object is set only in the case
* that the particle has undergone a shower emission. This is similar to
* the case of the decay of a normal Particle where
* the pointer to a Decayer object is set only in the case
* that the particle has undergone to a decay.
* In the case of particle connected directly to the hard subprocess,
* there is no pointer to the hard subprocess, but there is a method
* isFromHardSubprocess() which returns true only in this case.
*
* @see Particle
* @see ShowerConfig
* @see ShowerKinematics
*/
class ShowerParticle: public Particle {
public:
/**
* Struct for all the info on an evolution partner
*/
struct EvolutionPartner {
/**
* Constructor
*/
EvolutionPartner(tShowerParticlePtr p,double w, ShowerPartnerType::Type t,
Energy s) : partner(p), weight(w), type(t), scale(s)
{}
/**
* The partner
*/
tShowerParticlePtr partner;
/**
* Weight
*/
double weight;
/**
* Type
*/
ShowerPartnerType::Type type;
/**
* The assoicated evolution scale
*/
Energy scale;
};
/**
* Struct to store the evolution scales
*/
struct EvolutionScales {
/**
* Constructor
*/
EvolutionScales() : QED(),QCD_c(),QCD_ac(),
QED_noAO(),QCD_c_noAO(),QCD_ac_noAO(),EW(),
Max_Q2(Constants::MaxEnergy2)
{}
/**
* QED scale
*/
Energy QED;
/**
* QCD colour scale
*/
Energy QCD_c;
/**
* QCD anticolour scale
*/
Energy QCD_ac;
/**
* QED scale
*/
Energy QED_noAO;
/**
* QCD colour scale
*/
Energy QCD_c_noAO;
/**
* QCD anticolour scale
*/
Energy QCD_ac_noAO;
/**
* EW scale
*/
Energy EW;
/**
* Maximum allowed virtuality of the particle
*/
Energy2 Max_Q2;
};
/** @name Construction and descruction functions. */
//@{
/**
* Standard Constructor. Note that the default constructor is
* private - there is no particle without a pointer to a
* ParticleData object.
* @param x the ParticleData object
* @param fs Whether or not the particle is an inital or final-state particle
* @param tls Whether or not the particle initiates a time-like shower
*/
ShowerParticle(tcEventPDPtr x, bool fs, bool tls=false)
: Particle(x), _isFinalState(fs),
_perturbative(0), _initiatesTLS(tls), _x(1.0), _showerKinematics(),
_vMass(ZERO), _thePEGBase() {}
/**
* Copy constructor from a ThePEG Particle
* @param x ThePEG particle
* @param pert Where the particle came from
* @param fs Whether or not the particle is an inital or final-state particle
* @param tls Whether or not the particle initiates a time-like shower
*/
ShowerParticle(const Particle & x, unsigned int pert, bool fs, bool tls=false)
: Particle(x), _isFinalState(fs),
_perturbative(pert), _initiatesTLS(tls), _x(1.0), _showerKinematics(),
_vMass(ZERO), _thePEGBase(&x) {}
//@}
public:
/**
* Access/Set various flags about the state of the particle
*/
//@{
/**
* Access the flag that tells if the particle is final state
* or initial state.
*/
bool isFinalState() const { return _isFinalState; }
/**
* Access the flag that tells if the particle is initiating a
* time like shower when it has been emitted in an initial state shower.
*/
bool initiatesTLS() const { return _initiatesTLS; }
/**
* Access the flag which tells us where the particle came from
* This is 0 for a particle produced in the shower, 1 if the particle came
* from the hard sub-process and 2 is it came from a decay.
*/
unsigned int perturbative() const { return _perturbative; }
//@}
/**
* Set/Get the momentum fraction for initial-state particles
*/
//@{
/**
* For an initial state particle get the fraction of the beam momentum
*/
void x(double x) { _x = x; }
/**
* For an initial state particle set the fraction of the beam momentum
*/
double x() const { return _x; }
//@}
/**
* Set/Get methods for the ShowerKinematics objects
*/
//@{
/**
* Access/ the ShowerKinematics object.
*/
const ShoKinPtr & showerKinematics() const { return _showerKinematics; }
/**
* Set the ShowerKinematics object.
*/
void showerKinematics(const ShoKinPtr in) { _showerKinematics = in; }
//@}
/**
* Members relating to the initial evolution scale and partner for the particle
*/
//@{
/**
* Veto emission at a given scale
*/
void vetoEmission(ShowerPartnerType::Type type, Energy scale);
/**
* Access to the evolution scales
*/
const EvolutionScales & scales() const {return scales_;}
/**
* Access to the evolution scales
*/
EvolutionScales & scales() {return scales_;}
/**
* Return the virtual mass\f$
*/
Energy virtualMass() const { return _vMass; }
/**
* Set the virtual mass
*/
void virtualMass(Energy mass) { _vMass = mass; }
/**
* Return the partner
*/
tShowerParticlePtr partner() const { return _partner; }
/**
* Set the partner
*/
void partner(const tShowerParticlePtr partner) { _partner = partner; }
/**
* Get the possible partners
*/
vector<EvolutionPartner> & partners() { return partners_; }
/**
* Add a possible partners
*/
void addPartner(EvolutionPartner in );
/**
* Clear the evolution partners
*/
void clearPartners() { partners_.clear(); }
/**
* Return the progenitor of the shower
*/
tShowerParticlePtr progenitor() const { return _progenitor; }
/**
* Set the progenitor of the shower
*/
void progenitor(const tShowerParticlePtr progenitor) { _progenitor = progenitor; }
//@}
/**
* Members to store and provide access to variables for a specific
* shower evolution scheme
*/
//@{
struct Parameters {
Parameters() : alpha(1.), beta(), ptx(), pty(), pt() {}
double alpha;
double beta;
Energy ptx;
Energy pty;
Energy pt;
};
/**
* Set the vector containing dimensionless variables
*/
Parameters & showerParameters() { return _parameters; }
//@}
/**
* If this particle came from the hard process get a pointer to ThePEG particle
* it came from
*/
const tcPPtr thePEGBase() const { return _thePEGBase; }
+public:
+
+ /**
+ * Extract the rho matrix including mapping needed in the shower
+ */
+ RhoDMatrix extractRhoMatrix(ShoKinPtr kinematics, bool forward);
+
+protected:
+
+ /**
+ * For a particle which came from the hard process get the spin density and
+ * the mapping required to the basis used in the Shower
+ * @param rho The \f$\rho\f$ matrix
+ * @param mapping The mapping
+ * @param showerkin The ShowerKinematics object
+ */
+ bool getMapping(SpinPtr &, RhoDMatrix & map, ShoKinPtr showerkin);
+
+
protected:
/**
* Standard clone function.
*/
virtual PPtr clone() const;
/**
* Standard clone function.
*/
virtual PPtr fullclone() const;
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<ShowerParticle> initShowerParticle;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
ShowerParticle & operator=(const ShowerParticle &);
private:
/**
* Whether the particle is in the final or initial state
*/
bool _isFinalState;
/**
* Whether the particle came from
*/
unsigned int _perturbative;
/**
* Does a particle produced in the backward shower initiate a time-like shower
*/
bool _initiatesTLS;
/**
* Dimensionless parameters
*/
Parameters _parameters;
/**
* The beam energy fraction for particle's in the initial state
*/
double _x;
/**
* The shower kinematics for the particle
*/
ShoKinPtr _showerKinematics;
/**
* Storage of the evolution scales
*/
EvolutionScales scales_;
/**
* Virtual mass
*/
Energy _vMass;
/**
* Partners
*/
tShowerParticlePtr _partner;
/**
* Pointer to ThePEG Particle this ShowerParticle was created from
*/
const tcPPtr _thePEGBase;
/**
* Progenitor
*/
tShowerParticlePtr _progenitor;
/**
* Partners
*/
vector<EvolutionPartner> partners_;
};
inline ostream & operator<<(ostream & os, const ShowerParticle::EvolutionScales & es) {
os << "Scales: QED=" << es.QED / GeV
<< " QCD_c=" << es.QCD_c / GeV
<< " QCD_ac=" << es.QCD_ac / GeV
<< " QED_noAO=" << es.QED_noAO / GeV
<< " QCD_c_noAO=" << es.QCD_c_noAO / GeV
<< " QCD_ac_noAO=" << es.QCD_ac_noAO / GeV
<< '\n';
return os;
}
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of ShowerParticle. */
template <>
struct BaseClassTrait<Herwig::ShowerParticle,1> {
/** Typedef of the first base class of ShowerParticle. */
typedef Particle NthBase;
};
/** This template specialization informs ThePEG about the name of
* the ShowerParticle class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::ShowerParticle>
: public ClassTraitsBase<Herwig::ShowerParticle> {
/** Return a platform-independent class name */
static string className() { return "Herwig::ShowerParticle"; }
/** Create a Event object. */
static TPtr create() { return TPtr::Create(Herwig::ShowerParticle(tcEventPDPtr(),true)); }
};
/** @endcond */
}
#endif /* HERWIG_ShowerParticle_H */
diff --git a/Shower/Base/SudakovFormFactor.cc b/Shower/Base/SudakovFormFactor.cc
--- a/Shower/Base/SudakovFormFactor.cc
+++ b/Shower/Base/SudakovFormFactor.cc
@@ -1,714 +1,326 @@
// -*- C++ -*-
//
// SudakovFormFactor.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 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 SudakovFormFactor class.
//
#include "SudakovFormFactor.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ShowerKinematics.h"
#include "ShowerParticle.h"
-#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
-#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
-#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
-#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
-#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "Herwig/Shower/ShowerHandler.h"
using namespace Herwig;
DescribeAbstractClass<SudakovFormFactor,Interfaced>
describeSudakovFormFactor ("Herwig::SudakovFormFactor","");
void SudakovFormFactor::persistentOutput(PersistentOStream & os) const {
os << splittingFn_ << alpha_ << pdfmax_ << particles_ << pdffactor_
<< a_ << b_ << ounit(c_,GeV) << ounit(kinCutoffScale_,GeV) << cutOffOption_
<< ounit(vgcut_,GeV) << ounit(vqcut_,GeV)
<< ounit(pTmin_,GeV) << ounit(pT2min_,GeV2)
<< theFactorizationScaleFactor << theRenormalizationScaleFactor;
}
void SudakovFormFactor::persistentInput(PersistentIStream & is, int) {
is >> splittingFn_ >> alpha_ >> pdfmax_ >> particles_ >> pdffactor_
>> a_ >> b_ >> iunit(c_,GeV) >> iunit(kinCutoffScale_,GeV) >> cutOffOption_
>> iunit(vgcut_,GeV) >> iunit(vqcut_,GeV)
>> iunit(pTmin_,GeV) >> iunit(pT2min_,GeV2)
>> theFactorizationScaleFactor >> theRenormalizationScaleFactor;
}
void SudakovFormFactor::Init() {
static ClassDocumentation<SudakovFormFactor> documentation
("The SudakovFormFactor class is the base class for the implementation of Sudakov"
" form factors in Herwig");
static Reference<SudakovFormFactor,SplittingFunction>
interfaceSplittingFunction("SplittingFunction",
"A reference to the SplittingFunction object",
&Herwig::SudakovFormFactor::splittingFn_,
false, false, true, false);
static Reference<SudakovFormFactor,ShowerAlpha>
interfaceAlpha("Alpha",
"A reference to the Alpha object",
&Herwig::SudakovFormFactor::alpha_,
false, false, true, false);
static Parameter<SudakovFormFactor,double> interfacePDFmax
("PDFmax",
"Maximum value of PDF weight. ",
&SudakovFormFactor::pdfmax_, 35.0, 1.0, 100000.0,
false, false, Interface::limited);
static Switch<SudakovFormFactor,unsigned int> interfacePDFFactor
("PDFFactor",
"Include additional factors in the overestimate for the PDFs",
&SudakovFormFactor::pdffactor_, 0, false, false);
static SwitchOption interfacePDFFactorOff
(interfacePDFFactor,
"Off",
"Don't include any factors",
0);
static SwitchOption interfacePDFFactorOverZ
(interfacePDFFactor,
"OverZ",
"Include an additional factor of 1/z",
1);
static SwitchOption interfacePDFFactorOverOneMinusZ
(interfacePDFFactor,
"OverOneMinusZ",
"Include an additional factor of 1/(1-z)",
2);
static SwitchOption interfacePDFFactorOverZOneMinusZ
(interfacePDFFactor,
"OverZOneMinusZ",
"Include an additional factor of 1/z/(1-z)",
3);
static Switch<SudakovFormFactor,unsigned int> interfaceCutOffOption
("CutOffOption",
"The type of cut-off to use to end the shower",
&SudakovFormFactor::cutOffOption_, 0, false, false);
static SwitchOption interfaceCutOffOptionDefault
(interfaceCutOffOption,
"Default",
"Use the standard Herwig cut-off on virtualities with the minimum"
" virtuality depending on the mass of the branching particle",
0);
static SwitchOption interfaceCutOffOptionFORTRAN
(interfaceCutOffOption,
"FORTRAN",
"Use a FORTRAN-like cut-off on virtualities",
1);
static SwitchOption interfaceCutOffOptionpT
(interfaceCutOffOption,
"pT",
"Use a cut on the minimum allowed pT",
2);
static Parameter<SudakovFormFactor,double> interfaceaParameter
("aParameter",
"The a parameter for the kinematic cut-off",
&SudakovFormFactor::a_, 0.3, -10.0, 10.0,
false, false, Interface::limited);
static Parameter<SudakovFormFactor,double> interfacebParameter
("bParameter",
"The b parameter for the kinematic cut-off",
&SudakovFormFactor::b_, 2.3, -10.0, 10.0,
false, false, Interface::limited);
static Parameter<SudakovFormFactor,Energy> interfacecParameter
("cParameter",
"The c parameter for the kinematic cut-off",
&SudakovFormFactor::c_, GeV, 0.3*GeV, 0.1*GeV, 10.0*GeV,
false, false, Interface::limited);
static Parameter<SudakovFormFactor,Energy>
interfaceKinScale ("cutoffKinScale",
"kinematic cutoff scale for the parton shower phase"
" space (unit [GeV])",
&SudakovFormFactor::kinCutoffScale_, GeV,
2.3*GeV, 0.001*GeV, 10.0*GeV,false,false,false);
static Parameter<SudakovFormFactor,Energy> interfaceGluonVirtualityCut
("GluonVirtualityCut",
"For the FORTRAN cut-off option the minimum virtuality of the gluon",
&SudakovFormFactor::vgcut_, GeV, 0.85*GeV, 0.1*GeV, 10.0*GeV,
false, false, Interface::limited);
static Parameter<SudakovFormFactor,Energy> interfaceQuarkVirtualityCut
("QuarkVirtualityCut",
"For the FORTRAN cut-off option the minimum virtuality added to"
" the mass for particles other than the gluon",
&SudakovFormFactor::vqcut_, GeV, 0.85*GeV, 0.1*GeV, 10.0*GeV,
false, false, Interface::limited);
static Parameter<SudakovFormFactor,Energy> interfacepTmin
("pTmin",
"The minimum pT if using a cut-off on the pT",
&SudakovFormFactor::pTmin_, GeV, 1.0*GeV, ZERO, 10.0*GeV,
false, false, Interface::limited);
}
bool SudakovFormFactor::alphaSVeto(Energy2 pt2) const {
pt2 *= sqr(renormalizationScaleFactor());
return UseRandom::rnd() > ThePEG::Math::powi(alpha_->ratio(pt2),
splittingFn_->interactionOrder());
}
bool SudakovFormFactor::
PDFVeto(const Energy2 t, const double x,
const tcPDPtr parton0, const tcPDPtr parton1,
Ptr<BeamParticleData>::transient_const_pointer beam) const {
assert(pdf_);
Energy2 theScale = t * sqr(factorizationScaleFactor());
if (theScale < sqr(freeze_)) theScale = sqr(freeze_);
double newpdf(0.0), oldpdf(0.0);
//different treatment of MPI ISR is done via CascadeHandler::resetPDFs()
newpdf=pdf_->xfx(beam,parton0,theScale,x/z());
oldpdf=pdf_->xfx(beam,parton1,theScale,x);
if(newpdf<=0.) return true;
if(oldpdf<=0.) return false;
double ratio = newpdf/oldpdf;
double maxpdf(pdfmax_);
switch (pdffactor_) {
case 1:
maxpdf /= z();
break;
case 2:
maxpdf /= 1.-z();
break;
case 3:
maxpdf /= (z()*(1.-z()));
break;
}
// ratio / PDFMax must be a probability <= 1.0
if (ratio > maxpdf) {
generator()->log() << "PDFVeto warning: Ratio > " << name()
<< ":PDFmax (by a factor of "
<< ratio/maxpdf <<") for "
<< parton0->PDGName() << " to "
<< parton1->PDGName() << "\n";
}
return ratio < UseRandom::rnd()*maxpdf;
}
void SudakovFormFactor::addSplitting(const IdList & in) {
bool add=true;
for(unsigned int ix=0;ix<particles_.size();++ix) {
if(particles_[ix].size()==in.size()) {
bool match=true;
for(unsigned int iy=0;iy<in.size();++iy) {
if(particles_[ix][iy]!=in[iy]) {
match=false;
break;
}
}
if(match) {
add=false;
break;
}
}
}
if(add) particles_.push_back(in);
}
-namespace {
-
-LorentzRotation boostToShower(const vector<Lorentz5Momentum> & basis,
- ShowerKinematics::Frame frame,
- Lorentz5Momentum & porig) {
- LorentzRotation output;
- if(frame==ShowerKinematics::BackToBack) {
- // we are doing the evolution in the back-to-back frame
- // work out the boostvector
- Boost boostv(-(basis[0]+basis[1]).boostVector());
- // momentum of the parton
- Lorentz5Momentum ptest(basis[0]);
- // construct the Lorentz boost
- output = LorentzRotation(boostv);
- ptest *= output;
- Axis axis(ptest.vect().unit());
- // now rotate so along the z axis as needed for the splitting functions
- if(axis.perp2()>1e-10) {
- double sinth(sqrt(1.-sqr(axis.z())));
- output.rotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
- }
- else if(axis.z()<0.) {
- output.rotate(Constants::pi,Axis(1.,0.,0.));
- }
- porig = output*basis[0];
- porig.setX(ZERO);
- porig.setY(ZERO);
- }
- else {
- output = LorentzRotation(-basis[0].boostVector());
- porig = output*basis[0];
- porig.setX(ZERO);
- porig.setY(ZERO);
- porig.setZ(ZERO);
- }
- return output;
-}
-
-RhoDMatrix bosonMapping(ShowerParticle & particle,
- const Lorentz5Momentum & porig,
- VectorSpinPtr vspin,
- const LorentzRotation & rot) {
- // rotate the original basis
- vector<LorentzPolarizationVector> sbasis;
- for(unsigned int ix=0;ix<3;++ix) {
- sbasis.push_back(vspin->getProductionBasisState(ix));
- sbasis.back().transform(rot);
- }
- // splitting basis
- vector<LorentzPolarizationVector> fbasis;
- bool massless(particle.id()==ParticleID::g||particle.id()==ParticleID::gamma);
- VectorWaveFunction wave(porig,particle.dataPtr(),outgoing);
- for(unsigned int ix=0;ix<3;++ix) {
- if(massless&&ix==1) {
- fbasis.push_back(LorentzPolarizationVector());
- }
- else {
- wave.reset(ix);
- fbasis.push_back(wave.wave());
- }
- }
- // work out the mapping
- RhoDMatrix mapping=RhoDMatrix(PDT::Spin1,false);
- for(unsigned int ix=0;ix<3;++ix) {
- for(unsigned int iy=0;iy<3;++iy) {
- mapping(ix,iy)= sbasis[iy].dot(fbasis[ix].conjugate());
- if(particle.id()<0)
- mapping(ix,iy)=conj(mapping(ix,iy));
- }
- }
- // \todo need to fix this
- mapping = RhoDMatrix(PDT::Spin1,false);
- if(massless) {
- mapping(0,0) = 1.;
- mapping(2,2) = 1.;
- }
- else {
- mapping(0,0) = 1.;
- mapping(1,1) = 1.;
- mapping(2,2) = 1.;
- }
- return mapping;
-}
-
-RhoDMatrix fermionMapping(ShowerParticle & particle,
- const Lorentz5Momentum & porig,
- FermionSpinPtr fspin,
- const LorentzRotation & rot) {
- // extract the original basis states
- vector<LorentzSpinor<SqrtEnergy> > sbasis;
- for(unsigned int ix=0;ix<2;++ix) {
- sbasis.push_back(fspin->getProductionBasisState(ix));
- sbasis.back().transform(rot);
- }
- // calculate the states in the splitting basis
- vector<LorentzSpinor<SqrtEnergy> > fbasis;
- SpinorWaveFunction wave(porig,particle.dataPtr(),
- particle.id()>0 ? incoming : outgoing);
- for(unsigned int ix=0;ix<2;++ix) {
- wave.reset(ix);
- fbasis.push_back(wave.dimensionedWave());
- }
- RhoDMatrix mapping=RhoDMatrix(PDT::Spin1Half,false);
- for(unsigned int ix=0;ix<2;++ix) {
- if(fbasis[0].s2()==complex<SqrtEnergy>()) {
- mapping(ix,0) = sbasis[ix].s3()/fbasis[0].s3();
- mapping(ix,1) = sbasis[ix].s2()/fbasis[1].s2();
- }
- else {
- mapping(ix,0) = sbasis[ix].s2()/fbasis[0].s2();
- mapping(ix,1) = sbasis[ix].s3()/fbasis[1].s3();
- }
- }
- return mapping;
-}
-
-FermionSpinPtr createFermionSpinInfo(ShowerParticle & particle,
- const Lorentz5Momentum & porig,
- const LorentzRotation & rot,
- Helicity::Direction dir) {
- // calculate the splitting basis for the branching
- // and rotate back to construct the basis states
- LorentzRotation rinv = rot.inverse();
- SpinorWaveFunction wave;
- if(particle.id()>0)
- wave=SpinorWaveFunction(porig,particle.dataPtr(),incoming);
- else
- wave=SpinorWaveFunction(porig,particle.dataPtr(),outgoing);
- FermionSpinPtr fspin = new_ptr(FermionSpinInfo(particle.momentum(),dir==outgoing));
- for(unsigned int ix=0;ix<2;++ix) {
- wave.reset(ix);
- LorentzSpinor<SqrtEnergy> basis = wave.dimensionedWave();
- basis.transform(rinv);
- fspin->setBasisState(ix,basis);
- fspin->setDecayState(ix,basis);
- }
- particle.spinInfo(fspin);
- return fspin;
-}
-
-VectorSpinPtr createVectorSpinInfo(ShowerParticle & particle,
- const Lorentz5Momentum & porig,
- const LorentzRotation & rot,
- Helicity::Direction dir) {
- // calculate the splitting basis for the branching
- // and rotate back to construct the basis states
- LorentzRotation rinv = rot.inverse();
- bool massless(particle.id()==ParticleID::g||particle.id()==ParticleID::gamma);
- VectorWaveFunction wave(porig,particle.dataPtr(),dir);
- VectorSpinPtr vspin = new_ptr(VectorSpinInfo(particle.momentum(),dir==outgoing));
- for(unsigned int ix=0;ix<3;++ix) {
- LorentzPolarizationVector basis;
- if(massless&&ix==1) {
- basis = LorentzPolarizationVector();
- }
- else {
- wave.reset(ix);
- basis = wave.wave();
- }
- basis *= rinv;
- vspin->setBasisState(ix,basis);
- vspin->setDecayState(ix,basis);
- }
- particle.spinInfo(vspin);
- vspin-> DMatrix() = RhoDMatrix(PDT::Spin1);
- vspin->rhoMatrix() = RhoDMatrix(PDT::Spin1);
- if(massless) {
- vspin-> DMatrix()(0,0) = 0.5;
- vspin->rhoMatrix()(0,0) = 0.5;
- vspin-> DMatrix()(2,2) = 0.5;
- vspin->rhoMatrix()(2,2) = 0.5;
- }
- return vspin;
-}
-}
-
-bool SudakovFormFactor::getMapping(SpinPtr & output, RhoDMatrix & mapping,
- ShowerParticle & particle,ShoKinPtr showerkin) {
- // if the particle is not from the hard process
- if(!particle.perturbative()) {
- // mapping is the identity
- output=particle.spinInfo();
- mapping=RhoDMatrix(particle.dataPtr()->iSpin());
- if(output) {
- return false;
- }
- else {
- Lorentz5Momentum porig;
- LorentzRotation rot = boostToShower(showerkin->getBasis(),showerkin->frame(),porig);
- Helicity::Direction dir = particle.isFinalState() ? outgoing : incoming;
- if(particle.dataPtr()->iSpin()==PDT::Spin0) {
- assert(false);
- }
- else if(particle.dataPtr()->iSpin()==PDT::Spin1Half) {
- output = createFermionSpinInfo(particle,porig,rot,dir);
- }
- else if(particle.dataPtr()->iSpin()==PDT::Spin1) {
- output = createVectorSpinInfo(particle,porig,rot,dir);
- }
- else {
- assert(false);
- }
- return false;
- }
- }
- // if particle is final-state and is from the hard process
- else if(particle.isFinalState()) {
- assert(particle.perturbative()==1 || particle.perturbative()==2);
- // get transform to shower frame
- Lorentz5Momentum porig;
- LorentzRotation rot = boostToShower(showerkin->getBasis(),showerkin->frame(),porig);
- // the rest depends on the spin of the particle
- PDT::Spin spin(particle.dataPtr()->iSpin());
- mapping=RhoDMatrix(spin,false);
- // do the spin dependent bit
- if(spin==PDT::Spin0) {
- ScalarSpinPtr sspin=dynamic_ptr_cast<ScalarSpinPtr>(particle.spinInfo());
- if(!sspin) {
- ScalarWaveFunction::constructSpinInfo(&particle,outgoing,true);
- }
- output=particle.spinInfo();
- return false;
- }
- else if(spin==PDT::Spin1Half) {
- FermionSpinPtr fspin=dynamic_ptr_cast<FermionSpinPtr>(particle.spinInfo());
- // spin info exists get information from it
- if(fspin) {
- output=fspin;
- mapping = fermionMapping(particle,porig,fspin,rot);
- return true;
- }
- // spin info does not exist create it
- else {
- output = createFermionSpinInfo(particle,porig,rot,outgoing);
- return false;
- }
- }
- else if(spin==PDT::Spin1) {
- VectorSpinPtr vspin=dynamic_ptr_cast<VectorSpinPtr>(particle.spinInfo());
- // spin info exists get information from it
- if(vspin) {
- output=vspin;
- mapping = bosonMapping(particle,porig,vspin,rot);
- return true;
- }
- else {
- output = createVectorSpinInfo(particle,porig,rot,outgoing);
- return false;
- }
- }
- // not scalar/fermion/vector
- else
- assert(false);
- }
- // incoming to hard process
- else if(particle.perturbative()==1 && !particle.isFinalState()) {
- // get the basis vectors
- // get transform to shower frame
- Lorentz5Momentum porig;
- LorentzRotation rot = boostToShower(showerkin->getBasis(),showerkin->frame(),porig);
- porig *= particle.x();
- // the rest depends on the spin of the particle
- PDT::Spin spin(particle.dataPtr()->iSpin());
- mapping=RhoDMatrix(spin);
- // do the spin dependent bit
- if(spin==PDT::Spin0) {
- cerr << "testing spin 0 not yet implemented " << endl;
- assert(false);
- }
- // spin-1/2
- else if(spin==PDT::Spin1Half) {
- FermionSpinPtr fspin=dynamic_ptr_cast<FermionSpinPtr>(particle.spinInfo());
- // spin info exists get information from it
- if(fspin) {
- output=fspin;
- mapping = fermionMapping(particle,porig,fspin,rot);
- return true;
- }
- // spin info does not exist create it
- else {
- output = createFermionSpinInfo(particle,porig,rot,incoming);
- return false;
- }
- }
- // spin-1
- else if(spin==PDT::Spin1) {
- VectorSpinPtr vspin=dynamic_ptr_cast<VectorSpinPtr>(particle.spinInfo());
- // spinInfo exists map it
- if(vspin) {
- output=vspin;
- mapping = bosonMapping(particle,porig,vspin,rot);
- return true;
- }
- // create the spininfo
- else {
- output = createVectorSpinInfo(particle,porig,rot,incoming);
- return false;
- }
- }
- assert(false);
- }
- // incoming to decay
- else if(particle.perturbative() == 2 && !particle.isFinalState()) {
- // get the basis vectors
- Lorentz5Momentum porig;
- LorentzRotation rot=boostToShower(showerkin->getBasis(),
- showerkin->frame(),porig);
- // the rest depends on the spin of the particle
- PDT::Spin spin(particle.dataPtr()->iSpin());
- mapping=RhoDMatrix(spin);
- // do the spin dependent bit
- if(spin==PDT::Spin0) {
- cerr << "testing spin 0 not yet implemented " << endl;
- assert(false);
- }
- // spin-1/2
- else if(spin==PDT::Spin1Half) {
- // FermionSpinPtr fspin=dynamic_ptr_cast<FermionSpinPtr>(particle.spinInfo());
- // // spin info exists get information from it
- // if(fspin) {
- // output=fspin;
- // mapping = fermionMapping(particle,porig,fspin,rot);
- // return true;
- // // spin info does not exist create it
- // else {
- // output = createFermionSpinInfo(particle,porig,rot,incoming);
- // return false;
- // }
- // }
- assert(false);
- }
- // // spin-1
- // else if(spin==PDT::Spin1) {
- // VectorSpinPtr vspin=dynamic_ptr_cast<VectorSpinPtr>(particle.spinInfo());
- // // spinInfo exists map it
- // if(vspin) {
- // output=vspin;
- // mapping = bosonMapping(particle,porig,vspin,rot);
- // return true;
- // }
- // // create the spininfo
- // else {
- // output = createVectorSpinInfo(particle,porig,rot,incoming);
- // return false;
- // }
- // }
- // assert(false);
- assert(false);
- }
- else
- assert(false);
- return true;
-}
-
void SudakovFormFactor::removeSplitting(const IdList & in) {
for(vector<IdList>::iterator it=particles_.begin();
it!=particles_.end();++it) {
if(it->size()==in.size()) {
bool match=true;
for(unsigned int iy=0;iy<in.size();++iy) {
if((*it)[iy]!=in[iy]) {
match=false;
break;
}
}
if(match) {
vector<IdList>::iterator itemp=it;
--itemp;
particles_.erase(it);
it = itemp;
}
}
}
}
Energy2 SudakovFormFactor::guesst(Energy2 t1,unsigned int iopt,
const IdList &ids,
double enhance,bool ident) const {
unsigned int pdfopt = iopt!=1 ? 0 : pdffactor_;
double c =
1./((splittingFn_->integOverP(zlimits_.second,ids,pdfopt) -
splittingFn_->integOverP(zlimits_.first ,ids,pdfopt))*
alpha_->overestimateValue()/Constants::twopi*enhance);
assert(iopt<=2);
if(iopt==1) {
c/=pdfmax_;
//symmetry of FS gluon splitting
if(ident) c*=0.5;
}
else if(iopt==2) c*=-1.;
if(splittingFn_->interactionOrder()==1) {
double r = UseRandom::rnd();
if(iopt!=2 || c*log(r)<log(Constants::MaxEnergy2/t1)) {
return t1*pow(r,c);
}
else
return Constants::MaxEnergy2;
}
else {
assert(false && "Units are dubious here.");
int nm(splittingFn()->interactionOrder()-1);
c/=Math::powi(alpha_->overestimateValue()/Constants::twopi,nm);
return t1 / pow (1. - nm*c*log(UseRandom::rnd())
* Math::powi(t1*UnitRemoval::InvE2,nm)
,1./double(nm));
}
}
double SudakovFormFactor::guessz (unsigned int iopt, const IdList &ids) const {
unsigned int pdfopt = iopt!=1 ? 0 : pdffactor_;
double lower = splittingFn_->integOverP(zlimits_.first,ids,pdfopt);
return splittingFn_->invIntegOverP
(lower + UseRandom::rnd()*(splittingFn_->integOverP(zlimits_.second,ids,pdfopt) -
lower),ids,pdfopt);
}
void SudakovFormFactor::doinit() {
Interfaced::doinit();
pT2min_ = cutOffOption()==2 ? sqr(pTmin_) : ZERO;
}
const vector<Energy> & SudakovFormFactor::virtualMasses(const IdList & ids) {
static vector<Energy> output;
output.clear();
if(cutOffOption() == 0) {
for(unsigned int ix=0;ix<ids.size();++ix)
output.push_back(getParticleData(ids[ix])->mass());
Energy kinCutoff=
kinematicCutOff(kinScale(),*std::max_element(output.begin(),output.end()));
for(unsigned int ix=0;ix<output.size();++ix)
output[ix]=max(kinCutoff,output[ix]);
}
else if(cutOffOption() == 1) {
for(unsigned int ix=0;ix<ids.size();++ix) {
output.push_back(getParticleData(ids[ix])->mass());
output.back() += ids[ix]==ParticleID::g ? vgCut() : vqCut();
}
}
else if(cutOffOption() == 2) {
for(unsigned int ix=0;ix<ids.size();++ix)
output.push_back(getParticleData(ids[ix])->mass());
}
else {
throw Exception() << "Unknown option for the cut-off"
<< " in SudakovFormFactor::virtualMasses()"
<< Exception::runerror;
}
return output;
}
-
-RhoDMatrix SudakovFormFactor::extractRhoMatrix(ShowerParticle & particle,
- ShoKinPtr kinematics,bool forward) {
- // get the spin density matrix and the mapping
- RhoDMatrix mapping;
- SpinPtr inspin;
- bool needMapping = getMapping(inspin,mapping,particle,kinematics);
- // set the decayed flag
- inspin->decay();
- // get the spin density matrix
- RhoDMatrix rho = forward ? inspin->rhoMatrix() : inspin->DMatrix();
- // map to the shower basis if needed
- if(needMapping) {
- RhoDMatrix rhop(rho.iSpin(),false);
- for(int ixa=0;ixa<rho.iSpin();++ixa) {
- for(int ixb=0;ixb<rho.iSpin();++ixb) {
- for(int iya=0;iya<rho.iSpin();++iya) {
- for(int iyb=0;iyb<rho.iSpin();++iyb) {
- rhop(ixa,ixb) += rho(iya,iyb)*mapping(iya,ixa)*conj(mapping(iyb,ixb));
- }
- }
- }
- }
- rhop.normalize();
- rho = rhop;
- }
- return rho;
-}
-
diff --git a/Shower/Base/SudakovFormFactor.h b/Shower/Base/SudakovFormFactor.h
--- a/Shower/Base/SudakovFormFactor.h
+++ b/Shower/Base/SudakovFormFactor.h
@@ -1,704 +1,690 @@
// -*- C++ -*-
//
// SudakovFormFactor.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 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_SudakovFormFactor_H
#define HERWIG_SudakovFormFactor_H
//
// This is the declaration of the SudakovFormFactor class.
//
#include "ThePEG/Interface/Interfaced.h"
#include "Herwig/Shower/SplittingFunctions/SplittingFunction.h"
#include "Herwig/Shower/Couplings/ShowerAlpha.h"
#include "Herwig/Shower/SplittingFunctions/SplittingGenerator.fh"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/PDF/BeamParticleData.h"
#include "ThePEG/EventRecord/RhoDMatrix.h"
#include "ThePEG/EventRecord/SpinInfo.h"
#include "ShowerKinematics.fh"
#include "SudakovFormFactor.fh"
namespace Herwig {
using namespace ThePEG;
/**
* A typedef for the BeamParticleData
*/
typedef Ptr<BeamParticleData>::transient_const_pointer tcBeamPtr;
/** \ingroup Shower
*
* This is the definition of the Sudakov form factor class. In general this
* is the base class for the implementation of Sudakov form factors in Herwig.
* The methods generateNextTimeBranching(), generateNextDecayBranching() and
* generateNextSpaceBranching need to be implemented in classes inheriting from this
* one.
*
* In addition a number of methods are implemented to assist with the calculation
* of the form factor using the veto algorithm in classes inheriting from this one.
*
* In general the Sudakov form-factor, for final-state radiation, is given
* by
* \f[\Delta_{ba}(\tilde{q}_{i+1},\tilde{q}_i)=
* \exp\left\{
* -\int^{\tilde{q}^2_i}_{\tilde{q}^2_{i+1}}
* \frac{{\rm d}\tilde{q}^2}{\tilde{q}^2}
* \int\frac{\alpha_S(z,\tilde{q})}{2\pi}
* P_{ba}(z,\tilde{q})\Theta(p_T)
* \right\}.
* \f]
* We can solve this to obtain the next value of the scale \f$\tilde{q}_{i+1}\f$
* given the previous value \f$\tilde{q}_i\f$
* in the following way. First we obtain a simplified form of the integrand
* which is greater than or equal to the true integrand for all values of
* \f$\tilde{q}\f$.
*
* In practice it is easiest to obtain this over estimate in pieces. The ShowerAlpha
* object contains an over estimate for \f$\alpha_S\f$, the splitting function
* contains both an over estimate of the spltting function and its integral
* which is needed to compute the over estimate of the \f$\tilde{q}\f$ integrand,
* together with an over estimate of the limit of the \f$z\f$ integral.
*
* This gives an overestimate of the integrand
* \f[g(\tilde{q}^2) = \frac{c}{\tilde{q}^2}, \f]
* where because the over estimates are chosen to be independent of \f$\tilde{q}\f$ the
* parameter
* \f[c = \frac{\alpha_{\rm over}}{2\pi}\int^{z_1}_{z_0}P_{\rm over}(z),\f]
* is a constant independent of \f$\tilde{q}\f$.
*
* The guesst() member can then be used to generate generate the value of
* \f$\tilde{q}^2\f$ according to this result. This is done by solving the Sudakov
* form factor, with the over estimates, is equal to a random number
* \f$r\f$ in the interval \f$[0,1]\f$. This gives
* \f[\tilde{q}^2_{i+1}=G^{-1}\left[G(\tilde{q}^2_i)+\ln r\right],\f]
* where \f$G(\tilde{q}^2)=c\ln(\tilde{q}^2)\f$ is the infinite integral
* of \f$g(\tilde{q}^2)\f$ and \f$G^{-1}(x)=\exp\left(\frac{x}c\right)\f$
* is its inverse.
* It this case we therefore obtain
* \f[\tilde{q}^2_{i+1}=\tilde{q}^2_ir^{\frac1c}.\f]
* The value of \f$z\f$ can then be calculated in a similar way
* \f[z = I^{-1}\left[I(z_0)+r\left(I(z_1)-I(z_0)\right)\right],\f]
* using the guessz() member,
* where \f$I=\int P(z){\rm d}z\f$ and \f$I^{-1}\f$ is its inverse.
*
* The veto algorithm then uses rejection using the ratio of the
* true value to the overestimated one to obtain the original distribution.
* This is accomplished using the
* - alphaSVeto() member for the \f$\alpha_S\f$ veto
* - SplittingFnVeto() member for the veto on the value of the splitting function.
* in general there must also be a chech that the emission is in the allowed
* phase space but this is left to the inheriting classes as it will depend
* on the ordering variable.
*
* The Sudakov form factor for the initial-scale shower is different because
* it must include the PDF which guides the backward evolution.
* It is given by
* \f[\Delta_{ba}(\tilde{q}_{i+1},\tilde{q}_i)=
* \exp\left\{
* -\int^{\tilde{q}^2_i}_{\tilde{q}^2_{i+1}}
* \frac{{\rm d}\tilde{q}^2}{\tilde{q}^2}
* \int\frac{\alpha_S(z,\tilde{q})}{2\pi}
* P_{ba}(z,\tilde{q})\frac{x'f_a(\frac{x}z,\tilde{q}^2)}{xf_b(x,\tilde{q^2})}
* \right\},
* \f]
* where \f$x\f$ is the fraction of the beam momentum the parton \f$b\f$ had before
* the backward evolution.
* This can be solve in the same way as for the final-state branching but the constant
* becomes
* \f[c = \frac{\alpha_{\rm over}}{2\pi}\int^{z_1}_{z_0}P_{\rm over}(z)PDF_{\rm max},\f]
* where
* \f[PDF_{\rm max}=\max\frac{x'f_a(\frac{x}z,\tilde{q}^2)}{xf_b(x,\tilde{q^2})},\f]
* which can be set using an interface.
* In addition the PDFVeto() member then is needed to implement the relevant veto.
*
* @see SplittingGenerator
* @see SplittingFunction
* @see ShowerAlpha
* @see \ref SudakovFormFactorInterfaces "The interfaces"
* defined for SudakovFormFactor.
*/
class SudakovFormFactor: public Interfaced {
/**
* The SplittingGenerator is a friend to insert the particles in the
* branchings at initialisation
*/
friend class SplittingGenerator;
public:
/**
* The default constructor.
*/
SudakovFormFactor() : pdfmax_(35.0), pdffactor_(0),
cutOffOption_(0), a_(0.3), b_(2.3), c_(0.3*GeV),
kinCutoffScale_( 2.3*GeV ), vgcut_(0.85*GeV),
vqcut_(0.85*GeV), pTmin_(1.*GeV), pT2min_(ZERO),
z_( 0.0 ),phi_(0.0), pT_(),
theFactorizationScaleFactor(1.0),
theRenormalizationScaleFactor(1.0) {}
/**
* Members to generate the scale of the next branching
*/
//@{
/**
* Return the scale of the next time-like branching. If there is no
* branching then it returns ZERO.
* @param startingScale starting scale for the evolution
* @param ids The PDG codes of the particles in the splitting
* @param cc Whether this is the charge conjugate of the branching
* @param enhance The radiation enhancement factor
* @param maxQ2 The maximum \f$Q^2\f$ for the emission
*/
virtual ShoKinPtr generateNextTimeBranching(const Energy startingScale,
const IdList &ids,const bool cc,
double enhance, Energy2 maxQ2)=0;
/**
* Return the scale of the next space-like decay branching. If there is no
* branching then it returns ZERO.
* @param startingScale starting scale for the evolution
* @param stoppingScale stopping scale for the evolution
* @param minmass The minimum mass allowed for the spake-like particle.
* @param ids The PDG codes of the particles in the splitting
* @param cc Whether this is the charge conjugate of the branching
* defined.
* @param enhance The radiation enhancement factor
*/
virtual ShoKinPtr generateNextDecayBranching(const Energy startingScale,
const Energy stoppingScale,
const Energy minmass,
const IdList &ids,
const bool cc,
double enhance)=0;
/**
* Return the scale of the next space-like branching. If there is no
* branching then it returns ZERO.
* @param startingScale starting scale for the evolution
* @param ids The PDG codes of the particles in the splitting
* @param x The fraction of the beam momentum
* @param cc Whether this is the charge conjugate of the branching
* defined.
* @param beam The beam particle
* @param enhance The radiation enhancement factor
*/
virtual ShoKinPtr generateNextSpaceBranching(const Energy startingScale,
const IdList &ids,double x,
const bool cc,double enhance,
tcBeamPtr beam)=0;
//@}
/**
* Generate the azimuthal angle of the branching for forward evolution
* @param particle The branching particle
* @param ids The PDG codes of the particles in the branchings
* @param The Shower kinematics
*/
virtual double generatePhiForward(ShowerParticle & particle,const IdList & ids,
ShoKinPtr kinematics)=0;
/**
* Generate the azimuthal angle of the branching for backward evolution
* @param particle The branching particle
* @param ids The PDG codes of the particles in the branchings
* @param The Shower kinematics
*/
virtual double generatePhiBackward(ShowerParticle & particle,const IdList & ids,
ShoKinPtr kinematics)=0;
/**
* Generate the azimuthal angle of the branching for ISR in decays
* @param particle The branching particle
* @param ids The PDG codes of the particles in the branchings
* @param The Shower kinematics
*/
virtual double generatePhiDecay(ShowerParticle & particle,const IdList & ids,
ShoKinPtr kinematics)=0;
/**
* Methods to provide public access to the private member variables
*/
//@{
/**
* Return the pointer to the SplittingFunction object.
*/
tSplittingFnPtr splittingFn() const { return splittingFn_; }
/**
* Return the pointer to the ShowerAlpha object.
*/
tShowerAlphaPtr alpha() const { return alpha_; }
/**
* The type of interaction
*/
inline ShowerInteraction::Type interactionType() const
{return splittingFn_->interactionType();}
//@}
public:
/**
* Methods to access the kinematic variables for the branching
*/
//@{
/**
* The energy fraction
*/
double z() const { return z_; }
/**
* The azimuthal angle
*/
double phi() const { return phi_; }
/**
* The transverse momentum
*/
Energy pT() const { return pT_; }
//@}
/**
* Access the maximum weight for the PDF veto
*/
double pdfMax() const { return pdfmax_;}
/**
* Method to return the evolution scale given the
* transverse momentum, \f$p_T\f$ and \f$z\f$.
*/
virtual Energy calculateScale(double z, Energy pt, IdList ids,unsigned int iopt)=0;
/**
* Method to create the ShowerKinematics object for a final-state branching
*/
virtual ShoKinPtr createFinalStateBranching(Energy scale,double z,
double phi, Energy pt)=0;
/**
* Method to create the ShowerKinematics object for an initial-state branching
*/
virtual ShoKinPtr createInitialStateBranching(Energy scale,double z,
double phi, Energy pt)=0;
/**
* Method to create the ShowerKinematics object for a decay branching
*/
virtual ShoKinPtr createDecayBranching(Energy scale,double z,
double phi, Energy pt)=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 Standard Interfaced functions. */
//@{
/**
* 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();
//@}
protected:
/**
* Methods to implement the veto algorithm to generate the scale of
* the next branching
*/
//@{
/**
* Value of the energy fraction for the veto algorithm
* @param iopt The option for calculating z
* @param ids The PDG codes of the particles in the splitting
* - 0 is final-state
* - 1 is initial-state for the hard process
* - 2 is initial-state for particle decays
*/
double guessz (unsigned int iopt, const IdList &ids) const;
/**
* Value of the scale for the veto algorithm
* @param t1 The starting valoe of the scale
* @param iopt The option for calculating t
* @param ids The PDG codes of the particles in the splitting
* - 0 is final-state
* - 1 is initial-state for the hard process
* - 2 is initial-state for particle decays
* @param enhance The radiation enhancement factor
* @param identical Whether or not the outgoing particles are identical
*/
Energy2 guesst (Energy2 t1,unsigned int iopt, const IdList &ids,
double enhance, bool identical) const;
/**
* Veto on the PDF for the initial-state shower
* @param t The scale
* @param x The fraction of the beam momentum
* @param parton0 Pointer to the particleData for the
* new parent (this is the particle we evolved back to)
* @param parton1 Pointer to the particleData for the
* original particle
* @param beam The BeamParticleData object
*/
bool PDFVeto(const Energy2 t, const double x,
const tcPDPtr parton0, const tcPDPtr parton1,
tcBeamPtr beam) const;
/**
* The veto on the splitting function.
* @param t The scale
* @param ids The PDG codes of the particles in the splitting
* @param mass Whether or not to use the massive splitting functions
* @return true if vetoed
*/
bool SplittingFnVeto(const Energy2 t,
const IdList &ids,
const bool mass) const {
return UseRandom::rnd()>splittingFn_->ratioP(z_, t, ids,mass);
}
/**
* The veto on the coupling constant
* @param pt2 The value of ther transverse momentum squared, \f$p_T^2\f$.
* @return true if vetoed
*/
bool alphaSVeto(Energy2 pt2) const;
//@}
/**
* Methods to set the kinematic variables for the branching
*/
//@{
/**
* The energy fraction
*/
void z(double in) { z_=in; }
/**
* The azimuthal angle
*/
void phi(double in) { phi_=in; }
/**
* The transverse momentum
*/
void pT(Energy in) { pT_=in; }
//@}
/**
* Set/Get the limits on the energy fraction for the splitting
*/
//@{
/**
* Get the limits
*/
pair<double,double> zLimits() const { return zlimits_;}
/**
* Set the limits
*/
void zLimits(pair<double,double> in) { zlimits_=in; }
//@}
/**
* Set the particles in the splittings
*/
void addSplitting(const IdList &);
/**
* Delete the particles in the splittings
*/
void removeSplitting(const IdList &);
/**
* Access the potential branchings
*/
const vector<IdList> & particles() const { return particles_; }
- /**
- * For a particle which came from the hard process get the spin density and
- * the mapping required to the basis used in the Shower
- * @param rho The \f$\rho\f$ matrix
- * @param mapping The mapping
- * @param particle The particle
- * @param showerkin The ShowerKinematics object
- */
- bool getMapping(SpinPtr &, RhoDMatrix & map,
- ShowerParticle & particle,ShoKinPtr showerkin);
-
- RhoDMatrix extractRhoMatrix(ShowerParticle & particle,
- ShoKinPtr kinematics, bool forward);
-
public:
/**
* @name Methods for the cut-off
*/
//@{
/**
* The option being used
*/
unsigned int cutOffOption() const { return cutOffOption_; }
/**
* The kinematic scale
*/
Energy kinScale() const {return kinCutoffScale_;}
/**
* The virtuality cut-off on the gluon \f$Q_g=\frac{\delta-am_q}{b}\f$
* @param scale The scale \f$\delta\f$
* @param mq The quark mass \f$m_q\f$.
*/
Energy kinematicCutOff(Energy scale, Energy mq) const
{return max((scale -a_*mq)/b_,c_);}
/**
* The virtualilty cut-off for gluons
*/
Energy vgCut() const { return vgcut_; }
/**
* The virtuality cut-off for everything else
*/
Energy vqCut() const { return vqcut_; }
/**
* The minimum \f$p_T\f$ for the branching
*/
Energy pTmin() const { return pTmin_; }
/**
* The square of the minimum \f$p_T\f$
*/
Energy2 pT2min() const { return pT2min_; }
/**
* Calculate the virtual masses for a branchings
*/
const vector<Energy> & virtualMasses(const IdList & ids);
//@}
/**
* Set the PDF
*/
void setPDF(tcPDFPtr pdf, Energy scale) {
pdf_ = pdf;
freeze_ = scale;
}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
SudakovFormFactor & operator=(const SudakovFormFactor &);
private:
/**
* Pointer to the splitting function for this Sudakov form factor
*/
SplittingFnPtr splittingFn_;
/**
* Pointer to the coupling for this Sudakov form factor
*/
ShowerAlphaPtr alpha_;
/**
* Maximum value of the PDF weight
*/
double pdfmax_;
/**
* List of the particles this Sudakov is used for to aid in setting up
* interpolation tables if needed
*/
vector<IdList> particles_;
/**
* Option for the inclusion of a factor \f$1/(1-z)\f$ in the PDF estimate
*/
unsigned pdffactor_;
private:
/**
* Option for the type of cut-off to be applied
*/
unsigned int cutOffOption_;
/**
* Parameters for the default Herwig cut-off option, i.e. the parameters for
* the \f$Q_g=\max(\frac{\delta-am_q}{b},c)\f$ kinematic cut-off
*/
//@{
/**
* The \f$a\f$ parameter
*/
double a_;
/**
* The \f$b\f$ parameter
*/
double b_;
/**
* The \f$c\f$ parameter
*/
Energy c_;
/**
* Kinematic cutoff used in the parton shower phase space.
*/
Energy kinCutoffScale_;
//@}
/**
* Parameters for the FORTRAN-like cut-off
*/
//@{
/**
* The virtualilty cut-off for gluons
*/
Energy vgcut_;
/**
* The virtuality cut-off for everything else
*/
Energy vqcut_;
//@}
/**
* Parameters for the \f$p_T\f$ cut-off
*/
//@{
/**
* The minimum \f$p_T\f$ for the branching
*/
Energy pTmin_;
/**
* The square of the minimum \f$p_T\f$
*/
Energy2 pT2min_;
//@}
private:
/**
* Member variables to keep the shower kinematics information
* generated by a call to generateNextTimeBranching or generateNextSpaceBranching
*/
//@{
/**
* The energy fraction
*/
double z_;
/**
* The azimuthal angle
*/
double phi_;
/**
* The transverse momentum
*/
Energy pT_;
//@}
/**
* The limits of \f$z\f$ in the splitting
*/
pair<double,double> zlimits_;
/**
* Stuff for the PDFs
*/
//@{
/**
* PDf
*/
tcPDFPtr pdf_;
/**
* Freezing scale
*/
Energy freeze_;
//@}
public:
/**
* Get the factorization scale factor
*/
double factorizationScaleFactor() const { return theFactorizationScaleFactor; }
/**
* Set the factorization scale factor
*/
void factorizationScaleFactor(double f) { theFactorizationScaleFactor = f; }
/**
* Get the renormalization scale factor
*/
double renormalizationScaleFactor() const { return theRenormalizationScaleFactor; }
/**
* Set the renormalization scale factor
*/
void renormalizationScaleFactor(double f) { theRenormalizationScaleFactor = f; }
private:
/**
* The factorization scale factor.
*/
double theFactorizationScaleFactor;
/**
* The renormalization scale factor.
*/
double theRenormalizationScaleFactor;
};
}
#endif /* HERWIG_SudakovFormFactor_H */
diff --git a/Shower/Default/QTildeSudakov.cc b/Shower/Default/QTildeSudakov.cc
--- a/Shower/Default/QTildeSudakov.cc
+++ b/Shower/Default/QTildeSudakov.cc
@@ -1,939 +1,939 @@
// -*- C++ -*-
//
// QTildeSudakov.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 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 QTildeSudakov class.
//
#include "QTildeSudakov.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/PDT/ParticleData.h"
#include "ThePEG/EventRecord/Event.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Repository/CurrentGenerator.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "Herwig/Shower/Default/FS_QTildeShowerKinematics1to2.h"
#include "Herwig/Shower/Default/IS_QTildeShowerKinematics1to2.h"
#include "Herwig/Shower/Default/Decay_QTildeShowerKinematics1to2.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "Herwig/Shower/Base/ShowerVertex.h"
#include "Herwig/Shower/Base/ShowerParticle.h"
#include "Herwig/Shower/ShowerHandler.h"
#include "Herwig/Shower/Base/Evolver.h"
#include "Herwig/Shower/Base/PartnerFinder.h"
#include "Herwig/Shower/Base/ShowerModel.h"
#include "Herwig/Shower/Base/KinematicsReconstructor.h"
using namespace Herwig;
DescribeNoPIOClass<QTildeSudakov,Herwig::SudakovFormFactor>
describeQTildeSudakov ("Herwig::QTildeSudakov","HwShower.so");
void QTildeSudakov::Init() {
static ClassDocumentation<QTildeSudakov> documentation
("The QTildeSudakov class implements the Sudakov form factor for ordering it"
" qtilde");
}
bool QTildeSudakov::guessTimeLike(Energy2 &t,Energy2 tmin,double enhance) {
Energy2 told = t;
// calculate limits on z and if lower>upper return
if(!computeTimeLikeLimits(t)) return false;
// guess values of t and z
t = guesst(told,0,ids_,enhance,ids_[1]==ids_[2]);
z(guessz(0,ids_));
// actual values for z-limits
if(!computeTimeLikeLimits(t)) return false;
if(t<tmin) {
t=-1.0*GeV2;
return false;
}
else
return true;
}
bool QTildeSudakov::guessSpaceLike(Energy2 &t, Energy2 tmin, const double x,
double enhance) {
Energy2 told = t;
// calculate limits on z if lower>upper return
if(!computeSpaceLikeLimits(t,x)) return false;
// guess values of t and z
t = guesst(told,1,ids_,enhance,ids_[1]==ids_[2]);
z(guessz(1,ids_));
// actual values for z-limits
if(!computeSpaceLikeLimits(t,x)) return false;
if(t<tmin) {
t=-1.0*GeV2;
return false;
}
else
return true;
}
bool QTildeSudakov::PSVeto(const Energy2 t,
const Energy2 maxQ2) {
// still inside PS, return true if outside
// check vs overestimated limits
if(z() < zLimits().first || z() > zLimits().second) return true;
Energy2 q2 = z()*(1.-z())*t;
if(ids_[0]!=ParticleID::g &&
ids_[0]!=ParticleID::gamma ) q2 += masssquared_[0];
if(q2>maxQ2) return true;
// compute the pts
Energy2 pt2 = z()*(1.-z())*q2 - masssquared_[1]*(1.-z()) - masssquared_[2]*z();
// if pt2<0 veto
if(pt2<pT2min()) return true;
// otherwise calculate pt and return
pT(sqrt(pt2));
return false;
}
ShoKinPtr QTildeSudakov::generateNextTimeBranching(const Energy startingScale,
const IdList &ids,const bool cc,
double enhance, Energy2 maxQ2) {
// First reset the internal kinematics variables that can
// have been eventually set in the previous call to the method.
q_ = ZERO;
z(0.);
phi(0.);
// perform initialization
Energy2 tmax(sqr(startingScale)),tmin;
initialize(ids,tmin,cc);
// check max > min
if(tmax<=tmin) return ShoKinPtr();
// calculate next value of t using veto algorithm
Energy2 t(tmax);
do {
if(!guessTimeLike(t,tmin,enhance)) break;
}
while(PSVeto(t,maxQ2) || SplittingFnVeto(z()*(1.-z())*t,ids,true) ||
alphaSVeto(splittingFn()->angularOrdered() ? sqr(z()*(1.-z()))*t : z()*(1.-z())*t));
q_ = t > ZERO ? Energy(sqrt(t)) : -1.*MeV;
if(q_ < ZERO) return ShoKinPtr();
// return the ShowerKinematics object
return createFinalStateBranching(q_,z(),phi(),pT());
}
ShoKinPtr QTildeSudakov::
generateNextSpaceBranching(const Energy startingQ,
const IdList &ids,
double x,bool cc,
double enhance,
Ptr<BeamParticleData>::transient_const_pointer beam) {
// First reset the internal kinematics variables that can
// have been eventually set in the previous call to the method.
q_ = ZERO;
z(0.);
phi(0.);
// perform the initialization
Energy2 tmax(sqr(startingQ)),tmin;
initialize(ids,tmin,cc);
// check max > min
if(tmax<=tmin) return ShoKinPtr();
// extract the partons which are needed for the PDF veto
// Different order, incoming parton is id = 1, outgoing are id=0,2
tcPDPtr parton0 = getParticleData(ids[0]);
tcPDPtr parton1 = getParticleData(ids[1]);
if(cc) {
if(parton0->CC()) parton0 = parton0->CC();
if(parton1->CC()) parton1 = parton1->CC();
}
// calculate next value of t using veto algorithm
Energy2 t(tmax),pt2(ZERO);
do {
if(!guessSpaceLike(t,tmin,x,enhance)) break;
pt2=sqr(1.-z())*t-z()*masssquared_[2];
}
while(z() > zLimits().second ||
SplittingFnVeto((1.-z())*t/z(),ids,true) ||
alphaSVeto(splittingFn()->angularOrdered() ? sqr(1.-z())*t : (1.-z())*t) ||
PDFVeto(t,x,parton0,parton1,beam) || pt2 < pT2min() );
if(t > ZERO && zLimits().first < zLimits().second) q_ = sqrt(t);
else return ShoKinPtr();
pT(sqrt(pt2));
// create the ShowerKinematics and return it
return createInitialStateBranching(q_,z(),phi(),pT());
}
void QTildeSudakov::initialize(const IdList & ids, Energy2 & tmin,const bool cc) {
ids_=ids;
if(cc) {
for(unsigned int ix=0;ix<ids.size();++ix) {
if(getParticleData(ids[ix])->CC()) ids_[ix]*=-1;
}
}
tmin = cutOffOption() != 2 ? ZERO : 4.*pT2min();
masses_ = virtualMasses(ids);
masssquared_.clear();
for(unsigned int ix=0;ix<masses_.size();++ix) {
masssquared_.push_back(sqr(masses_[ix]));
if(ix>0) tmin=max(masssquared_[ix],tmin);
}
}
ShoKinPtr QTildeSudakov::generateNextDecayBranching(const Energy startingScale,
const Energy stoppingScale,
const Energy minmass,
const IdList &ids,
const bool cc,
double enhance) {
// First reset the internal kinematics variables that can
// have been eventually set in the previous call to this method.
q_ = Constants::MaxEnergy;
z(0.);
phi(0.);
// perform initialisation
Energy2 tmax(sqr(stoppingScale)),tmin;
initialize(ids,tmin,cc);
tmin=sqr(startingScale);
// check some branching possible
if(tmax<=tmin) return ShoKinPtr();
// perform the evolution
Energy2 t(tmin),pt2(-MeV2);
do {
if(!guessDecay(t,tmax,minmass,enhance)) break;
pt2 = sqr(1.-z())*(t-masssquared_[0])-z()*masssquared_[2];
}
while(SplittingFnVeto((1.-z())*t/z(),ids,true)||
alphaSVeto(splittingFn()->angularOrdered() ? sqr(1.-z())*t : (1.-z())*t ) ||
pt2<pT2min() ||
t*(1.-z())>masssquared_[0]-sqr(minmass));
if(t > ZERO) {
q_ = sqrt(t);
pT(sqrt(pt2));
}
else return ShoKinPtr();
phi(0.);
// create the ShowerKinematics object
return createDecayBranching(q_,z(),phi(),pT());
}
bool QTildeSudakov::guessDecay(Energy2 &t,Energy2 tmax, Energy minmass,
double enhance) {
// previous scale
Energy2 told = t;
// overestimated limits on z
if(tmax<masssquared_[0]) {
t=-1.0*GeV2;
return false;
}
Energy2 tm2 = tmax-masssquared_[0];
Energy tm = sqrt(tm2);
pair<double,double> limits=make_pair(sqr(minmass/masses_[0]),
1.-sqrt(masssquared_[2]+pT2min()+
0.25*sqr(masssquared_[2])/tm2)/tm
+0.5*masssquared_[2]/tm2);
zLimits(limits);
if(zLimits().second<zLimits().first) {
t=-1.0*GeV2;
return false;
}
// guess values of t and z
t = guesst(told,2,ids_,enhance,ids_[1]==ids_[2]);
z(guessz(2,ids_));
// actual values for z-limits
if(t<masssquared_[0]) {
t=-1.0*GeV2;
return false;
}
tm2 = t-masssquared_[0];
tm = sqrt(tm2);
limits=make_pair(sqr(minmass/masses_[0]),
1.-sqrt(masssquared_[2]+pT2min()+
0.25*sqr(masssquared_[2])/tm2)/tm
+0.5*masssquared_[2]/tm2);
zLimits(limits);
if(t>tmax||zLimits().second<zLimits().first) {
t=-1.0*GeV2;
return false;
}
else
return true;
}
bool QTildeSudakov::computeTimeLikeLimits(Energy2 & t) {
if (t < 1e-20 * GeV2) {
t=-1.*GeV2;
return false;
}
// special case for gluon radiating
pair<double,double> limits;
if(ids_[0]==ParticleID::g||ids_[0]==ParticleID::gamma) {
// no emission possible
if(t<16.*(masssquared_[1]+pT2min())) {
t=-1.*GeV2;
return false;
}
// overestimate of the limits
limits.first = 0.5*(1.-sqrt(1.-4.*sqrt((masssquared_[1]+pT2min())/t)));
limits.second = 1.-limits.first;
}
// special case for radiated particle is gluon
else if(ids_[2]==ParticleID::g||ids_[2]==ParticleID::gamma) {
limits.first = sqrt((masssquared_[1]+pT2min())/t);
limits.second = 1.-sqrt((masssquared_[2]+pT2min())/t);
}
else if(ids_[1]==ParticleID::g||ids_[1]==ParticleID::gamma) {
limits.second = sqrt((masssquared_[2]+pT2min())/t);
limits.first = 1.-sqrt((masssquared_[1]+pT2min())/t);
}
else {
limits.first = (masssquared_[1]+pT2min())/t;
limits.second = 1.-(masssquared_[2]+pT2min())/t;
}
if(limits.first>=limits.second) {
t=-1.*GeV2;
return false;
}
zLimits(limits);
return true;
}
bool QTildeSudakov::computeSpaceLikeLimits(Energy2 & t, double x) {
if (t < 1e-20 * GeV2) {
t=-1.*GeV2;
return false;
}
pair<double,double> limits;
// compute the limits
limits.first = x;
double yy = 1.+0.5*masssquared_[2]/t;
limits.second = yy - sqrt(sqr(yy)-1.+pT2min()/t);
// return false if lower>upper
zLimits(limits);
if(limits.second<limits.first) {
t=-1.*GeV2;
return false;
}
else
return true;
}
namespace {
tShowerParticlePtr findCorrelationPartner(ShowerParticle & particle,
bool forward,
ShowerInteraction::Type inter) {
tPPtr child = &particle;
tShowerParticlePtr mother;
if(forward) {
mother = !particle.parents().empty() ?
dynamic_ptr_cast<tShowerParticlePtr>(particle.parents()[0]) : tShowerParticlePtr();
}
else {
mother = particle.children().size()==2 ?
dynamic_ptr_cast<tShowerParticlePtr>(&particle) : tShowerParticlePtr();
}
tShowerParticlePtr partner;
while(mother) {
tPPtr otherChild;
if(forward) {
for (unsigned int ix=0;ix<mother->children().size();++ix) {
if(mother->children()[ix]!=child) {
otherChild = mother->children()[ix];
break;
}
}
}
else {
otherChild = mother->children()[1];
}
tShowerParticlePtr other = dynamic_ptr_cast<tShowerParticlePtr>(otherChild);
if((inter==ShowerInteraction::QCD && otherChild->dataPtr()->coloured()) ||
(inter==ShowerInteraction::QED && otherChild->dataPtr()->charged())) {
partner = other;
break;
}
if(forward && !other->isFinalState()) {
partner = dynamic_ptr_cast<tShowerParticlePtr>(mother);
break;
}
child = mother;
if(forward) {
mother = ! mother->parents().empty() ?
dynamic_ptr_cast<tShowerParticlePtr>(mother->parents()[0]) : tShowerParticlePtr();
}
else {
if(mother->children()[0]->children().size()!=2)
break;
tShowerParticlePtr mtemp =
dynamic_ptr_cast<tShowerParticlePtr>(mother->children()[0]);
if(!mtemp)
break;
else
mother=mtemp;
}
}
if(!partner) {
if(forward) {
partner = dynamic_ptr_cast<tShowerParticlePtr>( child)->partner();
}
else {
if(mother) {
tShowerParticlePtr parent;
if(!mother->children().empty()) {
parent = dynamic_ptr_cast<tShowerParticlePtr>(mother->children()[0]);
}
if(!parent) {
parent = dynamic_ptr_cast<tShowerParticlePtr>(mother);
}
partner = parent->partner();
}
else {
partner = dynamic_ptr_cast<tShowerParticlePtr>(&particle)->partner();
}
}
}
return partner;
}
pair<double,double> softPhiMin(double phi0, double phi1, double A, double B, double C, double D) {
double c01 = cos(phi0 - phi1);
double s01 = sin(phi0 - phi1);
double s012(sqr(s01)), c012(sqr(c01));
double A2(A*A), B2(B*B), C2(C*C), D2(D*D);
if(abs(B/A)<1e-10 && abs(D/C)<1e-10) return make_pair(phi0,phi0+Constants::pi);
double root = sqr(B2)*C2*D2*sqr(s012) + 2.*A*B2*B*C2*C*D*c01*s012 + 2.*A*B2*B*C*D2*D*c01*s012
+ 4.*A2*B2*C2*D2*c012 - A2*B2*C2*D2*s012 - A2*B2*sqr(D2)*s012 - sqr(B2)*sqr(C2)*s012
- sqr(B2)*C2*D2*s012 - 4.*A2*A*B*C*D2*D*c01 - 4.*A*B2*B*C2*C*D*c01 + sqr(A2)*sqr(D2)
+ 2.*A2*B2*C2*D2 + sqr(B2)*sqr(C2);
if(root<0.) return make_pair(phi0,phi0+Constants::pi);
root = sqrt(root);
double denom = (-2.*A*B*C*D*c01 + A2*D2 + B2*C2);
double denom2 = (-B*C*c01 + A*D);
double num = B2*C*D*s012;
return make_pair(atan2(B*s01*(-C*(num + root) / denom + D) / denom2, -(num + root ) / denom) + phi0,
atan2(B*s01*(-C*(num - root) / denom + D) / denom2, -(num - root ) / denom) + phi0);
}
}
double QTildeSudakov::generatePhiForward(ShowerParticle & particle,
const IdList & ids,
ShoKinPtr kinematics) {
// no correlations, return flat phi
if(! ShowerHandler::currentHandler()->evolver()->correlations())
return Constants::twopi*UseRandom::rnd();
// get the kinematic variables
double z = kinematics->z();
Energy2 t = z*(1.-z)*sqr(kinematics->scale());
Energy pT = kinematics->pT();
// if soft correlations
Energy2 pipj,pik;
bool canBeSoft[2] = {ids[1]==ParticleID::g || ids[1]==ParticleID::gamma,
ids[2]==ParticleID::g || ids[2]==ParticleID::gamma };
vector<Energy2> pjk(3,ZERO);
vector<Energy> Ek(3,ZERO);
Energy Ei,Ej;
Energy2 m12(ZERO),m22(ZERO);
InvEnergy2 aziMax(ZERO);
bool softAllowed = ShowerHandler::currentHandler()->evolver()->softCorrelations()&&
(canBeSoft[0] || canBeSoft[1]);
if(softAllowed) {
// find the partner for the soft correlations
tShowerParticlePtr partner=findCorrelationPartner(particle,true,splittingFn()->interactionType());
// remember we want the softer gluon
bool swapOrder = !canBeSoft[1] || (canBeSoft[0] && canBeSoft[1] && z < 0.5);
double zFact = !swapOrder ? (1.-z) : z;
// compute the transforms to the shower reference frame
// first the boost
vector<Lorentz5Momentum> basis = kinematics->getBasis();
Lorentz5Momentum pVect = basis[0], nVect = basis[1];
Boost beta_bb;
if(kinematics->frame()==ShowerKinematics::BackToBack) {
beta_bb = -(pVect + nVect).boostVector();
}
else if(kinematics->frame()==ShowerKinematics::Rest) {
beta_bb = -pVect.boostVector();
}
else
assert(false);
pVect.boost(beta_bb);
nVect.boost(beta_bb);
Axis axis;
if(kinematics->frame()==ShowerKinematics::BackToBack) {
axis = pVect.vect().unit();
}
else if(kinematics->frame()==ShowerKinematics::Rest) {
axis = nVect.vect().unit();
}
else
assert(false);
// and then the rotation
LorentzRotation rot;
if(axis.perp2()>0.) {
double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
rot.rotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
}
else if(axis.z()<0.) {
rot.rotate(Constants::pi,Axis(1.,0.,0.));
}
rot.invert();
pVect *= rot;
nVect *= rot;
// shower parameters
Energy2 pn = pVect*nVect, m2 = pVect.m2();
double alpha0 = particle.showerParameters().alpha;
double beta0 = 0.5/alpha0/pn*
(sqr(particle.dataPtr()->mass())-sqr(alpha0)*m2+sqr(particle.showerParameters().pt));
Lorentz5Momentum qperp0(particle.showerParameters().ptx,
particle.showerParameters().pty,ZERO,ZERO);
assert(partner);
Lorentz5Momentum pj = partner->momentum();
pj.boost(beta_bb);
pj *= rot;
// compute the two phi independent dot products
pik = 0.5*zFact*(sqr(alpha0)*m2 - sqr(particle.showerParameters().pt) + 2.*alpha0*beta0*pn )
+0.5*sqr(pT)/zFact;
Energy2 dot1 = pj*pVect;
Energy2 dot2 = pj*nVect;
Energy2 dot3 = pj*qperp0;
pipj = alpha0*dot1+beta0*dot2+dot3;
// compute the constants for the phi dependent dot product
pjk[0] = zFact*(alpha0*dot1+dot3-0.5*dot2/pn*(alpha0*m2-sqr(particle.showerParameters().pt)/alpha0))
+0.5*sqr(pT)*dot2/pn/zFact/alpha0;
pjk[1] = (pj.x() - dot2/alpha0/pn*qperp0.x())*pT;
pjk[2] = (pj.y() - dot2/alpha0/pn*qperp0.y())*pT;
m12 = sqr(particle.dataPtr()->mass());
m22 = sqr(partner->dataPtr()->mass());
if(swapOrder) {
pjk[1] *= -1.;
pjk[2] *= -1.;
}
Ek[0] = zFact*(alpha0*pVect.t()-0.5*nVect.t()/pn*(alpha0*m2-sqr(particle.showerParameters().pt)/alpha0))
+0.5*sqr(pT)*nVect.t()/pn/zFact/alpha0;
Ek[1] = -nVect.t()/alpha0/pn*qperp0.x()*pT;
Ek[2] = -nVect.t()/alpha0/pn*qperp0.y()*pT;
if(swapOrder) {
Ek[1] *= -1.;
Ek[2] *= -1.;
}
Energy mag2=sqrt(sqr(Ek[1])+sqr(Ek[2]));
Ei = alpha0*pVect.t()+beta0*nVect.t();
Ej = pj.t();
double phi0 = atan2(-pjk[2],-pjk[1]);
if(phi0<0.) phi0 += Constants::twopi;
double phi1 = atan2(-Ek[2],-Ek[1]);
if(phi1<0.) phi1 += Constants::twopi;
double xi_min = pik/Ei/(Ek[0]+mag2), xi_max = pik/Ei/(Ek[0]-mag2), xi_ij = pipj/Ei/Ej;
if(xi_min>xi_max) swap(xi_min,xi_max);
if(xi_min>xi_ij) softAllowed = false;
Energy2 mag = sqrt(sqr(pjk[1])+sqr(pjk[2]));
if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==1) {
aziMax = -m12/sqr(pik) -m22/sqr(pjk[0]+mag) +2.*pipj/pik/(pjk[0]-mag);
}
else if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==2) {
double A = (pipj*Ek[0]- Ej*pik)/Ej/sqr(Ej);
double B = -sqrt(sqr(pipj)*(sqr(Ek[1])+sqr(Ek[2])))/Ej/sqr(Ej);
double C = pjk[0]/sqr(Ej);
double D = -sqrt(sqr(pjk[1])+sqr(pjk[2]))/sqr(Ej);
pair<double,double> minima = softPhiMin(phi0,phi1,A,B,C,D);
aziMax = 0.5/pik/(Ek[0]-mag2)*(Ei-m12*(Ek[0]-mag2)/pik + max(Ej*(A+B*cos(minima.first -phi1))/(C+D*cos(minima.first -phi0)),
Ej*(A+B*cos(minima.second-phi1))/(C+D*cos(minima.second-phi0))));
}
else
assert(false);
}
// if spin correlations
vector<pair<int,Complex> > wgts;
if(ShowerHandler::currentHandler()->evolver()->spinCorrelations()) {
- RhoDMatrix rho = extractRhoMatrix(particle,kinematics,true);
+ RhoDMatrix rho = particle.extractRhoMatrix(kinematics,true);
// calculate the weights
wgts = splittingFn()->generatePhiForward(z,t,ids,rho);
}
else {
wgts = vector<pair<int,Complex> >(1,make_pair(0,1.));
}
// generate the azimuthal angle
double phi,wgt;
static const Complex ii(0.,1.);
unsigned int ntry(0);
double phiMax(0.),wgtMax(0.);
do {
phi = Constants::twopi*UseRandom::rnd();
// first the spin correlations bit (gives 1 if correlations off)
Complex spinWgt = 0.;
for(unsigned int ix=0;ix<wgts.size();++ix) {
if(wgts[ix].first==0)
spinWgt += wgts[ix].second;
else
spinWgt += exp(double(wgts[ix].first)*ii*phi)*wgts[ix].second;
}
wgt = spinWgt.real();
if(wgt-1.>1e-10) {
generator()->log() << "Forward spin weight problem " << wgt << " " << wgt-1.
<< " " << ids[0] << " " << ids[1] << " " << ids[2] << " " << " " << phi << "\n";
generator()->log() << "Weights \n";
for(unsigned int ix=0;ix<wgts.size();++ix)
generator()->log() << wgts[ix].first << " " << wgts[ix].second << "\n";
}
// soft correlations bit
double aziWgt = 1.;
if(softAllowed) {
Energy2 dot = pjk[0]+pjk[1]*cos(phi)+pjk[2]*sin(phi);
Energy Eg = Ek[0]+Ek[1]*cos(phi)+Ek[2]*sin(phi);
if(pipj*Eg>pik*Ej) {
if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==1) {
aziWgt = (-m12/sqr(pik) -m22/sqr(dot) +2.*pipj/pik/dot)/aziMax;
}
else if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==2) {
aziWgt = max(ZERO,0.5/pik/Eg*(Ei-m12*Eg/pik + (pipj*Eg - Ej*pik)/dot)/aziMax);
}
if(aziWgt-1.>1e-10||aziWgt<-1e-10) {
generator()->log() << "Forward soft weight problem " << aziWgt << " " << aziWgt-1.
<< " " << ids[0] << " " << ids[1] << " " << ids[2] << " " << " " << phi << "\n";
}
}
else {
aziWgt = 0.;
}
}
wgt *= aziWgt;
if(wgt>wgtMax) {
phiMax = phi;
wgtMax = wgt;
}
++ntry;
}
while(wgt<UseRandom::rnd()&&ntry<10000);
if(ntry==10000) {
generator()->log() << "Too many tries to generate phi in forward evolution\n";
phi = phiMax;
}
// return the azimuthal angle
return phi;
}
double QTildeSudakov::generatePhiBackward(ShowerParticle & particle,
const IdList & ids,
ShoKinPtr kinematics) {
// no correlations, return flat phi
if(! ShowerHandler::currentHandler()->evolver()->correlations())
return Constants::twopi*UseRandom::rnd();
// get the kinematic variables
double z = kinematics->z();
Energy2 t = (1.-z)*sqr(kinematics->scale())/z;
Energy pT = kinematics->pT();
// if soft correlations
bool softAllowed = ShowerHandler::currentHandler()->evolver()->softCorrelations() &&
(ids[2]==ParticleID::g || ids[2]==ParticleID::gamma);
Energy2 pipj,pik,m12(ZERO),m22(ZERO);
vector<Energy2> pjk(3,ZERO);
Energy Ei,Ej,Ek;
InvEnergy2 aziMax(ZERO);
if(softAllowed) {
// find the partner for the soft correlations
tShowerParticlePtr partner=findCorrelationPartner(particle,false,splittingFn()->interactionType());
double zFact = (1.-z);
// compute the transforms to the shower reference frame
// first the boost
vector<Lorentz5Momentum> basis = kinematics->getBasis();
Lorentz5Momentum pVect = basis[0];
Lorentz5Momentum nVect = basis[1];
assert(kinematics->frame()==ShowerKinematics::BackToBack);
Boost beta_bb = -(pVect + nVect).boostVector();
pVect.boost(beta_bb);
nVect.boost(beta_bb);
Axis axis = pVect.vect().unit();
// and then the rotation
LorentzRotation rot;
if(axis.perp2()>0.) {
double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
rot.rotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
}
else if(axis.z()<0.) {
rot.rotate(Constants::pi,Axis(1.,0.,0.));
}
rot.invert();
pVect *= rot;
nVect *= rot;
// shower parameters
Energy2 pn = pVect*nVect;
Energy2 m2 = pVect.m2();
double alpha0 = particle.x();
double beta0 = -0.5/alpha0/pn*sqr(alpha0)*m2;
Lorentz5Momentum pj = partner->momentum();
pj.boost(beta_bb);
pj *= rot;
double beta2 = 0.5*(1.-zFact)*(sqr(alpha0*zFact/(1.-zFact))*m2+sqr(pT))/alpha0/zFact/pn;
// compute the two phi independent dot products
Energy2 dot1 = pj*pVect;
Energy2 dot2 = pj*nVect;
pipj = alpha0*dot1+beta0*dot2;
pik = alpha0*(alpha0*zFact/(1.-zFact)*m2+pn*(beta2+zFact/(1.-zFact)*beta0));
// compute the constants for the phi dependent dot product
pjk[0] = alpha0*zFact/(1.-zFact)*dot1+beta2*dot2;
pjk[1] = pj.x()*pT;
pjk[2] = pj.y()*pT;
m12 = ZERO;
m22 = sqr(partner->dataPtr()->mass());
Energy2 mag = sqrt(sqr(pjk[1])+sqr(pjk[2]));
if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==1) {
aziMax = -m12/sqr(pik) -m22/sqr(pjk[0]+mag) +2.*pipj/pik/(pjk[0]-mag);
}
else if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==2) {
Ek = alpha0*zFact/(1.-zFact)*pVect.t()+beta2*nVect.t();
Ei = alpha0*pVect.t()+beta0*nVect.t();
Ej = pj.t();
if(pipj*Ek> Ej*pik) {
aziMax = 0.5/pik/Ek*(Ei-m12*Ek/pik + (pipj*Ek- Ej*pik)/(pjk[0]-mag));
}
else {
aziMax = 0.5/pik/Ek*(Ei-m12*Ek/pik);
}
}
else {
assert(ShowerHandler::currentHandler()->evolver()->softCorrelations()==0);
}
}
// if spin correlations
vector<pair<int,Complex> > wgts;
if(ShowerHandler::currentHandler()->evolver()->spinCorrelations()) {
// get the spin density matrix and the mapping
- RhoDMatrix rho = extractRhoMatrix(particle,kinematics,false);
+ RhoDMatrix rho = particle.extractRhoMatrix(kinematics,false);
// get the weights
wgts = splittingFn()->generatePhiBackward(z,t,ids,rho);
}
else {
wgts = vector<pair<int,Complex> >(1,make_pair(0,1.));
}
// generate the azimuthal angle
double phi,wgt;
static const Complex ii(0.,1.);
unsigned int ntry(0);
double phiMax(0.),wgtMax(0.);
do {
phi = Constants::twopi*UseRandom::rnd();
Complex spinWgt = 0.;
for(unsigned int ix=0;ix<wgts.size();++ix) {
if(wgts[ix].first==0)
spinWgt += wgts[ix].second;
else
spinWgt += exp(double(wgts[ix].first)*ii*phi)*wgts[ix].second;
}
wgt = spinWgt.real();
if(wgt-1.>1e-10) {
generator()->log() << "Backward weight problem " << wgt << " " << wgt-1.
<< " " << ids[0] << " " << ids[1] << " " << ids[2] << " " << " " << z << " " << phi << "\n";
generator()->log() << "Weights \n";
for(unsigned int ix=0;ix<wgts.size();++ix)
generator()->log() << wgts[ix].first << " " << wgts[ix].second << "\n";
}
// soft correlations bit
double aziWgt = 1.;
if(softAllowed) {
Energy2 dot = pjk[0]+pjk[1]*cos(phi)+pjk[2]*sin(phi);
if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==1) {
aziWgt = (-m12/sqr(pik) -m22/sqr(dot) +2.*pipj/pik/dot)/aziMax;
}
else if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==2) {
aziWgt = max(ZERO,0.5/pik/Ek*(Ei-m12*Ek/pik + pipj*Ek/dot - Ej*pik/dot)/aziMax);
}
if(aziWgt-1.>1e-10||aziWgt<-1e-10) {
generator()->log() << "Backward soft weight problem " << aziWgt << " " << aziWgt-1.
<< " " << ids[0] << " " << ids[1] << " " << ids[2] << " " << " " << phi << "\n";
}
}
wgt *= aziWgt;
if(wgt>wgtMax) {
phiMax = phi;
wgtMax = wgt;
}
++ntry;
}
while(wgt<UseRandom::rnd()&&ntry<10000);
if(ntry==10000) {
generator()->log() << "Too many tries to generate phi in backward evolution\n";
phi = phiMax;
}
// return the azimuthal angle
return phi;
}
double QTildeSudakov::generatePhiDecay(ShowerParticle & particle,
const IdList & ids,
ShoKinPtr kinematics) {
// only soft correlations in this case
// no correlations, return flat phi
if( !(ShowerHandler::currentHandler()->evolver()->softCorrelations() &&
(ids[2]==ParticleID::g || ids[2]==ParticleID::gamma )))
return Constants::twopi*UseRandom::rnd();
// get the kinematic variables
double z = kinematics->z();
Energy pT = kinematics->pT();
// if soft correlations
// find the partner for the soft correlations
tShowerParticlePtr partner = findCorrelationPartner(particle,true,splittingFn()->interactionType());
double zFact(1.-z);
vector<Lorentz5Momentum> basis = kinematics->getBasis();
// compute the transforms to the shower reference frame
// first the boost
Lorentz5Momentum pVect = basis[0];
Lorentz5Momentum nVect = basis[1];
assert(kinematics->frame()==ShowerKinematics::Rest);
Boost beta_bb = -pVect.boostVector();
pVect.boost(beta_bb);
nVect.boost(beta_bb);
Axis axis = nVect.vect().unit();
// and then the rotation
LorentzRotation rot;
if(axis.perp2()>0.) {
double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
rot.rotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
}
else if(axis.z()<0.) {
rot.rotate(Constants::pi,Axis(1.,0.,0.));
}
rot.invert();
pVect *= rot;
nVect *= rot;
// shower parameters
Energy2 pn = pVect*nVect;
Energy2 m2 = pVect.m2();
double alpha0 = particle.showerParameters().alpha;
double beta0 = 0.5/alpha0/pn*
(sqr(particle.dataPtr()->mass())-sqr(alpha0)*m2+sqr(particle.showerParameters().pt));
Lorentz5Momentum qperp0(particle.showerParameters().ptx,
particle.showerParameters().pty,ZERO,ZERO);
Lorentz5Momentum pj = partner->momentum();
pj.boost(beta_bb);
pj *= rot;
// compute the two phi independent dot products
Energy2 pik = 0.5*zFact*(sqr(alpha0)*m2 - sqr(particle.showerParameters().pt) + 2.*alpha0*beta0*pn )
+0.5*sqr(pT)/zFact;
Energy2 dot1 = pj*pVect;
Energy2 dot2 = pj*nVect;
Energy2 dot3 = pj*qperp0;
Energy2 pipj = alpha0*dot1+beta0*dot2+dot3;
// compute the constants for the phi dependent dot product
vector<Energy2> pjk(3,ZERO);
pjk[0] = zFact*(alpha0*dot1+dot3-0.5*dot2/pn*(alpha0*m2-sqr(particle.showerParameters().pt)/alpha0))
+0.5*sqr(pT)*dot2/pn/zFact/alpha0;
pjk[1] = (pj.x() - dot2/alpha0/pn*qperp0.x())*pT;
pjk[2] = (pj.y() - dot2/alpha0/pn*qperp0.y())*pT;
Energy2 m12 = sqr(particle.dataPtr()->mass());
Energy2 m22 = sqr(partner->dataPtr()->mass());
Energy2 mag = sqrt(sqr(pjk[1])+sqr(pjk[2]));
InvEnergy2 aziMax;
vector<Energy> Ek(3,ZERO);
Energy Ei,Ej;
if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==1) {
aziMax = -m12/sqr(pik) -m22/sqr(pjk[0]+mag) +2.*pipj/pik/(pjk[0]-mag);
}
else if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==2) {
Ek[0] = zFact*(alpha0*pVect.t()+-0.5*nVect.t()/pn*(alpha0*m2-sqr(particle.showerParameters().pt)/alpha0))
+0.5*sqr(pT)*nVect.t()/pn/zFact/alpha0;
Ek[1] = -nVect.t()/alpha0/pn*qperp0.x()*pT;
Ek[2] = -nVect.t()/alpha0/pn*qperp0.y()*pT;
Energy mag2=sqrt(sqr(Ek[1])+sqr(Ek[2]));
Ei = alpha0*pVect.t()+beta0*nVect.t();
Ej = pj.t();
aziMax = 0.5/pik/(Ek[0]-mag2)*(Ei-m12*(Ek[0]-mag2)/pik + pipj*(Ek[0]+mag2)/(pjk[0]-mag) - Ej*pik/(pjk[0]-mag) );
}
else
assert(ShowerHandler::currentHandler()->evolver()->softCorrelations()==0);
// generate the azimuthal angle
double phi,wgt(0.);
unsigned int ntry(0);
double phiMax(0.),wgtMax(0.);
do {
phi = Constants::twopi*UseRandom::rnd();
Energy2 dot = pjk[0]+pjk[1]*cos(phi)+pjk[2]*sin(phi);
if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==1) {
wgt = (-m12/sqr(pik) -m22/sqr(dot) +2.*pipj/pik/dot)/aziMax;
}
else if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==2) {
if(qperp0.m2()==ZERO) {
wgt = 1.;
}
else {
Energy Eg = Ek[0]+Ek[1]*cos(phi)+Ek[2]*sin(phi);
wgt = max(ZERO,0.5/pik/Eg*(Ei-m12*Eg/pik + (pipj*Eg - Ej*pik)/dot)/aziMax);
}
}
if(wgt-1.>1e-10||wgt<-1e-10) {
generator()->log() << "Decay soft weight problem " << wgt << " " << wgt-1.
<< " " << ids[0] << " " << ids[1] << " " << ids[2] << " " << " " << phi << "\n";
}
if(wgt>wgtMax) {
phiMax = phi;
wgtMax = wgt;
}
++ntry;
}
while(wgt<UseRandom::rnd()&&ntry<10000);
if(ntry==10000) {
phi = phiMax;
generator()->log() << "Too many tries to generate phi\n";
}
// return the azimuthal angle
return phi;
}
Energy QTildeSudakov::calculateScale(double zin, Energy pt, IdList ids,
unsigned int iopt) {
Energy2 tmin;
initialize(ids,tmin,false);
// final-state branching
if(iopt==0) {
Energy2 scale=(sqr(pt)+masssquared_[1]*(1.-zin)+masssquared_[2]*zin);
if(ids[0]!=ParticleID::g) scale -= zin*(1.-zin)*masssquared_[0];
scale /= sqr(zin*(1-zin));
return scale<=ZERO ? sqrt(tmin) : sqrt(scale);
}
else if(iopt==1) {
Energy2 scale=(sqr(pt)+zin*masssquared_[2])/sqr(1.-zin);
return scale<=ZERO ? sqrt(tmin) : sqrt(scale);
}
else if(iopt==2) {
Energy2 scale = (sqr(pt)+zin*masssquared_[2])/sqr(1.-zin)+masssquared_[0];
return scale<=ZERO ? sqrt(tmin) : sqrt(scale);
}
else {
throw Exception() << "Unknown option in QTildeSudakov::calculateScale() "
<< "iopt = " << iopt << Exception::runerror;
}
}
ShoKinPtr QTildeSudakov::createFinalStateBranching(Energy scale,double z,
double phi, Energy pt) {
ShoKinPtr showerKin = new_ptr(FS_QTildeShowerKinematics1to2());
showerKin->scale(scale);
showerKin->z(z);
showerKin->phi(phi);
showerKin->pT(pt);
showerKin->SudakovFormFactor(this);
return showerKin;
}
ShoKinPtr QTildeSudakov::createInitialStateBranching(Energy scale,double z,
double phi, Energy pt) {
ShoKinPtr showerKin = new_ptr(IS_QTildeShowerKinematics1to2());
showerKin->scale(scale);
showerKin->z(z);
showerKin->phi(phi);
showerKin->pT(pt);
showerKin->SudakovFormFactor(this);
return showerKin;
}
ShoKinPtr QTildeSudakov::createDecayBranching(Energy scale,double z,
double phi, Energy pt) {
ShoKinPtr showerKin = new_ptr(Decay_QTildeShowerKinematics1to2());
showerKin->scale(scale);
showerKin->z(z);
showerKin->phi(phi);
showerKin->pT(pt);
showerKin->SudakovFormFactor(this);
return showerKin;
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 5:52 PM (9 h, 26 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4015453
Default Alt Text
(102 KB)

Event Timeline