diff --git a/Hadronization/Cluster.h b/Hadronization/Cluster.h --- a/Hadronization/Cluster.h +++ b/Hadronization/Cluster.h @@ -1,346 +1,341 @@ // -*- C++ -*- // // Cluster.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_Cluster_H #define HERWIG_Cluster_H #include #include "Herwig/Utilities/EnumParticles.h" #include "CluHadConfig.h" #include "ClusterHadronizationHandler.fh" #include "Cluster.fh" #include "ThePEG/Utilities/ClassDescription.h" namespace Herwig { using namespace ThePEG; /** \ingroup Hadronization * \class Cluster * \brief This class describes a cluster object. * \author Philip Stephens * \author Alberto Ribon * * This class represents a cluster, which is a colour singlet made usually * of two components (quark-antiquark, quark-diquark, antiquark-antidiquark) * or rarely by three components (quark-quark-quark, antiquark-antiquark- * antiquark). A reference to the container with the pointers to its * Components is provided. * * The class provides access to the pointers which point to: * * - The cluster parent. In the case that the cluster it is a fission * product of a heavy cluster the parent is a cluster. If the cluster * is formed from the perturbative partons then the parents will be * the colour connected partons that formed the cluster. * - The children (usually two). In the case the cluster is a * heavy cluster that undergoes fission the children are clusters. * Occasionally the cluster has been "redefined" (re-interpreted). For * example in the case that three quark or anti-quark components * have been redefined as two components (quark+diquark, or antiquark+ * antidiquark). * - The (eventual) reshuffling partner, necessary for energy-momentum * conservation when light clusters are decayed into single hadron. Not * all clusters will have a reshuffling partner. * * Notice that in order to determine the cluster position from the positions * of the components, the Cluster class needs some parameters. * Because the Cluster class is neither interfaced nor persistent, * a static pointer to the ClusterHadronizationHandler class instance, * where the parameters are, is used. This static pointer is * set via the method setPointerClusterHadHandler(), during the * run initialization, doinitrun() of ClusterHadronizationHandler. * * @see ClusterHadronizationHandler */ class Cluster : public Particle { protected: /** @name Standard constructors and destructors. */ //@{ /** * Default constructor. Only used in PersistentIStream */ Cluster(); /** * The ClassTraits class must be a friend to be able to * use the private default constructor. */ friend struct ClassTraits; public: /** * Constructor with a particleData pointer */ Cluster(tcEventPDPtr); /** * This creates a cluster from 2 (or 3) partons. */ Cluster(tPPtr part1, tPPtr part2, tPPtr part3 = tPPtr()); /** * Also a constructor where a particle is given not a cluster. */ Cluster(const Particle &); //@} /** * Number of quark (diquark) constituents (normally two). */ unsigned int numComponents() const { return _numComp; } /** * Sum of the constituent masses of the components of the cluster. */ Energy sumConstituentMasses() const; /** * Returns the ith constituent. */ tPPtr particle(int i) const; /** * Returns the original constituent carrying colour */ tPPtr colParticle(bool anti = false) const; /** * Returns the original constituent carrying anticolour */ tPPtr antiColParticle() const; /** * Returns whether the ith constituent is from a perturbative process. */ bool isPerturbative(int) const; /** * Indicates whether the ith constituent is a beam remnant. */ bool isBeamRemnant(int) const; /** * Sets whether the ith constituent is a beam remnant. */ void setBeamRemnant(int,bool); /** * Returns the clusters id, not the same as the PDG id. */ int clusterId() const { return _id; } - /** - * Returns the pre-colour-reconnected cluster ancestors - */ - ClusterPair ancestors() const; - public: /** * Returns true when a constituent is a beam remnant. */ bool isBeamCluster() const; /** * Set the pointer to the reshuffling partner cluster. */ void flagAsReshuffled() { _hasReshuffled = true; } /** * Sets the component (if any) that points to "part" as a beam remnant. */ void isBeamCluster(tPPtr part); /** * Returns true if this cluster is to be handled by the hadronization. */ bool isAvailable() const { return _isAvailable; } /** * Sets the value of availability. */ void isAvailable(bool inputAvailable) { _isAvailable = inputAvailable; } /** * Return true if the cluster does not have cluster parent. */ bool isStatusInitial() const { return parents().empty(); } /** * Return true if the cluster does not have cluster children and * it is not already decayed (i.e. it does not have hadron children) * (to be used only after the fission of heavy clusters). */ bool isReadyToDecay() const { return children().empty(); } /** * Return true if the cluster has one and only one cluster children * and no hadron children: that means either that its three quarks or * anti-quarks components have been redefined as two components * (quark+diquark, or antiquark+antidiquark), or that the cluster * has been used as a partner for the momentum reshuffling necessary * to conserve energy-momentum when a light cluster is decayed into * a single hadron (notice that this latter light cluster has * isRedefined() false, because it has an hadron child). * In both cases, the unique cluster children is the new redefined * cluster. The two cases can be distinguish by the next method. */ bool isRedefined() const { return ( children().size() == 1 && children()[0]->id() == ParticleID::Cluster ); } /** * Return true when it has a reshuffling partner. * Notice that a cluster can have hasBeenReshuffled() true but * isRedefined() false: this is the case of a light cluster * that decays into a single hadron. */ bool hasBeenReshuffled() const { return _hasReshuffled; } /** * Return true if the cluster has hadron children. */ bool isStatusFinal() const; public: /** * Calculate the 4-position vector of the cluster * Method made static so can be used in other places * Displacement of the ith constituent given by momentum \f$p_i\f$ * vertex \f$x_i\f$ and mass \f$m_i\f$ is * \f[ D_i = -C \log(r) \frac{p_i}{\sqrt{(p_i^2 - m_i^2)^2 + v^4}} \f] * where \f$r\f$ is a random number [0,1], * \f$v\f$ is the minimum virtuality and \f$C\f$ is * a conversion factor from GeV to millimeters. We can then find the * difference in \f$s\f$ factors as * \f[ (s_1-s_2) = \frac{(\vec{p}_1 + \vec{p}_2) \cdot (\vec{x}_2 - * \vec{x}_1)}{(\vec{p}_1 + \vec{p}_2) \cdot \vec{D}_1}. * \f] * if \f$s_2>s_1\f$ then \f$s_1 = 1.0\f$ otherwise \f$s_2 = 1.0\f$. * These are then used to determine the value of the clusters vertex as * \f[ X = \frac{1}{2} ( x_1 +x_2 + s_1 D_1 + s_2 D_2). \f] */ static LorentzPoint calculateX(tPPtr q1, tPPtr q2); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual PPtr 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 PPtr fullclone() const; //@} public: /** @name Input and output functions. */ //@{ /** * Standard function for writing to a persistent stream. */ void persistentOutput(PersistentOStream &) const; /** * Standard function for reading from a persistent stream. */ void persistentInput(PersistentIStream &, int); private: /** * Private and non-existent assignment operator. */ Cluster & operator=(const Cluster &) = delete; /** * Calculate the 5-momentum vector of the cluster * The 5-momentum of the cluster is given by * \f[ P = \sum_i p_i \f] * and the mass of the cluster is \f$m^2 = P^2\f$ */ void calculateP(); /** * Determines whether constituent p is perturbative or not. */ bool initPerturbative(tPPtr p); /** * Describe an abstract base class with persistent data. */ static ClassDescription initCluster; bool _isAvailable; //!< Whether the cluster is hadronizing bool _hasReshuffled; //!< Whether the cluster has been reshuffled ParticleVector _component; //!< The constituent partons tParticleVector _original; //!< The original components vector _isBeamRemnant; //!< Whether a parton is a beam remnant vector _isPerturbative; //!< Whether a parton is perturbative unsigned int _numComp; //!< The number of constituents long _id; //!< The id of this cluster }; } // end namespace Herwig #include "ThePEG/Utilities/ClassTraits.h" namespace ThePEG { /** @cond TRAITSPECIALIZATIONS */ /** * The following template specialization informs ThePEG about the * base class of Cluster. */ template <> struct BaseClassTrait { /** Typedef of the base class of Cluster. */ typedef Particle NthBase; }; /** * The following template specialization informs ThePEG about the * name of this class and the shared object where it is defined. */ template <> struct ClassTraits: public ClassTraitsBase { /** Return the class name. */ static string className() { return "Herwig::Cluster"; } /** Create a Particle object. */ static TPtr create() { return TPtr::Create(Herwig::Cluster()); } }; /** @endcond */ } #endif // HERWIG_Cluster_H diff --git a/Hadronization/ClusterDecayer.cc b/Hadronization/ClusterDecayer.cc --- a/Hadronization/ClusterDecayer.cc +++ b/Hadronization/ClusterDecayer.cc @@ -1,473 +1,990 @@ // -*- C++ -*- // // ClusterDecayer.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 ClusterDecayer class. // #include "ClusterDecayer.h" #include #include #include #include #include #include #include #include #include "Herwig/Utilities/Kinematics.h" #include "CheckId.h" #include "Cluster.h" #include #include using namespace Herwig; +using CluVecIt = ClusterDecayer::CluVecIt; +using Constants::pi; +using Constants::twopi; DescribeClass describeClusterDecayer("Herwig::ClusterDecayer","Herwig.so"); ClusterDecayer::ClusterDecayer() : _clDirLight(1), _clDirBottom(1), _clDirCharm(1), _clDirExotic(1), _clSmrLight(0.0), _clSmrBottom(0.0), _clSmrCharm(0.0), _clSmrExotic(0.0), _onshell(false), - _masstry(20) + _masstry(20), + _swapConstituents(false), + _octetOption(0) {} IBPtr ClusterDecayer::clone() const { return new_ptr(*this); } IBPtr ClusterDecayer::fullclone() const { return new_ptr(*this); } void ClusterDecayer::persistentOutput(PersistentOStream & os) const { os << _hadronsSelector << _clDirLight << _clDirBottom << _clDirCharm << _clDirExotic << _clSmrLight << _clSmrBottom - << _clSmrCharm << _clSmrExotic << _onshell << _masstry; + << _clSmrCharm << _clSmrExotic << _onshell << _masstry + << _swapConstituents << _octetOption; } void ClusterDecayer::persistentInput(PersistentIStream & is, int) { is >> _hadronsSelector >> _clDirLight >> _clDirBottom >> _clDirCharm >> _clDirExotic >> _clSmrLight >> _clSmrBottom - >> _clSmrCharm >> _clSmrExotic >> _onshell >> _masstry; + >> _clSmrCharm >> _clSmrExotic >> _onshell >> _masstry + >> _swapConstituents >> _octetOption; } void ClusterDecayer::Init() { static ClassDocumentation documentation ("This class is responsible for the two-body decays of normal clusters"); static Reference interfaceHadronSelector("HadronSelector", "A reference to the HadronSelector object", &Herwig::ClusterDecayer::_hadronsSelector, false, false, true, false); //ClDir for light, Bottom, Charm and exotic (e.g Susy) quarks static Switch interfaceClDirLight ("ClDirLight", "Cluster direction for light quarks", &ClusterDecayer::_clDirLight, true, false, false); static SwitchOption interfaceClDirLightPreserve (interfaceClDirLight, "Preserve", "Preserve the direction of the quark from the perturbative stage" " as the direction of the meson produced from it", true); static SwitchOption interfaceClDirLightIsotropic (interfaceClDirLight, "Isotropic", "Assign the direction of the meson randomly", false); static Switch interfaceClDirBottom ("ClDirBottom", "Cluster direction for bottom quarks", &ClusterDecayer::_clDirBottom, true, false, false); static SwitchOption interfaceClDirBottomPreserve (interfaceClDirBottom, "Preserve", "Preserve the direction of the quark from the perturbative stage" " as the direction of the meson produced from it", true); static SwitchOption interfaceClDirBottomIsotropic (interfaceClDirBottom, "Isotropic", "Assign the direction of the meson randomly", false); static Switch interfaceClDirCharm ("ClDirCharm", "Cluster direction for charm quarks", &ClusterDecayer::_clDirCharm, true, false, false); static SwitchOption interfaceClDirCharmPreserve (interfaceClDirCharm, "Preserve", "Preserve the direction of the quark from the perturbative stage" " as the direction of the meson produced from it", true); static SwitchOption interfaceClDirCharmIsotropic (interfaceClDirCharm, "Isotropic", "Assign the direction of the meson randomly", false); static Switch interfaceClDirExotic ("ClDirExotic", "Cluster direction for exotic quarks", &ClusterDecayer::_clDirExotic, true, false, false); static SwitchOption interfaceClDirExoticPreserve (interfaceClDirExotic, "Preserve", "Preserve the direction of the quark from the perturbative stage" " as the direction of the meson produced from it", true); static SwitchOption interfaceClDirExoticIsotropic (interfaceClDirExotic, "Isotropic", "Assign the direction of the meson randomly", false); // ClSmr for ligth, Bottom, Charm and exotic (e.g. Susy) quarks static Parameter interfaceClSmrLight ("ClSmrLight", "cluster direction Gaussian smearing for light quark", &ClusterDecayer::_clSmrLight, 0, 0.0 , 0.0 , 2.0,false,false,false); static Parameter interfaceClSmrBottom ("ClSmrBottom", "cluster direction Gaussian smearing for b quark", &ClusterDecayer::_clSmrBottom, 0, 0.0 , 0.0 , 2.0,false,false,false); static Parameter interfaceClSmrCharm ("ClSmrCharm", "cluster direction Gaussian smearing for c quark", &ClusterDecayer::_clSmrCharm, 0, 0.0 , 0.0 , 2.0,false,false,false); static Parameter interfaceClSmrExotic ("ClSmrExotic", "cluster direction Gaussian smearing for exotic quark", &ClusterDecayer::_clSmrExotic, 0, 0.0 , 0.0 , 2.0,false,false,false); static Switch interfaceOnShell ("OnShell", "Whether or not the hadrons produced should by on shell or generated using the" " mass generator.", &ClusterDecayer::_onshell, false, false, false); static SwitchOption interfaceOnShellOnShell (interfaceOnShell, "Yes", "Produce the hadrons on shell", true); static SwitchOption interfaceOnShellOffShell (interfaceOnShell, "No", "Generate the masses using the mass generator.", false); static Parameter interfaceMassTry ("MassTry", "The number attempts to generate the masses of the hadrons produced" " in the cluster decay.", &ClusterDecayer::_masstry, 20, 1, 50, false, false, Interface::limited); + static Switch interfaceSwapConstituents + ("SwapConstituents", + "Choose whether to allow the clusters to swap constituents after " + "the pair production, but before finally creating hadrons", + &ClusterDecayer::_swapConstituents, false, false, false); + static SwitchOption interfaceSwapConstituentsOff + (interfaceSwapConstituents, + "No", + "Ordinary cluster decay of Herwig, i.e. no constituent swapping.", + false); + static SwitchOption interfaceSwapConstituentsOn + (interfaceSwapConstituents, + "Yes", + "Allow clusters to swap constituents.", + true); + + static Switch interfaceOctetTreatment + ("OctetTreatment", + "Which octets are not allowed to be reconnected", + &ClusterDecayer::_octetOption, 0, false, false); + static SwitchOption interfaceOctetTreatmentFinal + (interfaceOctetTreatment, + "Final", + "Only prevent for the final (usuaslly non-perturbative) g -> q qbar splitting", + 0); + static SwitchOption interfaceOctetTreatmentAll + (interfaceOctetTreatment, + "All", + "Prevent for all octets", + 1); + } +// Functions taken from ColourReconnector + +namespace { + inline bool hasDiquark(ClusterVector::const_iterator cit) { + for(int i = 0; i<(*cit)->numComponents(); i++) { + if (DiquarkMatcher::Check(*((*cit)->particle(i)->dataPtr()))) + return true; + } + return false; + } +} + +namespace { + +double calculateRapidityRF(const Lorentz5Momentum & q1, + const Lorentz5Momentum & p2) { + //calculate rapidity wrt the direction of q1 + //angle between the particles in the RF of cluster of q1 + + // calculate the z component of p2 w.r.t the direction of q1 + if(q1.rho2()==ZERO) return 0.; + const Energy pz = p2.vect() * q1.vect().unit(); + if ( pz == ZERO ) return 0.; + + // Transverse momentum of p2 w.r.t the direction of q1 + const Energy pt = sqrt(p2.vect().mag2() - sqr(pz)); + + // Transverse mass pf p2 w.r.t to the direction of q1 + const Energy mtrans = sqrt(p2.mass()*p2.mass() + (pt*pt)); + + // Correct formula + const double y2 = log((p2.t() + abs(pz))/mtrans); + + return ( pz < ZERO ) ? -y2 : y2; +} + +} + +bool ClusterDecayer::_isColour8(tcPPtr p, tcPPtr q) const { + bool octet = false; + + // make sure we have a triplet and an anti-triplet + if ( ( p->hasColour() && q->hasAntiColour() ) || + ( p->hasAntiColour() && q->hasColour() ) ) { + + // true if p and q are originated from a colour octet + if ( !p->parents().empty() && !q->parents().empty() ) { + octet = ( p->parents()[0] == q->parents()[0] ) && + ( p->parents()[0]->data().iColour() == PDT::Colour8 ); + } + + // (Final) option: check if same colour8 parent + // or already found an octet. + if(_octetOption==0||octet) return octet; + + // (All) option handling more octets + // by browsing particle history/colour lines. + tColinePtr cline,aline; + + // Get colourlines form final states. + if(p->hasColour() && q->hasAntiColour()) { + cline = p-> colourLine(); + aline = q->antiColourLine(); + } + else { + cline = q-> colourLine(); + aline = p->antiColourLine(); + } + + // Follow the colourline of p. + if ( !p->parents().empty() ) { + tPPtr parent = p->parents()[0]; + while (parent) { + if(parent->data().iColour() == PDT::Colour8) { + // Coulour8 particles should have a colour + // and an anticolour line. Currently the + // remnant has none of those. Since the children + // of the remnant are not allowed to emit currently, + // the colour octet remnant is handled by the return + // statement above. The assert also catches other + // colour octets without clines. If the children of + // a remnant should be allowed to emit, the remnant + // should get appropriate colour lines and + // colour states. + // See Ticket: #407 + // assert(parent->colourLine()&&parent->antiColourLine()); + octet = (parent-> colourLine()==cline && + parent->antiColourLine()==aline); + } + if(octet||parent->parents().empty()) break; + parent = parent->parents()[0]; + } + } + } + + return octet; +} + +// TODO: make the 'normal' decay mechanism into a function to reduce the +// amount of code void ClusterDecayer::decay(const ClusterVector & clusters, tPVector & finalhadrons) { // Loop over all clusters, and if they are not too heavy (that is // intermediate clusters that have undergone to fission) or not // too light (that is final clusters that have been already decayed // into single hadron) then decay them into two hadrons. - - for (ClusterVector::const_iterator it = clusters.begin(); - it != clusters.end(); ++it) { - if ((*it)->isAvailable() && !(*it)->isStatusFinal() - && (*it)->isReadyToDecay()) { - pair prod = decayIntoTwoHadrons(*it); - if(!prod.first) - throw Exception() << "Can't perform decay of cluster " << **it - << "in ClusterDecayer::decay()" - << Exception::eventerror; - (*it)->addChild(prod.first); - (*it)->addChild(prod.second); - finalhadrons.push_back(prod.first); - finalhadrons.push_back(prod.second); + if (!_swapConstituents){ + for (ClusterVector::const_iterator it = clusters.begin(); + it != clusters.end(); ++it) { + if ((*it)->isAvailable() && !(*it)->isStatusFinal() + && (*it)->isReadyToDecay()) { + pair prod = decayIntoTwoHadrons(*it); + if(!prod.first) + throw Exception() << "Can't perform decay of cluster " << **it + << "in ClusterDecayer::decay()" + << Exception::eventerror; + (*it)->addChild(prod.first); + (*it)->addChild(prod.second); + finalhadrons.push_back(prod.first); + finalhadrons.push_back(prod.second); + } } } + else { + // Copy of list + ClusterVector cluCopy = clusters; + // List of clusters that have already been considered + ClusterVector deleted; + deleted.reserve(clusters.size()); + // Loop over the clusters but do not yet make the hadrons + for (ClusterVector::iterator it = cluCopy.begin(); + it != cluCopy.end(); ++it) { + // Check if the cluster would be ready to undergo decay + if ((*it)->isAvailable() && (*it)->isStatusFinal() + && (*it)->isReadyToDecay()) { + // As a first model, ignore clusters that already contain diquarks + if (hasDiquark(it)){ + pair prod = decayIntoTwoHadrons(*it); + if(!prod.first) + throw Exception() << "Can't perform decay of cluster " << **it + << "in ClusterDecayer::decay()" + << Exception::eventerror; + (*it)->addChild(prod.first); + (*it)->addChild(prod.second); + finalhadrons.push_back(prod.first); + finalhadrons.push_back(prod.second); + } + // Skip any clusters that we've already considered + if (find(deleted.begin(), deleted.end(), *it) != deleted.end()) + continue; + // Find the 'nearest partner' in delta eta and delta phi space + // This algorithm is the same as the mesonic part of the + // baryonic colour reconnection algorithm + CluVecIt candidate = _findPartnerCluster(it, cluCopy, deleted); + + // If there is no suitable candidate, decay normally + if (candidate == it){ + pair prod = decayIntoTwoHadrons(*it); + if(!prod.first) + throw Exception() << "Can't perform decay of cluster " << **it + << "in ClusterDecayer::decay()" + << Exception::eventerror; + (*it)->addChild(prod.first); + (*it)->addChild(prod.second); + finalhadrons.push_back(prod.first); + finalhadrons.push_back(prod.second); + } + // Else, we allow the two to 'decay together' + ParticleVector prods = twoClusterDecay(*it, *candidate); + + // Push back the decayed products + + } + + } // End of loop over cluster + + } +} + +ParticleVector ClusterDecayer::twoClusterDecay(tClusterPtr clA, + tClusterPtr clB) { + ParticleVector hadronOutput = ParticleVector(); + + // This is similar to the decayIntoTwoHadrons function, though + // instead of forming the two (or in this case four) hadrons, + // we will allow the constituents to be swapped. + // We need to require (at least at the moment, maybe in the future we + // could change it) that the clusters have exactly two components. + // If this is not the case, then send a warning because it is not suppose + // to happen, and then return. + if ( clA->numComponents() != 2 || clB->numComponents() != 2) { + generator()->logWarning( Exception("ClusterDecayer::twoClusterDecay " + "***Still clusters with not exactly 2 components*** ", + Exception::warning) ); + return ParticleVector(); + } + // Extract the id and particle pointer of the components of the clusters. + tPPtr ptrA1 = clA->particle(0); + tPPtr ptrA2 = clA->particle(1); + tPPtr ptrB1 = clB->particle(0); + tPPtr ptrB2 = clB->particle(1); + tcPDPtr ptrA1data = ptrA1->dataPtr(); + tcPDPtr ptrA2data = ptrA2->dataPtr(); + tcPDPtr ptrB1data = ptrB1->dataPtr(); + tcPDPtr ptrB2data = ptrB2->dataPtr(); + + bool isHadA1FlavSpecial = false; + bool cluDirHadA1 = _clDirLight; + double cluSmearHadA1 = _clSmrLight; + bool isHadA2FlavSpecial = false; + bool cluDirHadA2 = _clDirLight; + double cluSmearHadA2 = _clSmrLight; + + bool isHadB1FlavSpecial = false; + bool cluDirHadB1 = _clDirLight; + double cluSmearHadB1 = _clSmrLight; + bool isHadB2FlavSpecial = false; + bool cluDirHadB2 = _clDirLight; + double cluSmearHadB2 = _clSmrLight; + + + if (CheckId::isExotic(ptrA1data)) { + isHadA1FlavSpecial = true; + cluDirHadA1 = _clDirExotic; + cluSmearHadA1 = _clSmrExotic; + } + else if (CheckId::hasBottom(ptrA1data)) { + isHadA1FlavSpecial = true; + cluDirHadA1 = _clDirBottom; + cluSmearHadA1 = _clSmrBottom; + } + else if (CheckId::hasCharm(ptrA1data)) { + isHadA1FlavSpecial = true; + cluDirHadA1 = _clDirCharm; + cluSmearHadA1 = _clSmrCharm; + } + + if (CheckId::isExotic(ptrA2data)) { + isHadA2FlavSpecial = true; + cluDirHadA2 = _clDirExotic; + cluSmearHadA2 = _clSmrExotic; + } + else if (CheckId::hasBottom(ptrA2data)) { + isHadA2FlavSpecial = true; + cluDirHadA2 = _clDirBottom; + cluSmearHadA2 = _clSmrBottom; + } + else if (CheckId::hasCharm(ptrA2data)) { + isHadA2FlavSpecial = true; + cluDirHadA2 = _clDirCharm; + cluSmearHadA2 = _clSmrCharm; + } + + if (CheckId::isExotic(ptrB1data)) { + isHadB1FlavSpecial = true; + cluDirHadB1 = _clDirExotic; + cluSmearHadB1 = _clSmrExotic; + } + else if (CheckId::hasBottom(ptrB1data)) { + isHadB1FlavSpecial = true; + cluDirHadB1 = _clDirBottom; + cluSmearHadB1 = _clSmrBottom; + } + else if (CheckId::hasCharm(ptrB1data)) { + isHadB1FlavSpecial = true; + cluDirHadB1 = _clDirCharm; + cluSmearHadB1 = _clSmrCharm; + } + + if (CheckId::isExotic(ptrB2data)) { + isHadB2FlavSpecial = true; + cluDirHadB2 = _clDirExotic; + cluSmearHadB2 = _clSmrExotic; + } + else if (CheckId::hasBottom(ptrB2data)) { + isHadB2FlavSpecial = true; + cluDirHadB2 = _clDirBottom; + cluSmearHadB2 = _clSmrBottom; + } + else if (CheckId::hasCharm(ptrB2data)) { + isHadB2FlavSpecial = true; + cluDirHadB2 = _clDirCharm; + cluSmearHadB2 = _clSmrCharm; + } + + + bool isOriginA1Perturbative = clA->isPerturbative(0); + bool isOriginA2Perturbative = clA->isPerturbative(1); + bool isOriginB1Perturbative = clB->isPerturbative(0); + bool isOriginB2Perturbative = clB->isPerturbative(1); + + // We have to decide which, if any, of the two hadrons will have + // the momentum, in the cluster parent frame, smeared around the + // direction of its constituent (for Had1 is the one pointed by + // ptr1, and for Had2 is the one pointed by ptr2). + // This happens only if the flag _ClDirX is 1 and the constituent is + // perturbative (that is not coming from nonperturbative gluon splitting + // or cluster fission). In the case that both the hadrons satisfy this + // two requirements (of course only one must be treated, because the other + // one will have the momentum automatically fixed by the momentum + // conservation) then more priority is given in the case of a b-hadron. + // Finally, in the case that the two hadrons have same priority, then + // we choose randomly, with equal probability, one of the two. + + + int priorityHadA1 = 0; + if ( cluDirHadA1 && isOriginA1Perturbative ) { + priorityHadA1 = isHadA1FlavSpecial ? 2 : 1; + } + int priorityHadA2 = 0; + if ( cluDirHadA2 && isOriginA2Perturbative ) { + priorityHadA2 = isHadA2FlavSpecial ? 2 : 1; + } + if ( priorityHadA2 && priorityHadA1 == priorityHadA2 && UseRandom::rndbool() ) { + priorityHadA2 = 0; + } + + int priorityHadB1 = 0; + if ( cluDirHadB1 && isOriginB1Perturbative ) { + priorityHadB1 = isHadB1FlavSpecial ? 2 : 1; + } + int priorityHadB2 = 0; + if ( cluDirHadB2 && isOriginB2Perturbative ) { + priorityHadB2 = isHadB2FlavSpecial ? 2 : 1; + } + if ( priorityHadB2 && priorityHadB1 == priorityHadB2 && UseRandom::rndbool() ) { + priorityHadB2 = 0; + } + + + Lorentz5Momentum pClu = clA->momentum(); + bool secondHad = false; + Axis uSmear_v3; + if ( priorityHadA1 || priorityHadA2 ) { + + double cluSmear; + Lorentz5Momentum pQ; + if ( priorityHadA1 > priorityHadA2 ) { + pQ = ptrA1->momentum(); + cluSmear = cluSmearHadA1; + } else { + pQ = ptrA2->momentum(); + cluSmear = cluSmearHadA2; + secondHad = true; + } + + // To find the momenta of the two hadrons in the parent cluster frame + // we proceed as follows. First, we determine the unit vector parallel + // to the direction of the constituent in the cluster frame. Then we + // have to smear such direction using the following prescription: + // --- in theta angle w.r.t. such direction (not the z-axis), + // the drawing of the cosine of such angle is done via: + // 1.0 + cluSmear*log( rnd() ) + // (repeated until it gives a value above -1.0) + // --- in phi angle w.r.t. such direction, the drawing is simply flat. + // Then, given the direction in the parent cluster frame of one of the + // two hadrons, it is straighforward to get the momenta of both hadrons + // (always in the same parent cluster frame). + + pQ.boost( -pClu.boostVector() ); // boost from Lab to Cluster frame + uSmear_v3 = pQ.vect().unit(); + // skip if cluSmear is too small + if ( cluSmear > 0.001 ) { + // generate the smearing angle + double cosSmear; + do cosSmear = 1.0 + cluSmear*log( UseRandom::rnd() ); + while ( cosSmear < -1.0 ); + double sinSmear = sqrt( 1.0 - sqr(cosSmear) ); + // calculate rotation to z axis + LorentzRotation rot; + double sinth(sqrt(1.-sqr(uSmear_v3.z()))); + if(abs(uSmear_v3.perp2()/uSmear_v3.z())>1e-10) + rot.setRotate(-acos(uSmear_v3.z()), + Axis(-uSmear_v3.y()/sinth,uSmear_v3.x()/sinth,0.)); + // + random azimuthal rotation + rot.rotateZ(UseRandom::rnd()*twopi); + // set direction in rotated frame + Lorentz5Vector ltemp(0.,sinSmear,cosSmear,0.); + // rotate back + rot.invert(); + ltemp *= rot; + uSmear_v3 = ltemp.vect(); + } + } + else { + + // Isotropic decay: flat in cosTheta and phi. + uSmear_v3 = Axis(1.0, 0.0, 0.0); // just to set the rho to 1 + uSmear_v3.setTheta( acos( UseRandom::rnd( -1.0 , 1.0 ) ) ); + uSmear_v3.setPhi( UseRandom::rnd( -pi , pi ) ); + + } + + // Choose the produced quark-antiquark pair for the two clusters + pair dataPairA + = _hadronsSelector->chooseHadronPair(clA->mass(), + ptrA1data, ptrA2data); + if(dataPairA.first == tcPDPtr() || + dataPairA.second == tcPDPtr()) return ParticleVector(); + // Get the produced quark species + long qA = _hadronsSelector->getQ(); + pair dataPairB + = _hadronsSelector->chooseHadronPair(clB->mass(), + ptrB1data, ptrB2data); + if(dataPairB.first == tcPDPtr() || + dataPairB.second == tcPDPtr()) return ParticleVector(); + // Get the produced quark specie/s + long qB = _hadronsSelector->getQ(); + + return hadronOutput; + +} + + + +CluVecIt ClusterDecayer::_findPartnerCluster( + CluVecIt cl, ClusterVector & cv, + const ClusterVector& deleted) const { + + using Constants::pi; + using Constants::twopi; + + // Returns a candidate for possible reconnection + CluVecIt candidate = cl; + + + double maxrap = 0.0; + double minrap = 0.0; + double maxrapNormal = 0.0; + double minrapNormal = 0.0; + double maxsumnormal = 0.0; + + + // boost into RF of cl + Lorentz5Momentum cl1 = (*cl)->momentum(); + const Boost boostv(-cl1.boostVector()); + cl1.boost(boostv); + // boost constituents of cl into RF of cl + Lorentz5Momentum p1col = (*cl)->colParticle()->momentum(); + Lorentz5Momentum p1anticol = (*cl)->antiColParticle()->momentum(); + p1col.boost(boostv); + p1anticol.boost(boostv); + + + for (CluVecIt cit=cv.begin(); cit != cv.end(); ++cit) { + //avoid looping over clusters containing diquarks + if ( hasDiquark(cit) ) continue; + // Clusters do not have 3 components, but if they do, this line + // will need to be fixed + //if ( (*cit)->numComponents()==3 ) continue; + // Skip the cluster itself + if ( cit==cl ) continue; + + //skip the cluster to be deleted later 3->2 cluster + if ( find(deleted.begin(), deleted.end(), *cit) != deleted.end() ) + continue; + + if ( (*cl)->isBeamCluster() && (*cit)->isBeamCluster() ) + continue; + + // stop it putting far apart clusters together + // TODO: what does this line really do? + //if ( ( (**cl).vertex()-(**cit).vertex() ).m() >_maxDistance ) + // continue; + + const bool Colour8 = + _isColour8( (*cl)->colParticle(), (*cit)->antiColParticle() ) + || + _isColour8( (*cit)->colParticle(), (*cl)->antiColParticle() ) ; + if ( Colour8 ) continue; + + + // boost constituents of cit into RF of cl + Lorentz5Momentum p2col = (*cit)->colParticle()->momentum(); + Lorentz5Momentum p2anticol = (*cit)->antiColParticle()->momentum(); + + p2col.boost(boostv); + p2anticol.boost(boostv); + + // calculate the rapidity of the other constituents of the clusters + // w.r.t axis of p1anticol.vect.unit + const double rapq = calculateRapidityRF(p1anticol,p2col); + const double rapqbar = calculateRapidityRF(p1anticol,p2anticol); + + // configuration for normal CR + if ( rapq > 0.0 && rapqbar < 0.0 + && rapq > maxrap + && rapqbar < minrap ) { + maxrap = rapq; + minrap = rapqbar; + //sum of rapidities of quarks + const double normalsum = abs(rapq) + abs(rapqbar); + if ( normalsum > maxsumnormal ) { + maxsumnormal = normalsum; + maxrapNormal = rapq; + minrapNormal = rapqbar; + candidate = cit; + } + } + + } + + return candidate; } pair ClusterDecayer::decayIntoTwoHadrons(tClusterPtr ptr) { using Constants::pi; using Constants::twopi; // To decay the cluster into two hadrons one must distinguish between // constituent quarks (or diquarks) that originate from perturbative // processes (hard process or parton shower) from those that are generated // by the non-perturbative gluon splitting or from fission of heavy clusters. // In the latter case the two body decay is assumed to be *isotropic*. // In the former case instead, if proper flag are activated, the two body // decay is assumed to "remember" the direction of the constituents inside // the cluster, in the cluster frame. The actual smearing of the hadron // directions around the direction of the constituents, in the cluster // frame, can be different between non-b hadrons and b-hadrons, but is given // by the same functional form: // cosThetaSmearing = 1 + smearFactor * log( rnd() ) // (repeated until cosThetaSmearing > -1) // where the factor smearFactor is different between b- and non-b hadrons. // // We need to require (at least at the moment, maybe in the future we // could change it) that the cluster has exactly two components. // If this is not the case, then send a warning because it is not suppose // to happen, and then return. if ( ptr->numComponents() != 2 ) { generator()->logWarning( Exception("ClusterDecayer::decayIntoTwoHadrons " "***Still cluster with not exactly 2 components*** ", Exception::warning) ); return pair(); } // Extract the id and particle pointer of the two components of the cluster. tPPtr ptr1 = ptr->particle(0); tPPtr ptr2 = ptr->particle(1); tcPDPtr ptr1data = ptr1->dataPtr(); tcPDPtr ptr2data = ptr2->dataPtr(); bool isHad1FlavSpecial = false; bool cluDirHad1 = _clDirLight; double cluSmearHad1 = _clSmrLight; bool isHad2FlavSpecial = false; bool cluDirHad2 = _clDirLight; double cluSmearHad2 = _clSmrLight; if (CheckId::isExotic(ptr1data)) { isHad1FlavSpecial = true; cluDirHad1 = _clDirExotic; cluSmearHad1 = _clSmrExotic; } else if (CheckId::hasBottom(ptr1data)) { isHad1FlavSpecial = true; cluDirHad1 = _clDirBottom; cluSmearHad1 = _clSmrBottom; } else if (CheckId::hasCharm(ptr1data)) { isHad1FlavSpecial = true; cluDirHad1 = _clDirCharm; cluSmearHad1 = _clSmrCharm; } if (CheckId::isExotic(ptr2data)) { isHad2FlavSpecial = true; cluDirHad2 = _clDirExotic; cluSmearHad2 = _clSmrExotic; } else if (CheckId::hasBottom(ptr2data)) { isHad2FlavSpecial = true; cluDirHad2 = _clDirBottom; cluSmearHad2 = _clSmrBottom; } else if (CheckId::hasCharm(ptr2data)) { isHad2FlavSpecial = true; cluDirHad2 = _clDirCharm; cluSmearHad2 = _clSmrCharm; } bool isOrigin1Perturbative = ptr->isPerturbative(0); bool isOrigin2Perturbative = ptr->isPerturbative(1); // We have to decide which, if any, of the two hadrons will have // the momentum, in the cluster parent frame, smeared around the // direction of its constituent (for Had1 is the one pointed by // ptr1, and for Had2 is the one pointed by ptr2). // This happens only if the flag _ClDirX is 1 and the constituent is // perturbative (that is not coming from nonperturbative gluon splitting // or cluster fission). In the case that both the hadrons satisfy this // two requirements (of course only one must be treated, because the other // one will have the momentum automatically fixed by the momentum // conservation) then more priority is given in the case of a b-hadron. // Finally, in the case that the two hadrons have same priority, then // we choose randomly, with equal probability, one of the two. int priorityHad1 = 0; if ( cluDirHad1 && isOrigin1Perturbative ) { priorityHad1 = isHad1FlavSpecial ? 2 : 1; } int priorityHad2 = 0; if ( cluDirHad2 && isOrigin2Perturbative ) { priorityHad2 = isHad2FlavSpecial ? 2 : 1; } if ( priorityHad2 && priorityHad1 == priorityHad2 && UseRandom::rndbool() ) { priorityHad2 = 0; } Lorentz5Momentum pClu = ptr->momentum(); bool secondHad = false; Axis uSmear_v3; if ( priorityHad1 || priorityHad2 ) { double cluSmear; Lorentz5Momentum pQ; if ( priorityHad1 > priorityHad2 ) { pQ = ptr1->momentum(); cluSmear = cluSmearHad1; } else { pQ = ptr2->momentum(); cluSmear = cluSmearHad2; secondHad = true; } // To find the momenta of the two hadrons in the parent cluster frame // we proceed as follows. First, we determine the unit vector parallel // to the direction of the constituent in the cluster frame. Then we // have to smear such direction using the following prescription: // --- in theta angle w.r.t. such direction (not the z-axis), // the drawing of the cosine of such angle is done via: // 1.0 + cluSmear*log( rnd() ) // (repeated until it gives a value above -1.0) // --- in phi angle w.r.t. such direction, the drawing is simply flat. // Then, given the direction in the parent cluster frame of one of the // two hadrons, it is straighforward to get the momenta of both hadrons // (always in the same parent cluster frame). pQ.boost( -pClu.boostVector() ); // boost from Lab to Cluster frame uSmear_v3 = pQ.vect().unit(); // skip if cluSmear is too small if ( cluSmear > 0.001 ) { // generate the smearing angle double cosSmear; do cosSmear = 1.0 + cluSmear*log( UseRandom::rnd() ); while ( cosSmear < -1.0 ); double sinSmear = sqrt( 1.0 - sqr(cosSmear) ); // calculate rotation to z axis LorentzRotation rot; double sinth(sqrt(1.-sqr(uSmear_v3.z()))); if(abs(uSmear_v3.perp2()/uSmear_v3.z())>1e-10) rot.setRotate(-acos(uSmear_v3.z()), Axis(-uSmear_v3.y()/sinth,uSmear_v3.x()/sinth,0.)); // + random azimuthal rotation rot.rotateZ(UseRandom::rnd()*twopi); // set direction in rotated frame Lorentz5Vector ltemp(0.,sinSmear,cosSmear,0.); // rotate back rot.invert(); ltemp *= rot; uSmear_v3 = ltemp.vect(); } } else { // Isotropic decay: flat in cosTheta and phi. uSmear_v3 = Axis(1.0, 0.0, 0.0); // just to set the rho to 1 uSmear_v3.setTheta( acos( UseRandom::rnd( -1.0 , 1.0 ) ) ); uSmear_v3.setPhi( UseRandom::rnd( -pi , pi ) ); } - pair dataPair = _hadronsSelector->chooseHadronPair(ptr->mass(), ptr1data, ptr2data); if(dataPair.first == tcPDPtr() || dataPair.second == tcPDPtr()) return pair(); - // Create the two hadron particle objects with the specified id. PPtr ptrHad1,ptrHad2; // produce the hadrons on mass shell if(_onshell) { ptrHad1 = dataPair.first ->produceParticle(dataPair.first ->mass()); ptrHad2 = dataPair.second->produceParticle(dataPair.second->mass()); } // produce the hadrons with mass given by the mass generator else { unsigned int ntry(0); do { ptrHad1 = dataPair.first ->produceParticle(); ptrHad2 = dataPair.second->produceParticle(); ++ntry; } while(ntry<_masstry&&ptrHad1->mass()+ptrHad2->mass()>ptr->mass()); // if fails produce on shell and issue warning (should never happen??) if( ptrHad1->mass() + ptrHad2->mass() > ptr->mass() ) { generator()->log() << "Failed to produce off-shell hadrons in " << "ClusterDecayer::decayIntoTwoHadrons producing hadrons " << "on-shell" << endl; ptrHad1 = dataPair.first ->produceParticle(dataPair.first ->mass()); ptrHad2 = dataPair.second->produceParticle(dataPair.second->mass()); } } if (!ptrHad1 || !ptrHad2) { ostringstream s; s << "ClusterDecayer::decayIntoTwoHadrons ***Cannot create the two hadrons***" << dataPair.first ->PDGName() << " and " << dataPair.second->PDGName(); cerr << s.str() << endl; generator()->logWarning( Exception(s.str(), Exception::warning) ); } else { Lorentz5Momentum pHad1, pHad2; // 5-momentum vectors that we have to determine if ( secondHad ) uSmear_v3 *= -1.0; if (pClu.m() < ptrHad1->mass()+ptrHad2->mass() ) { throw Exception() << "Impossible Kinematics in ClusterDecayer::decayIntoTwoHadrons()" << Exception::eventerror; } Kinematics::twoBodyDecay(pClu,ptrHad1->mass(),ptrHad2->mass(),uSmear_v3, pHad1,pHad2); ptrHad1->set5Momentum(pHad1); ptrHad2->set5Momentum(pHad2); // Determine the positions of the two children clusters. LorentzPoint positionHad1 = LorentzPoint(); LorentzPoint positionHad2 = LorentzPoint(); calculatePositions(pClu, ptr->vertex(), pHad1, pHad2, positionHad1, positionHad2); ptrHad1->setVertex(positionHad1); ptrHad2->setVertex(positionHad2); } return pair(ptrHad1,ptrHad2); } void ClusterDecayer:: calculatePositions(const Lorentz5Momentum &pClu, const LorentzPoint &positionClu, const Lorentz5Momentum &, const Lorentz5Momentum &, LorentzPoint &positionHad1, LorentzPoint &positionHad2 ) const { // First, determine the relative positions of the children hadrons // with respect to their parent cluster, in the cluster reference frame, // assuming gaussian smearing with width inversely proportional to the // parent cluster mass. Length smearingWidth = hbarc / pClu.m(); LorentzDistance distanceHad[2]; for ( int i = 0; i < 2; i++ ) { // i==0 is the first hadron; i==1 is the second one Length delta[4]={ZERO,ZERO,ZERO,ZERO}; // smearing of the four components of the LorentzDistance, two at the same time to improve speed for ( int j = 0; j < 3; j += 2 ) { delta[j] = UseRandom::rndGauss(smearingWidth, Length(ZERO)); delta[j+1] = UseRandom::rndGauss(smearingWidth, Length(ZERO)); } // set the distance delta[0] = abs(delta[0]) +sqrt(sqr(delta[1])+sqr(delta[2])+sqr(delta[3])); distanceHad[i] = LorentzDistance(delta[1],delta[2],delta[3],delta[0]); // Boost such relative positions of the children hadrons, // with respect to their parent cluster, // from the cluster reference frame to the Lab frame. distanceHad[i].boost(pClu.boostVector()); } // Finally, determine the absolute positions of the children hadrons // in the Lab frame. positionHad1 = distanceHad[0] + positionClu; positionHad2 = distanceHad[1] + positionClu; } diff --git a/Hadronization/ClusterDecayer.h b/Hadronization/ClusterDecayer.h --- a/Hadronization/ClusterDecayer.h +++ b/Hadronization/ClusterDecayer.h @@ -1,167 +1,195 @@ // -*- C++ -*- // // ClusterDecayer.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_ClusterDecayer_H #define HERWIG_ClusterDecayer_H #include #include #include "CluHadConfig.h" #include "HadronSelector.h" #include "ClusterDecayer.fh" namespace Herwig { using namespace ThePEG; /** \ingroup Hadronization * \class ClusterDecayer * \brief This class decays the "normal" clusters * \author Philip Stephens * \author Alberto Ribon * * This class decays the "normal" clusters, e.g. ones that are not heavy * enough for fission, and not too light to decay into one hadron. * * This class is directs the production of hadrons via 2-body cluster decays. * The selection of the hadron flavours is given by Herwig::HadronSelector. * * @see HadronSelector * @see \ref ClusterDecayerInterfaces "The interfaces" * defined for ClusterDecayer. */ class ClusterDecayer: public Interfaced { public: /** @name Standard constructors and destructors. */ //@{ /** * Default constructor. */ ClusterDecayer(); //@} /** Decays all remaining clusters into hadrons. * This routine decays the clusters that are left after * Herwig::ClusterFissioner::fission and * Herwig::LightClusterDecayer::decay have been called. These are all * the "normal" clusters which are not forced into hadrons by * the other functions. */ void decay(const ClusterVector & clusters, tPVector & finalhadrons) ; + using CluVecIt = ClusterVector::iterator; + public: /** * Standard ThePEG function for writing a persistent stream. */ void persistentOutput(PersistentOStream &) const; /** * Standard ThePEG function for reading from a persistent stream. */ void persistentInput(PersistentIStream &, int); /** * Standard Init function used to initialize the interfaces. */ 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; //@} private: /** * Private and non-existent assignment operator. */ ClusterDecayer & operator=(const ClusterDecayer &) = delete; public: /** Decays the cluster into two hadrons. * * This routine is used to take a given cluster and decay it into * two hadrons which are returned. If one of the constituents is from * the perturbative regime then the direction of the perturbative parton * is remembered and the decay is preferentially in that direction. The * direction of the decay is given by * \f[ \cos \theta = 1 + S \log r_1 \f] * where \f$ S \f$ is a parameter of the model and \f$ r_1 \f$ is a random * number [0,1]. */ pair decayIntoTwoHadrons(tClusterPtr ptr); + CluVecIt _findPartnerCluster(CluVecIt cl, ClusterVector & cv, + const ClusterVector& a) const; + + + + /** + * @return true, if the two partons are splitting products of the same + * gluon + */ + bool _isColour8(tcPPtr p, tcPPtr q) const; + + ParticleVector twoClusterDecay(tClusterPtr clA, + tClusterPtr clB); + private: /** Compute the positions of the new hadrons based on the clusters position. * * This method calculates the positions of the children hadrons by a * call to ThePEG::RandomGenerator::rndGaussTwoNumbers with width inversely * proportional to the cluster mass, around the parent cluster position. */ void calculatePositions( const Lorentz5Momentum &, const LorentzPoint &, const Lorentz5Momentum &, const Lorentz5Momentum &, LorentzPoint &, LorentzPoint &) const; /** * Pointer to a Herwig::HadronSelector for choosing decay types */ Ptr::pointer _hadronsSelector; //@{ /** * Whether a cluster decays along the perturbative parton direction. */ bool _clDirLight; bool _clDirBottom; bool _clDirCharm; bool _clDirExotic; /** * The S parameter from decayIntoTwoHadrons */ double _clSmrLight; double _clSmrBottom; double _clSmrCharm; double _clSmrExotic; //@} /** * Whether or not the hadrons produced should be on-shell * or generated used the MassGenerator */ bool _onshell; /** * Number of tries to generate the masses of the decay products */ unsigned int _masstry; + /** + * Flag to choose whether or not to allow clusters to swap + * constituents. + */ + bool _swapConstituents; + + + /** + * Option for handling octets + */ + unsigned int _octetOption; + }; } #endif /* HERWIG_ClusterDecayer_H */ diff --git a/Hadronization/HadronSelector.cc b/Hadronization/HadronSelector.cc --- a/Hadronization/HadronSelector.cc +++ b/Hadronization/HadronSelector.cc @@ -1,998 +1,999 @@ // -*- C++ -*- // // HadronSelector.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 HadronSelector class. // #include "HadronSelector.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/ParVector.h" #include "ThePEG/Interface/RefVector.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include #include #include #include #include "CheckId.h" #include using namespace Herwig; DescribeAbstractClass describeHadronSelector("Herwig::HadronSelector","Herwig.so"); namespace { // debug helper void dumpTable(const HadronSelector::HadronTable & tbl) { typedef HadronSelector::HadronTable::const_iterator TableIter; for (TableIter it = tbl.begin(); it != tbl.end(); ++it) { - cerr << it->first.first << ' ' + cerr << it->first.first << ' ' << it->first.second << '\n'; for (HadronSelector::KupcoData::const_iterator jt = it->second.begin(); jt != it->second.end(); ++jt) { cerr << '\t' << *jt << '\n'; } } } bool weightIsLess (pair a, pair b) { return a.second < b.second; } } -HadronSelector::HadronSelector(unsigned int opt) +HadronSelector::HadronSelector(unsigned int opt) : _pwtDquark( 1.0 ),_pwtUquark( 1.0 ),_pwtSquark( 1.0 ),_pwtCquark( 0.0 ), _pwtBquark( 0.0 ),_pwtDIquark( 1.0 ), _weight1S0(Nmax,1.),_weight3S1(Nmax,1.),_weight1P1(Nmax,1.),_weight3P0(Nmax,1.), _weight3P1(Nmax,1.),_weight3P2(Nmax,1.),_weight1D2(Nmax,1.),_weight3D1(Nmax,1.), _weight3D2(Nmax,1.),_weight3D3(Nmax,1.), _repwt(Lmax,vector >(Jmax,vector(Nmax))), _sngWt( 1.0 ),_decWt( 1.0 ), - _topt(opt),_trial(0), - _limBottom(), _limCharm(), _limExotic(), belowThreshold_(0) + _topt(opt),_trial(0), + _limBottom(), _limCharm(), _limExotic(), belowThreshold_(0) { // The mixing angles // the ideal mixing angle const double idealAngleMix = atan( sqrt(0.5) ) * 180.0 / Constants::pi; // \eta-\eta' mixing angle _etamix = -23.0; // phi-omega mixing angle _phimix = +36.0; // h_1'-h_1 mixing angle _h1mix = idealAngleMix; // f_0(1710)-f_0(1370) mixing angle _f0mix = idealAngleMix; // f_1(1420)-f_1(1285)\f$ mixing angle _f1mix = idealAngleMix; // f'_2-f_2\f$ mixing angle _f2mix = +26.0; // eta_2(1870)-eta_2(1645) mixing angle _eta2mix = idealAngleMix; // phi(???)-omega(1650) mixing angle _omhmix = idealAngleMix; // phi_3-omega_3 mixing angle _ph3mix = +28.0; // eta(1475)-eta(1295) mixing angle _eta2Smix = idealAngleMix; // phi(1680)-omega(1420) mixing angle _phi2Smix = idealAngleMix; } void HadronSelector::persistentOutput(PersistentOStream & os) const { - os << _partons << _pwtDquark << _pwtUquark << _pwtSquark + os << _partons << _pwtDquark << _pwtUquark << _pwtSquark << _pwtCquark << _pwtBquark << _pwtDIquark - << _etamix << _phimix << _h1mix << _f0mix << _f1mix << _f2mix - << _eta2mix << _omhmix << _ph3mix << _eta2Smix << _phi2Smix - << _weight1S0 << _weight3S1 << _weight1P1 << _weight3P0 << _weight3P1 + << _etamix << _phimix << _h1mix << _f0mix << _f1mix << _f2mix + << _eta2mix << _omhmix << _ph3mix << _eta2Smix << _phi2Smix + << _weight1S0 << _weight3S1 << _weight1P1 << _weight3P0 << _weight3P1 << _weight3P2 << _weight1D2 << _weight3D1 << _weight3D2 << _weight3D3 << _forbidden << _sngWt << _decWt << _repwt << _pwt << _limBottom << _limCharm << _limExotic << belowThreshold_ << _table; } void HadronSelector::persistentInput(PersistentIStream & is, int) { - is >> _partons >> _pwtDquark >> _pwtUquark >> _pwtSquark - >> _pwtCquark >> _pwtBquark - >> _pwtDIquark>> _etamix >> _phimix >> _h1mix >> _f0mix >> _f1mix >> _f2mix - >> _eta2mix >> _omhmix >> _ph3mix >> _eta2Smix >> _phi2Smix - >> _weight1S0 >> _weight3S1 >> _weight1P1 >> _weight3P0 >> _weight3P1 + is >> _partons >> _pwtDquark >> _pwtUquark >> _pwtSquark + >> _pwtCquark >> _pwtBquark + >> _pwtDIquark>> _etamix >> _phimix >> _h1mix >> _f0mix >> _f1mix >> _f2mix + >> _eta2mix >> _omhmix >> _ph3mix >> _eta2Smix >> _phi2Smix + >> _weight1S0 >> _weight3S1 >> _weight1P1 >> _weight3P0 >> _weight3P1 >> _weight3P2 >> _weight1D2 >> _weight3D1 >> _weight3D2 >> _weight3D3 >> _forbidden >> _sngWt >> _decWt >> _repwt >> _pwt >> _limBottom >> _limCharm >> _limExotic >> belowThreshold_ >> _table; } void HadronSelector::Init() { static ClassDocumentation documentation ("There is no documentation for the HadronSelector class"); static Parameter interfacePwtDquark("PwtDquark","Weight for choosing a quark D", &HadronSelector::_pwtDquark, 0, 1.0, 0.0, 10.0, false,false,false); static Parameter interfacePwtUquark("PwtUquark","Weight for choosing a quark U", &HadronSelector::_pwtUquark, 0, 1.0, 0.0, 10.0, false,false,false); static Parameter interfacePwtSquark("PwtSquark","Weight for choosing a quark S", &HadronSelector::_pwtSquark, 0, 1.0, 0.0, 10.0, false,false,false); static Parameter interfacePwtCquark("PwtCquark","Weight for choosing a quark C", &HadronSelector::_pwtCquark, 0, 0.0, 0.0, 10.0, false,false,false); static Parameter interfacePwtBquark("PwtBquark","Weight for choosing a quark B", &HadronSelector::_pwtBquark, 0, 0.0, 0.0, 10.0, false,false,false); static Parameter interfacePwtDIquark("PwtDIquark","Weight for choosing a DIquark", &HadronSelector::_pwtDIquark, 0, 1.0, 0.0, 100.0, false,false,false); static Parameter interfaceSngWt("SngWt","Weight for singlet baryons", &HadronSelector::_sngWt, 0, 1.0, 0.0, 10.0, false,false,false); static Parameter interfaceDecWt("DecWt","Weight for decuplet baryons", &HadronSelector::_decWt, 0, 1.0, 0.0, 10.0, false,false,false); static RefVector interfacePartons ("Partons", "The partons which are to be considered as the consistuents of the hadrons.", &HadronSelector::_partons, -1, false, false, true, false, false); static RefVector interfaceForbidden ("Forbidden", "The PDG codes of the particles which cannot be produced in the hadronization.", &HadronSelector::_forbidden, -1, false, false, true, false, false); // // mixing angles // // the ideal mixing angle const double idealAngleMix = atan( sqrt(0.5) ) * 180.0 / Constants::pi; static Parameter interface11S0Mixing ("11S0Mixing", "The mixing angle for the I=0 mesons from the 1 1S0 multiplet," " i.e. eta and etaprime.", &HadronSelector::_etamix, -23., -180., 180., false, false, Interface::limited); static Parameter interface13S1Mixing ("13S1Mixing", "The mixing angle for the I=0 mesons from the 1 3S1 multiplet," " i.e. phi and omega.", &HadronSelector::_phimix, +36., -180., 180., false, false, Interface::limited); static Parameter interface11P1Mixing ("11P1Mixing", "The mixing angle for the I=0 mesons from the 1 1P1 multiplet," " i.e. h_1' and h_1.", &HadronSelector::_h1mix, idealAngleMix, -180., 180., false, false, Interface::limited); static Parameter interface13P0Mixing ("13P0Mixing", "The mixing angle for the I=0 mesons from the 1 3P0 multiplet," " i.e. f_0(1710) and f_0(1370).", &HadronSelector::_f0mix, idealAngleMix, -180., 180., false, false, Interface::limited); static Parameter interface13P1Mixing ("13P1Mixing", "The mixing angle for the I=0 mesons from the 1 3P1 multiplet," " i.e. f_1(1420) and f_1(1285).", &HadronSelector::_f1mix, idealAngleMix, -180., 180., false, false, Interface::limited); static Parameter interface13P2Mixing ("13P2Mixing", "The mixing angle for the I=0 mesons from the 1 3P2 multiplet," " i.e. f'_2 and f_2.", &HadronSelector::_f2mix, 26.0, -180., 180., false, false, Interface::limited); static Parameter interface11D2Mixing ("11D2Mixing", "The mixing angle for the I=0 mesons from the 1 1D2 multiplet," " i.e. eta_2(1870) and eta_2(1645).", &HadronSelector::_eta2mix, idealAngleMix, -180., 180., false, false, Interface::limited); static Parameter interface13D0Mixing ("13D0Mixing", "The mixing angle for the I=0 mesons from the 1 3D0 multiplet," " i.e. eta_2(1870) phi(?) and omega(1650).", &HadronSelector::_omhmix, idealAngleMix, -180., 180., false, false, Interface::limited); static Parameter interface13D1Mixing ("13D1Mixing", "The mixing angle for the I=0 mesons from the 1 3D1 multiplet," " i.e. phi_3 and omega_3.", &HadronSelector::_ph3mix, 28.0, -180., 180., false, false, Interface::limited); static Parameter interface21S0Mixing ("21S0Mixing", "The mixing angle for the I=0 mesons from the 2 1S0 multiplet," " i.e. eta(1475) and eta(1295).", &HadronSelector::_eta2Smix, idealAngleMix, -180., 180., false, false, Interface::limited); static Parameter interface23S1Mixing ("23S1Mixing", "The mixing angle for the I=0 mesons from the 1 3S1 multiplet," " i.e. phi(1680) and omega(1420).", &HadronSelector::_phi2Smix, idealAngleMix, -180., 180., false, false, Interface::limited); // // the meson weights // static ParVector interface1S0Weights ("1S0Weights", "The weights for the 1S0 multiplets start with n=1.", &HadronSelector::_weight1S0, Nmax, 1.0, 0.0, 100.0, false, false, Interface::limited); static ParVector interface3S1Weights ("3S1Weights", "The weights for the 3S1 multiplets start with n=1.", &HadronSelector::_weight3S1, Nmax, 1.0, 0.0, 100.0, false, false, Interface::limited); static ParVector interface1P1Weights ("1P1Weights", "The weights for the 1P1 multiplets start with n=1.", &HadronSelector::_weight1P1, Nmax, 1.0, 0.0, 100.0, false, false, Interface::limited); static ParVector interface3P0Weights ("3P0Weights", "The weights for the 3P0 multiplets start with n=1.", &HadronSelector::_weight3P0, Nmax, 1.0, 0.0, 100.0, false, false, Interface::limited); static ParVector interface3P1Weights ("3P1Weights", "The weights for the 3P1 multiplets start with n=1.", &HadronSelector::_weight3P1, Nmax, 1.0, 0.0, 100.0, false, false, Interface::limited); static ParVector interface3P2Weights ("3P2Weights", "The weights for the 3P2 multiplets start with n=1.", &HadronSelector::_weight3P2, Nmax, 1.0, 0.0, 100.0, false, false, Interface::limited); static ParVector interface1D2Weights ("1D2Weights", "The weights for the 1D2 multiplets start with n=1.", &HadronSelector::_weight1D2, Nmax, 1.0, 0.0, 100.0, false, false, Interface::limited); static ParVector interface3D1Weights ("3D1Weights", "The weights for the 3D1 multiplets start with n=1.", &HadronSelector::_weight3D1, Nmax, 1.0, 0.0, 100.0, false, false, Interface::limited); static ParVector interface3D2Weights ("3D2Weights", "The weights for the 3D2 multiplets start with n=1.", &HadronSelector::_weight3D2, Nmax, 1.0, 0.0, 100.0, false, false, Interface::limited); static ParVector interface3D3Weights ("3D3Weights", "The weights for the 3D3 multiplets start with n=1.", &HadronSelector::_weight3D3, Nmax, 1.0, 0.0, 100.0, false, false, Interface::limited); static Switch interfaceTrial ("Trial", "A Debugging option to only produce certain types of hadrons", &HadronSelector::_trial, 0, false, false); static SwitchOption interfaceTrialAll (interfaceTrial, "All", "Produce all the hadrons", 0); static SwitchOption interfaceTrialPions (interfaceTrial, "Pions", "Only produce pions", 1); static SwitchOption interfaceTrialSpin2 (interfaceTrial, "Spin2", "Only mesons with spin less than or equal to two are produced", 2); static SwitchOption interfaceTrialSpin3 (interfaceTrial, "Spin3", "Only hadrons with spin less than or equal to three are produced", 3); static Parameter interfaceSingleHadronLimitBottom ("SingleHadronLimitBottom", "Threshold for one-hadron decay of b-cluster", &HadronSelector::_limBottom, 0, 0.0, 0.0, 100.0,false,false,false); static Parameter interfaceSingleHadronLimitCharm ("SingleHadronLimitCharm", "threshold for one-hadron decay of c-cluster", &HadronSelector::_limCharm, 0, 0.0, 0.0, 100.0,false,false,false); static Parameter interfaceSingleHadronLimitExotic ("SingleHadronLimitExotic", "threshold for one-hadron decay of exotic cluster", &HadronSelector::_limExotic, 0, 0.0, 0.0, 100.0,false,false,false); static Switch interfaceBelowThreshold ("BelowThreshold", "Option fo the selection of the hadrons if the cluster is below the pair threshold", &HadronSelector::belowThreshold_, 0, false, false); static SwitchOption interfaceBelowThresholdLightest (interfaceBelowThreshold, "Lightest", "Force cluster to decay to the lightest hadron with the appropriate flavours", 0); static SwitchOption interfaceBelowThresholdAll (interfaceBelowThreshold, "All", "Select from all the hadrons below the two hadron threshold according to their spin weights", 1); } double HadronSelector::mixingStateWeight(long id) const { switch(id) { case ParticleID::eta: return 0.5*probabilityMixing(_etamix ,1); case ParticleID::etaprime: return 0.5*probabilityMixing(_etamix ,2); case ParticleID::phi: return 0.5*probabilityMixing(_phimix ,1); case ParticleID::omega: return 0.5*probabilityMixing(_phimix ,2); case ParticleID::hprime_1: return 0.5*probabilityMixing(_h1mix ,1); case ParticleID::h_1: return 0.5*probabilityMixing(_h1mix ,2); case 10331: return 0.5*probabilityMixing(_f0mix ,1); case 10221: return 0.5*probabilityMixing(_f0mix ,2); case ParticleID::fprime_1: return 0.5*probabilityMixing(_f1mix ,1); case ParticleID::f_1: return 0.5*probabilityMixing(_f1mix ,2); case ParticleID::fprime_2: return 0.5*probabilityMixing(_f2mix ,1); case ParticleID::f_2: return 0.5*probabilityMixing(_f2mix ,2); case 10335: return 0.5*probabilityMixing(_eta2mix ,1); case 10225: return 0.5*probabilityMixing(_eta2mix ,2); // missing phi member of 13D1 should be here case 30223: return 0.5*probabilityMixing(_omhmix ,2); case 337: return 0.5*probabilityMixing(_ph3mix ,1); case 227: return 0.5*probabilityMixing(_ph3mix ,2); case 100331: return 0.5*probabilityMixing(_eta2mix ,1); case 100221: return 0.5*probabilityMixing(_eta2mix ,2); case 100333: return 0.5*probabilityMixing(_phi2Smix,1); case 100223: return 0.5*probabilityMixing(_phi2Smix,2); default: return 1./3.; } } void HadronSelector::doinit() { Interfaced::doinit(); // the default partons allowed // the quarks for ( int ix=1; ix<=5; ++ix ) { _partons.push_back(getParticleData(ix)); } // the diquarks for(unsigned int ix=1;ix<=5;++ix) { for(unsigned int iy=1; iy<=ix;++iy) { _partons.push_back(getParticleData(CheckId::makeDiquarkID(ix,iy))); } } // set the weights for the various excited mesons // set all to one to start with for (int l = 0; l < Lmax; ++l ) { for (int j = 0; j < Jmax; ++j) { for (int n = 0; n < Nmax; ++n) { _repwt[l][j][n] = 1.0; } } } // set the others from the relevant vectors for( int ix=0;ixid()]=1.; } _pwt[1] = _pwtDquark; _pwt[2] = _pwtUquark; _pwt[3] = _pwtSquark; _pwt[4] = _pwtCquark; _pwt[5] = _pwtBquark; _pwt[1103] = _pwtDIquark * _pwtDquark * _pwtDquark; _pwt[2101] = 0.5 * _pwtDIquark * _pwtUquark * _pwtDquark; _pwt[2203] = _pwtDIquark * _pwtUquark * _pwtUquark; _pwt[3101] = 0.5 * _pwtDIquark * _pwtSquark * _pwtDquark; _pwt[3201] = 0.5 * _pwtDIquark * _pwtSquark * _pwtUquark; _pwt[3303] = _pwtDIquark * _pwtSquark * _pwtSquark; // Commenting out heavy di-quark weights _pwt[4101] = 0.0; _pwt[4201] = 0.0; _pwt[4301] = 0.0; _pwt[4403] = 0.0; _pwt[5101] = 0.0; _pwt[5201] = 0.0; _pwt[5301] = 0.0; _pwt[5401] = 0.0; _pwt[5503] = 0.0; // find the maximum map::iterator pit = - max_element(_pwt.begin(),_pwt.end(),weightIsLess); + max_element(_pwt.begin(),_pwt.end(),weightIsLess); const double pmax = pit->second; for(pit=_pwt.begin(); pit!=_pwt.end(); ++pit) { pit->second/=pmax; } // construct the hadron tables constructHadronTable(); // for debugging // dumpTable(table()); } void HadronSelector::constructHadronTable() { // initialise the table _table.clear(); for(unsigned int ix=0; ix<_partons.size(); ++ix) { for(unsigned int iy=0; iy<_partons.size(); ++iy) { - if (!(DiquarkMatcher::Check(_partons[ix]->id()) + if (!(DiquarkMatcher::Check(_partons[ix]->id()) && DiquarkMatcher::Check(_partons[iy]->id()))) _table[make_pair(_partons[ix]->id(),_partons[iy]->id())] = KupcoData(); } } // get the particles from the event generator ParticleMap particles = generator()->particles(); // loop over the particles double maxdd(0.),maxss(0.),maxrest(0.); - for(ParticleMap::iterator it=particles.begin(); + for(ParticleMap::iterator it=particles.begin(); it!=particles.end(); ++it) { long pid = it->first; tPDPtr particle = it->second; int pspin = particle->iSpin(); // Don't include hadrons which are explicitly forbidden - if(find(_forbidden.begin(),_forbidden.end(),particle)!=_forbidden.end()) + if(find(_forbidden.begin(),_forbidden.end(),particle)!=_forbidden.end()) continue; // Don't include non-hadrons or antiparticles if(pid < 100) continue; // remove diffractive particles if(pspin == 0) continue; // K_0S and K_0L not made make K0 and Kbar0 if(pid==ParticleID::K_S0||pid==ParticleID::K_L0) continue; // Debugging options // Only include those with 2J+1 less than...5 if(_trial==2 && pspin >= 5) continue; // Only include those with 2J+1 less than...7 if(_trial==3 && pspin >= 7) continue; // Only include pions if(_trial==1 && pid!=111 && pid!=211) continue; // shouldn't be coloured if(particle->coloured()) continue; // Get the flavours - const int x4 = (pid/1000)%10; + const int x4 = (pid/1000)%10; const int x3 = (pid/100 )%10; const int x2 = (pid/10 )%10; const int x7 = (pid/1000000)%10; const bool wantSusy = x7 == 1 || x7 == 2; int flav1; int flav2; // Skip non-hadrons (susy particles, etc...) if(x3 == 0 || x2 == 0) continue; else if(x4 == 0) { // meson - flav1 = x2; - flav2 = x3; - } + flav1 = x2; + flav2 = x3; + } else { // baryon flav1 = CheckId::makeDiquarkID(x2,x3); flav2 = x4; } if (wantSusy) flav2 += 1000000 * x7; HadronInfo a(pid, particle, specialWeight(pid), particle->mass()); // set the weight to the number of spin states a.overallWeight = pspin; // identical light flavours if(flav1 == flav2 && flav1<=3) { // ddbar> uubar> admixture states if(flav1==1) { if(_topt != 0) a.overallWeight *= 0.5*a.swtef; _table[make_pair(1,1)].insert(a); _table[make_pair(2,2)].insert(a); if(_topt == 0 && a.overallWeight > maxdd) maxdd = a.overallWeight; } // load up ssbar> uubar> ddbar> admixture states else { a.wt = mixingStateWeight(pid); a.overallWeight *= a.wt; if(_topt != 0) a.overallWeight *= a.swtef; _table[make_pair(1,1)].insert(a); _table[make_pair(2,2)].insert(a); if(_topt == 0 && a.overallWeight > maxdd) maxdd = a.overallWeight; a.wt = (_topt != 0) ? 1.- 2.*a.wt : 1 - a.wt; if(a.wt > 0) { a.overallWeight = a.wt * a.swtef * pspin; _table[make_pair(3,3)].insert(a); if(_topt == 0 && a.overallWeight > maxss) maxss = a.overallWeight; } } } // light baryons with all quarks identical else if((flav1 == 1 && flav2 == 1103) || (flav1 == 1103 && flav2 == 1) || (flav1 == 2 && flav2 == 2203) || (flav1 == 2203 && flav2 == 2) || (flav1 == 3 && flav2 == 3303) || (flav1 == 3303 && flav2 == 3)) { if(_topt != 0) a.overallWeight *= 1.5*a.swtef; _table[make_pair(flav1,flav2)].insert(a); _table[make_pair(flav2,flav1)].insert(a); if(_topt == 0 && a.overallWeight > maxrest) maxrest = a.overallWeight; } // all other cases else { if(_topt != 0) a.overallWeight *=a.swtef; _table[make_pair(flav1,flav2)].insert(a); if(flav1 != flav2) _table[make_pair(flav2,flav1)].insert(a); if(_topt == 0 && a.overallWeight > maxrest) maxrest = a.overallWeight; } } // Account for identical combos of diquark/quarks and symmetrical elements // e.g. U UD = D UU - HadronTable::iterator tit; + HadronTable::iterator tit; for(tit=_table.begin();tit!=_table.end();++tit) { if(tit->first.first>ParticleID::c) continue; if(!DiquarkMatcher::Check(tit->first.second)) continue; long k, l, sub; if(tit->first.second>=ParticleID::bd_0) { k = ParticleID::b; sub = ParticleID::bd_0/100; } else if(tit->first.second>=ParticleID::cd_0) { k = ParticleID::c; sub = ParticleID::cd_0/100; } else if(tit->first.second>=ParticleID::sd_0) { k = ParticleID::s; sub = ParticleID::sd_0/100; } else if(tit->first.second>=ParticleID::ud_0) { k = ParticleID::u; sub = ParticleID::ud_0/100; } else if(tit->first.second==ParticleID::dd_1) { k = ParticleID::d; sub = ParticleID::dd_1/100; } else continue; sub=tit->first.second/100-sub+1; if(sub > tit->first.first) { l = 1000*sub+100*tit->first.first+1; } else if(sub==tit->first.first) { l = 1000*sub+ 100*tit->first.first+3; } else { l = 100*sub +1000*tit->first.first+1; } if(tit->second.empty()) { pair newpair(k,l); tit->second=_table[newpair]; newpair=make_pair(tit->first.second,tit->first.first); _table[newpair]=tit->second; }; } // normalise weights to one for first option if(_topt == 0) { HadronTable::const_iterator tit; KupcoData::iterator it; for(tit=_table.begin();tit!=_table.end();++tit) { double weight; if(tit->first.first==tit->first.second) { if(tit->first.first==1||tit->first.first==2) weight=1./maxdd; else if (tit->first.first==3) weight=1./maxss; else weight=1./maxrest; } else weight=1./maxrest; for(it = tit->second.begin(); it!=tit->second.end(); ++it) { it->rescale(weight); } } } } double HadronSelector::specialWeight(long id) const { - const int pspin = id % 10; + const int pspin = id % 10; // Only K0L and K0S have pspin == 0, should // not get them until Decay step assert( pspin != 0 ); // Baryon : J = 1/2 or 3/2 if(pspin == 2) { // Singlet (Lambda-like) baryon - if( (id/100)%10 < (id/10 )%10 ) return sqr(_sngWt); + if( (id/100)%10 < (id/10 )%10 ) return sqr(_sngWt); // octet else return 1.; - } + } // Decuplet baryon else if (pspin == 4) { return sqr(_decWt); - } + } // Meson else if(pspin % 2 == 1) { // Total angular momentum int j = (pspin - 1) / 2; // related to Orbital angular momentum l int nl = (id/10000 )%10; - int l = -999; + int l = -999; int n = (id/100000)%10; // Radial excitation if(j == 0) l = nl; else if(nl == 0) l = j - 1; else if(nl == 1 || nl == 2) l = j; else if(nl == 3) l = j + 1; // Angular or Radial excited meson if((l||j||n) && l>=0 && l= 5/2 (ispin >= 6), haven't got those return 1.0; } -int HadronSelector::signHadron(tcPDPtr idQ1, tcPDPtr idQ2, +int HadronSelector::signHadron(tcPDPtr idQ1, tcPDPtr idQ2, tcPDPtr hadron) const { // This method receives in input three PDG ids, whose the - // first two have proper signs (corresponding to particles, id > 0, + // first two have proper signs (corresponding to particles, id > 0, // or antiparticles, id < 0 ), whereas the third one must // be always positive (particle not antiparticle), // corresponding to: // --- quark-antiquark, or antiquark-quark, or // quark-diquark, or diquark-quark, or // antiquark-antidiquark, or antidiquark-antiquark // for the first two input (idQ1, idQ2); - // --- meson or baryon for the third input (idHad): + // --- meson or baryon for the third input (idHad): // The method returns: // --- + 1 if the two partons (idQ1, idQ2) are exactly // the constituents for the hadron idHad; // --- - 1 if the two partons (idQ1, idQ2) are exactly // the constituents for the anti-hadron -idHad; // --- + 0 otherwise. // The method it is therefore useful to decide the - // sign of the id of the produced hadron as appeared - // in the vector _vecHad (where only hadron idHad > 0 are present) + // sign of the id of the produced hadron as appeared + // in the vector _vecHad (where only hadron idHad > 0 are present) // given the two constituent partons. int sign = 0; long idHad = hadron->id(); assert(idHad > 0); int chargeIn = idQ1->iCharge() + idQ2->iCharge(); int chargeOut = hadron->iCharge(); // same charge if( chargeIn == chargeOut && chargeIn !=0 ) sign = +1; else if(chargeIn == -chargeOut && chargeIn !=0 ) sign = -1; - else if(chargeIn == 0 && chargeOut == 0 ) { + else if(chargeIn == 0 && chargeOut == 0 ) { // In the case of same null charge, there are four cases: // i) K0-like mesons, B0-like mesons, Bs-like mesons - // the PDG convention is to consider them "antiparticle" (idHad < 0) + // the PDG convention is to consider them "antiparticle" (idHad < 0) // if the "dominant" (heavier) flavour (respectively, s, b) // is a quark (idQ > 0): for instance, B0s = (b, sbar) has id < 0 // Remember that there is an important exception for K0L (id=130) and // K0S (id=310): they don't have antiparticles, therefore idHad > 0 // always. We use below the fact that K0L and K0S are the unique // hadrons having 0 the first (less significant) digit of their id. // 2) D0-like mesons: the PDG convention is to consider them "particle" // (idHad > 0) if the charm flavour is carried by a c: (c,ubar) has id>0 // 3) the remaining mesons should not have antiparticle, therefore their // sign is always positive. - // 4) for baryons, that is when one of idQ1 and idQ2 is a (anti-) quark and + // 4) for baryons, that is when one of idQ1 and idQ2 is a (anti-) quark and // the other one is a (anti-) diquark the sign is negative when both // constituents are "anti", that is both with id < 0; positive otherwise. // meson if(abs(int(idQ1->iColour()))== 3 && abs(int(idQ2->iColour())) == 3 && !DiquarkMatcher::Check(idQ1->id()) && !DiquarkMatcher::Check(idQ2->id())) { int idQa = abs(idQ1->id()); - int idQb = abs(idQ2->id()); + int idQb = abs(idQ2->id()); int dominant = idQ2->id(); if(idQa > idQb) { swap(idQa,idQb); dominant = idQ1->id(); } if((idQa==ParticleID::d && idQb==ParticleID::s) || (idQa==ParticleID::d && idQb==ParticleID::b) || - (idQa==ParticleID::s && idQb==ParticleID::b)) { + (idQa==ParticleID::s && idQb==ParticleID::b)) { // idHad%10 is zero for K0L,K0S if (dominant < 0 || idHad%10 == 0) sign = +1; else if(dominant > 0) sign = -1; - } + } else if((idQa==ParticleID::u && idQb==ParticleID::c) || (idQa==ParticleID::u && idQb==ParticleID::t) || (idQa==ParticleID::c && idQb==ParticleID::t)) { if (dominant > 0) sign = +1; else if(dominant < 0) sign = -1; - } + } else if(idQa==idQb) sign = +1; // sets sign for Susy particles else sign = (dominant > 0) ? +1 : -1; } // baryon else if(DiquarkMatcher::Check(idQ1->id()) || DiquarkMatcher::Check(idQ2->id())) { if (idQ1->id() > 0 && idQ2->id() > 0) sign = +1; else if(idQ1->id() < 0 && idQ2->id() < 0) sign = -1; } } if (sign == 0) { - cerr << "Could not work out sign for " - << idQ1->PDGName() << ' ' - << idQ2->PDGName() << " => " + cerr << "Could not work out sign for " + << idQ1->PDGName() << ' ' + << idQ2->PDGName() << " => " << hadron->PDGName() << '\n'; assert(false); } return sign; } pair HadronSelector::lightestHadronPair(tcPDPtr ptr1, tcPDPtr ptr2, tcPDPtr ptr3) const { // throw exception of id3!=0 as doesn't work - if ( ptr3 ) throw Exception() + if ( ptr3 ) throw Exception() << "ptr3!=0 not yet implemented in HadronSelector::lightestHadronPair" << Exception::abortnow; // charge int totalcharge = ptr1->iCharge() + ptr2->iCharge(); if ( ptr3 ) totalcharge += ptr3->iCharge(); tcPDPtr vIdHad1[2]={tcPDPtr(),tcPDPtr()},vIdHad2[2]={tcPDPtr(),tcPDPtr()}; bool vOk[2] = {false, false}; Energy vMassPair[2] = { ZERO, ZERO }; for (int i = 0; i < 2; i++) { tcPDPtr idPartner = i==0 ? getParticleData(ParticleID::d) : getParticleData(ParticleID::u); // Change sign to idPartner (transform it into a anti-quark) if it is not // possible to form a meson or a baryon. assert (ptr1 && idPartner); if (!CheckId::canBeHadron(ptr1, idPartner)) idPartner = idPartner->CC(); vIdHad1[i] = lightestHadron(ptr1, idPartner); vIdHad2[i] = lightestHadron(ptr2, idPartner->CC()); - if ( vIdHad1[i] && vIdHad2[i] && + if ( vIdHad1[i] && vIdHad2[i] && vIdHad1[i]->iCharge() + vIdHad2[i]->iCharge() == totalcharge ) { vOk[i] = true; vMassPair[i] = vIdHad1[i]->mass() + vIdHad2[i]->mass(); - } + } } // Take the lightest pair compatible with charge conservation. - if ( vOk[0] && ( ! vOk[1] || vMassPair[0] <= vMassPair[1] ) ) { + if ( vOk[0] && ( ! vOk[1] || vMassPair[0] <= vMassPair[1] ) ) { return make_pair(vIdHad1[0],vIdHad2[0]); - } - else if ( vOk[1] && ( ! vOk[0] || vMassPair[1] < vMassPair[0] ) ) { + } + else if ( vOk[1] && ( ! vOk[0] || vMassPair[1] < vMassPair[0] ) ) { return make_pair(vIdHad1[1],vIdHad2[1]); } else { return make_pair(tcPDPtr(),tcPDPtr()); - } + } } Energy HadronSelector::massLightestBaryonPair(tcPDPtr ptr1, tcPDPtr ptr2) const { // Make sure that we don't have any diquarks as input, return arbitrarily // large value if we do - Energy currentSum = Constants::MaxEnergy; + Energy currentSum = Constants::MaxEnergy; for(unsigned int ix=0; ix<_partons.size(); ++ix) { if(!DiquarkMatcher::Check(_partons[ix]->id())) continue; - HadronTable::const_iterator + HadronTable::const_iterator tit1=_table.find(make_pair(abs(ptr1->id()),_partons[ix]->id())), tit2=_table.find(make_pair(_partons[ix]->id(),abs(ptr2->id()))); if( tit1==_table.end() || tit2==_table.end()) continue; if(tit1->second.empty()||tit2->second.empty()) continue; Energy s = tit1->second.begin()->mass + tit2->second.begin()->mass; if(currentSum > s) currentSum = s; } return currentSum; } tcPDPtr HadronSelector::lightestHadron(tcPDPtr ptr1, tcPDPtr ptr2, #ifndef NDEBUG tcPDPtr ptr3) const { #else tcPDPtr ) const { #endif // The method assumes ptr3 == 0 rest not implemented assert(ptr1 && ptr2 && !ptr3); // find entry in the table pair ids = make_pair(abs(ptr1->id()),abs(ptr2->id())); HadronTable::const_iterator tit=_table.find(ids); // throw exception if flavours wrong - if (tit==_table.end()) - throw Exception() << "Could not find " - << ids.first << ' ' << ids.second + if (tit==_table.end()) + throw Exception() << "Could not find " + << ids.first << ' ' << ids.second << " in _table. " << "In HadronSelector::lightestHadron()" << Exception::eventerror; if(tit->second.empty()) throw Exception() << "HadronSelector::lightestHadron " - << "could not find any hadrons containing " + << "could not find any hadrons containing " << ptr1->id() << ' ' << ptr2->id() << '\n' - << tit->first.first << ' ' + << tit->first.first << ' ' << tit->first.second << Exception::eventerror; // find the lightest hadron int sign = signHadron(ptr1,ptr2,tit->second.begin()->ptrData); - tcPDPtr candidate = sign > 0 ? + tcPDPtr candidate = sign > 0 ? tit->second.begin()->ptrData : tit->second.begin()->ptrData->CC(); // \todo 20 GeV limit is temporary fudge to let SM particles go through. // \todo Use isExotic instead? - if (candidate->mass() > 20*GeV + if (candidate->mass() > 20*GeV && candidate->mass() < ptr1->constituentMass() + ptr2->constituentMass()) { generator()->log() << "HadronSelector::lightestHadron: " - << "chosen candidate " << candidate->PDGName() + << "chosen candidate " << candidate->PDGName() << " is lighter than its constituents " << ptr1->PDGName() << ", " << ptr2->PDGName() << '\n' << candidate->mass()/GeV << " < " << ptr1->constituentMass()/GeV << " + " << ptr2->constituentMass()/GeV << '\n' << "Check your particle data tables.\n"; assert(false); } return candidate; } -vector > +vector > HadronSelector::hadronsBelowThreshold(Energy threshold, tcPDPtr ptr1, tcPDPtr ptr2, #ifndef NDEBUG tcPDPtr ptr3) const { #else tcPDPtr ) const { #endif // The method assumes ptr3 == 0 rest not implemented assert(ptr1 && ptr2 && !ptr3); // find entry in the table pair ids = make_pair(abs(ptr1->id()),abs(ptr2->id())); HadronTable::const_iterator tit=_table.find(ids); // throw exception if flavours wrong - if (tit==_table.end()) - throw Exception() << "Could not find " - << ids.first << ' ' << ids.second + if (tit==_table.end()) + throw Exception() << "Could not find " + << ids.first << ' ' << ids.second << " in _table. " << "In HadronSelector::hadronsBelowThreshold()" << Exception::eventerror; if(tit->second.empty()) throw Exception() << "HadronSelector::hadronsBelowThreshold() " - << "could not find any hadrons containing " + << "could not find any hadrons containing " << ptr1->id() << ' ' << ptr2->id() << '\n' - << tit->first.first << ' ' + << tit->first.first << ' ' << tit->first.second << Exception::eventerror; vector > candidates; KupcoData::const_iterator hit = tit->second.begin(); // find the hadrons while(hit!=tit->second.end()&&hit->massptrData); tcPDPtr candidate = sign > 0 ? hit->ptrData : hit->ptrData->CC(); // \todo 20 GeV limit is temporary fudge to let SM particles go through. // \todo Use isExotic instead? - if (candidate->mass() > 20*GeV + if (candidate->mass() > 20*GeV && candidate->mass() < ptr1->constituentMass() + ptr2->constituentMass()) { generator()->log() << "HadronSelector::hadronsBelowTheshold: " - << "chosen candidate " << candidate->PDGName() + << "chosen candidate " << candidate->PDGName() << " is lighter than its constituents " << ptr1->PDGName() << ", " << ptr2->PDGName() << '\n' << candidate->mass()/GeV << " < " << ptr1->constituentMass()/GeV << " + " << ptr2->constituentMass()/GeV << '\n' << "Check your particle data tables.\n"; assert(false); - } + } candidates.push_back(make_pair(candidate,hit->overallWeight)); ++hit; } return candidates; } tcPDPtr HadronSelector::chooseSingleHadron(tcPDPtr par1, tcPDPtr par2, Energy mass) const { // Determine the sum of the nominal masses of the two lightest hadrons // with the right flavour numbers as the cluster under consideration. - // Notice that we don't need real masses (drawn by a Breit-Wigner + // Notice that we don't need real masses (drawn by a Breit-Wigner // distribution) because the lightest pair of hadrons does not involve // any broad resonance. Energy threshold = massLightestHadronPair(par1,par2); // Special: it allows one-hadron decays also above threshold. - if (CheckId::isExotic(par1,par2)) + if (CheckId::isExotic(par1,par2)) threshold *= (1.0 + UseRandom::rnd()*_limExotic); - else if (CheckId::hasBottom(par1,par2)) + else if (CheckId::hasBottom(par1,par2)) threshold *= (1.0 + UseRandom::rnd()*_limBottom); - else if (CheckId::hasCharm(par1,par2)) + else if (CheckId::hasCharm(par1,par2)) threshold *= (1.0 + UseRandom::rnd()*_limCharm); - + // only do one hadron decay is mass less than the threshold if(mass>=threshold) return tcPDPtr(); // select the hadron tcPDPtr hadron; // old option pick the lightest hadron if(belowThreshold_ == 0) { hadron= lightestHadron(par1,par2); + cout<<"below threshold"< > hadrons = + vector > hadrons = hadronsBelowThreshold(threshold,par1,par2); if(hadrons.size()==1) { hadron = hadrons[0].first; } else if(hadrons.empty()) { hadron= lightestHadron(par1,par2); } else { double totalWeight=0.; for(unsigned int ix=0;ix #include #include #include #include #include "HadronSelector.fh" namespace Herwig { using namespace ThePEG; /**\ingroup Hadronization * \class HadronSelector * \brief This class selects the hadron flavours of a cluster decay. * \author Philip Stephens * \author Alberto Ribon * \author Peter Richardson * * This is the base class for the selection of either a pair of hadrons, or - * in some cases a single hadron. The different approaches which were - * previously implemented in this class are now implemented in the + * in some cases a single hadron. The different approaches which were + * previously implemented in this class are now implemented in the * HwppSelector and Hw64Selector which inherit from this class. * * This class implements a number of methods which are needed by all models * and in addition contains the weights for the different meson multiplets and * mixing of the light \f$I=0\f$ mesons. * * @see \ref HadronSelectorInterfaces "The interfaces" * defined for HadronSelector. * @see HwppSelector * @see Hw64Selector */ class HadronSelector: public Interfaced { public: /** \ingroup Hadronization - * \class HadronInfo + * \class HadronInfo * \brief Class used to store all the hadron information for easy access. * \author Philip Stephens - * + * * Note that: - * - the hadrons in _table can be filled in any ordered + * - the hadrons in _table can be filled in any ordered * w.r.t. the mass value, and flavours for different * groups (for instance, (u,s) hadrons don't need to - * be placed after (d,s) or any other flavour), but + * be placed after (d,s) or any other flavour), but * all hadrons with the same flavours must be consecutive * ( for instance you cannot alternate hadrons of type * (d,s) with those of flavour (u,s) ). * Furthermore, it is assumed that particle and antiparticle * have the same weights, and therefore only one of them * must be entered in the table: we have chosen to refer * to the particle, defined as PDG id > 0, although if * an anti-particle is provided in input it is automatically * transform to its particle, simply by taking the modulus * of its id. */ - class HadronInfo { + class HadronInfo { public: - + /** * Constructor * @param idin The PDG code of the hadron * @param datain The pointer to the ParticleData object * @param swtin The singlet/decuplet/orbital factor * @param massin The mass of the hadron */ HadronInfo(long idin=0, tPDPtr datain=tPDPtr(), double swtin=1., Energy massin=ZERO) : id(idin), ptrData(datain), swtef(swtin), wt(1.0), overallWeight(0.0), mass(massin) {} - + /** * Comparision operator on mass */ bool operator<(const HadronInfo &x) const { if(mass!=x.mass) return mass < x.mass; else return id < x.id; } - + /** * The hadrons id. */ long id; - + /** * pointer to ParticleData, to get the spin, etc... */ tPDPtr ptrData; - + /** - * singlet/decuplet/orbital factor + * singlet/decuplet/orbital factor */ double swtef; - + /** * mixing factor */ - double wt; - + double wt; + /** * (2*J+1)*wt*swtef */ double overallWeight; - + /** * The hadrons mass */ Energy mass; - + /** * Rescale the weight for a given hadron */ - void rescale(double x) const { - const_cast(this)->overallWeight *= x; + void rescale(double x) const { + const_cast(this)->overallWeight *= x; } - + /** * Friend method used to print the value of a table element. */ - friend PersistentOStream & operator<< (PersistentOStream & os, + friend PersistentOStream & operator<< (PersistentOStream & os, const HadronInfo & hi ) { - os << hi.id << hi.ptrData << hi.swtef << hi.wt + os << hi.id << hi.ptrData << hi.swtef << hi.wt << hi.overallWeight << ounit(hi.mass,GeV); return os; } - + /** * debug output */ friend ostream & operator<< (ostream & os, const HadronInfo & hi ) { os << std::scientific << std::showpoint << std::setprecision(4) << setw(2) << hi.id << '\t' << hi.swtef << '\t' << hi.wt << '\t' << hi.overallWeight << '\t' << ounit(hi.mass,GeV); return os; } - + /** * Friend method used to read in the value of a table element. */ - friend PersistentIStream & operator>> (PersistentIStream & is, + friend PersistentIStream & operator>> (PersistentIStream & is, HadronInfo & hi ) { - is >> hi.id >> hi.ptrData >> hi.swtef >> hi.wt + is >> hi.id >> hi.ptrData >> hi.swtef >> hi.wt >> hi.overallWeight >> iunit(hi.mass,GeV); return is; } }; - + /** \ingroup Hadronization * \class Kupco * \brief Class designed to make STL routines easy to use. * \author Philip Stephens * - * This class is used to generate a list of the hadron pairs which can + * This class is used to generate a list of the hadron pairs which can * be produced that allows easy traversal and quick access. */ class Kupco { - + public: - + /** * Constructor * @param inidQ PDG code of the quark drawn from the vacuum. * @param inhad1 ParticleData for the first hadron produced. * @param inhad2 ParticleData for the second hadron produced. - * @param inwgt The weight for the hadron pair + * @param inwgt The weight for the hadron pair */ Kupco(tcPDPtr inidQ,tcPDPtr inhad1,tcPDPtr inhad2, Energy inwgt) : idQ(inidQ),hadron1(inhad1),hadron2(inhad2),weight(inwgt) {} - + /** * id of the quark drawn from the vacuum. */ tcPDPtr idQ; - + /** * The ParticleData object for the first hadron produced. */ tcPDPtr hadron1; - + /** * The ParticleData object for the second hadron produced. */ tcPDPtr hadron2; - + /** * Weight factor of this componation. */ Energy weight; }; public: /** * The helper classes */ //@{ /** * The type is used to contain all the hadrons info of a given flavour. */ typedef set KupcoData; //@} /** * The hadron table type. */ typedef map,KupcoData> HadronTable; public: /** * The default constructor. */ HadronSelector(unsigned int); /** * Method to return a pair of hadrons given the PDG codes of * two or three constituents * @param cluMass The mass of the cluster * @param par1 The first constituent * @param par2 The second constituent * @param par3 The third constituent */ - virtual pair chooseHadronPair(const Energy cluMass, tcPDPtr par1, + virtual pair chooseHadronPair(const Energy cluMass, tcPDPtr par1, tcPDPtr par2,tcPDPtr par3 = PDPtr()) const = 0; /** * Select the single hadron for a cluster decay * return null pointer if not a single hadron decay * @param par1 1st constituent * @param par2 2nd constituent * @param mass Mass of the cluster */ tcPDPtr chooseSingleHadron(tcPDPtr par1, tcPDPtr par2, Energy mass) const; /** * This returns the lightest pair of hadrons given by the flavours. * * Given the two (or three) constituents of a cluster, it returns * the two lightest hadrons with proper flavour numbers. * Furthermore, the first of the two hadrons must have the constituent with - * par1, and the second must have the constituent with par2. + * par1, and the second must have the constituent with par2. * \todo At the moment it does *nothing* in the case that also par3 is present. * - * The method is implemented by calling twice lightestHadron, + * The method is implemented by calling twice lightestHadron, * once with (par1,quarktopick->CC()) ,and once with (par2,quarktopick) * where quarktopick is either the pointer to - * d or u quarks . In fact, the idea is that whatever the flavour of par1 + * d or u quarks . In fact, the idea is that whatever the flavour of par1 * and par2, no matter if (anti-)quark or (anti-)diquark, the lightest - * pair of hadrons containing flavour par1 and par2 will have either + * pair of hadrons containing flavour par1 and par2 will have either * flavour d or u, being the lightest quarks. - * The method returns the pair (PDPtr(),PDPtr()) if anything goes wrong. + * The method returns the pair (PDPtr(),PDPtr()) if anything goes wrong. * - * \todo The method assumes par3 == PDPtr() (otherwise we don't know how to proceed: a - * possible, trivial way would be to randomly select two of the three + * \todo The method assumes par3 == PDPtr() (otherwise we don't know how to proceed: a + * possible, trivial way would be to randomly select two of the three * (anti-)quarks and treat them as a (anti-)diquark, reducing the problem * to two components as treated below. * In the normal (two components) situation, the strategy is the following: - * treat in the same way the two possibilities: (d dbar) (i=0) and + * treat in the same way the two possibilities: (d dbar) (i=0) and * (u ubar) (i=1) as the pair quark-antiquark necessary to form a * pair of hadrons containing the input flavour par1 and par2; finally, * select the one that produces the lightest pair of hadrons, compatible * with the charge conservation constraint. */ pair lightestHadronPair(tcPDPtr ptr1, tcPDPtr ptr2, tcPDPtr ptr3 = PDPtr ()) const; /** * Returns the mass of the lightest pair of hadrons with the given particles * @param ptr1 is the first constituent - * @param ptr2 is the second constituent + * @param ptr2 is the second constituent * @param ptr3 is the third constituent */ Energy massLightestHadronPair(tcPDPtr ptr1, tcPDPtr ptr2, tcPDPtr ptr3 = PDPtr ()) const { pair pairData = lightestHadronPair(ptr1, ptr2, ptr3); if ( ! pairData.first || ! pairData.second ) return ZERO; - return ( pairData.first->mass() + pairData.second->mass() ); + return ( pairData.first->mass() + pairData.second->mass() ); } /** * Returns the lightest hadron formed by the given particles. * * Given the id of two (or three) constituents of a cluster, it returns * the lightest hadron with proper flavour numbers. * At the moment it does *nothing* in the case that also 'ptr3' present. * @param ptr1 is the first constituent - * @param ptr2 is the second constituent + * @param ptr2 is the second constituent * @param ptr3 is the third constituent */ tcPDPtr lightestHadron(tcPDPtr ptr1, tcPDPtr ptr2, tcPDPtr ptr3 = PDPtr ()) const; /** * Returns the hadrons below the constituent mass threshold formed by the given particles, * together with their total weight * * Given the id of two (or three) constituents of a cluster, it returns * the lightest hadron with proper flavour numbers. * At the moment it does *nothing* in the case that also 'ptr3' present. * @param threshold The theshold * @param ptr1 is the first constituent - * @param ptr2 is the second constituent + * @param ptr2 is the second constituent * @param ptr3 is the third constituent */ - vector > + vector > hadronsBelowThreshold(Energy threshold, tcPDPtr ptr1, tcPDPtr ptr2, tcPDPtr ptr3 = PDPtr ()) const; /** + * Functions to update the species of quark produced from the vacuum + * and return the value to ClusterDecayer + */ + virtual long getQ() = 0; + + /** * Return the nominal mass of the hadron returned by lightestHadron() * @param ptr1 is the first constituent - * @param ptr2 is the second constituent - * @param ptr3 is the third constituent + * @param ptr2 is the second constituent + * @param ptr3 is the third constituent */ Energy massLightestHadron(tcPDPtr ptr1, tcPDPtr ptr2, #ifndef NDEBUG tcPDPtr ptr3 = PDPtr ()) const { #else tcPDPtr = PDPtr ()) const { #endif - // The method assumes ptr3 == empty + // The method assumes ptr3 == empty assert(!(ptr3)); // find entry in the table pair ids(abs(ptr1->id()),abs(ptr2->id())); HadronTable::const_iterator tit=_table.find(ids); // throw exception if flavours wrong if(tit==_table.end()||tit->second.empty()) throw Exception() << "HadronSelector::massLightestHadron " - << "failed for particle" << ptr1->id() << " " - << ptr2->id() + << "failed for particle" << ptr1->id() << " " + << ptr2->id() << Exception::eventerror; // return the mass return tit->second.begin()->mass; } /** * Returns the mass of the lightest pair of baryons. * @param ptr1 is the first constituent - * @param ptr2 is the second constituent + * @param ptr2 is the second constituent */ Energy massLightestBaryonPair(tcPDPtr ptr1, tcPDPtr ptr2) const; - + /** * Return the weights for the different quarks and diquarks */ //@{ /** * The down quark weight. */ double pwtDquark() const { return _pwtDquark; - } + } /** * The up quark weight. */ - double pwtUquark() const { + double pwtUquark() const { return _pwtUquark; } /** * The strange quark weight. */ - double pwtSquark() const { + double pwtSquark() const { return _pwtSquark; } /** * The charm quark weight. */ - double pwtCquark() const { + double pwtCquark() const { return _pwtCquark; } /** * The bottom quark weight. */ - double pwtBquark() const { + double pwtBquark() const { return _pwtBquark; - } - + } + /** * The diquark weight. */ double pwtDIquark() const { return _pwtDIquark; } //@} public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * * The array _repwt is initialized using the interfaces to set different * weights for different meson multiplets and the constructHadronTable() * method called to complete the construction of the hadron tables. * * @throws InitException if object could not be initialized properly. */ virtual void doinit(); //@} protected: /** * Construct the table of hadron data * This is the main method to initialize the hadron data (mainly the - * weights associated to each hadron, taking into account its spin, + * weights associated to each hadron, taking into account its spin, * eventual isoscalar-octect mixing, singlet-decuplet factor). This is * the method that one should update when new or updated hadron data is - * available. + * available. * - * This class implements the construction of the basic table but can be + * This class implements the construction of the basic table but can be * overridden if needed in inheriting classes. * - * The rationale for factors used for diquarks involving different quarks can + * The rationale for factors used for diquarks involving different quarks can * be can be explained by taking a prototype example that in the exact SU(2) limit, * in which: - * \f[m_u=m_d\f] + * \f[m_u=m_d\f] * \f[M_p=M_n=M_\Delta\f] * and we will have equal numbers of u and d quarks produced. - * Suppose that we weight 1 the diquarks made of the same + * Suppose that we weight 1 the diquarks made of the same * quark and 1/2 those made of different quarks, the fractions * of u and d baryons (p, n, Delta) we get are the following: * - \f$\Delta^{++}\f$: 1 possibility only u uu with weight 1 * - \f$\Delta^- \f$: 1 possibility only d dd with weight 1 * - \f$p,\Delta^+ \f$: 2 possibilities u ud with weight 1/2 * d uu with weight 1 * - \f$n,\Delta^0 \f$: 2 possibilities d ud with weight 1/2 * u dd with weight 1 - * In the latter two cases, we have to take into account the + * In the latter two cases, we have to take into account the * fact that p and n have spin 1/2 whereas Delta+ and Delta0 - * have spin 3/2 therefore from phase space we get a double weight + * have spin 3/2 therefore from phase space we get a double weight * for Delta+ and Delta0 relative to p and n respectively. * Therefore the relative amount of these baryons that is * produced is the following: * # p = # n = ( 1/2 + 1 ) * 1/3 = 1/2 * # Delta++ = # Delta- = 1 = ( 1/2 + 1) * 2/3 # Delta+ = # Delta0 * which is correct, and therefore the weight 1/2 for the * diquarks of different types of quarks is justified (at least * in this limit of exact SU(2) ). */ virtual void constructHadronTable(); /** * Access to the table of hadrons */ const HadronTable & table() const { return _table; } - + /** * Access to the list of partons */ const vector & partons() const { return _partons; } /** * Access the parton weights */ double pwt(long pid) const { map::const_iterator it = _pwt.find(abs(pid)); assert( it != _pwt.end() ); return it->second; } /** * Methods for the mixing of \f$I=0\f$ mesons */ //@{ /** * Return the probability of mixing for Octet-Singlet isoscalar mixing, - * the probability of the + * the probability of the * \f$\frac1{\sqrt{2}}(|u\bar{u}\rangle + |d\bar{d}\rangle)\f$ component * is returned. * @param angleMix The mixing angle in degrees (not radians) * @param order is 0 for no mixing, 1 for the first resonance of a pair, * 2 for the second one. * The mixing is defined so that for example with \f$\eta-\eta'\f$ mixing where * the mixing angle is \f$\theta=-23^0$ with $\eta\f$ as the first particle * and \f$\eta'\f$ the second one. - * The convention used is + * The convention used is * \f[\eta = \cos\theta|\eta_{\rm octet }\rangle * -\sin\theta|\eta_{\rm singlet}\rangle\f] * \f[\eta' = \sin\theta|\eta_{\rm octet }\rangle * -\cos\theta|\eta_{\rm singlet}\rangle\f] - * with + * with * \f[|\eta_{\rm singlet}\rangle = \frac1{\sqrt{3}} * \left[|u\bar{u}\rangle + |d\bar{d}\rangle + |s\bar{s}\rangle\right]\f] * \f[|\eta_{\rm octet }\rangle = \frac1{\sqrt{6}} * \left[|u\bar{u}\rangle + |d\bar{d}\rangle - 2|s\bar{s}\rangle\right]\f] */ double probabilityMixing(const double angleMix, const int order) const { static double convert=Constants::pi/180.0; - if (order == 1) + if (order == 1) return sqr( cos( angleMix*convert + atan( sqrt(2.0) ) ) ); - else if (order == 2) + else if (order == 2) return sqr( sin( angleMix*convert + atan( sqrt(2.0) ) ) ); - else + else return 1.; } /** * Returns the weight of given mixing state. * @param id The PDG code of the meson */ - virtual double mixingStateWeight(long id) const; + virtual double mixingStateWeight(long id) const; //@} /** * Calculates a special weight specific to a given hadron. * @param id The PDG code of the hadron */ double specialWeight(long id) const; /** * This method returns the proper sign ( > 0 hadron; < 0 anti-hadron ) * for the input PDG id idHad > 0, suppose to be made by the * two constituent particle pointers: par1 and par2 (both with proper sign). */ int signHadron(tcPDPtr ptr1, tcPDPtr ptr2, tcPDPtr hadron) const; private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ HadronSelector & operator=(const HadronSelector &) = delete; private: /** * The PDG codes of the constituent particles allowed */ vector _partons; /** * The PDG codes of the hadrons which cannot be produced in the hadronization */ vector _forbidden; /** * The weights for the different quarks and diquarks */ //@{ /** * The probability of producting a down quark. */ double _pwtDquark; /** * The probability of producting an up quark. */ double _pwtUquark; /** * The probability of producting a strange quark. */ double _pwtSquark; /** * The probability of producting a charm quark. */ double _pwtCquark; /** * The probability of producting a bottom quark. */ double _pwtBquark; /** * The probability of producting a diquark. */ double _pwtDIquark; /** * Weights for quarks and diquarks. */ map _pwt; //@} /** * The mixing angles for the \f$I=0\f$ mesons containing light quarks */ //@{ /** - * The \f$\eta-\eta'\f$ mixing angle + * The \f$\eta-\eta'\f$ mixing angle */ double _etamix; /** * The \f$\phi-\omega\f$ mixing angle */ double _phimix; /** * The \f$h_1'-h_1\f$ mixing angle */ double _h1mix; /** * The \f$f_0(1710)-f_0(1370)\f$ mixing angle */ double _f0mix; /** * The \f$f_1(1420)-f_1(1285)\f$ mixing angle */ double _f1mix; /** * The \f$f'_2-f_2\f$ mixing angle */ double _f2mix; /** * The \f$\eta_2(1870)-\eta_2(1645)\f$ mixing angle */ double _eta2mix; /** * The \f$\phi(???)-\omega(1650)\f$ mixing angle */ double _omhmix; /** * The \f$\phi_3-\omega_3\f$ mixing angle */ double _ph3mix; /** * The \f$\eta(1475)-\eta(1295)\f$ mixing angle */ double _eta2Smix; /** * The \f$\phi(1680)-\omega(1420)\f$ mixing angle */ double _phi2Smix; //@} /** * The weights for the various meson multiplets to be used to supress the * production of particular states */ //@{ /** * The weights for the \f$\phantom{1}^1S_0\f$ multiplets */ vector _weight1S0; /** * The weights for the \f$\phantom{1}^3S_1\f$ multiplets */ vector _weight3S1; /** * The weights for the \f$\phantom{1}^1P_1\f$ multiplets */ vector _weight1P1; /** * The weights for the \f$\phantom{1}^3P_0\f$ multiplets */ vector _weight3P0; /** * The weights for the \f$\phantom{1}^3P_1\f$ multiplets */ vector _weight3P1; /** * The weights for the \f$\phantom{1}^3P_2\f$ multiplets */ vector _weight3P2; /** * The weights for the \f$\phantom{1}^1D_2\f$ multiplets */ vector _weight1D2; /** * The weights for the \f$\phantom{1}^3D_1\f$ multiplets */ vector _weight3D1; /** * The weights for the \f$\phantom{1}^3D_2\f$ multiplets */ vector _weight3D2; /** * The weights for the \f$\phantom{1}^3D_3\f$ multiplets */ vector _weight3D3; //@} /** * The weights for the excited meson multiplets */ vector > > _repwt; /** * Singlet and Decuplet weights */ //@{ /** * The singlet weight */ - double _sngWt; + double _sngWt; /** * The decuplet weight */ - double _decWt; + double _decWt; //@} /** * The table of hadron data */ HadronTable _table; /** * Enums so arrays can be statically allocated */ //@{ /** * Defines values for array sizes. L,J,N max values for excited mesons. */ - enum MesonMultiplets { Lmax = 3, Jmax = 4, Nmax = 4}; + enum MesonMultiplets { Lmax = 3, Jmax = 4, Nmax = 4}; //@} /** * Option for the construction of the tables - */ + */ unsigned int _topt; /** * Which particles to produce for debugging purposes */ unsigned int _trial; /** * @name A parameter used for determining when clusters are too light. * * This parameter is used for setting the lower threshold, \f$ t \f$ as * \f[ t' = t(1 + r B^1_{\rm lim}) \f] * where \f$ r \f$ is a random number [0,1]. */ //@{ double _limBottom; double _limCharm; double _limExotic; //@} /** * Option for the selection of hadrons below the pair threshold */ unsigned int belowThreshold_; }; } #endif /* HERWIG_HadronSelector_H */ diff --git a/Hadronization/Hw64Selector.h b/Hadronization/Hw64Selector.h --- a/Hadronization/Hw64Selector.h +++ b/Hadronization/Hw64Selector.h @@ -1,111 +1,123 @@ // -*- C++ -*- // // Hw64Selector.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_Hw64Selector_H #define HERWIG_Hw64Selector_H // // This is the declaration of the Hw64Selector class. // #include "HadronSelector.h" #include "Hw64Selector.fh" namespace Herwig { using namespace ThePEG; /** \ingroup hadronization * The Hw64Selector class selects the hadrons produced in cluster decay using * the FORTRAN HERWIG variant of the cluster model. * * @see \ref Hw64SelectorInterfaces "The interfaces" * defined for Hw64Selector. */ class Hw64Selector: public HadronSelector { public: /** * The default constructor. */ Hw64Selector() : HadronSelector(0) {} /** * Method to return a pair of hadrons given the PDG codes of * two or three constituents * @param cluMass The mass of the cluster * @param par1 The particle pointer of the first constituent * @param par2 The particle pointer of the second constituent * @param par3 The particle pointer of the third constituent */ - virtual pair chooseHadronPair(const Energy cluMass,tcPDPtr par1, + virtual pair chooseHadronPair(const Energy cluMass,tcPDPtr par1, tcPDPtr par2,tcPDPtr par3 = PDPtr()) const ; + /** + * species picked + */ + long getQ(){ return _qPicked; } + 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 assignment operator is private and must never be called. * In fact, it should not even be implemented. */ Hw64Selector & operator=(const Hw64Selector &) = delete; + /** + * Internal variable to keep tracke of which quark species was pair + * produced + */ + mutable long _qPicked; + + }; } #endif /* HERWIG_Hw64Selector_H */ diff --git a/Hadronization/HwppSelector.cc b/Hadronization/HwppSelector.cc --- a/Hadronization/HwppSelector.cc +++ b/Hadronization/HwppSelector.cc @@ -1,245 +1,248 @@ // -*- 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; } void HwppSelector::persistentInput(PersistentIStream & is, int) { is >> _mode >> _enhanceSProb >> iunit(_m0Decay,GeV) >> _massMeasure; } 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 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.+pwtDIquark()) &&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 (quarktopick->id() == 3) { // Scaling strangeness enhancement 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); + // Update the picked quark species + // This is a cheeky work aroun to update the quark species + _qPicked = hadrons[ix].idQ->id(); 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,168 +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) {} /** * * 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(); + /** + * species picked + */ + long getQ(){ return _qPicked; } + 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.; + /** + * Internal variable to keep tracke of which quark species was pair + * produced + */ + mutable long _qPicked; + }; } #endif /* HERWIG_HwppSelector_H */