diff --git a/Shower/QTilde/SplittingFunctions/SudakovFormFactor.cc b/Shower/QTilde/SplittingFunctions/SudakovFormFactor.cc --- a/Shower/QTilde/SplittingFunctions/SudakovFormFactor.cc +++ b/Shower/QTilde/SplittingFunctions/SudakovFormFactor.cc @@ -1,1231 +1,1232 @@ // -*- C++ -*- // // SudakovFormFactor.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 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 "Herwig/Shower/QTilde/Kinematics/ShowerKinematics.h" #include "Herwig/Shower/QTilde/Base/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 "Herwig/Shower/QTilde/Kinematics/KinematicHelpers.h" #include "SudakovCutOff.h" #include using std::array; using namespace Herwig; DescribeClass describeSudakovFormFactor ("Herwig::SudakovFormFactor",""); void SudakovFormFactor::persistentOutput(PersistentOStream & os) const { os << splittingFn_ << alpha_ << pdfmax_ << particles_ << pdffactor_ << cutoff_; } void SudakovFormFactor::persistentInput(PersistentIStream & is, int) { is >> splittingFn_ >> alpha_ >> pdfmax_ >> particles_ >> pdffactor_ >> cutoff_; } 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 Reference interfaceCutoff("Cutoff", "A reference to the SudakovCutOff object", &Herwig::SudakovFormFactor::cutoff_, 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); } 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(); + if(ShowerHandler::currentHandler()) + factor *= ShowerHandler::currentHandler()->renormalizationScaleFactor(); return alpha_->showerRatio(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_); const double newpdf=pdf_->xfx(beam,parton0,theScale,x/z()); if(newpdf<=0.) return 0.; const double oldpdf=pdf_->xfx(beam,parton1,theScale,x); if(oldpdf<=0.) return 1.; const 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; } } } } void SudakovFormFactor::guesstz(Energy2 t1,unsigned int iopt, const IdList &ids, double enhance,bool ident, double detune, Energy2 &t_main, double &z_main){ unsigned int pdfopt = iopt!=1 ? 0 : pdffactor_; double lower = splittingFn_->integOverP(zlimits_.first ,ids,pdfopt); double upper = splittingFn_->integOverP(zlimits_.second,ids,pdfopt); double c = 1./((upper - lower) * alpha_->showerOverestimateValue()/Constants::twopi*enhance*detune); double r = UseRandom::rnd(); assert(iopt<=2); if(iopt==1) { c/=pdfmax_; //symmetry of FS gluon splitting if(ident) c*=0.5; } else if(iopt==2) c*=-1.; // guessing t if(iopt!=2 || c*log(r)invIntegOverP(lower + UseRandom::rnd() *(upper - lower),ids,pdfopt); } 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 guesstz(told,0,ids_,enhance,ids_[1]==ids_[2],detune,t,z_); // 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 guesstz(told,1,ids_,enhance,ids_[1]==ids_[2],detune,t,z_); // actual values for z-limits if(!computeSpaceLikeLimits(t,x)) return false; if(t zlimits_.second) return true; Energy2 m02 = (ids_[0]->id()!=ParticleID::g && ids_[0]->id()!=ParticleID::gamma) ? masssquared_[0] : Energy2(); Energy2 pt2 = QTildeKinematics::pT2_FSR(t,z(),m02,masssquared_[1],masssquared_[2], masssquared_[1],masssquared_[2]); // if pt2<0 veto if (pt2pT2min()) return true; // otherwise calculate pt and return pT_ = sqrt(pt2); return false; } ShoKinPtr SudakovFormFactor::generateNextTimeBranching(const Energy startingScale, 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 the method. q_ = ZERO; z_ = 0.; phi_ = 0.; // perform initialization Energy2 tmax(sqr(startingScale)),tmin; initialize(ids,tmin); // check max > min if(tmax<=tmin) return ShoKinPtr(); // calculate next value of t using veto algorithm Energy2 t(tmax); // no shower variations to calculate if(!ShowerHandler::currentHandlerIsSet() || 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) || 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); 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 new_ptr(FS_QTildeShowerKinematics1to2(q_,z(),phi(),pT(),this)); } 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 = QTildeKinematics::pT2_ISR(t,z(),masssquared_[2]); } while(pt2 < cutoff_->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 = QTildeKinematics::pT2_ISR(t,z(),masssquared_[2]); ptRew=pt2 < cutoff_->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 new_ptr(IS_QTildeShowerKinematics1to2(q_,z(),phi(),pT(),this)); } void SudakovFormFactor::initialize(const IdList & ids, Energy2 & tmin) { ids_=ids; tmin = 4.*cutoff_->pT2min(); masses_ = cutoff_->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 = QTildeKinematics::pT2_Decay(t,z(),masssquared_[0],masssquared_[2]); } while(SplittingFnVeto((1.-z())*t/z(),ids,true,rho,detuning)|| alphaSVeto(splittingFn()->pTScale() ? sqr(1.-z())*t : (1.-z())*t ) || pt2pT2min() || t*(1.-z())>masssquared_[0]-sqr(minmass)); if(t > ZERO) { q_ = sqrt(t); pT_ = sqrt(pt2); } else return ShoKinPtr(); phi_ = 0.; // create the ShowerKinematics object return new_ptr(Decay_QTildeShowerKinematics1to2(q_,z(),phi(),pT(),this)); } bool SudakovFormFactor::guessDecay(Energy2 &t,Energy2 tmax, Energy minmass, double enhance, double detune) { minmass = max(minmass,GeV); // previous scale Energy2 told = t; // overestimated limits on z if(tmaxpT2min()+ 0.25*sqr(masssquared_[2])/tm2)/tm +0.5*masssquared_[2]/tm2); if(zlimits_.secondpT2min()+ 0.25*sqr(masssquared_[2])/tm2)/tm +0.5*masssquared_[2]/tm2); if(t>tmax||zlimits_.secondpT2min(); // special case for gluon radiating if(ids_[0]->id()==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.+cutoff_->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); double num = B2*C*D*s012; double y1 = B*s01*(-C*(num + root) + D*denom) / denom2; double y2 = B*s01*(-C*(num - root) + D*denom) / denom2; double x1 = -(num + root ); double x2 = -(num - root ); if(denom<0.) { y1*=-1.; y2*=-1.; x1*=-1.; x2*=-1.; } return make_pair(atan2(y1,x1) + phi0,atan2(y2,x2) + phi0); } } double SudakovFormFactor::generatePhiForward(ShowerParticle & particle, const IdList & ids, ShoKinPtr kinematics, const RhoDMatrix & rho) { // no correlations, return flat phi if(ShowerHandler::currentHandlerIsSet() && ! 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 = (!ShowerHandler::currentHandlerIsSet() || 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(!ShowerHandler::currentHandlerIsSet() || 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 if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { aziMax = -m12/sqr(pik) -m22/sqr(pjk[0]+mag) +2.*pipj/pik/(pjk[0]-mag); } else assert(false); } // if spin correlations vector > wgts; if(!ShowerHandler::currentHandlerIsSet() || dynamic_ptr_cast(ShowerHandler::currentHandler())->spinCorrelations()) { // calculate the weights wgts = splittingFn()->generatePhiForward(z,t,ids,rho); } else { wgts = {{ {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(!ShowerHandler::currentHandlerIsSet() || dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==2) { aziWgt = max(ZERO,0.5/pik/Eg*(Ei-m12*Eg/pik + (pipj*Eg - Ej*pik)/dot)/aziMax); } else if(dynamic_ptr_cast(ShowerHandler::currentHandler())->softCorrelations()==1) { aziWgt = (-m12/sqr(pik) -m22/sqr(dot) +2.*pipj/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 = {{ {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; } } diff --git a/Shower/ShowerHandler.h b/Shower/ShowerHandler.h --- a/Shower/ShowerHandler.h +++ b/Shower/ShowerHandler.h @@ -1,908 +1,910 @@ // -*- C++ -*- // // ShowerHandler.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 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_ShowerHandler_H #define HERWIG_ShowerHandler_H // // This is the declaration of the ShowerHandler class. // #include "ThePEG/Handlers/EventHandler.h" #include "ThePEG/Handlers/CascadeHandler.h" #include "ShowerVariation.h" #include "Herwig/PDF/HwRemDecayer.fh" #include "ThePEG/EventRecord/RemnantParticle.fh" #include "UEBase.h" #include "PerturbativeProcess.h" #include "Herwig/MatrixElement/Matchbox/Matching/HardScaleProfile.h" #include "ShowerHandler.fh" namespace Herwig { using namespace ThePEG; /** \ingroup Shower * - * This class is the main driver of the shower: it is responsible for + * This class is the main driver of the shower: it is responsible for * the proper handling of all other specific collaborating classes * and for the storing of the produced particles in the event record. - * + * * @see \ref ShowerHandlerInterfaces "The interfaces" * * @see ThePEG::CascadeHandler * @see MPIHandler * @see HwRemDecayer */ class ShowerHandler: public CascadeHandler { public: /** * Typedef for a pair of ThePEG::RemnantParticle pointers. */ typedef pair RemPair; public: /** * Default constructor */ ShowerHandler(); /** * Destructor */ virtual ~ShowerHandler(); public: /** * The main method which manages the multiple interactions and starts * the shower by calling cascade(sub, lastXC). */ virtual void cascade(); /** * pointer to "this", the current ShowerHandler. */ static const tShowerHandlerPtr currentHandler() { - assert(currentHandler_); - return currentHandler_; + if(currentHandler_) + return currentHandler_; + else + return tShowerHandlerPtr(); } /** * pointer to "this", the current ShowerHandler. */ static bool currentHandlerIsSet() { return currentHandler_; } public: /** * Hook to allow vetoing of event after showering hard sub-process * as in e.g. MLM merging. */ virtual bool showerHardProcessVeto() const { return false; } /** * Return true, if this cascade handler will perform reshuffling from hard * process masses. */ virtual bool isReshuffling() const { return true; } - + /** * Return true, if this cascade handler will put the final state * particles to their constituent mass. If false the nominal mass is used. */ virtual bool retConstituentMasses() const { return useConstituentMasses_; } - - + + /** - * Return true, if the shower handler can generate a truncated + * Return true, if the shower handler can generate a truncated * shower for POWHEG style events generated using Matchbox */ virtual bool canHandleMatchboxTrunc() const { return false; } /** * Get the PDF freezing scale */ Energy pdfFreezingScale() const { return pdfFreezingScale_; } - + /** * Get the local PDFs. */ PDFPtr getPDFA() const {return PDFA_;} - + /** * Get the local PDFs. */ PDFPtr getPDFB() const {return PDFB_;} - + /** * Return true if currently the primary subprocess is showered. */ bool firstInteraction() const { if (!eventHandler()->currentCollision())return true; - return ( subProcess_ == + return ( subProcess_ == eventHandler()->currentCollision()->primarySubProcess() ); } /** * Return the remnant decayer. */ tHwRemDecPtr remnantDecayer() const { return remDec_; } /** * Split the hard process into production and decays * @param tagged The tagged particles from the StepHandler * @param hard The hard perturbative process * @param decay The decay particles */ void splitHardProcess(tPVector tagged, PerturbativeProcessPtr & hard, DecayProcessMap & decay) const; /** - * Information if the Showerhandler splits the hard process. + * Information if the Showerhandler splits the hard process. */ bool doesSplitHardProcess()const {return splitHardProcess_;} /** * Decay a particle. * radPhotons switches the generation of photon * radiation on/off. * Required for Dipole Shower but not QTilde Shower. */ tDMPtr decay(PerturbativeProcessPtr, DecayProcessMap & decay, bool radPhotons = false) const; - + /** * Cached lookup of decay modes. * Generator::findDecayMode() is not efficient. */ tDMPtr findDecayMode(const string & tag) const; /** * A struct to order the particles in the same way as in the DecayMode's */ - + struct ParticleOrdering { bool operator() (tcPDPtr p1, tcPDPtr p2) const; }; - + /** * A container for ordered particles required * for constructing tags for decay mode lookup. */ typedef multiset OrderedParticles; public: /** * @name Switches for initial- and final-state radiation */ //@{ /** * Switch for any radiation */ bool doRadiation() const {return doFSR_ || doISR_;} /** * Switch on or off final state radiation. */ bool doFSR() const { return doFSR_;} - + /** * Switch on or off initial state radiation. */ bool doISR() const { return doISR_;} //@} public: /** * @name Switches for scales */ //@{ /** * Return true if maximum pt should be deduced from the factorization scale */ bool hardScaleIsMuF() const { return maxPtIsMuF_; } /** * The factorization scale factor. */ double factorizationScaleFactor() const { return factorizationScaleFactor_; } - + /** * The renormalization scale factor. */ double renFac() const { return renormalizationScaleFactor_; } - + /** * The factorization scale factor. */ double facFac() const { return factorizationScaleFactor_; } - + /** * The renormalization scale factor. */ double renormalizationScaleFactor() const { return renormalizationScaleFactor_; } /** * The scale factor for the hard scale */ double hardScaleFactor() const { return hardScaleFactor_; } /** * Return true, if the phase space restrictions of the dipole shower should * be applied. */ bool restrictPhasespace() const { return restrictPhasespace_; } /** * Return profile scales */ Ptr::tptr profileScales() const { return hardScaleProfile_; } /** * Return the relevant hard scale to be used in the profile scales */ virtual Energy hardScale() const; /** * Return information about shower phase space choices */ virtual int showerPhaseSpaceOption() const { assert(false && "not implemented in general"); return -1; } //@} public: /** * Access the shower variations */ map& showerVariations() { return showerVariations_; } /** * Return the shower variations */ const map& showerVariations() const { return showerVariations_; } /** * Access the current Weights */ map& currentWeights() { return currentWeights_; } /** * Return the current Weights */ const map& currentWeights() const { return currentWeights_; } /** * Change the current reweighting factor */ void reweight(double w) { reweight_ = w; } /** * Return the current reweighting factor */ double reweight() const { return reweight_; } public : - + /** * Access to switches for spin correlations */ //@{ /** * Spin Correlations */ unsigned int spinCorrelations() const { return spinOpt_; } - + /** * Any correlations */ virtual bool correlations() const { return spinOpt_!=0; } //@} public: - /** + /** * struct that is used to catch exceptions which are thrown * due to energy conservation issues of additional scatters */ struct ExtraScatterVeto {}; - /** + /** * struct that is used to catch exceptions which are thrown * due to fact that the Shower has been invoked more than * a defined threshold on a certain configuration */ struct ShowerTriesVeto { /** variable to store the number of attempts */ const int tries; /** constructor */ ShowerTriesVeto(int t) : tries(t) {} }; 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 Functions to perform the cascade */ //@{ /** * The main method which manages the showering of a subprocess. */ virtual tPPair cascade(tSubProPtr sub, XCPtr xcomb); /** * Set up for the cascade */ - void prepareCascade(tSubProPtr sub) { - current_ = currentStep(); + void prepareCascade(tSubProPtr sub) { + current_ = currentStep(); subProcess_ = sub; - } + } /** * Boost all the particles in the collision so that the collision always occurs * in the rest frame with the incoming particles along the z axis */ void boostCollision(bool boost); //@} protected: /** * Set/unset the current shower handler */ //@{ /** * Set the current handler */ void setCurrentHandler() { currentHandler_ = tShowerHandlerPtr(this); } /** * Unset the current handler */ void unSetCurrentHandler() { currentHandler_ = tShowerHandlerPtr(); } //@} protected: /** * @name Members relating to the underlying event and MPI */ //@{ /** - * Return true if multiple parton interactions are switched on + * Return true if multiple parton interactions are switched on * and can be used for this beam setup. */ bool isMPIOn() const { return MPIHandler_ && MPIHandler_->beamOK(); } /** * Access function for the MPIHandler, it should only be called after * checking with isMPIOn. */ - tUEBasePtr getMPIHandler() const { + tUEBasePtr getMPIHandler() const { assert(MPIHandler_); - return MPIHandler_; + return MPIHandler_; } /** * Is a beam particle where hadronic structure is resolved */ bool isResolvedHadron(tPPtr); /** - * Get the remnants from the ThePEG::PartonBinInstance es and + * Get the remnants from the ThePEG::PartonBinInstance es and * do some checks. */ RemPair getRemnants(PBIPair incbins); /** * Reset the PDF's after the hard collision has been showered */ void setMPIPDFs(); //@} public: /** * Check if a particle decays in the shower * @param id The PDG code for the particle */ bool decaysInShower(long id) const { - return ( particlesDecayInShower_.find( abs(id) ) != + return ( particlesDecayInShower_.find( abs(id) ) != particlesDecayInShower_.end() ); } protected: /** * Members to handle splitting up of hard process and decays */ //@{ /** * Find decay products from the hard process and create decay processes * @param parent The parent particle * @param hard The hard process * @param decay The decay processes */ void findDecayProducts(PPtr parent, PerturbativeProcessPtr hard, DecayProcessMap & decay) const; /** * Find decay products from the hard process and create decay processes * @param parent The parent particle * @param hard The parent hard process * @param decay The decay processes */ void createDecayProcess(PPtr parent,PerturbativeProcessPtr hard, DecayProcessMap & decay) const; //@} /** * @name Functions to return information relevant to the process being showered */ //@{ /** * Return the currently used SubProcess. */ tSubProPtr currentSubProcess() const { assert(subProcess_); return subProcess_; } /** * Access to the incoming beam particles */ tPPair incomingBeams() const { return incoming_; } //@} protected: /** * Weight handling for shower variations */ //@ /** * Combine the variation weights which have been encountered */ void combineWeights(); /** * Initialise the weights in currentEvent() */ void initializeWeights(); /** * Reset the current weights */ void resetWeights(); //@} protected: /** * Return the maximum number of attempts for showering * a given subprocess. */ unsigned int maxtry() const { return maxtry_; } protected: - + /** * Parameters for the space-time model */ //@{ /** * Whether or not to include spa-cetime distances in the shower */ bool includeSpaceTime() const {return includeSpaceTime_;} /** * The minimum virtuality for the space-time model */ Energy2 vMin() const {return vMin_;} //@} protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); /** * Initialize this object. Called in the run phase just before * a run begins. */ virtual void doinitrun(); /** * Finalize this object. Called in the run phase just after a * run has ended. Used eg. to write out statistics. */ virtual void dofinish(); //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ ShowerHandler & operator=(const ShowerHandler &) = delete; private: /** * pointer to "this", the current ShowerHandler. */ static tShowerHandlerPtr currentHandler_; /** - * a MPIHandler to administer the creation of several (semihard) + * a MPIHandler to administer the creation of several (semihard) * partonic interactions. */ UEBasePtr MPIHandler_; /** * Pointer to the HwRemDecayer */ HwRemDecPtr remDec_; private: /** * Maximum tries for various stages of the showering process */ //@{ /** * Maximum number of attempts for the * main showering loop */ unsigned int maxtry_; /** * Maximum number of attempts for the regeneration of an additional * scattering, before the number of scatters is reduced. */ unsigned int maxtryMPI_; /** * Maximum number of attempts for the regeneration of an additional * hard scattering, before this event is vetoed. */ unsigned int maxtryDP_; /** * Maximum number of attempts to generate a decay */ unsigned int maxtryDecay_; //@} private: /** * Factors for the various scales */ //@{ /** * The factorization scale factor. */ double factorizationScaleFactor_; /** * The renormalization scale factor. */ double renormalizationScaleFactor_; /** * The scale factor for the hard scale */ double hardScaleFactor_; /** * True, if the phase space restrictions of the dipole shower should * be applied. */ bool restrictPhasespace_; /** * True if maximum pt should be deduced from the factorization scale */ bool maxPtIsMuF_; /** * The profile scales */ Ptr::ptr hardScaleProfile_; //@} /** * Option to include spin correlations */ unsigned int spinOpt_; private: /** * Storage of information about the current event */ //@{ /** * The incoming beam particles for the current collision */ tPPair incoming_; /** * Boost to get back to the lab */ LorentzRotation boost_; /** * Const pointer to the currently handeled ThePEG::SubProcess */ tSubProPtr subProcess_; /** * Const pointer to the current step */ tcStepPtr current_; //@} private: /** * PDFs to be used for the various stages and related parameters */ //@{ /** * The PDF freezing scale */ Energy pdfFreezingScale_; /** * PDFs to be used for the various stages and related parameters */ //@{ /** * The PDF for beam particle A. Overrides the particle's own PDF setting. */ PDFPtr PDFA_; /** * The PDF for beam particle B. Overrides the particle's own PDF setting. */ PDFPtr PDFB_; /** * The PDF for beam particle A for remnant splitting. Overrides the particle's own PDF setting. */ PDFPtr PDFARemnant_; /** * The PDF for beam particle B for remnant splitting. Overrides the particle's own PDF setting. */ PDFPtr PDFBRemnant_; /** * The MPI PDF's to be used for secondary scatters. */ pair mpipdfs_; /** * The MPI PDF's to be used for secondary scatters. */ pair rempdfs_; /** * The MPI PDF's to be used for secondary scatters. */ pair remmpipdfs_; //@} private: /** * @name Parameters for initial- and final-state radiation */ //@{ /** * Switch on or off final state radiation. */ bool doFSR_; /** * Switch on or off initial state radiation. */ bool doISR_; //@} private: /** * @name Parameters for particle decays */ //@{ /** * Whether or not to split into hard and decay trees */ bool splitHardProcess_; /** * PDG codes of the particles which decay during showering * this is fast storage for use during running */ set particlesDecayInShower_; /** * PDG codes of the particles which decay during showering * this is a vector that is interfaced so they can be changed */ vector inputparticlesDecayInShower_; //@} private: /** * Parameters for the space-time model */ //@{ /** * Whether or not to include space-time distances in the shower */ bool includeSpaceTime_; /** * The minimum virtuality for the space-time model */ Energy2 vMin_; //@} private: - + /** * Parameters for the constituent mass treatment. */ //@{ /** * True if shower should return constituent masses. */ bool useConstituentMasses_ = true; /** * Status code to tag intermediate state as after the showering; if zero no tagging will be performed */ int tagIntermediates_ = 0; //@} private: /** * Parameters relevant for reweight and variations */ //@{ /** * The shower variations */ map showerVariations_; /** * Command to add a shower variation */ string doAddVariation(string); /** * A reweighting factor applied by the showering */ double reweight_; /** * The shower variation weights */ map currentWeights_; //@} }; } #endif /* HERWIG_ShowerHandler_H */