Page MenuHomeHEPForge
Files F8308742

No OneTemporary
Actions

View Options

diff --git a/Shower/Base/Evolver.cc b/Shower/Base/Evolver.cc
--- a/Shower/Base/Evolver.cc
+++ b/Shower/Base/Evolver.cc
@@ -1,1503 +1,1525 @@
// -*- C++ -*-
//
// Evolver.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the 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 "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;
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 << _limitEmissions
<< ounit(_iptrms,GeV) << _beta << ounit(_gamma,GeV) << ounit(_iptmax,GeV)
<< _vetoes << _hardonly << _trunc_Mode << _hardEmissionMode;
}
void Evolver::persistentInput(PersistentIStream & is, int) {
is >> _model >> _splittingGenerator >> _maxtry
>> _meCorrMode >> _hardVetoMode >> _hardVetoRead >> _limitEmissions
>> iunit(_iptrms,GeV) >> _beta >> iunit(_gamma,GeV) >> iunit(_iptmax,GeV)
>> _vetoes >> _hardonly >> _trunc_Mode >> _hardEmissionMode;
}
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 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 Switch<Evolver,bool> interfaceHardOnly
("HardOnly",
"Only generate the emission supplied by the hardest emission"
" generator, for testing only.",
&Evolver::_hardonly, false, false, false);
static SwitchOption interfaceHardOnlyNo
(interfaceHardOnly,
"No",
"Generate full shower",
false);
static SwitchOption interfaceHardOnlyYes
(interfaceHardOnly,
"Yes",
"Only the hardest emission",
true);
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);
}
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(ShowerTreePtr hard,
vector<ShowerProgenitorPtr> p) {
// find out if hard partonic subprocess.
bool isPartonic(false);
map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
cit = _currenttree->incomingLines().begin();
Lorentz5Momentum pcm;
for(; cit!=hard->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()) {
// scattering process
if(hard->isHard()) {
// coloured incoming particles
if (isPartonic) {
map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
cjt = hard->outgoingLines().begin();
for(; cjt!=hard->outgoingLines().end(); ++cjt) {
if (cjt->first->progenitor()->coloured())
ptmax = max(ptmax,cjt->first->progenitor()->momentum().mt());
}
}
if (ptmax < ZERO) ptmax = pcm.m();
}
// decay, incoming() is the decaying particle.
else {
ptmax = hard->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);
// zero number of emissions
_nis = _nfs = 0;
// extract particles to shower
vector<ShowerProgenitorPtr> particlesToShower=setupShower(true);
// setup the maximum scales for the shower, given by the hard process
if (hardVetoOn()) setupMaximumScales(currentTree(), particlesToShower);
// generate the intrinsic p_T once and for all
generateIntrinsicpT(particlesToShower);
// main shower loop
unsigned int ntry(0);
do {
// clear results of last attempt if needed
if(ntry!=0) {
currentTree()->clear();
setEvolutionPartners(true,ShowerInteraction::QCD);
_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;
}
// set the progenitor
_progenitor=particlesToShower[ix];
// initial-state
if(!_progenitor->progenitor()->isFinalState()) {
if(!isISRadiationON()) continue;
// 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,
ShowerInteraction::QCD));
}
// final-state
else {
if(!isFSRadiationON()) continue;
// perform shower
progenitor()->hasEmitted(startTimeLikeShower(ShowerInteraction::QCD));
}
}
}
while(!showerModel()->kinematicsReconstructor()->
reconstructHardJets(hard,intrinsicpT())&&
maximumTries()>++ntry);
_hardme=HwMEBasePtr();
_hardtree=HardTreePtr();
if(_maxtry==ntry) throw ShowerHandler::ShowerTriesVeto(ntry);
// the tree has now showered
_currenttree->hasShowered(true);
}
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) {
+ ShowerInteraction::Type type,
+ bool first) {
// don't do anything if not needed
if(_limitEmissions == 1 || _limitEmissions == 3 ||
( _limitEmissions == 2 && _nfs != 0) ) return false;
- // 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)) break;
- // otherwise reset scale and continue - SO IS involved in veto algorithm
- particle->setEvolutionScale(fb.kinematics->scale());
- }
- // has emitted
- // Assign the shower kinematics to the emitting particle.
- particle->setShowerKinematics(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;
+ ShowerParticleVector theChildren;
+ do {
+ // 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)) break;
+ // otherwise reset scale and continue - SO IS involved in veto algorithm
+ particle->setEvolutionScale(fb.kinematics->scale());
+ }
+ // has emitted
+ // Assign the shower kinematics to the emitting particle.
+ particle->setShowerKinematics(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,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);
+ // branching has happened
+ particle->showerKinematics()->updateParent(particle, theChildren,true);
+ // clean up the vetoed emission
+ if(particle->virtualMass()==ZERO) {
+ particle->setShowerKinematics(ShoKinPtr());
+ for(unsigned int ix=0;ix<theChildren.size();++ix)
+ particle->abandonChild(theChildren[ix]);
+ theChildren.clear();
}
}
- ShowerParticleVector theChildren;
- 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,true);
- // update the history if needed
- if(particle==_currenttree->getFinalStateShowerProduct(_progenitor))
- _currenttree->updateFinalStateShowerProduct(_progenitor,
- particle,theChildren);
- _currenttree->addFinalStateBranching(particle,theChildren);
- // update number of emissions
- ++_nfs;
- if(_limitEmissions!=0) return true;
- // shower the first particle
- timeLikeShower(theChildren[0],type);
- // shower the second particle
- timeLikeShower(theChildren[1],type);
- // branching has happened
+ while(particle->virtualMass()==ZERO);
+ if(first) {
+// cerr << "testing before reset " << *particle << "\n";
+// cerr << "testing before reset " << particle->showerKinematics() << "\n";
+ particle->showerKinematics()->resetChildren(particle,theChildren);
+// cerr << "testing after reset " << *particle << "\n";
+ }
+ // needs sorting out
+// // update the history if needed
+// if(particle==_currenttree->getFinalStateShowerProduct(_progenitor))
+// _currenttree->updateFinalStateShowerProduct(_progenitor,
+// particle,theChildren);
+// _currenttree->addFinalStateBranching(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 || _limitEmissions == 3 ||
( _limitEmissions == 1 && _nis != 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)) break;
// otherwise reset scale and continue
particle->setEvolutionScale(bb.kinematics->scale());
}
// assign the splitting function and shower kinematics
particle->setShowerKinematics(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,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,true);
if(_limitEmissions!=0) return true;
// perform the shower of the final-state particle
- timeLikeShower(otherChild,type);
+ 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();
// extract particles to be shower, set scales and
// perform hard matrix element correction
vector<ShowerProgenitorPtr> particlesToShower=setupShower(false);
setupMaximumScales(currentTree(), particlesToShower);
// compute the minimum mass of the final-state
Energy minmass(ZERO);
for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
if(particlesToShower[ix]->progenitor()->isFinalState())
minmass+=particlesToShower[ix]->progenitor()->mass();
}
// main showering loop
unsigned int ntry(0);
do {
// clear results of last attempt
if(ntry!=0) {
currentTree()->clear();
setEvolutionPartners(false,ShowerInteraction::QCD);
}
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(ShowerInteraction::QCD));
}
// initial-state radiation
else {
if(!isISRadiationON()) continue;
// perform shower
// set the scales correctly. The current scale is the maximum scale for
// emission not the starting scale
Energy maxscale=progenitor()->progenitor()->evolutionScale();
Energy startScale=progenitor()->progenitor()->mass();
progenitor()->progenitor()->setEvolutionScale(startScale);
// perform the shower
progenitor()->hasEmitted(startSpaceLikeDecayShower(maxscale,minmass,
ShowerInteraction::QCD));
}
}
}
while(!showerModel()->kinematicsReconstructor()->reconstructDecayJets(decay)&&
maximumTries()>++ntry);
_decayme = HwDecayerBasePtr();
if(maximumTries()==ntry)
throw Exception() << "Failed to generate the shower after "
<< ntry << " attempts in Evolver::showerDecay()"
<< Exception::eventerror;
// tree has now showered
_currenttree->hasShowered(true);
}
bool Evolver::spaceLikeDecayShower(tShowerParticlePtr particle,
Energy maxscale,
Energy minmass,ShowerInteraction::Type type) {
Branching fb;
while (true) {
fb=_splittingGenerator->chooseDecayBranching(*particle,maxscale,minmass,
_initialenhance,type);
// return if no radiation
if(!fb.kinematics) return false;
// if not vetoed break
if(!spaceLikeDecayVetoed(fb,particle)) break;
// otherwise reset scale and continue
particle->setEvolutionScale(fb.kinematics->scale());
}
// has emitted
// Assign the shower kinematics to the emitting particle.
particle->setShowerKinematics(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,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],maxscale,minmass,type);
// shower the second particle
- timeLikeShower(theChildren[1],type);
+ 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,ShowerInteraction::QCD);
// get the particles to be showered
map<ShowerProgenitorPtr, ShowerParticlePtr>::const_iterator cit;
map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator cjt;
// generate hard me if needed
if(_hardEmissionMode==0) hardMatrixElementCorrection(hard);
// get the particles to be showered
vector<ShowerProgenitorPtr> particlesToShower;
// incoming particles
for(cit=currentTree()->incomingLines().begin();
cit!=currentTree()->incomingLines().end();++cit)
particlesToShower.push_back((*cit).first);
assert((particlesToShower.size()==1&&!hard)||
(particlesToShower.size()==2&&hard));
// outgoing particles
for(cjt=currentTree()->outgoingLines().begin();
cjt!=currentTree()->outgoingLines().end();++cjt)
particlesToShower.push_back(((*cjt).first));
// remake the colour partners if needed
if(_hardEmissionMode==0 && _currenttree->hardMatrixElementCorrection()) {
setEvolutionPartners(hard,ShowerInteraction::QCD);
_currenttree->resetShowerProducts();
}
return particlesToShower;
}
void Evolver::setEvolutionPartners(bool hard,ShowerInteraction::Type ) {
map<ShowerProgenitorPtr, ShowerParticlePtr>::const_iterator cit;
map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator cjt;
vector<ShowerParticlePtr> particles;
// match the particles in the ShowerTree and hardTree
if(hardTree() && !hardTree()->connect(currentTree()))
throw Exception() << "Can't match trees in "
<< "Evolver::setEvolutionPartners()"
<< Exception::eventerror;
// sort out the colour partners
for(cit=currentTree()->incomingLines().begin();
cit!=currentTree()->incomingLines().end();++cit)
particles.push_back(cit->first->progenitor());
assert((particles.size()==1&&!hard)||(particles.size()==2&&hard));
// outgoing particles
for(cjt=currentTree()->outgoingLines().begin();
cjt!=currentTree()->outgoingLines().end();++cjt)
particles.push_back(cjt->first->progenitor());
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]->setPartner(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,ShowerInteraction::QCD,!_hardtree);
}
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() ) {
return truncatedTimeLikeShower(progenitor()->progenitor(), mit->second ,type);
}
}
return hardOnly() ? false :
- timeLikeShower(progenitor()->progenitor() ,type) ;
+ timeLikeShower(progenitor()->progenitor() ,type,true) ;
}
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(Energy maxscale,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(),maxscale,
minimumMass, branch ,type);
}
}
return hardOnly() ? false :
spaceLikeDecayShower(progenitor()->progenitor(),maxscale,minimumMass,type);
}
bool Evolver::timeLikeVetoed(const Branching & fb,
ShowerParticlePtr particle) {
// 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()) 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) {
// 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()) 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 ) {
// 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()) 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)
_hardtree = _hardme->generateHardest( currentTree() );
else
_hardtree = _decayme->generateHardest( currentTree() );
if(!_hardtree) return;
// join up the two tree
connectTrees(currentTree(),_hardtree,hard);
}
else {
_hardtree = ShowerHandler::currentHandler()->generateCKKW(currentTree());
}
}
bool Evolver::truncatedTimeLikeShower(tShowerParticlePtr particle,
HardBranchingPtr branch,
ShowerInteraction::Type type) {
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->setEvolutionScale(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->setEvolutionScale(fb.kinematics->scale());
continue;
}
}
// pt veto
if(fb.kinematics->pT() > progenitor()->maximumpT()) {
particle->setEvolutionScale(fb.kinematics->scale());
continue;
}
// should do base class vetos as well
if(timeLikeVetoed(fb,particle)) {
particle->setEvolutionScale(fb.kinematics->scale());
continue;
}
break;
}
// if no branching set decay matrix and return
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->setEvolutionScale(branch->scale() );
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() );
// Assign the shower kinematics to the emitting particle.
particle->setShowerKinematics( 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,type==branch->sudakov()->interactionType());
// update the history if needed
if(particle==currentTree()->getFinalStateShowerProduct(progenitor()))
currentTree()->updateFinalStateShowerProduct(progenitor(),
particle,theChildren);
currentTree()->addFinalStateBranching(particle,theChildren);
// shower the first particle
if( branch->children()[0]->children().empty() ) {
if( ! hardOnly() )
- timeLikeShower(theChildren[0],type);
+ 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);
+ timeLikeShower( theChildren[1] , type,false);
}
else {
truncatedTimeLikeShower( theChildren[1],branch->children()[1] ,type);
}
return true;
}
// has emitted
// Assign the shower kinematics to the emitting particle.
particle->setShowerKinematics(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 , true);
// update the history if needed
if( particle == currentTree()->getFinalStateShowerProduct( progenitor() ) )
currentTree()->updateFinalStateShowerProduct( progenitor(),
particle, theChildren );
currentTree()->addFinalStateBranching( particle, theChildren );
// shower the first particle
if( iout == 1 ) truncatedTimeLikeShower( theChildren[0], branch , type );
- else timeLikeShower( theChildren[0] , type);
+ else timeLikeShower( theChildren[0] , type,false);
// shower the second particle
if( iout == 2 ) truncatedTimeLikeShower( theChildren[1], branch , type );
- else timeLikeShower( theChildren[1] , type);
+ else timeLikeShower( theChildren[1] , type,false);
// branching has happened
return true;
}
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;
// generate branching
tcPDPtr part[2];
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->setEvolutionScale(bb.kinematics->scale() );
continue;
}
// others
if( part[0]->id() != particle->id() || // if particle changes type
bb.kinematics->pT() > progenitor()->maximumpT() || // pt veto
bb.kinematics->scale() < branch->scale()) { // angular ordering veto
particle->setEvolutionScale(bb.kinematics->scale() );
continue;
}
// and those from the base class
if(spaceLikeVetoed(bb,particle)) {
particle->setEvolutionScale(bb.kinematics->scale() );
continue;
}
break;
}
if( !bb.kinematics ) {
//do the hard emission
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();
}
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->setShowerKinematics( 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, 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,
type==branch->sudakov()->interactionType() );
if(hardOnly()) return true;
// perform the shower of the final-state particle
if( timelike->children().empty() ) {
- timeLikeShower( otherChild , type);
+ timeLikeShower( otherChild , type,true);
}
else {
truncatedTimeLikeShower( otherChild, timelike , type);
}
// return the emitted
return true;
}
// assign the splitting function and shower kinematics
particle->setShowerKinematics( 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 , 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 , true);
// perform the shower of the final-state particle
- timeLikeShower( otherChild , type);
+ timeLikeShower( otherChild , type,true);
// return the emitted
return true;
}
bool Evolver::
truncatedSpaceLikeDecayShower(tShowerParticlePtr particle, Energy maxscale,
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,maxscale,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->setEvolutionScale(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->setEvolutionScale(fb.kinematics->scale());
continue;
}
}
// pt veto
if(fb.kinematics->pT() > progenitor()->maximumpT()) {
particle->setEvolutionScale(fb.kinematics->scale());
continue;
}
// should do base class vetos as well
// if not vetoed break
if(!spaceLikeDecayVetoed(fb,particle)) break;
// otherwise reset scale and continue
particle->setEvolutionScale(fb.kinematics->scale());
}
// if no branching set decay matrix and return
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->setEvolutionScale(branch->scale() );
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() );
// Assign the shower kinematics to the emitting particle.
particle->setShowerKinematics( 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,
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],maxscale,minmass,type);
}
else {
truncatedSpaceLikeDecayShower( theChildren[0],maxscale,minmass,
branch->children()[0],type);
}
// shower the second particle
if( branch->children()[1]->children().empty() ) {
- if( ! hardOnly() ) timeLikeShower( theChildren[1] , type);
+ 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],maxscale,minmass,type);
}
else {
truncatedSpaceLikeDecayShower( theChildren[1],maxscale,minmass,
branch->children()[1],type);
}
// shower the second particle
if( branch->children()[0]->children().empty() ) {
- if( ! hardOnly() ) timeLikeShower( theChildren[0] , type);
+ 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->setShowerKinematics(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,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],maxscale,minmass,branch,type);
// shower the second particle
- timeLikeShower(theChildren[1],type);
+ timeLikeShower(theChildren[1],type,true);
// branching has happened
return true;
}
void Evolver::connectTrees(ShowerTreePtr showerTree, HardTreePtr hardTree, bool hard )const {
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) {
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()) {
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(),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);
}
}
diff --git a/Shower/Base/Evolver.h b/Shower/Base/Evolver.h
--- a/Shower/Base/Evolver.h
+++ b/Shower/Base/Evolver.h
@@ -1,675 +1,676 @@
// -*- C++ -*-
//
// Evolver.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_Evolver_H
#define HERWIG_Evolver_H
//
// This is the declaration of the Evolver class.
//
#include "ThePEG/Interface/Interfaced.h"
#include "Herwig++/Shower/SplittingFunctions/SplittingGenerator.h"
#include "ShowerModel.h"
#include "ThePEG/PDF/BeamParticleData.h"
#include "ShowerTree.h"
#include "ShowerProgenitor.fh"
#include "Herwig++/Shower/ShowerHandler.fh"
#include "Branching.h"
#include "ShowerVeto.h"
#include "HardTree.h"
#include "ThePEG/Handlers/XComb.h"
#include "Evolver.fh"
#include "Herwig++/MatrixElement/HwMEBase.h"
#include "Herwig++/Decay/HwDecayerBase.h"
namespace Herwig {
using namespace ThePEG;
/** \ingroup Shower
* The Evolver class class performs the sohwer evolution of hard scattering
* and decay processes in Herwig++.
*
* @see \ref EvolverInterfaces "The interfaces"
* defined for Evolver.
*/
class Evolver: public Interfaced {
/**
* The ShowerHandler is a friend to set some parameters at initialisation
*/
friend class ShowerHandler;
public:
/**
* Pointer to an XComb object
*/
typedef Ptr<XComb>::pointer XCPtr;
public:
/**
* Default Constructor
*/
Evolver() : _maxtry(100), _meCorrMode(1), _hardVetoMode(1),
_hardVetoRead(0),
_iptrms(ZERO), _beta(0.), _gamma(ZERO), _iptmax(),
_limitEmissions(0), _initialenhance(1.), _finalenhance(1.),
_hardonly(false), _trunc_Mode(true), _hardEmissionMode(0)
{}
/**
* Members to perform the shower
*/
//@{
/**
* Perform the shower of the hard process
*/
virtual void showerHardProcess(ShowerTreePtr,XCPtr);
/**
* Perform the shower of a decay
*/
virtual void showerDecay(ShowerTreePtr);
//@}
/**
* Access to the flags and shower variables
*/
//@{
/**
* Is there any showering switched on
*/
bool showeringON() const { return isISRadiationON() || isFSRadiationON(); }
/**
* It returns true/false if the initial-state radiation is on/off.
*/
bool isISRadiationON() const { return _splittingGenerator->isISRadiationON(); }
/**
* It returns true/false if the final-state radiation is on/off.
*/
bool isFSRadiationON() const { return _splittingGenerator->isFSRadiationON(); }
/**
* Get the ShowerModel
*/
ShowerModelPtr showerModel() const {return _model;}
/**
* Get the SplittingGenerator
*/
tSplittingGeneratorPtr splittingGenerator() const { return _splittingGenerator; }
//@}
/**
* Connect the Hard and Shower trees
*/
virtual void connectTrees(ShowerTreePtr showerTree, HardTreePtr hardTree, bool hard )const;
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:
/**
* Generate the hard matrix element correction
*/
virtual void hardMatrixElementCorrection(bool);
/**
* Generate the hardest emission
*/
virtual void hardestEmission(bool hard);
/**
* Extract the particles to be showered, set the evolution scales
* and apply the hard matrix element correction
* @param hard Whether this is a hard process or decay
* @return The particles to be showered
*/
virtual vector<ShowerProgenitorPtr> setupShower(bool hard);
/**
* set the colour partners
*/
virtual void setEvolutionPartners(bool hard,ShowerInteraction::Type);
/**
* Methods to perform the evolution of an individual particle, including
* recursive calling on the products
*/
//@{
/**
* It does the forward evolution of the time-like input particle
* (and recursively for all its radiation products).
* accepting only emissions which conforms to the showerVariables
* and soft matrix element correction.
* If at least one emission has occurred then the method returns true.
* @param particle The particle to be showered
*/
- virtual bool timeLikeShower(tShowerParticlePtr particle, ShowerInteraction::Type);
+ virtual bool timeLikeShower(tShowerParticlePtr particle, ShowerInteraction::Type,
+ bool first);
/**
* It does the backward evolution of the space-like input particle
* (and recursively for all its time-like radiation products).
* accepting only emissions which conforms to the showerVariables.
* If at least one emission has occurred then the method returns true
* @param particle The particle to be showered
* @param beam The beam particle
*/
virtual bool spaceLikeShower(tShowerParticlePtr particle,PPtr beam,
ShowerInteraction::Type);
/**
* If does the forward evolution of the input on-shell particle
* involved in a decay
* (and recursively for all its time-like radiation products).
* accepting only emissions which conforms to the showerVariables.
* @param particle The particle to be showered
* @param maxscale The maximum scale for the shower.
* @param minimumMass The minimum mass of the final-state system
*/
virtual bool spaceLikeDecayShower(tShowerParticlePtr particle,
Energy maxscale,
Energy minimumMass,
ShowerInteraction::Type);
/**
* Truncated shower from a time-like particle
*/
virtual bool truncatedTimeLikeShower(tShowerParticlePtr particle,
HardBranchingPtr branch,
ShowerInteraction::Type type);
/**
* Truncated shower from a space-like particle
*/
virtual bool truncatedSpaceLikeShower(tShowerParticlePtr particle,PPtr beam,
HardBranchingPtr branch,
ShowerInteraction::Type type);
/**
* Truncated shower from a time-like particle
*/
virtual bool truncatedSpaceLikeDecayShower(tShowerParticlePtr particle,
Energy maxscale, Energy minimumMass,
HardBranchingPtr branch,
ShowerInteraction::Type type);
//@}
/**
* Switches for matrix element corrections
*/
//@{
/**
* Any ME correction?
*/
bool MECOn() const {
return _meCorrMode > 0 && _hardEmissionMode==0;
}
/**
* Any hard ME correction?
*/
bool hardMEC() const {
return (_meCorrMode == 1 || _meCorrMode == 2) && _hardEmissionMode==0;
}
/**
* Any soft ME correction?
*/
bool softMEC() const {
return (_meCorrMode == 1 || _meCorrMode > 2) && _hardEmissionMode==0;
}
//@}
/**
* Is the truncated shower on?
*/
bool isTruncatedShowerON() const {return _trunc_Mode;}
/**
* Switch for intrinsic pT
*/
//@{
/**
* Any intrinsic pT?
*/
bool ipTon() const {
return _iptrms != ZERO || ( _beta == 1.0 && _gamma != ZERO && _iptmax !=ZERO );
}
//@}
/**@name Additional shower vetoes */
//@{
/**
* Insert a veto.
*/
void addVeto (ShowerVetoPtr v) { _vetoes.push_back(v); }
/**
* Remove a veto.
*/
void removeVeto (ShowerVetoPtr v) {
vector<ShowerVetoPtr>::iterator vit = find(_vetoes.begin(),_vetoes.end(),v);
if (vit != _vetoes.end())
_vetoes.erase(vit);
}
//@}
/**
* Switches for vetoing hard emissions
*/
//@{
/**
* Vetos on?
*/
bool hardVetoOn() const { return _hardVetoMode > 0; }
/**
* veto hard emissions in IS shower?
*/
bool hardVetoIS() const { return _hardVetoMode == 1 || _hardVetoMode == 2; }
/**
* veto hard emissions in FS shower?
*/
bool hardVetoFS() const { return _hardVetoMode == 1 || _hardVetoMode > 2; }
/**
* veto hard emissions according to lastScale from XComb?
*/
bool hardVetoXComb() const {return (_hardVetoRead == 1);}
//@}
/**
* Enhancement factors for radiation needed to generate the soft matrix
* element correction.
*/
//@{
/**
* Access the enhancement factor for initial-state radiation
*/
double initialStateRadiationEnhancementFactor() const { return _initialenhance; }
/**
* Access the enhancement factor for final-state radiation
*/
double finalStateRadiationEnhancementFactor() const { return _finalenhance; }
/**
* Set the enhancement factor for initial-state radiation
*/
void initialStateRadiationEnhancementFactor(double in) { _initialenhance=in; }
/**
* Set the enhancement factor for final-state radiation
*/
void finalStateRadiationEnhancementFactor(double in) { _finalenhance=in; }
//@}
/**
* Access to set/get the HardTree currently beinging showered
*/
//@{
/**
* The HardTree currently being showered
*/
tHardTreePtr hardTree() {return _hardtree;}
/**
* The HardTree currently being showered
*/
void hardTree(tHardTreePtr in) {_hardtree = in;}
//@}
/**
* Access/set the beam particle for the current initial-state shower
*/
//@{
/**
* Get the beam particle data
*/
Ptr<BeamParticleData>::const_pointer beamParticle() const { return _beam; }
/**
* Set the beam particle data
*/
void setBeamParticle(Ptr<BeamParticleData>::const_pointer in) { _beam=in; }
//@}
/**
* Set/Get the current tree being evolverd for inheriting classes
*/
//@{
/**
* Get the tree
*/
tShowerTreePtr currentTree() { return _currenttree; }
/**
* Set the tree
*/
void currentTree(tShowerTreePtr tree) { _currenttree=tree; }
//@}
/**
* Access the maximum number of attempts to generate the shower
*/
unsigned int maximumTries() const { return _maxtry; }
/**
* Set/Get the ShowerProgenitor for the current shower
*/
//@{
/**
* Access the progenitor
*/
ShowerProgenitorPtr progenitor() { return _progenitor; }
/**
* Set the progenitor
*/
void progenitor(ShowerProgenitorPtr in) { _progenitor=in; }
//@}
/**
* Calculate the intrinsic \f$p_T\f$.
*/
virtual void generateIntrinsicpT(vector<ShowerProgenitorPtr>);
/**
* Access to the intrinsic \f$p_T\f$ for inheriting classes
*/
map<tShowerProgenitorPtr,pair<Energy,double> > & intrinsicpT() { return _intrinsic; }
/**
* find the maximally allowed pt acc to the hard process.
*/
void setupMaximumScales(ShowerTreePtr, vector<ShowerProgenitorPtr>);
protected:
/**
* Start the shower of a timelike particle
*/
virtual bool startTimeLikeShower(ShowerInteraction::Type);
/**
* Start the shower of a spacelike particle
*/
virtual bool startSpaceLikeShower(PPtr,ShowerInteraction::Type);
/**
* Start the shower of a spacelike particle
*/
virtual bool startSpaceLikeDecayShower(Energy maxscale,Energy minimumMass,
ShowerInteraction::Type);
/**
* Vetos for the timelike shower
*/
virtual bool timeLikeVetoed(const Branching &,ShowerParticlePtr);
/**
* Vetos for the spacelike shower
*/
virtual bool spaceLikeVetoed(const Branching &,ShowerParticlePtr);
/**
* Vetos for the spacelike shower
*/
virtual bool spaceLikeDecayVetoed(const Branching &,ShowerParticlePtr);
/**
* Only generate the hard emission, for testing only.
*/
bool hardOnly() const {return _hardonly;}
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<Evolver> initEvolver;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
Evolver & operator=(const Evolver &);
private:
/**
* Pointer to the model for the shower evolution model
*/
ShowerModelPtr _model;
/**
* Pointer to the splitting generator
*/
SplittingGeneratorPtr _splittingGenerator;
/**
* Maximum number of tries to generate the shower of a particular tree
*/
unsigned int _maxtry;
/**
* Matrix element correction switch
*/
unsigned int _meCorrMode;
/**
* Hard emission veto switch
*/
unsigned int _hardVetoMode;
/**
* Hard veto to be read switch
*/
unsigned int _hardVetoRead;
/**
* rms intrinsic pT of Gaussian distribution
*/
Energy _iptrms;
/**
* Proportion of inverse quadratic intrinsic pT distribution
*/
double _beta;
/**
* Parameter for inverse quadratic: 2*Beta*Gamma/(sqr(Gamma)+sqr(intrinsicpT))
*/
Energy _gamma;
/**
* Upper bound on intrinsic pT for inverse quadratic
*/
Energy _iptmax;
/**
* Limit the number of emissions for testing
*/
unsigned int _limitEmissions;
/**
* The progenitor of the current shower
*/
ShowerProgenitorPtr _progenitor;
/**
* Matrix element
*/
HwMEBasePtr _hardme;
/**
* Decayer
*/
HwDecayerBasePtr _decayme;
/**
* The ShowerTree currently being showered
*/
ShowerTreePtr _currenttree;
/**
* The HardTree currently being showered
*/
HardTreePtr _hardtree;
/**
* Radiation enhancement factors for use with the veto algorithm
* if needed by the soft matrix element correction
*/
//@{
/**
* Enhancement factor for initial-state radiation
*/
double _initialenhance;
/**
* Enhancement factor for final-state radiation
*/
double _finalenhance;
//@}
/**
* The beam particle data for the current initial-state shower
*/
Ptr<BeamParticleData>::const_pointer _beam;
/**
* Storage of the intrinsic \f$p_t\f$ of the particles
*/
map<tShowerProgenitorPtr,pair<Energy,double> > _intrinsic;
/**
* Vetoes
*/
vector<ShowerVetoPtr> _vetoes;
/**
* number of IS emissions
*/
unsigned int _nis;
/**
* Number of FS emissions
*/
unsigned int _nfs;
/**
* Only generate the emission from the hardest emission
* generate for testing only
*/
bool _hardonly;
/**
* Truncated shower switch
*/
bool _trunc_Mode;
/**
* Count of the number of truncated emissions
*/
unsigned int _truncEmissions;
/**
* Mode for the hard emissions
*/
unsigned int _hardEmissionMode;
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of Evolver. */
template <>
struct BaseClassTrait<Herwig::Evolver,1> {
/** Typedef of the first base class of Evolver. */
typedef Interfaced NthBase;
};
/** This template specialization informs ThePEG about the name of
* the Evolver class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::Evolver>
: public ClassTraitsBase<Herwig::Evolver> {
/** Return a platform-independent class name */
static string className() { return "Herwig::Evolver"; }
/**
* The name of a file containing the dynamic library where the class
* Evolver is implemented. It may also include several, space-separated,
* libraries if the class Evolver depends on other classes (base classes
* excepted). In this case the listed libraries will be dynamically
* linked in the order they are specified.
*/
static string library() { return "HwShower.so"; }
};
/** @endcond */
}
#endif /* HERWIG_Evolver_H */
diff --git a/Shower/Base/ShowerKinematics.cc b/Shower/Base/ShowerKinematics.cc
--- a/Shower/Base/ShowerKinematics.cc
+++ b/Shower/Base/ShowerKinematics.cc
@@ -1,64 +1,70 @@
// -*- C++ -*-
//
// ShowerKinematics.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the ShowerKinematics class.
//
#include "ShowerKinematics.h"
using namespace Herwig;
void ShowerKinematics::updateChildren(const tShowerParticlePtr,
const ShowerParticleVector &,
bool ) const {
throw Exception() << "Base class ShowerKinematics::updateChildren called,"
<< " should have been overriden in an inheriting class"
<< Exception::runerror;
}
+void ShowerKinematics::resetChildren(const tShowerParticlePtr,
+ const ShowerParticleVector &) const {
+ throw Exception() << "Base class ShowerKinematics::resetChildren called,"
+ << " should have been overriden in an inheriting class"
+ << Exception::runerror;
+}
void ShowerKinematics::updateParent(const tShowerParticlePtr,
const ShowerParticleVector &,
bool) const {
throw Exception() << "Base class ShowerKinematics::updateParent called,"
<< " should have been overriden in an inheriting class"
<< Exception::runerror;
}
void ShowerKinematics::reconstructChildren(const tShowerParticlePtr,
const ShowerParticleVector &) const {
throw Exception() << "Base class ShowerKinematics::reconstructChildren called,"
<< " should have been overriden in an inheriting class"
<< Exception::runerror;
}
void ShowerKinematics::reconstructParent(const tShowerParticlePtr,
const ParticleVector &) const {
throw Exception() << "Base class ShowerKinematics::reconstructParent called,"
<< " should have been overriden in an inheriting class"
<< Exception::runerror;
}
void ShowerKinematics::reconstructLast(const tShowerParticlePtr,
unsigned int,Energy) const {
throw Exception() << "Base class ShowerKinematics::reconstructLast called,"
<< " should have been overriden in an inheriting class"
<< Exception::runerror;
}
void ShowerKinematics::updateLast(const tShowerParticlePtr,
Energy,Energy) const {
throw Exception() << "Base class ShowerKinematics::updatetLast called,"
<< " should have been overriden in an inheriting class"
<< Exception::runerror;
}
void ShowerKinematics::initialize(ShowerParticle &,PPtr) {
throw Exception() << "Base class ShowerKinematics::initialize called "
<< Exception::runerror;
}
diff --git a/Shower/Base/ShowerKinematics.h b/Shower/Base/ShowerKinematics.h
--- a/Shower/Base/ShowerKinematics.h
+++ b/Shower/Base/ShowerKinematics.h
@@ -1,317 +1,320 @@
// -*- C++ -*-
//
// ShowerKinematics.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_ShowerKinematics_H
#define HERWIG_ShowerKinematics_H
//
// This is the declaration of the ShowerKinematics class.
//
#include "Herwig++/Shower/ShowerConfig.h"
#include "ThePEG/Config/ThePEG.h"
#include "Herwig++/Shower/Base/SudakovFormFactor.h"
#include "ShowerKinematics.fh"
namespace Herwig {
using namespace ThePEG;
/**\ingroup Shower
*
* This is the abstract base class from which all other shower
* kinematics classes derive. The main purpose of the
* shower kinematics classes is to allow the reconstruction
* of jet masses, at the end of the showering (indeed, for
* multi-scale showering, at the end of each scale-range evolution).
* This is necessary for the kinematics reshuffling
* in order to compensate the recoil of the emissions.
* The KinematicsReconstructor class is in
* charge of this job, and which is the main "user" of
* ShowerKinematics and its derived classes.
* How this is done depends on the choice of kinematics variables
* and whether the jet is time-like (forward evolved) or
* space-like (backward evolved), whereas the class ShowerKinematics
* describes only the common features which are independent by them.
*
* In general there are a number of methods specific to a shower approach
*
* @see KinematicsReconstructor
*/
class ShowerKinematics: public Base {
public:
/**
* The default constructor.
*/
ShowerKinematics() : Base(), _isTheJetStartingPoint( false ),
_scale(), _z( 0.0 ), _phi( 0.0 ), _pt(),
_sudakov() {}
/**
* The updateChildren and updateParent
* members to update the values of the \f$\alpha\f$ and
* \f$p_\perp\f$ variables during the shower evolution.
*/
//@{
/**
* Along with the showering evolution --- going forward for
* time-like (forward) evolution, and going backward for space-like
* (backward) evolution --- the kinematical variables of the
* branching products are calculated and updated from the knowledge
* of the parent kinematics.
* @param theParent The parent
* @param theChildren The children
* @param angularOrder Whether or not to apply angular ordering
*/
virtual void updateChildren(const tShowerParticlePtr theParent,
const ShowerParticleVector & theChildren,
bool angularOrder) const;
+ virtual void resetChildren( const tShowerParticlePtr theParent,
+ const ShowerParticleVector & theChildren) const;
+
/**
* Update the parent Kinematics from the knowledge of the kinematics
* of the children. This method will be used by the KinematicsReconstructor.
* @param theParent The parent
* @param theChildren The children
* @param angularOrder Whether or not to apply angular ordering
*/
virtual void updateParent(const tShowerParticlePtr theParent,
const ShowerParticleVector & theChildren,
bool angularOrder) const;
/**
* Update the kinematical data of a particle when a reconstruction
* fixpoint was found. This will highly depend on the kind of
* kinematics chosen and will be defined in the inherited concrete
* classes. This method will be used by the KinematicsReconstructor.
* @param theLast The particle.
* @param px The \f$x\f$ component of the \f$p_T\f$.
* @param py The \f$y\f$ component of the \f$p_T\f$.
*/
virtual void updateLast(const tShowerParticlePtr theLast,
Energy px, Energy py) const;
//@}
/**
* The reconstructLast, reconstructChildren and reconstructParent members
* are used during the reconstruction
*/
//@{
/**
* Along with the showering evolution --- going forward for
* time-like (forward) evolution, and going backward for space-like
* (backward) evolution --- the kinematical variables of the
* branching products are calculated and updated from the knowledge
* of the parent kinematics.
* @param theParent The parent
* @param theChildren The children
*/
virtual void reconstructChildren(const tShowerParticlePtr theParent,
const ShowerParticleVector & theChildren) const;
/**
* Reconstruct the parent Kinematics from the knowledge of the kinematics
* of the children. This method will be used by the KinematicsReconstructor.
* @param theParent The parent
* @param theChildren The children
*/
virtual void reconstructParent(const tShowerParticlePtr theParent,
const ParticleVector & theChildren) const;
/**
* Update the kinematical data of a particle when a reconstruction
* fixpoint was found. This will highly depend on the kind of
* kinematics chosen and will be defined in the inherited concrete
* classes. This method will be used by the KinematicsReconstructor.
* @param theLast The particle.
* @param iopt The option for the momentum reconstruction
* - 0 is in the rest frame of the pair of reference vectors
* - 1 is in the rest frame of the p vector
* @param mass The mass to be used, if less than zero on-shell
*/
virtual void reconstructLast(const tShowerParticlePtr theLast,
unsigned int iopt, Energy mass=-1.*GeV) const;
/**
* Perform any initial calculations needed after the branching has been selected
* @param particle The branching particle
* @param parent The bema particle for the jet if needed
*/
virtual void initialize(ShowerParticle & particle,PPtr parent);
//@}
public:
/**
* Set/access the flag that tells whether or not this ShowerKinematics
* object is associated to the starting particle of the jet: only in this
* case it is sensible to use the two main virtual methods below.
*/
//@{
/**
* Set the starting point flag
*/
void isTheJetStartingPoint(const bool );
/**
* Get the starting point flag
*/
bool isTheJetStartingPoint() const;
//@}
/**
* Virtual function to return a set of basis vectors, specific to
* the type of evolution. This function will be used by the
* ForwardShowerEvolver in order to access \f$p\f$ and \f$n\f$,
* which in turn are members of the concrete class QTildeShowerKinematics1to2.
*/
virtual vector<Lorentz5Momentum> getBasis() const = 0;
/**
* Set/Get methods for the kinematic variables
*/
//@{
/**
* Access the scale of the splitting.
*/
Energy scale() const { return _scale; }
/**
* Set the scale of the splitting.
*/
void scale(const Energy in) { _scale=in; }
/**
* Access the energy fraction, \f$z\f$.
*/
double z() const { return _z; }
/**
* Set the energy fraction, \f$z\f$.
*/
void z(const double in) { _z=in; }
/**
* Access the azimuthal angle, \f$\phi\f$.
*/
double phi() const { return _phi; }
/**
* Set the azimuthal angle, \f$\phi\f$.
*/
void phi(const double in) { _phi=in; }
/**
* Access the relative \f$p_T\f$ for the branching
*/
Energy pT() const { return _pt; }
/**
* Set the relative \f$p_T\f$ for the branching
*/
- void pT(const Energy in) { _pt=in; }
+ void pT(const Energy in) const { _pt=in; }
//@}
/**
* Set and get methods for the SplittingFunction object
*/
//@{
/**
* Access the SplittingFunction object responsible of the
* eventual branching of this particle.
*/
tSplittingFnPtr splittingFn() const { return _sudakov-> splittingFn(); }
//@}
/**
* Set and get methods for the SudakovFormFactor object
*/
/**
* Access the SudakovFormFactor object responsible of the
* eventual branching of this particle.
*/
tSudakovPtr SudakovFormFactor() const { return _sudakov; }
/**
* Set the SudakovFormFactor object responsible of the
* eventual branching of this particle.
*/
void SudakovFormFactor(const tSudakovPtr sud) { _sudakov=sud; }
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
ShowerKinematics & operator=(const ShowerKinematics &);
private:
/**
* Is this the starting point of the jet
*/
bool _isTheJetStartingPoint;
/**
* The \f$\tilde{q}\f$ evolution variable.
*/
Energy _scale;
/**
* The energy fraction, \f$z\f$
*/
double _z;
/**
* The azimuthal angle, \f$\phi\f$.
*/
double _phi;
/**
* The relative \f$p_T\f$
*/
- Energy _pt;
+ mutable Energy _pt;
/**
* The splitting function for the branching of the particle
*/
tSudakovPtr _sudakov;
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of ShowerKinematics. */
template <>
struct BaseClassTrait<Herwig::ShowerKinematics,1> {
/** Typedef of the first base class of ShowerKinematics. */
typedef Base NthBase;
};
/** This template specialization informs ThePEG about the name of
* the ShowerKinematics class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::ShowerKinematics>
: public ClassTraitsBase<Herwig::ShowerKinematics> {
/** Return a platform-independent class name */
static string className() { return "Herwig::ShowerKinematics"; }
};
/** @endcond */
}
#endif /* HERWIG_ShowerKinematics_H */
diff --git a/Shower/Base/ShowerParticle.h b/Shower/Base/ShowerParticle.h
--- a/Shower/Base/ShowerParticle.h
+++ b/Shower/Base/ShowerParticle.h
@@ -1,337 +1,351 @@
// -*- C++ -*-
//
// ShowerParticle.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_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:
/** @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), _reconstructionFixedPoint( false ),
_perturbative(0), _initiatesTLS(tls), _x(1.0), _showerKinematics(),
- _scale(ZERO), _thePEGBase() {}
+ _scale(ZERO), _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), _reconstructionFixedPoint( false ),
_perturbative(pert), _initiatesTLS(tls), _x(1.0), _showerKinematics(),
- _scale(ZERO), _thePEGBase(&x) {}
+ _scale(ZERO), _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 setShowerKinematics(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$
*/
Energy evolutionScale() const { return _scale; }
/**
* Set the evolution \f$\tilde{q}\f$ scale
*/
void setEvolutionScale(Energy scale) { _scale = scale; }
+ /**
+ * Return the virtual mass\f$
+ */
+ Energy virtualMass() const { return _vMass; }
+
+ /**
+ * Set the virtual mass
+ */
+ void setVirtualMass(Energy mass) { _vMass = mass; }
/**
* Return the partner
*/
tShowerParticlePtr partner() const { return _partner; }
/**
* Set the partner
*/
void setPartner(const tShowerParticlePtr partner) { _partner = partner; }
//@}
/**
* Access/Set the flag that tells if the particle should be
* treated in a special way during the kinematics reconstruction
* (see KinematicsReconstructor class).
* In practice, it returns true when either the particle is childless,
* or is a on-shell decaying particle (in which case we have to set the flag to
* true before the showering of this particle: it is not enough to check
* if decayer() is not null, because if it emits radiation
* the decays products will be "transferred" to the particle
* instance after the showering).
*/
//@{
/**
* Get the flag
*/
bool isReconstructionFixedPoint() const { return _reconstructionFixedPoint || children().empty(); }
/**
* Set the flag
*/
void setReconstructionFixedPoint(const bool in) { _reconstructionFixedPoint = in; }
//@}
/**
* 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 getThePEGBase() 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 is a reconstruction fixed point
*/
bool _reconstructionFixedPoint;
/**
* 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;
/**
* Evolution scales
*/
Energy _scale;
/**
+ * Virtual mass
+ */
+ Energy _vMass;
+
+ /**
* Partners
*/
tShowerParticlePtr _partner;
/**
* Pointer to ThePEG Particle this ShowerParticle was created from
*/
const tcPPtr _thePEGBase;
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of ShowerParticle. */
template <>
struct BaseClassTrait<Herwig::ShowerParticle,1> {
/** Typedef of the first base class of ShowerParticle. */
typedef Particle NthBase;
};
/** This template specialization informs ThePEG about the name of
* the ShowerParticle class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::ShowerParticle>
: public ClassTraitsBase<Herwig::ShowerParticle> {
/** Return a platform-independent class name */
static string className() { return "Herwig::ShowerParticle"; }
/** Create a Event object. */
static TPtr create() { return TPtr::Create(Herwig::ShowerParticle(tcEventPDPtr(),true)); }
};
/** @endcond */
}
#endif /* HERWIG_ShowerParticle_H */
diff --git a/Shower/Base/ShowerTree.cc b/Shower/Base/ShowerTree.cc
--- a/Shower/Base/ShowerTree.cc
+++ b/Shower/Base/ShowerTree.cc
@@ -1,914 +1,914 @@
// -*- C++ -*-
//
// ShowerTree.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#include "ShowerProgenitor.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ShowerTree.h"
#include "Herwig++/Shower/Base/ShowerParticle.h"
#include "ThePEG/PDT/DecayMode.h"
#include "ThePEG/Handlers/EventHandler.h"
#include "ThePEG/Handlers/XComb.h"
#include "KinematicsReconstructor.h"
#include <cassert>
#include "ThePEG/Repository/CurrentGenerator.h"
using namespace Herwig;
using namespace ThePEG;
set<long> ShowerTree::_decayInShower = set<long>();
namespace {
void findBeam(tPPtr & beam, PPtr incoming) {
while(!beam->children().empty()) {
bool found=false;
for(unsigned int ix=0;ix<beam->children().size();++ix) {
if(beam->children()[ix]==incoming) {
found = true;
break;
}
}
if(found) break;
beam = beam->children()[0];
}
}
}
// constructor from hard process
ShowerTree::ShowerTree(const PPair incoming, const ParticleVector & out,
ShowerDecayMap& decay)
: _hardMECorrection(false), _wasHard(true),
_parent(), _hasShowered(false) {
tPPair beam = CurrentGenerator::current().currentEvent()->incoming();
findBeam(beam.first ,incoming.first );
findBeam(beam.second,incoming.second);
_incoming = incoming;
double x1(_incoming.first ->momentum().rho()/beam.first ->momentum().rho());
double x2(_incoming.second->momentum().rho()/beam.second->momentum().rho());
// must have two incoming particles
assert(_incoming.first && _incoming.second);
// set the parent tree
_parent=ShowerTreePtr();
// temporary vectors to contain all the particles before insertion into
// the data structure
vector<PPtr> original,copy;
vector<ShowerParticlePtr> shower;
// create copies of ThePEG particles for the incoming particles
original.push_back(_incoming.first);
copy.push_back(new_ptr(Particle(*_incoming.first)));
original.push_back(_incoming.second);
copy.push_back(new_ptr(Particle(*_incoming.second)));
// and same for outgoing
map<PPtr,ShowerTreePtr> trees;
for (ParticleVector::const_iterator it = out.begin();
it != out.end(); ++it) {
// if decayed or should be decayed in shower make the tree
PPtr orig = *it;
if(!orig->children().empty()||
(decaysInShower(orig->id())&&!orig->dataPtr()->stable())) {
ShowerTreePtr newtree=new_ptr(ShowerTree(orig,decay));
newtree->setParents();
trees.insert(make_pair(orig,newtree));
Energy width=orig->dataPtr()->generateWidth(orig->mass());
decay.insert(make_pair(width,newtree));
}
original.push_back(orig);
copy.push_back(new_ptr(Particle(*orig)));
}
// colour isolate the hard process
colourIsolate(original,copy);
// now create the Shower particles
// create ShowerParticles for the incoming particles
assert(original.size() == copy.size());
for(unsigned int ix=0;ix<original.size();++ix) {
ShowerParticlePtr temp=new_ptr(ShowerParticle(*copy[ix],1,ix>=2));
fixColour(temp);
// incoming
if(ix<2) {
temp->x(ix==0 ? x1 : x2);
_incomingLines.insert(make_pair(new_ptr(ShowerProgenitor(original[ix],
copy[ix],temp)),
temp));
_backward.insert(temp);
}
// outgoing
else {
_outgoingLines.insert(make_pair(new_ptr(ShowerProgenitor(original[ix],
copy[ix],temp)),
temp));
_forward.insert(temp);
}
}
// set up the map of daughter trees
map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator mit;
for(mit=_outgoingLines.begin();mit!=_outgoingLines.end();++mit) {
map<PPtr,ShowerTreePtr>::const_iterator tit=trees.find(mit->first->original());
if(tit!=trees.end())
_treelinks.insert(make_pair(tit->second,
make_pair(mit->first,mit->first->progenitor())));
}
}
ShowerTree::ShowerTree(PPtr in,
ShowerDecayMap& decay)
: _hardMECorrection(false), _wasHard(false), _hasShowered(false) {
// there must be an incoming particle
assert(in);
// temporary vectors to contain all the particles before insertion into
// the data structure
vector<PPtr> original,copy;
// insert place holder for incoming particle
original.push_back(in);
copy.push_back(PPtr());
// we need to deal with the decay products if decayed
map<PPtr,ShowerTreePtr> trees;
if(!in->children().empty()) {
ParticleVector children=in->children();
for(unsigned int ix=0;ix<children.size();++ix) {
// if decayed or should be decayed in shower make the tree
PPtr orig=children[ix];
in->abandonChild(orig);
if(!orig->children().empty()||
(decaysInShower(orig->id())&&!orig->dataPtr()->stable())) {
ShowerTreePtr newtree=new_ptr(ShowerTree(orig,decay));
trees.insert(make_pair(orig,newtree));
Energy width=orig->dataPtr()->generateWidth(orig->mass());
decay.insert(make_pair(width,newtree));
newtree->setParents();
newtree->_parent=this;
}
original.push_back(orig);
copy.push_back(new_ptr(Particle(*orig)));
}
}
// create the incoming particle
copy[0] = new_ptr(Particle(*in));
// isolate the colour
colourIsolate(original,copy);
// create the parent
ShowerParticlePtr sparent(new_ptr(ShowerParticle(*copy[0],2,false)));
fixColour(sparent);
_incomingLines.insert(make_pair(new_ptr(ShowerProgenitor(original[0],copy[0],sparent))
,sparent));
// return if not decayed
if(original.size()==1) return;
// create the children
assert(copy.size() == original.size());
for (unsigned int ix=1;ix<original.size();++ix) {
ShowerParticlePtr stemp= new_ptr(ShowerParticle(*copy[ix],2,true));
fixColour(stemp);
_outgoingLines.insert(make_pair(new_ptr(ShowerProgenitor(original[ix],copy[ix],
stemp)),
stemp));
_forward.insert(stemp);
}
// set up the map of daughter trees
map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator mit;
for(mit=_outgoingLines.begin();mit!=_outgoingLines.end();++mit) {
map<PPtr,ShowerTreePtr>::const_iterator tit=trees.find(mit->first->original());
if(tit!=trees.end())
_treelinks.insert(make_pair(tit->second,
make_pair(mit->first,mit->first->progenitor())));
}
}
void ShowerTree::updateFinalStateShowerProduct(ShowerProgenitorPtr progenitor,
ShowerParticlePtr parent,
const ShowerParticleVector & children) {
assert(children.size()==2);
bool matches[2];
for(unsigned int ix=0;ix<2;++ix) {
matches[ix] = children[ix]->id()==progenitor->id();
}
ShowerParticlePtr newpart;
if(matches[0]&&matches[1]) {
if(parent->showerKinematics()->z()>0.5) newpart=children[0];
else newpart=children[1];
}
else if(matches[0]) newpart=children[0];
else if(matches[1]) newpart=children[1];
_outgoingLines[progenitor]=newpart;
}
void ShowerTree::updateInitialStateShowerProduct(ShowerProgenitorPtr progenitor,
ShowerParticlePtr newParent) {
_incomingLines[progenitor]=newParent;
}
void ShowerTree::colourIsolate(const vector<PPtr> & original,
const vector<PPtr> & copy) {
// vectors must have same size
assert(original.size()==copy.size());
// create a temporary map with all the particles to make looping easier
vector<PPair> particles;
particles.reserve(original.size());
for(unsigned int ix=0;ix<original.size();++ix)
particles.push_back(make_pair(copy[ix],original[ix]));
// reset the colour of the copies
vector<PPair>::const_iterator cit,cjt;
for(cit=particles.begin();cit!=particles.end();++cit)
if((*cit).first->colourInfo()) (*cit).first->colourInfo(new_ptr(ColourBase()));
map<tColinePtr,tColinePtr> cmap;
// make the colour connections of the copies
for(cit=particles.begin();cit!=particles.end();++cit) {
ColinePtr c1,newline;
// if particle has a colour line
if((*cit).second->colourLine()&&!(*cit).first->colourLine()) {
c1=(*cit).second->colourLine();
newline=ColourLine::create((*cit).first);
cmap[c1]=newline;
for(cjt=particles.begin();cjt!=particles.end();++cjt) {
if(cjt==cit) continue;
if((*cjt).second->colourLine()==c1)
newline->addColoured((*cjt).first);
else if((*cjt).second->antiColourLine()==c1)
newline->addColoured((*cjt).first,true);
}
}
// if anticolour line
if((*cit).second->antiColourLine()&&!(*cit).first->antiColourLine()) {
c1=(*cit).second->antiColourLine();
newline=ColourLine::create((*cit).first,true);
cmap[c1]=newline;
for(cjt=particles.begin();cjt!=particles.end();++cjt) {
if(cjt==cit) continue;
if((*cjt).second->colourLine()==c1)
newline->addColoured((*cjt).first);
else if((*cjt).second->antiColourLine()==c1)
newline->addColoured((*cjt).first,true);
}
}
}
// sort out sinks and sources
for(cit=particles.begin();cit!=particles.end();++cit) {
tColinePtr cline[2];
tColinePair cpair;
for(unsigned int ix=0;ix<4;++ix) {
cline[0] = ix<2 ? cit->second->colourLine() : cit->second->antiColourLine();
cline[1] = ix<2 ? cit->first ->colourLine() : cit->first ->antiColourLine();
if(cline[0]) {
switch (ix) {
case 0: case 2:
cpair = cline[0]->sinkNeighbours();
break;
case 1: case 3:
cpair = cline[0]->sourceNeighbours();
break;
};
}
else {
cpair = make_pair(tColinePtr(),tColinePtr());
}
if(cline[0]&&cpair.first) {
map<tColinePtr,tColinePtr>::const_iterator
mit[2] = {cmap.find(cpair.first),cmap.find(cpair.second)};
if(mit[0]!=cmap.end()&&mit[1]!=cmap.end()) {
if(ix==0||ix==2) {
cline[1]->setSinkNeighbours(mit[0]->second,mit[1]->second);
}
else {
cline[1]->setSourceNeighbours(mit[0]->second,mit[1]->second);
}
}
}
}
}
}
void ShowerTree::insertHard(StepPtr pstep, bool ISR, bool) {
assert(_incomingLines.size()==2);
_colour.clear();
map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator cit;
// construct the map of colour lines for hard process
for(cit=_incomingLines.begin();cit!=_incomingLines.end();++cit) {
if(!cit->first->perturbative()) continue;
if((*cit).first->copy()->colourLine())
_colour.insert(make_pair((*cit).first->copy()->colourLine(),
(*cit).first->original()->colourLine()));
if((*cit).first->copy()->antiColourLine())
_colour.insert(make_pair((*cit).first->copy()->antiColourLine(),
(*cit).first->original()->antiColourLine()));
}
map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator cjt;
for(cjt=_outgoingLines.begin();cjt!=_outgoingLines.end();++cjt) {
if(!cjt->first->perturbative()) continue;
if((*cjt).first->copy()->colourLine())
_colour.insert(make_pair((*cjt).first->copy()->colourLine(),
(*cjt).first->original()->colourLine()));
if((*cjt).first->copy()->antiColourLine())
_colour.insert(make_pair((*cjt).first->copy()->antiColourLine(),
(*cjt).first->original()->antiColourLine()));
}
// initial-state radiation
if(ISR) {
for(cit=incomingLines().begin();cit!=incomingLines().end();++cit) {
ShowerParticlePtr init=(*cit).first->progenitor();
assert(init->getThePEGBase());
PPtr original = (*cit).first->original();
if(original->parents().empty()) continue;
PPtr hadron= original->parents()[0];
assert(!original->children().empty());
PPtr copy=cit->first->copy();
ParticleVector intermediates=original->children();
for(unsigned int ix=0;ix<intermediates.size();++ix) {
init->abandonChild(intermediates[ix]);
copy->abandonChild(intermediates[ix]);
}
// if not from a matrix element correction
if(cit->first->perturbative()) {
// break mother/daugther relations
init->addChild(original);
hadron->abandonChild(original);
// if particle showers add shower
if(cit->first->hasEmitted()) {
addInitialStateShower(init,hadron,pstep,false);
}
// no showering for this particle
else {
updateColour(init);
hadron->addChild(init);
pstep->addIntermediate(init);
}
}
// from matrix element correction
else {
// break mother/daugther relations
hadron->abandonChild(original);
copy->addChild(original);
updateColour(copy);
init->addChild(copy);
pstep->addIntermediate(copy);
// if particle showers add shower
if(cit->first->hasEmitted()) {
addInitialStateShower(init,hadron,pstep,false);
}
// no showering for this particle
else {
updateColour(init);
hadron->addChild(init);
pstep->addIntermediate(init);
}
}
}
}
else {
for(cit=incomingLines().begin();cit!=incomingLines().end();++cit) {
ShowerParticlePtr init=(*cit).first->progenitor();
assert(init->getThePEGBase());
PPtr original = (*cit).first->original();
if(original->parents().empty()) continue;
PPtr hadron= original->parents()[0];
assert(!original->children().empty());
PPtr copy=cit->first->copy();
ParticleVector intermediates=original->children();
for(unsigned int ix=0;ix<intermediates.size();++ix) {
init->abandonChild(intermediates[ix]);
copy->abandonChild(intermediates[ix]);
}
// break mother/daugther relations
init->addChild(original);
hadron->abandonChild(original);
// no showering for this particle
updateColour(init);
hadron->addChild(init);
pstep->addIntermediate(init);
}
}
// final-state radiation
for(cjt=outgoingLines().begin();cjt!=outgoingLines().end();++cjt) {
ShowerParticlePtr init=(*cjt).first->progenitor();
assert(init->getThePEGBase());
// if not from a matrix element correction
if(cjt->first->perturbative()) {
// register the shower particle as a
// copy of the one from the hard process
tParticleVector parents=init->parents();
for(unsigned int ix=0;ix<parents.size();++ix)
parents[ix]->abandonChild(init);
(*cjt).first->original()->addChild(init);
pstep->addDecayProduct(init);
}
// from a matrix element correction
else {
if(cjt->first->original()==_incoming.first||
cjt->first->original()==_incoming.second) {
updateColour((*cjt).first->copy());
(*cjt).first->original()->parents()[0]->
addChild((*cjt).first->copy());
pstep->addDecayProduct((*cjt).first->copy());
(*cjt).first->copy()->addChild(init);
pstep->addDecayProduct(init);
}
else {
updateColour((*cjt).first->copy());
(*cjt).first->original()->addChild((*cjt).first->copy());
pstep->addDecayProduct((*cjt).first->copy());
(*cjt).first->copy()->addChild(init);
pstep->addDecayProduct(init);
}
}
updateColour(init);
// insert shower products
addFinalStateShower(init,pstep);
}
_colour.clear();
}
void ShowerTree::addFinalStateShower(PPtr p, StepPtr s) {
if(p->children().empty()) return;
ParticleVector::const_iterator child;
for(child=p->children().begin(); child != p->children().end(); ++child) {
updateColour(*child);
s->addDecayProduct(*child);
addFinalStateShower(*child,s);
}
}
void ShowerTree::updateColour(PPtr particle) {
// if attached to a colour line
if(particle->colourLine()) {
bool reset=false;
// if colour line from hard process reconnect
if(_colour.find(particle->colourLine())!=_colour.end()) {
ColinePtr c1=particle->colourLine();
c1->removeColoured(particle);
_colour[c1]->addColoured(particle);
reset=true;
}
// ensure properly connected to the line
if(!reset) {
ColinePtr c1=particle->colourLine();
c1->removeColoured(particle);
c1->addColoured(particle);
}
}
// if attached to an anticolour line
if(particle->antiColourLine()) {
bool reset=false;
// if anti colour line from hard process reconnect
if(_colour.find(particle->antiColourLine())!=_colour.end()) {
ColinePtr c1=particle->antiColourLine();
c1->removeColoured(particle,true);
_colour[c1]->addColoured(particle,true);
reset=true;
}
if(!reset) {
ColinePtr c1=particle->antiColourLine();
c1->removeColoured(particle,true);
c1->addColoured(particle,true);
}
}
}
void ShowerTree::addInitialStateShower(PPtr p, PPtr hadron,
StepPtr s, bool addchildren) {
// Each parton here should only have one parent
if(!p->parents().empty()) {
if(p->parents().size()!=1)
throw Exception() << "Particle must only have one parent in ShowerTree"
<< "::addInitialStateShower" << Exception::runerror;
addInitialStateShower(p->parents()[0],hadron,s);
}
else {
hadron->addChild(p);
s->addIntermediate(p);
}
updateColour(p);
ParticleVector::const_iterator child;
// if not adding children return
if(!addchildren) return;
// add children
for(child = p->children().begin(); child != p->children().end(); ++child) {
// if a final-state particle update the colour
ShowerParticlePtr schild =
dynamic_ptr_cast<ShowerParticlePtr>(*child);
if(schild && schild->isFinalState()) updateColour(*child);
// if there are grandchildren of p
if(!(*child)->children().empty()) {
// Add child as intermediate
s->addIntermediate(*child);
// If child is shower particle and final-state, add children
if(schild && schild->isFinalState()) addFinalStateShower(schild,s);
}
else
s->addDecayProduct(*child);
}
}
void ShowerTree::decay(ShowerDecayMap & decay) {
// must be one incoming particle
assert(_incomingLines.size()==1);
// if already decayed return
if(!_outgoingLines.empty()) return;
// otherwise decay it
// now we need to replace the particle with a new copy after the shower
// find particle after the shower
ShowerParticlePtr newparent=_parent->_treelinks[this].second;
// now make the new progenitor
vector<PPtr> original,copy;
original.push_back(newparent);
copy.push_back(new_ptr(Particle(*newparent)));
// reisolate the colour
colourIsolate(original,copy);
// make the new progenitor
ShowerParticlePtr stemp=new_ptr(ShowerParticle(*copy[0],2,false));
fixColour(stemp);
ShowerProgenitorPtr newprog=new_ptr(ShowerProgenitor(original[0],copy[0],stemp));
_incomingLines.clear();
_incomingLines.insert(make_pair(newprog,stemp));
// now we need to decay the copy
PPtr parent=copy[0];
unsigned int ntry = 0;
while (true) {
// exit if fails
if (++ntry>=200)
throw Exception() << "Failed to perform decay in ShowerTree::decay()"
<< " after " << 200
<< " attempts for " << parent->PDGName()
<< Exception::eventerror;
// select decay mode
tDMPtr dm(parent->data().selectMode(*parent));
if(!dm)
throw Exception() << "Failed to select decay mode in ShowerTree::decay()"
<< "for " << newparent->PDGName()
<< Exception::eventerror;
if(!dm->decayer())
throw Exception() << "No Decayer for selected decay mode "
<< " in ShowerTree::decay()"
<< Exception::runerror;
// start of try block
try {
ParticleVector children = dm->decayer()->decay(*dm, *parent);
// if no children have another go
if(children.empty()) continue;
// set up parent
parent->decayMode(dm);
// add children
for (unsigned int i = 0, N = children.size(); i < N; ++i ) {
children[i]->setLabVertex(parent->labDecayVertex());
parent->addChild(children[i]);
parent->scale(ZERO);
}
// if succeeded break out of loop
break;
}
catch(KinematicsReconstructionVeto) {}
}
// insert the trees from the children
ParticleVector children=parent->children();
map<PPtr,ShowerTreePtr> trees;
for(unsigned int ix=0;ix<children.size();++ix) {
PPtr orig=children[ix];
parent->abandonChild(orig);
// if particle has children or decays in shower
if(!orig->children().empty()||
(decaysInShower(orig->id())&&!orig->dataPtr()->stable())) {
ShowerTreePtr newtree=new_ptr(ShowerTree(orig,decay));
trees.insert(make_pair(orig,newtree));
Energy width=orig->dataPtr()->generateWidth(orig->mass());
decay.insert(make_pair(width,newtree));
}
// now create the shower progenitors
PPtr ncopy=new_ptr(Particle(*orig));
//copy[0]->addChild(ncopy);
ShowerParticlePtr nshow=new_ptr(ShowerParticle(*ncopy,2,true));
fixColour(nshow);
ShowerProgenitorPtr prog=new_ptr(ShowerProgenitor(children[ix],
ncopy,nshow));
_outgoingLines.insert(make_pair(prog,nshow));
}
// set up the map of daughter trees
map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator mit;
for(mit=_outgoingLines.begin();mit!=_outgoingLines.end();++mit) {
map<PPtr,ShowerTreePtr>::const_iterator tit=trees.find(mit->first->original());
if(tit!=trees.end()) {
_treelinks.insert(make_pair(tit->second,
make_pair(mit->first,
mit->first->progenitor())));
tit->second->_parent=this;
}
}
}
void ShowerTree::insertDecay(StepPtr pstep,bool ISR, bool) {
assert(_incomingLines.size()==1);
_colour.clear();
// find final particle from previous tree
PPtr final;
if(_parent&&!_parent->_treelinks.empty())
final = _parent->_treelinks[this].second;
else
final=_incomingLines.begin()->first->original();
// construct the map of colour lines
PPtr copy=_incomingLines.begin()->first->copy();
if(copy->colourLine())
_colour.insert(make_pair(copy->colourLine(),final->colourLine()));
if(copy->antiColourLine())
_colour.insert(make_pair(copy->antiColourLine(),final->antiColourLine()));
// initial-state radiation
if(ISR&&!_incomingLines.begin()->first->progenitor()->children().empty()) {
ShowerParticlePtr init=_incomingLines.begin()->first->progenitor();
updateColour(init);
final->addChild(init);
pstep->addDecayProduct(init);
// insert shower products
addFinalStateShower(init,pstep);
// sort out colour
final=_incomingLines.begin()->second;
_colour.clear();
if(copy->colourLine())
_colour.insert(make_pair(copy->colourLine(),final->colourLine()));
if(copy->antiColourLine())
_colour.insert(make_pair(copy->antiColourLine(),final->antiColourLine()));
}
// get the decaying particles
// make the copy
tColinePair cline=make_pair(copy->colourLine(),copy->antiColourLine());
updateColour(copy);
// sort out sinks and sources if needed
if(cline.first) {
if(cline.first->sourceNeighbours().first) {
copy->colourLine()->setSourceNeighbours(cline.first->sourceNeighbours().first,
cline.first->sourceNeighbours().second);
}
else if (cline.first->sinkNeighbours().first) {
copy->colourLine()->setSinkNeighbours(cline.first->sinkNeighbours().first,
cline.first->sinkNeighbours().second);
}
}
if(cline.second) {
if(cline.second->sourceNeighbours().first) {
copy->antiColourLine()->setSourceNeighbours(cline.second->sourceNeighbours().first,
cline.second->sourceNeighbours().second);
}
else if (cline.second->sinkNeighbours().first) {
copy->antiColourLine()->setSinkNeighbours(cline.second->sinkNeighbours().first,
cline.second->sinkNeighbours().second);
}
}
// copy of the one from the hard process
tParticleVector dpar=copy->parents();
for(unsigned int ix=0;ix<dpar.size();++ix) dpar[ix]->abandonChild(copy);
final->addChild(copy);
pstep->addDecayProduct(copy);
// final-state radiation
map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator cit;
for(cit=outgoingLines().begin();cit!=outgoingLines().end();++cit) {
ShowerParticlePtr init=cit->first->progenitor();
if(!init->getThePEGBase())
throw Exception() << "Final-state particle must have a ThePEGBase"
<< " in ShowerTree::fillEventRecord()"
<< Exception::runerror;
// if not from matrix element correction
if(cit->first->perturbative()) {
// add the child
updateColour(cit->first->copy());
PPtr orig=cit->first->original();
copy->addChild(orig);
pstep->addDecayProduct(orig);
orig->addChild(cit->first->copy());
pstep->addDecayProduct(cit->first->copy());
// register the shower particle as a
// copy of the one from the hard process
tParticleVector parents=init->parents();
for(unsigned int ix=0;ix<parents.size();++ix)
{parents[ix]->abandonChild(init);}
(*cit).first->copy()->addChild(init);
pstep->addDecayProduct(init);
updateColour(init);
}
// from a matrix element correction
else {
if(copy->children().end()==
find(copy->children().begin(),copy->children().end(),
cit->first->original())) {
updateColour(cit->first->original());
copy->addChild(cit->first->original());
pstep->addDecayProduct(cit->first->original());
}
updateColour(cit->first->copy());
cit->first->original()->addChild(cit->first->copy());
pstep->addDecayProduct(cit->first->copy());
// register the shower particle as a
// copy of the one from the hard process
tParticleVector parents=init->parents();
for(unsigned int ix=0;ix<parents.size();++ix)
{parents[ix]->abandonChild(init);}
(*cit).first->copy()->addChild(init);
pstep->addDecayProduct(init);
updateColour(init);
}
// insert shower products
addFinalStateShower(init,pstep);
}
_colour.clear();
}
void ShowerTree::clear() {
// reset the has showered flag
_hasShowered=false;
// clear the colour map
_colour.clear();
map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator cit;
map<ShowerProgenitorPtr, ShowerParticlePtr>::const_iterator cjt;
// abandon the children of the outgoing particles
for(cit=_outgoingLines.begin();cit!=_outgoingLines.end();++cit) {
ShowerParticlePtr orig=cit->first->progenitor();
orig->set5Momentum(cit->first->copy()->momentum());
ParticleVector children=orig->children();
for(unsigned int ix=0;ix<children.size();++ix) orig->abandonChild(children[ix]);
_outgoingLines[cit->first]=orig;
cit->first->hasEmitted(false);
}
// forward products
_forward.clear();
for(cit=_outgoingLines.begin();cit!=_outgoingLines.end();++cit)
_forward.insert(cit->first->progenitor());
// if a decay
if(!_wasHard) {
ShowerParticlePtr orig=_incomingLines.begin()->first->progenitor();
orig->set5Momentum(_incomingLines.begin()->first->copy()->momentum());
ParticleVector children=orig->children();
for(unsigned int ix=0;ix<children.size();++ix) orig->abandonChild(children[ix]);
}
// if a hard process
else {
for(cjt=_incomingLines.begin();cjt!=_incomingLines.end();++cjt) {
tPPtr parent = cjt->first->original()->parents().empty() ?
tPPtr() : cjt->first->original()->parents()[0];
ShowerParticlePtr temp=
new_ptr(ShowerParticle(*cjt->first->copy(),
cjt->first->progenitor()->perturbative(),
cjt->first->progenitor()->isFinalState()));
fixColour(temp);
temp->x(cjt->first->progenitor()->x());
cjt->first->hasEmitted(false);
if(!(cjt->first->progenitor()==cjt->second)&&cjt->second&&parent)
parent->abandonChild(cjt->second);
cjt->first->progenitor(temp);
_incomingLines[cjt->first]=temp;
}
}
// reset the particles at the end of the shower
_backward.clear();
// if hard process backward products
if(_wasHard)
for(cjt=_incomingLines.begin();cjt!=_incomingLines.end();++cjt)
_backward.insert(cjt->first->progenitor());
clearTransforms();
}
void ShowerTree::resetShowerProducts() {
map<ShowerProgenitorPtr, ShowerParticlePtr>::const_iterator cit;
map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator cjt;
_backward.clear();
_forward.clear();
for(cit=_incomingLines.begin();cit!=_incomingLines.end();++cit)
_backward.insert(cit->second);
for(cjt=_outgoingLines.begin();cjt!=_outgoingLines.end();++cjt)
_forward.insert(cjt->second);
}
void ShowerTree::updateAfterShower(ShowerDecayMap & decay) {
// update the links
map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator mit;
map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::iterator tit;
for(tit=_treelinks.begin();tit!=_treelinks.end();++tit) {
if(tit->second.first) {
mit=_outgoingLines.find(tit->second.first);
if(mit!=_outgoingLines.end()) tit->second.second=mit->second;
}
}
// get the particles coming from those in the hard process
set<tShowerParticlePtr> hard;
for(mit=_outgoingLines.begin();mit!=_outgoingLines.end();++mit)
hard.insert(mit->second);
// find the shower particles which should be decayed in the
// shower but didn't come from the hard process
set<tShowerParticlePtr>::const_iterator cit;
for(cit=_forward.begin();cit!=_forward.end();++cit) {
if(decaysInShower((**cit).id())&&
hard.find(*cit)==hard.end()) {
ShowerTreePtr newtree=new_ptr(ShowerTree(*cit,decay));
newtree->setParents();
newtree->_parent=this;
Energy width=(**cit).dataPtr()->generateWidth((**cit).mass());
decay.insert(make_pair(width,newtree));
_treelinks.insert(make_pair(newtree,
make_pair(tShowerProgenitorPtr(),*cit)));
}
}
}
void ShowerTree::addFinalStateBranching(ShowerParticlePtr parent,
- const ShowerParticleVector & children) {
+ const ShowerParticleVector & children) {
assert(children.size()==2);
_forward.erase(parent);
for(unsigned int ix=0; ix<children.size(); ++ix) {
_forward.insert(children[ix]);
}
}
void ShowerTree::addInitialStateBranching(ShowerParticlePtr oldParent,
ShowerParticlePtr newParent,
ShowerParticlePtr otherChild) {
_backward.erase(oldParent);
_backward.insert(newParent);
_forward.insert(otherChild);
}
void ShowerTree::setParents() {
// set the parent tree of the children
map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator tit;
for(tit=_treelinks.begin();tit!=_treelinks.end();++tit)
tit->first->_parent=this;
}
void ShowerTree::fixColour(tShowerParticlePtr part) {
ColinePtr line=part->colourLine();
if(line) {
line->removeColoured(part);
line->addColoured(part);
}
line=part->antiColourLine();
if(line) {
line->removeAntiColoured(part);
line->addAntiColoured(part);
}
}
vector<ShowerProgenitorPtr> ShowerTree::extractProgenitors() {
// extract the particles from the ShowerTree
map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator mit;
vector<ShowerProgenitorPtr> ShowerHardJets;
for(mit=incomingLines().begin();mit!=incomingLines().end();++mit)
ShowerHardJets.push_back((*mit).first);
map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator mjt;
for(mjt=outgoingLines().begin();mjt!=outgoingLines().end();++mjt)
ShowerHardJets.push_back((*mjt).first);
return ShowerHardJets;
}
void ShowerTree::transform(const LorentzRotation & boost, bool applyNow) {
if(applyNow) {
// now boost all the particles
map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator cit;
// incoming
for(cit=_incomingLines.begin();cit!=_incomingLines.end();++cit) {
cit->first->progenitor()->deepTransform(boost);
cit->first->copy()->deepTransform(boost);
}
// outgoing
map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator cjt;
for(cjt=_outgoingLines.begin();cjt!=_outgoingLines.end();++cjt) {
cjt->first->progenitor()->deepTransform(boost);
cjt->first->copy()->deepTransform(boost);
}
}
else {
Lorentz5Momentum ptemp1 = _incomingLines.begin()->first->progenitor()->momentum();
Lorentz5Momentum ptemp2 = ptemp1;
ptemp1 *= _transforms;
ptemp1 *= boost;
_transforms.transform(boost);
ptemp2 *= _transforms;
}
// child trees
for(map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator
tit=_treelinks.begin();tit!=_treelinks.end();++tit)
tit->first->transform(boost,applyNow);
}
void ShowerTree::applyTransforms() {
// now boost all the particles
map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator cit;
// incoming
for(cit=_incomingLines.begin();cit!=_incomingLines.end();++cit) {
cit->first->progenitor()->deepTransform(_transforms);
cit->first->copy()->deepTransform(_transforms);
}
// outgoing
map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator cjt;
for(cjt=_outgoingLines.begin();cjt!=_outgoingLines.end();++cjt) {
cjt->first->progenitor()->deepTransform(_transforms);
cjt->first->copy()->deepTransform(_transforms);
}
// child trees
for(map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator
tit=_treelinks.begin();tit!=_treelinks.end();++tit)
tit->first->applyTransforms();
_transforms = LorentzRotation();
}
void ShowerTree::clearTransforms() {
_transforms = LorentzRotation();
// child trees
for(map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator
tit=_treelinks.begin();tit!=_treelinks.end();++tit)
tit->first->clearTransforms();
}
diff --git a/Shower/Base/SudakovFormFactor.cc b/Shower/Base/SudakovFormFactor.cc
--- a/Shower/Base/SudakovFormFactor.cc
+++ b/Shower/Base/SudakovFormFactor.cc
@@ -1,281 +1,309 @@
// -*- C++ -*-
//
// SudakovFormFactor.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the SudakovFormFactor class.
//
#include "SudakovFormFactor.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ShowerKinematics.h"
using namespace Herwig;
void SudakovFormFactor::persistentOutput(PersistentOStream & os) const {
os << splittingFn_ << alpha_ << pdfmax_ << particles_ << pdffactor_
<< a_ << b_ << ounit(c_,GeV) << ounit(kinCutoffScale_,GeV) << cutOffOption_
<< ounit(vgcut_,GeV) << ounit(vqcut_,GeV)
<< ounit(pTmin_,GeV) << ounit(pT2min_,GeV2);
}
void SudakovFormFactor::persistentInput(PersistentIStream & is, int) {
is >> splittingFn_ >> alpha_ >> pdfmax_ >> particles_ >> pdffactor_
>> a_ >> b_ >> iunit(c_,GeV) >> iunit(kinCutoffScale_,GeV) >> cutOffOption_
>> iunit(vgcut_,GeV) >> iunit(vqcut_,GeV)
>> iunit(pTmin_,GeV) >> iunit(pT2min_,GeV2);
}
AbstractClassDescription<SudakovFormFactor> SudakovFormFactor::initSudakovFormFactor;
// Definition of the static class description member.
void SudakovFormFactor::Init() {
static ClassDocumentation<SudakovFormFactor> documentation
("The SudakovFormFactor class is the base class for the implementation of Sudakov"
" form factors in Herwig++");
static Reference<SudakovFormFactor,SplittingFunction>
interfaceSplittingFunction("SplittingFunction",
"A reference to the SplittingFunction object",
&Herwig::SudakovFormFactor::splittingFn_,
false, false, true, false);
static Reference<SudakovFormFactor,ShowerAlpha>
interfaceAlpha("Alpha",
"A reference to the Alpha object",
&Herwig::SudakovFormFactor::alpha_,
false, false, true, false);
static Parameter<SudakovFormFactor,double> interfacePDFmax
("PDFmax",
"Maximum value of PDF weight. ",
&SudakovFormFactor::pdfmax_, 35.0, 1.0, 100000.0,
false, false, Interface::limited);
static Switch<SudakovFormFactor,unsigned int> interfacePDFFactor
("PDFFactor",
"Include additional factors in the overestimate for the PDFs",
&SudakovFormFactor::pdffactor_, 0, false, false);
static SwitchOption interfacePDFFactorOff
(interfacePDFFactor,
"Off",
"Don't include any factors",
0);
static SwitchOption interfacePDFFactorOverZ
(interfacePDFFactor,
"OverZ",
"Include an additional factor of 1/z",
1);
static SwitchOption interfacePDFFactorOverOneMinusZ
(interfacePDFFactor,
"OverOneMinusZ",
"Include an additional factor of 1/(1-z)",
2);
static SwitchOption interfacePDFFactorOverZOneMinusZ
(interfacePDFFactor,
"OverZOneMinusZ",
"Include an additional factor of 1/z/(1-z)",
3);
static Switch<SudakovFormFactor,unsigned int> interfaceCutOffOption
("CutOffOption",
"The type of cut-off to use to end the shower",
&SudakovFormFactor::cutOffOption_, 0, false, false);
static SwitchOption interfaceCutOffOptionDefault
(interfaceCutOffOption,
"Default",
"Use the standard Herwig++ cut-off on virtualities with the minimum"
" virtuality depending on the mass of the branching particle",
0);
static SwitchOption interfaceCutOffOptionFORTRAN
(interfaceCutOffOption,
"FORTRAN",
"Use a FORTRAN-like cut-off on virtualities",
1);
static SwitchOption interfaceCutOffOptionpT
(interfaceCutOffOption,
"pT",
"Use a cut on the minimum allowed pT",
2);
static Parameter<SudakovFormFactor,double> interfaceaParameter
("aParameter",
"The a parameter for the kinematic cut-off",
&SudakovFormFactor::a_, 0.3, -10.0, 10.0,
false, false, Interface::limited);
static Parameter<SudakovFormFactor,double> interfacebParameter
("bParameter",
"The b parameter for the kinematic cut-off",
&SudakovFormFactor::b_, 2.3, -10.0, 10.0,
false, false, Interface::limited);
static Parameter<SudakovFormFactor,Energy> interfacecParameter
("cParameter",
"The c parameter for the kinematic cut-off",
&SudakovFormFactor::c_, GeV, 0.3*GeV, 0.1*GeV, 10.0*GeV,
false, false, Interface::limited);
static Parameter<SudakovFormFactor,Energy>
interfaceKinScale ("cutoffKinScale",
"kinematic cutoff scale for the parton shower phase"
" space (unit [GeV])",
&SudakovFormFactor::kinCutoffScale_, GeV,
2.3*GeV, 0.001*GeV, 10.0*GeV,false,false,false);
static Parameter<SudakovFormFactor,Energy> interfaceGluonVirtualityCut
("GluonVirtualityCut",
"For the FORTRAN cut-off option the minimum virtuality of the gluon",
&SudakovFormFactor::vgcut_, GeV, 0.85*GeV, 0.1*GeV, 10.0*GeV,
false, false, Interface::limited);
static Parameter<SudakovFormFactor,Energy> interfaceQuarkVirtualityCut
("QuarkVirtualityCut",
"For the FORTRAN cut-off option the minimum virtuality added to"
" the mass for particles other than the gluon",
&SudakovFormFactor::vqcut_, GeV, 0.85*GeV, 0.1*GeV, 10.0*GeV,
false, false, Interface::limited);
static Parameter<SudakovFormFactor,Energy> interfacepTmin
("pTmin",
"The minimum pT if using a cut-off on the pT",
&SudakovFormFactor::pTmin_, GeV, 1.0*GeV, ZERO, 10.0*GeV,
false, false, Interface::limited);
}
bool SudakovFormFactor::
PDFVeto(const Energy2 t, const double x,
const tcPDPtr parton0, const tcPDPtr parton1,
Ptr<BeamParticleData>::transient_const_pointer beam) const {
assert(pdf_);
Energy2 theScale = t;
if (theScale < sqr(freeze_)) theScale = sqr(freeze_);
double newpdf(0.0), oldpdf(0.0);
//different treatment of MPI ISR is done via CascadeHandler::resetPDFs()
newpdf=pdf_->xfx(beam,parton0,theScale,x/z());
oldpdf=pdf_->xfx(beam,parton1,theScale,x);
if(newpdf<=0.) return true;
if(oldpdf<=0.) return false;
double ratio = newpdf/oldpdf;
double maxpdf(pdfmax_);
switch (pdffactor_) {
case 1:
maxpdf /= z();
break;
case 2:
maxpdf /= 1.-z();
break;
case 3:
maxpdf /= (z()*(1.-z()));
break;
}
// ratio / PDFMax must be a probability <= 1.0
if (ratio > maxpdf) {
generator()->log() << "PDFVeto warning: Ratio > " << name()
<< ":PDFmax (by a factor of "
<< ratio/maxpdf <<") for "
<< parton0->PDGName() << " to "
<< parton1->PDGName() << "\n";
}
return ratio < UseRandom::rnd()*maxpdf;
}
void SudakovFormFactor::addSplitting(const IdList & in) {
bool add=true;
for(unsigned int ix=0;ix<particles_.size();++ix) {
if(particles_[ix].size()==in.size()) {
bool match=true;
for(unsigned int iy=0;iy<in.size();++iy) {
if(particles_[ix][iy]!=in[iy]) {
match=false;
break;
}
}
if(match) {
add=false;
break;
}
}
}
if(add) particles_.push_back(in);
}
void SudakovFormFactor::removeSplitting(const IdList & in) {
for(vector<IdList>::iterator it=particles_.begin();
it!=particles_.end();++it) {
if(it->size()==in.size()) {
bool match=true;
for(unsigned int iy=0;iy<in.size();++iy) {
if((*it)[iy]!=in[iy]) {
match=false;
break;
}
}
if(match) {
vector<IdList>::iterator itemp=it;
--itemp;
particles_.erase(it);
it = itemp;
}
}
}
}
Energy2 SudakovFormFactor::guesst(Energy2 t1,unsigned int iopt,
const IdList &ids,
double enhance,bool ident) const {
unsigned int pdfopt = iopt!=1 ? 0 : pdffactor_;
double c =
1./((splittingFn_->integOverP(zlimits_.second,ids,pdfopt) -
splittingFn_->integOverP(zlimits_.first ,ids,pdfopt))*
alpha_->overestimateValue()/Constants::twopi*enhance);
assert(iopt<=2);
if(iopt==1) {
c/=pdfmax_;
if(ident) c*=0.5;
}
else if(iopt==2) c*=-1.;
if(splittingFn_->interactionOrder()==1) {
return t1*pow(UseRandom::rnd(),c);
}
else {
assert(false && "Units are dubious here.");
int nm(splittingFn()->interactionOrder()-1);
c/=Math::powi(alpha_->overestimateValue()/Constants::twopi,nm);
return t1 / pow (1. - nm*c*log(UseRandom::rnd())
* Math::powi(t1*UnitRemoval::InvE2,nm)
,1./double(nm));
}
}
double SudakovFormFactor::guessz (unsigned int iopt, const IdList &ids) const {
unsigned int pdfopt = iopt!=1 ? 0 : pdffactor_;
double lower = splittingFn_->integOverP(zlimits_.first,ids,pdfopt);
return splittingFn_->invIntegOverP
(lower + UseRandom::rnd()*(splittingFn_->integOverP(zlimits_.second,ids,pdfopt) -
lower),ids,pdfopt);
}
void SudakovFormFactor::doinit() {
Interfaced::doinit();
pT2min_ = cutOffOption()==2 ? sqr(pTmin_) : ZERO;
}
+
+vector<Energy> SudakovFormFactor::virtualMasses(const IdList & ids) {
+ vector<Energy> output;
+ if(cutOffOption() == 0) {
+ for(unsigned int ix=0;ix<ids.size();++ix)
+ output.push_back(getParticleData(ids[ix])->mass());
+ Energy kinCutoff=
+ kinematicCutOff(kinScale(),*std::max_element(output.begin(),output.end()));
+ for(unsigned int ix=0;ix<output.size();++ix)
+ output[ix]=max(kinCutoff,output[ix]);
+ }
+ else if(cutOffOption() == 1) {
+ for(unsigned int ix=0;ix<ids.size();++ix) {
+ output.push_back(getParticleData(ids[ix])->mass());
+ output.back() += ids[ix]==ParticleID::g ? vgCut() : vqCut();
+ }
+ }
+ else if(cutOffOption() == 2) {
+ for(unsigned int ix=0;ix<ids.size();++ix)
+ output.push_back(getParticleData(ids[ix])->mass());
+ }
+ else {
+ throw Exception() << "Unknown option for the cut-off"
+ << " in SudakovFormFactor::virtualMasses()"
+ << Exception::runerror;
+ }
+ return output;
+}
diff --git a/Shower/Base/SudakovFormFactor.h b/Shower/Base/SudakovFormFactor.h
--- a/Shower/Base/SudakovFormFactor.h
+++ b/Shower/Base/SudakovFormFactor.h
@@ -1,679 +1,684 @@
// -*- C++ -*-
//
// SudakovFormFactor.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_SudakovFormFactor_H
#define HERWIG_SudakovFormFactor_H
//
// This is the declaration of the SudakovFormFactor class.
//
#include "ThePEG/Interface/Interfaced.h"
#include "Herwig++/Shower/SplittingFunctions/SplittingFunction.h"
#include "Herwig++/Shower/Couplings/ShowerAlpha.h"
#include "Herwig++/Shower/SplittingFunctions/SplittingGenerator.fh"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/PDF/BeamParticleData.h"
#include <cassert>
#include "ShowerKinematics.fh"
#include "SudakovFormFactor.fh"
namespace Herwig {
using namespace ThePEG;
/**
* A typedef for the BeamParticleData
*/
typedef Ptr<BeamParticleData>::transient_const_pointer tcBeamPtr;
/** \ingroup Shower
*
* This is the definition of the Sudakov form factor class. In general this
* is the base class for the implementation of Sudakov form factors in Herwig++.
* The methods generateNextTimeBranching(), generateNextDecayBranching() and
* generateNextSpaceBranching need to be implemented in classes inheriting from this
* one.
*
* In addition a number of methods are implemented to assist with the calculation
* of the form factor using the veto algorithm in classes inheriting from this one.
*
* In general the Sudakov form-factor, for final-state radiation, is given
* by
* \f[\Delta_{ba}(\tilde{q}_{i+1},\tilde{q}_i)=
* \exp\left\{
* -\int^{\tilde{q}^2_i}_{\tilde{q}^2_{i+1}}
* \frac{{\rm d}\tilde{q}^2}{\tilde{q}^2}
* \int\frac{\alpha_S(z,\tilde{q})}{2\pi}
* P_{ba}(z,\tilde{q})\Theta(p_T)
* \right\}.
* \f]
* We can solve this to obtain the next value of the scale \f$\tilde{q}_{i+1}\f$
* given the previous value \f$\tilde{q}_i\f$
* in the following way. First we obtain a simplified form of the integrand
* which is greater than or equal to the true integrand for all values of
* \f$\tilde{q}\f$.
*
* In practice it is easiest to obtain this over estimate in pieces. The ShowerAlpha
* object contains an over estimate for \f$\alpha_S\f$, the splitting function
* contains both an over estimate of the spltting function and its integral
* which is needed to compute the over estimate of the \f$\tilde{q}\f$ integrand,
* together with an over estimate of the limit of the \f$z\f$ integral.
*
* This gives an overestimate of the integrand
* \f[g(\tilde{q}^2) = \frac{c}{\tilde{q}^2}, \f]
* where because the over estimates are chosen to be independent of \f$\tilde{q}\f$ the
* parameter
* \f[c = \frac{\alpha_{\rm over}}{2\pi}\int^{z_1}_{z_0}P_{\rm over}(z),\f]
* is a constant independent of \f$\tilde{q}\f$.
*
* The guesst() member can then be used to generate generate the value of
* \f$\tilde{q}^2\f$ according to this result. This is done by solving the Sudakov
* form factor, with the over estimates, is equal to a random number
* \f$r\f$ in the interval \f$[0,1]\f$. This gives
* \f[\tilde{q}^2_{i+1}=G^{-1}\left[G(\tilde{q}^2_i)+\ln r\right],\f]
* where \f$G(\tilde{q}^2)=c\ln(\tilde{q}^2)\f$ is the infinite integral
* of \f$g(\tilde{q}^2)\f$ and \f$G^{-1}(x)=\exp\left(\frac{x}c\right)\f$
* is its inverse.
* It this case we therefore obtain
* \f[\tilde{q}^2_{i+1}=\tilde{q}^2_ir^{\frac1c}.\f]
* The value of \f$z\f$ can then be calculated in a similar way
* \f[z = I^{-1}\left[I(z_0)+r\left(I(z_1)-I(z_0)\right)\right],\f]
* using the guessz() member,
* where \f$I=\int P(z){\rm d}z\f$ and \f$I^{-1}\f$ is its inverse.
*
* The veto algorithm then uses rejection using the ratio of the
* true value to the overestimated one to obtain the original distribution.
* This is accomplished using the
* - alphaSVeto() member for the \f$\alpha_S\f$ veto
* - SplittingFnVeto() member for the veto on the value of the splitting function.
* in general there must also be a chech that the emission is in the allowed
* phase space but this is left to the inheriting classes as it will depend
* on the ordering variable.
*
* The Sudakov form factor for the initial-scale shower is different because
* it must include the PDF which guides the backward evolution.
* It is given by
* \f[\Delta_{ba}(\tilde{q}_{i+1},\tilde{q}_i)=
* \exp\left\{
* -\int^{\tilde{q}^2_i}_{\tilde{q}^2_{i+1}}
* \frac{{\rm d}\tilde{q}^2}{\tilde{q}^2}
* \int\frac{\alpha_S(z,\tilde{q})}{2\pi}
* P_{ba}(z,\tilde{q})\frac{x'f_a(\frac{x}z,\tilde{q}^2)}{xf_b(x,\tilde{q^2})}
* \right\},
* \f]
* where \f$x\f$ is the fraction of the beam momentum the parton \f$b\f$ had before
* the backward evolution.
* This can be solve in the same way as for the final-state branching but the constant
* becomes
* \f[c = \frac{\alpha_{\rm over}}{2\pi}\int^{z_1}_{z_0}P_{\rm over}(z)PDF_{\rm max},\f]
* where
* \f[PDF_{\rm max}=\max\frac{x'f_a(\frac{x}z,\tilde{q}^2)}{xf_b(x,\tilde{q^2})},\f]
* which can be set using an interface.
* In addition the PDFVeto() member then is needed to implement the relevant veto.
*
* @see SplittingGenerator
* @see SplittingFunction
* @see ShowerAlpha
* @see \ref SudakovFormFactorInterfaces "The interfaces"
* defined for SudakovFormFactor.
*/
class SudakovFormFactor: public Interfaced {
/**
* The SplittingGenerator is a friend to insert the particles in the
* branchings at initialisation
*/
friend class SplittingGenerator;
public:
/**
* The default constructor.
*/
SudakovFormFactor() : pdfmax_(35.0), pdffactor_(0),
cutOffOption_(0), a_(0.3), b_(2.3), c_(0.3*GeV),
kinCutoffScale_( 2.3*GeV ), vgcut_(0.85*GeV),
vqcut_(0.85*GeV), pTmin_(1.*GeV), pT2min_(ZERO),
z_( 0.0 ),phi_(0.0), pT_() {}
/**
* Members to generate the scale of the next branching
*/
//@{
/**
* Return the scale of the next time-like branching. If there is no
* branching then it returns ZERO.
* @param startingScale starting scale for the evolution
* @param ids The PDG codes of the particles in the splitting
* @param cc Whether this is the charge conjugate of the branching
* @param enhance The radiation enhancement factor
* defined.
*/
virtual ShoKinPtr generateNextTimeBranching(const Energy startingScale,
const IdList &ids,const bool cc,
double enhance)=0;
/**
* Return the scale of the next space-like decay branching. If there is no
* branching then it returns ZERO.
* @param startingScale starting scale for the evolution
* @param stoppingScale stopping scale for the evolution
* @param minmass The minimum mass allowed for the spake-like particle.
* @param ids The PDG codes of the particles in the splitting
* @param cc Whether this is the charge conjugate of the branching
* defined.
* @param enhance The radiation enhancement factor
*/
virtual ShoKinPtr generateNextDecayBranching(const Energy startingScale,
const Energy stoppingScale,
const Energy minmass,
const IdList &ids,
const bool cc,
double enhance)=0;
/**
* Return the scale of the next space-like branching. If there is no
* branching then it returns ZERO.
* @param startingScale starting scale for the evolution
* @param ids The PDG codes of the particles in the splitting
* @param x The fraction of the beam momentum
* @param cc Whether this is the charge conjugate of the branching
* defined.
* @param beam The beam particle
* @param enhance The radiation enhancement factor
*/
virtual ShoKinPtr generateNextSpaceBranching(const Energy startingScale,
const IdList &ids,double x,
const bool cc,double enhance,
tcBeamPtr beam)=0;
//@}
/**
* Methods to provide public access to the private member variables
*/
//@{
/**
* Return the pointer to the SplittingFunction object.
*/
tSplittingFnPtr splittingFn() const { return splittingFn_; }
/**
* Return the pointer to the ShowerAlpha object.
*/
tShowerAlphaPtr alpha() const { return alpha_; }
/**
* The type of interaction
*/
inline ShowerInteraction::Type interactionType() const
{return splittingFn_->interactionType();}
//@}
public:
/**
* Methods to access the kinematic variables for the branching
*/
//@{
/**
* The energy fraction
*/
double z() const { return z_; }
/**
* The azimuthal angle
*/
double phi() const { return phi_; }
/**
* The transverse momentum
*/
Energy pT() const { return pT_; }
//@}
/**
* Access the maximum weight for the PDF veto
*/
double pdfMax() const { return pdfmax_;}
/**
* Method to return the evolution scale given the
* transverse momentum, \f$p_T\f$ and \f$z\f$.
*/
virtual Energy calculateScale(double z, Energy pt, IdList ids,unsigned int iopt)=0;
/**
* Method to create the ShowerKinematics object for a final-state branching
*/
virtual ShoKinPtr createFinalStateBranching(Energy scale,double z,
double phi, Energy pt)=0;
/**
* Method to create the ShowerKinematics object for an initial-state branching
*/
virtual ShoKinPtr createInitialStateBranching(Energy scale,double z,
double phi, Energy pt)=0;
/**
* Method to create the ShowerKinematics object for a decay branching
*/
virtual ShoKinPtr createDecayBranching(Energy scale,double z,
double phi, Energy pt)=0;
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
//@}
protected:
/**
* Methods to implement the veto algorithm to generate the scale of
* the next branching
*/
//@{
/**
* Value of the energy fraction for the veto algorithm
* @param iopt The option for calculating z
* @param ids The PDG codes of the particles in the splitting
* - 0 is final-state
* - 1 is initial-state for the hard process
* - 2 is initial-state for particle decays
*/
double guessz (unsigned int iopt, const IdList &ids) const;
/**
* Value of the scale for the veto algorithm
* @param t1 The starting valoe of the scale
* @param iopt The option for calculating t
* @param ids The PDG codes of the particles in the splitting
* - 0 is final-state
* - 1 is initial-state for the hard process
* - 2 is initial-state for particle decays
* @param enhance The radiation enhancement factor
* @param identical Whether or not the outgoing particles are identical
*/
Energy2 guesst (Energy2 t1,unsigned int iopt, const IdList &ids,
double enhance, bool identical) const;
/**
* Veto on the PDF for the initial-state shower
* @param t The scale
* @param x The fraction of the beam momentum
* @param parton0 Pointer to the particleData for the
* new parent (this is the particle we evolved back to)
* @param parton1 Pointer to the particleData for the
* original particle
* @param beam The BeamParticleData object
*/
bool PDFVeto(const Energy2 t, const double x,
const tcPDPtr parton0, const tcPDPtr parton1,
tcBeamPtr beam) const;
/**
* The veto on the splitting function.
* @param t The scale
* @param ids The PDG codes of the particles in the splitting
* @param mass Whether or not to use the massive splitting functions
* @return true if vetoed
*/
bool SplittingFnVeto(const Energy2 t,
const IdList &ids,
const bool mass) const
{ return UseRandom::rnd()>splittingFn_->ratioP(z_, t, ids,mass); }
/**
* The veto on the coupling constant
* @param pt2 The value of ther transverse momentum squared, \f$p_T^2\f$.
* @return true if vetoed
*/
bool alphaSVeto(const Energy2 pt2) const
{return UseRandom::rnd() > ThePEG::Math::powi(alpha_->ratio(pt2),
splittingFn_->interactionOrder());}
//@}
/**
* Methods to set the kinematic variables for the branching
*/
//@{
/**
* The energy fraction
*/
void z(double in) { z_=in; }
/**
* The azimuthal angle
*/
void phi(double in) { phi_=in; }
/**
* The transverse momentum
*/
void pT(Energy in) { pT_=in; }
//@}
/**
* Set/Get the limits on the energy fraction for the splitting
*/
//@{
/**
* Get the limits
*/
pair<double,double> zLimits() const { return zlimits_;}
/**
* Set the limits
*/
void zLimits(pair<double,double> in) { zlimits_=in; }
//@}
/**
* Set the particles in the splittings
*/
void addSplitting(const IdList &);
/**
* Delete the particles in the splittings
*/
void removeSplitting(const IdList &);
/**
* Access the potential branchings
*/
vector<IdList> particles() const { return particles_; }
/**
* Methods to set the member variables for inheriting classes
*/
//@{
/**
* Method to set the SplittingFunction
*/
void splittingFn(tSplittingFnPtr in) { splittingFn_ = in;}
/**
* Method to set the coupling
*/
void alpha(tShowerAlphaPtr in) { alpha_ = in; }
/**
* Method to set the maximum PDF weight
*/
void pdfMax(double in) { pdfmax_ = in;}
/**
* Get the option for the PDF factor
*/
unsigned int PDFFactor() const { return pdffactor_; }
//@}
public:
/**
* @name Methods for the cut-off
*/
//@{
/**
* The option being used
*/
unsigned int cutOffOption() const { return cutOffOption_; }
/**
* The kinematic scale
*/
Energy kinScale() const {return kinCutoffScale_;}
/**
* The virtuality cut-off on the gluon \f$Q_g=\frac{\delta-am_q}{b}\f$
* @param scale The scale \f$\delta\f$
* @param mq The quark mass \f$m_q\f$.
*/
Energy kinematicCutOff(Energy scale, Energy mq) const
{return max((scale -a_*mq)/b_,c_);}
/**
* The virtualilty cut-off for gluons
*/
Energy vgCut() const { return vgcut_; }
/**
* The virtuality cut-off for everything else
*/
Energy vqCut() const { return vqcut_; }
/**
* The minimum \f$p_T\f$ for the branching
*/
Energy pTmin() const { return pTmin_; }
/**
* The square of the minimum \f$p_T\f$
*/
Energy2 pT2min() const { return pT2min_; }
+
+ /**
+ * Calculate the virtual masses for a branchings
+ */
+ vector<Energy> virtualMasses(const IdList & ids);
//@}
/**
* Set the PDF
*/
void setPDF(tcPDFPtr pdf, Energy scale) {
pdf_ = pdf;
freeze_ = scale;
}
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is an abstract class with persistent data.
*/
static AbstractClassDescription<SudakovFormFactor> initSudakovFormFactor;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
SudakovFormFactor & operator=(const SudakovFormFactor &);
private:
/**
* Pointer to the splitting function for this Sudakov form factor
*/
SplittingFnPtr splittingFn_;
/**
* Pointer to the coupling for this Sudakov form factor
*/
ShowerAlphaPtr alpha_;
/**
* Maximum value of the PDF weight
*/
double pdfmax_;
/**
* List of the particles this Sudakov is used for to aid in setting up
* interpolation tables if needed
*/
vector<IdList> particles_;
/**
* Option for the inclusion of a factor \f$1/(1-z)\f$ in the PDF estimate
*/
unsigned pdffactor_;
private:
/**
* Option for the type of cut-off to be applied
*/
unsigned int cutOffOption_;
/**
* Parameters for the default Herwig++ cut-off option, i.e. the parameters for
* the \f$Q_g=\max(\frac{\delta-am_q}{b},c)\f$ kinematic cut-off
*/
//@{
/**
* The \f$a\f$ parameter
*/
double a_;
/**
* The \f$b\f$ parameter
*/
double b_;
/**
* The \f$c\f$ parameter
*/
Energy c_;
/**
* Kinematic cutoff used in the parton shower phase space.
*/
Energy kinCutoffScale_;
//@}
/**
* Parameters for the FORTRAN-like cut-off
*/
//@{
/**
* The virtualilty cut-off for gluons
*/
Energy vgcut_;
/**
* The virtuality cut-off for everything else
*/
Energy vqcut_;
//@}
/**
* Parameters for the \f$p_T\f$ cut-off
*/
//@{
/**
* The minimum \f$p_T\f$ for the branching
*/
Energy pTmin_;
/**
* The square of the minimum \f$p_T\f$
*/
Energy2 pT2min_;
//@}
private:
/**
* Member variables to keep the shower kinematics information
* generated by a call to generateNextTimeBranching or generateNextSpaceBranching
*/
//@{
/**
* The energy fraction
*/
double z_;
/**
* The azimuthal angle
*/
double phi_;
/**
* The transverse momentum
*/
Energy pT_;
//@}
/**
* The limits of \f$z\f$ in the splitting
*/
pair<double,double> zlimits_;
/**
* Stuff for the PDFs
*/
//@{
/**
* PDf
*/
tcPDFPtr pdf_;
/**
* Freezing scale
*/
Energy freeze_;
//@}
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of SudakovFormFactor. */
template <>
struct BaseClassTrait<Herwig::SudakovFormFactor,1> {
/** Typedef of the first base class of SudakovFormFactor. */
typedef Interfaced NthBase;
};
/** This template specialization informs ThePEG about the name of
* the SudakovFormFactor class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::SudakovFormFactor>
: public ClassTraitsBase<Herwig::SudakovFormFactor> {
/** Return a platform-independent class name */
static string className() { return "Herwig::SudakovFormFactor"; }
};
/** @endcond */
}
#endif /* HERWIG_SudakovFormFactor_H */
diff --git a/Shower/Default/FS_QTildeShowerKinematics1to2.cc b/Shower/Default/FS_QTildeShowerKinematics1to2.cc
--- a/Shower/Default/FS_QTildeShowerKinematics1to2.cc
+++ b/Shower/Default/FS_QTildeShowerKinematics1to2.cc
@@ -1,154 +1,219 @@
// -*- C++ -*-
//
// FS_QTildeShowerKinematics1to2.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the FS_QTildeShowerKinematics1to2 class.
//
#include "FS_QTildeShowerKinematics1to2.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "Herwig++/Shower/SplittingFunctions/SplittingFunction.h"
#include "Herwig++/Shower/Base/ShowerParticle.h"
using namespace Herwig;
void FS_QTildeShowerKinematics1to2::
updateChildren(const tShowerParticlePtr theParent,
const ShowerParticleVector & theChildren,
bool angularOrder) const {
if(theChildren.size() != 2)
throw Exception() << "FS_QTildeShowerKinematics1to2::updateChildren() "
<< "Warning! too many children!" << Exception::eventerror;
// copy scales etc
Energy dqtilde = scale();
- double dz = z();
- double dphi = phi();
// resize the parameter vectors
if(theParent->showerVariables().empty()) {
theParent->showerVariables().resize(3);
theParent->showerParameters().resize(2);
theParent->showerParameters()[0]=1.;
}
theChildren[0]->showerVariables() .resize(3);
theChildren[0]->showerParameters().resize(2);
theChildren[1]->showerVariables() .resize(3);
theChildren[1]->showerParameters().resize(2);
// note that 1st child gets z, 2nd gets (1-z) by our convention.
if(angularOrder) {
- theChildren[0]->setEvolutionScale(dz*dqtilde);
- theChildren[1]->setEvolutionScale((1.-dz)*dqtilde);
+ theChildren[0]->setEvolutionScale(z()*dqtilde);
+ theChildren[1]->setEvolutionScale((1.-z())*dqtilde);
}
else {
theChildren[0]->setEvolutionScale(dqtilde);
theChildren[1]->setEvolutionScale(dqtilde);
}
// determine alphas of children according to interpretation of z
- theChildren[0]->showerParameters()[0]= dz *theParent->showerParameters()[0];
- theChildren[1]->showerParameters()[0]= (1.-dz)*theParent->showerParameters()[0];
+ theChildren[0]->showerParameters()[0]= z() *theParent->showerParameters()[0];
+ theChildren[1]->showerParameters()[0]= (1.-z())*theParent->showerParameters()[0];
// set the values
theChildren[0]->showerVariables()[0]=
- pT()*cos(dphi) + dz *theParent->showerVariables()[0];
+ pT()*cos(phi()) + z() *theParent->showerVariables()[0];
theChildren[0]->showerVariables()[1]=
- pT()*sin(dphi) + dz *theParent->showerVariables()[1];
+ pT()*sin(phi()) + z() *theParent->showerVariables()[1];
theChildren[1]->showerVariables()[0]=
- - pT()*cos(dphi) + (1.-dz)*theParent->showerVariables()[0];
+ - pT()*cos(phi()) + (1.-z())*theParent->showerVariables()[0];
theChildren[1]->showerVariables()[1]=
- - pT()*sin(dphi) + (1.-dz)*theParent->showerVariables()[1];
+ - pT()*sin(phi()) + (1.-z())*theParent->showerVariables()[1];
for(unsigned int ix=0;ix<2;++ix)
theChildren[ix]->showerVariables()[2]=
sqrt(sqr(theChildren[ix]->showerVariables()[0])+
sqr(theChildren[ix]->showerVariables()[1]));
// set up the colour connections
splittingFn()->colourConnection(theParent,theChildren[0],theChildren[1],false);
// make the products children of the parent
theParent->addChild(theChildren[0]);
theParent->addChild(theChildren[1]);
}
void FS_QTildeShowerKinematics1to2::
reconstructParent(const tShowerParticlePtr theParent,
const ParticleVector & theChildren ) const {
if(theChildren.size() != 2)
throw Exception() << "FS_QTildeShowerKinematics1to2::updateParent() "
<< "Warning! too many children!"
<< Exception::eventerror;
ShowerParticlePtr c1 = dynamic_ptr_cast<ShowerParticlePtr>(theChildren[0]);
ShowerParticlePtr c2 = dynamic_ptr_cast<ShowerParticlePtr>(theChildren[1]);
theParent->showerParameters()[1]=
c1->showerParameters()[1] + c2->showerParameters()[1];
theParent->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,
+ 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]->setVirtualMass(virtualMasses[1]);
+ if(children[1]->children().empty()) children[1]->setVirtualMass(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]);
+
+// cerr << "Result = " << pt2/GeV2 << "\n";
+ Energy2 q2 =
+ sqr(children[0]->virtualMass())/z() +
+ sqr(children[1]->virtualMass())/(1.-z()) +
+ pt2/z()/(1.-z());
+ if(pt2<ZERO) {
+ parent->setVirtualMass(ZERO);
+ }
+ else {
+ parent->setVirtualMass(sqrt(q2));
+ pT(sqrt(pt2));
+ }
+// cerr << "NEGATIVE!!!\n";
+
+// cerr << "testing branching " << pt2/GeV2 << " " << sqr(pT()/GeV) << " " << q2/GeV2
+// << "\n";
+
+// cerr << "testing got here ??? PARENT: "
+// << *parent << " " << virtualMasses[0]/GeV << "\n"
+// << "CHILD : " << *children[0] << " " << children[0]->virtualMass()/GeV << "\n"
+// << "CHILD : " << *children[1] << " " << children[1]->virtualMass()/GeV << "\n";
+
+
+}
+
+void FS_QTildeShowerKinematics1to2::
+resetChildren(const tShowerParticlePtr theParent,
+ const ShowerParticleVector & theChildren) const {
+ // set the values
+ theChildren[0]->showerVariables()[0]=
+ pT()*cos(phi()) + z() *theParent->showerVariables()[0];
+ theChildren[0]->showerVariables()[1]=
+ pT()*sin(phi()) + z() *theParent->showerVariables()[1];
+ theChildren[1]->showerVariables()[0]=
+ - pT()*cos(phi()) + (1.-z())*theParent->showerVariables()[0];
+ theChildren[1]->showerVariables()[1]=
+ - pT()*sin(phi()) + (1.-z())*theParent->showerVariables()[1];
+ for(unsigned int ix=0;ix<2;++ix)
+ theChildren[ix]->showerVariables()[2]=
+ sqrt(sqr(theChildren[ix]->showerVariables()[0])+
+ sqr(theChildren[ix]->showerVariables()[1]));
+ for(unsigned int ix=0;ix<theChildren.size();++ix) {
+ if(theChildren[ix]->children().empty()) continue;
+ ShowerParticleVector children;
+ for(unsigned int iy=0;iy<theChildren[ix]->children().size();++iy)
+ children.push_back(dynamic_ptr_cast<ShowerParticlePtr>
+ (theChildren[ix]->children()[iy]));
+ theChildren[ix]->showerKinematics()->resetChildren(theChildren[ix],children);
+ }
+}
diff --git a/Shower/Default/FS_QTildeShowerKinematics1to2.h b/Shower/Default/FS_QTildeShowerKinematics1to2.h
--- a/Shower/Default/FS_QTildeShowerKinematics1to2.h
+++ b/Shower/Default/FS_QTildeShowerKinematics1to2.h
@@ -1,105 +1,120 @@
// -*- C++ -*-
//
// FS_QTildeShowerKinematics1to2.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_FS_QTildeShowerKinematics1to2_H
#define HERWIG_FS_QTildeShowerKinematics1to2_H
//
// This is the declaration of the FS_QTildeShowerKinematics1to2 class.
//
#include "QTildeShowerKinematics1to2.h"
namespace Herwig {
using namespace ThePEG;
/** \ingroup Shower
*
* This (concrete) class provides the specific Final State shower
* kinematics information.
*
* @see QTildeShowerKinematics1to2
* @see IS_QTildeShowerKinematics1to2
* @see Decay_QTildeShowerKinematics1to2
* @see KinematicsReconstructor
*/
class FS_QTildeShowerKinematics1to2: public QTildeShowerKinematics1to2 {
public:
/**
* Default constructor
*/
inline FS_QTildeShowerKinematics1to2() {}
/**
* The updateChildren, updateParent and updateLast
* members to update the values of the \f$\alpha\f$ and
* \f$p_\perp\f$ variables during the shower evolution.
*/
//@{
/**
* Along with the showering evolution --- going forward for
* time-like (forward) evolution, and going backward for space-like
* (backward) evolution --- the kinematical variables of the
* branching products are calculated and updated from the knowledge
* of the parent kinematics. This method is used by the
* ForwardShowerEvolver.
* ***ACHTUNG*** Might be extended to update colour connections as well.
* @param theParent The branching particle
* @param theChildren The particles produced in the branching
* @param angularOrder Whether or not to apply angular ordering
*/
virtual void updateChildren( const tShowerParticlePtr theParent,
const ShowerParticleVector & theChildren,
bool angularOrder) const;
+ virtual void resetChildren( const tShowerParticlePtr theParent,
+ const ShowerParticleVector & theChildren) const;
+
+
+ /**
+ * Update the parent Kinematics from the knowledge of the kinematics
+ * of the children. This method will be used by the KinematicsReconstructor.
+ * @param theParent The parent
+ * @param theChildren The children
+ * @param angularOrder Whether or not to apply angular ordering
+ */
+ virtual void updateParent(const tShowerParticlePtr theParent,
+ const ShowerParticleVector & theChildren,
+ bool angularOrder) const;
+
/**
* Update the parent Kinematics from the knowledge of the kinematics
* of the children. This method will be used by the
* KinematicsReconstructor.
*/
virtual void reconstructParent( const tShowerParticlePtr theParent,
const ParticleVector & theChildren ) const;
/**
* Update the kinematical data of a particle when a reconstruction
* fixpoint was found. This will highly depend on the kind of
* kinematics chosen and will be defined in the inherited concrete
* classes. This method will be used by the KinematicsReconstructor.
* @param theLast The particle to update
* @param iopt The option for the momentum reconstruction
* - 0 is in the rest frame of the pair of reference vectors
* - 1 is in the rest frame of the p vector
* @param mass The mass to be used, if less than zero on-shell
*/
virtual void reconstructLast(const tShowerParticlePtr theLast,
unsigned int iopt, Energy mass=-1.*GeV) const;
/**
* Perform any initial calculations needed after the branching has been selected
* @param particle The branching particle
* @param parent The bema particle for the jet if needed
*/
virtual void initialize(ShowerParticle & particle,PPtr parent);
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
FS_QTildeShowerKinematics1to2 & operator=(const FS_QTildeShowerKinematics1to2 &);
};
}
#endif /* HERWIG_FS_QTildeShowerKinematics1to2_H */
diff --git a/Shower/Default/QTildeSudakov.cc b/Shower/Default/QTildeSudakov.cc
--- a/Shower/Default/QTildeSudakov.cc
+++ b/Shower/Default/QTildeSudakov.cc
@@ -1,383 +1,359 @@
// -*- C++ -*-
//
// QTildeSudakov.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the QTildeSudakov class.
//
#include "QTildeSudakov.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/PDT/ParticleData.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "Herwig++/Shower/Default/FS_QTildeShowerKinematics1to2.h"
#include "Herwig++/Shower/Default/IS_QTildeShowerKinematics1to2.h"
#include "Herwig++/Shower/Default/Decay_QTildeShowerKinematics1to2.h"
using namespace Herwig;
NoPIOClassDescription<QTildeSudakov> QTildeSudakov::initQTildeSudakov;
// Definition of the static class description member.
void QTildeSudakov::Init() {
static ClassDocumentation<QTildeSudakov> documentation
("The QTildeSudakov class implements the Sudakov form factor for ordering it"
" qtilde");
}
bool QTildeSudakov::guessTimeLike(Energy2 &t,Energy2 tmin,double enhance) {
Energy2 told = t;
// calculate limits on z and if lower>upper return
if(!computeTimeLikeLimits(t)) return false;
// guess values of t and z
t = guesst(told,0,ids_,enhance,ids_[1]==ids_[2]);
z(guessz(0,ids_));
// actual values for z-limits
if(!computeTimeLikeLimits(t)) return false;
if(t<tmin) {
t=-1.0*GeV2;
return false;
}
else
return true;
}
bool QTildeSudakov::guessSpaceLike(Energy2 &t, Energy2 tmin, const double x,
double enhance) {
Energy2 told = t;
// calculate limits on z if lower>upper return
if(!computeSpaceLikeLimits(t,x)) return false;
// guess values of t and z
t = guesst(told,1,ids_,enhance,ids_[1]==ids_[2]);
z(guessz(1,ids_));
// actual values for z-limits
if(!computeSpaceLikeLimits(t,x)) return false;
if(t<tmin) {
t=-1.0*GeV2;
return false;
}
else
return true;
}
bool QTildeSudakov::PSVeto(const Energy2 t) {
// still inside PS, return true if outside
// check vs overestimated limits
if(z() < zLimits().first || z() > zLimits().second) return true;
// compute the pts
Energy2 pt2=sqr(z()*(1.-z()))*t-masssquared_[1]*(1.-z())-masssquared_[2]*z();
if(ids_[0]!=ParticleID::g) pt2+=z()*(1.-z())*masssquared_[0];
// if pt2<0 veto
if(pt2<pT2min()) return true;
// otherwise calculate pt and return
pT(sqrt(pt2));
return false;
}
ShoKinPtr QTildeSudakov::generateNextTimeBranching(const Energy startingScale,
const IdList &ids,const bool cc,
double enhance) {
// First reset the internal kinematics variables that can
// have been eventually set in the previous call to the method.
q_ = ZERO;
z(0.);
phi(0.);
// perform initialization
Energy2 tmax(sqr(startingScale)),tmin;
initialize(ids,tmin,cc);
// check max > min
if(tmax<=tmin) return ShoKinPtr();
// calculate next value of t using veto algorithm
Energy2 t(tmax);
do {
if(!guessTimeLike(t,tmin,enhance)) break;
}
while(PSVeto(t) || SplittingFnVeto(z()*(1.-z())*t,ids,true) ||
alphaSVeto(sqr(z()*(1.-z()))*t));
if(t > ZERO) q_ = sqrt(t);
else q_ = -1.*MeV;
phi(Constants::twopi*UseRandom::rnd());
if(q_ < ZERO) return ShoKinPtr();
// return the ShowerKinematics object
return createFinalStateBranching(q_,z(),phi(),pT());
}
ShoKinPtr QTildeSudakov::
generateNextSpaceBranching(const Energy startingQ,
const IdList &ids,
double x,bool cc,
double enhance,
Ptr<BeamParticleData>::transient_const_pointer beam) {
// First reset the internal kinematics variables that can
// have been eventually set in the previous call to the method.
q_ = ZERO;
z(0.);
phi(0.);
// perform the initialization
Energy2 tmax(sqr(startingQ)),tmin;
initialize(ids,tmin,cc);
// check max > min
if(tmax<=tmin) return ShoKinPtr();
// extract the partons which are needed for the PDF veto
// Different order, incoming parton is id = 1, outgoing are id=0,2
tcPDPtr parton0 = getParticleData(ids[0]);
tcPDPtr parton1 = getParticleData(ids[1]);
if(cc) {
if(parton0->CC()) parton0 = parton0->CC();
if(parton1->CC()) parton1 = parton1->CC();
}
// calculate next value of t using veto algorithm
Energy2 t(tmax),pt2(ZERO);
do {
if(!guessSpaceLike(t,tmin,x,enhance)) break;
pt2=sqr(1.-z())*t-z()*masssquared_[2];
}
while(z() > zLimits().second ||
SplittingFnVeto((1.-z())*t/z(),ids,true) ||
alphaSVeto(sqr(1.-z())*t) ||
PDFVeto(t,x,parton0,parton1,beam) || pt2 < pT2min() );
if(t > ZERO && zLimits().first < zLimits().second) q_ = sqrt(t);
else return ShoKinPtr();
phi(Constants::twopi*UseRandom::rnd());
pT(sqrt(pt2));
// create the ShowerKinematics and return it
return createInitialStateBranching(q_,z(),phi(),pT());
}
void QTildeSudakov::initialize(const IdList & ids, Energy2 & tmin,const bool cc) {
ids_=ids;
if(cc) {
for(unsigned int ix=0;ix<ids.size();++ix) {
if(getParticleData(ids[ix])->CC()) ids_[ix]*=-1;
}
}
- masses_.clear();
+ tmin = cutOffOption() != 2 ? ZERO : 4.*pT2min();
+ masses_ = virtualMasses(ids);
masssquared_.clear();
- tmin=ZERO;
- if(cutOffOption() == 0) {
- for(unsigned int ix=0;ix<ids_.size();++ix)
- masses_.push_back(getParticleData(ids_[ix])->mass());
- Energy kinCutoff=
- kinematicCutOff(kinScale(),*std::max_element(masses_.begin(),masses_.end()));
- for(unsigned int ix=0;ix<masses_.size();++ix)
- masses_[ix]=max(kinCutoff,masses_[ix]);
- }
- else if(cutOffOption() == 1) {
- for(unsigned int ix=0;ix<ids_.size();++ix) {
- masses_.push_back(getParticleData(ids_[ix])->mass());
- masses_.back() += ids_[ix]==ParticleID::g ? vgCut() : vqCut();
- }
- }
- else if(cutOffOption() == 2) {
- for(unsigned int ix=0;ix<ids_.size();++ix)
- masses_.push_back(getParticleData(ids_[ix])->mass());
- tmin = 4.*pT2min();
- }
- else {
- throw Exception() << "Unknown option for the cut-off"
- << " in QTildeSudakov::initialize()"
- << Exception::runerror;
- }
for(unsigned int ix=0;ix<masses_.size();++ix) {
masssquared_.push_back(sqr(masses_[ix]));
if(ix>0) tmin=max(masssquared_[ix],tmin);
}
}
ShoKinPtr QTildeSudakov::generateNextDecayBranching(const Energy startingScale,
const Energy stoppingScale,
const Energy minmass,
const IdList &ids,
const bool cc,
double enhance) {
// First reset the internal kinematics variables that can
// have been eventually set in the previous call to this method.
q_ = Constants::MaxEnergy;
z(0.);
phi(0.);
// perform initialisation
Energy2 tmax(sqr(stoppingScale)),tmin;
initialize(ids,tmin,cc);
tmin=sqr(startingScale);
// check some branching possible
if(tmax<=tmin) return ShoKinPtr();
// perform the evolution
Energy2 t(tmin),pt2(-MeV2);
do {
if(!guessDecay(t,tmax,minmass,enhance)) break;
pt2 = sqr(1.-z())*(t-masssquared_[0])-z()*masssquared_[2];
}
while(SplittingFnVeto((1.-z())*t/z(),ids,true)||
alphaSVeto(sqr(1.-z())*t) ||
pt2<pT2min() ||
t*(1.-z())>masssquared_[0]-sqr(minmass));
if(t > ZERO) {
q_ = sqrt(t);
pT(sqrt(pt2));
}
else return ShoKinPtr();
phi(Constants::twopi*UseRandom::rnd());
// create the ShowerKinematics object
return createDecayBranching(q_,z(),phi(),pT());
}
bool QTildeSudakov::guessDecay(Energy2 &t,Energy2 tmax, Energy minmass,
double enhance) {
// previous scale
Energy2 told = t;
// overestimated limits on z
if(tmax<masssquared_[0]) {
t=-1.0*GeV2;
return false;
}
pair<double,double> limits=make_pair(sqr(minmass/masses_[0]),
1.-masses_[2]/sqrt(tmax-masssquared_[0])
+0.5*masssquared_[2]/(tmax-masssquared_[0]));
zLimits(limits);
if(zLimits().second<zLimits().first) {
t=-1.0*GeV2;
return false;
}
// guess values of t and z
t = guesst(told,2,ids_,enhance,ids_[1]==ids_[2]);
z(guessz(2,ids_));
// actual values for z-limits
if(t<masssquared_[0]) {
t=-1.0*GeV2;
return false;
}
limits=make_pair(sqr(minmass/masses_[0]),
1.-masses_[2]/sqrt(t-masssquared_[0])
+0.5*masssquared_[2]/(t-masssquared_[0]));
zLimits(limits);
if(t>tmax||zLimits().second<zLimits().first) {
t=-1.0*GeV2;
return false;
}
else
return true;
}
bool QTildeSudakov::computeTimeLikeLimits(Energy2 & t) {
if (t == ZERO) {
t=-1.*GeV2;
return false;
}
// special case for gluon radiating
pair<double,double> limits;
if(ids_[0]==ParticleID::g) {
// no emission possible
if(t<16.*masssquared_[1]) {
t=-1.*GeV2;
return false;
}
// overestimate of the limits
limits.first = 0.5*(1.-sqrt(1.-4.*sqrt((masssquared_[1]+pT2min())/t)));
limits.second = 1.-limits.first;
}
// special case for radiated particle is gluon
else if(ids_[2]==ParticleID::g) {
limits.first = sqrt((masssquared_[1]+pT2min())/t);
limits.second = 1.-sqrt((masssquared_[2]+pT2min())/t);
}
else if(ids_[1]==ParticleID::g) {
limits.second = sqrt((masssquared_[2]+pT2min())/t);
limits.first = 1.-sqrt((masssquared_[1]+pT2min())/t);
}
else {
limits.first = (masssquared_[1]+pT2min())/t;
limits.second = 1.-(masssquared_[2]+pT2min())/t;
}
if(limits.first>=limits.second) {
t=-1.*GeV2;
return false;
}
zLimits(limits);
return true;
}
bool QTildeSudakov::computeSpaceLikeLimits(Energy2 & t, double x) {
if (t == ZERO) {
t=-1.*GeV2;
return false;
}
pair<double,double> limits;
// compute the limits
limits.first = x;
double yy = 1.+0.5*masssquared_[2]/t;
limits.second = yy - sqrt(sqr(yy)-1.+pT2min()/t);
// return false if lower>upper
zLimits(limits);
if(limits.second<limits.first) {
t=-1.*GeV2;
return false;
}
else
return true;
}
Energy QTildeSudakov::calculateScale(double zin, Energy pt, IdList ids,
unsigned int iopt) {
Energy2 tmin;
initialize(ids,tmin,false);
// final-state branching
if(iopt==0) {
Energy2 scale=(sqr(pt)+masssquared_[1]*(1.-zin)+masssquared_[2]*zin);
if(ids[0]!=ParticleID::g) scale -= zin*(1.-zin)*masssquared_[0];
scale /= sqr(zin*(1-zin));
return scale<=ZERO ? sqrt(tmin) : sqrt(scale);
}
else if(iopt==1) {
Energy2 scale=(sqr(pt)+zin*masssquared_[2])/sqr(1.-zin);
return scale<=ZERO ? sqrt(tmin) : sqrt(scale);
}
else if(iopt==2) {
Energy2 scale = (sqr(pt)+zin*masssquared_[2])/sqr(1.-zin)+masssquared_[0];
return scale<=ZERO ? sqrt(tmin) : sqrt(scale);
}
else {
throw Exception() << "Unknown option in QTildeSudakov::calculateScale() "
<< "iopt = " << iopt << Exception::runerror;
}
}
ShoKinPtr QTildeSudakov::createFinalStateBranching(Energy scale,double z,
double phi, Energy pt) {
ShoKinPtr showerKin = new_ptr(FS_QTildeShowerKinematics1to2());
showerKin->scale(scale);
showerKin->z(z);
showerKin->phi(phi);
showerKin->pT(pt);
showerKin->SudakovFormFactor(this);
return showerKin;
}
ShoKinPtr QTildeSudakov::createInitialStateBranching(Energy scale,double z,
double phi, Energy pt) {
ShoKinPtr showerKin = new_ptr(IS_QTildeShowerKinematics1to2());
showerKin->scale(scale);
showerKin->z(z);
showerKin->phi(phi);
showerKin->pT(pt);
showerKin->SudakovFormFactor(this);
return showerKin;
}
ShoKinPtr QTildeSudakov::createDecayBranching(Energy scale,double z,
double phi, Energy pt) {
ShoKinPtr showerKin = new_ptr(Decay_QTildeShowerKinematics1to2());
showerKin->scale(scale);
showerKin->z(z);
showerKin->phi(phi);
showerKin->pT(pt);
showerKin->SudakovFormFactor(this);
return showerKin;
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 1:01 PM (1 d, 7 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4022874
Default Alt Text
(190 KB)

Event Timeline