diff --git a/Shower/Base/Evolver.cc b/Shower/Base/Evolver.cc --- a/Shower/Base/Evolver.cc +++ b/Shower/Base/Evolver.cc @@ -1,3227 +1,3222 @@ // -*- 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 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(); } progenitor()->progenitor()->scales().EW = progenitor()->progenitor()->mass(); progenitor()->progenitor()->scales().EW_noAO = progenitor()->progenitor()->mass(); // perform the shower progenitor()->hasEmitted(startSpaceLikeDecayShower(maxScales,minmass, interaction_)); } } } // do the kinematic reconstruction, checking if it worked reconstructed = hard ? showerModel()->kinematicsReconstructor()-> reconstructHardJets (currentTree(),intrinsicpT(),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/SudakovFormFactor.h b/Shower/Base/SudakovFormFactor.h --- a/Shower/Base/SudakovFormFactor.h +++ b/Shower/Base/SudakovFormFactor.h @@ -1,690 +1,698 @@ // -*- C++ -*- // // SudakovFormFactor.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_SudakovFormFactor_H #define HERWIG_SudakovFormFactor_H // // This is the declaration of the SudakovFormFactor class. // #include "ThePEG/Interface/Interfaced.h" #include "Herwig/Shower/SplittingFunctions/SplittingFunction.h" #include "Herwig/Shower/Couplings/ShowerAlpha.h" #include "Herwig/Shower/SplittingFunctions/SplittingGenerator.fh" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/PDF/BeamParticleData.h" #include "ThePEG/EventRecord/RhoDMatrix.h" #include "ThePEG/EventRecord/SpinInfo.h" #include "ShowerKinematics.fh" #include "SudakovFormFactor.fh" namespace Herwig { using namespace ThePEG; /** * A typedef for the BeamParticleData */ typedef Ptr::transient_const_pointer tcBeamPtr; /** \ingroup Shower * * This is the definition of the Sudakov form factor class. In general this * is the base class for the implementation of Sudakov form factors in Herwig. * The methods generateNextTimeBranching(), generateNextDecayBranching() and * generateNextSpaceBranching need to be implemented in classes inheriting from this * one. * * In addition a number of methods are implemented to assist with the calculation * of the form factor using the veto algorithm in classes inheriting from this one. * * In general the Sudakov form-factor, for final-state radiation, is given * by * \f[\Delta_{ba}(\tilde{q}_{i+1},\tilde{q}_i)= * \exp\left\{ * -\int^{\tilde{q}^2_i}_{\tilde{q}^2_{i+1}} * \frac{{\rm d}\tilde{q}^2}{\tilde{q}^2} * \int\frac{\alpha_S(z,\tilde{q})}{2\pi} * P_{ba}(z,\tilde{q})\Theta(p_T) * \right\}. * \f] * We can solve this to obtain the next value of the scale \f$\tilde{q}_{i+1}\f$ * given the previous value \f$\tilde{q}_i\f$ * in the following way. First we obtain a simplified form of the integrand * which is greater than or equal to the true integrand for all values of * \f$\tilde{q}\f$. * * In practice it is easiest to obtain this over estimate in pieces. The ShowerAlpha * object contains an over estimate for \f$\alpha_S\f$, the splitting function * contains both an over estimate of the spltting function and its integral * which is needed to compute the over estimate of the \f$\tilde{q}\f$ integrand, * together with an over estimate of the limit of the \f$z\f$ integral. * * This gives an overestimate of the integrand * \f[g(\tilde{q}^2) = \frac{c}{\tilde{q}^2}, \f] * where because the over estimates are chosen to be independent of \f$\tilde{q}\f$ the * parameter * \f[c = \frac{\alpha_{\rm over}}{2\pi}\int^{z_1}_{z_0}P_{\rm over}(z),\f] * is a constant independent of \f$\tilde{q}\f$. * * The guesst() member can then be used to generate generate the value of * \f$\tilde{q}^2\f$ according to this result. This is done by solving the Sudakov * form factor, with the over estimates, is equal to a random number * \f$r\f$ in the interval \f$[0,1]\f$. This gives * \f[\tilde{q}^2_{i+1}=G^{-1}\left[G(\tilde{q}^2_i)+\ln r\right],\f] * where \f$G(\tilde{q}^2)=c\ln(\tilde{q}^2)\f$ is the infinite integral * of \f$g(\tilde{q}^2)\f$ and \f$G^{-1}(x)=\exp\left(\frac{x}c\right)\f$ * is its inverse. * It this case we therefore obtain * \f[\tilde{q}^2_{i+1}=\tilde{q}^2_ir^{\frac1c}.\f] * The value of \f$z\f$ can then be calculated in a similar way * \f[z = I^{-1}\left[I(z_0)+r\left(I(z_1)-I(z_0)\right)\right],\f] * using the guessz() member, * where \f$I=\int P(z){\rm d}z\f$ and \f$I^{-1}\f$ is its inverse. * * The veto algorithm then uses rejection using the ratio of the * true value to the overestimated one to obtain the original distribution. * This is accomplished using the * - alphaSVeto() member for the \f$\alpha_S\f$ veto * - SplittingFnVeto() member for the veto on the value of the splitting function. * in general there must also be a chech that the emission is in the allowed * phase space but this is left to the inheriting classes as it will depend * on the ordering variable. * * The Sudakov form factor for the initial-scale shower is different because * it must include the PDF which guides the backward evolution. * It is given by * \f[\Delta_{ba}(\tilde{q}_{i+1},\tilde{q}_i)= * \exp\left\{ * -\int^{\tilde{q}^2_i}_{\tilde{q}^2_{i+1}} * \frac{{\rm d}\tilde{q}^2}{\tilde{q}^2} * \int\frac{\alpha_S(z,\tilde{q})}{2\pi} * P_{ba}(z,\tilde{q})\frac{x'f_a(\frac{x}z,\tilde{q}^2)}{xf_b(x,\tilde{q^2})} * \right\}, * \f] * where \f$x\f$ is the fraction of the beam momentum the parton \f$b\f$ had before * the backward evolution. * This can be solve in the same way as for the final-state branching but the constant * becomes * \f[c = \frac{\alpha_{\rm over}}{2\pi}\int^{z_1}_{z_0}P_{\rm over}(z)PDF_{\rm max},\f] * where * \f[PDF_{\rm max}=\max\frac{x'f_a(\frac{x}z,\tilde{q}^2)}{xf_b(x,\tilde{q^2})},\f] * which can be set using an interface. * In addition the PDFVeto() member then is needed to implement the relevant veto. * * @see SplittingGenerator * @see SplittingFunction * @see ShowerAlpha * @see \ref SudakovFormFactorInterfaces "The interfaces" * defined for SudakovFormFactor. */ class SudakovFormFactor: public Interfaced { /** * The SplittingGenerator is a friend to insert the particles in the * branchings at initialisation */ friend class SplittingGenerator; public: /** * The default constructor. */ SudakovFormFactor() : pdfmax_(35.0), pdffactor_(0), cutOffOption_(0), a_(0.3), b_(2.3), c_(0.3*GeV), kinCutoffScale_( 2.3*GeV ), vgcut_(0.85*GeV), vqcut_(0.85*GeV), pTmin_(1.*GeV), pT2min_(ZERO), z_( 0.0 ),phi_(0.0), pT_(), theFactorizationScaleFactor(1.0), theRenormalizationScaleFactor(1.0) {} /** * Members to generate the scale of the next branching */ //@{ /** * Return the scale of the next time-like branching. If there is no * branching then it returns ZERO. * @param startingScale starting scale for the evolution * @param ids The PDG codes of the particles in the splitting * @param cc Whether this is the charge conjugate of the branching * @param enhance The radiation enhancement factor * @param maxQ2 The maximum \f$Q^2\f$ for the emission */ virtual ShoKinPtr generateNextTimeBranching(const Energy startingScale, const IdList &ids,const bool cc, + const RhoDMatrix & rho, double enhance, Energy2 maxQ2)=0; /** * Return the scale of the next space-like decay branching. If there is no * branching then it returns ZERO. * @param startingScale starting scale for the evolution * @param stoppingScale stopping scale for the evolution * @param minmass The minimum mass allowed for the spake-like particle. * @param ids The PDG codes of the particles in the splitting * @param cc Whether this is the charge conjugate of the branching * defined. * @param enhance The radiation enhancement factor */ virtual ShoKinPtr generateNextDecayBranching(const Energy startingScale, const Energy stoppingScale, const Energy minmass, const IdList &ids, const bool cc, + const RhoDMatrix & rho, double enhance)=0; /** * Return the scale of the next space-like branching. If there is no * branching then it returns ZERO. * @param startingScale starting scale for the evolution * @param ids The PDG codes of the particles in the splitting * @param x The fraction of the beam momentum * @param cc Whether this is the charge conjugate of the branching * defined. * @param beam The beam particle * @param enhance The radiation enhancement factor */ virtual ShoKinPtr generateNextSpaceBranching(const Energy startingScale, const IdList &ids,double x, - const bool cc,double enhance, + const bool cc, + const RhoDMatrix & rho, + double enhance, tcBeamPtr beam)=0; //@} /** * Generate the azimuthal angle of the branching for forward evolution * @param particle The branching particle * @param ids The PDG codes of the particles in the branchings * @param The Shower kinematics */ virtual double generatePhiForward(ShowerParticle & particle,const IdList & ids, - ShoKinPtr kinematics)=0; + ShoKinPtr kinematics, + const RhoDMatrix & rho)=0; /** * Generate the azimuthal angle of the branching for backward evolution * @param particle The branching particle * @param ids The PDG codes of the particles in the branchings * @param The Shower kinematics */ virtual double generatePhiBackward(ShowerParticle & particle,const IdList & ids, - ShoKinPtr kinematics)=0; + ShoKinPtr kinematics, + const RhoDMatrix & rho)=0; /** * Generate the azimuthal angle of the branching for ISR in decays * @param particle The branching particle * @param ids The PDG codes of the particles in the branchings * @param The Shower kinematics */ virtual double generatePhiDecay(ShowerParticle & particle,const IdList & ids, - ShoKinPtr kinematics)=0; + ShoKinPtr kinematics, + const RhoDMatrix & rho)=0; /** * Methods to provide public access to the private member variables */ //@{ /** * Return the pointer to the SplittingFunction object. */ tSplittingFnPtr splittingFn() const { return splittingFn_; } /** * Return the pointer to the ShowerAlpha object. */ tShowerAlphaPtr alpha() const { return alpha_; } /** * The type of interaction */ inline ShowerInteraction::Type interactionType() const {return splittingFn_->interactionType();} //@} public: /** * Methods to access the kinematic variables for the branching */ //@{ /** * The energy fraction */ double z() const { return z_; } /** * The azimuthal angle */ double phi() const { return phi_; } /** * The transverse momentum */ Energy pT() const { return pT_; } //@} /** * Access the maximum weight for the PDF veto */ double pdfMax() const { return pdfmax_;} /** * Method to return the evolution scale given the * transverse momentum, \f$p_T\f$ and \f$z\f$. */ virtual Energy calculateScale(double z, Energy pt, IdList ids,unsigned int iopt)=0; /** * Method to create the ShowerKinematics object for a final-state branching */ virtual ShoKinPtr createFinalStateBranching(Energy scale,double z, double phi, Energy pt)=0; /** * Method to create the ShowerKinematics object for an initial-state branching */ virtual ShoKinPtr createInitialStateBranching(Energy scale,double z, double phi, Energy pt)=0; /** * Method to create the ShowerKinematics object for a decay branching */ virtual ShoKinPtr createDecayBranching(Energy scale,double z, double phi, Energy pt)=0; public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); //@} protected: /** * Methods to implement the veto algorithm to generate the scale of * the next branching */ //@{ /** * Value of the energy fraction for the veto algorithm * @param iopt The option for calculating z * @param ids The PDG codes of the particles in the splitting * - 0 is final-state * - 1 is initial-state for the hard process * - 2 is initial-state for particle decays */ double guessz (unsigned int iopt, const IdList &ids) const; /** * Value of the scale for the veto algorithm * @param t1 The starting valoe of the scale * @param iopt The option for calculating t * @param ids The PDG codes of the particles in the splitting * - 0 is final-state * - 1 is initial-state for the hard process * - 2 is initial-state for particle decays * @param enhance The radiation enhancement factor * @param identical Whether or not the outgoing particles are identical */ Energy2 guesst (Energy2 t1,unsigned int iopt, const IdList &ids, double enhance, bool identical) const; /** * Veto on the PDF for the initial-state shower * @param t The scale * @param x The fraction of the beam momentum * @param parton0 Pointer to the particleData for the * new parent (this is the particle we evolved back to) * @param parton1 Pointer to the particleData for the * original particle * @param beam The BeamParticleData object */ bool PDFVeto(const Energy2 t, const double x, const tcPDPtr parton0, const tcPDPtr parton1, tcBeamPtr beam) const; /** * The veto on the splitting function. * @param t The scale * @param ids The PDG codes of the particles in the splitting * @param mass Whether or not to use the massive splitting functions * @return true if vetoed */ bool SplittingFnVeto(const Energy2 t, const IdList &ids, - const bool mass) const { - return UseRandom::rnd()>splittingFn_->ratioP(z_, t, ids,mass); + const bool mass, + const RhoDMatrix & rho) const { + return UseRandom::rnd()>splittingFn_->ratioP(z_, t, ids,mass,rho); } /** * The veto on the coupling constant * @param pt2 The value of ther transverse momentum squared, \f$p_T^2\f$. * @return true if vetoed */ bool alphaSVeto(Energy2 pt2) const; //@} /** * Methods to set the kinematic variables for the branching */ //@{ /** * The energy fraction */ void z(double in) { z_=in; } /** * The azimuthal angle */ void phi(double in) { phi_=in; } /** * The transverse momentum */ void pT(Energy in) { pT_=in; } //@} /** * Set/Get the limits on the energy fraction for the splitting */ //@{ /** * Get the limits */ pair zLimits() const { return zlimits_;} /** * Set the limits */ void zLimits(pair in) { zlimits_=in; } //@} /** * Set the particles in the splittings */ void addSplitting(const IdList &); /** * Delete the particles in the splittings */ void removeSplitting(const IdList &); /** * Access the potential branchings */ const vector & particles() const { return particles_; } public: /** * @name Methods for the cut-off */ //@{ /** * The option being used */ unsigned int cutOffOption() const { return cutOffOption_; } /** * The kinematic scale */ Energy kinScale() const {return kinCutoffScale_;} /** * The virtuality cut-off on the gluon \f$Q_g=\frac{\delta-am_q}{b}\f$ * @param scale The scale \f$\delta\f$ * @param mq The quark mass \f$m_q\f$. */ Energy kinematicCutOff(Energy scale, Energy mq) const {return max((scale -a_*mq)/b_,c_);} /** * The virtualilty cut-off for gluons */ Energy vgCut() const { return vgcut_; } /** * The virtuality cut-off for everything else */ Energy vqCut() const { return vqcut_; } /** * The minimum \f$p_T\f$ for the branching */ Energy pTmin() const { return pTmin_; } /** * The square of the minimum \f$p_T\f$ */ Energy2 pT2min() const { return pT2min_; } /** * Calculate the virtual masses for a branchings */ const vector & virtualMasses(const IdList & ids); //@} /** * Set the PDF */ void setPDF(tcPDFPtr pdf, Energy scale) { pdf_ = pdf; freeze_ = scale; } private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ SudakovFormFactor & operator=(const SudakovFormFactor &); private: /** * Pointer to the splitting function for this Sudakov form factor */ SplittingFnPtr splittingFn_; /** * Pointer to the coupling for this Sudakov form factor */ ShowerAlphaPtr alpha_; /** * Maximum value of the PDF weight */ double pdfmax_; /** * List of the particles this Sudakov is used for to aid in setting up * interpolation tables if needed */ vector particles_; /** * Option for the inclusion of a factor \f$1/(1-z)\f$ in the PDF estimate */ unsigned pdffactor_; private: /** * Option for the type of cut-off to be applied */ unsigned int cutOffOption_; /** * Parameters for the default Herwig cut-off option, i.e. the parameters for * the \f$Q_g=\max(\frac{\delta-am_q}{b},c)\f$ kinematic cut-off */ //@{ /** * The \f$a\f$ parameter */ double a_; /** * The \f$b\f$ parameter */ double b_; /** * The \f$c\f$ parameter */ Energy c_; /** * Kinematic cutoff used in the parton shower phase space. */ Energy kinCutoffScale_; //@} /** * Parameters for the FORTRAN-like cut-off */ //@{ /** * The virtualilty cut-off for gluons */ Energy vgcut_; /** * The virtuality cut-off for everything else */ Energy vqcut_; //@} /** * Parameters for the \f$p_T\f$ cut-off */ //@{ /** * The minimum \f$p_T\f$ for the branching */ Energy pTmin_; /** * The square of the minimum \f$p_T\f$ */ Energy2 pT2min_; //@} private: /** * Member variables to keep the shower kinematics information * generated by a call to generateNextTimeBranching or generateNextSpaceBranching */ //@{ /** * The energy fraction */ double z_; /** * The azimuthal angle */ double phi_; /** * The transverse momentum */ Energy pT_; //@} /** * The limits of \f$z\f$ in the splitting */ pair zlimits_; /** * Stuff for the PDFs */ //@{ /** * PDf */ tcPDFPtr pdf_; /** * Freezing scale */ Energy freeze_; //@} public: /** * Get the factorization scale factor */ double factorizationScaleFactor() const { return theFactorizationScaleFactor; } /** * Set the factorization scale factor */ void factorizationScaleFactor(double f) { theFactorizationScaleFactor = f; } /** * Get the renormalization scale factor */ double renormalizationScaleFactor() const { return theRenormalizationScaleFactor; } /** * Set the renormalization scale factor */ void renormalizationScaleFactor(double f) { theRenormalizationScaleFactor = f; } private: /** * The factorization scale factor. */ double theFactorizationScaleFactor; /** * The renormalization scale factor. */ double theRenormalizationScaleFactor; }; } #endif /* HERWIG_SudakovFormFactor_H */ diff --git a/Shower/Default/QTildeSudakov.cc b/Shower/Default/QTildeSudakov.cc --- a/Shower/Default/QTildeSudakov.cc +++ b/Shower/Default/QTildeSudakov.cc @@ -1,937 +1,941 @@ // -*- C++ -*- // // QTildeSudakov.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 QTildeSudakov class. // #include "QTildeSudakov.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/PDT/ParticleData.h" #include "ThePEG/EventRecord/Event.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Repository/CurrentGenerator.h" #include "ThePEG/PDT/EnumParticles.h" #include "Herwig/Shower/Default/FS_QTildeShowerKinematics1to2.h" #include "Herwig/Shower/Default/IS_QTildeShowerKinematics1to2.h" #include "Herwig/Shower/Default/Decay_QTildeShowerKinematics1to2.h" #include "ThePEG/Utilities/DescribeClass.h" #include "Herwig/Shower/Base/ShowerVertex.h" #include "Herwig/Shower/Base/ShowerParticle.h" #include "Herwig/Shower/ShowerHandler.h" #include "Herwig/Shower/Base/Evolver.h" #include "Herwig/Shower/Base/PartnerFinder.h" #include "Herwig/Shower/Base/ShowerModel.h" #include "Herwig/Shower/Base/KinematicsReconstructor.h" using namespace Herwig; DescribeNoPIOClass describeQTildeSudakov ("Herwig::QTildeSudakov","HwShower.so"); void QTildeSudakov::Init() { static ClassDocumentation documentation ("The QTildeSudakov class implements the Sudakov form factor for ordering it" " qtilde"); } bool QTildeSudakov::guessTimeLike(Energy2 &t,Energy2 tmin,double enhance) { Energy2 told = t; // calculate limits on z and if lower>upper return if(!computeTimeLikeLimits(t)) return false; // guess values of t and z t = guesst(told,0,ids_,enhance,ids_[1]==ids_[2]); z(guessz(0,ids_)); // actual values for z-limits if(!computeTimeLikeLimits(t)) return false; if(tupper return if(!computeSpaceLikeLimits(t,x)) return false; // guess values of t and z t = guesst(told,1,ids_,enhance,ids_[1]==ids_[2]); z(guessz(1,ids_)); // actual values for z-limits if(!computeSpaceLikeLimits(t,x)) return false; if(t zLimits().second) return true; Energy2 q2 = z()*(1.-z())*t; if(ids_[0]!=ParticleID::g && ids_[0]!=ParticleID::gamma ) q2 += masssquared_[0]; if(q2>maxQ2) return true; // compute the pts Energy2 pt2 = z()*(1.-z())*q2 - masssquared_[1]*(1.-z()) - masssquared_[2]*z(); // if pt2<0 veto if(pt2 min if(tmax<=tmin) return ShoKinPtr(); // calculate next value of t using veto algorithm Energy2 t(tmax); do { if(!guessTimeLike(t,tmin,enhance)) break; } - while(PSVeto(t,maxQ2) || SplittingFnVeto(z()*(1.-z())*t,ids,true) || + while(PSVeto(t,maxQ2) || SplittingFnVeto(z()*(1.-z())*t,ids,true,rho) || alphaSVeto(splittingFn()->angularOrdered() ? sqr(z()*(1.-z()))*t : z()*(1.-z())*t)); q_ = t > ZERO ? Energy(sqrt(t)) : -1.*MeV; if(q_ < ZERO) return ShoKinPtr(); // return the ShowerKinematics object return createFinalStateBranching(q_,z(),phi(),pT()); } ShoKinPtr QTildeSudakov:: generateNextSpaceBranching(const Energy startingQ, const IdList &ids, double x,bool cc, + const RhoDMatrix & rho, double enhance, Ptr::transient_const_pointer beam) { // First reset the internal kinematics variables that can // have been eventually set in the previous call to the method. q_ = ZERO; z(0.); phi(0.); // perform the initialization Energy2 tmax(sqr(startingQ)),tmin; initialize(ids,tmin,cc); // check max > min if(tmax<=tmin) return ShoKinPtr(); // extract the partons which are needed for the PDF veto // Different order, incoming parton is id = 1, outgoing are id=0,2 tcPDPtr parton0 = getParticleData(ids[0]); tcPDPtr parton1 = getParticleData(ids[1]); if(cc) { if(parton0->CC()) parton0 = parton0->CC(); if(parton1->CC()) parton1 = parton1->CC(); } // calculate next value of t using veto algorithm Energy2 t(tmax),pt2(ZERO); do { if(!guessSpaceLike(t,tmin,x,enhance)) break; pt2=sqr(1.-z())*t-z()*masssquared_[2]; } while(z() > zLimits().second || - SplittingFnVeto((1.-z())*t/z(),ids,true) || + SplittingFnVeto((1.-z())*t/z(),ids,true,rho) || alphaSVeto(splittingFn()->angularOrdered() ? sqr(1.-z())*t : (1.-z())*t) || PDFVeto(t,x,parton0,parton1,beam) || pt2 < pT2min() ); if(t > ZERO && zLimits().first < zLimits().second) q_ = sqrt(t); else return ShoKinPtr(); pT(sqrt(pt2)); // create the ShowerKinematics and return it return createInitialStateBranching(q_,z(),phi(),pT()); } void QTildeSudakov::initialize(const IdList & ids, Energy2 & tmin,const bool cc) { ids_=ids; if(cc) { for(unsigned int ix=0;ixCC()) ids_[ix]*=-1; } } tmin = cutOffOption() != 2 ? ZERO : 4.*pT2min(); masses_ = virtualMasses(ids); masssquared_.clear(); for(unsigned int ix=0;ix0) tmin=max(masssquared_[ix],tmin); } } ShoKinPtr QTildeSudakov::generateNextDecayBranching(const Energy startingScale, - const Energy stoppingScale, - const Energy minmass, - const IdList &ids, - const bool cc, - double enhance) { + const Energy stoppingScale, + const Energy minmass, + const IdList &ids, + const bool cc, + const RhoDMatrix & rho, + double enhance) { // First reset the internal kinematics variables that can // have been eventually set in the previous call to this method. q_ = Constants::MaxEnergy; z(0.); phi(0.); // perform initialisation Energy2 tmax(sqr(stoppingScale)),tmin; initialize(ids,tmin,cc); tmin=sqr(startingScale); // check some branching possible if(tmax<=tmin) return ShoKinPtr(); // perform the evolution Energy2 t(tmin),pt2(-MeV2); do { if(!guessDecay(t,tmax,minmass,enhance)) break; pt2 = sqr(1.-z())*(t-masssquared_[0])-z()*masssquared_[2]; } - while(SplittingFnVeto((1.-z())*t/z(),ids,true)|| + while(SplittingFnVeto((1.-z())*t/z(),ids,true,rho)|| alphaSVeto(splittingFn()->angularOrdered() ? sqr(1.-z())*t : (1.-z())*t ) || pt2masssquared_[0]-sqr(minmass)); if(t > ZERO) { q_ = sqrt(t); pT(sqrt(pt2)); } else return ShoKinPtr(); phi(0.); // create the ShowerKinematics object return createDecayBranching(q_,z(),phi(),pT()); } bool QTildeSudakov::guessDecay(Energy2 &t,Energy2 tmax, Energy minmass, double enhance) { // previous scale Energy2 told = t; // overestimated limits on z if(tmax limits=make_pair(sqr(minmass/masses_[0]), 1.-sqrt(masssquared_[2]+pT2min()+ 0.25*sqr(masssquared_[2])/tm2)/tm +0.5*masssquared_[2]/tm2); zLimits(limits); if(zLimits().secondtmax||zLimits().second limits; if(ids_[0]==ParticleID::g||ids_[0]==ParticleID::gamma) { // no emission possible if(t<16.*(masssquared_[1]+pT2min())) { t=-1.*GeV2; return false; } // overestimate of the limits limits.first = 0.5*(1.-sqrt(1.-4.*sqrt((masssquared_[1]+pT2min())/t))); limits.second = 1.-limits.first; } // special case for radiated particle is gluon else if(ids_[2]==ParticleID::g||ids_[2]==ParticleID::gamma) { limits.first = sqrt((masssquared_[1]+pT2min())/t); limits.second = 1.-sqrt((masssquared_[2]+pT2min())/t); } else if(ids_[1]==ParticleID::g||ids_[1]==ParticleID::gamma) { limits.second = sqrt((masssquared_[2]+pT2min())/t); limits.first = 1.-sqrt((masssquared_[1]+pT2min())/t); } else { limits.first = (masssquared_[1]+pT2min())/t; limits.second = 1.-(masssquared_[2]+pT2min())/t; } if(limits.first>=limits.second) { t=-1.*GeV2; return false; } zLimits(limits); return true; } bool QTildeSudakov::computeSpaceLikeLimits(Energy2 & t, double x) { if (t < 1e-20 * GeV2) { t=-1.*GeV2; return false; } pair limits; // compute the limits limits.first = x; double yy = 1.+0.5*masssquared_[2]/t; limits.second = yy - sqrt(sqr(yy)-1.+pT2min()/t); // return false if lower>upper zLimits(limits); if(limits.second(particle.parents()[0]) : tShowerParticlePtr(); } else { mother = particle.children().size()==2 ? dynamic_ptr_cast(&particle) : tShowerParticlePtr(); } tShowerParticlePtr partner; while(mother) { tPPtr otherChild; if(forward) { for (unsigned int ix=0;ixchildren().size();++ix) { if(mother->children()[ix]!=child) { otherChild = mother->children()[ix]; break; } } } else { otherChild = mother->children()[1]; } tShowerParticlePtr other = dynamic_ptr_cast(otherChild); if((inter==ShowerInteraction::QCD && otherChild->dataPtr()->coloured()) || (inter==ShowerInteraction::QED && otherChild->dataPtr()->charged())) { partner = other; break; } if(forward && !other->isFinalState()) { partner = dynamic_ptr_cast(mother); break; } child = mother; if(forward) { mother = ! mother->parents().empty() ? dynamic_ptr_cast(mother->parents()[0]) : tShowerParticlePtr(); } else { if(mother->children()[0]->children().size()!=2) break; tShowerParticlePtr mtemp = dynamic_ptr_cast(mother->children()[0]); if(!mtemp) break; else mother=mtemp; } } if(!partner) { if(forward) { partner = dynamic_ptr_cast( child)->partner(); } else { if(mother) { tShowerParticlePtr parent; if(!mother->children().empty()) { parent = dynamic_ptr_cast(mother->children()[0]); } if(!parent) { parent = dynamic_ptr_cast(mother); } partner = parent->partner(); } else { partner = dynamic_ptr_cast(&particle)->partner(); } } } return partner; } pair softPhiMin(double phi0, double phi1, double A, double B, double C, double D) { double c01 = cos(phi0 - phi1); double s01 = sin(phi0 - phi1); double s012(sqr(s01)), c012(sqr(c01)); double A2(A*A), B2(B*B), C2(C*C), D2(D*D); if(abs(B/A)<1e-10 && abs(D/C)<1e-10) return make_pair(phi0,phi0+Constants::pi); double root = sqr(B2)*C2*D2*sqr(s012) + 2.*A*B2*B*C2*C*D*c01*s012 + 2.*A*B2*B*C*D2*D*c01*s012 + 4.*A2*B2*C2*D2*c012 - A2*B2*C2*D2*s012 - A2*B2*sqr(D2)*s012 - sqr(B2)*sqr(C2)*s012 - sqr(B2)*C2*D2*s012 - 4.*A2*A*B*C*D2*D*c01 - 4.*A*B2*B*C2*C*D*c01 + sqr(A2)*sqr(D2) + 2.*A2*B2*C2*D2 + sqr(B2)*sqr(C2); if(root<0.) return make_pair(phi0,phi0+Constants::pi); root = sqrt(root); double denom = (-2.*A*B*C*D*c01 + A2*D2 + B2*C2); double denom2 = (-B*C*c01 + A*D); double num = B2*C*D*s012; return make_pair(atan2(B*s01*(-C*(num + root) / denom + D) / denom2, -(num + root ) / denom) + phi0, atan2(B*s01*(-C*(num - root) / denom + D) / denom2, -(num - root ) / denom) + phi0); } } double QTildeSudakov::generatePhiForward(ShowerParticle & particle, const IdList & ids, - ShoKinPtr kinematics) { + ShoKinPtr kinematics, + const RhoDMatrix & rho) { // no correlations, return flat phi if(! ShowerHandler::currentHandler()->evolver()->correlations()) return Constants::twopi*UseRandom::rnd(); // get the kinematic variables double z = kinematics->z(); Energy2 t = z*(1.-z)*sqr(kinematics->scale()); Energy pT = kinematics->pT(); // if soft correlations Energy2 pipj,pik; bool canBeSoft[2] = {ids[1]==ParticleID::g || ids[1]==ParticleID::gamma, ids[2]==ParticleID::g || ids[2]==ParticleID::gamma }; vector pjk(3,ZERO); vector Ek(3,ZERO); Energy Ei,Ej; Energy2 m12(ZERO),m22(ZERO); InvEnergy2 aziMax(ZERO); bool softAllowed = ShowerHandler::currentHandler()->evolver()->softCorrelations()&& (canBeSoft[0] || canBeSoft[1]); if(softAllowed) { // find the partner for the soft correlations tShowerParticlePtr partner=findCorrelationPartner(particle,true,splittingFn()->interactionType()); // remember we want the softer gluon bool swapOrder = !canBeSoft[1] || (canBeSoft[0] && canBeSoft[1] && z < 0.5); double zFact = !swapOrder ? (1.-z) : z; // compute the transforms to the shower reference frame // first the boost Lorentz5Momentum pVect = particle.showerBasis()->pVector(); Lorentz5Momentum nVect = particle.showerBasis()->nVector(); Boost beta_bb; if(particle.showerBasis()->frame()==ShowerBasis::BackToBack) { beta_bb = -(pVect + nVect).boostVector(); } else if(particle.showerBasis()->frame()==ShowerBasis::Rest) { beta_bb = -pVect.boostVector(); } else assert(false); pVect.boost(beta_bb); nVect.boost(beta_bb); Axis axis; if(particle.showerBasis()->frame()==ShowerBasis::BackToBack) { axis = pVect.vect().unit(); } else if(particle.showerBasis()->frame()==ShowerBasis::Rest) { axis = nVect.vect().unit(); } else assert(false); // and then the rotation LorentzRotation rot; if(axis.perp2()>0.) { double sinth(sqrt(sqr(axis.x())+sqr(axis.y()))); rot.rotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.)); } else if(axis.z()<0.) { rot.rotate(Constants::pi,Axis(1.,0.,0.)); } rot.invert(); pVect *= rot; nVect *= rot; // shower parameters Energy2 pn = pVect*nVect, m2 = pVect.m2(); double alpha0 = particle.showerParameters().alpha; double beta0 = 0.5/alpha0/pn* (sqr(particle.dataPtr()->mass())-sqr(alpha0)*m2+sqr(particle.showerParameters().pt)); Lorentz5Momentum qperp0(particle.showerParameters().ptx, particle.showerParameters().pty,ZERO,ZERO); assert(partner); Lorentz5Momentum pj = partner->momentum(); pj.boost(beta_bb); pj *= rot; // compute the two phi independent dot products pik = 0.5*zFact*(sqr(alpha0)*m2 - sqr(particle.showerParameters().pt) + 2.*alpha0*beta0*pn ) +0.5*sqr(pT)/zFact; Energy2 dot1 = pj*pVect; Energy2 dot2 = pj*nVect; Energy2 dot3 = pj*qperp0; pipj = alpha0*dot1+beta0*dot2+dot3; // compute the constants for the phi dependent dot product pjk[0] = zFact*(alpha0*dot1+dot3-0.5*dot2/pn*(alpha0*m2-sqr(particle.showerParameters().pt)/alpha0)) +0.5*sqr(pT)*dot2/pn/zFact/alpha0; pjk[1] = (pj.x() - dot2/alpha0/pn*qperp0.x())*pT; pjk[2] = (pj.y() - dot2/alpha0/pn*qperp0.y())*pT; m12 = sqr(particle.dataPtr()->mass()); m22 = sqr(partner->dataPtr()->mass()); if(swapOrder) { pjk[1] *= -1.; pjk[2] *= -1.; } Ek[0] = zFact*(alpha0*pVect.t()-0.5*nVect.t()/pn*(alpha0*m2-sqr(particle.showerParameters().pt)/alpha0)) +0.5*sqr(pT)*nVect.t()/pn/zFact/alpha0; Ek[1] = -nVect.t()/alpha0/pn*qperp0.x()*pT; Ek[2] = -nVect.t()/alpha0/pn*qperp0.y()*pT; if(swapOrder) { Ek[1] *= -1.; Ek[2] *= -1.; } Energy mag2=sqrt(sqr(Ek[1])+sqr(Ek[2])); Ei = alpha0*pVect.t()+beta0*nVect.t(); Ej = pj.t(); double phi0 = atan2(-pjk[2],-pjk[1]); if(phi0<0.) phi0 += Constants::twopi; double phi1 = atan2(-Ek[2],-Ek[1]); if(phi1<0.) phi1 += Constants::twopi; double xi_min = pik/Ei/(Ek[0]+mag2), xi_max = pik/Ei/(Ek[0]-mag2), xi_ij = pipj/Ei/Ej; if(xi_min>xi_max) swap(xi_min,xi_max); if(xi_min>xi_ij) softAllowed = false; Energy2 mag = sqrt(sqr(pjk[1])+sqr(pjk[2])); if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==1) { aziMax = -m12/sqr(pik) -m22/sqr(pjk[0]+mag) +2.*pipj/pik/(pjk[0]-mag); } else if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==2) { double A = (pipj*Ek[0]- Ej*pik)/Ej/sqr(Ej); double B = -sqrt(sqr(pipj)*(sqr(Ek[1])+sqr(Ek[2])))/Ej/sqr(Ej); double C = pjk[0]/sqr(Ej); double D = -sqrt(sqr(pjk[1])+sqr(pjk[2]))/sqr(Ej); pair minima = softPhiMin(phi0,phi1,A,B,C,D); aziMax = 0.5/pik/(Ek[0]-mag2)*(Ei-m12*(Ek[0]-mag2)/pik + max(Ej*(A+B*cos(minima.first -phi1))/(C+D*cos(minima.first -phi0)), Ej*(A+B*cos(minima.second-phi1))/(C+D*cos(minima.second-phi0)))); } else assert(false); } // if spin correlations vector > wgts; if(ShowerHandler::currentHandler()->evolver()->spinCorrelations()) { - RhoDMatrix rho = particle.extractRhoMatrix(true); // calculate the weights wgts = splittingFn()->generatePhiForward(z,t,ids,rho); } else { wgts = vector >(1,make_pair(0,1.)); } // generate the azimuthal angle double phi,wgt; static const Complex ii(0.,1.); unsigned int ntry(0); double phiMax(0.),wgtMax(0.); do { phi = Constants::twopi*UseRandom::rnd(); // first the spin correlations bit (gives 1 if correlations off) Complex spinWgt = 0.; for(unsigned int ix=0;ix1e-10) { generator()->log() << "Forward spin weight problem " << wgt << " " << wgt-1. << " " << ids[0] << " " << ids[1] << " " << ids[2] << " " << " " << phi << "\n"; generator()->log() << "Weights \n"; for(unsigned int ix=0;ixlog() << wgts[ix].first << " " << wgts[ix].second << "\n"; } // soft correlations bit double aziWgt = 1.; if(softAllowed) { Energy2 dot = pjk[0]+pjk[1]*cos(phi)+pjk[2]*sin(phi); Energy Eg = Ek[0]+Ek[1]*cos(phi)+Ek[2]*sin(phi); if(pipj*Eg>pik*Ej) { if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==1) { aziWgt = (-m12/sqr(pik) -m22/sqr(dot) +2.*pipj/pik/dot)/aziMax; } else if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==2) { aziWgt = max(ZERO,0.5/pik/Eg*(Ei-m12*Eg/pik + (pipj*Eg - Ej*pik)/dot)/aziMax); } if(aziWgt-1.>1e-10||aziWgt<-1e-10) { generator()->log() << "Forward soft weight problem " << aziWgt << " " << aziWgt-1. << " " << ids[0] << " " << ids[1] << " " << ids[2] << " " << " " << phi << "\n"; } } else { aziWgt = 0.; } } wgt *= aziWgt; if(wgt>wgtMax) { phiMax = phi; wgtMax = wgt; } ++ntry; } while(wgtlog() << "Too many tries to generate phi in forward evolution\n"; phi = phiMax; } // return the azimuthal angle return phi; } double QTildeSudakov::generatePhiBackward(ShowerParticle & particle, const IdList & ids, - ShoKinPtr kinematics) { + ShoKinPtr kinematics, + const RhoDMatrix & rho) { // no correlations, return flat phi if(! ShowerHandler::currentHandler()->evolver()->correlations()) return Constants::twopi*UseRandom::rnd(); // get the kinematic variables double z = kinematics->z(); Energy2 t = (1.-z)*sqr(kinematics->scale())/z; Energy pT = kinematics->pT(); // if soft correlations bool softAllowed = ShowerHandler::currentHandler()->evolver()->softCorrelations() && (ids[2]==ParticleID::g || ids[2]==ParticleID::gamma); Energy2 pipj,pik,m12(ZERO),m22(ZERO); vector pjk(3,ZERO); Energy Ei,Ej,Ek; InvEnergy2 aziMax(ZERO); if(softAllowed) { // find the partner for the soft correlations tShowerParticlePtr partner=findCorrelationPartner(particle,false,splittingFn()->interactionType()); double zFact = (1.-z); // compute the transforms to the shower reference frame // first the boost Lorentz5Momentum pVect = particle.showerBasis()->pVector(); Lorentz5Momentum nVect = particle.showerBasis()->nVector(); assert(particle.showerBasis()->frame()==ShowerBasis::BackToBack); Boost beta_bb = -(pVect + nVect).boostVector(); pVect.boost(beta_bb); nVect.boost(beta_bb); Axis axis = pVect.vect().unit(); // and then the rotation LorentzRotation rot; if(axis.perp2()>0.) { double sinth(sqrt(sqr(axis.x())+sqr(axis.y()))); rot.rotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.)); } else if(axis.z()<0.) { rot.rotate(Constants::pi,Axis(1.,0.,0.)); } rot.invert(); pVect *= rot; nVect *= rot; // shower parameters Energy2 pn = pVect*nVect; Energy2 m2 = pVect.m2(); double alpha0 = particle.x(); double beta0 = -0.5/alpha0/pn*sqr(alpha0)*m2; Lorentz5Momentum pj = partner->momentum(); pj.boost(beta_bb); pj *= rot; double beta2 = 0.5*(1.-zFact)*(sqr(alpha0*zFact/(1.-zFact))*m2+sqr(pT))/alpha0/zFact/pn; // compute the two phi independent dot products Energy2 dot1 = pj*pVect; Energy2 dot2 = pj*nVect; pipj = alpha0*dot1+beta0*dot2; pik = alpha0*(alpha0*zFact/(1.-zFact)*m2+pn*(beta2+zFact/(1.-zFact)*beta0)); // compute the constants for the phi dependent dot product pjk[0] = alpha0*zFact/(1.-zFact)*dot1+beta2*dot2; pjk[1] = pj.x()*pT; pjk[2] = pj.y()*pT; m12 = ZERO; m22 = sqr(partner->dataPtr()->mass()); Energy2 mag = sqrt(sqr(pjk[1])+sqr(pjk[2])); if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==1) { aziMax = -m12/sqr(pik) -m22/sqr(pjk[0]+mag) +2.*pipj/pik/(pjk[0]-mag); } else if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==2) { Ek = alpha0*zFact/(1.-zFact)*pVect.t()+beta2*nVect.t(); Ei = alpha0*pVect.t()+beta0*nVect.t(); Ej = pj.t(); if(pipj*Ek> Ej*pik) { aziMax = 0.5/pik/Ek*(Ei-m12*Ek/pik + (pipj*Ek- Ej*pik)/(pjk[0]-mag)); } else { aziMax = 0.5/pik/Ek*(Ei-m12*Ek/pik); } } else { assert(ShowerHandler::currentHandler()->evolver()->softCorrelations()==0); } } // if spin correlations vector > wgts; if(ShowerHandler::currentHandler()->evolver()->spinCorrelations()) { // get the spin density matrix and the mapping - RhoDMatrix rho = particle.extractRhoMatrix(false); // get the weights wgts = splittingFn()->generatePhiBackward(z,t,ids,rho); } else { wgts = vector >(1,make_pair(0,1.)); } // generate the azimuthal angle double phi,wgt; static const Complex ii(0.,1.); unsigned int ntry(0); double phiMax(0.),wgtMax(0.); do { phi = Constants::twopi*UseRandom::rnd(); Complex spinWgt = 0.; for(unsigned int ix=0;ix1e-10) { generator()->log() << "Backward weight problem " << wgt << " " << wgt-1. << " " << ids[0] << " " << ids[1] << " " << ids[2] << " " << " " << z << " " << phi << "\n"; generator()->log() << "Weights \n"; for(unsigned int ix=0;ixlog() << wgts[ix].first << " " << wgts[ix].second << "\n"; } // soft correlations bit double aziWgt = 1.; if(softAllowed) { Energy2 dot = pjk[0]+pjk[1]*cos(phi)+pjk[2]*sin(phi); if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==1) { aziWgt = (-m12/sqr(pik) -m22/sqr(dot) +2.*pipj/pik/dot)/aziMax; } else if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==2) { aziWgt = max(ZERO,0.5/pik/Ek*(Ei-m12*Ek/pik + pipj*Ek/dot - Ej*pik/dot)/aziMax); } if(aziWgt-1.>1e-10||aziWgt<-1e-10) { generator()->log() << "Backward soft weight problem " << aziWgt << " " << aziWgt-1. << " " << ids[0] << " " << ids[1] << " " << ids[2] << " " << " " << phi << "\n"; } } wgt *= aziWgt; if(wgt>wgtMax) { phiMax = phi; wgtMax = wgt; } ++ntry; } while(wgtlog() << "Too many tries to generate phi in backward evolution\n"; phi = phiMax; } // return the azimuthal angle return phi; } double QTildeSudakov::generatePhiDecay(ShowerParticle & particle, const IdList & ids, - ShoKinPtr kinematics) { + ShoKinPtr kinematics, + const RhoDMatrix &) { // only soft correlations in this case // no correlations, return flat phi if( !(ShowerHandler::currentHandler()->evolver()->softCorrelations() && (ids[2]==ParticleID::g || ids[2]==ParticleID::gamma ))) return Constants::twopi*UseRandom::rnd(); // get the kinematic variables double z = kinematics->z(); Energy pT = kinematics->pT(); // if soft correlations // find the partner for the soft correlations tShowerParticlePtr partner = findCorrelationPartner(particle,true,splittingFn()->interactionType()); double zFact(1.-z); // compute the transforms to the shower reference frame // first the boost Lorentz5Momentum pVect = particle.showerBasis()->pVector(); Lorentz5Momentum nVect = particle.showerBasis()->nVector(); assert(particle.showerBasis()->frame()==ShowerBasis::Rest); Boost beta_bb = -pVect.boostVector(); pVect.boost(beta_bb); nVect.boost(beta_bb); Axis axis = nVect.vect().unit(); // and then the rotation LorentzRotation rot; if(axis.perp2()>0.) { double sinth(sqrt(sqr(axis.x())+sqr(axis.y()))); rot.rotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.)); } else if(axis.z()<0.) { rot.rotate(Constants::pi,Axis(1.,0.,0.)); } rot.invert(); pVect *= rot; nVect *= rot; // shower parameters Energy2 pn = pVect*nVect; Energy2 m2 = pVect.m2(); double alpha0 = particle.showerParameters().alpha; double beta0 = 0.5/alpha0/pn* (sqr(particle.dataPtr()->mass())-sqr(alpha0)*m2+sqr(particle.showerParameters().pt)); Lorentz5Momentum qperp0(particle.showerParameters().ptx, particle.showerParameters().pty,ZERO,ZERO); Lorentz5Momentum pj = partner->momentum(); pj.boost(beta_bb); pj *= rot; // compute the two phi independent dot products Energy2 pik = 0.5*zFact*(sqr(alpha0)*m2 - sqr(particle.showerParameters().pt) + 2.*alpha0*beta0*pn ) +0.5*sqr(pT)/zFact; Energy2 dot1 = pj*pVect; Energy2 dot2 = pj*nVect; Energy2 dot3 = pj*qperp0; Energy2 pipj = alpha0*dot1+beta0*dot2+dot3; // compute the constants for the phi dependent dot product vector pjk(3,ZERO); pjk[0] = zFact*(alpha0*dot1+dot3-0.5*dot2/pn*(alpha0*m2-sqr(particle.showerParameters().pt)/alpha0)) +0.5*sqr(pT)*dot2/pn/zFact/alpha0; pjk[1] = (pj.x() - dot2/alpha0/pn*qperp0.x())*pT; pjk[2] = (pj.y() - dot2/alpha0/pn*qperp0.y())*pT; Energy2 m12 = sqr(particle.dataPtr()->mass()); Energy2 m22 = sqr(partner->dataPtr()->mass()); Energy2 mag = sqrt(sqr(pjk[1])+sqr(pjk[2])); InvEnergy2 aziMax; vector Ek(3,ZERO); Energy Ei,Ej; if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==1) { aziMax = -m12/sqr(pik) -m22/sqr(pjk[0]+mag) +2.*pipj/pik/(pjk[0]-mag); } else if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==2) { Ek[0] = zFact*(alpha0*pVect.t()+-0.5*nVect.t()/pn*(alpha0*m2-sqr(particle.showerParameters().pt)/alpha0)) +0.5*sqr(pT)*nVect.t()/pn/zFact/alpha0; Ek[1] = -nVect.t()/alpha0/pn*qperp0.x()*pT; Ek[2] = -nVect.t()/alpha0/pn*qperp0.y()*pT; Energy mag2=sqrt(sqr(Ek[1])+sqr(Ek[2])); Ei = alpha0*pVect.t()+beta0*nVect.t(); Ej = pj.t(); aziMax = 0.5/pik/(Ek[0]-mag2)*(Ei-m12*(Ek[0]-mag2)/pik + pipj*(Ek[0]+mag2)/(pjk[0]-mag) - Ej*pik/(pjk[0]-mag) ); } else assert(ShowerHandler::currentHandler()->evolver()->softCorrelations()==0); // generate the azimuthal angle double phi,wgt(0.); unsigned int ntry(0); double phiMax(0.),wgtMax(0.); do { phi = Constants::twopi*UseRandom::rnd(); Energy2 dot = pjk[0]+pjk[1]*cos(phi)+pjk[2]*sin(phi); if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==1) { wgt = (-m12/sqr(pik) -m22/sqr(dot) +2.*pipj/pik/dot)/aziMax; } else if(ShowerHandler::currentHandler()->evolver()->softCorrelations()==2) { if(qperp0.m2()==ZERO) { wgt = 1.; } else { Energy Eg = Ek[0]+Ek[1]*cos(phi)+Ek[2]*sin(phi); wgt = max(ZERO,0.5/pik/Eg*(Ei-m12*Eg/pik + (pipj*Eg - Ej*pik)/dot)/aziMax); } } if(wgt-1.>1e-10||wgt<-1e-10) { generator()->log() << "Decay soft weight problem " << wgt << " " << wgt-1. << " " << ids[0] << " " << ids[1] << " " << ids[2] << " " << " " << phi << "\n"; } if(wgt>wgtMax) { phiMax = phi; wgtMax = wgt; } ++ntry; } while(wgtlog() << "Too many tries to generate phi\n"; } // return the azimuthal angle return phi; } Energy QTildeSudakov::calculateScale(double zin, Energy pt, IdList ids, unsigned int iopt) { Energy2 tmin; initialize(ids,tmin,false); // final-state branching if(iopt==0) { Energy2 scale=(sqr(pt)+masssquared_[1]*(1.-zin)+masssquared_[2]*zin); if(ids[0]!=ParticleID::g) scale -= zin*(1.-zin)*masssquared_[0]; scale /= sqr(zin*(1-zin)); return scale<=ZERO ? sqrt(tmin) : sqrt(scale); } else if(iopt==1) { Energy2 scale=(sqr(pt)+zin*masssquared_[2])/sqr(1.-zin); return scale<=ZERO ? sqrt(tmin) : sqrt(scale); } else if(iopt==2) { Energy2 scale = (sqr(pt)+zin*masssquared_[2])/sqr(1.-zin)+masssquared_[0]; return scale<=ZERO ? sqrt(tmin) : sqrt(scale); } else { throw Exception() << "Unknown option in QTildeSudakov::calculateScale() " << "iopt = " << iopt << Exception::runerror; } } ShoKinPtr QTildeSudakov::createFinalStateBranching(Energy scale,double z, double phi, Energy pt) { ShoKinPtr showerKin = new_ptr(FS_QTildeShowerKinematics1to2()); showerKin->scale(scale); showerKin->z(z); showerKin->phi(phi); showerKin->pT(pt); showerKin->SudakovFormFactor(this); return showerKin; } ShoKinPtr QTildeSudakov::createInitialStateBranching(Energy scale,double z, double phi, Energy pt) { ShoKinPtr showerKin = new_ptr(IS_QTildeShowerKinematics1to2()); showerKin->scale(scale); showerKin->z(z); showerKin->phi(phi); showerKin->pT(pt); showerKin->SudakovFormFactor(this); return showerKin; } ShoKinPtr QTildeSudakov::createDecayBranching(Energy scale,double z, double phi, Energy pt) { ShoKinPtr showerKin = new_ptr(Decay_QTildeShowerKinematics1to2()); showerKin->scale(scale); showerKin->z(z); showerKin->phi(phi); showerKin->pT(pt); showerKin->SudakovFormFactor(this); return showerKin; } diff --git a/Shower/Default/QTildeSudakov.h b/Shower/Default/QTildeSudakov.h --- a/Shower/Default/QTildeSudakov.h +++ b/Shower/Default/QTildeSudakov.h @@ -1,287 +1,294 @@ // -*- C++ -*- // // QTildeSudakov.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_QTildeSudakov_H #define HERWIG_QTildeSudakov_H // // This is the declaration of the QTildeSudakov class. // #include "Herwig/Shower/Base/SudakovFormFactor.h" namespace Herwig { using namespace ThePEG; /** \ingroup Shower * * The QTildeSudakov class implements the Sudakov form factor for evolution in * \f$\tilde{q}^2\f$ using the veto algorithm. * * @see \ref QTildeSudakovInterfaces "The interfaces" * defined for QTildeSudakov. */ class QTildeSudakov: public SudakovFormFactor { public: /** * The default constructor. */ inline QTildeSudakov() {} /** * Members to generate the scale of the next branching */ //@{ /** * Return the scale of the next time-like branching. If there is no * branching then it returns ZERO. * @param startingScale starting scale for the evolution * @param ids The PDG codes of the particles in the splitting * @param cc Whether this is the charge conjugate of the branching * defined. * @param enhance The radiation enhancement factor * @param maxQ2 The maximum \f$Q^2\f$ for the emission */ virtual ShoKinPtr generateNextTimeBranching(const Energy startingScale, const IdList &ids,const bool cc, + const RhoDMatrix & rho, double enhance, Energy2 maxQ2); /** * Return the scale of the next space-like decay branching. If there is no * branching then it returns ZERO. * @param startingScale starting scale for the evolution * @param stoppingScale stopping scale for the evolution * @param minmass The minimum mass allowed for the spake-like particle. * @param ids The PDG codes of the particles in the splitting * @param cc Whether this is the charge conjugate of the branching * defined. * @param enhance The radiation enhancement factor */ virtual ShoKinPtr generateNextDecayBranching(const Energy startingScale, - const Energy stoppingScale, - const Energy minmass, - const IdList &ids, - const bool cc, - double enhance); + const Energy stoppingScale, + const Energy minmass, + const IdList &ids, + const bool cc, + const RhoDMatrix & rho, + double enhance); /** * Return the scale of the next space-like branching. If there is no * branching then it returns ZERO. * @param startingScale starting scale for the evolution * @param ids The PDG codes of the particles in the splitting * @param x The fraction of the beam momentum * @param cc Whether this is the charge conjugate of the branching * defined. * @param enhance The radiation enhancement factor * @param beam The beam particle */ virtual ShoKinPtr generateNextSpaceBranching(const Energy startingScale, const IdList &ids,double x, - const bool cc, double enhance, + const bool cc, + const RhoDMatrix & rho, + double enhance, tcBeamPtr beam); //@} /** * Generate the azimuthal angle of the branching for forward branching * @param particle The branching particle * @param ids The PDG codes of the particles in the branchings * @param The Shower kinematics */ virtual double generatePhiForward(ShowerParticle & particle,const IdList & ids, - ShoKinPtr kinematics); + ShoKinPtr kinematics, + const RhoDMatrix & rho); /** * Generate the azimuthal angle of the branching for backward branching * @param particle The branching particle * @param ids The PDG codes of the particles in the branchings * @param The Shower kinematics */ virtual double generatePhiBackward(ShowerParticle & particle,const IdList & ids, - ShoKinPtr kinematics); + ShoKinPtr kinematics, + const RhoDMatrix & rho); /** * Generate the azimuthal angle of the branching for ISR in decays * @param particle The branching particle * @param ids The PDG codes of the particles in the branchings * @param The Shower kinematics */ virtual double generatePhiDecay(ShowerParticle & particle,const IdList & ids, - ShoKinPtr kinematics); + ShoKinPtr kinematics, + const RhoDMatrix & rho); /** * Method to return the evolution scale given the * transverse momentum, \f$p_T\f$ and \f$z\f$. */ virtual Energy calculateScale(double z, Energy pt, IdList ids,unsigned int iopt); /** * Method to create the ShowerKinematics object for a final-state branching */ virtual ShoKinPtr createFinalStateBranching(Energy scale,double z, double phi, Energy pt); /** * Method to create the ShowerKinematics object for an initial-state branching */ virtual ShoKinPtr createInitialStateBranching(Energy scale,double z, double phi, Energy pt); /** * Method to create the ShowerKinematics object for a decay branching */ virtual ShoKinPtr createDecayBranching(Energy scale,double z, double phi, Energy pt); 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: /** * Methods to provide the next value of the scale before the vetos * are applied. */ //@{ /** * Value of the energy fraction and scale for time-like branching * @param t The scale * @param tmin The minimum scale * @param enhance The radiation enhancement factor * @return False if scale less than minimum, true otherwise */ bool guessTimeLike(Energy2 &t, Energy2 tmin, double enhance); /** * Value of the energy fraction and scale for time-like branching * @param t The scale * @param tmax The maximum scale * @param minmass The minimum mass of the particle after the branching * @param enhance The radiation enhancement factor */ bool guessDecay(Energy2 &t, Energy2 tmax,Energy minmass, double enhance); /** * Value of the energy fraction and scale for space-like branching * @param t The scale * @param tmin The minimum scale * @param x Fraction of the beam momentum. * @param enhance The radiation enhancement factor */ bool guessSpaceLike(Energy2 &t, Energy2 tmin, const double x, double enhance); //@} /** * Initialize the values of the cut-offs and scales * @param tmin The minimum scale * @param ids The ids of the partics in the branching * @param cc Whether this is the charge conjugate of the branching */ void initialize(const IdList & ids,Energy2 &tmin, const bool cc); /** * Phase Space veto member to implement the \f$\Theta\f$ function as a veto * so that the emission is within the allowed phase space. * @param t The scale * @param maxQ2 The maximum virtuality * @return true if vetoed */ bool PSVeto(const Energy2 t,const Energy2 maxQ2); /** * Compute the limits on \f$z\f$ for time-like branching * @param scale The scale of the particle * @return True if lower limit less than upper, otherwise false */ bool computeTimeLikeLimits(Energy2 & scale); /** * Compute the limits on \f$z\f$ for space-like branching * @param scale The scale of the particle * @param x The energy fraction of the parton * @return True if lower limit less than upper, otherwise false */ bool computeSpaceLikeLimits(Energy2 & scale, double x); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ inline virtual IBPtr clone() const {return new_ptr(*this);} /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ inline virtual IBPtr fullclone() const {return new_ptr(*this);} //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ QTildeSudakov & operator=(const QTildeSudakov &); private: /** * The evolution scale, \f$\tilde{q}\f$. */ Energy q_; /** * The Ids of the particles in the current branching */ IdList ids_; /** * The masses of the particles in the current branching */ vector masses_; /** * The mass squared of the particles in the current branching */ vector masssquared_; }; } #endif /* HERWIG_QTildeSudakov_H */ diff --git a/Shower/SplittingFunctions/HalfHalfOneEWSplitFn.cc b/Shower/SplittingFunctions/HalfHalfOneEWSplitFn.cc --- a/Shower/SplittingFunctions/HalfHalfOneEWSplitFn.cc +++ b/Shower/SplittingFunctions/HalfHalfOneEWSplitFn.cc @@ -1,214 +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() { SplittingFunction::doinit(); 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 { + const IdList &ids, const bool mass, const RhoDMatrix & rho) 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)); + val *= (sqr(gL)*norm(rho(0,0))+sqr(gR)*norm(rho(1,1))); 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 { + const IdList & ids, const bool mass, const RhoDMatrix & rho) 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)); + val *= (sqr(gL)*norm(rho(0,0))+sqr(gR)*norm(rho(1,1))); 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/HalfHalfOneEWSplitFn.h b/Shower/SplittingFunctions/HalfHalfOneEWSplitFn.h --- a/Shower/SplittingFunctions/HalfHalfOneEWSplitFn.h +++ b/Shower/SplittingFunctions/HalfHalfOneEWSplitFn.h @@ -1,215 +1,217 @@ // -*- C++ -*- #ifndef Herwig_HalfHalfOneEWSplitFn_H #define Herwig_HalfHalfOneEWSplitFn_H // // This is the declaration of the HalfHalfOneEWSplitFn class. // #include "Herwig/Shower/SplittingFunctions/SplittingFunction.h" namespace Herwig { using namespace ThePEG; /** * The HalfHalfOneEWSplitFn class implements the splitting function for * \f$\frac12\to q\frac12 1\f$ where the spin-1 particle is a massive electroweak gauge boson. * * @see \ref HalfHalfOneEWSplitFnInterfaces "The interfaces" * defined for HalfHalfOneEWSplitFn. */ class HalfHalfOneEWSplitFn: public SplittingFunction { public: /** * The default constructor. */ HalfHalfOneEWSplitFn() : SplittingFunction(1) {} /** * Concrete implementation of the method to determine whether this splitting * function can be used for a given set of particles. * @param ids The PDG codes for the particles in the splitting. */ virtual bool accept(const IdList & ids) const; /** * Methods to return the splitting function. */ //@{ /** * The concrete implementation of the splitting function, \f$P(z,t)\f$. * @param z The energy fraction. * @param t The scale. * @param ids The PDG codes for the particles in the splitting. * @param mass Whether or not to include the mass dependent terms + * @param rho The spin density matrix */ virtual double P(const double z, const Energy2 t, const IdList & ids, - const bool mass) const; + const bool mass, const RhoDMatrix & rho) const; /** * The concrete implementation of the overestimate of the splitting function, * \f$P_{\rm over}\f$. * @param z The energy fraction. * @param ids The PDG codes for the particles in the splitting. */ virtual double overestimateP(const double z, const IdList & ids) const; /** * The concrete implementation of the * the ratio of the splitting function to the overestimate, i.e. * \f$P(z,t)/P_{\rm over}(z)\f$. * @param z The energy fraction. * @param t The scale. * @param ids The PDG codes for the particles in the splitting. * @param mass Whether or not to include the mass dependent terms + * @param rho The spin density matrix */ virtual double ratioP(const double z, const Energy2 t, const IdList & ids, - const bool mass) const; + const bool mass, const RhoDMatrix & rho) const; /** * The concrete implementation of the indefinite integral of the * overestimated splitting function, \f$P_{\rm over}\f$. * @param z The energy fraction. * @param ids The PDG codes for the particles in the splitting. * @param PDFfactor Which additional factor to include for the PDF * 0 is no additional factor, * 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$ */ virtual double integOverP(const double z, const IdList & ids, unsigned int PDFfactor=0) const; /** * The concrete implementation of the inverse of the indefinite integral. * @param r Value of the splitting function to be inverted * @param ids The PDG codes for the particles in the splitting. * @param PDFfactor Which additional factor to include for the PDF * 0 is no additional factor, * 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$ */ virtual double invIntegOverP(const double r, const IdList & ids, unsigned int PDFfactor=0) const; //@} /** * Method to calculate the azimuthal angle * @param z The energy fraction * @param t The scale \f$t=2p_j\cdot p_k\f$. * @param ids The PDG codes for the particles in the splitting. * @param The azimuthal angle, \f$\phi\f$. * @return The weight */ virtual vector > generatePhiForward(const double z, const Energy2 t, const IdList & ids, const RhoDMatrix &); /** * Method to calculate the azimuthal angle for backward evolution * @param z The energy fraction * @param t The scale \f$t=2p_j\cdot p_k\f$. * @param ids The PDG codes for the particles in the splitting. * @param The azimuthal angle, \f$\phi\f$. * @return The weight */ virtual vector > generatePhiBackward(const double z, const Energy2 t, const IdList & ids, const RhoDMatrix &); /** * Calculate the matrix element for the splitting * @param z The energy fraction * @param t The scale \f$t=2p_j\cdot p_k\f$. * @param ids The PDG codes for the particles in the splitting. * @param The azimuthal angle, \f$\phi\f$. */ virtual DecayMEPtr matrixElement(const double z, const Energy2 t, const IdList & ids, const double phi, bool timeLike); protected: /** * Get the couplings */ void getCouplings(double & gL, double & gR, const IdList & ids) const; public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ HalfHalfOneEWSplitFn & operator=(const HalfHalfOneEWSplitFn &); private: /** * Z couplings */ map > gZ_; /** * W couplings */ double gWL_; }; } #endif /* Herwig_HalfHalfOneEWSplitFn_H */ diff --git a/Shower/SplittingFunctions/HalfHalfOneSplitFn.cc b/Shower/SplittingFunctions/HalfHalfOneSplitFn.cc --- a/Shower/SplittingFunctions/HalfHalfOneSplitFn.cc +++ b/Shower/SplittingFunctions/HalfHalfOneSplitFn.cc @@ -1,136 +1,136 @@ // -*- C++ -*- // // HalfHalfOneSplitFn.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 HalfHalfOneSplitFn class. // #include "HalfHalfOneSplitFn.h" #include "ThePEG/PDT/ParticleData.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Utilities/DescribeClass.h" #include "Herwig/Decay/TwoBodyDecayMatrixElement.h" using namespace Herwig; DescribeNoPIOClass describeHalfHalfOneSplitFn ("Herwig::HalfHalfOneSplitFn","HwShower.so"); void HalfHalfOneSplitFn::Init() { static ClassDocumentation documentation ("The HalfHalfOneSplitFn class implements the q -> qg splitting function"); } double HalfHalfOneSplitFn::P(const double z, const Energy2 t, - const IdList &ids, const bool mass) const { + const IdList &ids, const bool mass, const RhoDMatrix &) const { double val = (1. + sqr(z))/(1.-z); if(mass) { Energy m = getParticleData(ids[0])->mass(); val -= 2.*sqr(m)/t; } return colourFactor(ids)*val; } double HalfHalfOneSplitFn::overestimateP(const double z, const IdList & ids) const { return 2.*colourFactor(ids)/(1.-z); } double HalfHalfOneSplitFn::ratioP(const double z, const Energy2 t, - const IdList & ids, const bool mass) const { + const IdList & ids, const bool mass, const RhoDMatrix &) const { double val = 1. + sqr(z); if(mass) { Energy m = getParticleData(ids[0])->mass(); val -= 2.*sqr(m)*(1.-z)/t; } return 0.5*val; } double HalfHalfOneSplitFn::integOverP(const double z, const IdList & ids, unsigned int PDFfactor) const { switch (PDFfactor) { case 0: return -2.*colourFactor(ids)*Math::log1m(z); case 1: return 2.*colourFactor(ids)*log(z/(1.-z)); case 2: return 2.*colourFactor(ids)/(1.-z); case 3: default: throw Exception() << "HalfHalfOneSplitFn::integOverP() invalid PDFfactor = " << PDFfactor << Exception::runerror; } } double HalfHalfOneSplitFn::invIntegOverP(const double r, const IdList & ids, unsigned int PDFfactor) const { switch (PDFfactor) { case 0: return 1. - exp(- 0.5*r/colourFactor(ids)); case 1: return 1./(1.-exp(-0.5*r/colourFactor(ids))); case 2: return 1.-2.*colourFactor(ids)/r; case 3: default: throw Exception() << "HalfHalfOneSplitFn::invIntegOverP() invalid PDFfactor = " << PDFfactor << Exception::runerror; } } bool HalfHalfOneSplitFn::accept(const IdList &ids) const { // 3 particles and in and out fermion same if(ids.size()!=3 || ids[0]!=ids[1]) return false; tcPDPtr q=getParticleData(ids[0]); tcPDPtr g=getParticleData(ids[2]); if(q->iSpin()!=PDT::Spin1Half || g->iSpin()!=PDT::Spin1) return false; return checkColours(ids); } vector > HalfHalfOneSplitFn::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 > HalfHalfOneSplitFn::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 HalfHalfOneSplitFn::matrixElement(const double z, const Energy2 t, const IdList & ids, const double phi, bool timeLike) { // calculate the kernal DecayMEPtr kernal(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin1Half,PDT::Spin1Half,PDT::Spin1))); Energy m = !timeLike ? ZERO : getParticleData(ids[0])->mass(); double mt = m/sqrt(t); double root = sqrt(1.-(1.-z)*sqr(m)/z/t); double romz = sqrt(1.-z); double rz = sqrt(z); Complex phase = exp(Complex(0.,1.)*phi); (*kernal)(0,0,0) = -root/romz*phase; (*kernal)(1,1,2) = -conj((*kernal)(0,0,0)); (*kernal)(0,0,2) = root/romz*z/phase; (*kernal)(1,1,0) = -conj((*kernal)(0,0,2)); (*kernal)(1,0,2) = mt*(1.-z)/rz; (*kernal)(0,1,0) = conj((*kernal)(1,0,2)); (*kernal)(0,1,2) = 0.; (*kernal)(1,0,0) = 0.; return kernal; } diff --git a/Shower/SplittingFunctions/HalfHalfOneSplitFn.h b/Shower/SplittingFunctions/HalfHalfOneSplitFn.h --- a/Shower/SplittingFunctions/HalfHalfOneSplitFn.h +++ b/Shower/SplittingFunctions/HalfHalfOneSplitFn.h @@ -1,186 +1,188 @@ // -*- C++ -*- // // HalfHalfOneSplitFn.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_HalfHalfOneSplitFn_H #define HERWIG_HalfHalfOneSplitFn_H // // This is the declaration of the HalfHalfOneSplitFn class. // #include "SplittingFunction.h" namespace Herwig { using namespace ThePEG; /**\ingroup Shower * * This class provides the concrete implementation of the exact leading-order * splitting function for \f$\frac12\to q\frac12 1\f$. * * In this case the splitting function is given by * \f[P(z,t) =C\left(\frac{1+z^2}{1-z}-2\frac{m^2_q}{t}\right),\f] * where \f$C\f$ is the corresponding colour factor. * Our choice for the overestimate is * \f[P_{\rm over}(z) = \frac{2C}{1-z},\f] * therefore the integral is * \f[\int P_{\rm over}(z) {\rm d}z = -2C\ln(1-z),\f] * and its inverse is * \f[1-\exp\left(\frac{r}{2C}\right).\f] * * @see \ref HalfHalfOneSplitFnInterfaces "The interfaces" * defined for HalfHalfOneSplitFn. */ class HalfHalfOneSplitFn: public SplittingFunction { public: /** * The default constructor. */ HalfHalfOneSplitFn() : SplittingFunction(1) {} /** * Concrete implementation of the method to determine whether this splitting * function can be used for a given set of particles. * @param ids The PDG codes for the particles in the splitting. */ virtual bool accept(const IdList & ids) const; /** * Methods to return the splitting function. */ //@{ /** * The concrete implementation of the splitting function, \f$P(z,t)\f$. * @param z The energy fraction. * @param t The scale. * @param ids The PDG codes for the particles in the splitting. * @param mass Whether or not to include the mass dependent terms + * @param rho The spin density matrix */ virtual double P(const double z, const Energy2 t, const IdList & ids, - const bool mass) const; + const bool mass, const RhoDMatrix & rho) const; /** * The concrete implementation of the overestimate of the splitting function, * \f$P_{\rm over}\f$. * @param z The energy fraction. * @param ids The PDG codes for the particles in the splitting. */ virtual double overestimateP(const double z, const IdList & ids) const; /** * The concrete implementation of the * the ratio of the splitting function to the overestimate, i.e. * \f$P(z,t)/P_{\rm over}(z)\f$. * @param z The energy fraction. * @param t The scale. * @param ids The PDG codes for the particles in the splitting. * @param mass Whether or not to include the mass dependent terms + * @param rho The spin density matrix */ virtual double ratioP(const double z, const Energy2 t, const IdList & ids, - const bool mass) const; + const bool mass, const RhoDMatrix & rho) const; /** * The concrete implementation of the indefinite integral of the * overestimated splitting function, \f$P_{\rm over}\f$. * @param z The energy fraction. * @param ids The PDG codes for the particles in the splitting. * @param PDFfactor Which additional factor to include for the PDF * 0 is no additional factor, * 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$ */ virtual double integOverP(const double z, const IdList & ids, unsigned int PDFfactor=0) const; /** * The concrete implementation of the inverse of the indefinite integral. * @param r Value of the splitting function to be inverted * @param ids The PDG codes for the particles in the splitting. * @param PDFfactor Which additional factor to include for the PDF * 0 is no additional factor, * 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$ */ virtual double invIntegOverP(const double r, const IdList & ids, unsigned int PDFfactor=0) const; //@} /** * Method to calculate the azimuthal angle * @param z The energy fraction * @param t The scale \f$t=2p_j\cdot p_k\f$. * @param ids The PDG codes for the particles in the splitting. * @param The azimuthal angle, \f$\phi\f$. * @return The weight */ virtual vector > generatePhiForward(const double z, const Energy2 t, const IdList & ids, const RhoDMatrix &); /** * Method to calculate the azimuthal angle for backward evolution * @param z The energy fraction * @param t The scale \f$t=2p_j\cdot p_k\f$. * @param ids The PDG codes for the particles in the splitting. * @param The azimuthal angle, \f$\phi\f$. * @return The weight */ virtual vector > generatePhiBackward(const double z, const Energy2 t, const IdList & ids, const RhoDMatrix &); /** * Calculate the matrix element for the splitting * @param z The energy fraction * @param t The scale \f$t=2p_j\cdot p_k\f$. * @param ids The PDG codes for the particles in the splitting. * @param The azimuthal angle, \f$\phi\f$. */ virtual DecayMEPtr matrixElement(const double z, const Energy2 t, const IdList & ids, const double phi, bool timeLike); public: /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const {return new_ptr(*this);} /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const {return new_ptr(*this);} //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ HalfHalfOneSplitFn & operator=(const HalfHalfOneSplitFn &); }; } #endif /* HERWIG_HalfHalfOneSplitFn_H */ diff --git a/Shower/SplittingFunctions/HalfOneHalfSplitFn.cc b/Shower/SplittingFunctions/HalfOneHalfSplitFn.cc --- a/Shower/SplittingFunctions/HalfOneHalfSplitFn.cc +++ b/Shower/SplittingFunctions/HalfOneHalfSplitFn.cc @@ -1,145 +1,145 @@ // -*- C++ -*- // // HalfOneHalfSplitFn.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 HalfOneHalfSplitFn class. // #include "HalfOneHalfSplitFn.h" #include "ThePEG/PDT/ParticleData.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Utilities/DescribeClass.h" #include "Herwig/Decay/TwoBodyDecayMatrixElement.h" using namespace Herwig; DescribeNoPIOClass describeHalfOneHalfSplitFn ("Herwig::HalfOneHalfSplitFn","HwShower.so"); void HalfOneHalfSplitFn::Init() { static ClassDocumentation documentation ("The HalfOneHalfSplitFn class implements the splitting " "function for q -> g q"); } double HalfOneHalfSplitFn::P(const double z, const Energy2 t, - const IdList &ids, const bool mass) const { + const IdList &ids, const bool mass, const RhoDMatrix & ) const { double val=(2.*(1.-z)+sqr(z))/z; if(mass) { Energy m = getParticleData(ids[0])->mass(); val-=2.*sqr(m)/t; } return colourFactor(ids)*val; } double HalfOneHalfSplitFn::overestimateP(const double z, const IdList &ids) const { return 2.*colourFactor(ids)/z; } double HalfOneHalfSplitFn::ratioP(const double z, const Energy2 t, - const IdList &ids,const bool mass) const { + const IdList &ids,const bool mass, const RhoDMatrix & ) const { double val=2.*(1.-z)+sqr(z); if(mass) { Energy m=getParticleData(ids[0])->mass(); val -=2.*sqr(m)*z/t; } return 0.5*val; } double HalfOneHalfSplitFn::integOverP(const double z, const IdList & ids, unsigned int PDFfactor) const { switch(PDFfactor) { case 0: return 2.*colourFactor(ids)*log(z); case 1: return -2.*colourFactor(ids)/z; case 2: return 2.*colourFactor(ids)*log(z/(1.-z)); case 3: default: throw Exception() << "HalfOneHalfSplitFn::integOverP() invalid PDFfactor = " << PDFfactor << Exception::runerror; } } double HalfOneHalfSplitFn::invIntegOverP(const double r, const IdList & ids, unsigned int PDFfactor) const { switch(PDFfactor) { case 0: return exp(0.5*r/colourFactor(ids)); case 1: return -2.*colourFactor(ids)/r; case 2: return 1./(1.+exp(-0.5*r/colourFactor(ids))); case 3: default: throw Exception() << "HalfOneHalfSplitFn::integOverP() invalid PDFfactor = " << PDFfactor << Exception::runerror; } } bool HalfOneHalfSplitFn::accept(const IdList &ids) const { // 3 particles and in and out fermion same if(ids.size()!=3 || ids[0]!=ids[2]) return false; tcPDPtr q=getParticleData(ids[0]); tcPDPtr g=getParticleData(ids[1]); if(q->iSpin()!=PDT::Spin1Half || g->iSpin()!=PDT::Spin1) return false; return checkColours(ids); } vector > HalfOneHalfSplitFn::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 > HalfOneHalfSplitFn::generatePhiBackward(const double z, const Energy2 t, const IdList & ids, const RhoDMatrix & rho) { assert(rho.iSpin()==PDT::Spin1); double mt = sqr(getParticleData(ids[0])->mass())/t; double diag = (1.+sqr(1.-z))/z - 2.*mt; double off = 2.*(1.-z)/z*(1.-mt*z/(1.-z)); double max = diag+2.*abs(rho(0,2))*off; vector > output; output.push_back(make_pair( 0, (rho(0,0)+rho(2,2))*diag/max)); output.push_back(make_pair( 2, -rho(0,2) * off/max)); output.push_back(make_pair(-2, -rho(2,0) * off/max)); return output; } DecayMEPtr HalfOneHalfSplitFn::matrixElement(const double z, const Energy2 t, const IdList & ids, const double phi, bool timeLike) { // calculate the kernal DecayMEPtr kernal(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin1Half,PDT::Spin1,PDT::Spin1Half))); Energy m = !timeLike ? ZERO : getParticleData(ids[0])->mass(); double mt = m/sqrt(t); double root = sqrt(1.-z*sqr(m)/(1.-z)/t); double romz = sqrt(1.-z); double rz = sqrt(z); Complex phase = exp(Complex(0.,1.)*phi); (*kernal)(0,0,0) = -root/rz/phase; (*kernal)(1,2,1) = -conj((*kernal)(0,0,0)); (*kernal)(0,2,0) = root/rz*(1.-z)*phase; (*kernal)(1,0,1) = -conj((*kernal)(0,2,0)); (*kernal)(1,2,0) = mt*z/romz; (*kernal)(0,0,1) = conj((*kernal)(1,2,0)); (*kernal)(0,2,1) = 0.; (*kernal)(1,0,0) = 0.; return kernal; } diff --git a/Shower/SplittingFunctions/HalfOneHalfSplitFn.h b/Shower/SplittingFunctions/HalfOneHalfSplitFn.h --- a/Shower/SplittingFunctions/HalfOneHalfSplitFn.h +++ b/Shower/SplittingFunctions/HalfOneHalfSplitFn.h @@ -1,186 +1,188 @@ // -*- C++ -*- // // HalfOneHalfSplitFn.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_HalfOneHalfSplitFn_H #define HERWIG_HalfOneHalfSplitFn_H // // This is the declaration of the HalfOneHalfSplitFn class. // #include "SplittingFunction.h" namespace Herwig { using namespace ThePEG; /** \ingroup Shower * * This classs provides the concrete implementation of the exact leading-order * splitting function for \f$\frac12\to 1\frac12\f$. * * In this case the splitting function is given by * \f[P(z,t) = C\left(\frac{2(1-z)+z^2}{z}-2\frac{m^2_q}t\right),\f] * where \f$C\f$ is the corresponding colour factor. * Our choice for the overestimate is * \f[P_{\rm over}(z) = 2C\frac1z,\f] * therefore the integral is * \f[\int P_{\rm over}(z) {\rm d}z = 2C\ln z,\f] * and its inverse is * \f[\exp\left(\frac{r}{2C}\right).\f] * * @see SplittingFunction */ class HalfOneHalfSplitFn: public SplittingFunction { public: /** * The default constructor. */ HalfOneHalfSplitFn() : SplittingFunction(1) {} /** * Concrete implementation of the method to determine whether this splitting * function can be used for a given set of particles. * @param ids The PDG codes for the particles in the splitting. */ virtual bool accept(const IdList & ids) const; /** * Methods to return the splitting function. */ //@{ /** * The concrete implementation of the splitting function, \f$P(z,t)\f$. * @param z The energy fraction. * @param t The scale. * @param ids The PDG codes for the particles in the splitting. * @param mass Whether or not to include the mass dependent terms + * @param rho The spin density matrix */ virtual double P(const double z, const Energy2 t, const IdList & ids, - const bool mass) const; + const bool mass, const RhoDMatrix & rho) const; /** * The concrete implementation of the overestimate of the splitting function, * \f$P_{\rm over}\f$. * @param z The energy fraction. * @param ids The PDG codes for the particles in the splitting. */ virtual double overestimateP(const double z, const IdList & ids) const; /** * The concrete implementation of the * the ratio of the splitting function to the overestimate, i.e. * \f$P(z,t)/P_{\rm over}(z)\f$. * @param z The energy fraction. * @param t The scale. * @param ids The PDG codes for the particles in the splitting. * @param mass Whether or not to include the mass dependent terms + * @param rho The spin density matrix */ virtual double ratioP(const double z, const Energy2 t, const IdList & ids, - const bool mass) const; + const bool mass, const RhoDMatrix & rho) const; /** * The concrete implementation of the indefinite integral of the * overestimated splitting function, \f$P_{\rm over}\f$. * @param z The energy fraction. * @param ids The PDG codes for the particles in the splitting. * @param PDFfactor Which additional factor to include for the PDF * 0 is no additional factor, * 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$ */ virtual double integOverP(const double z, const IdList & ids, unsigned int PDFfactor=0) const; /** * The concrete implementation of the inverse of the indefinite integral. * @param r Value of the splitting function to be inverted * @param ids The PDG codes for the particles in the splitting. * @param PDFfactor Which additional factor to include for the PDF * 0 is no additional factor, * 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$ */ virtual double invIntegOverP(const double r, const IdList & ids, unsigned int PDFfactor=0) const; //@} /** * Method to calculate the azimuthal angle for forward evolution * @param z The energy fraction * @param t The scale \f$t=2p_j\cdot p_k\f$. * @param ids The PDG codes for the particles in the splitting. * @param The azimuthal angle, \f$\phi\f$. * @return The weight */ virtual vector > generatePhiForward(const double z, const Energy2 t, const IdList & ids, const RhoDMatrix &); /** * Method to calculate the azimuthal angle for backward evolution * @param z The energy fraction * @param t The scale \f$t=2p_j\cdot p_k\f$. * @param ids The PDG codes for the particles in the splitting. * @param The azimuthal angle, \f$\phi\f$. * @return The weight */ virtual vector > generatePhiBackward(const double z, const Energy2 t, const IdList & ids, const RhoDMatrix &); /** * Calculate the matrix element for the splitting * @param z The energy fraction * @param t The scale \f$t=2p_j\cdot p_k\f$. * @param ids The PDG codes for the particles in the splitting. * @param The azimuthal angle, \f$\phi\f$. */ virtual DecayMEPtr matrixElement(const double z, const Energy2 t, const IdList & ids, const double phi, bool timeLike); public: /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const {return new_ptr(*this);} /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const {return new_ptr(*this);} //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ HalfOneHalfSplitFn & operator=(const HalfOneHalfSplitFn &); }; } #endif /* HERWIG_HalfOneHalfSplitFn_H */ diff --git a/Shower/SplittingFunctions/OneHalfHalfSplitFn.cc b/Shower/SplittingFunctions/OneHalfHalfSplitFn.cc --- a/Shower/SplittingFunctions/OneHalfHalfSplitFn.cc +++ b/Shower/SplittingFunctions/OneHalfHalfSplitFn.cc @@ -1,144 +1,144 @@ // -*- C++ -*- // // OneHalfHalfSplitFn.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 OneHalfHalfSplitFn class. // #include "OneHalfHalfSplitFn.h" #include "ThePEG/PDT/ParticleData.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Utilities/DescribeClass.h" #include "Herwig/Decay/TwoBodyDecayMatrixElement.h" using namespace Herwig; DescribeNoPIOClass describeOneHalfHalfSplitFn ("Herwig::OneHalfHalfSplitFn","HwShower.so"); void OneHalfHalfSplitFn::Init() { static ClassDocumentation documentation ("The OneHalfHalfSplitFn class implements the splitting function for g->q qbar"); } double OneHalfHalfSplitFn::P(const double z, const Energy2 t, - const IdList &ids, const bool mass) const { + const IdList &ids, const bool mass, const RhoDMatrix &) const { double zz = z*(1.-z); double val=1.-2.*zz; if(mass) { Energy m = getParticleData(ids[1])->mass(); val +=2.*sqr(m)/t; } return colourFactor(ids)*val; } double OneHalfHalfSplitFn::overestimateP(const double, const IdList &ids) const { return colourFactor(ids); } double OneHalfHalfSplitFn::ratioP(const double z, const Energy2 t, - const IdList &ids, const bool mass) const { + const IdList &ids, const bool mass, const RhoDMatrix &) const { double zz = z*(1.-z); double val = 1.-2.*zz; if(mass) { Energy m = getParticleData(ids[1])->mass(); val+= 2.*sqr(m)/t; } return val; } double OneHalfHalfSplitFn::integOverP(const double z, const IdList & ids, unsigned int PDFfactor) const { switch(PDFfactor) { case 0: return colourFactor(ids)*z; case 1: return colourFactor(ids)*log(z); case 2: return -colourFactor(ids)*log(1.-z); case 3: return colourFactor(ids)*log(z/(1.-z)); default: throw Exception() << "OneHalfHalfSplitFn::integOverP() invalid PDFfactor = " << PDFfactor << Exception::runerror; } } double OneHalfHalfSplitFn::invIntegOverP(const double r, const IdList & ids, unsigned int PDFfactor) const { switch(PDFfactor) { case 0: return r/colourFactor(ids); case 1: return exp(r/colourFactor(ids)); case 2: return 1.-exp(-r/colourFactor(ids)); case 3: return 1./(1.+exp(-r/colourFactor(ids))); default: throw Exception() << "OneHalfHalfSplitFn::integOverP() invalid PDFfactor = " << PDFfactor << Exception::runerror; } } bool OneHalfHalfSplitFn::accept(const IdList &ids) const { if(ids.size()!=3) return false; if(ids[1]!=-ids[2]) return false; tcPDPtr q=getParticleData(ids[1]); if(q->iSpin()!=PDT::Spin1Half) return false; tcPDPtr g=getParticleData(ids[0]); if(g->iSpin()!=PDT::Spin1) return false; return checkColours(ids); } vector > OneHalfHalfSplitFn::generatePhiForward(const double z, const Energy2 t, const IdList & ids, const RhoDMatrix & rho) { assert(rho.iSpin()==PDT::Spin1); double modRho = abs(rho(0,2)); Energy mq = getParticleData(ids[1])->mass(); Energy2 mq2 = sqr(mq); double fact = z*(1.-z)-mq2/t; double max = 1.+2.*fact*(-1.+2.*modRho); vector > output; output.push_back(make_pair( 0,(rho(0,0)+rho(2,2))*(1.-2.*fact)/max)); output.push_back(make_pair(-2,2.*fact*rho(0,2)/max)); output.push_back(make_pair( 2,2.*fact*rho(2,0)/max)); return output; } vector > OneHalfHalfSplitFn::generatePhiBackward(const double, const Energy2, const IdList &, const RhoDMatrix & ) { // no dependance return vector >(1,make_pair(0,1.)); } DecayMEPtr OneHalfHalfSplitFn::matrixElement(const double z, const Energy2 t, const IdList & ids, const double phi, bool timeLike) { static const Complex ii(0.,1.); // calculate the kernal DecayMEPtr kernal(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin1,PDT::Spin1Half,PDT::Spin1Half))); double mt = !timeLike ? ZERO : getParticleData(ids[1])->mass()/sqrt(t); double root =sqrt(1.-sqr(mt)/z/(1.-z)); (*kernal)(0,0,0) = mt/sqrt(z*(1.-z)); (*kernal)(2,1,1) = (*kernal)(0,0,0); (*kernal)(0,0,1) = -z*root*exp(-ii*phi); (*kernal)(2,1,0) = -conj((*kernal)(0,0,1)); (*kernal)(0,1,0) = (1.-z)*exp(-ii*phi)*root; (*kernal)(2,0,1) = -conj((*kernal)(0,1,0)); (*kernal)(0,1,1) = 0.; (*kernal)(2,0,0) = 0.; return kernal; } diff --git a/Shower/SplittingFunctions/OneHalfHalfSplitFn.h b/Shower/SplittingFunctions/OneHalfHalfSplitFn.h --- a/Shower/SplittingFunctions/OneHalfHalfSplitFn.h +++ b/Shower/SplittingFunctions/OneHalfHalfSplitFn.h @@ -1,191 +1,193 @@ // -*- C++ -*- // // OneHalfHalfSplitFn.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_OneHalfHalfSplitFn_H #define HERWIG_OneHalfHalfSplitFn_H // // This is the declaration of the OneHalfHalfSplitFn class. // #include "SplittingFunction.h" namespace Herwig { using namespace ThePEG; /**\ingroup Shower * * This class provides the concrete implementation of the exact leading-order * splitting function for \f$1\to \frac12\frac12\f$. * * In this case the splitting function is given by * \f[P(z,t) =C\left(1-2z(1-z)+2\frac{m_q^2}{t}\right),\f] * where \f$C\f$ is the corresponding colour factor * Our choice for the overestimate is * \f[P_{\rm over}(z) = C,\f] * therefore the integral is * \f[\int P_{\rm over}(z) {\rm d}z =Cz,\f] * and its inverse is * \f[\frac{r}{C}\f] * * @see \ref OneHalfHalfSplitFnInterfaces "The interfaces" * defined for OneHalfHalfSplitFn. */ class OneHalfHalfSplitFn: public SplittingFunction { public: /** * The default constructor. */ OneHalfHalfSplitFn() : SplittingFunction(1) {} /** * Concrete implementation of the method to determine whether this splitting * function can be used for a given set of particles. * @param ids The PDG codes for the particles in the splitting. */ virtual bool accept(const IdList & ids) const; /** * Methods to return the splitting function. */ //@{ /** * The concrete implementation of the splitting function, \f$P\f$. * @param z The energy fraction. * @param t The scale. * @param ids The PDG codes for the particles in the splitting. * @param mass Whether or not to include the mass dependent terms + * @param rho The spin density matrix */ virtual double P(const double z, const Energy2 t, const IdList & ids, - const bool mass) const; + const bool mass, const RhoDMatrix & rho) const; /** * The concrete implementation of the overestimate of the splitting function, * \f$P_{\rm over}\f$. * @param z The energy fraction. * @param ids The PDG codes for the particles in the splitting. */ virtual double overestimateP(const double z, const IdList & ids) const; /** * The concrete implementation of the * the ratio of the splitting function to the overestimate, i.e. * \f$P(z,\tilde{q}^2)/P_{\rm over}(z)\f$. * @param z The energy fraction. * @param t The scale. * @param ids The PDG codes for the particles in the splitting. * @param mass Whether or not to include the mass dependent terms + * @param rho The spin density matrix */ virtual double ratioP(const double z, const Energy2 t, const IdList & ids, - const bool mass) const; + const bool mass, const RhoDMatrix & rho) const; /** * The concrete implementation of the indefinite integral of the * overestimated splitting function, \f$P_{\rm over}\f$. * @param z The energy fraction. * @param ids The PDG codes for the particles in the splitting. * @param PDFfactor Which additional factor to include for the PDF * 0 is no additional factor, * 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$ */ virtual double integOverP(const double z, const IdList & ids, unsigned int PDFfactor=0) const; /** * The concrete implementation of the inverse of the indefinite integral. * @param r Value of the splitting function to be inverted * @param ids The PDG codes for the particles in the splitting. * @param PDFfactor Which additional factor to include for the PDF * 0 is no additional factor, * 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$ */ virtual double invIntegOverP(const double r, const IdList & ids, unsigned int PDFfactor=0) const; //@} /** * Method to calculate the azimuthal angle * @param particle The particle which is branching * @param showerkin The ShowerKinematics object * @param z The energy fraction * @param t The scale \f$t=2p_j\cdot p_k\f$. * @param ids The PDG codes for the particles in the splitting. * @param The azimuthal angle, \f$\phi\f$. * @return The weight */ virtual vector > generatePhiForward(const double z, const Energy2 t, const IdList & ids, const RhoDMatrix &); /** * Method to calculate the azimuthal angle * @param particle The particle which is branching * @param showerkin The ShowerKinematics object * @param z The energy fraction * @param t The scale \f$t=2p_j\cdot p_k\f$. * @param ids The PDG codes for the particles in the splitting. * @param The azimuthal angle, \f$\phi\f$. * @return The weight */ virtual vector > generatePhiBackward(const double z, const Energy2 t, const IdList & ids, const RhoDMatrix &); /** * Calculate the matrix element for the splitting * @param z The energy fraction * @param t The scale \f$t=2p_j\cdot p_k\f$. * @param ids The PDG codes for the particles in the splitting. * @param The azimuthal angle, \f$\phi\f$. */ virtual DecayMEPtr matrixElement(const double z, const Energy2 t, const IdList & ids, const double phi, bool timeLike); public: /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const {return new_ptr(*this);} /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const {return new_ptr(*this);} //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ OneHalfHalfSplitFn & operator=(const OneHalfHalfSplitFn &); }; } #endif /* HERWIG_OneHalfHalfSplitFn_H */ diff --git a/Shower/SplittingFunctions/OneOneOneSplitFn.cc b/Shower/SplittingFunctions/OneOneOneSplitFn.cc --- a/Shower/SplittingFunctions/OneOneOneSplitFn.cc +++ b/Shower/SplittingFunctions/OneOneOneSplitFn.cc @@ -1,134 +1,134 @@ // -*- C++ -*- // // OneOneOneSplitFn.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 OneOneOneSplitFn class. // #include "OneOneOneSplitFn.h" #include "ThePEG/PDT/ParticleData.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Utilities/DescribeClass.h" #include "Herwig/Decay/TwoBodyDecayMatrixElement.h" using namespace Herwig; DescribeNoPIOClass describeOneOneOneSplitFn ("Herwig::OneOneOneSplitFn","HwShower.so"); void OneOneOneSplitFn::Init() { static ClassDocumentation documentation ("The OneOneOneSplitFn class implements the g -> gg splitting function"); } double OneOneOneSplitFn::P(const double z, const Energy2, - const IdList & ids, const bool)const { + const IdList & ids, const bool, const RhoDMatrix &)const { // (this is historically important! the first physics - two years // after the birth of the project - in the Herwig shower! Alberto // & Stefan, 25/04/2002). return colourFactor(ids)*sqr(1.-z*(1.-z))/(z*(1.-z)); } double OneOneOneSplitFn::overestimateP(const double z, const IdList & ids) const { return colourFactor(ids)*(1/z + 1/(1.-z)); } double OneOneOneSplitFn::ratioP(const double z, const Energy2, - const IdList & , const bool) const { + const IdList & , const bool, const RhoDMatrix &) const { return sqr(1.-z*(1.-z)); } double OneOneOneSplitFn::invIntegOverP(const double r, const IdList & ids, unsigned int PDFfactor) const { switch(PDFfactor) { case 0: return 1./(1.+exp(-r/colourFactor(ids))); case 1: case 2: case 3: default: throw Exception() << "OneOneOneSplitFn::integOverP() invalid PDFfactor = " << PDFfactor << Exception::runerror; } } double OneOneOneSplitFn::integOverP(const double z, const IdList & ids, unsigned int PDFfactor) const { switch(PDFfactor) { case 0: assert(z>0.&&z<1.); return colourFactor(ids)*log(z/(1.-z)); case 1: case 2: case 3: default: throw Exception() << "OneOneOneSplitFn::integOverP() invalid PDFfactor = " << PDFfactor << Exception::runerror; } } bool OneOneOneSplitFn::accept(const IdList & ids) const { if(ids.size()!=3) return false; for(unsigned int ix=0;ixiSpin()!=PDT::Spin1) return false; } return checkColours(ids); } vector > OneOneOneSplitFn::generatePhiForward(const double z, const Energy2, const IdList &, const RhoDMatrix & rho) { assert(rho.iSpin()==PDT::Spin1); double modRho = abs(rho(0,2)); double max = 2.*z*modRho*(1.-z)+sqr(1.-(1.-z)*z)/(z*(1.-z)); vector > output; output.push_back(make_pair( 0,(rho(0,0)+rho(2,2))*sqr(1.-(1.-z)*z)/(z*(1.-z))/max)); output.push_back(make_pair(-2,-rho(0,2)*z*(1.-z)/max)); output.push_back(make_pair( 2,-rho(2,0)*z*(1.-z)/max)); return output; } vector > OneOneOneSplitFn::generatePhiBackward(const double z, const Energy2, const IdList &, const RhoDMatrix & rho) { assert(rho.iSpin()==PDT::Spin1); double diag = sqr(1 - (1 - z)*z)/(1 - z)/z; double off = (1.-z)/z; double max = 2.*abs(rho(0,2))*off+diag; vector > output; output.push_back(make_pair( 0, (rho(0,0)+rho(2,2))*diag/max)); output.push_back(make_pair( 2,-rho(0,2) * off/max)); output.push_back(make_pair(-2,-rho(2,0) * off/max)); return output; } DecayMEPtr OneOneOneSplitFn::matrixElement(const double z, const Energy2, const IdList &, const double phi, bool) { // calculate the kernal DecayMEPtr kernal(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin1,PDT::Spin1,PDT::Spin1))); double omz = 1.-z; double root = sqrt(z*omz); Complex phase = exp(Complex(0.,1.)*phi); (*kernal)(0,0,0) = phase/root; (*kernal)(2,2,2) = -conj((*kernal)(0,0,0)); (*kernal)(0,0,2) = -sqr(z)/root/phase; (*kernal)(2,2,0) = -conj((*kernal)(0,0,2)); (*kernal)(0,2,0) = -sqr(omz)/root/phase; (*kernal)(2,0,2) = -conj((*kernal)(0,2,0)); (*kernal)(0,2,2) = 0.; (*kernal)(2,0,0) = 0.; return kernal; } diff --git a/Shower/SplittingFunctions/OneOneOneSplitFn.h b/Shower/SplittingFunctions/OneOneOneSplitFn.h --- a/Shower/SplittingFunctions/OneOneOneSplitFn.h +++ b/Shower/SplittingFunctions/OneOneOneSplitFn.h @@ -1,187 +1,189 @@ // -*- C++ -*- // // OneOneOneSplitFn.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_OneOneOneSplitFn_H #define HERWIG_OneOneOneSplitFn_H // // This is the declaration of the OneOneOneSplitFn class. // #include "SplittingFunction.h" namespace Herwig { using namespace ThePEG; /** \ingroup Shower * * This class provides the concrete implementation of the exact leading-order * splitting function for \f$1\to 11\f$. * * In this case the splitting function is given by * \f[P(z) =2C*\left(\frac{z}{1-z}+\frac{1-z}{z}+z(1-z)\right),\f] * where \f$C=\f$ is the corresponding colour factor. * Our choice for the overestimate is * \f[P_{\rm over}(z) = 2C\left(\frac1z+\frac1{1-z}\right),\f] * therefore the integral is * \f[\int P_{\rm over}(z) {\rm d}z =2C\ln\left(\frac{z}{1-z}\right),\f] * and its inverse is * \f[\frac{\exp\left(\frac{r}{2C}\right)}{(1+\exp\left(\frac{r}{2C}\right)}\f] * * * @see \ref OneOneOneSplitFnInterfaces "The interfaces" * defined for OneOneOneSplitFn. */ class OneOneOneSplitFn: public SplittingFunction { public: /** * The default constructor. */ OneOneOneSplitFn() : SplittingFunction(1) {} /** * Concrete implementation of the method to determine whether this splitting * function can be used for a given set of particles. * @param ids The PDG codes for the particles in the splitting. */ virtual bool accept(const IdList & ids) const; /** * Methods to return the splitting function. */ //@{ /** * The concrete implementation of the splitting function, \f$P(z,t)\f$. * @param z The energy fraction. * @param t The scale. * @param ids The PDG codes for the particles in the splitting. * @param mass Whether or not to include the mass dependent terms + * @param rho The spin density matrix */ virtual double P(const double z, const Energy2 t, const IdList & ids, - const bool mass) const; + const bool mass, const RhoDMatrix & rho) const; /** * The concrete implementation of the overestimate of the splitting function, * \f$P_{\rm over}\f$. * @param z The energy fraction. * @param ids The PDG codes for the particles in the splitting. */ virtual double overestimateP(const double z, const IdList & ids) const; /** * The concrete implementation of the * the ratio of the splitting function to the overestimate, i.e. * \f$P(z,t)/P_{\rm over}(z)\f$. * @param z The energy fraction. * @param t The scale. * @param ids The PDG codes for the particles in the splitting. * @param mass Whether or not to include the mass dependent terms + * @param rho The spin density matrix */ virtual double ratioP(const double z, const Energy2 t, const IdList & ids, - const bool mass) const; + const bool mass, const RhoDMatrix & rho) const; /** * The concrete implementation of the indefinite integral of the * overestimated splitting function, \f$P_{\rm over}\f$. * @param z The energy fraction. * @param ids The PDG codes for the particles in the splitting. * @param PDFfactor Which additional factor to include for the PDF * 0 is no additional factor, * 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$ */ virtual double integOverP(const double z, const IdList & ids, unsigned int PDFfactor=0) const; /** * The concrete implementation of the inverse of the indefinite integral. * @param r Value of the splitting function to be inverted * @param ids The PDG codes for the particles in the splitting. * @param PDFfactor Which additional factor to include for the PDF * 0 is no additional factor, * 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$ */ virtual double invIntegOverP(const double r, const IdList & ids, unsigned int PDFfactor=0) const; //@} /** * Method to calculate the azimuthal angle for forward evolution * @param z The energy fraction * @param t The scale \f$t=2p_j\cdot p_k\f$. * @param ids The PDG codes for the particles in the splitting. * @param The azimuthal angle, \f$\phi\f$. * @return The weight */ virtual vector > generatePhiForward(const double z, const Energy2 t, const IdList & ids, const RhoDMatrix &); /** * Method to calculate the azimuthal angle for backward evolution * @param z The energy fraction * @param t The scale \f$t=2p_j\cdot p_k\f$. * @param ids The PDG codes for the particles in the splitting. * @param The azimuthal angle, \f$\phi\f$. * @return The weight */ virtual vector > generatePhiBackward(const double z, const Energy2 t, const IdList & ids, const RhoDMatrix &); /** * Calculate the matrix element for the splitting * @param z The energy fraction * @param t The scale \f$t=2p_j\cdot p_k\f$. * @param ids The PDG codes for the particles in the splitting. * @param The azimuthal angle, \f$\phi\f$. */ virtual DecayMEPtr matrixElement(const double z, const Energy2 t, const IdList & ids, const double phi, bool timeLike); public: /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const {return new_ptr(*this);} /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const {return new_ptr(*this);} //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ OneOneOneSplitFn & operator=(const OneOneOneSplitFn &); }; } #endif /* HERWIG_OneOneOneSplitFn_H */ diff --git a/Shower/SplittingFunctions/SplittingFunction.h b/Shower/SplittingFunctions/SplittingFunction.h --- a/Shower/SplittingFunctions/SplittingFunction.h +++ b/Shower/SplittingFunctions/SplittingFunction.h @@ -1,372 +1,374 @@ // -*- C++ -*- // // SplittingFunction.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_SplittingFunction_H #define HERWIG_SplittingFunction_H // // This is the declaration of the SplittingFunction class. // #include "ThePEG/Interface/Interfaced.h" #include "Herwig/Shower/ShowerConfig.h" #include "ThePEG/EventRecord/RhoDMatrix.h" #include "Herwig/Decay/DecayMatrixElement.h" #include "Herwig/Shower/Base/ShowerKinematics.fh" #include "ThePEG/EventRecord/ColourLine.h" #include "ThePEG/PDT/ParticleData.h" #include "SplittingFunction.fh" namespace Herwig { using namespace ThePEG; /** \ingroup Shower * Enum to define the possible types of colour structure which can occur in * the branching. */ enum ColourStructure {Undefined=0, TripletTripletOctet = 1,OctetOctetOctet =2, OctetTripletTriplet = 3,TripletOctetTriplet=4, SextetSextetOctet = 5, ChargedChargedNeutral=-1,ChargedNeutralCharged=-2, NeutralChargedCharged=-3,EW=-4}; /** \ingroup Shower * * This is an abstract class which defines the common interface * for all \f$1\to2\f$ splitting functions, for both initial-state * and final-state radiation. * * The SplittingFunction class contains a number of purely virtual members * which must be implemented in the inheriting classes. The class also stores * the interaction type of the spltting function. * * The inheriting classes need to specific the splitting function * \f$P(z,2p_j\cdot p_k)\f$, in terms of the energy fraction \f$z\f$ and * the evolution scale. In order to allow the splitting functions to be used * with different choices of evolution functions the scale is given by * \f[2p_j\cdot p_k=(p_j+p_k)^2-m_{jk}^2=Q^2-(p_j+p_k)^2=z(1-z)\tilde{q}^2= * \frac{p_T^2}{z(1-z)}-m_{jk}^2+\frac{m_j^2}{z}+\frac{m_k^2}{1-z},\f] * where \f$Q^2\f$ is the virtuality of the branching particle, * $p_T$ is the relative transverse momentum of the branching products and * \f$\tilde{q}^2\f$ is the angular variable described in hep-ph/0310083. * * In addition an overestimate of the * splitting function, \f$P_{\rm over}(z)\f$ which only depends upon \f$z\f$, * the integral, inverse of the integral for this overestimate and * ratio of the true splitting function to the overestimate must be provided * as they are necessary for the veto alogrithm used to implement the evolution. * * @see \ref SplittingFunctionInterfaces "The interfaces" * defined for SplittingFunction. */ class SplittingFunction: public Interfaced { public: /** * The default constructor. * @param b All splitting functions must have an interaction order */ SplittingFunction(unsigned int b) : Interfaced(), _interactionType(ShowerInteraction::UNDEFINED), _interactionOrder(b), _colourStructure(Undefined), _colourFactor(-1.), angularOrdered_(true) {} public: /** * Methods to return the interaction type and order for the splitting function */ //@{ /** * Return the type of the interaction */ ShowerInteraction::Type interactionType() const {return _interactionType;} /** * Return the order of the splitting function in the interaction */ unsigned int interactionOrder() const {return _interactionOrder;} /** * Return the colour structure */ ColourStructure colourStructure() const {return _colourStructure;} /** * Return the colour factor */ double colourFactor(const IdList &ids) const { if(_colourStructure>0) return _colourFactor; else if(_colourStructure<0) { if(_colourStructure==ChargedChargedNeutral || _colourStructure==ChargedNeutralCharged) { tPDPtr part=getParticleData(ids[0]); return sqr(double(part->iCharge())/3.); } else { tPDPtr part=getParticleData(ids[1]); return sqr(double(part->iCharge())/3.); } } else assert(false); } //@} /** * Purely virtual method which should determine whether this splitting * function can be used for a given set of particles. * @param ids The PDG codes for the particles in the splitting. */ virtual bool accept(const IdList & ids) const = 0; /** * Method to check the colours are correct */ virtual bool checkColours(const IdList & ids) const; /** * Methods to return the splitting function. */ //@{ /** * Purely virtual method which should return the exact value of the splitting function, * \f$P\f$ evaluated in terms of the energy fraction, \f$z\f$, and the evolution scale \f$\tilde{q}^2\f$. * @param z The energy fraction. * @param t The scale \f$t=2p_j\cdot p_k\f$. * @param ids The PDG codes for the particles in the splitting. * @param mass Whether or not to include the mass dependent terms + * @param rho The spin density matrix */ virtual double P(const double z, const Energy2 t, const IdList & ids, - const bool mass) const = 0; + const bool mass, const RhoDMatrix & rho) const = 0; /** * Purely virtual method which should return * an overestimate of the splitting function, * \f$P_{\rm over}\f$ such that the result \f$P_{\rm over}\geq P\f$. This function * should be simple enough that it does not depend on the evolution scale. * @param z The energy fraction. * @param ids The PDG codes for the particles in the splitting. */ virtual double overestimateP(const double z, const IdList & ids) const = 0; /** * Purely virtual method which should return * the ratio of the splitting function to the overestimate, i.e. * \f$P(z,\tilde{q}^2)/P_{\rm over}(z)\f$. * @param z The energy fraction. * @param t The scale \f$t=2p_j\cdot p_k\f$. * @param ids The PDG codes for the particles in the splitting. * @param mass Whether or not to include the mass dependent terms + * @param rho The spin density matrix */ virtual double ratioP(const double z, const Energy2 t, const IdList & ids, - const bool mass) const = 0; + const bool mass, const RhoDMatrix & rho) const = 0; /** * Purely virtual method which should return the indefinite integral of the * overestimated splitting function, \f$P_{\rm over}\f$. * @param z The energy fraction. * @param ids The PDG codes for the particles in the splitting. * @param PDFfactor Which additional factor to include for the PDF * 0 is no additional factor, * 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$ * */ virtual double integOverP(const double z, const IdList & ids, unsigned int PDFfactor=0) const = 0; /** * Purely virtual method which should return the inverse of the * indefinite integral of the * overestimated splitting function, \f$P_{\rm over}\f$ which is used to * generate the value of \f$z\f$. * @param r Value of the splitting function to be inverted * @param ids The PDG codes for the particles in the splitting. * @param PDFfactor Which additional factor to include for the PDF * 0 is no additional factor, * 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$ */ virtual double invIntegOverP(const double r, const IdList & ids, unsigned int PDFfactor=0) const = 0; //@} /** * Purely virtual method which should make the proper colour connection * between the emitting parent and the branching products. * @param parent The parent for the branching * @param first The first branching product * @param second The second branching product * @param partnerType The type of evolution partner * @param back Whether this is foward or backward evolution. */ virtual void colourConnection(tShowerParticlePtr parent, tShowerParticlePtr first, tShowerParticlePtr second, ShowerPartnerType::Type partnerType, const bool back) const; /** * Method to calculate the azimuthal angle for forward evolution * @param z The energy fraction * @param t The scale \f$t=2p_j\cdot p_k\f$. * @param ids The PDG codes for the particles in the splitting. * @param The azimuthal angle, \f$\phi\f$. * @return The weight */ virtual vector > generatePhiForward(const double z, const Energy2 t, const IdList & ids, const RhoDMatrix &) = 0; /** * Method to calculate the azimuthal angle for backward evolution * @param z The energy fraction * @param t The scale \f$t=2p_j\cdot p_k\f$. * @param ids The PDG codes for the particles in the splitting. * @return The weight */ virtual vector > generatePhiBackward(const double z, const Energy2 t, const IdList & ids, const RhoDMatrix &) = 0; /** * Calculate the matrix element for the splitting * @param z The energy fraction * @param t The scale \f$t=2p_j\cdot p_k\f$. * @param ids The PDG codes for the particles in the splitting. * @param phi The azimuthal angle, \f$\phi\f$. * @param timeLike Whether timelike or spacelike, affects inclusive of mass terms */ virtual DecayMEPtr matrixElement(const double z, const Energy2 t, const IdList & ids, const double phi, bool timeLike) = 0; /** * Whether or not the interaction is angular ordered */ bool angularOrdered() const {return angularOrdered_;} /** * Functions to state scales after branching happens */ //@{ /** * Sort out scales for final-state emission */ void evaluateFinalStateScales(ShowerPartnerType::Type type, Energy scale, double z, tShowerParticlePtr parent, tShowerParticlePtr first, tShowerParticlePtr second); /** * Sort out scales for initial-state emission */ void evaluateInitialStateScales(ShowerPartnerType::Type type, Energy scale, double z, tShowerParticlePtr parent, tShowerParticlePtr first, tShowerParticlePtr second); /** * Sort out scales for decay emission */ void evaluateDecayScales(ShowerPartnerType::Type type, Energy scale, double z, tShowerParticlePtr parent, tShowerParticlePtr first, tShowerParticlePtr second); //@} public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); //@} protected: /** * Set the colour factor */ void colourFactor(double in) {_colourFactor=in;} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ SplittingFunction & operator=(const SplittingFunction &); private: /** * The interaction type for the splitting function. */ ShowerInteraction::Type _interactionType; /** * The order of the splitting function in the coupling */ unsigned int _interactionOrder; /** * The colour structure */ ColourStructure _colourStructure; /** * The colour factor */ double _colourFactor; /** * Whether or not this interaction is angular-ordered */ bool angularOrdered_; }; } #endif /* HERWIG_SplittingFunction_H */ diff --git a/Shower/SplittingFunctions/SplittingGenerator.cc b/Shower/SplittingFunctions/SplittingGenerator.cc --- a/Shower/SplittingFunctions/SplittingGenerator.cc +++ b/Shower/SplittingFunctions/SplittingGenerator.cc @@ -1,611 +1,617 @@ // -*- C++ -*- // // SplittingGenerator.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 SplittingGenerator class. // #include "SplittingGenerator.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Command.h" #include "ThePEG/Utilities/StringUtils.h" #include "ThePEG/Repository/Repository.h" #include "Herwig/Shower/Base/ShowerParticle.h" #include "ThePEG/Utilities/Rebinder.h" #include #include "ThePEG/Utilities/DescribeClass.h" using namespace Herwig; namespace { bool checkInteraction(ShowerInteraction::Type allowed, ShowerInteraction::Type splitting) { if(allowed==ShowerInteraction::ALL) return true; else if(allowed==ShowerInteraction::QEDQCD && (splitting==ShowerInteraction::QED || splitting==ShowerInteraction::QCD )) return true; else if(allowed == splitting) return true; else return false; } } DescribeClass describeSplittingGenerator ("Herwig::SplittingGenerator",""); IBPtr SplittingGenerator::clone() const { return new_ptr(*this); } IBPtr SplittingGenerator::fullclone() const { return new_ptr(*this); } void SplittingGenerator::persistentOutput(PersistentOStream & os) const { os << _isr_Mode << _fsr_Mode << _bbranchings << _fbranchings; } void SplittingGenerator::persistentInput(PersistentIStream & is, int) { is >> _isr_Mode >> _fsr_Mode >> _bbranchings >> _fbranchings; } void SplittingGenerator::Init() { static ClassDocumentation documentation ("There class is responsible for initializing the Sudakov form factors ", "and generating splittings."); static Switch interfaceISRMode ("ISR", "Include initial-state radiation?", &SplittingGenerator::_isr_Mode, 1, false, false); static SwitchOption interfaceISRMode0 (interfaceISRMode,"No","ISR (Initial State Radiation) is OFF", 0); static SwitchOption interfaceISRMode1 (interfaceISRMode,"Yes","ISR (Initial State Radiation) is ON", 1); static Switch interfaceFSRMode ("FSR", "Include final-state radiation?", &SplittingGenerator::_fsr_Mode, 1, false, false); static SwitchOption interfaceFSRMode0 (interfaceFSRMode,"No","FSR (Final State Radiation) is OFF", 0); static SwitchOption interfaceFSRMode1 (interfaceFSRMode,"Yes","FSR (Final State Radiation) is ON", 1); static Command interfaceAddSplitting ("AddFinalSplitting", "Adds another splitting to the list of splittings considered " "in the shower. Command is a->b,c; Sudakov", &SplittingGenerator::addFinalSplitting); static Command interfaceAddInitialSplitting ("AddInitialSplitting", "Adds another splitting to the list of initial splittings to consider " "in the shower. Command is a->b,c; Sudakov. Here the particle a is the " "particle that is PRODUCED by the splitting. b is the initial state " "particle that is splitting in the shower.", &SplittingGenerator::addInitialSplitting); static Command interfaceDeleteSplitting ("DeleteFinalSplitting", "Deletes a splitting from the list of splittings considered " "in the shower. Command is a->b,c; Sudakov", &SplittingGenerator::deleteFinalSplitting); static Command interfaceDeleteInitialSplitting ("DeleteInitialSplitting", "Deletes a splitting from the list of initial splittings to consider " "in the shower. Command is a->b,c; Sudakov. Here the particle a is the " "particle that is PRODUCED by the splitting. b is the initial state " "particle that is splitting in the shower.", &SplittingGenerator::deleteInitialSplitting); } string SplittingGenerator::addSplitting(string arg, bool final) { string partons = StringUtils::car(arg); string sudakov = StringUtils::cdr(arg); vector products; string::size_type next = partons.find("->"); if(next == string::npos) return "Error: Invalid string for splitting " + arg; if(partons.find(';') == string::npos) return "Error: Invalid string for splitting " + arg; tPDPtr parent = Repository::findParticle(partons.substr(0,next)); partons = partons.substr(next+2); do { next = min(partons.find(','), partons.find(';')); tPDPtr pdp = Repository::findParticle(partons.substr(0,next)); partons = partons.substr(next+1); if(pdp) products.push_back(pdp); else return "Error: Could not create splitting from " + arg; } while(partons[0] != ';' && partons.size()); SudakovPtr s; s = dynamic_ptr_cast(Repository::TraceObject(sudakov)); if(!s) return "Error: Could not load Sudakov " + sudakov + '\n'; IdList ids; ids.push_back(parent->id()); for(vector::iterator it = products.begin(); it!=products.end(); ++it) ids.push_back((*it)->id()); // check splitting can handle this if(!s->splittingFn()->accept(ids)) return "Error: Sudakov " + sudakov + "can't handle particles\n"; // add to map addToMap(ids,s,final); return ""; } string SplittingGenerator::deleteSplitting(string arg, bool final) { string partons = StringUtils::car(arg); string sudakov = StringUtils::cdr(arg); vector products; string::size_type next = partons.find("->"); if(next == string::npos) return "Error: Invalid string for splitting " + arg; if(partons.find(';') == string::npos) return "Error: Invalid string for splitting " + arg; tPDPtr parent = Repository::findParticle(partons.substr(0,next)); partons = partons.substr(next+2); do { next = min(partons.find(','), partons.find(';')); tPDPtr pdp = Repository::findParticle(partons.substr(0,next)); partons = partons.substr(next+1); if(pdp) products.push_back(pdp); else return "Error: Could not create splitting from " + arg; } while(partons[0] != ';' && partons.size()); SudakovPtr s; s = dynamic_ptr_cast(Repository::TraceObject(sudakov)); if(!s) return "Error: Could not load Sudakov " + sudakov + '\n'; IdList ids; ids.push_back(parent->id()); for(vector::iterator it = products.begin(); it!=products.end(); ++it) ids.push_back((*it)->id()); // check splitting can handle this if(!s->splittingFn()->accept(ids)) return "Error: Sudakov " + sudakov + "can't handle particles\n"; // delete from map deleteFromMap(ids,s,final); return ""; } void SplittingGenerator::addToMap(const IdList &ids, const SudakovPtr &s, bool final) { if(isISRadiationON() && !final) { _bbranchings.insert(BranchingInsert(ids[1],BranchingElement(s,ids))); s->addSplitting(ids); } if(isFSRadiationON() && final) { _fbranchings.insert(BranchingInsert(ids[0],BranchingElement(s,ids))); s->addSplitting(ids); } } void SplittingGenerator::deleteFromMap(const IdList &ids, const SudakovPtr &s, bool final) { if(isISRadiationON() && !final) { pair range = _bbranchings.equal_range(ids[1]); for(BranchingList::iterator it=range.first; it!=range.second&&it!=_bbranchings.end()&&it->first==ids[1];++it) { if(it->second.first==s&&it->second.second==ids) { BranchingList::iterator it2=it; --it; _bbranchings.erase(it2); } } s->removeSplitting(ids); } if(isFSRadiationON() && final) { pair range = _fbranchings.equal_range(ids[0]); for(BranchingList::iterator it=range.first; it!=range.second&&it!=_fbranchings.end()&&it->first==ids[0];++it) { if(it->second.first==s&&it->second.second==ids) { BranchingList::iterator it2 = it; --it; _fbranchings.erase(it2); } } s->removeSplitting(ids); } } Branching SplittingGenerator::chooseForwardBranching(ShowerParticle &particle, double enhance, ShowerInteraction::Type type) const { + RhoDMatrix rho = particle.extractRhoMatrix(true); Energy newQ = ZERO; ShoKinPtr kinematics = ShoKinPtr(); ShowerPartnerType::Type partnerType(ShowerPartnerType::Undefined); SudakovPtr sudakov = SudakovPtr(); IdList ids; // First, find the eventual branching, corresponding to the highest scale. long index = abs(particle.data().id()); // if no branchings return empty branching struct if( _fbranchings.find(index) == _fbranchings.end() ) return Branching(ShoKinPtr(), IdList(),SudakovPtr(),ShowerPartnerType::Undefined); // otherwise select branching for(BranchingList::const_iterator cit = _fbranchings.lower_bound(index); cit != _fbranchings.upper_bound(index); ++cit) { // check either right interaction or doing both if(!checkInteraction(type,cit->second.first->interactionType())) continue; // whether or not this interaction should be angular ordered bool angularOrdered = cit->second.first->splittingFn()->angularOrdered(); ShoKinPtr newKin; ShowerPartnerType::Type type; // work out which starting scale we need if(cit->second.first->interactionType()==ShowerInteraction::QED) { type = ShowerPartnerType::QED; Energy startingScale = angularOrdered ? particle.scales().QED : particle.scales().QED_noAO; newKin = cit->second.first-> generateNextTimeBranching(startingScale,cit->second.second, - particle.id()!=cit->first,enhance, + particle.id()!=cit->first,rho,enhance, particle.scales().Max_Q2); } else if(cit->second.first->interactionType()==ShowerInteraction::QCD) { // special for octets if(particle.dataPtr()->iColour()==PDT::Colour8) { // octet -> octet octet if(cit->second.first->splittingFn()->colourStructure()==OctetOctetOctet) { type = ShowerPartnerType::QCDColourLine; Energy startingScale = angularOrdered ? particle.scales().QCD_c : particle.scales().QCD_c_noAO; newKin= cit->second.first-> generateNextTimeBranching(startingScale,cit->second.second, - particle.id()!=cit->first,0.5*enhance, + particle.id()!=cit->first,rho,0.5*enhance, particle.scales().Max_Q2); startingScale = angularOrdered ? particle.scales().QCD_ac : particle.scales().QCD_ac_noAO; ShoKinPtr newKin2 = cit->second.first-> generateNextTimeBranching(startingScale,cit->second.second, - particle.id()!=cit->first,0.5*enhance, + particle.id()!=cit->first,rho,0.5*enhance, particle.scales().Max_Q2); // pick the one with the highest scale if( ( newKin && newKin2 && newKin2->scale() > newKin->scale()) || (!newKin && newKin2) ) { newKin = newKin2; type = ShowerPartnerType::QCDAntiColourLine; } } // other g -> q qbar else { Energy startingScale = angularOrdered ? max(particle.scales().QCD_c , particle.scales().QCD_ac ) : max(particle.scales().QCD_c_noAO, particle.scales().QCD_ac_noAO); newKin= cit->second.first-> generateNextTimeBranching(startingScale, cit->second.second, - particle.id()!=cit->first,enhance, + particle.id()!=cit->first,rho,enhance, particle.scales().Max_Q2); type = UseRandom::rndbool() ? ShowerPartnerType::QCDColourLine : ShowerPartnerType::QCDAntiColourLine; } } // everything else q-> qg etc else { Energy startingScale; if(particle.hasColour()) { type = ShowerPartnerType::QCDColourLine; startingScale = angularOrdered ? particle.scales().QCD_c : particle.scales().QCD_c_noAO; } else { type = ShowerPartnerType::QCDAntiColourLine; startingScale = angularOrdered ? particle.scales().QCD_ac : particle.scales().QCD_ac_noAO; } newKin= cit->second.first-> generateNextTimeBranching(startingScale,cit->second.second, - particle.id()!=cit->first,enhance, + particle.id()!=cit->first,rho,enhance, particle.scales().Max_Q2); } } else if(cit->second.first->interactionType()==ShowerInteraction::EW) { type = ShowerPartnerType::EW; Energy startingScale = particle.scales().EW; newKin = cit->second.first-> generateNextTimeBranching(startingScale,cit->second.second, - particle.id()!=cit->first,enhance, + particle.id()!=cit->first,rho,enhance, particle.scales().Max_Q2); } // shouldn't be anything else else assert(false); // if no kinematics contine if(!newKin) continue; // select highest scale if( newKin->scale() > newQ ) { kinematics = newKin; newQ = newKin->scale(); ids = cit->second.second; sudakov = cit->second.first; partnerType = type; } } // return empty branching if nothing happened if(!kinematics) return Branching(ShoKinPtr(), IdList(),SudakovPtr(), ShowerPartnerType::Undefined); + // and generate phi + // if not hard generate phi + kinematics->phi(sudakov->generatePhiForward(particle,ids,kinematics,rho)); // and return it return Branching(kinematics, ids,sudakov,partnerType); } Branching SplittingGenerator:: chooseDecayBranching(ShowerParticle &particle, const ShowerParticle::EvolutionScales & stoppingScales, Energy minmass, double enhance, ShowerInteraction::Type interaction) const { + RhoDMatrix rho(particle.dataPtr()->iSpin()); Energy newQ = Constants::MaxEnergy; ShoKinPtr kinematics; SudakovPtr sudakov; ShowerPartnerType::Type partnerType(ShowerPartnerType::Undefined); IdList ids; // First, find the eventual branching, corresponding to the lowest scale. long index = abs(particle.data().id()); // if no branchings return empty branching struct if(_fbranchings.find(index) == _fbranchings.end()) return Branching(ShoKinPtr(), IdList(),SudakovPtr(),ShowerPartnerType::Undefined); // otherwise select branching for(BranchingList::const_iterator cit = _fbranchings.lower_bound(index); cit != _fbranchings.upper_bound(index); ++cit) { // check interaction doesn't change flavour if(cit->second.second[1]!=index&&cit->second.second[2]!=index) continue; // check either right interaction or doing both if(!checkInteraction(interaction,cit->second.first->interactionType())) continue; // whether or not this interaction should be angular ordered bool angularOrdered = cit->second.first->splittingFn()->angularOrdered(); ShoKinPtr newKin; ShowerPartnerType::Type type; // work out which starting scale we need if(cit->second.first->interactionType()==ShowerInteraction::QED) { type = ShowerPartnerType::QED; Energy stoppingScale = angularOrdered ? stoppingScales.QED : stoppingScales.QED_noAO; Energy startingScale = angularOrdered ? particle.scales().QED : particle.scales().QED_noAO; if(startingScale < stoppingScale ) { newKin = cit->second.first-> generateNextDecayBranching(startingScale,stoppingScale,minmass,cit->second.second, - particle.id()!=cit->first,enhance); + particle.id()!=cit->first,rho,enhance); } } else if(cit->second.first->interactionType()==ShowerInteraction::QCD) { // special for octets if(particle.dataPtr()->iColour()==PDT::Colour8) { // octet -> octet octet if(cit->second.first->splittingFn()->colourStructure()==OctetOctetOctet) { Energy stoppingColour = angularOrdered ? stoppingScales.QCD_c : stoppingScales.QCD_c_noAO; Energy stoppingAnti = angularOrdered ? stoppingScales.QCD_ac : stoppingScales.QCD_ac_noAO; Energy startingColour = angularOrdered ? particle.scales().QCD_c : particle.scales().QCD_c_noAO; Energy startingAnti = angularOrdered ? particle.scales().QCD_ac : particle.scales().QCD_ac_noAO; type = ShowerPartnerType::QCDColourLine; if(startingColoursecond.first-> generateNextDecayBranching(startingColour,stoppingColour,minmass, cit->second.second, - particle.id()!=cit->first,0.5*enhance); + particle.id()!=cit->first,rho,0.5*enhance); } ShoKinPtr newKin2; if(startingAntisecond.first-> generateNextDecayBranching(startingAnti,stoppingAnti,minmass, cit->second.second, - particle.id()!=cit->first,0.5*enhance); + particle.id()!=cit->first,rho,0.5*enhance); } // pick the one with the lowest scale if( (newKin&&newKin2&&newKin2->scale()scale()) || (!newKin&&newKin2) ) { newKin = newKin2; type = ShowerPartnerType::QCDAntiColourLine; } } // other else { assert(false); } } // everything else else { Energy startingScale,stoppingScale; if(particle.hasColour()) { type = ShowerPartnerType::QCDColourLine; stoppingScale = angularOrdered ? stoppingScales.QCD_c : stoppingScales.QCD_c_noAO; startingScale = angularOrdered ? particle.scales().QCD_c : particle.scales().QCD_c_noAO; } else { type = ShowerPartnerType::QCDAntiColourLine; stoppingScale = angularOrdered ? stoppingScales.QCD_ac : stoppingScales.QCD_ac_noAO; startingScale = angularOrdered ? particle.scales().QCD_ac : particle.scales().QCD_ac_noAO; } if(startingScale < stoppingScale ) { newKin = cit->second.first-> generateNextDecayBranching(startingScale,stoppingScale,minmass,cit->second.second, - particle.id()!=cit->first,enhance); + particle.id()!=cit->first,rho,enhance); } } } else if(cit->second.first->interactionType()==ShowerInteraction::EW) { type = ShowerPartnerType::EW; Energy stoppingScale = stoppingScales.EW; Energy startingScale = particle.scales().EW; if(startingScale < stoppingScale ) { newKin = cit->second.first-> generateNextDecayBranching(startingScale,stoppingScale,minmass,cit->second.second, - particle.id()!=cit->first,enhance); + particle.id()!=cit->first,rho,enhance); } } // shouldn't be anything else else assert(false); if(!newKin) continue; // select highest scale if(newKin->scale() < newQ ) { newQ = newKin->scale(); ids = cit->second.second; kinematics=newKin; sudakov=cit->second.first; partnerType = type; } } // return empty branching if nothing happened if(!kinematics) return Branching(ShoKinPtr(), IdList(),SudakovPtr(), ShowerPartnerType::Undefined); // and generate phi - kinematics->phi(sudakov->generatePhiDecay(particle,ids,kinematics)); + kinematics->phi(sudakov->generatePhiDecay(particle,ids,kinematics,rho)); // and return it return Branching(kinematics, ids,sudakov,partnerType); } Branching SplittingGenerator:: chooseBackwardBranching(ShowerParticle &particle,PPtr, double enhance, Ptr::transient_const_pointer beam, ShowerInteraction::Type type, tcPDFPtr pdf, Energy freeze) const { + RhoDMatrix rho = particle.extractRhoMatrix(false); Energy newQ=ZERO; ShoKinPtr kinematics=ShoKinPtr(); ShowerPartnerType::Type partnerType(ShowerPartnerType::Undefined); SudakovPtr sudakov; IdList ids; // First, find the eventual branching, corresponding to the highest scale. long index = abs(particle.id()); // if no possible branching return if(_bbranchings.find(index) == _bbranchings.end()) return Branching(ShoKinPtr(), IdList(),SudakovPtr(),ShowerPartnerType::Undefined); // otherwise select branching for(BranchingList::const_iterator cit = _bbranchings.lower_bound(index); cit != _bbranchings.upper_bound(index); ++cit ) { // check either right interaction or doing both if(!checkInteraction(type,cit->second.first->interactionType())) continue; // setup the PDF cit->second.first->setPDF(pdf,freeze); // whether or not this interaction should be angular ordered bool angularOrdered = cit->second.first->splittingFn()->angularOrdered(); ShoKinPtr newKin; ShowerPartnerType::Type type; if(cit->second.first->interactionType()==ShowerInteraction::QED) { type = ShowerPartnerType::QED; Energy startingScale = angularOrdered ? particle.scales().QED : particle.scales().QED_noAO; newKin=cit->second.first-> generateNextSpaceBranching(startingScale,cit->second.second, particle.x(), - particle.id()!=cit->first,enhance,beam); + particle.id()!=cit->first,rho,enhance,beam); } else if(cit->second.first->interactionType()==ShowerInteraction::QCD) { // special for octets if(particle.dataPtr()->iColour()==PDT::Colour8) { // octet -> octet octet if(cit->second.first->splittingFn()->colourStructure()==OctetOctetOctet) { type = ShowerPartnerType::QCDColourLine; Energy startingScale = angularOrdered ? particle.scales().QCD_c : particle.scales().QCD_c_noAO; newKin = cit->second.first-> generateNextSpaceBranching(startingScale,cit->second.second, particle.x(), - particle.id()!=cit->first,0.5*enhance,beam); + particle.id()!=cit->first,rho,0.5*enhance,beam); startingScale = angularOrdered ? particle.scales().QCD_ac : particle.scales().QCD_ac_noAO; ShoKinPtr newKin2 = cit->second.first-> generateNextSpaceBranching(startingScale,cit->second.second, particle.x(), - particle.id()!=cit->first,0.5*enhance,beam); + particle.id()!=cit->first,rho,0.5*enhance,beam); // pick the one with the highest scale if( (newKin&&newKin2&&newKin2->scale()>newKin->scale()) || (!newKin&&newKin2) ) { newKin = newKin2; type = ShowerPartnerType::QCDAntiColourLine; } } else { Energy startingScale = angularOrdered ? max(particle.scales().QCD_c , particle.scales().QCD_ac ) : max(particle.scales().QCD_c_noAO, particle.scales().QCD_ac_noAO); type = UseRandom::rndbool() ? ShowerPartnerType::QCDColourLine : ShowerPartnerType::QCDAntiColourLine; newKin=cit->second.first-> generateNextSpaceBranching(startingScale,cit->second.second, particle.x(), - particle.id()!=cit->first,enhance,beam); + particle.id()!=cit->first,rho,enhance,beam); } } // everything else else { Energy startingScale; if(particle.hasColour()) { type = ShowerPartnerType::QCDColourLine; startingScale = angularOrdered ? particle.scales().QCD_c : particle.scales().QCD_c_noAO; } else { type = ShowerPartnerType::QCDAntiColourLine; startingScale = angularOrdered ? particle.scales().QCD_ac : particle.scales().QCD_ac_noAO; } newKin=cit->second.first-> generateNextSpaceBranching(startingScale,cit->second.second, particle.x(), - particle.id()!=cit->first,enhance,beam); + particle.id()!=cit->first,rho,enhance,beam); } } else if(cit->second.first->interactionType()==ShowerInteraction::EW) { type = ShowerPartnerType::EW; Energy startingScale = particle.scales().EW; newKin=cit->second.first-> generateNextSpaceBranching(startingScale,cit->second.second, particle.x(), - particle.id()!=cit->first,enhance,beam); + particle.id()!=cit->first,rho,enhance,beam); } // shouldn't be anything else else assert(false); // if no kinematics contine if(!newKin) continue; // select highest scale if(newKin->scale() > newQ) { newQ = newKin->scale(); kinematics=newKin; ids = cit->second.second; sudakov=cit->second.first; partnerType = type; } } // return empty branching if nothing happened if(!kinematics) return Branching(ShoKinPtr(), IdList(),SudakovPtr(), ShowerPartnerType::Undefined); // initialize the ShowerKinematics // and return it // and generate phi - kinematics->phi(sudakov->generatePhiBackward(particle,ids,kinematics)); + kinematics->phi(sudakov->generatePhiBackward(particle,ids,kinematics,rho)); // return the answer return Branching(kinematics, ids,sudakov,partnerType); } void SplittingGenerator::rebind(const TranslationMap & trans) { BranchingList::iterator cit; for(cit=_fbranchings.begin();cit!=_fbranchings.end();++cit) {(cit->second).first=trans.translate((cit->second).first);} for(cit=_bbranchings.begin();cit!=_bbranchings.end();++cit) {(cit->second).first=trans.translate((cit->second).first);} Interfaced::rebind(trans); } IVector SplittingGenerator::getReferences() { IVector ret = Interfaced::getReferences(); BranchingList::iterator cit; for(cit=_fbranchings.begin();cit!=_fbranchings.end();++cit) {ret.push_back((cit->second).first);} for(cit=_bbranchings.begin();cit!=_bbranchings.end();++cit) {ret.push_back((cit->second).first);} return ret; } void SplittingGenerator::factorizationScaleFactor(double f) { BranchingList::iterator cit; for(cit=_fbranchings.begin();cit!=_fbranchings.end();++cit) {(cit->second).first->factorizationScaleFactor(f);} for(cit=_bbranchings.begin();cit!=_bbranchings.end();++cit) {(cit->second).first->factorizationScaleFactor(f);} } void SplittingGenerator::renormalizationScaleFactor(double f) { BranchingList::iterator cit; for(cit=_fbranchings.begin();cit!=_fbranchings.end();++cit) {(cit->second).first->renormalizationScaleFactor(f);} for(cit=_bbranchings.begin();cit!=_bbranchings.end();++cit) {(cit->second).first->renormalizationScaleFactor(f);} } diff --git a/Shower/SplittingFunctions/ZeroZeroOneSplitFn.cc b/Shower/SplittingFunctions/ZeroZeroOneSplitFn.cc --- a/Shower/SplittingFunctions/ZeroZeroOneSplitFn.cc +++ b/Shower/SplittingFunctions/ZeroZeroOneSplitFn.cc @@ -1,120 +1,120 @@ // -*- C++ -*- // // PhitoPhiGSplitFn.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 ZeroZeroOneSplitFn class. // #include "ZeroZeroOneSplitFn.h" #include "ThePEG/PDT/ParticleData.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Utilities/DescribeClass.h" #include "Herwig/Decay/TwoBodyDecayMatrixElement.h" using namespace Herwig; DescribeNoPIOClass describeZeroZeroOneSplitFn ("Herwig::ZeroZeroOneSplitFn","HwShower.so"); void ZeroZeroOneSplitFn::Init() { static ClassDocumentation documentation ("The ZeroZeroOneSplitFn class implements the splitting function for the " "radiation of a gluon by a scalar coloured particle"); } double ZeroZeroOneSplitFn::P(const double z, const Energy2 t, - const IdList &ids, const bool mass) const { + const IdList &ids, const bool mass, const RhoDMatrix &) const { double val = z/(1.-z); if(mass) { Energy m = getParticleData(ids[0])->mass(); val-= sqr(m)/t; } return 2.*colourFactor(ids)*val; } double ZeroZeroOneSplitFn::overestimateP(const double z, const IdList &ids) const { return 2.*colourFactor(ids)/(1.-z); } double ZeroZeroOneSplitFn::ratioP(const double z, const Energy2 t, - const IdList &ids,const bool mass) const { + const IdList &ids,const bool mass, const RhoDMatrix &) const { double val = z; if(mass) { Energy m = getParticleData(ids[0])->mass(); val-=sqr(m)*(1.-z)/t; } return val; } double ZeroZeroOneSplitFn::integOverP(const double z, const IdList & ids, unsigned int PDFfactor) const { switch(PDFfactor) { case 0: return -2.*colourFactor(ids)*log(1.-z); case 1: case 2: case 3: default: throw Exception() << "ZeroZeroOneSplitFn::integOverP() invalid PDFfactor = " << PDFfactor << Exception::runerror; } } double ZeroZeroOneSplitFn::invIntegOverP(const double r, const IdList & ids, unsigned int PDFfactor) const { switch(PDFfactor) { case 0: return 1. - exp(- 0.5*r/colourFactor(ids)); case 1: case 2: case 3: default: throw Exception() << "ZeroZeroOneSplitFn::integOverP() invalid PDFfactor = " << PDFfactor << Exception::runerror; } } bool ZeroZeroOneSplitFn::accept(const IdList &ids) const { if(ids.size()!=3) return false; if(ids[0]!=ids[1]) return false; tcPDPtr q=getParticleData(ids[0]); tcPDPtr g=getParticleData(ids[2]); if(q->iSpin()!=PDT::Spin0 || g->iSpin()!=PDT::Spin1) return false; return checkColours(ids); } vector > ZeroZeroOneSplitFn::generatePhiForward(const double, const Energy2, const IdList &, const RhoDMatrix &) { // scalar so no dependence return vector >(1,make_pair(0,1.)); } vector > ZeroZeroOneSplitFn::generatePhiBackward(const double, const Energy2, const IdList &, const RhoDMatrix &) { // scalar so no dependence assert(false); return vector >(1,make_pair(0,1.)); } DecayMEPtr ZeroZeroOneSplitFn::matrixElement(const double z, const Energy2 t, const IdList & ids, const double phi, bool timeLike) { // calculate the kernal DecayMEPtr kernal(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin0,PDT::Spin0,PDT::Spin1))); Energy m = timeLike ? getParticleData(ids[0])->mass() : ZERO; (*kernal)(0,0,0) = -exp(Complex(0.,1.)*phi)*sqrt(1.-(1.-z)*sqr(m)/z/t)*sqrt(z/(1.-z)); (*kernal)(0,0,2) = -conj((*kernal)(0,0,0)); return kernal; } diff --git a/Shower/SplittingFunctions/ZeroZeroOneSplitFn.h b/Shower/SplittingFunctions/ZeroZeroOneSplitFn.h --- a/Shower/SplittingFunctions/ZeroZeroOneSplitFn.h +++ b/Shower/SplittingFunctions/ZeroZeroOneSplitFn.h @@ -1,191 +1,193 @@ // -*- C++ -*- // // ZeroZeroOneSplitFn.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_ZeroZeroOneSplitFn_H #define HERWIG_ZeroZeroOneSplitFn_H // // This is the declaration of the ZeroZeroOneSplitFn class. // #include "SplittingFunction.h" namespace Herwig { using namespace ThePEG; /** \ingroup Shower * This class provides the concrete implementation of the exact leading-order * splitting function for \f$\phi\to \phi g\f$. * * In this case the splitting function is given by * \f[P(z,t) = 2C\left(\frac{z}{1-z}-\frac{m^2_\phi}{t}\right),\f] * where \f$C\f$ is the corresponding colour factor. * Our choice for the overestimate is * \f[P_{\rm over}(z) = \frac{2C}{1-z},\f] * therefore the integral is * \f[\int P_{\rm over}(z) {\rm d}z = -2C\ln(1-z),\f] * and its inverse is * \f[1-\exp\left(\frac{r}{2C}\right).\f] * * @see \ref ZeroZeroOneSplitFnInterfaces "The interfaces" * defined for ZeroZeroOneSplitFn. */ class ZeroZeroOneSplitFn: public SplittingFunction { public: /** * The default constructor. */ ZeroZeroOneSplitFn() : SplittingFunction(1) {} /** * Concrete implementation of the method to determine whether this splitting * function can be used for a given set of particles. * @param ids The PDG codes for the particles in the splitting. */ virtual bool accept(const IdList & ids) const; /** * Methods to return the splitting function. */ //@{ /** * The concrete implementation of the splitting function, \f$P\f$. * @param z The energy fraction. * @param t The scale. * @param ids The PDG codes for the particles in the splitting. * @param mass Whether or not to include the mass dependent terms + * @param rho The spin density matrix */ virtual double P(const double z, const Energy2 t, const IdList & ids, - bool mass) const; + bool mass, const RhoDMatrix & rho) const; /** * The concrete implementation of the overestimate of the splitting function, * \f$P_{\rm over}\f$. * @param z The energy fraction. * @param ids The PDG codes for the particles in the splitting. */ virtual double overestimateP(const double z, const IdList & ids) const; /** * The concrete implementation of the * the ratio of the splitting function to the overestimate, i.e. * \f$P(z,\tilde{q}^2)/P_{\rm over}(z)\f$. * @param z The energy fraction. * @param t The scale. * @param ids The PDG codes for the particles in the splitting. * @param mass Whether or not to include the mass dependent terms + * @param rho The spin density matrix */ virtual double ratioP(const double z, const Energy2 t, const IdList & ids, - bool mass) const; + bool mass, const RhoDMatrix & rho) const; /** * The concrete implementation of the indefinite integral of the * overestimated splitting function, \f$P_{\rm over}\f$. * @param z The energy fraction. * @param ids The PDG codes for the particles in the splitting. * @param PDFfactor Which additional factor to include for the PDF * 0 is no additional factor, * 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$ */ virtual double integOverP(const double z, const IdList & ids, unsigned int PDFfactor=0) const; /** * The concrete implementation of the inverse of the indefinite integral. * @param r Value of the splitting function to be inverted * @param ids The PDG codes for the particles in the splitting. * @param PDFfactor Which additional factor to include for the PDF * 0 is no additional factor, * 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$ */ virtual double invIntegOverP(const double r, const IdList & ids, unsigned int PDFfactor=0) const; //@} /** * Method to calculate the azimuthal angle for forward evolution * @param particle The particle which is branching * @param showerkin The ShowerKinematics object * @param z The energy fraction * @param t The scale \f$t=2p_j\cdot p_k\f$. * @param ids The PDG codes for the particles in the splitting. * @param The azimuthal angle, \f$\phi\f$. * @return The weight */ virtual vector > generatePhiForward(const double z, const Energy2 t, const IdList & ids, const RhoDMatrix &); /** * Method to calculate the azimuthal angle for backward * Shouldn't be needed and NOT IMPLEMENTED * @param particle The particle which is branching * @param showerkin The ShowerKinematics object * @param z The energy fraction * @param t The scale \f$t=2p_j\cdot p_k\f$. * @param ids The PDG codes for the particles in the splitting. * @param The azimuthal angle, \f$\phi\f$. * @return The weight */ virtual vector > generatePhiBackward(const double z, const Energy2 t, const IdList & ids, const RhoDMatrix &); /** * Calculate the matrix element for the splitting * @param particle The particle which is branching * @param showerkin The ShowerKinematics object * @param z The energy fraction * @param t The scale \f$t=2p_j\cdot p_k\f$. * @param ids The PDG codes for the particles in the splitting. * @param The azimuthal angle, \f$\phi\f$. */ virtual DecayMEPtr matrixElement(const double z, const Energy2 t, const IdList & ids, const double phi, bool timeLike); public: /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const {return new_ptr(*this);} /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const {return new_ptr(*this);} //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ ZeroZeroOneSplitFn & operator=(const ZeroZeroOneSplitFn &); }; } #endif /* HERWIG_ZeroZeroOneSplitFn_H */