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 using namespace Herwig; PPtr ShowerParticle::clone() const { return new_ptr(*this); } PPtr ShowerParticle::fullclone() const { return new_ptr(*this); } ClassDescription 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 & 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 sbasis; + for(unsigned int ix=0;ix<3;++ix) { + sbasis.push_back(vspin->getProductionBasisState(ix)); + sbasis.back().transform(rot); + } + // splitting basis + vector 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 > 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 > 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()) { + 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 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;ixaperturbative()) { + // 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(this->spinInfo()); + if(!sspin) { + ScalarWaveFunction::constructSpinInfo(this,outgoing,true); + } + output=this->spinInfo(); + return false; + } + else if(spin==PDT::Spin1Half) { + FermionSpinPtr fspin=dynamic_ptr_cast(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(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(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(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(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(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 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 & 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 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 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 { /** 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 : public ClassTraitsBase { /** 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 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 documentation ("The SudakovFormFactor class is the base class for the implementation of Sudakov" " form factors in Herwig"); static Reference interfaceSplittingFunction("SplittingFunction", "A reference to the SplittingFunction object", &Herwig::SudakovFormFactor::splittingFn_, false, false, true, false); static Reference interfaceAlpha("Alpha", "A reference to the Alpha object", &Herwig::SudakovFormFactor::alpha_, false, false, true, false); static Parameter interfacePDFmax ("PDFmax", "Maximum value of PDF weight. ", &SudakovFormFactor::pdfmax_, 35.0, 1.0, 100000.0, false, false, Interface::limited); static Switch 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 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 interfaceaParameter ("aParameter", "The a parameter for the kinematic cut-off", &SudakovFormFactor::a_, 0.3, -10.0, 10.0, false, false, Interface::limited); static Parameter interfacebParameter ("bParameter", "The b parameter for the kinematic cut-off", &SudakovFormFactor::b_, 2.3, -10.0, 10.0, false, false, Interface::limited); static Parameter 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 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 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 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 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::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 & 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 sbasis; - for(unsigned int ix=0;ix<3;++ix) { - sbasis.push_back(vspin->getProductionBasisState(ix)); - sbasis.back().transform(rot); - } - // splitting basis - vector 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 > 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 > 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()) { - 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 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(particle.spinInfo()); - if(!sspin) { - ScalarWaveFunction::constructSpinInfo(&particle,outgoing,true); - } - output=particle.spinInfo(); - return false; - } - else if(spin==PDT::Spin1Half) { - FermionSpinPtr fspin=dynamic_ptr_cast(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(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(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(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(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(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::iterator it=particles_.begin(); it!=particles_.end();++it) { if(it->size()==in.size()) { bool match=true; for(unsigned int iy=0;iy::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)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 & SudakovFormFactor::virtualMasses(const IdList & ids) { static vector output; output.clear(); if(cutOffOption() == 0) { for(unsigned int ix=0;ixmass()); Energy kinCutoff= kinematicCutOff(kinScale(),*std::max_element(output.begin(),output.end())); for(unsigned int ix=0;ixmass()); output.back() += ids[ix]==ParticleID::g ? vgCut() : vqCut(); } } else if(cutOffOption() == 2) { for(unsigned int ix=0;ixmass()); } 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::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 zLimits() const { return zlimits_;} /** * Set the limits */ void zLimits(pair 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 & 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 & 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 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 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 describeQTildeSudakov ("Herwig::QTildeSudakov","HwShower.so"); void QTildeSudakov::Init() { static ClassDocumentation 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(tupper 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 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 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::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;ixCC()) ids_[ix]*=-1; } } tmin = cutOffOption() != 2 ? ZERO : 4.*pT2min(); masses_ = virtualMasses(ids); masssquared_.clear(); for(unsigned int ix=0;ix0) 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 ) || pt2masssquared_[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 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().secondtmax||zLimits().second 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 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(particle.parents()[0]) : tShowerParticlePtr(); } else { mother = particle.children().size()==2 ? dynamic_ptr_cast(&particle) : tShowerParticlePtr(); } tShowerParticlePtr partner; while(mother) { tPPtr otherChild; if(forward) { for (unsigned int ix=0;ixchildren().size();++ix) { if(mother->children()[ix]!=child) { otherChild = mother->children()[ix]; break; } } } else { otherChild = mother->children()[1]; } tShowerParticlePtr other = dynamic_ptr_cast(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(mother); break; } child = mother; if(forward) { mother = ! mother->parents().empty() ? dynamic_ptr_cast(mother->parents()[0]) : tShowerParticlePtr(); } else { if(mother->children()[0]->children().size()!=2) break; tShowerParticlePtr mtemp = dynamic_ptr_cast(mother->children()[0]); if(!mtemp) break; else mother=mtemp; } } if(!partner) { if(forward) { partner = dynamic_ptr_cast( child)->partner(); } else { if(mother) { tShowerParticlePtr parent; if(!mother->children().empty()) { parent = dynamic_ptr_cast(mother->children()[0]); } if(!parent) { parent = dynamic_ptr_cast(mother); } partner = parent->partner(); } else { partner = dynamic_ptr_cast(&particle)->partner(); } } } return partner; } pair 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 pjk(3,ZERO); vector 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 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 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 > 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 >(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;ix1e-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;ixlog() << 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(wgtlog() << "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 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 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 > 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 >(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;ix1e-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;ixlog() << 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(wgtlog() << "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 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 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 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(wgtlog() << "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; }