Page MenuHomeHEPForge
Files F8724420

No OneTemporary
Actions

View Options

diff --git a/Shower/Base/SudakovFormFactor.cc b/Shower/Base/SudakovFormFactor.cc
--- a/Shower/Base/SudakovFormFactor.cc
+++ b/Shower/Base/SudakovFormFactor.cc
@@ -1,595 +1,670 @@
// -*- 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/Utilities/DescribeClass.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);
}
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);
}
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::
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;
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) {
- // 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 porig(basis[0]);
- // construct the Lorentz boost
- LorentzRotation output(boostv);
- porig *= output;
- Axis axis(porig.vect().unit());
- // now rotate so along the z axis as needed for the splitting functions
- if(axis.perp2()>0.) {
- double sinth(sqrt(1.-sqr(axis.z())));
- output.rotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
+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()>0.) {
+ 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 if(axis.z()<0.) {
- output.rotate(Constants::pi,Axis(1.,0.,0.));
+ 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()==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());
- output=particle.spinInfo();
- assert(output);
- return false;
+ if(output) {
+ return false;
+ }
+ else {
+ assert(particle.parents().size()==1);
+ Lorentz5Momentum porig;
+ LorentzRotation rot = boostToShower(showerkin->getBasis(),showerkin->frame(),porig);
+ if(particle.dataPtr()->iSpin()==PDT::Spin0) {
+ assert(false);
+ }
+ else if(particle.dataPtr()->iSpin()==PDT::Spin1Half) {
+ output = createFermionSpinInfo(particle,porig,rot,outgoing);
+ }
+ else if(particle.dataPtr()->iSpin()==PDT::Spin1) {
+ output = createVectorSpinInfo(particle,porig,rot,outgoing);
+ }
+ 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 the basis vectors
- vector<Lorentz5Momentum> basis=showerkin->getBasis();
// get transform to shower frame
- LorentzRotation rot = boostToShower(basis);
- Lorentz5Momentum porig = rot*basis[0];
- porig.setX(ZERO);
- porig.setY(ZERO);
+ 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) {
cerr << "testing spin 0 not yet implemented " << endl;
assert(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);
}
- else if(particle.perturbative() && !particle.isFinalState()) {
- assert(particle.perturbative()==1);
+ // incoming to hard process
+ else if(particle.perturbative()==1 && !particle.isFinalState()) {
// get the basis vectors
- vector<Lorentz5Momentum> basis=showerkin->getBasis();
// get transform to shower frame
- LorentzRotation rot = boostToShower(basis);
- Lorentz5Momentum porig = basis[0]*particle.x();
- porig *= rot;
- porig.setX(ZERO);
- porig.setY(ZERO);
+ 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_;
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;
}
diff --git a/Shower/Base/SudakovFormFactor.h b/Shower/Base/SudakovFormFactor.h
--- a/Shower/Base/SudakovFormFactor.h
+++ b/Shower/Base/SudakovFormFactor.h
@@ -1,658 +1,667 @@
// -*- 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_() {}
/**
* 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
* defined.
*/
virtual ShoKinPtr generateNextTimeBranching(const Energy startingScale,
const IdList &ids,const bool cc,
double enhance)=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(const Energy2 pt2) const {
return UseRandom::rnd() > ThePEG::Math::powi(alpha_->ratio(pt2),
splittingFn_->interactionOrder());
}
//@}
/**
* 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);
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_;
//@}
};
}
#endif /* HERWIG_SudakovFormFactor_H */
diff --git a/Shower/Default/FS_QTildeShowerKinematics1to2.cc b/Shower/Default/FS_QTildeShowerKinematics1to2.cc
--- a/Shower/Default/FS_QTildeShowerKinematics1to2.cc
+++ b/Shower/Default/FS_QTildeShowerKinematics1to2.cc
@@ -1,204 +1,215 @@
// -*- C++ -*-
//
// FS_QTildeShowerKinematics1to2.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 FS_QTildeShowerKinematics1to2 class.
//
#include "FS_QTildeShowerKinematics1to2.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "Herwig++/Shower/SplittingFunctions/SplittingFunction.h"
#include "Herwig++/Shower/Base/ShowerParticle.h"
#include "ThePEG/Utilities/Debug.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"
+#include "Herwig++/Shower/Base/ShowerVertex.h"
using namespace Herwig;
void FS_QTildeShowerKinematics1to2::
updateParameters(tShowerParticlePtr theParent,
tShowerParticlePtr theChild0,
tShowerParticlePtr theChild1,
bool setAlpha) const {
const ShowerParticle::Parameters & parent = theParent->showerParameters();
ShowerParticle::Parameters & child0 = theChild0->showerParameters();
ShowerParticle::Parameters & child1 = theChild1->showerParameters();
// determine alphas of children according to interpretation of z
if ( setAlpha ) {
child0.alpha = z() * parent.alpha;
child1.alpha = (1.-z()) * parent.alpha;
}
// set the values
double cphi = cos(phi());
double sphi = sin(phi());
child0.ptx = pT() * cphi + z() * parent.ptx;
child0.pty = pT() * sphi + z() * parent.pty;
child0.pt = sqrt( sqr(child0.ptx) + sqr(child0.pty) );
child1.ptx = -pT() * cphi + (1.-z())* parent.ptx;
child1.pty = -pT() * sphi + (1.-z())* parent.pty;
child1.pt = sqrt( sqr(child1.ptx) + sqr(child1.pty) );
}
void FS_QTildeShowerKinematics1to2::
updateChildren(const tShowerParticlePtr parent,
const ShowerParticleVector & children,
ShowerPartnerType::Type partnerType) const {
assert(children.size()==2);
// calculate the scales
splittingFn()->evaluateFinalStateScales(partnerType,scale(),z(),parent,
children[0],children[1]);
// update the parameters
updateParameters(parent, children[0], children[1], true);
// set up the colour connections
splittingFn()->colourConnection(parent,children[0],children[1],partnerType,false);
// make the products children of the parent
parent->addChild(children[0]);
parent->addChild(children[1]);
// sort out the helicity stuff
if(! ShowerHandler::currentHandler()->evolver()->correlations()) return;
SpinPtr pspin(parent->spinInfo());
if(!pspin) return;
- // get the vertex
- VertexPtr vertex(const_ptr_cast<VertexPtr>(pspin->decayVertex()));
- if(!vertex) return;
+ Energy2 t = sqr(scale())*z()*(1.-z());
+ IdList ids;
+ ids.push_back(parent->id());
+ ids.push_back(children[0]->id());
+ ids.push_back(children[1]->id());
+ // compute the matrix element for spin correlations
+ DecayMatrixElement me(splittingFn()->matrixElement(z(),t,ids,phi()));
+ // create the vertex
+ SVertexPtr vertex(new_ptr(ShowerVertex()));
+ // set the matrix element
+ vertex->ME().reset(me);
+ // set the incoming particle for the vertex
+ parent->spinInfo()->decayVertex(vertex);
ShowerParticleVector::const_iterator pit;
for(pit=children.begin();pit!=children.end();++pit) {
// construct the spin info for the children
constructSpinInfo(*pit,true);
// connect the spinInfo object to the vertex
(*pit)->spinInfo()->productionVertex(vertex);
}
}
void FS_QTildeShowerKinematics1to2::
reconstructParent(const tShowerParticlePtr parent,
const ParticleVector & children ) const {
assert(children.size() == 2);
ShowerParticlePtr c1 = dynamic_ptr_cast<ShowerParticlePtr>(children[0]);
ShowerParticlePtr c2 = dynamic_ptr_cast<ShowerParticlePtr>(children[1]);
parent->showerParameters().beta=
c1->showerParameters().beta + c2->showerParameters().beta;
parent->set5Momentum( c1->momentum() + c2->momentum() );
}
void FS_QTildeShowerKinematics1to2::reconstructLast(const tShowerParticlePtr theLast,
Energy mass) const {
// set beta component and consequently all missing data from that,
// using the nominal (i.e. PDT) mass.
Energy theMass = mass > ZERO ? mass : theLast->data().constituentMass();
ShowerParticle::Parameters & last = theLast->showerParameters();
last.beta = ( sqr(theMass) + sqr(last.pt) - sqr(last.alpha) * pVector().m2() )
/ ( 2. * last.alpha * p_dot_n() );
// set that new momentum
theLast->set5Momentum(sudakov2Momentum( last.alpha, last.beta,
last.ptx, last.pty) );
}
void FS_QTildeShowerKinematics1to2::initialize(ShowerParticle & particle,PPtr) {
// set the basis vectors
Lorentz5Momentum p,n;
Frame frame;
if(particle.perturbative()!=0) {
// find the partner and its momentum
ShowerParticlePtr partner=particle.partner();
Lorentz5Momentum ppartner(partner->momentum());
// momentum of the emitting particle
p = particle.momentum();
Lorentz5Momentum pcm;
// if the partner is a final-state particle then the reference
// vector is along the partner in the rest frame of the pair
if(partner->isFinalState()) {
Boost boost=(p + ppartner).findBoostToCM();
pcm = ppartner;
pcm.boost(boost);
n = Lorentz5Momentum(ZERO,pcm.vect());
n.boost( -boost);
}
else if(!partner->isFinalState()) {
// if the partner is an initial-state particle then the reference
// vector is along the partner which should be massless
if(particle.perturbative()==1)
{n = Lorentz5Momentum(ZERO,ppartner.vect());}
// if the partner is an initial-state decaying particle then the reference
// vector is along the backwards direction in rest frame of decaying particle
else {
Boost boost=ppartner.findBoostToCM();
pcm = p;
pcm.boost(boost);
n = Lorentz5Momentum( ZERO, -pcm.vect());
n.boost( -boost);
}
}
frame = BackToBack;
}
else if(particle.initiatesTLS()) {
tShoKinPtr kin=dynamic_ptr_cast<ShowerParticlePtr>
(particle.parents()[0]->children()[0])->showerKinematics();
p = kin->getBasis()[0];
n = kin->getBasis()[1];
frame = kin->frame();
}
else {
tShoKinPtr kin=dynamic_ptr_cast<ShowerParticlePtr>(particle.parents()[0])
->showerKinematics();
p = kin->getBasis()[0];
n = kin->getBasis()[1];
frame = kin->frame();
}
// set the basis vectors
setBasis(p,n,frame);
}
void FS_QTildeShowerKinematics1to2::updateParent(const tShowerParticlePtr parent,
const ShowerParticleVector & children,
ShowerPartnerType::Type) const {
IdList ids(3);
ids[0] = parent->id();
ids[1] = children[0]->id();
ids[2] = children[1]->id();
const vector<Energy> & virtualMasses = SudakovFormFactor()->virtualMasses(ids);
if(children[0]->children().empty()) children[0]->virtualMass(virtualMasses[1]);
if(children[1]->children().empty()) children[1]->virtualMass(virtualMasses[2]);
// compute the new pT of the branching
Energy2 pt2=sqr(z()*(1.-z()))*sqr(scale())
- sqr(children[0]->virtualMass())*(1.-z())
- sqr(children[1]->virtualMass())* z() ;
if(ids[0]!=ParticleID::g) pt2 += z()*(1.-z())*sqr(virtualMasses[0]);
Energy2 q2 =
sqr(children[0]->virtualMass())/z() +
sqr(children[1]->virtualMass())/(1.-z()) +
pt2/z()/(1.-z());
if(pt2<ZERO) {
parent->virtualMass(ZERO);
}
else {
parent->virtualMass(sqrt(q2));
pT(sqrt(pt2));
}
}
void FS_QTildeShowerKinematics1to2::
resetChildren(const tShowerParticlePtr parent,
const ShowerParticleVector & children) const {
updateParameters(parent, children[0], children[1], false);
for(unsigned int ix=0;ix<children.size();++ix) {
if(children[ix]->children().empty()) continue;
ShowerParticleVector newChildren;
for(unsigned int iy=0;iy<children[ix]->children().size();++iy)
newChildren.push_back(dynamic_ptr_cast<ShowerParticlePtr>
(children[ix]->children()[iy]));
children[ix]->showerKinematics()->resetChildren(children[ix],newChildren);
}
}
diff --git a/Shower/Default/IS_QTildeShowerKinematics1to2.cc b/Shower/Default/IS_QTildeShowerKinematics1to2.cc
--- a/Shower/Default/IS_QTildeShowerKinematics1to2.cc
+++ b/Shower/Default/IS_QTildeShowerKinematics1to2.cc
@@ -1,179 +1,191 @@
// -*- C++ -*-
//
// IS_QTildeShowerKinematics1to2.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 IS_QTildeShowerKinematics1to2 class.
//
#include "IS_QTildeShowerKinematics1to2.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "Herwig++/Shower/Base/ShowerParticle.h"
#include "ThePEG/Utilities/Debug.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"
+#include "Herwig++/Shower/Base/ShowerVertex.h"
#include <cassert>
using namespace Herwig;
void IS_QTildeShowerKinematics1to2::
updateChildren( const tShowerParticlePtr theParent,
const ShowerParticleVector & theChildren,
ShowerPartnerType::Type) const {
const ShowerParticle::Parameters & parent = theParent->showerParameters();
ShowerParticle::Parameters & child0 = theChildren[0]->showerParameters();
ShowerParticle::Parameters & child1 = theChildren[1]->showerParameters();
double cphi = cos(phi());
double sphi = sin(phi());
child1.alpha = (1.-z()) * parent.alpha;
child1.ptx = (1.-z()) * parent.ptx - cphi * pT();
child1.pty = (1.-z()) * parent.pty - sphi * pT();
child1.pt = sqrt( sqr(child1.ptx) + sqr(child1.pty) );
// space-like child
child0.alpha = parent.alpha - child1.alpha;
child0.beta = parent.beta - child1.beta;
child0.ptx = parent.ptx - child1.ptx;
child0.pty = parent.pty - child1.pty;
}
void IS_QTildeShowerKinematics1to2::
updateParent(const tShowerParticlePtr parent,
const ShowerParticleVector & children,
ShowerPartnerType::Type partnerType) const {
// calculate the scales
splittingFn()->evaluateInitialStateScales(partnerType,scale(),z(),parent,
children[0],children[1]);
// set proper colour connections
splittingFn()->colourConnection(parent,children[0],children[1],
partnerType,true);
// set proper parent/child relationships
parent->addChild(children[0]);
parent->addChild(children[1]);
parent->x(children[0]->x()/z());
// sort out the helicity stuff
if(! ShowerHandler::currentHandler()->evolver()->correlations()) return;
SpinPtr pspin(children[0]->spinInfo());
if(!pspin) return;
- // get the vertex
- VertexPtr vertex(const_ptr_cast<VertexPtr>(pspin->decayVertex()));
- if(!vertex) return;
+ // compute the matrix element for spin correlations
+ IdList ids;
+ ids.push_back(parent->id());
+ ids.push_back(children[0]->id());
+ ids.push_back(children[1]->id());
+ Energy2 t = (1.-z())*sqr(scale())/z();
+ DecayMatrixElement me(splittingFn()->matrixElement(z(),t,ids,phi()));
+ // create the vertex
+ SVertexPtr vertex(new_ptr(ShowerVertex()));
+ // set the matrix element
+ vertex->ME().reset(me);
+ // set the incoming particle for the vertex
+ // (in reality the first child as going backwards)
+ pspin->decayVertex(vertex);
// construct the spin info for parent and timelike child
// temporary assignment of shower parameters to calculate correlations
parent->showerParameters().alpha = parent->x();
children[1]->showerParameters().alpha = (1.-z()) * parent->x();
children[1]->showerParameters().ptx = - cos(phi()) * pT();
children[1]->showerParameters().pty = - sin(phi()) * pT();
children[1]->showerParameters().pt = pT();
// construct the spin infos
constructSpinInfo(parent,false);
constructSpinInfo(children[1],true);
// connect the spinInfo objects to the vertex
parent ->spinInfo()->productionVertex(vertex);
children[1]->spinInfo()->productionVertex(vertex);
}
void IS_QTildeShowerKinematics1to2::
reconstructParent(const tShowerParticlePtr theParent,
const ParticleVector & theChildren ) const {
PPtr c1 = theChildren[0];
ShowerParticlePtr c2 = dynamic_ptr_cast<ShowerParticlePtr>(theChildren[1]);
ShowerParticle::Parameters & c2param = c2->showerParameters();
// get shower variables from 1st child in order to keep notation
// parent->(c1, c2) clean even though the splitting was initiated
// from c1. The name updateParent is still referring to the
// timelike branching though.
// on-shell child
c2param.beta = 0.5*( sqr(c2->data().constituentMass()) + sqr(c2param.pt) )
/ ( c2param.alpha * p_dot_n() );
c2->set5Momentum( sudakov2Momentum(c2param.alpha, c2param.beta,
c2param.ptx , c2param.pty) );
// spacelike child
Lorentz5Momentum pc1(theParent->momentum() - c2->momentum());
pc1.rescaleMass();
c1->set5Momentum(pc1);
}
void IS_QTildeShowerKinematics1to2::
updateLast( const tShowerParticlePtr theLast,Energy px,Energy py) const {
if(theLast->isFinalState()) return;
ShowerParticle::Parameters & last = theLast->showerParameters();
Energy2 pt2 = sqr(px) + sqr(py);
last.alpha = theLast->x();
last.beta = 0.5 * pt2 / last.alpha / p_dot_n();
last.ptx = ZERO;
last.pty = ZERO;
last.pt = ZERO;
// momentum
Lorentz5Momentum ntemp = Lorentz5Momentum(ZERO,-pVector().vect());
double beta = 0.5 * pt2 / last.alpha / (pVector() * ntemp);
Lorentz5Momentum plast =
Lorentz5Momentum( (pVector().z()>ZERO ? px : -px), py, ZERO, ZERO)
+ theLast->x() * pVector() + beta * ntemp;
plast.rescaleMass();
theLast->set5Momentum(plast);
}
void IS_QTildeShowerKinematics1to2::initialize(ShowerParticle & particle, PPtr parent) {
// For the time being we are considering only 1->2 branching
Lorentz5Momentum p, n, pthis, pcm;
assert(particle.perturbative()!=2);
Frame frame;
if(particle.perturbative()==1) {
// find the partner and its momentum
ShowerParticlePtr partner=particle.partner();
assert(partner);
if(partner->isFinalState()) {
Lorentz5Momentum pa = -particle.momentum()+partner->momentum();
Lorentz5Momentum pb = particle.momentum();
Energy scale=parent->momentum().t();
Lorentz5Momentum pbasis(ZERO,parent->momentum().vect().unit()*scale);
Axis axis(pa.vect().unit());
LorentzRotation rot;
double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
if(axis.perp2()>1e-20) {
rot.setRotate(-acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
rot.rotateX(Constants::pi);
}
if(abs(1.-pa.e()/pa.vect().mag())>1e-6) rot.boostZ( pa.e()/pa.vect().mag());
pb *= rot;
if(pb.perp2()/GeV2>1e-20) {
Boost trans = -1./pb.e()*pb.vect();
trans.setZ(0.);
rot.boost(trans);
}
pbasis *=rot;
rot.invert();
n = rot*Lorentz5Momentum(ZERO,-pbasis.vect());
p = rot*Lorentz5Momentum(ZERO, pbasis.vect());
}
else {
pcm = parent->momentum();
p = Lorentz5Momentum(ZERO, pcm.vect());
n = Lorentz5Momentum(ZERO, -pcm.vect());
}
frame = BackToBack;
}
else {
p = dynamic_ptr_cast<ShowerParticlePtr>(particle.children()[0])
->showerKinematics()->getBasis()[0];
n = dynamic_ptr_cast<ShowerParticlePtr>(particle.children()[0])
->showerKinematics()->getBasis()[1];
frame = dynamic_ptr_cast<ShowerParticlePtr>(particle.children()[0])
->showerKinematics()->frame();
}
setBasis(p,n,frame);
}
diff --git a/Shower/Default/QTildeSudakov.cc b/Shower/Default/QTildeSudakov.cc
--- a/Shower/Default/QTildeSudakov.cc
+++ b/Shower/Default/QTildeSudakov.cc
@@ -1,509 +1,585 @@
// -*- 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/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) {
// still inside PS, return true if outside
// check vs overestimated limits
if(z() < zLimits().first || z() > zLimits().second) return true;
// compute the pts
Energy2 pt2=sqr(z()*(1.-z()))*t-masssquared_[1]*(1.-z())-masssquared_[2]*z();
if(ids_[0]!=ParticleID::g) pt2+=z()*(1.-z())*masssquared_[0];
// 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) {
// 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) || SplittingFnVeto(z()*(1.-z())*t,ids,true) ||
alphaSVeto(sqr(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(sqr(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(sqr(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(Constants::twopi*UseRandom::rnd());
+ 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]) {
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;
}
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 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=inspin->rhoMatrix();
// 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;
}
// get the kinematic variables
double z = kinematics->z();
Energy2 t = z*(1.-z)*sqr(kinematics->scale());
// generate the azimuthal angle
double phi;
Complex wgt;
vector<pair<int,Complex> >
wgts = splittingFn()->generatePhiForward(z,t,ids,rho);
static const Complex ii(0.,1.);
do {
phi = Constants::twopi*UseRandom::rnd();
wgt = 0.;
for(unsigned int ix=0;ix<wgts.size();++ix) {
if(wgts[ix].first==0)
wgt += wgts[ix].second;
else
wgt += exp(double(wgts[ix].first)*ii*phi)*wgts[ix].second;
}
if(wgt.real()-1.>1e-10) {
cerr << "Forward weight problem " << wgt << " " << wgt.real()-1.
<< " " << ids[0] << " " << ids[1] << " " << ids[2] << " " << " " << z << " " << phi << "\n";
cerr << "Weights \n";
for(unsigned int ix=0;ix<wgts.size();++ix)
cerr << wgts[ix].first << " " << wgts[ix].second << "\n";
}
}
while(wgt.real()<UseRandom::rnd());
- // compute the matrix element for spin correlations
- DecayMatrixElement me(splittingFn()->matrixElement(z,t,ids,phi));
- // create the vertex
- SVertexPtr Svertex(new_ptr(ShowerVertex()));
- // set the matrix element
- Svertex->ME().reset(me);
- // set the incoming particle for the vertex
- inspin->decayVertex(Svertex);
// 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 spin density matrix and the mapping
RhoDMatrix mapping;
SpinPtr inspin;
bool needMapping = getMapping(inspin,mapping,particle,kinematics);
// set the decayed flag (counterintuitive but going backward)
inspin->decay();
// get the spin density matrix
RhoDMatrix rho=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;
}
// get the kinematic variables
double z = kinematics->z();
Energy2 t = (1.-z)*sqr(kinematics->scale())/z;
// generate the azimuthal angle
double phi;
Complex wgt;
vector<pair<int,Complex> >
wgts = splittingFn()->generatePhiBackward(z,t,ids,rho);
static const Complex ii(0.,1.);
do {
phi = Constants::twopi*UseRandom::rnd();
wgt = 0.;
for(unsigned int ix=0;ix<wgts.size();++ix) {
if(wgts[ix].first==0)
wgt += wgts[ix].second;
else
wgt += exp(double(wgts[ix].first)*ii*phi)*wgts[ix].second;
}
if(wgt.real()-1.>1e-10) {
cerr << "Backward weight problem " << wgt << " " << wgt.real()-1.
<< " " << ids[0] << " " << ids[1] << " " << ids[2] << " " << " " << z << " " << phi << "\n";
cerr << "Weights \n";
for(unsigned int ix=0;ix<wgts.size();++ix)
cerr << wgts[ix].first << " " << wgts[ix].second << "\n";
}
}
while(wgt.real()<UseRandom::rnd());
- // compute the matrix element for spin correlations
- DecayMatrixElement me(splittingFn()->matrixElement(z,t,ids,phi));
- // create the vertex
- SVertexPtr Svertex(new_ptr(ShowerVertex()));
- // set the matrix element
- Svertex->ME().reset(me);
- // set the incoming particle for the vertex
- // (in reality the first child as going backwards)
- inspin->decayVertex(Svertex);
// return the azimuthal angle
return phi;
}
+double QTildeSudakov::generatePhiDecay(ShowerParticle & particle,
+ const IdList & ids,
+ ShoKinPtr kinematics) {
+ return Constants::twopi*UseRandom::rnd();
+ // cerr << particle.isFinalState() << " " << particle << "\n";
+ // cerr << particle.spinInfo() << "\n";
+ // if(particle.spinInfo()) {
+ // cerr << "testing spin info " << particle.spinInfo()->productionVertex() << " "
+ // << particle.spinInfo()->decayVertex() << " "
+ // << particle.spinInfo()->developed() << " " << particle.spinInfo()->decayed() << "\n";
+ // if(particle.spinInfo()->productionVertex()) {
+ // cerr << "production "
+ // << particle.spinInfo()->productionVertex()->incoming().size() << " "
+ // << particle.spinInfo()->productionVertex()->outgoing().size() << "\n";
+ // }
+ // if(particle.spinInfo()->decayVertex()) {
+ // cerr << "decay "
+ // << particle.spinInfo()->decayVertex()->incoming().size() << " "
+ // << particle.spinInfo()->decayVertex()->outgoing().size() << "\n";
+ // }
+ // }
+
+ // cerr << "testing variables " << kinematics->z() << " " << kinematics->scale()/GeV << " "
+ // << particle.mass()/GeV << "\n";;
+
+ // // no correlations, return flat phi
+ // if(! ShowerHandler::currentHandler()->evolver()->correlations())
+ // return Constants::twopi*UseRandom::rnd();
+ // // get the spin density matrix and the mapping
+ // RhoDMatrix mapping;
+ // SpinPtr inspin;
+ // bool needMapping = getMapping(inspin,mapping,particle,kinematics);
+ // assert(false);
+ //
+ // // set the decayed flag
+ // inspin->decay();
+ // // get the spin density matrix
+ // RhoDMatrix rho=inspin->rhoMatrix();
+ // // 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;
+ // }
+ // // get the kinematic variables
+ // double z = kinematics->z();
+ // Energy2 t = z*(1.-z)*sqr(kinematics->scale());
+ // // generate the azimuthal angle
+ // double phi;
+ // Complex wgt;
+ // vector<pair<int,Complex> >
+ // wgts = splittingFn()->generatePhiForward(z,t,ids,rho);
+ // static const Complex ii(0.,1.);
+ // do {
+ // phi = Constants::twopi*UseRandom::rnd();
+ // wgt = 0.;
+ // for(unsigned int ix=0;ix<wgts.size();++ix) {
+ // if(wgts[ix].first==0)
+ // wgt += wgts[ix].second;
+ // else
+ // wgt += exp(double(wgts[ix].first)*ii*phi)*wgts[ix].second;
+ // }
+ // if(wgt.real()-1.>1e-10) {
+ // cerr << "Forward weight problem " << wgt << " " << wgt.real()-1.
+ // << " " << ids[0] << " " << ids[1] << " " << ids[2] << " " << " " << z << " " << phi << "\n";
+ // cerr << "Weights \n";
+ // for(unsigned int ix=0;ix<wgts.size();++ix)
+ // cerr << wgts[ix].first << " " << wgts[ix].second << "\n";
+ // }
+ // }
+ // while(wgt.real()<UseRandom::rnd());
+ // // compute the matrix element for spin correlations
+ // DecayMatrixElement me(splittingFn()->matrixElement(z,t,ids,phi));
+ // // create the vertex
+ // SVertexPtr Svertex(new_ptr(ShowerVertex()));
+ // // set the matrix element
+ // Svertex->ME().reset(me);
+ // // set the incoming particle for the vertex
+ // inspin->decayVertex(Svertex);
+ // // 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;
}
diff --git a/Shower/Default/QTildeSudakov.h b/Shower/Default/QTildeSudakov.h
--- a/Shower/Default/QTildeSudakov.h
+++ b/Shower/Default/QTildeSudakov.h
@@ -1,275 +1,285 @@
// -*- C++ -*-
//
// QTildeSudakov.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_QTildeSudakov_H
#define HERWIG_QTildeSudakov_H
//
// This is the declaration of the QTildeSudakov class.
//
#include "Herwig++/Shower/Base/SudakovFormFactor.h"
namespace Herwig {
using namespace ThePEG;
/** \ingroup Shower
*
* The QTildeSudakov class implements the Sudakov form factor for evolution in
* \f$\tilde{q}^2\f$ using the veto algorithm.
*
* @see \ref QTildeSudakovInterfaces "The interfaces"
* defined for QTildeSudakov.
*/
class QTildeSudakov: public SudakovFormFactor {
public:
/**
* The default constructor.
*/
inline QTildeSudakov() {}
/**
* 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
* defined.
* @param enhance The radiation enhancement factor
*/
virtual ShoKinPtr generateNextTimeBranching(const Energy startingScale,
const IdList &ids,const bool cc,
double enhance);
/**
* 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);
/**
* 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 enhance The radiation enhancement factor
* @param beam The beam particle
*/
virtual ShoKinPtr generateNextSpaceBranching(const Energy startingScale,
const IdList &ids,double x,
const bool cc, double enhance,
tcBeamPtr beam);
//@}
/**
* Generate the azimuthal angle of the branching for forward branching
* @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);
/**
* Generate the azimuthal angle of the branching for backward branching
* @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);
+
+ /**
+ * 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);
+
/**
* 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);
/**
* Method to create the ShowerKinematics object for a final-state branching
*/
virtual ShoKinPtr createFinalStateBranching(Energy scale,double z,
double phi, Energy pt);
/**
* Method to create the ShowerKinematics object for an initial-state branching
*/
virtual ShoKinPtr createInitialStateBranching(Energy scale,double z,
double phi, Energy pt);
/**
* Method to create the ShowerKinematics object for a decay branching
*/
virtual ShoKinPtr createDecayBranching(Energy scale,double z,
double phi, Energy pt);
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:
/**
* Methods to provide the next value of the scale before the vetos
* are applied.
*/
//@{
/**
* Value of the energy fraction and scale for time-like branching
* @param t The scale
* @param tmin The minimum scale
* @param enhance The radiation enhancement factor
* @return False if scale less than minimum, true otherwise
*/
bool guessTimeLike(Energy2 &t, Energy2 tmin, double enhance);
/**
* Value of the energy fraction and scale for time-like branching
* @param t The scale
* @param tmax The maximum scale
* @param minmass The minimum mass of the particle after the branching
* @param enhance The radiation enhancement factor
*/
bool guessDecay(Energy2 &t, Energy2 tmax,Energy minmass,
double enhance);
/**
* Value of the energy fraction and scale for space-like branching
* @param t The scale
* @param tmin The minimum scale
* @param x Fraction of the beam momentum.
* @param enhance The radiation enhancement factor
*/
bool guessSpaceLike(Energy2 &t, Energy2 tmin, const double x,
double enhance);
//@}
/**
* Initialize the values of the cut-offs and scales
* @param tmin The minimum scale
* @param ids The ids of the partics in the branching
* @param cc Whether this is the charge conjugate of the branching
*/
void initialize(const IdList & ids,Energy2 &tmin, const bool cc);
/**
* Phase Space veto member to implement the \f$\Theta\f$ function as a veto
* so that the emission is within the allowed phase space.
* @param t The scale
* @return true if vetoed
*/
bool PSVeto(const Energy2 t);
/**
* Compute the limits on \f$z\f$ for time-like branching
* @param scale The scale of the particle
* @return True if lower limit less than upper, otherwise false
*/
bool computeTimeLikeLimits(Energy2 & scale);
/**
* Compute the limits on \f$z\f$ for space-like branching
* @param scale The scale of the particle
* @param x The energy fraction of the parton
* @return True if lower limit less than upper, otherwise false
*/
bool computeSpaceLikeLimits(Energy2 & scale, double x);
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
inline virtual IBPtr clone() const {return new_ptr(*this);}
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
inline virtual IBPtr fullclone() const {return new_ptr(*this);}
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
QTildeSudakov & operator=(const QTildeSudakov &);
private:
/**
* The evolution scale, \f$\tilde{q}\f$.
*/
Energy q_;
/**
* The Ids of the particles in the current branching
*/
IdList ids_;
/**
* The masses of the particles in the current branching
*/
vector<Energy> masses_;
/**
* The mass squared of the particles in the current branching
*/
vector<Energy2> masssquared_;
};
}
#endif /* HERWIG_QTildeSudakov_H */
diff --git a/Shower/SplittingFunctions/SplittingGenerator.cc b/Shower/SplittingFunctions/SplittingGenerator.cc
--- a/Shower/SplittingFunctions/SplittingGenerator.cc
+++ b/Shower/SplittingFunctions/SplittingGenerator.cc
@@ -1,549 +1,550 @@
// -*- C++ -*-
//
// SplittingGenerator.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 SplittingGenerator class.
//
#include "SplittingGenerator.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Command.h"
#include "ThePEG/Utilities/StringUtils.h"
#include "ThePEG/Repository/Repository.h"
#include "Herwig++/Shower/Base/ShowerParticle.h"
#include "ThePEG/Utilities/Rebinder.h"
#include <cassert>
#include "ThePEG/Utilities/DescribeClass.h"
using namespace Herwig;
DescribeClass<SplittingGenerator,Interfaced>
describeSplittingGenerator ("Herwig::SplittingGenerator","");
IBPtr SplittingGenerator::clone() const {
return new_ptr(*this);
}
IBPtr SplittingGenerator::fullclone() const {
return new_ptr(*this);
}
void SplittingGenerator::persistentOutput(PersistentOStream & os) const {
os << _isr_Mode << _fsr_Mode << _bbranchings << _fbranchings;
}
void SplittingGenerator::persistentInput(PersistentIStream & is, int) {
is >> _isr_Mode >> _fsr_Mode >> _bbranchings >> _fbranchings;
}
void SplittingGenerator::Init() {
static ClassDocumentation<SplittingGenerator> documentation
("There class is responsible for initializing the Sudakov form factors ",
"and generating splittings.");
static Switch<SplittingGenerator, bool> interfaceISRMode
("ISR",
"Include initial-state radiation?",
&SplittingGenerator::_isr_Mode, 1, false, false);
static SwitchOption interfaceISRMode0
(interfaceISRMode,"No","ISR (Initial State Radiation) is OFF", 0);
static SwitchOption interfaceISRMode1
(interfaceISRMode,"Yes","ISR (Initial State Radiation) is ON", 1);
static Switch<SplittingGenerator, bool> interfaceFSRMode
("FSR",
"Include final-state radiation?",
&SplittingGenerator::_fsr_Mode, 1, false, false);
static SwitchOption interfaceFSRMode0
(interfaceFSRMode,"No","FSR (Final State Radiation) is OFF", 0);
static SwitchOption interfaceFSRMode1
(interfaceFSRMode,"Yes","FSR (Final State Radiation) is ON", 1);
static Command<SplittingGenerator> interfaceAddSplitting
("AddFinalSplitting",
"Adds another splitting to the list of splittings considered "
"in the shower. Command is a->b,c; Sudakov",
&SplittingGenerator::addFinalSplitting);
static Command<SplittingGenerator> interfaceAddInitialSplitting
("AddInitialSplitting",
"Adds another splitting to the list of initial splittings to consider "
"in the shower. Command is a->b,c; Sudakov. Here the particle a is the "
"particle that is PRODUCED by the splitting. b is the initial state "
"particle that is splitting in the shower.",
&SplittingGenerator::addInitialSplitting);
static Command<SplittingGenerator> interfaceDeleteSplitting
("DeleteFinalSplitting",
"Deletes a splitting from the list of splittings considered "
"in the shower. Command is a->b,c; Sudakov",
&SplittingGenerator::deleteFinalSplitting);
static Command<SplittingGenerator> interfaceDeleteInitialSplitting
("DeleteInitialSplitting",
"Deletes a splitting from the list of initial splittings to consider "
"in the shower. Command is a->b,c; Sudakov. Here the particle a is the "
"particle that is PRODUCED by the splitting. b is the initial state "
"particle that is splitting in the shower.",
&SplittingGenerator::deleteInitialSplitting);
}
string SplittingGenerator::addSplitting(string arg, bool final) {
string partons = StringUtils::car(arg);
string sudakov = StringUtils::cdr(arg);
vector<tPDPtr> products;
string::size_type next = partons.find("->");
if(next == string::npos)
return "Error: Invalid string for splitting " + arg;
if(partons.find(';') == string::npos)
return "Error: Invalid string for splitting " + arg;
tPDPtr parent = Repository::findParticle(partons.substr(0,next));
partons = partons.substr(next+2);
do {
next = min(partons.find(','), partons.find(';'));
tPDPtr pdp = Repository::findParticle(partons.substr(0,next));
partons = partons.substr(next+1);
if(pdp) products.push_back(pdp);
else return "Error: Could not create splitting from " + arg;
} while(partons[0] != ';' && partons.size());
SudakovPtr s;
s = dynamic_ptr_cast<SudakovPtr>(Repository::TraceObject(sudakov));
if(!s) return "Error: Could not load Sudakov " + sudakov + '\n';
IdList ids;
ids.push_back(parent->id());
for(vector<tPDPtr>::iterator it = products.begin(); it!=products.end(); ++it)
ids.push_back((*it)->id());
// check splitting can handle this
if(!s->splittingFn()->accept(ids))
return "Error: Sudakov " + sudakov + "can't handle particles\n";
// add to map
addToMap(ids,s,final);
return "";
}
string SplittingGenerator::deleteSplitting(string arg, bool final) {
string partons = StringUtils::car(arg);
string sudakov = StringUtils::cdr(arg);
vector<tPDPtr> products;
string::size_type next = partons.find("->");
if(next == string::npos)
return "Error: Invalid string for splitting " + arg;
if(partons.find(';') == string::npos)
return "Error: Invalid string for splitting " + arg;
tPDPtr parent = Repository::findParticle(partons.substr(0,next));
partons = partons.substr(next+2);
do {
next = min(partons.find(','), partons.find(';'));
tPDPtr pdp = Repository::findParticle(partons.substr(0,next));
partons = partons.substr(next+1);
if(pdp) products.push_back(pdp);
else return "Error: Could not create splitting from " + arg;
} while(partons[0] != ';' && partons.size());
SudakovPtr s;
s = dynamic_ptr_cast<SudakovPtr>(Repository::TraceObject(sudakov));
if(!s) return "Error: Could not load Sudakov " + sudakov + '\n';
IdList ids;
ids.push_back(parent->id());
for(vector<tPDPtr>::iterator it = products.begin(); it!=products.end(); ++it)
ids.push_back((*it)->id());
// check splitting can handle this
if(!s->splittingFn()->accept(ids))
return "Error: Sudakov " + sudakov + "can't handle particles\n";
// delete from map
deleteFromMap(ids,s,final);
return "";
}
void SplittingGenerator::addToMap(const IdList &ids, const SudakovPtr &s, bool final) {
if(isISRadiationON() && !final) {
_bbranchings.insert(BranchingInsert(ids[1],BranchingElement(s,ids)));
s->addSplitting(ids);
}
if(isFSRadiationON() && final) {
_fbranchings.insert(BranchingInsert(ids[0],BranchingElement(s,ids)));
s->addSplitting(ids);
}
}
void SplittingGenerator::deleteFromMap(const IdList &ids,
const SudakovPtr &s, bool final) {
if(isISRadiationON() && !final) {
pair<BranchingList::iterator,BranchingList::iterator>
range = _bbranchings.equal_range(ids[1]);
for(BranchingList::iterator it=range.first;it!=range.second&&it->first==ids[1];++it) {
if(it->second.first==s&&it->second.second==ids)
_bbranchings.erase(it);
}
s->removeSplitting(ids);
}
if(isFSRadiationON() && final) {
pair<BranchingList::iterator,BranchingList::iterator>
range = _fbranchings.equal_range(ids[0]);
for(BranchingList::iterator it=range.first;it!=range.second&&it->first==ids[0];++it) {
if(it->second.first==s&&it->second.second==ids)
_fbranchings.erase(it);
}
s->removeSplitting(ids);
}
}
Branching SplittingGenerator::chooseForwardBranching(ShowerParticle &particle,
double enhance,
ShowerInteraction::Type type) const {
Energy newQ = ZERO;
ShoKinPtr kinematics = ShoKinPtr();
ShowerPartnerType::Type partnerType(ShowerPartnerType::Undefined);
SudakovPtr sudakov = SudakovPtr();
IdList ids;
// First, find the eventual branching, corresponding to the highest scale.
long index = abs(particle.data().id());
// if no branchings return empty branching struct
if( _fbranchings.find(index) == _fbranchings.end() )
return Branching(ShoKinPtr(), IdList(),SudakovPtr(),ShowerPartnerType::Undefined);
// otherwise select branching
for(BranchingList::const_iterator cit = _fbranchings.lower_bound(index);
cit != _fbranchings.upper_bound(index); ++cit) {
// check either right interaction or doing both
if(type != cit->second.first->interactionType() &&
type != ShowerInteraction::Both ) continue;
// whether or not this interaction should be angular ordered
bool angularOrdered = cit->second.first->splittingFn()->angularOrdered();
ShoKinPtr newKin;
ShowerPartnerType::Type type;
// work out which starting scale we need
if(cit->second.first->interactionType()==ShowerInteraction::QED) {
type = ShowerPartnerType::QED;
Energy startingScale = angularOrdered ? particle.scales().QED : particle.scales().QED_noAO;
newKin = cit->second.first->
generateNextTimeBranching(startingScale,cit->second.second,
particle.id()!=cit->first,enhance);
}
else if(cit->second.first->interactionType()==ShowerInteraction::QCD) {
// special for octets
if(particle.dataPtr()->iColour()==PDT::Colour8) {
// octet -> octet octet
if(cit->second.first->splittingFn()->colourStructure()==OctetOctetOctet) {
type = ShowerPartnerType::QCDColourLine;
Energy startingScale = angularOrdered ? particle.scales().QCD_c : particle.scales().QCD_c_noAO;
newKin= cit->second.first->
generateNextTimeBranching(startingScale,cit->second.second,
particle.id()!=cit->first,0.5*enhance);
startingScale = angularOrdered ? particle.scales().QCD_ac : particle.scales().QCD_ac_noAO;
ShoKinPtr newKin2 = cit->second.first->
generateNextTimeBranching(startingScale,cit->second.second,
particle.id()!=cit->first,0.5*enhance);
// pick the one with the highest scale
if( ( newKin && newKin2 && newKin2->scale() > newKin->scale()) ||
(!newKin && newKin2) ) {
newKin = newKin2;
type = ShowerPartnerType::QCDAntiColourLine;
}
}
// other g -> q qbar
else {
Energy startingScale = angularOrdered ?
max(particle.scales().QCD_c , particle.scales().QCD_ac ) :
max(particle.scales().QCD_c_noAO, particle.scales().QCD_ac_noAO);
newKin= cit->second.first->
generateNextTimeBranching(startingScale, cit->second.second,
particle.id()!=cit->first,enhance);
type = UseRandom::rndbool() ?
ShowerPartnerType::QCDColourLine : ShowerPartnerType::QCDAntiColourLine;
}
}
// everything else q-> qg etc
else {
Energy startingScale;
if(particle.hasColour()) {
type = ShowerPartnerType::QCDColourLine;
startingScale = angularOrdered ? particle.scales().QCD_c : particle.scales().QCD_c_noAO;
}
else {
type = ShowerPartnerType::QCDAntiColourLine;
startingScale = angularOrdered ? particle.scales().QCD_ac : particle.scales().QCD_ac_noAO;
}
newKin= cit->second.first->
generateNextTimeBranching(startingScale,cit->second.second,
particle.id()!=cit->first,enhance);
}
}
// shouldn't be anything else
else
assert(false);
// if no kinematics contine
if(!newKin) continue;
// select highest scale
if( newKin->scale() > newQ ) {
kinematics = newKin;
newQ = newKin->scale();
ids = cit->second.second;
sudakov = cit->second.first;
partnerType = type;
}
}
// return empty branching if nothing happened
if(!kinematics) return Branching(ShoKinPtr(), IdList(),SudakovPtr(),
ShowerPartnerType::Undefined);
// If a branching has been selected initialize it
kinematics->initialize(particle,PPtr());
// and generate phi
kinematics->phi(sudakov->generatePhiForward(particle,ids,kinematics));
// and return it
return Branching(kinematics, ids,sudakov,partnerType);
}
Branching SplittingGenerator::
chooseDecayBranching(ShowerParticle &particle,
const ShowerParticle::EvolutionScales & stoppingScales,
Energy minmass, double enhance,
ShowerInteraction::Type interaction) const {
Energy newQ = Constants::MaxEnergy;
ShoKinPtr kinematics;
SudakovPtr sudakov;
ShowerPartnerType::Type partnerType(ShowerPartnerType::Undefined);
IdList ids;
// First, find the eventual branching, corresponding to the lowest scale.
long index = abs(particle.data().id());
// if no branchings return empty branching struct
if(_fbranchings.find(index) == _fbranchings.end())
return Branching(ShoKinPtr(), IdList(),SudakovPtr(),ShowerPartnerType::Undefined);
// otherwise select branching
for(BranchingList::const_iterator cit = _fbranchings.lower_bound(index);
cit != _fbranchings.upper_bound(index); ++cit) {
// check interaction doesn't change flavour
if(cit->second.second[1]!=index&&cit->second.second[2]!=index) continue;
// check either right interaction or doing both
if(interaction != cit->second.first->interactionType() &&
interaction != ShowerInteraction::Both ) continue;
// whether or not this interaction should be angular ordered
bool angularOrdered = cit->second.first->splittingFn()->angularOrdered();
ShoKinPtr newKin;
ShowerPartnerType::Type type;
// work out which starting scale we need
if(cit->second.first->interactionType()==ShowerInteraction::QED) {
type = ShowerPartnerType::QED;
Energy stoppingScale = angularOrdered ? stoppingScales.QED : stoppingScales.QED_noAO;
Energy startingScale = angularOrdered ? particle.scales().QED : particle.scales().QED_noAO;
if(startingScale < stoppingScale ) {
newKin = cit->second.first->
generateNextDecayBranching(startingScale,stoppingScale,minmass,cit->second.second,
particle.id()!=cit->first,enhance);
}
}
else if(cit->second.first->interactionType()==ShowerInteraction::QCD) {
// special for octets
if(particle.dataPtr()->iColour()==PDT::Colour8) {
// octet -> octet octet
if(cit->second.first->splittingFn()->colourStructure()==OctetOctetOctet) {
Energy stoppingColour = angularOrdered ? stoppingScales.QCD_c : stoppingScales.QCD_c_noAO;
Energy stoppingAnti = angularOrdered ? stoppingScales.QCD_ac : stoppingScales.QCD_ac_noAO;
Energy startingColour = angularOrdered ? particle.scales().QCD_c : particle.scales().QCD_c_noAO;
Energy startingAnti = angularOrdered ? particle.scales().QCD_ac : particle.scales().QCD_ac_noAO;
type = ShowerPartnerType::QCDColourLine;
if(startingColour<stoppingColour) {
newKin= cit->second.first->
generateNextDecayBranching(startingColour,stoppingColour,minmass,
cit->second.second,
particle.id()!=cit->first,0.5*enhance);
}
ShoKinPtr newKin2;
if(startingAnti<stoppingAnti) {
newKin2 = cit->second.first->
generateNextDecayBranching(startingAnti,stoppingAnti,minmass,
cit->second.second,
particle.id()!=cit->first,0.5*enhance);
}
// pick the one with the lowest scale
if( (newKin&&newKin2&&newKin2->scale()<newKin->scale()) ||
(!newKin&&newKin2) ) {
newKin = newKin2;
type = ShowerPartnerType::QCDAntiColourLine;
}
}
// other
else {
assert(false);
}
}
// everything else
else {
Energy startingScale,stoppingScale;
if(particle.hasColour()) {
type = ShowerPartnerType::QCDColourLine;
stoppingScale = angularOrdered ? stoppingScales.QCD_c : stoppingScales.QCD_c_noAO;
startingScale = angularOrdered ? particle.scales().QCD_c : particle.scales().QCD_c_noAO;
}
else {
type = ShowerPartnerType::QCDAntiColourLine;
stoppingScale = angularOrdered ? stoppingScales.QCD_ac : stoppingScales.QCD_ac_noAO;
startingScale = angularOrdered ? particle.scales().QCD_ac : particle.scales().QCD_ac_noAO;
}
if(startingScale < stoppingScale ) {
newKin = cit->second.first->
generateNextDecayBranching(startingScale,stoppingScale,minmass,cit->second.second,
particle.id()!=cit->first,enhance);
}
}
}
// shouldn't be anything else
else
assert(false);
if(!newKin) continue;
// select highest scale
if(newKin->scale() < newQ ) {
newQ = newKin->scale();
ids = cit->second.second;
kinematics=newKin;
sudakov=cit->second.first;
partnerType = type;
}
}
// return empty branching if nothing happened
if(!kinematics) return Branching(ShoKinPtr(), IdList(),SudakovPtr(),
ShowerPartnerType::Undefined);
// initialize the branching
kinematics->initialize(particle,PPtr());
+ // and generate phi
+ kinematics->phi(sudakov->generatePhiDecay(particle,ids,kinematics));
// and return it
return Branching(kinematics, ids,sudakov,partnerType);
}
Branching SplittingGenerator::
chooseBackwardBranching(ShowerParticle &particle,PPtr beamparticle,
double enhance,
Ptr<BeamParticleData>::transient_const_pointer beam,
ShowerInteraction::Type type,
tcPDFPtr pdf, Energy freeze) const {
Energy newQ=ZERO;
ShoKinPtr kinematics=ShoKinPtr();
ShowerPartnerType::Type partnerType(ShowerPartnerType::Undefined);
SudakovPtr sudakov;
IdList ids;
// First, find the eventual branching, corresponding to the highest scale.
long index = abs(particle.id());
// if no possible branching return
if(_bbranchings.find(index) == _bbranchings.end())
return Branching(ShoKinPtr(), IdList(),SudakovPtr(),ShowerPartnerType::Undefined);
// otherwise select branching
for(BranchingList::const_iterator cit = _bbranchings.lower_bound(index);
cit != _bbranchings.upper_bound(index); ++cit ) {
// check either right interaction or doing both
if(type != cit->second.first->interactionType() &&
type != ShowerInteraction::Both ) continue;
// setup the PDF
cit->second.first->setPDF(pdf,freeze);
// whether or not this interaction should be angular ordered
bool angularOrdered = cit->second.first->splittingFn()->angularOrdered();
ShoKinPtr newKin;
ShowerPartnerType::Type type;
if(cit->second.first->interactionType()==ShowerInteraction::QED) {
type = ShowerPartnerType::QED;
Energy startingScale = angularOrdered ? particle.scales().QED : particle.scales().QED_noAO;
newKin=cit->second.first->
generateNextSpaceBranching(startingScale,cit->second.second, particle.x(),
particle.id()!=cit->first,enhance,beam);
}
else if(cit->second.first->interactionType()==ShowerInteraction::QCD) {
// special for octets
if(particle.dataPtr()->iColour()==PDT::Colour8) {
// octet -> octet octet
if(cit->second.first->splittingFn()->colourStructure()==OctetOctetOctet) {
type = ShowerPartnerType::QCDColourLine;
Energy startingScale = angularOrdered ? particle.scales().QCD_c : particle.scales().QCD_c_noAO;
newKin = cit->second.first->
generateNextSpaceBranching(startingScale,cit->second.second, particle.x(),
particle.id()!=cit->first,0.5*enhance,beam);
startingScale = angularOrdered ? particle.scales().QCD_ac : particle.scales().QCD_ac_noAO;
ShoKinPtr newKin2 = cit->second.first->
generateNextSpaceBranching(startingScale,cit->second.second, particle.x(),
particle.id()!=cit->first,0.5*enhance,beam);
// pick the one with the highest scale
if( (newKin&&newKin2&&newKin2->scale()>newKin->scale()) ||
(!newKin&&newKin2) ) {
newKin = newKin2;
type = ShowerPartnerType::QCDAntiColourLine;
}
}
else {
Energy startingScale = angularOrdered ?
max(particle.scales().QCD_c , particle.scales().QCD_ac ) :
max(particle.scales().QCD_c_noAO, particle.scales().QCD_ac_noAO);
type = UseRandom::rndbool() ?
ShowerPartnerType::QCDColourLine : ShowerPartnerType::QCDAntiColourLine;
newKin=cit->second.first->
generateNextSpaceBranching(startingScale,cit->second.second, particle.x(),
particle.id()!=cit->first,enhance,beam);
}
}
// everything else
else {
Energy startingScale;
if(particle.hasColour()) {
type = ShowerPartnerType::QCDColourLine;
startingScale = angularOrdered ? particle.scales().QCD_c : particle.scales().QCD_c_noAO;
}
else {
type = ShowerPartnerType::QCDAntiColourLine;
startingScale = angularOrdered ? particle.scales().QCD_ac : particle.scales().QCD_ac_noAO;
}
newKin=cit->second.first->
generateNextSpaceBranching(startingScale,cit->second.second, particle.x(),
particle.id()!=cit->first,enhance,beam);
}
}
// shouldn't be anything else
else
assert(false);
// if no kinematics contine
if(!newKin) continue;
// select highest scale
if(newKin->scale() > newQ) {
newQ = newKin->scale();
kinematics=newKin;
ids = cit->second.second;
sudakov=cit->second.first;
partnerType = type;
}
}
// return empty branching if nothing happened
if(!kinematics) return Branching(ShoKinPtr(), IdList(),SudakovPtr(),
ShowerPartnerType::Undefined);
// initialize the ShowerKinematics
// and return it
kinematics->initialize(particle,beamparticle);
// and generate phi
kinematics->phi(sudakov->generatePhiBackward(particle,ids,kinematics));
// return the answer
return Branching(kinematics, ids,sudakov,partnerType);
}
-void SplittingGenerator::rebind(const TranslationMap & trans)
- {
+void SplittingGenerator::rebind(const TranslationMap & trans) {
BranchingList::iterator cit;
for(cit=_fbranchings.begin();cit!=_fbranchings.end();++cit)
{(cit->second).first=trans.translate((cit->second).first);}
for(cit=_bbranchings.begin();cit!=_bbranchings.end();++cit)
{(cit->second).first=trans.translate((cit->second).first);}
Interfaced::rebind(trans);
}
IVector SplittingGenerator::getReferences() {
IVector ret = Interfaced::getReferences();
BranchingList::iterator cit;
for(cit=_fbranchings.begin();cit!=_fbranchings.end();++cit)
{ret.push_back((cit->second).first);}
for(cit=_bbranchings.begin();cit!=_bbranchings.end();++cit)
{ret.push_back((cit->second).first);}
return ret;
}

File Metadata

Mime Type
text/x-diff
Expires
Mon, Jan 20, 11:28 PM (1 d, 9 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4242849
Default Alt Text
(110 KB)

Event Timeline