diff --git a/Hadronization/HwppSelector.cc b/Hadronization/HwppSelector.cc --- a/Hadronization/HwppSelector.cc +++ b/Hadronization/HwppSelector.cc @@ -1,273 +1,273 @@ // -*- C++ -*- // // HwppSelector.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 HwppSelector class. // #include "HwppSelector.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/Utilities/Kinematics.h" #include "ThePEG/Utilities/Selector.h" #include "ThePEG/Repository/UseRandom.h" #include "CheckId.h" #include #include using namespace Herwig; DescribeClass describeHwppSelector("Herwig::HwppSelector","Herwig.so"); IBPtr HwppSelector::clone() const { return new_ptr(*this); } IBPtr HwppSelector::fullclone() const { return new_ptr(*this); } void HwppSelector::doinit() { HadronSelector::doinit(); } void HwppSelector::persistentOutput(PersistentOStream & os) const { os << _mode << _enhanceSProb << ounit(_m0Decay,GeV) << _massMeasure - << _charmDIQuarkWtFactor << _bottomDIQuarkWtFactor; + << _scHadronWtFactor << _sbHadronWtFactor; } void HwppSelector::persistentInput(PersistentIStream & is, int) { is >> _mode >> _enhanceSProb >> iunit(_m0Decay,GeV) >> _massMeasure - >> _charmDIQuarkWtFactor >> _bottomDIQuarkWtFactor; + >> _scHadronWtFactor >> _sbHadronWtFactor; } void HwppSelector::Init() { static ClassDocumentation documentation ("The HwppSelector class implements the Herwig algorithm for selecting" " the hadrons", "The hadronization used the selection algorithm described in \\cite{Kupco:1998fx}.", "%\\cite{Kupco:1998fx}\n" "\\bibitem{Kupco:1998fx}\n" " A.~Kupco,\n" " ``Cluster hadronization in HERWIG 5.9,''\n" " arXiv:hep-ph/9906412.\n" " %%CITATION = HEP-PH/9906412;%%\n" ); // put useMe() only in correct place! static Switch interfaceMode ("Mode", "Which algorithm to use", &HwppSelector::_mode, 1, false, false); static SwitchOption interfaceModeKupco (interfaceMode, "Kupco", "Use the Kupco approach", 0); static SwitchOption interfaceModeHwpp (interfaceMode, "Hwpp", "Use the Herwig approach", 1); static Switch interfaceEnhanceSProb ("EnhanceSProb", "Option for enhancing strangeness", &HwppSelector::_enhanceSProb, 0, false, false); static SwitchOption interfaceEnhanceSProbNo (interfaceEnhanceSProb, "No", "No strangeness enhancement.", 0); static SwitchOption interfaceEnhanceSProbScaled (interfaceEnhanceSProb, "Scaled", "Scaled strangeness enhancement", 1); static SwitchOption interfaceEnhanceSProbExponential (interfaceEnhanceSProb, "Exponential", "Exponential strangeness enhancement", 2); static Switch interfaceMassMeasure ("MassMeasure", "Option to use different mass measures", &HwppSelector::_massMeasure,0,false,false); static SwitchOption interfaceMassMeasureMass (interfaceMassMeasure, "Mass", "Mass Measure", 0); static SwitchOption interfaceMassMeasureLambda (interfaceMassMeasure, "Lambda", "Lambda Measure", 1); - static Parameter interfaceCharmDIQuarkWtFactor - ("charmDIQuarkWtFactor", - "Charmed diquark wight coefficients", - &HwppSelector::_charmDIQuarkWtFactor, 1., 0., 10., + static Parameter interfacescHadronWtFactor + ("scHadronWtFactor", + "Wight factor for strenge-charm heavy hadrns", + &HwppSelector::_scHadronWtFactor, 1., 0., 10., false, false, Interface::limited); - static Parameter interfaceBottomDIQuarkWtFactor - ("bottomDIQuarkWtFactor", - "Bottom diquark wight coefficients", - &HwppSelector::_bottomDIQuarkWtFactor, 1., 0., 10., + static Parameter interfacesbHadronWtFactor + ("sbHadronWtFactor", + "Wight factor for strenge-bottom heavy hadrns", + &HwppSelector::_sbHadronWtFactor, 1., 0., 10., false, false, Interface::limited); static Parameter interfaceDecayMassScale ("DecayMassScale", "Cluster decay mass scale", &HwppSelector::_m0Decay, GeV, 1.0*GeV, 0.1*GeV, 50.*GeV, false, false, Interface::limited); } pair HwppSelector::chooseHadronPair(const Energy cluMass,tcPDPtr par1, tcPDPtr par2,tcPDPtr ) const { // if either of the input partons is a diquark don't allow diquarks to be // produced bool diquark = !(DiquarkMatcher::Check(par1->id()) || DiquarkMatcher::Check(par2->id())); bool quark = true; // if the Herwig algorithm if(_mode ==1) { if(UseRandom::rnd() > 1./(1.+pwtDIquarkS0()+pwtDIquarkS1()) && cluMass > massLightestBaryonPair(par1,par2)) { diquark = true; quark = false; } else { useMe(); diquark = false; quark = true; } } // weights for the different possibilities Energy weight, wgtsum(ZERO); // loop over all hadron pairs with the allowed flavours static vector hadrons; hadrons.clear(); for(unsigned int ix=0;ixiColour())) == 3 && !DiquarkMatcher::Check(quarktopick->id())) continue; if(!diquark && abs(int(quarktopick->iColour())) == 3 && DiquarkMatcher::Check(quarktopick->id())) continue; HadronTable::const_iterator tit1 = table().find(make_pair(abs(par1->id()),quarktopick->id())); HadronTable::const_iterator tit2 = table().find(make_pair(quarktopick->id(),abs(par2->id()))); // If not in table skip if(tit1 == table().end()||tit2==table().end()) continue; // tables empty skip const KupcoData & T1 = tit1->second; const KupcoData & T2 = tit2->second; if(T1.empty()||T2.empty()) continue; // if too massive skip if(cluMass <= T1.begin()->mass + T2.begin()->mass) continue; // quark weight double quarkWeight = pwt(quarktopick->id()); if(abs(quarktopick->id()) == 3) { // Decoupling the weight of heavy strenge hadrons if(_enhanceSProb == 0 && abs(par1->id()) == 4) { - quarkWeight = pwt(quarktopick->id())*_charmDIQuarkWtFactor; + quarkWeight = pwt(quarktopick->id())*_scHadronWtFactor; } else if(_enhanceSProb == 0 && abs(par1->id()) == 5) { - quarkWeight = pwt(quarktopick->id())*_bottomDIQuarkWtFactor; + quarkWeight = pwt(quarktopick->id())*_sbHadronWtFactor; } // Scaling strangeness enhancement else if(_enhanceSProb == 1) { double scale = double(sqr(_m0Decay/cluMass)); quarkWeight = (_maxScale < scale) ? 0. : pow(quarkWeight,scale); } // Exponential strangeness enhancement else if(_enhanceSProb == 2) { Energy2 mass2; Energy endpointmass = par1->mass() + par2->mass(); // Choose to use either the cluster mass // or to use the lambda measure mass2 = (_massMeasure == 0) ? sqr(cluMass) : sqr(cluMass) - sqr(endpointmass); double scale = double(sqr(_m0Decay)/mass2); quarkWeight = (_maxScale < scale) ? 0. : exp(-scale); } } // loop over the hadrons KupcoData::const_iterator H1,H2; for(H1 = T1.begin();H1 != T1.end(); ++H1) { for(H2 = T2.begin();H2 != T2.end(); ++H2) { // break if cluster too light if(cluMass < H1->mass + H2->mass) break; weight = quarkWeight * H1->overallWeight * H2->overallWeight * Kinematics::pstarTwoBodyDecay(cluMass, H1->mass, H2->mass); int signQ = 0; assert (par1 && quarktopick); assert (par2); assert(quarktopick->CC()); if(CheckId::canBeHadron(par1, quarktopick->CC()) && CheckId::canBeHadron(quarktopick, par2)) signQ = +1; else if(CheckId::canBeHadron(par1, quarktopick) && CheckId::canBeHadron(quarktopick->CC(), par2)) signQ = -1; else { cerr << "Could not make sign for" << par1->id()<< " " << quarktopick->id() << " " << par2->id() << "\n"; assert(false); } if (signQ == -1) quarktopick = quarktopick->CC(); // construct the object with the info Kupco a(quarktopick, H1->ptrData, H2->ptrData, weight); hadrons.push_back(a); wgtsum += weight; } } } if (hadrons.empty()) return make_pair(tcPDPtr(),tcPDPtr()); // select the hadron wgtsum *= UseRandom::rnd(); unsigned int ix=0; do { wgtsum-= hadrons[ix].weight; ++ix; } while(wgtsum > ZERO && ix < hadrons.size()); if(ix == hadrons.size() && wgtsum > ZERO) return make_pair(tcPDPtr(),tcPDPtr()); --ix; assert(hadrons[ix].idQ); int signHad1 = signHadron(par1, hadrons[ix].idQ->CC(), hadrons[ix].hadron1); int signHad2 = signHadron(par2, hadrons[ix].idQ, hadrons[ix].hadron2); assert( signHad1 != 0 && signHad2 != 0 ); return make_pair ( signHad1 > 0 ? hadrons[ix].hadron1 : tcPDPtr(hadrons[ix].hadron1->CC()), signHad2 > 0 ? hadrons[ix].hadron2 : tcPDPtr(hadrons[ix].hadron2->CC())); } diff --git a/Hadronization/HwppSelector.h b/Hadronization/HwppSelector.h --- a/Hadronization/HwppSelector.h +++ b/Hadronization/HwppSelector.h @@ -1,176 +1,179 @@ // -*- C++ -*- // // HwppSelector.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_HwppSelector_H #define HERWIG_HwppSelector_H // // This is the declaration of the HwppSelector class. // #include "HadronSelector.h" #include "HwppSelector.fh" namespace Herwig { using namespace ThePEG; /** \ingroup hadronization * The HwppSelector class selects the hadrons produced in cluster decay using * the Herwig variant of the cluster model. * * @see \ref HwppSelectorInterfaces "The interfaces" * defined for HwppSelector. */ class HwppSelector: public HadronSelector { public: /** * The default constructor. */ HwppSelector() : HadronSelector(1), _mode(1), _enhanceSProb(0), _m0Decay(1.*GeV), - _charmDIQuarkWtFactor(1.), _bottomDIQuarkWtFactor(1.) + _scHadronWtFactor(1.), _sbHadronWtFactor(1.) {} /** * * This method is used to choose a pair of hadrons. * * Given the mass of a cluster and the particle pointers of its * two (or three) constituents, this returns the pair of particle pointers of * the two hadrons with proper flavour numbers. * Furthermore, the first of the two hadron must have the * constituent with par1, and the second must have the constituent with par2. * At the moment it does *nothing* in the case that also par3 is present. * * Kupco's method is used, rather than one used in FORTRAN HERWIG * The idea is to build on the fly a table of all possible pairs * of hadrons (Had1,Had2) (that we can call "cluster decay channels") * which are kinematically above threshold and have flavour * Had1=(par1,quarktopick->CC()), Had2=(quarktopick,par2), where quarktopick * is the poniter of: * --- d, u, s, c, b * if either par1 or par2 is a diquark; * --- d, u, s, c, b, dd, ud, uu, sd, su, ss, * cd, cu, cs, cc, bd, bu, bs, bc, bb * if both par1 and par2 are quarks. * The weight associated with each channel is given by the product * of: the phase space available including the spin factor 2*J+1, * the constant weight factor for chosen idQ, * the octet-singlet isoscalar mixing factor, and finally * the singlet-decuplet weight factor. */ pair chooseHadronPair(const Energy cluMass,tcPDPtr par1, tcPDPtr par2,tcPDPtr par3 = PDPtr()) const ; public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ HwppSelector & operator=(const HwppSelector &) = delete; private: /** * Which algorithm to use */ unsigned int _mode; /** * Flag that switches between no strangeness enhancement, scaling enhancement, * and exponential enhancement (in numerical order) */ int _enhanceSProb; /** * Parameter that governs the strangeness enhancement scaling */ Energy _m0Decay; /** * Flag that switches between mass measures used in strangeness enhancement: * cluster mass, or the lambda measure - ( m_{clu}^2 - (m_q + m_{qbar})^2 ) */ int _massMeasure; /** * Constant variable that stops the scale in strangeness enhancement from * becoming too large */ const double _maxScale = 20.; /** - * Heavy diquark wight coefficients + * Heavy strange-charm hadron wight coefficient */ - double _charmDIQuarkWtFactor; + double _scHadronWtFactor; - double _bottomDIQuarkWtFactor; + /** + * Heavy strange-bottom hadron wight coefficient + */ + double _sbHadronWtFactor; }; } #endif /* HERWIG_HwppSelector_H */ diff --git a/src/defaults/Hadronization.in b/src/defaults/Hadronization.in --- a/src/defaults/Hadronization.in +++ b/src/defaults/Hadronization.in @@ -1,101 +1,101 @@ # -*- ThePEG-repository -*- ############################################################ # Setup of default hadronization # # There are no user servicable parts inside. # # Anything that follows below should only be touched if you # know what you're doing. ############################################################# cd /Herwig/Particles create ThePEG::ParticleData Cluster setup Cluster 81 Cluster 0.00990 0.0 0.0 0.0 0 0 0 1 create ThePEG::ParticleData Remnant setup Remnant 82 Remnant 0.00990 0.0 0.0 0.0 0 0 0 1 mkdir /Herwig/Hadronization cd /Herwig/Hadronization create Herwig::ClusterHadronizationHandler ClusterHadHandler create Herwig::PartonSplitter PartonSplitter create Herwig::ClusterFinder ClusterFinder create Herwig::ColourReconnector ColourReconnector create Herwig::ClusterFissioner ClusterFissioner create Herwig::LightClusterDecayer LightClusterDecayer create Herwig::ClusterDecayer ClusterDecayer create Herwig::HwppSelector HadronSelector newdef ClusterHadHandler:PartonSplitter PartonSplitter newdef ClusterHadHandler:ClusterFinder ClusterFinder newdef ClusterHadHandler:ColourReconnector ColourReconnector newdef ClusterHadHandler:ClusterFissioner ClusterFissioner newdef ClusterHadHandler:LightClusterDecayer LightClusterDecayer newdef ClusterHadHandler:ClusterDecayer ClusterDecayer newdef ClusterHadHandler:MinVirtuality2 0.1*GeV2 newdef ClusterHadHandler:MaxDisplacement 1.0e-10*millimeter newdef ClusterHadHandler:UnderlyingEventHandler NULL newdef ClusterFissioner:HadronSelector HadronSelector newdef LightClusterDecayer:HadronSelector HadronSelector newdef ClusterDecayer:HadronSelector HadronSelector newdef ColourReconnector:ColourReconnection Yes newdef ColourReconnector:ReconnectionProbabilityBaryonic 0.7 newdef ColourReconnector:ReconnectionProbability 0.95 newdef ColourReconnector:Algorithm Baryonic newdef ColourReconnector:OctetTreatment All # Clustering parameters for light quarks newdef ClusterFissioner:ClMaxLight 3.649 newdef ClusterFissioner:ClPowLight 2.780 newdef ClusterFissioner:PSplitLight 0.899 newdef ClusterDecayer:ClDirLight 1 newdef ClusterDecayer:ClSmrLight 0.78 # Clustering parameters for b-quarks newdef ClusterFissioner:ClMaxBottom 3.757 newdef ClusterFissioner:ClPowBottom 0.547 newdef ClusterFissioner:PSplitBottom 0.625 newdef ClusterDecayer:ClDirBottom 1 newdef ClusterDecayer:ClSmrBottom 0.078 newdef HadronSelector:SingleHadronLimitBottom 0.000 # Clustering parameters for c-quarks newdef ClusterFissioner:ClMaxCharm 3.950 newdef ClusterFissioner:ClPowCharm 2.559 newdef ClusterFissioner:PSplitCharm 0.994 newdef ClusterDecayer:ClDirCharm 1 newdef ClusterDecayer:ClSmrCharm 0.163 newdef HadronSelector:SingleHadronLimitCharm 0.000 # Clustering parameters for exotic quarks # (e.g. hadronizing Susy particles) newdef ClusterFissioner:ClMaxExotic 2.7*GeV newdef ClusterFissioner:ClPowExotic 1.46 newdef ClusterFissioner:PSplitExotic 1.00 newdef ClusterDecayer:ClDirExotic 1 newdef ClusterDecayer:ClSmrExotic 0. newdef HadronSelector:SingleHadronLimitExotic 0. # newdef PartonSplitter:SplitPwtSquark 0.824135 newdef PartonSplitter:Split uds # newdef HadronSelector:PwtDquark 1.0 newdef HadronSelector:PwtUquark 1.0 newdef HadronSelector:PwtSquark 0.291717 newdef HadronSelector:PwtCquark 0.0 newdef HadronSelector:PwtBquark 0.0 newdef HadronSelector:PwtDIquarkS0 0.2 newdef HadronSelector:PwtDIquarkS1 0.1 newdef HadronSelector:SngWt 0.74 newdef HadronSelector:DecWt 0.62 newdef HadronSelector:Mode 1 newdef HadronSelector:BelowThreshold All -newdef HadronSelector:charmDIQuarkWtFactor 2.5 -newdef HadronSelector:bottomDIQuarkWtFactor 2.5 +newdef HadronSelector:scHadronWtFactor 2.5 +newdef HadronSelector:sbHadronWtFactor 2.5 create Herwig::SpinHadronizer SpinHadronizer