Page MenuHomeHEPForge

No OneTemporary

This file is larger than 256 KB, so syntax highlighting was skipped.
diff --git a/Decay/Perturbative/SMTopDecayer.cc b/Decay/Perturbative/SMTopDecayer.cc
--- a/Decay/Perturbative/SMTopDecayer.cc
+++ b/Decay/Perturbative/SMTopDecayer.cc
@@ -1,1201 +1,1202 @@
// -*- C++ -*-
//
// SMTopDecayer.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 SMTopDecayer class.
//
#include "SMTopDecayer.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/DecayMode.h"
#include "Herwig++/Decay/DecayVertex.h"
#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
#include "Herwig++/PDT/ThreeBodyAllOn1IntegralCalculator.h"
using namespace Herwig;
using namespace ThePEG::Helicity;
SMTopDecayer::SMTopDecayer()
: _wquarkwgt(6,0.),_wleptonwgt(3,0.), _xg_sampling(1.5),
_initialenhance(1.), _finalenhance(2.3), _useMEforT2(true) {
_wleptonwgt[0] = 0.302583;
_wleptonwgt[1] = 0.301024;
_wleptonwgt[2] = 0.299548;
_wquarkwgt[0] = 0.851719;
_wquarkwgt[1] = 0.0450162;
_wquarkwgt[2] = 0.0456962;
_wquarkwgt[3] = 0.859839;
_wquarkwgt[4] = 3.9704e-06;
_wquarkwgt[5] = 0.000489657;
generateIntermediates(true);
}
bool SMTopDecayer::accept(tcPDPtr parent, const tPDVector & children) const {
if(abs(parent->id()) != ParticleID::t) return false;
int id0(0),id1(0),id2(0);
for(tPDVector::const_iterator it = children.begin();
it != children.end();++it) {
int id=(**it).id(),absid(abs(id));
if(absid==ParticleID::b&&double(id)/double(parent->id())>0) {
id0=id;
}
else {
switch (absid) {
case ParticleID::nu_e:
case ParticleID::nu_mu:
case ParticleID::nu_tau:
id1 = id;
break;
case ParticleID::eminus:
case ParticleID::muminus:
case ParticleID::tauminus:
id2 = id;
break;
case ParticleID::b:
case ParticleID::d:
case ParticleID::s:
id1 = id;
break;
case ParticleID::u:
case ParticleID::c:
id2=id;
break;
default :
break;
}
}
}
if(id0==0||id1==0||id2==0) return false;
if(double(id1)/double(id2)>0) return false;
return true;
}
ParticleVector SMTopDecayer::decay(const Particle & parent,
const tPDVector & children) const {
int id1(0),id2(0);
for(tPDVector::const_iterator it = children.begin();
it != children.end();++it) {
int id=(**it).id(),absid=abs(id);
if(absid == ParticleID::b && double(id)/double(parent.id())>0) continue;
//leptons
if(absid > 10 && absid%2==0) id1=absid;
if(absid > 10 && absid%2==1) id2=absid;
//quarks
if(absid < 10 && absid%2==0) id2=absid;
if(absid < 10 && absid%2==1) id1=absid;
}
unsigned int imode(0);
if(id2 >=11 && id2<=16) imode = (id1-12)/2;
else imode = id1+1+id2/2;
bool cc = parent.id() == ParticleID::tbar;
ParticleVector out(generate(true,cc,imode,parent));
//arrange colour flow
PPtr pparent=const_ptr_cast<PPtr>(&parent);
out[1]->incomingColour(pparent,out[1]->id()<0);
ParticleVector products = out[0]->children();
if(products[0]->hasColour())
products[0]->colourNeighbour(products[1],true);
else if(products[0]->hasAntiColour())
products[0]->colourNeighbour(products[1],false);
return out;
}
void SMTopDecayer::persistentOutput(PersistentOStream & os) const {
os << _wvertex << _wquarkwgt << _wleptonwgt << _wplus << _alpha
<< _initialenhance << _finalenhance << _xg_sampling << _useMEforT2;
}
void SMTopDecayer::persistentInput(PersistentIStream & is, int) {
is >> _wvertex >> _wquarkwgt >> _wleptonwgt >> _wplus >> _alpha
>> _initialenhance >> _finalenhance >> _xg_sampling >> _useMEforT2;
}
ClassDescription<SMTopDecayer> SMTopDecayer::initSMTopDecayer;
// Definition of the static class description member.
void SMTopDecayer::Init() {
static ClassDocumentation<SMTopDecayer> documentation
("This is the implementation of the SMTopDecayer which "
"decays top quarks into bottom quarks and either leptons "
"or quark-antiquark pairs including the matrix element for top decay",
"The matrix element correction for top decay \\cite{Hamilton:2006ms}.",
"%\\cite{Hamilton:2006ms}\n"
"\\bibitem{Hamilton:2006ms}\n"
" K.~Hamilton and P.~Richardson,\n"
" ``A simulation of QCD radiation in top quark decays,''\n"
" JHEP {\\bf 0702}, 069 (2007)\n"
" [arXiv:hep-ph/0612236].\n"
" %%CITATION = JHEPA,0702,069;%%\n");
static ParVector<SMTopDecayer,double> interfaceQuarkWeights
("QuarkWeights",
"Maximum weights for the hadronic decays",
&SMTopDecayer::_wquarkwgt, 6, 1.0, 0.0, 10.0,
false, false, Interface::limited);
static ParVector<SMTopDecayer,double> interfaceLeptonWeights
("LeptonWeights",
"Maximum weights for the semi-leptonic decays",
&SMTopDecayer::_wleptonwgt, 3, 1.0, 0.0, 10.0,
false, false, Interface::limited);
static Parameter<SMTopDecayer,double> interfaceEnhancementFactor
("InitialEnhancementFactor",
"The enhancement factor for initial-state radiation in the shower to ensure"
" the weight for the matrix element correction is less than one.",
&SMTopDecayer::_initialenhance, 1.0, 1.0, 10000.0,
false, false, Interface::limited);
static Parameter<SMTopDecayer,double> interfaceFinalEnhancementFactor
("FinalEnhancementFactor",
"The enhancement factor for final-state radiation in the shower to ensure"
" the weight for the matrix element correction is less than one",
&SMTopDecayer::_finalenhance, 1.6, 1.0, 1000.0,
false, false, Interface::limited);
static Parameter<SMTopDecayer,double> interfaceSamplingTopHardMEC
("SamplingTopHardMEC",
"The importance sampling power for choosing an initial xg, "
"to sample xg according to xg^-_xg_sampling",
&SMTopDecayer::_xg_sampling, 1.5, 1.2, 2.0,
false, false, Interface::limited);
static Switch<SMTopDecayer,bool> interfaceUseMEForT2
("UseMEForT2",
"Use the matrix element correction, if available to fill the T2"
" region for the decay shower and don't fill using the shower",
&SMTopDecayer::_useMEforT2, true, false, false);
static SwitchOption interfaceUseMEForT2Shower
(interfaceUseMEForT2,
"Shower",
"Use the shower to fill the T2 region",
false);
static SwitchOption interfaceUseMEForT2ME
(interfaceUseMEForT2,
"ME",
"Use the Matrix element to fill the T2 region",
true);
static Reference<SMTopDecayer,ShowerAlpha> interfaceCoupling
("Coupling",
"Pointer to the object to calculate the coupling for the correction",
&SMTopDecayer::_alpha, false, false, true, false, false);
}
double SMTopDecayer::me2(const int, const Particle & inpart,
const ParticleVector & decay,
MEOption meopt) const {
// spinors etc for the decaying particle
if(meopt==Initialize) {
// spinors and rho
if(inpart.id()>0)
SpinorWaveFunction ::calculateWaveFunctions(_inHalf,_rho,
const_ptr_cast<tPPtr>(&inpart),
incoming);
else
SpinorBarWaveFunction::calculateWaveFunctions(_inHalfBar,_rho,
const_ptr_cast<tPPtr>(&inpart),
incoming);
ME(DecayMatrixElement(PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1Half,PDT::Spin1Half));
}
// setup spin info when needed
if(meopt==Terminate) {
// for the decaying particle
if(inpart.id()>0) {
SpinorWaveFunction::
constructSpinInfo(_inHalf,const_ptr_cast<tPPtr>(&inpart),incoming,true);
SpinorBarWaveFunction::constructSpinInfo(_inHalfBar,decay[0],outgoing,true);
SpinorWaveFunction ::constructSpinInfo(_outHalf ,decay[1],outgoing,true);
SpinorBarWaveFunction::constructSpinInfo(_outHalfBar,decay[2],outgoing,true);
}
else {
SpinorBarWaveFunction::
constructSpinInfo(_inHalfBar,const_ptr_cast<tPPtr>(&inpart),incoming,true);
SpinorWaveFunction::constructSpinInfo(_inHalf,decay[0],outgoing,true);
SpinorBarWaveFunction::constructSpinInfo(_outHalfBar,decay[1],outgoing,true);
SpinorWaveFunction ::constructSpinInfo(_outHalf ,decay[2],outgoing,true);
}
}
if ( ( decay[1]->momentum() + decay[2]->momentum() ).m()
< decay[1]->data().constituentMass() + decay[2]->data().constituentMass() )
return 0.0;
// spinors for the decay product
if(inpart.id()>0) {
SpinorBarWaveFunction::calculateWaveFunctions(_inHalfBar ,decay[0],outgoing);
SpinorWaveFunction ::calculateWaveFunctions(_outHalf ,decay[1],outgoing);
SpinorBarWaveFunction::calculateWaveFunctions(_outHalfBar,decay[2],outgoing);
}
else {
SpinorWaveFunction ::calculateWaveFunctions(_inHalf ,decay[0],outgoing);
SpinorBarWaveFunction::calculateWaveFunctions(_outHalfBar,decay[1],outgoing);
SpinorWaveFunction ::calculateWaveFunctions(_outHalf ,decay[2],outgoing);
}
Energy2 scale(sqr(inpart.mass()));
if(inpart.id() == ParticleID::t) {
//Define intermediate vector wave-function for Wplus
tcPDPtr Wplus(getParticleData(ParticleID::Wplus));
VectorWaveFunction inter;
unsigned int thel,bhel,fhel,afhel;
for(thel = 0;thel<2;++thel){
for(bhel = 0;bhel<2;++bhel){
inter = _wvertex->evaluate(scale,1,Wplus,_inHalf[thel],
_inHalfBar[bhel]);
for(afhel=0;afhel<2;++afhel){
for(fhel=0;fhel<2;++fhel){
ME()(thel,bhel,afhel,fhel) =
_wvertex->evaluate(scale,_outHalf[afhel],
_outHalfBar[fhel],inter);
}
}
}
}
}
else if(inpart.id() == ParticleID::tbar) {
VectorWaveFunction inter;
tcPDPtr Wminus(getParticleData(ParticleID::Wminus));
unsigned int tbhel,bbhel,afhel,fhel;
for(tbhel = 0;tbhel<2;++tbhel){
for(bbhel = 0;bbhel<2;++bbhel){
inter = _wvertex->
evaluate(scale,1,Wminus,_inHalf[bbhel],_inHalfBar[tbhel]);
for(afhel=0;afhel<2;++afhel){
for(fhel=0;fhel<2;++fhel){
ME()(tbhel,bbhel,fhel,afhel) =
_wvertex->evaluate(scale,_outHalf[afhel],
_outHalfBar[fhel],inter);
}
}
}
}
}
double output = (ME().contract(_rho)).real();
if(abs(decay[1]->id())<=6) output *=3.;
return output;
}
void SMTopDecayer::doinit() {
DecayIntegrator::doinit();
//get vertices from SM object
tcHwSMPtr hwsm = dynamic_ptr_cast<tcHwSMPtr>(standardModel());
if(!hwsm) throw InitException() << "Must have Herwig::StandardModel in "
<< "SMTopDecayer::doinit()";
_wvertex = hwsm->vertexFFW();
//initialise
_wvertex->init();
//set up decay modes
_wplus = getParticleData(ParticleID::Wplus);
DecayPhaseSpaceModePtr mode;
DecayPhaseSpaceChannelPtr Wchannel;
tPDVector extpart(4);
vector<double> wgt(1,1.0);
extpart[0] = getParticleData(ParticleID::t);
extpart[1] = getParticleData(ParticleID::b);
//lepton modes
for(int i=11; i<17;i+=2) {
extpart[2] = getParticleData(-i);
extpart[3] = getParticleData(i+1);
mode = new_ptr(DecayPhaseSpaceMode(extpart,this));
Wchannel = new_ptr(DecayPhaseSpaceChannel(mode));
Wchannel->addIntermediate(extpart[0],0,0.0,-1,1);
Wchannel->addIntermediate(_wplus,0,0.0,2,3);
Wchannel->init();
mode->addChannel(Wchannel);
addMode(mode,_wleptonwgt[(i-11)/2],wgt);
}
//quark modes
unsigned int iz=0;
for(int ix=1;ix<6;ix+=2) {
for(int iy=2;iy<6;iy+=2) {
// check that the combination of particles is allowed
if(_wvertex->allowed(-ix,iy,ParticleID::Wminus)) {
extpart[2] = getParticleData(-ix);
extpart[3] = getParticleData( iy);
mode = new_ptr(DecayPhaseSpaceMode(extpart,this));
Wchannel = new_ptr(DecayPhaseSpaceChannel(mode));
Wchannel->addIntermediate(extpart[0],0,0.0,-1,1);
Wchannel->addIntermediate(_wplus,0,0.0,2,3);
Wchannel->init();
mode->addChannel(Wchannel);
addMode(mode,_wquarkwgt[iz],wgt);
++iz;
}
else {
throw InitException() << "SMTopDecayer::doinit() the W vertex"
<< "cannot handle all the quark modes"
<< Exception::abortnow;
}
}
}
}
void SMTopDecayer::dataBaseOutput(ofstream & os,bool header) const {
if(header) os << "update decayers set parameters=\"";
// parameters for the DecayIntegrator base class
for(unsigned int ix=0;ix<_wquarkwgt.size();++ix) {
os << "newdef " << name() << ":QuarkWeights " << ix << " "
<< _wquarkwgt[ix] << "\n";
}
for(unsigned int ix=0;ix<_wleptonwgt.size();++ix) {
os << "newdef " << name() << ":LeptonWeights " << ix << " "
<< _wleptonwgt[ix] << "\n";
}
DecayIntegrator::dataBaseOutput(os,false);
if(header) os << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl;
}
void SMTopDecayer::doinitrun() {
DecayIntegrator::doinitrun();
if(initialize()) {
for(unsigned int ix=0;ix<numberModes();++ix) {
if(ix<3) _wleptonwgt[ix ] = mode(ix)->maxWeight();
else _wquarkwgt [ix-3] = mode(ix)->maxWeight();
}
}
}
WidthCalculatorBasePtr SMTopDecayer::threeBodyMEIntegrator(const DecayMode & dm) const {
// identify W decay products
int sign = dm.parent()->id() > 0 ? 1 : -1;
int iferm(0),ianti(0);
for(ParticleMSet::const_iterator pit=dm.products().begin();
pit!=dm.products().end();++pit) {
int id = (**pit).id();
if(id*sign != ParticleID::b) {
if (id*sign > 0 ) iferm = id*sign;
else ianti = id*sign;
}
}
assert(iferm!=0&&ianti!=0);
// work out which mode we are doing
int imode(-1);
for(unsigned int ix=0;ix<numberModes();++ix) {
if(mode(ix)->externalParticles(2)->id() == ianti &&
mode(ix)->externalParticles(3)->id() == iferm ) {
imode = ix;
break;
}
}
assert(imode>=0);
// get the masses we need
Energy m[3] = {mode(imode)->externalParticles(1)->mass(),
mode(imode)->externalParticles(3)->mass(),
mode(imode)->externalParticles(2)->mass()};
return
new_ptr(ThreeBodyAllOn1IntegralCalculator<SMTopDecayer>
(3,_wplus->mass(),_wplus->width(),0.0,*this,imode,m[0],m[1],m[2]));
}
InvEnergy SMTopDecayer::threeBodydGammads(const int imode, const Energy2 mt2,
const Energy2 mffb2, const Energy mb,
const Energy mf, const Energy mfb) const {
Energy mffb(sqrt(mffb2));
Energy mw(_wplus->mass());
Energy2 mw2(sqr(mw)),gw2(sqr(_wplus->width()));
Energy mt(sqrt(mt2));
Energy Eb = 0.5*(mt2-mffb2-sqr(mb))/mffb;
Energy Ef = 0.5*(mffb2-sqr(mfb)+sqr(mf))/mffb;
Energy Ebm = sqrt(sqr(Eb)-sqr(mb));
Energy Efm = sqrt(sqr(Ef)-sqr(mf));
Energy2 upp = sqr(Eb+Ef)-sqr(Ebm-Efm);
Energy2 low = sqr(Eb+Ef)-sqr(Ebm+Efm);
InvEnergy width=(dGammaIntegrand(mffb2,upp,mt,mb,mf,mfb,mw)-
dGammaIntegrand(mffb2,low,mt,mb,mf,mfb,mw))
/32./mt2/mt/8/pow(Constants::pi,3)/(sqr(mffb2-mw2)+mw2*gw2);
// couplings
width *= 0.25*sqr(4.*Constants::pi*generator()->standardModel()->alphaEM(mt2)/
generator()->standardModel()->sin2ThetaW());
width *= generator()->standardModel()->CKM(*mode(imode)->externalParticles(0),
*mode(imode)->externalParticles(1));
if(abs(mode(imode)->externalParticles(2)->id())<=6) {
width *=3.;
if(abs(mode(imode)->externalParticles(2)->id())%2==0)
width *=generator()->standardModel()->CKM(*mode(imode)->externalParticles(2),
*mode(imode)->externalParticles(3));
else
width *=generator()->standardModel()->CKM(*mode(imode)->externalParticles(3),
*mode(imode)->externalParticles(2));
}
// final spin average
assert(!isnan(width*GeV));
return 0.5*width;
}
Energy6 SMTopDecayer::dGammaIntegrand(Energy2 mffb2, Energy2 mbf2, Energy mt,
Energy mb, Energy mf, Energy mfb, Energy mw) const {
Energy2 mt2(sqr(mt)) ,mb2(sqr(mb)) ,mf2(sqr(mf )),mfb2(sqr(mfb )),mw2(sqr(mw ));
Energy4 mt4(sqr(mt2)),mb4(sqr(mb2)),mf4(sqr(mf2)),mfb4(sqr(mfb2)),mw4(sqr(mw2));
return -mbf2 * ( + 6 * mb2 * mf2 * mfb2 * mffb2 + 6 * mb2 * mt2 * mfb2 * mffb2
+ 6 * mb2 * mt2 * mf2 * mffb2 + 12 * mb2 * mt2 * mf2 * mfb2
- 3 * mb2 * mfb4 * mffb2 + 3 * mb2 * mf2 * mffb2 * mffb2
- 3 * mb2 * mf4 * mffb2 - 6 * mb2 * mt2 * mfb4
- 6 * mb2 * mt2 * mf4 - 3 * mb4 * mfb2 * mffb2
- 3 * mb4 * mf2 * mffb2 - 6 * mb4 * mf2 * mfb2
+ 3 * mt4 * mf4 + 3 * mb4 * mfb4
+ 3 * mb4 * mf4 + 3 * mt4 * mfb4
+ 3 * mb2 * mfb2 * mffb2 * mffb2 + 3 * mt2 * mfb2 * mffb2 * mffb2
- 3 * mt2 * mfb4 * mffb2 + 3 * mt2 * mf2 * mffb2 * mffb2
- 3 * mt2 * mf4 * mffb2 - 3 * mt4 * mfb2 * mffb2
- 3 * mt4 * mf2 * mffb2 - 6 * mt4 * mf2 * mfb2
+ 6 * mt2 * mf2 * mfb2 * mffb2 + 12 * mt2 * mf2 * mw4
+ 12 * mb2 * mfb2 * mw4 + 12 * mb2 * mt2 * mw4
+ 6 * mw2 * mt2 * mfb2 * mbf2 - 12 * mw2 * mt2 * mf2 * mffb2
- 6 * mw2 * mt2 * mf2 * mbf2 - 12 * mw2 * mt2 * mf2 * mfb2
- 12 * mw2 * mb2 * mfb2 * mffb2 - 6 * mw2 * mb2 * mfb2 * mbf2
+ 6 * mw2 * mb2 * mf2 * mbf2 - 12 * mw2 * mb2 * mf2 * mfb2
- 12 * mw2 * mb2 * mt2 * mfb2 - 12 * mw2 * mb2 * mt2 * mf2
+ 12 * mf2 * mfb2 * mw4 + 4 * mbf2 * mbf2 * mw4
- 6 * mfb2 * mbf2 * mw4 - 6 * mf2 * mbf2 * mw4
- 6 * mt2 * mbf2 * mw4 - 6 * mb2 * mbf2 * mw4
+ 12 * mw2 * mt2 * mf4 + 12 * mw2 * mt4 * mf2
+ 12 * mw2 * mb2 * mfb4 + 12 * mw2 * mb4 * mfb2) /mw4 / 3.;
}
void SMTopDecayer::initializeMECorrection(ShowerTreePtr tree, double & initial,
double & final) {
// check the outgoing particles
map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator cit;
ShowerParticlePtr part[2];
unsigned int ix(0);
for(cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit) {
part[ix]=cit->first->progenitor();
++ix;
}
// check the final-state particles and get the masses
if(abs(part[0]->id())==ParticleID::Wplus&&abs(part[1]->id())==ParticleID::b) {
_ma=part[0]->mass();
_mc=part[1]->mass();
}
else if(abs(part[1]->id())==ParticleID::Wplus&&abs(part[0]->id())==ParticleID::b) {
_ma=part[1]->mass();
_mc=part[0]->mass();
}
else {
return;
}
// set the top mass
_mt=tree->incomingLines().begin()->first->progenitor()->mass();
// set the gluon mass
_mg=getParticleData(ParticleID::g)->constituentMass();
// set the radiation enhancement factors
initial = _initialenhance;
final = _finalenhance;
// reduced mass parameters
_a=sqr(_ma/_mt);
_g=sqr(_mg/_mt);
_c=sqr(_mc/_mt);
useMe();
}
void SMTopDecayer::applyHardMatrixElementCorrection(ShowerTreePtr tree) {
// Get b and a and put them in particle vector ba in that order...
ParticleVector ba;
map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator cit;
for(cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit)
ba.push_back(cit->first->copy());
PPtr temp;
if(abs(ba[0]->id())!=5) swap(ba[0],ba[1]);
// Get the starting scales for the showers $\tilde{\kappa}_{b}$
// $\tilde{\kappa}_{c}$. These are needed in order to know the boundary
// of the dead zone.
map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator cjt;
- for(cjt = tree->incomingLines().begin();
- cjt!= tree->incomingLines().end();++cjt) {
- if(abs(cjt->first->progenitor()->id())==6)
- _ktb=sqr(cjt->first->progenitor()->
- evolutionScale(true,cjt->first->progenitor()->id()>0 ?
- ShowerPartnerType::QCDColourLine :
- ShowerPartnerType::QCDAntiColourLine)/_mt);
- }
- for(cit = tree->outgoingLines().begin();
- cit!= tree->outgoingLines().end();++cit) {
- if(abs(cit->first->progenitor()->id())==5)
- _ktc=sqr(cit->first->progenitor()->
- evolutionScale(true,cit->first->progenitor()->id()>0 ?
- ShowerPartnerType::QCDColourLine :
- ShowerPartnerType::QCDAntiColourLine)/_mt);
- }
+ // for(cjt = tree->incomingLines().begin();
+ // cjt!= tree->incomingLines().end();++cjt) {
+ // if(abs(cjt->first->progenitor()->id())==6)
+ // _ktb=sqr(cjt->first->progenitor()->
+ // evolutionScale(true,cjt->first->progenitor()->id()>0 ?
+ // ShowerPartnerType::QCDColourLine :
+ // ShowerPartnerType::QCDAntiColourLine)/_mt);
+ // }
+ // for(cit = tree->outgoingLines().begin();
+ // cit!= tree->outgoingLines().end();++cit) {
+ // if(abs(cit->first->progenitor()->id())==5)
+ // _ktc=sqr(cit->first->progenitor()->
+ // evolutionScale(true,cit->first->progenitor()->id()>0 ?
+ // ShowerPartnerType::QCDColourLine :
+ // ShowerPartnerType::QCDAntiColourLine)/_mt);
+ // }
+ assert(false);
if (_ktb<=0.||_ktc<=0.) {
throw Exception()
<< "SMTopDecayer::applyHardMatrixElementCorrection()"
<< " did not set ktb,ktc"
<< Exception::abortnow;
}
// Now decide if we get an emission into the dead region.
// If there is an emission newfs stores momenta for a,c,g
// according to NLO decay matrix element.
vector<Lorentz5Momentum> newfs = applyHard(ba,_ktb,_ktc);
// If there was no gluon emitted return.
if(newfs.size()!=3) return;
// Sanity checks to ensure energy greater than mass etc :)
bool check = true;
tcPDPtr gluondata=getParticleData(ParticleID::g);
if (newfs[0].e()<ba[0]->data().constituentMass()) check = false;
if (newfs[1].e()<ba[1]->mass()) check = false;
if (newfs[2].e()<gluondata->constituentMass()) check = false;
// Return if insane:
if (!check) return;
// Set masses in 5-vectors:
newfs[0].setMass(ba[0]->mass());
newfs[1].setMass(ba[1]->mass());
newfs[2].setMass(ZERO);
// The next part of this routine sets the colour structure.
// To do this for decays we assume that the gluon comes from c!
// First create new particle objects for c, a and gluon:
PPtr newg = gluondata->produceParticle(newfs[2]);
PPtr newc,newa;
newc = ba[0]->data().produceParticle(newfs[0]);
newa = ba[1];
newa->set5Momentum(newfs[1]);
// set the colour lines
ColinePtr col;
if(ba[0]->id()>0) {
col=ba[0]->colourLine();
col->addColoured(newg);
newg->colourNeighbour(newc);
}
else {
col=ba[0]->antiColourLine();
col->addAntiColoured(newg);
newg->antiColourNeighbour(newc);
}
// change the existing quark and antiquark
PPtr orig;
for(cit=tree->outgoingLines().begin();cit!=tree->outgoingLines().end();++cit) {
if(cit->first->progenitor()->id()==newc->id()) {
// remove old particles from colour line
if(newc->id()>0) {
col->removeColoured(cit->first->copy());
col->removeColoured(cit->first->progenitor());
}
else {
col->removeAntiColoured(cit->first->copy());
col->removeAntiColoured(cit->first->progenitor());
}
// insert new particles
cit->first->copy(newc);
ShowerParticlePtr sp(new_ptr(ShowerParticle(*newc,2,true)));
cit->first->progenitor(sp);
tree->outgoingLines()[cit->first]=sp;
cit->first->perturbative(false);
orig=cit->first->original();
}
else {
cit->first->copy(newa);
ShowerParticlePtr sp(new_ptr(ShowerParticle(*newa,2,true)));
map<tShowerTreePtr,pair<tShowerProgenitorPtr,
tShowerParticlePtr> >::const_iterator tit;
for(tit = tree->treelinks().begin();
tit != tree->treelinks().end();++tit) {
if(tit->second.first && tit->second.second==cit->first->progenitor())
break;
}
cit->first->progenitor(sp);
if(tit!=tree->treelinks().end())
tree->updateLink(tit->first,make_pair(cit->first,sp));
tree->outgoingLines()[cit->first]=sp;
cit->first->perturbative(true);
}
}
// Add the gluon to the shower:
ShowerParticlePtr sg =new_ptr(ShowerParticle(*newg,2,true));
ShowerProgenitorPtr gluon=new_ptr(ShowerProgenitor(orig,newg,sg));
gluon->perturbative(false);
tree->outgoingLines().insert(make_pair(gluon,sg));
if(!inTheDeadRegion(_xg,_xa,_ktb,_ktc)) {
generator()->log()
<< "SMTopDecayer::applyHardMatrixElementCorrection()\n"
<< "Just found a point that escaped from the dead region!\n"
<< " _xg: " << _xg << " _xa: " << _xa
<< " newfs.size(): " << newfs.size() << endl;
}
tree->hardMatrixElementCorrection(true);
}
vector<Lorentz5Momentum> SMTopDecayer::
applyHard(const ParticleVector &p,double ktb, double ktc) {
// ********************************* //
// First we see if we get a dead //
// region event: _xa,_xg //
// ********************************* //
vector<Lorentz5Momentum> fs;
// Return if there is no (NLO) gluon emission:
double weight = getHard(ktb,ktc);
if(weight>1.) {
generator()->log() << "Weight greater than 1 for hard emission in "
<< "SMTopDecayer::applyHard xg = " << _xg
<< " xa = " << _xa << "\n";
weight=1.;
}
// Accept/Reject
if (weight<UseRandom::rnd()||p.size()!= 2) return fs;
// Drop events if getHard returned a negative weight
// as in events that, somehow have escaped from the dead region
// or, worse, the allowed region.
if(weight<0.) return fs;
// Calculate xc by momentum conservation:
_xc = 2.-_xa-_xg;
// ************************************ //
// Now we get the boosts & rotations to //
// go from lab to top rest frame with //
// a in the +z direction. //
// ************************************ //
Lorentz5Momentum pa_lab,pb_lab,pc_lab,pg_lab;
// Calculate momentum of b:
pb_lab = p[0]->momentum() + p[1]->momentum();
// Define/assign momenta of c,a and the gluon:
if(abs(p[0]->id())==5) {
pc_lab = p[0]->momentum();
pa_lab = p[1]->momentum();
} else {
pc_lab = p[1]->momentum();
pa_lab = p[0]->momentum();
}
// Calculate the boost to the b rest frame:
SpinOneLorentzRotation rot0(pb_lab.findBoostToCM());
// Calculate the rotation matrix to position a along the +z direction
// in the rest frame of b and does a random rotation about z:
SpinOneLorentzRotation rot1 = rotateToZ(rot0*pa_lab);
// Calculate the boost from the b rest frame back to the lab:
// and the inverse of the random rotation about the z-axis and the
// rotation required to align a with +z:
SpinOneLorentzRotation invrot = rot0.inverse()*rot1.inverse();
// ************************************ //
// Now we construct the momenta in the //
// b rest frame using _xa,_xg. //
// First we construct b, then c and g, //
// finally we generate a by momentum //
// conservation. //
// ************************************ //
Lorentz5Momentum pa_brf, pb_brf(_mt), pc_brf, pg_brf;
// First we set the top quark to being on-shell and at rest.
// Second we set the energies of c and g,
pc_brf.setE(0.5*_mt*(2.-_xa-_xg));
pg_brf.setE(0.5*_mt*_xg);
// then their masses,
pc_brf.setMass(_mc);
pg_brf.setMass(ZERO);
// Now set the z-component of c and g. For pg we simply start from
// _xa and _xg, while for pc we assume it is equal to minus the sum
// of the z-components of a (assumed to point in the +z direction) and g.
double root=sqrt(_xa*_xa-4.*_a);
pg_brf.setZ(_mt*(1.-_xa-_xg+0.5*_xa*_xg-_c+_a)/root);
pc_brf.setZ(-1.*( pg_brf.z()+_mt*0.5*root));
// Now set the y-component of c and g's momenta
pc_brf.setY(ZERO);
pg_brf.setY(ZERO);
// Now set the x-component of c and g's momenta
pg_brf.setX(sqrt(sqr(pg_brf.t())-sqr(pg_brf.z())));
pc_brf.setX(-pg_brf.x());
// Momenta b,c,g are now set. Now we obtain a from momentum conservation,
pa_brf = pb_brf-pc_brf-pg_brf;
pa_brf.setMass(pa_brf.m());
pa_brf.rescaleEnergy();
// ************************************ //
// Now we orient the momenta and boost //
// them back to the original lab frame. //
// ************************************ //
// As in herwig6507 we assume that, in the rest frame
// of b, we have aligned the W boson momentum in the
// +Z direction by rot1*rot0*pa_lab, therefore
// we obtain the new pa_lab by applying:
// invrot*pa_brf.
pa_lab = invrot*pa_brf;
pb_lab = invrot*pb_brf;
pc_lab = invrot*pc_brf;
pg_lab = invrot*pg_brf;
fs.push_back(pc_lab);
fs.push_back(pa_lab);
fs.push_back(pg_lab);
return fs;
}
double SMTopDecayer::getHard(double ktb, double ktc) {
// zero the variables
_xg = 0.;
_xa = 0.;
_xc = 0.;
// Get a phase space point in the dead region:
double volume_factor = deadRegionxgxa(ktb,ktc);
// if outside region return -1
if(volume_factor<0) return volume_factor;
// Compute the weight for this phase space point:
double weight = volume_factor*me(_xa,_xg)*(1.+_a-_c-_xa);
// Alpha_S and colour factors - this hard wired Alpha_S needs removing.
weight *= (4./3.)/Constants::pi
*(_alpha->value(_mt*_mt*_xg*(1.-_xa+_a-_c)
/(2.-_xg-_xa-_c)));
return weight;
}
bool SMTopDecayer::softMatrixElementVeto(ShowerProgenitorPtr initial,
ShowerParticlePtr parent,Branching br) {
// check if we need to apply the full correction
long id[2]={abs(initial->progenitor()->id()),abs(parent->id())};
// the initial-state correction
if(id[0]==ParticleID::t&&id[1]==ParticleID::t)
{
Energy pt=br.kinematics->pT();
// check if hardest so far
// if not just need to remove effect of enhancement
bool veto(false);
// if not hardest so far
if(pt<initial->highestpT())
veto=!UseRandom::rndbool(1./_initialenhance);
// if hardest so far do calculation
else
{
// values of kappa and z
double z(br.kinematics->z()),kappa(sqr(br.kinematics->scale()/_mt));
// parameters for the translation
double w(1.-(1.-z)*(kappa-1.)),u(1.+_a-_c-(1.-z)*kappa),v(sqr(u)-4.*_a*w*z);
// veto if outside phase space
if(v<0.)
veto=true;
// otherwise calculate the weight
else
{
v = sqrt(v);
double xa((0.5*(u+v)/w+0.5*(u-v)/z)),xg((1.-z)*kappa);
double f(me(xa,xg)),
J(0.5*(u+v)/sqr(w)-0.5*(u-v)/sqr(z)+_a*sqr(w-z)/(v*w*z));
double wgt(f*J*2./kappa/(1.+sqr(z)-2.*z/kappa)/_initialenhance);
// This next `if' prevents the hardest emission from the
// top shower ever entering the so-called T2 region of the
// phase space if that region is to be populated by the hard MEC.
if(_useMEforT2&&xg>xgbcut(_ktb)) wgt = 0.;
if(wgt>1.) {
generator()->log() << "Violation of maximum for initial-state "
<< " soft veto in "
<< "SMTopDecayer::softMatrixElementVeto"
<< "xg = " << xg << " xa = " << xa
<< "weight = " << wgt << "\n";
wgt=1.;
}
// compute veto from weight
veto = !UseRandom::rndbool(wgt);
}
// if not vetoed reset max
if(!veto) initial->highestpT(pt);
}
// if vetoing reset the scale
if(veto) parent->vetoEmission(br.type,br.kinematics->scale());
// return the veto
return veto;
}
// final-state correction
else if(id[0]==ParticleID::b&&id[1]==ParticleID::b)
{
Energy pt=br.kinematics->pT();
// check if hardest so far
// if not just need to remove effect of enhancement
bool veto(false);
// if not hardest so far
if(pt<initial->highestpT()) return !UseRandom::rndbool(1./_finalenhance);
// if hardest so far do calculation
// values of kappa and z
double z(br.kinematics->z()),kappa(sqr(br.kinematics->scale()/_mt));
// momentum fractions
double xa(1.+_a-_c-z*(1.-z)*kappa),r(0.5*(1.+_c/(1.+_a-xa))),root(sqr(xa)-4.*_a);
if(root<0.) {
generator()->log() << "Imaginary root for final-state veto in "
<< "SMTopDecayer::softMatrixElementVeto"
<< "\nz = " << z << "\nkappa = " << kappa
<< "\nxa = " << xa
<< "\nroot^2= " << root;
parent->vetoEmission(br.type,br.kinematics->scale());
return true;
}
root=sqrt(root);
double xg((2.-xa)*(1.-r)-(z-r)*root);
// xfact (below) is supposed to equal xg/(1-z).
double xfact(z*kappa/2./(z*(1.-z)*kappa+_c)*(2.-xa-root)+root);
// calculate the full result
double f(me(xa,xg));
// jacobian
double J(z*root);
double wgt(f*J*2.*kappa/(1.+sqr(z)-2.*_c/kappa/z)/sqr(xfact)/_finalenhance);
if(wgt>1.) {
generator()->log() << "Violation of maximum for final-state soft veto in "
<< "SMTopDecayer::softMatrixElementVeto"
<< "xg = " << xg << " xa = " << xa
<< "weight = " << wgt << "\n";
wgt=1.;
}
// compute veto from weight
veto = !UseRandom::rndbool(wgt);
// if not vetoed reset max
if(!veto) initial->highestpT(pt);
// if vetoing reset the scale
if(veto) parent->vetoEmission(br.type,br.kinematics->scale());
// return the veto
return veto;
}
// otherwise don't veto
else return !UseRandom::rndbool(1./_finalenhance);
}
double SMTopDecayer::me(double xw,double xg) {
double prop(1.+_a-_c-xw),xg2(sqr(xg));
double lambda=sqrt(1.+_a*_a+_c*_c-2.*_a-2.*_c-2.*_a*_c);
double denom=(1.-2*_a*_a+_a+_c*_a+_c*_c-2.*_c);
double wgt=-_c*xg2/prop+(1.-_a+_c)*xg-(prop*(1 - xg)+xg2)
+(0.5*(1.+2.*_a+_c)*sqr(prop-xg)*xg+2.*_a*prop*xg2)/denom;
return wgt/(lambda*prop);
}
// This function is auxiliary to the xab function.
double SMTopDecayer::xgbr(int toggle) {
return 1.+toggle*sqrt(_a)-_c*(1.-toggle*sqrt(_a))/(1.-_a);
}
// This function is auxiliary to the xab function.
double SMTopDecayer::ktr(double xgb, int toggle) {
return 2.*xgb/
(xgb+toggle*sqrt((1.-1./_a)
*(xgb-xgbr( 1))
*(xgb-xgbr(-1))));
}
// Function xab determines xa (2*W energy fraction) for a given value
// of xg (2*gluon energy fraction) and kappa tilde (q tilde squared over
// m_top squared). Hence this function allows you to draw 1: the total
// phase space volume in the xa vs xg plane 2: for a given value of
// kappa tilde (i.e. starting evolution scale) the associated contour
// in the xa vs xg plane (and hence the regions that either shower can
// populate). This calculation is done assuming the emission came from
// the top quark i.e. kappa tilde here is the q tilde squared of the TOP
// quark divided by m_top squared.
double SMTopDecayer::xab(double xgb, double kt, int toggle) {
double xab;
if(toggle==2) {
// This applies for g==0.&&kt==ktr(a,c,0.,xgb,1).
xab = -2.*_a*(xgb-2.)/(1.+_a-_c-xgb);
} else if(toggle==1) {
// This applies for kt==1&&g==0.
double lambda = sqrt(sqr(xgb-1.+_a+_c)-4.*_a*_c);
xab = (0.5/(kt-xgb))*(kt*(1.+_a-_c-xgb)-lambda)
+ (0.5/(kt+xgb*(1.-kt)))*(kt*(1.+_a-_c-xgb)+lambda);
} else {
// This is the form of xab FOR _g=0.
double ktmktrpktmktrm = kt*kt - 4.*_a*(kt-1.)*xgb*xgb
/ (sqr(1.-_a-_c-xgb)-4.*_a*_c);
if(fabs(kt-(2.*xgb-2.*_g)/(xgb-sqrt(xgb*xgb-4.*_g)))/kt>1.e-6) {
double lambda = sqrt((sqr(1.-_a-_c-xgb)-4.*_a*_c)*ktmktrpktmktrm);
xab = (0.5/(kt-xgb))*(kt*(1.+_a-_c-xgb)-lambda)
+ (0.5/(kt+xgb*(1.-kt)))*(kt*(1.+_a-_c-xgb)+lambda);
}
else {
// This is the value of xa as a function of xb when kt->infinity.
// Where we take any kt > (2.*xgb-2.*_g)/(xgb-sqrt(xgb*xgb-4.*_g))
// as being effectively infinite. This kt value is actually the
// maximum allowed value kt can have if the phase space is calculated
// without the approximation of _g=0 (massless gluon). This formula
// for xab below is then valid for _g=0 AND kt=infinity only.
xab = ( 2.*_c+_a*(xgb-2.)
+ 3.*xgb
- xgb*(_c+xgb+sqrt(_a*_a-2.*(_c-xgb+1.)*_a+sqr(_c+xgb-1.)))
- 2.
)/2./(xgb-1.);
}
}
if(isnan(xab)) {
double ktmktrpktmktrm = ( sqr(xgb*kt-2.*(xgb-_g))
-kt*kt*(1.-1./_a)*(xgb-xgbr( 1)-_g/(1.+sqrt(_a)))
*(xgb-xgbr(-1)-_g/(1.-sqrt(_a)))
)/
(xgb*xgb-(1.-1./_a)*(xgb-xgbr( 1)-_g/(1.+sqrt(_a)))
*(xgb-xgbr(-1)-_g/(1.-sqrt(_a)))
);
double lambda = sqrt((xgb-1.+sqr(sqrt(_a)+sqrt(_c-_g)))
*(xgb-1.+sqr(sqrt(_a)-sqrt(_c-_g)))*
ktmktrpktmktrm);
xab = (0.5/(kt-xgb+_g))*(kt*(1.+_a-_c+_g-xgb)-lambda)
+ (0.5/(kt+xgb*(1.-kt)-_g))*(kt*(1.+_a-_c+_g-xgb)+lambda);
if(isnan(xab))
throw Exception() << "TopMECorrection::xab complex x_a value.\n"
<< " xgb = " << xgb << "\n"
<< " xab = " << xab << "\n"
<< " toggle = " << toggle << "\n"
<< " ktmktrpktmktrm = " << ktmktrpktmktrm
<< Exception::eventerror;
}
return xab;
}
// xgbcut is the point along the xg axis where the upper bound on the
// top quark (i.e. b) emission phase space goes back on itself in the
// xa vs xg plane i.e. roughly mid-way along the xg axis in
// the xa vs xg Dalitz plot.
double SMTopDecayer::xgbcut(double kt) {
double lambda2 = 1.+_a*_a+_c*_c-2.*_a-2.*_c-2.*_a*_c;
double num1 = kt*kt*(1.-_a-_c);
double num2 = 2.*kt*sqrt(_a*(kt*kt*_c+lambda2*(kt-1.)));
return (num1-num2)/(kt*kt-4.*_a*(kt-1.));
}
double SMTopDecayer::xaccut(double kt) {
return 1.+_a-_c-0.25*kt;
}
double SMTopDecayer::z(double xac, double kt,
int toggle1, int toggle2) {
double z = -1.0;
if(toggle2==0) {
z = (kt+toggle1*sqrt(kt*(kt-4.*(1.+_a-_c-xac))))/(2.*kt);
} else if(toggle2==1) {
z = ((1.+_a+_c-xac)+toggle1*(1.+_a-_c-xac))
/(2.*(1.+_a-xac));
} else if(toggle2==2) {
z = 0.5;
} else {
throw Exception() << "Cannot determine z in SMTopDecayer::z()"
<< Exception::eventerror;
}
return z;
}
double SMTopDecayer::xgc(double xac, double kt,
int toggle1, int toggle2) {
double tiny(1.e-6);
double xaToMinBoundary(xac*xac-4.*_a);
if(xaToMinBoundary<0) {
if(fabs(xaToMinBoundary/(1.-_a)/(1.-_a))<tiny)
xaToMinBoundary *= -1.;
else
throw Exception() << "SMTopDecayer::xgc xa not in phase space!"
<< Exception::eventerror;
}
return (2.-xac)*(1.-0.5*(1.+_c/(1.+_a-xac)))
-(z(xac,kt,toggle1,toggle2)-0.5*(1.+_c/(1.+_a-xac)))
*sqrt(xaToMinBoundary);
}
double SMTopDecayer::xginvc0(double xg , double kt) {
// The function xg(kappa_tilde_c,xa) surely, enough, draws a
// line of constant kappa_tilde_c in the xg, xa Dalitz plot.
// Such a function can therefore draw the upper and lower
// edges of the phase space for emission from c (the b-quark).
// However, to sample the soft part of the dead zone effectively
// we want to generate a value of xg first and THEN distribute
// xa in the associated allowed part of the dead zone. Hence, the
// function we want, to define the dead zone in xa for a given
// xg, is the inverse of xg(kappa_tilde_c,xa). The full expression
// for xg(kappa_tilde_c,xa) is complicated and, sure enough,
// does not invert. Therefore we try to overestimate the size
// of the dead zone initially, rejecting events which do not
// fall exactly inside it afterwards, with the immediate aim
// of getting an approximate version of xg(kappa_tilde_c,xa)
// that can be inverted. We do this by simply setting c=0 i.e.
// the b-quark mass to zero (and the gluon mass of course), in
// the full expression xg(...). The result of inverting this
// function is the output of this routine (a value of xa) hence
// the name xginvc0. xginvc0 is calculated to be,
// xginvc0 = (1./3.)*(1.+a+pow((U+sqrt(4.*V*V*V+U*U))/2.,1./3.)
// -V*pow(2./(U+sqrt(4.*V*V*V+U*U)),1./3.)
// )
// U = 2.*a*a*a - 66.*a*a + 9.*a*kt*xg + 18.*a*kt
// - 66.*a + 27.*kt*xg*xg - 45.*kt*xg +18.*kt +2. ;
// V = -1.-a*a-14.*a-3.kt*xg+3.*kt;
// This function, as with many functions in this ME correction,
// is plagued by cuts that have to handled carefully in numerical
// implementation. We have analysed the cuts and hence we implement
// it in the following way, with a series of 'if' statements.
//
// A useful -definition- to know in deriving the v<0 terms is
// that tanh^-1(z) = 0.5*(log(1.+z)-log(1.-z)).
double u,v,output;
u = 2.*_a*_a*_a-66.*_a*_a
+9.*xg*kt*_a+18.*kt*_a
-66.*_a+27.*xg*xg*kt
-45.*xg*kt+18.*kt+2.;
v = -_a*_a-14.*_a-3.*xg*kt+3.*kt-1.;
double u2=u*u,v3=v*v*v;
if(v<0.) {
if(u>0.&&(4.*v3+u2)<0.) output = cos( atan(sqrt(-4.*v3-u2)/u)/3.);
else if(u>0.&&(4.*v3+u2)>0.) output = cosh(atanh(sqrt( 4.*v3+u2)/u)/3.);
else output = cos(( atan(sqrt(-4.*v3-u2)/u)
+Constants::pi)/3.);
output *= 2.*sqrt(-v);
} else {
output = sinh(log((u+sqrt(4.*v3+u2))/(2.*sqrt(v3)))/3.);
output *= 2.*sqrt(v);
}
if(isnan(output)||isinf(output)) {
throw Exception() << "TopMECorrection::xginvc0:\n"
<< "possible numerical instability detected.\n"
<< "\n v = " << v << " u = " << u << "\n4.*v3+u2 = " << 4.*v3+u2
<< "\n_a = " << _a << " ma = " << sqrt(_a*_mt*_mt/GeV2)
<< "\n_c = " << _c << " mc = " << sqrt(_c*_mt*_mt/GeV2)
<< "\n_g = " << _g << " mg = " << sqrt(_g*_mt*_mt/GeV2)
<< Exception::eventerror;
}
return ( 1.+_a +output)/3.;
}
double SMTopDecayer::approxDeadMaxxa(double xg,double ktb,double ktc) {
double maxxa(0.);
double x = min(xginvc0(xg,ktc),
xab(xg,(2.*xg-2.*_g)/(xg-sqrt(xg*xg-4.*_g)),0));
double y(-9999999999.);
if(xg>2.*sqrt(_g)&&xg<=xgbcut(ktb)) {
y = max(xab(xg,ktb,0),xab(xg,1.,1));
} else if(xg>=xgbcut(ktb)&&xg<=1.-sqr(sqrt(_a)+sqrt(_c))) {
y = max(xab(xg,ktr(xg,1),2),xab(xg,1.,1));
}
if(xg>2.*sqrt(_g)&&xg<=1.-sqr(sqrt(_a)+sqrt(_c))) {
if(x>=y) { maxxa = x ; }
else { maxxa = -9999999.; }
} else {
maxxa = -9999999.;
}
return maxxa;
}
double SMTopDecayer::approxDeadMinxa(double xg,double ktb,double ktc) {
double minxa(0.);
double x = min(xginvc0(xg,ktc),
xab(xg,(2.*xg-2.*_g)/(xg-sqrt(xg*xg-4.*_g)),0));
double y(-9999999999.);
if(xg>2.*sqrt(_g)&&xg<=xgbcut(ktb)) {
y = max(xab(xg,ktb,0),xab(xg,1.,1));
} else if(xg>=xgbcut(ktb)&&xg<=1.-sqr(sqrt(_a)+sqrt(_c))) {
if(_useMEforT2) y = xab(xg,1.,1);
else y = max(xab(xg,ktr(xg,1),2),xab(xg,1.,1));
}
if(xg>2.*sqrt(_g)&&xg<=1.-sqr(sqrt(_a)+sqrt(_c))) {
if(x>=y) { minxa = y ; }
else { minxa = 9999999.; }
} else {
minxa = 9999999.;
}
return minxa;
}
// This function returns true if the phase space point (xg,xa) is in the
// kinematically allowed phase space.
bool SMTopDecayer::inTheAllowedRegion(double xg , double xa) {
bool output(true);
if(xg<2.*sqrt(_g)||xg>1.-sqr(sqrt(_a)+sqrt(_c))) output = false;
if(xa<xab(xg,1.,1)) output = false;
if(xa>xab(xg,(2.*xg-2.*_g)/(xg-sqrt(xg*xg-4.*_g)),0)) output = false;
return output;
}
// This function returns true if the phase space point (xg,xa) is in the
// approximate (overestimated) dead region.
bool SMTopDecayer::inTheApproxDeadRegion(double xg , double xa,
double ktb, double ktc) {
bool output(true);
if(!inTheAllowedRegion(xg,xa)) output = false;
if(xa<approxDeadMinxa(xg,ktb,ktc)) output = false;
if(xa>approxDeadMaxxa(xg,ktb,ktc)) output = false;
return output;
}
// This function returns true if the phase space point (xg,xa) is in the
// dead region.
bool SMTopDecayer::inTheDeadRegion(double xg , double xa,
double ktb, double ktc) {
bool output(true);
if(!inTheApproxDeadRegion(xg,xa,ktb,ktc)) output = false;
if(xa>xaccut(ktc)) {
if(xg<xgc(max(xaccut(ktc),2.*sqrt(_a)),ktc, 1,2)&&
xg>xgc(xa,ktc, 1,0)) { output = false; }
if(xg>xgc(max(xaccut(ktc),2.*sqrt(_a)),ktc,-1,2)&&
xg<xgc(xa,ktc,-1,0)) { output = false; }
}
return output;
}
// This function attempts to generate a phase space point in the dead
// region and returns the associated phase space volume factor needed for
// the associated event weight.
double SMTopDecayer::deadRegionxgxa(double ktb,double ktc) {
_xg=0.;
_xa=0.;
// Here we set limits on xg and generate a value inside the bounds.
double xgmin(2.*sqrt(_g)),xgmax(1.-sqr(sqrt(_a)+sqrt(_c)));
// Generate _xg.
if(_xg_sampling==2.) {
_xg=xgmin*xgmax/(xgmin+UseRandom::rnd()*(xgmax-xgmin));
} else {
_xg=xgmin*xgmax/pow(( pow(xgmin,_xg_sampling-1.)
+ UseRandom::rnd()*(pow(xgmax,_xg_sampling-1.)
-pow(xgmin,_xg_sampling-1.))
),1./(_xg_sampling-1.));
}
// Here we set the bounds on _xa for given _xg.
if(_xg<xgmin||xgmin>xgmax)
throw Exception() << "TopMECorrection::deadRegionxgxa:\n"
<< "upper xg bound is less than the lower xg bound.\n"
<< "\n_xg = " << _xg
<< "\n2.*sqrt(_g) = " << 2.*sqrt(_g)
<< "\n_a = " << _a << " ma = " << sqrt(_a*_mt*_mt/GeV2)
<< "\n_c = " << _c << " mc = " << sqrt(_c*_mt*_mt/GeV2)
<< "\n_g = " << _g << " mg = " << sqrt(_g*_mt*_mt/GeV2)
<< Exception::eventerror;
double xamin(approxDeadMinxa(_xg,ktb,ktc));
double xamax(approxDeadMaxxa(_xg,ktb,ktc));
// Are the bounds sensible? If not return.
if(xamax<=xamin) return -1.;
_xa=1.+_a-(1.+_a-xamax)*pow((1.+_a-xamin)/(1.+_a-xamax),UseRandom::rnd());
// If outside the allowed region return -1.
if(!inTheDeadRegion(_xg,_xa,ktb,ktc)) return -1.;
// The integration volume for the weight
double xg_vol,xa_vol;
if(_xg_sampling==2.) {
xg_vol = (xgmax-xgmin)
/ (xgmax*xgmin);
} else {
xg_vol = (pow(xgmax,_xg_sampling-1.)-pow(xgmin,_xg_sampling-1.))
/ ((_xg_sampling-1.)*pow(xgmax*xgmin,_xg_sampling-1.));
}
xa_vol = log((1.+_a-xamin)/(1.+_a-xamax));
// Here we return the integral volume factor multiplied by the part of the
// weight left over which is not included in the BRACES function, i.e.
// the part of _xg^-2 which is not absorbed in the integration measure.
return xg_vol*xa_vol*pow(_xg,_xg_sampling-2.);
}
LorentzRotation SMTopDecayer::rotateToZ(Lorentz5Momentum v) {
// compute the rotation matrix
LorentzRotation trans;
// rotate so in z-y plane
trans.rotateZ(-atan2(v.y(),v.x()));
// rotate so along Z
trans.rotateY(-acos(v.z()/v.vect().mag()));
// generate random rotation
double c,s,cs;
do
{
c = 2.*UseRandom::rnd()-1.;
s = 2.*UseRandom::rnd()-1.;
cs = c*c+s*s;
}
while(cs>1.||cs==0.);
double cost=(c*c-s*s)/cs,sint=2.*c*s/cs;
// apply random azimuthal rotation
trans.rotateZ(atan2(sint,cost));
return trans;
}
diff --git a/Shower/Base/Evolver.cc b/Shower/Base/Evolver.cc
--- a/Shower/Base/Evolver.cc
+++ b/Shower/Base/Evolver.cc
@@ -1,2047 +1,2051 @@
// -*- C++ -*-
//
// Evolver.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 Evolver class.
//
#include "Evolver.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/Shower/Base/ShowerParticle.h"
#include "ThePEG/Utilities/EnumIO.h"
#include "ShowerKinematics.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Handlers/EventHandler.h"
#include "ThePEG/Utilities/Throw.h"
#include "ShowerTree.h"
#include "ShowerProgenitor.h"
#include "KinematicsReconstructor.h"
#include "PartnerFinder.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include "ThePEG/PDT/DecayMode.h"
#include "Herwig++/Shower/ShowerHandler.h"
using namespace Herwig;
namespace {
void findChildren(tShowerParticlePtr parent,set<ShowerParticlePtr> & fs) {
for(unsigned int ix=0;ix<parent->children().size();++ix) {
tShowerParticlePtr child=
dynamic_ptr_cast<tShowerParticlePtr>(parent->children()[ix]);
if(child) findChildren(child,fs);
}
if(parent->children().empty()) {
if(parent->isFinalState()) fs.insert(parent);
}
}
}
IBPtr Evolver::clone() const {
return new_ptr(*this);
}
IBPtr Evolver::fullclone() const {
return new_ptr(*this);
}
void Evolver::persistentOutput(PersistentOStream & os) const {
os << _model << _splittingGenerator << _maxtry
<< _meCorrMode << _hardVetoMode << _hardVetoRead << _hardVetoReadOption
<< _limitEmissions
<< ounit(_iptrms,GeV) << _beta << ounit(_gamma,GeV) << ounit(_iptmax,GeV)
<< _vetoes << _trunc_Mode << _hardEmissionMode
<< _colourEvolutionMethod << _reconOpt
<< interaction_<< interactions_.size();
for(unsigned int ix=0;ix<interactions_.size();++ix)
os << oenum(interactions_[ix]);
}
void Evolver::persistentInput(PersistentIStream & is, int) {
unsigned int isize;
is >> _model >> _splittingGenerator >> _maxtry
>> _meCorrMode >> _hardVetoMode >> _hardVetoRead >> _hardVetoReadOption
>> _limitEmissions
>> iunit(_iptrms,GeV) >> _beta >> iunit(_gamma,GeV) >> iunit(_iptmax,GeV)
>> _vetoes >> _trunc_Mode >> _hardEmissionMode
>> _colourEvolutionMethod >> _reconOpt
>> interaction_ >> isize;
interactions_.resize(isize);
for(unsigned int ix=0;ix<interactions_.size();++ix)
is >> ienum(interactions_[ix]);
}
void Evolver::doinit() {
Interfaced::doinit();
if(interaction_==0) {
interactions_.push_back(ShowerInteraction::QCD);
interactions_.push_back(ShowerInteraction::QED);
}
else if(interaction_==1) {
interactions_.push_back(ShowerInteraction::QCD);
}
else if(interaction_==2) {
interactions_.push_back(ShowerInteraction::QED);
interactions_.push_back(ShowerInteraction::QCD);
}
else if(interaction_==3) {
interactions_.push_back(ShowerInteraction::QED);
}
else if(interaction_==4) {
interactions_.push_back(ShowerInteraction::Both);
}
}
ClassDescription<Evolver> Evolver::initEvolver;
// Definition of the static class description member.
void Evolver::Init() {
static ClassDocumentation<Evolver> documentation
("This class is responsible for carrying out the showering,",
"including the kinematics reconstruction, in a given scale range,"
"including the option of the POWHEG approach to simulated next-to-leading order"
" radiation\\cite{Nason:2004rx}.",
"%\\cite{Nason:2004rx}\n"
"\\bibitem{Nason:2004rx}\n"
" P.~Nason,\n"
" ``A new method for combining NLO QCD with shower Monte Carlo algorithms,''\n"
" JHEP {\\bf 0411} (2004) 040\n"
" [arXiv:hep-ph/0409146].\n"
" %%CITATION = JHEPA,0411,040;%%\n");
static Reference<Evolver,SplittingGenerator>
interfaceSplitGen("SplittingGenerator",
"A reference to the SplittingGenerator object",
&Herwig::Evolver::_splittingGenerator,
false, false, true, false);
static Reference<Evolver,ShowerModel> interfaceShowerModel
("ShowerModel",
"The pointer to the object which defines the shower evolution model.",
&Evolver::_model, false, false, true, false, false);
static Parameter<Evolver,unsigned int> interfaceMaxTry
("MaxTry",
"The maximum number of attempts to generate the shower from a"
" particular ShowerTree",
&Evolver::_maxtry, 100, 1, 1000,
false, false, Interface::limited);
static Switch<Evolver, unsigned int> ifaceMECorrMode
("MECorrMode",
"Choice of the ME Correction Mode",
&Evolver::_meCorrMode, 1, false, false);
static SwitchOption off
(ifaceMECorrMode,"No","MECorrections off", 0);
static SwitchOption on
(ifaceMECorrMode,"Yes","hard+soft on", 1);
static SwitchOption hard
(ifaceMECorrMode,"Hard","only hard on", 2);
static SwitchOption soft
(ifaceMECorrMode,"Soft","only soft on", 3);
static Switch<Evolver, unsigned int> ifaceHardVetoMode
("HardVetoMode",
"Choice of the Hard Veto Mode",
&Evolver::_hardVetoMode, 1, false, false);
static SwitchOption HVoff
(ifaceHardVetoMode,"No","hard vetos off", 0);
static SwitchOption HVon
(ifaceHardVetoMode,"Yes","hard vetos on", 1);
static SwitchOption HVIS
(ifaceHardVetoMode,"Initial", "only IS emissions vetoed", 2);
static SwitchOption HVFS
(ifaceHardVetoMode,"Final","only FS emissions vetoed", 3);
static Switch<Evolver, unsigned int> ifaceHardVetoRead
("HardVetoScaleSource",
"If hard veto scale is to be read",
&Evolver::_hardVetoRead, 0, false, false);
static SwitchOption HVRcalc
(ifaceHardVetoRead,"Calculate","Calculate from hard process", 0);
static SwitchOption HVRread
(ifaceHardVetoRead,"Read","Read from XComb->lastScale", 1);
static Switch<Evolver, bool> ifaceHardVetoReadOption
("HardVetoReadOption",
"Apply read-in scale veto to all collisions or just the primary one?",
&Evolver::_hardVetoReadOption, false, false, false);
static SwitchOption AllCollisions
(ifaceHardVetoReadOption,
"AllCollisions",
"Read-in pT veto applied to primary and secondary collisions.",
false);
static SwitchOption PrimaryCollision
(ifaceHardVetoReadOption,
"PrimaryCollision",
"Read-in pT veto applied to primary but not secondary collisions.",
true);
static Parameter<Evolver, Energy> ifaceiptrms
("IntrinsicPtGaussian",
"RMS of intrinsic pT of Gaussian distribution:\n"
"2*(1-Beta)*exp(-sqr(intrinsicpT/RMS))/sqr(RMS)",
&Evolver::_iptrms, GeV, ZERO, ZERO, 1000000.0*GeV,
false, false, Interface::limited);
static Parameter<Evolver, double> ifacebeta
("IntrinsicPtBeta",
"Proportion of inverse quadratic distribution in generating intrinsic pT.\n"
"(1-Beta) is the proportion of Gaussian distribution",
&Evolver::_beta, 0, 0, 1,
false, false, Interface::limited);
static Parameter<Evolver, Energy> ifacegamma
("IntrinsicPtGamma",
"Parameter for inverse quadratic:\n"
"2*Beta*Gamma/(sqr(Gamma)+sqr(intrinsicpT))",
&Evolver::_gamma,GeV, ZERO, ZERO, 100000.0*GeV,
false, false, Interface::limited);
static Parameter<Evolver, Energy> ifaceiptmax
("IntrinsicPtIptmax",
"Upper bound on intrinsic pT for inverse quadratic",
&Evolver::_iptmax,GeV, ZERO, ZERO, 100000.0*GeV,
false, false, Interface::limited);
static RefVector<Evolver,ShowerVeto> ifaceVetoes
("Vetoes",
"The vetoes to be checked during showering",
&Evolver::_vetoes, -1,
false,false,true,true,false);
static Switch<Evolver,unsigned int> interfaceLimitEmissions
("LimitEmissions",
"Limit the number and type of emissions for testing",
&Evolver::_limitEmissions, 0, false, false);
static SwitchOption interfaceLimitEmissionsNoLimit
(interfaceLimitEmissions,
"NoLimit",
"Allow an arbitrary number of emissions",
0);
static SwitchOption interfaceLimitEmissionsOneInitialStateEmission
(interfaceLimitEmissions,
"OneInitialStateEmission",
"Allow one emission in the initial state and none in the final state",
1);
static SwitchOption interfaceLimitEmissionsOneFinalStateEmission
(interfaceLimitEmissions,
"OneFinalStateEmission",
"Allow one emission in the final state and none in the initial state",
2);
static SwitchOption interfaceLimitEmissionsHardOnly
(interfaceLimitEmissions,
"HardOnly",
"Only allow radiation from the hard ME correction",
3);
static SwitchOption interfaceLimitEmissionsOneEmission
(interfaceLimitEmissions,
"OneEmission",
"Allow one emission in either the final state or initial state, but not both",
4);
static Switch<Evolver,bool> interfaceTruncMode
("TruncatedShower", "Include the truncated shower?",
&Evolver::_trunc_Mode, 1, false, false);
static SwitchOption interfaceTruncMode0
(interfaceTruncMode,"No","Truncated Shower is OFF", 0);
static SwitchOption interfaceTruncMode1
(interfaceTruncMode,"Yes","Truncated Shower is ON", 1);
static Switch<Evolver,unsigned int> interfaceHardEmissionMode
("HardEmissionMode",
"Whether to use ME corrections or POWHEG for the hardest emission",
&Evolver::_hardEmissionMode, 0, false, false);
static SwitchOption interfaceHardEmissionModeMECorrection
(interfaceHardEmissionMode,
"MECorrection",
"Old fashioned ME correction",
0);
static SwitchOption interfaceHardEmissionModePOWHEG
(interfaceHardEmissionMode,
"POWHEG",
"Powheg style hard emission",
1);
static Switch<Evolver,int> interfaceColourEvolutionMethod
("ColourEvolutionMethod",
"Choice of method for choosing the colour factor in gluon evolution",
&Evolver::_colourEvolutionMethod, 0, false, false);
static SwitchOption interfaceColourEvolutionMethodDefault
(interfaceColourEvolutionMethod,
"Default",
"Colour factor is CA for all scales",
0);
static SwitchOption interfaceColourEvolutionMethodHalfCA
(interfaceColourEvolutionMethod,
"HalfCA",
"Only use half the normal radiation until second scale is reached",
1);
static Switch<Evolver,unsigned int > interfaceInteractions
("Interactions",
"The interactions to be used in the shower",
&Evolver::interaction_, 1, false, false);
static SwitchOption interfaceInteractionsQCDFirst
(interfaceInteractions,
"QCDFirst",
"QCD first then QED",
0);
static SwitchOption interfaceInteractionsQCDOnly
(interfaceInteractions,
"QCDOnly",
"Only QCD",
1);
static SwitchOption interfaceInteractionsQEDFirst
(interfaceInteractions,
"QEDFirst",
"QED first then QCD",
2);
static SwitchOption interfaceInteractionsQEDOnly
(interfaceInteractions,
"QEDOnly",
"Only QED",
3);
static SwitchOption interfaceInteractionsBothAtOnce
(interfaceInteractions,
"BothAtOnce",
"Generate both at the same time",
4);
static Switch<Evolver,unsigned int> interfaceReconstructionOption
("ReconstructionOption",
"Treatment of the reconstruction of the transverse momentum of "
"a branching from the evolution scale.",
&Evolver::_reconOpt, 0, false, false);
static SwitchOption interfaceReconstructionOptionCutOff
(interfaceReconstructionOption,
"CutOff",
"Use the cut-off masses in the calculation",
0);
static SwitchOption interfaceReconstructionOptionOffShell
(interfaceReconstructionOption,
"OffShell",
"Use the off-shell masses in the calculation",
1);
}
void Evolver::generateIntrinsicpT(vector<ShowerProgenitorPtr> particlesToShower) {
_intrinsic.clear();
if ( !ipTon() || !isISRadiationON() ) return;
// don't do anything for the moment for secondary scatters
if( !ShowerHandler::currentHandler()->firstInteraction() ) return;
// generate intrinsic pT
for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
// only consider initial-state particles
if(particlesToShower[ix]->progenitor()->isFinalState()) continue;
if(!particlesToShower[ix]->progenitor()->dataPtr()->coloured()) continue;
Energy ipt;
if(UseRandom::rnd() > _beta) {
ipt=_iptrms*sqrt(-log(UseRandom::rnd()));
}
else {
ipt=_gamma*sqrt(pow(1.+sqr(_iptmax/_gamma), UseRandom::rnd())-1.);
}
pair<Energy,double> pt = make_pair(ipt,UseRandom::rnd(Constants::twopi));
_intrinsic[particlesToShower[ix]] = pt;
}
}
void Evolver::setupMaximumScales(vector<ShowerProgenitorPtr> & p) {
// return if no vetos
if (!hardVetoOn()) return;
// find out if hard partonic subprocess.
bool isPartonic(false);
map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
cit = _currenttree->incomingLines().begin();
Lorentz5Momentum pcm;
for(; cit!=currentTree()->incomingLines().end(); ++cit) {
pcm += cit->first->progenitor()->momentum();
isPartonic |= cit->first->progenitor()->coloured();
}
// find maximum pt from hard process, the maximum pt from all outgoing
// coloured lines (this is simpler and more general than
// 2stu/(s^2+t^2+u^2)). Maximum scale for scattering processes will
// be transverse mass.
Energy ptmax = -1.0*GeV;
// general case calculate the scale
if (!hardVetoXComb()||
(hardVetoReadOption()&&
!ShowerHandler::currentHandler()->firstInteraction())) {
// scattering process
if(currentTree()->isHard()) {
// coloured incoming particles
if (isPartonic) {
map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
cjt = currentTree()->outgoingLines().begin();
for(; cjt!=currentTree()->outgoingLines().end(); ++cjt) {
if (cjt->first->progenitor()->coloured())
ptmax = max(ptmax,cjt->first->progenitor()->momentum().mt());
}
}
if (ptmax < ZERO) ptmax = pcm.m();
if(hardVetoXComb()&&hardVetoReadOption()&&
!ShowerHandler::currentHandler()->firstInteraction()) {
ptmax=min(ptmax,sqrt(ShowerHandler::currentHandler()
->lastXCombPtr()->lastScale()));
}
}
// decay, incoming() is the decaying particle.
else {
ptmax = currentTree()->incomingLines().begin()->first
->progenitor()->momentum().mass();
}
}
// hepeup.SCALUP is written into the lastXComb by the
// LesHouchesReader itself - use this by user's choice.
// Can be more general than this.
else {
ptmax = sqrt( ShowerHandler::currentHandler()
->lastXCombPtr()->lastScale() );
}
// set maxHardPt for all progenitors. For partonic processes this
// is now the max pt in the FS, for non-partonic processes or
// processes with no coloured FS the invariant mass of the IS
vector<ShowerProgenitorPtr>::const_iterator ckt = p.begin();
for (; ckt != p.end(); ckt++) (*ckt)->maxHardPt(ptmax);
}
void Evolver::showerHardProcess(ShowerTreePtr hard, XCPtr xcomb) {
_hardme = HwMEBasePtr();
// extract the matrix element
tStdXCombPtr lastXC = dynamic_ptr_cast<tStdXCombPtr>(xcomb);
if(lastXC) {
_hardme = dynamic_ptr_cast<HwMEBasePtr>(lastXC->matrixElement());
}
_decayme = HwDecayerBasePtr();
// set the current tree
currentTree(hard);
hardTree(HardTreePtr());
// number of attempts if more than one interaction switched on
unsigned int interactionTry=0;
do {
try {
// generate the showering
doShowering(true);
// if no vetos return
return;
}
catch (InteractionVeto) {
currentTree()->clear();
++interactionTry;
}
}
while(interactionTry<=5);
throw Exception() << "Too many tries for shower in "
<< "Evolver::showerHardProcess()"
<< Exception::eventerror;
}
void Evolver::hardMatrixElementCorrection(bool hard) {
// set the initial enhancement factors for the soft correction
_initialenhance = 1.;
_finalenhance = 1.;
// if hard matrix element switched off return
if(!MECOn()) return;
// see if we can get the correction from the matrix element
// or decayer
if(hard) {
if(_hardme&&_hardme->hasMECorrection()) {
_hardme->initializeMECorrection(_currenttree,
_initialenhance,_finalenhance);
if(hardMEC())
_hardme->applyHardMatrixElementCorrection(_currenttree);
}
}
else {
if(_decayme&&_decayme->hasMECorrection()) {
_decayme->initializeMECorrection(_currenttree,
_initialenhance,_finalenhance);
if(hardMEC())
_decayme->applyHardMatrixElementCorrection(_currenttree);
}
}
}
bool Evolver::timeLikeShower(tShowerParticlePtr particle,
ShowerInteraction::Type type,
bool first) {
// don't do anything if not needed
if(_limitEmissions == 1 || hardOnly() ||
( _limitEmissions == 2 && _nfs != 0) ||
( _limitEmissions == 4 && _nfs + _nis != 0) ) return false;
ShowerParticleVector theChildren;
int ntry=0;
do {
++ntry;
// generate the emission
Branching fb;
while (true) {
fb=_splittingGenerator->chooseForwardBranching(*particle,_finalenhance,type);
// no emission return
if(!fb.kinematics) return false;
// if emission OK break
if(!timeLikeVetoed(fb,particle,type)) break;
// otherwise reset scale and continue - SO IS involved in veto algorithm
particle->vetoEmission(fb.type,fb.kinematics->scale());
}
// has emitted
// Assign the shower kinematics to the emitting particle.
particle->showerKinematics(fb.kinematics);
// Assign the splitting function to the emitting particle.
// For the time being we are considering only 1->2 branching
// Create the ShowerParticle objects for the two children of
// the emitting particle; set the parent/child relationship
// if same as definition create particles, otherwise create cc
tcPDPtr pdata[2];
for(unsigned int ix=0;ix<2;++ix) pdata[ix]=getParticleData(fb.ids[ix+1]);
if(particle->id()!=fb.ids[0]) {
for(unsigned int ix=0;ix<2;++ix) {
tPDPtr cc(pdata[ix]->CC());
if(cc) pdata[ix]=cc;
}
}
theChildren.push_back(new_ptr(ShowerParticle(pdata[0],true)));
theChildren.push_back(new_ptr(ShowerParticle(pdata[1],true)));
// update the children
particle->showerKinematics()->updateChildren(particle, theChildren,
fb.type,true);
// update number of emissions
++_nfs;
if(_limitEmissions!=0) return true;
// shower the first particle
timeLikeShower(theChildren[0],type,false);
// shower the second particle
timeLikeShower(theChildren[1],type,false);
// that's if for old approach
if(_reconOpt==0) break;
// branching has happened
particle->showerKinematics()->updateParent(particle, theChildren,
fb.type,true);
// clean up the vetoed emission
if(particle->virtualMass()==ZERO) {
particle->showerKinematics(ShoKinPtr());
for(unsigned int ix=0;ix<theChildren.size();++ix)
particle->abandonChild(theChildren[ix]);
theChildren.clear();
}
}
while(particle->virtualMass()==ZERO&&ntry<50);
if(first)
particle->showerKinematics()->resetChildren(particle,theChildren);
return true;
}
bool
Evolver::spaceLikeShower(tShowerParticlePtr particle, PPtr beam,
ShowerInteraction::Type type) {
//using the pdf's associated with the ShowerHandler assures, that
//modified pdf's are used for the secondary interactions via
//CascadeHandler::resetPDFs(...)
tcPDFPtr pdf;
if(ShowerHandler::currentHandler()->firstPDF().particle() == _beam)
pdf = ShowerHandler::currentHandler()->firstPDF().pdf();
if(ShowerHandler::currentHandler()->secondPDF().particle() == _beam)
pdf = ShowerHandler::currentHandler()->secondPDF().pdf();
Energy freeze = ShowerHandler::currentHandler()->pdfFreezingScale();
// don't do anything if not needed
if(_limitEmissions == 2 || hardOnly() ||
( _limitEmissions == 1 && _nis != 0 ) ||
( _limitEmissions == 4 && _nis + _nfs != 0 ) ) return false;
Branching bb;
// generate branching
while (true) {
bb=_splittingGenerator->chooseBackwardBranching(*particle,beam,
_initialenhance,
_beam,type,
pdf,freeze);
// return if no emission
if(!bb.kinematics) return false;
// if not vetoed break
if(!spaceLikeVetoed(bb,particle,type)) break;
// otherwise reset scale and continue
particle->vetoEmission(bb.type,bb.kinematics->scale());
}
// assign the splitting function and shower kinematics
particle->showerKinematics(bb.kinematics);
// For the time being we are considering only 1->2 branching
// particles as in Sudakov form factor
tcPDPtr part[2]={getParticleData(bb.ids[0]),
getParticleData(bb.ids[2])};
if(particle->id()!=bb.ids[1]) {
if(part[0]->CC()) part[0]=part[0]->CC();
if(part[1]->CC()) part[1]=part[1]->CC();
}
// Now create the actual particles, make the otherChild a final state
// particle, while the newParent is not
ShowerParticlePtr newParent=new_ptr(ShowerParticle(part[0],false));
ShowerParticlePtr otherChild = new_ptr(ShowerParticle(part[1],true,true));
ShowerParticleVector theChildren;
theChildren.push_back(particle);
theChildren.push_back(otherChild);
//this updates the evolution scale
particle->showerKinematics()->updateParent(newParent, theChildren,
bb.type,true);
// update the history if needed
_currenttree->updateInitialStateShowerProduct(_progenitor,newParent);
_currenttree->addInitialStateBranching(particle,newParent,otherChild);
// for the reconstruction of kinematics, parent/child
// relationships are according to the branching process:
// now continue the shower
++_nis;
bool emitted = _limitEmissions==0 ?
spaceLikeShower(newParent,beam,type) : false;
// now reconstruct the momentum
if(!emitted) {
if(_intrinsic.find(_progenitor)==_intrinsic.end()) {
bb.kinematics->updateLast(newParent,ZERO,ZERO);
}
else {
pair<Energy,double> kt=_intrinsic[_progenitor];
bb.kinematics->updateLast(newParent,
kt.first*cos(kt.second),
kt.first*sin(kt.second));
}
}
particle->showerKinematics()->updateChildren(newParent, theChildren,
bb.type,true);
if(_limitEmissions!=0) return true;
// perform the shower of the final-state particle
timeLikeShower(otherChild,type,true);
// return the emitted
return true;
}
void Evolver::showerDecay(ShowerTreePtr decay) {
_decayme = HwDecayerBasePtr();
_hardme = HwMEBasePtr();
// find the decayer
// try the normal way if possible
tDMPtr dm = decay->incomingLines().begin()->first->original() ->decayMode();
if(!dm) dm = decay->incomingLines().begin()->first->copy() ->decayMode();
if(!dm) dm = decay->incomingLines().begin()->first->progenitor()->decayMode();
// otherwise make a string and look it up
if(!dm) {
string tag = decay->incomingLines().begin()->first->original()->dataPtr()->name()
+ "->";
for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
it=decay->outgoingLines().begin();it!=decay->outgoingLines().end();++it) {
if(it!=decay->outgoingLines().begin()) tag += ",";
tag += it->first->original()->dataPtr()->name();
}
tag += ";";
dm = generator()->findDecayMode(tag);
}
if(dm) _decayme = dynamic_ptr_cast<HwDecayerBasePtr>(dm->decayer());
// set the ShowerTree to be showered
currentTree(decay);
decay->applyTransforms();
hardTree(HardTreePtr());
unsigned int interactionTry=0;
do {
try {
// generate the showering
doShowering(false);
// if no vetos return
return;
}
catch (InteractionVeto) {
currentTree()->clear();
++interactionTry;
}
}
while(interactionTry<=5);
throw Exception() << "Too many tries for QED shower in Evolver::showerDecay()"
<< Exception::eventerror;
}
bool Evolver::spaceLikeDecayShower(tShowerParticlePtr particle,
const map<ShowerPartnerType::Type,pair<Energy,Energy> > & maxScales,
Energy minmass,ShowerInteraction::Type type) {
Branching fb;
while (true) {
fb=_splittingGenerator->chooseDecayBranching(*particle,maxScales,minmass,
_initialenhance,type);
// return if no radiation
if(!fb.kinematics) return false;
// if not vetoed break
if(!spaceLikeDecayVetoed(fb,particle,type)) break;
// otherwise reset scale and continue
particle->vetoEmission(fb.type,fb.kinematics->scale());
}
// has emitted
// Assign the shower kinematics to the emitting particle.
particle->showerKinematics(fb.kinematics);
// For the time being we are considering only 1->2 branching
// Create the ShowerParticle objects for the two children of
// the emitting particle; set the parent/child relationship
// if same as definition create particles, otherwise create cc
tcPDPtr pdata[2];
for(unsigned int ix=0;ix<2;++ix) pdata[ix]=getParticleData(fb.ids[ix+1]);
if(particle->id()!=fb.ids[0]) {
for(unsigned int ix=0;ix<2;++ix) {
tPDPtr cc(pdata[ix]->CC());
if(cc) pdata[ix]=cc;
}
}
ShowerParticleVector theChildren;
theChildren.push_back(new_ptr(ShowerParticle(pdata[0],true)));
theChildren.push_back(new_ptr(ShowerParticle(pdata[1],true)));
// some code moved to updateChildren
particle->showerKinematics()->updateChildren(particle, theChildren,
fb.type,true);
// In the case of splittings which involves coloured particles,
// set properly the colour flow of the branching.
// update the history if needed
_currenttree->updateInitialStateShowerProduct(_progenitor,theChildren[0]);
_currenttree->addInitialStateBranching(particle,theChildren[0],theChildren[1]);
// shower the first particle
spaceLikeDecayShower(theChildren[0],maxScales,minmass,type);
// shower the second particle
timeLikeShower(theChildren[1],type,true);
// branching has happened
return true;
}
vector<ShowerProgenitorPtr> Evolver::setupShower(bool hard) {
// generate POWHEG hard emission if needed
if(_hardEmissionMode==1) hardestEmission(hard);
// set the initial colour partners
setEvolutionPartners(hard,_hardtree ?
_hardtree->interaction() : interactions_[0]);
// generate hard me if needed
if(_hardEmissionMode==0) hardMatrixElementCorrection(hard);
// get the particles to be showered
vector<ShowerProgenitorPtr> particlesToShower =
currentTree()->extractProgenitors();
// remake the colour partners if needed
if(_hardEmissionMode==0 && _currenttree->hardMatrixElementCorrection()) {
setEvolutionPartners(hard,interactions_[0]);
_currenttree->resetShowerProducts();
}
// return the answer
return particlesToShower;
}
void Evolver::setEvolutionPartners(bool hard,ShowerInteraction::Type type) {
// match the particles in the ShowerTree and hardTree
if(hardTree() && !hardTree()->connect(currentTree()))
throw Exception() << "Can't match trees in "
<< "Evolver::setEvolutionPartners()"
<< Exception::eventerror;
// extract the progenitors
vector<ShowerParticlePtr> particles =
currentTree()->extractProgenitorParticles();
// sort out the colour partners
if(hardTree()) {
// find the partner
for(unsigned int ix=0;ix<particles.size();++ix) {
tHardBranchingPtr partner =
hardTree()->particles()[particles[ix]]->colourPartner();
if(!partner) continue;
for(map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
it=hardTree()->particles().begin();
it!=hardTree()->particles().end();++it) {
if(it->second==partner) particles[ix]->partner(it->first);
}
if(!particles[ix]->partner())
throw Exception() << "Can't match partners in "
<< "Evolver::setEvolutionPartners()"
<< Exception::eventerror;
}
}
// Set the initial evolution scales
showerModel()->partnerFinder()->
setInitialEvolutionScales(particles,!hard,type,!_hardtree);
}
void Evolver::updateHistory(tShowerParticlePtr particle) {
if(!particle->children().empty()) {
ShowerParticleVector theChildren;
for(unsigned int ix=0;ix<particle->children().size();++ix) {
ShowerParticlePtr part = dynamic_ptr_cast<ShowerParticlePtr>
(particle->children()[ix]);
theChildren.push_back(part);
}
// update the history if needed
if(particle==_currenttree->getFinalStateShowerProduct(_progenitor))
_currenttree->updateFinalStateShowerProduct(_progenitor,
particle,theChildren);
_currenttree->addFinalStateBranching(particle,theChildren);
for(unsigned int ix=0;ix<theChildren.size();++ix)
updateHistory(theChildren[ix]);
}
}
bool Evolver::startTimeLikeShower(ShowerInteraction::Type type) {
if(hardTree()) {
map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
eit=hardTree()->particles().end(),
mit = hardTree()->particles().find(progenitor()->progenitor());
if( mit != eit && !mit->second->children().empty() ) {
bool output=truncatedTimeLikeShower(progenitor()->progenitor(),
mit->second ,type);
if(output) updateHistory(progenitor()->progenitor());
return output;
}
}
bool output = hardOnly() ? false :
timeLikeShower(progenitor()->progenitor() ,type,true) ;
if(output) updateHistory(progenitor()->progenitor());
return output;
}
bool Evolver::startSpaceLikeShower(PPtr parent, ShowerInteraction::Type type) {
if(hardTree()) {
map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
eit =hardTree()->particles().end(),
mit = hardTree()->particles().find(progenitor()->progenitor());
if( mit != eit && mit->second->parent() ) {
return truncatedSpaceLikeShower( progenitor()->progenitor(),
parent, mit->second->parent(), type );
}
}
return hardOnly() ? false :
spaceLikeShower(progenitor()->progenitor(),parent,type);
}
bool Evolver::
startSpaceLikeDecayShower(const map<ShowerPartnerType::Type,pair<Energy,Energy> > & maxScales,
Energy minimumMass,ShowerInteraction::Type type) {
if(hardTree()) {
map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
eit =hardTree()->particles().end(),
mit = hardTree()->particles().find(progenitor()->progenitor());
if( mit != eit && mit->second->parent() ) {
HardBranchingPtr branch=mit->second;
while(branch->parent()) branch=branch->parent();
return truncatedSpaceLikeDecayShower(progenitor()->progenitor(),maxScales,
minimumMass, branch ,type);
}
}
return hardOnly() ? false :
spaceLikeDecayShower(progenitor()->progenitor(),maxScales,minimumMass,type);
}
bool Evolver::timeLikeVetoed(const Branching & fb,
ShowerParticlePtr particle,
ShowerInteraction::Type type) {
// check whether emission was harder than largest pt of hard subprocess
if ( hardVetoFS() && fb.kinematics->pT() > _progenitor->maxHardPt() )
return true;
// soft matrix element correction veto
if( softMEC()) {
if(_hardme && _hardme->hasMECorrection()) {
if(_hardme->softMatrixElementVeto(_progenitor,particle,fb))
return true;
}
else if(_decayme && _decayme->hasMECorrection()) {
if(_decayme->softMatrixElementVeto(_progenitor,particle,fb))
return true;
}
}
// veto on maximum pt
if(fb.kinematics->pT()>_progenitor->maximumpT(type)) return true;
// general vetos
if (fb.kinematics && !_vetoes.empty()) {
bool vetoed=false;
for (vector<ShowerVetoPtr>::iterator v = _vetoes.begin();
v != _vetoes.end(); ++v) {
bool test = (**v).vetoTimeLike(_progenitor,particle,fb);
switch((**v).vetoType()) {
case ShowerVeto::Emission:
vetoed |= test;
break;
case ShowerVeto::Shower:
if(test) throw VetoShower();
break;
case ShowerVeto::Event:
if(test) throw Veto();
break;
}
}
if(vetoed) return true;
}
return false;
}
bool Evolver::spaceLikeVetoed(const Branching & bb,ShowerParticlePtr particle,
ShowerInteraction::Type type) {
// check whether emission was harder than largest pt of hard subprocess
if (hardVetoIS() && bb.kinematics->pT() > _progenitor->maxHardPt())
return true;
// apply the soft correction
if( softMEC() && _hardme && _hardme->hasMECorrection() ) {
if(_hardme->softMatrixElementVeto(_progenitor,particle,bb))
return true;
}
// the more general vetos
// check vs max pt for the shower
if(bb.kinematics->pT()>_progenitor->maximumpT(type)) return true;
if (!_vetoes.empty()) {
bool vetoed=false;
for (vector<ShowerVetoPtr>::iterator v = _vetoes.begin();
v != _vetoes.end(); ++v) {
bool test = (**v).vetoSpaceLike(_progenitor,particle,bb);
switch ((**v).vetoType()) {
case ShowerVeto::Emission:
vetoed |= test;
break;
case ShowerVeto::Shower:
if(test) throw VetoShower();
break;
case ShowerVeto::Event:
if(test) throw Veto();
break;
}
}
if (vetoed) return true;
}
return false;
}
bool Evolver::spaceLikeDecayVetoed( const Branching & fb,
ShowerParticlePtr particle,
ShowerInteraction::Type type ) {
// apply the soft correction
if( softMEC() && _decayme && _decayme->hasMECorrection() ) {
if(_decayme->softMatrixElementVeto(_progenitor,particle,fb))
return true;
}
// veto on hardest pt in the shower
if(fb.kinematics->pT()> _progenitor->maximumpT(type)) return true;
// general vetos
if (!_vetoes.empty()) {
bool vetoed=false;
for (vector<ShowerVetoPtr>::iterator v = _vetoes.begin();
v != _vetoes.end(); ++v) {
bool test = (**v).vetoSpaceLike(_progenitor,particle,fb);
switch((**v).vetoType()) {
case ShowerVeto::Emission:
vetoed |= test;
break;
case ShowerVeto::Shower:
if(test) throw VetoShower();
break;
case ShowerVeto::Event:
if(test) throw Veto();
break;
}
if (vetoed) return true;
}
}
return false;
}
void Evolver::hardestEmission(bool hard) {
if( ( _hardme && _hardme->hasPOWHEGCorrection()) ||
(_decayme && _decayme->hasPOWHEGCorrection())) {
if(_hardme) {
assert(hard);
_hardtree = _hardme->generateHardest( currentTree(),interactions_ );
}
else {
assert(!hard);
_hardtree = _decayme->generateHardest( currentTree() );
}
if(!_hardtree) return;
// join up the two trees
connectTrees(currentTree(),_hardtree,hard);
}
else {
_hardtree = ShowerHandler::currentHandler()->generateCKKW(currentTree());
}
}
bool Evolver::truncatedTimeLikeShower(tShowerParticlePtr particle,
HardBranchingPtr branch,
ShowerInteraction::Type type) {
int ntry=0;
do {
++ntry;
Branching fb;
unsigned int iout=0;
tcPDPtr pdata[2];
while (true) {
// no truncated shower break
if(!isTruncatedShowerON()||hardOnly()) break;
// generate emission
fb=splittingGenerator()->chooseForwardBranching(*particle,1.,type);
// no emission break
if(!fb.kinematics) break;
// check haven't evolved too far
if(fb.kinematics->scale() < branch->scale()) {
fb=Branching();
break;
}
// get the particle data objects
for(unsigned int ix=0;ix<2;++ix) pdata[ix]=getParticleData(fb.ids[ix+1]);
if(particle->id()!=fb.ids[0]) {
for(unsigned int ix=0;ix<2;++ix) {
tPDPtr cc(pdata[ix]->CC());
if(cc) pdata[ix]=cc;
}
}
// find the truncated line
iout=0;
if(pdata[0]->id()!=pdata[1]->id()) {
if(pdata[0]->id()==particle->id()) iout=1;
else if (pdata[1]->id()==particle->id()) iout=2;
}
else if(pdata[0]->id()==particle->id()) {
if(fb.kinematics->z()>0.5) iout=1;
else iout=2;
}
// apply the vetos for the truncated shower
// no flavour changing branchings
if(iout==0) {
particle->vetoEmission(fb.type,fb.kinematics->scale());
continue;
}
double zsplit = iout==1 ? fb.kinematics->z() : 1-fb.kinematics->z();
// only if same interaction for forced branching
// and evolution
if(type==branch->sudakov()->interactionType()) {
if(zsplit < 0.5 || // hardest line veto
fb.kinematics->scale()*zsplit < branch->scale() ) { // angular ordering veto
particle->vetoEmission(fb.type,fb.kinematics->scale());
continue;
}
}
// pt veto
if(fb.kinematics->pT() > progenitor()->maximumpT(type)) {
particle->vetoEmission(fb.type,fb.kinematics->scale());
continue;
}
// should do base class vetos as well
if(timeLikeVetoed(fb,particle,type)) {
particle->vetoEmission(fb.type,fb.kinematics->scale());
continue;
}
break;
}
// if no branching force trunctaed emission
if(!fb.kinematics) {
// construct the kinematics for the hard emission
ShoKinPtr showerKin=
branch->sudakov()->createFinalStateBranching(branch->scale(),
branch->children()[0]->z(),
branch->phi(),
branch->children()[0]->pT());
- particle->evolutionScale(branch->type(),make_pair(branch->scale(),branch->scale()));
+ // particle->evolutionScale(branch->type(),make_pair(branch->scale(),branch->scale()));
+ assert(false);
showerKin->initialize( *particle,PPtr() );
IdList idlist(3);
idlist[0] = particle->id();
idlist[1] = branch->children()[0]->branchingParticle()->id();
idlist[2] = branch->children()[1]->branchingParticle()->id();
fb = Branching( showerKin, idlist, branch->sudakov(),branch->type() );
// Assign the shower kinematics to the emitting particle.
particle->showerKinematics( fb.kinematics );
// Assign the splitting function to the emitting particle.
// For the time being we are considering only 1->2 branching
// Create the ShowerParticle objects for the two children of
// the emitting particle; set the parent/child relationship
// if same as definition create particles, otherwise create cc
ShowerParticleVector theChildren;
theChildren.push_back(new_ptr(ShowerParticle(branch->children()[0]->
branchingParticle()->dataPtr(),true)));
theChildren.push_back(new_ptr(ShowerParticle(branch->children()[1]->
branchingParticle()->dataPtr(),true)));
particle->showerKinematics()->
updateChildren(particle, theChildren,fb.type,
type==branch->sudakov()->interactionType());
// shower the first particle
if( branch->children()[0]->children().empty() ) {
if( ! hardOnly() )
timeLikeShower(theChildren[0],type,false);
}
else {
truncatedTimeLikeShower( theChildren[0],branch->children()[0],type);
}
// shower the second particle
if( branch->children()[1]->children().empty() ) {
if( ! hardOnly() )
timeLikeShower( theChildren[1] , type,false);
}
else {
truncatedTimeLikeShower( theChildren[1],branch->children()[1] ,type);
}
// that's if for old approach
if(_reconOpt==0) return true;
// branching has happened
particle->showerKinematics()->updateParent(particle, theChildren,
fb.type,true);
// clean up the vetoed emission
if(particle->virtualMass()==ZERO) {
particle->showerKinematics(ShoKinPtr());
for(unsigned int ix=0;ix<theChildren.size();++ix)
particle->abandonChild(theChildren[ix]);
theChildren.clear();
}
else return true;
}
// has emitted
// Assign the shower kinematics to the emitting particle.
particle->showerKinematics(fb.kinematics);
// Assign the splitting function to the emitting particle.
// For the time being we are considering only 1->2 branching
// Create the ShowerParticle objects for the two children of
// the emitting particle; set the parent/child relationship
// if same as definition create particles, otherwise create cc
ShowerParticleVector theChildren;
theChildren.push_back( new_ptr( ShowerParticle( pdata[0], true ) ) );
theChildren.push_back( new_ptr( ShowerParticle( pdata[1], true ) ) );
particle->showerKinematics()->
updateChildren( particle, theChildren , fb.type, true);
// shower the first particle
if( iout == 1 ) truncatedTimeLikeShower( theChildren[0], branch , type );
else timeLikeShower( theChildren[0] , type,false);
// shower the second particle
if( iout == 2 ) truncatedTimeLikeShower( theChildren[1], branch , type );
else timeLikeShower( theChildren[1] , type,false);
// that's if for old approach
if(_reconOpt==0) return true;
// branching has happened
particle->showerKinematics()->updateParent(particle, theChildren,fb.type,true);
// clean up the vetoed emission
if(particle->virtualMass()==ZERO) {
particle->showerKinematics(ShoKinPtr());
for(unsigned int ix=0;ix<theChildren.size();++ix)
particle->abandonChild(theChildren[ix]);
theChildren.clear();
}
else return true;
}
while(ntry<50);
return false;
}
bool Evolver::truncatedSpaceLikeShower(tShowerParticlePtr particle, PPtr beam,
HardBranchingPtr branch,
ShowerInteraction::Type type) {
tcPDFPtr pdf;
if(ShowerHandler::currentHandler()->firstPDF().particle() == beamParticle())
pdf = ShowerHandler::currentHandler()->firstPDF().pdf();
if(ShowerHandler::currentHandler()->secondPDF().particle() == beamParticle())
pdf = ShowerHandler::currentHandler()->secondPDF().pdf();
Energy freeze = ShowerHandler::currentHandler()->pdfFreezingScale();
Branching bb;
// parameters of the force branching
double z(0.);
HardBranchingPtr timelike;
for( unsigned int ix = 0; ix < branch->children().size(); ++ix ) {
if( branch->children()[ix]->status() ==HardBranching::Outgoing) {
timelike = branch->children()[ix];
}
if( branch->children()[ix]->status() ==HardBranching::Incoming )
z = branch->children()[ix]->z();
}
// generate truncated branching
tcPDPtr part[2];
if(z>=0.&&z<=1.) {
while (true) {
if( !isTruncatedShowerON() || hardOnly() ) break;
bb = splittingGenerator()->chooseBackwardBranching( *particle,
beam, 1., beamParticle(),
type , pdf,freeze);
if( !bb.kinematics || bb.kinematics->scale() < branch->scale() ) {
bb = Branching();
break;
}
// particles as in Sudakov form factor
part[0] = getParticleData( bb.ids[0] );
part[1] = getParticleData( bb.ids[2] );
//is emitter anti-particle
if( particle->id() != bb.ids[1]) {
if( part[0]->CC() ) part[0] = part[0]->CC();
if( part[1]->CC() ) part[1] = part[1]->CC();
}
double zsplit = bb.kinematics->z();
// apply the vetos for the truncated shower
// if doesn't carry most of momentum
if(type==branch->sudakov()->interactionType() &&
zsplit < 0.5) {
particle->vetoEmission(bb.type,bb.kinematics->scale());
continue;
}
// others
if( part[0]->id() != particle->id() || // if particle changes type
bb.kinematics->pT() > progenitor()->maximumpT(type) || // pt veto
bb.kinematics->scale() < branch->scale()) { // angular ordering veto
particle->vetoEmission(bb.type,bb.kinematics->scale());
continue;
}
// and those from the base class
if(spaceLikeVetoed(bb,particle,type)) {
particle->vetoEmission(bb.type,bb.kinematics->scale());
continue;
}
break;
}
}
if( !bb.kinematics ) {
//do the hard emission
ShoKinPtr kinematics =
branch->sudakov()->createInitialStateBranching( branch->scale(), z, branch->phi(),
branch->children()[0]->pT() );
kinematics->initialize( *particle, beam );
// assign the splitting function and shower kinematics
particle->showerKinematics( kinematics );
// For the time being we are considering only 1->2 branching
// Now create the actual particles, make the otherChild a final state
// particle, while the newParent is not
ShowerParticlePtr newParent =
new_ptr( ShowerParticle( branch->branchingParticle()->dataPtr(), false ) );
ShowerParticlePtr otherChild =
new_ptr( ShowerParticle( timelike->branchingParticle()->dataPtr(),
true, true ) );
ShowerParticleVector theChildren;
theChildren.push_back( particle );
theChildren.push_back( otherChild );
particle->showerKinematics()->
updateParent( newParent, theChildren, branch->type(),
type==branch->sudakov()->interactionType() );
// update the history if needed
currentTree()->updateInitialStateShowerProduct( progenitor(), newParent );
currentTree()->addInitialStateBranching( particle, newParent, otherChild );
// for the reconstruction of kinematics, parent/child
// relationships are according to the branching process:
// now continue the shower
bool emitted=false;
if(!hardOnly()) {
if( branch->parent() ) {
emitted = truncatedSpaceLikeShower( newParent, beam, branch->parent() , type);
}
else {
emitted = spaceLikeShower( newParent, beam , type);
}
}
if( !emitted ) {
if( intrinsicpT().find( progenitor() ) == intrinsicpT().end() ) {
kinematics->updateLast( newParent, ZERO, ZERO );
}
else {
pair<Energy,double> kt = intrinsicpT()[progenitor()];
kinematics->updateLast( newParent,
kt.first*cos( kt.second ),
kt.first*sin( kt.second ) );
}
}
particle->showerKinematics()->
updateChildren( newParent, theChildren,bb.type,
type==branch->sudakov()->interactionType() );
if(hardOnly()) return true;
// perform the shower of the final-state particle
if( timelike->children().empty() ) {
timeLikeShower( otherChild , type,true);
}
else {
truncatedTimeLikeShower( otherChild, timelike , type);
}
// return the emitted
return true;
}
// assign the splitting function and shower kinematics
particle->showerKinematics( bb.kinematics );
// For the time being we are considering only 1->2 branching
// Now create the actual particles, make the otherChild a final state
// particle, while the newParent is not
ShowerParticlePtr newParent = new_ptr( ShowerParticle( part[0], false ) );
ShowerParticlePtr otherChild = new_ptr( ShowerParticle( part[1], true, true ) );
ShowerParticleVector theChildren;
theChildren.push_back( particle );
theChildren.push_back( otherChild );
particle->showerKinematics()->updateParent( newParent, theChildren ,
bb.type, true);
// update the history if needed
currentTree()->updateInitialStateShowerProduct( progenitor(), newParent );
currentTree()->addInitialStateBranching( particle, newParent, otherChild );
// for the reconstruction of kinematics, parent/child
// relationships are according to the branching process:
// now continue the shower
bool emitted = truncatedSpaceLikeShower( newParent, beam, branch,type);
// now reconstruct the momentum
if( !emitted ) {
if( intrinsicpT().find( progenitor() ) == intrinsicpT().end() ) {
bb.kinematics->updateLast( newParent, ZERO, ZERO );
}
else {
pair<Energy,double> kt = intrinsicpT()[ progenitor() ];
bb.kinematics->updateLast( newParent,
kt.first*cos( kt.second ),
kt.first*sin( kt.second ) );
}
}
particle->showerKinematics()->updateChildren( newParent, theChildren ,
bb.type, true);
// perform the shower of the final-state particle
timeLikeShower( otherChild , type,true);
// return the emitted
return true;
}
bool Evolver::
truncatedSpaceLikeDecayShower(tShowerParticlePtr particle,
const map<ShowerPartnerType::Type,pair<Energy,Energy> > & maxScales,
Energy minmass, HardBranchingPtr branch,
ShowerInteraction::Type type) {
Branching fb;
unsigned int iout=0;
tcPDPtr pdata[2];
while (true) {
// no truncated shower break
if(!isTruncatedShowerON()||hardOnly()) break;
fb=splittingGenerator()->chooseDecayBranching(*particle,maxScales,minmass,1.,type);
// return if no radiation
if(!fb.kinematics) break;
// check haven't evolved too far
if(fb.kinematics->scale() < branch->scale()) {
fb=Branching();
break;
}
// get the particle data objects
for(unsigned int ix=0;ix<2;++ix) pdata[ix]=getParticleData(fb.ids[ix+1]);
if(particle->id()!=fb.ids[0]) {
for(unsigned int ix=0;ix<2;++ix) {
tPDPtr cc(pdata[ix]->CC());
if(cc) pdata[ix]=cc;
}
}
// find the truncated line
iout=0;
if(pdata[0]->id()!=pdata[1]->id()) {
if(pdata[0]->id()==particle->id()) iout=1;
else if (pdata[1]->id()==particle->id()) iout=2;
}
else if(pdata[0]->id()==particle->id()) {
if(fb.kinematics->z()>0.5) iout=1;
else iout=2;
}
// apply the vetos for the truncated shower
// no flavour changing branchings
if(iout==0) {
particle->vetoEmission(fb.type,fb.kinematics->scale());
continue;
}
double zsplit = iout==1 ? fb.kinematics->z() : 1-fb.kinematics->z();
if(type==branch->sudakov()->interactionType()) {
if(zsplit < 0.5 || // hardest line veto
fb.kinematics->scale()*zsplit < branch->scale() ) { // angular ordering veto
particle->vetoEmission(fb.type,fb.kinematics->scale());
continue;
}
}
// pt veto
if(fb.kinematics->pT() > progenitor()->maximumpT(type)) {
particle->vetoEmission(fb.type,fb.kinematics->scale());
continue;
}
// should do base class vetos as well
// if not vetoed break
if(!spaceLikeDecayVetoed(fb,particle,type)) break;
// otherwise reset scale and continue
particle->vetoEmission(fb.type,fb.kinematics->scale());
}
// if no branching insert hard emission and continue
if(!fb.kinematics) {
// construct the kinematics for the hard emission
ShoKinPtr showerKin=
branch->sudakov()->createDecayBranching(branch->scale(),
branch->children()[0]->z(),
branch->phi(),
branch->children()[0]->pT());
- particle->evolutionScale(branch->type(),make_pair(branch->scale(),branch->scale()));
+ // particle->evolutionScale(branch->type(),make_pair(branch->scale(),branch->scale()));
+ assert(false);
showerKin->initialize( *particle,PPtr() );
IdList idlist(3);
idlist[0] = particle->id();
idlist[1] = branch->children()[0]->branchingParticle()->id();
idlist[2] = branch->children()[1]->branchingParticle()->id();
assert(false);
fb = Branching( showerKin, idlist, branch->sudakov(),ShowerPartnerType::QCDColourLine );
// Assign the shower kinematics to the emitting particle.
particle->showerKinematics( fb.kinematics );
// Assign the splitting function to the emitting particle.
// For the time being we are considering only 1->2 branching
// Create the ShowerParticle objects for the two children of
// the emitting particle; set the parent/child relationship
// if same as definition create particles, otherwise create cc
ShowerParticleVector theChildren;
theChildren.push_back(new_ptr(ShowerParticle(branch->children()[0]->
branchingParticle()->dataPtr(),true)));
theChildren.push_back(new_ptr(ShowerParticle(branch->children()[1]->
branchingParticle()->dataPtr(),true)));
particle->showerKinematics()->
updateChildren(particle, theChildren,fb.type,
type==branch->sudakov()->interactionType());
if(theChildren[0]->id()==particle->id()) {
// update the history if needed
currentTree()->updateInitialStateShowerProduct(progenitor(),theChildren[0]);
currentTree()->addInitialStateBranching(particle,theChildren[0],theChildren[1]);
// shower the space-like particle
if( branch->children()[0]->children().empty() ) {
if( ! hardOnly() ) spaceLikeDecayShower(theChildren[0],maxScales,minmass,type);
}
else {
truncatedSpaceLikeDecayShower( theChildren[0],maxScales,minmass,
branch->children()[0],type);
}
// shower the second particle
if( branch->children()[1]->children().empty() ) {
if( ! hardOnly() ) timeLikeShower( theChildren[1] , type,true);
}
else {
truncatedTimeLikeShower( theChildren[1],branch->children()[1] ,type);
}
}
else {
// update the history if needed
currentTree()->updateInitialStateShowerProduct(progenitor(),theChildren[1]);
currentTree()->addInitialStateBranching(particle,theChildren[0],theChildren[1]);
// shower the space-like particle
if( branch->children()[1]->children().empty() ) {
if( ! hardOnly() ) spaceLikeDecayShower(theChildren[1],maxScales,minmass,type);
}
else {
truncatedSpaceLikeDecayShower( theChildren[1],maxScales,minmass,
branch->children()[1],type);
}
// shower the second particle
if( branch->children()[0]->children().empty() ) {
if( ! hardOnly() ) timeLikeShower( theChildren[0] , type,true);
}
else {
truncatedTimeLikeShower( theChildren[0],branch->children()[0] ,type);
}
}
return true;
}
// has emitted
// Assign the shower kinematics to the emitting particle.
particle->showerKinematics(fb.kinematics);
// For the time being we are considering only 1->2 branching
// Create the ShowerParticle objects for the two children of
// the emitting particle; set the parent/child relationship
// if same as definition create particles, otherwise create cc
ShowerParticleVector theChildren;
theChildren.push_back(new_ptr(ShowerParticle(pdata[0],true)));
theChildren.push_back(new_ptr(ShowerParticle(pdata[1],true)));
particle->showerKinematics()->updateChildren(particle, theChildren,fb.type,true);
// In the case of splittings which involves coloured particles,
// set properly the colour flow of the branching.
// update the history if needed
currentTree()->updateInitialStateShowerProduct(progenitor(),theChildren[0]);
currentTree()->addInitialStateBranching(particle,theChildren[0],theChildren[1]);
// shower the first particle
truncatedSpaceLikeDecayShower(theChildren[0],maxScales,minmass,branch,type);
// shower the second particle
timeLikeShower(theChildren[1],type,true);
// branching has happened
return true;
}
bool Evolver::constructDecayTree(vector<ShowerProgenitorPtr> & particlesToShower,
ShowerInteraction::Type inter) {
Energy ptmax(-GeV);
// get the maximum pt is all ready a hard tree
if(hardTree()) {
for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
if(particlesToShower[ix]->maximumpT(inter)>ptmax&&
particlesToShower[ix]->progenitor()->isFinalState())
ptmax = particlesToShower[ix]->maximumpT(inter);
}
}
vector<HardBranchingPtr> spaceBranchings,allBranchings;
for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
if(particlesToShower[ix]->progenitor()->isFinalState()) {
HardBranchingPtr newBranch;
if(particlesToShower[ix]->hasEmitted()) {
newBranch =
new_ptr(HardBranching(particlesToShower[ix]->progenitor(),
particlesToShower[ix]->progenitor()->
showerKinematics()->SudakovFormFactor(),
HardBranchingPtr(),HardBranching::Outgoing));
constructTimeLikeLine(newBranch,particlesToShower[ix]->progenitor());
}
else {
newBranch =
new_ptr(HardBranching(particlesToShower[ix]->progenitor(),
SudakovPtr(),HardBranchingPtr(),
HardBranching::Outgoing));
}
allBranchings.push_back(newBranch);
}
else {
HardBranchingPtr newBranch;
if(particlesToShower[ix]->hasEmitted()) {
newBranch =
new_ptr(HardBranching(particlesToShower[ix]->progenitor(),
particlesToShower[ix]->progenitor()->
showerKinematics()->SudakovFormFactor(),
HardBranchingPtr(),HardBranching::Decay));
constructTimeLikeLine(newBranch,particlesToShower[ix]->progenitor());
HardBranchingPtr last=newBranch;
do {
for(unsigned int ix=0;ix<last->children().size();++ix) {
if(last->children()[ix]->branchingParticle()->id()==
particlesToShower[ix]->id()) {
last = last->children()[ix];
continue;
}
}
}
while(!last->children().empty());
last->status(HardBranching::Incoming);
spaceBranchings.push_back(newBranch);
allBranchings .push_back(last);
}
else {
newBranch =
new_ptr(HardBranching(particlesToShower[ix]->progenitor(),
SudakovPtr(),HardBranchingPtr(),
HardBranching::Incoming));
spaceBranchings.push_back(newBranch);
allBranchings .push_back(newBranch);
}
}
}
HardTreePtr QCDTree = new_ptr(HardTree(allBranchings,spaceBranchings,inter));
// set the charge partners
ShowerParticleVector particles;
particles.push_back(spaceBranchings.back()->branchingParticle());
for(set<HardBranchingPtr>::iterator cit=QCDTree->branchings().begin();
cit!=QCDTree->branchings().end();++cit) {
if((*cit)->status()==HardBranching::Outgoing)
particles.push_back((*cit)->branchingParticle());
}
// get the partners
showerModel()->partnerFinder()->setInitialEvolutionScales(particles,true,inter,true);
// do the inverse recon
if(!showerModel()->kinematicsReconstructor()->
deconstructDecayJets(QCDTree,this,inter)) {
return false;
}
// clear the old shower
currentTree()->clear();
// set the hard tree
hardTree(QCDTree);
// set the charge partners
setEvolutionPartners(false,inter);
// get the particles to be showered
map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator cit;
map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator cjt;
particlesToShower.clear();
// incoming particles
for(cit=currentTree()->incomingLines().begin();
cit!=currentTree()->incomingLines().end();++cit)
particlesToShower.push_back(((*cit).first));
assert(particlesToShower.size()==1);
// outgoing particles
for(cjt=currentTree()->outgoingLines().begin();
cjt!=currentTree()->outgoingLines().end();++cjt) {
particlesToShower.push_back(((*cjt).first));
if(ptmax>ZERO) particlesToShower.back()->maximumpT(ptmax,inter);
}
for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
eit=hardTree()->particles().end(),
mit = hardTree()->particles().find(particlesToShower[ix]->progenitor());
if( mit != eit) {
if(mit->second->status()==HardBranching::Outgoing)
particlesToShower[ix]->progenitor()->set5Momentum(mit->second->pVector());
}
}
return true;
}
bool Evolver::constructHardTree(vector<ShowerProgenitorPtr> & particlesToShower,
ShowerInteraction::Type inter) {
bool noEmission = true;
vector<HardBranchingPtr> spaceBranchings,allBranchings;
for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
if(particlesToShower[ix]->progenitor()->isFinalState()) {
HardBranchingPtr newBranch;
if(particlesToShower[ix]->hasEmitted()) {
noEmission = false;
newBranch =
new_ptr(HardBranching(particlesToShower[ix]->progenitor(),
particlesToShower[ix]->progenitor()->
showerKinematics()->SudakovFormFactor(),
HardBranchingPtr(),HardBranching::Outgoing));
constructTimeLikeLine(newBranch,particlesToShower[ix]->progenitor());
}
else {
newBranch =
new_ptr(HardBranching(particlesToShower[ix]->progenitor(),
SudakovPtr(),HardBranchingPtr(),
HardBranching::Outgoing));
}
allBranchings.push_back(newBranch);
}
else {
HardBranchingPtr first,last;
if(!particlesToShower[ix]->progenitor()->parents().empty()) {
noEmission = false;
constructSpaceLikeLine(particlesToShower[ix]->progenitor(),
first,last,SudakovPtr(),
particlesToShower[ix]->original()->parents()[0]);
}
else {
first = new_ptr(HardBranching(particlesToShower[ix]->progenitor(),
SudakovPtr(),HardBranchingPtr(),
HardBranching::Incoming));
if(particlesToShower[ix]->original()->parents().empty())
first->beam(particlesToShower[ix]->original());
else
first->beam(particlesToShower[ix]->original()->parents()[0]);
last = first;
}
spaceBranchings.push_back(first);
allBranchings.push_back(last);
}
}
if(!noEmission) {
HardTreePtr QCDTree = new_ptr(HardTree(allBranchings,spaceBranchings,
inter));
// set the charge partners
ShowerParticleVector particles;
for(set<HardBranchingPtr>::iterator cit=QCDTree->branchings().begin();
cit!=QCDTree->branchings().end();++cit) {
particles.push_back((*cit)->branchingParticle());
}
// get the partners
showerModel()->partnerFinder()->setInitialEvolutionScales(particles,false,
inter,true);
// do the inverse recon
if(!showerModel()->kinematicsReconstructor()->
deconstructHardJets(QCDTree,this,inter))
throw Exception() << "Can't to shower deconstruction for QED shower in"
<< "QEDEvolver::showerHard" << Exception::eventerror;
// set the hard tree
hardTree(QCDTree);
}
// clear the old shower
currentTree()->clear();
// set the charge partners
setEvolutionPartners(true,inter);
// get the particles to be showered
particlesToShower = currentTree()->extractProgenitors();
// reset momenta
if(hardTree()) {
for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
eit=hardTree()->particles().end(),
mit = hardTree()->particles().find(particlesToShower[ix]->progenitor());
if( mit != eit) {
particlesToShower[ix]->progenitor()->set5Momentum(mit->second->showerMomentum());
}
}
}
return true;
}
void Evolver::constructTimeLikeLine(tHardBranchingPtr branch,
tShowerParticlePtr particle) {
for(unsigned int ix=0;ix<particle->children().size();++ix) {
HardBranching::Status status = branch->status();
tShowerParticlePtr child =
dynamic_ptr_cast<ShowerParticlePtr>(particle->children()[ix]);
if(child->children().empty()) {
HardBranchingPtr newBranch =
new_ptr(HardBranching(child,SudakovPtr(),branch,status));
branch->addChild(newBranch);
}
else {
HardBranchingPtr newBranch =
new_ptr(HardBranching(child,child->showerKinematics()->SudakovFormFactor(),
branch,status));
constructTimeLikeLine(newBranch,child);
branch->addChild(newBranch);
}
}
// sort out the type of interaction
if(!branch->children().empty()) {
if(branch->branchingParticle()->id()==ParticleID::gamma ||
branch->children()[0]->branchingParticle()->id()==ParticleID::gamma ||
branch->children()[1]->branchingParticle()->id()==ParticleID::gamma)
branch->type(ShowerPartnerType::QED);
else {
if(branch->branchingParticle()->id()==
branch->children()[0]->branchingParticle()->id()) {
if(branch->branchingParticle()->dataPtr()->iColour()==PDT::Colour8) {
tShowerParticlePtr emittor =
branch->branchingParticle()->showerKinematics()->z()>0.5 ?
branch->children()[0]->branchingParticle() :
branch->children()[1]->branchingParticle();
if(branch->branchingParticle()->colourLine()==emittor->colourLine())
branch->type(ShowerPartnerType::QCDAntiColourLine);
else if(branch->branchingParticle()->antiColourLine()==emittor->antiColourLine())
branch->type(ShowerPartnerType::QCDColourLine);
else
assert(false);
}
else if(branch->branchingParticle()->colourLine()) {
branch->type(ShowerPartnerType::QCDColourLine);
}
else if(branch->branchingParticle()->antiColourLine()) {
branch->type(ShowerPartnerType::QCDAntiColourLine);
}
else
assert(false);
}
else if(branch->branchingParticle()->id()==ParticleID::g &&
branch->children()[0]->branchingParticle()->id()==
-branch->children()[1]->branchingParticle()->id()) {
if(branch->branchingParticle()->showerKinematics()->z()>0.5)
branch->type(ShowerPartnerType::QCDAntiColourLine);
else
branch->type(ShowerPartnerType::QCDColourLine);
}
else
assert(false);
}
}
}
void Evolver::constructSpaceLikeLine(tShowerParticlePtr particle,
HardBranchingPtr & first,
HardBranchingPtr & last,
SudakovPtr sud,PPtr beam) {
if(!particle) return;
if(!particle->parents().empty()) {
tShowerParticlePtr parent =
dynamic_ptr_cast<ShowerParticlePtr>(particle->parents()[0]);
SudakovPtr newSud=particle->showerKinematics()->SudakovFormFactor();
constructSpaceLikeLine(parent,first,last,newSud,beam);
}
HardBranchingPtr newBranch =
new_ptr(HardBranching(particle,sud,last,HardBranching::Incoming));
newBranch->beam(beam);
if(!first) {
first=newBranch;
last =newBranch;
return;
}
last->addChild(newBranch);
tShowerParticlePtr timeChild =
dynamic_ptr_cast<ShowerParticlePtr>(particle->parents()[0]->children()[1]);
HardBranchingPtr timeBranch;
if(!timeChild->children().empty()) {
timeBranch =
new_ptr(HardBranching(timeChild,
timeChild->showerKinematics()->SudakovFormFactor(),
last,HardBranching::Outgoing));
constructTimeLikeLine(timeBranch,timeChild);
}
else {
timeBranch =
new_ptr(HardBranching(timeChild,SudakovPtr(),last,HardBranching::Outgoing));
}
last->addChild(timeBranch);
last=newBranch;
}
void Evolver::connectTrees(ShowerTreePtr showerTree,
HardTreePtr hardTree, bool hard ) {
ShowerParticleVector particles;
// find the Sudakovs
for(set<HardBranchingPtr>::iterator cit=hardTree->branchings().begin();
cit!=hardTree->branchings().end();++cit) {
// Sudakovs for ISR
if((**cit).parent()&&(**cit).status()==HardBranching::Incoming) {
++_nis;
IdList br(3);
br[0] = (**cit).parent()->branchingParticle()->id();
br[1] = (**cit). branchingParticle()->id();
br[2] = (**cit).parent()->children()[0]==*cit ?
(**cit).parent()->children()[1]->branchingParticle()->id() :
(**cit).parent()->children()[0]->branchingParticle()->id();
BranchingList branchings = splittingGenerator()->initialStateBranchings();
if(br[1]<0&&br[0]==br[1]) {
br[0] = abs(br[0]);
br[1] = abs(br[1]);
}
else if(br[1]<0) {
br[1] = -br[1];
br[2] = -br[2];
}
long index = abs(br[1]);
SudakovPtr sudakov;
for(BranchingList::const_iterator cjt = branchings.lower_bound(index);
cjt != branchings.upper_bound(index); ++cjt ) {
IdList ids = cjt->second.second;
if(ids[0]==br[0]&&ids[1]==br[1]&&ids[2]==br[2]) {
sudakov=cjt->second.first;
break;
}
}
if(!sudakov) throw Exception() << "Can't find Sudakov for the hard emission in "
<< "Evolver::connectTrees() for ISR"
<< Exception::runerror;
(**cit).parent()->sudakov(sudakov);
}
// Sudakovs for FSR
else if(!(**cit).children().empty()) {
++_nfs;
IdList br(3);
br[0] = (**cit) .branchingParticle()->id();
br[1] = (**cit).children()[0]->branchingParticle()->id();
br[2] = (**cit).children()[1]->branchingParticle()->id();
BranchingList branchings = splittingGenerator()->finalStateBranchings();
if(br[0]<0) {
br[0] = abs(br[0]);
br[1] = abs(br[1]);
br[2] = abs(br[2]);
}
long index = br[0];
SudakovPtr sudakov;
for(BranchingList::const_iterator cjt = branchings.lower_bound(index);
cjt != branchings.upper_bound(index); ++cjt ) {
IdList ids = cjt->second.second;
if(ids[0]==br[0]&&ids[1]==br[1]&&ids[2]==br[2]) {
sudakov=cjt->second.first;
break;
}
}
if(!sudakov) throw Exception() << "Can't find Sudakov for the hard emission in "
<< "Evolver::connectTrees()"
<< Exception::runerror;
(**cit).sudakov(sudakov);
}
}
// calculate the evolution scale
for(set<HardBranchingPtr>::iterator cit=hardTree->branchings().begin();
cit!=hardTree->branchings().end();++cit) {
particles.push_back((*cit)->branchingParticle());
}
showerModel()->partnerFinder()->
setInitialEvolutionScales(particles,!hard,hardTree->interaction(),
!hardTree->partnersSet());
hardTree->partnersSet(true);
// inverse reconstruction
if(hard)
showerModel()->kinematicsReconstructor()->
deconstructHardJets(hardTree,ShowerHandler::currentHandler()->evolver(),
hardTree->interaction());
else
showerModel()->kinematicsReconstructor()->
deconstructDecayJets(hardTree,ShowerHandler::currentHandler()->evolver(),
hardTree->interaction());
// now reset the momenta of the showering particles
vector<ShowerProgenitorPtr> particlesToShower;
for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
cit=showerTree->incomingLines().begin();
cit!=showerTree->incomingLines().end();++cit )
particlesToShower.push_back(cit->first);
// extract the showering particles
for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
cit=showerTree->outgoingLines().begin();
cit!=showerTree->outgoingLines().end();++cit )
particlesToShower.push_back(cit->first);
// match them
vector<bool> matched(particlesToShower.size(),false);
for(set<HardBranchingPtr>::const_iterator cit=hardTree->branchings().begin();
cit!=hardTree->branchings().end();++cit) {
Energy2 dmin( 1e30*GeV2 );
int iloc(-1);
for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
if(matched[ix]) continue;
if( (**cit).branchingParticle()->id() != particlesToShower[ix]->progenitor()->id() ) continue;
if( (**cit).branchingParticle()->isFinalState() !=
particlesToShower[ix]->progenitor()->isFinalState() ) continue;
Energy2 dtest =
sqr( particlesToShower[ix]->progenitor()->momentum().x() - (**cit).showerMomentum().x() ) +
sqr( particlesToShower[ix]->progenitor()->momentum().y() - (**cit).showerMomentum().y() ) +
sqr( particlesToShower[ix]->progenitor()->momentum().z() - (**cit).showerMomentum().z() ) +
sqr( particlesToShower[ix]->progenitor()->momentum().t() - (**cit).showerMomentum().t() );
// add mass difference for identical particles (e.g. Z0 Z0 production)
dtest += 1e10*sqr(particlesToShower[ix]->progenitor()->momentum().m()-
(**cit).showerMomentum().m());
if( dtest < dmin ) {
iloc = ix;
dmin = dtest;
}
}
if(iloc<0) throw Exception() << "Failed to match shower and hard trees in Evolver::hardestEmission"
<< Exception::eventerror;
particlesToShower[iloc]->progenitor()->set5Momentum((**cit).showerMomentum());
matched[iloc] = true;
}
// correction boosts for daughter trees
for(map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator
tit = showerTree->treelinks().begin();
tit != showerTree->treelinks().end();++tit) {
ShowerTreePtr decayTree = tit->first;
map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
cit = decayTree->incomingLines().begin();
// reset the momentum of the decay particle
Lorentz5Momentum oldMomentum = cit->first->progenitor()->momentum();
Lorentz5Momentum newMomentum = tit->second.second->momentum();
LorentzRotation boost( oldMomentum.findBoostToCM(),oldMomentum.e()/oldMomentum.mass());
boost.boost (-newMomentum.findBoostToCM(),newMomentum.e()/newMomentum.mass());
decayTree->transform(boost,true);
}
}
void Evolver::doShowering(bool hard) {
// order of the interactions
bool showerOrder(true);
// zero number of emissions
_nis = _nfs = 0;
// extract particles to shower
vector<ShowerProgenitorPtr> particlesToShower(setupShower(hard));
// setup the maximum scales for the shower
setupMaximumScales(particlesToShower);
// specifc stuff for hard processes and decays
Energy minmass(ZERO), mIn(ZERO);
// hard process generate the intrinsic p_T once and for all
if(hard) {
generateIntrinsicpT(particlesToShower);
}
// decay compute the minimum mass of the final-state
else {
for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
if(particlesToShower[ix]->progenitor()->isFinalState()) {
minmass += max( particlesToShower[ix]->progenitor()->mass(),
particlesToShower[ix]->progenitor()->dataPtr()->constituentMass() );
}
else {
mIn = particlesToShower[ix]->progenitor()->mass();
}
}
// throw exception if decay can'ty happen
if ( minmass > mIn ) {
throw Exception() << "Evolver.cc: Mass of decaying particle is "
<< "below constituent masses of decay products."
<< Exception::eventerror;
}
}
// check if interactions in right order
if(hardTree() && hardTree()->interaction()!=interactions_[0]) {
showerOrder = false;
swap(interactions_[0],interactions_[1]);
}
// loop over possible interactions
for(unsigned int inter=0;inter<interactions_.size();++inter) {
// set up for second pass if required
if(inter!=0) {
// zero intrinsic pt so only added first time round
intrinsicpT().clear();
// construct the tree and throw veto if not possible
if(!(hard ?
constructHardTree (particlesToShower,interactions_[inter]) :
constructDecayTree(particlesToShower,interactions_[inter])))
throw InteractionVeto();
}
// main shower loop
unsigned int ntry(0);
bool reconstructed = false;
do {
// clear results of last attempt if needed
if(ntry!=0) {
currentTree()->clear();
setEvolutionPartners(hard,interactions_[inter]);
_nis = _nfs = 0;
}
// generate the shower
// pick random starting point
unsigned int istart=UseRandom::irnd(particlesToShower.size());
unsigned int istop = particlesToShower.size();
// loop over particles with random starting point
for(unsigned int ix=istart;ix<=istop;++ix) {
if(ix==particlesToShower.size()) {
if(istart!=0) {
istop = istart-1;
ix=0;
}
else break;
}
// extract the progenitor
progenitor(particlesToShower[ix]);
// final-state radiation
if(progenitor()->progenitor()->isFinalState()) {
if(!isFSRadiationON()) continue;
// perform shower
progenitor()->hasEmitted(startTimeLikeShower(interactions_[inter]));
}
// initial-state radiation
else {
if(!isISRadiationON()) continue;
// hard process
if(hard) {
// get the PDF
setBeamParticle(_progenitor->beam());
assert(beamParticle());
// perform the shower
// set the beam particle
tPPtr beamparticle=progenitor()->original();
if(!beamparticle->parents().empty())
beamparticle=beamparticle->parents()[0];
// generate the shower
progenitor()->hasEmitted(startSpaceLikeShower(beamparticle,
interactions_[inter]));
}
// decay
else {
// skip colour and electrically neutral particles
if(!progenitor()->progenitor()->dataPtr()->coloured() &&
!progenitor()->progenitor()->dataPtr()->charged()) {
progenitor()->hasEmitted(false);
continue;
}
// perform shower
// set the scales correctly. The current scale is the maximum scale for
// emission not the starting scale
- map<ShowerPartnerType::Type,pair<Energy,Energy> >
- maxScales = progenitor()->progenitor()->evolutionScales();
- // starting scale is mass
- Energy startScale=progenitor()->progenitor()->mass();
- for(map<ShowerPartnerType::Type,pair<Energy,Energy> >::iterator it=maxScales.begin();
- it!=maxScales.end();++it)
- progenitor()->progenitor()->evolutionScale(it->first,make_pair(startScale,startScale));
+ map<ShowerPartnerType::Type,pair<Energy,Energy> > maxScales;
+ // map<ShowerPartnerType::Type,pair<Energy,Energy> >
+ // maxScales = progenitor()->progenitor()->evolutionScales();
+ // // starting scale is mass
+ // Energy startScale=progenitor()->progenitor()->mass();
+ // for(map<ShowerPartnerType::Type,pair<Energy,Energy> >::iterator it=maxScales.begin();
+ // it!=maxScales.end();++it)
+ // progenitor()->progenitor()->evolutionScale(it->first,make_pair(startScale,startScale));
+ assert(false);
// perform the shower
progenitor()->hasEmitted(startSpaceLikeDecayShower(maxScales,minmass,
interactions_[inter]));
}
}
}
// do the kinematic reconstruction, checking if it worked
reconstructed = hard ?
showerModel()->kinematicsReconstructor()->
reconstructHardJets (currentTree(),intrinsicpT(),interactions_[inter]) :
showerModel()->kinematicsReconstructor()->
reconstructDecayJets(currentTree(),interactions_[inter]);
}
while(!reconstructed&&maximumTries()>++ntry);
// check if failed to generate the shower
if(ntry==maximumTries()) {
if(hard)
throw ShowerHandler::ShowerTriesVeto(ntry);
else
throw Exception() << "Failed to generate the shower after "
<< ntry << " attempts in Evolver::showerDecay()"
<< Exception::eventerror;
}
}
// tree has now showered
_currenttree->hasShowered(true);
if(!showerOrder) swap(interactions_[0],interactions_[1]);
hardTree(HardTreePtr());
}
diff --git a/Shower/Base/PartnerFinder.cc b/Shower/Base/PartnerFinder.cc
--- a/Shower/Base/PartnerFinder.cc
+++ b/Shower/Base/PartnerFinder.cc
@@ -1,523 +1,523 @@
// -*- C++ -*-
//
// PartnerFinder.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 PartnerFinder class.
//
#include "PartnerFinder.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ShowerParticle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Utilities/Debug.h"
using namespace Herwig;
// some useful functions to avoid using #define
namespace {
// return bool if final-state particle
inline bool FS(const tShowerParticlePtr a) {
return a->isFinalState();
}
// return colour line pointer
inline Ptr<ThePEG::ColourLine>::transient_pointer
CL(const tShowerParticlePtr a, unsigned int index=0) {
return a->colourInfo()->colourLines().empty() ? ThePEG::tColinePtr() :
const_ptr_cast<ThePEG::tColinePtr>(a->colourInfo()->colourLines()[index]);
}
// return colour line size
inline unsigned int
CLSIZE(const tShowerParticlePtr a) {
return a->colourInfo()->colourLines().size();
}
inline Ptr<ThePEG::ColourLine>::transient_pointer
ACL(const tShowerParticlePtr a, unsigned int index=0) {
return a->colourInfo()->antiColourLines().empty() ? ThePEG::tColinePtr() :
const_ptr_cast<ThePEG::tColinePtr>(a->colourInfo()->antiColourLines()[index]);
}
inline unsigned int
ACLSIZE(const tShowerParticlePtr a) {
return a->colourInfo()->antiColourLines().size();
}
}
void PartnerFinder::persistentOutput(PersistentOStream & os) const {
os << partnerMethod_ << QEDPartner_ << scaleChoice_;
}
void PartnerFinder::persistentInput(PersistentIStream & is, int) {
is >> partnerMethod_ >> QEDPartner_ >> scaleChoice_;
}
AbstractClassDescription<PartnerFinder> PartnerFinder::initPartnerFinder;
// Definition of the static class description member.
void PartnerFinder::Init() {
static ClassDocumentation<PartnerFinder> documentation
("This class is responsible for finding the partners for each interaction types ",
"and within the evolution scale range specified by the ShowerVariables ",
"then to determine the initial evolution scales for each pair of partners.");
static Switch<PartnerFinder,int> interfacePartnerMethod
("PartnerMethod",
"Choice of partner finding method for gluon evolution.",
&PartnerFinder::partnerMethod_, 0, false, false);
static SwitchOption interfacePartnerMethodRandom
(interfacePartnerMethod,
"Random",
"Choose partners of a gluon randomly.",
0);
static SwitchOption interfacePartnerMethodMaximum
(interfacePartnerMethod,
"Maximum",
"Choose partner of gluon with largest angle.",
1);
static Switch<PartnerFinder,int> interfaceQEDPartner
("QEDPartner",
"Control of which particles to use as the partner for QED radiation",
&PartnerFinder::QEDPartner_, 0, false, false);
static SwitchOption interfaceQEDPartnerAll
(interfaceQEDPartner,
"All",
"Consider all possible choices which give a positive contribution"
" in the soft limit.",
0);
static SwitchOption interfaceQEDPartnerIIandFF
(interfaceQEDPartner,
"IIandFF",
"Only allow initial-initial or final-final combinations",
1);
static SwitchOption interfaceQEDPartnerIF
(interfaceQEDPartner,
"IF",
"Only allow initial-final combinations",
2);
static Switch<PartnerFinder,int> interfaceScaleChoice
("ScaleChoice",
"The choice of the evolution scales",
&PartnerFinder::scaleChoice_, 0, false, false);
static SwitchOption interfaceScaleChoicePartner
(interfaceScaleChoice,
"Partner",
"Scale of all interactions is that of the evolution partner",
0);
static SwitchOption interfaceScaleChoiceDifferent
(interfaceScaleChoice,
"Different",
"Allow each interaction to have different scales",
1);
}
void PartnerFinder::setInitialEvolutionScales(const ShowerParticleVector &particles,
const bool isDecayCase,
ShowerInteraction::Type type,
const bool setPartners) {
// clear the existing partners
for(ShowerParticleVector::const_iterator cit = particles.begin();
cit != particles.end(); ++cit) (*cit)->clearPartners();
// set them
if(type==ShowerInteraction::QCD) {
setInitialQCDEvolutionScales(particles,isDecayCase,setPartners);
}
else if(type==ShowerInteraction::QED) {
setInitialQEDEvolutionScales(particles,isDecayCase,setPartners);
}
else {
setInitialQCDEvolutionScales(particles,isDecayCase,setPartners);
setInitialQEDEvolutionScales(particles,isDecayCase,false);
}
// print out for debugging
if(Debug::level>=10) {
for(ShowerParticleVector::const_iterator cit = particles.begin();
cit != particles.end(); ++cit) {
generator()->log() << "Particle: " << **cit << "\n";
if(!(**cit).partner()) continue;
generator()->log() << "Primary partner: " << *(**cit).partner() << "\n";
for(vector<ShowerParticle::EvolutionPartner>::const_iterator it= (**cit).partners().begin();
it!=(**cit).partners().end();++it) {
generator()->log() << it->type << " " << it->weight << " "
<< it->scale/GeV << " "
<< *(it->partner)
<< "\n";
}
}
generator()->log() << flush;
}
}
void PartnerFinder::setInitialQCDEvolutionScales(const ShowerParticleVector &particles,
const bool isDecayCase,
const bool setPartners) {
// set the partners and the scales
if(setPartners) {
// Loop over particles and consider only coloured particles which don't
// have already their colour partner fixed and that don't have children
// (the latter requirement is relaxed in the case isDecayCase is true).
// Build a map which has as key one of these particles (i.e. a pointer
// to a ShowerParticle object) and as a corresponding value the vector
// of all its possible *normal* candidate colour partners, defined as follows:
// --- have colour, and no children (this is not required in the case
// isDecayCase is true);
// --- if both are initial (incoming) state particles, then the (non-null) colourLine()
// of one of them must match the (non-null) antiColourLine() of the other.
// --- if one is an initial (incoming) state particle and the other is
// a final (outgoing) state particle, then both must have the
// same (non-null) colourLine() or the same (non-null) antiColourLine();
// Notice that this definition exclude the special case of baryon-violating
// processes (as in R-parity Susy), which will show up as particles
// without candidate colour partners, and that we will be treated a part later
// (this means that no modifications of the following loop is needed!)
ShowerParticleVector::const_iterator cit, cjt;
for(cit = particles.begin(); cit != particles.end(); ++cit) {
// Skip colourless particles
if(!(*cit)->data().coloured()) continue;
// find the partners
vector< pair<ShowerPartnerType::Type, tShowerParticlePtr> > partners =
findQCDPartners(*cit,particles);
// must have a partner
if(partners.empty()) {
throw Exception() << "`Failed to make colour connections in "
<< "PartnerFinder::setQCDInitialEvolutionScales"
<< (**cit)
<< Exception::eventerror;
}
// Calculate the evolution scales for all possible pairs of of particles
vector<pair<Energy,Energy> > scales;
for(unsigned int ix=0;ix< partners.size();++ix) {
scales.push_back(calculateInitialEvolutionScales(ShowerPPair(*cit,partners[ix].second),
isDecayCase));
}
// In the case of more than one candidate colour partners,
// there are now two approaches to choosing the partner. The
// first method is based on two assumptions:
// 1) the choice of which is the colour partner is done
// *randomly* between the available candidates.
// 2) the choice of which is the colour partner is done
// *independently* from each particle: in other words,
// if for a particle "i" its selected colour partner is
// the particle "j", then the colour partner of "j"
// does not have to be necessarily "i".
// The second method always chooses the furthest partner
// for hard gluons and gluinos.
// store the choice
int position(-1);
// random choice
if( partnerMethod_ == 0 ) {
// random choice of partner
position = UseRandom::irnd(partners.size());
}
// take the one with largest angle
else if (partnerMethod_ == 1 ) {
if ((*cit)->perturbative() == 1 &&
(*cit)->dataPtr()->iColour()==PDT::Colour8 ) {
assert(partners.size()==2);
// Determine largest angle
double maxAngle(0.);
for(unsigned int ix=0;ix<partners.size();++ix) {
double angle = (*cit)->momentum().vect().
angle(partners[ix].second->momentum().vect());
if(angle>maxAngle) {
maxAngle = angle;
position = ix;
}
}
}
}
else
assert(false);
// set the evolution partner
(*cit)->partner(partners[position].second);
for(unsigned int ix=0;ix<partners.size();++ix) {
(**cit).addPartner(ShowerParticle::EvolutionPartner(partners[ix].second,
1.,partners[ix].first,
scales[ix].first));
}
// now scales
// set scales for all interactions to that of the partner, default
- if(scaleChoice_==0) {
- pair<Energy,Energy> scalePair(scales[position].first,scales[position].first);
- for(unsigned int ix=0;ix<partners.size();++ix) {
- (**cit).evolutionScale(partners[ix].first,scalePair);
- }
- }
- // leave them as is
- else {
- for(unsigned int ix=0;ix<partners.size();++ix) {
-// pair<Energy,Energy> scale = (**cit).evolutionScale(partners[ix].first);
-// if(scales[ix].first<=scale.first)
-// continue;
- (**cit).evolutionScale(partners[ix].first,make_pair(scales[ix].first,
- scales[ix].first));
- }
- }
+ assert(false);
+ // if(scaleChoice_==0) {
+ // pair<Energy,Energy> scalePair(scales[position].first,scales[position].first);
+ // for(unsigned int ix=0;ix<partners.size();++ix) {
+ // (**cit).evolutionScale(partners[ix].first,scalePair);
+ // }
+ // }
+ // // leave them as is
+ // else {
+ // for(unsigned int ix=0;ix<partners.size();++ix) {
+ // (**cit).evolutionScale(partners[ix].first,make_pair(scales[ix].first,
+ // scales[ix].first));
+ // }
+ // }
}
}
// primary partner set, set the others and do the scale
else {
for(ShowerParticleVector::const_iterator cit = particles.begin();
cit != particles.end(); ++cit) {
// Skip colourless particles
if(!(*cit)->data().coloured()) continue;
// find the partners
vector< pair<ShowerPartnerType::Type, tShowerParticlePtr> > partners =
findQCDPartners(*cit,particles);
// must have a partner
if(partners.empty()) {
throw Exception() << "`Failed to make colour connections in "
<< "PartnerFinder::setQCDInitialEvolutionScales"
<< (**cit)
<< Exception::eventerror;
}
// Calculate the evolution scales for all possible pairs of of particles
vector<pair<Energy,Energy> > scales;
int position(-1);
for(unsigned int ix=0;ix< partners.size();++ix) {
if(partners[ix].second) position = ix;
scales.push_back(calculateInitialEvolutionScales(ShowerPPair(*cit,partners[ix].second),
isDecayCase));
}
assert(position>=0);
for(unsigned int ix=0;ix<partners.size();++ix) {
(**cit).addPartner(ShowerParticle::EvolutionPartner(partners[ix].second,
1.,partners[ix].first,
scales[ix].first));
}
// now scales
- // set scales for all interactions to that of the partner, default
- if(scaleChoice_==0) {
- pair<Energy,Energy> scalePair(scales[position].first,scales[position].first);
- for(unsigned int ix=0;ix<partners.size();++ix) {
- (**cit).evolutionScale(partners[ix].first,scalePair);
- }
- }
- // leave them as is
- else {
- for(unsigned int ix=0;ix<partners.size();++ix) {
- pair<Energy,Energy> scale = (**cit).evolutionScale(partners[ix].first);
- if(scales[ix].first<=scale.first)
- continue;
- (**cit).evolutionScale(partners[ix].first,make_pair(scales[ix].first,
- scales[ix].first));
- }
- }
+ // // set scales for all interactions to that of the partner, default
+ // if(scaleChoice_==0) {
+ // pair<Energy,Energy> scalePair(scales[position].first,scales[position].first);
+ // for(unsigned int ix=0;ix<partners.size();++ix) {
+ // (**cit).evolutionScale(partners[ix].first,scalePair);
+ // }
+ // }
+ // // leave them as is
+ // else {
+ // for(unsigned int ix=0;ix<partners.size();++ix) {
+ // pair<Energy,Energy> scale = (**cit).evolutionScale(partners[ix].first);
+ // if(scales[ix].first<=scale.first)
+ // continue;
+ // (**cit).evolutionScale(partners[ix].first,make_pair(scales[ix].first,
+ // scales[ix].first));
+ // }
+ // }
+ assert(false);
}
}
}
void PartnerFinder::setInitialQEDEvolutionScales(const ShowerParticleVector &particles,
const bool isDecayCase,
const bool setPartners) {
// loop over all the particles
for(ShowerParticleVector::const_iterator cit = particles.begin();
cit != particles.end(); ++cit) {
// not charge continue
if(!(**cit).dataPtr()->charged()) continue;
// find the potential partners
vector<pair<double,tShowerParticlePtr> > partners = findQEDPartners(*cit,particles);
if(partners.empty()) {
throw Exception() << "Failed to partner in "
<< "PartnerFinder::setQEDInitialEvolutionScales"
<< (**cit) << Exception::eventerror;
}
// calculate the probabilities
double prob(0.);
for(unsigned int ix=0;ix<partners.size();++ix) prob += partners[ix].first;
// normalise
for(unsigned int ix=0;ix<partners.size();++ix) partners[ix].first /= prob;
// set the partner if required
int position(-1);
if(setPartners||!(*cit)->partner()) {
prob = UseRandom::rnd();
for(unsigned int ix=0;ix<partners.size();++ix) {
if(partners[ix].first>prob) {
position = ix;
break;
}
prob -= partners[ix].first;
}
if(position>=0) (*cit)->partner(partners[position].second);
}
// partner set, find it
else {
for(unsigned int ix=0;ix<partners.size();++ix) {
if((*cit)->partner()==partners[ix].second) {
position = ix;
break;
}
}
}
// must have a partner
if(position<0) throw Exception() << "Failed to partner in "
<< "PartnerFinder::setQEDInitialEvolutionScales"
<< (**cit) << Exception::eventerror;
// Calculate the evolution scales for all possible pairs of of particles
vector<pair<Energy,Energy> > scales;
for(unsigned int ix=0;ix< partners.size();++ix) {
scales.push_back(calculateInitialEvolutionScales(ShowerPPair(*cit,partners[ix].second),
isDecayCase));
}
// store all the possible partners
for(unsigned int ix=0;ix<partners.size();++ix) {
(**cit).addPartner(ShowerParticle::EvolutionPartner(partners[ix].second,
partners[ix].first,
ShowerPartnerType::QED,
scales[ix].first));
}
// set scales
- pair<Energy,Energy> scalePair(scales[position].first,scales[position].first);
- (**cit).evolutionScale(ShowerPartnerType::QED,scalePair);
+ // pair<Energy,Energy> scalePair(scales[position].first,scales[position].first);
+ // (**cit).evolutionScale(ShowerPartnerType::QED,scalePair);
+ assert(false);
}
}
pair<Energy,Energy> PartnerFinder::
calculateInitialEvolutionScales(const ShowerPPair &particlePair,
const bool isDecayCase) {
bool FS1=FS(particlePair.first),FS2= FS(particlePair.second);
if(FS1 && FS2)
return calculateFinalFinalScales(particlePair);
else if(FS1 && !FS2) {
ShowerPPair a(particlePair.second, particlePair.first);
pair<Energy,Energy> rval = calculateInitialFinalScales(a,isDecayCase);
return pair<Energy,Energy>(rval.second,rval.first);
}
else if(!FS1 &&FS2)
return calculateInitialFinalScales(particlePair,isDecayCase);
else
return calculateInitialInitialScales(particlePair);
}
vector< pair<ShowerPartnerType::Type, tShowerParticlePtr> >
PartnerFinder::findQCDPartners(tShowerParticlePtr particle,
const ShowerParticleVector &particles) {
vector< pair<ShowerPartnerType::Type, tShowerParticlePtr> > partners;
ShowerParticleVector::const_iterator cjt;
for(cjt = particles.begin(); cjt != particles.end(); ++cjt) {
if(!(*cjt)->data().coloured() || particle==*cjt) continue;
// one initial-state and one final-state particle
if(FS(particle) != FS(*cjt)) {
// loop over all the colours of both particles
for(unsigned int ix=0; ix<CLSIZE(particle); ++ix) {
for(unsigned int jx=0; jx<CLSIZE(*cjt); ++jx) {
if((CL(particle,ix) && CL(particle,ix)==CL(*cjt,jx))) {
partners.push_back(make_pair(ShowerPartnerType:: QCDColourLine,*cjt));
}
}
}
//loop over all the anti-colours of both particles
for(unsigned int ix=0; ix<ACLSIZE(particle); ++ix) {
for(unsigned int jx=0; jx<ACLSIZE(*cjt); ++jx) {
if((ACL(particle,ix) && ACL(particle,ix)==ACL(*cjt,jx))) {
partners.push_back(make_pair(ShowerPartnerType::QCDAntiColourLine,*cjt));
}
}
}
}
// two initial-state or two final-state particles
else {
//loop over the colours of the first particle and the anti-colours of the other
for(unsigned int ix=0; ix<CLSIZE(particle); ++ix){
for(unsigned int jx=0; jx<ACLSIZE(*cjt); ++jx){
if(CL(particle,ix) && CL(particle,ix)==ACL(*cjt,jx)) {
partners.push_back(make_pair(ShowerPartnerType:: QCDColourLine,*cjt));
}
}
}
//loop over the anti-colours of the first particle and the colours of the other
for(unsigned int ix=0; ix<ACLSIZE(particle); ++ix){
for(unsigned int jx=0; jx<CLSIZE(*cjt); jx++){
if(ACL(particle,ix) && ACL(particle,ix)==CL(*cjt,jx)) {
partners.push_back(make_pair(ShowerPartnerType::QCDAntiColourLine,*cjt));
}
}
}
}
}
// if we haven't found any partners look for RPV
if (partners.empty()) {
// special for RPV
tColinePtr col = CL(particle);
if(FS(particle)&&col&&col->sourceNeighbours().first) {
tColinePair cpair = col->sourceNeighbours();
for(cjt=particles.begin();cjt!=particles.end();++cjt) {
if(( FS(*cjt) && ( CL(*cjt) == cpair.first || CL(*cjt) == cpair.second))||
(!FS(*cjt) && (ACL(*cjt) == cpair.first || ACL(*cjt) == cpair.second ))) {
partners.push_back(make_pair(ShowerPartnerType:: QCDColourLine,*cjt));
}
}
}
else if(col&&col->sinkNeighbours().first) {
tColinePair cpair = col->sinkNeighbours();
for(cjt=particles.begin();cjt!=particles.end();++cjt) {
if(( FS(*cjt) && (ACL(*cjt) == cpair.first || ACL(*cjt) == cpair.second))||
(!FS(*cjt) && ( CL(*cjt) == cpair.first || CL(*cjt) == cpair.second))) {
partners.push_back(make_pair(ShowerPartnerType:: QCDColourLine,*cjt));
}
}
}
col = ACL(particle);
if(FS(particle)&&col&&col->sinkNeighbours().first) {
tColinePair cpair = col->sinkNeighbours();
for(cjt=particles.begin();cjt!=particles.end();++cjt) {
if(( FS(*cjt) && (ACL(*cjt) == cpair.first || ACL(*cjt) == cpair.second))||
(!FS(*cjt) && ( CL(*cjt) == cpair.first || CL(*cjt) == cpair.second ))) {
partners.push_back(make_pair(ShowerPartnerType::QCDAntiColourLine,*cjt));
}
}
}
else if(col&&col->sourceNeighbours().first) {
tColinePair cpair = col->sourceNeighbours();
for(cjt=particles.begin();cjt!=particles.end();++cjt) {
if(( FS(*cjt) && ( CL(*cjt) == cpair.first || CL(*cjt) == cpair.second))||
(!FS(*cjt) && (ACL(*cjt) == cpair.first ||ACL(*cjt) == cpair.second))) {
partners.push_back(make_pair(ShowerPartnerType::QCDAntiColourLine,*cjt));
}
}
}
}
// return the partners
return partners;
}
vector< pair<double, tShowerParticlePtr> >
PartnerFinder::findQEDPartners(tShowerParticlePtr particle,
const ShowerParticleVector &particles) {
vector< pair<double, tShowerParticlePtr> > partners;
ShowerParticleVector::const_iterator cjt;
for(cjt = particles.begin(); cjt != particles.end(); ++cjt) {
if(!(*cjt)->data().charged() || particle == *cjt) continue;
double charge = double(particle->data().iCharge()*(*cjt)->data().iCharge());
if( FS(particle) != FS(*cjt) ) charge *=-1.;
if( QEDPartner_ != 0 ) {
// only include II and FF as requested
if( QEDPartner_ == 1 && FS(particle) != FS(*cjt) )
continue;
// ony include IF is requested
else if(QEDPartner_ == 2 && FS(particle) == FS(*cjt) )
continue;
}
// only keep positive dipoles
if(charge<0.) partners.push_back(make_pair(-charge,*cjt));
}
return partners;
}
diff --git a/Shower/Base/ShowerParticle.cc b/Shower/Base/ShowerParticle.cc
--- a/Shower/Base/ShowerParticle.cc
+++ b/Shower/Base/ShowerParticle.cc
@@ -1,48 +1,49 @@
// -*- C++ -*-
//
// ShowerParticle.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the ShowerParticle class.
//
#include "ShowerParticle.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Repository/CurrentGenerator.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
PPtr ShowerParticle::clone() const {
return new_ptr(*this);
}
PPtr ShowerParticle::fullclone() const {
return new_ptr(*this);
}
ClassDescription<ShowerParticle> ShowerParticle::initShowerParticle;
// Definition of the static class description member.
void ShowerParticle::vetoEmission(ShowerPartnerType::Type type, Energy scale) {
- for(map<ShowerPartnerType::Type ,pair<Energy,Energy> >::iterator
- it=evolutionScales_.begin();it!=evolutionScales_.end();++it) {
- if(it->first==type) {
- it->second = make_pair(scale,scale);
- }
- else {
- if(it->second.first >scale) it->second.first = scale;
- if(it->second.second>scale) it->second.second = scale;
- }
- }
+ // for(map<ShowerPartnerType::Type ,pair<Energy,Energy> >::iterator
+ // it=evolutionScales_.begin();it!=evolutionScales_.end();++it) {
+ // if(it->first==type) {
+ // it->second = make_pair(scale,scale);
+ // }
+ // else {
+ // if(it->second.first >scale) it->second.first = scale;
+ // if(it->second.second>scale) it->second.second = scale;
+ // }
+ // }
+ assert(false);
}
void ShowerParticle::addPartner(EvolutionPartner in ) {
partners_.push_back(in);
}
diff --git a/Shower/Base/ShowerParticle.h b/Shower/Base/ShowerParticle.h
--- a/Shower/Base/ShowerParticle.h
+++ b/Shower/Base/ShowerParticle.h
@@ -1,425 +1,440 @@
// -*- C++ -*-
//
// ShowerParticle.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_ShowerParticle_H
#define HERWIG_ShowerParticle_H
//
// This is the declaration of the ShowerParticle class.
//
#include "ThePEG/EventRecord/Particle.h"
#include "Herwig++/Shower/SplittingFunctions/SplittingFunction.fh"
#include "Herwig++/Shower/ShowerConfig.h"
#include "ShowerKinematics.h"
#include "ShowerParticle.fh"
namespace Herwig {
using namespace ThePEG;
/** \ingroup Shower
* This class represents a particle in the showering process.
* It inherits from the Particle class of ThePEG and has some
* specifics information useful only during the showering process.
*
* Notice that:
* - for forward evolution, it is clear what is meant by parent/child;
* for backward evolution, however, it depends whether we want
* to keep a physical picture or a Monte-Carlo effective one.
* In the former case, an incoming particle (emitting particle)
* splits into an emitted particle and the emitting particle after
* the emission: the latter two are then children of the
* emitting particle, the parent. In the Monte-Carlo effective
* picture, we have that the particle close to the hard subprocess,
* with higher (space-like) virtuality, splits into an emitted particle
* and the emitting particle at lower virtuality: the latter two are,
* in this case, the children of the first one, the parent. However we
* choose a more physical picture where the new emitting particle is the
* parented of the emitted final-state particle and the original emitting
* particle.
* - the pointer to a SplitFun object is set only in the case
* that the particle has undergone a shower emission. This is similar to
* the case of the decay of a normal Particle where
* the pointer to a Decayer object is set only in the case
* that the particle has undergone to a decay.
* In the case of particle connected directly to the hard subprocess,
* there is no pointer to the hard subprocess, but there is a method
* isFromHardSubprocess() which returns true only in this case.
*
* @see Particle
* @see ShowerConfig
* @see ShowerKinematics
*/
class ShowerParticle: public Particle {
public:
/**
* Struct for all the info on an evolution partner
*/
struct EvolutionPartner {
/**
* Constructor
*/
EvolutionPartner(tShowerParticlePtr p,double w, ShowerPartnerType::Type t,
Energy s) : partner(p), weight(w), type(t), scale(s)
{}
/**
* The partner
*/
tShowerParticlePtr partner;
/**
* Weight
*/
double weight;
/**
* Type
*/
ShowerPartnerType::Type type;
/**
* The assoicated evolution scale
*/
Energy scale;
};
+ /**
+ * Struct to store the evolution scales
+ */
+ struct EvolutionScales {
+
+ /**
+ * Constructor
+ */
+ EvolutionScales() : QED(),QCD_c(),QCD_ac(),
+ QED_noAO(),QCD_c_noAO(),QCD_ac_noAO()
+ {}
+
+ /**
+ * QED scale
+ */
+ Energy QED;
+
+ /**
+ * QCD colour scale
+ */
+ Energy QCD_c;
+
+ /**
+ * QCD anticolour scale
+ */
+ Energy QCD_ac;
+
+ /**
+ * QED scale
+ */
+ Energy QED_noAO;
+
+ /**
+ * QCD colour scale
+ */
+ Energy QCD_c_noAO;
+
+ /**
+ * QCD anticolour scale
+ */
+ Energy QCD_ac_noAO;
+
+ };
+
+
/** @name Construction and descruction functions. */
//@{
/**
* Standard Constructor. Note that the default constructor is
* private - there is no particle without a pointer to a
* ParticleData object.
* @param x the ParticleData object
* @param fs Whether or not the particle is an inital or final-state particle
* @param tls Whether or not the particle initiates a time-like shower
*/
ShowerParticle(tcEventPDPtr x, bool fs, bool tls=false)
: Particle(x), _isFinalState(fs),
_perturbative(0), _initiatesTLS(tls), _x(1.0), _showerKinematics(),
_vMass(ZERO), _thePEGBase() {}
/**
* Copy constructor from a ThePEG Particle
* @param x ThePEG particle
* @param pert Where the particle came from
* @param fs Whether or not the particle is an inital or final-state particle
* @param tls Whether or not the particle initiates a time-like shower
*/
ShowerParticle(const Particle & x, unsigned int pert, bool fs, bool tls=false)
: Particle(x), _isFinalState(fs),
_perturbative(pert), _initiatesTLS(tls), _x(1.0), _showerKinematics(),
_vMass(ZERO), _thePEGBase(&x) {}
//@}
public:
/**
* Access/Set various flags about the state of the particle
*/
//@{
/**
* Access the flag that tells if the particle is final state
* or initial state.
*/
bool isFinalState() const { return _isFinalState; }
/**
* Access the flag that tells if the particle is initiating a
* time like shower when it has been emitted in an initial state shower.
*/
bool initiatesTLS() const { return _initiatesTLS; }
/**
* Access the flag which tells us where the particle came from
* This is 0 for a particle produced in the shower, 1 if the particle came
* from the hard sub-process and 2 is it came from a decay.
*/
unsigned int perturbative() const { return _perturbative; }
//@}
/**
* Set/Get the momentum fraction for initial-state particles
*/
//@{
/**
* For an initial state particle get the fraction of the beam momentum
*/
void x(double x) { _x = x; }
/**
* For an initial state particle set the fraction of the beam momentum
*/
double x() const { return _x; }
//@}
/**
* Set/Get methods for the ShowerKinematics objects
*/
//@{
/**
* Access/ the ShowerKinematics object.
*/
const ShoKinPtr & showerKinematics() const { return _showerKinematics; }
/**
* Set the ShowerKinematics object.
*/
void showerKinematics(const ShoKinPtr in) { _showerKinematics = in; }
//@}
/**
* Members relating to the initial evolution scale and partner for the particle
*/
//@{
/**
- * Return the evolution scale \f$\tilde{q}\f$, the first value
- * includes angular ordering, while the second does not.
- */
- pair<Energy,Energy> evolutionScale(ShowerPartnerType::Type type) const {
- map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
- it = evolutionScales_.find(type);
- if(it!=evolutionScales_.end())
- return it->second;
- else
- return make_pair(ZERO,ZERO);
- }
-
- /**
- * Return the evolution scale \f$\tilde{q}\f$
- */
- Energy evolutionScale(bool AO, ShowerPartnerType::Type type) const {
- map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
- it = evolutionScales_.find(type);
- if(it!=evolutionScales_.end()) {
- return AO ? it->second.first : it->second.second;
- }
- else
- return ZERO;
- }
-
- /**
* Veto emission at a given scale
*/
void vetoEmission(ShowerPartnerType::Type type, Energy scale);
/**
- * Set the evolution \f$\tilde{q}\f$ scale
+ * Access to the evolution scales
*/
- void evolutionScale(ShowerPartnerType::Type type, pair<Energy,Energy> scale) {
- evolutionScales_[type] = scale;
- }
+ const EvolutionScales & scales() const {return scales_;}
/**
- * Access to the map of evolution scales
+ * Access to the evolution scales
*/
- map<ShowerPartnerType::Type,pair<Energy,Energy> > & evolutionScales() {
- return evolutionScales_;
- }
+ EvolutionScales & scales() {return scales_;}
/**
* Return the virtual mass\f$
*/
Energy virtualMass() const { return _vMass; }
/**
* Set the virtual mass
*/
void virtualMass(Energy mass) { _vMass = mass; }
/**
* Return the partner
*/
tShowerParticlePtr partner() const { return _partner; }
/**
* Set the partner
*/
void partner(const tShowerParticlePtr partner) { _partner = partner; }
/**
* Get the possible partners
*/
vector<EvolutionPartner> & partners() { return partners_; }
/**
* Add a possible partners
*/
void addPartner(EvolutionPartner in );
/**
* Clear the evolution partners
*/
void clearPartners() { partners_.clear(); }
/**
* Return the progenitor of the shower
*/
tShowerParticlePtr progenitor() const { return _progenitor; }
/**
* Set the progenitor of the shower
*/
void progenitor(const tShowerParticlePtr progenitor) { _progenitor = progenitor; }
//@}
/**
* Members to store and provide access to variables for a specific
* shower evolution scheme
*/
//@{
/**
* Set the vector containing dimensionless variables
*/
vector<double> & showerParameters() { return _parameters; }
/**
* Set the vector containing dimensionful variables
*/
vector<Energy> & showerVariables() { return _variables; }
//@}
/**
* If this particle came from the hard process get a pointer to ThePEG particle
* it came from
*/
const tcPPtr thePEGBase() const { return _thePEGBase; }
protected:
/**
* Standard clone function.
*/
virtual PPtr clone() const;
/**
* Standard clone function.
*/
virtual PPtr fullclone() const;
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<ShowerParticle> initShowerParticle;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
ShowerParticle & operator=(const ShowerParticle &);
private:
/**
* Whether the particle is in the final or initial state
*/
bool _isFinalState;
/**
* Whether the particle came from
*/
unsigned int _perturbative;
/**
* Does a particle produced in the backward shower initiate a time-like shower
*/
bool _initiatesTLS;
/**
* Dimensionless parameters
*/
vector<double> _parameters;
/**
* Dimensionful parameters
*/
vector<Energy> _variables;
/**
* The beam energy fraction for particle's in the initial state
*/
double _x;
/**
* The shower kinematics for the particle
*/
ShoKinPtr _showerKinematics;
/**
* Storage of the evolution scales
*/
- map<ShowerPartnerType::Type ,pair<Energy,Energy> > evolutionScales_;
+ EvolutionScales scales_;
/**
* Virtual mass
*/
Energy _vMass;
/**
* Partners
*/
tShowerParticlePtr _partner;
/**
* Pointer to ThePEG Particle this ShowerParticle was created from
*/
const tcPPtr _thePEGBase;
/**
* Progenitor
*/
tShowerParticlePtr _progenitor;
/**
* Partners
*/
vector<EvolutionPartner> partners_;
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of ShowerParticle. */
template <>
struct BaseClassTrait<Herwig::ShowerParticle,1> {
/** Typedef of the first base class of ShowerParticle. */
typedef Particle NthBase;
};
/** This template specialization informs ThePEG about the name of
* the ShowerParticle class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::ShowerParticle>
: public ClassTraitsBase<Herwig::ShowerParticle> {
/** Return a platform-independent class name */
static string className() { return "Herwig::ShowerParticle"; }
/** Create a Event object. */
static TPtr create() { return TPtr::Create(Herwig::ShowerParticle(tcEventPDPtr(),true)); }
};
/** @endcond */
}
#endif /* HERWIG_ShowerParticle_H */
diff --git a/Shower/Default/Decay_QTildeShowerKinematics1to2.cc b/Shower/Default/Decay_QTildeShowerKinematics1to2.cc
--- a/Shower/Default/Decay_QTildeShowerKinematics1to2.cc
+++ b/Shower/Default/Decay_QTildeShowerKinematics1to2.cc
@@ -1,143 +1,109 @@
// -*- C++ -*-
//
// Decay_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 Decay_QTildeShowerKinematics1to2 class.
//
#include "Decay_QTildeShowerKinematics1to2.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "Herwig++/Shower/SplittingFunctions/SplittingFunction.h"
#include "Herwig++/Shower/Base/ShowerParticle.h"
#include <cassert>
using namespace Herwig;
void Decay_QTildeShowerKinematics1to2::
updateChildren(const tShowerParticlePtr parent,
const ShowerParticleVector & children,
ShowerPartnerType::Type partnerType,
bool angularOrder ) const {
assert(children.size() == 2);
// set the values
if(parent->showerVariables().empty()) {
parent->showerVariables().resize(3);
parent->showerParameters().resize(2);
parent->showerParameters()[0]=1.;
}
for(unsigned int ix=0;ix<2;++ix) {
children[ix]->showerVariables() .resize(3);
children[ix]->showerParameters().resize(2);
}
- // angular-ordered scale for 2nd child
- Energy AOScale = angularOrder ? (1.-z())*scale() : scale();
- // QED
- if(partnerType==ShowerPartnerType::QED) {
- for(map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
- it = parent->evolutionScales().begin();
- it!=parent->evolutionScales().end();++it) {
- children[0]->evolutionScale(it->first,make_pair(min(scale(),it->second.first ),
- min(scale(),it->second.second)));
- if(it->first==ShowerPartnerType::QED) {
- children[1]->evolutionScale(it->first,make_pair(min(AOScale,it->second.first ),
- min(scale(),it->second.second)));
- }
- }
- }
- // QCD
- else {
- // scales for the emittor
- for(map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
- it = parent->evolutionScales().begin();
- it!=parent->evolutionScales().end();++it) {
- children[0]->evolutionScale(it->first,make_pair(min(scale(),it->second.first ),
- min(scale(),it->second.second)));
- }
- PDT::Colour emittedColour = children[1]->dataPtr()->iColour();
- if(emittedColour==PDT::Colour3||emittedColour==PDT::Colour8||emittedColour==PDT::Colour6) {
- children[1]->evolutionScale(ShowerPartnerType::QCDColourLine,
- make_pair(AOScale,scale()));
- }
- if(emittedColour==PDT::Colour3bar||emittedColour==PDT::Colour8||emittedColour==PDT::Colour6bar) {
- children[1]->evolutionScale(ShowerPartnerType::QCDAntiColourLine,
- make_pair(AOScale,scale()));
- }
- if(children[1]->dataPtr()->charged()) {
- children[1]->evolutionScale(ShowerPartnerType::QED,make_pair(AOScale,scale()));
- }
- }
+ // calculate the scales
+ splittingFn()->evaluateDecayScales(partnerType,scale(),z(),parent,
+ children[0],children[1]);
// determine alphas of children according to interpretation of z
children[0]->showerParameters()[0]= z() *parent->showerParameters()[0];
children[1]->showerParameters()[0]=(1.-z())*parent->showerParameters()[0];
children[0]->showerVariables()[0]= pT()*cos(phi()) +
z() *parent->showerVariables()[0];
children[0]->showerVariables()[1]= pT()*sin(phi()) +
z() *parent->showerVariables()[1];
children[1]->showerVariables()[0]= - pT()*cos(phi()) +
(1.-z())*parent->showerVariables()[0];
children[1]->showerVariables()[1]= - pT()*sin(phi()) +
(1.-z())*parent->showerVariables()[1];
for(unsigned int ix=0;ix<2;++ix)
children[ix]->showerVariables()[2]=
sqrt(sqr(children[ix]->showerVariables()[0])+
sqr(children[ix]->showerVariables()[1]));
// 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]);
}
void Decay_QTildeShowerKinematics1to2::
reconstructParent( const tShowerParticlePtr, const ParticleVector &) const {
throw Exception() << "Decay_QTildeShowerKinematics1to2::updateParent not implemented"
<< Exception::abortnow;
}
void Decay_QTildeShowerKinematics1to2::
reconstructLast(const tShowerParticlePtr theLast,
unsigned int iopt,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();
theLast->showerParameters()[1]=
(sqr(theMass) + sqr(theLast->showerVariables()[2])
- sqr( theLast->showerParameters()[0] )*pVector().m2())
/ ( 2.*theLast->showerParameters()[0]*p_dot_n() );
// set that new momentum
theLast->set5Momentum( sudakov2Momentum( theLast->showerParameters()[0],
theLast->showerParameters()[1],
theLast->showerVariables()[0],
theLast->showerVariables()[1],iopt));
}
void Decay_QTildeShowerKinematics1to2::initialize(ShowerParticle & particle,PPtr) {
Lorentz5Momentum p, n, ppartner, pcm;
assert(particle.perturbative()!=1);
// this is for the initial decaying particle
if(particle.perturbative()==2) {
p = particle.momentum();
ShowerParticlePtr partner=particle.partner();
Lorentz5Momentum ppartner(partner->momentum());
// reomved to make inverse recon work properly
//if(partner->thePEGBase()) ppartner=partner->thePEGBase()->momentum();
pcm=ppartner;
Boost boost(p.findBoostToCM());
pcm.boost(boost);
n = Lorentz5Momentum( ZERO,0.5*p.mass()*pcm.vect().unit());
n.boost( -boost);
}
else {
tShoKinPtr kin=dynamic_ptr_cast<ShowerParticlePtr>(particle.parents()[0])
->showerKinematics();
p = kin->getBasis()[0];
n = kin->getBasis()[1];
}
setBasis(p,n);
}
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,314 +1,201 @@
// -*- 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"
using namespace Herwig;
void FS_QTildeShowerKinematics1to2::
updateChildren(const tShowerParticlePtr parent,
const ShowerParticleVector & children,
ShowerPartnerType::Type partnerType,
bool angularOrder) const {
assert(children.size()==2);
// resize the parameter vectors
if(parent->showerVariables().empty()) {
parent->showerVariables().resize(3);
parent->showerParameters().resize(2);
parent->showerParameters()[0]=1.;
}
children[0]->showerVariables() .resize(3);
children[0]->showerParameters().resize(2);
children[1]->showerVariables() .resize(3);
children[1]->showerParameters().resize(2);
- // identify emittor and spectator
- tShowerParticlePtr emittor = children[0];
- tShowerParticlePtr emitted = children[1];
- double zChild[2] = {z(),1.-z()};
- bool bosonSplitting(false);
- // special for g -> gg, particle highest z is emittor
- if(emittor->id()==emitted->id()&&emittor->id()==parent->id()) {
- if(zChild[1]>zChild[0]) {
- swap(zChild[0],zChild[1]);
- swap(emitted,emittor);
- }
- }
- // otherwise if particle ID same
- else if(emitted->id()==parent->id()) {
- swap(zChild[0],zChild[1]);
- swap(emitted,emittor);
- }
- // no real emittor/eemitted
- else if(emittor->id()!=parent->id()) {
- bosonSplitting = true;
- }
- // scales for angular ordering
- Energy AOScale[2];
- if(angularOrder) {
- for(unsigned int ix=0;ix<2;++ix) AOScale [ix] = zChild[ix]*scale();
- }
- else {
- for(unsigned int ix=0;ix<2;++ix) AOScale [ix] = scale();
- }
- // if not these cases doesn't matter
- // now the various scales
- if(partnerType==ShowerPartnerType::QED) {
- // normal case
- if(!bosonSplitting) {
- for(map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
- it = parent->evolutionScales().begin();
- it!=parent->evolutionScales().end();++it) {
- emittor->evolutionScale(it->first,make_pair(min(AOScale[0],it->second.first ),
- min(scale() ,it->second.second)));
- if(it->first==ShowerPartnerType::QED) {
- emitted->evolutionScale(it->first,make_pair(min(AOScale[1],it->second.first ),
- min(scale() ,it->second.second)));
- }
- }
- }
- // gamma -> f fbar
- else {
- // set QED scales as normal
- for(map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
- it = parent->evolutionScales().begin();
- it!=parent->evolutionScales().end();++it) {
- assert(it->first==ShowerPartnerType::QED);
- emittor->evolutionScale(it->first,make_pair(AOScale[0],scale()));
- emitted->evolutionScale(it->first,make_pair(AOScale[1],scale()));
- }
- // and any QCD scales needed
- PDT::Colour emittorColour = emittor->dataPtr()->iColour();
- if(emittorColour==PDT::Colour3||emittorColour==PDT::Colour8||emittorColour==PDT::Colour6) {
- emittor->evolutionScale(ShowerPartnerType::QCDColourLine,
- make_pair(AOScale[0],scale()));
- }
- if(emittorColour==PDT::Colour3bar||emittorColour==PDT::Colour8||emittorColour==PDT::Colour6bar) {
- emittor->evolutionScale(ShowerPartnerType::QCDAntiColourLine,
- make_pair(AOScale[0],scale()));
- }
- PDT::Colour emittedColour = emitted->dataPtr()->iColour();
- if(emittedColour==PDT::Colour3||emittedColour==PDT::Colour8||emittedColour==PDT::Colour6) {
- emitted->evolutionScale(ShowerPartnerType::QCDColourLine,
- make_pair(AOScale[1],scale()));
- }
- if(emittedColour==PDT::Colour3bar||emittedColour==PDT::Colour8||emittedColour==PDT::Colour6bar) {
- emitted->evolutionScale(ShowerPartnerType::QCDAntiColourLine,
- make_pair(AOScale[1],scale()));
- }
- }
- }
- else {
- // normal case
- if(!bosonSplitting) {
- // scales for the emittor
- for(map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
- it = parent->evolutionScales().begin();
- it!=parent->evolutionScales().end();++it) {
- emittor->evolutionScale(it->first,make_pair(min(AOScale[0],it->second.first ),
- min(scale() ,it->second.second)));
- }
- }
- else {
- // QCD scales needed
- PDT::Colour emittorColour = emittor->dataPtr()->iColour();
- if(emittorColour==PDT::Colour3||emittorColour==PDT::Colour8||emittorColour==PDT::Colour6) {
- emittor->evolutionScale(ShowerPartnerType::QCDColourLine,
- make_pair(AOScale[0],scale()));
- }
- if(emittorColour==PDT::Colour3bar||emittorColour==PDT::Colour8||emittorColour==PDT::Colour6bar) {
- emittor->evolutionScale(ShowerPartnerType::QCDAntiColourLine,
- make_pair(AOScale[0],scale()));
- }
- // QED scales
- if(emittor->dataPtr()->charged()) {
- emittor->evolutionScale(ShowerPartnerType::QED,make_pair(AOScale[0],scale()));
- }
- }
- PDT::Colour emittedColour = emitted->dataPtr()->iColour();
- if(emittedColour==PDT::Colour3||emittedColour==PDT::Colour8||emittedColour==PDT::Colour6) {
- emitted->evolutionScale(ShowerPartnerType::QCDColourLine,
- make_pair(AOScale[1],scale()));
- }
- if(emittedColour==PDT::Colour3bar||emittedColour==PDT::Colour8||emittedColour==PDT::Colour6bar) {
- emitted->evolutionScale(ShowerPartnerType::QCDAntiColourLine,
- make_pair(AOScale[1],scale()));
- }
- if(emitted->dataPtr()->charged()) {
- emitted->evolutionScale(ShowerPartnerType::QED,make_pair(AOScale[1],scale()));
- }
- }
+ // calculate the scales
+ splittingFn()->evaluateFinalStateScales(partnerType,scale(),z(),parent,
+ children[0],children[1]);
// debugging printout if needed
if(Debug::level >= 10 ) printScales(parent,children[0],children[1]);
// determine alphas of children according to interpretation of z
children[0]->showerParameters()[0]= z() *parent->showerParameters()[0];
children[1]->showerParameters()[0]= (1.-z())*parent->showerParameters()[0];
// set the values
children[0]->showerVariables()[0]=
pT()*cos(phi()) + z() *parent->showerVariables()[0];
children[0]->showerVariables()[1]=
pT()*sin(phi()) + z() *parent->showerVariables()[1];
children[1]->showerVariables()[0]=
- pT()*cos(phi()) + (1.-z())*parent->showerVariables()[0];
children[1]->showerVariables()[1]=
- pT()*sin(phi()) + (1.-z())*parent->showerVariables()[1];
for(unsigned int ix=0;ix<2;++ix)
children[ix]->showerVariables()[2]=
sqrt(sqr(children[ix]->showerVariables()[0])+
sqr(children[ix]->showerVariables()[1]));
// 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]);
}
void FS_QTildeShowerKinematics1to2::
reconstructParent(const tShowerParticlePtr parent,
const ParticleVector & children ) const {
if(children.size() != 2)
throw Exception() << "FS_QTildeShowerKinematics1to2::updateParent() "
<< "Warning! too many children!"
<< Exception::eventerror;
ShowerParticlePtr c1 = dynamic_ptr_cast<ShowerParticlePtr>(children[0]);
ShowerParticlePtr c2 = dynamic_ptr_cast<ShowerParticlePtr>(children[1]);
parent->showerParameters()[1]=
c1->showerParameters()[1] + c2->showerParameters()[1];
parent->set5Momentum( c1->momentum() + c2->momentum() );
}
void FS_QTildeShowerKinematics1to2::reconstructLast(const tShowerParticlePtr theLast,
unsigned int iopt,
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();
theLast->showerParameters()[1]=
(sqr(theMass) + sqr(theLast->showerVariables()[2])
- sqr( theLast->showerParameters()[0] )*pVector().m2())
/ ( 2.*theLast->showerParameters()[0]*p_dot_n() );
// set that new momentum
theLast->set5Momentum(sudakov2Momentum( theLast->showerParameters()[0],
theLast->showerParameters()[1],
theLast->showerVariables()[0],
theLast->showerVariables()[1],iopt));
}
void FS_QTildeShowerKinematics1to2::initialize(ShowerParticle & particle,PPtr) {
// set the basis vectors
Lorentz5Momentum p,n;
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);
}
}
}
else if(particle.initiatesTLS()) {
tShoKinPtr kin=dynamic_ptr_cast<ShowerParticlePtr>
(particle.parents()[0]->children()[0])->showerKinematics();
p = kin->getBasis()[0];
n = kin->getBasis()[1];
}
else {
tShoKinPtr kin=dynamic_ptr_cast<ShowerParticlePtr>(particle.parents()[0])
->showerKinematics();
p = kin->getBasis()[0];
n = kin->getBasis()[1];
}
// set the basis vectors
setBasis(p,n);
}
void FS_QTildeShowerKinematics1to2::updateParent(const tShowerParticlePtr parent,
const ShowerParticleVector & children,
ShowerPartnerType::Type,
bool) const {
IdList ids(3);
ids[0] = parent->id();
ids[1] = children[0]->id();
ids[2] = children[1]->id();
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 {
// set the values
children[0]->showerVariables()[0]=
pT()*cos(phi()) + z() *parent->showerVariables()[0];
children[0]->showerVariables()[1]=
pT()*sin(phi()) + z() *parent->showerVariables()[1];
children[1]->showerVariables()[0]=
- pT()*cos(phi()) + (1.-z())*parent->showerVariables()[0];
children[1]->showerVariables()[1]=
- pT()*sin(phi()) + (1.-z())*parent->showerVariables()[1];
for(unsigned int ix=0;ix<2;++ix)
children[ix]->showerVariables()[2]=
sqrt(sqr(children[ix]->showerVariables()[0])+
sqr(children[ix]->showerVariables()[1]));
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,264 +1,168 @@
// -*- 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 <cassert>
using namespace Herwig;
void IS_QTildeShowerKinematics1to2::
updateChildren( const tShowerParticlePtr theParent,
const ShowerParticleVector & theChildren,
ShowerPartnerType::Type,
bool ) const {
theChildren[1]->showerVariables().resize(3);
theChildren[1]->showerParameters().resize(2);
theChildren[1]->showerParameters()[0]=
(1.-z())*theParent->showerParameters()[0];
double cphi = cos(phi()),sphi = sin(phi());
theChildren[1]->showerVariables()[0] =
(1.-z())*theParent->showerVariables()[0] - cphi*pT();
theChildren[1]->showerVariables()[1]=
(1.-z())*theParent->showerVariables()[1] - sphi*pT();
theChildren[1]->showerVariables()[2]=
sqrt(sqr(theChildren[1]->showerVariables()[0])+
sqr(theChildren[1]->showerVariables()[1]));
// space-like child
theChildren[0]->showerParameters()[0] =
theParent->showerParameters()[0] - theChildren[1]->showerParameters()[0];
theChildren[0]->showerParameters()[1] =
theParent->showerParameters()[1] - theChildren[1]->showerParameters()[1];
theChildren[0]->showerVariables()[0]=
theParent->showerVariables()[0] - theChildren[1]->showerVariables()[0];
theChildren[0]->showerVariables()[1]=
theParent->showerVariables()[1] - theChildren[1]->showerVariables()[1];
}
void IS_QTildeShowerKinematics1to2::
updateParent(const tShowerParticlePtr parent,
const ShowerParticleVector & children,
ShowerPartnerType::Type partnerType,
bool angularOrder) const {
- // scale for time-like child
- Energy AOScale = (angularOrder ? (1.-z()) : 1. )*scale();
- // scales for parent if same as spacelike child
- if(parent->id()==children[0]->id()) {
- for(map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
- it = children[0]->evolutionScales().begin();
- it!=children[0]->evolutionScales().end();++it) {
- parent->evolutionScale(it->first,make_pair(min(scale(),it->second.first ),
- min(scale(),it->second.second)));
- }
- if(partnerType==ShowerPartnerType::QED) {
- children[1]->evolutionScale(partnerType,make_pair(AOScale,scale()));
- }
- PDT::Colour childColour = children[1]->dataPtr()->iColour();
- if(childColour==PDT::Colour3||childColour==PDT::Colour8||childColour==PDT::Colour6) {
- children[1]->evolutionScale(ShowerPartnerType::QCDColourLine,make_pair(AOScale,scale()));
- }
- if(childColour==PDT::Colour3bar||childColour==PDT::Colour8||childColour==PDT::Colour6bar) {
- children[1]->evolutionScale(ShowerPartnerType::QCDAntiColourLine,make_pair(AOScale,scale()));
- }
- }
- // scales if parent same as timelike child
- else if(parent->id()==children[1]->id()) {
- if(parent->dataPtr()->charged()) {
- parent ->evolutionScale(ShowerPartnerType::QED,make_pair(scale(),scale()));
- children[1]->evolutionScale(ShowerPartnerType::QED,make_pair(AOScale,scale()));
- }
- PDT::Colour childColour = children[1]->dataPtr()->iColour();
- if(partnerType==ShowerPartnerType::QED) {
- if(childColour==PDT::Colour3||childColour==PDT::Colour8||childColour==PDT::Colour6) {
- parent ->evolutionScale(ShowerPartnerType::QCDColourLine,make_pair(scale(),scale()));
- children[1]->evolutionScale(ShowerPartnerType::QCDColourLine,make_pair(AOScale,scale()));
- }
- if(childColour==PDT::Colour3bar||childColour==PDT::Colour8||childColour==PDT::Colour6bar) {
- parent ->evolutionScale(ShowerPartnerType::QCDAntiColourLine,make_pair(scale(),scale()));
- children[1]->evolutionScale(ShowerPartnerType::QCDAntiColourLine,make_pair(AOScale,scale()));
- }
- }
- else {
- if(childColour==PDT::Colour3||childColour==PDT::Colour8||childColour==PDT::Colour6) {
- pair<Energy,Energy> pScale = children[0]->evolutionScale(ShowerPartnerType::QCDColourLine);
- parent ->evolutionScale(ShowerPartnerType::QCDColourLine,
- make_pair(min(scale(),pScale.first ),
- min(scale(),pScale.second)));
- children[1]->evolutionScale(ShowerPartnerType::QCDColourLine,make_pair(AOScale,scale()));
- }
- if(childColour==PDT::Colour3bar||childColour==PDT::Colour8||childColour==PDT::Colour6bar) {
- pair<Energy,Energy> pScale = children[0]->evolutionScale(ShowerPartnerType::QCDAntiColourLine);
- parent ->evolutionScale(ShowerPartnerType::QCDAntiColourLine,
- make_pair(min(scale(),pScale.first ),
- min(scale(),pScale.second)));
- children[1]->evolutionScale(ShowerPartnerType::QCDAntiColourLine,make_pair(AOScale,scale()));
- }
- }
- }
- // g -> q qbar, or gamma -> e+e-
- else if(children[0]->id()==-children[1]->id()) {
- PDT::Colour childColour = children[1]->dataPtr()->iColour();
- if(partnerType==ShowerPartnerType::QED) {
- parent ->evolutionScale(ShowerPartnerType::QED,make_pair(scale(),scale()));
- children[1]->evolutionScale(ShowerPartnerType::QED,make_pair(AOScale,scale()));
- if(childColour==PDT::Colour3||childColour==PDT::Colour8||childColour==PDT::Colour6) {
- pair<Energy,Energy> pScale = children[0]->evolutionScale(ShowerPartnerType::QCDAntiColourLine);
- children[1]->evolutionScale(ShowerPartnerType::QCDColourLine,
- make_pair(min(AOScale,pScale.first ),
- min(AOScale,pScale.second)));
- }
- if(childColour==PDT::Colour3bar||childColour==PDT::Colour8||childColour==PDT::Colour6bar) {
- pair<Energy,Energy> pScale = children[0]->evolutionScale(ShowerPartnerType::QCDColourLine );
- children[1]->evolutionScale(ShowerPartnerType::QCDAntiColourLine,
- make_pair(min(AOScale,pScale.first ),
- min(AOScale,pScale.second)));
- }
- }
- else {
- pair<Energy,Energy> pScale = children[0]->evolutionScale(ShowerPartnerType::QED);
- children[1]->evolutionScale(ShowerPartnerType::QED,
- make_pair(min(AOScale,pScale.first ),
- min(AOScale,pScale.second)));
- if(childColour==PDT::Colour3||childColour==PDT::Colour8||childColour==PDT::Colour6) {
- pScale = children[0]->evolutionScale(ShowerPartnerType::QCDAntiColourLine);
- children[1]->evolutionScale(ShowerPartnerType::QCDColourLine,
- make_pair(min(AOScale,pScale.first ),
- min(AOScale,pScale.second)));
- }
- if(childColour==PDT::Colour3bar||childColour==PDT::Colour8||childColour==PDT::Colour6bar) {
- pScale = children[0]->evolutionScale(ShowerPartnerType::QCDColourLine );
- children[1]->evolutionScale(ShowerPartnerType::QCDAntiColourLine,
- make_pair(min(AOScale,pScale.first ),
- min(AOScale,pScale.second)));
- }
- parent->evolutionScale(ShowerPartnerType::QCDColourLine,
- make_pair(min(scale(),pScale.first ),
- min(scale(),pScale.second)));
- parent->evolutionScale(ShowerPartnerType::QCDAntiColourLine,
- make_pair(min(scale(),pScale.first ),
- min(scale(),pScale.second)));
- }
- }
+ // calculate the scales
+ splittingFn()->evaluateInitialStateScales(partnerType,scale(),z(),parent,
+ children[0],children[1]);
// debugging printout if needed
if(Debug::level >= 10 ) printScales(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());
// create the storage for the shower variables
parent->showerVariables().resize(3);
parent->showerParameters().resize(2);
if(children[0]->showerVariables().empty()) {
children[0]->showerVariables().resize(3);
children[0]->showerParameters().resize(2);
}
}
void IS_QTildeShowerKinematics1to2::
reconstructParent(const tShowerParticlePtr theParent,
const ParticleVector & theChildren ) const {
ShowerParticlePtr c1 = dynamic_ptr_cast<ShowerParticlePtr>(theChildren[0]);
ShowerParticlePtr c2 = dynamic_ptr_cast<ShowerParticlePtr>(theChildren[1]);
// 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
c2->showerParameters()[1]=
0.5*(sqr(c2->data().constituentMass())+sqr(c2->showerVariables()[2]))
/(c2->showerParameters()[0]*p_dot_n());
c2->set5Momentum(sudakov2Momentum(c2->showerParameters()[0],
c2->showerParameters()[1],
c2->showerVariables()[0],
c2->showerVariables()[1],0));
// 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;
Energy2 pt2=sqr(px)+sqr(py);
theLast->showerParameters()[0]=theLast->x();
theLast->showerParameters()[1]=0.5*pt2/
theLast->showerParameters()[0]/p_dot_n();
theLast->showerVariables().resize(3);
theLast->showerParameters().resize(2);
for(unsigned int ix=0;ix<3;++ix) theLast->showerVariables()[ix]=ZERO;
// momentum
Lorentz5Momentum ntemp=Lorentz5Momentum(ZERO,-pVector().vect());
double beta = 0.5*pt2/
theLast->showerParameters()[0]/(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);
if(particle.perturbative()==1) {
// find the partner and its momentum
ShowerParticlePtr partner=particle.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());
}
}
else {
p = dynamic_ptr_cast<ShowerParticlePtr>(particle.children()[0])
->showerKinematics()->getBasis()[0];
n = dynamic_ptr_cast<ShowerParticlePtr>(particle.children()[0])
->showerKinematics()->getBasis()[1];
}
setBasis(p,n);
}
diff --git a/Shower/Default/QTildeShowerKinematics1to2.cc b/Shower/Default/QTildeShowerKinematics1to2.cc
--- a/Shower/Default/QTildeShowerKinematics1to2.cc
+++ b/Shower/Default/QTildeShowerKinematics1to2.cc
@@ -1,132 +1,133 @@
// -*- C++ -*-
//
// 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 QTildeShowerKinematics1to2 class.
//
#include "QTildeShowerKinematics1to2.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "Herwig++/Shower/Base/ShowerParticle.h"
#include "ThePEG/Repository/CurrentGenerator.h"
using namespace Herwig;
vector<Lorentz5Momentum> QTildeShowerKinematics1to2::getBasis() const {
vector<Lorentz5Momentum> dum;
dum.push_back( _pVector );
dum.push_back( _nVector );
return dum;
}
void QTildeShowerKinematics1to2::setBasis(const Lorentz5Momentum &p,
const Lorentz5Momentum & n) {
_pVector=p;
_nVector=n;
}
Lorentz5Momentum QTildeShowerKinematics1to2::
sudakov2Momentum(double alpha, double beta, Energy px,
Energy py,unsigned int iopt) const {
if(isnan(beta)||isinf(beta))
throw Exception() << "beta infinite in "
<< "QTildeShowerKinematics1to2::sudakov2Momentum()"
<< Exception::eventerror;
Lorentz5Momentum dq;
if(iopt==0) {
const Boost beta_bb = -(_pVector + _nVector).boostVector();
Lorentz5Momentum p_bb = _pVector;
Lorentz5Momentum n_bb = _nVector;
p_bb.boost( beta_bb );
n_bb.boost( beta_bb );
// set first in b2b frame along z-axis (assuming that p and n are
// b2b as checked above)
dq=Lorentz5Momentum(ZERO, ZERO, (alpha - beta)*p_bb.vect().mag(),
alpha*p_bb.t() + beta*n_bb.t());
// add transverse components
dq.setX(px);
dq.setY(py);
// rotate to have z-axis parallel to p
// this rotation changed by PR to a different rotation with the same effect
// but different azimuthal angle to make implementing spin correlations easier
// dq.rotateUz( unitVector(p_bb.vect()) );
Axis axis(p_bb.vect().unit());
if(axis.perp2()>0.) {
LorentzRotation rot;
double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
rot.setRotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
dq.transform(rot);
}
else if(axis.z()<0.) {
dq.setZ(-dq.z());
}
// boost back
dq.boost( -beta_bb );
dq.rescaleMass();
// return the momentum
}
else {
const Boost beta_bb = -pVector().boostVector();
Lorentz5Momentum p_bb = pVector();
Lorentz5Momentum n_bb = nVector();
p_bb.boost( beta_bb );
n_bb.boost( beta_bb );
// set first in b2b frame along z-axis (assuming that p and n are
// b2b as checked above)
dq=Lorentz5Momentum (ZERO, ZERO, 0.5*beta*pVector().mass(),
alpha*pVector().mass() + 0.5*beta*pVector().mass());
// add transverse components
dq.setX(px);
dq.setY(py);
// changed to be same as other case
// dq.rotateUz( unitVector(n_bb.vect()) );
Axis axis(n_bb.vect().unit());
if(axis.perp2()>0.) {
LorentzRotation rot;
double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
rot.setRotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
dq.transform(rot);
}
else if(axis.z()<0.) {
dq.setZ(-dq.z());
}
// boost back
dq.boost( -beta_bb );
dq.rescaleMass();
}
return dq;
}
void QTildeShowerKinematics1to2::printScales(tShowerParticlePtr parent,
tShowerParticlePtr child1,
tShowerParticlePtr child2) const {
- CurrentGenerator::log() << *parent << "\n" << *child1 << " "
- << *child2 << "\n";
- CurrentGenerator::log() << "testing parent\n";
- for(map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
- it = parent->evolutionScales().begin();
- it!=parent->evolutionScales().end();++it) {
- CurrentGenerator::log() << it->first << " " << it->second.first/GeV << " "
- << it->second.second/GeV << "\n";
- }
- CurrentGenerator::log() << "testing child[0]\n";
- for(map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
- it = child1->evolutionScales().begin();
- it!=child1->evolutionScales().end();++it) {
- CurrentGenerator::log() << it->first << " " << it->second.first/GeV << " "
- << it->second.second/GeV << "\n";
- }
- CurrentGenerator::log() << "testing child[1]\n";
- for(map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
- it = child2->evolutionScales().begin();
- it!=child2->evolutionScales().end();++it) {
- CurrentGenerator::log() << it->first << " " << it->second.first/GeV << " "
- << it->second.second/GeV << "\n";
- }
+ assert(false);
+ // CurrentGenerator::log() << *parent << "\n" << *child1 << " "
+ // << *child2 << "\n";
+ // CurrentGenerator::log() << "testing parent\n";
+ // for(map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
+ // it = parent->evolutionScales().begin();
+ // it!=parent->evolutionScales().end();++it) {
+ // CurrentGenerator::log() << it->first << " " << it->second.first/GeV << " "
+ // << it->second.second/GeV << "\n";
+ // }
+ // CurrentGenerator::log() << "testing child[0]\n";
+ // for(map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
+ // it = child1->evolutionScales().begin();
+ // it!=child1->evolutionScales().end();++it) {
+ // CurrentGenerator::log() << it->first << " " << it->second.first/GeV << " "
+ // << it->second.second/GeV << "\n";
+ // }
+ // CurrentGenerator::log() << "testing child[1]\n";
+ // for(map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
+ // it = child2->evolutionScales().begin();
+ // it!=child2->evolutionScales().end();++it) {
+ // CurrentGenerator::log() << it->first << " " << it->second.first/GeV << " "
+ // << it->second.second/GeV << "\n";
+ // }
}
diff --git a/Shower/SplittingFunctions/SplittingFunction.cc b/Shower/SplittingFunctions/SplittingFunction.cc
--- a/Shower/SplittingFunctions/SplittingFunction.cc
+++ b/Shower/SplittingFunctions/SplittingFunction.cc
@@ -1,608 +1,884 @@
// -*- C++ -*-
//
// SplittingFunction.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 SplittingFunction class.
//
#include "SplittingFunction.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Utilities/EnumIO.h"
#include "Herwig++/Shower/Base/ShowerParticle.h"
using namespace Herwig;
// Static variable needed for the type description system in ThePEG.
DescribeAbstractClass<SplittingFunction,Interfaced>
describeThePEGSplittingFunction("Herwig::SplittingFunction", "Herwig.so");
void SplittingFunction::Init() {
static ClassDocumentation<SplittingFunction> documentation
("The SplittingFunction class is the based class for 1->2 splitting functions"
" in Herwig++");
static Switch<SplittingFunction,ColourStructure> interfaceColourStructure
("ColourStructure",
"The colour structure for the splitting function",
&SplittingFunction::_colourStructure, Undefined, false, false);
static SwitchOption interfaceColourStructureTripletTripletOctet
(interfaceColourStructure,
"TripletTripletOctet",
"3 -> 3 8",
TripletTripletOctet);
static SwitchOption interfaceColourStructureOctetOctetOctet
(interfaceColourStructure,
"OctetOctetOctet",
"8 -> 8 8",
OctetOctetOctet);
static SwitchOption interfaceColourStructureOctetTripletTriplet
(interfaceColourStructure,
"OctetTripletTriplet",
"8 -> 3 3bar",
OctetTripletTriplet);
static SwitchOption interfaceColourStructureTripletOctetTriplet
(interfaceColourStructure,
"TripletOctetTriplet",
"3 -> 8 3",
TripletOctetTriplet);
static SwitchOption interfaceColourStructureSextetSextetOctet
(interfaceColourStructure,
"SextetSextetOctet",
"6 -> 6 8",
SextetSextetOctet);
static SwitchOption interfaceColourStructureChargedChargedNeutral
(interfaceColourStructure,
"ChargedChargedNeutral",
"q -> q 0",
ChargedChargedNeutral);
static SwitchOption interfaceColourStructureNeutralChargedCharged
(interfaceColourStructure,
"NeutralChargedCharged",
"0 -> q qbar",
NeutralChargedCharged);
static SwitchOption interfaceColourStructureChargedNeutralCharged
(interfaceColourStructure,
"ChargedNeutralCharged",
"q -> 0 q",
ChargedNeutralCharged);
static Switch<SplittingFunction,ShowerInteraction::Type>
interfaceInteractionType
("InteractionType",
"Type of the interaction",
&SplittingFunction::_interactionType,
ShowerInteraction::UNDEFINED, false, false);
static SwitchOption interfaceInteractionTypeQCD
(interfaceInteractionType,
"QCD","QCD",ShowerInteraction::QCD);
static SwitchOption interfaceInteractionTypeQED
(interfaceInteractionType,
"QED","QED",ShowerInteraction::QED);
static Switch<SplittingFunction,int> interfaceSplittingColourMethod
("SplittingColourMethod",
"Choice of assigning colour in 8->88 splittings.",
&SplittingFunction::_splittingColourMethod, 0, false, false);
static SwitchOption interfaceSplittingColourMethodRandom
(interfaceSplittingColourMethod,
"Random",
"Choose colour assignments randomly.",
0);
static SwitchOption interfaceSplittingColourMethodCorrectLines
(interfaceSplittingColourMethod,
"CorrectLines",
"Choose correct lines for colour.",
1);
static SwitchOption interfaceSplittingColourMethodRandomRecord
(interfaceSplittingColourMethod,
"RandomRecord",
"Choose colour assignments randomly and record the result.",
2);
static Switch<SplittingFunction,bool> interfaceAngularOrdered
("AngularOrdered",
"Whether or not this interaction is angular ordered, "
"normally only g->q qbar and gamma-> f fbar are the only ones which aren't.",
&SplittingFunction::angularOrdered_, true, false, false);
static SwitchOption interfaceAngularOrderedYes
(interfaceAngularOrdered,
"Yes",
"Interaction is angular ordered",
true);
static SwitchOption interfaceAngularOrderedNo
(interfaceAngularOrdered,
"No",
"Interaction isn't angular ordered",
false);
}
void SplittingFunction::persistentOutput(PersistentOStream & os) const {
using namespace ShowerInteraction;
os << oenum(_interactionType) << _interactionOrder
<< oenum(_colourStructure) << _colourFactor
<< angularOrdered_ << _splittingColourMethod;
}
void SplittingFunction::persistentInput(PersistentIStream & is, int) {
using namespace ShowerInteraction;
is >> ienum(_interactionType) >> _interactionOrder
>> ienum(_colourStructure) >> _colourFactor
>> angularOrdered_ >> _splittingColourMethod;
}
void SplittingFunction::colourConnection(tShowerParticlePtr parent,
tShowerParticlePtr first,
tShowerParticlePtr second,
ShowerPartnerType::Type partnerType,
const bool back) const {
if(_colourStructure==TripletTripletOctet) {
if(!back) {
ColinePair cparent = ColinePair(parent->colourLine(),
parent->antiColourLine());
// ensure input consistency
assert(( cparent.first && !cparent.second &&
partnerType==ShowerPartnerType::QCDColourLine) ||
( !cparent.first && cparent.second &&
partnerType==ShowerPartnerType::QCDAntiColourLine));
// q -> q g
if(cparent.first) {
ColinePtr newline=new_ptr(ColourLine());
cparent.first->addColoured(second);
newline->addColoured ( first);
newline->addAntiColoured (second);
}
// qbar -> qbar g
else {
ColinePtr newline=new_ptr(ColourLine());
cparent.second->addAntiColoured(second);
newline->addColoured(second);
newline->addAntiColoured(first);
}
// Set progenitor
first->progenitor(parent->progenitor());
second->progenitor(parent->progenitor());
}
else {
ColinePair cfirst = ColinePair(first->colourLine(),
first->antiColourLine());
// ensure input consistency
assert(( cfirst.first && !cfirst.second &&
partnerType==ShowerPartnerType::QCDColourLine) ||
( !cfirst.first && cfirst.second &&
partnerType==ShowerPartnerType::QCDAntiColourLine));
// q -> q g
if(cfirst.first) {
ColinePtr newline=new_ptr(ColourLine());
cfirst.first->addAntiColoured(second);
newline->addColoured(second);
newline->addColoured(parent);
}
// qbar -> qbar g
else {
ColinePtr newline=new_ptr(ColourLine());
cfirst.second->addColoured(second);
newline->addAntiColoured(second);
newline->addAntiColoured(parent);
}
// Set progenitor
parent->progenitor(first->progenitor());
second->progenitor(first->progenitor());
}
}
else if(_colourStructure==OctetOctetOctet) {
if(!back) {
ColinePair cparent = ColinePair(parent->colourLine(),
parent->antiColourLine());
// ensure input consistency
assert(cparent.first&&cparent.second);
// ensure first gluon is hardest
if( first->id()==second->id() && parent->showerKinematics()->z()<0.5 )
swap(first,second);
// colour line radiates
if(partnerType==ShowerPartnerType::QCDColourLine) {
// The colour line is radiating
ColinePtr newline=new_ptr(ColourLine());
cparent.first->addColoured(second);
cparent.second->addAntiColoured(first);
newline->addColoured(first);
newline->addAntiColoured(second);
}
// anti colour line radiates
else if(partnerType==ShowerPartnerType::QCDAntiColourLine) {
ColinePtr newline=new_ptr(ColourLine());
cparent.first->addColoured(first);
cparent.second->addAntiColoured(second);
newline->addColoured(second);
newline->addAntiColoured(first);
}
else
assert(false);
}
else {
ColinePair cfirst = ColinePair(first->colourLine(),
first->antiColourLine());
// ensure input consistency
assert(cfirst.first&&cfirst.second);
// The colour line is radiating
if(partnerType==ShowerPartnerType::QCDColourLine) {
ColinePtr newline=new_ptr(ColourLine());
cfirst.first->addAntiColoured(second);
cfirst.second->addAntiColoured(parent);
newline->addColoured(parent);
newline->addColoured(second);
}
// anti colour line radiates
else if(partnerType==ShowerPartnerType::QCDAntiColourLine) {
ColinePtr newline=new_ptr(ColourLine());
cfirst.first->addColoured(parent);
cfirst.second->addColoured(second);
newline->addAntiColoured(second);
newline->addAntiColoured(parent);
}
else
assert(false);
}
}
else if(_colourStructure == OctetTripletTriplet) {
if(!back) {
ColinePair cparent = ColinePair(parent->colourLine(),
parent->antiColourLine());
// ensure input consistency
assert(cparent.first&&cparent.second);
cparent.first ->addColoured ( first);
cparent.second->addAntiColoured(second);
// Set progenitor
first->progenitor(parent->progenitor());
second->progenitor(parent->progenitor());
}
else {
ColinePair cfirst = ColinePair(first->colourLine(),
first->antiColourLine());
// ensure input consistency
assert(( cfirst.first && !cfirst.second) ||
(!cfirst.first && cfirst.second));
// g -> q qbar
if(cfirst.first) {
ColinePtr newline=new_ptr(ColourLine());
cfirst.first->addColoured(parent);
newline->addAntiColoured(second);
newline->addAntiColoured(parent);
}
// g -> qbar q
else {
ColinePtr newline=new_ptr(ColourLine());
cfirst.second->addAntiColoured(parent);
newline->addColoured(second);
newline->addColoured(parent);
}
// Set progenitor
parent->progenitor(first->progenitor());
second->progenitor(first->progenitor());
}
}
else if(_colourStructure == TripletOctetTriplet) {
if(!back) {
ColinePair cparent = ColinePair(parent->colourLine(),
parent->antiColourLine());
// ensure input consistency
assert(( cparent.first && !cparent.second) ||
(!cparent.first && cparent.second));
// q -> g q
if(cparent.first) {
ColinePtr newline=new_ptr(ColourLine());
cparent.first->addColoured(first);
newline->addColoured (second);
newline->addAntiColoured( first);
}
// qbar -> g qbar
else {
ColinePtr newline=new_ptr(ColourLine());
cparent.second->addAntiColoured(first);
newline->addColoured ( first);
newline->addAntiColoured(second);
}
// Set progenitor
first->progenitor(parent->progenitor());
second->progenitor(parent->progenitor());
}
else {
ColinePair cfirst = ColinePair(first->colourLine(),
first->antiColourLine());
// ensure input consistency
assert(cfirst.first&&cfirst.second);
// q -> g q
if(parent->id()>0) {
cfirst.first ->addColoured(parent);
cfirst.second->addColoured(second);
}
else {
cfirst.first ->addAntiColoured(second);
cfirst.second->addAntiColoured(parent);
}
// Set progenitor
parent->progenitor(first->progenitor());
second->progenitor(first->progenitor());
}
}
else if(_colourStructure==SextetSextetOctet) {
//make sure we're not doing backward evolution
assert(!back);
//make sure something sensible
assert(parent->colourLine() || parent->antiColourLine());
//get the colour lines or anti-colour lines
bool isAntiColour=true;
ColinePair cparent;
if(parent->colourLine()) {
cparent = ColinePair(const_ptr_cast<tColinePtr>(parent->colourInfo()->colourLines()[0]),
const_ptr_cast<tColinePtr>(parent->colourInfo()->colourLines()[1]));
isAntiColour=false;
}
else {
cparent = ColinePair(const_ptr_cast<tColinePtr>(parent->colourInfo()->antiColourLines()[0]),
const_ptr_cast<tColinePtr>(parent->colourInfo()->antiColourLines()[1]));
}
//check for sensible input
// assert(cparent.first && cparent.second);
// sextet has 2 colour lines
if(!isAntiColour) {
//pick at random which of the colour topolgies to take
double topology = UseRandom::rnd();
if(topology < 0.25) {
ColinePtr newline=new_ptr(ColourLine());
cparent.first->addColoured(second);
cparent.second->addColoured(first);
newline->addColoured(first);
newline->addAntiColoured(second);
}
else if(topology >=0.25 && topology < 0.5) {
ColinePtr newline=new_ptr(ColourLine());
cparent.first->addColoured(first);
cparent.second->addColoured(second);
newline->addColoured(first);
newline->addAntiColoured(second);
}
else if(topology >= 0.5 && topology < 0.75) {
ColinePtr newline=new_ptr(ColourLine());
cparent.first->addColoured(second);
cparent.second->addColoured(first);
newline->addColoured(first);
newline->addAntiColoured(second);
}
else {
ColinePtr newline=new_ptr(ColourLine());
cparent.first->addColoured(first);
cparent.second->addColoured(second);
newline->addColoured(first);
newline->addAntiColoured(second);
}
}
// sextet has 2 anti-colour lines
else {
double topology = UseRandom::rnd();
if(topology < 0.25){
ColinePtr newline=new_ptr(ColourLine());
cparent.first->addAntiColoured(second);
cparent.second->addAntiColoured(first);
newline->addAntiColoured(first);
newline->addColoured(second);
}
else if(topology >=0.25 && topology < 0.5) {
ColinePtr newline=new_ptr(ColourLine());
cparent.first->addAntiColoured(first);
cparent.second->addAntiColoured(second);
newline->addAntiColoured(first);
newline->addColoured(second);
}
else if(topology >= 0.5 && topology < 0.75) {
ColinePtr newline=new_ptr(ColourLine());
cparent.first->addAntiColoured(second);
cparent.second->addAntiColoured(first);
newline->addAntiColoured(first);
newline->addColoured(second);
}
else {
ColinePtr newline=new_ptr(ColourLine());
cparent.first->addAntiColoured(first);
cparent.second->addAntiColoured(second);
newline->addAntiColoured(first);
newline->addColoured(second);
}
}
}
else if(_colourStructure == ChargedChargedNeutral) {
if(!parent->data().coloured()) return;
if(!back) {
ColinePair cparent = ColinePair(parent->colourLine(),
parent->antiColourLine());
// q -> q g
if(cparent.first) {
cparent.first->addColoured(first);
}
// qbar -> qbar g
if(cparent.second) {
cparent.second->addAntiColoured(first);
}
}
else {
ColinePair cfirst = ColinePair(first->colourLine(),
first->antiColourLine());
// q -> q g
if(cfirst.first) {
cfirst.first->addColoured(parent);
}
// qbar -> qbar g
if(cfirst.second) {
cfirst.second->addAntiColoured(parent);
}
}
}
else if(_colourStructure == ChargedNeutralCharged) {
if(!parent->data().coloured()) return;
if(!back) {
ColinePair cparent = ColinePair(parent->colourLine(),
parent->antiColourLine());
// q -> q g
if(cparent.first) {
cparent.first->addColoured(second);
}
// qbar -> qbar g
if(cparent.second) {
cparent.second->addAntiColoured(second);
}
}
else {
if (second->dataPtr()->iColour()==PDT::Colour3 ) {
ColinePtr newline=new_ptr(ColourLine());
newline->addColoured(second);
newline->addColoured(parent);
}
else if (second->dataPtr()->iColour()==PDT::Colour3bar ) {
ColinePtr newline=new_ptr(ColourLine());
newline->addAntiColoured(second);
newline->addAntiColoured(parent);
}
}
}
else if(_colourStructure == NeutralChargedCharged ) {
if(!back) {
if(first->dataPtr()->coloured()) {
ColinePtr newline=new_ptr(ColourLine());
if(first->dataPtr()->iColour()==PDT::Colour3) {
newline->addColoured (first );
newline->addAntiColoured(second);
}
else if (first->dataPtr()->iColour()==PDT::Colour3bar) {
newline->addColoured (second);
newline->addAntiColoured(first );
}
else
assert(false);
}
}
else {
ColinePair cfirst = ColinePair(first->colourLine(),
first->antiColourLine());
// gamma -> q qbar
if(cfirst.first) {
cfirst.first->addAntiColoured(second);
}
// gamma -> qbar q
else if(cfirst.second) {
cfirst.second->addColoured(second);
}
else
assert(false);
}
}
else {
assert(false);
}
}
void SplittingFunction::doinit() {
Interfaced::doinit();
assert(_interactionType!=ShowerInteraction::UNDEFINED);
assert((_colourStructure>0&&_interactionType==ShowerInteraction::QCD) ||
(_colourStructure<0&&_interactionType==ShowerInteraction::QED) );
if(_colourFactor>0.) return;
// compute the colour factors if need
if(_colourStructure==TripletTripletOctet) {
_colourFactor = 4./3.;
}
else if(_colourStructure==OctetOctetOctet) {
_colourFactor = 3.;
}
else if(_colourStructure==OctetTripletTriplet) {
_colourFactor = 0.5;
}
else if(_colourStructure==TripletOctetTriplet) {
_colourFactor = 4./3.;
}
else if(_colourStructure==SextetSextetOctet) {
_colourFactor = 10./3.;
}
else if(_colourStructure<0) {
_colourFactor = 1.;
}
else {
assert(false);
}
}
bool SplittingFunction::checkColours(const IdList & ids) const {
tcPDPtr pd[3]={getParticleData(ids[0]),
getParticleData(ids[1]),
getParticleData(ids[2])};
if(_colourStructure==TripletTripletOctet) {
if(ids[0]!=ids[1]) return false;
if((pd[0]->iColour()==PDT::Colour3||pd[0]->iColour()==PDT::Colour3bar) &&
pd[2]->iColour()==PDT::Colour8) return true;
return false;
}
else if(_colourStructure==OctetOctetOctet) {
for(unsigned int ix=0;ix<3;++ix) {
if(pd[ix]->iColour()!=PDT::Colour8) return false;
}
return true;
}
else if(_colourStructure==OctetTripletTriplet) {
if(pd[0]->iColour()!=PDT::Colour8) return false;
if(pd[1]->iColour()==PDT::Colour3&&pd[2]->iColour()==PDT::Colour3bar)
return true;
if(pd[1]->iColour()==PDT::Colour3bar&&pd[2]->iColour()==PDT::Colour3)
return true;
return false;
}
else if(_colourStructure==TripletOctetTriplet) {
if(ids[0]!=ids[2]) return false;
if((pd[0]->iColour()==PDT::Colour3||pd[0]->iColour()==PDT::Colour3bar) &&
pd[1]->iColour()==PDT::Colour8) return true;
return false;
}
else if(_colourStructure==SextetSextetOctet) {
if(ids[0]!=ids[1]) return false;
if((pd[0]->iColour()==PDT::Colour6 || pd[0]->iColour()==PDT::Colour6bar) &&
pd[2]->iColour()==PDT::Colour8) return true;
return false;
}
else if(_colourStructure==ChargedChargedNeutral) {
if(ids[0]!=ids[1]) return false;
if(pd[2]->iCharge()!=0) return false;
if(pd[0]->iCharge()==pd[1]->iCharge()) return true;
return false;
}
else if(_colourStructure==ChargedNeutralCharged) {
if(ids[0]!=ids[2]) return false;
if(pd[1]->iCharge()!=0) return false;
if(pd[0]->iCharge()==pd[2]->iCharge()) return true;
return false;
}
else if(_colourStructure==NeutralChargedCharged) {
if(ids[1]!=-ids[2]) return false;
if(pd[0]->iCharge()!=0) return false;
if(pd[1]->iCharge()==-pd[2]->iCharge()) return true;
return false;
}
else {
assert(false);
}
return false;
}
+
+void SplittingFunction::evaluateFinalStateScales(ShowerPartnerType::Type type,
+ Energy scale, double z,
+ tShowerParticlePtr parent,
+ tShowerParticlePtr first,
+ tShowerParticlePtr second) {
+ // // identify emittor and spectator
+ // tShowerParticlePtr emittor = children[0];
+ // tShowerParticlePtr emitted = children[1];
+ // double zChild[2] = {z(),1.-z()};
+ // bool bosonSplitting(false);
+ // // special for g -> gg, particle highest z is emittor
+ // if(emittor->id()==emitted->id()&&emittor->id()==parent->id()) {
+ // if(zChild[1]>zChild[0]) {
+ // swap(zChild[0],zChild[1]);
+ // swap(emitted,emittor);
+ // }
+ // }
+ // // otherwise if particle ID same
+ // else if(emitted->id()==parent->id()) {
+ // swap(zChild[0],zChild[1]);
+ // swap(emitted,emittor);
+ // }
+ // // no real emittor/eemitted
+ // else if(emittor->id()!=parent->id()) {
+ // bosonSplitting = true;
+ // }
+ // // scales for angular ordering
+ // Energy AOScale[2];
+ // if(angularOrder) {
+ // for(unsigned int ix=0;ix<2;++ix) AOScale [ix] = zChild[ix]*scale();
+ // }
+ // else {
+ // for(unsigned int ix=0;ix<2;++ix) AOScale [ix] = scale();
+ // }
+ // // if not these cases doesn't matter
+ // // now the various scales
+ // if(partnerType==ShowerPartnerType::QED) {
+ // // normal case
+ // if(!bosonSplitting) {
+ // for(map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
+ // it = parent->evolutionScales().begin();
+ // it!=parent->evolutionScales().end();++it) {
+ // emittor->evolutionScale(it->first,make_pair(min(AOScale[0],it->second.first ),
+ // min(scale() ,it->second.second)));
+ // if(it->first==ShowerPartnerType::QED) {
+ // emitted->evolutionScale(it->first,make_pair(min(AOScale[1],it->second.first ),
+ // min(scale() ,it->second.second)));
+ // }
+ // }
+ // }
+ // // gamma -> f fbar
+ // else {
+ // // set QED scales as normal
+ // for(map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
+ // it = parent->evolutionScales().begin();
+ // it!=parent->evolutionScales().end();++it) {
+ // assert(it->first==ShowerPartnerType::QED);
+ // emittor->evolutionScale(it->first,make_pair(AOScale[0],scale()));
+ // emitted->evolutionScale(it->first,make_pair(AOScale[1],scale()));
+ // }
+ // // and any QCD scales needed
+ // PDT::Colour emittorColour = emittor->dataPtr()->iColour();
+ // if(emittorColour==PDT::Colour3||emittorColour==PDT::Colour8||emittorColour==PDT::Colour6) {
+ // emittor->evolutionScale(ShowerPartnerType::QCDColourLine,
+ // make_pair(AOScale[0],scale()));
+ // }
+ // if(emittorColour==PDT::Colour3bar||emittorColour==PDT::Colour8||emittorColour==PDT::Colour6bar) {
+ // emittor->evolutionScale(ShowerPartnerType::QCDAntiColourLine,
+ // make_pair(AOScale[0],scale()));
+ // }
+ // PDT::Colour emittedColour = emitted->dataPtr()->iColour();
+ // if(emittedColour==PDT::Colour3||emittedColour==PDT::Colour8||emittedColour==PDT::Colour6) {
+ // emitted->evolutionScale(ShowerPartnerType::QCDColourLine,
+ // make_pair(AOScale[1],scale()));
+ // }
+ // if(emittedColour==PDT::Colour3bar||emittedColour==PDT::Colour8||emittedColour==PDT::Colour6bar) {
+ // emitted->evolutionScale(ShowerPartnerType::QCDAntiColourLine,
+ // make_pair(AOScale[1],scale()));
+ // }
+ // }
+ // }
+ // else {
+ // // normal case
+ // if(!bosonSplitting) {
+ // // scales for the emittor
+ // for(map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
+ // it = parent->evolutionScales().begin();
+ // it!=parent->evolutionScales().end();++it) {
+ // emittor->evolutionScale(it->first,make_pair(min(AOScale[0],it->second.first ),
+ // min(scale() ,it->second.second)));
+ // }
+ // }
+ // else {
+ // // QCD scales needed
+ // PDT::Colour emittorColour = emittor->dataPtr()->iColour();
+ // if(emittorColour==PDT::Colour3||emittorColour==PDT::Colour8||emittorColour==PDT::Colour6) {
+ // emittor->evolutionScale(ShowerPartnerType::QCDColourLine,
+ // make_pair(AOScale[0],scale()));
+ // }
+ // if(emittorColour==PDT::Colour3bar||emittorColour==PDT::Colour8||emittorColour==PDT::Colour6bar) {
+ // emittor->evolutionScale(ShowerPartnerType::QCDAntiColourLine,
+ // make_pair(AOScale[0],scale()));
+ // }
+ // // QED scales
+ // if(emittor->dataPtr()->charged()) {
+ // emittor->evolutionScale(ShowerPartnerType::QED,make_pair(AOScale[0],scale()));
+ // }
+ // }
+ // PDT::Colour emittedColour = emitted->dataPtr()->iColour();
+ // if(emittedColour==PDT::Colour3||emittedColour==PDT::Colour8||emittedColour==PDT::Colour6) {
+ // emitted->evolutionScale(ShowerPartnerType::QCDColourLine,
+ // make_pair(AOScale[1],scale()));
+ // }
+ // if(emittedColour==PDT::Colour3bar||emittedColour==PDT::Colour8||emittedColour==PDT::Colour6bar) {
+ // emitted->evolutionScale(ShowerPartnerType::QCDAntiColourLine,
+ // make_pair(AOScale[1],scale()));
+ // }
+ // if(emitted->dataPtr()->charged()) {
+ // emitted->evolutionScale(ShowerPartnerType::QED,make_pair(AOScale[1],scale()));
+ // }
+ // }
+ assert(false);
+}
+
+void SplittingFunction::evaluateInitialStateScales(ShowerPartnerType::Type type,
+ Energy scale, double z,
+ tShowerParticlePtr parent,
+ tShowerParticlePtr first,
+ tShowerParticlePtr second) {
+ // // scale for time-like child
+ // Energy AOScale = (angularOrder ? (1.-z()) : 1. )*scale();
+ // // scales for parent if same as spacelike child
+ // if(parent->id()==children[0]->id()) {
+ // for(map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
+ // it = children[0]->evolutionScales().begin();
+ // it!=children[0]->evolutionScales().end();++it) {
+ // parent->evolutionScale(it->first,make_pair(min(scale(),it->second.first ),
+ // min(scale(),it->second.second)));
+ // }
+ // if(partnerType==ShowerPartnerType::QED) {
+ // children[1]->evolutionScale(partnerType,make_pair(AOScale,scale()));
+ // }
+ // PDT::Colour childColour = children[1]->dataPtr()->iColour();
+ // if(childColour==PDT::Colour3||childColour==PDT::Colour8||childColour==PDT::Colour6) {
+ // children[1]->evolutionScale(ShowerPartnerType::QCDColourLine,make_pair(AOScale,scale()));
+ // }
+ // if(childColour==PDT::Colour3bar||childColour==PDT::Colour8||childColour==PDT::Colour6bar) {
+ // children[1]->evolutionScale(ShowerPartnerType::QCDAntiColourLine,make_pair(AOScale,scale()));
+ // }
+ // }
+ // // scales if parent same as timelike child
+ // else if(parent->id()==children[1]->id()) {
+ // if(parent->dataPtr()->charged()) {
+ // parent ->evolutionScale(ShowerPartnerType::QED,make_pair(scale(),scale()));
+ // children[1]->evolutionScale(ShowerPartnerType::QED,make_pair(AOScale,scale()));
+ // }
+ // PDT::Colour childColour = children[1]->dataPtr()->iColour();
+ // if(partnerType==ShowerPartnerType::QED) {
+ // if(childColour==PDT::Colour3||childColour==PDT::Colour8||childColour==PDT::Colour6) {
+ // parent ->evolutionScale(ShowerPartnerType::QCDColourLine,make_pair(scale(),scale()));
+ // children[1]->evolutionScale(ShowerPartnerType::QCDColourLine,make_pair(AOScale,scale()));
+ // }
+ // if(childColour==PDT::Colour3bar||childColour==PDT::Colour8||childColour==PDT::Colour6bar) {
+ // parent ->evolutionScale(ShowerPartnerType::QCDAntiColourLine,make_pair(scale(),scale()));
+ // children[1]->evolutionScale(ShowerPartnerType::QCDAntiColourLine,make_pair(AOScale,scale()));
+ // }
+ // }
+ // else {
+ // if(childColour==PDT::Colour3||childColour==PDT::Colour8||childColour==PDT::Colour6) {
+ // pair<Energy,Energy> pScale = children[0]->evolutionScale(ShowerPartnerType::QCDColourLine);
+ // parent ->evolutionScale(ShowerPartnerType::QCDColourLine,
+ // make_pair(min(scale(),pScale.first ),
+ // min(scale(),pScale.second)));
+ // children[1]->evolutionScale(ShowerPartnerType::QCDColourLine,make_pair(AOScale,scale()));
+ // }
+ // if(childColour==PDT::Colour3bar||childColour==PDT::Colour8||childColour==PDT::Colour6bar) {
+ // pair<Energy,Energy> pScale = children[0]->evolutionScale(ShowerPartnerType::QCDAntiColourLine);
+ // parent ->evolutionScale(ShowerPartnerType::QCDAntiColourLine,
+ // make_pair(min(scale(),pScale.first ),
+ // min(scale(),pScale.second)));
+ // children[1]->evolutionScale(ShowerPartnerType::QCDAntiColourLine,make_pair(AOScale,scale()));
+ // }
+ // }
+ // }
+ // // g -> q qbar, or gamma -> e+e-
+ // else if(children[0]->id()==-children[1]->id()) {
+ // PDT::Colour childColour = children[1]->dataPtr()->iColour();
+ // if(partnerType==ShowerPartnerType::QED) {
+ // parent ->evolutionScale(ShowerPartnerType::QED,make_pair(scale(),scale()));
+ // children[1]->evolutionScale(ShowerPartnerType::QED,make_pair(AOScale,scale()));
+ // if(childColour==PDT::Colour3||childColour==PDT::Colour8||childColour==PDT::Colour6) {
+ // pair<Energy,Energy> pScale = children[0]->evolutionScale(ShowerPartnerType::QCDAntiColourLine);
+ // children[1]->evolutionScale(ShowerPartnerType::QCDColourLine,
+ // make_pair(min(AOScale,pScale.first ),
+ // min(AOScale,pScale.second)));
+ // }
+ // if(childColour==PDT::Colour3bar||childColour==PDT::Colour8||childColour==PDT::Colour6bar) {
+ // pair<Energy,Energy> pScale = children[0]->evolutionScale(ShowerPartnerType::QCDColourLine );
+ // children[1]->evolutionScale(ShowerPartnerType::QCDAntiColourLine,
+ // make_pair(min(AOScale,pScale.first ),
+ // min(AOScale,pScale.second)));
+ // }
+ // }
+ // else {
+ // pair<Energy,Energy> pScale = children[0]->evolutionScale(ShowerPartnerType::QED);
+ // children[1]->evolutionScale(ShowerPartnerType::QED,
+ // make_pair(min(AOScale,pScale.first ),
+ // min(AOScale,pScale.second)));
+ // if(childColour==PDT::Colour3||childColour==PDT::Colour8||childColour==PDT::Colour6) {
+ // pScale = children[0]->evolutionScale(ShowerPartnerType::QCDAntiColourLine);
+ // children[1]->evolutionScale(ShowerPartnerType::QCDColourLine,
+ // make_pair(min(AOScale,pScale.first ),
+ // min(AOScale,pScale.second)));
+ // }
+ // if(childColour==PDT::Colour3bar||childColour==PDT::Colour8||childColour==PDT::Colour6bar) {
+ // pScale = children[0]->evolutionScale(ShowerPartnerType::QCDColourLine );
+ // children[1]->evolutionScale(ShowerPartnerType::QCDAntiColourLine,
+ // make_pair(min(AOScale,pScale.first ),
+ // min(AOScale,pScale.second)));
+ // }
+ // parent->evolutionScale(ShowerPartnerType::QCDColourLine,
+ // make_pair(min(scale(),pScale.first ),
+ // min(scale(),pScale.second)));
+ // parent->evolutionScale(ShowerPartnerType::QCDAntiColourLine,
+ // make_pair(min(scale(),pScale.first ),
+ // min(scale(),pScale.second)));
+ // }
+ // }
+ assert(false);
+}
+
+void SplittingFunction::evaluateDecayScales(ShowerPartnerType::Type type,
+ Energy scale, double z,
+ tShowerParticlePtr parent,
+ tShowerParticlePtr first,
+ tShowerParticlePtr second) {
+ // // angular-ordered scale for 2nd child
+ // Energy AOScale = angularOrder ? (1.-z())*scale() : scale();
+ // // QED
+ // if(partnerType==ShowerPartnerType::QED) {
+ // for(map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
+ // it = parent->evolutionScales().begin();
+ // it!=parent->evolutionScales().end();++it) {
+ // children[0]->evolutionScale(it->first,make_pair(min(scale(),it->second.first ),
+ // min(scale(),it->second.second)));
+ // if(it->first==ShowerPartnerType::QED) {
+ // children[1]->evolutionScale(it->first,make_pair(min(AOScale,it->second.first ),
+ // min(scale(),it->second.second)));
+ // }
+ // }
+ // }
+ // // QCD
+ // else {
+ // // scales for the emittor
+ // for(map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
+ // it = parent->evolutionScales().begin();
+ // it!=parent->evolutionScales().end();++it) {
+ // children[0]->evolutionScale(it->first,make_pair(min(scale(),it->second.first ),
+ // min(scale(),it->second.second)));
+ // }
+ // PDT::Colour emittedColour = children[1]->dataPtr()->iColour();
+ // if(emittedColour==PDT::Colour3||emittedColour==PDT::Colour8||emittedColour==PDT::Colour6) {
+ // children[1]->evolutionScale(ShowerPartnerType::QCDColourLine,
+ // make_pair(AOScale,scale()));
+ // }
+ // if(emittedColour==PDT::Colour3bar||emittedColour==PDT::Colour8||emittedColour==PDT::Colour6bar) {
+ // children[1]->evolutionScale(ShowerPartnerType::QCDAntiColourLine,
+ // make_pair(AOScale,scale()));
+ // }
+ // if(children[1]->dataPtr()->charged()) {
+ // children[1]->evolutionScale(ShowerPartnerType::QED,make_pair(AOScale,scale()));
+ // }
+ // }
+ assert(false);
+}
diff --git a/Shower/SplittingFunctions/SplittingFunction.h b/Shower/SplittingFunctions/SplittingFunction.h
--- a/Shower/SplittingFunctions/SplittingFunction.h
+++ b/Shower/SplittingFunctions/SplittingFunction.h
@@ -1,314 +1,346 @@
// -*- C++ -*-
//
// SplittingFunction.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_SplittingFunction_H
#define HERWIG_SplittingFunction_H
//
// This is the declaration of the SplittingFunction class.
//
#include "ThePEG/Interface/Interfaced.h"
#include "Herwig++/Shower/ShowerConfig.h"
#include "ThePEG/EventRecord/ColourLine.h"
#include "ThePEG/PDT/ParticleData.h"
#include "SplittingFunction.fh"
namespace Herwig {
using namespace ThePEG;
/** \ingroup Shower
* Enum to define the possible types of colour structure which can occur in
* the branching.
*/
enum ColourStructure {Undefined=0,
TripletTripletOctet = 1,OctetOctetOctet =2,
OctetTripletTriplet = 3,TripletOctetTriplet=4,
SextetSextetOctet = 5,
ChargedChargedNeutral=-1,ChargedNeutralCharged=-2,
NeutralChargedCharged=-3};
/** \ingroup Shower
*
* This is an abstract class which defines the common interface
* for all \f$1\to2\f$ splitting functions, for both initial-state
* and final-state radiation.
*
* The SplittingFunction class contains a number of purely virtual members
* which must be implemented in the inheriting classes. The class also stores
* the interaction type of the spltting function.
*
* The inheriting classes need to specific the splitting function
* \f$P(z,2p_j\cdot p_k)\f$, in terms of the energy fraction \f$z\f$ and
* the evolution scale. In order to allow the splitting functions to be used
* with different choices of evolution functions the scale is given by
* \f[2p_j\cdot p_k=(p_j+p_k)^2-m_{jk}^2=Q^2-(p_j+p_k)^2=z(1-z)\tilde{q}^2=
* \frac{p_T^2}{z(1-z)}-m_{jk}^2+\frac{m_j^2}{z}+\frac{m_k^2}{1-z},\f]
* where \f$Q^2\f$ is the virtuality of the branching particle,
* $p_T$ is the relative transverse momentum of the branching products and
* \f$\tilde{q}^2\f$ is the angular variable described in hep-ph/0310083.
*
* In addition an overestimate of the
* splitting function, \f$P_{\rm over}(z)\f$ which only depends upon \f$z\f$,
* the integral, inverse of the integral for this overestimate and
* ratio of the true splitting function to the overestimate must be provided
* as they are necessary for the veto alogrithm used to implement the evolution.
*
* @see \ref SplittingFunctionInterfaces "The interfaces"
* defined for SplittingFunction.
*/
class SplittingFunction: public Interfaced {
public:
/**
* The default constructor.
* @param b All splitting functions must have an interaction order
*/
SplittingFunction(unsigned int b)
: Interfaced(), _interactionType(ShowerInteraction::UNDEFINED),
_interactionOrder(b),
_colourStructure(Undefined), _colourFactor(-1.),
angularOrdered_(true),
_splittingColourMethod(0) {}
public:
/**
* Methods to return the interaction type and order for the splitting function
*/
//@{
/**
* Return the type of the interaction
*/
ShowerInteraction::Type interactionType() const {return _interactionType;}
/**
* Return the order of the splitting function in the interaction
*/
unsigned int interactionOrder() const {return _interactionOrder;}
/**
* Return the colour structure
*/
ColourStructure colourStructure() const {return _colourStructure;}
/**
* Return the colour factor
*/
double colourFactor(const IdList &ids) const {
if(_colourStructure>0)
return _colourFactor;
else if(_colourStructure<0) {
if(_colourStructure==ChargedChargedNeutral ||
_colourStructure==ChargedNeutralCharged) {
tPDPtr part=getParticleData(ids[0]);
return sqr(double(part->iCharge())/3.);
}
else {
tPDPtr part=getParticleData(ids[1]);
return sqr(double(part->iCharge())/3.);
}
}
else
assert(false);
}
//@}
/**
* Purely virtual method which should determine whether this splitting
* function can be used for a given set of particles.
* @param ids The PDG codes for the particles in the splitting.
*/
virtual bool accept(const IdList & ids) const = 0;
/**
* Method to check the colours are correct
*/
virtual bool checkColours(const IdList & ids) const;
/**
* Methods to return the splitting function.
*/
//@{
/**
* Purely virtual method which should return the exact value of the splitting function,
* \f$P\f$ evaluated in terms of the energy fraction, \f$z\f$, and the evolution scale
\f$\tilde{q}^2\f$.
* @param z The energy fraction.
* @param t The scale \f$t=2p_j\cdot p_k\f$.
* @param ids The PDG codes for the particles in the splitting.
* @param mass Whether or not to include the mass dependent terms
*/
virtual double P(const double z, const Energy2 t, const IdList & ids,
const bool mass) const = 0;
/**
* Purely virtual method which should return
* an overestimate of the splitting function,
* \f$P_{\rm over}\f$ such that the result \f$P_{\rm over}\geq P\f$. This function
* should be simple enough that it does not depend on the evolution scale.
* @param z The energy fraction.
* @param ids The PDG codes for the particles in the splitting.
*/
virtual double overestimateP(const double z, const IdList & ids) const = 0;
/**
* Purely virtual method which should return
* the ratio of the splitting function to the overestimate, i.e.
* \f$P(z,\tilde{q}^2)/P_{\rm over}(z)\f$.
* @param z The energy fraction.
* @param t The scale \f$t=2p_j\cdot p_k\f$.
* @param ids The PDG codes for the particles in the splitting.
* @param mass Whether or not to include the mass dependent terms
*/
virtual double ratioP(const double z, const Energy2 t, const IdList & ids,
const bool mass) const = 0;
/**
* Purely virtual method which should return the indefinite integral of the
* overestimated splitting function, \f$P_{\rm over}\f$.
* @param z The energy fraction.
* @param ids The PDG codes for the particles in the splitting.
* @param PDFfactor Which additional factor to include for the PDF
* 0 is no additional factor,
* 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$
*
*/
virtual double integOverP(const double z, const IdList & ids,
unsigned int PDFfactor=0) const = 0;
/**
* Purely virtual method which should return the inverse of the
* indefinite integral of the
* overestimated splitting function, \f$P_{\rm over}\f$ which is used to
* generate the value of \f$z\f$.
* @param r Value of the splitting function to be inverted
* @param ids The PDG codes for the particles in the splitting.
* @param PDFfactor Which additional factor to include for the PDF
* 0 is no additional factor,
* 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$
*/
virtual double invIntegOverP(const double r, const IdList & ids,
unsigned int PDFfactor=0) const = 0;
//@}
/**
* Purely virtual method which should make the proper colour connection
* between the emitting parent and the branching products.
* @param parent The parent for the branching
* @param first The first branching product
* @param second The second branching product
* @param partnerType The type of evolution partner
* @param back Whether this is foward or backward evolution.
*/
virtual void colourConnection(tShowerParticlePtr parent,
tShowerParticlePtr first,
tShowerParticlePtr second,
ShowerPartnerType::Type partnerType,
const bool back) const;
/**
* Whether or not the interaction is angular ordered
*/
bool angularOrdered() const {return angularOrdered_;}
+ /**
+ * Functions to state scales after branching happens
+ */
+ //@{
+ /**
+ * Sort out scales for final-state emission
+ */
+ void evaluateFinalStateScales(ShowerPartnerType::Type type,
+ Energy scale, double z,
+ tShowerParticlePtr parent,
+ tShowerParticlePtr first,
+ tShowerParticlePtr second);
+
+ /**
+ * Sort out scales for initial-state emission
+ */
+ void evaluateInitialStateScales(ShowerPartnerType::Type type,
+ Energy scale, double z,
+ tShowerParticlePtr parent,
+ tShowerParticlePtr first,
+ tShowerParticlePtr second);
+
+ /**
+ * Sort out scales for decay emission
+ */
+ void evaluateDecayScales(ShowerPartnerType::Type type,
+ Energy scale, double z,
+ tShowerParticlePtr parent,
+ tShowerParticlePtr first,
+ tShowerParticlePtr second);
+ //@}
+
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:
/**
* Set the colour factor
*/
void colourFactor(double in) {_colourFactor=in;}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
SplittingFunction & operator=(const SplittingFunction &);
private:
/**
* The interaction type for the splitting function.
*/
ShowerInteraction::Type _interactionType;
/**
* The order of the splitting function in the coupling
*/
unsigned int _interactionOrder;
/**
* The colour structure
*/
ColourStructure _colourStructure;
/**
* The colour factor
*/
double _colourFactor;
/**
* Whether or not this interaction is angular-ordered
*/
bool angularOrdered_;
/**
* The method for assigning colour
* The default, 0, will assign colour lines for octets
* randomly without keeping a record of which lines radiate.
* For option 1 only the "correct" lines will radiate until
* the lowest scale is reached.
* For option 2 there will be random radiation, but the
* line which radiates is recorded
*/
int _splittingColourMethod;
};
}
#endif /* HERWIG_SplittingFunction_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,552 @@
// -*- 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/Reference.h"
#include "ThePEG/Interface/Command.h"
#include "ThePEG/Utilities/StringUtils.h"
#include "ThePEG/Repository/Repository.h"
#include "Herwig++/Shower/Base/ShowerParticle.h"
#include "Herwig++/Shower/ShowerHandler.h"
#include "ThePEG/Utilities/Rebinder.h"
#include <cassert>
using namespace Herwig;
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;
}
ClassDescription<SplittingGenerator> SplittingGenerator::initSplittingGenerator;
// Definition of the static class description member.
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;
- newKin = cit->second.first->
- generateNextTimeBranching(particle.evolutionScale(angularOrdered,type),
- 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;
- newKin= cit->second.first->
- generateNextTimeBranching(particle.evolutionScale(angularOrdered,type),
- cit->second.second,particle.id()!=cit->first,
- 0.5*enhance);
- ShoKinPtr newKin2 = cit->second.first->
- generateNextTimeBranching(particle.evolutionScale(angularOrdered,
- ShowerPartnerType::QCDAntiColourLine),
- 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
- else {
- Energy startingScale =
- max(particle.evolutionScale(angularOrdered,ShowerPartnerType::QCDColourLine),
- particle.evolutionScale(angularOrdered,ShowerPartnerType::QCDAntiColourLine));
- newKin= cit->second.first->
- generateNextTimeBranching(startingScale,
- cit->second.second,particle.id()!=cit->first,
- enhance);
- type = cit->second.second[0]<0 ?
- ShowerPartnerType::QCDColourLine : ShowerPartnerType::QCDAntiColourLine;
- }
- }
- // everything else
- else {
- type = particle.hasColour() ?
- ShowerPartnerType::QCDColourLine : ShowerPartnerType::QCDAntiColourLine;
- newKin= cit->second.first->
- generateNextTimeBranching(particle.evolutionScale(angularOrdered,type),
- cit->second.second,
- particle.id()!=cit->first,enhance);
- }
- }
- // shouldn't be anything else
- else
- assert(false);
+ // // work out which starting scale we need
+ // if(cit->second.first->interactionType()==ShowerInteraction::QED) {
+ // type = ShowerPartnerType::QED;
+ // newKin = cit->second.first->
+ // generateNextTimeBranching(particle.evolutionScale(angularOrdered,type),
+ // 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;
+ // newKin= cit->second.first->
+ // generateNextTimeBranching(particle.evolutionScale(angularOrdered,type),
+ // cit->second.second,particle.id()!=cit->first,
+ // 0.5*enhance);
+ // ShoKinPtr newKin2 = cit->second.first->
+ // generateNextTimeBranching(particle.evolutionScale(angularOrdered,
+ // ShowerPartnerType::QCDAntiColourLine),
+ // 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
+ // else {
+ // Energy startingScale =
+ // max(particle.evolutionScale(angularOrdered,ShowerPartnerType::QCDColourLine),
+ // particle.evolutionScale(angularOrdered,ShowerPartnerType::QCDAntiColourLine));
+ // newKin= cit->second.first->
+ // generateNextTimeBranching(startingScale,
+ // cit->second.second,particle.id()!=cit->first,
+ // enhance);
+ // type = cit->second.second[0]<0 ?
+ // ShowerPartnerType::QCDColourLine : ShowerPartnerType::QCDAntiColourLine;
+ // }
+ // }
+ // // everything else
+ // else {
+ // type = particle.hasColour() ?
+ // ShowerPartnerType::QCDColourLine : ShowerPartnerType::QCDAntiColourLine;
+ // newKin= cit->second.first->
+ // generateNextTimeBranching(particle.evolutionScale(angularOrdered,type),
+ // cit->second.second,
+ // particle.id()!=cit->first,enhance);
+ // }
+ // }
+ // // shouldn't be anything else
+ // else
+ // assert(false);
+ 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 return it
return Branching(kinematics, ids,sudakov,partnerType);
}
Branching SplittingGenerator::
chooseDecayBranching(ShowerParticle &particle,
const map<ShowerPartnerType::Type,pair<Energy,Energy> > & 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 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;
- map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
- it=stoppingScales.find(type);
- if(it==stoppingScales.end()) continue;
- Energy stoppingScale = angularOrdered ? it->second.first : it->second.second;
- if(particle.evolutionScale(angularOrdered,type) < stoppingScale ) {
- newKin = cit->second.first->
- generateNextDecayBranching(particle.evolutionScale(angularOrdered,type),
- stoppingScale,minmass,cit->second.second,
- particle.id()!=cit->first,enhance);
- cerr << "testing in QED " << particle << "\n";
- assert(false);
- }
- }
- else if(cit->second.first->interactionType()==ShowerInteraction::QCD) {
- // special for octets
- if(particle.dataPtr()->iColour()==PDT::Colour8) {
- map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
- it=stoppingScales.find(ShowerPartnerType::QCDColourLine);
- if(it==stoppingScales.end()) continue;
- Energy stoppingColour = angularOrdered ? it->second.first : it->second.second;
- it=stoppingScales.find(ShowerPartnerType::QCDAntiColourLine);
- Energy stoppingAnti = angularOrdered ? it->second.first : it->second.second;
- // octet -> octet octet
- if(cit->second.first->splittingFn()->colourStructure()==OctetOctetOctet) {
- type = ShowerPartnerType::QCDColourLine;
- newKin= cit->second.first->
- generateNextDecayBranching(particle.evolutionScale(angularOrdered,type),
- stoppingColour,minmass,
- cit->second.second,
- particle.id()!=cit->first,0.5*enhance);
- ShoKinPtr newKin2 = cit->second.first->
- generateNextDecayBranching(particle.evolutionScale(angularOrdered,
- ShowerPartnerType::QCDAntiColourLine),
- 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 {
- Energy startingScale =
- min(particle.evolutionScale(angularOrdered,ShowerPartnerType::QCDColourLine),
- particle.evolutionScale(angularOrdered,ShowerPartnerType::QCDAntiColourLine));
- newKin = cit->second.first->
- generateNextDecayBranching(startingScale,
- max(stoppingColour,stoppingAnti),minmass,
- cit->second.second,
- particle.id()!=cit->first,enhance);
- type = cit->second.second[0]<0 ?
- ShowerPartnerType::QCDColourLine : ShowerPartnerType::QCDAntiColourLine;
- }
- }
- // everything else
- else {
- type = particle.hasColour() ?
- ShowerPartnerType::QCDColourLine : ShowerPartnerType::QCDAntiColourLine;
- map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
- it=stoppingScales.find(type);
- if(it==stoppingScales.end()) continue;
- Energy stoppingScale = angularOrdered ? it->second.first : it->second.second;
- if(particle.evolutionScale(angularOrdered,type) < stoppingScale ) {
- newKin = cit->second.first->
- generateNextDecayBranching(particle.evolutionScale(angularOrdered,type),
- stoppingScale,minmass,cit->second.second,
- particle.id()!=cit->first,enhance);
- }
- }
- }
- // shouldn't be anything else
- else
- assert(false);
+ // if(cit->second.first->interactionType()==ShowerInteraction::QED) {
+ // type = ShowerPartnerType::QED;
+ // map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
+ // it=stoppingScales.find(type);
+ // if(it==stoppingScales.end()) continue;
+ // Energy stoppingScale = angularOrdered ? it->second.first : it->second.second;
+ // if(particle.evolutionScale(angularOrdered,type) < stoppingScale ) {
+ // newKin = cit->second.first->
+ // generateNextDecayBranching(particle.evolutionScale(angularOrdered,type),
+ // stoppingScale,minmass,cit->second.second,
+ // particle.id()!=cit->first,enhance);
+ // cerr << "testing in QED " << particle << "\n";
+ // assert(false);
+ // }
+ // }
+ // else if(cit->second.first->interactionType()==ShowerInteraction::QCD) {
+ // // special for octets
+ // if(particle.dataPtr()->iColour()==PDT::Colour8) {
+ // map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
+ // it=stoppingScales.find(ShowerPartnerType::QCDColourLine);
+ // if(it==stoppingScales.end()) continue;
+ // Energy stoppingColour = angularOrdered ? it->second.first : it->second.second;
+ // it=stoppingScales.find(ShowerPartnerType::QCDAntiColourLine);
+ // Energy stoppingAnti = angularOrdered ? it->second.first : it->second.second;
+ // // octet -> octet octet
+ // if(cit->second.first->splittingFn()->colourStructure()==OctetOctetOctet) {
+ // type = ShowerPartnerType::QCDColourLine;
+ // newKin= cit->second.first->
+ // generateNextDecayBranching(particle.evolutionScale(angularOrdered,type),
+ // stoppingColour,minmass,
+ // cit->second.second,
+ // particle.id()!=cit->first,0.5*enhance);
+ // ShoKinPtr newKin2 = cit->second.first->
+ // generateNextDecayBranching(particle.evolutionScale(angularOrdered,
+ // ShowerPartnerType::QCDAntiColourLine),
+ // 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 {
+ // Energy startingScale =
+ // min(particle.evolutionScale(angularOrdered,ShowerPartnerType::QCDColourLine),
+ // particle.evolutionScale(angularOrdered,ShowerPartnerType::QCDAntiColourLine));
+ // newKin = cit->second.first->
+ // generateNextDecayBranching(startingScale,
+ // max(stoppingColour,stoppingAnti),minmass,
+ // cit->second.second,
+ // particle.id()!=cit->first,enhance);
+ // type = cit->second.second[0]<0 ?
+ // ShowerPartnerType::QCDColourLine : ShowerPartnerType::QCDAntiColourLine;
+ // }
+ // }
+ // // everything else
+ // else {
+ // type = particle.hasColour() ?
+ // ShowerPartnerType::QCDColourLine : ShowerPartnerType::QCDAntiColourLine;
+ // map<ShowerPartnerType::Type,pair<Energy,Energy> >::const_iterator
+ // it=stoppingScales.find(type);
+ // if(it==stoppingScales.end()) continue;
+ // Energy stoppingScale = angularOrdered ? it->second.first : it->second.second;
+ // if(particle.evolutionScale(angularOrdered,type) < stoppingScale ) {
+ // newKin = cit->second.first->
+ // generateNextDecayBranching(particle.evolutionScale(angularOrdered,type),
+ // stoppingScale,minmass,cit->second.second,
+ // particle.id()!=cit->first,enhance);
+ // }
+ // }
+ // }
+ // // shouldn't be anything else
+ // else
+ // assert(false);
+ 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 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;
- newKin=cit->second.first->
- generateNextSpaceBranching(particle.evolutionScale(angularOrdered,type),
- cit->second.second, particle.x(),
- particle.id()!=cit->first,enhance,beam);
- cerr << "testing in QED " << particle << "\n";
- assert(false);
- }
- 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;
- newKin=cit->second.first->
- generateNextSpaceBranching(particle.evolutionScale(angularOrdered,type),
- cit->second.second, particle.x(),
- particle.id()!=cit->first,0.5*enhance,beam);
- ShoKinPtr newKin2 = cit->second.first->
- generateNextSpaceBranching(particle.evolutionScale(angularOrdered,
- ShowerPartnerType::QCDAntiColourLine),
- 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;
- }
- }
- // other
- else {
- Energy startingScale =
- max(particle.evolutionScale(angularOrdered,ShowerPartnerType::QCDColourLine),
- particle.evolutionScale(angularOrdered,ShowerPartnerType::QCDAntiColourLine));
- type = UseRandom::rndbool() ?
- ShowerPartnerType::QCDColourLine : ShowerPartnerType::QCDAntiColourLine;
- newKin=cit->second.first->
- generateNextSpaceBranching(particle.evolutionScale(angularOrdered,type),
- cit->second.second, particle.x(),
- particle.id()!=cit->first,enhance,beam);
- }
- }
- // everything else
- else {
- type = particle.hasColour() ?
- ShowerPartnerType::QCDColourLine : ShowerPartnerType::QCDAntiColourLine;
- newKin=cit->second.first->
- generateNextSpaceBranching(particle.evolutionScale(angularOrdered,type),
- cit->second.second, particle.x(),
- particle.id()!=cit->first,enhance,beam);
- }
- }
- // shouldn't be anything else
- else
- assert(false);
+ // if(cit->second.first->interactionType()==ShowerInteraction::QED) {
+ // type = ShowerPartnerType::QED;
+ // newKin=cit->second.first->
+ // generateNextSpaceBranching(particle.evolutionScale(angularOrdered,type),
+ // cit->second.second, particle.x(),
+ // particle.id()!=cit->first,enhance,beam);
+ // cerr << "testing in QED " << particle << "\n";
+ // assert(false);
+ // }
+ // 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;
+ // newKin=cit->second.first->
+ // generateNextSpaceBranching(particle.evolutionScale(angularOrdered,type),
+ // cit->second.second, particle.x(),
+ // particle.id()!=cit->first,0.5*enhance,beam);
+ // ShoKinPtr newKin2 = cit->second.first->
+ // generateNextSpaceBranching(particle.evolutionScale(angularOrdered,
+ // ShowerPartnerType::QCDAntiColourLine),
+ // 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;
+ // }
+ // }
+ // // other
+ // else {
+ // Energy startingScale =
+ // max(particle.evolutionScale(angularOrdered,ShowerPartnerType::QCDColourLine),
+ // particle.evolutionScale(angularOrdered,ShowerPartnerType::QCDAntiColourLine));
+ // type = UseRandom::rndbool() ?
+ // ShowerPartnerType::QCDColourLine : ShowerPartnerType::QCDAntiColourLine;
+ // newKin=cit->second.first->
+ // generateNextSpaceBranching(particle.evolutionScale(angularOrdered,type),
+ // cit->second.second, particle.x(),
+ // particle.id()!=cit->first,enhance,beam);
+ // }
+ // }
+ // // everything else
+ // else {
+ // type = particle.hasColour() ?
+ // ShowerPartnerType::QCDColourLine : ShowerPartnerType::QCDAntiColourLine;
+ // newKin=cit->second.first->
+ // generateNextSpaceBranching(particle.evolutionScale(angularOrdered,type),
+ // cit->second.second, particle.x(),
+ // particle.id()!=cit->first,enhance,beam);
+ // }
+ // }
+ // // shouldn't be anything else
+ // else
+ // assert(false);
+ 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);
// return the answer
return Branching(kinematics, ids,sudakov,partnerType);
}
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, 10:33 PM (1 d, 14 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4230067
Default Alt Text
(274 KB)

Event Timeline