Page MenuHomeHEPForge
Files F7878650

No OneTemporary
Actions

View Options

diff --git a/Shower/Base/Evolver.cc b/Shower/Base/Evolver.cc
deleted file mode 100644
--- a/Shower/Base/Evolver.cc
+++ /dev/null
@@ -1,3280 +0,0 @@
-// -*- C++ -*-
-//
-// Evolver.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
-// Copyright (C) 2002-2011 The Herwig Collaboration
-//
-// Herwig is licenced under version 2 of the GPL, see COPYING for details.
-// Please respect the MCnet academic guidelines, see GUIDELINES for details.
-//
-//
-// This is the implementation of the non-inlined, non-templated member
-// functions of the Evolver class.
-//
-#include "Evolver.h"
-#include "ThePEG/Interface/ClassDocumentation.h"
-#include "ThePEG/Interface/Reference.h"
-#include "ThePEG/Interface/RefVector.h"
-#include "ThePEG/Interface/Switch.h"
-#include "ThePEG/Interface/Parameter.h"
-#include "ThePEG/Persistency/PersistentOStream.h"
-#include "ThePEG/Persistency/PersistentIStream.h"
-#include "Herwig/Shower/Base/ShowerParticle.h"
-#include "ThePEG/Utilities/EnumIO.h"
-#include "ShowerKinematics.h"
-#include "ThePEG/PDT/EnumParticles.h"
-#include "ThePEG/Repository/EventGenerator.h"
-#include "ThePEG/Handlers/EventHandler.h"
-#include "ThePEG/Handlers/StandardEventHandler.h"
-#include "ThePEG/Utilities/Throw.h"
-#include "ShowerTree.h"
-#include "ShowerProgenitor.h"
-#include "KinematicsReconstructor.h"
-#include "PartnerFinder.h"
-#include "ThePEG/Handlers/StandardXComb.h"
-#include "ThePEG/PDT/DecayMode.h"
-#include "Herwig/Shower/ShowerHandler.h"
-#include "ThePEG/Utilities/DescribeClass.h"
-#include "ShowerVertex.h"
-#include "ThePEG/Repository/CurrentGenerator.h"
-#include "Herwig/MatrixElement/Matchbox/Base/SubtractedME.h"
-#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.h"
-#include "ThePEG/Handlers/StandardXComb.h"
-
-using namespace Herwig;
-
-namespace {
-
- /**
- * A struct to order the particles in the same way as in the DecayMode's
- */
- struct ParticleOrdering {
- /**
- * Operator for the ordering
- * @param p1 The first ParticleData object
- * @param p2 The second ParticleData object
- */
- bool operator() (tcPDPtr p1, tcPDPtr p2) {
- return abs(p1->id()) > abs(p2->id()) ||
- ( abs(p1->id()) == abs(p2->id()) && p1->id() > p2->id() ) ||
- ( p1->id() == p2->id() && p1->fullName() > p2->fullName() );
- }
- };
- typedef multiset<tcPDPtr,ParticleOrdering> OrderedParticles;
-
- /**
- * Cached lookup of decay modes.
- * Generator::findDecayMode() is not efficient.
- */
- tDMPtr findDecayMode(const string & tag) {
- static map<string,DMPtr> cache;
- map<string,DMPtr>::const_iterator pos = cache.find(tag);
-
- if ( pos != cache.end() )
- return pos->second;
-
- tDMPtr dm = CurrentGenerator::current().findDecayMode(tag);
- cache[tag] = dm;
- return dm;
- }
-}
-
-DescribeClass<Evolver,Interfaced>
-describeEvolver ("Herwig::Evolver","HwShower.so");
-
-bool Evolver::_hardEmissionModeWarn = true;
-bool Evolver::_missingTruncWarn = true;
-
-IBPtr Evolver::clone() const {
- return new_ptr(*this);
-}
-
-IBPtr Evolver::fullclone() const {
- return new_ptr(*this);
-}
-
-void Evolver::persistentOutput(PersistentOStream & os) const {
- os << _model << _splittingGenerator << _maxtry
- << _meCorrMode << _hardVetoMode << _hardVetoRead << _hardVetoReadOption
- << _limitEmissions << _spinOpt << _softOpt << _hardPOWHEG
- << ounit(_iptrms,GeV) << _beta << ounit(_gamma,GeV) << ounit(_iptmax,GeV)
- << _vetoes << _fullShowerVetoes << _nReWeight << _reWeight
- << _trunc_Mode << _hardEmissionMode << _reconOpt
- << isMCatNLOSEvent << isMCatNLOHEvent
- << isPowhegSEvent << isPowhegHEvent
- << theFactorizationScaleFactor << theRenormalizationScaleFactor << ounit(muPt,GeV)
- << oenum(interaction_) << _maxTryFSR << _maxFailFSR << _fracFSR;
-}
-
-void Evolver::persistentInput(PersistentIStream & is, int) {
- is >> _model >> _splittingGenerator >> _maxtry
- >> _meCorrMode >> _hardVetoMode >> _hardVetoRead >> _hardVetoReadOption
- >> _limitEmissions >> _spinOpt >> _softOpt >> _hardPOWHEG
- >> iunit(_iptrms,GeV) >> _beta >> iunit(_gamma,GeV) >> iunit(_iptmax,GeV)
- >> _vetoes >> _fullShowerVetoes >> _nReWeight >> _reWeight
- >> _trunc_Mode >> _hardEmissionMode >> _reconOpt
- >> isMCatNLOSEvent >> isMCatNLOHEvent
- >> isPowhegSEvent >> isPowhegHEvent
- >> theFactorizationScaleFactor >> theRenormalizationScaleFactor >> iunit(muPt,GeV)
- >> ienum(interaction_) >> _maxTryFSR >> _maxFailFSR >> _fracFSR;
-}
-
-void Evolver::doinit() {
- Interfaced::doinit();
- // calculate max no of FSR vetos
- _maxFailFSR = max(int(_maxFailFSR), int(_fracFSR*double(generator()->N())));
- // 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;
- }
-}
-
-void Evolver::Init() {
-
- static ClassDocumentation<Evolver> documentation
- ("This class is responsible for carrying out the showering,",
- "including the kinematics reconstruction, in a given scale range,"
- "including the option of the POWHEG approach to simulated next-to-leading order"
- " radiation\\cite{Nason:2004rx}.",
- "%\\cite{Nason:2004rx}\n"
- "\\bibitem{Nason:2004rx}\n"
- " P.~Nason,\n"
- " ``A new method for combining NLO QCD with shower Monte Carlo algorithms,''\n"
- " JHEP {\\bf 0411} (2004) 040\n"
- " [arXiv:hep-ph/0409146].\n"
- " %%CITATION = JHEPA,0411,040;%%\n");
-
- static Reference<Evolver,SplittingGenerator>
- interfaceSplitGen("SplittingGenerator",
- "A reference to the SplittingGenerator object",
- &Herwig::Evolver::_splittingGenerator,
- false, false, true, false);
-
- static Reference<Evolver,ShowerModel> interfaceShowerModel
- ("ShowerModel",
- "The pointer to the object which defines the shower evolution model.",
- &Evolver::_model, false, false, true, false, false);
-
- static Parameter<Evolver,unsigned int> interfaceMaxTry
- ("MaxTry",
- "The maximum number of attempts to generate the shower from a"
- " particular ShowerTree",
- &Evolver::_maxtry, 100, 1, 100000,
- false, false, Interface::limited);
-
- static Parameter<Evolver,unsigned int> interfaceNReWeight
- ("NReWeight",
- "The number of attempts for the shower when reweighting",
- &Evolver::_nReWeight, 100, 10, 10000,
- false, false, Interface::limited);
-
- static Switch<Evolver, unsigned int> ifaceMECorrMode
- ("MECorrMode",
- "Choice of the ME Correction Mode",
- &Evolver::_meCorrMode, 1, false, false);
- static SwitchOption off
- (ifaceMECorrMode,"No","MECorrections off", 0);
- static SwitchOption on
- (ifaceMECorrMode,"Yes","hard+soft on", 1);
- static SwitchOption hard
- (ifaceMECorrMode,"Hard","only hard on", 2);
- static SwitchOption soft
- (ifaceMECorrMode,"Soft","only soft on", 3);
-
- static Switch<Evolver, unsigned int> ifaceHardVetoMode
- ("HardVetoMode",
- "Choice of the Hard Veto Mode",
- &Evolver::_hardVetoMode, 1, false, false);
- static SwitchOption HVoff
- (ifaceHardVetoMode,"No","hard vetos off", 0);
- static SwitchOption HVon
- (ifaceHardVetoMode,"Yes","hard vetos on", 1);
- static SwitchOption HVIS
- (ifaceHardVetoMode,"Initial", "only IS emissions vetoed", 2);
- static SwitchOption HVFS
- (ifaceHardVetoMode,"Final","only FS emissions vetoed", 3);
-
- static Switch<Evolver, unsigned int> ifaceHardVetoRead
- ("HardVetoScaleSource",
- "If hard veto scale is to be read",
- &Evolver::_hardVetoRead, 0, false, false);
- static SwitchOption HVRcalc
- (ifaceHardVetoRead,"Calculate","Calculate from hard process", 0);
- static SwitchOption HVRread
- (ifaceHardVetoRead,"Read","Read from XComb->lastScale", 1);
-
- static Switch<Evolver, bool> ifaceHardVetoReadOption
- ("HardVetoReadOption",
- "Apply read-in scale veto to all collisions or just the primary one?",
- &Evolver::_hardVetoReadOption, false, false, false);
- static SwitchOption AllCollisions
- (ifaceHardVetoReadOption,
- "AllCollisions",
- "Read-in pT veto applied to primary and secondary collisions.",
- false);
- static SwitchOption PrimaryCollision
- (ifaceHardVetoReadOption,
- "PrimaryCollision",
- "Read-in pT veto applied to primary but not secondary collisions.",
- true);
-
- static Parameter<Evolver, Energy> ifaceiptrms
- ("IntrinsicPtGaussian",
- "RMS of intrinsic pT of Gaussian distribution:\n"
- "2*(1-Beta)*exp(-sqr(intrinsicpT/RMS))/sqr(RMS)",
- &Evolver::_iptrms, GeV, ZERO, ZERO, 1000000.0*GeV,
- false, false, Interface::limited);
-
- static Parameter<Evolver, double> ifacebeta
- ("IntrinsicPtBeta",
- "Proportion of inverse quadratic distribution in generating intrinsic pT.\n"
- "(1-Beta) is the proportion of Gaussian distribution",
- &Evolver::_beta, 0, 0, 1,
- false, false, Interface::limited);
-
- static Parameter<Evolver, Energy> ifacegamma
- ("IntrinsicPtGamma",
- "Parameter for inverse quadratic:\n"
- "2*Beta*Gamma/(sqr(Gamma)+sqr(intrinsicpT))",
- &Evolver::_gamma,GeV, ZERO, ZERO, 100000.0*GeV,
- false, false, Interface::limited);
-
- static Parameter<Evolver, Energy> ifaceiptmax
- ("IntrinsicPtIptmax",
- "Upper bound on intrinsic pT for inverse quadratic",
- &Evolver::_iptmax,GeV, ZERO, ZERO, 100000.0*GeV,
- false, false, Interface::limited);
-
- static RefVector<Evolver,ShowerVeto> ifaceVetoes
- ("Vetoes",
- "The vetoes to be checked during showering",
- &Evolver::_vetoes, -1,
- false,false,true,true,false);
-
- static RefVector<Evolver,FullShowerVeto> interfaceFullShowerVetoes
- ("FullShowerVetoes",
- "The vetos to be appliede on the full final state of the shower",
- &Evolver::_fullShowerVetoes, -1, false, false, true, false, false);
-
- static Switch<Evolver,unsigned int> interfaceLimitEmissions
- ("LimitEmissions",
- "Limit the number and type of emissions for testing",
- &Evolver::_limitEmissions, 0, false, false);
- static SwitchOption interfaceLimitEmissionsNoLimit
- (interfaceLimitEmissions,
- "NoLimit",
- "Allow an arbitrary number of emissions",
- 0);
- static SwitchOption interfaceLimitEmissionsOneInitialStateEmission
- (interfaceLimitEmissions,
- "OneInitialStateEmission",
- "Allow one emission in the initial state and none in the final state",
- 1);
- static SwitchOption interfaceLimitEmissionsOneFinalStateEmission
- (interfaceLimitEmissions,
- "OneFinalStateEmission",
- "Allow one emission in the final state and none in the initial state",
- 2);
- static SwitchOption interfaceLimitEmissionsHardOnly
- (interfaceLimitEmissions,
- "HardOnly",
- "Only allow radiation from the hard ME correction",
- 3);
- static SwitchOption interfaceLimitEmissionsOneEmission
- (interfaceLimitEmissions,
- "OneEmission",
- "Allow one emission in either the final state or initial state, but not both",
- 4);
-
- static Switch<Evolver,bool> interfaceTruncMode
- ("TruncatedShower", "Include the truncated shower?",
- &Evolver::_trunc_Mode, 1, false, false);
- static SwitchOption interfaceTruncMode0
- (interfaceTruncMode,"No","Truncated Shower is OFF", 0);
- static SwitchOption interfaceTruncMode1
- (interfaceTruncMode,"Yes","Truncated Shower is ON", 1);
-
- static Switch<Evolver,int> interfaceHardEmissionMode
- ("HardEmissionMode",
- "Whether to use ME corrections or POWHEG for the hardest emission",
- &Evolver::_hardEmissionMode, 0, false, false);
- static SwitchOption interfaceHardEmissionModeDecayMECorrection
- (interfaceHardEmissionMode,
- "DecayMECorrection",
- "Old fashioned ME correction for decays only",
- -1);
- static SwitchOption interfaceHardEmissionModeMECorrection
- (interfaceHardEmissionMode,
- "MECorrection",
- "Old fashioned ME correction",
- 0);
- static SwitchOption interfaceHardEmissionModePOWHEG
- (interfaceHardEmissionMode,
- "POWHEG",
- "Powheg style hard emission using internal matrix elements",
- 1);
- static SwitchOption interfaceHardEmissionModeMatchboxPOWHEG
- (interfaceHardEmissionMode,
- "MatchboxPOWHEG",
- "Powheg style emission for the hard process using Matchbox",
- 2);
- static SwitchOption interfaceHardEmissionModeFullPOWHEG
- (interfaceHardEmissionMode,
- "FullPOWHEG",
- "Powheg style emission for the hard process using Matchbox "
- "and decays using internal matrix elements",
- 3);
-
- static Switch<Evolver,ShowerInteraction::Type > interfaceInteraction
- ("Interaction",
- "The interactions to be used in the shower",
- &Evolver::interaction_, ShowerInteraction::QEDQCD, false, false);
- static SwitchOption interfaceInteractionQCD
- (interfaceInteraction,
- "QCDOnly",
- "Only QCD",
- ShowerInteraction::QCD);
- static SwitchOption interfaceInteractionQEDOnly
- (interfaceInteraction,
- "QEDOnly",
- "Only QED",
- ShowerInteraction::QED);
- static SwitchOption interfaceInteractionEWOnly
- (interfaceInteraction,
- "EWOnly",
- "Only EW",
- ShowerInteraction::EW);
- static SwitchOption interfaceInteractionQEDQCD
- (interfaceInteraction,
- "QEDQCD",
- "QED and QCD",
- ShowerInteraction::QEDQCD);
- static SwitchOption interfaceInteractionALL
- (interfaceInteraction,
- "ALL",
- "QED, QCD and EW",
- ShowerInteraction::ALL);
-
- static Switch<Evolver,unsigned int> interfaceReconstructionOption
- ("ReconstructionOption",
- "Treatment of the reconstruction of the transverse momentum of "
- "a branching from the evolution scale.",
- &Evolver::_reconOpt, 0, false, false);
- static SwitchOption interfaceReconstructionOptionCutOff
- (interfaceReconstructionOption,
- "CutOff",
- "Use the cut-off masses in the calculation",
- 0);
- static SwitchOption interfaceReconstructionOptionOffShell
- (interfaceReconstructionOption,
- "OffShell",
- "Use the off-shell masses in the calculation veto the emission of the parent,"
- " no veto in generation of emissions from children",
- 1);
- static SwitchOption interfaceReconstructionOptionOffShell2
- (interfaceReconstructionOption,
- "OffShell2",
- "Use the off-shell masses in the calculation veto the emissions from the children."
- " no veto in generation of emissions from children",
- 2);
- static SwitchOption interfaceReconstructionOptionOffShell3
- (interfaceReconstructionOption,
- "OffShell3",
- "Use the off-shell masses in the calculation veto the emissions from the children."
- " veto in generation of emissions from children using cut-off for second parton",
- 3);
- static SwitchOption interfaceReconstructionOptionOffShell4
- (interfaceReconstructionOption,
- "OffShell4",
- "Ass OffShell3 but with a restriction on the mass of final-state"
- " jets produced via backward evolution.",
- 4);
-
- static Switch<Evolver,unsigned int> interfaceSpinCorrelations
- ("SpinCorrelations",
- "Treatment of spin correlations in the parton shower",
- &Evolver::_spinOpt, 1, false, false);
- static SwitchOption interfaceSpinCorrelationsOff
- (interfaceSpinCorrelations,
- "No",
- "No spin correlations",
- 0);
- static SwitchOption interfaceSpinCorrelationsSpin
- (interfaceSpinCorrelations,
- "Yes",
- "Include the azimuthal spin correlations only",
- 1);
-
- static Switch<Evolver,unsigned int> interfaceSoftCorrelations
- ("SoftCorrelations",
- "Option for the treatment of soft correlations in the parton shower",
- &Evolver::_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<Evolver,bool> interfaceHardPOWHEG
- ("HardPOWHEG",
- "Treatment of powheg emissions which are too hard to have a shower interpretation",
- &Evolver::_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 Parameter<Evolver,unsigned int> interfaceMaxTryFSR
- ("MaxTryFSR",
- "The maximum number of attempted FSR emissions in"
- " the generation of the FSR",
- &Evolver::_maxTryFSR, 100000, 10, 100000000,
- false, false, Interface::limited);
-
- static Parameter<Evolver,unsigned int> interfaceMaxFailFSR
- ("MaxFailFSR",
- "Maximum number of failures generating the FSR",
- &Evolver::_maxFailFSR, 100, 1, 100000000,
- false, false, Interface::limited);
-
-
- static Parameter<Evolver,double> interfaceFSRFailureFraction
- ("FSRFailureFraction",
- "Maximum fraction of events allowed to fail due to too many FSR emissions",
- &Evolver::_fracFSR, 0.001, 1e-10, 1,
- false, false, Interface::limited);
-}
-
-void Evolver::generateIntrinsicpT(vector<ShowerProgenitorPtr> particlesToShower) {
- _intrinsic.clear();
- if ( !ipTon() || !isISRadiationON() ) return;
- // don't do anything for the moment for secondary scatters
- if( !ShowerHandler::currentHandler()->firstInteraction() ) return;
- // generate intrinsic pT
- for(unsigned int ix=0;ix<particlesToShower.size();++ix) {
- // only consider initial-state particles
- if(particlesToShower[ix]->progenitor()->isFinalState()) continue;
- if(!particlesToShower[ix]->progenitor()->dataPtr()->coloured()) continue;
- Energy ipt;
- if(UseRandom::rnd() > _beta) {
- ipt=_iptrms*sqrt(-log(UseRandom::rnd()));
- }
- else {
- ipt=_gamma*sqrt(pow(1.+sqr(_iptmax/_gamma), UseRandom::rnd())-1.);
- }
- pair<Energy,double> pt = make_pair(ipt,UseRandom::rnd(Constants::twopi));
- _intrinsic[particlesToShower[ix]] = pt;
- }
-}
-
-void Evolver::setupMaximumScales(const vector<ShowerProgenitorPtr> & p,
- XCPtr xcomb) {
- // let POWHEG events radiate freely
- if(_hardEmissionMode>0&&hardTree()) {
- vector<ShowerProgenitorPtr>::const_iterator ckt = p.begin();
- for (; ckt != p.end(); ckt++) (*ckt)->maxHardPt(Constants::MaxEnergy);
- return;
- }
- // return if no vetos
- if (!hardVetoOn()) return;
- // find out if hard partonic subprocess.
- bool isPartonic(false);
- map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
- cit = _currenttree->incomingLines().begin();
- Lorentz5Momentum pcm;
- for(; cit!=currentTree()->incomingLines().end(); ++cit) {
- pcm += cit->first->progenitor()->momentum();
- isPartonic |= cit->first->progenitor()->coloured();
- }
- // find 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 (!hardVetoXComb()||
- (hardVetoReadOption()&&
- !ShowerHandler::currentHandler()->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(hardVetoXComb()&&hardVetoReadOption()&&
- !ShowerHandler::currentHandler()->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 *= ShowerHandler::currentHandler()->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 Evolver::setupHardScales(const vector<ShowerProgenitorPtr> & p,
- XCPtr xcomb) {
- if ( hardVetoXComb() &&
- (!hardVetoReadOption() ||
- ShowerHandler::currentHandler()->firstInteraction()) ) {
- Energy hardScale = ZERO;
- if(currentTree()->isHard()) {
- assert(xcomb);
- hardScale = sqrt( xcomb->lastShowerScale() );
- }
- else {
- hardScale = currentTree()->incomingLines().begin()->first
- ->progenitor()->momentum().mass();
- }
- hardScale *= ShowerHandler::currentHandler()->hardScaleFactor();
- vector<ShowerProgenitorPtr>::const_iterator ckt = p.begin();
- for (; ckt != p.end(); ckt++) (*ckt)->hardScale(hardScale);
- muPt = hardScale;
- }
-}
-
-void Evolver::showerHardProcess(ShowerTreePtr hard, XCPtr xcomb) {
-
-
- isMCatNLOSEvent = false;
- isMCatNLOHEvent = false;
- isPowhegSEvent = false;
- isPowhegHEvent = false;
-
- Ptr<SubtractedME>::tptr subme;
- Ptr<MatchboxMEBase>::tptr me;
- Ptr<SubtractionDipole>::tptr dipme;
-
- Ptr<StandardXComb>::ptr sxc = dynamic_ptr_cast<Ptr<StandardXComb>::ptr>(xcomb);
-
- if ( sxc ) {
- subme = dynamic_ptr_cast<Ptr<SubtractedME>::tptr>(sxc->matrixElement());
- me = dynamic_ptr_cast<Ptr<MatchboxMEBase>::tptr>(sxc->matrixElement());
- dipme = dynamic_ptr_cast<Ptr<SubtractionDipole>::tptr>(sxc->matrixElement());
- }
-
- if ( subme ) {
- if ( subme->showerApproximation() ) {
- theShowerApproximation = subme->showerApproximation();
- // separate MCatNLO and POWHEG-type corrections
- if ( !subme->showerApproximation()->needsSplittingGenerator() ) {
- if ( subme->realShowerSubtraction() )
- isMCatNLOHEvent = true;
- else if ( subme->virtualShowerSubtraction() )
- isMCatNLOSEvent = true;
- }
- else {
- if ( subme->realShowerSubtraction() )
- isPowhegHEvent = true;
- else if ( subme->virtualShowerSubtraction() || subme->loopSimSubtraction() )
- isPowhegSEvent = true;
- }
- }
- } else if ( me ) {
- if ( me->factory()->showerApproximation() ) {
- theShowerApproximation = me->factory()->showerApproximation();
- if ( !me->factory()->showerApproximation()->needsSplittingGenerator() )
- isMCatNLOSEvent = true;
- else
- isPowhegSEvent = true;
- }
- }
-
- string error = "Inconsistent hard emission set-up in Evolver::showerHardProcess(). ";
- if ( ( isMCatNLOSEvent || isMCatNLOHEvent ) ){
- if (_hardEmissionMode > 1)
- throw Exception() << error
- << "Cannot generate POWHEG matching with MC@NLO shower "
- << "approximation. Add 'set Evolver:HardEmissionMode 0' to input file."
- << Exception::runerror;
- if ( ShowerHandler::currentHandler()->canHandleMatchboxTrunc())
- throw Exception() << error
- << "Cannot use truncated qtilde shower with MC@NLO shower "
- << "approximation. Set LHCGenerator:EventHandler"
- << ":CascadeHandler to '/Herwig/Shower/ShowerHandler' or "
- << "'/Herwig/DipoleShower/DipoleShowerHandler'."
- << Exception::runerror;
- }
- else if ( ((isPowhegSEvent || isPowhegHEvent) || dipme) &&
- _hardEmissionMode < 2){
- if ( ShowerHandler::currentHandler()->canHandleMatchboxTrunc())
- throw Exception() << error
- << "Unmatched events requested for POWHEG shower "
- << "approximation. Set Evolver:HardEmissionMode to "
- << "'MatchboxPOWHEG' or 'FullPOWHEG'."
- << Exception::runerror;
- else if (_hardEmissionModeWarn){
- _hardEmissionModeWarn = false;
- _hardEmissionMode+=2;
- throw Exception() << error
- << "Unmatched events requested for POWHEG shower "
- << "approximation. Changing Evolver:HardEmissionMode from "
- << _hardEmissionMode-2 << " to " << _hardEmissionMode
- << Exception::warning;
- }
- }
-
- if ( isPowhegSEvent || isPowhegHEvent) {
- if (theShowerApproximation->needsTruncatedShower() &&
- !ShowerHandler::currentHandler()->canHandleMatchboxTrunc() )
- throw Exception() << error
- << "Current shower handler cannot generate truncated shower. "
- << "Set Generator:EventHandler:CascadeHandler to "
- << "'/Herwig/Shower/PowhegShowerHandler'."
- << Exception::runerror;
- }
- else if ( dipme && _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 &&
- ShowerHandler::currentHandler()->firstInteraction())
- throw Exception() << error
- << "POWHEG matching requested for LO events. Include "
- << "'set Factory:ShowerApproximation MEMatching' in input file."
- << Exception::runerror;
-
- _hardme = HwMEBasePtr();
- // extract the matrix element
- tStdXCombPtr lastXC = dynamic_ptr_cast<tStdXCombPtr>(xcomb);
- if(lastXC) {
- _hardme = dynamic_ptr_cast<HwMEBasePtr>(lastXC->matrixElement());
- }
- _decayme = HwDecayerBasePtr();
- // set the current tree
- currentTree(hard);
- hardTree(HardTreePtr());
- // number of attempts if more than one interaction switched on
- unsigned int interactionTry=0;
- do {
- try {
- // generate the showering
- doShowering(true,xcomb);
- // if no vetos return
- return;
- }
- catch (InteractionVeto) {
- currentTree()->clear();
- ++interactionTry;
- }
- }
- while(interactionTry<=5);
- throw Exception() << "Too many tries for shower in "
- << "Evolver::showerHardProcess()"
- << Exception::eventerror;
-}
-
-void Evolver::hardMatrixElementCorrection(bool hard) {
- // set the initial enhancement factors for the soft correction
- _initialenhance = 1.;
- _finalenhance = 1.;
- // if hard matrix element switched off return
- if(!MECOn(hard)) return;
- // see if we can get the correction from the matrix element
- // or decayer
- if(hard) {
- if(_hardme&&_hardme->hasMECorrection()) {
- _hardme->initializeMECorrection(_currenttree,
- _initialenhance,_finalenhance);
- if(hardMEC(hard))
- _hardme->applyHardMatrixElementCorrection(_currenttree);
- }
- }
- else {
- if(_decayme&&_decayme->hasMECorrection()) {
- _decayme->initializeMECorrection(_currenttree,
- _initialenhance,_finalenhance);
- if(hardMEC(hard))
- _decayme->applyHardMatrixElementCorrection(_currenttree);
- }
- }
-}
-
-ShowerParticleVector Evolver::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())
- children[ix]->set5Momentum(Lorentz5Momentum(_progenitor->progenitor()->mass()));
- else
- children[ix]->set5Momentum(Lorentz5Momentum(ids[ix+1]->mass()));
- }
- return children;
-}
-
-bool Evolver::timeLikeShower(tShowerParticlePtr particle,
- ShowerInteraction::Type 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;
- }
- // too many tries
- if(_nFSR>=_maxTryFSR) {
- ++_nFailedFSR;
- // too many failed events
- if(_nFailedFSR>=_maxFailFSR)
- throw Exception() << "Too many events have failed due to too many shower emissions, in\n"
- << "Evolver::timeLikeShower(). Terminating run\n"
- << Exception::runerror;
- throw Exception() << "Too many attempted emissions in Evolver::timeLikeShower()\n"
- << Exception::eventerror;
- }
- // generate the emission
- ShowerParticleVector children;
- int ntry=0;
- // 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];
- bool setupChildren = true;
- while (ntry<50) {
- fc[0] = Branching();
- fc[1] = Branching();
- ++ntry;
- assert(fb.kinematics);
- // has emitted
- // Assign the shower kinematics to the emitting particle.
- if(setupChildren) {
- ++_nFSR;
- 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,_reconOpt>=3);
- // 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;
- }
- setupChildren = false;
- }
- // select branchings for children
- fc[0] = selectTimeLikeBranching(children[0],type,HardBranchingPtr());
- fc[1] = selectTimeLikeBranching(children[1],type,HardBranchingPtr());
- // old default
- if(_reconOpt==0) {
- // 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();
- break;
- }
- // Herwig default
- else if(_reconOpt==1) {
- // 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();
- // branching has happened
- particle->showerKinematics()->updateParent(particle, children,fb.type);
- // clean up the vetoed emission
- if(particle->virtualMass()==ZERO) {
- particle->showerKinematics(ShoKinPtr());
- for(unsigned int ix=0;ix<children.size();++ix)
- particle->abandonChild(children[ix]);
- children.clear();
- if(particle->spinInfo()) particle->spinInfo()->decayVertex(VertexPtr());
- particle->vetoEmission(fb.type,fb.kinematics->scale());
- // generate the new emission
- fb = selectTimeLikeBranching(particle,type,HardBranchingPtr());
- // no emission, return
- if(!fb.kinematics) {
- if(particle->spinInfo()) particle->spinInfo()->develop();
- return false;
- }
- setupChildren = true;
- continue;
- }
- else
- break;
- }
- // veto children
- else if(_reconOpt>=2) {
- // cut-off masses for the branching
- const vector<Energy> & virtualMasses = fb.sudakov->virtualMasses(fb.ids);
- // compute the masses of the children
- Energy masses[3];
- for(unsigned int ix=0;ix<2;++ix) {
- if(fc[ix].kinematics) {
- const vector<Energy> & vm = fc[ix].sudakov->virtualMasses(fc[ix].ids);
- Energy2 q2 =
- fc[ix].kinematics->z()*(1.-fc[ix].kinematics->z())*sqr(fc[ix].kinematics->scale());
- if(fc[ix].ids[0]->id()!=ParticleID::g) q2 += sqr(vm[0]);
- masses[ix+1] = sqrt(q2);
- }
- else {
- masses[ix+1] = virtualMasses[ix+1];
- }
- }
- masses[0] = fb.ids[0]->id()!=ParticleID::g ? virtualMasses[0] : ZERO;
- double z = fb.kinematics->z();
- Energy2 pt2 = z*(1.-z)*(z*(1.-z)*sqr(fb.kinematics->scale()) + sqr(masses[0]))
- - sqr(masses[1])*(1.-z) - sqr(masses[2])*z;
- if(pt2>=ZERO) {
- break;
- }
- else {
- // reset the scales for the children
- for(unsigned int ix=0;ix<2;++ix) {
- if(fc[ix].kinematics)
- children[ix]->vetoEmission(fc[ix].type,fc[ix].kinematics->scale());
- else
- children[ix]->vetoEmission(ShowerPartnerType::QCDColourLine,ZERO);
- children[ix]->virtualMass(ZERO);
- }
- }
- }
- };
- if(_reconOpt>=2) {
- // 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();
- // branching has happened
- particle->showerKinematics()->updateParent(particle, children,fb.type);
- }
- if(first&&!children.empty())
- particle->showerKinematics()->resetChildren(particle,children);
- if(particle->spinInfo()) particle->spinInfo()->develop();
- return true;
-}
-
-bool
-Evolver::spaceLikeShower(tShowerParticlePtr particle, PPtr beam,
- ShowerInteraction::Type type) {
- //using the pdf's associated with the ShowerHandler assures, that
- //modified pdf's are used for the secondary interactions via
- //CascadeHandler::resetPDFs(...)
- tcPDFPtr pdf;
- if(ShowerHandler::currentHandler()->firstPDF().particle() == _beam)
- pdf = ShowerHandler::currentHandler()->firstPDF().pdf();
- if(ShowerHandler::currentHandler()->secondPDF().particle() == _beam)
- pdf = ShowerHandler::currentHandler()->secondPDF().pdf();
- Energy freeze = ShowerHandler::currentHandler()->pdfFreezingScale();
- // don't do anything if not needed
- if(_limitEmissions == 2 || hardOnly() ||
- ( _limitEmissions == 1 && _nis != 0 ) ||
- ( _limitEmissions == 4 && _nis + _nfs != 0 ) ) {
- 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,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,_reconOpt>=4);
- 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 Evolver::showerDecay(ShowerTreePtr decay) {
- _decayme = HwDecayerBasePtr();
- _hardme = HwMEBasePtr();
- // find the decayer
- // try the normal way if possible
- tDMPtr dm = decay->incomingLines().begin()->first->original() ->decayMode();
- if(!dm) dm = decay->incomingLines().begin()->first->copy() ->decayMode();
- if(!dm) dm = decay->incomingLines().begin()->first->progenitor()->decayMode();
- // otherwise make a string and look it up
- if(!dm) {
- string tag = decay->incomingLines().begin()->first->original()->dataPtr()->name()
- + "->";
- 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());
- unsigned int interactionTry=0;
- do {
- try {
- // 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();
- }
- // and then return
- return;
- }
- catch (InteractionVeto) {
- currentTree()->clear();
- ++interactionTry;
- }
- }
- while(interactionTry<=5);
- throw Exception() << "Too many tries for QED shower in Evolver::showerDecay()"
- << Exception::eventerror;
-}
-
-bool Evolver::spaceLikeDecayShower(tShowerParticlePtr particle,
- const ShowerParticle::EvolutionScales & maxScales,
- Energy minmass,ShowerInteraction::Type type,
- Branching fb) {
- // too many tries
- if(_nFSR>=_maxTryFSR) {
- ++_nFailedFSR;
- // too many failed events
- if(_nFailedFSR>=_maxFailFSR)
- throw Exception() << "Too many events have failed due to too many shower emissions, in\n"
- << "Evolver::timeLikeShower(). Terminating run\n"
- << Exception::runerror;
- throw Exception() << "Too many attempted emissions in Evolver::timeLikeShower()\n"
- << Exception::eventerror;
- }
- // generate the emission
- ShowerParticleVector children;
- int ntry=0;
- // generate the emission
- if(!fb.kinematics)
- fb = selectSpaceLikeDecayBranching(particle,maxScales,minmass,type,
- HardBranchingPtr());
- // no emission, return
- if(!fb.kinematics) return false;
- Branching fc[2];
- bool setupChildren = true;
- while (ntry<50) {
- if(particle->virtualMass()==ZERO)
- particle->virtualMass(_progenitor->progenitor()->mass());
- fc[0] = Branching();
- fc[1] = Branching();
- ++ntry;
- assert(fb.kinematics);
- // has emitted
- // Assign the shower kinematics to the emitting particle.
- if(setupChildren) {
- ++_nFSR;
- // 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,_reconOpt>=3);
- setupChildren = false;
- }
- // select branchings for children
- fc[0] = selectSpaceLikeDecayBranching(children[0],maxScales,minmass,
- type,HardBranchingPtr());
- fc[1] = selectTimeLikeBranching (children[1],type,HardBranchingPtr());
- // old default
- if(_reconOpt==0) {
- // 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]);
- // branching has happened
- break;
- }
- // Herwig default
- else if(_reconOpt==1) {
- // 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]);
- // branching has happened
- particle->showerKinematics()->updateParent(particle, children,fb.type);
- // clean up the vetoed emission
- if(particle->virtualMass()==ZERO) {
- particle->showerKinematics(ShoKinPtr());
- for(unsigned int ix=0;ix<children.size();++ix)
- particle->abandonChild(children[ix]);
- children.clear();
- particle->vetoEmission(fb.type,fb.kinematics->scale());
- // generate the new emission
- fb = selectSpaceLikeDecayBranching(particle,maxScales,minmass,type,
- HardBranchingPtr());
- // no emission, return
- if(!fb.kinematics) {
- return false;
- }
- setupChildren = true;
- continue;
- }
- else
- break;
- }
- else if(_reconOpt>=2) {
- // cut-off masses for the branching
- const vector<Energy> & virtualMasses = fb.sudakov->virtualMasses(fb.ids);
- // compute the masses of the children
- Energy masses[3];
- // space-like children
- masses[1] = children[0]->virtualMass();
- // time-like child
- if(fc[1].kinematics) {
- const vector<Energy> & vm = fc[1].sudakov->virtualMasses(fc[1].ids);
- Energy2 q2 =
- fc[1].kinematics->z()*(1.-fc[1].kinematics->z())*sqr(fc[1].kinematics->scale());
- if(fc[1].ids[0]->id()!=ParticleID::g) q2 += sqr(vm[0]);
- masses[2] = sqrt(q2);
- }
- else {
- masses[2] = virtualMasses[2];
- }
- masses[0]=particle->virtualMass();
- double z = fb.kinematics->z();
- Energy2 pt2 = (1.-z)*(z*sqr(masses[0])-sqr(masses[1])-z/(1.-z)*sqr(masses[2]));
- if(pt2>=ZERO) {
- break;
- }
- else {
- // reset the scales for the children
- for(unsigned int ix=0;ix<2;++ix) {
- if(fc[ix].kinematics)
- children[ix]->vetoEmission(fc[ix].type,fc[ix].kinematics->scale());
- else {
- if(ix==0)
- children[ix]->vetoEmission(ShowerPartnerType::QCDColourLine,Constants::MaxEnergy);
- else
- children[ix]->vetoEmission(ShowerPartnerType::QCDColourLine,ZERO);
- }
- }
- children[0]->virtualMass(_progenitor->progenitor()->mass());
- children[1]->virtualMass(ZERO);
- }
- }
- };
- if(_reconOpt>=2) {
- // In the case of splittings which involves coloured particles,
- // set properly the colour flow of the branching.
- // update the history if needed
- _currenttree->updateInitialStateShowerProduct(_progenitor,children[0]);
- _currenttree->addInitialStateBranching(particle,children[0],children[1]);
- // shower the first particle
- 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]);
- // branching has happened
- particle->showerKinematics()->updateParent(particle, children,fb.type);
- }
- // branching has happened
- return true;
-}
-
-vector<ShowerProgenitorPtr> Evolver::setupShower(bool hard) {
- // generate POWHEG hard emission if needed
- if(_hardEmissionMode>0) hardestEmission(hard);
- ShowerInteraction::Type inter = interaction_;
- if(_hardtree&&inter!=ShowerInteraction::QEDQCD && inter!=ShowerInteraction::ALL) {
- inter = _hardtree->interaction();
- }
- // set the initial colour partners
- setEvolutionPartners(hard,inter,false);
- // generate hard me if needed
- if(_hardEmissionMode==0 ||
- (!hard && _hardEmissionMode==-1)) hardMatrixElementCorrection(hard);
- // get the particles to be showered
- vector<ShowerProgenitorPtr> particlesToShower =
- currentTree()->extractProgenitors();
- // remake the colour partners if needed
- if(_currenttree->hardMatrixElementCorrection()) {
- setEvolutionPartners(hard,interaction_,true);
- _currenttree->resetShowerProducts();
- }
- // return the answer
- return particlesToShower;
-}
-
-void Evolver::setEvolutionPartners(bool hard,ShowerInteraction::Type type,
- bool clear) {
- // match the particles in the ShowerTree and hardTree
- if(hardTree() && !hardTree()->connect(currentTree()))
- throw Exception() << "Can't match trees in "
- << "Evolver::setEvolutionPartners()"
- << Exception::eventerror;
- // extract the progenitors
- vector<ShowerParticlePtr> particles =
- currentTree()->extractProgenitorParticles();
- // 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) {
- tHardBranchingPtr partner =
- hardTree()->particles()[particles[ix]]->colourPartner();
- if(!partner) continue;
- for(map<ShowerParticlePtr,tHardBranchingPtr>::const_iterator
- it=hardTree()->particles().begin();
- it!=hardTree()->particles().end();++it) {
- if(it->second==partner) particles[ix]->partner(it->first);
- }
- if(!particles[ix]->partner())
- throw Exception() << "Can't match partners in "
- << "Evolver::setEvolutionPartners()"
- << Exception::eventerror;
- }
- }
- // Set the initial evolution scales
- showerModel()->partnerFinder()->
- setInitialEvolutionScales(particles,!hard,type,!_hardtree);
- 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 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;
- }
- else if(type==ShowerPartnerType::EW) {
- tooHard |= particles[ix]->scales().EW<hardScale;
- }
- }
- }
- if(tooHard) convertHardTree(hard,type);
- }
-}
-
-void Evolver::updateHistory(tShowerParticlePtr particle) {
- if(!particle->children().empty()) {
- ShowerParticleVector theChildren;
- for(unsigned int ix=0;ix<particle->children().size();++ix) {
- ShowerParticlePtr part = dynamic_ptr_cast<ShowerParticlePtr>
- (particle->children()[ix]);
- theChildren.push_back(part);
- }
- // update the history if needed
- if(particle==_currenttree->getFinalStateShowerProduct(_progenitor))
- _currenttree->updateFinalStateShowerProduct(_progenitor,
- particle,theChildren);
- _currenttree->addFinalStateBranching(particle,theChildren);
- for(unsigned int ix=0;ix<theChildren.size();++ix)
- updateHistory(theChildren[ix]);
- }
-}
-
-bool Evolver::startTimeLikeShower(ShowerInteraction::Type type) {
- _nFSR = 0;
- // initialize basis vectors etc
- 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 Evolver::startSpaceLikeShower(PPtr parent, ShowerInteraction::Type type) {
- // initialise the basis vectors
- 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 Evolver::
-startSpaceLikeDecayShower(const ShowerParticle::EvolutionScales & maxScales,
- Energy minimumMass,ShowerInteraction::Type type) {
- _nFSR = 0;
- // set up the particle basis vectors
- 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 Evolver::timeLikeVetoed(const Branching & fb,
- ShowerParticlePtr particle) {
- // work out type of interaction
- ShowerInteraction::Type type = convertInteraction(fb.type);
- // check whether emission was harder than largest pt of hard subprocess
- if ( hardVetoFS() && fb.kinematics->pT() > _progenitor->maxHardPt() )
- return true;
- // soft matrix element correction veto
- if( softMEC()) {
- if(_hardme && _hardme->hasMECorrection()) {
- if(_hardme->softMatrixElementVeto(_progenitor,particle,fb))
- return true;
- }
- else if(_decayme && _decayme->hasMECorrection()) {
- if(_decayme->softMatrixElementVeto(_progenitor,particle,fb))
- return true;
- }
- }
- // veto on maximum pt
- if(fb.kinematics->pT()>_progenitor->maximumpT(type)) return true;
- // general vetos
- if (fb.kinematics && !_vetoes.empty()) {
- bool vetoed=false;
- for (vector<ShowerVetoPtr>::iterator v = _vetoes.begin();
- v != _vetoes.end(); ++v) {
- bool test = (**v).vetoTimeLike(_progenitor,particle,fb);
- switch((**v).vetoType()) {
- case ShowerVeto::Emission:
- vetoed |= test;
- break;
- case ShowerVeto::Shower:
- if(test) throw VetoShower();
- break;
- case ShowerVeto::Event:
- if(test) throw Veto();
- break;
- }
- }
- if(vetoed) return true;
- }
- if ( ShowerHandler::currentHandler()->firstInteraction() &&
- ShowerHandler::currentHandler()->profileScales() ) {
- double weight =
- ShowerHandler::currentHandler()->profileScales()->
- hardScaleProfile(_progenitor->hardScale(),fb.kinematics->pT());
- if ( UseRandom::rnd() > weight )
- return true;
- }
- return false;
-}
-
-bool Evolver::spaceLikeVetoed(const Branching & bb,
- ShowerParticlePtr particle) {
- // work out type of interaction
- ShowerInteraction::Type type = convertInteraction(bb.type);
- // check whether emission was harder than largest pt of hard subprocess
- if (hardVetoIS() && bb.kinematics->pT() > _progenitor->maxHardPt())
- return true;
- // apply the soft correction
- if( softMEC() && _hardme && _hardme->hasMECorrection() ) {
- if(_hardme->softMatrixElementVeto(_progenitor,particle,bb))
- return true;
- }
- // the more general vetos
-
- // check vs max pt for the shower
- if(bb.kinematics->pT()>_progenitor->maximumpT(type)) return true;
-
- if (!_vetoes.empty()) {
- bool vetoed=false;
- for (vector<ShowerVetoPtr>::iterator v = _vetoes.begin();
- v != _vetoes.end(); ++v) {
- bool test = (**v).vetoSpaceLike(_progenitor,particle,bb);
- switch ((**v).vetoType()) {
- case ShowerVeto::Emission:
- vetoed |= test;
- break;
- case ShowerVeto::Shower:
- if(test) throw VetoShower();
- break;
- case ShowerVeto::Event:
- if(test) throw Veto();
- break;
- }
- }
- if (vetoed) return true;
- }
- if ( ShowerHandler::currentHandler()->firstInteraction() &&
- ShowerHandler::currentHandler()->profileScales() ) {
- double weight =
- ShowerHandler::currentHandler()->profileScales()->
- hardScaleProfile(_progenitor->hardScale(),bb.kinematics->pT());
- if ( UseRandom::rnd() > weight )
- return true;
- }
- return false;
-}
-
-bool Evolver::spaceLikeDecayVetoed( const Branching & fb,
- ShowerParticlePtr particle) {
- // work out type of interaction
- ShowerInteraction::Type type = convertInteraction(fb.type);
- // apply the soft correction
- if( softMEC() && _decayme && _decayme->hasMECorrection() ) {
- if(_decayme->softMatrixElementVeto(_progenitor,particle,fb))
- return true;
- }
- // veto on hardest pt in the shower
- if(fb.kinematics->pT()> _progenitor->maximumpT(type)) return true;
- // general vetos
- if (!_vetoes.empty()) {
- bool vetoed=false;
- for (vector<ShowerVetoPtr>::iterator v = _vetoes.begin();
- v != _vetoes.end(); ++v) {
- bool test = (**v).vetoSpaceLike(_progenitor,particle,fb);
- switch((**v).vetoType()) {
- case ShowerVeto::Emission:
- vetoed |= test;
- break;
- case ShowerVeto::Shower:
- if(test) throw VetoShower();
- break;
- case ShowerVeto::Event:
- if(test) throw Veto();
- break;
- }
- if (vetoed) return true;
- }
- }
- return false;
-}
-
-void Evolver::hardestEmission(bool hard) {
- HardTreePtr ISRTree;
- if( ( _hardme && _hardme->hasPOWHEGCorrection()!=0 && _hardEmissionMode< 2) ||
- ( _decayme && _decayme->hasPOWHEGCorrection()!=0 && _hardEmissionMode!=2) ) {
- if(_hardme) {
- assert(hard);
- _hardtree = _hardme->generateHardest( currentTree(),interaction_);
- }
- else {
- assert(!hard);
- _hardtree = _decayme->generateHardest( currentTree() );
- }
- // store initial state POWHEG radiation
- if(_hardtree && _hardme && _hardme->hasPOWHEGCorrection()==1)
- ISRTree=_hardtree;
- }
-
- else if (_hardEmissionMode>1 && hard) {
- // Get minimum pT cutoff used in shower approximation
- Energy maxpt = 1.*GeV;
- 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+=1;
- }
- for( map< ShowerProgenitorPtr, ShowerParticlePtr >::iterator it
- = currentTree()->incomingLines().begin();
- it != currentTree()->incomingLines().end(); ++it ) {
- if( ! it->second->coloured() ) colouredIn+=1;
- }
-
- if ( theShowerApproximation ){
- if ( theShowerApproximation->ffPtCut() == theShowerApproximation->fiPtCut() &&
- theShowerApproximation->ffPtCut() == theShowerApproximation->iiPtCut() )
- maxpt = theShowerApproximation->ffPtCut();
- else if ( colouredIn == 2 && colouredOut == 0 )
- maxpt = theShowerApproximation->iiPtCut();
- else if ( colouredIn == 0 && colouredOut > 1 )
- maxpt = theShowerApproximation->ffPtCut();
- else if ( colouredIn == 2 && colouredOut == 1 )
- maxpt = min(theShowerApproximation->iiPtCut(), theShowerApproximation->fiPtCut());
- else if ( colouredIn == 1 && colouredOut > 1 )
- maxpt = min(theShowerApproximation->ffPtCut(), theShowerApproximation->fiPtCut());
- else
- maxpt = min(min(theShowerApproximation->iiPtCut(), theShowerApproximation->fiPtCut()),
- theShowerApproximation->ffPtCut());
- }
-
- // Generate hardtree from born and real emission subprocesses
- _hardtree = ShowerHandler::currentHandler()->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(ShowerHandler::currentHandler()->lastXCombPtr()->lastShowerScale()),maxpt);
-
- // Set maxpt to pT of emission when showering POWHEG real-emission subprocesses
- if (!isPowhegSEvent && !isPowhegHEvent){
- vector<int> outGluon;
- vector<int> outQuark;
- map< ShowerProgenitorPtr, tShowerParticlePtr >::iterator it;
- for( it = currentTree()->outgoingLines().begin();
- it != currentTree()->outgoingLines().end(); ++it ) {
- if ( abs(it->second->id())< 6) outQuark.push_back(it->second->id());
- if ( it->second->id()==21 ) outGluon.push_back(it->second->id());
- }
- if (outGluon.size() + outQuark.size() == 1){
- for( it = currentTree()->outgoingLines().begin();
- it != currentTree()->outgoingLines().end(); ++it ) {
- if ( abs(it->second->id())< 6 || it->second->id()==21 )
- maxpt = it->second->momentum().perp();
- }
- }
- else if (outGluon.size() + outQuark.size() > 1){
- // assume qqbar pair from a Z/gamma
- if (outGluon.size()==1 && outQuark.size() == 2 && outQuark[0]==-outQuark[1]){
- for( it = currentTree()->outgoingLines().begin();
- it != currentTree()->outgoingLines().end(); ++it ) {
- if ( it->second->id()==21 )
- maxpt = it->second->momentum().perp();
- }
- }
- // otherwise take the lowest pT avoiding born DY events
- else {
- maxpt = generator()->maximumCMEnergy();
- for( it = currentTree()->outgoingLines().begin();
- it != currentTree()->outgoingLines().end(); ++it ) {
- if ( abs(it->second->id())< 6 || it->second->id()==21 )
- maxpt = min(maxpt,it->second->momentum().perp());
- }
- }
- }
- }
-
- // set maximum pT for subsequent emissions from S events
- if ( isPowhegSEvent || (!isPowhegSEvent && !isPowhegHEvent)){
- 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 = ShowerHandler::currentHandler()->generateCKKW(currentTree());
-
- // if hard me doesn't have a FSR powheg
- // correction use decay powheg correction
- if (_hardme && _hardme->hasPOWHEGCorrection()<2) {
- // 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){
- if(_hardtree) connectTrees(currentTree(),_hardtree,hard);
- return;
- }
-
- map<ShowerProgenitorPtr, tShowerParticlePtr > out = currentTree()->outgoingLines();
- // ignore cases where outgoing particles are not coloured
- if (out.size()!=2 ||
- out. begin()->second->dataPtr()->iColour()==PDT::Colour0 ||
- out.rbegin()->second->dataPtr()->iColour()==PDT::Colour0) {
- if(_hardtree) connectTrees(currentTree(),_hardtree,hard);
- 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){
- if(_hardtree) connectTrees(currentTree(),_hardtree,hard);
- return;
- }
-
- // generate the hardest emission
- ShowerDecayMap decay;
- PPtr in = new_ptr(*inter[0]);
- ShowerTreePtr decayTree = new_ptr(ShowerTree(in, decay));
- HardTreePtr FSRTree = decayer->generateHardest(decayTree);
- if (!FSRTree) {
- if(_hardtree) connectTrees(currentTree(),_hardtree,hard);
- 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;
- }
- }
- if(_hardtree){
- connectTrees(currentTree(),_hardtree,hard);
- }
-}
-
-bool Evolver::truncatedTimeLikeShower(tShowerParticlePtr particle,
- HardBranchingPtr branch,
- ShowerInteraction::Type 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;
- int ntry=0;
- Branching fc[2];
- bool setupChildren = true;
- while (ntry<50) {
- if(!fc[0].hard) fc[0] = Branching();
- if(!fc[1].hard) fc[1] = Branching();
- ++ntry;
- // Assign the shower kinematics to the emitting particle.
- if(setupChildren) {
- ++_nFSR;
- // 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,_reconOpt>=3);
- setupChildren = false;
- }
- // 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());
- }
- // old default
- if(_reconOpt==0 || (_reconOpt==1 && fb.hard) ) {
- // 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,fb.type);
- break;
- }
- // H7 default
- else if(_reconOpt==1) {
- // 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);
- 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);
- else
- timeLikeShower(children[1],type,fc[1],false);
- }
- if(children[1]->spinInfo()) children[1]->spinInfo()->develop();
- // branching has happened
- particle->showerKinematics()->updateParent(particle, children,fb.type);
- // clean up the vetoed emission
- if(particle->virtualMass()==ZERO) {
- particle->showerKinematics(ShoKinPtr());
- for(unsigned int ix=0;ix<children.size();++ix)
- particle->abandonChild(children[ix]);
- children.clear();
- if(particle->spinInfo()) particle->spinInfo()->decayVertex(VertexPtr());
- particle->vetoEmission(fb.type,fb.kinematics->scale());
- // generate the new emission
- fb = selectTimeLikeBranching(particle,type,branch);
- // must be at least hard emission
- assert(fb.kinematics);
- setupChildren = true;
- continue;
- }
- else
- break;
- }
- else if(_reconOpt>=2) {
- // cut-off masses for the branching
- const vector<Energy> & virtualMasses = fb.sudakov->virtualMasses(fb.ids);
- // compute the masses of the children
- Energy masses[3];
- for(unsigned int ix=0;ix<2;++ix) {
- if(fc[ix].kinematics) {
- const vector<Energy> & vm = fc[ix].sudakov->virtualMasses(fc[ix].ids);
- Energy2 q2 =
- fc[ix].kinematics->z()*(1.-fc[ix].kinematics->z())*sqr(fc[ix].kinematics->scale());
- if(fc[ix].ids[0]->id()!=ParticleID::g) q2 += sqr(vm[0]);
- masses[ix+1] = sqrt(q2);
- }
- else {
- masses[ix+1] = virtualMasses[ix+1];
- }
- }
- masses[0] = fb.ids[0]->id()!=ParticleID::g ? virtualMasses[0] : ZERO;
- double z = fb.kinematics->z();
- Energy2 pt2 = z*(1.-z)*(z*(1.-z)*sqr(fb.kinematics->scale()) + sqr(masses[0]))
- - sqr(masses[1])*(1.-z) - sqr(masses[2])*z;
- if(pt2>=ZERO) {
- break;
- }
- // if only the hard emission have to accept it
- else if ((fc[0].hard && !fc[1].kinematics) ||
- (fc[1].hard && !fc[0].kinematics) ) {
- break;
- }
- else {
- // reset the scales for the children
- for(unsigned int ix=0;ix<2;++ix) {
- if(fc[ix].hard) continue;
- if(fc[ix].kinematics && ! fc[ix].hard )
- children[ix]->vetoEmission(fc[ix].type,fc[ix].kinematics->scale());
- else
- children[ix]->vetoEmission(ShowerPartnerType::QCDColourLine,ZERO);
- children[ix]->virtualMass(ZERO);
- }
- }
- }
- };
- if(_reconOpt>=2) {
- // 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,fb.type);
- }
- if(first&&!children.empty())
- particle->showerKinematics()->resetChildren(particle,children);
- if(particle->spinInfo()) particle->spinInfo()->develop();
- return true;
-}
-
-bool Evolver::truncatedSpaceLikeShower(tShowerParticlePtr particle, PPtr beam,
- HardBranchingPtr branch,
- ShowerInteraction::Type type) {
- tcPDFPtr pdf;
- if(ShowerHandler::currentHandler()->firstPDF().particle() == beamParticle())
- pdf = ShowerHandler::currentHandler()->firstPDF().pdf();
- if(ShowerHandler::currentHandler()->secondPDF().particle() == beamParticle())
- pdf = ShowerHandler::currentHandler()->secondPDF().pdf();
- Energy freeze = ShowerHandler::currentHandler()->pdfFreezingScale();
- Branching bb;
- // 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::Type 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 =
- branch->sudakov()->createInitialStateBranching( branch->scale(), z, branch->phi(),
- branch->children()[0]->pT() );
- // 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, 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,false);
- 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, 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,false);
- // perform the shower of the final-state particle
- timeLikeShower( otherChild , type,Branching(),true);
- updateHistory(otherChild);
- // return the emitted
- return true;
-}
-
-bool Evolver::
-truncatedSpaceLikeDecayShower(tShowerParticlePtr particle,
- const ShowerParticle::EvolutionScales & maxScales,
- Energy minmass, HardBranchingPtr branch,
- ShowerInteraction::Type 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;
- int ntry=0;
- Branching fc[2];
- bool setupChildren = true;
- while (ntry<50) {
- if(!fc[0].hard) fc[0] = Branching();
- if(!fc[1].hard) fc[1] = Branching();
- ++ntry;
- if(setupChildren) {
- ++_nFSR;
- // 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,_reconOpt>=3);
- setupChildren = false;
- }
- // 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
- if(_reconOpt==0 || (_reconOpt==1 && fb.hard) ) {
- // 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]);
- // branching has happened
- break;
- }
- // H7 default
- else if(_reconOpt==1) {
- // 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);
- }
- }
- // clean up the vetoed emission
- if(particle->virtualMass()==ZERO) {
- particle->showerKinematics(ShoKinPtr());
- for(unsigned int ix=0;ix<children.size();++ix)
- particle->abandonChild(children[ix]);
- children.clear();
- particle->vetoEmission(fb.type,fb.kinematics->scale());
- // generate the new emission
- fb = selectSpaceLikeDecayBranching(particle,maxScales,minmass,type,branch);
- // must be at least hard emission
- assert(fb.kinematics);
- setupChildren = true;
- continue;
- }
- else {
- updateHistory(children[1]);
- break;
- }
- }
- else if(_reconOpt>=2) {
- // cut-off masses for the branching
- const vector<Energy> & virtualMasses = fb.sudakov->virtualMasses(fb.ids);
- // compute the masses of the children
- Energy masses[3];
- // space-like children
- masses[1] = children[0]->virtualMass();
- // time-like child
- if(fc[1].kinematics) {
- const vector<Energy> & vm = fc[1].sudakov->virtualMasses(fc[1].ids);
- Energy2 q2 =
- fc[1].kinematics->z()*(1.-fc[1].kinematics->z())*sqr(fc[1].kinematics->scale());
- if(fc[1].ids[0]->id()!=ParticleID::g) q2 += sqr(vm[0]);
- masses[2] = sqrt(q2);
- }
- else {
- masses[2] = virtualMasses[2];
- }
- masses[0]=particle->virtualMass();
- double z = fb.kinematics->z();
- Energy2 pt2 = (1.-z)*(z*sqr(masses[0])-sqr(masses[1])-z/(1.-z)*sqr(masses[2]));
- if(pt2>=ZERO) {
- break;
- }
- else {
- // reset the scales for the children
- for(unsigned int ix=0;ix<2;++ix) {
- if(fc[ix].kinematics)
- children[ix]->vetoEmission(fc[ix].type,fc[ix].kinematics->scale());
- else {
- if(ix==0)
- children[ix]->vetoEmission(ShowerPartnerType::QCDColourLine,Constants::MaxEnergy);
- else
- children[ix]->vetoEmission(ShowerPartnerType::QCDColourLine,ZERO);
- }
- }
- children[0]->virtualMass(_progenitor->progenitor()->mass());
- children[1]->virtualMass(ZERO);
- }
- }
- };
- if(_reconOpt>=2) {
- // 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]);
- }
- return true;
-}
-
-void Evolver::connectTrees(ShowerTreePtr showerTree,
- HardTreePtr hardTree, bool hard ) {
- ShowerParticleVector particles;
- // find the Sudakovs
- for(set<HardBranchingPtr>::iterator cit=hardTree->branchings().begin();
- cit!=hardTree->branchings().end();++cit) {
- // Sudakovs for ISR
- if((**cit).parent()&&(**cit).status()==HardBranching::Incoming) {
- ++_nis;
- vector<long> br(3);
- br[0] = (**cit).parent()->branchingParticle()->id();
- br[1] = (**cit). branchingParticle()->id();
- br[2] = (**cit).parent()->children()[0]==*cit ?
- (**cit).parent()->children()[1]->branchingParticle()->id() :
- (**cit).parent()->children()[0]->branchingParticle()->id();
- BranchingList branchings = splittingGenerator()->initialStateBranchings();
- if(br[1]<0&&br[0]==br[1]) {
- br[0] = abs(br[0]);
- br[1] = abs(br[1]);
- }
- else if(br[1]<0) {
- br[1] = -br[1];
- br[2] = -br[2];
- }
- long index = abs(br[1]);
- SudakovPtr sudakov;
- for(BranchingList::const_iterator cjt = branchings.lower_bound(index);
- cjt != branchings.upper_bound(index); ++cjt ) {
- IdList ids = cjt->second.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 "
- << "Evolver::connectTrees() for ISR"
- << Exception::runerror;
- (**cit).parent()->sudakov(sudakov);
- }
- // Sudakovs for FSR
- else if(!(**cit).children().empty()) {
- ++_nfs;
- vector<long> br(3);
- br[0] = (**cit) .branchingParticle()->id();
- br[1] = (**cit).children()[0]->branchingParticle()->id();
- br[2] = (**cit).children()[1]->branchingParticle()->id();
- BranchingList branchings = splittingGenerator()->finalStateBranchings();
- if(br[0]<0) {
- br[0] = abs(br[0]);
- br[1] = abs(br[1]);
- br[2] = abs(br[2]);
- }
- long index = br[0];
- SudakovPtr sudakov;
- for(BranchingList::const_iterator cjt = branchings.lower_bound(index);
- cjt != branchings.upper_bound(index); ++cjt ) {
- IdList ids = cjt->second.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 "
- << "Evolver::connectTrees()"
- << Exception::runerror;
- (**cit).sudakov(sudakov);
- }
- }
- // calculate the evolution scale
- for(set<HardBranchingPtr>::iterator cit=hardTree->branchings().begin();
- cit!=hardTree->branchings().end();++cit) {
- particles.push_back((*cit)->branchingParticle());
- }
- showerModel()->partnerFinder()->
- setInitialEvolutionScales(particles,!hard,hardTree->interaction(),
- !hardTree->partnersSet());
- hardTree->partnersSet(true);
- // inverse reconstruction
- if(hard) {
- showerModel()->kinematicsReconstructor()->
- deconstructHardJets(hardTree,ShowerHandler::currentHandler()->evolver(),
- hardTree->interaction());
- }
- else
- showerModel()->kinematicsReconstructor()->
- deconstructDecayJets(hardTree,ShowerHandler::currentHandler()->evolver(),
- hardTree->interaction());
- // now reset the momenta of the showering particles
- vector<ShowerProgenitorPtr> particlesToShower;
- for(map<ShowerProgenitorPtr,ShowerParticlePtr>::const_iterator
- cit=showerTree->incomingLines().begin();
- cit!=showerTree->incomingLines().end();++cit )
- particlesToShower.push_back(cit->first);
- // extract the showering particles
- for(map<ShowerProgenitorPtr,tShowerParticlePtr>::const_iterator
- cit=showerTree->outgoingLines().begin();
- cit!=showerTree->outgoingLines().end();++cit )
- particlesToShower.push_back(cit->first);
- // match them
- 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 Evolver::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 Evolver::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 ( isMCatNLOHEvent && _limitEmissions != 0 ) {
- _nis = _nfs = 1;
- }
- // extract particles to shower
- vector<ShowerProgenitorPtr> particlesToShower(setupShower(hard));
- // setup the maximum scales for the shower
- if (hardVetoOn()) 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())
- minmass += particlesToShower[ix]->progenitor()->dataPtr()->constituentMass();
- else
- minmass += particlesToShower[ix]->progenitor()->mass();
- }
- else {
- mIn = particlesToShower[ix]->progenitor()->mass();
- }
- }
- // throw exception if decay can't happen
- if ( minmass > mIn ) {
- throw Exception() << "Evolver.cc: Mass of decaying particle is "
- << "below constituent masses of decay products."
- << Exception::eventerror;
- }
- }
- bool reWeighting = _reWeight && hard && ShowerHandler::currentHandler()->firstInteraction();
- double eventWeight=0.;
- unsigned int nTryReWeight(0);
- // create random particle vector
- 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 ( 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());
- }
- }
- }
- // 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(!isFSRadiationON()) continue;
- // perform shower
- progenitor()->hasEmitted(startTimeLikeShower(interaction_));
- }
- // initial-state radiation
- else {
- if(!isISRadiationON()) continue;
- // hard process
- if(hard) {
- // get the PDF
- setBeamParticle(_progenitor->beam());
- assert(beamParticle());
- // perform the shower
- // set the beam particle
- tPPtr beamparticle=progenitor()->original();
- if(!beamparticle->parents().empty())
- beamparticle=beamparticle->parents()[0];
- // generate the shower
- progenitor()->hasEmitted(startSpaceLikeShower(beamparticle,
- 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();
- }
- progenitor()->progenitor()->scales().EW = progenitor()->progenitor()->mass();
- progenitor()->progenitor()->scales().EW_noAO = progenitor()->progenitor()->mass();
- // perform the shower
- progenitor()->hasEmitted(startSpaceLikeDecayShower(maxScales,minmass,
- interaction_));
- }
- }
- // do the kinematic reconstruction, checking if it worked
- reconstructed = hard ?
- showerModel()->kinematicsReconstructor()->
- reconstructHardJets (currentTree(),intrinsicpT(),interactions_[inter],
- switchRecon && ntry>maximumTries()/2) :
- showerModel()->kinematicsReconstructor()->
- reconstructDecayJets(currentTree(),interactions_[inter]);
- 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();
- }
- }
- }
- // do the kinematic reconstruction, checking if it worked
- reconstructed = hard ?
- showerModel()->kinematicsReconstructor()->
- reconstructHardJets (currentTree(),intrinsicpT(),interaction_,
- switchRecon && ntry>maximumTries()/2) :
- showerModel()->kinematicsReconstructor()->
- reconstructDecayJets(currentTree(),interaction_);
- }
- while(!reconstructed&&maximumTries()>++ntry);
- // check if failed to generate the shower
- if(ntry==maximumTries()) {
- if(hard)
- throw ShowerHandler::ShowerTriesVeto(ntry);
- else
- throw Exception() << "Failed to generate the shower after "
- << ntry << " attempts in Evolver::showerDecay()"
- << Exception::eventerror;
- }
- // 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 Evolver:: convertHardTree(bool hard,ShowerInteraction::Type 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
- showerModel()->partnerFinder()->
- setInitialEvolutionScales(particles,!hard,type,!_hardtree);
-}
-
-Branching Evolver::selectTimeLikeBranching(tShowerParticlePtr particle,
- ShowerInteraction::Type 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::Type 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;
- }
- 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=
- branch->sudakov()->createFinalStateBranching(branch->scale(),
- branch->children()[0]->z(),
- branch->phi(),
- branch->children()[0]->pT());
- 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 Evolver::selectSpaceLikeDecayBranching(tShowerParticlePtr particle,
- const ShowerParticle::EvolutionScales & maxScales,
- Energy minmass,ShowerInteraction::Type 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::Type 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=
- branch->sudakov()->createDecayBranching(branch->scale(),
- branch->children()[0]->z(),
- branch->phi(),
- branch->children()[0]->pT());
- 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;
-}
diff --git a/Shower/Base/Evolver.h b/Shower/Base/Evolver.h
deleted file mode 100644
--- a/Shower/Base/Evolver.h
+++ /dev/null
@@ -1,935 +0,0 @@
-// -*- C++ -*-
-//
-// Evolver.h is a part of Herwig - A multi-purpose Monte Carlo event generator
-// Copyright (C) 2002-2011 The Herwig Collaboration
-//
-// Herwig is licenced under version 2 of the GPL, see COPYING for details.
-// Please respect the MCnet academic guidelines, see GUIDELINES for details.
-//
-#ifndef HERWIG_Evolver_H
-#define HERWIG_Evolver_H
-//
-// This is the declaration of the Evolver class.
-//
-
-#include "ThePEG/Interface/Interfaced.h"
-#include "Herwig/Shower/SplittingFunctions/SplittingGenerator.h"
-#include "ShowerModel.h"
-#include "ThePEG/PDF/BeamParticleData.h"
-#include "ShowerTree.h"
-#include "ShowerProgenitor.fh"
-#include "Herwig/Shower/ShowerHandler.fh"
-#include "Branching.h"
-#include "ShowerVeto.h"
-#include "FullShowerVeto.h"
-#include "HardTree.h"
-#include "ThePEG/Handlers/XComb.h"
-#include "Evolver.fh"
-#include "Herwig/MatrixElement/HwMEBase.h"
-#include "Herwig/Decay/HwDecayerBase.h"
-#include "Herwig/MatrixElement/Matchbox/Matching/ShowerApproximation.h"
-#include "Herwig/Utilities/Statistic.h"
-
-namespace Herwig {
-
-using namespace ThePEG;
-
-/**\ingroup Shower
- * Exception class
- * used to communicate failure of QED shower
- */
-struct InteractionVeto {};
-
-/** \ingroup Shower
- * The Evolver class class performs the sohwer evolution of hard scattering
- * and decay processes in Herwig.
- *
- * @see \ref EvolverInterfaces "The interfaces"
- * defined for Evolver.
- */
-class Evolver: public Interfaced {
-
-/**
- * The ShowerHandler is a friend to set some parameters at initialisation
- */
-friend class ShowerHandler;
-
-public:
-
- /**
- * Pointer to an XComb object
- */
- typedef Ptr<XComb>::pointer XCPtr;
-
-public:
-
- /**
- * Default Constructor
- */
- Evolver() : _maxtry(100), _meCorrMode(1), _hardVetoMode(1),
- _hardVetoRead(0), _reconOpt(0),
- _hardVetoReadOption(false),
- _iptrms(ZERO), _beta(0.), _gamma(ZERO), _iptmax(),
- _limitEmissions(0), _initialenhance(1.), _finalenhance(1.),
- _nReWeight(100), _reWeight(false),
- interaction_(ShowerInteraction::QEDQCD), _trunc_Mode(true), _hardEmissionMode(0),
- _spinOpt(1), _softOpt(2), _hardPOWHEG(false),
- theFactorizationScaleFactor(1.0),
- theRenormalizationScaleFactor(1.0), muPt(ZERO),
- _maxTryFSR(100000),_maxFailFSR(100),_fracFSR(0.001),
- _nFSR(0), _nFailedFSR(0)
- {}
-
- /**
- * Members to perform the shower
- */
- //@{
- /**
- * Perform the shower of the hard process
- */
- virtual void showerHardProcess(ShowerTreePtr,XCPtr);
-
- /**
- * Perform the shower of a decay
- */
- virtual void showerDecay(ShowerTreePtr);
- //@}
-
- /**
- * Access to the flags and shower variables
- */
- //@{
- /**
- * Is there any showering switched on
- */
- bool showeringON() const { return isISRadiationON() || isFSRadiationON(); }
-
- /**
- * It returns true/false if the initial-state radiation is on/off.
- */
- bool isISRadiationON() const { return _splittingGenerator->isISRadiationON(); }
-
- /**
- * It returns true/false if the final-state radiation is on/off.
- */
- bool isFSRadiationON() const { return _splittingGenerator->isFSRadiationON(); }
-
- /**
- * Get the ShowerModel
- */
- ShowerModelPtr showerModel() const {return _model;}
-
- /**
- * Get the SplittingGenerator
- */
- tSplittingGeneratorPtr splittingGenerator() const { return _splittingGenerator; }
-
- /**
- * Mode for hard emissions
- */
- int hardEmissionMode() const {return _hardEmissionMode;}
-
- /**
- * Switch on or off hard vetoes
- */
- void restrictPhasespace(bool yes) {
- if ( yes )
- _hardVetoMode = 1;
- else
- _hardVetoMode = 0;
- }
-
- /**
- * Switch on or off hard veto scale from muF
- */
- void hardScaleIsMuF(bool yes) {
- if ( yes )
- _hardVetoRead = 1;
- else
- _hardVetoRead = 0;
- }
- //@}
-
- /**
- * Connect the Hard and Shower trees
- */
- virtual void connectTrees(ShowerTreePtr showerTree, HardTreePtr hardTree, bool hard );
-
- /**
- * Access to switches for spin correlations
- */
- //@{
- /**
- * Spin Correlations
- */
- unsigned int spinCorrelations() const {
- return _spinOpt;
- }
-
- /**
- * Soft correlations
- */
- unsigned int softCorrelations() const {
- return _softOpt;
- }
-
- /**
- * Any correlations
- */
- bool correlations() const {
- return _spinOpt!=0||_softOpt!=0;
- }
- //@}
-
-
- /**
- * Set the factorization scale factor
- */
- void factorizationScaleFactor(double f) {
- if ( f == theFactorizationScaleFactor )
- return;
- theFactorizationScaleFactor = f;
- splittingGenerator()->factorizationScaleFactor(f);
- }
-
- /**
- * Set the renormalization scale factor
- */
- void renormalizationScaleFactor(double f) {
- if ( f == theRenormalizationScaleFactor )
- return;
- theRenormalizationScaleFactor = f;
- splittingGenerator()->renormalizationScaleFactor(f);
- }
-
-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:
-
- /**
- * Perform the shower
- */
- void doShowering(bool hard,XCPtr);
-
- /**
- * Generate the hard matrix element correction
- */
- virtual void hardMatrixElementCorrection(bool);
-
- /**
- * Generate the hardest emission
- */
- virtual void hardestEmission(bool hard);
-
- /**
- * Extract the particles to be showered, set the evolution scales
- * and apply the hard matrix element correction
- * @param hard Whether this is a hard process or decay
- * @return The particles to be showered
- */
- virtual vector<ShowerProgenitorPtr> setupShower(bool hard);
-
- /**
- * set the colour partners
- */
- virtual void setEvolutionPartners(bool hard,ShowerInteraction::Type,
- bool clear);
-
- /**
- * Methods to perform the evolution of an individual particle, including
- * recursive calling on the products
- */
- //@{
- /**
- * It does the forward evolution of the time-like input particle
- * (and recursively for all its radiation products).
- * accepting only emissions which conforms to the showerVariables
- * and soft matrix element correction.
- * If at least one emission has occurred then the method returns true.
- * @param particle The particle to be showered
- */
- virtual bool timeLikeShower(tShowerParticlePtr particle, ShowerInteraction::Type,
- Branching fb, bool first);
-
- /**
- * It does the backward evolution of the space-like input particle
- * (and recursively for all its time-like radiation products).
- * accepting only emissions which conforms to the showerVariables.
- * If at least one emission has occurred then the method returns true
- * @param particle The particle to be showered
- * @param beam The beam particle
- */
- virtual bool spaceLikeShower(tShowerParticlePtr particle,PPtr beam,
- ShowerInteraction::Type);
-
- /**
- * If does the forward evolution of the input on-shell particle
- * involved in a decay
- * (and recursively for all its time-like radiation products).
- * accepting only emissions which conforms to the showerVariables.
- * @param particle The particle to be showered
- * @param maxscale The maximum scale for the shower.
- * @param minimumMass The minimum mass of the final-state system
- */
- virtual bool
- spaceLikeDecayShower(tShowerParticlePtr particle,
- const ShowerParticle::EvolutionScales & maxScales,
- Energy minimumMass,ShowerInteraction::Type,
- Branching fb);
-
- /**
- * Truncated shower from a time-like particle
- */
- virtual bool truncatedTimeLikeShower(tShowerParticlePtr particle,
- HardBranchingPtr branch,
- ShowerInteraction::Type type,
- Branching fb, bool first);
-
- /**
- * Truncated shower from a space-like particle
- */
- virtual bool truncatedSpaceLikeShower(tShowerParticlePtr particle,PPtr beam,
- HardBranchingPtr branch,
- ShowerInteraction::Type type);
-
- /**
- * Truncated shower from a time-like particle
- */
- virtual bool truncatedSpaceLikeDecayShower(tShowerParticlePtr particle,
- const ShowerParticle::EvolutionScales & maxScales,
- Energy minimumMass, HardBranchingPtr branch,
- ShowerInteraction::Type type, Branching fb);
- //@}
-
- /**
- * Switches for matrix element corrections
- */
- //@{
- /**
- * Any ME correction?
- */
- bool MECOn(bool hard) const {
- return ( _hardEmissionMode == 0 ||
- (!hard && _hardEmissionMode ==-1) ) &&
- _meCorrMode > 0;
- }
-
- /**
- * Any hard ME correction?
- */
- bool hardMEC(bool hard) const {
- return ( _hardEmissionMode == 0 ||
- (!hard && _hardEmissionMode ==-1) ) &&
- (_meCorrMode == 1 || _meCorrMode == 2);
- }
-
- /**
- * Any soft ME correction?
- */
- bool softMEC() const {
- return ( _hardEmissionMode == 0 ||
- (_currenttree->isDecay() && _hardEmissionMode ==-1) ) &&
- (_meCorrMode == 1 || _meCorrMode > 2);
- }
- //@}
-
- /**
- * Is the truncated shower on?
- */
- bool isTruncatedShowerON() const {return _trunc_Mode;}
-
- /**
- * Switch for intrinsic pT
- */
- //@{
- /**
- * Any intrinsic pT?
- */
- bool ipTon() const {
- return _iptrms != ZERO || ( _beta == 1.0 && _gamma != ZERO && _iptmax !=ZERO );
- }
- //@}
-
- /**@name Additional shower vetoes */
- //@{
- /**
- * Insert a veto.
- */
- void addVeto (ShowerVetoPtr v) { _vetoes.push_back(v); }
-
- /**
- * Remove a veto.
- */
- void removeVeto (ShowerVetoPtr v) {
- vector<ShowerVetoPtr>::iterator vit = find(_vetoes.begin(),_vetoes.end(),v);
- if (vit != _vetoes.end())
- _vetoes.erase(vit);
- }
-
- //@}
-
- /**
- * Switches for vetoing hard emissions
- */
- //@{
- /**
- * Vetos on?
- */
- bool hardVetoOn() const { return _hardVetoMode > 0; }
-
- /**
- * veto hard emissions in IS shower?
- */
- bool hardVetoIS() const { return _hardVetoMode == 1 || _hardVetoMode == 2; }
-
- /**
- * veto hard emissions in FS shower?
- */
- bool hardVetoFS() const { return _hardVetoMode == 1 || _hardVetoMode > 2; }
-
- /**
- * veto hard emissions according to lastScale from XComb?
- */
- bool hardVetoXComb() const {return (_hardVetoRead == 1);}
-
- /**
- * Returns true if the hard veto read-in is to be applied to only
- * the primary collision and false otherwise.
- */
- bool hardVetoReadOption() const {return _hardVetoReadOption;}
- //@}
-
- /**
- * Enhancement factors for radiation needed to generate the soft matrix
- * element correction.
- */
- //@{
- /**
- * Access the enhancement factor for initial-state radiation
- */
- double initialStateRadiationEnhancementFactor() const { return _initialenhance; }
-
- /**
- * Access the enhancement factor for final-state radiation
- */
- double finalStateRadiationEnhancementFactor() const { return _finalenhance; }
-
- /**
- * Set the enhancement factor for initial-state radiation
- */
- void initialStateRadiationEnhancementFactor(double in) { _initialenhance=in; }
-
- /**
- * Set the enhancement factor for final-state radiation
- */
- void finalStateRadiationEnhancementFactor(double in) { _finalenhance=in; }
- //@}
-
- /**
- * Access to set/get the HardTree currently beinging showered
- */
- //@{
- /**
- * The HardTree currently being showered
- */
- tHardTreePtr hardTree() {return _hardtree;}
-
- /**
- * The HardTree currently being showered
- */
- void hardTree(tHardTreePtr in) {_hardtree = in;}
- //@}
-
- /**
- * Access/set the beam particle for the current initial-state shower
- */
- //@{
- /**
- * Get the beam particle data
- */
- Ptr<BeamParticleData>::const_pointer beamParticle() const { return _beam; }
-
- /**
- * Set the beam particle data
- */
- void setBeamParticle(Ptr<BeamParticleData>::const_pointer in) { _beam=in; }
- //@}
-
- /**
- * Set/Get the current tree being evolverd for inheriting classes
- */
- //@{
- /**
- * Get the tree
- */
- tShowerTreePtr currentTree() { return _currenttree; }
-
- /**
- * Set the tree
- */
- void currentTree(tShowerTreePtr tree) { _currenttree=tree; }
-
- //@}
-
- /**
- * Access the maximum number of attempts to generate the shower
- */
- unsigned int maximumTries() const { return _maxtry; }
-
- /**
- * Set/Get the ShowerProgenitor for the current shower
- */
- //@{
- /**
- * Access the progenitor
- */
- ShowerProgenitorPtr progenitor() { return _progenitor; }
-
- /**
- * Set the progenitor
- */
- void progenitor(ShowerProgenitorPtr in) { _progenitor=in; }
- //@}
-
- /**
- * Calculate the intrinsic \f$p_T\f$.
- */
- virtual void generateIntrinsicpT(vector<ShowerProgenitorPtr>);
-
- /**
- * Access to the intrinsic \f$p_T\f$ for inheriting classes
- */
- map<tShowerProgenitorPtr,pair<Energy,double> > & intrinsicpT() { return _intrinsic; }
-
- /**
- * find the maximally allowed pt acc to the hard process.
- */
- void setupMaximumScales(const vector<ShowerProgenitorPtr> &,XCPtr);
-
- /**
- * find the relevant hard scales for profile scales.
- */
- void setupHardScales(const vector<ShowerProgenitorPtr> &,XCPtr);
-
- /**
- * Return the relevant hard scale to be used in the profile scales
- */
- Energy hardScale() const {
- return muPt;
- }
-
- /**
- * Convert the HardTree into an extra shower emission
- */
- void convertHardTree(bool hard,ShowerInteraction::Type type);
-
-protected:
-
- /**
- * Start the shower of a timelike particle
- */
- virtual bool startTimeLikeShower(ShowerInteraction::Type);
-
- /**
- * Update of the time-like stuff
- */
- void updateHistory(tShowerParticlePtr particle);
-
- /**
- * Start the shower of a spacelike particle
- */
- virtual bool startSpaceLikeShower(PPtr,ShowerInteraction::Type);
-
- /**
- * Start the shower of a spacelike particle
- */
- virtual bool
- startSpaceLikeDecayShower(const ShowerParticle::EvolutionScales & maxScales,
- Energy minimumMass,ShowerInteraction::Type);
-
- /**
- * Select the branching for the next time-like emission
- */
- Branching selectTimeLikeBranching(tShowerParticlePtr particle,
- ShowerInteraction::Type type,
- HardBranchingPtr branch);
-
- /**
- * Select the branching for the next space-like emission in a decay
- */
- Branching selectSpaceLikeDecayBranching(tShowerParticlePtr particle,
- const ShowerParticle::EvolutionScales & maxScales,
- Energy minmass,ShowerInteraction::Type type,
- HardBranchingPtr branch);
- /**
- * Create the timelike child of a branching
- */
- ShowerParticleVector createTimeLikeChildren(tShowerParticlePtr particle,
- IdList ids);
-
- /**
- * Vetos for the timelike shower
- */
- virtual bool timeLikeVetoed(const Branching &,ShowerParticlePtr);
-
- /**
- * Vetos for the spacelike shower
- */
- virtual bool spaceLikeVetoed(const Branching &,ShowerParticlePtr);
-
- /**
- * Vetos for the spacelike shower
- */
- virtual bool spaceLikeDecayVetoed(const Branching &,ShowerParticlePtr);
-
- /**
- * Only generate the hard emission, for testing only.
- */
- bool hardOnly() const {return _limitEmissions==3;}
-
-public:
-
- /** @name MC@NLO diagnostics */
- //@{
-
- /**
- * True, if Matchbox MC@NLO S-event
- */
- bool wasMCatNLOSEvent() const { return isMCatNLOSEvent; }
-
- /**
- * True, if matchbox MC@NLO H-event
- */
- bool wasMCatNLOHEvent() const { return isMCatNLOHEvent; }
-
- //@}
-
-protected:
-
- /** @name Clone Methods. */
- //@{
- /**
- * Make a simple clone of this object.
- * @return a pointer to the new object.
- */
- virtual IBPtr clone() const;
-
- /** Make a clone of this object, possibly modifying the cloned object
- * to make it sane.
- * @return a pointer to the new object.
- */
- virtual IBPtr fullclone() const;
- //@}
-
-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();
- //@}
-
-private:
-
- /**
- * The assignment operator is private and must never be called.
- * In fact, it should not even be implemented.
- */
- Evolver & operator=(const Evolver &);
-
-private:
-
- /**
- * Pointer to the model for the shower evolution model
- */
- ShowerModelPtr _model;
-
- /**
- * Pointer to the splitting generator
- */
- SplittingGeneratorPtr _splittingGenerator;
-
- /**
- * Maximum number of tries to generate the shower of a particular tree
- */
- unsigned int _maxtry;
-
- /**
- * Matrix element correction switch
- */
- unsigned int _meCorrMode;
-
- /**
- * Hard emission veto switch
- */
- unsigned int _hardVetoMode;
-
- /**
- * Hard veto to be read switch
- */
- unsigned int _hardVetoRead;
-
- /**
- * Control of the reconstruction option
- */
- unsigned int _reconOpt;
-
- /**
- * If hard veto pT scale is being read-in this determines
- * whether the read-in value is applied to primary and
- * secondary (MPI) scatters or just the primary one, with
- * the usual computation of the veto being performed for
- * the secondary (MPI) scatters.
- */
- bool _hardVetoReadOption;
-
- /**
- * rms intrinsic pT of Gaussian distribution
- */
- Energy _iptrms;
-
- /**
- * Proportion of inverse quadratic intrinsic pT distribution
- */
- double _beta;
-
- /**
- * Parameter for inverse quadratic: 2*Beta*Gamma/(sqr(Gamma)+sqr(intrinsicpT))
- */
- Energy _gamma;
-
- /**
- * Upper bound on intrinsic pT for inverse quadratic
- */
- Energy _iptmax;
-
- /**
- * Limit the number of emissions for testing
- */
- unsigned int _limitEmissions;
-
- /**
- * The progenitor of the current shower
- */
- ShowerProgenitorPtr _progenitor;
-
- /**
- * Matrix element
- */
- HwMEBasePtr _hardme;
-
- /**
- * Decayer
- */
- HwDecayerBasePtr _decayme;
-
- /**
- * The ShowerTree currently being showered
- */
- ShowerTreePtr _currenttree;
-
- /**
- * The HardTree currently being showered
- */
- HardTreePtr _hardtree;
-
- /**
- * Radiation enhancement factors for use with the veto algorithm
- * if needed by the soft matrix element correction
- */
- //@{
- /**
- * Enhancement factor for initial-state radiation
- */
- double _initialenhance;
-
- /**
- * Enhancement factor for final-state radiation
- */
- double _finalenhance;
- //@}
-
- /**
- * The beam particle data for the current initial-state shower
- */
- Ptr<BeamParticleData>::const_pointer _beam;
-
- /**
- * Storage of the intrinsic \f$p_t\f$ of the particles
- */
- map<tShowerProgenitorPtr,pair<Energy,double> > _intrinsic;
-
- /**
- * Vetoes
- */
- vector<ShowerVetoPtr> _vetoes;
-
- /**
- * Full Shower Vetoes
- */
- vector<FullShowerVetoPtr> _fullShowerVetoes;
-
- /**
- * Number of iterations for reweighting
- */
- unsigned int _nReWeight;
-
- /**
- * Whether or not we are reweighting
- */
- bool _reWeight;
-
- /**
- * number of IS emissions
- */
- unsigned int _nis;
-
- /**
- * Number of FS emissions
- */
- unsigned int _nfs;
-
- /**
- * The option for wqhich interactions to use
- */
- ShowerInteraction::Type interaction_;
-
- /**
- * Truncated shower switch
- */
- bool _trunc_Mode;
-
- /**
- * Count of the number of truncated emissions
- */
- unsigned int _truncEmissions;
-
- /**
- * Mode for the hard emissions
- */
- int _hardEmissionMode;
-
- /**
- * Option to include spin correlations
- */
- unsigned int _spinOpt;
-
- /**
- * Option for the kernal for soft correlations
- */
- unsigned int _softOpt;
-
- /**
- * Option for hard radiation in POWHEG events
- */
- bool _hardPOWHEG;
-
- /**
- * True, if Matchbox MC@NLO S-event
- */
- bool isMCatNLOSEvent;
-
- /**
- * True, if matchbox MC@NLO H-event
- */
- bool isMCatNLOHEvent;
-
- /**
- * True, if Matchbox Powheg S-event
- */
- bool isPowhegSEvent;
-
- /**
- * True, if matchbox Powheg H-event
- */
- bool isPowhegHEvent;
-
- /**
- * The shower approximation to provide the hard scale profile
- */
- Ptr<ShowerApproximation>::tptr theShowerApproximation;
-
- /**
- * The factorization scale factor.
- */
- double theFactorizationScaleFactor;
-
- /**
- * The renormalization scale factor.
- */
- double theRenormalizationScaleFactor;
-
- /**
- * True if no warnings about incorrect hard emission
- * mode setting have been issued yet
- */
- static bool _hardEmissionModeWarn;
-
- /**
- * True if no warnings about missing truncated shower
- * have been issued yet
- */
- static bool _missingTruncWarn;
-
- /**
- * The relevant hard scale to be used in the profile scales
- */
- Energy muPt;
-
- /**
- * Maximum number of emission attempts for FSR
- */
- unsigned int _maxTryFSR;
-
- /**
- * Maximum number of failures for FSR generation
- */
- unsigned int _maxFailFSR;
-
- /**
- * Failure fraction for FSR generation
- */
- double _fracFSR;
-
- /**
- * Counter for number of FSR emissions
- */
- unsigned int _nFSR;
-
- /**
- * Counter for the number of failed events due to FSR emissions
- */
- unsigned int _nFailedFSR;
-};
-
-}
-
-#endif /* HERWIG_Evolver_H */
diff --git a/Shower/QTilde/Base/Evolver.fh b/Shower/QTilde/Base/Evolver.fh
deleted file mode 100644
--- a/Shower/QTilde/Base/Evolver.fh
+++ /dev/null
@@ -1,22 +0,0 @@
-// -*- C++ -*-
-//
-// This is the forward declaration of the Evolver class.
-//
-#ifndef HERWIG_Evolver_FH
-#define HERWIG_Evolver_FH
-
-#include "ThePEG/Config/Pointers.h"
-
-namespace Herwig {
-
-class Evolver;
-
-}
-
-namespace ThePEG {
-
-ThePEG_DECLARE_POINTERS(Herwig::Evolver,EvolverPtr);
-
-}
-
-#endif
diff --git a/Shower/QTilde/Base/PartnerFinder.h b/Shower/QTilde/Base/PartnerFinder.h
--- a/Shower/QTilde/Base/PartnerFinder.h
+++ b/Shower/QTilde/Base/PartnerFinder.h
@@ -1,235 +1,234 @@
// -*- C++ -*-
//
// PartnerFinder.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 The Herwig Collaboration
//
// Herwig is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_PartnerFinder_H
#define HERWIG_PartnerFinder_H
//
// This is the declaration of the PartnerFinder class.
//
#include "Herwig/Shower/QTilde/ShowerConfig.h"
#include "ThePEG/Interface/Interfaced.h"
-#include "Evolver.fh"
#include "PartnerFinder.fh"
namespace Herwig {
using namespace ThePEG;
/**
* typedef of a pair of particle for calculating the evolution scales
*/
typedef pair<tShowerParticlePtr,tShowerParticlePtr> ShowerPPair;
/** \ingroup Shower
*
* This class is responsible of two related tasks:
* - it finds the partners
* - for each pair of partners (and interaction therefore)
* it sets the initial evolution scales of both of them.
*
* In general the finding of the partners is performed by this class but
* the calculation of the initial evolution scales should be implemented
* for different shower evolution models in classes inheriting from this one.
* Notice that a given particle has, in general, a different partner
* for each different interaction; however, given a partner, its
* initial evolution scale, Q, is purely a kinematical relationship
* between the pair, without dependence on the dynamics (i.e. type of interaction).
*
* @see \ref PartnerFinderInterfaces "The interfaces"
* defined for PartnerFinder.
*/
class PartnerFinder: public Interfaced {
public:
/**
* The default constructor.
*/
PartnerFinder() : partnerMethod_(0), QEDPartner_(0), scaleChoice_(0) {}
/**
* Given in input a collection of particles (ShowerParticle objects),
* each of these methods set the initial evolution scales of those particles,
* between the ones given in input, that do not have yet their
* evolution scale set.
* The input collection of particles can be either the full collection of
* showering particles (kept in the main class ShowerHandler,
* in the case isDecayCase is false, or simply, in the case isDecayCase
* is true, the decaying particle and its decay products.
* The methods returns true, unless something wrong (inconsistencies,
* or undefined values) happens.
*
* These methods are virtual but in most cases inheriting classes should not
* need to overide them as they simply find the relevant partner and call
* one of the calculateScale members to calculate the scale.
*/
//@{
/**
* Set the initial scales
* @param particles The particles to be considered
* @param isDecayCase Whether or not this is a decay
* @param setPartners Whether to set the colour partners or just the scales
*/
virtual void setInitialEvolutionScales(const ShowerParticleVector &particles,
const bool isDecayCase,
ShowerInteraction::Type,
const bool setPartners=true);
//@}
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:
/**
* Members to set the scales for different interactions
*/
//@{
/**
* Set initial scales for a QCD interaction
*/
virtual void setInitialQCDEvolutionScales(const ShowerParticleVector &particles,
const bool isDecayCase,
const bool setPartners=true);
/**
* Set initial scales for a QED interaction
*/
virtual void setInitialQEDEvolutionScales(const ShowerParticleVector &particles,
const bool isDecayCase,
const bool setPartners=true);
/**
* Set initial scales for a EW interaction
*/
virtual void setInitialEWEvolutionScales(const ShowerParticleVector &particles,
const bool isDecayCase,
const bool setPartners=true);
//@}
/**
* Find the QCD partners
* @param particle The particle to find the partners for
* @param particles The full set of particles to search
*/
vector< pair<ShowerPartnerType::Type, tShowerParticlePtr> >
findQCDPartners(tShowerParticlePtr particle, const ShowerParticleVector &particles);
/**
* Find the QED partners
* @param particle The particle to find the partners for
* @param particles The full set of particles to search
*/
vector< pair<double, tShowerParticlePtr> >
findQEDPartners(tShowerParticlePtr particle, const ShowerParticleVector &particles,
const bool isDecayCase);
/**
* Find the EW partners
* @param particle The particle to find the partners for
* @param particles The full set of particles to search
*/
vector< pair<double, tShowerParticlePtr> >
findEWPartners(tShowerParticlePtr particle, const ShowerParticleVector &particles,
const bool isDecayCase);
/**
* Given a pair of particles, supposedly partners w.r.t. an interaction,
* this method returns their initial evolution scales as a pair.
* If something wrong happens, it returns the null (ZERO,ZERO) pair.
* This method is used by the above setXXXInitialEvolutionScales
* methods.
* These methods must be overiden in inheriting classes
*/
//@{
/**
* General method to calculate the initial evolution scales
*/
virtual pair<Energy,Energy> calculateInitialEvolutionScales(const ShowerPPair &,
const bool isDecayCase);
/**
* Calculate the initial evolution scales for two final-state particles
*/
virtual pair<Energy,Energy> calculateFinalFinalScales(const ShowerPPair &)=0;
/**
* Calculate the initial evolution scales for two initial-state particles
*/
virtual pair<Energy,Energy> calculateInitialInitialScales(const ShowerPPair &)=0;
/**
* Calculate the initial evolution scales for one initial
* and one final-state particles
*/
virtual pair<Energy,Energy> calculateInitialFinalScales(const ShowerPPair &,
const bool isDecayCase)=0;
//@}
protected:
/**
* Find weakling interacting particles
*/
bool weaklyInteracting(tcPDPtr pd) {
long id = abs(pd->id());
return ( id==ParticleID::Wplus || id ==ParticleID::Z0 || (id>=1 && id<=6 ) || (id>=11 && id<=16));
}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
PartnerFinder & operator=(const PartnerFinder &);
private:
/**
* Method for choosing colour partner
*/
int partnerMethod_;
/**
* Choice for the QED radiation partner
*/
int QEDPartner_;
/**
* Choice of the scale
*/
int scaleChoice_;
};
}
#endif /* HERWIG_PartnerFinder_H */

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 6:33 PM (1 d, 20 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805612
Default Alt Text
(153 KB)

Event Timeline