diff --git a/Shower/Base/Evolver.cc b/Shower/Base/Evolver.cc --- a/Shower/Base/Evolver.cc +++ b/Shower/Base/Evolver.cc @@ -1,3220 +1,3225 @@ // -*- C++ -*- // // Evolver.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2011 The Herwig Collaboration // // Herwig is licenced under version 2 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the Evolver class. // #include "Evolver.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/RefVector.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/Shower/Base/ShowerParticle.h" #include "ThePEG/Utilities/EnumIO.h" #include "ShowerKinematics.h" #include "ThePEG/PDT/EnumParticles.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Handlers/EventHandler.h" #include "ThePEG/Utilities/Throw.h" #include "ShowerTree.h" #include "ShowerProgenitor.h" #include "KinematicsReconstructor.h" #include "PartnerFinder.h" #include "ThePEG/Handlers/StandardXComb.h" #include "ThePEG/PDT/DecayMode.h" #include "Herwig/Shower/ShowerHandler.h" #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 << _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 >> _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()))); } 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, 1000, 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 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 interfaceInteractionBothAtOnce + 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 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 particle, 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 tcPDPtr pdata[2]; for(unsigned int ix=0;ix<2;++ix) pdata[ix]=getParticleData(ids[ix+1]); if(particle->id()!=ids[0]) { for(unsigned int ix=0;ix<2;++ix) { tPDPtr cc(pdata[ix]->CC()); if(cc) pdata[ix]=cc; } } ShowerParticleVector children; for(unsigned int ix=0;ix<2;++ix) { children.push_back(new_ptr(ShowerParticle(pdata[ix],true))); if(children[ix]->id()==_progenitor->id()&&!pdata[ix]->stable()) children[ix]->set5Momentum(Lorentz5Momentum(_progenitor->progenitor()->mass())); else children[ix]->set5Momentum(Lorentz5Momentum(pdata[ix]->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); // generate phi fb.kinematics->phi(fb.sudakov->generatePhiForward(*particle,fb.ids,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]!=ParticleID::g) q2 += sqr(vm[0]); masses[ix+1] = sqrt(q2); } else { masses[ix+1] = virtualMasses[ix+1]; } } masses[0] = fb.ids[0]!=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]={getParticleData(bb.ids[0]), getParticleData(bb.ids[2])}; if(particle->id()!=bb.ids[1]) { if(part[0]->CC()) part[0]=part[0]->CC(); if(part[1]->CC()) part[1]=part[1]->CC(); } // Now create the actual particles, make the otherChild a final state // particle, while the newParent is not ShowerParticlePtr newParent=new_ptr(ShowerParticle(part[0],false)); ShowerParticlePtr otherChild = new_ptr(ShowerParticle(part[1],true,true)); ShowerParticleVector theChildren; theChildren.push_back(particle); theChildren.push_back(otherChild); //this updates the evolution scale particle->showerKinematics()-> updateParent(newParent, theChildren,bb.type); // 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,false); 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]!=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()); // if not hard generate phi if(!fb.hard) fb.kinematics->phi(fb.sudakov->generatePhiForward(*particle,fb.ids,fb.kinematics)); // 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]!=ParticleID::g) q2 += sqr(vm[0]); masses[ix+1] = sqrt(q2); } else { masses[ix+1] = virtualMasses[ix+1]; } } masses[0] = fb.ids[0]!=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] = getParticleData( bb.ids[0] ); part[1] = getParticleData( bb.ids[2] ); //is emitter anti-particle if( particle->id() != bb.ids[1]) { if( part[0]->CC() ) part[0] = part[0]->CC(); if( part[1]->CC() ) part[1] = part[1]->CC(); } double zsplit = bb.kinematics->z(); // apply the vetos for the truncated shower // if doesn't carry most of momentum 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]!=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; IdList br(3); br[0] = (**cit).parent()->branchingParticle()->id(); br[1] = (**cit). branchingParticle()->id(); br[2] = (**cit).parent()->children()[0]==*cit ? (**cit).parent()->children()[1]->branchingParticle()->id() : (**cit).parent()->children()[0]->branchingParticle()->id(); BranchingList branchings = splittingGenerator()->initialStateBranchings(); if(br[1]<0&&br[0]==br[1]) { br[0] = abs(br[0]); br[1] = abs(br[1]); } else if(br[1]<0) { br[1] = -br[1]; br[2] = -br[2]; } long index = abs(br[1]); SudakovPtr sudakov; for(BranchingList::const_iterator cjt = branchings.lower_bound(index); cjt != branchings.upper_bound(index); ++cjt ) { IdList ids = cjt->second.second; if(ids[0]==br[0]&&ids[1]==br[1]&&ids[2]==br[2]) { sudakov=cjt->second.first; break; } } if(!sudakov) throw Exception() << "Can't find Sudakov for the hard emission in " << "Evolver::connectTrees() for ISR" << Exception::runerror; (**cit).parent()->sudakov(sudakov); } // Sudakovs for FSR else if(!(**cit).children().empty()) { ++_nfs; IdList br(3); br[0] = (**cit) .branchingParticle()->id(); br[1] = (**cit).children()[0]->branchingParticle()->id(); br[2] = (**cit).children()[1]->branchingParticle()->id(); BranchingList branchings = splittingGenerator()->finalStateBranchings(); if(br[0]<0) { br[0] = abs(br[0]); br[1] = abs(br[1]); br[2] = abs(br[2]); } long index = br[0]; SudakovPtr sudakov; for(BranchingList::const_iterator cjt = branchings.lower_bound(index); cjt != branchings.upper_bound(index); ++cjt ) { IdList ids = cjt->second.second; if(ids[0]==br[0]&&ids[1]==br[1]&&ids[2]==br[2]) { sudakov=cjt->second.first; break; } } if(!sudakov) throw Exception() << "Can't find Sudakov for the hard emission in " << "Evolver::connectTrees()" << Exception::runerror; (**cit).sudakov(sudakov); } } // calculate the evolution scale for(set::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; } } // 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(); } // perform the shower progenitor()->hasEmitted(startSpaceLikeDecayShower(maxScales,minmass, interaction_)); } } } // 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; } // 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; tcPDPtr pdata[2]; 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; } // get the particle data objects for(unsigned int ix=0;ix<2;++ix) pdata[ix]=getParticleData(fb.ids[ix+1]); if(particle->id()!=fb.ids[0]) { for(unsigned int ix=0;ix<2;++ix) { tPDPtr cc(pdata[ix]->CC()); if(cc) pdata[ix]=cc; } } // find the truncated line iout=0; if(pdata[0]->id()!=pdata[1]->id()) { if(pdata[0]->id()==particle->id()) iout=1; else if (pdata[1]->id()==particle->id()) iout=2; } else if(pdata[0]->id()==particle->id()) { if(fb.kinematics->z()>0.5) iout=1; else iout=2; } // apply the vetos for the truncated shower // no flavour changing branchings if(iout==0) { particle->vetoEmission(fb.type,fb.kinematics->scale()); continue; } double zsplit = iout==1 ? fb.kinematics->z() : 1-fb.kinematics->z(); // only if same interaction for forced branching 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->id(); idlist[1] = branch->children()[0]->branchingParticle()->id(); idlist[2] = branch->children()[1]->branchingParticle()->id(); 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; tcPDPtr pdata[2]; 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; } // get the particle data objects for(unsigned int ix=0;ix<2;++ix) pdata[ix]=getParticleData(fb.ids[ix+1]); if(particle->id()!=fb.ids[0]) { for(unsigned int ix=0;ix<2;++ix) { tPDPtr cc(pdata[ix]->CC()); if(cc) pdata[ix]=cc; } } // find the truncated line iout=0; if(pdata[0]->id()!=pdata[1]->id()) { if(pdata[0]->id()==particle->id()) iout=1; else if (pdata[1]->id()==particle->id()) iout=2; } else if(pdata[0]->id()==particle->id()) { if(fb.kinematics->z()>0.5) iout=1; else iout=2; } // apply the vetos for the truncated shower // no flavour changing branchings if(iout==0) { particle->vetoEmission(fb.type,fb.kinematics->scale()); continue; } 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->id(); idlist[1] = branch->children()[0]->branchingParticle()->id(); idlist[2] = branch->children()[1]->branchingParticle()->id(); // 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/PartnerFinder.cc b/Shower/Base/PartnerFinder.cc --- a/Shower/Base/PartnerFinder.cc +++ b/Shower/Base/PartnerFinder.cc @@ -1,528 +1,623 @@ // -*- C++ -*- // // PartnerFinder.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2011 The Herwig Collaboration // // Herwig is licenced under version 2 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the PartnerFinder class. // #include "PartnerFinder.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ShowerParticle.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Utilities/Debug.h" #include "ThePEG/Utilities/DescribeClass.h" using namespace Herwig; DescribeAbstractClass describePartnerFinder ("Herwig::PartnerFinder","HwShower.so"); // some useful functions to avoid using #define namespace { // return bool if final-state particle inline bool FS(const tShowerParticlePtr a) { return a->isFinalState(); } // return colour line pointer inline Ptr::transient_pointer CL(const tShowerParticlePtr a, unsigned int index=0) { return a->colourInfo()->colourLines().empty() ? ThePEG::tColinePtr() : const_ptr_cast(a->colourInfo()->colourLines()[index]); } // return colour line size inline unsigned int CLSIZE(const tShowerParticlePtr a) { return a->colourInfo()->colourLines().size(); } inline Ptr::transient_pointer ACL(const tShowerParticlePtr a, unsigned int index=0) { return a->colourInfo()->antiColourLines().empty() ? ThePEG::tColinePtr() : const_ptr_cast(a->colourInfo()->antiColourLines()[index]); } inline unsigned int ACLSIZE(const tShowerParticlePtr a) { return a->colourInfo()->antiColourLines().size(); } } void PartnerFinder::persistentOutput(PersistentOStream & os) const { os << partnerMethod_ << QEDPartner_ << scaleChoice_; } void PartnerFinder::persistentInput(PersistentIStream & is, int) { is >> partnerMethod_ >> QEDPartner_ >> scaleChoice_; } void PartnerFinder::Init() { static ClassDocumentation documentation ("This class is responsible for finding the partners for each interaction types ", "and within the evolution scale range specified by the ShowerVariables ", "then to determine the initial evolution scales for each pair of partners."); static Switch interfacePartnerMethod ("PartnerMethod", "Choice of partner finding method for gluon evolution.", &PartnerFinder::partnerMethod_, 0, false, false); static SwitchOption interfacePartnerMethodRandom (interfacePartnerMethod, "Random", "Choose partners of a gluon randomly.", 0); static SwitchOption interfacePartnerMethodMaximum (interfacePartnerMethod, "Maximum", "Choose partner of gluon with largest angle.", 1); static Switch interfaceQEDPartner ("QEDPartner", "Control of which particles to use as the partner for QED radiation", &PartnerFinder::QEDPartner_, 0, false, false); static SwitchOption interfaceQEDPartnerAll (interfaceQEDPartner, "All", "Consider all possible choices which give a positive contribution" " in the soft limit.", 0); static SwitchOption interfaceQEDPartnerIIandFF (interfaceQEDPartner, "IIandFF", "Only allow initial-initial or final-final combinations", 1); static SwitchOption interfaceQEDPartnerIF (interfaceQEDPartner, "IF", "Only allow initial-final combinations", 2); static Switch interfaceScaleChoice ("ScaleChoice", "The choice of the evolution scales", &PartnerFinder::scaleChoice_, 0, false, false); static SwitchOption interfaceScaleChoicePartner (interfaceScaleChoice, "Partner", "Scale of all interactions is that of the evolution partner", 0); static SwitchOption interfaceScaleChoiceDifferent (interfaceScaleChoice, "Different", "Allow each interaction to have different scales", 1); } void PartnerFinder::setInitialEvolutionScales(const ShowerParticleVector &particles, const bool isDecayCase, ShowerInteraction::Type type, const bool setPartners) { // clear the existing partners for(ShowerParticleVector::const_iterator cit = particles.begin(); cit != particles.end(); ++cit) (*cit)->clearPartners(); // set them if(type==ShowerInteraction::QCD) { setInitialQCDEvolutionScales(particles,isDecayCase,setPartners); } else if(type==ShowerInteraction::QED) { setInitialQEDEvolutionScales(particles,isDecayCase,setPartners); } - else { + else if(type==ShowerInteraction::EW) { + setInitialEWEvolutionScales(particles,isDecayCase,false); + } + else if(type==ShowerInteraction::QEDQCD) { setInitialQCDEvolutionScales(particles,isDecayCase,setPartners); setInitialQEDEvolutionScales(particles,isDecayCase,false); } + else if(type==ShowerInteraction::ALL) { + setInitialQCDEvolutionScales(particles,isDecayCase,setPartners); + setInitialQEDEvolutionScales(particles,isDecayCase,false); + setInitialEWEvolutionScales(particles,isDecayCase,false); + } + else + assert(false); // \todo EW scales here // print out for debugging if(Debug::level>=10) { for(ShowerParticleVector::const_iterator cit = particles.begin(); cit != particles.end(); ++cit) { generator()->log() << "Particle: " << **cit << "\n"; if(!(**cit).partner()) continue; generator()->log() << "Primary partner: " << *(**cit).partner() << "\n"; for(vector::const_iterator it= (**cit).partners().begin(); it!=(**cit).partners().end();++it) { generator()->log() << it->type << " " << it->weight << " " << it->scale/GeV << " " << *(it->partner) << "\n"; } } generator()->log() << flush; } } void PartnerFinder::setInitialQCDEvolutionScales(const ShowerParticleVector &particles, const bool isDecayCase, const bool setPartners) { // set the partners and the scales if(setPartners) { // Loop over particles and consider only coloured particles which don't // have already their colour partner fixed and that don't have children // (the latter requirement is relaxed in the case isDecayCase is true). // Build a map which has as key one of these particles (i.e. a pointer // to a ShowerParticle object) and as a corresponding value the vector // of all its possible *normal* candidate colour partners, defined as follows: // --- have colour, and no children (this is not required in the case // isDecayCase is true); // --- if both are initial (incoming) state particles, then the (non-null) colourLine() // of one of them must match the (non-null) antiColourLine() of the other. // --- if one is an initial (incoming) state particle and the other is // a final (outgoing) state particle, then both must have the // same (non-null) colourLine() or the same (non-null) antiColourLine(); // Notice that this definition exclude the special case of baryon-violating // processes (as in R-parity Susy), which will show up as particles // without candidate colour partners, and that we will be treated a part later // (this means that no modifications of the following loop is needed!) ShowerParticleVector::const_iterator cit, cjt; for(cit = particles.begin(); cit != particles.end(); ++cit) { // Skip colourless particles if(!(*cit)->data().coloured()) continue; // find the partners vector< pair > partners = findQCDPartners(*cit,particles); // must have a partner if(partners.empty()) { throw Exception() << "`Failed to make colour connections in " << "PartnerFinder::setQCDInitialEvolutionScales" << (**cit) << Exception::eventerror; } // Calculate the evolution scales for all possible pairs of of particles vector > scales; for(unsigned int ix=0;ix< partners.size();++ix) { scales.push_back(calculateInitialEvolutionScales(ShowerPPair(*cit,partners[ix].second), isDecayCase)); } // In the case of more than one candidate colour partners, // there are now two approaches to choosing the partner. The // first method is based on two assumptions: // 1) the choice of which is the colour partner is done // *randomly* between the available candidates. // 2) the choice of which is the colour partner is done // *independently* from each particle: in other words, // if for a particle "i" its selected colour partner is // the particle "j", then the colour partner of "j" // does not have to be necessarily "i". // The second method always chooses the furthest partner // for hard gluons and gluinos. // store the choice int position(-1); // random choice if( partnerMethod_ == 0 ) { // random choice of partner position = UseRandom::irnd(partners.size()); } // take the one with largest angle else if (partnerMethod_ == 1 ) { if ((*cit)->perturbative() == 1 && (*cit)->dataPtr()->iColour()==PDT::Colour8 ) { assert(partners.size()==2); // Determine largest angle double maxAngle(0.); for(unsigned int ix=0;ixmomentum().vect(). angle(partners[ix].second->momentum().vect()); if(angle>maxAngle) { maxAngle = angle; position = ix; } } } else position = UseRandom::irnd(partners.size()); } else assert(false); // set the evolution partner (*cit)->partner(partners[position].second); for(unsigned int ix=0;ixdata().coloured()) continue; // find the partners vector< pair > partners = findQCDPartners(*cit,particles); // must have a partner if(partners.empty()) { throw Exception() << "`Failed to make colour connections in " << "PartnerFinder::setQCDInitialEvolutionScales" << (**cit) << Exception::eventerror; } // Calculate the evolution scales for all possible pairs of of particles vector > scales; int position(-1); for(unsigned int ix=0;ix< partners.size();++ix) { if(partners[ix].second) position = ix; scales.push_back(calculateInitialEvolutionScales(ShowerPPair(*cit,partners[ix].second), isDecayCase)); } assert(position>=0); for(unsigned int ix=0;ixcharged()) continue; // find the potential partners vector > partners = findQEDPartners(*cit,particles,isDecayCase); if(partners.empty()) { - throw Exception() << "Failed to partner in " + throw Exception() << "Failed to find partner in " << "PartnerFinder::setQEDInitialEvolutionScales" << (**cit) << Exception::eventerror; } // calculate the probabilities double prob(0.); for(unsigned int ix=0;ixpartner()) { for(unsigned int ix=0;ixpartner()==partners[ix].second) { position = ix; break; } } } // set the partner if(setPartners||!(*cit)->partner()||position<0) { prob = UseRandom::rnd(); for(unsigned int ix=0;ixprob) { position = ix; break; } prob -= partners[ix].first; } if(position>=0&&(setPartners||!(*cit)->partner())) { (*cit)->partner(partners[position].second); } } // must have a partner - if(position<0) throw Exception() << "Failed to partner in " + if(position<0) throw Exception() << "Failed to find partner in " << "PartnerFinder::setQEDInitialEvolutionScales" << (**cit) << Exception::eventerror; // Calculate the evolution scales for all possible pairs of of particles vector > scales; for(unsigned int ix=0;ix< partners.size();++ix) { scales.push_back(calculateInitialEvolutionScales(ShowerPPair(*cit,partners[ix].second), isDecayCase)); } // store all the possible partners for(unsigned int ix=0;ix PartnerFinder:: calculateInitialEvolutionScales(const ShowerPPair &particlePair, const bool isDecayCase) { bool FS1=FS(particlePair.first),FS2= FS(particlePair.second); if(FS1 && FS2) return calculateFinalFinalScales(particlePair); else if(FS1 && !FS2) { ShowerPPair a(particlePair.second, particlePair.first); pair rval = calculateInitialFinalScales(a,isDecayCase); return pair(rval.second,rval.first); } else if(!FS1 &&FS2) return calculateInitialFinalScales(particlePair,isDecayCase); else return calculateInitialInitialScales(particlePair); } -vector< pair > +vector< pair > PartnerFinder::findQCDPartners(tShowerParticlePtr particle, const ShowerParticleVector &particles) { vector< pair > partners; ShowerParticleVector::const_iterator cjt; for(cjt = particles.begin(); cjt != particles.end(); ++cjt) { if(!(*cjt)->data().coloured() || particle==*cjt) continue; // one initial-state and one final-state particle if(FS(particle) != FS(*cjt)) { // loop over all the colours of both particles for(unsigned int ix=0; ixsourceNeighbours().first) { tColinePair cpair = col->sourceNeighbours(); for(cjt=particles.begin();cjt!=particles.end();++cjt) { if(( FS(*cjt) && ( CL(*cjt) == cpair.first || CL(*cjt) == cpair.second))|| (!FS(*cjt) && (ACL(*cjt) == cpair.first || ACL(*cjt) == cpair.second ))) { partners.push_back(make_pair(ShowerPartnerType:: QCDColourLine,*cjt)); } } } else if(col&&col->sinkNeighbours().first) { tColinePair cpair = col->sinkNeighbours(); for(cjt=particles.begin();cjt!=particles.end();++cjt) { if(( FS(*cjt) && (ACL(*cjt) == cpair.first || ACL(*cjt) == cpair.second))|| (!FS(*cjt) && ( CL(*cjt) == cpair.first || CL(*cjt) == cpair.second))) { partners.push_back(make_pair(ShowerPartnerType:: QCDColourLine,*cjt)); } } } col = ACL(particle); if(FS(particle)&&col&&col->sinkNeighbours().first) { tColinePair cpair = col->sinkNeighbours(); for(cjt=particles.begin();cjt!=particles.end();++cjt) { if(( FS(*cjt) && (ACL(*cjt) == cpair.first || ACL(*cjt) == cpair.second))|| (!FS(*cjt) && ( CL(*cjt) == cpair.first || CL(*cjt) == cpair.second ))) { partners.push_back(make_pair(ShowerPartnerType::QCDAntiColourLine,*cjt)); } } } else if(col&&col->sourceNeighbours().first) { tColinePair cpair = col->sourceNeighbours(); for(cjt=particles.begin();cjt!=particles.end();++cjt) { if(( FS(*cjt) && ( CL(*cjt) == cpair.first || CL(*cjt) == cpair.second))|| (!FS(*cjt) && (ACL(*cjt) == cpair.first ||ACL(*cjt) == cpair.second))) { partners.push_back(make_pair(ShowerPartnerType::QCDAntiColourLine,*cjt)); } } } } // return the partners return partners; } -vector< pair > +vector< pair > PartnerFinder::findQEDPartners(tShowerParticlePtr particle, const ShowerParticleVector &particles, const bool isDecayCase) { vector< pair > partners; ShowerParticleVector::const_iterator cjt; for(cjt = particles.begin(); cjt != particles.end(); ++cjt) { if(!(*cjt)->data().charged() || particle == *cjt) continue; double charge = double(particle->data().iCharge()*(*cjt)->data().iCharge()); if( FS(particle) != FS(*cjt) ) charge *=-1.; if( QEDPartner_ != 0 && !isDecayCase ) { // only include II and FF as requested if( QEDPartner_ == 1 && FS(particle) != FS(*cjt) ) continue; // only include IF is requested else if(QEDPartner_ == 2 && FS(particle) == FS(*cjt) ) continue; } // only keep positive dipoles if(charge<0.) partners.push_back(make_pair(-charge,*cjt)); } return partners; } + +void PartnerFinder::setInitialEWEvolutionScales(const ShowerParticleVector &particles, + const bool isDecayCase, + const bool setPartners) { + // loop over all the particles + for(ShowerParticleVector::const_iterator cit = particles.begin(); + cit != particles.end(); ++cit) { + // if not weakly interacting continue + if( !weaklyInteracting( (**cit).dataPtr())) + continue; + // find the potential partners + vector > partners = findEWPartners(*cit,particles,isDecayCase); + if(partners.empty()) { + throw Exception() << "Failed to find partner in " + << "PartnerFinder::setEWInitialEvolutionScales" + << (**cit) << Exception::eventerror; + } + // calculate the probabilities + double prob(0.); + for(unsigned int ix=0;ixpartner()) { + for(unsigned int ix=0;ixpartner()==partners[ix].second) { + position = ix; + break; + } + } + } + // set the partner + if(setPartners||!(*cit)->partner()||position<0) { + prob = UseRandom::rnd(); + for(unsigned int ix=0;ixprob) { + position = ix; + break; + } + prob -= partners[ix].first; + } + if(position>=0&&(setPartners||!(*cit)->partner())) { + (*cit)->partner(partners[position].second); + } + } + // must have a partner + if(position<0) throw Exception() << "Failed to find partner in " + << "PartnerFinder::setEWInitialEvolutionScales" + << (**cit) << Exception::eventerror; + // Calculate the evolution scales for all possible pairs of of particles + vector > scales; + for(unsigned int ix=0;ix< partners.size();++ix) { + scales.push_back(calculateInitialEvolutionScales(ShowerPPair(*cit,partners[ix].second), + isDecayCase)); + } + // store all the possible partners + for(unsigned int ix=0;ix > +PartnerFinder::findEWPartners(tShowerParticlePtr particle, + const ShowerParticleVector &particles, + const bool ) { + vector< pair > partners; + ShowerParticleVector::const_iterator cjt; + for(cjt = particles.begin(); cjt != particles.end(); ++cjt) { + if(!weaklyInteracting((*cjt)->dataPtr()) || + particle == *cjt) continue; + double charge = 1.; + // only keep positive dipoles + if(charge>0.) partners.push_back(make_pair(charge,*cjt)); + } + return partners; +} diff --git a/Shower/Base/PartnerFinder.h b/Shower/Base/PartnerFinder.h --- a/Shower/Base/PartnerFinder.h +++ b/Shower/Base/PartnerFinder.h @@ -1,209 +1,236 @@ // -*- 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/ShowerConfig.h" +#include "ThePEG/PDT/ParticleData.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 > + 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 > + 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 */ diff --git a/Shower/Base/ShowerParticle.h b/Shower/Base/ShowerParticle.h --- a/Shower/Base/ShowerParticle.h +++ b/Shower/Base/ShowerParticle.h @@ -1,536 +1,541 @@ // -*- C++ -*- // // ShowerParticle.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2011 The Herwig Collaboration // // Herwig is licenced under version 2 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_ShowerParticle_H #define HERWIG_ShowerParticle_H // // This is the declaration of the ShowerParticle class. // #include "ThePEG/EventRecord/Particle.h" #include "Herwig/Shower/SplittingFunctions/SplittingFunction.fh" #include "Herwig/Shower/ShowerConfig.h" #include "ShowerBasis.h" #include "ShowerKinematics.h" #include "ShowerParticle.fh" #include namespace Herwig { using namespace ThePEG; /** \ingroup Shower * This class represents a particle in the showering process. * It inherits from the Particle class of ThePEG and has some * specifics information useful only during the showering process. * * Notice that: * - for forward evolution, it is clear what is meant by parent/child; * for backward evolution, however, it depends whether we want * to keep a physical picture or a Monte-Carlo effective one. * In the former case, an incoming particle (emitting particle) * splits into an emitted particle and the emitting particle after * the emission: the latter two are then children of the * emitting particle, the parent. In the Monte-Carlo effective * picture, we have that the particle close to the hard subprocess, * with higher (space-like) virtuality, splits into an emitted particle * and the emitting particle at lower virtuality: the latter two are, * in this case, the children of the first one, the parent. However we * choose a more physical picture where the new emitting particle is the * parented of the emitted final-state particle and the original emitting * particle. * - the pointer to a SplitFun object is set only in the case * that the particle has undergone a shower emission. This is similar to * the case of the decay of a normal Particle where * the pointer to a Decayer object is set only in the case * that the particle has undergone to a decay. * In the case of particle connected directly to the hard subprocess, * there is no pointer to the hard subprocess, but there is a method * isFromHardSubprocess() which returns true only in this case. * * @see Particle * @see ShowerConfig * @see ShowerKinematics */ class ShowerParticle: public Particle { public: /** * Struct for all the info on an evolution partner */ struct EvolutionPartner { /** * Constructor */ EvolutionPartner(tShowerParticlePtr p,double w, ShowerPartnerType::Type t, Energy s) : partner(p), weight(w), type(t), scale(s) {} /** * The partner */ tShowerParticlePtr partner; /** * Weight */ double weight; /** * Type */ ShowerPartnerType::Type type; /** * The assoicated evolution scale */ Energy scale; }; /** * Struct to store the evolution scales */ struct EvolutionScales { /** * Constructor */ EvolutionScales() : QED(),QCD_c(),QCD_ac(), QED_noAO(),QCD_c_noAO(),QCD_ac_noAO(),EW(), Max_Q2(Constants::MaxEnergy2) {} /** * QED scale */ Energy QED; /** * QCD colour scale */ Energy QCD_c; /** * QCD anticolour scale */ Energy QCD_ac; /** * QED scale */ Energy QED_noAO; /** * QCD colour scale */ Energy QCD_c_noAO; /** * QCD anticolour scale */ Energy QCD_ac_noAO; /** * EW scale */ Energy EW; /** + * EW scale + */ + Energy EW_noAO; + + /** * Maximum allowed virtuality of the particle */ Energy2 Max_Q2; }; /** @name Construction and descruction functions. */ //@{ /** * Standard Constructor. Note that the default constructor is * private - there is no particle without a pointer to a * ParticleData object. * @param x the ParticleData object * @param fs Whether or not the particle is an inital or final-state particle * @param tls Whether or not the particle initiates a time-like shower */ ShowerParticle(tcEventPDPtr x, bool fs, bool tls=false) : Particle(x), _isFinalState(fs), _perturbative(0), _initiatesTLS(tls), _x(1.0), _showerKinematics(), _vMass(ZERO), _thePEGBase() {} /** * Copy constructor from a ThePEG Particle * @param x ThePEG particle * @param pert Where the particle came from * @param fs Whether or not the particle is an inital or final-state particle * @param tls Whether or not the particle initiates a time-like shower */ ShowerParticle(const Particle & x, unsigned int pert, bool fs, bool tls=false) : Particle(x), _isFinalState(fs), _perturbative(pert), _initiatesTLS(tls), _x(1.0), _showerKinematics(), _vMass(ZERO), _thePEGBase(&x) {} //@} public: /** * Set a preliminary momentum for the particle */ void setShowerMomentum(bool timelike); /** * Construct the spin info object for a shower particle */ void constructSpinInfo(bool timelike); /** * Perform any initial calculations needed after the branching has been selected */ void initializeDecay(); /** * Perform any initial calculations needed after the branching has been selected * @param parent The beam particle */ void initializeInitialState(PPtr parent); /** * Perform any initial calculations needed after the branching has been selected */ void initializeFinalState(); /** * Access/Set various flags about the state of the particle */ //@{ /** * Access the flag that tells if the particle is final state * or initial state. */ bool isFinalState() const { return _isFinalState; } /** * Access the flag that tells if the particle is initiating a * time like shower when it has been emitted in an initial state shower. */ bool initiatesTLS() const { return _initiatesTLS; } /** * Access the flag which tells us where the particle came from * This is 0 for a particle produced in the shower, 1 if the particle came * from the hard sub-process and 2 is it came from a decay. */ unsigned int perturbative() const { return _perturbative; } //@} /** * Set/Get the momentum fraction for initial-state particles */ //@{ /** * For an initial state particle get the fraction of the beam momentum */ void x(double x) { _x = x; } /** * For an initial state particle set the fraction of the beam momentum */ double x() const { return _x; } //@} /** * Set/Get methods for the ShowerKinematics objects */ //@{ /** * Access/ the ShowerKinematics object. */ const ShoKinPtr & showerKinematics() const { return _showerKinematics; } /** * Set the ShowerKinematics object. */ void showerKinematics(const ShoKinPtr in) { _showerKinematics = in; } //@} /** * Set/Get methods for the ShowerBasis objects */ //@{ /** * Access/ the ShowerBasis object. */ const ShowerBasisPtr & showerBasis() const { return _showerBasis; } /** * Set the ShowerBasis object. */ void showerBasis(const ShowerBasisPtr in, bool copy) { if(!copy) _showerBasis = in; else { _showerBasis = new_ptr(ShowerBasis()); _showerBasis->setBasis(in->pVector(),in->nVector(),in->frame()); } } //@} /** * Members relating to the initial evolution scale and partner for the particle */ //@{ /** * Veto emission at a given scale */ void vetoEmission(ShowerPartnerType::Type type, Energy scale); /** * Access to the evolution scales */ const EvolutionScales & scales() const {return scales_;} /** * Access to the evolution scales */ EvolutionScales & scales() {return scales_;} /** * Return the virtual mass\f$ */ Energy virtualMass() const { return _vMass; } /** * Set the virtual mass */ void virtualMass(Energy mass) { _vMass = mass; } /** * Return the partner */ tShowerParticlePtr partner() const { return _partner; } /** * Set the partner */ void partner(const tShowerParticlePtr partner) { _partner = partner; } /** * Get the possible partners */ vector & partners() { return partners_; } /** * Add a possible partners */ void addPartner(EvolutionPartner in ); /** * Clear the evolution partners */ void clearPartners() { partners_.clear(); } /** * Return the progenitor of the shower */ tShowerParticlePtr progenitor() const { return _progenitor; } /** * Set the progenitor of the shower */ void progenitor(const tShowerParticlePtr progenitor) { _progenitor = progenitor; } //@} /** * Members to store and provide access to variables for a specific * shower evolution scheme */ //@{ struct Parameters { Parameters() : alpha(1.), beta(), ptx(), pty(), pt() {} double alpha; double beta; Energy ptx; Energy pty; Energy pt; }; /** * Set the vector containing dimensionless variables */ Parameters & showerParameters() { return _parameters; } //@} /** * If this particle came from the hard process get a pointer to ThePEG particle * it came from */ const tcPPtr thePEGBase() const { return _thePEGBase; } public: /** * Extract the rho matrix including mapping needed in the shower */ RhoDMatrix extractRhoMatrix(bool forward); protected: /** * For a particle which came from the hard process get the spin density and * the mapping required to the basis used in the Shower * @param rho The \f$\rho\f$ matrix * @param mapping The mapping * @param showerkin The ShowerKinematics object */ bool getMapping(SpinPtr &, RhoDMatrix & map); protected: /** * Standard clone function. */ virtual PPtr clone() const; /** * Standard clone function. */ virtual PPtr fullclone() const; private: /** * The static object used to initialize the description of this class. * Indicates that this is a concrete class with persistent data. */ static ClassDescription initShowerParticle; /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ ShowerParticle & operator=(const ShowerParticle &); private: /** * Whether the particle is in the final or initial state */ bool _isFinalState; /** * Whether the particle came from */ unsigned int _perturbative; /** * Does a particle produced in the backward shower initiate a time-like shower */ bool _initiatesTLS; /** * Dimensionless parameters */ Parameters _parameters; /** * The beam energy fraction for particle's in the initial state */ double _x; /** * The shower kinematics for the particle */ ShoKinPtr _showerKinematics; /** * The shower basis for the particle */ ShowerBasisPtr _showerBasis; /** * Storage of the evolution scales */ EvolutionScales scales_; /** * Virtual mass */ Energy _vMass; /** * Partners */ tShowerParticlePtr _partner; /** * Pointer to ThePEG Particle this ShowerParticle was created from */ const tcPPtr _thePEGBase; /** * Progenitor */ tShowerParticlePtr _progenitor; /** * Partners */ vector partners_; }; inline ostream & operator<<(ostream & os, const ShowerParticle::EvolutionScales & es) { os << "Scales: QED=" << es.QED / GeV << " QCD_c=" << es.QCD_c / GeV << " QCD_ac=" << es.QCD_ac / GeV << " QED_noAO=" << es.QED_noAO / GeV << " QCD_c_noAO=" << es.QCD_c_noAO / GeV << " QCD_ac_noAO=" << es.QCD_ac_noAO / GeV << '\n'; return os; } } #include "ThePEG/Utilities/ClassTraits.h" namespace ThePEG { /** @cond TRAITSPECIALIZATIONS */ /** This template specialization informs ThePEG about the * base classes of ShowerParticle. */ template <> struct BaseClassTrait { /** Typedef of the first base class of ShowerParticle. */ typedef Particle NthBase; }; /** This template specialization informs ThePEG about the name of * the ShowerParticle class and the shared object where it is defined. */ template <> struct ClassTraits : public ClassTraitsBase { /** Return a platform-independent class name */ static string className() { return "Herwig::ShowerParticle"; } /** Create a Event object. */ static TPtr create() { return TPtr::Create(Herwig::ShowerParticle(tcEventPDPtr(),true)); } }; /** @endcond */ } #endif /* HERWIG_ShowerParticle_H */ diff --git a/Shower/SplittingFunctions/HalfHalfOneEWSplitFn.cc b/Shower/SplittingFunctions/HalfHalfOneEWSplitFn.cc --- a/Shower/SplittingFunctions/HalfHalfOneEWSplitFn.cc +++ b/Shower/SplittingFunctions/HalfHalfOneEWSplitFn.cc @@ -1,216 +1,214 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the HalfHalfOneEWSplitFn class. // #include "HalfHalfOneEWSplitFn.h" #include "ThePEG/StandardModel/StandardModelBase.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/PDT/ParticleData.h" #include "Herwig/Decay/TwoBodyDecayMatrixElement.h" using namespace Herwig; IBPtr HalfHalfOneEWSplitFn::clone() const { return new_ptr(*this); } IBPtr HalfHalfOneEWSplitFn::fullclone() const { return new_ptr(*this); } void HalfHalfOneEWSplitFn::persistentOutput(PersistentOStream & os) const { os << gZ_ << gWL_; } void HalfHalfOneEWSplitFn::persistentInput(PersistentIStream & is, int) { is >> gZ_ >> gWL_; } // The following static variable is needed for the type description system in ThePEG. DescribeClass describeHerwigHalfHalfOneEWSplitFn("Herwig::HalfHalfOneEWSplitFn", "HwShower.so"); void HalfHalfOneEWSplitFn::Init() { static ClassDocumentation documentation ("The HalfHalfOneEWSplitFn class implements the splitting q->qWand q->qZ"); } void HalfHalfOneEWSplitFn::doinit() { - cerr << "testing in do init\n"; SplittingFunction::doinit(); - cerr << "testing in do init\n"; tcSMPtr sm = generator()->standardModel(); double sw2 = sm->sin2ThetaW(); // left-handled W coupling gWL_ = 1./sqrt(2.*sw2); // Z couplings double fact = 0.25/sqrt(sw2*(1.-sw2)); for(int ix=1;ix<4;++ix) { gZ_[2*ix-1] = make_pair(fact*(sm->vd() + sm->ad()), fact*(sm->vd() - sm->ad() )); gZ_[2*ix ] = make_pair(fact*(sm->vu() + sm->au() ), fact*(sm->vu() - sm->au() )); gZ_[2*ix+9 ] = make_pair(fact*(sm->ve() + sm->ae() ), fact*(sm->ve() - sm->ae() )); gZ_[2*ix+10] = make_pair(fact*(sm->vnu() + sm->anu()), fact*(sm->vnu() - sm->anu())); } } void HalfHalfOneEWSplitFn::getCouplings(double & gL, double & gR, const IdList & ids) const { if(ids[2]==ParticleID::Z0) { map >::const_iterator it = gZ_.find(abs(ids[0])); assert(it!=gZ_.end()); gL = it->second.first ; gR = it->second.second; } else if(abs(ids[2])==ParticleID::Wplus) { gL = gWL_; } else assert(false); } double HalfHalfOneEWSplitFn::P(const double z, const Energy2 t, const IdList &ids, const bool mass) const { double gL(0.),gR(0.); getCouplings(gL,gR,ids); double val = (1. + sqr(z))/(1.-z); if(mass) { Energy m = getParticleData(ids[2])->mass(); val -= sqr(m)/t; } val *= (sqr(gL)+sqr(gR)); return colourFactor(ids)*val; } double HalfHalfOneEWSplitFn::overestimateP(const double z, const IdList & ids) const { double gL(0.),gR(0.); getCouplings(gL,gR,ids); return 2.*sqr(max(gL,gR))*colourFactor(ids)/(1.-z); } double HalfHalfOneEWSplitFn::ratioP(const double z, const Energy2 t, const IdList & ids, const bool mass) const { double gL(0.),gR(0.); getCouplings(gL,gR,ids); double val = 1. + sqr(z); if(mass) { Energy m = getParticleData(ids[2])->mass(); val -= (1.-z)*sqr(m)/t; } val *= 0.5*(sqr(gL)+sqr(gR))/sqr(max(gL,gR)); return val; } double HalfHalfOneEWSplitFn::integOverP(const double z, const IdList & ids, unsigned int PDFfactor) const { double gL(0.),gR(0.); getCouplings(gL,gR,ids); double pre = colourFactor(ids)*sqr(max(gL,gR)); switch (PDFfactor) { case 0: return -2.*pre*Math::log1m(z); case 1: return 2.*pre*log(z/(1.-z)); case 2: return 2.*pre/(1.-z); case 3: default: throw Exception() << "HalfHalfOneEWSplitFn::integOverP() invalid PDFfactor = " << PDFfactor << Exception::runerror; } } double HalfHalfOneEWSplitFn::invIntegOverP(const double r, const IdList & ids, unsigned int PDFfactor) const { double gL(0.),gR(0.); getCouplings(gL,gR,ids); double pre = colourFactor(ids)*sqr(max(gL,gR)); switch (PDFfactor) { case 0: return 1. - exp(- 0.5*r/pre); case 1: return 1./(1.-exp(-0.5*r/pre)); case 2: return 1.-2.*pre/r; case 3: default: throw Exception() << "HalfHalfOneEWSplitFn::invIntegOverP() invalid PDFfactor = " << PDFfactor << Exception::runerror; } } bool HalfHalfOneEWSplitFn::accept(const IdList &ids) const { if(ids.size()!=3) return false; if(ids[2]==ParticleID::Z0) { if(ids[0]==ids[1] && ((ids[0]>=1 && ids[0]<=6) || (ids[0]>=11&&ids[0]<=16) )) return true; } else if(abs(ids[2])==ParticleID::Wplus) { if(!((ids[0]>=1 && ids[0]<=6) || (ids[0]>=11&&ids[0]<=16) )) return false; if(!((ids[1]>=1 && ids[1]<=6) || (ids[1]>=11&&ids[1]<=16) )) return false; if(ids[0]+1!=ids[1] && ids[0]-1!=ids[1]) return false; int out = getParticleData(ids[1])->iCharge()+getParticleData(ids[2])->iCharge(); if(getParticleData(ids[0])->iCharge()==out) return true; } return false; } vector > HalfHalfOneEWSplitFn::generatePhiForward(const double, const Energy2, const IdList & , const RhoDMatrix &) { // no dependence on the spin density matrix, dependence on off-diagonal terms cancels // and rest = splitting function for Tr(rho)=1 as required by defn return vector >(1,make_pair(0,1.)); } vector > HalfHalfOneEWSplitFn::generatePhiBackward(const double, const Energy2, const IdList & , const RhoDMatrix &) { // no dependence on the spin density matrix, dependence on off-diagonal terms cancels // and rest = splitting function for Tr(rho)=1 as required by defn return vector >(1,make_pair(0,1.)); } DecayMEPtr HalfHalfOneEWSplitFn::matrixElement(const double z, const Energy2 t, const IdList & ids, const double phi, bool) { // calculate the kernal DecayMEPtr kernal(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin1Half,PDT::Spin1Half,PDT::Spin1))); Energy m = getParticleData(ids[2])->mass(); double gL(0.),gR(0.); getCouplings(gL,gR,ids); double mt = m/sqrt(t); double root = sqrt(1.-sqr(m)/t/(1-z)); double romz = sqrt(1.-z); double rz = sqrt(z); double r2 = sqrt(2.); Complex phase = exp(Complex(0.,1.)*phi); Complex cphase = conj(phase); (*kernal)(0,0,0) = -phase*root*gL/romz; (*kernal)(1,1,2) = cphase*root*gR/romz; (*kernal)(0,0,2) = cphase*z*root*gL/romz; (*kernal)(1,1,0) = -phase*z*root*gR/romz; // long terms (*kernal)(1,1,1) =-gR*mt*r2*rz/(1-z); (*kernal)(0,0,1) =-gL*mt*r2*rz/(1-z); // +- -+ terms zero due quark mass for(unsigned int ix=0;ix<3;++ix) { (*kernal)(1,0,ix) = 0.; (*kernal)(0,1,ix) = 0.; } // return the answer return kernal; } diff --git a/Shower/SplittingFunctions/SplittingFunction.cc b/Shower/SplittingFunctions/SplittingFunction.cc --- a/Shower/SplittingFunctions/SplittingFunction.cc +++ b/Shower/SplittingFunctions/SplittingFunction.cc @@ -1,940 +1,968 @@ // -*- C++ -*- // // SplittingFunction.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2011 The Herwig Collaboration // // Herwig is licenced under version 2 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the SplittingFunction class. // #include "SplittingFunction.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Utilities/EnumIO.h" #include "Herwig/Shower/Base/ShowerParticle.h" #include "ThePEG/Utilities/DescribeClass.h" using namespace Herwig; DescribeAbstractClass describeSplittingFunction ("Herwig::SplittingFunction",""); void SplittingFunction::Init() { static ClassDocumentation documentation ("The SplittingFunction class is the based class for 1->2 splitting functions" " in Herwig"); static Switch interfaceColourStructure ("ColourStructure", "The colour structure for the splitting function", &SplittingFunction::_colourStructure, Undefined, false, false); static SwitchOption interfaceColourStructureTripletTripletOctet (interfaceColourStructure, "TripletTripletOctet", "3 -> 3 8", TripletTripletOctet); static SwitchOption interfaceColourStructureOctetOctetOctet (interfaceColourStructure, "OctetOctetOctet", "8 -> 8 8", OctetOctetOctet); static SwitchOption interfaceColourStructureOctetTripletTriplet (interfaceColourStructure, "OctetTripletTriplet", "8 -> 3 3bar", OctetTripletTriplet); static SwitchOption interfaceColourStructureTripletOctetTriplet (interfaceColourStructure, "TripletOctetTriplet", "3 -> 8 3", TripletOctetTriplet); static SwitchOption interfaceColourStructureSextetSextetOctet (interfaceColourStructure, "SextetSextetOctet", "6 -> 6 8", SextetSextetOctet); static SwitchOption interfaceColourStructureChargedChargedNeutral (interfaceColourStructure, "ChargedChargedNeutral", "q -> q 0", ChargedChargedNeutral); static SwitchOption interfaceColourStructureNeutralChargedCharged (interfaceColourStructure, "NeutralChargedCharged", "0 -> q qbar", NeutralChargedCharged); static SwitchOption interfaceColourStructureChargedNeutralCharged (interfaceColourStructure, "ChargedNeutralCharged", "q -> 0 q", ChargedNeutralCharged); static SwitchOption interfaceColourStructureEW (interfaceColourStructure, "EW", "q -> q W/Z", EW); static Switch interfaceInteractionType ("InteractionType", "Type of the interaction", &SplittingFunction::_interactionType, ShowerInteraction::UNDEFINED, false, false); static SwitchOption interfaceInteractionTypeQCD (interfaceInteractionType, "QCD","QCD",ShowerInteraction::QCD); static SwitchOption interfaceInteractionTypeQED (interfaceInteractionType, "QED","QED",ShowerInteraction::QED); static SwitchOption interfaceInteractionTypeEW (interfaceInteractionType, "EW","EW",ShowerInteraction::EW); static Switch interfaceAngularOrdered ("AngularOrdered", "Whether or not this interaction is angular ordered, " "normally only g->q qbar and gamma-> f fbar are the only ones which aren't.", &SplittingFunction::angularOrdered_, true, false, false); static SwitchOption interfaceAngularOrderedYes (interfaceAngularOrdered, "Yes", "Interaction is angular ordered", true); static SwitchOption interfaceAngularOrderedNo (interfaceAngularOrdered, "No", "Interaction isn't angular ordered", false); } void SplittingFunction::persistentOutput(PersistentOStream & os) const { using namespace ShowerInteraction; os << oenum(_interactionType) << _interactionOrder << oenum(_colourStructure) << _colourFactor << angularOrdered_; } void SplittingFunction::persistentInput(PersistentIStream & is, int) { using namespace ShowerInteraction; is >> ienum(_interactionType) >> _interactionOrder >> ienum(_colourStructure) >> _colourFactor >> angularOrdered_; } void SplittingFunction::colourConnection(tShowerParticlePtr parent, tShowerParticlePtr first, tShowerParticlePtr second, ShowerPartnerType::Type partnerType, const bool back) const { if(_colourStructure==TripletTripletOctet) { if(!back) { ColinePair cparent = ColinePair(parent->colourLine(), parent->antiColourLine()); // ensure input consistency assert(( cparent.first && !cparent.second && partnerType==ShowerPartnerType::QCDColourLine) || ( !cparent.first && cparent.second && partnerType==ShowerPartnerType::QCDAntiColourLine)); // q -> q g if(cparent.first) { ColinePtr newline=new_ptr(ColourLine()); cparent.first->addColoured(second); newline->addColoured ( first); newline->addAntiColoured (second); } // qbar -> qbar g else { ColinePtr newline=new_ptr(ColourLine()); cparent.second->addAntiColoured(second); newline->addColoured(second); newline->addAntiColoured(first); } // Set progenitor first->progenitor(parent->progenitor()); second->progenitor(parent->progenitor()); } else { ColinePair cfirst = ColinePair(first->colourLine(), first->antiColourLine()); // ensure input consistency assert(( cfirst.first && !cfirst.second && partnerType==ShowerPartnerType::QCDColourLine) || ( !cfirst.first && cfirst.second && partnerType==ShowerPartnerType::QCDAntiColourLine)); // q -> q g if(cfirst.first) { ColinePtr newline=new_ptr(ColourLine()); cfirst.first->addAntiColoured(second); newline->addColoured(second); newline->addColoured(parent); } // qbar -> qbar g else { ColinePtr newline=new_ptr(ColourLine()); cfirst.second->addColoured(second); newline->addAntiColoured(second); newline->addAntiColoured(parent); } // Set progenitor parent->progenitor(first->progenitor()); second->progenitor(first->progenitor()); } } else if(_colourStructure==OctetOctetOctet) { if(!back) { ColinePair cparent = ColinePair(parent->colourLine(), parent->antiColourLine()); // ensure input consistency assert(cparent.first&&cparent.second); // ensure first gluon is hardest if( first->id()==second->id() && parent->showerKinematics()->z()<0.5 ) swap(first,second); // colour line radiates if(partnerType==ShowerPartnerType::QCDColourLine) { // The colour line is radiating ColinePtr newline=new_ptr(ColourLine()); cparent.first->addColoured(second); cparent.second->addAntiColoured(first); newline->addColoured(first); newline->addAntiColoured(second); } // anti colour line radiates else if(partnerType==ShowerPartnerType::QCDAntiColourLine) { ColinePtr newline=new_ptr(ColourLine()); cparent.first->addColoured(first); cparent.second->addAntiColoured(second); newline->addColoured(second); newline->addAntiColoured(first); } else assert(false); } else { ColinePair cfirst = ColinePair(first->colourLine(), first->antiColourLine()); // ensure input consistency assert(cfirst.first&&cfirst.second); // The colour line is radiating if(partnerType==ShowerPartnerType::QCDColourLine) { ColinePtr newline=new_ptr(ColourLine()); cfirst.first->addAntiColoured(second); cfirst.second->addAntiColoured(parent); newline->addColoured(parent); newline->addColoured(second); } // anti colour line radiates else if(partnerType==ShowerPartnerType::QCDAntiColourLine) { ColinePtr newline=new_ptr(ColourLine()); cfirst.first->addColoured(parent); cfirst.second->addColoured(second); newline->addAntiColoured(second); newline->addAntiColoured(parent); } else assert(false); } } else if(_colourStructure == OctetTripletTriplet) { if(!back) { ColinePair cparent = ColinePair(parent->colourLine(), parent->antiColourLine()); // ensure input consistency assert(cparent.first&&cparent.second); cparent.first ->addColoured ( first); cparent.second->addAntiColoured(second); // Set progenitor first->progenitor(parent->progenitor()); second->progenitor(parent->progenitor()); } else { ColinePair cfirst = ColinePair(first->colourLine(), first->antiColourLine()); // ensure input consistency assert(( cfirst.first && !cfirst.second) || (!cfirst.first && cfirst.second)); // g -> q qbar if(cfirst.first) { ColinePtr newline=new_ptr(ColourLine()); cfirst.first->addColoured(parent); newline->addAntiColoured(second); newline->addAntiColoured(parent); } // g -> qbar q else { ColinePtr newline=new_ptr(ColourLine()); cfirst.second->addAntiColoured(parent); newline->addColoured(second); newline->addColoured(parent); } // Set progenitor parent->progenitor(first->progenitor()); second->progenitor(first->progenitor()); } } else if(_colourStructure == TripletOctetTriplet) { if(!back) { ColinePair cparent = ColinePair(parent->colourLine(), parent->antiColourLine()); // ensure input consistency assert(( cparent.first && !cparent.second) || (!cparent.first && cparent.second)); // q -> g q if(cparent.first) { ColinePtr newline=new_ptr(ColourLine()); cparent.first->addColoured(first); newline->addColoured (second); newline->addAntiColoured( first); } // qbar -> g qbar else { ColinePtr newline=new_ptr(ColourLine()); cparent.second->addAntiColoured(first); newline->addColoured ( first); newline->addAntiColoured(second); } // Set progenitor first->progenitor(parent->progenitor()); second->progenitor(parent->progenitor()); } else { ColinePair cfirst = ColinePair(first->colourLine(), first->antiColourLine()); // ensure input consistency assert(cfirst.first&&cfirst.second); // q -> g q if(parent->id()>0) { cfirst.first ->addColoured(parent); cfirst.second->addColoured(second); } else { cfirst.first ->addAntiColoured(second); cfirst.second->addAntiColoured(parent); } // Set progenitor parent->progenitor(first->progenitor()); second->progenitor(first->progenitor()); } } else if(_colourStructure==SextetSextetOctet) { //make sure we're not doing backward evolution assert(!back); //make sure something sensible assert(parent->colourLine() || parent->antiColourLine()); //get the colour lines or anti-colour lines bool isAntiColour=true; ColinePair cparent; if(parent->colourLine()) { cparent = ColinePair(const_ptr_cast(parent->colourInfo()->colourLines()[0]), const_ptr_cast(parent->colourInfo()->colourLines()[1])); isAntiColour=false; } else { cparent = ColinePair(const_ptr_cast(parent->colourInfo()->antiColourLines()[0]), const_ptr_cast(parent->colourInfo()->antiColourLines()[1])); } //check for sensible input // assert(cparent.first && cparent.second); // sextet has 2 colour lines if(!isAntiColour) { //pick at random which of the colour topolgies to take double topology = UseRandom::rnd(); if(topology < 0.25) { ColinePtr newline=new_ptr(ColourLine()); cparent.first->addColoured(second); cparent.second->addColoured(first); newline->addColoured(first); newline->addAntiColoured(second); } else if(topology >=0.25 && topology < 0.5) { ColinePtr newline=new_ptr(ColourLine()); cparent.first->addColoured(first); cparent.second->addColoured(second); newline->addColoured(first); newline->addAntiColoured(second); } else if(topology >= 0.5 && topology < 0.75) { ColinePtr newline=new_ptr(ColourLine()); cparent.first->addColoured(second); cparent.second->addColoured(first); newline->addColoured(first); newline->addAntiColoured(second); } else { ColinePtr newline=new_ptr(ColourLine()); cparent.first->addColoured(first); cparent.second->addColoured(second); newline->addColoured(first); newline->addAntiColoured(second); } } // sextet has 2 anti-colour lines else { double topology = UseRandom::rnd(); if(topology < 0.25){ ColinePtr newline=new_ptr(ColourLine()); cparent.first->addAntiColoured(second); cparent.second->addAntiColoured(first); newline->addAntiColoured(first); newline->addColoured(second); } else if(topology >=0.25 && topology < 0.5) { ColinePtr newline=new_ptr(ColourLine()); cparent.first->addAntiColoured(first); cparent.second->addAntiColoured(second); newline->addAntiColoured(first); newline->addColoured(second); } else if(topology >= 0.5 && topology < 0.75) { ColinePtr newline=new_ptr(ColourLine()); cparent.first->addAntiColoured(second); cparent.second->addAntiColoured(first); newline->addAntiColoured(first); newline->addColoured(second); } else { ColinePtr newline=new_ptr(ColourLine()); cparent.first->addAntiColoured(first); cparent.second->addAntiColoured(second); newline->addAntiColoured(first); newline->addColoured(second); } } } else if(_colourStructure == ChargedChargedNeutral) { if(!parent->data().coloured()) return; if(!back) { ColinePair cparent = ColinePair(parent->colourLine(), parent->antiColourLine()); // q -> q g if(cparent.first) { cparent.first->addColoured(first); } // qbar -> qbar g if(cparent.second) { cparent.second->addAntiColoured(first); } } else { ColinePair cfirst = ColinePair(first->colourLine(), first->antiColourLine()); // q -> q g if(cfirst.first) { cfirst.first->addColoured(parent); } // qbar -> qbar g if(cfirst.second) { cfirst.second->addAntiColoured(parent); } } } else if(_colourStructure == ChargedNeutralCharged) { if(!parent->data().coloured()) return; if(!back) { ColinePair cparent = ColinePair(parent->colourLine(), parent->antiColourLine()); // q -> q g if(cparent.first) { cparent.first->addColoured(second); } // qbar -> qbar g if(cparent.second) { cparent.second->addAntiColoured(second); } } else { if (second->dataPtr()->iColour()==PDT::Colour3 ) { ColinePtr newline=new_ptr(ColourLine()); newline->addColoured(second); newline->addColoured(parent); } else if (second->dataPtr()->iColour()==PDT::Colour3bar ) { ColinePtr newline=new_ptr(ColourLine()); newline->addAntiColoured(second); newline->addAntiColoured(parent); } } } else if(_colourStructure == NeutralChargedCharged ) { if(!back) { if(first->dataPtr()->coloured()) { ColinePtr newline=new_ptr(ColourLine()); if(first->dataPtr()->iColour()==PDT::Colour3) { newline->addColoured (first ); newline->addAntiColoured(second); } else if (first->dataPtr()->iColour()==PDT::Colour3bar) { newline->addColoured (second); newline->addAntiColoured(first ); } else assert(false); } } else { ColinePair cfirst = ColinePair(first->colourLine(), first->antiColourLine()); // gamma -> q qbar if(cfirst.first) { cfirst.first->addAntiColoured(second); } // gamma -> qbar q else if(cfirst.second) { cfirst.second->addColoured(second); } else assert(false); } } else if(_colourStructure == EW) { if(!parent->data().coloured()) return; if(!back) { ColinePair cparent = ColinePair(parent->colourLine(), parent->antiColourLine()); // q -> q g if(cparent.first) { cparent.first->addColoured(first); } // qbar -> qbar g if(cparent.second) { cparent.second->addAntiColoured(first); } } else { ColinePair cfirst = ColinePair(first->colourLine(), first->antiColourLine()); // q -> q g if(cfirst.first) { cfirst.first->addColoured(parent); } // qbar -> qbar g if(cfirst.second) { cfirst.second->addAntiColoured(parent); } } } else { assert(false); } } void SplittingFunction::doinit() { Interfaced::doinit(); assert(_interactionType!=ShowerInteraction::UNDEFINED); assert((_colourStructure>0&&_interactionType==ShowerInteraction::QCD) || (_colourStructure<0&&(_interactionType==ShowerInteraction::QED || _interactionType==ShowerInteraction::EW)) ); if(_colourFactor>0.) return; // compute the colour factors if need if(_colourStructure==TripletTripletOctet) { _colourFactor = 4./3.; } else if(_colourStructure==OctetOctetOctet) { _colourFactor = 3.; } else if(_colourStructure==OctetTripletTriplet) { _colourFactor = 0.5; } else if(_colourStructure==TripletOctetTriplet) { _colourFactor = 4./3.; } else if(_colourStructure==SextetSextetOctet) { _colourFactor = 10./3.; } else if(_colourStructure<0) { _colourFactor = 1.; } else { assert(false); } } bool SplittingFunction::checkColours(const IdList & ids) const { tcPDPtr pd[3]={getParticleData(ids[0]), getParticleData(ids[1]), getParticleData(ids[2])}; if(_colourStructure==TripletTripletOctet) { if(ids[0]!=ids[1]) return false; if((pd[0]->iColour()==PDT::Colour3||pd[0]->iColour()==PDT::Colour3bar) && pd[2]->iColour()==PDT::Colour8) return true; return false; } else if(_colourStructure==OctetOctetOctet) { for(unsigned int ix=0;ix<3;++ix) { if(pd[ix]->iColour()!=PDT::Colour8) return false; } return true; } else if(_colourStructure==OctetTripletTriplet) { if(pd[0]->iColour()!=PDT::Colour8) return false; if(pd[1]->iColour()==PDT::Colour3&&pd[2]->iColour()==PDT::Colour3bar) return true; if(pd[1]->iColour()==PDT::Colour3bar&&pd[2]->iColour()==PDT::Colour3) return true; return false; } else if(_colourStructure==TripletOctetTriplet) { if(ids[0]!=ids[2]) return false; if((pd[0]->iColour()==PDT::Colour3||pd[0]->iColour()==PDT::Colour3bar) && pd[1]->iColour()==PDT::Colour8) return true; return false; } else if(_colourStructure==SextetSextetOctet) { if(ids[0]!=ids[1]) return false; if((pd[0]->iColour()==PDT::Colour6 || pd[0]->iColour()==PDT::Colour6bar) && pd[2]->iColour()==PDT::Colour8) return true; return false; } else if(_colourStructure==ChargedChargedNeutral) { if(ids[0]!=ids[1]) return false; if(pd[2]->iCharge()!=0) return false; if(pd[0]->iCharge()==pd[1]->iCharge()) return true; return false; } else if(_colourStructure==ChargedNeutralCharged) { if(ids[0]!=ids[2]) return false; if(pd[1]->iCharge()!=0) return false; if(pd[0]->iCharge()==pd[2]->iCharge()) return true; return false; } else if(_colourStructure==NeutralChargedCharged) { if(ids[1]!=-ids[2]) return false; if(pd[0]->iCharge()!=0) return false; if(pd[1]->iCharge()==-pd[2]->iCharge()) return true; return false; } else { assert(false); } return false; } namespace { bool hasColour(tPPtr p) { PDT::Colour colour = p->dataPtr()->iColour(); return colour==PDT::Colour3 || colour==PDT::Colour8 || colour == PDT::Colour6; } bool hasAntiColour(tPPtr p) { PDT::Colour colour = p->dataPtr()->iColour(); return colour==PDT::Colour3bar || colour==PDT::Colour8 || colour == PDT::Colour6bar; } } void SplittingFunction::evaluateFinalStateScales(ShowerPartnerType::Type partnerType, Energy scale, double z, tShowerParticlePtr parent, tShowerParticlePtr emitter, tShowerParticlePtr emitted) { // identify emitter and emitted double zEmitter = z, zEmitted = 1.-z; bool bosonSplitting(false); // special for g -> gg, particle highest z is emitter if(emitter->id() == emitted->id() && emitter->id() == parent->id() && zEmitted > zEmitter) { swap(zEmitted,zEmitter); swap( emitted, emitter); } // otherwise if particle ID same else if(emitted->id()==parent->id()) { swap(zEmitted,zEmitter); swap( emitted, emitter); } // no real emitter/emitted else if(emitter->id()!=parent->id()) { bosonSplitting = true; } // may need to add angularOrder flag here // now the various scales // QED if(partnerType==ShowerPartnerType::QED) { assert(colourStructure()==ChargedChargedNeutral || colourStructure()==ChargedNeutralCharged || colourStructure()==NeutralChargedCharged ); // normal case if(!bosonSplitting) { assert(colourStructure()==ChargedChargedNeutral || colourStructure()==ChargedNeutralCharged ); // set the scales // emitter emitter->scales().QED = zEmitter*scale; emitter->scales().QED_noAO = scale; emitter->scales().QCD_c = min(scale,parent->scales().QCD_c ); emitter->scales().QCD_c_noAO = min(scale,parent->scales().QCD_c_noAO ); emitter->scales().QCD_ac = min(scale,parent->scales().QCD_ac ); emitter->scales().QCD_ac_noAO = min(scale,parent->scales().QCD_ac_noAO); + emitter->scales().EW = min(scale,parent->scales().EW ); + emitter->scales().EW_noAO = min(scale,parent->scales().EW_noAO ); // emitted emitted->scales().QED = zEmitted*scale; emitted->scales().QED_noAO = scale; emitted->scales().QCD_c = ZERO; emitted->scales().QCD_c_noAO = ZERO; emitted->scales().QCD_ac = ZERO; emitted->scales().QCD_ac_noAO = ZERO; + emitted->scales().EW = min(scale,parent->scales().EW ); + emitted->scales().EW_noAO = min(scale,parent->scales().EW_noAO ); } // gamma -> f fbar else { assert(colourStructure()==NeutralChargedCharged ); // emitter - emitter->scales().QED = zEmitter*scale; - emitter->scales().QED_noAO = scale; + emitter->scales().QED = zEmitter*scale; + emitter->scales().QED_noAO = scale; if(hasColour(emitter)) { emitter->scales().QCD_c = zEmitter*scale; emitter->scales().QCD_c_noAO = scale; } if(hasAntiColour(emitter)) { emitter->scales().QCD_ac = zEmitter*scale; emitter->scales().QCD_ac_noAO = scale; } + emitter->scales().EW = zEmitter*scale; + emitter->scales().EW_noAO = scale; // emitted - emitted->scales().QED = zEmitted*scale; - emitted->scales().QED_noAO = scale; + emitted->scales().QED = zEmitted*scale; + emitted->scales().QED_noAO = scale; if(hasColour(emitted)) { emitted->scales().QCD_c = zEmitted*scale; emitted->scales().QCD_c_noAO = scale; } if(hasAntiColour(emitted)) { emitted->scales().QCD_ac = zEmitted*scale; emitted->scales().QCD_ac_noAO = scale; } + emitted->scales().EW = zEmitted*scale; + emitted->scales().EW_noAO = scale; } } // QCD else if (partnerType==ShowerPartnerType::QCDColourLine || partnerType==ShowerPartnerType::QCDAntiColourLine) { // normal case eg q -> q g and g -> g g if(!bosonSplitting) { emitter->scales().QED = min(scale,parent->scales().QED ); emitter->scales().QED_noAO = min(scale,parent->scales().QED_noAO); + emitter->scales().EW = min(scale,parent->scales().EW ); + emitter->scales().EW_noAO = min(scale,parent->scales().EW_noAO); if(partnerType==ShowerPartnerType::QCDColourLine) { emitter->scales().QCD_c = zEmitter*scale; emitter->scales().QCD_c_noAO = scale; emitter->scales().QCD_ac = min(zEmitter*scale,parent->scales().QCD_ac ); emitter->scales().QCD_ac_noAO = min( scale,parent->scales().QCD_ac_noAO); } else { emitter->scales().QCD_c = min(zEmitter*scale,parent->scales().QCD_c ); emitter->scales().QCD_c_noAO = min( scale,parent->scales().QCD_c_noAO ); emitter->scales().QCD_ac = zEmitter*scale; emitter->scales().QCD_ac_noAO = scale; } // emitted emitted->scales().QED = ZERO; emitted->scales().QED_noAO = ZERO; emitted->scales().QCD_c = zEmitted*scale; emitted->scales().QCD_c_noAO = scale; emitted->scales().QCD_ac = zEmitted*scale; emitted->scales().QCD_ac_noAO = scale; + emitted->scales().EW = min(scale,parent->scales().EW ); + emitted->scales().EW_noAO = min(scale,parent->scales().EW_noAO); } // g -> q qbar else { // emitter if(emitter->dataPtr()->charged()) { emitter->scales().QED = zEmitter*scale; emitter->scales().QED_noAO = scale; } - emitter->scales().QCD_c = zEmitter*scale; - emitter->scales().QCD_c_noAO = scale; - emitter->scales().QCD_ac = zEmitter*scale; - emitter->scales().QCD_ac_noAO = scale; + emitter->scales().EW = zEmitter*scale; + emitter->scales().EW_noAO = scale; + emitter->scales().QCD_c = zEmitter*scale; + emitter->scales().QCD_c_noAO = scale; + emitter->scales().QCD_ac = zEmitter*scale; + emitter->scales().QCD_ac_noAO = scale; // emitted if(emitted->dataPtr()->charged()) { emitted->scales().QED = zEmitted*scale; emitted->scales().QED_noAO = scale; } - emitted->scales().QCD_c = zEmitted*scale; - emitted->scales().QCD_c_noAO = scale; - emitted->scales().QCD_ac = zEmitted*scale; - emitted->scales().QCD_ac_noAO = scale; + emitted->scales().EW = zEmitted*scale; + emitted->scales().EW_noAO = scale; + emitted->scales().QCD_c = zEmitted*scale; + emitted->scales().QCD_c_noAO = scale; + emitted->scales().QCD_ac = zEmitted*scale; + emitted->scales().QCD_ac_noAO = scale; } } else if(partnerType==ShowerPartnerType::EW) { - assert(false); - // // emitter - // emitter->scales().QED = min(scale,parent->scales().QED ); - // emitter->scales().QED_noAO = min(scale,parent->scales().QED_noAO); - // emitter->scales().QCD_c = min(scale,parent->scales().QCD_c ); - // emitter->scales().QCD_c_noAO = min(scale,parent->scales().QCD_c_noAO ); - // emitter->scales().QCD_ac = min(scale,parent->scales().QCD_ac ); - // emitter->scales().QCD_ac_noAO = min(scale,parent->scales().QCD_ac_noAO); - // // emitted - // emitter->scales().QED = zEmitter*scale; - // emitter->scales().QED_noAO = scale; - // emitted->scales().QED = zEmitted*scale; - // emitted->scales().QED_noAO = scale; - // emitted->scales().QCD_c = ZERO; - // emitted->scales().QCD_c_noAO = ZERO; - // emitted->scales().QCD_ac = ZERO; - // emitted->scales().QCD_ac_noAO = ZERO; + // EW + emitter->scales().EW = zEmitter*scale; + emitter->scales().EW_noAO = scale; + emitted->scales().EW = zEmitted*scale; + emitted->scales().EW_noAO = scale; + // QED + // W radiation AO + if(emitted->dataPtr()->charged()) { + emitter->scales().QED = zEmitter*scale; + emitter->scales().QED_noAO = scale; + emitted->scales().QED = zEmitted*scale; + emitted->scales().QED_noAO = scale; + } + // Z don't + else { + emitter->scales().QED = min(scale,parent->scales().QED ); + emitter->scales().QED_noAO = min(scale,parent->scales().QED_noAO); + emitted->scales().QED = ZERO; + emitted->scales().QED_noAO = ZERO; + } + // QCD + emitter->scales().QCD_c = min(scale,parent->scales().QCD_c ); + emitter->scales().QCD_c_noAO = min(scale,parent->scales().QCD_c_noAO ); + emitter->scales().QCD_ac = min(scale,parent->scales().QCD_ac ); + emitter->scales().QCD_ac_noAO = min(scale,parent->scales().QCD_ac_noAO); + emitted->scales().QCD_c = ZERO; + emitted->scales().QCD_c_noAO = ZERO; + emitted->scales().QCD_ac = ZERO; + emitted->scales().QCD_ac_noAO = ZERO; } else assert(false); } void SplittingFunction::evaluateInitialStateScales(ShowerPartnerType::Type partnerType, Energy scale, double z, tShowerParticlePtr parent, tShowerParticlePtr spacelike, tShowerParticlePtr timelike) { // scale for time-like child Energy AOScale = (1.-z)*scale; // QED if(partnerType==ShowerPartnerType::QED) { if(parent->id()==spacelike->id()) { // parent parent ->scales().QED = scale; parent ->scales().QED_noAO = scale; parent ->scales().QCD_c = min(scale,spacelike->scales().QCD_c ); parent ->scales().QCD_c_noAO = min(scale,spacelike->scales().QCD_c_noAO ); parent ->scales().QCD_ac = min(scale,spacelike->scales().QCD_ac ); parent ->scales().QCD_ac_noAO = min(scale,spacelike->scales().QCD_ac_noAO); // timelike timelike->scales().QED = AOScale; timelike->scales().QED_noAO = scale; timelike->scales().QCD_c = ZERO; timelike->scales().QCD_c_noAO = ZERO; timelike->scales().QCD_ac = ZERO; timelike->scales().QCD_ac_noAO = ZERO; } else if(parent->id()==timelike->id()) { parent ->scales().QED = scale; parent ->scales().QED_noAO = scale; if(hasColour(parent)) { parent ->scales().QCD_c = scale; parent ->scales().QCD_c_noAO = scale; } if(hasAntiColour(parent)) { parent ->scales().QCD_ac = scale; parent ->scales().QCD_ac_noAO = scale; } // timelike timelike->scales().QED = AOScale; timelike->scales().QED_noAO = scale; if(hasColour(timelike)) { timelike->scales().QCD_c = AOScale; timelike->scales().QCD_c_noAO = scale; } if(hasAntiColour(timelike)) { timelike->scales().QCD_ac = AOScale; timelike->scales().QCD_ac_noAO = scale; } } else { parent ->scales().QED = scale; parent ->scales().QED_noAO = scale; parent ->scales().QCD_c = ZERO ; parent ->scales().QCD_c_noAO = ZERO ; parent ->scales().QCD_ac = ZERO ; parent ->scales().QCD_ac_noAO = ZERO ; // timelike timelike->scales().QED = AOScale; timelike->scales().QED_noAO = scale; if(hasColour(timelike)) { timelike->scales().QCD_c = min(AOScale,spacelike->scales().QCD_ac ); timelike->scales().QCD_c_noAO = min( scale,spacelike->scales().QCD_ac_noAO); } if(hasAntiColour(timelike)) { timelike->scales().QCD_ac = min(AOScale,spacelike->scales().QCD_c ); timelike->scales().QCD_ac_noAO = min( scale,spacelike->scales().QCD_c_noAO ); } } } // QCD else if (partnerType==ShowerPartnerType::QCDColourLine || partnerType==ShowerPartnerType::QCDAntiColourLine) { // timelike if(timelike->dataPtr()->charged()) { timelike->scales().QED = AOScale; timelike->scales().QED_noAO = scale; } if(hasColour(timelike)) { timelike->scales().QCD_c = AOScale; timelike->scales().QCD_c_noAO = scale; } if(hasAntiColour(timelike)) { timelike->scales().QCD_ac = AOScale; timelike->scales().QCD_ac_noAO = scale; } if(parent->id()==spacelike->id()) { parent ->scales().QED = min(scale,spacelike->scales().QED ); parent ->scales().QED_noAO = min(scale,spacelike->scales().QED_noAO ); parent ->scales().QCD_c = min(scale,spacelike->scales().QCD_c ); parent ->scales().QCD_c_noAO = min(scale,spacelike->scales().QCD_c_noAO ); parent ->scales().QCD_ac = min(scale,spacelike->scales().QCD_ac ); parent ->scales().QCD_ac_noAO = min(scale,spacelike->scales().QCD_ac_noAO); } else { if(parent->dataPtr()->charged()) { parent ->scales().QED = scale; parent ->scales().QED_noAO = scale; } if(hasColour(parent)) { parent ->scales().QCD_c = scale; parent ->scales().QCD_c_noAO = scale; } if(hasAntiColour(parent)) { parent ->scales().QCD_ac = scale; parent ->scales().QCD_ac_noAO = scale; } } } else assert(false); } void SplittingFunction::evaluateDecayScales(ShowerPartnerType::Type partnerType, Energy scale, double z, tShowerParticlePtr parent, tShowerParticlePtr spacelike, tShowerParticlePtr timelike) { assert(parent->id()==spacelike->id()); // angular-ordered scale for 2nd child Energy AOScale = (1.-z)*scale; // QED if(partnerType==ShowerPartnerType::QED) { // timelike timelike->scales().QED = AOScale; timelike->scales().QED_noAO = scale; timelike->scales().QCD_c = ZERO; timelike->scales().QCD_c_noAO = ZERO; timelike->scales().QCD_ac = ZERO; timelike->scales().QCD_ac_noAO = ZERO; // spacelike spacelike->scales().QED = scale; spacelike->scales().QED_noAO = scale; } // QCD else if(partnerType==ShowerPartnerType::QCDColourLine || partnerType==ShowerPartnerType::QCDAntiColourLine) { // timelike timelike->scales().QED = ZERO; timelike->scales().QED_noAO = ZERO; timelike->scales().QCD_c = AOScale; timelike->scales().QCD_c_noAO = scale; timelike->scales().QCD_ac = AOScale; timelike->scales().QCD_ac_noAO = scale; // spacelike spacelike->scales().QED = max(scale,parent->scales().QED ); spacelike->scales().QED_noAO = max(scale,parent->scales().QED_noAO ); } else assert(false); spacelike->scales().QCD_c = max(scale,parent->scales().QCD_c ); spacelike->scales().QCD_c_noAO = max(scale,parent->scales().QCD_c_noAO ); spacelike->scales().QCD_ac = max(scale,parent->scales().QCD_ac ); spacelike->scales().QCD_ac_noAO = max(scale,parent->scales().QCD_ac_noAO); }