diff --git a/Shower/Dipole/AlphaS/alpha_s.cc b/Shower/Dipole/AlphaS/alpha_s.cc --- a/Shower/Dipole/AlphaS/alpha_s.cc +++ b/Shower/Dipole/AlphaS/alpha_s.cc @@ -1,236 +1,234 @@ // -*- C++ -*- // couplings/alpha_s.cc is part of matchbox // (C) 2008 Simon Platzer -- sp@particle.uni-karlsruhe.de #include "alpha_s.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Command.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Repository/Repository.h" #include "ThePEG/PDT/ParticleData.h" using namespace matchbox; alpha_s::alpha_s() : AlphaSBase(), min_active_flavours_(3), max_active_flavours_(6), matched_(false), scale_factor_(1.), quark_masses_squared_(), lambda_squared_(), alpha_s_in_(.1176), scale_in_(91.1876*GeV), lambda_range_(1.*MeV2,1.e6*MeV2), fixed_(false) { } -alpha_s::~alpha_s() {} - void alpha_s::persistentOutput(PersistentOStream & os) const { os << min_active_flavours_ << max_active_flavours_ << matched_ << scale_factor_; for (size_t f = 0; f < 7; ++f) os << ounit(quark_masses_squared_[f],MeV2) << ounit(lambda_squared_[f],MeV2); for (size_t f = 0; f < 6; ++f) os << ounit(nfvector[f],MeV2); os << alpha_s_in_ << ounit(scale_in_,GeV) << ounit(lambda_range_.first,MeV2) << ounit(lambda_range_.second,MeV2) << fixed_; } void alpha_s::persistentInput(PersistentIStream & is, int) { is >> min_active_flavours_ >> max_active_flavours_ >> matched_ >> scale_factor_; for (size_t f = 0; f < 7; ++f) is >> iunit(quark_masses_squared_[f],MeV2) >> iunit(lambda_squared_[f],MeV2); for (size_t f = 0; f < 6; ++f) is >> iunit(nfvector[f],MeV2); is >> alpha_s_in_ >> iunit(scale_in_,GeV) >> iunit(lambda_range_.first,MeV2) >> iunit(lambda_range_.second,MeV2) >> fixed_; } AbstractClassDescription alpha_s::initalpha_s; // Definition of the static class description member. void alpha_s::Init() { static ClassDocumentation documentation ("Base class for strong coupoling as used in matchbox"); static Parameter interfacemin_active_flavours ("min_active_flavours", "Minimum number of active flavours", &alpha_s::min_active_flavours_, 3, 0, 6, true, false, Interface::limited); static Parameter interfacemax_active_flavours ("max_active_flavours", "Maximum number of active flavours", &alpha_s::max_active_flavours_, 6, 0, 6, true, false, Interface::limited); static Parameter interfaceinput_alpha_s ("input_alpha_s", "alpha_s value at input scale", &alpha_s::alpha_s_in_, .1176, 0.0, 1.0, true, false, Interface::limited); static Parameter interfaceinput_scale ("input_scale", "Input scale for alpha_s value", &alpha_s::scale_in_, GeV, 91.1876*GeV, 0.*GeV, 0.*GeV, true, false, Interface::lowerlim); static Command interfacecheck ("check", "check", &alpha_s::check, false); static Parameter interfacescale_factor ("scale_factor", "scale factor for argument", &alpha_s::scale_factor_, 1., 0.0, 100.0, true, false, Interface::limited); static Switch interfacefixed ("fixed", "", &alpha_s::fixed_, false, false, false); static SwitchOption interfacefixedYes (interfacefixed, "Yes", "", true); static SwitchOption interfacefixedNo (interfacefixed, "No", "", false); } string alpha_s::check (string args) { istringstream argin(args); double Q_low, Q_high; long n_steps; argin >> Q_low >> Q_high >> n_steps; string fname; argin >> fname; Repository::clog() << "checking alpha_s in range [" << Q_low << "," << Q_high << "] GeV in " << n_steps << " steps.\nResults are written to " << fname << "\n"; double step_width = (Q_high-Q_low)/n_steps; match_thresholds(); Repository::clog() << "threshold matching results:\n" << "(m_Q^2 -> Lambda^2) / GeV^2 for dynamic flavours in range [" << min_active_flavours_ << "," << max_active_flavours_ << "]\n"; for (size_t f = 0; f < 7; ++f) { Repository::clog() << (quark_masses_squared_[f]/GeV2) << " " << (lambda_squared_[f]/GeV2) << "\n"; } ofstream out (fname.c_str()); for (long k = 0; k <= n_steps; ++k) { Energy Q = Q_low*GeV + k*step_width*GeV; out << (Q/GeV) << " " << (operator () (Q*Q)) << "\n"; } return "alpha_s check finished"; } void alpha_s::match_thresholds () { if (matched_) return; // get the quark masses quark_masses_squared_[0] = 0.*MeV2; for (long f = 1; f < 7; ++f) { if ( quarkMasses().empty() ) quark_masses_squared_[static_cast(f)] = sqr(getParticleData(f)->mass()); else quark_masses_squared_[static_cast(f)] = sqr(quarkMasses()[static_cast(f-1)]); } if ( quark_masses_squared_[1] > quark_masses_squared_[2] ) swap(quark_masses_squared_[1],quark_masses_squared_[2]); unsigned int active_at_input = active_flavours(sqr(scale_in_)); // solve for input lambda solve_input_lambda input_equation (this,active_at_input,alpha_s_in_,sqr(scale_in_)); gsl::bisection_root_solver,100> input_solver(input_equation); lambda_squared_[active_at_input] = MeV2 * input_solver.solve({lambda_range_.first/MeV2,lambda_range_.second/MeV2}); // get lambdas down to min active flavours unsigned int below = active_at_input; while (below > min_active_flavours_) { solve_lambda_below match_equation (this,below, lambda_squared_[below], quark_masses_squared_[below]); gsl::bisection_root_solver,100> match_solver(match_equation); lambda_squared_[below-1] = MeV2 * match_solver.solve({lambda_range_.first/MeV2,lambda_range_.second/MeV2}); --below; } // get lambdas up to max active flavours unsigned int above = active_at_input; while (above < max_active_flavours_) { solve_lambda_above match_equation (this,above, lambda_squared_[above], quark_masses_squared_[above+1]); gsl::bisection_root_solver,100> match_solver(match_equation); lambda_squared_[above+1] = MeV2 *match_solver.solve({lambda_range_.first/MeV2,lambda_range_.second/MeV2}); ++above; } if (min_active_flavours_ > 0) { for (size_t f = 0; f < min_active_flavours_; ++f) { lambda_squared_[f] = lambda_squared_[min_active_flavours_]; } } if (max_active_flavours_ < 6) { for (size_t f = max_active_flavours_+1; f < 7; ++f) { lambda_squared_[f] = lambda_squared_[max_active_flavours_]; } } matched_ = true; return; } diff --git a/Shower/Dipole/AlphaS/alpha_s.h b/Shower/Dipole/AlphaS/alpha_s.h --- a/Shower/Dipole/AlphaS/alpha_s.h +++ b/Shower/Dipole/AlphaS/alpha_s.h @@ -1,352 +1,344 @@ // -*- C++ -*- // couplings/alpha_s.h is part of matchbox // (C) 2008 Simon Platzer -- sp@particle.uni-karlsruhe.de #ifndef matchbox_couplings_alpha_s_h #define matchbox_couplings_alpha_s_h #include #include #include "ThePEG/Interface/Interfaced.h" #include "ThePEG/StandardModel/AlphaSBase.h" #include "gsl.h" namespace matchbox { using namespace ThePEG; template struct solve_lambda_below { typedef AlphaS alpha_s; inline solve_lambda_below (alpha_s* a, unsigned int n, Energy2 lambda2n, Energy2 mass2) : alpha(a), nf_in(n), lambda2_nf_in(lambda2n), threshold(mass2) {} alpha_s * alpha; unsigned int nf_in; Energy2 lambda2_nf_in; Energy2 threshold; inline double operator () (double lambda2) { return ((*alpha)(threshold,lambda2_nf_in,nf_in) - (*alpha)(threshold,lambda2*MeV2,nf_in-1)); } }; template struct solve_lambda_above { typedef AlphaS alpha_s; inline solve_lambda_above (alpha_s * a, unsigned int n, Energy2 lambda2n, Energy2 mass2) : alpha(a), nf_in(n), lambda2_nf_in(lambda2n), threshold(mass2) {} alpha_s * alpha; unsigned int nf_in; Energy2 lambda2_nf_in; Energy2 threshold; inline double operator () (double lambda2) { return ((*alpha)(threshold,lambda2_nf_in,nf_in) - (*alpha)(threshold,lambda2*MeV2,nf_in+1)); } }; template struct solve_input_lambda { typedef AlphaS alpha_s; inline solve_input_lambda (alpha_s * a, unsigned int n, double inalpha, Energy2 inscale) : alpha(a), nf_in(n), alpha_in(inalpha), scale_in(inscale) {} alpha_s * alpha; unsigned int nf_in; double alpha_in; Energy2 scale_in; inline double operator () (double lambda2) { return ((*alpha)(scale_in,lambda2*MeV2,nf_in) - alpha_in); } }; /** * Base class for the strong coupling. * * @see \ref alpha_sInterfaces "The interfaces" * defined for alpha_s. */ class alpha_s : public AlphaSBase { public: - /** @name Standard constructors and destructors. */ - //@{ /** * The default constructor. */ alpha_s(); - - /** - * The destructor. - */ - virtual ~alpha_s(); - //@} public: /** @name Virtual functions as required by AlphaSBase. */ //@{ /** * The \f$\alpha_S\f$. Return the QCD coupling for a given \a scale * using the given standard model object \a sm. */ virtual inline double value(Energy2 scale, const StandardModelBase &) const { return operator() (scale); } /** * Return the flavour thresholds used. The returned vector contains * (in position i) the scales when the active number of * flavours changes from i to i+1. */ virtual inline vector flavourThresholds() const { assert(!nfvector.empty()); return nfvector; } /** * Return the \f$\Lambda_{QCD}\f$ used for different numbers of * active flavours. */ virtual inline vector LambdaQCDs() const { vector res; for (size_t k = 0; k < 7; ++k) res.push_back(sqrt(lambda_squared_[k])); return res; } //@} public: /// return alpha_s as function of scale inline double operator () (Energy2 scale) const { if ( fixed_ ) return alpha_s_in_; assert(matched()); unsigned int active = active_flavours(scale_factor_*scale); return operator () (scale_factor_*scale,lambda_squared_[active],active); } /// return alpha_s as function of scale, QCD scale /// and number of active flavours virtual double operator () (Energy2 scale, Energy2 lambda2, unsigned int nf) const = 0; /// match thresholds and write alpha_s /// to specified file; arguments are /// Q_low/GeV Q_high/GeV n_steps filename string check (string args); public: /// return minimum number of active flavours inline unsigned int min_active_flavours () const { return min_active_flavours_; } /// set minimum number of active flavours inline void min_active_flavours (unsigned int nf) { min_active_flavours_ = nf; } /// return maximum number of active flavours inline unsigned int max_active_flavours () const { return max_active_flavours_; } /// set maximum number of active flavours inline void max_active_flavours (unsigned int nf) { max_active_flavours_ = nf; } /// return the number of active flavours at the given scale inline unsigned int active_flavours (Energy2 scale) const { unsigned int active = 0; if (scale > 0.*GeV2) { while(quark_mass_squared(active) < scale) { if (++active == max_active_flavours_+1) break; } active -= 1; } else { active = 0; } return active; } /// return the lambda squared for the given number of flavours inline Energy2 lambda_squared (unsigned int f) const { assert(f < 7); return lambda_squared_[f]; } /// return the mass squared for given flavour inline Energy2 quark_mass_squared (unsigned int f) const { assert(f < 7); return quark_masses_squared_[f]; } /// set the mass squared for given flavour inline void quark_mass_squared (unsigned int f, Energy2 m2) { assert(f < 7); quark_masses_squared_[f] = m2; matched_ = false; } public: /// perform the threshold matching /// given alpha_s value at reference scale void match_thresholds (); /// return true, if threshold matching has been /// performed inline bool matched () const { return matched_; } 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 inline void doinit() { match_thresholds(); copy(quark_masses_squared_.begin()+1, quark_masses_squared_.end(),nfvector.begin()); AlphaSBase::doinit(); } //@} /// return the scale factor double scale_factor () const { return scale_factor_; } public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @name os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @name is the persistent input stream read from. * @name 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(); private: /** * The static object used to initialize the description of this class. * Indicates that this is an abstract class with persistent data. */ static AbstractClassDescription initalpha_s; /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ alpha_s & operator=(const alpha_s &) = delete; private: unsigned int min_active_flavours_; unsigned int max_active_flavours_; bool matched_; double scale_factor_; std::array quark_masses_squared_; std::array lambda_squared_; vector nfvector=vector(6); double alpha_s_in_; Energy scale_in_; pair lambda_range_; bool fixed_; }; } #include "ThePEG/Utilities/ClassTraits.h" namespace ThePEG { /** @cond TRAITSPECIALIZATIONS */ /** This template specialization informs ThePEG about the * base classes of alpha_s. */ template <> struct BaseClassTrait { /** Typedef of the first base class of alpha_s. */ typedef AlphaSBase NthBase; }; /** This template specialization informs ThePEG about the name of * the alpha_s class and the shared object where it is defined. */ template <> struct ClassTraits : public ClassTraitsBase { /** Return a platform-independent class name */ static string className() { return "matchbox::alpha_s"; } /** * The name of a file containing the dynamic library where the class * alpha_s is implemented. It may also include several, space-separated, * libraries if the class alpha_s depends on other classes (base classes * excepted). In this case the listed libraries will be dynamically * linked in the order they are specified. */ static string library() { return "HwDipoleShowerAlphaS.so"; } }; /** @endcond */ } #endif /* matchbox_couplings_alpha_s_h */ diff --git a/Shower/Dipole/AlphaS/lo_alpha_s.cc b/Shower/Dipole/AlphaS/lo_alpha_s.cc --- a/Shower/Dipole/AlphaS/lo_alpha_s.cc +++ b/Shower/Dipole/AlphaS/lo_alpha_s.cc @@ -1,67 +1,65 @@ // -*- C++ -*- // couplings/lo_alpha_s.cc is part of matchbox // (C) 2008 Simon Platzer -- sp@particle.uni-karlsruhe.de #include "lo_alpha_s.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" using namespace matchbox; lo_alpha_s::lo_alpha_s() : alpha_s(), freezing_scale_(1.*GeV) {} -lo_alpha_s::~lo_alpha_s() {} - IBPtr lo_alpha_s::clone() const { return new_ptr(*this); } IBPtr lo_alpha_s::fullclone() const { return new_ptr(*this); } void lo_alpha_s::persistentOutput(PersistentOStream & os) const { os << ounit(freezing_scale_,GeV); } void lo_alpha_s::persistentInput(PersistentIStream & is, int) { is >> iunit(freezing_scale_,GeV); } ClassDescription lo_alpha_s::initlo_alpha_s; // Definition of the static class description member. void lo_alpha_s::Init() { static ClassDocumentation documentation ("LO running alpha_s"); static Parameter interfacefreezing_scale ("freezing_scale", "Freeze alpha_s below given scale", &lo_alpha_s::freezing_scale_, GeV, 1.0*GeV, 0.0*GeV, 0*GeV, true, false, Interface::lowerlim); } double lo_alpha_s::operator () (Energy2 scale, Energy2 lambda2, unsigned int nf) const { if (scale < sqr(freezing_scale_)) { scale = sqr(freezing_scale_); nf = active_flavours(scale); lambda2 = lambda_squared(nf); } double beta0 = (33.-2.*nf)/(12.*Constants::pi); return 1./(beta0*log(scale/lambda2)); } diff --git a/Shower/Dipole/AlphaS/lo_alpha_s.h b/Shower/Dipole/AlphaS/lo_alpha_s.h --- a/Shower/Dipole/AlphaS/lo_alpha_s.h +++ b/Shower/Dipole/AlphaS/lo_alpha_s.h @@ -1,167 +1,159 @@ // -*- C++ -*- // couplings/lo_alpha_s.h is part of matchbox // (C) 2008 Simon Platzer -- sp@particle.uni-karlsruhe.de #ifndef matchbox_couplings_lo_alpha_s_h #define matchbox_couplings_lo_alpha_s_h #include "alpha_s.h" namespace matchbox { using namespace ThePEG; /** * LO running alpha_s * * @see \ref lo_alpha_sInterfaces "The interfaces" * defined for lo_alpha_s. */ class lo_alpha_s : public alpha_s { public: - /** @name Standard constructors and destructors. */ - //@{ /** * The default constructor. */ lo_alpha_s(); - /** - * The destructor. - */ - virtual ~lo_alpha_s(); - //@} - public: /// return alpha_s as function of scale, QCD scale /// and number of active flavours virtual double operator () (Energy2 scale, Energy2 lambda2, unsigned int nf) const; /// return the number of loops which determine this running virtual unsigned int nloops () const { return 1; } public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @name os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @name is the persistent input stream read from. * @name 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 inline void doinit() { freezing_scale_ *= scale_factor(); alpha_s::doinit(); } //@} 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; //@} private: /** * The static object used to initialize the description of this class. * Indicates that this is an abstract class with persistent data. */ static ClassDescription initlo_alpha_s; /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ lo_alpha_s & operator=(const lo_alpha_s &) = delete; private: Energy freezing_scale_; }; } #include "ThePEG/Utilities/ClassTraits.h" namespace ThePEG { /** @cond TRAITSPECIALIZATIONS */ /** This template specialization informs ThePEG about the * base classes of lo_alpha_s. */ template <> struct BaseClassTrait { /** Typedef of the first base class of lo_alpha_s. */ typedef matchbox::alpha_s NthBase; }; /** This template specialization informs ThePEG about the name of * the lo_alpha_s class and the shared object where it is defined. */ template <> struct ClassTraits : public ClassTraitsBase { /** Return a platform-independent class name */ static string className() { return "matchbox::lo_alpha_s"; } /** * The name of a file containing the dynamic library where the class * lo_alpha_s is implemented. It may also include several, space-separated, * libraries if the class lo_alpha_s depends on other classes (base classes * excepted). In this case the listed libraries will be dynamically * linked in the order they are specified. */ static string library() { return "HwDipoleShowerAlphaS.so"; } }; /** @endcond */ } #endif /* matchbox_couplings_lo_alpha_s_h */ diff --git a/Shower/Dipole/AlphaS/nlo_alpha_s.cc b/Shower/Dipole/AlphaS/nlo_alpha_s.cc --- a/Shower/Dipole/AlphaS/nlo_alpha_s.cc +++ b/Shower/Dipole/AlphaS/nlo_alpha_s.cc @@ -1,135 +1,133 @@ // -*- C++ -*- // couplings/nlo_alpha_s.cc is part of matchbox // (C) 2008 Simon Platzer -- sp@particle.uni-karlsruhe.de #include "nlo_alpha_s.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" using namespace matchbox; nlo_alpha_s::nlo_alpha_s() : alpha_s(), freezing_scale_(1.*GeV), exact_evaluation_(true), two_largeq_terms_(true) { } -nlo_alpha_s::~nlo_alpha_s() {} - IBPtr nlo_alpha_s::clone() const { return new_ptr(*this); } IBPtr nlo_alpha_s::fullclone() const { return new_ptr(*this); } void nlo_alpha_s::persistentOutput(PersistentOStream & os) const { os << ounit(freezing_scale_,GeV) << exact_evaluation_ << two_largeq_terms_; } void nlo_alpha_s::persistentInput(PersistentIStream & is, int) { is >> iunit(freezing_scale_,GeV) >> exact_evaluation_ >> two_largeq_terms_; } ClassDescription nlo_alpha_s::initnlo_alpha_s; // Definition of the static class description member. void nlo_alpha_s::Init() { static ClassDocumentation documentation ("NLO running alpha_s"); static Parameter interfacefreezing_scale ("freezing_scale", "Freeze alpha_s below given scale", &nlo_alpha_s::freezing_scale_, GeV, 1.0*GeV, 0.0*GeV, 0*GeV, true, false, Interface::lowerlim); static Switch interfaceexact_evaluation ("exact_evaluation", "Wether to exactly evaluate the running or use running for large scales", &nlo_alpha_s::exact_evaluation_, true, true, false); static SwitchOption interfaceexact_evaluationexact (interfaceexact_evaluation, "exact", "Perform exact evaluation", true); static SwitchOption interfaceexact_evaluationlarge_scale (interfaceexact_evaluation, "large_scale", "Perform approximate evaluation for large scales", false); static Switch interfacetwo_largeq_terms ("two_largeq_terms", "Include two terms in the large q expansion.", &nlo_alpha_s::two_largeq_terms_, true, false, false); static SwitchOption interfacetwo_largeq_termsYes (interfacetwo_largeq_terms, "Yes", "Include two terms.", true); static SwitchOption interfacetwo_largeq_termsNo (interfacetwo_largeq_terms, "No", "Only include one term.", false); } double nlo_alpha_s::operator () (Energy2 scale, Energy2 lambda2, unsigned int nf) const { if (scale < sqr(freezing_scale_)) { scale = sqr(freezing_scale_); nf = active_flavours(scale); lambda2 = lambda_squared(nf); } double beta0 = (33.-2.*nf)/(12.*Constants::pi); double beta1 = (153.-19.*nf)/(24.*sqr(Constants::pi)); if (exact_evaluation_) { rg_solver().f.slog = log(scale/lambda2); rg_solver().f.nf = nf; double slog = rg_solver().f.slog; double center = (1./(beta0*slog))* (1. - (beta1/sqr(beta0)) * log(slog)/slog + sqr(beta1/(sqr(beta0)*slog)) * (sqr(log(slog)-.5) - 5./4.)); return rg_solver().solve({.5*center,1.5*center}); } else { double slog = log(scale/lambda2); double res = (1./(beta0*slog))* (1. - (beta1/sqr(beta0)) * log(slog)/slog); if ( two_largeq_terms_ ) res += (1./(beta0*slog))* (sqr(beta1/(sqr(beta0)*slog)) * (sqr(log(slog)-.5) - 5./4.)); return res; } return 0.; } diff --git a/Shower/Dipole/AlphaS/nlo_alpha_s.h b/Shower/Dipole/AlphaS/nlo_alpha_s.h --- a/Shower/Dipole/AlphaS/nlo_alpha_s.h +++ b/Shower/Dipole/AlphaS/nlo_alpha_s.h @@ -1,197 +1,189 @@ // -*- C++ -*- // couplings/nlo_alpha_s.h is part of matchbox // (C) 2008 Simon Platzer -- sp@particle.uni-karlsruhe.de #ifndef matchbox_couplings_nlo_alpha_s_h #define matchbox_couplings_nlo_alpha_s_h #include "alpha_s.h" namespace matchbox { using namespace ThePEG; /** * NLO running alpha_s * * @see \ref nlo_alpha_sInterfaces "The interfaces" * defined for nlo_alpha_s. */ class nlo_alpha_s : public alpha_s { public: - /** @name Standard constructors and destructors. */ - //@{ /** * The default constructor. */ nlo_alpha_s(); - /** - * The destructor. - */ - virtual ~nlo_alpha_s(); - //@} - public: /// return alpha_s as function of scale, QCD scale /// and number of active flavours virtual double operator () (Energy2 scale, Energy2 lambda2, unsigned int nf) const; /// return the number of loops which determine this running virtual unsigned int nloops () const { return 2; } public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @name os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @name is the persistent input stream read from. * @name 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 inline void doinit() { freezing_scale_ *= scale_factor(); alpha_s::doinit(); } //@} 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; //@} private: /** * The static object used to initialize the description of this class. * Indicates that this is an abstract class with persistent data. */ static ClassDescription initnlo_alpha_s; /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ nlo_alpha_s & operator=(const nlo_alpha_s &) = delete; private: struct rg_solution { inline double operator () (double alpha) { double beta0 = (33.-2.*nf)/(12.*Constants::pi); double beta1 = (153.-19.*nf)/(24.*sqr(Constants::pi)); return ((1./alpha)+(beta1/beta0)*log(alpha/(beta0+beta1*alpha))- beta0*slog); } double slog; unsigned int nf; }; Energy freezing_scale_; bool exact_evaluation_; static rg_solution& rg () { static rg_solution rg_; return rg_; } static gsl::bisection_root_solver& rg_solver () { static gsl::bisection_root_solver rg_solver_(rg()); return rg_solver_; } bool two_largeq_terms_; }; } #include "ThePEG/Utilities/ClassTraits.h" namespace ThePEG { /** @cond TRAITSPECIALIZATIONS */ /** This template specialization informs ThePEG about the * base classes of nlo_alpha_s. */ template <> struct BaseClassTrait { /** Typedef of the first base class of nlo_alpha_s. */ typedef matchbox::alpha_s NthBase; }; /** This template specialization informs ThePEG about the name of * the nlo_alpha_s class and the shared object where it is defined. */ template <> struct ClassTraits : public ClassTraitsBase { /** Return a platform-independent class name */ static string className() { return "matchbox::nlo_alpha_s"; } /** * The name of a file containing the dynamic library where the class * nlo_alpha_s is implemented. It may also include several, space-separated, * libraries if the class nlo_alpha_s depends on other classes (base classes * excepted). In this case the listed libraries will be dynamically * linked in the order they are specified. */ static string library() { return "HwDipoleShowerAlphaS.so"; } }; /** @endcond */ } #endif /* matchbox_couplings_nlo_alpha_s_h */ diff --git a/Shower/Dipole/Utility/ConstituentReshuffler.cc b/Shower/Dipole/Utility/ConstituentReshuffler.cc --- a/Shower/Dipole/Utility/ConstituentReshuffler.cc +++ b/Shower/Dipole/Utility/ConstituentReshuffler.cc @@ -1,622 +1,617 @@ // -*- C++ -*- // // ConstituentReshuffler.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. // // // This is the implementation of the non-inlined, non-templated member // functions of the ConstituentReshuffler class. // #include #include "ConstituentReshuffler.h" #include "ThePEG/Interface/ClassDocumentation.h" #include #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "DipolePartonSplitter.h" #include "Herwig/Utilities/GSLBisection.h" #include "Herwig/Shower/Dipole/DipoleShowerHandler.h" #include "Herwig/Shower/ShowerHandler.h" using namespace Herwig; -ConstituentReshuffler::ConstituentReshuffler() - : HandlerBase() {} - -ConstituentReshuffler::~ConstituentReshuffler() {} - IBPtr ConstituentReshuffler::clone() const { return new_ptr(*this); } IBPtr ConstituentReshuffler::fullclone() const { return new_ptr(*this); } void ConstituentReshuffler::reshuffle(PList& out, PPair& in, PList& intermediates, const bool decay, PList& decayPartons, PList& decayRecoilers) { assert(ShowerHandler::currentHandler()->retConstituentMasses()); if ( !decay ) { if (out.size() == 0) return; if (out.size() == 1) { PPtr recoiler; PPtr parton = out.front(); if (DipolePartonSplitter::colourConnected(parton,in.first) && DipolePartonSplitter::colourConnected(parton,in.second)) { if (UseRandom::rnd() < .5) recoiler = in.first; else recoiler = in.second; } else if (DipolePartonSplitter::colourConnected(parton,in.first)) { recoiler = in.first; } else if (DipolePartonSplitter::colourConnected(parton,in.second)) { recoiler = in.second; } else assert(false); assert(abs(recoiler->momentum().vect().perp2()/GeV2) < 1e-6); double sign = recoiler->momentum().z() < 0.*GeV ? -1. : 1.; Energy2 qperp2 = parton->momentum().perp2(); if (qperp2/GeV2 < Constants::epsilon) { // no emission off a 2 -> singlet process which // needed a single forced splitting: should never happen (?) assert(false); throw Veto(); } Energy2 m2 = sqr(parton->dataPtr()->constituentMass()); Energy abs_q = parton->momentum().vect().mag(); Energy qz = parton->momentum().z(); Energy abs_pz = recoiler->momentum().t(); assert(abs_pz > 0.*GeV); Energy xi_pz = sign*(2.*qperp2*abs_pz + m2*(abs_q + sign*qz))/(2.*qperp2); Energy x_qz = (2.*qperp2*qz + m2*(qz+sign*abs_q))/(2.*qperp2); Lorentz5Momentum recoiler_momentum (0.*GeV,0.*GeV,xi_pz,xi_pz < 0.*GeV ? - xi_pz : xi_pz); recoiler_momentum.rescaleMass(); Lorentz5Momentum parton_momentum (parton->momentum().x(),parton->momentum().y(),x_qz,sqrt(m2+qperp2+x_qz*x_qz)); parton_momentum.rescaleMass(); PPtr n_parton = new_ptr(Particle(parton->dataPtr())); n_parton->set5Momentum(parton_momentum); DipolePartonSplitter::change(parton,n_parton,false); out.pop_front(); intermediates.push_back(parton); out.push_back(n_parton); PPtr n_recoiler = new_ptr(Particle(recoiler->dataPtr())); n_recoiler->set5Momentum(recoiler_momentum); DipolePartonSplitter::change(recoiler,n_recoiler,true); intermediates.push_back(recoiler); if (recoiler == in.first) { in.first = n_recoiler; } if (recoiler == in.second) { in.second = n_recoiler; } return; } } Energy zero (0.*GeV); Lorentz5Momentum Q (zero,zero,zero,zero); for (PList::iterator p = out.begin(); p != out.end(); ++p) { Q += (**p).momentum(); } Boost beta = Q.findBoostToCM(); list mbackup; bool need_boost = (beta.mag2() > Constants::epsilon); if (need_boost) { for (PList::iterator p = out.begin(); p != out.end(); ++p) { Lorentz5Momentum mom = (**p).momentum(); mbackup.push_back(mom); (**p).set5Momentum(mom.boost(beta)); } } double xi; // Only partons if ( decayRecoilers.size()==0 ) { list masses; for ( auto p : out ) masses.push_back(p->dataPtr()->constituentMass()); ReshuffleEquation::const_iterator> solve (Q.m(),out.begin(),out.end(), masses.begin(),masses.end()); GSLBisection solver(1e-10,1e-8,10000); try { xi = solver.value(solve,0.0,1.1); } catch (GSLBisection::GSLerror) { throw DipoleShowerHandler::RedoShower(); } catch (GSLBisection::IntervalError) { throw DipoleShowerHandler::RedoShower(); } } // Partons and decaying recoilers else { DecayReshuffleEquation solve (Q.m(),decayPartons.begin(),decayPartons.end(),decayRecoilers.begin(),decayRecoilers.end()); GSLBisection solver(1e-10,1e-8,10000); try { xi = solver.value(solve,0.0,1.1); } catch (GSLBisection::GSLerror) { throw DipoleShowerHandler::RedoShower(); } catch (GSLBisection::IntervalError) { throw DipoleShowerHandler::RedoShower(); } } PList reshuffled; list::const_iterator backup_it; if (need_boost) backup_it = mbackup.begin(); // Reshuffling of non-decaying partons only if ( decayRecoilers.size()==0 ) { for (PList::iterator p = out.begin(); p != out.end(); ++p) { PPtr rp = new_ptr(Particle((**p).dataPtr())); DipolePartonSplitter::change(*p,rp,false); Lorentz5Momentum rm; rm = Lorentz5Momentum (xi*(**p).momentum().x(), xi*(**p).momentum().y(), xi*(**p).momentum().z(), sqrt(sqr((**p).dataPtr()->constituentMass()) + xi*xi*(sqr((**p).momentum().t())-sqr((**p).dataPtr()->mass())))); rm.rescaleMass(); if (need_boost) { (**p).set5Momentum(*backup_it); ++backup_it; rm.boost(-beta); } rp->set5Momentum(rm); intermediates.push_back(*p); reshuffled.push_back(rp); } } // For the case of a decay process with non-partonic recoilers else { assert ( decay ); for (PList::iterator p = out.begin(); p != out.end(); ++p) { // Flag to update spinInfo bool updateSpin = false; PPtr rp = new_ptr(Particle((**p).dataPtr())); DipolePartonSplitter::change(*p,rp,false); Lorentz5Momentum rm; // If the particle is a parton and not a recoiler if ( find( decayRecoilers.begin(), decayRecoilers.end(), *p ) == decayRecoilers.end() ) { rm = Lorentz5Momentum (xi*(**p).momentum().x(), xi*(**p).momentum().y(), xi*(**p).momentum().z(), sqrt(sqr((**p).dataPtr()->constituentMass()) + xi*xi*(sqr((**p).momentum().t())-sqr((**p).dataPtr()->mass())))); } // Otherwise the parton is a recoiler // and its invariant mass must be preserved else { if ( (*p)-> spinInfo() ) updateSpin = true; rm = Lorentz5Momentum (xi*(**p).momentum().x(), xi*(**p).momentum().y(), xi*(**p).momentum().z(), sqrt(sqr((**p).momentum().m()) + xi*xi*(sqr((**p).momentum().t())-sqr((**p).momentum().m())))); } rm.rescaleMass(); if (need_boost) { (**p).set5Momentum(*backup_it); ++backup_it; rm.boost(-beta); } rp->set5Momentum(rm); // Update SpinInfo if required if ( updateSpin ) updateSpinInfo(*p, rp); intermediates.push_back(*p); reshuffled.push_back(rp); } } out.clear(); out.splice(out.end(),reshuffled); } void ConstituentReshuffler::hardProcDecayReshuffle(PList& decaying, PList& eventOutgoing, PList& eventHard, PPair& eventIncoming, PList& eventIntermediates) { // Note, when this function is called, the particle pointers // in theDecays/decaying are those prior to the showering. // Here we find the newest pointers in the outgoing. // The update of the PPtrs in theDecays is done in DipoleShowerHandler::constituentReshuffle() // as this needs to be done if ConstituentReshuffling is switched off. //Make sure the shower should return constituent masses: assert(ShowerHandler::currentHandler()->retConstituentMasses()); // Find the outgoing decaying particles PList recoilers; for ( PList::iterator decIt = decaying.begin(); decIt != decaying.end(); ++decIt) { // First find the particles in the intermediates PList::iterator pos = find(eventIntermediates.begin(),eventIntermediates.end(), *decIt); // Colourless particle or coloured particle that did not radiate. if(pos==eventIntermediates.end()) { // Check that this is not a particle from a subsequent decay. // e.g. the W from a top decay from an LHE file. if ( find( eventHard.begin(), eventHard.end(), *decIt ) == eventHard.end() && find( eventOutgoing.begin(), eventOutgoing.end(), *decIt ) == eventOutgoing.end() ) continue; else recoilers.push_back( *decIt ); } // Coloured decaying particle that radiated else { PPtr unstable = *pos; while(!unstable->children().empty()) { unstable = unstable->children()[0]; } assert( find( eventOutgoing.begin(),eventOutgoing.end(), unstable ) != eventOutgoing.end() ); recoilers.push_back( unstable ); } } // Make a list of partons PList partons; for ( PList::iterator outPos = eventOutgoing.begin(); outPos != eventOutgoing.end(); ++outPos ) { if ( find (recoilers.begin(), recoilers.end(), *outPos ) == recoilers.end() ) { partons.push_back( *outPos ); } } // If no outgoing partons, do nothing if ( partons.size() == 0 ){ return; } // Otherwise reshuffling needs to be done. // If there is only one parton, attempt to reshuffle with // the incoming to be consistent with the reshuffle for a // hard process with no decays. else if ( partons.size() == 1 && ( DipolePartonSplitter::colourConnected(partons.front(),eventIncoming.first) || DipolePartonSplitter::colourConnected(partons.front(),eventIncoming.second) ) ) { // Erase the parton from the event outgoing eventOutgoing.erase( find( eventOutgoing.begin(), eventOutgoing.end(), partons.front() ) ); // Perform the reshuffle, this update the intermediates and the incoming reshuffle(partons, eventIncoming, eventIntermediates); // Update the outgoing eventOutgoing.push_back(partons.front()); return; } // If reshuffling amongst the incoming is not possible // or if we have multiple outgoing partons. else { // Create a complete list of the outgoing from the process PList out; // Make an empty list for storing the new intermediates PList intermediates; // Empty incoming particles pair PPair in; // A single parton which cannot be reshuffled // with the incoming. if ( partons.size() == 1 ) { // Populate the out for the reshuffling out.insert(out.end(),partons.begin(),partons.end()); out.insert(out.end(),recoilers.begin(),recoilers.end()); assert( out.size() > 1 ); // Perform the reshuffle with the temporary particle lists reshuffle(out, in, intermediates, true, partons, recoilers); } // If there is more than one parton, reshuffle only // amongst the partons else { assert(partons.size() > 1); // Populate the out for the reshuffling out.insert(out.end(),partons.begin(),partons.end()); assert( out.size() > 1 ); // Perform the reshuffle with the temporary particle lists reshuffle(out, in, intermediates, true); } // Update the dipole event record updateEvent(intermediates, eventIntermediates, out, eventOutgoing, eventHard ); return; } } void ConstituentReshuffler::decayReshuffle(PerturbativeProcessPtr& decayProc, PList& eventOutgoing, PList& eventHard, PList& eventIntermediates ) { // Separate particles into those to be assigned constituent masses // i.e. non-decaying coloured partons // and those which must only absorb recoil // i.e. non-coloured and decaying particles PList partons; PList recoilers; //Make sure the shower should return constituent masses: assert(ShowerHandler::currentHandler()->retConstituentMasses()); // Populate the particle lists from the outgoing of the decay process for( unsigned int ix = 0; ixoutgoing().size(); ++ix) { // Identify recoilers if ( !decayProc->outgoing()[ix].first->coloured() || ShowerHandler::currentHandler()->decaysInShower(decayProc->outgoing()[ix].first->id() ) ) recoilers.push_back(decayProc->outgoing()[ix].first); else partons.push_back(decayProc->outgoing()[ix].first); } // If there are no outgoing partons, then no reshuffling // needs to be done if ( partons.size() == 0 ) return; // Reshuffling needs to be done: else { // Create a complete list of the outgoing from the process PList out; // Make an empty list for storing the new intermediates PList intermediates; // Empty incoming particles pair PPair in; // SW - 15/06/2018, 31/01/2019 - Always include 'recoilers' in // reshuffling, regardless of the number of partons to be put on their // constituent mass shell. This is because reshuffling between 2 partons // frequently leads to a redoShower exception. This treatment is // consistent with the AO shower // Populate the out for the reshuffling out.insert(out.end(),partons.begin(),partons.end()); out.insert(out.end(),recoilers.begin(),recoilers.end()); assert( out.size() > 1 ); // Perform the reshuffle with the temporary particle lists reshuffle(out, in, intermediates, true, partons, recoilers); // Update the dipole event record and the decay process updateEvent(intermediates, eventIntermediates, out, eventOutgoing, eventHard, decayProc ); return; } } void ConstituentReshuffler::updateEvent( PList& intermediates, PList& eventIntermediates, #ifndef NDEBUG PList& out, #else PList&, #endif PList& eventOutgoing, PList& eventHard, PerturbativeProcessPtr decayProc ) { // Loop over the new intermediates following the reshuffling for (PList::iterator p = intermediates.begin(); p != intermediates.end(); ++p) { // Update the event record intermediates eventIntermediates.push_back(*p); // Identify the reshuffled particle assert( (*p)->children().size()==1 ); PPtr reshuffled = (*p)->children()[0]; assert( find(out.begin(), out.end(), reshuffled) != out.end() ); // Update the event record outgoing PList::iterator posOut = find(eventOutgoing.begin(), eventOutgoing.end(), *p); if ( posOut != eventOutgoing.end() ) { eventOutgoing.erase(posOut); eventOutgoing.push_back(reshuffled); } else { PList::iterator posHard = find(eventHard.begin(), eventHard.end(), *p); assert( posHard != eventHard.end() ); eventHard.erase(posHard); eventHard.push_back(reshuffled); } // Replace the particle in the the decay process outgoing if ( decayProc ) { vector >::iterator decayOutIt = decayProc->outgoing().end(); for ( decayOutIt = decayProc->outgoing().begin(); decayOutIt!= decayProc->outgoing().end(); ++decayOutIt ) { if ( decayOutIt->first == *p ){ break; } } assert( decayOutIt != decayProc->outgoing().end() ); decayOutIt->first = reshuffled; } } } void ConstituentReshuffler::updateSpinInfo( PPtr& oldPart, PPtr& newPart ) { const Lorentz5Momentum& oldMom = oldPart->momentum(); const Lorentz5Momentum& newMom = newPart->momentum(); // Rotation from old momentum to +ve z-axis LorentzRotation oldToZAxis; Axis axisOld(oldMom.vect().unit()); if( axisOld.perp2() > 1e-12 ) { double sinth(sqrt(1.-sqr(axisOld.z()))); oldToZAxis.rotate( -acos(axisOld.z()),Axis(-axisOld.y()/sinth,axisOld.x()/sinth,0.)); } // Rotation from new momentum to +ve z-axis LorentzRotation newToZAxis; Axis axisNew(newMom.vect().unit()); if( axisNew.perp2() > 1e-12 ) { double sinth(sqrt(1.-sqr(axisNew.z()))); newToZAxis.rotate( -acos(axisNew.z()),Axis(-axisNew.y()/sinth,axisNew.x()/sinth,0.)); } // Boost from old momentum to new momentum along z-axis Lorentz5Momentum momOldRotated = oldToZAxis*Lorentz5Momentum(oldMom); Lorentz5Momentum momNewRotated = newToZAxis*Lorentz5Momentum(newMom); Energy2 a = sqr(momOldRotated.z()) + sqr(momNewRotated.t()); Energy2 b = 2.*momOldRotated.t()*momOldRotated.z(); Energy2 c = sqr(momOldRotated.t()) - sqr(momNewRotated.t()); double beta; // The rotated momentum should always lie along the +ve z-axis if ( momOldRotated.z() > ZERO ) beta = (-b + sqrt(sqr(b)-4.*a*c)) / 2. / a; else beta = (-b - sqrt(sqr(b)-4.*a*c)) / 2. / a; LorentzRotation boostOldToNew(0., 0., beta); // Total transform LorentzRotation transform = (newToZAxis.inverse())*boostOldToNew*oldToZAxis; // Assign the same spin info to the old and new particles newPart->spinInfo(oldPart->spinInfo()); newPart->spinInfo()->transform(oldMom, transform); } // If needed, insert default implementations of virtual function defined // in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs). void ConstituentReshuffler::persistentOutput(PersistentOStream &) const { } void ConstituentReshuffler::persistentInput(PersistentIStream &, int) { } ClassDescription ConstituentReshuffler::initConstituentReshuffler; // Definition of the static class description member. void ConstituentReshuffler::Init() { static ClassDocumentation documentation ("The ConstituentReshuffler class implements reshuffling " "of partons on their nominal mass shell to their constituent " "mass shells."); } diff --git a/Shower/Dipole/Utility/ConstituentReshuffler.h b/Shower/Dipole/Utility/ConstituentReshuffler.h --- a/Shower/Dipole/Utility/ConstituentReshuffler.h +++ b/Shower/Dipole/Utility/ConstituentReshuffler.h @@ -1,274 +1,259 @@ // -*- C++ -*- // // ConstituentReshuffler.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_ConstituentReshuffler_H #define HERWIG_ConstituentReshuffler_H // // This is the declaration of the ConstituentReshuffler class. // #include "ThePEG/Handlers/HandlerBase.h" #include "ThePEG/Utilities/Exception.h" #include "Herwig/Shower/PerturbativeProcess.h" #include "Herwig/Utilities/Reshuffler.h" namespace Herwig { using namespace ThePEG; /** * \ingroup DipoleShower * \author Simon Platzer, Stephen Webster * * \brief The ConstituentReshuffler class implements reshuffling * of partons on their nominal mass shell to their constituent * mass shells. * */ class ConstituentReshuffler: public HandlerBase, public Reshuffler { public: - /** @name Standard constructors and destructors. */ - //@{ - /** - * The default constructor. - */ - ConstituentReshuffler(); - - /** - * The destructor. - */ - virtual ~ConstituentReshuffler(); - //@} - -public: - /** * Reshuffle the outgoing partons to constituent * masses. Optionally, incoming partons are given * to absorb recoils. Add the non-reshuffled partons * to the intermediates list. Throw ConstituentReshufflerProblem * if a numerical problem prevents the solution of * the reshuffling equation. */ void reshuffle(PList& out, PPair& in, PList& intermediates, const bool decay, PList& decayPartons, PList& decayRecoilers); /** * Reshuffle the outgoing partons to constituent * masses. Optionally, incoming partons are given * to absorb recoils. Add the non-reshuffled partons * to the intermediates list. Throw ConstituentReshufflerProblem * if a numerical problem prevents the solution of * the reshuffling equation. */ void reshuffle(PList& out, PPair& in, PList& intermediates, const bool decay=false) { PList decayPartons; PList decayRecoilers; reshuffle(out, in, intermediates, decay, decayPartons, decayRecoilers); } /** * Reshuffle the outgoing partons following the showering * of the initial hard interaction to constituent masses, * for the case of outgoing decaying particles. * Throw ConstituentReshufflerProblem * if a numerical problem prevents the solution of * the reshuffling equation. */ void hardProcDecayReshuffle(PList& decaying, PList& eventOutgoing, PList& eventHard, PPair& eventIncoming, PList& eventIntermediates) ; /** * Reshuffle the outgoing partons following the showering * of a particle decay to constituent masses. * Throw ConstituentReshufflerProblem * if a numerical problem prevents the solution of * the reshuffling equation. */ void decayReshuffle(PerturbativeProcessPtr& decayProc, PList& eventOutgoing, PList& eventHard, PList& eventIntermediates) ; /** * Update the dipole event record and, if appropriate, * the relevant decay process. **/ void updateEvent( PList& intermediates, PList& eventIntermediates, PList& out, PList& eventOutgoing, PList& eventHard, PerturbativeProcessPtr decayProc = PerturbativeProcessPtr() ) ; /** * Update the spinInfo of a particle following reshuffling * to take account of the change in momentum. * Used only for unstable particles that need to be dealt with. **/ void updateSpinInfo( PPtr& oldPart, PPtr& newPart ) ; protected: /** * The function object defining the equation * to be solved in the case of separate recoilers * TODO - refine the whole implementation of separate partons and recoilers */ struct DecayReshuffleEquation { DecayReshuffleEquation (Energy q, PList::iterator m_begin, PList::iterator m_end, PList::iterator n_begin, PList::iterator n_end) : w(q), p_begin(m_begin), p_end(m_end), r_begin(n_begin), r_end(n_end) {} typedef double ArgType; typedef double ValType; static double aUnit(); static double vUnit(); double operator() (double xi) const; Energy w; PList::iterator p_begin; PList::iterator p_end; PList::iterator r_begin; PList::iterator r_end; }; public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} // If needed, insert declarations of virtual function defined in the // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs). private: /** * The static object used to initialize the description of this class. * Indicates that this is a concrete class with persistent data. */ static ClassDescription initConstituentReshuffler; /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ ConstituentReshuffler & operator=(const ConstituentReshuffler &) = delete; }; } #include "ThePEG/Utilities/ClassTraits.h" namespace ThePEG { /** @cond TRAITSPECIALIZATIONS */ /** This template specialization informs ThePEG about the * base classes of ConstituentReshuffler. */ template <> struct BaseClassTrait { /** Typedef of the first base class of ConstituentReshuffler. */ typedef HandlerBase NthBase; }; /** This template specialization informs ThePEG about the name of * the ConstituentReshuffler class and the shared object where it is defined. */ template <> struct ClassTraits : public ClassTraitsBase { /** Return a platform-independent class name */ static string className() { return "Herwig::ConstituentReshuffler"; } /** * The name of a file containing the dynamic library where the class * ConstituentReshuffler is implemented. It may also include several, space-separated, * libraries if the class ConstituentReshuffler depends on other classes (base classes * excepted). In this case the listed libraries will be dynamically * linked in the order they are specified. */ static string library() { return "HwDipoleShower.so"; } }; /** @endcond */ } #endif /* HERWIG_ConstituentReshuffler_H */ diff --git a/Shower/Dipole/Utility/DipoleMCCheck.cc b/Shower/Dipole/Utility/DipoleMCCheck.cc --- a/Shower/Dipole/Utility/DipoleMCCheck.cc +++ b/Shower/Dipole/Utility/DipoleMCCheck.cc @@ -1,283 +1,281 @@ // -*- C++ -*- // // DipoleMCCheck.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 DipoleMCCheck class. // #include "DipoleMCCheck.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" using namespace Herwig; DipoleMCCheck::DipoleMCCheck() : HandlerBase(), theHardPtBins(10), theEmitterXBins(5), theSpectatorXBins(5), thePtBins(100), theZBins(100) { } -DipoleMCCheck::~DipoleMCCheck() {} - IBPtr DipoleMCCheck::clone() const { return new_ptr(*this); } IBPtr DipoleMCCheck::fullclone() const { return new_ptr(*this); } vector DipoleMCCheck::makeLogBins(double xlow, double xup, unsigned int n) const { vector res; double c = log10(xup/xlow) / (n-1.); for ( unsigned int k = 0; k < n; ++k ) res.push_back(xlow*pow(10.0,k*c)); return res; } // If needed, insert default implementations of virtual function defined // in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs). void DipoleMCCheck::doinitrun() { HandlerBase::doinitrun(); vector ptbins; double w = 0.5/theHardPtBins; for ( unsigned int k = 1; k <= theHardPtBins; ++k ) ptbins.push_back(k*w); vector xebins; if ( theEmitterXBins > 1 ) xebins = makeLogBins(1e-7,1.0,theEmitterXBins); else xebins.push_back(1.0); vector xsbins; if ( theSpectatorXBins > 1 ) xsbins = makeLogBins(1e-7,1.0,theSpectatorXBins); else xsbins.push_back(1.0); for ( vector::const_iterator xeit = xebins.begin(); xeit != xebins.end(); ++xeit ) { map::ptr,Ptr::ptr> > > xebin; for ( vector::const_iterator xsit = xsbins.begin(); xsit != xsbins.end(); ++xsit ) { map::ptr,Ptr::ptr> > xsbin; for ( vector::const_iterator ptit = ptbins.begin(); ptit != ptbins.end(); ++ptit ) { pair::ptr,Ptr::ptr> ptbin (new_ptr(Histogram(0.0,0.5,thePtBins)), new_ptr(Histogram(0.0,1.0,theZBins))); xsbin[*ptit] = ptbin; } xebin[*xsit] = xsbin; } histoMap[*xeit] = xebin; } } void DipoleMCCheck::dofinish() { HandlerBase::dofinish(); map::ptr,Ptr::ptr> > > >::iterator xeit; map::ptr,Ptr::ptr> > >::iterator xsit; map::ptr,Ptr::ptr> >::iterator ptit; double xelow = 0.0; double xeup = 0.0; for ( xeit = histoMap.begin(); xeit != histoMap.end(); ++xeit ) { xeup = xelow + xeit->first; double xslow = 0.0; double xsup = 0.0; for ( xsit = xeit->second.begin(); xsit != xeit->second.end(); ++xsit ) { xsup = xslow + xsit->first; // open files here ostringstream ptFileName; ptFileName << name() << "_pt_" << xelow << "_" << xeup << "_" << xslow << "_" << xsup << ".dat"; ofstream ptFile(ptFileName.str().c_str()); ostringstream zFileName; zFileName << name() << "_z_" << xelow << "_" << xeup << "_" << xslow << "_" << xsup << ".dat"; ofstream zFile(zFileName.str().c_str()); double ptlow = 0.0; double ptup = 0.0; for ( ptit = xsit->second.begin(); ptit != xsit->second.end(); ++ptit ) { ptup = ptlow + ptit->first; // dump histos here ptFile << "#\n# " << ptlow << " < \\kappa < " << ptup << "\n#\n"; zFile << "#\n# " << ptlow << " < \\kappa < " << ptup << "\n#\n"; ptit->second.first->simpleOutput(ptFile,false); ptit->second.second->simpleOutput(zFile,false); ptFile << "\n\n\n"; zFile << "\n\n\n"; ptlow = ptup; } xslow = xsup; } xelow = xeup; } } void DipoleMCCheck::book(double xe,double xs, Energy dScale, Energy hardPt, Energy pt, double z, double weight) { map::ptr,Ptr::ptr> > > >::iterator xeit; map::ptr,Ptr::ptr> > >::iterator xsit; map::ptr,Ptr::ptr> >::iterator ptit; if ( theEmitterXBins == 1 || xe >= 1. ) xeit = --histoMap.end(); else xeit = histoMap.upper_bound(xe); if ( theSpectatorXBins == 1 || xs >= 1. ) xsit = --(xeit->second.end()); else xsit = xeit->second.upper_bound(xs); if ( theHardPtBins == 1 || hardPt/dScale >= 0.5 ) ptit = --(xsit->second.end()); else ptit = xsit->second.upper_bound(hardPt/dScale); ptit->second.first->addWeighted(pt/dScale,weight); ptit->second.second->addWeighted(z,weight); } void DipoleMCCheck::persistentOutput(PersistentOStream & os) const { os << theHardPtBins << theEmitterXBins << theSpectatorXBins << thePtBins << theZBins; } void DipoleMCCheck::persistentInput(PersistentIStream & is, int) { is >> theHardPtBins >> theEmitterXBins >> theSpectatorXBins >> thePtBins >> theZBins; } ClassDescription DipoleMCCheck::initDipoleMCCheck; // Definition of the static class description member. void DipoleMCCheck::Init() { static ClassDocumentation documentation ("DipoleMCCheck is used to perform checks for " "the dipole shower."); static Parameter interfaceHardPtBins ("HardPtBins", "HardPtBins", &DipoleMCCheck::theHardPtBins, 10, 1, 0, false, false, Interface::lowerlim); static Parameter interfaceEmitterXBins ("EmitterXBins", "EmitterXBins", &DipoleMCCheck::theEmitterXBins, 5, 1, 0, false, false, Interface::lowerlim); static Parameter interfaceSpectatorXBins ("SpectatorXBins", "SpectatorXBins", &DipoleMCCheck::theSpectatorXBins, 5, 1, 0, false, false, Interface::lowerlim); static Parameter interfacePtBins ("PtBins", "PtBins", &DipoleMCCheck::thePtBins, 100, 1, 0, false, false, Interface::lowerlim); static Parameter interfaceZBins ("ZBins", "ZBins", &DipoleMCCheck::theZBins, 100, 1, 0, false, false, Interface::lowerlim); } diff --git a/Shower/Dipole/Utility/DipoleMCCheck.h b/Shower/Dipole/Utility/DipoleMCCheck.h --- a/Shower/Dipole/Utility/DipoleMCCheck.h +++ b/Shower/Dipole/Utility/DipoleMCCheck.h @@ -1,231 +1,223 @@ // -*- C++ -*- // // DipoleMCCheck.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_DipoleMCCheck_H #define HERWIG_DipoleMCCheck_H // // This is the declaration of the DipoleMCCheck class. // #include "ThePEG/Handlers/HandlerBase.h" #include "Herwig/Utilities/Histogram.h" namespace Herwig { using namespace ThePEG; /** * \ingroup DipoleShower * \author Simon Platzer * * \brief DipoleMCCheck is used to perform checks for * the dipole shower. * * @see \ref DipoleMCCheckInterfaces "The interfaces" * defined for DipoleMCCheck. */ class DipoleMCCheck: public HandlerBase { public: - /** @name Standard constructors and destructors. */ - //@{ /** * The default constructor. */ DipoleMCCheck(); - /** - * The destructor. - */ - virtual ~DipoleMCCheck(); - //@} - public: /** * Book a point. */ void book(double xe,double xs, Energy dScale, Energy hardPt, Energy pt, double z, double weight); public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} protected: /** * 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(); // If needed, insert declarations of virtual function defined in the // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs). private: /** * The number of bins in the starting scale * divided by the dipole scale, the upper bound * is 1/2. */ unsigned int theHardPtBins; /** * The number of bins in the emitter fraction. * The lenght of the zero bin is taken to be * 10^(-7). */ unsigned int theEmitterXBins; /** * The number of bins in the spectator fraction. * The lenght of the zero bin is taken to be * 10^(-7). */ unsigned int theSpectatorXBins; /** * The number of bins in pt dicided by the * dipole scale; the upper bound is 1/2 */ unsigned int thePtBins; /** * The number of bins in z */ unsigned int theZBins; /** * The recursive map structure: xe, xs, hard pt / GeV * to histograms for pt and z; * output is done such that there's one file for each * xe,xs bin, including all histograms for the * hard pt bins. */ map::ptr,Ptr::ptr> > > > histoMap; /** * Helper to make logarithmic bins. */ vector makeLogBins(double xlow, double xup, unsigned int n) const; private: /** * The static object used to initialize the description of this class. * Indicates that this is a concrete class with persistent data. */ static ClassDescription initDipoleMCCheck; /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ DipoleMCCheck & operator=(const DipoleMCCheck &) = delete; }; } #include "ThePEG/Utilities/ClassTraits.h" namespace ThePEG { /** @cond TRAITSPECIALIZATIONS */ /** This template specialization informs ThePEG about the * base classes of DipoleMCCheck. */ template <> struct BaseClassTrait { /** Typedef of the first base class of DipoleMCCheck. */ typedef HandlerBase NthBase; }; /** This template specialization informs ThePEG about the name of * the DipoleMCCheck class and the shared object where it is defined. */ template <> struct ClassTraits : public ClassTraitsBase { /** Return a platform-independent class name */ static string className() { return "Herwig::DipoleMCCheck"; } /** * The name of a file containing the dynamic library where the class * DipoleMCCheck is implemented. It may also include several, space-separated, * libraries if the class DipoleMCCheck depends on other classes (base classes * excepted). In this case the listed libraries will be dynamically * linked in the order they are specified. */ static string library() { return "HwDipoleShower.so"; } }; /** @endcond */ } #endif /* HERWIG_DipoleMCCheck_H */ diff --git a/Shower/Dipole/Utility/IntrinsicPtGenerator.cc b/Shower/Dipole/Utility/IntrinsicPtGenerator.cc --- a/Shower/Dipole/Utility/IntrinsicPtGenerator.cc +++ b/Shower/Dipole/Utility/IntrinsicPtGenerator.cc @@ -1,198 +1,196 @@ // -*- C++ -*- // // IntrinsicPtGenerator.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 IntrinsicPtGenerator class. // #include "IntrinsicPtGenerator.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Config/Constants.h" #include "DipolePartonSplitter.h" #include "Herwig/Shower/ShowerHandler.h" #include "Herwig/PDF/HwRemDecayer.h" using namespace Herwig; IntrinsicPtGenerator::IntrinsicPtGenerator() : HandlerBase(), theValenceIntrinsicPtScale(1.0*GeV), theSeaIntrinsicPtScale(1.0*GeV) {} -IntrinsicPtGenerator::~IntrinsicPtGenerator() {} - IBPtr IntrinsicPtGenerator::clone() const { return new_ptr(*this); } IBPtr IntrinsicPtGenerator::fullclone() const { return new_ptr(*this); } LorentzRotation IntrinsicPtGenerator::kick(PPair& in, PList& intermediates) { if ( theValenceIntrinsicPtScale == 0.0*GeV && theSeaIntrinsicPtScale == 0.0*GeV ) return LorentzRotation(); assert(ShowerHandler::currentHandler()); tHwRemDecPtr remDec = ShowerHandler::currentHandler()->remnantDecayer(); assert(remDec); Lorentz5Momentum Q = in.first->momentum() + in.second->momentum(); // first parton if (in.first->coloured()) { Axis perp; Energy pt = 0.*GeV; double phi = 2.*Constants::pi*UseRandom::rnd(); Axis dir = in.first->momentum().vect().unit(); perp = dir.orthogonal(); perp.rotate(phi,dir); double r = sqrt(-log(1.-UseRandom::rnd())); if ( remDec->content().first.isValenceQuark(in.first) ) pt = sqrt(2.) * theValenceIntrinsicPtScale * r; else { assert(in.first->id() == ParticleID::g || remDec->content().first.isSeaQuark(in.first)); pt = sqrt(2.) * theSeaIntrinsicPtScale * r; } PPtr nin = new_ptr(Particle(in.first->dataPtr())); DipolePartonSplitter::change(in.first,nin,true); nin->set5Momentum(Lorentz5Momentum(0.*GeV,in.first->momentum().vect() + pt * perp)); intermediates.push_back(in.first); in.first = nin; } // second parton if (in.second->coloured()) { Axis perp; Energy pt = 0.*GeV; double phi = 2.*Constants::pi*UseRandom::rnd(); Axis dir = in.second->momentum().vect().unit(); perp = dir.orthogonal(); perp.rotate(phi,dir); double r = sqrt(-log(1.-UseRandom::rnd())); if ( remDec->content().second.isValenceQuark(in.second) ) pt = sqrt(2.) * theValenceIntrinsicPtScale * r; else { assert(in.second->id() == ParticleID::g || remDec->content().second.isSeaQuark(in.second)); pt = sqrt(2.) * theSeaIntrinsicPtScale * r; } PPtr nin = new_ptr(Particle(in.second->dataPtr())); nin->colourInfo(new_ptr(ColourBase())); DipolePartonSplitter::change(in.second,nin,true); nin->set5Momentum(Lorentz5Momentum(0.*GeV,in.second->momentum().vect() + pt * perp)); intermediates.push_back(in.second); in.second = nin; } // restore mometum conservation Lorentz5Momentum nQ = in.first->momentum() + in.second->momentum(); double x = Q.m()/nQ.m(); nQ *= x; // work around for constructor problems // in Lorentz5Vector constructor, mass is not guaranteed // to be zero, event it was set as such before, // since gets recalculated in the LorentzVector Lorentz5Momentum scaled = x * in.first->momentum(); scaled.setMass(ZERO); scaled.rescaleEnergy(); in.first->set5Momentum(scaled); scaled = x * in.second->momentum(); scaled.setMass(ZERO); scaled.rescaleEnergy(); in.second->set5Momentum(scaled); // apparently, this is more stable than boosting directly with // n_Q.boostVector()-Q.boostVector() Boost beta1 = -Q.boostVector(); Boost beta2 = nQ.boostVector(); LorentzRotation transform (beta1); transform.boost(beta2); return transform; } // If needed, insert default implementations of virtual function defined // in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs). void IntrinsicPtGenerator::persistentOutput(PersistentOStream & os) const { os << ounit(theValenceIntrinsicPtScale,GeV) << ounit(theSeaIntrinsicPtScale,GeV); } void IntrinsicPtGenerator::persistentInput(PersistentIStream & is, int) { is >> iunit(theValenceIntrinsicPtScale,GeV) >> iunit(theSeaIntrinsicPtScale,GeV); } ClassDescription IntrinsicPtGenerator::initIntrinsicPtGenerator; // Definition of the static class description member. void IntrinsicPtGenerator::Init() { static ClassDocumentation documentation ("IntrinsicPtGenerator generates intrinsic pt for massless " "incoming partons in a shower independent way."); static Parameter interfaceValenceIntrinsicPtScale ("ValenceIntrinsicPtScale", "The width of the intrinsic pt Gaussian distribution for valence partons.", &IntrinsicPtGenerator::theValenceIntrinsicPtScale, GeV, 1.0*GeV, 0.0*GeV, 0*GeV, false, false, Interface::lowerlim); static Parameter interfaceSeaIntrinsicPtScale ("SeaIntrinsicPtScale", "The width of the intrinsic pt Gaussian distribution for sea partons.", &IntrinsicPtGenerator::theSeaIntrinsicPtScale, GeV, 1.0*GeV, 0.0*GeV, 0*GeV, false, false, Interface::lowerlim); } diff --git a/Shower/Dipole/Utility/IntrinsicPtGenerator.h b/Shower/Dipole/Utility/IntrinsicPtGenerator.h --- a/Shower/Dipole/Utility/IntrinsicPtGenerator.h +++ b/Shower/Dipole/Utility/IntrinsicPtGenerator.h @@ -1,174 +1,166 @@ // -*- C++ -*- // // IntrinsicPtGenerator.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_IntrinsicPtGenerator_H #define HERWIG_IntrinsicPtGenerator_H // // This is the declaration of the IntrinsicPtGenerator class. // #include "ThePEG/Handlers/HandlerBase.h" #include "ThePEG/Vectors/LorentzRotation.h" namespace Herwig { using namespace ThePEG; /** * \ingroup DipoleShower * \author Simon Platzer * * \brief IntrinsicPtGenerator generates intrinsic pt for massless * incoming partons in a shower independent way. * * @see \ref IntrinsicPtGeneratorInterfaces "The interfaces" * defined for IntrinsicPtGenerator. */ class IntrinsicPtGenerator: public HandlerBase { public: - /** @name Standard constructors and destructors. */ - //@{ /** * The default constructor. */ IntrinsicPtGenerator(); - /** - * The destructor. - */ - virtual ~IntrinsicPtGenerator(); - //@} - public: /** * Generate intrinsic pt for the given incoming * partons and return the transformation to be * applied on the final state particles. Add the * old incoming partons to the given list. */ LorentzRotation kick(PPair& in, PList& intermediates); public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} // If needed, insert declarations of virtual function defined in the // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs). private: /** * The mean of the Gaussian distribution for * the intrinsic pt of valence partons. */ Energy theValenceIntrinsicPtScale; /** * The mean of the Gaussian distribution for * the intrinsic pt of sea partons. */ Energy theSeaIntrinsicPtScale; private: /** * The static object used to initialize the description of this class. * Indicates that this is a concrete class with persistent data. */ static ClassDescription initIntrinsicPtGenerator; /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ IntrinsicPtGenerator & operator=(const IntrinsicPtGenerator &) = delete; }; } #include "ThePEG/Utilities/ClassTraits.h" namespace ThePEG { /** @cond TRAITSPECIALIZATIONS */ /** This template specialization informs ThePEG about the * base classes of IntrinsicPtGenerator. */ template <> struct BaseClassTrait { /** Typedef of the first base class of IntrinsicPtGenerator. */ typedef HandlerBase NthBase; }; /** This template specialization informs ThePEG about the name of * the IntrinsicPtGenerator class and the shared object where it is defined. */ template <> struct ClassTraits : public ClassTraitsBase { /** Return a platform-independent class name */ static string className() { return "Herwig::IntrinsicPtGenerator"; } /** * The name of a file containing the dynamic library where the class * IntrinsicPtGenerator is implemented. It may also include several, space-separated, * libraries if the class IntrinsicPtGenerator depends on other classes (base classes * excepted). In this case the listed libraries will be dynamically * linked in the order they are specified. */ static string library() { return "HwDipoleShower.so"; } }; /** @endcond */ } #endif /* HERWIG_IntrinsicPtGenerator_H */ diff --git a/Shower/Dipole/Utility/PDFRatio.cc b/Shower/Dipole/Utility/PDFRatio.cc --- a/Shower/Dipole/Utility/PDFRatio.cc +++ b/Shower/Dipole/Utility/PDFRatio.cc @@ -1,144 +1,142 @@ // -*- C++ -*- // // PDFRatio.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 PDFRatio class. // #include "PDFRatio.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/Shower/ShowerHandler.h" #include "Herwig/PDF/HwRemDecayer.h" using namespace Herwig; PDFRatio::PDFRatio() : HandlerBase(), theValenceExtrapolation(0.7), theSeaExtrapolation(0.6), theFreezingScale(1.0*GeV) {} -PDFRatio::~PDFRatio() {} - IBPtr PDFRatio::clone() const { return new_ptr(*this); } IBPtr PDFRatio::fullclone() const { return new_ptr(*this); } double PDFRatio::operator() (const PDF& pdf, Energy2 scale, tcPDPtr from, tcPDPtr to, double x, double z) const { if ( x/z > 1.0 ) return 0.0; if(z==1)return 1.0; if ( scale < sqr(theFreezingScale) ) scale = sqr(theFreezingScale); tcHwRemDecPtr remDec = ShowerHandler::currentHandler()->remnantDecayer(); const HwRemDecayer::HadronContent* theContent = 0; if ( remDec->content().first.hadron == pdf.particle() ) { theContent = &(remDec->content().first); } if ( remDec->content().second.hadron == pdf.particle() ) { theContent = &(remDec->content().second); } assert(theContent); double exFrom = theContent->isValenceQuarkData(from) ? theValenceExtrapolation : theSeaExtrapolation; double exTo = theContent->isValenceQuarkData(to) ? theValenceExtrapolation : theSeaExtrapolation; double xTo = x/z; double fromPDF = x < exFrom ? pdf.xfx(from,scale,x) : ((1.-x)/(1.-exFrom)) * pdf.xfx(from,scale,exFrom); if ( fromPDF < 1e-8 ) fromPDF = 0.0; double toPDF = xTo < exTo ? pdf.xfx(to,scale,xTo) : ((1.-xTo)/(1.-exTo)) * pdf.xfx(to,scale,exTo); if ( toPDF < 1e-8 ) toPDF = 0.0; if ( toPDF == 0.0 || fromPDF == 0.0 ) return 0.0; return toPDF / fromPDF; } // If needed, insert default implementations of virtual function defined // in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs). void PDFRatio::persistentOutput(PersistentOStream & os) const { os << theValenceExtrapolation << theSeaExtrapolation << ounit(theFreezingScale,GeV); } void PDFRatio::persistentInput(PersistentIStream & is, int) { is >> theValenceExtrapolation >> theSeaExtrapolation >> iunit(theFreezingScale,GeV); } ClassDescription PDFRatio::initPDFRatio; // Definition of the static class description member. void PDFRatio::Init() { static ClassDocumentation documentation ("PDFRatio implements numerically stable pdf ratios."); static Parameter interfaceValenceExtrapolation ("ValenceExtrapolation", "The x from which on extrapolation should be done for valence partons.", &PDFRatio::theValenceExtrapolation, 0.7, 0.0, 1.0, false, false, Interface::limited); static Parameter interfaceSeaExtrapolation ("SeaExtrapolation", "The x from which on extrapolation should be done for valence partons.", &PDFRatio::theSeaExtrapolation, 0.6, 0.0, 1.0, false, false, Interface::limited); static Parameter interfaceFreezingScale ("FreezingScale", "The scale below which the PDFs are frozen.", &PDFRatio::theFreezingScale, GeV, 1.0*GeV, 0.0*GeV, 0*GeV, false, false, Interface::lowerlim); } diff --git a/Shower/Dipole/Utility/PDFRatio.h b/Shower/Dipole/Utility/PDFRatio.h --- a/Shower/Dipole/Utility/PDFRatio.h +++ b/Shower/Dipole/Utility/PDFRatio.h @@ -1,178 +1,170 @@ // -*- C++ -*- // // PDFRatio.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_PDFRatio_H #define HERWIG_PDFRatio_H // // This is the declaration of the PDFRatio class. // #include "ThePEG/Handlers/HandlerBase.h" #include "ThePEG/PDF/PDF.h" namespace Herwig { using namespace ThePEG; /** * \ingroup DipoleShower * \author Simon Platzer * * \brief PDFRatio implements numerically stable PDF ratios. * * @see \ref PDFRatioInterfaces "The interfaces" * defined for PDFRatio. */ class PDFRatio: public HandlerBase { public: - /** @name Standard constructors and destructors. */ - //@{ /** * The default constructor. */ PDFRatio(); - /** - * The destructor. - */ - virtual ~PDFRatio(); - //@} - public: /** * For the given PDF, scale and partons from and to and * x,z values return the ratio xf_to(x/z) / xf_from(x) */ double operator() (const PDF& pdf, Energy2 scale, tcPDPtr from, tcPDPtr to, double x, double z) const; public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} // If needed, insert declarations of virtual function defined in the // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs). private: /** * The x from which on extrapolation should * be done for valence partons. */ double theValenceExtrapolation; /** * The x from which on extrapolation should * be done for sea partons. */ double theSeaExtrapolation; /** * The scale below which the PDF will be frozen */ Energy theFreezingScale; private: /** * The static object used to initialize the description of this class. * Indicates that this is a concrete class with persistent data. */ static ClassDescription initPDFRatio; /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ PDFRatio & operator=(const PDFRatio &) = delete; }; } #include "ThePEG/Utilities/ClassTraits.h" namespace ThePEG { /** @cond TRAITSPECIALIZATIONS */ /** This template specialization informs ThePEG about the * base classes of PDFRatio. */ template <> struct BaseClassTrait { /** Typedef of the first base class of PDFRatio. */ typedef HandlerBase NthBase; }; /** This template specialization informs ThePEG about the name of * the PDFRatio class and the shared object where it is defined. */ template <> struct ClassTraits : public ClassTraitsBase { /** Return a platform-independent class name */ static string className() { return "Herwig::PDFRatio"; } /** * The name of a file containing the dynamic library where the class * PDFRatio is implemented. It may also include several, space-separated, * libraries if the class PDFRatio depends on other classes (base classes * excepted). In this case the listed libraries will be dynamically * linked in the order they are specified. */ static string library() { return "HwDipoleShower.so"; } }; /** @endcond */ } #endif /* HERWIG_PDFRatio_H */