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 OrderedParticles; - - /** - * Cached lookup of decay modes. - * Generator::findDecayMode() is not efficient. - */ - tDMPtr findDecayMode(const string & tag) { - static map cache; - map::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 -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 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 - interfaceSplitGen("SplittingGenerator", - "A reference to the SplittingGenerator object", - &Herwig::Evolver::_splittingGenerator, - false, false, true, false); - - static Reference interfaceShowerModel - ("ShowerModel", - "The pointer to the object which defines the shower evolution model.", - &Evolver::_model, false, false, true, false, false); - - static Parameter 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 interfaceNReWeight - ("NReWeight", - "The number of attempts for the shower when reweighting", - &Evolver::_nReWeight, 100, 10, 10000, - false, false, Interface::limited); - - static Switch 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 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 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 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 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 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 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 ifaceiptmax - ("IntrinsicPtIptmax", - "Upper bound on intrinsic pT for inverse quadratic", - &Evolver::_iptmax,GeV, ZERO, ZERO, 100000.0*GeV, - false, false, Interface::limited); - - static RefVector ifaceVetoes - ("Vetoes", - "The vetoes to be checked during showering", - &Evolver::_vetoes, -1, - false,false,true,true,false); - - static RefVector 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 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 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 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 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 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 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 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 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 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 interfaceMaxFailFSR - ("MaxFailFSR", - "Maximum number of failures generating the FSR", - &Evolver::_maxFailFSR, 100, 1, 100000000, - false, false, Interface::limited); - - - static Parameter 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 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;ixprogenitor()->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 pt = make_pair(ipt,UseRandom::rnd(Constants::twopi)); - _intrinsic[particlesToShower[ix]] = pt; - } -} - -void Evolver::setupMaximumScales(const vector & p, - XCPtr xcomb) { - // let POWHEG events radiate freely - if(_hardEmissionMode>0&&hardTree()) { - vector::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::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::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::const_iterator ckt = p.begin(); - for (; ckt != p.end(); ckt++) (*ckt)->maxHardPt(ptmax); -} - -void Evolver::setupHardScales(const vector & 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::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::tptr subme; - Ptr::tptr me; - Ptr::tptr dipme; - - Ptr::ptr sxc = dynamic_ptr_cast::ptr>(xcomb); - - if ( sxc ) { - subme = dynamic_ptr_cast::tptr>(sxc->matrixElement()); - me = dynamic_ptr_cast::tptr>(sxc->matrixElement()); - dipme = dynamic_ptr_cast::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(xcomb); - if(lastXC) { - _hardme = dynamic_ptr_cast(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;ixabandonChild(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 & 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 & 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 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::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::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(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;ixabandonChild(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 & 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 & 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 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 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 particles = - currentTree()->extractProgenitorParticles(); - // clear the partners if needed - if(clear) { - for(unsigned int ix=0;ixpartner(ShowerParticlePtr()); - particles[ix]->clearPartners(); - } - } - // sort out the colour partners - if(hardTree()) { - // find the partner - for(unsigned int ix=0;ixparticles()[particles[ix]]->colourPartner(); - if(!partner) continue; - for(map::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::const_iterator - eit=hardTree()->particles().end(); - for(unsigned int ix=0;ix::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_noAOscales().QCD_c_noAOscales().QCD_ac_noAOscales().EWchildren().empty()) { - ShowerParticleVector theChildren; - for(unsigned int ix=0;ixchildren().size();++ix) { - ShowerParticlePtr part = dynamic_ptr_cast - (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;ixprogenitor()->initializeFinalState(); - if(hardTree()) { - map::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::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::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::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::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::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::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 outGluon; - vector 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 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 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(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 inBranch,hardBranch; - for(map::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::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 FSROut, ISROut; - set::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;ixabandonChild(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 & 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 & 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 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 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;ixabandonChild(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 & 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 & 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::iterator cit=hardTree->branchings().begin(); - cit!=hardTree->branchings().end();++cit) { - // Sudakovs for ISR - if((**cit).parent()&&(**cit).status()==HardBranching::Incoming) { - ++_nis; - vector 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 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::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 particlesToShower; - for(map::const_iterator - cit=showerTree->incomingLines().begin(); - cit!=showerTree->incomingLines().end();++cit ) - particlesToShower.push_back(cit->first); - // extract the showering particles - for(map::const_iterator - cit=showerTree->outgoingLines().begin(); - cit!=showerTree->outgoingLines().end();++cit ) - particlesToShower.push_back(cit->first); - // match them - map partners; - for(set::const_iterator bit=hardTree->branchings().begin(); - bit!=hardTree->branchings().end();++bit) { - Energy2 dmin( 1e30*GeV2 ); - ShowerProgenitorPtr partner; - for(vector::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::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 >::const_iterator - tit = showerTree->treelinks().begin(); - tit != showerTree->treelinks().end();++tit) { - ShowerTreePtr decayTree = tit->first; - map::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 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;ixprogenitor()->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 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;ixprogenitor()->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; ixprogenitor()->spinInfo(); - if(spin && spin->decayVertex() && - dynamic_ptr_cast(spin->decayVertex())) { - spin->decayVertex(VertexPtr()); - } - } - } - // loop over particles - for(unsigned int ix=0;ixprogenitor()->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::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(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 cmap; - // incoming particles - for(map::const_iterator - cit=currentTree()->incomingLines().begin();cit!=currentTree()->incomingLines().end();++cit) { - map::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::const_iterator - cit=currentTree()->outgoingLines().begin();cit!=currentTree()->outgoingLines().end();++cit) { - map >::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::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 particles = - currentTree()->extractProgenitorParticles(); - // clear the partners - for(unsigned int ix=0;ixpartner(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::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 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::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::const_pointer beamParticle() const { return _beam; } - - /** - * Set the beam particle data - */ - void setBeamParticle(Ptr::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); - - /** - * Access to the intrinsic \f$p_T\f$ for inheriting classes - */ - map > & intrinsicpT() { return _intrinsic; } - - /** - * find the maximally allowed pt acc to the hard process. - */ - void setupMaximumScales(const vector &,XCPtr); - - /** - * find the relevant hard scales for profile scales. - */ - void setupHardScales(const vector &,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::const_pointer _beam; - - /** - * Storage of the intrinsic \f$p_t\f$ of the particles - */ - map > _intrinsic; - - /** - * Vetoes - */ - vector _vetoes; - - /** - * Full Shower Vetoes - */ - vector _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::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 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 > 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 > 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 > 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 calculateInitialEvolutionScales(const ShowerPPair &, const bool isDecayCase); /** * Calculate the initial evolution scales for two final-state particles */ virtual pair calculateFinalFinalScales(const ShowerPPair &)=0; /** * Calculate the initial evolution scales for two initial-state particles */ virtual pair calculateInitialInitialScales(const ShowerPPair &)=0; /** * Calculate the initial evolution scales for one initial * and one final-state particles */ virtual pair 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 */