diff --git a/Shower/QTilde/Base/SudakovFormFactor.cc b/Shower/QTilde/Base/SudakovFormFactor.cc --- a/Shower/QTilde/Base/SudakovFormFactor.cc +++ b/Shower/QTilde/Base/SudakovFormFactor.cc @@ -1,1395 +1,1384 @@ // -*- C++ -*- // // SudakovFormFactor.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the SudakovFormFactor class. // #include "SudakovFormFactor.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Parameter.h" #include "ShowerKinematics.h" #include "ShowerParticle.h" #include "ThePEG/Utilities/DescribeClass.h" #include "Herwig/Shower/QTilde/QTildeShowerHandler.h" #include "Herwig/Shower/QTilde/Kinematics/FS_QTildeShowerKinematics1to2.h" #include "Herwig/Shower/QTilde/Kinematics/IS_QTildeShowerKinematics1to2.h" #include "Herwig/Shower/QTilde/Kinematics/Decay_QTildeShowerKinematics1to2.h" #include using std::array; using namespace Herwig; DescribeClass describeSudakovFormFactor ("Herwig::SudakovFormFactor",""); void SudakovFormFactor::persistentOutput(PersistentOStream & os) const { os << splittingFn_ << alpha_ << pdfmax_ << particles_ << pdffactor_ << a_ << b_ << ounit(c_,GeV) << ounit(kinCutoffScale_,GeV) << cutOffOption_ << ounit(vgcut_,GeV) << ounit(vqcut_,GeV) << ounit(pTmin_,GeV) << ounit(pT2min_,GeV2); } void SudakovFormFactor::persistentInput(PersistentIStream & is, int) { is >> splittingFn_ >> alpha_ >> pdfmax_ >> particles_ >> pdffactor_ >> a_ >> b_ >> iunit(c_,GeV) >> iunit(kinCutoffScale_,GeV) >> cutOffOption_ >> iunit(vgcut_,GeV) >> iunit(vqcut_,GeV) >> iunit(pTmin_,GeV) >> iunit(pT2min_,GeV2); } void SudakovFormFactor::Init() { static ClassDocumentation documentation ("The SudakovFormFactor class is the base class for the implementation of Sudakov" " form factors in Herwig"); static Reference interfaceSplittingFunction("SplittingFunction", "A reference to the SplittingFunction object", &Herwig::SudakovFormFactor::splittingFn_, false, false, true, false); static Reference interfaceAlpha("Alpha", "A reference to the Alpha object", &Herwig::SudakovFormFactor::alpha_, false, false, true, false); static Parameter interfacePDFmax ("PDFmax", "Maximum value of PDF weight. ", &SudakovFormFactor::pdfmax_, 35.0, 1.0, 1000000.0, false, false, Interface::limited); static Switch interfacePDFFactor ("PDFFactor", "Include additional factors in the overestimate for the PDFs", &SudakovFormFactor::pdffactor_, 0, false, false); static SwitchOption interfacePDFFactorNo (interfacePDFFactor, "No", "Don't include any factors", 0); static SwitchOption interfacePDFFactorOverZ (interfacePDFFactor, "OverZ", "Include an additional factor of 1/z", 1); static SwitchOption interfacePDFFactorOverOneMinusZ (interfacePDFFactor, "OverOneMinusZ", "Include an additional factor of 1/(1-z)", 2); static SwitchOption interfacePDFFactorOverZOneMinusZ (interfacePDFFactor, "OverZOneMinusZ", "Include an additional factor of 1/z/(1-z)", 3); static SwitchOption interfacePDFFactorOverRootZ (interfacePDFFactor, "OverRootZ", "Include an additional factor of 1/sqrt(z)", 4); static SwitchOption interfacePDFFactorRootZ (interfacePDFFactor, "RootZ", "Include an additional factor of sqrt(z)", 5); static Switch interfaceCutOffOption ("CutOffOption", "The type of cut-off to use to end the shower", &SudakovFormFactor::cutOffOption_, 0, false, false); static SwitchOption interfaceCutOffOptionDefault (interfaceCutOffOption, "Default", "Use the standard Herwig cut-off on virtualities with the minimum" " virtuality depending on the mass of the branching particle", 0); static SwitchOption interfaceCutOffOptionFORTRAN (interfaceCutOffOption, "FORTRAN", "Use a FORTRAN-like cut-off on virtualities", 1); static SwitchOption interfaceCutOffOptionpT (interfaceCutOffOption, "pT", "Use a cut on the minimum allowed pT", 2); static Parameter interfaceaParameter ("aParameter", "The a parameter for the kinematic cut-off", &SudakovFormFactor::a_, 0.3, -10.0, 10.0, false, false, Interface::limited); static Parameter interfacebParameter ("bParameter", "The b parameter for the kinematic cut-off", &SudakovFormFactor::b_, 2.3, -10.0, 10.0, false, false, Interface::limited); static Parameter interfacecParameter ("cParameter", "The c parameter for the kinematic cut-off", &SudakovFormFactor::c_, GeV, 0.3*GeV, 0.1*GeV, 10.0*GeV, false, false, Interface::limited); static Parameter interfaceKinScale ("cutoffKinScale", "kinematic cutoff scale for the parton shower phase" " space (unit [GeV])", &SudakovFormFactor::kinCutoffScale_, GeV, 2.3*GeV, 0.001*GeV, 10.0*GeV,false,false,false); static Parameter interfaceGluonVirtualityCut ("GluonVirtualityCut", "For the FORTRAN cut-off option the minimum virtuality of the gluon", &SudakovFormFactor::vgcut_, GeV, 0.85*GeV, 0.1*GeV, 10.0*GeV, false, false, Interface::limited); static Parameter interfaceQuarkVirtualityCut ("QuarkVirtualityCut", "For the FORTRAN cut-off option the minimum virtuality added to" " the mass for particles other than the gluon", &SudakovFormFactor::vqcut_, GeV, 0.85*GeV, 0.1*GeV, 10.0*GeV, false, false, Interface::limited); static Parameter interfacepTmin ("pTmin", "The minimum pT if using a cut-off on the pT", &SudakovFormFactor::pTmin_, GeV, 1.0*GeV, ZERO, 10.0*GeV, false, false, Interface::limited); } bool SudakovFormFactor::alphaSVeto(Energy2 pt2) const { double ratio=alphaSVetoRatio(pt2,1.); return UseRandom::rnd() > ratio; } double SudakovFormFactor::alphaSVetoRatio(Energy2 pt2, double factor) const { factor *= ShowerHandler::currentHandler()->renormalizationScaleFactor(); - return ThePEG::Math::powi(alpha_->ratio(pt2, factor), - splittingFn_->interactionOrder()); + return alpha_->ratio(pt2, factor); } bool SudakovFormFactor::PDFVeto(const Energy2 t, const double x, const tcPDPtr parton0, const tcPDPtr parton1, Ptr::transient_const_pointer beam) const { double ratio=PDFVetoRatio(t,x,parton0,parton1,beam,1.); return UseRandom::rnd() > ratio; } double SudakovFormFactor::PDFVetoRatio(const Energy2 t, const double x, const tcPDPtr parton0, const tcPDPtr parton1, Ptr::transient_const_pointer beam,double factor) const { assert(pdf_); Energy2 theScale = t * sqr(ShowerHandler::currentHandler()->factorizationScaleFactor()*factor); if (theScale < sqr(freeze_)) theScale = sqr(freeze_); double newpdf(0.0), oldpdf(0.0); newpdf=pdf_->xfx(beam,parton0,theScale,x/z()); oldpdf=pdf_->xfx(beam,parton1,theScale,x); if(newpdf<=0.) return 0.; if(oldpdf<=0.) return 1.; double ratio = newpdf/oldpdf; double maxpdf = pdfmax_; switch (pdffactor_) { case 0: break; case 1: maxpdf /= z(); break; case 2: maxpdf /= 1.-z(); break; case 3: maxpdf /= (z()*(1.-z())); break; case 4: maxpdf /= sqrt(z()); break; case 5: maxpdf *= sqrt(z()); break; default : throw Exception() << "SudakovFormFactor::PDFVetoRatio invalid PDFfactor = " << pdffactor_ << Exception::runerror; } if (ratio > maxpdf) { generator()->log() << "PDFVeto warning: Ratio > " << name() << ":PDFmax (by a factor of " << ratio/maxpdf <<") for " << parton0->PDGName() << " to " << parton1->PDGName() << "\n"; } return ratio/maxpdf ; } void SudakovFormFactor::addSplitting(const IdList & in) { bool add=true; for(unsigned int ix=0;ix::iterator it=particles_.begin(); it!=particles_.end();++it) { if(it->size()==in.size()) { bool match=true; for(unsigned int iy=0;iy::iterator itemp=it; --itemp; particles_.erase(it); it = itemp; } } } } Energy2 SudakovFormFactor::guesst(Energy2 t1,unsigned int iopt, const IdList &ids, double enhance,bool ident, double detune) const { unsigned int pdfopt = iopt!=1 ? 0 : pdffactor_; double c = 1./((splittingFn_->integOverP(zlimits_.second,ids,pdfopt) - splittingFn_->integOverP(zlimits_.first ,ids,pdfopt))* alpha_->overestimateValue()/Constants::twopi*enhance*detune); assert(iopt<=2); if(iopt==1) { c/=pdfmax_; //symmetry of FS gluon splitting if(ident) c*=0.5; } else if(iopt==2) c*=-1.; - if(splittingFn_->interactionOrder()==1) { - double r = UseRandom::rnd(); - if(iopt!=2 || c*log(r)interactionOrder()-1); - c/=Math::powi(alpha_->overestimateValue()/Constants::twopi,nm); - return t1 / pow (1. - nm*c*log(UseRandom::rnd()) - * Math::powi(t1*UnitRemoval::InvE2,nm) - ,1./double(nm)); - } + else + return Constants::MaxEnergy2; } double SudakovFormFactor::guessz (unsigned int iopt, const IdList &ids) const { unsigned int pdfopt = iopt!=1 ? 0 : pdffactor_; double lower = splittingFn_->integOverP(zlimits_.first,ids,pdfopt); return splittingFn_->invIntegOverP (lower + UseRandom::rnd()*(splittingFn_->integOverP(zlimits_.second,ids,pdfopt) - lower),ids,pdfopt); } void SudakovFormFactor::doinit() { Interfaced::doinit(); pT2min_ = cutOffOption()==2 ? sqr(pTmin_) : ZERO; } const vector & SudakovFormFactor::virtualMasses(const IdList & ids) { static vector output; output.clear(); if(cutOffOption() == 0) { for(unsigned int ix=0;ixmass()); Energy kinCutoff= kinematicCutOff(kinScale(),*std::max_element(output.begin(),output.end())); for(unsigned int ix=0;ixmass()); output.back() += ids[ix]->id()==ParticleID::g ? vgCut() : vqCut(); } } else if(cutOffOption() == 2) { for(unsigned int ix=0;ixmass()); } else { throw Exception() << "Unknown option for the cut-off" << " in SudakovFormFactor::virtualMasses()" << Exception::runerror; } return output; } bool SudakovFormFactor::guessTimeLike(Energy2 &t,Energy2 tmin,double enhance, double detune) { 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],detune); 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],detune); 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]->id()!=ParticleID::g && ids_[0]->id()!=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); // no shower variations to calculate if(ShowerHandler::currentHandler()->showerVariations().empty()){ // Without variations do the usual Veto algorithm // No need for more if-statements in this loop. do { if(!guessTimeLike(t,tmin,enhance,detuning)) break; } while(PSVeto(t,maxQ2) || SplittingFnVeto(z()*(1.-z())*t,ids,true,rho,detuning) || alphaSVeto(splittingFn()->pTScale() ? sqr(z()*(1.-z()))*t : z()*(1.-z())*t)); } else { bool alphaRew(true),PSRew(true),SplitRew(true); do { if(!guessTimeLike(t,tmin,enhance,detuning)) break; PSRew=PSVeto(t,maxQ2); if (PSRew) continue; SplitRew=SplittingFnVeto(z()*(1.-z())*t,ids,true,rho,detuning); alphaRew=alphaSVeto(splittingFn()->pTScale() ? sqr(z()*(1.-z()))*t : z()*(1.-z())*t); double factor=alphaSVetoRatio(splittingFn()->pTScale() ? sqr(z()*(1.-z()))*t : z()*(1.-z())*t,1.)* SplittingFnVetoRatio(z()*(1.-z())*t,ids,true,rho,detuning); tShowerHandlerPtr ch = ShowerHandler::currentHandler(); if( !(SplitRew || alphaRew) ) { //Emission q_ = t > ZERO ? Energy(sqrt(t)) : -1.*MeV; if (q_ <= ZERO) break; } for ( map::const_iterator var = ch->showerVariations().begin(); var != ch->showerVariations().end(); ++var ) { if ( ( ch->firstInteraction() && var->second.firstInteraction ) || ( !ch->firstInteraction() && var->second.secondaryInteractions ) ) { double newfactor = alphaSVetoRatio(splittingFn()->pTScale() ? sqr(z()*(1.-z()))*t : z()*(1.-z())*t,var->second.renormalizationScaleFactor) * SplittingFnVetoRatio(z()*(1.-z())*t,ids,true,rho,detuning); double varied; if ( SplitRew || alphaRew ) { // No Emission varied = (1. - newfactor) / (1. - factor); } else { // Emission varied = newfactor / factor; } map::iterator wi = ch->currentWeights().find(var->first); if ( wi != ch->currentWeights().end() ) wi->second *= varied; else { assert(false); //ch->currentWeights()[var->first] = varied; } } } } while(PSRew || SplitRew || alphaRew); } q_ = t > ZERO ? Energy(sqrt(t)) : -1.*MeV; if(q_ < ZERO) return ShoKinPtr(); // return the ShowerKinematics object return createFinalStateBranching(q_,z(),phi(),pT()); } ShoKinPtr SudakovFormFactor:: generateNextSpaceBranching(const Energy startingQ, const IdList &ids, double x, const RhoDMatrix & rho, double enhance, Ptr::transient_const_pointer beam, double detuning) { // 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); // check max > min if(tmax<=tmin) return ShoKinPtr(); // calculate next value of t using veto algorithm Energy2 t(tmax),pt2(ZERO); // no shower variations if(ShowerHandler::currentHandler()->showerVariations().empty()){ // Without variations do the usual Veto algorithm // No need for more if-statements in this loop. do { if(!guessSpaceLike(t,tmin,x,enhance,detuning)) break; pt2=sqr(1.-z())*t-z()*masssquared_[2]; } while(pt2 < pT2min()|| z() > zlimits_.second|| SplittingFnVeto((1.-z())*t/z(),ids,false,rho,detuning)|| alphaSVeto(splittingFn()->pTScale() ? sqr(1.-z())*t : (1.-z())*t)|| PDFVeto(t,x,ids[0],ids[1],beam)); } // shower variations else { bool alphaRew(true),PDFRew(true),ptRew(true),zRew(true),SplitRew(true); do { if(!guessSpaceLike(t,tmin,x,enhance,detuning)) break; pt2=sqr(1.-z())*t-z()*masssquared_[2]; ptRew=pt2 < pT2min(); zRew=z() > zlimits_.second; if (ptRew||zRew) continue; SplitRew=SplittingFnVeto((1.-z())*t/z(),ids,false,rho,detuning); alphaRew=alphaSVeto(splittingFn()->pTScale() ? sqr(1.-z())*t : (1.-z())*t); PDFRew=PDFVeto(t,x,ids[0],ids[1],beam); double factor=PDFVetoRatio(t,x,ids[0],ids[1],beam,1.)* alphaSVetoRatio(splittingFn()->pTScale() ? sqr(1.-z())*t : (1.-z())*t,1.)* SplittingFnVetoRatio((1.-z())*t/z(),ids,false,rho,detuning); tShowerHandlerPtr ch = ShowerHandler::currentHandler(); if( !(PDFRew || SplitRew || alphaRew) ) { //Emission q_ = t > ZERO ? Energy(sqrt(t)) : -1.*MeV; if (q_ <= ZERO) break; } for ( map::const_iterator var = ch->showerVariations().begin(); var != ch->showerVariations().end(); ++var ) { if ( ( ch->firstInteraction() && var->second.firstInteraction ) || ( !ch->firstInteraction() && var->second.secondaryInteractions ) ) { double newfactor = PDFVetoRatio(t,x,ids[0],ids[1],beam,var->second.factorizationScaleFactor)* alphaSVetoRatio(splittingFn()->pTScale() ? sqr(1.-z())*t : (1.-z())*t,var->second.renormalizationScaleFactor) *SplittingFnVetoRatio((1.-z())*t/z(),ids,false,rho,detuning); double varied; if( PDFRew || SplitRew || alphaRew) { // No Emission varied = (1. - newfactor) / (1. - factor); } else { // Emission varied = newfactor / factor; } map::iterator wi = ch->currentWeights().find(var->first); if ( wi != ch->currentWeights().end() ) wi->second *= varied; else { assert(false); //ch->currentWeights()[var->first] = varied; } } } } while( PDFRew || SplitRew || alphaRew); } 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 SudakovFormFactor::initialize(const IdList & ids, Energy2 & tmin) { ids_=ids; tmin = cutOffOption() != 2 ? ZERO : 4.*pT2min(); masses_ = virtualMasses(ids); masssquared_.clear(); for(unsigned int ix=0;ix0) tmin=max(masssquared_[ix],tmin); } } ShoKinPtr SudakovFormFactor::generateNextDecayBranching(const Energy startingScale, const Energy stoppingScale, const Energy minmass, const IdList &ids, const RhoDMatrix & rho, double enhance, double detuning) { // 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); 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,detuning)) break; pt2 = sqr(1.-z())*(t-masssquared_[0])-z()*masssquared_[2]; } while(SplittingFnVeto((1.-z())*t/z(),ids,true,rho,detuning)|| alphaSVeto(splittingFn()->pTScale() ? 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 SudakovFormFactor::guessDecay(Energy2 &t,Energy2 tmax, Energy minmass, double enhance, double detune) { // previous scale Energy2 told = t; // overestimated limits on z if(tmaxtmax||zlimits_.secondid()==ParticleID::g||ids_[0]->id()==ParticleID::gamma) { // no emission possible if(t<16.*(masssquared_[1]+pT2min())) { t=-1.*GeV2; return false; } // overestimate of the limits zlimits_.first = 0.5*(1.-sqrt(1.-4.*sqrt((masssquared_[1]+pT2min())/t))); zlimits_.second = 1.-zlimits_.first; } // special case for radiated particle is gluon else if(ids_[2]->id()==ParticleID::g||ids_[2]->id()==ParticleID::gamma) { zlimits_.first = sqrt((masssquared_[1]+pT2min())/t); zlimits_.second = 1.-sqrt((masssquared_[2]+pT2min())/t); } else if(ids_[1]->id()==ParticleID::g||ids_[1]->id()==ParticleID::gamma) { zlimits_.second = sqrt((masssquared_[2]+pT2min())/t); zlimits_.first = 1.-sqrt((masssquared_[1]+pT2min())/t); } else { zlimits_.first = (masssquared_[1]+pT2min())/t; zlimits_.second = 1.-(masssquared_[2]+pT2min())/t; } if(zlimits_.first>=zlimits_.second) { t=-1.*GeV2; return false; } return true; } bool SudakovFormFactor::computeSpaceLikeLimits(Energy2 & t, double x) { if (t < 1e-20 * GeV2) { t=-1.*GeV2; return false; } // compute the limits zlimits_.first = x; double yy = 1.+0.5*masssquared_[2]/t; zlimits_.second = yy - sqrt(sqr(yy)-1.+pT2min()/t); // return false if lower>upper if(zlimits_.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); if(denom==ZERO || denom2==0) return make_pair(phi0,phi0+Constants::pi); 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 SudakovFormFactor::generatePhiForward(ShowerParticle & particle, const IdList & ids, ShoKinPtr kinematics, const RhoDMatrix & rho) { // no correlations, return flat phi if(! dynamic_ptr_cast(ShowerHandler::currentHandler())->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]->id()==ParticleID::g || ids[1]->id()==ParticleID::gamma, ids[2]->id()==ParticleID::g || ids[2]->id()==ParticleID::gamma }; array pjk; array Ek; Energy Ei,Ej; Energy2 m12(ZERO),m22(ZERO); InvEnergy2 aziMax(ZERO); bool softAllowed = dynamic_ptr_cast(ShowerHandler::currentHandler())->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(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { aziMax = -m12/sqr(pik) -m22/sqr(pjk[0]+mag) +2.*pipj/pik/(pjk[0]-mag); } else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->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(dynamic_ptr_cast(ShowerHandler::currentHandler())->spinCorrelations()) { // 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]->id() << " " << ids[1]->id() << " " << ids[2]->id() << " " << " " << 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(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { aziWgt = (-m12/sqr(pik) -m22/sqr(dot) +2.*pipj/pik/dot)/aziMax; } else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->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]->id() << " " << ids[1]->id() << " " << ids[2]->id() << " " << " " << 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 SudakovFormFactor::generatePhiBackward(ShowerParticle & particle, const IdList & ids, ShoKinPtr kinematics, const RhoDMatrix & rho) { // no correlations, return flat phi if(! dynamic_ptr_cast(ShowerHandler::currentHandler())->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 = dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations() && (ids[2]->id()==ParticleID::g || ids[2]->id()==ParticleID::gamma); Energy2 pipj,pik,m12(ZERO),m22(ZERO); array pjk; 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(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { aziMax = -m12/sqr(pik) -m22/sqr(pjk[0]+mag) +2.*pipj/pik/(pjk[0]-mag); } else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->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(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==0); } } // if spin correlations vector > wgts; if(dynamic_ptr_cast(ShowerHandler::currentHandler())->spinCorrelations()) { // 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]->id() << " " << ids[1]->id() << " " << ids[2]->id() << " " << " " << 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(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { aziWgt = (-m12/sqr(pik) -m22/sqr(dot) +2.*pipj/pik/dot)/aziMax; } else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->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]->id() << " " << ids[1]->id() << " " << ids[2]->id() << " " << " " << 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 SudakovFormFactor::generatePhiDecay(ShowerParticle & particle, const IdList & ids, ShoKinPtr kinematics, const RhoDMatrix &) { // only soft correlations in this case // no correlations, return flat phi if( !(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations() && (ids[2]->id()==ParticleID::g || ids[2]->id()==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 array pjk; 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; array Ek; Energy Ei,Ej; if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { aziMax = -m12/sqr(pik) -m22/sqr(pjk[0]+mag) +2.*pipj/pik/(pjk[0]-mag); } else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->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(dynamic_ptr_cast(ShowerHandler::currentHandler())->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(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { wgt = (-m12/sqr(pik) -m22/sqr(dot) +2.*pipj/pik/dot)/aziMax; } else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->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]->id() << " " << ids[1]->id() << " " << ids[2]->id() << " " << " " << 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 SudakovFormFactor::calculateScale(double zin, Energy pt, IdList ids, unsigned int iopt) { Energy2 tmin; initialize(ids,tmin); // final-state branching if(iopt==0) { Energy2 scale=(sqr(pt)+masssquared_[1]*(1.-zin)+masssquared_[2]*zin); if(ids[0]->id()!=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 SudakovFormFactor::calculateScale() " << "iopt = " << iopt << Exception::runerror; } } ShoKinPtr SudakovFormFactor::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 SudakovFormFactor::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 SudakovFormFactor::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/QTilde/SplittingFunctions/CMWHalfHalfOneSplitFn.h b/Shower/QTilde/SplittingFunctions/CMWHalfHalfOneSplitFn.h --- a/Shower/QTilde/SplittingFunctions/CMWHalfHalfOneSplitFn.h +++ b/Shower/QTilde/SplittingFunctions/CMWHalfHalfOneSplitFn.h @@ -1,163 +1,158 @@ // -*- C++ -*- // // CMWHalfHalfOneSplitFn.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_CMWHalfHalfOneSplitFn_H #define HERWIG_CMWHalfHalfOneSplitFn_H // // This is the declaration of the CMWHalfHalfOneSplitFn class. // #include "HalfHalfOneSplitFn.h" #include "Herwig/Shower/ShowerAlpha.h" namespace Herwig { using namespace ThePEG; /** \ingroup Shower * * This class provides the concrete implementation * of the CMW enhanced expressions for the * splitting function for \f$\frac12\to q\frac12 1\f$. * * The kernel uses the same overestimate as the * corresponding HalfHalfOneSplitFn and thus only needs to * implement the spitting function and ratio to the overestimate. * * TODO: For a more efficient sampling one needs can rewrite the * overestimation to contain the alpha_max*Kgmax factors. * * @see \ref CMWHalfHalfOneSplitFnInterfaces "The interfaces" * defined for CMWHalfHalfOneSplitFn. */ class CMWHalfHalfOneSplitFn: public HalfHalfOneSplitFn { public: /** - * The default constructor. - */ - CMWHalfHalfOneSplitFn() : HalfHalfOneSplitFn() {} - - /** * Methods to return the splitting function. */ //@{ /** * Very similar to HalfHalfOneSplitFn. * Here the kernel only contains the soft part multiplied by the * alphas/2pi * Kg from * Nucl.Phys. B349 (1991) 635-654 * */ virtual double P(const double z, const Energy2 t, const IdList & ids, const bool mass, const RhoDMatrix & rho) const; /** * Very similar to HalfHalfOneSplitFn. * Since we use only the 1/1-z part for overestimating the kernel * in the first place we can keep the same overestimation related functions * for the CMW kernels. */ virtual double ratioP(const double z, const Energy2 t, const IdList & ids, const bool mass, const RhoDMatrix & rho) const; /** * Return the correction term from: * Nucl.Phys. B349 (1991) 635-654 */ double Kg(Energy2 )const{ //TODO: Should be t-dependent int Nf=5;//alpha_->Nf(t) return (3.*(67./18.-1./6.*sqr(Constants::pi))-5./9.*Nf); } //@} 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(){ HalfHalfOneSplitFn::doinit(); }; //@} private: // Pointer to the alpha_s object in use. ShowerAlphaPtr alpha_; // Provide information if the kernel is used for initial state. // as the pt definition contains an additional factor of z. bool isIS_=false; 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. */ CMWHalfHalfOneSplitFn & operator=(const CMWHalfHalfOneSplitFn &); }; } #endif /* HERWIG_CMWHalfHalfOneSplitFn_H */ diff --git a/Shower/QTilde/SplittingFunctions/CMWOneOneOneSplitFn.h b/Shower/QTilde/SplittingFunctions/CMWOneOneOneSplitFn.h --- a/Shower/QTilde/SplittingFunctions/CMWOneOneOneSplitFn.h +++ b/Shower/QTilde/SplittingFunctions/CMWOneOneOneSplitFn.h @@ -1,159 +1,154 @@ // -*- C++ -*- // // CMWOneOneOneSplitFn.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_CMWOneOneOneSplitFn_H #define HERWIG_CMWOneOneOneSplitFn_H // // This is the declaration of the CMWOneOneOneSplitFn class. // #include "OneOneOneSplitFn.h" #include "Herwig/Shower/ShowerAlpha.h" namespace Herwig { using namespace ThePEG; /** \ingroup Shower * * This class provides the concrete implementation * of the CMW enhanced expressions for the * splitting function for \f$1\to 11\f$. * * The kernel uses the same overestimate as the * corresponding OneOneOneSplitFn and thus only needs to * implement the spitting function and ratio to the overestimate. * * TODO: For a more efficient sampling one needs can rewrite the * overestimation to contain the alpha_max*Kgmax factors. * * @see \ref CMWOneOneOneSplitFnInterfaces "The interfaces" * defined for CMWOneOneOneSplitFn. */ class CMWOneOneOneSplitFn: public OneOneOneSplitFn { public: /** - * The default constructor. - */ - CMWOneOneOneSplitFn() : OneOneOneSplitFn() {} - - /** * Methods to return the splitting function. */ //@{ /** * Very similar to HalfHalfOneSplitFn. * Here the kernel only contains the soft part multiplied by the * alphas/2pi * Kg from * Nucl.Phys. B349 (1991) 635-654 * */ virtual double P(const double z, const Energy2 t, const IdList & ids, const bool mass, const RhoDMatrix & rho) const; /** * Very similar to HalfHalfOneSplitFn. * Since we use only the 1/1-z part for overestimating the kernel * in the first place we can keep the same overestimation related functions * for the CMW kernels. */ virtual double ratioP(const double z, const Energy2 t, const IdList & ids, const bool mass, const RhoDMatrix & rho) const; /** * Return the correction term from: * Nucl.Phys. B349 (1991) 635-654 */ double Kg(Energy2 )const{ //TODO: Might be t-dependent int Nf=5;//alpha_->Nf(t) return (3.*(67./18.-1./6.*sqr(Constants::pi))-5./9.*Nf); } //@} 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(){ OneOneOneSplitFn::doinit(); }; //@} private: // Pointer to the alpha_s object in use. ShowerAlphaPtr alpha_; // Provide information if the kernel is used for initial state. // as the pt definition contains an additional factor of z. bool isIS_=false; 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. */ CMWOneOneOneSplitFn & operator=(const CMWOneOneOneSplitFn &); }; } #endif /* HERWIG_CMWOneOneOneSplitFn_H */ diff --git a/Shower/QTilde/SplittingFunctions/HalfHalfOneSplitFn.h b/Shower/QTilde/SplittingFunctions/HalfHalfOneSplitFn.h --- a/Shower/QTilde/SplittingFunctions/HalfHalfOneSplitFn.h +++ b/Shower/QTilde/SplittingFunctions/HalfHalfOneSplitFn.h @@ -1,188 +1,183 @@ // -*- C++ -*- // // HalfHalfOneSplitFn.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_HalfHalfOneSplitFn_H #define HERWIG_HalfHalfOneSplitFn_H // // This is the declaration of the HalfHalfOneSplitFn class. // #include "Herwig/Shower/QTilde/SplittingFunctions/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 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 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/QTilde/SplittingFunctions/HalfOneHalfSplitFn.h b/Shower/QTilde/SplittingFunctions/HalfOneHalfSplitFn.h --- a/Shower/QTilde/SplittingFunctions/HalfOneHalfSplitFn.h +++ b/Shower/QTilde/SplittingFunctions/HalfOneHalfSplitFn.h @@ -1,188 +1,183 @@ // -*- C++ -*- // // HalfOneHalfSplitFn.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_HalfOneHalfSplitFn_H #define HERWIG_HalfOneHalfSplitFn_H // // This is the declaration of the HalfOneHalfSplitFn class. // #include "Herwig/Shower/QTilde/SplittingFunctions/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 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 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/QTilde/SplittingFunctions/OneHalfHalfSplitFn.h b/Shower/QTilde/SplittingFunctions/OneHalfHalfSplitFn.h --- a/Shower/QTilde/SplittingFunctions/OneHalfHalfSplitFn.h +++ b/Shower/QTilde/SplittingFunctions/OneHalfHalfSplitFn.h @@ -1,193 +1,188 @@ // -*- C++ -*- // // OneHalfHalfSplitFn.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_OneHalfHalfSplitFn_H #define HERWIG_OneHalfHalfSplitFn_H // // This is the declaration of the OneHalfHalfSplitFn class. // #include "Herwig/Shower/QTilde/SplittingFunctions/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 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 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/QTilde/SplittingFunctions/OneOneOneMassiveSplitFn.h b/Shower/QTilde/SplittingFunctions/OneOneOneMassiveSplitFn.h --- a/Shower/QTilde/SplittingFunctions/OneOneOneMassiveSplitFn.h +++ b/Shower/QTilde/SplittingFunctions/OneOneOneMassiveSplitFn.h @@ -1,189 +1,184 @@ // -*- C++ -*- // // OneOneOneSplitFn.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_OneOneOneMassiveSplitFn_H #define HERWIG_OneOneOneMassiveSplitFn_H // // This is the declaration of the OneOneOneMassiveSplitFn class. // #include "Herwig/Shower/QTilde/SplittingFunctions/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$ where the emitting particle is massi e * * In this case the splitting function is given by * \f[P(z) =2C\left( \frac{z}{1-z}-\frac{m^2}{q^2-m^2}) +2\rho_{00]\frac{(1-z)^2m^2}{q^2-m^2} + (1-\rho_{00})\left(\frac{1-z}{z}+z(1-z)-\frac{(1-z)^2m^2}{q^2-m^2}\right)\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 OneOneOneMassiveSplitFnInterfaces "The interfaces" * defined for OneOneOneMassiveSplitFn. */ class OneOneOneMassiveSplitFn: public SplittingFunction { public: /** - * The default constructor. - */ - OneOneOneMassiveSplitFn() : 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 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 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. */ OneOneOneMassiveSplitFn & operator=(const OneOneOneMassiveSplitFn &); }; } #endif /* HERWIG_OneOneOneMassiveSplitFn_H */ diff --git a/Shower/QTilde/SplittingFunctions/OneOneOneSplitFn.h b/Shower/QTilde/SplittingFunctions/OneOneOneSplitFn.h --- a/Shower/QTilde/SplittingFunctions/OneOneOneSplitFn.h +++ b/Shower/QTilde/SplittingFunctions/OneOneOneSplitFn.h @@ -1,189 +1,184 @@ // -*- C++ -*- // // OneOneOneSplitFn.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_OneOneOneSplitFn_H #define HERWIG_OneOneOneSplitFn_H // // This is the declaration of the OneOneOneSplitFn class. // #include "Herwig/Shower/QTilde/SplittingFunctions/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 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 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/QTilde/SplittingFunctions/SplittingFunction.cc b/Shower/QTilde/SplittingFunctions/SplittingFunction.cc --- a/Shower/QTilde/SplittingFunctions/SplittingFunction.cc +++ b/Shower/QTilde/SplittingFunctions/SplittingFunction.cc @@ -1,891 +1,891 @@ // -*- C++ -*- // // SplittingFunction.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the SplittingFunction class. // #include "SplittingFunction.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Utilities/EnumIO.h" #include "Herwig/Shower/QTilde/Base/ShowerParticle.h" #include "ThePEG/Utilities/DescribeClass.h" using namespace Herwig; DescribeAbstractClass describeSplittingFunction ("Herwig::SplittingFunction",""); void SplittingFunction::Init() { static ClassDocumentation documentation ("The SplittingFunction class is the based class for 1->2 splitting functions" " in Herwig"); static Switch interfaceColourStructure ("ColourStructure", "The colour structure for the splitting function", &SplittingFunction::_colourStructure, Undefined, false, false); static SwitchOption interfaceColourStructureTripletTripletOctet (interfaceColourStructure, "TripletTripletOctet", "3 -> 3 8", TripletTripletOctet); static SwitchOption interfaceColourStructureOctetOctetOctet (interfaceColourStructure, "OctetOctetOctet", "8 -> 8 8", OctetOctetOctet); static SwitchOption interfaceColourStructureOctetTripletTriplet (interfaceColourStructure, "OctetTripletTriplet", "8 -> 3 3bar", OctetTripletTriplet); static SwitchOption interfaceColourStructureTripletOctetTriplet (interfaceColourStructure, "TripletOctetTriplet", "3 -> 8 3", TripletOctetTriplet); static SwitchOption interfaceColourStructureSextetSextetOctet (interfaceColourStructure, "SextetSextetOctet", "6 -> 6 8", SextetSextetOctet); static SwitchOption interfaceColourStructureChargedChargedNeutral (interfaceColourStructure, "ChargedChargedNeutral", "q -> q 0", ChargedChargedNeutral); static SwitchOption interfaceColourStructureNeutralChargedCharged (interfaceColourStructure, "NeutralChargedCharged", "0 -> q qbar", NeutralChargedCharged); static SwitchOption interfaceColourStructureChargedNeutralCharged (interfaceColourStructure, "ChargedNeutralCharged", "q -> 0 q", ChargedNeutralCharged); static Switch interfaceInteractionType ("InteractionType", "Type of the interaction", &SplittingFunction::_interactionType, ShowerInteraction::UNDEFINED, false, false); static SwitchOption interfaceInteractionTypeQCD (interfaceInteractionType, "QCD","QCD",ShowerInteraction::QCD); static SwitchOption interfaceInteractionTypeQED (interfaceInteractionType, "QED","QED",ShowerInteraction::QED); static Switch interfaceAngularOrdered ("AngularOrdered", "Whether or not this interaction is angular ordered, " "normally only g->q qbar and gamma-> f fbar are the only ones which aren't.", &SplittingFunction::angularOrdered_, true, false, false); static SwitchOption interfaceAngularOrderedYes (interfaceAngularOrdered, "Yes", "Interaction is angular ordered", true); static SwitchOption interfaceAngularOrderedNo (interfaceAngularOrdered, "No", "Interaction isn't angular ordered", false); static Switch interfaceScaleChoice ("ScaleChoice", "The scale choice to be used", &SplittingFunction::scaleChoice_, 2, false, false); static SwitchOption interfaceScaleChoicepT (interfaceScaleChoice, "pT", "pT of the branching", 0); static SwitchOption interfaceScaleChoiceQ2 (interfaceScaleChoice, "Q2", "Q2 of the branching", 1); static SwitchOption interfaceScaleChoiceFromAngularOrdering (interfaceScaleChoice, "FromAngularOrdering", "If angular order use pT, otherwise Q2", 2); } void SplittingFunction::persistentOutput(PersistentOStream & os) const { - os << oenum(_interactionType) << _interactionOrder + os << oenum(_interactionType) << oenum(_colourStructure) << _colourFactor << angularOrdered_ << scaleChoice_; } void SplittingFunction::persistentInput(PersistentIStream & is, int) { - is >> ienum(_interactionType) >> _interactionOrder + is >> ienum(_interactionType) >> ienum(_colourStructure) >> _colourFactor >> angularOrdered_ >> scaleChoice_; } void SplittingFunction::colourConnection(tShowerParticlePtr parent, tShowerParticlePtr first, tShowerParticlePtr second, ShowerPartnerType partnerType, const bool back) const { if(_colourStructure==TripletTripletOctet) { if(!back) { ColinePair cparent = ColinePair(parent->colourLine(), parent->antiColourLine()); // ensure input consistency assert(( cparent.first && !cparent.second && partnerType==ShowerPartnerType::QCDColourLine) || ( !cparent.first && cparent.second && partnerType==ShowerPartnerType::QCDAntiColourLine)); // q -> q g if(cparent.first) { ColinePtr newline=new_ptr(ColourLine()); cparent.first->addColoured(second); newline->addColoured ( first); newline->addAntiColoured (second); } // qbar -> qbar g else { ColinePtr newline=new_ptr(ColourLine()); cparent.second->addAntiColoured(second); newline->addColoured(second); newline->addAntiColoured(first); } // Set progenitor first->progenitor(parent->progenitor()); second->progenitor(parent->progenitor()); } else { ColinePair cfirst = ColinePair(first->colourLine(), first->antiColourLine()); // ensure input consistency assert(( cfirst.first && !cfirst.second && partnerType==ShowerPartnerType::QCDColourLine) || ( !cfirst.first && cfirst.second && partnerType==ShowerPartnerType::QCDAntiColourLine)); // q -> q g if(cfirst.first) { ColinePtr newline=new_ptr(ColourLine()); cfirst.first->addAntiColoured(second); newline->addColoured(second); newline->addColoured(parent); } // qbar -> qbar g else { ColinePtr newline=new_ptr(ColourLine()); cfirst.second->addColoured(second); newline->addAntiColoured(second); newline->addAntiColoured(parent); } // Set progenitor parent->progenitor(first->progenitor()); second->progenitor(first->progenitor()); } } else if(_colourStructure==OctetOctetOctet) { if(!back) { ColinePair cparent = ColinePair(parent->colourLine(), parent->antiColourLine()); // ensure input consistency assert(cparent.first&&cparent.second); // ensure first gluon is hardest if( first->id()==second->id() && parent->showerKinematics()->z()<0.5 ) swap(first,second); // colour line radiates if(partnerType==ShowerPartnerType::QCDColourLine) { // The colour line is radiating ColinePtr newline=new_ptr(ColourLine()); cparent.first->addColoured(second); cparent.second->addAntiColoured(first); newline->addColoured(first); newline->addAntiColoured(second); } // anti colour line radiates else if(partnerType==ShowerPartnerType::QCDAntiColourLine) { ColinePtr newline=new_ptr(ColourLine()); cparent.first->addColoured(first); cparent.second->addAntiColoured(second); newline->addColoured(second); newline->addAntiColoured(first); } else assert(false); } else { ColinePair cfirst = ColinePair(first->colourLine(), first->antiColourLine()); // ensure input consistency assert(cfirst.first&&cfirst.second); // The colour line is radiating if(partnerType==ShowerPartnerType::QCDColourLine) { ColinePtr newline=new_ptr(ColourLine()); cfirst.first->addAntiColoured(second); cfirst.second->addAntiColoured(parent); newline->addColoured(parent); newline->addColoured(second); } // anti colour line radiates else if(partnerType==ShowerPartnerType::QCDAntiColourLine) { ColinePtr newline=new_ptr(ColourLine()); cfirst.first->addColoured(parent); cfirst.second->addColoured(second); newline->addAntiColoured(second); newline->addAntiColoured(parent); } else assert(false); } } else if(_colourStructure == OctetTripletTriplet) { if(!back) { ColinePair cparent = ColinePair(parent->colourLine(), parent->antiColourLine()); // ensure input consistency assert(cparent.first&&cparent.second); cparent.first ->addColoured ( first); cparent.second->addAntiColoured(second); // Set progenitor first->progenitor(parent->progenitor()); second->progenitor(parent->progenitor()); } else { ColinePair cfirst = ColinePair(first->colourLine(), first->antiColourLine()); // ensure input consistency assert(( cfirst.first && !cfirst.second) || (!cfirst.first && cfirst.second)); // g -> q qbar if(cfirst.first) { ColinePtr newline=new_ptr(ColourLine()); cfirst.first->addColoured(parent); newline->addAntiColoured(second); newline->addAntiColoured(parent); } // g -> qbar q else { ColinePtr newline=new_ptr(ColourLine()); cfirst.second->addAntiColoured(parent); newline->addColoured(second); newline->addColoured(parent); } // Set progenitor parent->progenitor(first->progenitor()); second->progenitor(first->progenitor()); } } else if(_colourStructure == TripletOctetTriplet) { if(!back) { ColinePair cparent = ColinePair(parent->colourLine(), parent->antiColourLine()); // ensure input consistency assert(( cparent.first && !cparent.second) || (!cparent.first && cparent.second)); // q -> g q if(cparent.first) { ColinePtr newline=new_ptr(ColourLine()); cparent.first->addColoured(first); newline->addColoured (second); newline->addAntiColoured( first); } // qbar -> g qbar else { ColinePtr newline=new_ptr(ColourLine()); cparent.second->addAntiColoured(first); newline->addColoured ( first); newline->addAntiColoured(second); } // Set progenitor first->progenitor(parent->progenitor()); second->progenitor(parent->progenitor()); } else { ColinePair cfirst = ColinePair(first->colourLine(), first->antiColourLine()); // ensure input consistency assert(cfirst.first&&cfirst.second); // q -> g q if(parent->id()>0) { cfirst.first ->addColoured(parent); cfirst.second->addColoured(second); } else { cfirst.first ->addAntiColoured(second); cfirst.second->addAntiColoured(parent); } // Set progenitor parent->progenitor(first->progenitor()); second->progenitor(first->progenitor()); } } else if(_colourStructure==SextetSextetOctet) { //make sure we're not doing backward evolution assert(!back); //make sure something sensible assert(parent->colourLine() || parent->antiColourLine()); //get the colour lines or anti-colour lines bool isAntiColour=true; ColinePair cparent; if(parent->colourLine()) { cparent = ColinePair(const_ptr_cast(parent->colourInfo()->colourLines()[0]), const_ptr_cast(parent->colourInfo()->colourLines()[1])); isAntiColour=false; } else { cparent = ColinePair(const_ptr_cast(parent->colourInfo()->antiColourLines()[0]), const_ptr_cast(parent->colourInfo()->antiColourLines()[1])); } //check for sensible input // assert(cparent.first && cparent.second); // sextet has 2 colour lines if(!isAntiColour) { //pick at random which of the colour topolgies to take double topology = UseRandom::rnd(); if(topology < 0.25) { ColinePtr newline=new_ptr(ColourLine()); cparent.first->addColoured(second); cparent.second->addColoured(first); newline->addColoured(first); newline->addAntiColoured(second); } else if(topology >=0.25 && topology < 0.5) { ColinePtr newline=new_ptr(ColourLine()); cparent.first->addColoured(first); cparent.second->addColoured(second); newline->addColoured(first); newline->addAntiColoured(second); } else if(topology >= 0.5 && topology < 0.75) { ColinePtr newline=new_ptr(ColourLine()); cparent.first->addColoured(second); cparent.second->addColoured(first); newline->addColoured(first); newline->addAntiColoured(second); } else { ColinePtr newline=new_ptr(ColourLine()); cparent.first->addColoured(first); cparent.second->addColoured(second); newline->addColoured(first); newline->addAntiColoured(second); } } // sextet has 2 anti-colour lines else { double topology = UseRandom::rnd(); if(topology < 0.25){ ColinePtr newline=new_ptr(ColourLine()); cparent.first->addAntiColoured(second); cparent.second->addAntiColoured(first); newline->addAntiColoured(first); newline->addColoured(second); } else if(topology >=0.25 && topology < 0.5) { ColinePtr newline=new_ptr(ColourLine()); cparent.first->addAntiColoured(first); cparent.second->addAntiColoured(second); newline->addAntiColoured(first); newline->addColoured(second); } else if(topology >= 0.5 && topology < 0.75) { ColinePtr newline=new_ptr(ColourLine()); cparent.first->addAntiColoured(second); cparent.second->addAntiColoured(first); newline->addAntiColoured(first); newline->addColoured(second); } else { ColinePtr newline=new_ptr(ColourLine()); cparent.first->addAntiColoured(first); cparent.second->addAntiColoured(second); newline->addAntiColoured(first); newline->addColoured(second); } } } else if(_colourStructure == ChargedChargedNeutral) { if(!parent->data().coloured()) return; if(!back) { ColinePair cparent = ColinePair(parent->colourLine(), parent->antiColourLine()); // q -> q g if(cparent.first) { cparent.first->addColoured(first); } // qbar -> qbar g if(cparent.second) { cparent.second->addAntiColoured(first); } } else { ColinePair cfirst = ColinePair(first->colourLine(), first->antiColourLine()); // q -> q g if(cfirst.first) { cfirst.first->addColoured(parent); } // qbar -> qbar g if(cfirst.second) { cfirst.second->addAntiColoured(parent); } } } else if(_colourStructure == ChargedNeutralCharged) { if(!parent->data().coloured()) return; if(!back) { ColinePair cparent = ColinePair(parent->colourLine(), parent->antiColourLine()); // q -> q g if(cparent.first) { cparent.first->addColoured(second); } // qbar -> qbar g if(cparent.second) { cparent.second->addAntiColoured(second); } } else { if (second->dataPtr()->iColour()==PDT::Colour3 ) { ColinePtr newline=new_ptr(ColourLine()); newline->addColoured(second); newline->addColoured(parent); } else if (second->dataPtr()->iColour()==PDT::Colour3bar ) { ColinePtr newline=new_ptr(ColourLine()); newline->addAntiColoured(second); newline->addAntiColoured(parent); } } } else if(_colourStructure == NeutralChargedCharged ) { if(!back) { if(first->dataPtr()->coloured()) { ColinePtr newline=new_ptr(ColourLine()); if(first->dataPtr()->iColour()==PDT::Colour3) { newline->addColoured (first ); newline->addAntiColoured(second); } else if (first->dataPtr()->iColour()==PDT::Colour3bar) { newline->addColoured (second); newline->addAntiColoured(first ); } else assert(false); } } else { ColinePair cfirst = ColinePair(first->colourLine(), first->antiColourLine()); // gamma -> q qbar if(cfirst.first) { cfirst.first->addAntiColoured(second); } // gamma -> qbar q else if(cfirst.second) { cfirst.second->addColoured(second); } else assert(false); } } else { assert(false); } } void SplittingFunction::doinit() { Interfaced::doinit(); assert(_interactionType!=ShowerInteraction::UNDEFINED); assert((_colourStructure>0&&_interactionType==ShowerInteraction::QCD) || (_colourStructure<0&&_interactionType==ShowerInteraction::QED) ); if(_colourFactor>0.) return; // compute the colour factors if need if(_colourStructure==TripletTripletOctet) { _colourFactor = 4./3.; } else if(_colourStructure==OctetOctetOctet) { _colourFactor = 3.; } else if(_colourStructure==OctetTripletTriplet) { _colourFactor = 0.5; } else if(_colourStructure==TripletOctetTriplet) { _colourFactor = 4./3.; } else if(_colourStructure==SextetSextetOctet) { _colourFactor = 10./3.; } else if(_colourStructure<0) { _colourFactor = 1.; } else { assert(false); } } bool SplittingFunction::checkColours(const IdList & ids) const { if(_colourStructure==TripletTripletOctet) { if(ids[0]!=ids[1]) return false; if((ids[0]->iColour()==PDT::Colour3||ids[0]->iColour()==PDT::Colour3bar) && ids[2]->iColour()==PDT::Colour8) return true; return false; } else if(_colourStructure==OctetOctetOctet) { for(unsigned int ix=0;ix<3;++ix) { if(ids[ix]->iColour()!=PDT::Colour8) return false; } return true; } else if(_colourStructure==OctetTripletTriplet) { if(ids[0]->iColour()!=PDT::Colour8) return false; if(ids[1]->iColour()==PDT::Colour3&&ids[2]->iColour()==PDT::Colour3bar) return true; if(ids[1]->iColour()==PDT::Colour3bar&&ids[2]->iColour()==PDT::Colour3) return true; return false; } else if(_colourStructure==TripletOctetTriplet) { if(ids[0]!=ids[2]) return false; if((ids[0]->iColour()==PDT::Colour3||ids[0]->iColour()==PDT::Colour3bar) && ids[1]->iColour()==PDT::Colour8) return true; return false; } else if(_colourStructure==SextetSextetOctet) { if(ids[0]!=ids[1]) return false; if((ids[0]->iColour()==PDT::Colour6 || ids[0]->iColour()==PDT::Colour6bar) && ids[2]->iColour()==PDT::Colour8) return true; return false; } else if(_colourStructure==ChargedChargedNeutral) { if(ids[0]!=ids[1]) return false; if(ids[2]->iCharge()!=0) return false; if(ids[0]->iCharge()==ids[1]->iCharge()) return true; return false; } else if(_colourStructure==ChargedNeutralCharged) { if(ids[0]!=ids[2]) return false; if(ids[1]->iCharge()!=0) return false; if(ids[0]->iCharge()==ids[2]->iCharge()) return true; return false; } else if(_colourStructure==NeutralChargedCharged) { if(ids[1]->id()!=-ids[2]->id()) return false; if(ids[0]->iCharge()!=0) return false; if(ids[1]->iCharge()==-ids[2]->iCharge()) return true; return false; } else { assert(false); } return false; } namespace { bool hasColour(tPPtr p) { PDT::Colour colour = p->dataPtr()->iColour(); return colour==PDT::Colour3 || colour==PDT::Colour8 || colour == PDT::Colour6; } bool hasAntiColour(tPPtr p) { PDT::Colour colour = p->dataPtr()->iColour(); return colour==PDT::Colour3bar || colour==PDT::Colour8 || colour == PDT::Colour6bar; } } void SplittingFunction::evaluateFinalStateScales(ShowerPartnerType partnerType, Energy scale, double z, tShowerParticlePtr parent, tShowerParticlePtr emitter, tShowerParticlePtr emitted) { // identify emitter and emitted double zEmitter = z, zEmitted = 1.-z; bool bosonSplitting(false); // special for g -> gg, particle highest z is emitter if(emitter->id() == emitted->id() && emitter->id() == parent->id() && zEmitted > zEmitter) { swap(zEmitted,zEmitter); swap( emitted, emitter); } // otherwise if particle ID same else if(emitted->id()==parent->id()) { swap(zEmitted,zEmitter); swap( emitted, emitter); } // no real emitter/emitted else if(emitter->id()!=parent->id()) { bosonSplitting = true; } // may need to add angularOrder flag here // now the various scales // QED if(partnerType==ShowerPartnerType::QED) { assert(colourStructure()==ChargedChargedNeutral || colourStructure()==ChargedNeutralCharged || colourStructure()==NeutralChargedCharged ); // normal case if(!bosonSplitting) { assert(colourStructure()==ChargedChargedNeutral || colourStructure()==ChargedNeutralCharged ); // set the scales // emitter emitter->scales().QED = zEmitter*scale; emitter->scales().QED_noAO = scale; emitter->scales().QCD_c = min(scale,parent->scales().QCD_c ); emitter->scales().QCD_c_noAO = min(scale,parent->scales().QCD_c_noAO ); emitter->scales().QCD_ac = min(scale,parent->scales().QCD_ac ); emitter->scales().QCD_ac_noAO = min(scale,parent->scales().QCD_ac_noAO); // emitted emitted->scales().QED = zEmitted*scale; emitted->scales().QED_noAO = scale; emitted->scales().QCD_c = ZERO; emitted->scales().QCD_c_noAO = ZERO; emitted->scales().QCD_ac = ZERO; emitted->scales().QCD_ac_noAO = ZERO; } // gamma -> f fbar else { assert(colourStructure()==NeutralChargedCharged ); // emitter emitter->scales().QED = zEmitter*scale; emitter->scales().QED_noAO = scale; if(hasColour(emitter)) { emitter->scales().QCD_c = zEmitter*scale; emitter->scales().QCD_c_noAO = scale; } if(hasAntiColour(emitter)) { emitter->scales().QCD_ac = zEmitter*scale; emitter->scales().QCD_ac_noAO = scale; } // emitted emitted->scales().QED = zEmitted*scale; emitted->scales().QED_noAO = scale; if(hasColour(emitted)) { emitted->scales().QCD_c = zEmitted*scale; emitted->scales().QCD_c_noAO = scale; } if(hasAntiColour(emitted)) { emitted->scales().QCD_ac = zEmitted*scale; emitted->scales().QCD_ac_noAO = scale; } } } // QCD else { // normal case eg q -> q g and g -> g g if(!bosonSplitting) { emitter->scales().QED = min(scale,parent->scales().QED ); emitter->scales().QED_noAO = min(scale,parent->scales().QED_noAO); if(partnerType==ShowerPartnerType::QCDColourLine) { emitter->scales().QCD_c = zEmitter*scale; emitter->scales().QCD_c_noAO = scale; emitter->scales().QCD_ac = min(zEmitter*scale,parent->scales().QCD_ac ); emitter->scales().QCD_ac_noAO = min( scale,parent->scales().QCD_ac_noAO); } else { emitter->scales().QCD_c = min(zEmitter*scale,parent->scales().QCD_c ); emitter->scales().QCD_c_noAO = min( scale,parent->scales().QCD_c_noAO ); emitter->scales().QCD_ac = zEmitter*scale; emitter->scales().QCD_ac_noAO = scale; } // emitted emitted->scales().QED = ZERO; emitted->scales().QED_noAO = ZERO; emitted->scales().QCD_c = zEmitted*scale; emitted->scales().QCD_c_noAO = scale; emitted->scales().QCD_ac = zEmitted*scale; emitted->scales().QCD_ac_noAO = scale; } // g -> q qbar else { // emitter if(emitter->dataPtr()->charged()) { emitter->scales().QED = zEmitter*scale; emitter->scales().QED_noAO = scale; } emitter->scales().QCD_c = zEmitter*scale; emitter->scales().QCD_c_noAO = scale; emitter->scales().QCD_ac = zEmitter*scale; emitter->scales().QCD_ac_noAO = scale; // emitted if(emitted->dataPtr()->charged()) { emitted->scales().QED = zEmitted*scale; emitted->scales().QED_noAO = scale; } emitted->scales().QCD_c = zEmitted*scale; emitted->scales().QCD_c_noAO = scale; emitted->scales().QCD_ac = zEmitted*scale; emitted->scales().QCD_ac_noAO = scale; } } } void SplittingFunction::evaluateInitialStateScales(ShowerPartnerType partnerType, Energy scale, double z, tShowerParticlePtr parent, tShowerParticlePtr spacelike, tShowerParticlePtr timelike) { // scale for time-like child Energy AOScale = (1.-z)*scale; // QED if(partnerType==ShowerPartnerType::QED) { if(parent->id()==spacelike->id()) { // parent parent ->scales().QED = scale; parent ->scales().QED_noAO = scale; parent ->scales().QCD_c = min(scale,spacelike->scales().QCD_c ); parent ->scales().QCD_c_noAO = min(scale,spacelike->scales().QCD_c_noAO ); parent ->scales().QCD_ac = min(scale,spacelike->scales().QCD_ac ); parent ->scales().QCD_ac_noAO = min(scale,spacelike->scales().QCD_ac_noAO); // timelike timelike->scales().QED = AOScale; timelike->scales().QED_noAO = scale; timelike->scales().QCD_c = ZERO; timelike->scales().QCD_c_noAO = ZERO; timelike->scales().QCD_ac = ZERO; timelike->scales().QCD_ac_noAO = ZERO; } else if(parent->id()==timelike->id()) { parent ->scales().QED = scale; parent ->scales().QED_noAO = scale; if(hasColour(parent)) { parent ->scales().QCD_c = scale; parent ->scales().QCD_c_noAO = scale; } if(hasAntiColour(parent)) { parent ->scales().QCD_ac = scale; parent ->scales().QCD_ac_noAO = scale; } // timelike timelike->scales().QED = AOScale; timelike->scales().QED_noAO = scale; if(hasColour(timelike)) { timelike->scales().QCD_c = AOScale; timelike->scales().QCD_c_noAO = scale; } if(hasAntiColour(timelike)) { timelike->scales().QCD_ac = AOScale; timelike->scales().QCD_ac_noAO = scale; } } else { parent ->scales().QED = scale; parent ->scales().QED_noAO = scale; parent ->scales().QCD_c = ZERO ; parent ->scales().QCD_c_noAO = ZERO ; parent ->scales().QCD_ac = ZERO ; parent ->scales().QCD_ac_noAO = ZERO ; // timelike timelike->scales().QED = AOScale; timelike->scales().QED_noAO = scale; if(hasColour(timelike)) { timelike->scales().QCD_c = min(AOScale,spacelike->scales().QCD_ac ); timelike->scales().QCD_c_noAO = min( scale,spacelike->scales().QCD_ac_noAO); } if(hasAntiColour(timelike)) { timelike->scales().QCD_ac = min(AOScale,spacelike->scales().QCD_c ); timelike->scales().QCD_ac_noAO = min( scale,spacelike->scales().QCD_c_noAO ); } } } // QCD else { // timelike if(timelike->dataPtr()->charged()) { timelike->scales().QED = AOScale; timelike->scales().QED_noAO = scale; } if(hasColour(timelike)) { timelike->scales().QCD_c = AOScale; timelike->scales().QCD_c_noAO = scale; } if(hasAntiColour(timelike)) { timelike->scales().QCD_ac = AOScale; timelike->scales().QCD_ac_noAO = scale; } if(parent->id()==spacelike->id()) { parent ->scales().QED = min(scale,spacelike->scales().QED ); parent ->scales().QED_noAO = min(scale,spacelike->scales().QED_noAO ); parent ->scales().QCD_c = min(scale,spacelike->scales().QCD_c ); parent ->scales().QCD_c_noAO = min(scale,spacelike->scales().QCD_c_noAO ); parent ->scales().QCD_ac = min(scale,spacelike->scales().QCD_ac ); parent ->scales().QCD_ac_noAO = min(scale,spacelike->scales().QCD_ac_noAO); } else { if(parent->dataPtr()->charged()) { parent ->scales().QED = scale; parent ->scales().QED_noAO = scale; } if(hasColour(parent)) { parent ->scales().QCD_c = scale; parent ->scales().QCD_c_noAO = scale; } if(hasAntiColour(parent)) { parent ->scales().QCD_ac = scale; parent ->scales().QCD_ac_noAO = scale; } } } } void SplittingFunction::evaluateDecayScales(ShowerPartnerType partnerType, Energy scale, double z, tShowerParticlePtr parent, tShowerParticlePtr spacelike, tShowerParticlePtr timelike) { assert(parent->id()==spacelike->id()); // angular-ordered scale for 2nd child Energy AOScale = (1.-z)*scale; // QED if(partnerType==ShowerPartnerType::QED) { // timelike timelike->scales().QED = AOScale; timelike->scales().QED_noAO = scale; timelike->scales().QCD_c = ZERO; timelike->scales().QCD_c_noAO = ZERO; timelike->scales().QCD_ac = ZERO; timelike->scales().QCD_ac_noAO = ZERO; // spacelike spacelike->scales().QED = scale; spacelike->scales().QED_noAO = scale; } // QCD else { // timelike timelike->scales().QED = ZERO; timelike->scales().QED_noAO = ZERO; timelike->scales().QCD_c = AOScale; timelike->scales().QCD_c_noAO = scale; timelike->scales().QCD_ac = AOScale; timelike->scales().QCD_ac_noAO = scale; // spacelike spacelike->scales().QED = max(scale,parent->scales().QED ); spacelike->scales().QED_noAO = max(scale,parent->scales().QED_noAO ); } spacelike->scales().QCD_c = max(scale,parent->scales().QCD_c ); spacelike->scales().QCD_c_noAO = max(scale,parent->scales().QCD_c_noAO ); spacelike->scales().QCD_ac = max(scale,parent->scales().QCD_ac ); spacelike->scales().QCD_ac_noAO = max(scale,parent->scales().QCD_ac_noAO); } diff --git a/Shower/QTilde/SplittingFunctions/SplittingFunction.h b/Shower/QTilde/SplittingFunctions/SplittingFunction.h --- a/Shower/QTilde/SplittingFunctions/SplittingFunction.h +++ b/Shower/QTilde/SplittingFunctions/SplittingFunction.h @@ -1,388 +1,377 @@ // -*- C++ -*- // // SplittingFunction.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_SplittingFunction_H #define HERWIG_SplittingFunction_H // // This is the declaration of the SplittingFunction class. // #include "ThePEG/Interface/Interfaced.h" #include "Herwig/Shower/QTilde/ShowerConfig.h" #include "ThePEG/EventRecord/RhoDMatrix.h" #include "Herwig/Decay/DecayMatrixElement.h" #include "Herwig/Shower/QTilde/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}; /** \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) + SplittingFunction() : Interfaced(), _interactionType(ShowerInteraction::UNDEFINED), - _interactionOrder(b), _colourStructure(Undefined), _colourFactor(-1.), angularOrdered_(true), scaleChoice_(2) {} public: /** * Methods to return the interaction type and order for the splitting function */ //@{ /** * Return the type of the interaction */ ShowerInteraction interactionType() const {return _interactionType;} /** - * Return the 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) { return sqr(double(ids[0]->iCharge())/3.); } else if(_colourStructure==NeutralChargedCharged) { double fact = sqr(double(ids[1]->iCharge())/3.); if(ids[1]->coloured()) fact *= abs(double(ids[1]->iColour())); return fact; } else assert(false); } else assert(false); } //@} /** * Purely virtual method which should determine whether this splitting * function can be used for a given set of particles. * @param ids The PDG codes for the particles in the splitting. */ virtual bool accept(const IdList & ids) const = 0; /** * Method to check the colours are correct */ virtual bool checkColours(const IdList & ids) const; /** * Methods to return the splitting function. */ //@{ /** * Purely virtual method which should return the exact value of the splitting function, * \f$P\f$ evaluated in terms of the energy fraction, \f$z\f$, and the evolution scale \f$\tilde{q}^2\f$. * @param z The energy fraction. * @param t The scale \f$t=2p_j\cdot p_k\f$. * @param ids The PDG codes for the particles in the splitting. * @param mass Whether or not to include the mass dependent terms */ virtual double P(const double z, const Energy2 t, const IdList & ids, const bool mass, const RhoDMatrix & rho) const = 0; /** * Purely virtual method which should return * an overestimate of the splitting function, * \f$P_{\rm over}\f$ such that the result \f$P_{\rm over}\geq P\f$. This function * should be simple enough that it does not depend on the evolution scale. * @param z The energy fraction. * @param ids The PDG codes for the particles in the splitting. */ virtual double overestimateP(const double z, const IdList & ids) const = 0; /** * Purely virtual method which should return * the ratio of the splitting function to the overestimate, i.e. * \f$P(z,\tilde{q}^2)/P_{\rm over}(z)\f$. * @param z The energy fraction. * @param t The scale \f$t=2p_j\cdot p_k\f$. * @param ids The PDG codes for the particles in the splitting. * @param mass Whether or not to include the mass dependent terms */ virtual double ratioP(const double z, const Energy2 t, const IdList & ids, const bool mass, const RhoDMatrix & rho) const = 0; /** * Purely virtual method which should return the indefinite integral of the * overestimated splitting function, \f$P_{\rm over}\f$. * @param z The energy fraction. * @param ids The PDG codes for the particles in the splitting. * @param PDFfactor Which additional factor to include for the PDF * 0 is no additional factor, * 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$ * */ virtual double integOverP(const double z, const IdList & ids, unsigned int PDFfactor=0) const = 0; /** * Purely virtual method which should return the inverse of the * indefinite integral of the * overestimated splitting function, \f$P_{\rm over}\f$ which is used to * generate the value of \f$z\f$. * @param r Value of the splitting function to be inverted * @param ids The PDG codes for the particles in the splitting. * @param PDFfactor Which additional factor to include for the PDF * 0 is no additional factor, * 1 is \f$1/z\f$, 2 is \f$1/(1-z)\f$ and 3 is \f$1/z/(1-z)\f$ */ virtual double invIntegOverP(const double r, const IdList & ids, unsigned int PDFfactor=0) const = 0; //@} /** * Purely virtual method which should make the proper colour connection * between the emitting parent and the branching products. * @param parent The parent for the branching * @param first The first branching product * @param second The second branching product * @param partnerType The type of evolution partner * @param back Whether this is foward or backward evolution. */ virtual void colourConnection(tShowerParticlePtr parent, tShowerParticlePtr first, tShowerParticlePtr second, ShowerPartnerType partnerType, const bool back) const; /** * Method to calculate the azimuthal angle for forward evolution * @param z The energy fraction * @param t The scale \f$t=2p_j\cdot p_k\f$. * @param ids The PDG codes for the particles in the splitting. * @param The azimuthal angle, \f$\phi\f$. * @return The weight */ virtual vector > 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_;} /** * Scale choice */ bool pTScale() const { return scaleChoice_ == 2 ? angularOrdered_ : scaleChoice_ == 0; } /** * Functions to state scales after branching happens */ //@{ /** * Sort out scales for final-state emission */ void evaluateFinalStateScales(ShowerPartnerType type, Energy scale, double z, tShowerParticlePtr parent, tShowerParticlePtr first, tShowerParticlePtr second); /** * Sort out scales for initial-state emission */ void evaluateInitialStateScales(ShowerPartnerType type, Energy scale, double z, tShowerParticlePtr parent, tShowerParticlePtr first, tShowerParticlePtr second); /** * Sort out scales for decay emission */ void evaluateDecayScales(ShowerPartnerType type, Energy scale, double z, tShowerParticlePtr parent, tShowerParticlePtr first, tShowerParticlePtr second); //@} public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); //@} protected: /** * Set the colour factor */ void colourFactor(double in) {_colourFactor=in;} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ SplittingFunction & operator=(const SplittingFunction &); private: /** * The interaction type for the splitting function. */ ShowerInteraction _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_; /** * The choice of scale */ unsigned int scaleChoice_; }; } #endif /* HERWIG_SplittingFunction_H */ diff --git a/Shower/QTilde/SplittingFunctions/ZeroZeroOneSplitFn.h b/Shower/QTilde/SplittingFunctions/ZeroZeroOneSplitFn.h --- a/Shower/QTilde/SplittingFunctions/ZeroZeroOneSplitFn.h +++ b/Shower/QTilde/SplittingFunctions/ZeroZeroOneSplitFn.h @@ -1,191 +1,186 @@ // -*- C++ -*- // // ZeroZeroOneSplitFn.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_ZeroZeroOneSplitFn_H #define HERWIG_ZeroZeroOneSplitFn_H // // This is the declaration of the ZeroZeroOneSplitFn class. // #include "Herwig/Shower/QTilde/SplittingFunctions/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 */ virtual double P(const double z, const Energy2 t, const IdList & ids, 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 */ virtual double ratioP(const double z, const Energy2 t, const IdList & ids, 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 */