Page MenuHomeHEPForge
Files F7878949

No OneTemporary
Actions

View Options

diff --git a/Shower/QTilde/QTildeShowerHandler.cc b/Shower/QTilde/QTildeShowerHandler.cc
--- a/Shower/QTilde/QTildeShowerHandler.cc
+++ b/Shower/QTilde/QTildeShowerHandler.cc
@@ -1,3216 +1,3216 @@
// -*- C++ -*-
//
// QTildeShowerHandler.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 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 QTildeShowerHandler class.
//
#include "QTildeShowerHandler.h"
#include "ThePEG/Interface/Deleted.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Utilities/EnumIO.h"
#include "Herwig/Shower/QTilde/Base/ShowerParticle.h"
#include "Herwig/PDF/MPIPDF.h"
#include "Herwig/PDF/MinBiasPDF.h"
#include "Herwig/Shower/QTilde/Base/ShowerTree.h"
#include "Herwig/Shower/QTilde/Base/HardTree.h"
#include "Herwig/Shower/QTilde/Kinematics/KinematicsReconstructor.h"
#include "Herwig/Shower/QTilde/Base/PartnerFinder.h"
#include "Herwig/PDF/HwRemDecayer.h"
#include "Herwig/Shower/QTilde/Base/ShowerVertex.h"
#include "ThePEG/Repository/CurrentGenerator.h"
#include "Herwig/MatrixElement/Matchbox/Base/SubtractedME.h"
#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"
#include "ThePEG/PDF/PartonExtractor.h"
#include "Herwig/Shower/RealEmissionProcess.h"
#include "Herwig/Shower/QTilde/Kinematics/FS_QTildeShowerKinematics1to2.h"
#include "Herwig/Shower/QTilde/Kinematics/IS_QTildeShowerKinematics1to2.h"
#include "Herwig/Shower/QTilde/Kinematics/Decay_QTildeShowerKinematics1to2.h"
using namespace Herwig;
bool QTildeShowerHandler::_hardEmissionWarn = true;
bool QTildeShowerHandler::_missingTruncWarn = true;
QTildeShowerHandler::QTildeShowerHandler() :
- _maxtry(100), _meCorrMode(1), _evolutionScheme(2),
+ _maxtry(100), _meCorrMode(1), _evolutionScheme(1),
_hardVetoReadOption(false),
_iptrms(ZERO), _beta(0.), _gamma(ZERO), _iptmax(),
_limitEmissions(0), _initialenhance(1.), _finalenhance(1.),
_nReWeight(100), _reWeight(false),
interaction_(ShowerInteraction::Both),
_trunc_Mode(true), _hardEmission(1),
_softOpt(2), _hardPOWHEG(false), muPt(ZERO)
{}
QTildeShowerHandler::~QTildeShowerHandler() {}
IBPtr QTildeShowerHandler::clone() const {
return new_ptr(*this);
}
IBPtr QTildeShowerHandler::fullclone() const {
return new_ptr(*this);
}
void QTildeShowerHandler::persistentOutput(PersistentOStream & os) const {
os << _splittingGenerator << _maxtry
<< _meCorrMode << _hardVetoReadOption
<< _limitEmissions << _softOpt << _hardPOWHEG
<< ounit(_iptrms,GeV) << _beta << ounit(_gamma,GeV) << ounit(_iptmax,GeV)
<< _vetoes << _fullShowerVetoes << _nReWeight << _reWeight
<< _trunc_Mode << _hardEmission << _evolutionScheme
<< ounit(muPt,GeV) << oenum(interaction_)
<< _reconstructor << _partnerfinder;
}
void QTildeShowerHandler::persistentInput(PersistentIStream & is, int) {
is >> _splittingGenerator >> _maxtry
>> _meCorrMode >> _hardVetoReadOption
>> _limitEmissions >> _softOpt >> _hardPOWHEG
>> iunit(_iptrms,GeV) >> _beta >> iunit(_gamma,GeV) >> iunit(_iptmax,GeV)
>> _vetoes >> _fullShowerVetoes >> _nReWeight >> _reWeight
>> _trunc_Mode >> _hardEmission >> _evolutionScheme
>> iunit(muPt,GeV) >> ienum(interaction_)
>> _reconstructor >> _partnerfinder;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<QTildeShowerHandler,ShowerHandler>
describeHerwigQTildeShowerHandler("Herwig::QTildeShowerHandler", "HwShower.so");
void QTildeShowerHandler::Init() {
static ClassDocumentation<QTildeShowerHandler> documentation
("TheQTildeShowerHandler class is the main class"
" for the angular-ordered parton shower",
"The Shower evolution was performed using an algorithm described in "
"\\cite{Marchesini:1983bm,Marchesini:1987cf,Gieseke:2003rz,Bahr:2008pv}.",
"%\\cite{Marchesini:1983bm}\n"
"\\bibitem{Marchesini:1983bm}\n"
" G.~Marchesini and B.~R.~Webber,\n"
" ``Simulation Of QCD Jets Including Soft Gluon Interference,''\n"
" Nucl.\\ Phys.\\ B {\\bf 238}, 1 (1984).\n"
" %%CITATION = NUPHA,B238,1;%%\n"
"%\\cite{Marchesini:1987cf}\n"
"\\bibitem{Marchesini:1987cf}\n"
" G.~Marchesini and B.~R.~Webber,\n"
" ``Monte Carlo Simulation of General Hard Processes with Coherent QCD\n"
" Radiation,''\n"
" Nucl.\\ Phys.\\ B {\\bf 310}, 461 (1988).\n"
" %%CITATION = NUPHA,B310,461;%%\n"
"%\\cite{Gieseke:2003rz}\n"
"\\bibitem{Gieseke:2003rz}\n"
" S.~Gieseke, P.~Stephens and B.~Webber,\n"
" ``New formalism for QCD parton showers,''\n"
" JHEP {\\bf 0312}, 045 (2003)\n"
" [arXiv:hep-ph/0310083].\n"
" %%CITATION = JHEPA,0312,045;%%\n"
);
static Reference<QTildeShowerHandler,SplittingGenerator>
interfaceSplitGen("SplittingGenerator",
"A reference to the SplittingGenerator object",
&Herwig::QTildeShowerHandler::_splittingGenerator,
false, false, true, false);
static Parameter<QTildeShowerHandler,unsigned int> interfaceMaxTry
("MaxTry",
"The maximum number of attempts to generate the shower from a"
" particular ShowerTree",
&QTildeShowerHandler::_maxtry, 100, 1, 100000,
false, false, Interface::limited);
static Parameter<QTildeShowerHandler,unsigned int> interfaceNReWeight
("NReWeight",
"The number of attempts for the shower when reweighting",
&QTildeShowerHandler::_nReWeight, 100, 10, 10000,
false, false, Interface::limited);
static Switch<QTildeShowerHandler, unsigned int> ifaceMECorrMode
("MECorrMode",
"Choice of the ME Correction Mode",
&QTildeShowerHandler::_meCorrMode, 1, false, false);
static SwitchOption on
(ifaceMECorrMode,"HardPlusSoft","hard+soft on", 1);
static SwitchOption hard
(ifaceMECorrMode,"Hard","only hard on", 2);
static SwitchOption soft
(ifaceMECorrMode,"Soft","only soft on", 3);
static Switch<QTildeShowerHandler, bool> ifaceHardVetoReadOption
("HardVetoReadOption",
"Apply read-in scale veto to all collisions or just the primary one?",
&QTildeShowerHandler::_hardVetoReadOption, false, false, false);
static SwitchOption AllCollisions
(ifaceHardVetoReadOption,
"AllCollisions",
"Read-in pT veto applied to primary and secondary collisions.",
false);
static SwitchOption PrimaryCollision
(ifaceHardVetoReadOption,
"PrimaryCollision",
"Read-in pT veto applied to primary but not secondary collisions.",
true);
static Parameter<QTildeShowerHandler, Energy> ifaceiptrms
("IntrinsicPtGaussian",
"RMS of intrinsic pT of Gaussian distribution:\n"
"2*(1-Beta)*exp(-sqr(intrinsicpT/RMS))/sqr(RMS)",
&QTildeShowerHandler::_iptrms, GeV, ZERO, ZERO, 1000000.0*GeV,
false, false, Interface::limited);
static Parameter<QTildeShowerHandler, double> ifacebeta
("IntrinsicPtBeta",
"Proportion of inverse quadratic distribution in generating intrinsic pT.\n"
"(1-Beta) is the proportion of Gaussian distribution",
&QTildeShowerHandler::_beta, 0, 0, 1,
false, false, Interface::limited);
static Parameter<QTildeShowerHandler, Energy> ifacegamma
("IntrinsicPtGamma",
"Parameter for inverse quadratic:\n"
"2*Beta*Gamma/(sqr(Gamma)+sqr(intrinsicpT))",
&QTildeShowerHandler::_gamma,GeV, ZERO, ZERO, 100000.0*GeV,
false, false, Interface::limited);
static Parameter<QTildeShowerHandler, Energy> ifaceiptmax
("IntrinsicPtIptmax",
"Upper bound on intrinsic pT for inverse quadratic",
&QTildeShowerHandler::_iptmax,GeV, ZERO, ZERO, 100000.0*GeV,
false, false, Interface::limited);
static RefVector<QTildeShowerHandler,ShowerVeto> ifaceVetoes
("Vetoes",
"The vetoes to be checked during showering",
&QTildeShowerHandler::_vetoes, -1,
false,false,true,true,false);
static RefVector<QTildeShowerHandler,FullShowerVeto> interfaceFullShowerVetoes
("FullShowerVetoes",
"The vetos to be appliede on the full final state of the shower",
&QTildeShowerHandler::_fullShowerVetoes, -1, false, false, true, false, false);
static Switch<QTildeShowerHandler,unsigned int> interfaceLimitEmissions
("LimitEmissions",
"Limit the number and type of emissions for testing",
&QTildeShowerHandler::_limitEmissions, 0, false, false);
static SwitchOption interfaceLimitEmissionsNoLimit
(interfaceLimitEmissions,
"NoLimit",
"Allow an arbitrary number of emissions",
0);
static SwitchOption interfaceLimitEmissionsOneInitialStateEmission
(interfaceLimitEmissions,
"OneInitialStateEmission",
"Allow one emission in the initial state and none in the final state",
1);
static SwitchOption interfaceLimitEmissionsOneFinalStateEmission
(interfaceLimitEmissions,
"OneFinalStateEmission",
"Allow one emission in the final state and none in the initial state",
2);
static SwitchOption interfaceLimitEmissionsHardOnly
(interfaceLimitEmissions,
"HardOnly",
"Only allow radiation from the hard ME correction",
3);
static SwitchOption interfaceLimitEmissionsOneEmission
(interfaceLimitEmissions,
"OneEmission",
"Allow one emission in either the final state or initial state, but not both",
4);
static Switch<QTildeShowerHandler,bool> interfaceTruncMode
("TruncatedShower", "Include the truncated shower?",
&QTildeShowerHandler::_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<QTildeShowerHandler,int> interfaceHardEmission
("HardEmission",
"Whether to use ME corrections or POWHEG for the hardest emission",
&QTildeShowerHandler::_hardEmission, 0, false, false);
static SwitchOption interfaceHardEmissionNone
(interfaceHardEmission,
"None",
"No Corrections",
0);
static SwitchOption interfaceHardEmissionMECorrection
(interfaceHardEmission,
"MECorrection",
"Old fashioned ME correction",
1);
static SwitchOption interfaceHardEmissionPOWHEG
(interfaceHardEmission,
"POWHEG",
"Powheg style hard emission",
2);
static Switch<QTildeShowerHandler,ShowerInteraction> interfaceInteractions
("Interactions",
"The interactions to be used in the shower",
&QTildeShowerHandler::interaction_, ShowerInteraction::Both, false, false);
static SwitchOption interfaceInteractionsQCD
(interfaceInteractions,
"QCD",
"Only QCD radiation",
ShowerInteraction::QCD);
static SwitchOption interfaceInteractionsQED
(interfaceInteractions,
"QED",
"Only QEd radiation",
ShowerInteraction::QED);
static SwitchOption interfaceInteractionsQCDandQED
(interfaceInteractions,
"QCDandQED",
"Both QED and QCD radiation",
ShowerInteraction::Both);
static Deleted<QTildeShowerHandler> delReconstructionOption
("ReconstructionOption", "The old reconstruction option switch has been replaced with"
" the new EvolutionScheme switch, see arXiv:1904.11866 for details");
static Switch<QTildeShowerHandler,unsigned int> interfaceEvolutionScheme
("EvolutionScheme",
"The scheme to interpret the evolution variable in the case of multple emission.",
&QTildeShowerHandler::_evolutionScheme, 1, false, false);
static SwitchOption interfaceEvolutionSchemepT
(interfaceEvolutionScheme,
"pT",
"pT scheme",
0);
static SwitchOption interfaceEvolutionSchemeDotProduct
(interfaceEvolutionScheme,
"DotProduct",
"Dot-product scheme",
1);
static SwitchOption interfaceEvolutionSchemeQ2
(interfaceEvolutionScheme,
"Q2",
"Q2 scheme",
2);
static Switch<QTildeShowerHandler,unsigned int> interfaceSoftCorrelations
("SoftCorrelations",
"Option for the treatment of soft correlations in the parton shower",
&QTildeShowerHandler::_softOpt, 2, false, false);
static SwitchOption interfaceSoftCorrelationsNone
(interfaceSoftCorrelations,
"No",
"No soft correlations",
0);
static SwitchOption interfaceSoftCorrelationsFull
(interfaceSoftCorrelations,
"Full",
"Use the full eikonal",
1);
static SwitchOption interfaceSoftCorrelationsSingular
(interfaceSoftCorrelations,
"Singular",
"Use original Webber-Marchisini form",
2);
static Switch<QTildeShowerHandler,bool> interfaceHardPOWHEG
("HardPOWHEG",
"Treatment of powheg emissions which are too hard to have a shower interpretation",
&QTildeShowerHandler::_hardPOWHEG, false, false, false);
static SwitchOption interfaceHardPOWHEGAsShower
(interfaceHardPOWHEG,
"AsShower",
"Still interpret as shower emissions",
false);
static SwitchOption interfaceHardPOWHEGRealEmission
(interfaceHardPOWHEG,
"RealEmission",
"Generate shower from the real emmission configuration",
true);
static Reference<QTildeShowerHandler,KinematicsReconstructor> interfaceKinematicsReconstructor
("KinematicsReconstructor",
"Reference to the KinematicsReconstructor object",
&QTildeShowerHandler::_reconstructor, false, false, true, false, false);
static Reference<QTildeShowerHandler,PartnerFinder> interfacePartnerFinder
("PartnerFinder",
"Reference to the PartnerFinder object",
&QTildeShowerHandler::_partnerfinder, false, false, true, false, false);
}
tPPair QTildeShowerHandler::cascade(tSubProPtr sub,
XCPtr xcomb) {
// use me for reference in tex file etc
useMe();
prepareCascade(sub);
// set things up in the base class
resetWeights();
hard_=ShowerTreePtr();
decay_.clear();
done_.clear();
// start of the try block for the whole showering process
unsigned int countFailures=0;
while (countFailures<maxtry()) {
try {
decay_.clear();
done_.clear();
PerturbativeProcessPtr hard;
DecayProcessMap decay;
splitHardProcess(firstInteraction() ? tagged() :
tPVector(currentSubProcess()->outgoing().begin(),
currentSubProcess()->outgoing().end()),
hard,decay);
ShowerTree::constructTrees(hard_,decay_,hard,decay);
// if no hard process
if(!hard_) throw Exception() << "Shower starting with a decay"
<< "is not implemented"
<< Exception::runerror;
// perform the shower for the hard process
showerHardProcess(hard_,xcomb);
done_.push_back(hard_);
hard_->updateAfterShower(decay_);
// if no decaying particles to shower break out of the loop
if(decay_.empty()) break;
// shower the decay products
while(!decay_.empty()) {
// find particle whose production process has been showered
ShowerDecayMap::iterator dit = decay_.begin();
while(!dit->second->parent()->hasShowered() && dit!=decay_.end()) ++dit;
assert(dit!=decay_.end());
// get the particle
ShowerTreePtr decayingTree = dit->second;
// remove it from the multimap
decay_.erase(dit);
// make sure the particle has been decayed
QTildeShowerHandler::decay(decayingTree,decay_);
// now shower the decay
showerDecay(decayingTree);
done_.push_back(decayingTree);
decayingTree->updateAfterShower(decay_);
}
// suceeded break out of the loop
break;
}
catch (KinematicsReconstructionVeto) {
resetWeights();
++countFailures;
}
catch ( ... ) {
hard_=ShowerTreePtr();
decay_.clear();
done_.clear();
throw;
}
}
// if loop exited because of too many tries, throw event away
if (countFailures >= maxtry()) {
resetWeights();
hard_=ShowerTreePtr();
decay_.clear();
done_.clear();
throw Exception() << "Too many tries for main while loop "
<< "in QTildeShowerHandler::cascade()."
<< Exception::eventerror;
}
//enter the particles in the event record
fillEventRecord();
// clear storage
hard_=ShowerTreePtr();
decay_.clear();
done_.clear();
// non hadronic case return
if (!isResolvedHadron(incomingBeams().first ) &&
!isResolvedHadron(incomingBeams().second) )
return incomingBeams();
// remake the remnants (needs to be after the colours are sorted
// out in the insertion into the event record)
if ( firstInteraction() ) return remakeRemnant(sub->incoming());
//Return the new pair of incoming partons. remakeRemnant is not
//necessary here, because the secondary interactions are not yet
//connected to the remnants.
return make_pair(findFirstParton(sub->incoming().first ),
findFirstParton(sub->incoming().second));
}
void QTildeShowerHandler::fillEventRecord() {
// create a new step
StepPtr pstep = newStep();
assert(!done_.empty());
assert(done_[0]->isHard());
// insert the steps
for(unsigned int ix=0;ix<done_.size();++ix) {
done_[ix]->fillEventRecord(pstep,doISR(),doFSR());
}
}
HardTreePtr QTildeShowerHandler::generateCKKW(ShowerTreePtr ) const {
return HardTreePtr();
}
void QTildeShowerHandler::doinit() {
ShowerHandler::doinit();
// check on the reweighting
for(unsigned int ix=0;ix<_fullShowerVetoes.size();++ix) {
if(_fullShowerVetoes[ix]->behaviour()==1) {
_reWeight = true;
break;
}
}
if(_reWeight && maximumTries()<_nReWeight) {
throw Exception() << "Reweight being performed in the shower but the number of attempts for the"
<< "shower is less than that for the reweighting.\n"
<< "Maximum number of attempt for the shower "
<< fullName() << ":MaxTry is " << maximumTries() << "\nand for reweighting is "
<< fullName() << ":NReWeight is " << _nReWeight << "\n"
<< "we recommend the number of attempts is 10 times the number for reweighting\n"
<< Exception::runerror;
}
ShowerTree::_vmin2 = vMin();
ShowerTree::_spaceTime = includeSpaceTime();
}
void QTildeShowerHandler::doinitrun() {
ShowerHandler::doinitrun();
ShowerTree::_vmin2 = vMin();
ShowerTree::_spaceTime = includeSpaceTime();
}
void QTildeShowerHandler::generateIntrinsicpT(vector<ShowerProgenitorPtr> particlesToShower) {
_intrinsic.clear();
if ( !ipTon() || !doISR() ) return;
// don't do anything for the moment for secondary scatters
if( !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 QTildeShowerHandler::setupMaximumScales(const vector<ShowerProgenitorPtr> & p,
XCPtr xcomb) {
// let POWHEG events radiate freely
if(_hardEmission==2&&hardTree()) {
vector<ShowerProgenitorPtr>::const_iterator ckt = p.begin();
for (; ckt != p.end(); ckt++) (*ckt)->maxHardPt(Constants::MaxEnergy);
return;
}
// return if no vetos
if (!restrictPhasespace()) return;
// find out if hard partonic subprocess.
bool isPartonic(false);
map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
cit = _currenttree->incomingLines().begin();
Lorentz5Momentum pcm;
for(; cit!=currentTree()->incomingLines().end(); ++cit) {
pcm += cit->first->progenitor()->momentum();
isPartonic |= cit->first->progenitor()->coloured();
}
// find minimum 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 = generator()->maximumCMEnergy();
// general case calculate the scale
if ( !hardScaleIsMuF() || (hardVetoReadOption()&&!firstInteraction()) ) {
// scattering process
if(currentTree()->isHard()) {
assert(xcomb);
// coloured incoming particles
if (isPartonic) {
map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
cjt = currentTree()->outgoingLines().begin();
for(; cjt!=currentTree()->outgoingLines().end(); ++cjt) {
if (cjt->first->progenitor()->coloured())
ptmax = min(ptmax,cjt->first->progenitor()->momentum().mt());
}
}
if (ptmax == generator()->maximumCMEnergy() ) ptmax = pcm.m();
if(hardScaleIsMuF()&&hardVetoReadOption()&&
!firstInteraction()) {
ptmax=min(ptmax,sqrt(xcomb->lastShowerScale()));
}
}
// decay, incoming() is the decaying particle.
else {
ptmax = currentTree()->incomingLines().begin()->first
->progenitor()->momentum().mass();
}
}
// hepeup.SCALUP is written into the lastXComb by the
// LesHouchesReader itself - use this by user's choice.
// Can be more general than this.
else {
if(currentTree()->isHard()) {
assert(xcomb);
ptmax = sqrt( xcomb->lastShowerScale() );
}
else {
ptmax = currentTree()->incomingLines().begin()->first
->progenitor()->momentum().mass();
}
}
ptmax *= hardScaleFactor();
// 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 QTildeShowerHandler::setupHardScales(const vector<ShowerProgenitorPtr> & p,
XCPtr xcomb) {
if ( hardScaleIsMuF() &&
(!hardVetoReadOption() || firstInteraction()) ) {
Energy hardScale = ZERO;
if(currentTree()->isHard()) {
assert(xcomb);
hardScale = sqrt( xcomb->lastShowerScale() );
}
else {
hardScale = currentTree()->incomingLines().begin()->first
->progenitor()->momentum().mass();
}
hardScale *= hardScaleFactor();
vector<ShowerProgenitorPtr>::const_iterator ckt = p.begin();
for (; ckt != p.end(); ckt++) (*ckt)->hardScale(hardScale);
muPt = hardScale;
}
}
void QTildeShowerHandler::showerHardProcess(ShowerTreePtr hard, XCPtr xcomb) {
_hardme = HwMEBasePtr();
// extract the matrix element
tStdXCombPtr lastXC = dynamic_ptr_cast<tStdXCombPtr>(xcomb);
if(lastXC) {
_hardme = dynamic_ptr_cast<HwMEBasePtr>(lastXC->matrixElement());
}
_decayme = HwDecayerBasePtr();
// set the current tree
currentTree(hard);
hardTree(HardTreePtr());
// work out the type of event
currentTree()->xcombPtr(dynamic_ptr_cast<StdXCombPtr>(xcomb));
currentTree()->identifyEventType();
checkFlags();
// generate the showering
doShowering(true,xcomb);
}
RealEmissionProcessPtr QTildeShowerHandler::hardMatrixElementCorrection(bool hard) {
// set the initial enhancement factors for the soft correction
_initialenhance = 1.;
_finalenhance = 1.;
// see if we can get the correction from the matrix element
// or decayer
RealEmissionProcessPtr real;
if(hard) {
if(_hardme&&_hardme->hasMECorrection()) {
_hardme->initializeMECorrection(_currenttree->perturbativeProcess(),
_initialenhance,_finalenhance);
if(hardMEC())
real =
_hardme->applyHardMatrixElementCorrection(_currenttree->perturbativeProcess());
}
}
else {
if(_decayme&&_decayme->hasMECorrection()) {
_decayme->initializeMECorrection(_currenttree->perturbativeProcess(),
_initialenhance,_finalenhance);
if(hardMEC())
real = _decayme->applyHardMatrixElementCorrection(_currenttree->perturbativeProcess());
}
}
return real;
}
ShowerParticleVector QTildeShowerHandler::createTimeLikeChildren(tShowerParticlePtr, IdList ids) {
// 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 children;
for(unsigned int ix=0;ix<2;++ix) {
children.push_back(new_ptr(ShowerParticle(ids[ix+1],true)));
if(children[ix]->id()==_progenitor->id()&&!ids[ix+1]->stable()&&abs(ids[ix+1]->id())!=ParticleID::tauminus)
children[ix]->set5Momentum(Lorentz5Momentum(_progenitor->progenitor()->mass()));
else
children[ix]->set5Momentum(Lorentz5Momentum(ids[ix+1]->mass()));
}
return children;
}
bool QTildeShowerHandler::timeLikeShower(tShowerParticlePtr particle,
ShowerInteraction type,
Branching fb, bool first) {
// don't do anything if not needed
if(_limitEmissions == 1 || hardOnly() ||
( _limitEmissions == 2 && _nfs != 0) ||
( _limitEmissions == 4 && _nfs + _nis != 0) ) {
if(particle->spinInfo()) particle->spinInfo()->develop();
return false;
}
// generate the emission
ShowerParticleVector children;
// generate the emission
if(!fb.kinematics)
fb = selectTimeLikeBranching(particle,type,HardBranchingPtr());
// no emission, return
if(!fb.kinematics) {
if(particle->spinInfo()) particle->spinInfo()->develop();
return false;
}
Branching fc[2] = {Branching(),Branching()};
assert(fb.kinematics);
// has emitted
// Assign the shower kinematics to the emitting particle.
particle->showerKinematics(fb.kinematics);
// check highest pT
if(fb.kinematics->pT()>progenitor()->highestpT())
progenitor()->highestpT(fb.kinematics->pT());
// create the children
children = createTimeLikeChildren(particle,fb.ids);
// update the children
particle->showerKinematics()->
updateChildren(particle, children,fb.type);
// update number of emissions
++_nfs;
if(_limitEmissions!=0) {
if(children[0]->spinInfo()) children[0]->spinInfo()->develop();
if(children[1]->spinInfo()) children[1]->spinInfo()->develop();
if(particle->spinInfo()) particle->spinInfo()->develop();
return true;
}
// select branchings for children
fc[0] = selectTimeLikeBranching(children[0],type,HardBranchingPtr());
fc[1] = selectTimeLikeBranching(children[1],type,HardBranchingPtr());
// shower the first particle
if(fc[0].kinematics) timeLikeShower(children[0],type,fc[0],false);
if(children[0]->spinInfo()) children[0]->spinInfo()->develop();
// shower the second particle
if(fc[1].kinematics) timeLikeShower(children[1],type,fc[1],false);
if(children[1]->spinInfo()) children[1]->spinInfo()->develop();
particle->showerKinematics()->updateParent(particle, children,_evolutionScheme,fb.type);
// branching has happened
if(first&&!children.empty())
particle->showerKinematics()->resetChildren(particle,children);
if(particle->spinInfo()) particle->spinInfo()->develop();
return true;
}
bool
QTildeShowerHandler::spaceLikeShower(tShowerParticlePtr particle, PPtr beam,
ShowerInteraction 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(firstPDF().particle() == _beam)
pdf = firstPDF().pdf();
if(secondPDF().particle() == _beam)
pdf = secondPDF().pdf();
Energy freeze = pdfFreezingScale();
// don't do anything if not needed
if(_limitEmissions == 2 || hardOnly() ||
( _limitEmissions == 1 && _nis != 0 ) ||
( _limitEmissions == 4 && _nis + _nfs != 0 ) ) {
if(particle->spinInfo()) particle->spinInfo()->develop();
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) {
if(particle->spinInfo()) particle->spinInfo()->develop();
return false;
}
// if not vetoed break
if(!spaceLikeVetoed(bb,particle)) break;
// otherwise reset scale and continue
particle->vetoEmission(bb.type,bb.kinematics->scale());
if(particle->spinInfo()) particle->spinInfo()->decayVertex(VertexPtr());
}
// assign the splitting function and shower kinematics
particle->showerKinematics(bb.kinematics);
if(bb.kinematics->pT()>progenitor()->highestpT())
progenitor()->highestpT(bb.kinematics->pT());
// For the time being we are considering only 1->2 branching
// particles as in Sudakov form factor
tcPDPtr part[2]={bb.ids[0],bb.ids[2]};
// 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,_evolutionScheme,bb.type);
// 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;
if(newParent->spinInfo()) newParent->spinInfo()->develop();
// now reconstruct the momentum
if(!emitted) {
if(_intrinsic.find(_progenitor)==_intrinsic.end()) {
bb.kinematics->updateLast(newParent,ZERO,ZERO);
}
else {
pair<Energy,double> kt=_intrinsic[_progenitor];
bb.kinematics->updateLast(newParent,
kt.first*cos(kt.second),
kt.first*sin(kt.second));
}
}
particle->showerKinematics()->
updateChildren(newParent, theChildren,bb.type);
if(_limitEmissions!=0) {
if(particle->spinInfo()) particle->spinInfo()->develop();
return true;
}
// perform the shower of the final-state particle
timeLikeShower(otherChild,type,Branching(),true);
updateHistory(otherChild);
if(theChildren[1]->spinInfo()) theChildren[1]->spinInfo()->develop();
// return the emitted
if(particle->spinInfo()) particle->spinInfo()->develop();
return true;
}
void QTildeShowerHandler::showerDecay(ShowerTreePtr decay) {
// work out the type of event
currentTree()->xcombPtr(StdXCombPtr());
currentTree()->identifyEventType();
_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()
+ "->";
OrderedParticles outgoing;
for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
it=decay->outgoingLines().begin();it!=decay->outgoingLines().end();++it) {
if(abs(decay->incomingLines().begin()->first->original()->id()) == ParticleID::t &&
abs(it->first->original()->id())==ParticleID::Wplus &&
decay->treelinks().size() == 1) {
ShowerTreePtr Wtree = decay->treelinks().begin()->first;
for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
it2=Wtree->outgoingLines().begin();it2!=Wtree->outgoingLines().end();++it2) {
outgoing.insert(it2->first->original()->dataPtr());
}
}
else {
outgoing.insert(it->first->original()->dataPtr());
}
}
for(OrderedParticles::const_iterator it=outgoing.begin(); it!=outgoing.end();++it) {
if(it!=outgoing.begin()) tag += ",";
tag +=(**it).name();
}
tag += ";";
dm = findDecayMode(tag);
}
if(dm) _decayme = dynamic_ptr_cast<HwDecayerBasePtr>(dm->decayer());
// set the ShowerTree to be showered
currentTree(decay);
decay->applyTransforms();
hardTree(HardTreePtr());
// generate the showering
doShowering(false,XCPtr());
// if no vetos
// force calculation of spin correlations
SpinPtr spInfo = decay->incomingLines().begin()->first->progenitor()->spinInfo();
if(spInfo) {
if(!spInfo->developed()) spInfo->needsUpdate();
spInfo->develop();
}
}
bool QTildeShowerHandler::spaceLikeDecayShower(tShowerParticlePtr particle,
const ShowerParticle::EvolutionScales & maxScales,
Energy minmass,ShowerInteraction type,
Branching fb) {
// don't do anything if not needed
if(_limitEmissions == 1 || hardOnly() ||
( _limitEmissions == 3 && _nis != 0) ||
( _limitEmissions == 4 && _nfs + _nis != 0) ) {
return false;
}
// generate the emission
ShowerParticleVector children;
// generate the emission
if(!fb.kinematics)
fb = selectSpaceLikeDecayBranching(particle,maxScales,minmass,type,
HardBranchingPtr());
// no emission, return
if(!fb.kinematics) return false;
Branching fc[2];
if(particle->virtualMass()==ZERO)
particle->virtualMass(_progenitor->progenitor()->mass());
fc[0] = Branching();
fc[1] = Branching();
assert(fb.kinematics);
// has emitted
// Assign the shower kinematics to the emitting particle.
particle->showerKinematics(fb.kinematics);
if(fb.kinematics->pT()>progenitor()->highestpT())
progenitor()->highestpT(fb.kinematics->pT());
// create the ShowerParticle objects for the two children
children = createTimeLikeChildren(particle,fb.ids);
// updateChildren the children
particle->showerKinematics()->
updateChildren(particle, children, fb.type);
// select branchings for children
fc[0] = selectSpaceLikeDecayBranching(children[0],maxScales,minmass,
type,HardBranchingPtr());
fc[1] = selectTimeLikeBranching (children[1],type,HardBranchingPtr());
// old default
++_nis;
// shower the first particle
_currenttree->updateInitialStateShowerProduct(_progenitor,children[0]);
_currenttree->addInitialStateBranching(particle,children[0],children[1]);
if(fc[0].kinematics) spaceLikeDecayShower(children[0],maxScales,minmass,type,Branching());
// shower the second particle
if(fc[1].kinematics) timeLikeShower(children[1],type,fc[1],true);
updateHistory(children[1]);
// TODO NEED AN UPDATE HERE FOR RECONOPT!=0
// branching has happened
return true;
}
vector<ShowerProgenitorPtr> QTildeShowerHandler::setupShower(bool hard) {
RealEmissionProcessPtr real;
// generate hard me if needed
if(_hardEmission==1) {
real = hardMatrixElementCorrection(hard);
if(real&&!real->outgoing().empty()) setupMECorrection(real);
}
// generate POWHEG hard emission if needed
else if(_hardEmission==2)
hardestEmission(hard);
// set the initial colour partners
setEvolutionPartners(hard,interaction_,false);
// get the particles to be showered
vector<ShowerProgenitorPtr> particlesToShower =
currentTree()->extractProgenitors();
// return the answer
return particlesToShower;
}
void QTildeShowerHandler::setEvolutionPartners(bool hard,ShowerInteraction type,
bool clear) {
// match the particles in the ShowerTree and hardTree
if(hardTree() && !hardTree()->connect(currentTree()))
throw Exception() << "Can't match trees in "
<< "QTildeShowerHandler::setEvolutionPartners()"
<< Exception::eventerror;
// extract the progenitors
vector<ShowerParticlePtr> particles =
currentTree()->extractProgenitorParticles();
// clear the partners if needed
if(clear) {
for(unsigned int ix=0;ix<particles.size();++ix) {
particles[ix]->partner(ShowerParticlePtr());
particles[ix]->clearPartners();
}
}
// sort out the colour partners
if(hardTree()) {
// find the partner
for(unsigned int ix=0;ix<particles.size();++ix) {
tShowerParticlePtr partner = hardTree()->particles()[particles[ix]]->branchingParticle()->partner();
if(!partner) continue;
for(map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
it=hardTree()->particles().begin();
it!=hardTree()->particles().end();++it) {
if(it->second->branchingParticle()==partner) {
particles[ix]->partner(it->first);
break;
}
}
if(!particles[ix]->partner())
throw Exception() << "Can't match partners in "
<< "QTildeShowerHandler::setEvolutionPartners()"
<< Exception::eventerror;
}
}
// Set the initial evolution scales
partnerFinder()->
setInitialEvolutionScales(particles,!hard,interaction_,!_hardtree);
if(hardTree() && _hardPOWHEG) {
bool tooHard=false;
map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
eit=hardTree()->particles().end();
for(unsigned int ix=0;ix<particles.size();++ix) {
map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
mit = hardTree()->particles().find(particles[ix]);
Energy hardScale(ZERO);
ShowerPartnerType type(ShowerPartnerType::Undefined);
// final-state
if(particles[ix]->isFinalState()) {
if(mit!= eit && !mit->second->children().empty()) {
hardScale = mit->second->scale();
type = mit->second->type();
}
}
// initial-state
else {
if(mit!= eit && mit->second->parent()) {
hardScale = mit->second->parent()->scale();
type = mit->second->parent()->type();
}
}
if(type!=ShowerPartnerType::Undefined) {
if(type==ShowerPartnerType::QED) {
tooHard |= particles[ix]->scales().QED_noAO<hardScale;
}
else if(type==ShowerPartnerType::QCDColourLine) {
tooHard |= particles[ix]->scales().QCD_c_noAO<hardScale;
}
else if(type==ShowerPartnerType::QCDAntiColourLine) {
tooHard |= particles[ix]->scales().QCD_ac_noAO<hardScale;
}
}
}
if(tooHard) convertHardTree(hard,type);
}
}
void QTildeShowerHandler::updateHistory(tShowerParticlePtr particle) {
if(!particle->children().empty()) {
ShowerParticleVector theChildren;
for(unsigned int ix=0;ix<particle->children().size();++ix) {
ShowerParticlePtr part = dynamic_ptr_cast<ShowerParticlePtr>
(particle->children()[ix]);
theChildren.push_back(part);
}
// update the history if needed
if(particle==_currenttree->getFinalStateShowerProduct(_progenitor))
_currenttree->updateFinalStateShowerProduct(_progenitor,
particle,theChildren);
_currenttree->addFinalStateBranching(particle,theChildren);
for(unsigned int ix=0;ix<theChildren.size();++ix)
updateHistory(theChildren[ix]);
}
}
bool QTildeShowerHandler::startTimeLikeShower(ShowerInteraction type) {
// initialize basis vectors etc
if(!progenitor()->progenitor()->partner()) return false;
progenitor()->progenitor()->initializeFinalState();
if(hardTree()) {
map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
eit=hardTree()->particles().end(),
mit = hardTree()->particles().find(progenitor()->progenitor());
if( mit != eit && !mit->second->children().empty() ) {
bool output=truncatedTimeLikeShower(progenitor()->progenitor(),
mit->second ,type,Branching(),true);
if(output) updateHistory(progenitor()->progenitor());
return output;
}
}
// do the shower
bool output = hardOnly() ? false :
timeLikeShower(progenitor()->progenitor() ,type,Branching(),true) ;
if(output) updateHistory(progenitor()->progenitor());
return output;
}
bool QTildeShowerHandler::startSpaceLikeShower(PPtr parent, ShowerInteraction type) {
// initialise the basis vectors
if(!progenitor()->progenitor()->partner()) return false;
progenitor()->progenitor()->initializeInitialState(parent);
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 );
}
}
// perform the shower
return hardOnly() ? false :
spaceLikeShower(progenitor()->progenitor(),parent,type);
}
bool QTildeShowerHandler::
startSpaceLikeDecayShower(const ShowerParticle::EvolutionScales & maxScales,
Energy minimumMass,ShowerInteraction type) {
// set up the particle basis vectors
if(!progenitor()->progenitor()->partner()) return false;
progenitor()->progenitor()->initializeDecay();
if(hardTree()) {
map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
eit =hardTree()->particles().end(),
mit = hardTree()->particles().find(progenitor()->progenitor());
if( mit != eit && mit->second->parent() ) {
HardBranchingPtr branch=mit->second;
while(branch->parent()) branch=branch->parent();
return truncatedSpaceLikeDecayShower(progenitor()->progenitor(),maxScales,
minimumMass, branch ,type, Branching());
}
}
// perform the shower
return hardOnly() ? false :
spaceLikeDecayShower(progenitor()->progenitor(),maxScales,minimumMass,type,Branching());
}
bool QTildeShowerHandler::timeLikeVetoed(const Branching & fb,
ShowerParticlePtr particle) {
// work out type of interaction
ShowerInteraction type = convertInteraction(fb.type);
// check whether emission was harder than largest pt of hard subprocess
if ( restrictPhasespace() && fb.kinematics->pT() > _progenitor->maxHardPt() )
return true;
// soft matrix element correction veto
if( softMEC()) {
if(_hardme && _hardme->hasMECorrection()) {
if(_hardme->softMatrixElementVeto(particle,
_progenitor->progenitor(),
particle->isFinalState(),
_progenitor->highestpT(),
fb.ids, fb.kinematics->z(),
fb.kinematics->scale(),
fb.kinematics->pT()))
return true;
}
else if(_decayme && _decayme->hasMECorrection()) {
if(_decayme->softMatrixElementVeto(particle,
_progenitor->progenitor(),
particle->isFinalState(),
_progenitor->highestpT(),
fb.ids, fb.kinematics->z(),
fb.kinematics->scale(),
fb.kinematics->pT()))
return true;
}
}
// veto on maximum pt
if(fb.kinematics->pT()>_progenitor->maximumpT(type)) return true;
// general vetos
if (fb.kinematics && !_vetoes.empty()) {
bool vetoed=false;
for (vector<ShowerVetoPtr>::iterator v = _vetoes.begin();
v != _vetoes.end(); ++v) {
bool test = (**v).vetoTimeLike(_progenitor,particle,fb,currentTree());
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;
}
if ( firstInteraction() &&
profileScales() ) {
double weight =
profileScales()->
hardScaleProfile(_progenitor->hardScale(),fb.kinematics->pT());
if ( UseRandom::rnd() > weight )
return true;
}
return false;
}
bool QTildeShowerHandler::spaceLikeVetoed(const Branching & bb,
ShowerParticlePtr particle) {
// work out type of interaction
ShowerInteraction type = convertInteraction(bb.type);
// check whether emission was harder than largest pt of hard subprocess
if (restrictPhasespace() && bb.kinematics->pT() > _progenitor->maxHardPt())
return true;
// apply the soft correction
if( softMEC() && _hardme && _hardme->hasMECorrection() ) {
if(_hardme->softMatrixElementVeto(particle,
_progenitor->progenitor(),
particle->isFinalState(),
_progenitor->highestpT(),
bb.ids, bb.kinematics->z(),
bb.kinematics->scale(),
bb.kinematics->pT()))
return true;
}
// the more general vetos
// check vs max pt for the shower
if(bb.kinematics->pT()>_progenitor->maximumpT(type)) return true;
if (!_vetoes.empty()) {
bool vetoed=false;
for (vector<ShowerVetoPtr>::iterator v = _vetoes.begin();
v != _vetoes.end(); ++v) {
bool test = (**v).vetoSpaceLike(_progenitor,particle,bb,currentTree());
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;
}
if ( firstInteraction() &&
profileScales() ) {
double weight =
profileScales()->
hardScaleProfile(_progenitor->hardScale(),bb.kinematics->pT());
if ( UseRandom::rnd() > weight )
return true;
}
return false;
}
bool QTildeShowerHandler::spaceLikeDecayVetoed( const Branching & fb,
ShowerParticlePtr particle) {
// work out type of interaction
ShowerInteraction type = convertInteraction(fb.type);
// apply the soft correction
if( softMEC() && _decayme && _decayme->hasMECorrection() ) {
if(_decayme->softMatrixElementVeto(particle,
_progenitor->progenitor(),
particle->isFinalState(),
_progenitor->highestpT(),
fb.ids, fb.kinematics->z(),
fb.kinematics->scale(),
fb.kinematics->pT()))
return true;
}
// veto on hardest pt in the shower
if(fb.kinematics->pT()> _progenitor->maximumpT(type)) return true;
// general vetos
if (!_vetoes.empty()) {
bool vetoed=false;
for (vector<ShowerVetoPtr>::iterator v = _vetoes.begin();
v != _vetoes.end(); ++v) {
bool test = (**v).vetoSpaceLike(_progenitor,particle,fb,currentTree());
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 QTildeShowerHandler::hardestEmission(bool hard) {
HardTreePtr ISRTree;
// internal POWHEG in production or decay
if( (( _hardme && _hardme->hasPOWHEGCorrection()!=0 ) ||
( _decayme && _decayme->hasPOWHEGCorrection()!=0 ) ) ) {
RealEmissionProcessPtr real;
unsigned int type(0);
// production
if(_hardme) {
assert(hard);
real = _hardme->generateHardest( currentTree()->perturbativeProcess(),
interaction_);
type = _hardme->hasPOWHEGCorrection();
}
// decay
else {
assert(!hard);
real = _decayme->generateHardest( currentTree()->perturbativeProcess() );
type = _decayme->hasPOWHEGCorrection();
}
if(real) {
// set up ther hard tree
if(!real->outgoing().empty()) _hardtree = new_ptr(HardTree(real));
// set up the vetos
currentTree()->setVetoes(real->pT(),type);
}
// store initial state POWHEG radiation
if(_hardtree && _hardme && _hardme->hasPOWHEGCorrection()==1)
ISRTree = _hardtree;
}
else if (hard) {
// Get minimum pT cutoff used in shower approximation
Energy maxpt = 1.*GeV;
if ( currentTree()->showerApproximation() ) {
int colouredIn = 0;
int colouredOut = 0;
for( map< ShowerProgenitorPtr, tShowerParticlePtr >::iterator it
= currentTree()->outgoingLines().begin();
it != currentTree()->outgoingLines().end(); ++it ) {
if( it->second->coloured() ) ++colouredOut;
}
for( map< ShowerProgenitorPtr, ShowerParticlePtr >::iterator it
= currentTree()->incomingLines().begin();
it != currentTree()->incomingLines().end(); ++it ) {
if( it->second->coloured() ) ++colouredIn;
}
if ( currentTree()->showerApproximation()->ffPtCut() == currentTree()->showerApproximation()->fiPtCut() &&
currentTree()->showerApproximation()->ffPtCut() == currentTree()->showerApproximation()->iiPtCut() )
maxpt = currentTree()->showerApproximation()->ffPtCut();
else if ( colouredIn == 2 && colouredOut == 0 )
maxpt = currentTree()->showerApproximation()->iiPtCut();
else if ( colouredIn == 0 && colouredOut > 1 )
maxpt = currentTree()->showerApproximation()->ffPtCut();
else if ( colouredIn == 2 && colouredOut == 1 )
maxpt = min(currentTree()->showerApproximation()->iiPtCut(), currentTree()->showerApproximation()->fiPtCut());
else if ( colouredIn == 1 && colouredOut > 1 )
maxpt = min(currentTree()->showerApproximation()->ffPtCut(), currentTree()->showerApproximation()->fiPtCut());
else
maxpt = min(min(currentTree()->showerApproximation()->iiPtCut(), currentTree()->showerApproximation()->fiPtCut()),
currentTree()->showerApproximation()->ffPtCut());
}
// Generate hardtree from born and real emission subprocesses
_hardtree = generateCKKW(currentTree());
// Find transverse momentum of hardest emission
if (_hardtree){
for(set<HardBranchingPtr>::iterator it=_hardtree->branchings().begin();
it!=_hardtree->branchings().end();++it) {
if ((*it)->parent() && (*it)->status()==HardBranching::Incoming)
maxpt=(*it)->branchingParticle()->momentum().perp();
if ((*it)->children().size()==2 && (*it)->status()==HardBranching::Outgoing){
if ((*it)->branchingParticle()->id()!=21 &&
abs((*it)->branchingParticle()->id())>5 ){
if ((*it)->children()[0]->branchingParticle()->id()==21 ||
abs((*it)->children()[0]->branchingParticle()->id())<6)
maxpt=(*it)->children()[0]->branchingParticle()->momentum().perp();
else if ((*it)->children()[1]->branchingParticle()->id()==21 ||
abs((*it)->children()[1]->branchingParticle()->id())<6)
maxpt=(*it)->children()[1]->branchingParticle()->momentum().perp();
}
else {
if ( abs((*it)->branchingParticle()->id())<6){
if (abs((*it)->children()[0]->branchingParticle()->id())<6)
maxpt = (*it)->children()[1]->branchingParticle()->momentum().perp();
else
maxpt = (*it)->children()[0]->branchingParticle()->momentum().perp();
}
else maxpt = (*it)->children()[1]->branchingParticle()->momentum().perp();
}
}
}
}
// Hardest (pt) emission should be the first powheg emission.
maxpt=min(sqrt(lastXCombPtr()->lastShowerScale()),maxpt);
// set maximum pT for subsequent emissions from S events
if ( currentTree()->isPowhegSEvent() ) {
for( map< ShowerProgenitorPtr, tShowerParticlePtr >::iterator it
= currentTree()->outgoingLines().begin();
it != currentTree()->outgoingLines().end(); ++it ) {
if( ! it->second->coloured() ) continue;
it->first->maximumpT(maxpt, ShowerInteraction::QCD );
}
for( map< ShowerProgenitorPtr, ShowerParticlePtr >::iterator it
= currentTree()->incomingLines().begin();
it != currentTree()->incomingLines().end(); ++it ) {
if( ! it->second->coloured() ) continue;
it->first->maximumpT(maxpt, ShowerInteraction::QCD );
}
}
}
else
_hardtree = generateCKKW(currentTree());
// if hard me doesn't have a FSR powheg
// correction use decay powheg correction
if (_hardme && _hardme->hasPOWHEGCorrection()<2) {
addFSRUsingDecayPOWHEG(ISRTree);
}
// connect the trees
if(_hardtree) {
connectTrees(currentTree(),_hardtree,hard);
}
}
void QTildeShowerHandler::addFSRUsingDecayPOWHEG(HardTreePtr ISRTree) {
// check for intermediate colour singlet resonance
const ParticleVector inter = _hardme->subProcess()->intermediates();
if (inter.size()!=1 || inter[0]->momentum().m2()/GeV2 < 0 ||
inter[0]->dataPtr()->iColour()!=PDT::Colour0) {
return;
}
// ignore cases where outgoing particles are not coloured
map<ShowerProgenitorPtr, tShowerParticlePtr > out = currentTree()->outgoingLines();
if (out.size() != 2 ||
out. begin()->second->dataPtr()->iColour()==PDT::Colour0 ||
out.rbegin()->second->dataPtr()->iColour()==PDT::Colour0) {
return;
}
// look up decay mode
tDMPtr dm;
string tag;
string inParticle = inter[0]->dataPtr()->name() + "->";
vector<string> outParticles;
outParticles.push_back(out.begin ()->first->progenitor()->dataPtr()->name());
outParticles.push_back(out.rbegin()->first->progenitor()->dataPtr()->name());
for (int it=0; it<2; ++it){
tag = inParticle + outParticles[it] + "," + outParticles[(it+1)%2] + ";";
dm = generator()->findDecayMode(tag);
if(dm) break;
}
// get the decayer
HwDecayerBasePtr decayer;
if(dm) decayer = dynamic_ptr_cast<HwDecayerBasePtr>(dm->decayer());
// check if decayer has a FSR POWHEG correction
if (!decayer || decayer->hasPOWHEGCorrection()<2) {
return;
}
// generate the hardest emission
// create RealEmissionProcess
PPtr in = new_ptr(*inter[0]);
RealEmissionProcessPtr newProcess(new_ptr(RealEmissionProcess()));
newProcess->bornIncoming().push_back(in);
newProcess->bornOutgoing().push_back(out.begin ()->first->progenitor());
newProcess->bornOutgoing().push_back(out.rbegin()->first->progenitor());
// generate the FSR
newProcess = decayer->generateHardest(newProcess);
HardTreePtr FSRTree;
if(newProcess) {
// set up ther hard tree
if(!newProcess->outgoing().empty()) FSRTree = new_ptr(HardTree(newProcess));
// set up the vetos
currentTree()->setVetoes(newProcess->pT(),2);
}
if(!FSRTree) return;
// if there is no ISRTree make _hardtree from FSRTree
if (!ISRTree){
vector<HardBranchingPtr> inBranch,hardBranch;
for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
cit =currentTree()->incomingLines().begin();
cit!=currentTree()->incomingLines().end();++cit ) {
inBranch.push_back(new_ptr(HardBranching(cit->second,SudakovPtr(),
HardBranchingPtr(),
HardBranching::Incoming)));
inBranch.back()->beam(cit->first->original()->parents()[0]);
hardBranch.push_back(inBranch.back());
}
if(inBranch[0]->branchingParticle()->dataPtr()->coloured()) {
inBranch[0]->colourPartner(inBranch[1]);
inBranch[1]->colourPartner(inBranch[0]);
}
for(set<HardBranchingPtr>::iterator it=FSRTree->branchings().begin();
it!=FSRTree->branchings().end();++it) {
if((**it).branchingParticle()->id()!=in->id())
hardBranch.push_back(*it);
}
hardBranch[2]->colourPartner(hardBranch[3]);
hardBranch[3]->colourPartner(hardBranch[2]);
HardTreePtr newTree = new_ptr(HardTree(hardBranch,inBranch,
ShowerInteraction::QCD));
_hardtree = newTree;
}
// Otherwise modify the ISRTree to include the emission in FSRTree
else {
vector<tShowerParticlePtr> FSROut, ISROut;
set<HardBranchingPtr>::iterator itFSR, itISR;
// get outgoing particles
for(itFSR =FSRTree->branchings().begin();
itFSR!=FSRTree->branchings().end();++itFSR){
if ((**itFSR).status()==HardBranching::Outgoing)
FSROut.push_back((*itFSR)->branchingParticle());
}
for(itISR =ISRTree->branchings().begin();
itISR!=ISRTree->branchings().end();++itISR){
if ((**itISR).status()==HardBranching::Outgoing)
ISROut.push_back((*itISR)->branchingParticle());
}
// find COM frame formed by outgoing particles
LorentzRotation eventFrameFSR, eventFrameISR;
eventFrameFSR = ((FSROut[0]->momentum()+FSROut[1]->momentum()).findBoostToCM());
eventFrameISR = ((ISROut[0]->momentum()+ISROut[1]->momentum()).findBoostToCM());
// find rotation between ISR and FSR frames
int j=0;
if (ISROut[0]->id()!=FSROut[0]->id()) j=1;
eventFrameISR.rotateZ( (eventFrameFSR*FSROut[0]->momentum()).phi()-
(eventFrameISR*ISROut[j]->momentum()).phi() );
eventFrameISR.rotateY( (eventFrameFSR*FSROut[0]->momentum()).theta()-
(eventFrameISR*ISROut[j]->momentum()).theta() );
eventFrameISR.invert();
for (itFSR=FSRTree->branchings().begin();
itFSR!=FSRTree->branchings().end();++itFSR){
if ((**itFSR).branchingParticle()->id()==in->id()) continue;
for (itISR =ISRTree->branchings().begin();
itISR!=ISRTree->branchings().end();++itISR){
if ((**itISR).status()==HardBranching::Incoming) continue;
if ((**itFSR).branchingParticle()->id()==
(**itISR).branchingParticle()->id()){
// rotate FSRTree particle to ISRTree event frame
(**itISR).branchingParticle()->setMomentum(eventFrameISR*
eventFrameFSR*
(**itFSR).branchingParticle()->momentum());
(**itISR).branchingParticle()->rescaleMass();
// add the children of the FSRTree particles to the ISRTree
if(!(**itFSR).children().empty()){
(**itISR).addChild((**itFSR).children()[0]);
(**itISR).addChild((**itFSR).children()[1]);
// rotate momenta to ISRTree event frame
(**itISR).children()[0]->branchingParticle()->setMomentum(eventFrameISR*
eventFrameFSR*
(**itFSR).children()[0]->branchingParticle()->momentum());
(**itISR).children()[1]->branchingParticle()->setMomentum(eventFrameISR*
eventFrameFSR*
(**itFSR).children()[1]->branchingParticle()->momentum());
}
}
}
}
_hardtree = ISRTree;
}
}
bool QTildeShowerHandler::truncatedTimeLikeShower(tShowerParticlePtr particle,
HardBranchingPtr branch,
ShowerInteraction type,
Branching fb, bool first) {
// select a branching if we don't have one
if(!fb.kinematics)
fb = selectTimeLikeBranching(particle,type,branch);
// must be an emission, the forced one it not a truncated one
assert(fb.kinematics);
ShowerParticleVector children;
Branching fc[2] = {Branching(),Branching()};
// Assign the shower kinematics to the emitting particle.
particle->showerKinematics(fb.kinematics);
if(fb.kinematics->pT()>progenitor()->highestpT())
progenitor()->highestpT(fb.kinematics->pT());
// create the children
children = createTimeLikeChildren(particle,fb.ids);
// update the children
particle->showerKinematics()->
updateChildren(particle, children,fb.type);
// select branchings for children
if(!fc[0].kinematics) {
// select branching for first particle
if(!fb.hard && fb.iout ==1 )
fc[0] = selectTimeLikeBranching(children[0],type,branch);
else if(fb.hard && !branch->children()[0]->children().empty() )
fc[0] = selectTimeLikeBranching(children[0],type,branch->children()[0]);
else
fc[0] = selectTimeLikeBranching(children[0],type,HardBranchingPtr());
}
// select branching for the second particle
if(!fc[1].kinematics) {
// select branching for first particle
if(!fb.hard && fb.iout ==2 )
fc[1] = selectTimeLikeBranching(children[1],type,branch);
else if(fb.hard && !branch->children()[1]->children().empty() )
fc[1] = selectTimeLikeBranching(children[1],type,branch->children()[1]);
else
fc[1] = selectTimeLikeBranching(children[1],type,HardBranchingPtr());
}
// shower the first particle
if(fc[0].kinematics) {
// the parent has truncated emission and following line
if(!fb.hard && fb.iout == 1)
truncatedTimeLikeShower(children[0],branch,type,fc[0],false);
// hard emission and subsquent hard emissions
else if(fb.hard && !branch->children()[0]->children().empty() )
truncatedTimeLikeShower(children[0],branch->children()[0],type,fc[0],false);
// normal shower
else
timeLikeShower(children[0],type,fc[0],false);
}
if(children[0]->spinInfo()) children[0]->spinInfo()->develop();
// shower the second particle
if(fc[1].kinematics) {
// the parent has truncated emission and following line
if(!fb.hard && fb.iout == 2)
truncatedTimeLikeShower(children[1],branch,type,fc[1],false);
// hard emission and subsquent hard emissions
else if(fb.hard && !branch->children()[1]->children().empty() )
truncatedTimeLikeShower(children[1],branch->children()[1],type,fc[1],false);
else
timeLikeShower(children[1],type,fc[1],false);
}
if(children[1]->spinInfo()) children[1]->spinInfo()->develop();
// branching has happened
particle->showerKinematics()->updateParent(particle, children,_evolutionScheme,fb.type);
if(first&&!children.empty())
particle->showerKinematics()->resetChildren(particle,children);
if(particle->spinInfo()) particle->spinInfo()->develop();
// TODO NEED AN UPDATE HERE FOR RECONOPT!=0 ?
return true;
}
bool QTildeShowerHandler::truncatedSpaceLikeShower(tShowerParticlePtr particle, PPtr beam,
HardBranchingPtr branch,
ShowerInteraction type) {
tcPDFPtr pdf;
if(firstPDF().particle() == beamParticle())
pdf = firstPDF().pdf();
if(secondPDF().particle() == beamParticle())
pdf = secondPDF().pdf();
Energy freeze = pdfFreezingScale();
Branching bb;
// parameters of the force branching
double z(0.);
HardBranchingPtr timelike;
for( unsigned int ix = 0; ix < branch->children().size(); ++ix ) {
if( branch->children()[ix]->status() ==HardBranching::Outgoing) {
timelike = branch->children()[ix];
}
if( branch->children()[ix]->status() ==HardBranching::Incoming )
z = branch->children()[ix]->z();
}
// generate truncated branching
tcPDPtr part[2];
if(z>=0.&&z<=1.) {
while (true) {
if( !isTruncatedShowerON() || hardOnly() ) break;
bb = splittingGenerator()->chooseBackwardBranching( *particle,
beam, 1., beamParticle(),
type , pdf,freeze);
if( !bb.kinematics || bb.kinematics->scale() < branch->scale() ) {
bb = Branching();
break;
}
// particles as in Sudakov form factor
part[0] = bb.ids[0];
part[1] = bb.ids[2];
double zsplit = bb.kinematics->z();
// apply the vetos for the truncated shower
// if doesn't carry most of momentum
ShowerInteraction type2 = convertInteraction(bb.type);
if(type2==branch->sudakov()->interactionType() &&
zsplit < 0.5) {
particle->vetoEmission(bb.type,bb.kinematics->scale());
continue;
}
// others
if( part[0]->id() != particle->id() || // if particle changes type
bb.kinematics->pT() > progenitor()->maximumpT(type2) || // pt veto
bb.kinematics->scale() < branch->scale()) { // angular ordering veto
particle->vetoEmission(bb.type,bb.kinematics->scale());
continue;
}
// and those from the base class
if(spaceLikeVetoed(bb,particle)) {
particle->vetoEmission(bb.type,bb.kinematics->scale());
continue;
}
break;
}
}
if( !bb.kinematics ) {
//do the hard emission
ShoKinPtr kinematics = new_ptr(IS_QTildeShowerKinematics1to2(
branch->scale(), z, branch->phi(),
branch->children()[0]->pT(), branch->sudakov()
));
// assign the splitting function and shower kinematics
particle->showerKinematics( kinematics );
if(kinematics->pT()>progenitor()->highestpT())
progenitor()->highestpT(kinematics->pT());
// 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,_evolutionScheme, branch->type());
// update the history if needed
currentTree()->updateInitialStateShowerProduct( progenitor(), newParent );
currentTree()->addInitialStateBranching( particle, newParent, otherChild );
// for the reconstruction of kinematics, parent/child
// relationships are according to the branching process:
// now continue the shower
bool emitted=false;
if(!hardOnly()) {
if( branch->parent() ) {
emitted = truncatedSpaceLikeShower( newParent, beam, branch->parent() , type);
}
else {
emitted = spaceLikeShower( newParent, beam , type);
}
}
if( !emitted ) {
if( intrinsicpT().find( progenitor() ) == intrinsicpT().end() ) {
kinematics->updateLast( newParent, ZERO, ZERO );
}
else {
pair<Energy,double> kt = intrinsicpT()[progenitor()];
kinematics->updateLast( newParent,
kt.first*cos( kt.second ),
kt.first*sin( kt.second ) );
}
}
particle->showerKinematics()->
updateChildren( newParent, theChildren,bb.type);
if(hardOnly()) return true;
// perform the shower of the final-state particle
if( timelike->children().empty() ) {
timeLikeShower( otherChild , type,Branching(),true);
}
else {
truncatedTimeLikeShower( otherChild, timelike , type,Branching(), true);
}
updateHistory(otherChild);
// return the emitted
return true;
}
// assign the splitting function and shower kinematics
particle->showerKinematics( bb.kinematics );
if(bb.kinematics->pT()>progenitor()->highestpT())
progenitor()->highestpT(bb.kinematics->pT());
// 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,_evolutionScheme, bb.type);
// update the history if needed
currentTree()->updateInitialStateShowerProduct( progenitor(), newParent );
currentTree()->addInitialStateBranching( particle, newParent, otherChild );
// for the reconstruction of kinematics, parent/child
// relationships are according to the branching process:
// now continue the shower
bool emitted = truncatedSpaceLikeShower( newParent, beam, branch,type);
// now reconstruct the momentum
if( !emitted ) {
if( intrinsicpT().find( progenitor() ) == intrinsicpT().end() ) {
bb.kinematics->updateLast( newParent, ZERO, ZERO );
}
else {
pair<Energy,double> kt = intrinsicpT()[ progenitor() ];
bb.kinematics->updateLast( newParent,
kt.first*cos( kt.second ),
kt.first*sin( kt.second ) );
}
}
particle->showerKinematics()->
updateChildren( newParent, theChildren, bb.type);
// perform the shower of the final-state particle
timeLikeShower( otherChild , type,Branching(),true);
updateHistory(otherChild);
// return the emitted
return true;
}
bool QTildeShowerHandler::
truncatedSpaceLikeDecayShower(tShowerParticlePtr particle,
const ShowerParticle::EvolutionScales & maxScales,
Energy minmass, HardBranchingPtr branch,
ShowerInteraction type, Branching fb) {
// select a branching if we don't have one
if(!fb.kinematics)
fb = selectSpaceLikeDecayBranching(particle,maxScales,minmass,type,branch);
// must be an emission, the forced one it not a truncated one
assert(fb.kinematics);
ShowerParticleVector children;
Branching fc[2]={Branching(),Branching()};
// Assign the shower kinematics to the emitting particle.
particle->showerKinematics(fb.kinematics);
if(fb.kinematics->pT()>progenitor()->highestpT())
progenitor()->highestpT(fb.kinematics->pT());
// create the ShowerParticle objects for the two children
children = createTimeLikeChildren(particle,fb.ids);
// updateChildren the children
particle->showerKinematics()->
updateChildren(particle, children, fb.type);
// select branchings for children
if(!fc[0].kinematics) {
if(children[0]->id()==particle->id()) {
// select branching for first particle
if(!fb.hard)
fc[0] = selectSpaceLikeDecayBranching(children[0],maxScales,minmass,type,branch);
else if(fb.hard && ! branch->children()[0]->children().empty() )
fc[0] = selectSpaceLikeDecayBranching(children[0],maxScales,minmass,type,
branch->children()[0]);
else
fc[0] = selectSpaceLikeDecayBranching(children[0],maxScales,minmass,type,
HardBranchingPtr());
}
else {
// select branching for first particle
if(fb.hard && !branch->children()[0]->children().empty() )
fc[0] = selectTimeLikeBranching(children[0],type,branch->children()[0]);
else
fc[0] = selectTimeLikeBranching(children[0],type,HardBranchingPtr());
}
}
// select branching for the second particle
if(!fc[1].kinematics) {
if(children[1]->id()==particle->id()) {
// select branching for first particle
if(!fb.hard)
fc[1] = selectSpaceLikeDecayBranching(children[1],maxScales,minmass,type,branch);
else if(fb.hard && ! branch->children()[1]->children().empty() )
fc[1] = selectSpaceLikeDecayBranching(children[1],maxScales,minmass,type,
branch->children()[1]);
else
fc[1] = selectSpaceLikeDecayBranching(children[1],maxScales,minmass,type,
HardBranchingPtr());
}
else {
if(fb.hard && !branch->children()[1]->children().empty() )
fc[1] = selectTimeLikeBranching(children[1],type,branch->children()[1]);
else
fc[1] = selectTimeLikeBranching(children[1],type,HardBranchingPtr());
}
}
// old default
// update the history if needed
currentTree()->updateInitialStateShowerProduct(progenitor(),children[0]);
currentTree()->addInitialStateBranching(particle,children[0],children[1]);
// shower the first particle
if(fc[0].kinematics) {
if(children[0]->id()==particle->id()) {
if(!fb.hard)
truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
branch,type,fc[0]);
else if(fb.hard && ! branch->children()[0]->children().empty() )
truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
branch->children()[0],type,fc[0]);
else
spaceLikeDecayShower( children[0],maxScales,minmass,type,fc[0]);
}
else {
if(fb.hard && !branch->children()[0]->children().empty() )
truncatedTimeLikeShower(children[0],branch->children()[0],type,fc[0],false);
// normal shower
else
timeLikeShower(children[0],type,fc[0],false);
}
}
// shower the second particle
if(fc[1].kinematics) {
if(children[0]->id()==particle->id()) {
if(!fb.hard)
truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
branch,type,fc[1]);
else if(fb.hard && ! branch->children()[0]->children().empty() )
truncatedSpaceLikeDecayShower( children[0],maxScales,minmass,
branch->children()[0],type,fc[1]);
else
spaceLikeDecayShower( children[0],maxScales,minmass,type,fc[1]);
}
else {
if(fb.hard && !branch->children()[0]->children().empty() )
truncatedTimeLikeShower(children[0],branch->children()[0],type,fc[1],false);
// normal shower
else
timeLikeShower(children[0],type,fc[1],false);
}
}
updateHistory(children[1]);
// TODO DO WE NEED A CHECK IF RECONPT !=0
return true;
}
void QTildeShowerHandler::connectTrees(ShowerTreePtr showerTree,
HardTreePtr hardTree, bool hard ) {
ShowerParticleVector particles;
// find the Sudakovs
for(set<HardBranchingPtr>::iterator cit=hardTree->branchings().begin();
cit!=hardTree->branchings().end();++cit) {
// Sudakovs for ISR
if((**cit).parent()&&(**cit).status()==HardBranching::Incoming) {
++_nis;
array<long,3> br;
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.particles;
if(ids[0]->id()==br[0]&&ids[1]->id()==br[1]&&ids[2]->id()==br[2]) {
sudakov=cjt->second.sudakov;
break;
}
}
if(!sudakov) throw Exception() << "Can't find Sudakov for the hard emission in "
<< "QTildeShowerHandler::connectTrees() for ISR"
<< Exception::runerror;
(**cit).parent()->sudakov(sudakov);
}
// Sudakovs for FSR
else if(!(**cit).children().empty()) {
++_nfs;
array<long,3> br;
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.particles;
if(ids[0]->id()==br[0]&&ids[1]->id()==br[1]&&ids[2]->id()==br[2]) {
sudakov=cjt->second.sudakov;
break;
}
}
if(!sudakov) {
throw Exception() << "Can't find Sudakov for the hard emission in "
<< "QTildeShowerHandler::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());
}
partnerFinder()->
setInitialEvolutionScales(particles,!hard,interaction_,true);
hardTree->partnersSet(true);
// inverse reconstruction
if(hard) {
kinematicsReconstructor()->
deconstructHardJets(hardTree,interaction_);
}
else
kinematicsReconstructor()->
deconstructDecayJets(hardTree,interaction_);
// now reset the momenta of the showering particles
vector<ShowerProgenitorPtr> particlesToShower=showerTree->extractProgenitors();
// match them
map<ShowerProgenitorPtr,HardBranchingPtr> partners;
for(set<HardBranchingPtr>::const_iterator bit=hardTree->branchings().begin();
bit!=hardTree->branchings().end();++bit) {
Energy2 dmin( 1e30*GeV2 );
ShowerProgenitorPtr partner;
for(vector<ShowerProgenitorPtr>::const_iterator pit=particlesToShower.begin();
pit!=particlesToShower.end();++pit) {
if(partners.find(*pit)!=partners.end()) continue;
if( (**bit).branchingParticle()->id() != (**pit).progenitor()->id() ) continue;
if( (**bit).branchingParticle()->isFinalState() !=
(**pit).progenitor()->isFinalState() ) continue;
if( (**pit).progenitor()->isFinalState() ) {
Energy2 dtest =
sqr( (**pit).progenitor()->momentum().x() - (**bit).showerMomentum().x() ) +
sqr( (**pit).progenitor()->momentum().y() - (**bit).showerMomentum().y() ) +
sqr( (**pit).progenitor()->momentum().z() - (**bit).showerMomentum().z() ) +
sqr( (**pit).progenitor()->momentum().t() - (**bit).showerMomentum().t() );
// add mass difference for identical particles (e.g. Z0 Z0 production)
dtest += 1e10*sqr((**pit).progenitor()->momentum().m()-(**bit).showerMomentum().m());
if( dtest < dmin ) {
partner = *pit;
dmin = dtest;
}
}
else {
// ensure directions are right
if((**pit).progenitor()->momentum().z()/(**bit).showerMomentum().z()>ZERO) {
partner = *pit;
break;
}
}
}
if(!partner) throw Exception() << "Failed to match shower and hard trees in QTildeShowerHandler::hardestEmission"
<< Exception::eventerror;
partners[partner] = *bit;
}
for(vector<ShowerProgenitorPtr>::const_iterator pit=particlesToShower.begin();
pit!=particlesToShower.end();++pit) {
HardBranchingPtr partner = partners[*pit];
if((**pit).progenitor()->dataPtr()->stable()) {
(**pit).progenitor()->set5Momentum(partner->showerMomentum());
(**pit).copy()->set5Momentum(partner->showerMomentum());
}
else {
Lorentz5Momentum oldMomentum = (**pit).progenitor()->momentum();
Lorentz5Momentum newMomentum = partner->showerMomentum();
LorentzRotation boost( oldMomentum.findBoostToCM(),oldMomentum.e()/oldMomentum.mass());
(**pit).progenitor()->transform(boost);
(**pit).copy() ->transform(boost);
boost = LorentzRotation(-newMomentum.findBoostToCM(),newMomentum.e()/newMomentum.mass());
(**pit).progenitor()->transform(boost);
(**pit).copy() ->transform(boost);
}
}
// 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());
decayTree->transform(boost,true);
boost = LorentzRotation(-newMomentum.findBoostToCM(),newMomentum.e()/newMomentum.mass());
decayTree->transform(boost,true);
}
}
void QTildeShowerHandler::doShowering(bool hard,XCPtr xcomb) {
// zero number of emissions
_nis = _nfs = 0;
// if MC@NLO H event and limited emissions
// indicate both final and initial state emission
if ( currentTree()->isMCatNLOHEvent() && _limitEmissions != 0 ) {
_nis = _nfs = 1;
}
// extract particles to shower
vector<ShowerProgenitorPtr> particlesToShower(setupShower(hard));
// check if we should shower
bool colCharge = false;
for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
if(particlesToShower[ix]->progenitor()->dataPtr()->coloured() ||
particlesToShower[ix]->progenitor()->dataPtr()->charged()) {
colCharge = true;
break;
}
}
if(!colCharge) {
_currenttree->hasShowered(true);
return;
}
// setup the maximum scales for the shower
if (restrictPhasespace()) setupMaximumScales(particlesToShower,xcomb);
// set the hard scales for the profiles
setupHardScales(particlesToShower,xcomb);
// specific stuff for hard processes and decays
Energy minmass(ZERO), mIn(ZERO);
// hard process generate the intrinsic p_T once and for all
if(hard) {
generateIntrinsicpT(particlesToShower);
}
// decay compute the minimum mass of the final-state
else {
for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
if(particlesToShower[ix]->progenitor()->isFinalState()) {
if(particlesToShower[ix]->progenitor()->dataPtr()->stable()){
auto dm= ShowerHandler::currentHandler()->retConstituentMasses()?
particlesToShower[ix]->progenitor()->dataPtr()->constituentMass():
particlesToShower[ix]->progenitor()->dataPtr()->mass();
minmass += dm;
}else
minmass += particlesToShower[ix]->progenitor()->mass();
}
else {
mIn = particlesToShower[ix]->progenitor()->mass();
}
}
// throw exception if decay can't happen
if ( minmass > mIn ) {
throw Exception() << "QTildeShowerHandler.cc: Mass of decaying particle is "
<< "below constituent masses of decay products."
<< Exception::eventerror;
}
}
// setup for reweighted
bool reWeighting = _reWeight && hard && ShowerHandler::currentHandler()->firstInteraction();
double eventWeight=0.;
unsigned int nTryReWeight(0);
// create random particle vector (only need to do once)
vector<ShowerProgenitorPtr> tmp;
unsigned int nColouredIncoming = 0;
while(particlesToShower.size()>0){
unsigned int xx=UseRandom::irnd(particlesToShower.size());
tmp.push_back(particlesToShower[xx]);
particlesToShower.erase(particlesToShower.begin()+xx);
}
particlesToShower=tmp;
for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
if(!particlesToShower[ix]->progenitor()->isFinalState() &&
particlesToShower[ix]->progenitor()->coloured()) ++nColouredIncoming;
}
bool switchRecon = hard && nColouredIncoming !=1;
// main shower loop
unsigned int ntry(0);
bool reconstructed = false;
do {
// clear results of last attempt if needed
if(ntry!=0) {
currentTree()->clear();
setEvolutionPartners(hard,interaction_,true);
_nis = _nfs = 0;
// if MC@NLO H event and limited emissions
// indicate both final and initial state emission
if ( currentTree()->isMCatNLOHEvent() && _limitEmissions != 0 ) {
_nis = _nfs = 1;
}
for(unsigned int ix=0; ix<particlesToShower.size();++ix) {
SpinPtr spin = particlesToShower[ix]->progenitor()->spinInfo();
if(spin && spin->decayVertex() &&
dynamic_ptr_cast<tcSVertexPtr>(spin->decayVertex())) {
spin->decayVertex(VertexPtr());
}
}
for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
if(particlesToShower[ix]->progenitor()->isFinalState() ||
(hard && !particlesToShower[ix]->progenitor()->isFinalState())) {
if(particlesToShower[ix]->progenitor()->spinInfo())
particlesToShower[ix]->progenitor()->spinInfo()->reset();
}
}
}
// loop over particles
for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
// extract the progenitor
progenitor(particlesToShower[ix]);
// final-state radiation
if(progenitor()->progenitor()->isFinalState()) {
if(!doFSR()) continue;
// perform shower
progenitor()->hasEmitted(startTimeLikeShower(interaction_));
}
// initial-state radiation
else {
if(!doISR()) continue;
// hard process
if(hard) {
// get the PDF
setBeamParticle(_progenitor->beam());
if(!beamParticle()) {
throw Exception() << "Incorrect type of beam particle in "
<< "QTildeShowerHandler::doShowering(). "
<< "This should not happen for conventional choices but may happen if you have used a"
<< " non-default choice and have not changed the create ParticleData line in the input files"
<< " for this particle to create BeamParticleData."
<< Exception::runerror;
}
// 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,
interaction_));
}
// decay
else {
// skip colour and electrically neutral particles
if(!progenitor()->progenitor()->dataPtr()->coloured() &&
!progenitor()->progenitor()->dataPtr()->charged()) {
progenitor()->hasEmitted(false);
continue;
}
// perform shower
// set the scales correctly. The current scale is the maximum scale for
// emission not the starting scale
ShowerParticle::EvolutionScales maxScales(progenitor()->progenitor()->scales());
progenitor()->progenitor()->scales() = ShowerParticle::EvolutionScales();
if(progenitor()->progenitor()->dataPtr()->charged()) {
progenitor()->progenitor()->scales().QED = progenitor()->progenitor()->mass();
progenitor()->progenitor()->scales().QED_noAO = progenitor()->progenitor()->mass();
}
if(progenitor()->progenitor()->hasColour()) {
progenitor()->progenitor()->scales().QCD_c = progenitor()->progenitor()->mass();
progenitor()->progenitor()->scales().QCD_c_noAO = progenitor()->progenitor()->mass();
}
if(progenitor()->progenitor()->hasAntiColour()) {
progenitor()->progenitor()->scales().QCD_ac = progenitor()->progenitor()->mass();
progenitor()->progenitor()->scales().QCD_ac_noAO = progenitor()->progenitor()->mass();
}
// perform the shower
progenitor()->hasEmitted(startSpaceLikeDecayShower(maxScales,minmass,
interaction_));
}
}
}
// do the kinematic reconstruction, checking if it worked
reconstructed = hard ?
kinematicsReconstructor()->
reconstructHardJets (currentTree(),intrinsicpT(),interaction_,
switchRecon && ntry>maximumTries()/2) :
kinematicsReconstructor()->
reconstructDecayJets(currentTree(),interaction_);
if(!reconstructed) continue;
// apply vetos on the full shower
for(vector<FullShowerVetoPtr>::const_iterator it=_fullShowerVetoes.begin();
it!=_fullShowerVetoes.end();++it) {
int veto = (**it).applyVeto(currentTree());
if(veto<0) continue;
// veto the shower
if(veto==0) {
reconstructed = false;
break;
}
// veto the shower and reweight
else if(veto==1) {
reconstructed = false;
break;
}
// veto the event
else if(veto==2) {
throw Veto();
}
}
if(reWeighting) {
if(reconstructed) eventWeight += 1.;
reconstructed=false;
++nTryReWeight;
if(nTryReWeight==_nReWeight) {
reWeighting = false;
if(eventWeight==0.) throw Veto();
}
}
}
while(!reconstructed&&maximumTries()>++ntry);
// check if failed to generate the shower
if(ntry==maximumTries()) {
if(hard)
throw ShowerHandler::ShowerTriesVeto(ntry);
else
throw Exception() << "Failed to generate the shower after "
<< ntry << " attempts in QTildeShowerHandler::showerDecay()"
<< Exception::eventerror;
}
// handle the weights and apply any reweighting required
if(nTryReWeight>0) {
tStdEHPtr seh = dynamic_ptr_cast<tStdEHPtr>(generator()->currentEventHandler());
static bool first = true;
if(seh) {
seh->reweight(eventWeight/double(nTryReWeight));
}
else if(first) {
generator()->log() << "Reweighting the shower only works with internal Herwig7 processes"
<< "Presumably you are showering Les Houches Events. These will not be"
<< "reweighted\n";
first = false;
}
}
// tree has now showered
_currenttree->hasShowered(true);
hardTree(HardTreePtr());
}
void QTildeShowerHandler:: convertHardTree(bool hard,ShowerInteraction type) {
map<ColinePtr,ColinePtr> cmap;
// incoming particles
for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
cit=currentTree()->incomingLines().begin();cit!=currentTree()->incomingLines().end();++cit) {
map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
mit = hardTree()->particles().find(cit->first->progenitor());
// put the colour lines in the map
ShowerParticlePtr oldParticle = cit->first->progenitor();
ShowerParticlePtr newParticle = mit->second->branchingParticle();
ColinePtr cLine = oldParticle-> colourLine();
ColinePtr aLine = oldParticle->antiColourLine();
if(newParticle->colourLine() &&
cmap.find(newParticle-> colourLine())==cmap.end())
cmap[newParticle-> colourLine()] = cLine;
if(newParticle->antiColourLine() &&
cmap.find(newParticle->antiColourLine())==cmap.end())
cmap[newParticle->antiColourLine()] = aLine;
// check whether or not particle emits
bool emission = mit->second->parent();
if(emission) {
if(newParticle->colourLine()) {
ColinePtr ctemp = newParticle-> colourLine();
ctemp->removeColoured(newParticle);
}
if(newParticle->antiColourLine()) {
ColinePtr ctemp = newParticle->antiColourLine();
ctemp->removeAntiColoured(newParticle);
}
newParticle = mit->second->parent()->branchingParticle();
}
// get the new colour lines
ColinePtr newCLine,newALine;
// sort out colour lines
if(newParticle->colourLine()) {
ColinePtr ctemp = newParticle-> colourLine();
ctemp->removeColoured(newParticle);
if(cmap.find(ctemp)!=cmap.end()) {
newCLine = cmap[ctemp];
}
else {
newCLine = new_ptr(ColourLine());
cmap[ctemp] = newCLine;
}
}
// and anticolour lines
if(newParticle->antiColourLine()) {
ColinePtr ctemp = newParticle->antiColourLine();
ctemp->removeAntiColoured(newParticle);
if(cmap.find(ctemp)!=cmap.end()) {
newALine = cmap[ctemp];
}
else {
newALine = new_ptr(ColourLine());
cmap[ctemp] = newALine;
}
}
// remove colour lines from old particle
if(aLine) {
aLine->removeAntiColoured(cit->first->copy());
aLine->removeAntiColoured(cit->first->progenitor());
}
if(cLine) {
cLine->removeColoured(cit->first->copy());
cLine->removeColoured(cit->first->progenitor());
}
// add particle to colour lines
if(newCLine) newCLine->addColoured (newParticle);
if(newALine) newALine->addAntiColoured(newParticle);
// insert new particles
cit->first->copy(newParticle);
ShowerParticlePtr sp(new_ptr(ShowerParticle(*newParticle,1,false)));
cit->first->progenitor(sp);
currentTree()->incomingLines()[cit->first]=sp;
cit->first->perturbative(!emission);
// and the emitted particle if needed
if(emission) {
ShowerParticlePtr newOut = mit->second->parent()->children()[1]->branchingParticle();
if(newOut->colourLine()) {
ColinePtr ctemp = newOut-> colourLine();
ctemp->removeColoured(newOut);
assert(cmap.find(ctemp)!=cmap.end());
cmap[ctemp]->addColoured (newOut);
}
if(newOut->antiColourLine()) {
ColinePtr ctemp = newOut->antiColourLine();
ctemp->removeAntiColoured(newOut);
assert(cmap.find(ctemp)!=cmap.end());
cmap[ctemp]->addAntiColoured(newOut);
}
ShowerParticlePtr sout=new_ptr(ShowerParticle(*newOut,1,true));
ShowerProgenitorPtr out=new_ptr(ShowerProgenitor(cit->first->original(),newOut,sout));
out->perturbative(false);
currentTree()->outgoingLines().insert(make_pair(out,sout));
}
if(hard) {
// sort out the value of x
if(mit->second->beam()->momentum().z()>ZERO) {
sp->x(newParticle->momentum(). plus()/mit->second->beam()->momentum(). plus());
}
else {
sp->x(newParticle->momentum().minus()/mit->second->beam()->momentum().minus());
}
}
}
// outgoing particles
for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
cit=currentTree()->outgoingLines().begin();cit!=currentTree()->outgoingLines().end();++cit) {
map<tShowerTreePtr,pair<tShowerProgenitorPtr,
tShowerParticlePtr> >::const_iterator tit;
for(tit = currentTree()->treelinks().begin();
tit != currentTree()->treelinks().end();++tit) {
if(tit->second.first && tit->second.second==cit->first->progenitor())
break;
}
map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
mit = hardTree()->particles().find(cit->first->progenitor());
if(mit==hardTree()->particles().end()) continue;
// put the colour lines in the map
ShowerParticlePtr oldParticle = cit->first->progenitor();
ShowerParticlePtr newParticle = mit->second->branchingParticle();
ShowerParticlePtr newOut;
ColinePtr cLine = oldParticle-> colourLine();
ColinePtr aLine = oldParticle->antiColourLine();
if(newParticle->colourLine() &&
cmap.find(newParticle-> colourLine())==cmap.end())
cmap[newParticle-> colourLine()] = cLine;
if(newParticle->antiColourLine() &&
cmap.find(newParticle->antiColourLine())==cmap.end())
cmap[newParticle->antiColourLine()] = aLine;
// check whether or not particle emits
bool emission = !mit->second->children().empty();
if(emission) {
if(newParticle->colourLine()) {
ColinePtr ctemp = newParticle-> colourLine();
ctemp->removeColoured(newParticle);
}
if(newParticle->antiColourLine()) {
ColinePtr ctemp = newParticle->antiColourLine();
ctemp->removeAntiColoured(newParticle);
}
newParticle = mit->second->children()[0]->branchingParticle();
newOut = mit->second->children()[1]->branchingParticle();
if(newParticle->id()!=oldParticle->id()&&newParticle->id()==newOut->id())
swap(newParticle,newOut);
}
// get the new colour lines
ColinePtr newCLine,newALine;
// sort out colour lines
if(newParticle->colourLine()) {
ColinePtr ctemp = newParticle-> colourLine();
ctemp->removeColoured(newParticle);
if(cmap.find(ctemp)!=cmap.end()) {
newCLine = cmap[ctemp];
}
else {
newCLine = new_ptr(ColourLine());
cmap[ctemp] = newCLine;
}
}
// and anticolour lines
if(newParticle->antiColourLine()) {
ColinePtr ctemp = newParticle->antiColourLine();
ctemp->removeAntiColoured(newParticle);
if(cmap.find(ctemp)!=cmap.end()) {
newALine = cmap[ctemp];
}
else {
newALine = new_ptr(ColourLine());
cmap[ctemp] = newALine;
}
}
// remove colour lines from old particle
if(aLine) {
aLine->removeAntiColoured(cit->first->copy());
aLine->removeAntiColoured(cit->first->progenitor());
}
if(cLine) {
cLine->removeColoured(cit->first->copy());
cLine->removeColoured(cit->first->progenitor());
}
// special for unstable particles
if(newParticle->id()==oldParticle->id() &&
(tit!=currentTree()->treelinks().end()||!oldParticle->dataPtr()->stable())) {
Lorentz5Momentum oldMomentum = oldParticle->momentum();
Lorentz5Momentum newMomentum = newParticle->momentum();
LorentzRotation boost( oldMomentum.findBoostToCM(),oldMomentum.e()/oldMomentum.mass());
if(tit!=currentTree()->treelinks().end()) tit->first->transform(boost,false);
oldParticle->transform(boost);
boost = LorentzRotation(-newMomentum.findBoostToCM(),newMomentum.e()/newMomentum.mass());
oldParticle->transform(boost);
if(tit!=currentTree()->treelinks().end()) tit->first->transform(boost,false);
newParticle=oldParticle;
}
// add particle to colour lines
if(newCLine) newCLine->addColoured (newParticle);
if(newALine) newALine->addAntiColoured(newParticle);
// insert new particles
cit->first->copy(newParticle);
ShowerParticlePtr sp(new_ptr(ShowerParticle(*newParticle,1,true)));
cit->first->progenitor(sp);
currentTree()->outgoingLines()[cit->first]=sp;
cit->first->perturbative(!emission);
// and the emitted particle if needed
if(emission) {
if(newOut->colourLine()) {
ColinePtr ctemp = newOut-> colourLine();
ctemp->removeColoured(newOut);
assert(cmap.find(ctemp)!=cmap.end());
cmap[ctemp]->addColoured (newOut);
}
if(newOut->antiColourLine()) {
ColinePtr ctemp = newOut->antiColourLine();
ctemp->removeAntiColoured(newOut);
assert(cmap.find(ctemp)!=cmap.end());
cmap[ctemp]->addAntiColoured(newOut);
}
ShowerParticlePtr sout=new_ptr(ShowerParticle(*newOut,1,true));
ShowerProgenitorPtr out=new_ptr(ShowerProgenitor(cit->first->original(),newOut,sout));
out->perturbative(false);
currentTree()->outgoingLines().insert(make_pair(out,sout));
}
// update any decay products
if(tit!=currentTree()->treelinks().end())
currentTree()->updateLink(tit->first,make_pair(cit->first,sp));
}
// reset the tree
currentTree()->resetShowerProducts();
// reextract the particles and set the colour partners
vector<ShowerParticlePtr> particles =
currentTree()->extractProgenitorParticles();
// clear the partners
for(unsigned int ix=0;ix<particles.size();++ix) {
particles[ix]->partner(ShowerParticlePtr());
particles[ix]->clearPartners();
}
// clear the tree
hardTree(HardTreePtr());
// Set the initial evolution scales
partnerFinder()->
setInitialEvolutionScales(particles,!hard,type,!_hardtree);
}
Branching QTildeShowerHandler::selectTimeLikeBranching(tShowerParticlePtr particle,
ShowerInteraction type,
HardBranchingPtr branch) {
Branching fb;
unsigned int iout=0;
while (true) {
// break if doing truncated shower and no truncated shower needed
if(branch && (!isTruncatedShowerON()||hardOnly())) break;
fb=_splittingGenerator->chooseForwardBranching(*particle,_finalenhance,type);
// no emission break
if(!fb.kinematics) break;
// special for truncated shower
if(branch) {
// check haven't evolved too far
if(fb.kinematics->scale() < branch->scale()) {
fb=Branching();
break;
}
// find the truncated line
iout=0;
if(fb.ids[1]->id()!=fb.ids[2]->id()) {
if(fb.ids[1]->id()==particle->id()) iout=1;
else if (fb.ids[2]->id()==particle->id()) iout=2;
}
else if(fb.ids[1]->id()==particle->id()) {
if(fb.kinematics->z()>0.5) iout=1;
else iout=2;
}
// apply the vetos for the truncated shower
// no flavour changing branchings
if(iout==0) {
particle->vetoEmission(fb.type,fb.kinematics->scale());
continue;
}
double zsplit = iout==1 ? fb.kinematics->z() : 1-fb.kinematics->z();
// only if same interaction for forced branching
ShowerInteraction type2 = convertInteraction(fb.type);
// and evolution
if(type2==branch->sudakov()->interactionType()) {
if(zsplit < 0.5 || // hardest line veto
fb.kinematics->scale()*zsplit < branch->scale() ) { // angular ordering veto
particle->vetoEmission(fb.type,fb.kinematics->scale());
continue;
}
}
// pt veto
if(fb.kinematics->pT() > progenitor()->maximumpT(type2)) {
particle->vetoEmission(fb.type,fb.kinematics->scale());
continue;
}
}
// standard vetos for all emissions
if(timeLikeVetoed(fb,particle)) {
particle->vetoEmission(fb.type,fb.kinematics->scale());
if(particle->spinInfo()) particle->spinInfo()->decayVertex(VertexPtr());
continue;
}
// special for already decayed particles
// don't allow flavour changing branchings
bool vetoDecay = false;
for(map<tShowerTreePtr,pair<tShowerProgenitorPtr,
tShowerParticlePtr> >::const_iterator tit = currentTree()->treelinks().begin();
tit != currentTree()->treelinks().end();++tit) {
if(tit->second.first == progenitor()) {
map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
it = currentTree()->outgoingLines().find(progenitor());
if(it!=currentTree()->outgoingLines().end() && particle == it->second &&
fb.ids[0]!=fb.ids[1] && fb.ids[1]!=fb.ids[2]) {
vetoDecay = true;
break;
}
}
}
if(vetoDecay) {
particle->vetoEmission(fb.type,fb.kinematics->scale());
if(particle->spinInfo()) particle->spinInfo()->decayVertex(VertexPtr());
continue;
}
break;
}
// normal case
if(!branch) {
if(fb.kinematics) fb.hard = false;
return fb;
}
// truncated emission
if(fb.kinematics) {
fb.hard = false;
fb.iout = iout;
return fb;
}
// otherwise need to return the hard emission
// construct the kinematics for the hard emission
ShoKinPtr showerKin = new_ptr(FS_QTildeShowerKinematics1to2(
branch->scale(),
branch->children()[0]->z(),
branch->phi(),
branch->children()[0]->pT(),
branch->sudakov()
));
IdList idlist(3);
idlist[0] = particle->dataPtr();
idlist[1] = branch->children()[0]->branchingParticle()->dataPtr();
idlist[2] = branch->children()[1]->branchingParticle()->dataPtr();
fb = Branching( showerKin, idlist, branch->sudakov(),branch->type() );
fb.hard = true;
fb.iout=0;
// return it
return fb;
}
Branching QTildeShowerHandler::selectSpaceLikeDecayBranching(tShowerParticlePtr particle,
const ShowerParticle::EvolutionScales & maxScales,
Energy minmass,ShowerInteraction type,
HardBranchingPtr branch) {
Branching fb;
unsigned int iout=0;
while (true) {
// break if doing truncated shower and no truncated shower needed
if(branch && (!isTruncatedShowerON()||hardOnly())) break;
// select branching
fb=_splittingGenerator->chooseDecayBranching(*particle,maxScales,minmass,
_initialenhance,type);
// return if no radiation
if(!fb.kinematics) break;
// special for truncated shower
if(branch) {
// check haven't evolved too far
if(fb.kinematics->scale() < branch->scale()) {
fb=Branching();
break;
}
// find the truncated line
iout=0;
if(fb.ids[1]->id()!=fb.ids[2]->id()) {
if(fb.ids[1]->id()==particle->id()) iout=1;
else if (fb.ids[2]->id()==particle->id()) iout=2;
}
else if(fb.ids[1]->id()==particle->id()) {
if(fb.kinematics->z()>0.5) iout=1;
else iout=2;
}
// apply the vetos for the truncated shower
// no flavour changing branchings
if(iout==0) {
particle->vetoEmission(fb.type,fb.kinematics->scale());
continue;
}
ShowerInteraction type2 = convertInteraction(fb.type);
double zsplit = iout==1 ? fb.kinematics->z() : 1-fb.kinematics->z();
if(type2==branch->sudakov()->interactionType()) {
if(zsplit < 0.5 || // hardest line veto
fb.kinematics->scale()*zsplit < branch->scale() ) { // angular ordering veto
particle->vetoEmission(fb.type,fb.kinematics->scale());
continue;
}
}
// pt veto
if(fb.kinematics->pT() > progenitor()->maximumpT(type2)) {
particle->vetoEmission(fb.type,fb.kinematics->scale());
continue;
}
}
// if not vetoed break
if(spaceLikeDecayVetoed(fb,particle)) {
// otherwise reset scale and continue
particle->vetoEmission(fb.type,fb.kinematics->scale());
continue;
}
break;
}
// normal case
if(!branch) {
if(fb.kinematics) fb.hard = false;
return fb;
}
// truncated emission
if(fb.kinematics) {
fb.hard = false;
fb.iout = iout;
return fb;
}
// otherwise need to return the hard emission
// construct the kinematics for the hard emission
ShoKinPtr showerKin = new_ptr(Decay_QTildeShowerKinematics1to2(
branch->scale(),
branch->children()[0]->z(),
branch->phi(),
branch->children()[0]->pT(),
branch->sudakov()));
IdList idlist(3);
idlist[0] = particle->dataPtr();
idlist[1] = branch->children()[0]->branchingParticle()->dataPtr();
idlist[2] = branch->children()[1]->branchingParticle()->dataPtr();
// create the branching
fb = Branching( showerKin, idlist, branch->sudakov(),ShowerPartnerType::QCDColourLine );
fb.hard=true;
fb.iout=0;
// return it
return fb;
}
void QTildeShowerHandler::checkFlags() {
string error = "Inconsistent hard emission set-up in QTildeShowerHandler::showerHardProcess(). ";
if ( ( currentTree()->isMCatNLOSEvent() || currentTree()->isMCatNLOHEvent() ) ) {
if (_hardEmission ==2 )
throw Exception() << error
<< "Cannot generate POWHEG matching with MC@NLO shower "
<< "approximation. Add 'set QTildeShowerHandler:HardEmission 0' to input file."
<< Exception::runerror;
if ( canHandleMatchboxTrunc() )
throw Exception() << error
<< "Cannot use truncated qtilde shower with MC@NLO shower "
<< "approximation. Set LHCGenerator:EventHandler"
<< ":CascadeHandler to '/Herwig/Shower/ShowerHandler' or "
<< "'/Herwig/Shower/Dipole/DipoleShowerHandler'."
<< Exception::runerror;
}
else if ( ((currentTree()->isPowhegSEvent() || currentTree()->isPowhegHEvent()) ) &&
_hardEmission != 2){
if ( canHandleMatchboxTrunc())
throw Exception() << error
<< "Unmatched events requested for POWHEG shower "
<< "approximation. Set QTildeShowerHandler:HardEmission to "
<< "'POWHEG'."
<< Exception::runerror;
else if (_hardEmissionWarn) {
_hardEmissionWarn = false;
_hardEmission=2;
throw Exception() << error
<< "Unmatched events requested for POWHEG shower "
<< "approximation. Changing QTildeShowerHandler:HardEmission from "
<< _hardEmission << " to 2"
<< Exception::warning;
}
}
if ( currentTree()->isPowhegSEvent() || currentTree()->isPowhegHEvent()) {
if (currentTree()->showerApproximation()->needsTruncatedShower() &&
!canHandleMatchboxTrunc() )
throw Exception() << error
<< "Current shower handler cannot generate truncated shower. "
<< "Set Generator:EventHandler:CascadeHandler to "
<< "'/Herwig/Shower/PowhegShowerHandler'."
<< Exception::runerror;
}
else if ( currentTree()->truncatedShower() && _missingTruncWarn) {
_missingTruncWarn=false;
throw Exception() << "Warning: POWHEG shower approximation used without "
<< "truncated shower. Set Generator:EventHandler:"
<< "CascadeHandler to '/Herwig/Shower/PowhegShowerHandler' and "
<< "'MEMatching:TruncatedShower Yes'."
<< Exception::warning;
}
// else if ( !dipme && _hardEmissionMode > 1 &&
// firstInteraction())
// throw Exception() << error
// << "POWHEG matching requested for LO events. Include "
// << "'set Factory:ShowerApproximation MEMatching' in input file."
// << Exception::runerror;
}
tPPair QTildeShowerHandler::remakeRemnant(tPPair oldp){
// get the parton extractor
PartonExtractor & pex = *lastExtractor();
// get the new partons
tPPair newp = make_pair(findFirstParton(oldp.first ),
findFirstParton(oldp.second));
// if the same do nothing
if(newp == oldp) return oldp;
// Creates the new remnants and returns the new PartonBinInstances
// ATTENTION Broken here for very strange configuration
PBIPair newbins = pex.newRemnants(oldp, newp, newStep());
newStep()->addIntermediate(newp.first);
newStep()->addIntermediate(newp.second);
// return the new partons
return newp;
}
PPtr QTildeShowerHandler::findFirstParton(tPPtr seed) const{
if(seed->parents().empty()) return seed;
tPPtr parent = seed->parents()[0];
//if no parent there this is a loose end which will
//be connected to the remnant soon.
if(!parent || parent == incomingBeams().first ||
parent == incomingBeams().second ) return seed;
else return findFirstParton(parent);
}
void QTildeShowerHandler::decay(ShowerTreePtr tree, ShowerDecayMap & decay) {
// must be one incoming particle
assert(tree->incomingLines().size()==1);
// apply any transforms
tree->applyTransforms();
// if already decayed return
if(!tree->outgoingLines().empty()) return;
// now we need to replace the particle with a new copy after the shower
// find particle after the shower
map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator
tit = tree->parent()->treelinks().find(tree);
assert(tit!=tree->parent()->treelinks().end());
ShowerParticlePtr newparent=tit->second.second;
PerturbativeProcessPtr newProcess = new_ptr(PerturbativeProcess());
newProcess->incoming().push_back(make_pair(newparent,PerturbativeProcessPtr()));
DecayProcessMap decayMap;
ShowerHandler::decay(newProcess,decayMap);
ShowerTree::constructTrees(tree,decay,newProcess,decayMap);
}
namespace {
ShowerProgenitorPtr
findFinalStateLine(ShowerTreePtr tree, long id, Lorentz5Momentum momentum) {
map<ShowerProgenitorPtr,tShowerParticlePtr>::iterator partner;
Energy2 dmin(1e30*GeV2);
for(map<ShowerProgenitorPtr,tShowerParticlePtr>::iterator
cit =tree->outgoingLines().begin(); cit!=tree->outgoingLines().end(); ++cit) {
if(cit->second->id()!=id) continue;
Energy2 test =
sqr(cit->second->momentum().x()-momentum.x())+
sqr(cit->second->momentum().y()-momentum.y())+
sqr(cit->second->momentum().z()-momentum.z())+
sqr(cit->second->momentum().t()-momentum.t());
if(test<dmin) {
dmin = test;
partner = cit;
}
}
return partner->first;
}
ShowerProgenitorPtr
findInitialStateLine(ShowerTreePtr tree, long id, Lorentz5Momentum momentum) {
map<ShowerProgenitorPtr,ShowerParticlePtr>::iterator partner;
Energy2 dmin(1e30*GeV2);
for(map<ShowerProgenitorPtr,ShowerParticlePtr>::iterator
cit =tree->incomingLines().begin(); cit!=tree->incomingLines().end(); ++cit) {
if(cit->second->id()!=id) continue;
Energy2 test =
sqr(cit->second->momentum().x()-momentum.x())+
sqr(cit->second->momentum().y()-momentum.y())+
sqr(cit->second->momentum().z()-momentum.z())+
sqr(cit->second->momentum().t()-momentum.t());
if(test<dmin) {
dmin = test;
partner = cit;
}
}
return partner->first;
}
void fixSpectatorColours(PPtr newSpect,ShowerProgenitorPtr oldSpect,
ColinePair & cline,ColinePair & aline, bool reconnect) {
cline.first = oldSpect->progenitor()->colourLine();
cline.second = newSpect->colourLine();
aline.first = oldSpect->progenitor()->antiColourLine();
aline.second = newSpect->antiColourLine();
if(!reconnect) return;
if(cline.first) {
cline.first ->removeColoured(oldSpect->copy());
cline.first ->removeColoured(oldSpect->progenitor());
cline.second->removeColoured(newSpect);
cline.first ->addColoured(newSpect);
}
if(aline.first) {
aline.first ->removeAntiColoured(oldSpect->copy());
aline.first ->removeAntiColoured(oldSpect->progenitor());
aline.second->removeAntiColoured(newSpect);
aline.first ->addAntiColoured(newSpect);
}
}
void fixInitialStateEmitter(ShowerTreePtr tree, PPtr newEmit,PPtr emitted, ShowerProgenitorPtr emitter,
ColinePair cline,ColinePair aline,double x) {
// sort out the colours
if(emitted->dataPtr()->iColour()==PDT::Colour8) {
// emitter
if(cline.first && cline.first == emitter->progenitor()->antiColourLine() &&
cline.second !=newEmit->antiColourLine()) {
// sort out not radiating line
ColinePtr col = emitter->progenitor()->colourLine();
if(col) {
col->removeColoured(emitter->copy());
col->removeColoured(emitter->progenitor());
newEmit->colourLine()->removeColoured(newEmit);
col->addColoured(newEmit);
}
}
else if(aline.first && aline.first == emitter->progenitor()->colourLine() &&
aline.second !=newEmit->colourLine()) {
// sort out not radiating line
ColinePtr anti = emitter->progenitor()->antiColourLine();
if(anti) {
anti->removeAntiColoured(emitter->copy());
anti->removeAntiColoured(emitter->progenitor());
newEmit->colourLine()->removeAntiColoured(newEmit);
anti->addAntiColoured(newEmit);
}
}
else
assert(false);
// emitted
if(cline.first && cline.second==emitted->colourLine()) {
cline.second->removeColoured(emitted);
cline.first->addColoured(emitted);
}
else if(aline.first && aline.second==emitted->antiColourLine()) {
aline.second->removeAntiColoured(emitted);
aline.first->addAntiColoured(emitted);
}
else
assert(false);
}
else {
if(emitter->progenitor()->antiColourLine() ) {
ColinePtr col = emitter->progenitor()->antiColourLine();
col->removeAntiColoured(emitter->copy());
col->removeAntiColoured(emitter->progenitor());
if(newEmit->antiColourLine()) {
newEmit->antiColourLine()->removeAntiColoured(newEmit);
col->addAntiColoured(newEmit);
}
else if (emitted->colourLine()) {
emitted->colourLine()->removeColoured(emitted);
col->addColoured(emitted);
}
else
assert(false);
}
if(emitter->progenitor()->colourLine() ) {
ColinePtr col = emitter->progenitor()->colourLine();
col->removeColoured(emitter->copy());
col->removeColoured(emitter->progenitor());
if(newEmit->colourLine()) {
newEmit->colourLine()->removeColoured(newEmit);
col->addColoured(newEmit);
}
else if (emitted->antiColourLine()) {
emitted->antiColourLine()->removeAntiColoured(emitted);
col->addAntiColoured(emitted);
}
else
assert(false);
}
}
// update the emitter
emitter->copy(newEmit);
ShowerParticlePtr sp = new_ptr(ShowerParticle(*newEmit,1,false));
sp->x(x);
emitter->progenitor(sp);
tree->incomingLines()[emitter]=sp;
emitter->perturbative(false);
// add emitted
sp=new_ptr(ShowerParticle(*emitted,1,true));
ShowerProgenitorPtr gluon=new_ptr(ShowerProgenitor(emitter->original(),emitted,sp));
gluon->perturbative(false);
tree->outgoingLines().insert(make_pair(gluon,sp));
}
void fixFinalStateEmitter(ShowerTreePtr tree, PPtr newEmit,PPtr emitted, ShowerProgenitorPtr emitter,
ColinePair cline,ColinePair aline) {
map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator tit;
// special case if decayed
for(tit = tree->treelinks().begin(); tit != tree->treelinks().end();++tit) {
if(tit->second.first && tit->second.second==emitter->progenitor())
break;
}
// sort out the colour lines
if(cline.first && cline.first == emitter->progenitor()->antiColourLine() &&
cline.second !=newEmit->antiColourLine()) {
// sort out not radiating line
ColinePtr col = emitter->progenitor()->colourLine();
if(col) {
col->removeColoured(emitter->copy());
col->removeColoured(emitter->progenitor());
newEmit->colourLine()->removeColoured(newEmit);
col->addColoured(newEmit);
}
}
else if(aline.first && aline.first == emitter->progenitor()->colourLine() &&
aline.second !=newEmit->colourLine()) {
// sort out not radiating line
ColinePtr anti = emitter->progenitor()->antiColourLine();
if(anti) {
anti->removeAntiColoured(emitter->copy());
anti->removeAntiColoured(emitter->progenitor());
newEmit->colourLine()->removeAntiColoured(newEmit);
anti->addAntiColoured(newEmit);
}
}
else
assert(false);
// update the emitter
emitter->copy(newEmit);
ShowerParticlePtr sp = new_ptr(ShowerParticle(*newEmit,1,true));
emitter->progenitor(sp);
tree->outgoingLines()[emitter]=sp;
emitter->perturbative(false);
// update for decaying particles
if(tit!=tree->treelinks().end())
tree->updateLink(tit->first,make_pair(emitter,sp));
// add the emitted particle
// sort out the colour
if(cline.first && cline.second==emitted->antiColourLine()) {
cline.second->removeAntiColoured(emitted);
cline.first->addAntiColoured(emitted);
}
else if(aline.first && aline.second==emitted->colourLine()) {
aline.second->removeColoured(emitted);
aline.first->addColoured(emitted);
}
else
assert(false);
sp=new_ptr(ShowerParticle(*emitted,1,true));
ShowerProgenitorPtr gluon=new_ptr(ShowerProgenitor(emitter->original(),
emitted,sp));
gluon->perturbative(false);
tree->outgoingLines().insert(make_pair(gluon,sp));
}
}
void QTildeShowerHandler::setupMECorrection(RealEmissionProcessPtr real) {
assert(real);
// II emission
if(real->emitter() < real->incoming().size() &&
real->spectator() < real->incoming().size()) {
// recoiling system
for( map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
cjt= currentTree()->outgoingLines().begin();
cjt != currentTree()->outgoingLines().end();++cjt ) {
cjt->first->progenitor()->transform(real->transformation());
cjt->first->copy()->transform(real->transformation());
}
// the the radiating system
ShowerProgenitorPtr emitter,spectator;
unsigned int iemit = real->emitter();
unsigned int ispect = real->spectator();
int ig = int(real->emitted())-int(real->incoming().size());
emitter = findInitialStateLine(currentTree(),
real->bornIncoming()[iemit]->id(),
real->bornIncoming()[iemit]->momentum());
spectator = findInitialStateLine(currentTree(),
real->bornIncoming()[ispect]->id(),
real->bornIncoming()[ispect]->momentum());
// sort out the colours
ColinePair cline,aline;
fixSpectatorColours(real->incoming()[ispect],spectator,cline,aline,true);
// update the spectator
spectator->copy(real->incoming()[ispect]);
ShowerParticlePtr sp(new_ptr(ShowerParticle(*real->incoming()[ispect],1,false)));
sp->x(ispect ==0 ? real->x().first :real->x().second);
spectator->progenitor(sp);
currentTree()->incomingLines()[spectator]=sp;
spectator->perturbative(true);
// now for the emitter
fixInitialStateEmitter(currentTree(),real->incoming()[iemit],real->outgoing()[ig],
emitter,cline,aline,iemit ==0 ? real->x().first :real->x().second);
}
// FF emission
else if(real->emitter() >= real->incoming().size() &&
real->spectator() >= real->incoming().size()) {
assert(real->outgoing()[real->emitted()-real->incoming().size()]->id()==ParticleID::g);
// find the emitter and spectator in the shower tree
ShowerProgenitorPtr emitter,spectator;
int iemit = int(real->emitter())-int(real->incoming().size());
emitter = findFinalStateLine(currentTree(),
real->bornOutgoing()[iemit]->id(),
real->bornOutgoing()[iemit]->momentum());
int ispect = int(real->spectator())-int(real->incoming().size());
spectator = findFinalStateLine(currentTree(),
real->bornOutgoing()[ispect]->id(),
real->bornOutgoing()[ispect]->momentum());
map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator tit;
// first the spectator
// special case if decayed
for(tit = currentTree()->treelinks().begin(); tit != currentTree()->treelinks().end();++tit) {
if(tit->second.first && tit->second.second==spectator->progenitor())
break;
}
// sort out the colours
ColinePair cline,aline;
fixSpectatorColours(real->outgoing()[ispect],spectator,cline,aline,true);
// update the spectator
spectator->copy(real->outgoing()[ispect]);
ShowerParticlePtr sp(new_ptr(ShowerParticle(*real->outgoing()[ispect],1,true)));
spectator->progenitor(sp);
currentTree()->outgoingLines()[spectator]=sp;
spectator->perturbative(true);
// update for decaying particles
if(tit!=currentTree()->treelinks().end())
currentTree()->updateLink(tit->first,make_pair(spectator,sp));
// now the emitting particle
int ig = int(real->emitted())-int(real->incoming().size());
fixFinalStateEmitter(currentTree(),real->outgoing()[iemit],
real->outgoing()[ig],
emitter,cline,aline);
}
// IF emission
else {
// scattering process
if(real->incoming().size()==2) {
ShowerProgenitorPtr emitter,spectator;
unsigned int iemit = real->emitter();
unsigned int ispect = real->spectator();
int ig = int(real->emitted())-int(real->incoming().size());
ColinePair cline,aline;
// incoming spectator
if(ispect<2) {
spectator = findInitialStateLine(currentTree(),
real->bornIncoming()[ispect]->id(),
real->bornIncoming()[ispect]->momentum());
fixSpectatorColours(real->incoming()[ispect],spectator,cline,aline,true);
// update the spectator
spectator->copy(real->incoming()[ispect]);
ShowerParticlePtr sp(new_ptr(ShowerParticle(*real->incoming()[ispect],1,false)));
sp->x(ispect ==0 ? real->x().first :real->x().second);
spectator->progenitor(sp);
currentTree()->incomingLines()[spectator]=sp;
spectator->perturbative(true);
}
// outgoing spectator
else {
spectator = findFinalStateLine(currentTree(),
real->bornOutgoing()[ispect-real->incoming().size()]->id(),
real->bornOutgoing()[ispect-real->incoming().size()]->momentum());
// special case if decayed
map<tShowerTreePtr,pair<tShowerProgenitorPtr,tShowerParticlePtr> >::const_iterator tit;
for(tit = currentTree()->treelinks().begin(); tit != currentTree()->treelinks().end();++tit) {
if(tit->second.first && tit->second.second==spectator->progenitor())
break;
}
fixSpectatorColours(real->outgoing()[ispect-real->incoming().size()],spectator,cline,aline,true);
// update the spectator
spectator->copy(real->outgoing()[ispect-real->incoming().size()]);
ShowerParticlePtr sp(new_ptr(ShowerParticle(*real->outgoing()[ispect-real->incoming().size()],1,true)));
spectator->progenitor(sp);
currentTree()->outgoingLines()[spectator]=sp;
spectator->perturbative(true);
// update for decaying particles
if(tit!=currentTree()->treelinks().end())
currentTree()->updateLink(tit->first,make_pair(spectator,sp));
}
// incoming emitter
if(iemit<2) {
emitter = findInitialStateLine(currentTree(),
real->bornIncoming()[iemit]->id(),
real->bornIncoming()[iemit]->momentum());
fixInitialStateEmitter(currentTree(),real->incoming()[iemit],real->outgoing()[ig],
emitter,aline,cline,iemit ==0 ? real->x().first :real->x().second);
}
// outgoing emitter
else {
emitter = findFinalStateLine(currentTree(),
real->bornOutgoing()[iemit-real->incoming().size()]->id(),
real->bornOutgoing()[iemit-real->incoming().size()]->momentum());
fixFinalStateEmitter(currentTree(),real->outgoing()[iemit-real->incoming().size()],
real->outgoing()[ig],emitter,aline,cline);
}
}
// decay process
else {
assert(real->spectator()==0);
unsigned int iemit = real->emitter()-real->incoming().size();
int ig = int(real->emitted())-int(real->incoming().size());
ColinePair cline,aline;
// incoming spectator
ShowerProgenitorPtr spectator = findInitialStateLine(currentTree(),
real->bornIncoming()[0]->id(),
real->bornIncoming()[0]->momentum());
fixSpectatorColours(real->incoming()[0],spectator,cline,aline,false);
// find the emitter
ShowerProgenitorPtr emitter =
findFinalStateLine(currentTree(),
real->bornOutgoing()[iemit]->id(),
real->bornOutgoing()[iemit]->momentum());
// recoiling system
for( map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
cjt= currentTree()->outgoingLines().begin();
cjt != currentTree()->outgoingLines().end();++cjt ) {
if(cjt->first==emitter) continue;
cjt->first->progenitor()->transform(real->transformation());
cjt->first->copy()->transform(real->transformation());
}
// sort out the emitter
fixFinalStateEmitter(currentTree(),real->outgoing()[iemit],
real->outgoing()[ig],emitter,aline,cline);
}
}
// clean up the shower tree
_currenttree->resetShowerProducts();
}
diff --git a/Shower/QTilde/SplittingFunctions/SplittingFunction.cc b/Shower/QTilde/SplittingFunctions/SplittingFunction.cc
--- a/Shower/QTilde/SplittingFunctions/SplittingFunction.cc
+++ b/Shower/QTilde/SplittingFunctions/SplittingFunction.cc
@@ -1,900 +1,915 @@
// -*- C++ -*-
//
// SplittingFunction.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the SplittingFunction class.
//
#include "SplittingFunction.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Utilities/EnumIO.h"
#include "Herwig/Shower/QTilde/Base/ShowerParticle.h"
#include "ThePEG/Utilities/DescribeClass.h"
using namespace Herwig;
DescribeAbstractClass<SplittingFunction,Interfaced>
describeSplittingFunction ("Herwig::SplittingFunction","");
void SplittingFunction::Init() {
static ClassDocumentation<SplittingFunction> documentation
("The SplittingFunction class is the based class for 1->2 splitting functions"
" in Herwig");
static Switch<SplittingFunction,ColourStructure> interfaceColourStructure
("ColourStructure",
"The colour structure for the splitting function",
&SplittingFunction::_colourStructure, Undefined, false, false);
static SwitchOption interfaceColourStructureTripletTripletOctet
(interfaceColourStructure,
"TripletTripletOctet",
"3 -> 3 8",
TripletTripletOctet);
static SwitchOption interfaceColourStructureOctetOctetOctet
(interfaceColourStructure,
"OctetOctetOctet",
"8 -> 8 8",
OctetOctetOctet);
static SwitchOption interfaceColourStructureOctetTripletTriplet
(interfaceColourStructure,
"OctetTripletTriplet",
"8 -> 3 3bar",
OctetTripletTriplet);
static SwitchOption interfaceColourStructureTripletOctetTriplet
(interfaceColourStructure,
"TripletOctetTriplet",
"3 -> 8 3",
TripletOctetTriplet);
static SwitchOption interfaceColourStructureSextetSextetOctet
(interfaceColourStructure,
"SextetSextetOctet",
"6 -> 6 8",
SextetSextetOctet);
static SwitchOption interfaceColourStructureChargedChargedNeutral
(interfaceColourStructure,
"ChargedChargedNeutral",
"q -> q 0",
ChargedChargedNeutral);
static SwitchOption interfaceColourStructureNeutralChargedCharged
(interfaceColourStructure,
"NeutralChargedCharged",
"0 -> q qbar",
NeutralChargedCharged);
static SwitchOption interfaceColourStructureChargedNeutralCharged
(interfaceColourStructure,
"ChargedNeutralCharged",
"q -> 0 q",
ChargedNeutralCharged);
static Switch<SplittingFunction,ShowerInteraction>
interfaceInteractionType
("InteractionType",
"Type of the interaction",
&SplittingFunction::_interactionType,
ShowerInteraction::UNDEFINED, false, false);
static SwitchOption interfaceInteractionTypeQCD
(interfaceInteractionType,
"QCD","QCD",ShowerInteraction::QCD);
static SwitchOption interfaceInteractionTypeQED
(interfaceInteractionType,
"QED","QED",ShowerInteraction::QED);
static Switch<SplittingFunction,bool> interfaceAngularOrdered
("AngularOrdered",
"Whether or not this interaction is angular ordered, "
"normally only g->q qbar and gamma-> f fbar are the only ones which aren't.",
&SplittingFunction::angularOrdered_, true, false, false);
static SwitchOption interfaceAngularOrderedYes
(interfaceAngularOrdered,
"Yes",
"Interaction is angular ordered",
true);
static SwitchOption interfaceAngularOrderedNo
(interfaceAngularOrdered,
"No",
"Interaction isn't angular ordered",
false);
static Switch<SplittingFunction,unsigned int> interfaceScaleChoice
("ScaleChoice",
"The scale choice to be used",
&SplittingFunction::scaleChoice_, 2, false, false);
static SwitchOption interfaceScaleChoicepT
(interfaceScaleChoice,
"pT",
"pT of the branching",
0);
static SwitchOption interfaceScaleChoiceQ2
(interfaceScaleChoice,
"Q2",
"Q2 of the branching",
1);
static SwitchOption interfaceScaleChoiceFromAngularOrdering
(interfaceScaleChoice,
"FromAngularOrdering",
"If angular order use pT, otherwise Q2",
2);
+ static Switch<SplittingFunction,bool> interfaceStrictAO
+ ("StrictAO",
+ "Whether or not to apply strict angular-ordering,"
+ " i.e. for QED even in QCD emission, and vice versa",
+ &SplittingFunction::strictAO_, true, false, false);
+ static SwitchOption interfaceStrictAOYes
+ (interfaceStrictAO,
+ "Yes",
+ "Apply strict ordering",
+ true);
+ static SwitchOption interfaceStrictAONo
+ (interfaceStrictAO,
+ "No",
+ "Don't apply strict ordering",
+ false);
}
void SplittingFunction::persistentOutput(PersistentOStream & os) const {
os << oenum(_interactionType)
<< oenum(_colourStructure) << _colourFactor
<< angularOrdered_ << scaleChoice_ << strictAO_;
}
void SplittingFunction::persistentInput(PersistentIStream & is, int) {
is >> ienum(_interactionType)
>> ienum(_colourStructure) >> _colourFactor
>> angularOrdered_ >> scaleChoice_ >> strictAO_;
}
void SplittingFunction::colourConnection(tShowerParticlePtr parent,
tShowerParticlePtr first,
tShowerParticlePtr second,
ShowerPartnerType partnerType,
const bool back) const {
if(_colourStructure==TripletTripletOctet) {
if(!back) {
ColinePair cparent = ColinePair(parent->colourLine(),
parent->antiColourLine());
// ensure input consistency
assert(( cparent.first && !cparent.second &&
partnerType==ShowerPartnerType::QCDColourLine) ||
( !cparent.first && cparent.second &&
partnerType==ShowerPartnerType::QCDAntiColourLine));
// q -> q g
if(cparent.first) {
ColinePtr newline=new_ptr(ColourLine());
cparent.first->addColoured(second);
newline->addColoured ( first);
newline->addAntiColoured (second);
}
// qbar -> qbar g
else {
ColinePtr newline=new_ptr(ColourLine());
cparent.second->addAntiColoured(second);
newline->addColoured(second);
newline->addAntiColoured(first);
}
// Set progenitor
first->progenitor(parent->progenitor());
second->progenitor(parent->progenitor());
}
else {
ColinePair cfirst = ColinePair(first->colourLine(),
first->antiColourLine());
// ensure input consistency
assert(( cfirst.first && !cfirst.second &&
partnerType==ShowerPartnerType::QCDColourLine) ||
( !cfirst.first && cfirst.second &&
partnerType==ShowerPartnerType::QCDAntiColourLine));
// q -> q g
if(cfirst.first) {
ColinePtr newline=new_ptr(ColourLine());
cfirst.first->addAntiColoured(second);
newline->addColoured(second);
newline->addColoured(parent);
}
// qbar -> qbar g
else {
ColinePtr newline=new_ptr(ColourLine());
cfirst.second->addColoured(second);
newline->addAntiColoured(second);
newline->addAntiColoured(parent);
}
// Set progenitor
parent->progenitor(first->progenitor());
second->progenitor(first->progenitor());
}
}
else if(_colourStructure==OctetOctetOctet) {
if(!back) {
ColinePair cparent = ColinePair(parent->colourLine(),
parent->antiColourLine());
// ensure input consistency
assert(cparent.first&&cparent.second);
// ensure first gluon is hardest
if( first->id()==second->id() && parent->showerKinematics()->z()<0.5 )
swap(first,second);
// colour line radiates
if(partnerType==ShowerPartnerType::QCDColourLine) {
// The colour line is radiating
ColinePtr newline=new_ptr(ColourLine());
cparent.first->addColoured(second);
cparent.second->addAntiColoured(first);
newline->addColoured(first);
newline->addAntiColoured(second);
}
// anti colour line radiates
else if(partnerType==ShowerPartnerType::QCDAntiColourLine) {
ColinePtr newline=new_ptr(ColourLine());
cparent.first->addColoured(first);
cparent.second->addAntiColoured(second);
newline->addColoured(second);
newline->addAntiColoured(first);
}
else
assert(false);
}
else {
ColinePair cfirst = ColinePair(first->colourLine(),
first->antiColourLine());
// ensure input consistency
assert(cfirst.first&&cfirst.second);
// The colour line is radiating
if(partnerType==ShowerPartnerType::QCDColourLine) {
ColinePtr newline=new_ptr(ColourLine());
cfirst.first->addAntiColoured(second);
cfirst.second->addAntiColoured(parent);
newline->addColoured(parent);
newline->addColoured(second);
}
// anti colour line radiates
else if(partnerType==ShowerPartnerType::QCDAntiColourLine) {
ColinePtr newline=new_ptr(ColourLine());
cfirst.first->addColoured(parent);
cfirst.second->addColoured(second);
newline->addAntiColoured(second);
newline->addAntiColoured(parent);
}
else
assert(false);
}
}
else if(_colourStructure == OctetTripletTriplet) {
if(!back) {
ColinePair cparent = ColinePair(parent->colourLine(),
parent->antiColourLine());
// ensure input consistency
assert(cparent.first&&cparent.second);
cparent.first ->addColoured ( first);
cparent.second->addAntiColoured(second);
// Set progenitor
first->progenitor(parent->progenitor());
second->progenitor(parent->progenitor());
}
else {
ColinePair cfirst = ColinePair(first->colourLine(),
first->antiColourLine());
// ensure input consistency
assert(( cfirst.first && !cfirst.second) ||
(!cfirst.first && cfirst.second));
// g -> q qbar
if(cfirst.first) {
ColinePtr newline=new_ptr(ColourLine());
cfirst.first->addColoured(parent);
newline->addAntiColoured(second);
newline->addAntiColoured(parent);
}
// g -> qbar q
else {
ColinePtr newline=new_ptr(ColourLine());
cfirst.second->addAntiColoured(parent);
newline->addColoured(second);
newline->addColoured(parent);
}
// Set progenitor
parent->progenitor(first->progenitor());
second->progenitor(first->progenitor());
}
}
else if(_colourStructure == TripletOctetTriplet) {
if(!back) {
ColinePair cparent = ColinePair(parent->colourLine(),
parent->antiColourLine());
// ensure input consistency
assert(( cparent.first && !cparent.second) ||
(!cparent.first && cparent.second));
// q -> g q
if(cparent.first) {
ColinePtr newline=new_ptr(ColourLine());
cparent.first->addColoured(first);
newline->addColoured (second);
newline->addAntiColoured( first);
}
// qbar -> g qbar
else {
ColinePtr newline=new_ptr(ColourLine());
cparent.second->addAntiColoured(first);
newline->addColoured ( first);
newline->addAntiColoured(second);
}
// Set progenitor
first->progenitor(parent->progenitor());
second->progenitor(parent->progenitor());
}
else {
ColinePair cfirst = ColinePair(first->colourLine(),
first->antiColourLine());
// ensure input consistency
assert(cfirst.first&&cfirst.second);
// q -> g q
if(parent->id()>0) {
cfirst.first ->addColoured(parent);
cfirst.second->addColoured(second);
}
else {
cfirst.first ->addAntiColoured(second);
cfirst.second->addAntiColoured(parent);
}
// Set progenitor
parent->progenitor(first->progenitor());
second->progenitor(first->progenitor());
}
}
else if(_colourStructure==SextetSextetOctet) {
//make sure we're not doing backward evolution
assert(!back);
//make sure something sensible
assert(parent->colourLine() || parent->antiColourLine());
//get the colour lines or anti-colour lines
bool isAntiColour=true;
ColinePair cparent;
if(parent->colourLine()) {
cparent = ColinePair(const_ptr_cast<tColinePtr>(parent->colourInfo()->colourLines()[0]),
const_ptr_cast<tColinePtr>(parent->colourInfo()->colourLines()[1]));
isAntiColour=false;
}
else {
cparent = ColinePair(const_ptr_cast<tColinePtr>(parent->colourInfo()->antiColourLines()[0]),
const_ptr_cast<tColinePtr>(parent->colourInfo()->antiColourLines()[1]));
}
//check for sensible input
// assert(cparent.first && cparent.second);
// sextet has 2 colour lines
if(!isAntiColour) {
//pick at random which of the colour topolgies to take
double topology = UseRandom::rnd();
if(topology < 0.25) {
ColinePtr newline=new_ptr(ColourLine());
cparent.first->addColoured(second);
cparent.second->addColoured(first);
newline->addColoured(first);
newline->addAntiColoured(second);
}
else if(topology >=0.25 && topology < 0.5) {
ColinePtr newline=new_ptr(ColourLine());
cparent.first->addColoured(first);
cparent.second->addColoured(second);
newline->addColoured(first);
newline->addAntiColoured(second);
}
else if(topology >= 0.5 && topology < 0.75) {
ColinePtr newline=new_ptr(ColourLine());
cparent.first->addColoured(second);
cparent.second->addColoured(first);
newline->addColoured(first);
newline->addAntiColoured(second);
}
else {
ColinePtr newline=new_ptr(ColourLine());
cparent.first->addColoured(first);
cparent.second->addColoured(second);
newline->addColoured(first);
newline->addAntiColoured(second);
}
}
// sextet has 2 anti-colour lines
else {
double topology = UseRandom::rnd();
if(topology < 0.25){
ColinePtr newline=new_ptr(ColourLine());
cparent.first->addAntiColoured(second);
cparent.second->addAntiColoured(first);
newline->addAntiColoured(first);
newline->addColoured(second);
}
else if(topology >=0.25 && topology < 0.5) {
ColinePtr newline=new_ptr(ColourLine());
cparent.first->addAntiColoured(first);
cparent.second->addAntiColoured(second);
newline->addAntiColoured(first);
newline->addColoured(second);
}
else if(topology >= 0.5 && topology < 0.75) {
ColinePtr newline=new_ptr(ColourLine());
cparent.first->addAntiColoured(second);
cparent.second->addAntiColoured(first);
newline->addAntiColoured(first);
newline->addColoured(second);
}
else {
ColinePtr newline=new_ptr(ColourLine());
cparent.first->addAntiColoured(first);
cparent.second->addAntiColoured(second);
newline->addAntiColoured(first);
newline->addColoured(second);
}
}
}
else if(_colourStructure == ChargedChargedNeutral) {
if(!parent->data().coloured()) return;
if(!back) {
ColinePair cparent = ColinePair(parent->colourLine(),
parent->antiColourLine());
// q -> q g
if(cparent.first) {
cparent.first->addColoured(first);
}
// qbar -> qbar g
if(cparent.second) {
cparent.second->addAntiColoured(first);
}
}
else {
ColinePair cfirst = ColinePair(first->colourLine(),
first->antiColourLine());
// q -> q g
if(cfirst.first) {
cfirst.first->addColoured(parent);
}
// qbar -> qbar g
if(cfirst.second) {
cfirst.second->addAntiColoured(parent);
}
}
}
else if(_colourStructure == ChargedNeutralCharged) {
if(!parent->data().coloured()) return;
if(!back) {
ColinePair cparent = ColinePair(parent->colourLine(),
parent->antiColourLine());
// q -> q g
if(cparent.first) {
cparent.first->addColoured(second);
}
// qbar -> qbar g
if(cparent.second) {
cparent.second->addAntiColoured(second);
}
}
else {
if (second->dataPtr()->iColour()==PDT::Colour3 ) {
ColinePtr newline=new_ptr(ColourLine());
newline->addColoured(second);
newline->addColoured(parent);
}
else if (second->dataPtr()->iColour()==PDT::Colour3bar ) {
ColinePtr newline=new_ptr(ColourLine());
newline->addAntiColoured(second);
newline->addAntiColoured(parent);
}
}
}
else if(_colourStructure == NeutralChargedCharged ) {
if(!back) {
if(first->dataPtr()->coloured()) {
ColinePtr newline=new_ptr(ColourLine());
if(first->dataPtr()->iColour()==PDT::Colour3) {
newline->addColoured (first );
newline->addAntiColoured(second);
}
else if (first->dataPtr()->iColour()==PDT::Colour3bar) {
newline->addColoured (second);
newline->addAntiColoured(first );
}
else
assert(false);
}
}
else {
ColinePair cfirst = ColinePair(first->colourLine(),
first->antiColourLine());
// gamma -> q qbar
if(cfirst.first) {
cfirst.first->addAntiColoured(second);
}
// gamma -> qbar q
else if(cfirst.second) {
cfirst.second->addColoured(second);
}
else
assert(false);
}
}
else {
assert(false);
}
}
void SplittingFunction::doinit() {
Interfaced::doinit();
assert(_interactionType!=ShowerInteraction::UNDEFINED);
assert((_colourStructure>0&&_interactionType==ShowerInteraction::QCD) ||
(_colourStructure<0&&_interactionType==ShowerInteraction::QED) );
if(_colourFactor>0.) return;
// compute the colour factors if need
if(_colourStructure==TripletTripletOctet) {
_colourFactor = 4./3.;
}
else if(_colourStructure==OctetOctetOctet) {
_colourFactor = 3.;
}
else if(_colourStructure==OctetTripletTriplet) {
_colourFactor = 0.5;
}
else if(_colourStructure==TripletOctetTriplet) {
_colourFactor = 4./3.;
}
else if(_colourStructure==SextetSextetOctet) {
_colourFactor = 10./3.;
}
else if(_colourStructure<0) {
_colourFactor = 1.;
}
else {
assert(false);
}
}
bool SplittingFunction::checkColours(const IdList & ids) const {
if(_colourStructure==TripletTripletOctet) {
if(ids[0]!=ids[1]) return false;
if((ids[0]->iColour()==PDT::Colour3||ids[0]->iColour()==PDT::Colour3bar) &&
ids[2]->iColour()==PDT::Colour8) return true;
return false;
}
else if(_colourStructure==OctetOctetOctet) {
for(unsigned int ix=0;ix<3;++ix) {
if(ids[ix]->iColour()!=PDT::Colour8) return false;
}
return true;
}
else if(_colourStructure==OctetTripletTriplet) {
if(ids[0]->iColour()!=PDT::Colour8) return false;
if(ids[1]->iColour()==PDT::Colour3&&ids[2]->iColour()==PDT::Colour3bar)
return true;
if(ids[1]->iColour()==PDT::Colour3bar&&ids[2]->iColour()==PDT::Colour3)
return true;
return false;
}
else if(_colourStructure==TripletOctetTriplet) {
if(ids[0]!=ids[2]) return false;
if((ids[0]->iColour()==PDT::Colour3||ids[0]->iColour()==PDT::Colour3bar) &&
ids[1]->iColour()==PDT::Colour8) return true;
return false;
}
else if(_colourStructure==SextetSextetOctet) {
if(ids[0]!=ids[1]) return false;
if((ids[0]->iColour()==PDT::Colour6 || ids[0]->iColour()==PDT::Colour6bar) &&
ids[2]->iColour()==PDT::Colour8) return true;
return false;
}
else if(_colourStructure==ChargedChargedNeutral) {
if(ids[0]!=ids[1]) return false;
if(ids[2]->iCharge()!=0) return false;
if(ids[0]->iCharge()==ids[1]->iCharge()) return true;
return false;
}
else if(_colourStructure==ChargedNeutralCharged) {
if(ids[0]!=ids[2]) return false;
if(ids[1]->iCharge()!=0) return false;
if(ids[0]->iCharge()==ids[2]->iCharge()) return true;
return false;
}
else if(_colourStructure==NeutralChargedCharged) {
if(ids[1]->id()!=-ids[2]->id()) return false;
if(ids[0]->iCharge()!=0) return false;
if(ids[1]->iCharge()==-ids[2]->iCharge()) return true;
return false;
}
else {
assert(false);
}
return false;
}
namespace {
bool hasColour(tPPtr p) {
PDT::Colour colour = p->dataPtr()->iColour();
return colour==PDT::Colour3 || colour==PDT::Colour8 || colour == PDT::Colour6;
}
bool hasAntiColour(tPPtr p) {
PDT::Colour colour = p->dataPtr()->iColour();
return colour==PDT::Colour3bar || colour==PDT::Colour8 || colour == PDT::Colour6bar;
}
}
void SplittingFunction::evaluateFinalStateScales(ShowerPartnerType partnerType,
Energy scale, double z,
tShowerParticlePtr parent,
tShowerParticlePtr emitter,
tShowerParticlePtr emitted) {
// identify emitter and emitted
double zEmitter = z, zEmitted = 1.-z;
bool bosonSplitting(false);
// special for g -> gg, particle highest z is emitter
if(emitter->id() == emitted->id() && emitter->id() == parent->id() &&
zEmitted > zEmitter) {
swap(zEmitted,zEmitter);
swap( emitted, emitter);
}
// otherwise if particle ID same
else if(emitted->id()==parent->id()) {
swap(zEmitted,zEmitter);
swap( emitted, emitter);
}
// no real emitter/emitted
else if(emitter->id()!=parent->id()) {
bosonSplitting = true;
}
// may need to add angularOrder flag here
// now the various scales
// QED
if(partnerType==ShowerPartnerType::QED) {
assert(colourStructure()==ChargedChargedNeutral ||
colourStructure()==ChargedNeutralCharged ||
colourStructure()==NeutralChargedCharged );
// normal case
if(!bosonSplitting) {
assert(colourStructure()==ChargedChargedNeutral ||
colourStructure()==ChargedNeutralCharged );
// set the scales
// emitter
emitter->scales().QED = zEmitter*scale;
emitter->scales().QED_noAO = scale;
if(strictAO_)
emitter->scales().QCD_c = min(zEmitter*scale,parent->scales().QCD_c );
else
emitter->scales().QCD_c = min( scale,parent->scales().QCD_c );
emitter->scales().QCD_c_noAO = min(scale,parent->scales().QCD_c_noAO );
if(strictAO_)
emitter->scales().QCD_ac = min(zEmitter*scale,parent->scales().QCD_ac );
else
emitter->scales().QCD_ac = min( scale,parent->scales().QCD_ac );
emitter->scales().QCD_ac_noAO = min(scale,parent->scales().QCD_ac_noAO);
// emitted
emitted->scales().QED = zEmitted*scale;
emitted->scales().QED_noAO = scale;
emitted->scales().QCD_c = ZERO;
emitted->scales().QCD_c_noAO = ZERO;
emitted->scales().QCD_ac = ZERO;
emitted->scales().QCD_ac_noAO = ZERO;
}
// gamma -> f fbar
else {
assert(colourStructure()==NeutralChargedCharged );
// emitter
emitter->scales().QED = zEmitter*scale;
emitter->scales().QED_noAO = scale;
if(hasColour(emitter)) {
emitter->scales().QCD_c = zEmitter*scale;
emitter->scales().QCD_c_noAO = scale;
}
if(hasAntiColour(emitter)) {
emitter->scales().QCD_ac = zEmitter*scale;
emitter->scales().QCD_ac_noAO = scale;
}
// emitted
emitted->scales().QED = zEmitted*scale;
emitted->scales().QED_noAO = scale;
if(hasColour(emitted)) {
emitted->scales().QCD_c = zEmitted*scale;
emitted->scales().QCD_c_noAO = scale;
}
if(hasAntiColour(emitted)) {
emitted->scales().QCD_ac = zEmitted*scale;
emitted->scales().QCD_ac_noAO = scale;
}
}
}
// QCD
else {
// normal case eg q -> q g and g -> g g
if(!bosonSplitting) {
if(strictAO_)
emitter->scales().QED = min(zEmitter*scale,parent->scales().QED );
else
emitter->scales().QED = min( scale,parent->scales().QED );
emitter->scales().QED_noAO = min(scale,parent->scales().QED_noAO);
if(partnerType==ShowerPartnerType::QCDColourLine) {
emitter->scales().QCD_c = zEmitter*scale;
emitter->scales().QCD_c_noAO = scale;
emitter->scales().QCD_ac = min(zEmitter*scale,parent->scales().QCD_ac );
emitter->scales().QCD_ac_noAO = min( scale,parent->scales().QCD_ac_noAO);
}
else {
emitter->scales().QCD_c = min(zEmitter*scale,parent->scales().QCD_c );
emitter->scales().QCD_c_noAO = min( scale,parent->scales().QCD_c_noAO );
emitter->scales().QCD_ac = zEmitter*scale;
emitter->scales().QCD_ac_noAO = scale;
}
// emitted
emitted->scales().QED = ZERO;
emitted->scales().QED_noAO = ZERO;
emitted->scales().QCD_c = zEmitted*scale;
emitted->scales().QCD_c_noAO = scale;
emitted->scales().QCD_ac = zEmitted*scale;
emitted->scales().QCD_ac_noAO = scale;
}
// g -> q qbar
else {
// emitter
if(emitter->dataPtr()->charged()) {
emitter->scales().QED = zEmitter*scale;
emitter->scales().QED_noAO = scale;
}
emitter->scales().QCD_c = zEmitter*scale;
emitter->scales().QCD_c_noAO = scale;
emitter->scales().QCD_ac = zEmitter*scale;
emitter->scales().QCD_ac_noAO = scale;
// emitted
if(emitted->dataPtr()->charged()) {
emitted->scales().QED = zEmitted*scale;
emitted->scales().QED_noAO = scale;
}
emitted->scales().QCD_c = zEmitted*scale;
emitted->scales().QCD_c_noAO = scale;
emitted->scales().QCD_ac = zEmitted*scale;
emitted->scales().QCD_ac_noAO = scale;
}
}
}
void SplittingFunction::evaluateInitialStateScales(ShowerPartnerType partnerType,
Energy scale, double z,
tShowerParticlePtr parent,
tShowerParticlePtr spacelike,
tShowerParticlePtr timelike) {
// scale for time-like child
Energy AOScale = (1.-z)*scale;
// QED
if(partnerType==ShowerPartnerType::QED) {
if(parent->id()==spacelike->id()) {
// parent
parent ->scales().QED = scale;
parent ->scales().QED_noAO = scale;
parent ->scales().QCD_c = min(scale,spacelike->scales().QCD_c );
parent ->scales().QCD_c_noAO = min(scale,spacelike->scales().QCD_c_noAO );
parent ->scales().QCD_ac = min(scale,spacelike->scales().QCD_ac );
parent ->scales().QCD_ac_noAO = min(scale,spacelike->scales().QCD_ac_noAO);
// timelike
timelike->scales().QED = AOScale;
timelike->scales().QED_noAO = scale;
timelike->scales().QCD_c = ZERO;
timelike->scales().QCD_c_noAO = ZERO;
timelike->scales().QCD_ac = ZERO;
timelike->scales().QCD_ac_noAO = ZERO;
}
else if(parent->id()==timelike->id()) {
parent ->scales().QED = scale;
parent ->scales().QED_noAO = scale;
if(hasColour(parent)) {
parent ->scales().QCD_c = scale;
parent ->scales().QCD_c_noAO = scale;
}
if(hasAntiColour(parent)) {
parent ->scales().QCD_ac = scale;
parent ->scales().QCD_ac_noAO = scale;
}
// timelike
timelike->scales().QED = AOScale;
timelike->scales().QED_noAO = scale;
if(hasColour(timelike)) {
timelike->scales().QCD_c = AOScale;
timelike->scales().QCD_c_noAO = scale;
}
if(hasAntiColour(timelike)) {
timelike->scales().QCD_ac = AOScale;
timelike->scales().QCD_ac_noAO = scale;
}
}
else {
parent ->scales().QED = scale;
parent ->scales().QED_noAO = scale;
parent ->scales().QCD_c = ZERO ;
parent ->scales().QCD_c_noAO = ZERO ;
parent ->scales().QCD_ac = ZERO ;
parent ->scales().QCD_ac_noAO = ZERO ;
// timelike
timelike->scales().QED = AOScale;
timelike->scales().QED_noAO = scale;
if(hasColour(timelike)) {
timelike->scales().QCD_c = min(AOScale,spacelike->scales().QCD_ac );
timelike->scales().QCD_c_noAO = min( scale,spacelike->scales().QCD_ac_noAO);
}
if(hasAntiColour(timelike)) {
timelike->scales().QCD_ac = min(AOScale,spacelike->scales().QCD_c );
timelike->scales().QCD_ac_noAO = min( scale,spacelike->scales().QCD_c_noAO );
}
}
}
// QCD
else {
// timelike
if(timelike->dataPtr()->charged()) {
timelike->scales().QED = AOScale;
timelike->scales().QED_noAO = scale;
}
if(hasColour(timelike)) {
timelike->scales().QCD_c = AOScale;
timelike->scales().QCD_c_noAO = scale;
}
if(hasAntiColour(timelike)) {
timelike->scales().QCD_ac = AOScale;
timelike->scales().QCD_ac_noAO = scale;
}
if(parent->id()==spacelike->id()) {
parent ->scales().QED = min(scale,spacelike->scales().QED );
parent ->scales().QED_noAO = min(scale,spacelike->scales().QED_noAO );
parent ->scales().QCD_c = min(scale,spacelike->scales().QCD_c );
parent ->scales().QCD_c_noAO = min(scale,spacelike->scales().QCD_c_noAO );
parent ->scales().QCD_ac = min(scale,spacelike->scales().QCD_ac );
parent ->scales().QCD_ac_noAO = min(scale,spacelike->scales().QCD_ac_noAO);
}
else {
if(parent->dataPtr()->charged()) {
parent ->scales().QED = scale;
parent ->scales().QED_noAO = scale;
}
if(hasColour(parent)) {
parent ->scales().QCD_c = scale;
parent ->scales().QCD_c_noAO = scale;
}
if(hasAntiColour(parent)) {
parent ->scales().QCD_ac = scale;
parent ->scales().QCD_ac_noAO = scale;
}
}
}
}
void SplittingFunction::evaluateDecayScales(ShowerPartnerType partnerType,
Energy scale, double z,
tShowerParticlePtr parent,
tShowerParticlePtr spacelike,
tShowerParticlePtr timelike) {
assert(parent->id()==spacelike->id());
// angular-ordered scale for 2nd child
Energy AOScale = (1.-z)*scale;
// QED
if(partnerType==ShowerPartnerType::QED) {
// timelike
timelike->scales().QED = AOScale;
timelike->scales().QED_noAO = scale;
timelike->scales().QCD_c = ZERO;
timelike->scales().QCD_c_noAO = ZERO;
timelike->scales().QCD_ac = ZERO;
timelike->scales().QCD_ac_noAO = ZERO;
// spacelike
spacelike->scales().QED = scale;
spacelike->scales().QED_noAO = scale;
}
// QCD
else {
// timelike
timelike->scales().QED = ZERO;
timelike->scales().QED_noAO = ZERO;
timelike->scales().QCD_c = AOScale;
timelike->scales().QCD_c_noAO = scale;
timelike->scales().QCD_ac = AOScale;
timelike->scales().QCD_ac_noAO = scale;
// spacelike
spacelike->scales().QED = max(scale,parent->scales().QED );
spacelike->scales().QED_noAO = max(scale,parent->scales().QED_noAO );
}
spacelike->scales().QCD_c = max(scale,parent->scales().QCD_c );
spacelike->scales().QCD_c_noAO = max(scale,parent->scales().QCD_c_noAO );
spacelike->scales().QCD_ac = max(scale,parent->scales().QCD_ac );
spacelike->scales().QCD_ac_noAO = max(scale,parent->scales().QCD_ac_noAO);
}
diff --git a/Shower/QTilde/SplittingFunctions/SplittingFunction.h b/Shower/QTilde/SplittingFunctions/SplittingFunction.h
--- a/Shower/QTilde/SplittingFunctions/SplittingFunction.h
+++ b/Shower/QTilde/SplittingFunctions/SplittingFunction.h
@@ -1,382 +1,382 @@
// -*- C++ -*-
//
// SplittingFunction.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_SplittingFunction_H
#define HERWIG_SplittingFunction_H
//
// This is the declaration of the SplittingFunction class.
//
#include "ThePEG/Interface/Interfaced.h"
#include "Herwig/Shower/QTilde/ShowerConfig.h"
#include "ThePEG/EventRecord/RhoDMatrix.h"
#include "Herwig/Decay/DecayMatrixElement.h"
#include "Herwig/Shower/QTilde/Kinematics/ShowerKinematics.fh"
#include "ThePEG/EventRecord/ColourLine.h"
#include "ThePEG/PDT/ParticleData.h"
#include "SplittingFunction.fh"
namespace Herwig {
using namespace ThePEG;
/** \ingroup Shower
* Enum to define the possible types of colour structure which can occur in
* the branching.
*/
enum ColourStructure {Undefined=0,
TripletTripletOctet = 1,OctetOctetOctet =2,
OctetTripletTriplet = 3,TripletOctetTriplet=4,
SextetSextetOctet = 5,
ChargedChargedNeutral=-1,ChargedNeutralCharged=-2,
NeutralChargedCharged=-3};
/** \ingroup Shower
*
* This is an abstract class which defines the common interface
* for all \f$1\to2\f$ splitting functions, for both initial-state
* and final-state radiation.
*
* The SplittingFunction class contains a number of purely virtual members
* which must be implemented in the inheriting classes. The class also stores
* the interaction type of the spltting function.
*
* The inheriting classes need to specific the splitting function
* \f$P(z,2p_j\cdot p_k)\f$, in terms of the energy fraction \f$z\f$ and
* the evolution scale. In order to allow the splitting functions to be used
* with different choices of evolution functions the scale is given by
* \f[2p_j\cdot p_k=(p_j+p_k)^2-m_{jk}^2=Q^2-(p_j+p_k)^2=z(1-z)\tilde{q}^2=
* \frac{p_T^2}{z(1-z)}-m_{jk}^2+\frac{m_j^2}{z}+\frac{m_k^2}{1-z},\f]
* where \f$Q^2\f$ is the virtuality of the branching particle,
* $p_T$ is the relative transverse momentum of the branching products and
* \f$\tilde{q}^2\f$ is the angular variable described in hep-ph/0310083.
*
* In addition an overestimate of the
* splitting function, \f$P_{\rm over}(z)\f$ which only depends upon \f$z\f$,
* the integral, inverse of the integral for this overestimate and
* ratio of the true splitting function to the overestimate must be provided
* as they are necessary for the veto alogrithm used to implement the evolution.
*
* @see \ref SplittingFunctionInterfaces "The interfaces"
* defined for SplittingFunction.
*/
class SplittingFunction: public Interfaced {
public:
/**
* The default constructor.
* @param b All splitting functions must have an interaction order
*/
SplittingFunction()
: Interfaced(), _interactionType(ShowerInteraction::UNDEFINED),
_colourStructure(Undefined), _colourFactor(-1.),
- angularOrdered_(true), scaleChoice_(2), strictAO_(false) {}
+ angularOrdered_(true), scaleChoice_(2), strictAO_(true) {}
public:
/**
* Methods to return the interaction type and order for the splitting function
*/
//@{
/**
* Return the type of the interaction
*/
ShowerInteraction interactionType() const {return _interactionType;}
/**
* Return the colour structure
*/
ColourStructure colourStructure() const {return _colourStructure;}
/**
* Return the colour factor
*/
double colourFactor(const IdList &ids) const {
if(_colourStructure>0)
return _colourFactor;
else if(_colourStructure<0) {
if(_colourStructure==ChargedChargedNeutral ||
_colourStructure==ChargedNeutralCharged) {
return sqr(double(ids[0]->iCharge())/3.);
}
else if(_colourStructure==NeutralChargedCharged) {
double fact = sqr(double(ids[1]->iCharge())/3.);
if(ids[1]->coloured())
fact *= abs(double(ids[1]->iColour()));
return fact;
}
else
assert(false);
}
else
assert(false);
}
//@}
/**
* Purely virtual method which should determine whether this splitting
* function can be used for a given set of particles.
* @param ids The PDG codes for the particles in the splitting.
*/
virtual bool accept(const IdList & ids) const = 0;
/**
* Method to check the colours are correct
*/
virtual bool checkColours(const IdList & ids) const;
/**
* Methods to return the splitting function.
*/
//@{
/**
* Purely virtual method which should return the exact value of the splitting function,
* \f$P\f$ evaluated in terms of the energy fraction, \f$z\f$, and the evolution scale
\f$\tilde{q}^2\f$.
* @param z The energy fraction.
* @param t The scale \f$t=2p_j\cdot p_k\f$.
* @param ids The PDG codes for the particles in the splitting.
* @param mass Whether or not to include the mass dependent terms
*/
virtual double P(const double z, const Energy2 t, const IdList & ids,
const bool mass, const RhoDMatrix & rho) const = 0;
/**
* Purely virtual method which should return
* an overestimate of the splitting function,
* \f$P_{\rm over}\f$ such that the result \f$P_{\rm over}\geq P\f$. This function
* should be simple enough that it does not depend on the evolution scale.
* @param z The energy fraction.
* @param ids The PDG codes for the particles in the splitting.
*/
virtual double overestimateP(const double z, const IdList & ids) const = 0;
/**
* Purely virtual method which should return
* the ratio of the splitting function to the overestimate, i.e.
* \f$P(z,\tilde{q}^2)/P_{\rm over}(z)\f$.
* @param z The energy fraction.
* @param t The scale \f$t=2p_j\cdot p_k\f$.
* @param ids The PDG codes for the particles in the splitting.
* @param mass Whether or not to include the mass dependent terms
*/
virtual double ratioP(const double z, const Energy2 t, const IdList & ids,
const bool mass, const RhoDMatrix & rho) const = 0;
/**
* Purely virtual method which should return the indefinite integral of the
* overestimated splitting function, \f$P_{\rm over}\f$.
* @param z The energy fraction.
* @param ids The PDG codes for the particles in the splitting.
* @param PDFfactor Which additional factor to include for the PDF
* 0 is no additional factor,
* 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$
*
*/
virtual double integOverP(const double z, const IdList & ids,
unsigned int PDFfactor=0) const = 0;
/**
* Purely virtual method which should return the inverse of the
* indefinite integral of the
* overestimated splitting function, \f$P_{\rm over}\f$ which is used to
* generate the value of \f$z\f$.
* @param r Value of the splitting function to be inverted
* @param ids The PDG codes for the particles in the splitting.
* @param PDFfactor Which additional factor to include for the PDF
* 0 is no additional factor,
* 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$
*/
virtual double invIntegOverP(const double r, const IdList & ids,
unsigned int PDFfactor=0) const = 0;
//@}
/**
* Purely virtual method which should make the proper colour connection
* between the emitting parent and the branching products.
* @param parent The parent for the branching
* @param first The first branching product
* @param second The second branching product
* @param partnerType The type of evolution partner
* @param back Whether this is foward or backward evolution.
*/
virtual void colourConnection(tShowerParticlePtr parent,
tShowerParticlePtr first,
tShowerParticlePtr second,
ShowerPartnerType partnerType,
const bool back) const;
/**
* Method to calculate the azimuthal angle for forward evolution
* @param z The energy fraction
* @param t The scale \f$t=2p_j\cdot p_k\f$.
* @param ids The PDG codes for the particles in the splitting.
* @param The azimuthal angle, \f$\phi\f$.
* @return The weight
*/
virtual vector<pair<int,Complex> >
generatePhiForward(const double z, const Energy2 t, const IdList & ids,
const RhoDMatrix &) = 0;
/**
* Method to calculate the azimuthal angle for backward evolution
* @param z The energy fraction
* @param t The scale \f$t=2p_j\cdot p_k\f$.
* @param ids The PDG codes for the particles in the splitting.
* @return The weight
*/
virtual vector<pair<int,Complex> >
generatePhiBackward(const double z, const Energy2 t, const IdList & ids,
const RhoDMatrix &) = 0;
/**
* Calculate the matrix element for the splitting
* @param z The energy fraction
* @param t The scale \f$t=2p_j\cdot p_k\f$.
* @param ids The PDG codes for the particles in the splitting.
* @param phi The azimuthal angle, \f$\phi\f$.
* @param timeLike Whether timelike or spacelike, affects inclusive of mass terms
*/
virtual DecayMEPtr matrixElement(const double z, const Energy2 t,
const IdList & ids, const double phi,
bool timeLike) = 0;
/**
* Whether or not the interaction is angular ordered
*/
bool angularOrdered() const {return angularOrdered_;}
/**
* Scale choice
*/
bool pTScale() const {
return scaleChoice_ == 2 ? angularOrdered_ : scaleChoice_ == 0;
}
/**
* Functions to state scales after branching happens
*/
//@{
/**
* Sort out scales for final-state emission
*/
void evaluateFinalStateScales(ShowerPartnerType type,
Energy scale, double z,
tShowerParticlePtr parent,
tShowerParticlePtr first,
tShowerParticlePtr second);
/**
* Sort out scales for initial-state emission
*/
void evaluateInitialStateScales(ShowerPartnerType type,
Energy scale, double z,
tShowerParticlePtr parent,
tShowerParticlePtr first,
tShowerParticlePtr second);
/**
* Sort out scales for decay emission
*/
void evaluateDecayScales(ShowerPartnerType type,
Energy scale, double z,
tShowerParticlePtr parent,
tShowerParticlePtr first,
tShowerParticlePtr second);
//@}
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
//@}
protected:
/**
* Set the colour factor
*/
void colourFactor(double in) {_colourFactor=in;}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
SplittingFunction & operator=(const SplittingFunction &) = delete;
private:
/**
* The interaction type for the splitting function.
*/
ShowerInteraction _interactionType;
/**
* The colour structure
*/
ColourStructure _colourStructure;
/**
* The colour factor
*/
double _colourFactor;
/**
* Whether or not this interaction is angular-ordered
*/
bool angularOrdered_;
/**
* The choice of scale
*/
unsigned int scaleChoice_;
/**
* Enforce strict AO
*/
bool strictAO_;
};
}
#endif /* HERWIG_SplittingFunction_H */
diff --git a/src/defaults/Shower.in b/src/defaults/Shower.in
--- a/src/defaults/Shower.in
+++ b/src/defaults/Shower.in
@@ -1,312 +1,320 @@
# -*- ThePEG-repository -*-
############################################################
# Setup of default parton shower
#
# Useful switches for users are marked near the top of
# this file.
#
# Don't edit this file directly, but reset the switches
# in your own input files!
############################################################
library HwMPI.so
library HwShower.so
library HwMatching.so
mkdir /Herwig/Shower
cd /Herwig/Shower
create Herwig::QTildeShowerHandler ShowerHandler
newdef ShowerHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler
newdef ShowerHandler:RemDecayer /Herwig/Partons/RemnantDecayer
# use LO PDFs for Shower, can be changed later
newdef ShowerHandler:PDFA /Herwig/Partons/ShowerLOPDF
newdef ShowerHandler:PDFB /Herwig/Partons/ShowerLOPDF
newdef ShowerHandler:PDFARemnant /Herwig/Partons/RemnantPDF
newdef ShowerHandler:PDFBRemnant /Herwig/Partons/RemnantPDF
#####################################
# initial setup, don't change these!
#####################################
create Herwig::SplittingGenerator SplittingGenerator
create Herwig::ShowerAlphaQCD AlphaQCD
create Herwig::ShowerAlphaQED AlphaQED
create Herwig::PartnerFinder PartnerFinder
newdef PartnerFinder:PartnerMethod 1
newdef PartnerFinder:ScaleChoice 1
create Herwig::KinematicsReconstructor KinematicsReconstructor
newdef KinematicsReconstructor:ReconstructionOption Colour3
newdef KinematicsReconstructor:InitialStateReconOption SofterFraction
newdef KinematicsReconstructor:InitialInitialBoostOption LongTransBoost
newdef /Herwig/Partons/RemnantDecayer:AlphaS AlphaQCD
newdef /Herwig/Partons/RemnantDecayer:AlphaEM AlphaQED
newdef ShowerHandler:PartnerFinder PartnerFinder
newdef ShowerHandler:KinematicsReconstructor KinematicsReconstructor
newdef ShowerHandler:SplittingGenerator SplittingGenerator
newdef ShowerHandler:SpinCorrelations Yes
newdef ShowerHandler:SoftCorrelations Singular
##################################################################
# Intrinsic pT
#
# Recommended:
# 1.9 GeV for Tevatron W/Z production.
# 2.1 GeV for LHC W/Z production at 10 TeV
# 2.2 GeV for LHC W/Z production at 14 TeV
#
# Set all parameters to 0 to disable
##################################################################
newdef ShowerHandler:IntrinsicPtGaussian 1.3*GeV
newdef ShowerHandler:IntrinsicPtBeta 0
newdef ShowerHandler:IntrinsicPtGamma 0*GeV
newdef ShowerHandler:IntrinsicPtIptmax 0*GeV
#############################################################
# Set up truncated shower handler.
#############################################################
create Herwig::PowhegShowerHandler PowhegShowerHandler
set PowhegShowerHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler
set PowhegShowerHandler:RemDecayer /Herwig/Partons/RemnantDecayer
newdef PowhegShowerHandler:PDFA /Herwig/Partons/ShowerLOPDF
newdef PowhegShowerHandler:PDFB /Herwig/Partons/ShowerLOPDF
newdef PowhegShowerHandler:PDFARemnant /Herwig/Partons/RemnantPDF
newdef PowhegShowerHandler:PDFBRemnant /Herwig/Partons/RemnantPDF
newdef PowhegShowerHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler
newdef PowhegShowerHandler:RemDecayer /Herwig/Partons/RemnantDecayer
newdef PowhegShowerHandler:PDFA /Herwig/Partons/ShowerLOPDF
newdef PowhegShowerHandler:PDFB /Herwig/Partons/ShowerLOPDF
newdef PowhegShowerHandler:PDFARemnant /Herwig/Partons/RemnantPDF
newdef PowhegShowerHandler:PDFBRemnant /Herwig/Partons/RemnantPDF
newdef PowhegShowerHandler:PartnerFinder PartnerFinder
newdef PowhegShowerHandler:KinematicsReconstructor KinematicsReconstructor
newdef PowhegShowerHandler:SplittingGenerator SplittingGenerator
newdef PowhegShowerHandler:Interactions QCDandQED
newdef PowhegShowerHandler:SpinCorrelations Yes
newdef PowhegShowerHandler:SoftCorrelations Singular
newdef PowhegShowerHandler:IntrinsicPtGaussian 1.3*GeV
newdef PowhegShowerHandler:IntrinsicPtBeta 0
newdef PowhegShowerHandler:IntrinsicPtGamma 0*GeV
newdef PowhegShowerHandler:IntrinsicPtIptmax 0*GeV
newdef PowhegShowerHandler:EvolutionScheme Q2
#############################################################
# End of interesting user servicable section.
#
# Anything that follows below should only be touched if you
# know what you're doing.
#
# Really.
#############################################################
#
# a few default values
newdef ShowerHandler:MECorrMode 1
newdef ShowerHandler:EvolutionScheme Q2
newdef AlphaQCD:ScaleFactor 1.0
newdef AlphaQCD:NPAlphaS 2
newdef AlphaQCD:Qmin 0.935
newdef AlphaQCD:NumberOfLoops 2
newdef AlphaQCD:AlphaIn 0.126234
#
#
# Lets set up all the splittings
create Herwig::HalfHalfOneSplitFn QtoQGammaSplitFn
set QtoQGammaSplitFn:InteractionType QED
set QtoQGammaSplitFn:ColourStructure ChargedChargedNeutral
set QtoQGammaSplitFn:AngularOrdered Yes
+set QtoQGammaSplitFn:StrictAO No
create Herwig::HalfHalfOneSplitFn QtoQGSplitFn
newdef QtoQGSplitFn:InteractionType QCD
newdef QtoQGSplitFn:ColourStructure TripletTripletOctet
set QtoQGSplitFn:AngularOrdered Yes
+set QtoQGSplitFn:StrictAO No
create Herwig::OneOneOneSplitFn GtoGGSplitFn
newdef GtoGGSplitFn:InteractionType QCD
newdef GtoGGSplitFn:ColourStructure OctetOctetOctet
set GtoGGSplitFn:AngularOrdered Yes
+set GtoGGSplitFn:StrictAO No
create Herwig::OneOneOneMassiveSplitFn WtoWGammaSplitFn
newdef WtoWGammaSplitFn:InteractionType QED
newdef WtoWGammaSplitFn:ColourStructure ChargedChargedNeutral
set WtoWGammaSplitFn:AngularOrdered Yes
+set WtoWGammaSplitFn:StrictAO No
create Herwig::OneHalfHalfSplitFn GtoQQbarSplitFn
newdef GtoQQbarSplitFn:InteractionType QCD
newdef GtoQQbarSplitFn:ColourStructure OctetTripletTriplet
set GtoQQbarSplitFn:AngularOrdered Yes
+set GtoQQbarSplitFn:StrictAO No
create Herwig::OneHalfHalfSplitFn GammatoQQbarSplitFn
newdef GammatoQQbarSplitFn:InteractionType QED
newdef GammatoQQbarSplitFn:ColourStructure NeutralChargedCharged
set GammatoQQbarSplitFn:AngularOrdered Yes
+set GammatoQQbarSplitFn:StrictAO No
create Herwig::HalfOneHalfSplitFn QtoGQSplitFn
newdef QtoGQSplitFn:InteractionType QCD
newdef QtoGQSplitFn:ColourStructure TripletOctetTriplet
set QtoGQSplitFn:AngularOrdered Yes
+set QtoGQSplitFn:StrictAO No
create Herwig::HalfOneHalfSplitFn QtoGammaQSplitFn
newdef QtoGammaQSplitFn:InteractionType QED
newdef QtoGammaQSplitFn:ColourStructure ChargedNeutralCharged
set QtoGammaQSplitFn:AngularOrdered Yes
+set QtoGammaQSplitFn:StrictAO No
#
# Now the Sudakovs
create Herwig::PTCutOff PTCutOff
newdef PTCutOff:pTmin 1.222798*GeV
create Herwig::SudakovFormFactor SudakovCommon
newdef SudakovCommon:Alpha AlphaQCD
newdef SudakovCommon:Cutoff PTCutOff
newdef SudakovCommon:PDFmax 1.0
cp SudakovCommon QtoQGSudakov
newdef QtoQGSudakov:SplittingFunction QtoQGSplitFn
newdef QtoQGSudakov:PDFmax 1.9
cp SudakovCommon QtoQGammaSudakov
set QtoQGammaSudakov:SplittingFunction QtoQGammaSplitFn
set QtoQGammaSudakov:Alpha AlphaQED
set QtoQGammaSudakov:PDFmax 1.9
cp QtoQGammaSudakov LtoLGammaSudakov
cp PTCutOff LtoLGammaPTCutOff
# Technical parameter to stop evolution.
set LtoLGammaPTCutOff:pTmin 0.000001
set LtoLGammaSudakov:Cutoff LtoLGammaPTCutOff
cp SudakovCommon GtoGGSudakov
newdef GtoGGSudakov:SplittingFunction GtoGGSplitFn
newdef GtoGGSudakov:PDFmax 2.0
cp SudakovCommon WtoWGammaSudakov
newdef WtoWGammaSudakov:SplittingFunction WtoWGammaSplitFn
set WtoWGammaSudakov:Alpha AlphaQED
cp SudakovCommon GtoQQbarSudakov
newdef GtoQQbarSudakov:SplittingFunction GtoQQbarSplitFn
newdef GtoQQbarSudakov:PDFmax 120.0
cp SudakovCommon GammatoQQbarSudakov
newdef GammatoQQbarSudakov:SplittingFunction GammatoQQbarSplitFn
set GammatoQQbarSudakov:Alpha AlphaQED
newdef GammatoQQbarSudakov:PDFmax 120.0
cp SudakovCommon GtobbbarSudakov
newdef GtobbbarSudakov:SplittingFunction GtoQQbarSplitFn
newdef GtobbbarSudakov:PDFmax 40000.0
cp SudakovCommon GtoccbarSudakov
newdef GtoccbarSudakov:SplittingFunction GtoQQbarSplitFn
newdef GtoccbarSudakov:PDFmax 2000.0
cp SudakovCommon QtoGQSudakov
newdef QtoGQSudakov:SplittingFunction QtoGQSplitFn
cp SudakovCommon QtoGammaQSudakov
newdef QtoGammaQSudakov:SplittingFunction QtoGammaQSplitFn
set QtoGammaQSudakov:Alpha AlphaQED
cp SudakovCommon utoGuSudakov
newdef utoGuSudakov:SplittingFunction QtoGQSplitFn
newdef utoGuSudakov:PDFFactor OverOneMinusZ
newdef utoGuSudakov:PDFmax 5.0
cp SudakovCommon dtoGdSudakov
newdef dtoGdSudakov:SplittingFunction QtoGQSplitFn
newdef dtoGdSudakov:PDFFactor OverOneMinusZ
#
# Now add the final splittings
#
do SplittingGenerator:AddFinalSplitting u->u,g; QtoQGSudakov
do SplittingGenerator:AddFinalSplitting d->d,g; QtoQGSudakov
do SplittingGenerator:AddFinalSplitting s->s,g; QtoQGSudakov
do SplittingGenerator:AddFinalSplitting c->c,g; QtoQGSudakov
do SplittingGenerator:AddFinalSplitting b->b,g; QtoQGSudakov
do SplittingGenerator:AddFinalSplitting t->t,g; QtoQGSudakov
#
do SplittingGenerator:AddFinalSplitting g->g,g; GtoGGSudakov
#
do SplittingGenerator:AddFinalSplitting g->u,ubar; GtoQQbarSudakov
do SplittingGenerator:AddFinalSplitting g->d,dbar; GtoQQbarSudakov
do SplittingGenerator:AddFinalSplitting g->s,sbar; GtoQQbarSudakov
do SplittingGenerator:AddFinalSplitting g->c,cbar; GtoccbarSudakov
do SplittingGenerator:AddFinalSplitting g->b,bbar; GtobbbarSudakov
do SplittingGenerator:AddFinalSplitting g->t,tbar; GtoQQbarSudakov
#
do SplittingGenerator:AddFinalSplitting gamma->u,ubar; GammatoQQbarSudakov
do SplittingGenerator:AddFinalSplitting gamma->d,dbar; GammatoQQbarSudakov
do SplittingGenerator:AddFinalSplitting gamma->s,sbar; GammatoQQbarSudakov
do SplittingGenerator:AddFinalSplitting gamma->c,cbar; GammatoQQbarSudakov
do SplittingGenerator:AddFinalSplitting gamma->b,bbar; GammatoQQbarSudakov
do SplittingGenerator:AddFinalSplitting gamma->t,tbar; GammatoQQbarSudakov
do SplittingGenerator:AddFinalSplitting gamma->e-,e+; GammatoQQbarSudakov
do SplittingGenerator:AddFinalSplitting gamma->mu-,mu+; GammatoQQbarSudakov
do SplittingGenerator:AddFinalSplitting gamma->tau-,tau+; GammatoQQbarSudakov
#
do SplittingGenerator:AddFinalSplitting u->u,gamma; QtoQGammaSudakov
do SplittingGenerator:AddFinalSplitting d->d,gamma; QtoQGammaSudakov
do SplittingGenerator:AddFinalSplitting s->s,gamma; QtoQGammaSudakov
do SplittingGenerator:AddFinalSplitting c->c,gamma; QtoQGammaSudakov
do SplittingGenerator:AddFinalSplitting b->b,gamma; QtoQGammaSudakov
do SplittingGenerator:AddFinalSplitting t->t,gamma; QtoQGammaSudakov
do SplittingGenerator:AddFinalSplitting e-->e-,gamma; LtoLGammaSudakov
do SplittingGenerator:AddFinalSplitting mu-->mu-,gamma; LtoLGammaSudakov
do SplittingGenerator:AddFinalSplitting tau-->tau-,gamma; LtoLGammaSudakov
do SplittingGenerator:AddFinalSplitting W+->W+,gamma; WtoWGammaSudakov
#
# Now lets add the initial splittings. Remember the form a->b,c; means
# that the current particle b is given and we backward branch to new
# particle a which is initial state and new particle c which is final state
#
do SplittingGenerator:AddInitialSplitting u->u,g; QtoQGSudakov
do SplittingGenerator:AddInitialSplitting d->d,g; QtoQGSudakov
do SplittingGenerator:AddInitialSplitting s->s,g; QtoQGSudakov
do SplittingGenerator:AddInitialSplitting c->c,g; QtoQGSudakov
do SplittingGenerator:AddInitialSplitting b->b,g; QtoQGSudakov
do SplittingGenerator:AddInitialSplitting u->u,gamma; QtoQGammaSudakov
do SplittingGenerator:AddInitialSplitting d->d,gamma; QtoQGammaSudakov
do SplittingGenerator:AddInitialSplitting s->s,gamma; QtoQGammaSudakov
do SplittingGenerator:AddInitialSplitting c->c,gamma; QtoQGammaSudakov
do SplittingGenerator:AddInitialSplitting b->b,gamma; QtoQGammaSudakov
do SplittingGenerator:AddInitialSplitting t->t,gamma; QtoQGammaSudakov
do SplittingGenerator:AddInitialSplitting g->g,g; GtoGGSudakov
#
do SplittingGenerator:AddInitialSplitting g->d,dbar; GtoQQbarSudakov
do SplittingGenerator:AddInitialSplitting g->u,ubar; GtoQQbarSudakov
do SplittingGenerator:AddInitialSplitting g->s,sbar; GtoQQbarSudakov
do SplittingGenerator:AddInitialSplitting g->c,cbar; GtoccbarSudakov
do SplittingGenerator:AddInitialSplitting g->b,bbar; GtobbbarSudakov
#
do SplittingGenerator:AddInitialSplitting gamma->d,dbar; GammatoQQbarSudakov
do SplittingGenerator:AddInitialSplitting gamma->u,ubar; GammatoQQbarSudakov
do SplittingGenerator:AddInitialSplitting gamma->s,sbar; GammatoQQbarSudakov
do SplittingGenerator:AddInitialSplitting gamma->c,cbar; GammatoQQbarSudakov
do SplittingGenerator:AddInitialSplitting gamma->b,bbar; GammatoQQbarSudakov
#
do SplittingGenerator:AddInitialSplitting d->g,d; dtoGdSudakov
do SplittingGenerator:AddInitialSplitting u->g,u; utoGuSudakov
do SplittingGenerator:AddInitialSplitting s->g,s; QtoGQSudakov
do SplittingGenerator:AddInitialSplitting c->g,c; QtoGQSudakov
do SplittingGenerator:AddInitialSplitting b->g,b; QtoGQSudakov
do SplittingGenerator:AddInitialSplitting dbar->g,dbar; dtoGdSudakov
do SplittingGenerator:AddInitialSplitting ubar->g,ubar; utoGuSudakov
do SplittingGenerator:AddInitialSplitting sbar->g,sbar; QtoGQSudakov
do SplittingGenerator:AddInitialSplitting cbar->g,cbar; QtoGQSudakov
do SplittingGenerator:AddInitialSplitting bbar->g,bbar; QtoGQSudakov
#
do SplittingGenerator:AddInitialSplitting d->gamma,d; QtoGammaQSudakov
do SplittingGenerator:AddInitialSplitting u->gamma,u; QtoGammaQSudakov
do SplittingGenerator:AddInitialSplitting s->gamma,s; QtoGammaQSudakov
do SplittingGenerator:AddInitialSplitting c->gamma,c; QtoGammaQSudakov
do SplittingGenerator:AddInitialSplitting b->gamma,b; QtoGammaQSudakov
do SplittingGenerator:AddInitialSplitting dbar->gamma,dbar; QtoGammaQSudakov
do SplittingGenerator:AddInitialSplitting ubar->gamma,ubar; QtoGammaQSudakov
do SplittingGenerator:AddInitialSplitting sbar->gamma,sbar; QtoGammaQSudakov
do SplittingGenerator:AddInitialSplitting cbar->gamma,cbar; QtoGammaQSudakov
do SplittingGenerator:AddInitialSplitting bbar->gamma,bbar; QtoGammaQSudakov

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 7:02 PM (1 d, 12 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805745
Default Alt Text
(182 KB)

Event Timeline