diff --git a/Hadronization/Cluster.cc b/Hadronization/Cluster.cc --- a/Hadronization/Cluster.cc +++ b/Hadronization/Cluster.cc @@ -1,271 +1,271 @@ // -*- C++ -*- // // Cluster.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the Cluster class. // #include "Cluster.h" #include #include #include #include "ClusterHadronizationHandler.h" using namespace Herwig; PPtr Cluster::clone() const { return new_ptr(*this); } PPtr Cluster::fullclone() const { return clone(); } ClassDescription Cluster::initCluster; Cluster::Cluster() : Particle(), _isAvailable(true), _hasReshuffled(false), _component(), _original(), _isBeamRemnant(), _isPerturbative(), _numComp(0), _id(0) {} void Cluster::persistentOutput(PersistentOStream & os) const { os << _isAvailable << _hasReshuffled << _component << _original << _isBeamRemnant << _isPerturbative << _numComp << _id; } void Cluster::persistentInput(PersistentIStream & is, int) { is >> _isAvailable >> _hasReshuffled >> _component >> _original >> _isBeamRemnant >> _isPerturbative >> _numComp >> _id; } Cluster::Cluster(tPPtr p1, tPPtr p2, tPPtr p3) : Particle(CurrentGenerator::current(). getParticleData(long(ParticleID::Cluster))), _isAvailable(true), _hasReshuffled(false) { if(!dataPtr()) { cerr << "Cluster Particle Data not defined. Cannot complete Hadronization " << "without ParticleData for id " << ParticleID::Cluster << '\n'; } _component.push_back(new_ptr(Particle(*p1))); _component.push_back(new_ptr(Particle(*p2))); if(p3) _component.push_back(new_ptr(Particle(*p3))); _original.push_back(p1); _original.push_back(p2); if(p3) _original.push_back(p3); _isPerturbative.push_back(initPerturbative(p1)); _isPerturbative.push_back(initPerturbative(p2)); if(p3) _isPerturbative.push_back(initPerturbative(p3)); else _isPerturbative.push_back(false); for(int i = 0; i<3; i++) _isBeamRemnant.push_back(false); if(p3) { _numComp = 3; _id = 100*abs(p1->id()) + 10*abs(p2->id()) + abs(p3->id()); } else { _numComp = 2; int i1,i2; if(p2->id() > 10) { i1 = abs(p2->id()/100); i2 = abs(p1->id()); } else if(p1->id() > 10) { i1 = abs(p1->id()/100); i2 = abs(p2->id()); } else { i1 = abs(p1->id()); i2 = abs(p2->id()); } if(i1>i2) swap (i1,i2); _id = 10*i1+i2; } // calculate the momentum calculateP(); // calculate the position // Only in the case of two components we have a definition of cluster // position in terms of the two components. if ( _numComp != 2 ) { // take the average setVertex(LorentzPoint()); } else { setVertex(calculateX(_component[0],_component[1])); } } Cluster::Cluster(tcEventPDPtr x) : Particle(x), _isAvailable(false), _hasReshuffled(false), _component(), _original(), _isBeamRemnant(), _isPerturbative(), _numComp(0), _id(0) {} Cluster::Cluster(const Particle &x) : Particle(x), _isAvailable(false), _hasReshuffled(false), _component(), _original(), _isBeamRemnant(), _isPerturbative(), _numComp(0), _id(0) {} Energy Cluster::sumConstituentMasses() const { if(_numComp == 3) { return _component[0]->mass() + _component[1]->mass() + _component[2]->mass(); } else if(_numComp == 2) return _component[0]->mass() + _component[1]->mass(); else return ZERO; } void Cluster::calculateP() { Lorentz5Momentum m; - for(int i = 0; i<_numComp; i++) + for(unsigned int i = 0; i<_numComp; i++) m += _component[i]->momentum(); m.rescaleMass(); set5Momentum(m); } LorentzPoint Cluster::calculateX(tPPtr q1, tPPtr q2) { // Get the needed parameters. Energy2 vmin2 = ClusterHadronizationHandler::currentHandler()->minVirtuality2(); Length dmax = ClusterHadronizationHandler::currentHandler()->maxDisplacement(); // Get the positions and displacements of the two components (Lab frame). LorentzPoint pos1 = q1->vertex(); Lorentz5Momentum p1 = q1->momentum(); LorentzDistance displace1 = -log( UseRandom::rnd() ) * hbarc * p1 * (1 / sqrt(sqr(p1.m2() - p1.mass2()) + sqr(vmin2))); if ( abs( displace1.m() ) > dmax ) { displace1 *= dmax / abs( displace1.m() ); } LorentzPoint pos2 = q2->vertex(); Lorentz5Momentum p2 = q2->momentum(); LorentzDistance displace2 = -log( UseRandom::rnd() ) * hbarc * p2 * (1 / sqrt(sqr(p2.m2() - p2.mass2()) + sqr(vmin2))); if ( abs( displace2.m() ) > dmax ) { displace2 *= dmax / abs( displace2.m() ); } double s1 = 0.0, s2 = 0.0; Lorentz5Momentum pcl = p1 + p2; if ( abs( pcl.vect().dot( displace1.vect() ) ) > 1.0e-20*MeV*mm && abs( pcl.vect().dot( displace2.vect() ) ) > 1.0e-20*MeV*mm ) { // The displacement with the smallest projection along pcl.vect() // is scaled up such that both displacements have equal projections // along pcl.vect(). double ratio = ( abs( pcl.vect().dot( displace1.vect() ) ) / abs( pcl.vect().dot( displace2.vect() ) ) ); if ( pcl.vect().dot(displace1.vect()) * pcl.vect().dot(displace2.vect()) < 0.0*sqr(MeV*mm) ) { ratio *= -1; } if ( abs( ratio ) > 1.0 ) { displace2 *= ratio; } else { displace1 *= ratio; } // Now determine the s1 and s2 values. double s1minusS2 = ( pcl.vect().dot( pos2.vect() - pos1.vect() ) / pcl.vect().dot( displace1.vect() ) ); if ( s1minusS2 < 0 ) { s1 = 1.0; s2 = s1 - s1minusS2; } else if ( s1minusS2 > 0 ) { s2 = 1; s1 = s2 + s1minusS2; } } // Now, finally, determine the cluster position LorentzPoint position = 0.5 * (pos1 + pos2 + s1*displace1 + s2*displace2); // set the decay vertex of the two particles via the lifeLength q1->setLifeLength(position-q1->vertex()); q2->setLifeLength(position-q2->vertex()); // return the answer return position; } bool Cluster::isBeamCluster() const { - for(int i = 0; i<_numComp; i++) + for(unsigned int i = 0; i<_numComp; i++) if(_isBeamRemnant[i]) return true; return false; } void Cluster::isBeamCluster(tPPtr part) { - for(int i = 0; i<_numComp; i++) { + for(unsigned int i = 0; i<_numComp; i++) { if(_original[i] == part) { _isBeamRemnant[i] = true; break; } } } bool Cluster::isStatusFinal() const { int s = children().size(); for(unsigned int i = 0; iPDGName() == "Cluster") s--; return ( s > 0); } tPPtr Cluster::particle(int i) const { return (i < _numComp) ? _component[i] : PPtr(); } tPPtr Cluster::colParticle(bool anti) const { if ( _numComp != 2 ) return PPtr(); if ( _original[0]->hasColour(anti) ) return _original[0]; else if ( _original[1]->hasColour(anti) ) return _original[1]; else return PPtr(); } tPPtr Cluster::antiColParticle() const { return colParticle(true); } bool Cluster::isPerturbative(int i) const { return _isPerturbative[i]; } bool Cluster::isBeamRemnant(int i) const { return _isBeamRemnant[i]; } void Cluster::setBeamRemnant(int i, bool b) { if(i < _numComp) _isBeamRemnant[i] = b; } bool Cluster::initPerturbative(tPPtr p) { Energy mg = CurrentGenerator::current().getParticleData(ParticleID::g)->constituentMass(); return p->scale() > sqr(mg); } diff --git a/Hadronization/Cluster.h b/Hadronization/Cluster.h --- a/Hadronization/Cluster.h +++ b/Hadronization/Cluster.h @@ -1,341 +1,341 @@ // -*- C++ -*- // // Cluster.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_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). */ - int numComponents() const + 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; } 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 - int _numComp; //!< The number of constituents + 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,473 @@ // -*- C++ -*- // // ClusterDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the 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; DescribeClass -describeClusterDecayer("Herwig::ClusterDecayer",""); +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) {} 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; } void ClusterDecayer::persistentInput(PersistentIStream & is, int) { is >> _hadronsSelector >> _clDirLight >> _clDirBottom >> _clDirCharm >> _clDirExotic >> _clSmrLight >> _clSmrBottom >> _clSmrCharm >> _clSmrExotic >> _onshell >> _masstry; } 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); } 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); } } } 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/ClusterFinder.cc b/Hadronization/ClusterFinder.cc --- a/Hadronization/ClusterFinder.cc +++ b/Hadronization/ClusterFinder.cc @@ -1,485 +1,485 @@ // -*- C++ -*- // // ClusterFinder.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the ClusterFinder class. // #include "ClusterFinder.h" #include #include #include #include #include #include #include #include #include "CheckId.h" #include "Herwig/Utilities/EnumParticles.h" #include "Herwig/Utilities/Kinematics.h" #include "Cluster.h" #include using namespace Herwig; DescribeClass -describeClusterFinder("Herwig::ClusterFinder",""); +describeClusterFinder("Herwig::ClusterFinder","Herwig.so"); IBPtr ClusterFinder::clone() const { return new_ptr(*this); } IBPtr ClusterFinder::fullclone() const { return new_ptr(*this); } void ClusterFinder::persistentOutput(PersistentOStream & os) const { os << heavyDiquarks_ << diQuarkSelection_ << diQuarkOnShell_; } void ClusterFinder::persistentInput(PersistentIStream & is, int) { is >> heavyDiquarks_ >> diQuarkSelection_ >> diQuarkOnShell_; } void ClusterFinder::Init() { static ClassDocumentation documentation ("This class is responsible of finding clusters."); static Switch interfaceHeavyDiquarks ("HeavyDiquarks", "How to treat heavy quarks in baryon number violating clusters", &ClusterFinder::heavyDiquarks_, 2, false, false); static SwitchOption interfaceHeavyDiquarksDefault (interfaceHeavyDiquarks, "Allow", "No special treatment, allow both heavy and doubly heavy diquarks", 0); static SwitchOption interfaceHeavyDiquarksNoDoublyHeavy (interfaceHeavyDiquarks, "NoDoublyHeavy", "Avoid having diquarks with twoo heavy quarks", 1); static SwitchOption interfaceHeavyDiquarksNoHeavy (interfaceHeavyDiquarks, "NoHeavy", "Try and avoid diquarks contain c and b altogether", 2); static Switch interfaceDiQuarkSelection ("DiQuarkSelection", "Option controlling the selection of quarks to merge into a diquark in baryon-number violating clusters", &ClusterFinder::diQuarkSelection_, 1, false, false); static SwitchOption interfaceDiQuarkSelectionRandom (interfaceDiQuarkSelection, "Random", "Randomly pick a pair to combine", 0); static SwitchOption interfaceDiQuarkSelectionLowestMass (interfaceDiQuarkSelection, "LowestMass", "Combine the lowest mass pair", 1); static Switch interfaceDiQuarkOnShell ("DiQuarkOnShell", "Force the diquark produced in baryon-number violating clusters to be on-shell", &ClusterFinder::diQuarkOnShell_, false, false, false); static SwitchOption interfaceDiQuarkOnShellYes (interfaceDiQuarkOnShell, "Yes", "Force to be on-shell", true); static SwitchOption interfaceDiQuarkOnShellNo (interfaceDiQuarkOnShell, "No", "Leave off-shell", false); } ClusterVector ClusterFinder::formClusters(const PVector & partons) { set examinedSet; // colour particles already included in a cluster map > quarkQuark; // quark quark map > aQuarkQuark; // anti quark anti quark ParticleSet inputParticles(partons.begin(),partons.end()); ClusterVector clusters; // Loop over all current particles. for(PVector::const_iterator pit=partons.begin();pit!=partons.end();++pit){ // Skip to the next particle if it is not coloured or already examined. assert(*pit); assert((*pit)->dataPtr()); if(!(**pit).data().coloured() || examinedSet.find(*pit) != examinedSet.end()) { continue; } // We assume that a cluster is made of, at most, 3 constituents although // in most cases the number will be 2; however, for baryon violating decays // (for example in Susy model without R parity conservation) we can have 3 // constituents. In the latter case, a quark (antiquark) do not have an // anticolour (colour) partner as usual, but its colour line either stems // from a colour source, or ends in a colour sink. In the case of double // baryon violating decays, but with overall baryon conservation // ( for instance: // tilde_u_R -> dbar_1 + dbar_2 // tilde_u_R_star -> d1 + d2 // where tilde_u_R and tilde_u_R_star are colour connected ) // a special treatment is needed, because first we have to process all // partons in the current step, and then for each left pair of quarks which // stem from a colour source we have to find the corresponding pair of // anti-quarks which ends in a colour sink and is connected with the // above colour source. These special pairs are kept into the maps: // spec/CluHadConfig.hialQuarkQuarkMap and specialAntiQuarkAntiQuarkMap. tParticleVector connected(3); int iElement = 0; connected[iElement++] = *pit; bool specialCase = false; if((*pit)->hasColour()) { tPPtr partner = (*pit)->colourLine()->getColouredParticle(partons.begin(), partons.end(), true); if(partner) { connected[iElement++]= partner; } // colour source : baryon-violating process else { if((*pit)->colourLine()->sourceNeighbours() != tColinePair()) { tColinePair sourcePair = (*pit)->colourLine()->sourceNeighbours(); tColinePtr intCL = tColinePtr(); for(int i = 0; i < 2; ++i) { tColinePtr pLine = i==0 ? sourcePair.first : sourcePair.second; int saveNumElements = iElement; for(tPVector::const_iterator cit = pLine->coloured().begin(); cit != pLine->coloured().end(); ++cit ) { ParticleSet::const_iterator cjt = inputParticles.find(*cit); if(cjt!=inputParticles.end()) connected[iElement++]= (*cit); } if(iElement == saveNumElements) intCL = pLine; } if(intCL && iElement == 2) { specialCase = true; pair qp=pair(connected[0],connected[1]); quarkQuark.insert(pair >(intCL,qp)); } else if(iElement != 3) { throw Exception() << "Colour connections fail in the hadronization for " << **pit << "in ClusterFinder::formClusters" << " for a coloured particle." << " Failed to find particles from a source" << Exception::runerror; } } else { throw Exception() << "Colour connections fail in the hadronization for " << **pit << "in ClusterFinder::formClusters for" << " a coloured particle" << Exception::runerror; } } } if((*pit)->hasAntiColour()) { tPPtr partner = (*pit)->antiColourLine()->getColouredParticle(partons.begin(), partons.end(), false); if(partner) { connected[iElement++]=partner; } // colour sink : baryon-violating process else { if((*pit)->antiColourLine()->sinkNeighbours() != tColinePair()) { tColinePair sinkPair = (*pit)->antiColourLine()->sinkNeighbours(); tColinePtr intCL = tColinePtr(); for(int i = 0; i < 2; ++i) { tColinePtr pLine = i==0 ? sinkPair.first : sinkPair.second; int saveNumElements = iElement; for(tPVector::const_iterator cit = pLine->antiColoured().begin(); cit != pLine->antiColoured().end(); ++cit ) { ParticleSet::const_iterator cjt = inputParticles.find(*cit); if(cjt!=inputParticles.end()) connected[iElement++]= (*cit); } if(iElement == saveNumElements) intCL = pLine; } if(intCL && iElement == 2) { specialCase = true; pair aqp=pair(connected[0],connected[1]); aQuarkQuark.insert(pair >(intCL,aqp)); } else if( iElement !=3) { throw Exception() << "Colour connections fail in the hadronization for " << **pit << "in ClusterFinder::formClusters for" << " an anti-coloured particle." << " Failed to find particles from a sink" << Exception::runerror; } } else { throw Exception() << "Colour connections fail in the hadronization for " << **pit << "in ClusterFinder::formClusters for" << " an anti-coloured particle" << Exception::runerror; } } } if(!specialCase) { // Tag the components of the found cluster as already examined. for (int i=0; iaddChild(cluPtr); connected[1]->addChild(cluPtr); if(connected[2]) connected[2]->addChild(cluPtr); clusters.push_back(cluPtr); // Check if any of the components is a beam remnant, and if this // is the case then inform the cluster. // this will only work for baryon collisions for (int i=0; iid()==ParticleID::Remnant) { fromRemnant = true; break; } parent = parent->parents().empty() ? tPPtr() : parent->parents()[0]; } if(fromRemnant&&DiquarkMatcher::Check(connected[i]->id())) cluPtr->isBeamCluster(connected[i]); } } } // Treat now the special cases, if any. The idea is to find for each pair // of quarks coming from a common colour source the corresponding pair of // antiquarks coming from a common colour sink, connected to the above // colour source via the same colour line. Then, randomly couple one of // the two quarks with one of the two antiquarks, and do the same with the // quark and antiquark left. for(map >::const_iterator cit = quarkQuark.begin(); cit != quarkQuark.end(); ++cit ) { tColinePtr coline = cit->first; pair quarkPair = cit->second; if(aQuarkQuark.find( coline ) != aQuarkQuark.end()) { pair antiQuarkPair = aQuarkQuark.find(coline)->second; ClusterPtr cluPtr1, cluPtr2; if ( UseRandom::rndbool() ) { cluPtr1 = new_ptr(Cluster(quarkPair.first , antiQuarkPair.first)); cluPtr2 = new_ptr(Cluster(quarkPair.second , antiQuarkPair.second)); quarkPair.first->addChild(cluPtr1); antiQuarkPair.first->addChild(cluPtr1); quarkPair.second->addChild(cluPtr2); antiQuarkPair.second->addChild(cluPtr2); } else { cluPtr1 = new_ptr(Cluster(quarkPair.first , antiQuarkPair.second)); cluPtr2 = new_ptr(Cluster(quarkPair.second , antiQuarkPair.first)); quarkPair.second->addChild(cluPtr2); antiQuarkPair.first->addChild(cluPtr2); quarkPair.first->addChild(cluPtr1); antiQuarkPair.second->addChild(cluPtr1); } clusters.push_back(cluPtr1); clusters.push_back(cluPtr2); } else { throw Exception() << "ClusterFinder::formClusters : " << "***Skip event: unable to match pairs in " << "Baryon-violating processes***" << Exception::eventerror; } } return clusters; } namespace { bool PartOrdering(tPPtr p1,tPPtr p2) { return abs(p1->id())id()); } } void ClusterFinder::reduceToTwoComponents(ClusterVector & clusters) { // In order to preserve all of the information, we do not modify the // directly the 3-component clusters, but instead we define new clusters, // which are related to the original ones by a child-parent relationship, // by considering two randomly chosen components as a diquark (or anti-diquark). // These new clusters are first added to the vector vecNewRedefinedCluPtr, // and at the end, when all input clusters have been examined, the elements of // this vector will be copied in collecCluPtr (the reason is that it is not // allowed to modify a STL container while iterating over it). vector redefinedClusters; for(ClusterVector::iterator cluIter = clusters.begin() ; cluIter != clusters.end() ; ++cluIter) { tParticleVector vec; if ( (*cluIter)->numComponents() != 3 || ! (*cluIter)->isAvailable() ) continue; tPPtr other; - for(int i = 0; i<(*cluIter)->numComponents(); i++) { + for(unsigned int i = 0; i<(*cluIter)->numComponents(); i++) { tPPtr part = (*cluIter)->particle(i); if(!DiquarkMatcher::Check(*(part->dataPtr()))) vec.push_back(part); else other = part; } if(vec.size()<2) { throw Exception() << "Could not make a diquark for a baryonic cluster decay from " << (*cluIter)->particle(0)->PDGName() << " " << (*cluIter)->particle(1)->PDGName() << " " << (*cluIter)->particle(2)->PDGName() << " " << " in ClusterFinder::reduceToTwoComponents()." << Exception::eventerror; } // order the vector so heaviest at the end std::sort(vec.begin(),vec.end(),PartOrdering); // Special treatment of heavy quarks // avoid doubly heavy diquarks if(heavyDiquarks_>=1 && vec.size()>2 && abs(vec[1]->id())>3 && abs(vec[0]->id())<=3) { if(UseRandom::rndbool()) swap(vec[1],vec[2]); other = vec[2]; vec.pop_back(); } // avoid singly heavy diquarks if(heavyDiquarks_==2 && vec.size()>2 && abs(vec[2]->id())>3 && abs(vec[1]->id())<=3) { other = vec[2]; vec.pop_back(); } // if there's a choice pick the pair to make a diquark from if(vec.size()>2) { unsigned int ichoice(0); // random choice if(diQuarkSelection_==0) { ichoice = UseRandom::rnd3(1.0, 1.0, 1.0); } // pick the lightest quark pair else if(diQuarkSelection_==1) { Energy m12 = (vec[0]->momentum()+vec[1]->momentum()).m(); Energy m13 = (vec[0]->momentum()+vec[2]->momentum()).m(); Energy m23 = (vec[1]->momentum()+vec[2]->momentum()).m(); if (m13<=m12&&m13<=m23) ichoice = 2; else if(m23<=m12&&m23<=m13) ichoice = 1; } else assert(false); // make the swaps so select pair first switch (ichoice) { case 0: break; case 1: swap(vec[2],vec[0]); break; case 2: swap(vec[2],vec[1]); break; } } // set up tcPDPtr temp1 = vec[0]->dataPtr(); tcPDPtr temp2 = vec[1]->dataPtr(); if(!other) other = vec[2]; tcPDPtr dataDiquark = CheckId::makeDiquark(temp1,temp2); if(!dataDiquark) throw Exception() << "Could not make a diquark from" << temp1->PDGName() << " and " << temp2->PDGName() << " in ClusterFinder::reduceToTwoComponents()" << Exception::eventerror; // Create the new cluster (with two components) and assign to it the same // momentum and position of the original (with three components) one. // Furthermore, assign to the diquark component a momentum given by the // sum of the two original components from which has been formed; for the // position, we are assuming, very simply, that the diquark position is // the average positions of the two original components. // Notice that the mass (5-th component of the 5-momentum) of the diquark // is set by hand to the constituent mass of the diquark (which is equal // to the sum of the constituent masses of the two quarks which form the // diquark) because the sum of 5-component vectors do add only the "normal" // 4-components, not the 5-th one. After that, the 5-momentum of the diquark // is in an inconsistent state, because the mass (5-th component) is not // equal to the invariant mass obtained from the 4-momemtum. This is not // unique to this kind of component (all perturbative components are in // a similar situation), but it is not harmful. // construct the diquark PPtr diquark = dataDiquark->produceParticle(); vec[0]->addChild(diquark); vec[1]->addChild(diquark); diquark->set5Momentum(Lorentz5Momentum(vec[0]->momentum() + vec[1]->momentum(), dataDiquark->constituentMass())); // use the same method as for cluster to determine the diquark position diquark->setVertex(Cluster::calculateX(vec[0],vec[1])); // put on-shell if required if(diQuarkOnShell_) { Lorentz5Momentum psum = diquark->momentum()+other->momentum(); psum.rescaleMass(); Boost boost = psum.boostVector(); Lorentz5Momentum pother = other->momentum(); Lorentz5Momentum pdiquark = diquark->momentum(); pother.boost(-boost); pdiquark.boost(-boost); Energy pcm = Kinematics::pstarTwoBodyDecay(psum.mass(), other->dataPtr()->constituentMass(), diquark->dataPtr()->constituentMass()); if(pcm>ZERO) { double fact = pcm/pother.vect().mag(); pother *= fact; pdiquark *= fact; pother .setMass(other->dataPtr()->constituentMass()); pdiquark.setMass(dataDiquark ->constituentMass()); pother .rescaleEnergy(); pdiquark.rescaleEnergy(); pother .boost(boost); pdiquark.boost(boost); other->set5Momentum(pother); diquark->set5Momentum(pdiquark); } } // make the new cluster ClusterPtr nclus = new_ptr(Cluster(other,diquark)); //vec[0]->addChild(nclus); //diquark->addChild(nclus); // Set the parent/children relationship between the original cluster // (the one with three components) with the new one (the one with two components) // and add the latter to the vector of new redefined clusters. (*cluIter)->addChild(nclus); redefinedClusters.push_back(nclus); } // Add to collecCluPtr all of the redefined new clusters (indeed the // pointers to them are added) contained in vecNewRedefinedCluPtr. /// \todo why do we keep the original of the redefined clusters? for (tClusterVector::const_iterator it = redefinedClusters.begin(); it != redefinedClusters.end(); ++it) { clusters.push_back(*it); } } diff --git a/Hadronization/ClusterFissioner.cc b/Hadronization/ClusterFissioner.cc --- a/Hadronization/ClusterFissioner.cc +++ b/Hadronization/ClusterFissioner.cc @@ -1,1121 +1,1121 @@ // -*- C++ -*- // // ClusterFissioner.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // Thisk is the implementation of the non-inlined, non-templated member // functions of the ClusterFissioner class. // #include "ClusterFissioner.h" #include #include #include #include #include #include #include #include "Herwig/Utilities/Kinematics.h" #include "CheckId.h" #include "Cluster.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include using namespace Herwig; DescribeClass -describeClusterFissioner("Herwig::ClusterFissioner",""); +describeClusterFissioner("Herwig::ClusterFissioner","Herwig.so"); ClusterFissioner::ClusterFissioner() : _clMaxLight(3.35*GeV), _clMaxBottom(3.35*GeV), _clMaxCharm(3.35*GeV), _clMaxExotic(3.35*GeV), _clPowLight(2.0), _clPowBottom(2.0), _clPowCharm(2.0), _clPowExotic(2.0), _pSplitLight(1.0), _pSplitBottom(1.0), _pSplitCharm(1.0), _pSplitExotic(1.0), _fissionCluster(0), _fissionPwtUquark(1), _fissionPwtDquark(1), _fissionPwtSquark(0.5), _btClM(1.0*GeV), _iopRem(1), _kappa(1.0e15*GeV/meter), _enhanceSProb(0), _m0Fission(2.*GeV), _massMeasure(0) {} IBPtr ClusterFissioner::clone() const { return new_ptr(*this); } IBPtr ClusterFissioner::fullclone() const { return new_ptr(*this); } void ClusterFissioner::persistentOutput(PersistentOStream & os) const { os << _hadronsSelector << ounit(_clMaxLight,GeV) << ounit(_clMaxBottom,GeV) << ounit(_clMaxCharm,GeV) << ounit(_clMaxExotic,GeV) << _clPowLight << _clPowBottom << _clPowCharm << _clPowExotic << _pSplitLight << _pSplitBottom << _pSplitCharm << _pSplitExotic << _fissionCluster << _fissionPwtUquark << _fissionPwtDquark << _fissionPwtSquark << ounit(_btClM,GeV) << _iopRem << ounit(_kappa, GeV/meter) << _enhanceSProb << ounit(_m0Fission,GeV) << _massMeasure; } void ClusterFissioner::persistentInput(PersistentIStream & is, int) { is >> _hadronsSelector >> iunit(_clMaxLight,GeV) >> iunit(_clMaxBottom,GeV) >> iunit(_clMaxCharm,GeV) >> iunit(_clMaxExotic,GeV) >> _clPowLight >> _clPowBottom >> _clPowCharm >> _clPowExotic >> _pSplitLight >> _pSplitBottom >> _pSplitCharm >> _pSplitExotic >> _fissionCluster >> _fissionPwtUquark >> _fissionPwtDquark >> _fissionPwtSquark >> iunit(_btClM,GeV) >> _iopRem >> iunit(_kappa, GeV/meter) >> _enhanceSProb >> iunit(_m0Fission,GeV) >> _massMeasure; } void ClusterFissioner::Init() { static ClassDocumentation documentation ("Class responsibles for chopping up the clusters"); static Reference interfaceHadronSelector("HadronSelector", "A reference to the HadronSelector object", &Herwig::ClusterFissioner::_hadronsSelector, false, false, true, false); // ClMax for light, Bottom, Charm and exotic (e.g. Susy) quarks static Parameter interfaceClMaxLight ("ClMaxLight","cluster max mass for light quarks (unit [GeV])", &ClusterFissioner::_clMaxLight, GeV, 3.35*GeV, ZERO, 10.0*GeV, false,false,false); static Parameter interfaceClMaxBottom ("ClMaxBottom","cluster max mass for b quarks (unit [GeV])", &ClusterFissioner::_clMaxBottom, GeV, 3.35*GeV, ZERO, 10.0*GeV, false,false,false); static Parameter interfaceClMaxCharm ("ClMaxCharm","cluster max mass for c quarks (unit [GeV])", &ClusterFissioner::_clMaxCharm, GeV, 3.35*GeV, ZERO, 10.0*GeV, false,false,false); static Parameter interfaceClMaxExotic ("ClMaxExotic","cluster max mass for exotic quarks (unit [GeV])", &ClusterFissioner::_clMaxExotic, GeV, 3.35*GeV, ZERO, 10.0*GeV, false,false,false); // ClPow for light, Bottom, Charm and exotic (e.g. Susy) quarks static Parameter interfaceClPowLight ("ClPowLight","cluster mass exponent for light quarks", &ClusterFissioner::_clPowLight, 0, 2.0, 0.0, 10.0,false,false,false); static Parameter interfaceClPowBottom ("ClPowBottom","cluster mass exponent for b quarks", &ClusterFissioner::_clPowBottom, 0, 2.0, 0.0, 10.0,false,false,false); static Parameter interfaceClPowCharm ("ClPowCharm","cluster mass exponent for c quarks", &ClusterFissioner::_clPowCharm, 0, 2.0, 0.0, 10.0,false,false,false); static Parameter interfaceClPowExotic ("ClPowExotic","cluster mass exponent for exotic quarks", &ClusterFissioner::_clPowExotic, 0, 2.0, 0.0, 10.0,false,false,false); // PSplit for light, Bottom, Charm and exotic (e.g. Susy) quarks static Parameter interfacePSplitLight ("PSplitLight","cluster mass splitting param for light quarks", &ClusterFissioner::_pSplitLight, 0, 1.0, 0.0, 10.0,false,false,false); static Parameter interfacePSplitBottom ("PSplitBottom","cluster mass splitting param for b quarks", &ClusterFissioner::_pSplitBottom, 0, 1.0, 0.0, 10.0,false,false,false); static Parameter interfacePSplitCharm ("PSplitCharm","cluster mass splitting param for c quarks", &ClusterFissioner::_pSplitCharm, 0, 1.0, 0.0, 10.0,false,false,false); static Parameter interfacePSplitExotic ("PSplitExotic","cluster mass splitting param for exotic quarks", &ClusterFissioner::_pSplitExotic, 0, 1.0, 0.0, 10.0,false,false,false); static Switch interfaceFission ("Fission", "Option for different Fission options", &ClusterFissioner::_fissionCluster, 1, false, false); static SwitchOption interfaceFissionDefault (interfaceFission, "default", "Normal cluster fission which depends on the hadron selector class.", 0); static SwitchOption interfaceFissionNew (interfaceFission, "new", "Alternative cluster fission which does not depend on the hadron selector class", 1); static Parameter interfaceFissionPwtUquark ("FissionPwtUquark", "Weight for fission in U quarks", &ClusterFissioner::_fissionPwtUquark, 1, 0.0, 1.0, false, false, Interface::limited); static Parameter interfaceFissionPwtDquark ("FissionPwtDquark", "Weight for fission in D quarks", &ClusterFissioner::_fissionPwtDquark, 1, 0.0, 1.0, false, false, Interface::limited); static Parameter interfaceFissionPwtSquark ("FissionPwtSquark", "Weight for fission in S quarks", &ClusterFissioner::_fissionPwtSquark, 0.5, 0.0, 1.0, false, false, Interface::limited); static Switch interfaceRemnantOption ("RemnantOption", "Option for the treatment of remnant clusters", &ClusterFissioner::_iopRem, 1, false, false); static SwitchOption interfaceRemnantOptionSoft (interfaceRemnantOption, "Soft", "Both clusters produced in the fission of the beam cluster" " are treated as soft clusters.", 0); static SwitchOption interfaceRemnantOptionHard (interfaceRemnantOption, "Hard", "Only the cluster containing the remnant is treated as a soft cluster.", 1); static SwitchOption interfaceRemnantOptionVeryHard (interfaceRemnantOption, "VeryHard", "Even remnant clusters are treated as hard, i.e. all clusters the same", 2); static Parameter interfaceBTCLM ("SoftClusterFactor", "Parameter for the mass spectrum of remnant clusters", &ClusterFissioner::_btClM, GeV, 1.*GeV, 0.1*GeV, 10.0*GeV, false, false, Interface::limited); static Parameter interfaceStringTension ("StringTension", "String tension used in vertex displacement calculation", &ClusterFissioner::_kappa, GeV/meter, 1.0e15*GeV/meter, ZERO, ZERO, false, false, Interface::lowerlim); static Switch interfaceEnhanceSProb ("EnhanceSProb", "Option for enhancing strangeness", &ClusterFissioner::_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", &ClusterFissioner::_massMeasure,0,false,false); static SwitchOption interfaceMassMeasureMass (interfaceMassMeasure, "Mass", "Mass Measure", 0); static SwitchOption interfaceMassMeasureLambda (interfaceMassMeasure, "Lambda", "Lambda Measure", 1); static Parameter interfaceFissionMassScale ("FissionMassScale", "Cluster fission mass scale", &ClusterFissioner::_m0Fission, GeV, 2.0*GeV, 0.1*GeV, 50.*GeV, false, false, Interface::limited); } tPVector ClusterFissioner::fission(ClusterVector & clusters, bool softUEisOn) { // return if no clusters if (clusters.empty()) return tPVector(); /***************** * Loop over the (input) collection of cluster pointers, and store in * the vector splitClusters all the clusters that need to be split * (these are beam clusters, if soft underlying event is off, and * heavy non-beam clusters). ********************/ stack splitClusters; for(ClusterVector::iterator it = clusters.begin() ; it != clusters.end() ; ++it) { /************** * Skip 3-component clusters that have been redefined (as 2-component * clusters) or not available clusters. The latter check is indeed * redundant now, but it is used for possible future extensions in which, * for some reasons, some of the clusters found by ClusterFinder are tagged * straight away as not available. **************/ if((*it)->isRedefined() || !(*it)->isAvailable()) continue; // if the cluster is a beam cluster add it to the vector of clusters // to be split or if it is heavy if((*it)->isBeamCluster() || isHeavy(*it)) splitClusters.push(*it); } tPVector finalhadrons; cut(splitClusters, clusters, finalhadrons, softUEisOn); return finalhadrons; } void ClusterFissioner::cut(stack & clusterStack, ClusterVector &clusters, tPVector & finalhadrons, bool softUEisOn) { /************************************************** * This method does the splitting of the cluster pointed by cluPtr * and "recursively" by all of its cluster children, if heavy. All of these * new children clusters are added (indeed the pointers to them) to the * collection of cluster pointers collecCluPtr. The method works as follows. * Initially the vector vecCluPtr contains just the input pointer to the * cluster to be split. Then it will be filled "recursively" by all * of the cluster's children that are heavy enough to require, in their turn, * to be split. In each loop, the last element of the vector vecCluPtr is * considered (only once because it is then removed from the vector). * This approach is conceptually recursive, but avoid the overhead of * a concrete recursive function. Furthermore it requires minimal changes * in the case that the fission of an heavy cluster could produce more * than two cluster children as assumed now. * * Draw the masses: for normal, non-beam clusters a power-like mass dist * is used, whereas for beam clusters a fast-decreasing exponential mass * dist is used instead (to avoid many iterative splitting which could * produce an unphysical large transverse energy from a supposed soft beam * remnant process). ****************************************/ // Here we recursively loop over clusters in the stack and cut them while (!clusterStack.empty()) { // take the last element of the vector ClusterPtr iCluster = clusterStack.top(); clusterStack.pop(); // split it cutType ct = iCluster->numComponents() == 2 ? cutTwo(iCluster, finalhadrons, softUEisOn) : cutThree(iCluster, finalhadrons, softUEisOn); // There are cases when we don't want to split, even if it fails mass test if(!ct.first.first || !ct.second.first) { // if an unsplit beam cluster leave if for the underlying event if(iCluster->isBeamCluster() && softUEisOn) iCluster->isAvailable(false); continue; } // check if clusters ClusterPtr one = dynamic_ptr_cast(ct.first.first); ClusterPtr two = dynamic_ptr_cast(ct.second.first); // is a beam cluster must be split into two clusters if(iCluster->isBeamCluster() && (!one||!two) && softUEisOn) { iCluster->isAvailable(false); continue; } // There should always be a intermediate quark(s) from the splitting assert(ct.first.second && ct.second.second); /// \todo sort out motherless quark pairs here. Watch out for 'quark in final state' errors iCluster->addChild(ct.first.first); // iCluster->addChild(ct.first.second); // ct.first.second->addChild(ct.first.first); iCluster->addChild(ct.second.first); // iCluster->addChild(ct.second.second); // ct.second.second->addChild(ct.second.first); // Sometimes the clusters decay C -> H + C' rather then C -> C' + C'' if(one) { clusters.push_back(one); if(one->isBeamCluster() && softUEisOn) one->isAvailable(false); if(isHeavy(one) && one->isAvailable()) clusterStack.push(one); } if(two) { clusters.push_back(two); if(two->isBeamCluster() && softUEisOn) two->isAvailable(false); if(isHeavy(two) && two->isAvailable()) clusterStack.push(two); } } } namespace { /** * Check if can't make a hadron from the partons */ bool cantMakeHadron(tcPPtr p1, tcPPtr p2) { return ! CheckId::canBeHadron(p1->dataPtr(), p2->dataPtr()); } /** * Check if can't make a diquark from the partons */ bool cantMakeDiQuark(tcPPtr p1, tcPPtr p2) { long id1 = p1->id(), id2 = p2->id(); return ! (QuarkMatcher::Check(id1) && QuarkMatcher::Check(id2) && id1*id2>0); } } ClusterFissioner::cutType ClusterFissioner::cutTwo(ClusterPtr & cluster, tPVector & finalhadrons, bool softUEisOn) { // need to make sure only 2-cpt clusters get here assert(cluster->numComponents() == 2); tPPtr ptrQ1 = cluster->particle(0); tPPtr ptrQ2 = cluster->particle(1); Energy Mc = cluster->mass(); assert(ptrQ1); assert(ptrQ2); // And check if those particles are from a beam remnant bool rem1 = cluster->isBeamRemnant(0); bool rem2 = cluster->isBeamRemnant(1); // workout which distribution to use bool soft1(false),soft2(false); switch (_iopRem) { case 0: soft1 = rem1 || rem2; soft2 = rem2 || rem1; break; case 1: soft1 = rem1; soft2 = rem2; break; } // Initialization for the exponential ("soft") mass distribution. static const int max_loop = 1000; int counter = 0; Energy Mc1 = ZERO, Mc2 = ZERO,m1=ZERO,m2=ZERO,m=ZERO; tcPDPtr toHadron1, toHadron2; PPtr newPtr1 = PPtr (); PPtr newPtr2 = PPtr (); bool succeeded = false; do { succeeded = false; ++counter; if (_enhanceSProb == 0){ drawNewFlavour(newPtr1,newPtr2); } else { Energy2 mass2 = clustermass(cluster); drawNewFlavourEnhanced(newPtr1,newPtr2,mass2); } // check for right ordering assert (ptrQ2); assert (newPtr2); assert (ptrQ2->dataPtr()); assert (newPtr2->dataPtr()); if(cantMakeHadron(ptrQ1, newPtr1) || cantMakeHadron(ptrQ2, newPtr2)) { swap(newPtr1, newPtr2); // check again if(cantMakeHadron(ptrQ1, newPtr1) || cantMakeHadron(ptrQ2, newPtr2)) { throw Exception() << "ClusterFissioner cannot split the cluster (" << ptrQ1->PDGName() << ' ' << ptrQ2->PDGName() << ") into hadrons.\n" << Exception::runerror; } } // Check that new clusters can produce particles and there is enough // phase space to choose the drawn flavour m1 = ptrQ1->data().constituentMass(); m2 = ptrQ2->data().constituentMass(); m = newPtr1->data().constituentMass(); // Do not split in the case there is no phase space available if(Mc < m1+m + m2+m) continue; // power for splitting double exp1=_pSplitLight; double exp2=_pSplitLight; if (CheckId::isExotic(ptrQ1->dataPtr())) exp1 = _pSplitExotic; else if(CheckId::hasBottom(ptrQ1->dataPtr()))exp1 = _pSplitBottom; else if(CheckId::hasCharm(ptrQ1->dataPtr())) exp1 = _pSplitCharm; if (CheckId::isExotic(ptrQ2->dataPtr())) exp2 = _pSplitExotic; else if(CheckId::hasBottom(ptrQ2->dataPtr())) exp2 = _pSplitBottom; else if(CheckId::hasCharm(ptrQ2->dataPtr())) exp2 = _pSplitCharm; // If, during the drawing of candidate masses, too many attempts fail // (because the phase space available is tiny) /// \todo run separate loop here? Mc1 = drawChildMass(Mc,m1,m2,m,exp1,soft1); Mc2 = drawChildMass(Mc,m2,m1,m,exp2,soft2); if(Mc1 < m1+m || Mc2 < m+m2 || Mc1+Mc2 > Mc) continue; /************************** * New (not present in Fortran Herwig): * check whether the fragment masses Mc1 and Mc2 are above the * threshold for the production of the lightest pair of hadrons with the * right flavours. If not, then set by hand the mass to the lightest * single hadron with the right flavours, in order to solve correctly * the kinematics, and (later in this method) create directly such hadron * and add it to the children hadrons of the cluster that undergoes the * fission (i.e. the one pointed by iCluPtr). Notice that in this special * case, the heavy cluster that undergoes the fission has one single * cluster child and one single hadron child. We prefer this approach, * rather than to create a light cluster, with the mass set equal to * the lightest hadron, and let then the class LightClusterDecayer to do * the job to decay it to that single hadron, for two reasons: * First, because the sum of the masses of the two constituents can be, * in this case, greater than the mass of that hadron, hence it would * be impossible to solve the kinematics for such two components, and * therefore we would have a cluster whose components are undefined. * Second, the algorithm is faster, because it avoids the reshuffling * procedure that would be necessary if we used LightClusterDecayer * to decay the light cluster to the lightest hadron. ****************************/ toHadron1 = _hadronsSelector->chooseSingleHadron(ptrQ1->dataPtr(), newPtr1->dataPtr(),Mc1); if(toHadron1) Mc1 = toHadron1->mass(); toHadron2 = _hadronsSelector->chooseSingleHadron(ptrQ2->dataPtr(), newPtr2->dataPtr(),Mc2); if(toHadron2) Mc2 = toHadron2->mass(); // if a beam cluster not allowed to decay to hadrons if(cluster->isBeamCluster() && (toHadron1||toHadron2) && softUEisOn) continue; // Check if the decay kinematics is still possible: if not then // force the one-hadron decay for the other cluster as well. if(Mc1 + Mc2 > Mc) { if(!toHadron1) { toHadron1 = _hadronsSelector->chooseSingleHadron(ptrQ1->dataPtr(), newPtr1->dataPtr(),Mc-Mc2); if(toHadron1) Mc1 = toHadron1->mass(); } else if(!toHadron2) { toHadron2 = _hadronsSelector->chooseSingleHadron(ptrQ2->dataPtr(), newPtr2->dataPtr(),Mc-Mc1); if(toHadron2) Mc2 = toHadron2->mass(); } } succeeded = (Mc >= Mc1+Mc2); } while (!succeeded && counter < max_loop); if(counter >= max_loop) { static const PPtr null = PPtr(); return cutType(PPair(null,null),PPair(null,null)); } // Determined the (5-components) momenta (all in the LAB frame) Lorentz5Momentum pClu = cluster->momentum(); // known Lorentz5Momentum p0Q1 = ptrQ1->momentum(); // known (mom Q1 before fission) Lorentz5Momentum pClu1, pClu2, pQ1, pQone, pQtwo, pQ2; //unknown pClu1.setMass(Mc1); pClu2.setMass(Mc2); pQ1.setMass(m1); pQ2.setMass(m2); pQone.setMass(m); pQtwo.setMass(m); calculateKinematics(pClu,p0Q1,toHadron1,toHadron2, pClu1,pClu2,pQ1,pQone,pQtwo,pQ2); // out /****************** * The previous methods have determined the kinematics and positions * of C -> C1 + C2. * In the case that one of the two product is light, that means either * decayOneHadronClu1 or decayOneHadronClu2 is true, then the momenta * of the components of that light product have not been determined, * and a (light) cluster will not be created: the heavy father cluster * decays, in this case, into a single (not-light) cluster and a * single hadron. In the other, "normal", cases the father cluster * decays into two clusters, each of which has well defined components. * Notice that, in the case of components which point to particles, the * momenta of the components is properly set to the new values, whereas * we do not change the momenta of the pointed particles, because we * want to keep all of the information (that is the new momentum of a * component after the splitting, which is contained in the _momentum * member of the Component class, and the (old) momentum of that component * before the splitting, which is contained in the momentum of the * pointed particle). Please not make confusion of this only apparent * inconsistency! ********************/ LorentzPoint posC,pos1,pos2; posC = cluster->vertex(); calculatePositions(pClu, posC, pClu1, pClu2, pos1, pos2); cutType rval; if(toHadron1) { rval.first = produceHadron(toHadron1, newPtr1, pClu1, pos1); finalhadrons.push_back(rval.first.first); } else { rval.first = produceCluster(ptrQ1, newPtr1, pClu1, pos1, pQ1, pQone, rem1); } if(toHadron2) { rval.second = produceHadron(toHadron2, newPtr2, pClu2, pos2); finalhadrons.push_back(rval.second.first); } else { rval.second = produceCluster(ptrQ2, newPtr2, pClu2, pos2, pQ2, pQtwo, rem2); } return rval; } ClusterFissioner::cutType ClusterFissioner::cutThree(ClusterPtr & cluster, tPVector & finalhadrons, bool softUEisOn) { // need to make sure only 3-cpt clusters get here assert(cluster->numComponents() == 3); // extract quarks tPPtr ptrQ[3] = {cluster->particle(0),cluster->particle(1),cluster->particle(2)}; assert( ptrQ[0] && ptrQ[1] && ptrQ[2] ); // find maximum mass pair Energy mmax(ZERO); Lorentz5Momentum pDiQuark; int iq1(-1),iq2(-1); Lorentz5Momentum psum; for(int q1=0;q1<3;++q1) { psum+= ptrQ[q1]->momentum(); for(int q2=q1+1;q2<3;++q2) { Lorentz5Momentum ptest = ptrQ[q1]->momentum()+ptrQ[q2]->momentum(); ptest.rescaleMass(); Energy mass = ptest.m(); if(mass>mmax) { mmax = mass; pDiQuark = ptest; iq1 = q1; iq2 = q2; } } } // and the spectators int iother(-1); for(int ix=0;ix<3;++ix) if(ix!=iq1&&ix!=iq2) iother=ix; assert(iq1>=0&&iq2>=0&&iother>=0); // And check if those particles are from a beam remnant bool rem1 = cluster->isBeamRemnant(iq1); bool rem2 = cluster->isBeamRemnant(iq2); // workout which distribution to use bool soft1(false),soft2(false); switch (_iopRem) { case 0: soft1 = rem1 || rem2; soft2 = rem2 || rem1; break; case 1: soft1 = rem1; soft2 = rem2; break; } // Initialization for the exponential ("soft") mass distribution. static const int max_loop = 1000; int counter = 0; Energy Mc1 = ZERO, Mc2 = ZERO, m1=ZERO, m2=ZERO, m=ZERO; tcPDPtr toHadron; bool toDiQuark(false); PPtr newPtr1 = PPtr(),newPtr2 = PPtr(); PDPtr diquark; bool succeeded = false; do { succeeded = false; ++counter; if (_enhanceSProb == 0) { drawNewFlavour(newPtr1,newPtr2); } else { Energy2 mass2 = clustermass(cluster); drawNewFlavourEnhanced(newPtr1,newPtr2, mass2); } // randomly pick which will be (anti)diquark and which a mesonic cluster if(UseRandom::rndbool()) { swap(iq1,iq2); swap(rem1,rem2); } // check first order if(cantMakeHadron(ptrQ[iq1], newPtr1) || cantMakeDiQuark(ptrQ[iq2], newPtr2)) { swap(newPtr1,newPtr2); } // check again if(cantMakeHadron(ptrQ[iq1], newPtr1) || cantMakeDiQuark(ptrQ[iq2], newPtr2)) { throw Exception() << "ClusterFissioner cannot split the cluster (" << ptrQ[iq1]->PDGName() << ' ' << ptrQ[iq2]->PDGName() << ") into a hadron and diquark.\n" << Exception::runerror; } // Check that new clusters can produce particles and there is enough // phase space to choose the drawn flavour m1 = ptrQ[iq1]->data().constituentMass(); m2 = ptrQ[iq2]->data().constituentMass(); m = newPtr1->data().constituentMass(); // Do not split in the case there is no phase space available if(mmax < m1+m + m2+m) continue; // power for splitting double exp1(_pSplitLight),exp2(_pSplitLight); if (CheckId::isExotic (ptrQ[iq1]->dataPtr())) exp1 = _pSplitExotic; else if(CheckId::hasBottom(ptrQ[iq1]->dataPtr())) exp1 = _pSplitBottom; else if(CheckId::hasCharm (ptrQ[iq1]->dataPtr())) exp1 = _pSplitCharm; if (CheckId::isExotic (ptrQ[iq2]->dataPtr())) exp2 = _pSplitExotic; else if(CheckId::hasBottom(ptrQ[iq2]->dataPtr())) exp2 = _pSplitBottom; else if(CheckId::hasCharm (ptrQ[iq2]->dataPtr())) exp2 = _pSplitCharm; // If, during the drawing of candidate masses, too many attempts fail // (because the phase space available is tiny) /// \todo run separate loop here? Mc1 = drawChildMass(mmax,m1,m2,m,exp1,soft1); Mc2 = drawChildMass(mmax,m2,m1,m,exp2,soft2); if(Mc1 < m1+m || Mc2 < m+m2 || Mc1+Mc2 > mmax) continue; // check if need to force meson clster to hadron toHadron = _hadronsSelector->chooseSingleHadron(ptrQ[iq1]->dataPtr(), newPtr1->dataPtr(),Mc1); if(toHadron) Mc1 = toHadron->mass(); // check if need to force diquark cluster to be on-shell toDiQuark = false; diquark = CheckId::makeDiquark(ptrQ[iq2]->dataPtr(), newPtr2->dataPtr()); if(Mc2 < diquark->constituentMass()) { Mc2 = diquark->constituentMass(); toDiQuark = true; } // if a beam cluster not allowed to decay to hadrons if(cluster->isBeamCluster() && toHadron && softUEisOn) continue; // Check if the decay kinematics is still possible: if not then // force the one-hadron decay for the other cluster as well. if(Mc1 + Mc2 > mmax) { if(!toHadron) { toHadron = _hadronsSelector->chooseSingleHadron(ptrQ[iq1]->dataPtr(), newPtr1->dataPtr(),mmax-Mc2); if(toHadron) Mc1 = toHadron->mass(); } else if(!toDiQuark) { Mc2 = _hadronsSelector->massLightestHadron(ptrQ[iq2]->dataPtr(), newPtr2->dataPtr()); toDiQuark = true; } } succeeded = (mmax >= Mc1+Mc2); } while (!succeeded && counter < max_loop); // check no of tries if(counter >= max_loop) return cutType(); // Determine the (5-components) momenta (all in the LAB frame) Lorentz5Momentum p0Q1 = ptrQ[iq1]->momentum(); // to be determined Lorentz5Momentum pClu1(Mc1), pClu2(Mc2), pQ1(m1), pQone(m), pQtwo(m), pQ2(m2); calculateKinematics(pDiQuark,p0Q1,toHadron,toDiQuark, pClu1,pClu2,pQ1,pQone,pQtwo,pQ2); // positions of the new clusters LorentzPoint pos1,pos2; Lorentz5Momentum pBaryon = pClu2+ptrQ[iother]->momentum(); calculatePositions(cluster->momentum(), cluster->vertex(), pClu1, pBaryon, pos1, pos2); // first the mesonic cluster/meson cutType rval; if(toHadron) { rval.first = produceHadron(toHadron, newPtr1, pClu1, pos1); finalhadrons.push_back(rval.first.first); } else { rval.first = produceCluster(ptrQ[iq1], newPtr1, pClu1, pos1, pQ1, pQone, rem1); } if(toDiQuark) { rem2 |= cluster->isBeamRemnant(iother); PPtr newDiQuark = diquark->produceParticle(pClu2); rval.second = produceCluster(newDiQuark, ptrQ[iother], pBaryon, pos2, pClu2, ptrQ[iother]->momentum(), rem2); } else { rval.second = produceCluster(ptrQ[iq2], newPtr2, pBaryon, pos2, pQ2, pQtwo, rem2, ptrQ[iother],cluster->isBeamRemnant(iother)); } cluster->isAvailable(false); return rval; } ClusterFissioner::PPair ClusterFissioner::produceHadron(tcPDPtr hadron, tPPtr newPtr, const Lorentz5Momentum &a, const LorentzPoint &b) const { PPair rval; if(hadron->coloured()) { rval.first = (_hadronsSelector->lightestHadron(hadron,newPtr->dataPtr()))->produceParticle(); } else rval.first = hadron->produceParticle(); rval.second = newPtr; rval.first->set5Momentum(a); rval.first->setVertex(b); return rval; } ClusterFissioner::PPair ClusterFissioner::produceCluster(tPPtr ptrQ, tPPtr newPtr, const Lorentz5Momentum & a, const LorentzPoint & b, const Lorentz5Momentum & c, const Lorentz5Momentum & d, bool isRem, tPPtr spect, bool remSpect) const { PPair rval; rval.second = newPtr; ClusterPtr cluster = !spect ? new_ptr(Cluster(ptrQ,rval.second)) : new_ptr(Cluster(ptrQ,rval.second,spect)); rval.first = cluster; cluster->set5Momentum(a); cluster->setVertex(b); assert(cluster->particle(0)->id() == ptrQ->id()); cluster->particle(0)->set5Momentum(c); cluster->particle(1)->set5Momentum(d); cluster->setBeamRemnant(0,isRem); if(remSpect) cluster->setBeamRemnant(2,remSpect); return rval; } void ClusterFissioner::drawNewFlavour(PPtr& newPtrPos,PPtr& newPtrNeg) const { // Flavour is assumed to be only u, d, s, with weights // (which are not normalized probabilities) given // by the same weights as used in HadronsSelector for // the decay of clusters into two hadrons. double prob_d; double prob_u; double prob_s; switch(_fissionCluster){ case 0: prob_d = _hadronsSelector->pwtDquark(); prob_u = _hadronsSelector->pwtUquark(); prob_s = _hadronsSelector->pwtSquark(); break; case 1: prob_d = _fissionPwtDquark; prob_u = _fissionPwtUquark; prob_s = _fissionPwtSquark; break; default : assert(false); } int choice = UseRandom::rnd3(prob_u, prob_d, prob_s); long idNew = 0; switch (choice) { case 0: idNew = ThePEG::ParticleID::u; break; case 1: idNew = ThePEG::ParticleID::d; break; case 2: idNew = ThePEG::ParticleID::s; break; } newPtrPos = getParticle(idNew); newPtrNeg = getParticle(-idNew); assert (newPtrPos); assert(newPtrNeg); assert (newPtrPos->dataPtr()); assert(newPtrNeg->dataPtr()); } void ClusterFissioner::drawNewFlavourEnhanced(PPtr& newPtrPos,PPtr& newPtrNeg, Energy2 mass2) const { // Flavour is assumed to be only u, d, s, with weights // (which are not normalized probabilities) given // by the same weights as used in HadronsSelector for // the decay of clusters into two hadrons. double prob_d; double prob_u; double prob_s = 0.; double scale = abs(double(sqr(_m0Fission)/mass2)); // Choose which splitting weights you wish to use switch(_fissionCluster){ // 0: ClusterFissioner and ClusterDecayer use the same weights case 0: prob_d = _hadronsSelector->pwtDquark(); prob_u = _hadronsSelector->pwtUquark(); /* Strangeness enhancement: Case 1: probability scaling Case 2: Exponential scaling */ if (_enhanceSProb == 1) prob_s = (_maxScale < scale) ? 0. : pow(_hadronsSelector->pwtSquark(),scale); else if (_enhanceSProb == 2) prob_s = (_maxScale < scale) ? 0. : exp(-scale); break; /* 1: ClusterFissioner uses its own unique set of weights, i.e. decoupled from ClusterDecayer */ case 1: prob_d = _fissionPwtDquark; prob_u = _fissionPwtUquark; if (_enhanceSProb == 1) prob_s = (_maxScale < scale) ? 0. : pow(_fissionPwtSquark,scale); else if (_enhanceSProb == 2) prob_s = (_maxScale < scale) ? 0. : exp(-scale); break; } int choice = UseRandom::rnd3(prob_u, prob_d, prob_s); long idNew = 0; switch (choice) { case 0: idNew = ThePEG::ParticleID::u; break; case 1: idNew = ThePEG::ParticleID::d; break; case 2: idNew = ThePEG::ParticleID::s; break; } newPtrPos = getParticle(idNew); newPtrNeg = getParticle(-idNew); assert (newPtrPos); assert(newPtrNeg); assert (newPtrPos->dataPtr()); assert(newPtrNeg->dataPtr()); } Energy2 ClusterFissioner::clustermass(const ClusterPtr & cluster){ Lorentz5Momentum pIn = cluster->momentum(); Energy2 endpointmass2 = sqr(cluster->particle(0)->mass() + cluster->particle(1)->mass()); Energy2 singletm2 = pIn.m2(); // Return either the cluster mass, or the lambda measure return (_massMeasure == 0) ? singletm2 : singletm2 - endpointmass2; } Energy ClusterFissioner::drawChildMass(const Energy M, const Energy m1, const Energy m2, const Energy m, const double expt, const bool soft) const { /*************************** * This method, given in input the cluster mass Mclu of an heavy cluster C, * made of consituents of masses m1 and m2, draws the masses Mclu1 and Mclu2 * of, respectively, the children cluster C1, made of constituent masses m1 * and m, and cluster C2, of mass Mclu2 and made of constituent masses m2 * and m. The mass is extracted from one of the two following mass * distributions: * --- power-like ("normal" distribution) * d(Prob) / d(M^exponent) = const * where the exponent can be different from the two children C1 (exp1) * and C2 (exponent2). * --- exponential ("soft" distribution) * d(Prob) / d(M^2) = exp(-b*M) * where b = 2.0 / average. * Such distributions are limited below by the masses of * the constituents quarks, and above from the mass of decaying cluster C. * The choice of which of the two mass distributions to use for each of the * two cluster children is dictated by iRemnant (see below). * If the number of attempts to extract a pair of mass values that are * kinematically acceptable is above some fixed number (max_loop, see below) * the method gives up and returns false; otherwise, when it succeeds, it * returns true. * * These distributions have been modified from HERWIG: * Before these were: * Mclu1 = m1 + (Mclu - m1 - m2)*pow( rnd(), 1.0/exponent1 ); * The new one coded here is a more efficient version, same density * but taking into account 'in phase space from' beforehand ***************************/ // hard cluster if(!soft) { return pow(UseRandom::rnd(pow((M-m1-m2-m)*UnitRemoval::InvE, expt), pow(m*UnitRemoval::InvE, expt)), 1./expt )*UnitRemoval::E + m1; } // Otherwise it uses a soft mass distribution else { static const InvEnergy b = 2.0 / _btClM; Energy max = M-m1-m2-2.0*m; double rmin = b*max; rmin = ( rmin < 50 ) ? exp(-rmin) : 0.; double r1; do { r1 = UseRandom::rnd(rmin, 1.0) * UseRandom::rnd(rmin, 1.0); } while (r1 < rmin); return m1 + m - log(r1)/b; } } void ClusterFissioner::calculateKinematics(const Lorentz5Momentum & pClu, const Lorentz5Momentum & p0Q1, const bool toHadron1, const bool toHadron2, Lorentz5Momentum & pClu1, Lorentz5Momentum & pClu2, Lorentz5Momentum & pQ1, Lorentz5Momentum & pQbar, Lorentz5Momentum & pQ, Lorentz5Momentum & pQ2bar) const { /****************** * This method solves the kinematics of the two body cluster decay: * C (Q1 Q2bar) ---> C1 (Q1 Qbar) + C2 (Q Q2bar) * In input we receive the momentum of C, pClu, and the momentum * of the quark Q1 (constituent of C), p0Q1, both in the LAB frame. * Furthermore, two boolean variables inform whether the two fission * products (C1, C2) decay immediately into a single hadron (in which * case the cluster itself is identify with that hadron) and we do * not have to solve the kinematics of the components (Q1,Qbar) for * C1 and (Q,Q2bar) for C2. * The output is given by the following momenta (all 5-components, * and all in the LAB frame): * pClu1 , pClu2 respectively of C1 , C2 * pQ1 , pQbar respectively of Q1 , Qbar in C1 * pQ , pQ2bar respectively of Q , Q2 in C2 * The assumption, suggested from the string model, is that, in C frame, * C1 and its constituents Q1 and Qbar are collinear, and collinear to * the direction of Q1 in C (that is before cluster decay); similarly, * (always in the C frame) C2 and its constituents Q and Q2bar are * collinear (and therefore anti-collinear with C1,Q1,Qbar). * The solution is then obtained by using Lorentz boosts, as follows. * The kinematics of C1 and C2 is solved in their parent C frame, * and then boosted back in the LAB. The kinematics of Q1 and Qbar * is solved in their parent C1 frame and then boosted back in the LAB; * similarly, the kinematics of Q and Q2bar is solved in their parent * C2 frame and then boosted back in the LAB. In each of the three * "two-body decay"-like cases, we use the fact that the direction * of the motion of the decay products is known in the rest frame of * their parent. This is obvious for the first case in which the * parent rest frame is C; but it is also true in the other two cases * where the rest frames are C1 and C2. This is because C1 and C2 * are boosted w.r.t. C in the same direction where their components, * respectively (Q1,Qbar) and (Q,Q2bar) move in C1 and C2 rest frame * respectively. * Of course, although the notation used assumed that C = (Q1 Q2bar) * where Q1 is a quark and Q2bar an antiquark, indeed everything remain * unchanged also in all following cases: * Q1 quark, Q2bar antiquark; --> Q quark; * Q1 antiquark , Q2bar quark; --> Q antiquark; * Q1 quark, Q2bar diquark; --> Q quark * Q1 antiquark, Q2bar anti-diquark; --> Q antiquark * Q1 diquark, Q2bar quark --> Q antiquark * Q1 anti-diquark, Q2bar antiquark; --> Q quark **************************/ // Calculate the unit three-vector, in the C frame, along which // all of the constituents and children clusters move. Lorentz5Momentum u(p0Q1); u.boost( -pClu.boostVector() ); // boost from LAB to C // the unit three-vector is then u.vect().unit() // Calculate the momenta of C1 and C2 in the (parent) C frame first, // where the direction of C1 is u.vect().unit(), and then boost back in the // LAB frame. if (pClu.m() < pClu1.mass() + pClu2.mass() ) { throw Exception() << "Impossible Kinematics in ClusterFissioner::calculateKinematics() (A)" << Exception::eventerror; } Kinematics::twoBodyDecay(pClu, pClu1.mass(), pClu2.mass(), u.vect().unit(), pClu1, pClu2); // In the case that cluster1 does not decay immediately into a single hadron, // calculate the momenta of Q1 (as constituent of C1) and Qbar in the // (parent) C1 frame first, where the direction of Q1 is u.vect().unit(), // and then boost back in the LAB frame. if(!toHadron1) { if (pClu1.m() < pQ1.mass() + pQbar.mass() ) { throw Exception() << "Impossible Kinematics in ClusterFissioner::calculateKinematics() (B)" << Exception::eventerror; } Kinematics::twoBodyDecay(pClu1, pQ1.mass(), pQbar.mass(), u.vect().unit(), pQ1, pQbar); } // In the case that cluster2 does not decay immediately into a single hadron, // Calculate the momenta of Q and Q2bar (as constituent of C2) in the // (parent) C2 frame first, where the direction of Q is u.vect().unit(), // and then boost back in the LAB frame. if(!toHadron2) { if (pClu2.m() < pQ.mass() + pQ2bar.mass() ) { throw Exception() << "Impossible Kinematics in ClusterFissioner::calculateKinematics() (C)" << Exception::eventerror; } Kinematics::twoBodyDecay(pClu2, pQ.mass(), pQ2bar.mass(), u.vect().unit(), pQ, pQ2bar); } } void ClusterFissioner::calculatePositions(const Lorentz5Momentum & pClu, const LorentzPoint & positionClu, const Lorentz5Momentum & pClu1, const Lorentz5Momentum & pClu2, LorentzPoint & positionClu1, LorentzPoint & positionClu2) const { // Determine positions of cluster children. // See Marc Smith's thesis, page 127, formulas (4.122) and (4.123). Energy Mclu = pClu.m(); Energy Mclu1 = pClu1.m(); Energy Mclu2 = pClu2.m(); // Calculate the unit three-vector, in the C frame, along which // children clusters move. Lorentz5Momentum u(pClu1); u.boost( -pClu.boostVector() ); // boost from LAB to C frame // the unit three-vector is then u.vect().unit() Energy pstarChild = Kinematics::pstarTwoBodyDecay(Mclu,Mclu1,Mclu2); // First, determine the relative positions of the children clusters // in the parent cluster reference frame. Length x1 = ( 0.25*Mclu + 0.5*( pstarChild + (sqr(Mclu2) - sqr(Mclu1))/(2.0*Mclu)))/_kappa; Length t1 = Mclu/_kappa - x1; LorentzDistance distanceClu1( x1 * u.vect().unit(), t1 ); Length x2 = (-0.25*Mclu + 0.5*(-pstarChild + (sqr(Mclu2) - sqr(Mclu1))/(2.0*Mclu)))/_kappa; Length t2 = Mclu/_kappa + x2; LorentzDistance distanceClu2( x2 * u.vect().unit(), t2 ); // Then, transform such relative positions from the parent cluster // reference frame to the Lab frame. distanceClu1.boost( pClu.boostVector() ); distanceClu2.boost( pClu.boostVector() ); // Finally, determine the absolute positions in the Lab frame. positionClu1 = positionClu + distanceClu1; positionClu2 = positionClu + distanceClu2; } bool ClusterFissioner::isHeavy(tcClusterPtr clu) { // default double clpow = _clPowLight; Energy clmax = _clMaxLight; // particle data for constituents tcPDPtr cptr[3]={tcPDPtr(),tcPDPtr(),tcPDPtr()}; for(int ix=0;ixnumComponents(),3);++ix) { cptr[ix]=clu->particle(ix)->dataPtr(); } // different parameters for exotic, bottom and charm clusters if(CheckId::isExotic(cptr[0],cptr[1],cptr[1])) { clpow = _clPowExotic; clmax = _clMaxExotic; } else if(CheckId::hasBottom(cptr[0],cptr[1],cptr[1])) { clpow = _clPowBottom; clmax = _clMaxBottom; } else if(CheckId::hasCharm(cptr[0],cptr[1],cptr[1])) { clpow = _clPowCharm; clmax = _clMaxCharm; } bool aboveCutoff = ( pow(clu->mass()*UnitRemoval::InvE , clpow) > pow(clmax*UnitRemoval::InvE, clpow) + pow(clu->sumConstituentMasses()*UnitRemoval::InvE, clpow) ); // required test for SUSY clusters, since aboveCutoff alone // cannot guarantee (Mc > m1 + m2 + 2*m) in cut() static const Energy minmass = getParticleData(ParticleID::d)->constituentMass(); bool canSplitMinimally = clu->mass() > clu->sumConstituentMasses() + 2.0 * minmass; return aboveCutoff && canSplitMinimally; } diff --git a/Hadronization/ClusterHadronizationHandler.cc b/Hadronization/ClusterHadronizationHandler.cc --- a/Hadronization/ClusterHadronizationHandler.cc +++ b/Hadronization/ClusterHadronizationHandler.cc @@ -1,307 +1,307 @@ // -*- C++ -*- // // ClusterHadronizationHandler.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the ClusterHadronizationHandler class. // #include "ClusterHadronizationHandler.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "Herwig/Utilities/EnumParticles.h" #include "CluHadConfig.h" #include "Cluster.h" #include using namespace Herwig; ClusterHadronizationHandler * ClusterHadronizationHandler::currentHandler_ = 0; DescribeClass -describeClusterHadronizationHandler("Herwig::ClusterHadronizationHandler",""); +describeClusterHadronizationHandler("Herwig::ClusterHadronizationHandler","Herwig.so"); IBPtr ClusterHadronizationHandler::clone() const { return new_ptr(*this); } IBPtr ClusterHadronizationHandler::fullclone() const { return new_ptr(*this); } void ClusterHadronizationHandler::persistentOutput(PersistentOStream & os) const { os << _partonSplitter << _clusterFinder << _colourReconnector << _clusterFissioner << _lightClusterDecayer << _clusterDecayer << ounit(_minVirtuality2,GeV2) << ounit(_maxDisplacement,mm) << _underlyingEventHandler << _reduceToTwoComponents; } void ClusterHadronizationHandler::persistentInput(PersistentIStream & is, int) { is >> _partonSplitter >> _clusterFinder >> _colourReconnector >> _clusterFissioner >> _lightClusterDecayer >> _clusterDecayer >> iunit(_minVirtuality2,GeV2) >> iunit(_maxDisplacement,mm) >> _underlyingEventHandler >> _reduceToTwoComponents; } void ClusterHadronizationHandler::Init() { static ClassDocumentation documentation ("This is the main handler class for the Cluster Hadronization", "The hadronization was performed using the cluster model of \\cite{Webber:1983if}.", "%\\cite{Webber:1983if}\n" "\\bibitem{Webber:1983if}\n" " B.~R.~Webber,\n" " ``A QCD Model For Jet Fragmentation Including Soft Gluon Interference,''\n" " Nucl.\\ Phys.\\ B {\\bf 238}, 492 (1984).\n" " %%CITATION = NUPHA,B238,492;%%\n" // main manual ); static Reference interfacePartonSplitter("PartonSplitter", "A reference to the PartonSplitter object", &Herwig::ClusterHadronizationHandler::_partonSplitter, false, false, true, false); static Reference interfaceClusterFinder("ClusterFinder", "A reference to the ClusterFinder object", &Herwig::ClusterHadronizationHandler::_clusterFinder, false, false, true, false); static Reference interfaceColourReconnector("ColourReconnector", "A reference to the ColourReconnector object", &Herwig::ClusterHadronizationHandler::_colourReconnector, false, false, true, false); static Reference interfaceClusterFissioner("ClusterFissioner", "A reference to the ClusterFissioner object", &Herwig::ClusterHadronizationHandler::_clusterFissioner, false, false, true, false); static Reference interfaceLightClusterDecayer("LightClusterDecayer", "A reference to the LightClusterDecayer object", &Herwig::ClusterHadronizationHandler::_lightClusterDecayer, false, false, true, false); static Reference interfaceClusterDecayer("ClusterDecayer", "A reference to the ClusterDecayer object", &Herwig::ClusterHadronizationHandler::_clusterDecayer, false, false, true, false); static Parameter interfaceMinVirtuality2 ("MinVirtuality2", "Minimum virtuality^2 of partons to use in calculating distances (unit [GeV2]).", &ClusterHadronizationHandler::_minVirtuality2, GeV2, 0.1*GeV2, ZERO, 10.0*GeV2,false,false,false); static Parameter interfaceMaxDisplacement ("MaxDisplacement", "Maximum displacement that is allowed for a particle (unit [millimeter]).", &ClusterHadronizationHandler::_maxDisplacement, mm, 1.0e-10*mm, 0.0*mm, 1.0e-9*mm,false,false,false); static Reference interfaceUnderlyingEventHandler ("UnderlyingEventHandler", "Pointer to the handler for the Underlying Event. " "Set to NULL to disable.", &ClusterHadronizationHandler::_underlyingEventHandler, false, false, true, true, false); static Switch interfaceReduceToTwoComponents ("ReduceToTwoComponents", "Whether or not to reduce three component baryon-number violating clusters to two components before cluster splitting or leave" " this till after the cluster splitting", &ClusterHadronizationHandler::_reduceToTwoComponents, true, false, false); static SwitchOption interfaceReduceToTwoComponentsYes (interfaceReduceToTwoComponents, "BeforeSplitting", "Reduce to two components", true); static SwitchOption interfaceReduceToTwoComponentsNo (interfaceReduceToTwoComponents, "AfterSplitting", "Treat as three components", false); } namespace { void extractChildren(tPPtr p, set & all) { if (p->children().empty()) return; for (PVector::const_iterator child = p->children().begin(); child != p->children().end(); ++child) { all.insert(*child); extractChildren(*child, all); } } } void ClusterHadronizationHandler:: handle(EventHandler & ch, const tPVector & tagged, const Hint &) { useMe(); currentHandler_ = this; PVector currentlist(tagged.begin(),tagged.end()); // set the scale for coloured particles to just above the gluon mass squared // if less than this so they are classed as perturbative Energy2 Q02 = 1.01*sqr(getParticleData(ParticleID::g)->constituentMass()); for(unsigned int ix=0;ixscale()scale(Q02); } // split the gluons _partonSplitter->split(currentlist); // form the clusters ClusterVector clusters = _clusterFinder->formClusters(currentlist); // reduce BV clusters to two components now if needed if(_reduceToTwoComponents) _clusterFinder->reduceToTwoComponents(clusters); // perform colour reconnection if needed and then // decay the clusters into one hadron bool lightOK = false; short tried = 0; const ClusterVector savedclusters = clusters; tPVector finalHadrons; // only needed for partonic decayer while (!lightOK && tried++ < 10) { // no colour reconnection with baryon-number-violating (BV) clusters ClusterVector CRclusters, BVclusters; CRclusters.reserve( clusters.size() ); BVclusters.reserve( clusters.size() ); for (size_t ic = 0; ic < clusters.size(); ++ic) { ClusterPtr cl = clusters.at(ic); bool hasClusterParent = false; for (unsigned int ix=0; ix < cl->parents().size(); ++ix) { if (cl->parents()[ix]->id() == ParticleID::Cluster) { hasClusterParent = true; break; } } if (cl->numComponents() > 2 || hasClusterParent) BVclusters.push_back(cl); else CRclusters.push_back(cl); } // colour reconnection _colourReconnector->rearrange(CRclusters); // tag new clusters as children of the partons to hadronize _setChildren(CRclusters); // forms diquarks _clusterFinder->reduceToTwoComponents(CRclusters); // recombine vectors of (possibly) reconnected and BV clusters clusters.clear(); clusters.insert( clusters.end(), CRclusters.begin(), CRclusters.end() ); clusters.insert( clusters.end(), BVclusters.begin(), BVclusters.end() ); // fission of heavy clusters // NB: during cluster fission, light hadrons might be produced straight away finalHadrons = _clusterFissioner->fission(clusters,isSoftUnderlyingEventON()); // if clusters not previously reduced to two components do it now if(!_reduceToTwoComponents) _clusterFinder->reduceToTwoComponents(clusters); lightOK = _lightClusterDecayer->decay(clusters,finalHadrons); // if the decay of the light clusters was not successful, undo the cluster // fission and decay steps and revert to the original state of the event // record if (!lightOK) { clusters = savedclusters; for_each(clusters.begin(), clusters.end(), mem_fun(&Particle::undecay)); } } if (!lightOK) { throw Exception("CluHad::handle(): tried LightClusterDecayer 10 times!", Exception::eventerror); } // decay the remaining clusters _clusterDecayer->decay(clusters,finalHadrons); // ***************************************** // ***************************************** // ***************************************** StepPtr pstep = newStep(); set allDecendants; for (tPVector::const_iterator it = tagged.begin(); it != tagged.end(); ++it) { extractChildren(*it, allDecendants); } for(set::const_iterator it = allDecendants.begin(); it != allDecendants.end(); ++it) { // this is a workaround because the set sometimes // re-orders parents after their children if ((*it)->children().empty()) pstep->addDecayProduct(*it); else { pstep->addDecayProduct(*it); pstep->addIntermediate(*it); } } // ***************************************** // ***************************************** // ***************************************** // soft underlying event if needed if (isSoftUnderlyingEventON()) { assert(_underlyingEventHandler); ch.performStep(_underlyingEventHandler,Hint::Default()); } } // Sets parent child relationship of all clusters with two components // Relationships for clusters with more than two components are set elsewhere in the Colour Reconnector void ClusterHadronizationHandler::_setChildren(const ClusterVector & clusters) const { // erase existing information about the partons' children tPVector partons; for ( const auto & cl : clusters ) { if ( cl->numComponents() > 2 ) continue; partons.push_back( cl->colParticle() ); partons.push_back( cl->antiColParticle() ); } // erase all previous information about parent child relationship for_each(partons.begin(), partons.end(), mem_fun(&Particle::undecay)); // give new parents to the clusters: their constituents for ( const auto & cl : clusters ) { if ( cl->numComponents() > 2 ) continue; cl->colParticle()->addChild(cl); cl->antiColParticle()->addChild(cl); } } diff --git a/Hadronization/ColourReconnector.cc b/Hadronization/ColourReconnector.cc --- a/Hadronization/ColourReconnector.cc +++ b/Hadronization/ColourReconnector.cc @@ -1,822 +1,822 @@ // -*- C++ -*- // // ColourReconnector.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the ColourReconnector class. // #include "ColourReconnector.h" #include "Cluster.h" #include #include #include #include #include #include #include #include "Herwig/Utilities/Maths.h" using namespace Herwig; using CluVecIt = ColourReconnector::CluVecIt; using Constants::pi; using Constants::twopi; DescribeClass -describeColourReconnector("Herwig::ColourReconnector",""); +describeColourReconnector("Herwig::ColourReconnector","Herwig.so"); IBPtr ColourReconnector::clone() const { return new_ptr(*this); } IBPtr ColourReconnector::fullclone() const { return new_ptr(*this); } void ColourReconnector::rearrange(ClusterVector & clusters) { if (_clreco == 0) return; // need at least two clusters if (clusters.size() < 2) return; // do the colour reconnection switch (_algorithm) { case 0: _doRecoPlain(clusters); break; case 1: _doRecoStatistical(clusters); break; case 2: _doRecoBaryonic(clusters); break; } } Energy2 ColourReconnector::_clusterMassSum(const PVector & q, const PVector & aq) const { const size_t nclusters = q.size(); assert (aq.size() == nclusters); Energy2 sum = ZERO; for (size_t i = 0; i < nclusters; i++) sum += ( q[i]->momentum() + aq[i]->momentum() ).m2(); return sum; } bool ColourReconnector::_containsColour8(const ClusterVector & cv, const vector & P) const { assert (P.size() == cv.size()); for (size_t i = 0; i < cv.size(); i++) { tcPPtr p = cv[i]->colParticle(); tcPPtr q = cv[P[i]]->antiColParticle(); if (_isColour8(p, q)) return true; } return false; } void ColourReconnector::_doRecoStatistical(ClusterVector & cv) const { const size_t nclusters = cv.size(); // initially, enumerate (anti)quarks as given in the cluster vector ParticleVector q, aq; for (size_t i = 0; i < nclusters; i++) { q.push_back( cv[i]->colParticle() ); aq.push_back( cv[i]->antiColParticle() ); } // annealing scheme Energy2 t, delta; Energy2 lambda = _clusterMassSum(q,aq); const unsigned _ntries = _triesPerStepFactor * nclusters; // find appropriate starting temperature by measuring the largest lambda // difference in some dry-run random rearrangements { vector typical; for (int i = 0; i < 10; i++) { const pair toswap = _shuffle(q,aq,5); ParticleVector newaq = aq; swap (newaq[toswap.first], newaq[toswap.second]); Energy2 newlambda = _clusterMassSum(q,newaq); typical.push_back( abs(newlambda - lambda) ); } t = _initTemp * Math::median(typical); } // anneal in up to _annealingSteps temperature steps for (unsigned step = 0; step < _annealingSteps; step++) { // For this temperature step, try to reconnect _ntries times. Stop the // algorithm if no successful reconnection happens. unsigned nSuccess = 0; for (unsigned it = 0; it < _ntries; it++) { // make a random rearrangement const unsigned maxtries = 10; const pair toswap = _shuffle(q,aq,maxtries); const int i = toswap.first; const int j = toswap.second; // stop here if we cannot find any allowed reconfiguration if (i == -1) break; // create a new antiquark vector with the two partons swapped ParticleVector newaq = aq; swap (newaq[i], newaq[j]); // Check if lambda would decrease. If yes, accept the reconnection. If no, // accept it only with a probability given by the current Boltzmann // factor. In the latter case we set p = 0 if the temperature is close to // 0, to avoid division by 0. Energy2 newlambda = _clusterMassSum(q,newaq); delta = newlambda - lambda; double prob = 1.0; if (delta > ZERO) prob = ( abs(t) < 1e-8*MeV2 ) ? 0.0 : exp(-delta/t); if (UseRandom::rnd() < prob) { lambda = newlambda; swap (newaq, aq); nSuccess++; } } if (nSuccess == 0) break; // reduce temperature t *= _annealingFactor; } // construct the new cluster vector ClusterVector newclusters; for (size_t i = 0; i < nclusters; i++) { ClusterPtr cl = new_ptr( Cluster( q[i], aq[i] ) ); newclusters.push_back(cl); } swap(newclusters,cv); return; } void ColourReconnector::_doRecoPlain(ClusterVector & cv) const { ClusterVector newcv = cv; // try to avoid systematic errors by randomising the reconnection order long (*p_irnd)(long) = UseRandom::irnd; random_shuffle( newcv.begin(), newcv.end(), p_irnd ); // iterate over all clusters for (CluVecIt cit = newcv.begin(); cit != newcv.end(); cit++) { // find the cluster which, if reconnected with *cit, would result in the // smallest sum of cluster masses // NB this method returns *cit if no reconnection partner can be found CluVecIt candidate = _findRecoPartner(cit, newcv); // skip this cluster if no possible reshuffling partner can be found if (candidate == cit) continue; // accept the reconnection with probability _preco. if (UseRandom::rnd() < _preco) { pair reconnected = _reconnect(*cit, *candidate); // Replace the clusters in the ClusterVector. The order of the // colour-triplet partons in the cluster vector is retained here. // replace *cit by reconnected.first *cit = reconnected.first; // replace candidate by reconnected.second *candidate = reconnected.second; } } swap(cv,newcv); return; } namespace { inline bool hasDiquark(CluVecIt cit) { for(int i = 0; i<(*cit)->numComponents(); i++) { if (DiquarkMatcher::Check(*((*cit)->particle(i)->dataPtr()))) return true; } return false; } } // Implementation of the baryonic reconnection algorithm void ColourReconnector::_doRecoBaryonic(ClusterVector & cv) const { ClusterVector newcv = cv; ClusterVector deleted; deleted.reserve(cv.size()); // try to avoid systematic errors by randomising the reconnection order long (*p_irnd)(long) = UseRandom::irnd; random_shuffle( newcv.begin(), newcv.end(), p_irnd ); // iterate over all clusters for (CluVecIt cit = newcv.begin(); cit != newcv.end(); ++cit) { //avoid clusters already containing diuarks if (hasDiquark(cit)) continue; //skip the cluster to be deleted later 3->2 cluster if (find(deleted.begin(), deleted.end(), *cit) != deleted.end()) continue; // Skip all found baryonic clusters, this biases the algorithm but implementing // something like re-reconnection is ongoing work if ((*cit)->numComponents()==3) continue; // Find a candidate suitable for reconnection CluVecIt baryonic1, baryonic2; bool isBaryonicCandidate = false; CluVecIt candidate = _findPartnerBaryonic(cit, newcv, isBaryonicCandidate, deleted, baryonic1, baryonic2); // skip this cluster if no possible reconnection partner can be found if ( !isBaryonicCandidate && candidate==cit ) continue; if ( isBaryonicCandidate && UseRandom::rnd() < _precoBaryonic ) { deleted.push_back(*baryonic2); // Function that does the reconnection from 3 -> 2 clusters ClusterPtr b1, b2; _makeBaryonicClusters(*cit,*baryonic1,*baryonic2, b1, b2); *cit = b1; *baryonic1 = b2; // Baryonic2 is easily skipped in the next loop } // Normal 2->2 Colour reconnection if ( !isBaryonicCandidate && UseRandom::rnd() < _preco ) { auto reconnected = _reconnectBaryonic(*cit, *candidate); *cit = reconnected.first; *candidate = reconnected.second; } } // create a new vector of clusters except for the ones which are "deleted" during // baryonic reconnection ClusterVector clustervector; for ( const auto & cluster : newcv ) if ( find(deleted.begin(), deleted.end(), cluster) == deleted.end() ) clustervector.push_back(cluster); swap(cv,clustervector); } 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 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; } } CluVecIt ColourReconnector::_findPartnerBaryonic( CluVecIt cl, ClusterVector & cv, bool & baryonicCand, const ClusterVector& deleted, CluVecIt &baryonic1, CluVecIt &baryonic2 ) const { using Constants::pi; using Constants::twopi; // Returns a candidate for possible reconnection CluVecIt candidate = cl; bool bcand = false; double maxrap = 0.0; double minrap = 0.0; double maxrapNormal = 0.0; double minrapNormal = 0.0; double maxsumnormal = 0.0; double maxsum = 0.0; double secondsum = 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; if ( (*cit)->numComponents()==3 ) continue; 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 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; bcand = false; candidate = cit; } } if ( rapq < 0.0 && rapqbar >0.0 && rapqbar > maxrapNormal && rapq < minrapNormal ) { maxrap = rapqbar; minrap = rapq; const double sumrap = abs(rapqbar) + abs(rapq); // first candidate gets here. If second baryonic candidate has higher Ysum than the first // one, the second candidate becomes the first one and the first the second. if (sumrap > maxsum) { if(maxsum != 0){ baryonic2 = baryonic1; baryonic1 = cit; bcand = true; } else { baryonic1 = cit; } maxsum = sumrap; } else { if (sumrap > secondsum && sumrap != maxsum) { secondsum = sumrap; bcand = true; baryonic2 = cit; } } } } if(bcand == true){ baryonicCand = true; } return candidate; } CluVecIt ColourReconnector::_findRecoPartner(CluVecIt cl, ClusterVector & cv) const { CluVecIt candidate = cl; Energy minMass = 1*TeV; for (CluVecIt cit=cv.begin(); cit != cv.end(); ++cit) { // don't even look at original cluster if(cit==cl) continue; // don't allow colour octet clusters if ( _isColour8( (*cl)->colParticle(), (*cit)->antiColParticle() ) || _isColour8( (*cit)->colParticle(), (*cl)->antiColParticle() ) ) { continue; } // stop it putting beam remnants together if((*cl)->isBeamCluster() && (*cit)->isBeamCluster()) continue; // stop it putting far apart clusters together if(((**cl).vertex()-(**cit).vertex()).m()>_maxDistance) continue; // momenta of the old clusters Lorentz5Momentum p1 = (*cl)->colParticle()->momentum() + (*cl)->antiColParticle()->momentum(); Lorentz5Momentum p2 = (*cit)->colParticle()->momentum() + (*cit)->antiColParticle()->momentum(); // momenta of the new clusters Lorentz5Momentum p3 = (*cl)->colParticle()->momentum() + (*cit)->antiColParticle()->momentum(); Lorentz5Momentum p4 = (*cit)->colParticle()->momentum() + (*cl)->antiColParticle()->momentum(); Energy oldMass = abs( p1.m() ) + abs( p2.m() ); Energy newMass = abs( p3.m() ) + abs( p4.m() ); if ( newMass < oldMass && newMass < minMass ) { minMass = newMass; candidate = cit; } } return candidate; } // forms two baryonic clusters from three clusters void ColourReconnector::_makeBaryonicClusters( ClusterPtr &c1, ClusterPtr &c2, ClusterPtr &c3, ClusterPtr &newcluster1, ClusterPtr &newcluster2) const{ //make sure they all have 2 components assert(c1->numComponents()==2); assert(c2->numComponents()==2); assert(c3->numComponents()==2); //abandon children c1->colParticle()->abandonChild(c1); c1->antiColParticle()->abandonChild(c1); c2->colParticle()->abandonChild(c2); c2->antiColParticle()->abandonChild(c2); c3->colParticle()->abandonChild(c3); c3->antiColParticle()->abandonChild(c3); newcluster1 = new_ptr(Cluster(c1->colParticle(),c2->colParticle(), c3->colParticle())); c1->colParticle()->addChild(newcluster1); c2->colParticle()->addChild(newcluster1); c3->colParticle()->addChild(newcluster1); newcluster1->setVertex(LorentzPoint()); newcluster2 = new_ptr(Cluster(c1->antiColParticle(), c2->antiColParticle(), c3->antiColParticle())); c1->antiColParticle()->addChild(newcluster2); c2->antiColParticle()->addChild(newcluster2); c3->antiColParticle()->addChild(newcluster2); newcluster2->setVertex(LorentzPoint()); } pair ColourReconnector::_reconnect(ClusterPtr &c1, ClusterPtr &c2) const { // choose the other possibility to form two clusters from the given // constituents assert(c1->numComponents()==2); assert(c2->numComponents()==2); int c1_col(-1),c1_anti(-1),c2_col(-1),c2_anti(-1); for(unsigned int ix=0;ix<2;++ix) { if (c1->particle(ix)->hasColour(false)) c1_col = ix; else if(c1->particle(ix)->hasColour(true )) c1_anti = ix; if (c2->particle(ix)->hasColour(false)) c2_col = ix; else if(c2->particle(ix)->hasColour(true )) c2_anti = ix; } assert(c1_col>=0&&c2_col>=0&&c1_anti>=0&&c2_anti>=0); ClusterPtr newCluster1 = new_ptr( Cluster( c1->colParticle(), c2->antiColParticle() ) ); newCluster1->setVertex(0.5*( c1->colParticle()->vertex() + c2->antiColParticle()->vertex() )); if(c1->isBeamRemnant(c1_col )) newCluster1->setBeamRemnant(0,true); if(c2->isBeamRemnant(c2_anti)) newCluster1->setBeamRemnant(1,true); ClusterPtr newCluster2 = new_ptr( Cluster( c2->colParticle(), c1->antiColParticle() ) ); newCluster2->setVertex(0.5*( c2->colParticle()->vertex() + c1->antiColParticle()->vertex() )); if(c2->isBeamRemnant(c2_col )) newCluster2->setBeamRemnant(0,true); if(c1->isBeamRemnant(c1_anti)) newCluster2->setBeamRemnant(1,true); return pair (newCluster1, newCluster2); } pair ColourReconnector::_reconnectBaryonic(ClusterPtr &c1, ClusterPtr &c2) const { // choose the other possibility to form two clusters from the given // constituents assert(c1->numComponents()==2); assert(c2->numComponents()==2); int c1_col(-1),c1_anti(-1),c2_col(-1),c2_anti(-1); for(unsigned int ix=0;ix<2;++ix) { if (c1->particle(ix)->hasColour(false)) c1_col = ix; else if(c1->particle(ix)->hasColour(true )) c1_anti = ix; if (c2->particle(ix)->hasColour(false)) c2_col = ix; else if(c2->particle(ix)->hasColour(true )) c2_anti = ix; } assert(c1_col>=0&&c2_col>=0&&c1_anti>=0&&c2_anti>=0); c1->colParticle()->abandonChild(c1); c2->antiColParticle()->abandonChild(c2); ClusterPtr newCluster1 = new_ptr( Cluster( c1->colParticle(), c2->antiColParticle() ) ); c1->colParticle()->addChild(newCluster1); c2->antiColParticle()->addChild(newCluster1); newCluster1->setVertex(0.5*( c1->colParticle()->vertex() + c2->antiColParticle()->vertex() )); if(c1->isBeamRemnant(c1_col )) newCluster1->setBeamRemnant(0,true); if(c2->isBeamRemnant(c2_anti)) newCluster1->setBeamRemnant(1,true); c1->antiColParticle()->abandonChild(c1); c2->colParticle()->abandonChild(c2); ClusterPtr newCluster2 = new_ptr( Cluster( c2->colParticle(), c1->antiColParticle() ) ); c1->antiColParticle()->addChild(newCluster2); c2->colParticle()->addChild(newCluster2); newCluster2->setVertex(0.5*( c2->colParticle()->vertex() + c1->antiColParticle()->vertex() )); if(c2->isBeamRemnant(c2_col )) newCluster2->setBeamRemnant(0,true); if(c1->isBeamRemnant(c1_anti)) newCluster2->setBeamRemnant(1,true); return pair (newCluster1, newCluster2); } pair ColourReconnector::_shuffle (const PVector & q, const PVector & aq, unsigned maxtries) const { const size_t nclusters = q.size(); assert (nclusters > 1); assert (aq.size() == nclusters); int i, j; unsigned tries = 0; bool octet; do { // find two different random integers in the range [0, nclusters) i = UseRandom::irnd( nclusters ); do { j = UseRandom::irnd( nclusters ); } while (i == j); // check if one of the two potential clusters would be a colour octet state octet = _isColour8( q[i], aq[j] ) || _isColour8( q[j], aq[i] ) ; tries++; } while (octet && tries < maxtries); if (octet) i = j = -1; return make_pair(i,j); } bool ColourReconnector::_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; } void ColourReconnector::persistentOutput(PersistentOStream & os) const { os << _clreco << _preco << _precoBaryonic << _algorithm << _initTemp << _annealingFactor << _annealingSteps << _triesPerStepFactor << ounit(_maxDistance,femtometer) << _octetOption; } void ColourReconnector::persistentInput(PersistentIStream & is, int) { is >> _clreco >> _preco >> _precoBaryonic >> _algorithm >> _initTemp >> _annealingFactor >> _annealingSteps >> _triesPerStepFactor >> iunit(_maxDistance,femtometer) >> _octetOption; } void ColourReconnector::Init() { static ClassDocumentation documentation ("This class is responsible of the colour reconnection."); static Switch interfaceColourReconnection ("ColourReconnection", "Colour reconnections", &ColourReconnector::_clreco, 0, true, false); static SwitchOption interfaceColourReconnectionNo (interfaceColourReconnection, "No", "Colour reconnections off", 0); static SwitchOption interfaceColourReconnectionYes (interfaceColourReconnection, "Yes", "Colour reconnections on", 1); static Parameter interfaceMtrpAnnealingFactor ("AnnealingFactor", "The annealing factor is the ratio of the temperatures in two successive " "temperature steps.", &ColourReconnector::_annealingFactor, 0.9, 0.0, 1.0, false, false, Interface::limited); static Parameter interfaceMtrpAnnealingSteps ("AnnealingSteps", "Number of temperature steps in the statistical annealing algorithm", &ColourReconnector::_annealingSteps, 50, 1, 10000, false, false, Interface::limited); static Parameter interfaceMtrpTriesPerStepFactor ("TriesPerStepFactor", "The number of reconnection tries per temperature steps is the number of " "clusters times this factor.", &ColourReconnector::_triesPerStepFactor, 5.0, 0.0, 100.0, false, false, Interface::limited); static Parameter interfaceMtrpInitialTemp ("InitialTemperature", "Factor used to determine the initial temperature from the median of the " "energy change in a few random rearrangements.", &ColourReconnector::_initTemp, 0.1, 0.00001, 100.0, false, false, Interface::limited); static Parameter interfaceRecoProb ("ReconnectionProbability", "Probability that a found reconnection possibility is actually accepted", &ColourReconnector::_preco, 0.5, 0.0, 1.0, false, false, Interface::limited); static Parameter interfaceRecoProbBaryonic ("ReconnectionProbabilityBaryonic", "Probability that a found reconnection possibility is actually accepted", &ColourReconnector::_precoBaryonic, 0.5, 0.0, 1.0, false, false, Interface::limited); static Switch interfaceAlgorithm ("Algorithm", "Specifies the colour reconnection algorithm", &ColourReconnector::_algorithm, 0, true, false); static SwitchOption interfaceAlgorithmPlain (interfaceAlgorithm, "Plain", "Plain colour reconnection as in Herwig 2.5.0", 0); static SwitchOption interfaceAlgorithmStatistical (interfaceAlgorithm, "Statistical", "Statistical colour reconnection using simulated annealing", 1); static SwitchOption interfaceAlgorithmBaryonic (interfaceAlgorithm, "Baryonic", "Baryonic cluster reconnection", 2); static Parameter interfaceMaxDistance ("MaxDistance", "Maximum distance between the clusters at which to consider rearrangement" " to avoid colour reconneections of displaced vertices", &ColourReconnector::_maxDistance, femtometer, 1000.*femtometer, 0.0*femtometer, 1e100*femtometer, false, false, Interface::limited); static Switch interfaceOctetTreatment ("OctetTreatment", "Which octets are not allowed to be reconnected", &ColourReconnector::_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); } diff --git a/Hadronization/HadronSelector.cc b/Hadronization/HadronSelector.cc --- a/Hadronization/HadronSelector.cc +++ b/Hadronization/HadronSelector.cc @@ -1,1003 +1,1003 @@ // -*- C++ -*- // // HadronSelector.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the 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",""); +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 << ' ' // << 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; } } ostream & operator<< (ostream & os, const HadronSelector::HadronInfo & hi ) { os << std::scientific << std::showpoint << std::setprecision(4) << setw(2) << hi.id << '\t' // << hi.ptrData << ' ' << hi.swtef << '\t' << hi.wt << '\t' << hi.overallWeight << '\t' << ounit(hi.mass,GeV); return os; } HadronSelector::HadronSelector(unsigned int opt) : _pwtDquark( 1.0 ),_pwtUquark( 1.0 ),_pwtSquark( 1.0 ),_pwtCquark( 1.0 ), _pwtBquark( 1.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) { // 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 << _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::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 >> _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, 1.0, 0.0, 10.0, false,false,false); static Parameter interfacePwtBquark("PwtBquark","Weight for choosing a quark B", &HadronSelector::_pwtBquark, 0, 1.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); 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()) && 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(); 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()) 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 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; } 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; 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; // 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); // 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 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, tcPDPtr hadron) const { // This method receives in input three PDG ids, whose the // 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): // 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) // 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 ) { // 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) // 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 // 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 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)) { // 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() << " => " << 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() << "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] && 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] ) ) { return make_pair(vIdHad1[0],vIdHad2[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; for(unsigned int ix=0; ix<_partons.size(); ++ix) { if(!DiquarkMatcher::Check(_partons[ix]->id())) continue; 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,tcPDPtr ptr3) const { // 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 << " in _table. " << "In HadronSelector::lightestHadron()" << Exception::eventerror; if(tit->second.empty()) throw Exception() << "HadronSelector::lightestHadron " << "could not find any hadrons containing " << ptr1->id() << ' ' << ptr2->id() << '\n' << 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 ? 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 && candidate->mass() < ptr1->constituentMass() + ptr2->constituentMass()) { generator()->log() << "HadronSelector::lightestHadron: " << "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 > HadronSelector::hadronsBelowThreshold(Energy threshold, tcPDPtr ptr1, tcPDPtr ptr2, tcPDPtr ptr3) const { // 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 << " in _table. " << "In HadronSelector::hadronsBelowThreshold()" << Exception::eventerror; if(tit->second.empty()) throw Exception() << "HadronSelector::hadronsBelowThreshold() " << "could not find any hadrons containing " << ptr1->id() << ' ' << ptr2->id() << '\n' << 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 && candidate->mass() < ptr1->constituentMass() + ptr2->constituentMass()) { generator()->log() << "HadronSelector::hadronsBelowTheshold: " << "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 // 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)) threshold *= (1.0 + UseRandom::rnd()*_limExotic); else if (CheckId::hasBottom(par1,par2)) threshold *= (1.0 + UseRandom::rnd()*_limBottom); 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); } // new option select from those available else if(belowThreshold_ == 1) { 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 using namespace Herwig; DescribeNoPIOClass -describeHw64Selector("Herwig::Hw64Selector",""); +describeHw64Selector("Herwig::Hw64Selector","Herwig.so"); IBPtr Hw64Selector::clone() const { return new_ptr(*this); } IBPtr Hw64Selector::fullclone() const { return new_ptr(*this); } void Hw64Selector::Init() { static ClassDocumentation documentation ("There is no documentation for the Hw64Selector class"); } pair Hw64Selector::chooseHadronPair(const Energy cluMass,tcPDPtr par1, tcPDPtr par2,tcPDPtr) const { bool diquark = !(DiquarkMatcher::Check(par1->id()) || DiquarkMatcher::Check(par2->id())); pair lighthad = lightestHadronPair(par1, par2); if(!lighthad.first || !lighthad.second) throw Exception() << "Hw64Selector::chooseHadronPair " << "We have 0's! First id = " << par1->id() << " second = " << par2->id() << ". This is probably a problem with either" << " undecayed heavy particles or colour connections" << Exception::eventerror; // calculate maximum momentum Energy PCMax = Kinematics::pstarTwoBodyDecay(cluMass,lighthad.first->mass(), lighthad.second->mass()); tcPDPtr had1 = tcPDPtr(); tcPDPtr had2 = tcPDPtr(); int ntry = 0; tcPDPtr quark = tcPDPtr(); const int nmax = 5000; Energy p; do { quark = partons()[UseRandom::irnd(partons().size())]; if(diquark && DiquarkMatcher::Check(quark->id())) continue; if(pwt(quark->id()) <= UseRandom::rnd()) continue; pair pid(abs(par1->id()),quark->id()); KupcoData::const_iterator it1,it2; const HadronTable::const_iterator tit = table().find(pid); assert(tit != table().end()); const KupcoData & hdata = tit->second; do { it1 = hdata.begin(); advance(it1,int(hdata.size()*UseRandom::rnd())); } while(it1 != hdata.end() && it1->overallWeight < UseRandom::rnd()); had1 = it1->ptrData; pid = make_pair(quark->id(),abs(par2->id())); do { it2 = hdata.begin(); advance(it2,int(hdata.size()*UseRandom::rnd())); } while(it2 != hdata.end() && it2->overallWeight < UseRandom::rnd()); had2 = it2->ptrData; if(had1 && had2) { p = Kinematics::pstarTwoBodyDecay(cluMass, it1->mass, it2->mass); if(p/PCMax < UseRandom::rnd()) { had1 = had2 = tcPDPtr(); ntry++; } } } while((!had1|| !had2) && ntry < nmax); if(ntry >= nmax) return lighthad; int signHad1 = 0; int signHad2 = 0; if(CheckId::canBeHadron(par1,quark->CC()) && CheckId::canBeHadron(quark,par2)) { signHad1 = signHadron(par1, quark->CC(), had1); signHad2 = signHadron(par2, quark, had2); } else if(CheckId::canBeHadron(par1,quark) && CheckId::canBeHadron(quark->CC(),par2)) { signHad1 = signHadron(par1, quark, had1); signHad2 = signHadron(par2, quark->CC(), had2); } else throw Exception() << "Hw64Selector::chooseHadronPair()" << Exception::abortnow; return make_pair( signHad1 > 0 ? had1 : tcPDPtr(had1->CC()), signHad2 > 0 ? had2 : tcPDPtr(had2->CC())); } diff --git a/Hadronization/HwppSelector.cc b/Hadronization/HwppSelector.cc --- a/Hadronization/HwppSelector.cc +++ b/Hadronization/HwppSelector.cc @@ -1,253 +1,253 @@ // -*- C++ -*- // // HwppSelector.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the 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",""); +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; // 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; // calculate the weight double pwtstrange; if (quarktopick->id() == 3){ // Strangeness weight takes the automatic flat weight pwtstrange = pwt(3); // Scaling strangeness enhancement if (_enhanceSProb == 1){ double scale = double(sqr(_m0Decay/cluMass)); pwtstrange = (_maxScale < scale) ? 0. : pow(pwtstrange,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); pwtstrange = (_maxScale < scale) ? 0. : exp(-scale); } weight = pwtstrange * H1->overallWeight * H2->overallWeight * Kinematics::pstarTwoBodyDecay(cluMass, H1->mass, H2->mass ); } else { weight = pwt(quarktopick->id()) * H1->overallWeight * H2->overallWeight * Kinematics::pstarTwoBodyDecay(cluMass, H1->mass, H2->mass ); } int signQ = 0; assert (par1 && quarktopick); assert (par2); assert(quarktopick->CC()); if(CheckId::canBeHadron(par1, quarktopick->CC()) && CheckId::canBeHadron(quarktopick, par2)) signQ = +1; else if(CheckId::canBeHadron(par1, quarktopick) && CheckId::canBeHadron(quarktopick->CC(), par2)) signQ = -1; else { cerr << "Could not make sign for" << par1->id()<< " " << quarktopick->id() << " " << par2->id() << "\n"; assert(false); } if (signQ == -1) quarktopick = quarktopick->CC(); // construct the object with the info Kupco a(quarktopick, H1->ptrData, H2->ptrData, weight); hadrons.push_back(a); wgtsum += weight; } } } if (hadrons.empty()) return make_pair(tcPDPtr(),tcPDPtr()); // select the hadron wgtsum *= UseRandom::rnd(); unsigned int ix=0; do { wgtsum-= hadrons[ix].weight; ++ix; } while(wgtsum > ZERO && ix < hadrons.size()); if(ix == hadrons.size() && wgtsum > ZERO) return make_pair(tcPDPtr(),tcPDPtr()); --ix; assert(hadrons[ix].idQ); int signHad1 = signHadron(par1, hadrons[ix].idQ->CC(), hadrons[ix].hadron1); int signHad2 = signHadron(par2, hadrons[ix].idQ, hadrons[ix].hadron2); assert( signHad1 != 0 && signHad2 != 0 ); return make_pair ( signHad1 > 0 ? hadrons[ix].hadron1 : tcPDPtr(hadrons[ix].hadron1->CC()), signHad2 > 0 ? hadrons[ix].hadron2 : tcPDPtr(hadrons[ix].hadron2->CC())); } diff --git a/Hadronization/LightClusterDecayer.cc b/Hadronization/LightClusterDecayer.cc --- a/Hadronization/LightClusterDecayer.cc +++ b/Hadronization/LightClusterDecayer.cc @@ -1,386 +1,386 @@ // -*- C++ -*- // // LightClusterDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the LightClusterDecayer class. // #include "LightClusterDecayer.h" #include #include #include #include #include #include #include #include #include "Cluster.h" #include "CheckId.h" #include "Herwig/Utilities/Kinematics.h" #include using namespace Herwig; DescribeClass -describeLightClusterDecayer("Herwig::LightClusterDecayer",""); +describeLightClusterDecayer("Herwig::LightClusterDecayer","Herwig.so"); IBPtr LightClusterDecayer::clone() const { return new_ptr(*this); } IBPtr LightClusterDecayer::fullclone() const { return new_ptr(*this); } void LightClusterDecayer::persistentOutput(PersistentOStream & os) const { os << _hadronSelector; } void LightClusterDecayer::persistentInput(PersistentIStream & is, int) { is >> _hadronSelector; } void LightClusterDecayer::Init() { static ClassDocumentation documentation ("There is the class responsible for the one-hadron decay of light clusters"); static Reference interfaceHadronSelector("HadronSelector", "A reference to the HadronSelector object", &Herwig::LightClusterDecayer::_hadronSelector, false, false, true, false); } bool LightClusterDecayer::decay(ClusterVector & clusters, tPVector & finalhadrons) { // Loop over all clusters, and for those that were not heavy enough // to undergo to fission, check if they are below the threshold // for normal two-hadron decays. If this is the case, then the cluster // should be decayed into a single hadron: this can happen only if // it is possible to reshuffle momenta between the cluster and // another one; in the rare occasions in which such exchange of momenta // is not possible (because all of the clusters are too light) then // the event is skipped. // Notice that, differently from what happens in Fortran Herwig, // light (that is below the threshold for the production of the lightest // pair of hadrons with the proper flavours) fission products, produced // by the fission of heavy clusters in class ClusterFissioner // have been already "decayed" into single hadron (the lightest one // with proper flavour) by the same latter class, without requiring // any reshuffling. Therefore the light clusters that are treated in // this LightClusterDecayer class are produced directly // (originally) by the class ClusterFinder. // To preserve all of the information, the cluster partner with which // the light cluster (that decays into a single hadron) exchanges // momentum in the reshuffling procedure is redefined and inserted // in the vector vecNewRedefinedCluPtr. Only at the end, when all // light clusters have been examined, the elements this vector will be // copied in collecCluPtr (the reason is that it is not allowed to // modify a STL container while iterating over it. At the same time, // this ensures that a cluster can be redefined only once, which seems // sensible although not strictly necessary). // Notice that the cluster reshuffling partner is normally redefined // and inserted in the vector vecNewRedefinedCluPtr, but not always: // in the case it is also light, then it is also decayed immediately // into a single hadron, without redefining it (the reason being that, // otherwise, the would-be redefined cluster could have undefined // components). vector redefinedClusters; for (ClusterVector::const_iterator it = clusters.begin(); it != clusters.end(); ++it) { // Skip the clusters that are not available or that are // heavy, intermediate, clusters that have undergone to fission, if ( ! (*it)->isAvailable() || ! (*it)->isReadyToDecay() ) continue; // We need to require (at least at the moment, maybe in the future we // could change it) that the cluster has exactly two components, // because otherwise we don't know how to deal with the kinematics. // If this is not the case, then send a warning because it is not suppose // to happen, and then do nothing with (ignore) such cluster. if ( (*it)->numComponents() != 2 ) { generator()->logWarning( Exception("LightClusterDecayer::decay " "***Still cluster with not exactly" " 2 components*** ", Exception::warning) ); continue; } // select the hadron for single hadron decay tcPDPtr hadron = _hadronSelector->chooseSingleHadron((*it)->particle(0)->dataPtr(), (*it)->particle(1)->dataPtr(), (**it).mass()); // if not single decay continue if(!hadron) continue; // We assume that the candidate reshuffling cluster partner, // with whom the light cluster can exchange momenta, // is chosen as the closest in space-time between the available // clusters. Notice that an alternative, sensible approach // could be to consider instead the "closeness" in the colour // structure... // Notice that nor a light cluster (which decays into a single hadron) // neither its cluster reshuffling partner (which either has a // redefined cluster or also decays into a single hadron) can be // a reshuffling partner of another light cluster. // This because we are requiring that the considered candidate cluster // reshuffling partner has the status "isAvailable && isReadyToDecay" true; // furthermore, the new redefined clusters are not added to the collection // of cluster before the end of the entire reshuffling procedure, avoiding // in this way that the redefined cluster of a cluster reshuffling partner // is used again later. Needless to say, this is just an assumption, // although reasonable, but nothing more than that! // Build a multimap of available reshuffling cluster partners, // with key given by the module of the invariant space-time distance // w.r.t. the light cluster, so that this new collection is automatically // ordered in increasing distance values. // We use a multimap, rather than a map, just for precaution against not properly // defined cluster positions which could produce all identical (null) distances. multimap candidates; for ( ClusterVector::iterator jt = clusters.begin(); jt != clusters.end(); ++jt ) { if ((*jt)->isAvailable() && (*jt)->isReadyToDecay() && jt != it) { Length distance = abs (((*it)->vertex() - (*jt)->vertex()).m()); candidates.insert(pair(distance,*jt)); } } // Loop sequentially the multimap. multimap::const_iterator mmapIt = candidates.begin(); bool found = false; while (!found && mmapIt != candidates.end()) { found = reshuffling(hadron, *it, (*mmapIt).second, redefinedClusters, finalhadrons); if (!found) ++mmapIt; } if (!found) return partonicReshuffle(hadron,*it,finalhadrons); } // end loop over collecCluPtr // Add to collecCluPtr all of the redefined new clusters (indeed the // pointers to them are added) contained in vecNewRedefinedCluPtr. for (tClusterVector::const_iterator it = redefinedClusters.begin(); it != redefinedClusters.end(); ++it) { clusters.push_back(*it); } return true; } bool LightClusterDecayer::reshuffling(const tcPDPtr pdata1, tClusterPtr cluPtr1, tClusterPtr cluPtr2, tClusterVector & redefinedClusters, tPVector & finalhadrons) { // don't reshuffle with beam clusters if(cluPtr2->isBeamCluster()) return false; // This method does the reshuffling of momenta between the cluster "1", // that must decay into a single hadron (with id equal to idhad1), and // the candidate cluster "2". It returns true if the reshuffling succeed, // false otherwise. PPtr ptrhad1 = pdata1->produceParticle(); if ( ! ptrhad1 ) { generator()->logWarning( Exception("LightClusterDecayer::reshuffling" "***Cannot create a particle with specified id***", Exception::warning) ); return false; } Energy mhad1 = ptrhad1->mass(); // Let's call "3" and "4" the two constituents of the second cluster tPPtr part3 = cluPtr2->particle(0); tPPtr part4 = cluPtr2->particle(1); // Check if the system of the two clusters can kinematically be replaced by // an hadron of mass mhad1 (which is the lightest single hadron with the // same flavour numbers as the first cluster) and the second cluster. // If not, then try to replace the second cluster with the lightest hadron // with the same flavour numbers; if it still fails, then give up! Lorentz5Momentum pSystem = cluPtr1->momentum() + cluPtr2->momentum(); pSystem.rescaleMass(); // set the mass as the invariant of the quadri-vector Energy mSystem = pSystem.mass(); Energy mclu2 = cluPtr2->mass(); bool singleHadron = false; Energy mLHP2 = _hadronSelector->massLightestHadronPair(part3->dataPtr(),part4->dataPtr()); Energy mLH2 = _hadronSelector->massLightestHadron(part3->dataPtr(),part4->dataPtr()); if(mSystem > mhad1 + mclu2 && mclu2 > mLHP2) { singleHadron = false; } else if(mSystem > mhad1 + mLH2) { singleHadron = true; mclu2 = mLH2; } else return false; // Let's call from now on "Sys" the system of the two clusters, and // had1 (of mass mhad1) the lightest hadron in which the first // cluster decays, and clu2 (of mass mclu2) either the second // cluster or the lightest hadron in which it decays (depending // which one is kinematically allowed, see above). // The idea behind the reshuffling is to replace the system of the // two clusters by the system of the hadron had1 and (cluster or hadron) clu2, // but leaving the overall system unchanged. Furthermore, the motion // of had1 and clu2 in the Sys frame is assumed to be parallel to, respectively, // those of the original cluster1 and cluster2 in the same Sys frame. // Calculate the unit three-vector, in the frame "Sys" along which the // two initial clusters move. Lorentz5Momentum u( cluPtr1->momentum() ); u.boost( - pSystem.boostVector() ); // boost from LAB to Sys // Calculate the momenta of had1 and clu2 in the Sys frame first, // and then boost back in the LAB frame. Lorentz5Momentum phad1, pclu2; if (pSystem.m() < mhad1 + mclu2 ) { throw Exception() << "Impossible Kinematics in LightClusterDecayer::reshuffling()" << Exception::eventerror; } Kinematics::twoBodyDecay(pSystem, mhad1, mclu2, u.vect().unit(), phad1, pclu2); ptrhad1->set5Momentum( phad1 ); // set momentum of first hadron. ptrhad1->setVertex(cluPtr1->vertex()); // set hadron vertex position to the // parent cluster position. cluPtr1->addChild(ptrhad1); finalhadrons.push_back(ptrhad1); cluPtr1->flagAsReshuffled(); cluPtr2->flagAsReshuffled(); if(singleHadron) { // In the case that also the cluster reshuffling partner is light // it is decayed into a single hadron, *without* creating the // redefined cluster (this choice is justified in order to avoid // clusters that could have undefined components). PPtr ptrhad2 = _hadronSelector->lightestHadron(part3->dataPtr(),part4->dataPtr()) ->produceParticle(); ptrhad2->set5Momentum( pclu2 ); ptrhad2->setVertex( cluPtr2->vertex() ); // set hadron vertex position to the // parent cluster position. cluPtr2->addChild(ptrhad2); finalhadrons.push_back(ptrhad2); } else { // Create the new cluster which is the redefinitions of the cluster // partner (cluster "2") used in the reshuffling procedure of the // light cluster (cluster "1"). // The rationale of this is to preserve completely all of the information. ClusterPtr cluPtr2new = ClusterPtr(); if(part3 && part4) cluPtr2new = new_ptr(Cluster(part3,part4)); cluPtr2new->set5Momentum( pclu2 ); cluPtr2new->setVertex( cluPtr2->vertex() ); cluPtr2->addChild( cluPtr2new ); redefinedClusters.push_back( cluPtr2new ); // Set consistently the momenta of the two components of the second cluster // after the reshuffling. To do that we first calculate the momenta of the // constituents in the initial cluster rest frame; then we boost them back // in the lab but using this time the new cluster rest frame. Finally we store // these information in the new cluster. Notice that we do *not* set // consistently also the momenta of the (eventual) particles pointed by the // two components: that's because we do not need to do so, being the momentum // an explicit private member of the class Component (which is set equal // to the momentum of the eventual particle pointed only in the constructor, // but then later should not necessary be the same), and furthermore it allows // us not to loose any information, in the sense that we can always, later on, // to find the original momenta of the two components before the reshuffling. Lorentz5Momentum p3 = part3->momentum(); //p3new->momentum(); p3.boost( - (cluPtr2->momentum()).boostVector() ); // from LAB to clu2 (old) frame p3.boost( pclu2.boostVector() ); // from clu2 (new) to LAB frame Lorentz5Momentum p4 = part4->momentum(); //p4new->momentum(); p4.boost( - (cluPtr2->momentum()).boostVector() ); // from LAB to clu2 (old) frame p4.boost( pclu2.boostVector() ); // from clu2 (new) to LAB frame cluPtr2new->particle(0)->set5Momentum(p3); cluPtr2new->particle(1)->set5Momentum(p4); } // end of if (singleHadron) return true; } bool LightClusterDecayer::partonicReshuffle(const tcPDPtr had, const PPtr cluster, tPVector & finalhadrons) { tPPtr meson(cluster); if(!meson->parents().empty()) meson=meson->parents()[0]; if(!meson->parents().empty()) meson=meson->parents()[0]; // check b/c hadron decay int ptype(abs(meson->id())%10000); bool heavy = (ptype/1000 == 5 || ptype/1000 ==4 ); heavy |= (ptype/100 == 5 || ptype/100 ==4 ); heavy |= (ptype/10 == 5 || ptype/10 ==4 ); if(!heavy) return false; // find the leptons tPVector leptons; for(unsigned int ix=0;ixchildren().size();++ix) { if(!(meson->children()[ix]->dataPtr()->coloured())) { leptons.push_back(meson->children()[ix]); } } if(leptons.size()==1) { tPPtr w=leptons[0]; leptons.pop_back(); for(unsigned int ix=0;ixchildren().size();++ix) { if(!w->children()[ix]->dataPtr()->coloured()) { leptons.push_back(w->children()[ix]); } } } if(leptons.size()!=2) return false; // get momentum of leptonic system and the its minimum possible mass Energy mmin(ZERO); Lorentz5Momentum pw; for(unsigned int ix=0;ixmomentum(); mmin+=leptons[ix]->mass(); } pw.rescaleMass(); // check we can do the reshuffling PPtr ptrhad = had->produceParticle(); // total momentum fo the system Lorentz5Momentum pSystem = pw + cluster->momentum(); pSystem.rescaleMass(); // normal case get additional energy by rescaling momentum in rest frame of // system if(pSystem.mass()>ptrhad->mass()+pw.mass()&&pw.mass()>mmin) { // Calculate the unit three-vector, in the frame "Sys" along which the // two initial clusters move. Lorentz5Momentum u(cluster->momentum()); u.boost( - pSystem.boostVector() ); // Calculate the momenta of had1 and clu2 in the Sys frame first, // and then boost back in the LAB frame. Lorentz5Momentum phad1, pclu2; Kinematics::twoBodyDecay(pSystem, ptrhad->mass(), pw.mass(), u.vect().unit(), phad1, pclu2); // set momentum of first hadron. ptrhad->set5Momentum( phad1 ); // set hadron vertex position to the parent cluster position. ptrhad->setLabVertex(cluster->vertex()); // add hadron cluster->addChild(ptrhad); finalhadrons.push_back(ptrhad); // reshuffle the leptons // boost the leptons to the rest frame of the system Boost boost1(-pw.boostVector()); Boost boost2( pclu2.boostVector()); for(unsigned int ix=0;ixdeepBoost(boost1); leptons[ix]->deepBoost(boost2); } return true; } else { return false; } } diff --git a/Hadronization/PartonSplitter.cc b/Hadronization/PartonSplitter.cc --- a/Hadronization/PartonSplitter.cc +++ b/Hadronization/PartonSplitter.cc @@ -1,555 +1,555 @@ // -*- C++ -*- // // PartonSplitter.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the PartonSplitter class. // #include "PartonSplitter.h" #include #include #include #include #include #include #include #include "ThePEG/Interface/Parameter.h" #include #include #include "ThePEG/Repository/UseRandom.h" #include "Herwig/Utilities/Kinematics.h" #include #include "ClusterHadronizationHandler.h" #include #include #include "CheckId.h" using namespace Herwig; IBPtr PartonSplitter::clone() const { return new_ptr(*this); } IBPtr PartonSplitter::fullclone() const { return new_ptr(*this); } void PartonSplitter::persistentOutput(PersistentOStream & os) const { os << _quarkSelector << ounit(_gluonDistance,femtometer) << _splitGluon << _splitPwtUquark << _splitPwtDquark << _splitPwtSquark << _enhanceSProb << ounit(_m0,GeV) << _massMeasure; } void PartonSplitter::persistentInput(PersistentIStream & is, int) { is >> _quarkSelector >> iunit(_gluonDistance,femtometer) >> _splitGluon >> _splitPwtUquark >> _splitPwtDquark >> _splitPwtSquark >>_enhanceSProb >> iunit(_m0,GeV) >> _massMeasure; } DescribeClass -describePartonSplitter("Herwig::PartonSplitter",""); +describePartonSplitter("Herwig::PartonSplitter","Herwig.so"); void PartonSplitter::Init() { static ClassDocumentation documentation ("This class is reponsible of the nonperturbative splitting of partons"); static Switch interfaceSplit ("Split", "Option for different splitting options", &PartonSplitter::_splitGluon, 0, false, false); static SwitchOption interfaceSplitDefault (interfaceSplit, "ud", "Normal cluster splitting where only u and d quarks are drawn is used.", 0); static SwitchOption interfaceSplitAll (interfaceSplit, "uds", "Alternative cluster splitting where all light quark pairs (u, d, s) can be drawn.", 1); static Parameter interfaceSplitPwtUquark ("SplitPwtUquark", "Weight for splitting in U quarks", &PartonSplitter::_splitPwtUquark, 1, 0.0, 1.0, false, false, Interface::limited); static Parameter interfaceSplitPwtDquark ("SplitPwtDquark", "Weight for splitting in D quarks", &PartonSplitter::_splitPwtDquark, 1, 0.0, 1.0, false, false, Interface::limited); static Parameter interfaceSplitPwtSquark ("SplitPwtSquark", "Weight for splitting in S quarks", &PartonSplitter::_splitPwtSquark, 0.5, 0.0, 1.0, false, false, Interface::limited); static Switch interfaceEnhanceSProb ("EnhanceSProb", "Option to enhance the strangeness weight (MassSplit Switch needs to be on)", &PartonSplitter::_enhanceSProb,0,false,false); static SwitchOption interfaceEnhanceSProbNo (interfaceEnhanceSProb, "No", "No scaling for strangeness", 0); static SwitchOption interfaceEnhanceSProbScale (interfaceEnhanceSProb, "Scaled", "Power-law scaling for strangeness", 1); static SwitchOption interfaceEnhanceSProbExp (interfaceEnhanceSProb, "Exponential", "Exponential suppression for strangeness", 2); static Switch interfaceMassMeasure ("MassMeasure", "Option to use different mass measures", &PartonSplitter::_massMeasure,0,false,false); static SwitchOption interfaceMassMeasureMass (interfaceMassMeasure, "Mass", "Mass Measure", 0); static SwitchOption interfaceMassMeasureLambda (interfaceMassMeasure, "Lambda", "Lambda Measure", 1); static Parameter interfaceMassScale ("MassScale", "Mass scale for g->qqb strangeness enhancement", &PartonSplitter::_m0, GeV, 20.*GeV, 1.*GeV, 1000.*GeV, false, false, Interface::limited); } void PartonSplitter::split(PVector & tagged) { // set the gluon c tau once and for all static bool first = true; if(first) { _gluonDistance = hbarc*getParticleData(ParticleID::g)->constituentMass()/ ClusterHadronizationHandler::currentHandler()->minVirtuality2(); first = false; } // Copy incoming for the (possible sorting and) splitting PVector particles = tagged; // Switch to enhance strangeness if (_enhanceSProb >= 1){ colourSorted(particles); } PVector newtag; Energy2 Q02 = 0.99*sqr(getParticleData(ParticleID::g)->constituentMass()); // Loop over all of the particles in the event. // Loop counter for colourSorted for(PVector::iterator pit = particles.begin(); pit!=particles.end(); ++pit) { // only considering gluons so add other particles to list of particles if( (**pit).data().id() != ParticleID::g ) { newtag.push_back(*pit); continue; } // should not have been called for massless or space-like gluons if((**pit).momentum().m2() <= 0.0*sqr(MeV) ) { throw Exception() << "Spacelike or massless gluon m2= " << (**pit).momentum().m2()/GeV2 << "GeV2 in PartonSplitter::split()" << Exception::eventerror; } // time like gluon gets split PPtr ptrQ = PPtr(); PPtr ptrQbar = PPtr(); if (_enhanceSProb == 0){ splitTimeLikeGluon(*pit,ptrQ,ptrQbar); } else { size_t i = pit - particles.begin(); massSplitTimeLikeGluon(*pit, ptrQ, ptrQbar, i); } ptrQ->scale(Q02); ptrQbar->scale(Q02); (*pit)->colourLine()->addColoured(ptrQ); (*pit)->addChild(ptrQ); newtag.push_back(ptrQ); (*pit)->antiColourLine()->addAntiColoured(ptrQbar); (*pit)->addChild(ptrQbar); newtag.push_back(ptrQbar); // set the life length of gluon Length distance = UseRandom::rndExp(_gluonDistance); (**pit).setLifeLength((distance/(**pit).mass())*(**pit).momentum()); // assume quarks same position as gluon ptrQ ->setVertex((**pit).decayVertex()); ptrQ ->setLifeLength(Lorentz5Distance()); ptrQbar->setVertex((**pit).decayVertex()); ptrQbar->setLifeLength(Lorentz5Distance()); } swap(tagged,newtag); } void PartonSplitter::splitTimeLikeGluon(tcPPtr ptrGluon, PPtr & ptrQ, PPtr & ptrQbar){ // select the quark flavour tPDPtr quark; long idNew=0; switch(_splitGluon){ case 0: quark = _quarkSelector.select(UseRandom::rnd()); break; case 1: if ( ptrGluon->momentum().m() < 2.0 *getParticle(ThePEG::ParticleID::s)->data().constituentMass() ) { throw Exception() << "Impossible Kinematics in PartonSplitter::splitTimeLikeGluon()" << Exception::runerror; } // Only allow light quarks u,d,s with the probabilities double prob_d = _splitPwtDquark; double prob_u = _splitPwtUquark; double prob_s = _splitPwtSquark; int choice = UseRandom::rnd3(prob_u, prob_d, prob_s); switch(choice) { case 0: idNew = ThePEG::ParticleID::u; break; case 1: idNew = ThePEG::ParticleID::d; break; case 2: idNew = ThePEG::ParticleID::s; break; } ptrQ = getParticle(idNew); ptrQbar = getParticle(-idNew); break; } // Solve the kinematics of the two body decay G --> Q + Qbar Lorentz5Momentum momentumQ; Lorentz5Momentum momentumQbar; double cosThetaStar = UseRandom::rnd( -1.0 , 1.0 ); using Constants::pi; double phiStar = UseRandom::rnd( -pi , pi ); Energy constituentQmass; if(_splitGluon==0) { constituentQmass = quark->constituentMass(); }else{ constituentQmass = ptrQ->data().constituentMass(); } if (ptrGluon->momentum().m() < 2.0*constituentQmass) { throw Exception() << "Impossible Kinematics in PartonSplitter::splitTimeLikeGluon()" << Exception::eventerror; } Kinematics::twoBodyDecay(ptrGluon->momentum(), constituentQmass, constituentQmass, cosThetaStar, phiStar, momentumQ, momentumQbar ); // Create quark and anti-quark particles of the chosen flavour // and set they 5-momentum (the mass is the constituent one). if(_splitGluon==0) { ptrQ = new_ptr(Particle(quark )); ptrQbar = new_ptr(Particle(quark->CC())); } ptrQ ->set5Momentum( momentumQ ); ptrQbar ->set5Momentum( momentumQbar ); } void PartonSplitter::doinit() { Interfaced::doinit(); // calculate the probabilties for the gluon to branch into each quark type // based on the available phase-space, as in fortran. Energy mg=getParticleData(ParticleID::g)->constituentMass(); for( int ix=1; ix<6; ++ix ) { PDPtr quark = getParticleData(ix); Energy pcm = Kinematics::pstarTwoBodyDecay(mg,quark->constituentMass(), quark->constituentMass()); if(pcm>ZERO) _quarkSelector.insert(pcm/GeV,quark); } if(_quarkSelector.empty()) throw InitException() << "At least one quark must have constituent mass less " << "then the constituent mass of the gluon in " << "PartonSplitter::doinit()" << Exception::runerror; } // Method to colour sort the event and calculate the masses of the // pre-clusters // Convention is to have the triplet of the colour singlet first, // then all gluons, then the antitriplet (and repeat for all the // colour-singlets in the event) void PartonSplitter::colourSorted(PVector& tagged) { // Set up the output PVector sorted; // Reset the storage variables for doing the mass-based strangeness // enhancement _colSingletSize.resize(0); _colSingletm2.resize(0); // Variable to exit out of gluon loops bool gluonLoop = false; // Partons left to consider PVector left = tagged; // Loop over all the tagged particles while (int(left.size()) > 0){ // Pick the first particle available PPtr p = left[0]; // Temporary holding pen for momenta Lorentz5Momentum tempMom(ZERO, ZERO, ZERO, ZERO); // Don't do anything if the particle has no colour // Simply add it back into the list of particles if ( !p->coloured() ) { sorted.push_back(p); // nparts is the index of the particle after adding it to the sorted list int nparts = 1; // Add on the last entry of colSingletSize if the vector is not empty // This is essentially the index the particle will have once it // Get placed into the new colour sorted event nparts += (_colSingletSize.empty()) ? 0 : _colSingletSize.back(); tempMom += p->momentum(); Energy2 singletm2 = tempMom.m2(); // Store the number of particles and mass of the colour-singlet _colSingletSize.push_back(nparts); _colSingletm2.push_back(singletm2); // Remove the particle from the list of particles left left.erase(remove(left.begin(),left.end(), p), left.end()); } // Temporary holding pen for partons PVector temp; // Variable to sum end-point masses i.e. triplets and anti-triplets Energy endPointMass = ZERO; // If the particle in question is a gluon, search up it's antiColourLine // until we get to the triplet. // Note there are situations where we have a gluon loop if ( p->hasColour() && p->hasAntiColour() ){ // Search up its anticolour line until we get to the start particle tPPtr partner = p->antiColourLine()->endParticle(); // Dummy index used to loop over the anticolour line trace tPPtr dummy = partner; // Store the partner particle temp.push_back(partner); tempMom += partner->momentum(); left.erase(remove(left.begin(),left.end(), partner), left.end()); // While loop continues while we still reach a particle with with // anti-colour, i.e. a gluon while ( dummy->hasAntiColour() ){ dummy = dummy->antiColourLine()->endParticle(); // Check that we haven't already included it via colour indices // If we have, it is a gluon loop. if ( find(left.begin(), left.end(), dummy) == left.end() ) { gluonLoop = true; break; } // Store the dummy partons in reverse temp.push_back(dummy); tempMom += dummy->momentum(); // Remove counted ones from the remaining list left.erase(remove(left.begin(),left.end(), dummy), left.end()); } // Number of particles in this colour singlets so far int nparts = int(temp.size()); // Insert the new particles in the reverse order sorted.insert(sorted.end(), temp.rbegin(), temp.rend()); endPointMass += ((temp.back())->mass()); // If it is a gluon loop we've already looped over the colour-singlet // in its entirety, so we need to end early if (gluonLoop){ // Store the index of the final entry nparts += (_colSingletSize.empty()) ? 0 : _colSingletSize.back(); // Insert the new particles in the correct order // i.e. triplet first, then the gluons we have seen so far Energy2 singletm2 = tempMom.m2(); _colSingletSize.push_back(nparts); _colSingletm2.push_back(singletm2); continue; } // If it is not a gluon loop, we now need to trace the colourLine // down, until we reach the triplet. // Works similarly to the first half // Reset the temp PVector temp.resize(0); // Push back the particle in question temp.push_back(p); // tempMom hasn't been reset, add the particle we started with tempMom += p->momentum(); left.erase(remove(left.begin(),left.end(), p), left.end()); // Search down its colour line until we get to the end particle dummy = p->colourLine()->startParticle(); temp.push_back(dummy); tempMom += dummy->momentum(); left.erase(remove(left.begin(),left.end(), dummy), left.end()); while ( dummy->hasColour() ){ dummy = dummy->colourLine()->startParticle(); temp.push_back(dummy); tempMom += dummy->momentum(); left.erase(remove(left.begin(),left.end(), dummy), left.end()); } endPointMass += ((temp.back())->mass()); // Update size of colour singlet nparts += int(temp.size()); nparts += (_colSingletSize.empty()) ? 0 : _colSingletSize.back(); // Insert the new particles in the correct order Energy2 singletm2 = tempMom.m2(); sorted.insert(sorted.end(), temp.begin(), temp.end()); endPointMass += ((sorted.back())->mass()); _colSingletSize.push_back(nparts); // Chooses whether to use invariant mass of singlet // or to use the lambda measure i.e. m^2 - (\sum m_i)^2 Energy2 m2 = (_massMeasure == 0) ? singletm2 : singletm2 - sqr(endPointMass); _colSingletm2.push_back(m2); } // Else if it's a quark else if ( p->hasColour() ) { // Search up its colour line until we get to the start particle // Works the same way as the second half of the gluon handling tPPtr partner = p->colourLine()->startParticle(); tPPtr dummy = partner; temp.push_back(p); endPointMass += ((temp.back())->mass()); temp.push_back(partner); tempMom += p->momentum(); tempMom += partner->momentum(); left.erase(remove(left.begin(),left.end(), p), left.end()); left.erase(remove(left.begin(),left.end(), partner), left.end()); while ( dummy->hasColour() ){ dummy = dummy->colourLine()->startParticle(); temp.push_back(dummy); tempMom += dummy->momentum(); left.erase(remove(left.begin(),left.end(), dummy), left.end()); } // Number of particles in this colour singlets int nparts = int(temp.size()); nparts += (_colSingletSize.empty()) ? 0 : _colSingletSize.back(); // Insert the new particles in the correct order Energy2 singletm2 = tempMom.m2(); sorted.insert(sorted.end(), temp.begin(), temp.end()); endPointMass += ((sorted.back())->mass()); _colSingletSize.push_back(nparts); Energy2 m2 = (_massMeasure == 0) ? singletm2 : singletm2 - sqr(endPointMass); _colSingletm2.push_back(m2); } // Else it's an antiquark else if ( p->hasAntiColour() ) { // Search along anti-colour line, storing particles, and reversing the order // at the end // Works in the same way as the first half of the gluon handling tPPtr partner = p->antiColourLine()->endParticle(); tPPtr dummy = partner; temp.push_back(p); endPointMass += ((temp.back())->mass()); temp.push_back(partner); tempMom += p->momentum(); tempMom += partner->momentum(); left.erase(remove(left.begin(),left.end(), p), left.end()); left.erase(remove(left.begin(),left.end(), partner), left.end()); while ( dummy->hasAntiColour() ){ dummy = dummy->antiColourLine()->endParticle(); temp.push_back(dummy); tempMom += dummy->momentum(); left.erase(remove(left.begin(),left.end(), dummy), left.end()); } // Number of particles in this colour singlets int nparts = int(temp.size()); nparts += (_colSingletSize.empty()) ? 0 : _colSingletSize.back(); // Insert the particles in the reverse order Energy2 singletm2 = tempMom.m2(); sorted.insert(sorted.end(), temp.rbegin(), temp.rend()); endPointMass += ((temp.back())->mass()); _colSingletSize.push_back(nparts); Energy2 m2 = (_massMeasure == 0) ? singletm2 : singletm2 - sqr(endPointMass); _colSingletm2.push_back(m2); } } // Check that the algorithm hasn't missed any particles. assert( sorted.size() == tagged.size() ); swap(sorted,tagged); } void PartonSplitter::massSplitTimeLikeGluon(tcPPtr ptrGluon, PPtr & ptrQ, PPtr & ptrQbar, size_t i){ // select the quark flavour tPDPtr quark; long idNew=0; if ( ptrGluon->momentum().m() < 2.0 *getParticle(ThePEG::ParticleID::s)->data().constituentMass() ) { throw Exception() << "Impossible Kinematics in PartonSplitter::massSplitTimeLikeGluon()" << Exception::runerror; } // Only allow light quarks u,d,s with the probabilities double prob_d = _splitPwtDquark; double prob_u = _splitPwtUquark; double prob_s = enhanceStrange(i); int choice = UseRandom::rnd3(prob_u, prob_d, prob_s); switch(choice) { case 0: idNew = ThePEG::ParticleID::u; break; case 1: idNew = ThePEG::ParticleID::d; break; case 2: idNew = ThePEG::ParticleID::s; break; } ptrQ = getParticle(idNew); ptrQbar = getParticle(-idNew); // Solve the kinematics of the two body decay G --> Q + Qbar Lorentz5Momentum momentumQ; Lorentz5Momentum momentumQbar; double cosThetaStar = UseRandom::rnd( -1.0 , 1.0 ); using Constants::pi; double phiStar = UseRandom::rnd( -pi , pi ); Energy constituentQmass; constituentQmass = ptrQ->data().constituentMass(); if (ptrGluon->momentum().m() < 2.0*constituentQmass) { throw Exception() << "Impossible Kinematics in PartonSplitter::massSplitTimeLikeGluon()" << Exception::eventerror; } Kinematics::twoBodyDecay(ptrGluon->momentum(), constituentQmass, constituentQmass, cosThetaStar, phiStar, momentumQ, momentumQbar ); // Create quark and anti-quark particles of the chosen flavour // and set they 5-momentum (the mass is the constituent one). ptrQ ->set5Momentum( momentumQ ); ptrQbar ->set5Momentum( momentumQbar ); } double PartonSplitter::enhanceStrange(size_t i){ // Get the m2 of the relevant colour singlet // First we need to get the index of the colour-singlet the indexed i // parton has auto const it = lower_bound(_colSingletSize.begin(), _colSingletSize.end(), i); // Get the index of the colourSinglet mass int indx = distance(_colSingletSize.begin(), it); Energy2 mass2 = _colSingletm2[indx]; Energy2 m2 = _m0*_m0; // Scaling strangeness enhancement if (_enhanceSProb == 1){ // If the mass is significantly smaller than the characteristic mass, // just set the prob to 0 double scale = double(m2/mass2); return (_maxScale < scale || scale < 0.) ? 0. : pow(_splitPwtSquark,scale); } // Exponential strangeness enhancement else if (_enhanceSProb == 2){ double scale = double(m2/mass2); return (_maxScale < scale || scale < 0.) ? 0. : exp(-scale); } else return _splitPwtSquark; }