diff --git a/Hadronization/Cluster.cc b/Hadronization/Cluster.cc --- a/Hadronization/Cluster.cc +++ b/Hadronization/Cluster.cc @@ -1,271 +1,275 @@ // -*- C++ -*- // // Cluster.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the 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(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(unsigned int i = 0; i<_numComp; i++) if(_isBeamRemnant[i]) return true; return false; } void Cluster::isBeamCluster(tPPtr part) { 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(); + return (i < (int)_numComp) ? _component[i] : PPtr(); +} + +tPPtr Cluster::particleB(int i) const{ + return (i < (int)_numComp) ? _original[i] : tPPtr(); } 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,346 @@ // -*- C++ -*- // // Cluster.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_Cluster_H #define HERWIG_Cluster_H #include #include "Herwig/Utilities/EnumParticles.h" #include "CluHadConfig.h" #include "ClusterHadronizationHandler.fh" #include "Cluster.fh" #include "ThePEG/Utilities/ClassDescription.h" namespace Herwig { using namespace ThePEG; /** \ingroup Hadronization * \class Cluster * \brief This class describes a cluster object. * \author Philip Stephens * \author Alberto Ribon * * This class represents a cluster, which is a colour singlet made usually * of two components (quark-antiquark, quark-diquark, antiquark-antidiquark) * or rarely by three components (quark-quark-quark, antiquark-antiquark- * antiquark). A reference to the container with the pointers to its * Components is provided. * * The class provides access to the pointers which point to: * * - The cluster parent. In the case that the cluster it is a fission * product of a heavy cluster the parent is a cluster. If the cluster * is formed from the perturbative partons then the parents will be * the colour connected partons that formed the cluster. * - The children (usually two). In the case the cluster is a * heavy cluster that undergoes fission the children are clusters. * Occasionally the cluster has been "redefined" (re-interpreted). For * example in the case that three quark or anti-quark components * have been redefined as two components (quark+diquark, or antiquark+ * antidiquark). * - The (eventual) reshuffling partner, necessary for energy-momentum * conservation when light clusters are decayed into single hadron. Not * all clusters will have a reshuffling partner. * * Notice that in order to determine the cluster position from the positions * of the components, the Cluster class needs some parameters. * Because the Cluster class is neither interfaced nor persistent, * a static pointer to the ClusterHadronizationHandler class instance, * where the parameters are, is used. This static pointer is * set via the method setPointerClusterHadHandler(), during the * run initialization, doinitrun() of ClusterHadronizationHandler. * * @see ClusterHadronizationHandler */ class Cluster : public Particle { protected: /** @name Standard constructors and destructors. */ //@{ /** * Default constructor. Only used in PersistentIStream */ Cluster(); /** * The ClassTraits class must be a friend to be able to * use the private default constructor. */ friend struct ClassTraits; public: /** * Constructor with a particleData pointer */ Cluster(tcEventPDPtr); /** * This creates a cluster from 2 (or 3) partons. */ Cluster(tPPtr part1, tPPtr part2, tPPtr part3 = tPPtr()); /** * Also a constructor where a particle is given not a cluster. */ Cluster(const Particle &); //@} /** * Number of quark (diquark) constituents (normally two). */ unsigned int numComponents() const { return _numComp; } /** * Sum of the constituent masses of the components of the cluster. */ Energy sumConstituentMasses() const; /** * Returns the ith constituent. */ tPPtr particle(int i) const; /** + * Returns the ith original constituent. + */ + tPPtr particleB(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 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/ColourReconnector.cc b/Hadronization/ColourReconnector.cc --- a/Hadronization/ColourReconnector.cc +++ b/Hadronization/ColourReconnector.cc @@ -1,826 +1,1713 @@ // -*- C++ -*- // // ColourReconnector.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the 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","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; + case 0: + _doRecoPlain(clusters); + break; + case 1: + if (_prePlainCR) _doRecoPlain(clusters); + _doRecoStatistical(clusters); + break; + case 2: + if (_prePlainCR) _doRecoPlain(clusters); + _doRecoBaryonic(clusters); + break; + case 3: + if (_prePlainCR) _doRecoPlain(clusters); + _doRecoBaryonicMesonic(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; } +double ColourReconnector::_displacement(tcPPtr p, tcPPtr q) const { + double deltaRap = (p->rapidity() - q->rapidity()); + double deltaPhi = (p->momentum().phi() - q->momentum().phi()); + return sqrt(deltaRap * deltaRap + deltaPhi * deltaPhi); +} + + +double ColourReconnector::_displacementBaryonic(tcPPtr q1, tcPPtr q2, tcPPtr q3) const { + if (_junctionMBCR) { + /** + * Junction-like option i.e. displacement + * from "junction centre" (mean rapidity/phi) + */ + double rap1=q1->rapidity(); + double rap2=q2->rapidity(); + double rap3=q3->rapidity(); + + double phi1=q1->momentum().phi(); + double phi2=q2->momentum().phi(); + double phi3=q3->momentum().phi(); + double meanRap=(rap1 + rap2 + rap3)/3.0; + double meanPhi=(phi1 + phi2 + phi3)/3.0; + double delR; + + delR = sqrt( (rap1-meanRap)*(rap1-meanRap) + (phi1-meanPhi)*(phi1-meanPhi) ); + delR += sqrt( (rap2-meanRap)*(rap2-meanRap) + (phi2-meanPhi)*(phi2-meanPhi) ); + delR += sqrt( (rap3-meanRap)*(rap3-meanRap) + (phi3-meanPhi)*(phi3-meanPhi) ); + return delR; + } else { + /* just summing up all possible 2 quark displacements */ + return _displacement(q1, q2) + _displacement(q1, q3) + _displacement(q2, q3); + } +} 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 ); + 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; +inline bool hasDiquark(CluVecIt cit) { + for (unsigned 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 newcv = cv; - ClusterVector deleted; deleted.reserve(cv.size()); + 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 ); + // 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 + // 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 = _reconnect(*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); + + + +} + + +bool ColourReconnector::_clustersFarApart( const std::vector & clu ) const { + int Ncl=clu.size(); + assert(Ncl<=3); + if (Ncl==1) { + return false; + } else if (Ncl==2) { + // veto if Clusters further apart than _maxDistance + if (_localCR && ((*clu[0])->vertex().vect()-(*clu[1])->vertex().vect()).mag() > _maxDistance) return true; + // veto if Clusters have negative spacetime difference + if (_causalCR && ((*clu[0])->vertex()-(*clu[1])->vertex()).m() < ZERO) return true; + } else if (Ncl==3) { + // veto if Clusters further apart than _maxDistance + if (_localCR && ((*clu[0])->vertex().vect()-(*clu[1])->vertex().vect()).mag() > _maxDistance) return true; + if (_localCR && ((*clu[1])->vertex().vect()-(*clu[2])->vertex().vect()).mag() > _maxDistance) return true; + if (_localCR && ((*clu[0])->vertex().vect()-(*clu[2])->vertex().vect()).mag() > _maxDistance) return true; + // veto if Clusters have negative spacetime difference + if (_causalCR && ((*clu[0])->vertex()-(*clu[1])->vertex()).m() < ZERO) return true; + if (_causalCR && ((*clu[1])->vertex()-(*clu[2])->vertex()).m() < ZERO) return true; + if (_causalCR && ((*clu[0])->vertex()-(*clu[2])->vertex()).m() < ZERO) return true; + } + + return false; +} + + + +void ColourReconnector::_doReco2BeamClusters(ClusterVector & cv) const { + // try other option + PPtr p1Di=(cv[0])->colParticle(); + PPtr p2Di=(cv[1])->colParticle(); + + PPtr p1Q=(cv[0])->antiColParticle(); + PPtr p2Q=(cv[1])->antiColParticle(); + + double min_dist=_displacement(p1Di,p1Q)+_displacement(p2Di,p2Q); + + if ((_displacement(p1Di,p2Q)+_displacement(p1Di,p2Q)) sel; + + unsigned number_of_tries = _stepFactor*cv.size()*cv.size(); + if (number_of_tries<1) number_of_tries=1; + + long (*p_irnd)(long) = UseRandom::irnd; + for (unsigned reconnections_tries = 0; reconnections_tries < number_of_tries; reconnections_tries++) { + num_meson = 0; + num_baryon = 0; + + // flag if we are able to find a suitable combinations of clusters + bool _found = false; + + // Shuffle list of clusters to avoid systematic bias in cluster selection + random_shuffle(newcv.begin(), newcv.end(), p_irnd); + + // loop over clustervector to find CR candidates for (CluVecIt cit = newcv.begin(); cit != newcv.end(); ++cit) { - //avoid clusters already containing diuarks + + // skip the clusters to be deleted later from 3->2 cluster CR + if (find(deleted.begin(), deleted.end(), *cit) != deleted.end()) continue; + + // 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()) + + // add to selection + sel.push_back(cit); + + if (_clustersFarApart(sel)) { + // reject far appart CR + // TODO: after discussion maybe to be omitted + sel.pop_back(); 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); + bool isMeson=((*cit)->numComponents() == 2); - // skip this cluster if no possible reconnection partner can be found - if ( !isBaryonicCandidate && candidate==cit ) + if ( isMeson && (num_meson ==0|| num_meson==1) && num_baryon ==0) { + num_meson++; + /** + * now we habe either 1 or 2 mesonic clusters and have to continue + */ 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; + } else if ( isMeson && (num_baryon == 1 || num_meson ==2)) { + num_meson++; + _found = true; + /** + * we have either 3 mesonic or 1 mesonic and 1 baryonic cluster + * and try to colour reconnect + */ + break; + } else if (num_baryon ==0 && num_meson==0) { + num_baryon++; + /** + * now we have 1 baryonic cluster and have to continue + */ + continue; + } else if (num_meson == 2) { + /** + * we already have 2 mesonic clusters and dont want a baryonic one + * since this is an invalid selection + */ + // remove previously added cluster + sel.pop_back(); + continue; + } else { + num_baryon++; + _found = true; + /** + * now we have either 2 baryonic clusters or 1 mesonic and 1 baryonic cluster + * and try to colour reconnect + */ + break; } } - // 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); + // added for more efficent rejection if some reco probabilities are 0 + if ( _found ) { - swap(cv,clustervector); + // reject MBtoMB candidates if _precoMB_MB=0 + if ( _precoMB_MB == 0 && (num_baryon == 1 && num_meson == 1) ) { + _found=false; + } + // reject BbarBto3M candidates if _precoBbarB_3M=0 + if ( _precoBbarB_3M== 0 && num_baryon == 2 ) { + bool isBbarBto3Mcandidate=( + (*sel[0])->particle(0)->hasColour() && (*sel[1])->particle(0)->hasColour(true) ) + || ( (*sel[0])->particle(0)->hasColour(true) && (*sel[1])->particle(0)->hasColour() ); + if ( isBbarBto3Mcandidate) _found=false; + } + // reject 2Bto2B candidates if _preco2B_2B=0 + if ( _preco2B_2B == 0 && num_baryon == 2 ) { + bool is2Bto2Bcandidate=( + (*sel[0])->particle(0)->hasColour() && (*sel[1])->particle(0)->hasColour() ) + || ( (*sel[0])->particle(0)->hasColour(true) && (*sel[1])->particle(0)->hasColour(true) ); + + if ( is2Bto2Bcandidate ) _found=false; + } + } + // were we able to find a combination? + if (_found==false) { + // clear the selection if we did not find a valid set of clusters + sel.erase(sel.begin(), sel.end()); + continue; + } + assert(sel.size()<4); + assert(sel.size()>1); + + string kind_of_reco = ""; + int reco_info[3]; + + // find best CR option for the selection + _findbestreconnectionoption(sel, num_baryon, kind_of_reco, reco_info); + + if (kind_of_reco == "") { + // no reconnection was found + sel.erase(sel.begin(), sel.end()); + continue; + } else if (kind_of_reco == "3Mto3M" && UseRandom::rnd() < _preco3M_3M) { + // 3Mto3M colour reconnection + auto reconnected = _reconnect3Mto3M(*sel[0], *sel[1], *sel[2], + reco_info); + (*sel[0]) = std::get<0>(reconnected); + (*sel[1]) = std::get<1>(reconnected); + (*sel[2]) = std::get<2>(reconnected); + } else if (kind_of_reco=="2Bto3M" && UseRandom::rnd() < _precoBbarB_3M) { + // antibaryonic and baryonic to 3 mesonic reconnecion + auto reconnected = _reconnectBbarBto3M(*sel[0], *sel[1], + reco_info[0], reco_info[1], reco_info[2]); + (*sel[0]) = std::get<0>(reconnected); + (*sel[1]) = std::get<1>(reconnected); + newcv.push_back(std::get<2>(reconnected)); + } else if (kind_of_reco=="3Mto2B" && UseRandom::rnd() < _preco3M_BBbar) { + // 3 mesonic to antibaryonic and baryonic reconnection + ClusterPtr b1, b2; + _makeBaryonicClusters(*sel[0], *sel[1], *sel[2], b1, b2); + (*sel[0]) = b1; + (*sel[1]) = b2; + deleted.push_back(*sel[2]); + } else if (kind_of_reco=="2Bto2B" && UseRandom::rnd() < _preco2B_2B) { + // 2 (anti)baryonic to 2 (anti)baryonic reconnection + auto reconnected = _reconnect2Bto2B(*sel[0], *sel[1], + reco_info[0], reco_info[1]); + (*sel[0]) = reconnected.first; + (*sel[1]) = reconnected.second; + } else if (kind_of_reco=="MBtoMB" && UseRandom::rnd() < _precoMB_MB) { + // (anti)baryonic and mesonic to (anti)baryonic and mesonic reconnection + auto reconnected = _reconnectMBtoMB(*sel[0], *sel[1], + reco_info[0]); + (*sel[0]) = reconnected.first; + (*sel[1]) = reconnected.second; + } + // erase the sel-vector + sel.erase(sel.begin(), sel.end()); + } + + // write to clustervector new CR'd clusters and deleting + // all deleted clusters + 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, +double calculateRapidityRF(const Lorentz5Momentum & q1, const Lorentz5Momentum & p2) { //calculate rapidity wrt the direction of q1 //angle between the particles in the RF of cluster of q1 // calculate the z component of p2 w.r.t the direction of q1 if(q1.rho2()==ZERO) return 0.; const Energy pz = p2.vect() * q1.vect().unit(); if ( pz == ZERO ) return 0.; - + // Transverse momentum of p2 w.r.t the direction of q1 const Energy pt = sqrt(p2.vect().mag2() - sqr(pz)); - + // Transverse mass pf p2 w.r.t to the direction of q1 const Energy mtrans = sqrt(p2.mass()*p2.mass() + (pt*pt)); // Correct formula const double y2 = log((p2.t() + abs(pz))/mtrans); return ( pz < ZERO ) ? -y2 : y2; } } -CluVecIt ColourReconnector::_findPartnerBaryonic( - CluVecIt cl, ClusterVector & cv, - bool & baryonicCand, - const ClusterVector& deleted, - CluVecIt &baryonic1, - CluVecIt &baryonic2 ) const { +void ColourReconnector::_findbestreconnectionoption(std::vector & cls, const unsigned & baryonic, + string & kind_of_reco, int (&reco_info)[3]) const { + double min_displacement; + if (baryonic==0) { + // case with 3 mesonic clusters + assert(cls.size()==3); - using Constants::pi; - using Constants::twopi; + // calculate the initial displacement sum + min_displacement = _mesonToBaryonFactor * _displacement((*cls[0])->particle(0), (*cls[0])->particle(1)); + min_displacement += _mesonToBaryonFactor * _displacement((*cls[1])->particle(0), (*cls[1])->particle(1)); + min_displacement += _mesonToBaryonFactor * _displacement((*cls[2])->particle(0), (*cls[2])->particle(1)); - // Returns a candidate for possible reconnection - CluVecIt candidate = cl; + // find best CR reco_info and kind_of_reco + _3MtoXreconnectionfinder(cls, + reco_info[0], reco_info[1], reco_info[2], min_displacement, kind_of_reco); + /** + * kind_of_reco either "3Mto3M" or "3Mto2B" (or "" if no better configuration is found) + * case 3Mto3M: the coloured particle of the i-th cluster forms a new cluster with the + * antiparticle of the reco_info[i]-th cluster + * case 3MtoBbarB: all 3 (anti)coloured particle form a new (anti)baryonic cluster + */ + } else if (baryonic == 1) { + // case 1 baryonic and 1 mesonic cluster + assert(cls.size()==2); - 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; - } - } - } - + // make mesonic cluster always the cls[0] + if ((*cls[0])->numComponents() == 3) { + ClusterPtr zw = *cls[0]; + *cls[0] = *cls[1]; + *cls[1] = zw; } - if(bcand == true){ - baryonicCand = true; + // calculate the initial displacement sum + min_displacement = _mesonToBaryonFactor *_displacement((*cls[0])->particle(0), (*cls[0])->particle(1)); + min_displacement += _displacementBaryonic((*cls[1])->particle(0), (*cls[1])->particle(1), (*cls[1])->particle(2)); + + // find best CR reco_info and kind_of_reco + _BMtoBMreconnectionfinder(*cls[0], *cls[1], + reco_info[0], min_displacement, kind_of_reco); + /** + * reco_info[0] is the index of the (anti)quarks of the baryonic cluster cls[1], which should + * be swapped with the (anti)quarks of the mesonic cluster cls[0] + */ + + } else { + assert(baryonic==2); + assert(cls.size()==2); + + // calculate the initial displacement sum + min_displacement = _displacementBaryonic((*cls[0])->particle(0), (*cls[0])->particle(1), (*cls[0])->particle(2)); + min_displacement += _displacementBaryonic((*cls[1])->particle(0), (*cls[1])->particle(1), (*cls[1])->particle(2)); + + // case 2 (anti)baryonic clusters to 2 other (anti)baryonic clusters + if ( ( (*cls[0])->particle(0)->hasColour() && (*cls[1])->particle(0)->hasColour() ) + || ( (*cls[0])->particle(0)->hasColour(true) && (*cls[1])->particle(0)->hasColour(true) ) ) { + // find best CR reco_info and kind_of_reco + _2Bto2BreconnectionFinder(*cls[0], *cls[1], + reco_info[0], reco_info[1], min_displacement, kind_of_reco); + /** + * swap the reco_info[0]-th particle of the first cluster in the vector with the + * reco_info[1]-th particle of the second cluster + */ + } else { + // case 1 baryonic and 1 antibaryonic cluster to 3 mesonic clusters + + // find best CR reco_info and kind_of_reco + _BbarBto3MreconnectionFinder(*cls[0], *cls[1], + reco_info[0], reco_info[1], reco_info[2], min_displacement, kind_of_reco); + /** + * the i-th particle of the first cluster form a new mesonic cluster with the + * reco_info[i]-th particle of the second cluster + */ + } + } + return; +} + + +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; + + // veto if Clusters further apart than _maxDistance + if (_localCR && ((**cl).vertex().vect()-(**cit).vertex().vect()).mag() > _maxDistance) continue; + // veto if Clusters have negative spacetime difference + if (_causalCR && ((**cl).vertex()-(**cit).vertex()).m() < ZERO) 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; + } } - return candidate; + 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 { + 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; + if (cit==cl) continue; // don't allow colour octet clusters if ( _isColour8( (*cl)->colParticle(), - (*cit)->antiColParticle() ) || + (*cit)->antiColParticle() ) || _isColour8( (*cit)->colParticle(), - (*cl)->antiColParticle() ) ) { + (*cl)->antiColParticle() ) ) { continue; } // stop it putting beam remnants together - if((*cl)->isBeamCluster() && (*cit)->isBeamCluster()) continue; + if ((*cl)->isBeamCluster() && (*cit)->isBeamCluster()) continue; - // stop it putting far apart clusters together - if(((**cl).vertex()-(**cit).vertex()).m()>_maxDistance) continue; + // veto if Clusters further apart than _maxDistance + if (_localCR && ((**cl).vertex().vect()-(**cit).vertex().vect()).mag() > _maxDistance) continue; + // veto if Clusters have negative spacetime difference + if (_causalCR && ((**cl).vertex()-(**cit).vertex()).m() < ZERO) continue; // momenta of the old clusters - Lorentz5Momentum p1 = (*cl)->colParticle()->momentum() + + Lorentz5Momentum p1 = (*cl)->colParticle()->momentum() + (*cl)->antiColParticle()->momentum(); - Lorentz5Momentum p2 = (*cit)->colParticle()->momentum() + + Lorentz5Momentum p2 = (*cit)->colParticle()->momentum() + (*cit)->antiColParticle()->momentum(); // momenta of the new clusters - Lorentz5Momentum p3 = (*cl)->colParticle()->momentum() + + Lorentz5Momentum p3 = (*cl)->colParticle()->momentum() + (*cit)->antiColParticle()->momentum(); - Lorentz5Momentum p4 = (*cit)->colParticle()->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{ + 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); + //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()); + 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()); + 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::_reconnect2Bto2B(ClusterPtr &c1, ClusterPtr &c2, const int s1, const int s2) const { + + // form the first new cluster + + // separate the quarks from their original cluster + c1->particleB((s1+1)%3)->abandonChild(c1); + c1->particleB((s1+2)%3)->abandonChild(c1); + c2->particleB(s2)->abandonChild(c2); + + // now the new cluster + ClusterPtr newCluster1 = new_ptr(Cluster(c1->particleB((s1+1)%3), c1->particleB((s1+2)%3), c2->particleB(s2))); + + c1->particleB((s1+1)%3)->addChild(newCluster1); + c1->particleB((s1+2)%3)->addChild(newCluster1); + c2->particleB(s2)->addChild(newCluster1); + + // set new vertex + newCluster1->setVertex(LorentzPoint()); + + // set beam remnants for new cluster + if (c1->isBeamRemnant((s1+1)%3)) newCluster1->setBeamRemnant(0, true); + if (c1->isBeamRemnant((s1+2)%3)) newCluster1->setBeamRemnant(1, true); + if (c2->isBeamRemnant(s2)) newCluster1->setBeamRemnant(2, true); + + // for the second cluster same procedure + c2->particleB((s2+1)%3)->abandonChild(c2); + c2->particleB((s2+2)%3)->abandonChild(c2); + c1->particleB(s1)->abandonChild(c1); + + ClusterPtr newCluster2 = new_ptr(Cluster(c2->particleB((s2+1)%3), c2->particleB((s2+2)%3), c1->particleB(s1))); + + c2->particleB((s2+1)%3)->addChild(newCluster2); + c2->particleB((s2+2)%3)->addChild(newCluster2); + c1->particleB(s1)->addChild(newCluster2); + + newCluster2->setVertex(LorentzPoint()); + + if (c2->isBeamRemnant((s2+1)%3)) newCluster2->setBeamRemnant(0, true); + if (c2->isBeamRemnant((s2+2)%3)) newCluster2->setBeamRemnant(1, true); + if (c1->isBeamRemnant(s1)) newCluster2->setBeamRemnant(2, true); + + return pair (newCluster1, newCluster2); +} + + +std::tuple +ColourReconnector::_reconnectBbarBto3M(ClusterPtr & c1, ClusterPtr & c2, const int s0, const int s1, const int s2) const { + // make sure they all have 3 components + assert(c1->numComponents()==3); + assert(c2->numComponents()==3); + + // first Cluster + c1->particleB(0)->abandonChild(c1); + c2->particleB(s0)->abandonChild(c2); + + ClusterPtr newCluster1 = new_ptr(Cluster(c1->particleB(0), c2->particleB(s0))); + + c1->particleB(0)->addChild(newCluster1); + c2->particleB(s0)->addChild(newCluster1); + + // set new vertex + newCluster1->setVertex(0.5*(c1->particleB(0)->vertex() + c2->particleB(s0)->vertex())); + + // set beam remnants for new cluster + if (c1->isBeamRemnant(0)) newCluster1->setBeamRemnant(0, true); + if (c2->isBeamRemnant(s0)) newCluster1->setBeamRemnant(1, true); + + // same for second cluster + c1->particleB(1)->abandonChild(c1); + c2->particleB(s1)->abandonChild(c2); + + ClusterPtr newCluster2 = new_ptr(Cluster(c1->particleB(1), c2->particleB(s1))); + + c1->particleB(1)->addChild(newCluster2); + c2->particleB(s1)->addChild(newCluster2); + + newCluster2->setVertex(0.5*(c1->particleB(1)->vertex() + c2->particleB(s1)->vertex())); + + if (c1->isBeamRemnant(1)) newCluster2->setBeamRemnant(0, true); + if (c2->isBeamRemnant(s1)) newCluster2->setBeamRemnant(1, true); + + // same for third cluster + c1->particleB(2)->abandonChild(c1); + c2->particleB(s2)->abandonChild(c2); + + ClusterPtr newCluster3 = new_ptr(Cluster(c1->particleB(2), c2->particleB(s2))); + + c1->particleB(2)->addChild(newCluster3); + c2->particleB(s2)->addChild(newCluster3); + + newCluster3->setVertex(0.5*(c1->particleB(2)->vertex() + c2->particleB(s2)->vertex())); + + if (c1->isBeamRemnant(2)) newCluster3->setBeamRemnant(0, true); + if (c2->isBeamRemnant(s2)) newCluster3->setBeamRemnant(1, true); + + return std::tuple (newCluster1, newCluster2, newCluster3); } 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) { + 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); + 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() )); + /* + * TODO: Questionable setting of the vertex + * */ + 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); + c1->antiColParticle()->addChild(newCluster2); + c2->colParticle()->addChild(newCluster2); - newCluster2->setVertex(0.5*( c2->colParticle()->vertex() + - c1->antiColParticle()->vertex() )); + /* + * TODO: Questionable setting of the vertex + * */ + 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); } +std::tuple +ColourReconnector::_reconnect3Mto3M(ClusterPtr & c1, ClusterPtr & c2, ClusterPtr & c3, const int infos [3]) const { + // check if mesonic clusters + assert(c1->numComponents()==2); + assert(c2->numComponents()==2); + assert(c3->numComponents()==2); + + ClusterVector oldclusters = {c1, c2, c3}; + ClusterVector newclusters; + + for (int i=0; i<3; i++) { + + int c1_col=-1; + int c2_anticol=-1; + + // get which index is coloured and which anticolour + for (unsigned int ix=0; ix<2; ++ix) { + if (oldclusters[i]->particle(ix)->hasColour(false)) c1_col = ix; + if (oldclusters[infos[i]]->particle(ix)->hasColour(true)) c2_anticol = ix; + } + + assert(c1_col>=0); + assert(c2_anticol>=0); + + oldclusters[i]->colParticle()->abandonChild(oldclusters[i]); + oldclusters[infos[i]]->antiColParticle()->abandonChild(oldclusters[infos[i]]); + + // form new cluster + ClusterPtr newCluster = new_ptr(Cluster(oldclusters[i]->colParticle(), oldclusters[infos[i]]->antiColParticle())); + + oldclusters[i]->colParticle()->addChild(newCluster); + oldclusters[infos[i]]->antiColParticle()->addChild(newCluster); + + // set new vertex + newCluster->setVertex(0.5*(oldclusters[i]->colParticle()->vertex() + + oldclusters[infos[i]]->antiColParticle()->vertex())); + + // set beam remnants for new cluster + if (oldclusters[i]->isBeamRemnant(c1_col)) newCluster->setBeamRemnant(0, true); + if (oldclusters[infos[i]]->isBeamRemnant(c2_anticol)) newCluster->setBeamRemnant(1, true); + newclusters.push_back(newCluster); + } + return std::tuple (newclusters[0], newclusters[1], newclusters[2]); +} + + +pair +ColourReconnector::_reconnectMBtoMB(ClusterPtr & c1, ClusterPtr & c2, const int s0) const { + // make c1 the mesonic cluster + if (c1->numComponents()==2) { + assert(c2->numComponents()==3); + } else { + return _reconnectMBtoMB(c2,c1,s0); + } + + int c1_col=-1; + int c1_anti=-1; + // get which index is coloured and which anticolour + 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; + + } + assert(c1_col>=0); + assert(c1_anti>=0); + + // pointers for the new clusters + ClusterPtr newCluster1; + ClusterPtr newCluster2; + if (c2->particle(0)->hasColour()==true) { + // first case: we have a baryonic clusters + + // first make the new mesonic cluster + c1->antiColParticle()->abandonChild(c1); + c2->particleB(s0)->abandonChild(c2); + + newCluster1 = new_ptr(Cluster(c1->antiColParticle(), c2->particleB(s0))); + + c1->antiColParticle()->addChild(newCluster1); + c2->particleB(s0)->addChild(newCluster1); + + // set new vertex + newCluster1->setVertex(0.5*(c1->antiColParticle()->vertex() + + c2->particleB(s0)->vertex())); + + // set beam remnants for new cluster + if (c1->isBeamRemnant(c1_anti)) newCluster1->setBeamRemnant(0, true); + if (c2->isBeamRemnant(s0)) newCluster1->setBeamRemnant(1, true); + + // then the baryonic one + c1->colParticle()->abandonChild(c1); + c2->particleB((s0+1)%3)->abandonChild(c2); + c2->particleB((s0+2)%3)->abandonChild(c2); + + newCluster2 = new_ptr(Cluster(c1->colParticle(), c2->particleB((s0+1)%3), c2->particleB((s0+2)%3))); + + c1->colParticle()->addChild(newCluster2); + c2->particleB((s0+1)%3)->addChild(newCluster2); + c2->particleB((s0+2)%3)->addChild(newCluster2); + + // set new vertex + newCluster2->setVertex(LorentzPoint()); + } else { + // second case we have an antibaryonic cluster + + // first make the new mesonic cluster + c1->colParticle()->abandonChild(c1); + c2->particleB(s0)->abandonChild(c2); + + newCluster1 = new_ptr(Cluster(c1->colParticle(), c2->particleB(s0))); + + c1->colParticle()->addChild(newCluster1); + c2->particleB(s0)->addChild(newCluster1); + + // set new vertex + newCluster1->setVertex(0.5*(c1->colParticle()->vertex() + + c2->particleB(s0)->vertex())); + + // set beam remnants for new cluster + if (c1->isBeamRemnant(c1_col)) newCluster1->setBeamRemnant(0, true); + if (c2->isBeamRemnant(s0)) newCluster1->setBeamRemnant(1, true); + + // then the baryonic one + c1->antiColParticle()->abandonChild(c1); + c2->particleB((s0+1)%3)->abandonChild(c2); + c2->particleB((s0+2)%3)->abandonChild(c2); + + newCluster2 = new_ptr(Cluster(c1->antiColParticle(), c2->particleB((s0+1)%3), c2->particleB((s0+2)%3))); + + c1->antiColParticle()->addChild(newCluster2); + c2->particleB((s0+1)%3)->addChild(newCluster2); + c2->particleB((s0+2)%3)->addChild(newCluster2); + + // set new vertex + newCluster2->setVertex(LorentzPoint()); + } + return pair (newCluster1, newCluster2); +} + +void ColourReconnector::_2Bto2BreconnectionFinder(ClusterPtr & c1, ClusterPtr & c2, + int & bswap1, int & bswap2, double min_displ_sum, string & kind_of_reco) const { + double tmp_delta; + for (int i=0; i<3; i++) { + for (int j=0; j<3; j++) { + // try swapping particle i of c1 with particle j of c2 + tmp_delta = _displacementBaryonic(c2->particle(j), c1->particle((i+1)%3), c1->particle((i+2)%3)); + tmp_delta += _displacementBaryonic(c1->particle(i), c2->particle((j+1)%3), c2->particle((j+2)%3)); + + if (tmp_delta < min_displ_sum) { + // if minimal displacement select the 2Bto2B CR option + min_displ_sum = tmp_delta; + bswap1 = i; + bswap2 = j; + kind_of_reco = "2Bto2B"; + } + } + } + +} + +void ColourReconnector::_BbarBto3MreconnectionFinder(ClusterPtr & c1, ClusterPtr & c2, int & mswap0, int & mswap1, int & mswap2, + double min_displ_sum, string & kind_of_reco) const { + double pre_tmp_delta; + double tmp_delta; + for (int p1=0; p1 <3; p1++) { + // make sure not to form a mesonic octet + if (_isColour8(c1->particle(0), c2->particle(p1))) continue; + + pre_tmp_delta = _displacement(c1->particle(0), c2->particle(p1)); + for (int p2=1; p2<3; p2++) { + + // make sure not to form a mesonic octet + if (_isColour8(c1->particle(1), c2->particle((p1+p2)%3))) continue; + if (_isColour8(c1->particle(2), c2->particle(3-p1-((p1+p2)%3)))) continue; + + tmp_delta = pre_tmp_delta + _displacement(c1->particle(1), c2->particle((p1+p2)%3)); + tmp_delta += _displacement(c1->particle(2), c2->particle(3-p1-((p1+p2)%3))); + + // factor _mesonToBaryonFactor to compare Baryonic an mesonic cluster + tmp_delta *=_mesonToBaryonFactor; + + if (tmp_delta < min_displ_sum) { + // if minimal displacement select the 2Bto3M CR option + min_displ_sum = tmp_delta; + mswap0 = p1; + mswap1 = (p1+p2)%3; + mswap2 = 3-p1-((p1+p2)%3); + kind_of_reco = "2Bto3M"; + + } + } + } +} + +void ColourReconnector::_BMtoBMreconnectionfinder(ClusterPtr & c1, ClusterPtr & c2, int & swap, double min_displ_sum, + string & kind_of_reco) const { + assert(c1->numComponents()==2); + assert(c2->numComponents()==3); + double tmp_displ = 0; + for (int i=0; i<3; i++) { + // Differ if the second cluster is baryonic or antibaryonic + if (c2->particle(0)->hasColour()) { + // c2 is baryonic + + // veto mesonic octets + if (_isColour8(c2->particle(i), c1->antiColParticle())) continue; + + // factor _mesonToBaryonFactor to compare Baryonic an mesonic cluster + tmp_displ = _mesonToBaryonFactor * _displacement(c2->particle(i), c1->antiColParticle()); + tmp_displ += _displacementBaryonic(c1->colParticle(), c2->particle((i+1)%3), c2->particle((i+2)%3)); + } else { + // c2 is antibaryonic + + // veto mesonic octets + if (_isColour8(c2->particle(i), c1->colParticle())) continue; + + // factor _mesonToBaryonFactor to compare Baryonic an mesonic cluster + tmp_displ = _mesonToBaryonFactor * _displacement(c2->particle(i), c1->colParticle()); + tmp_displ *= _displacementBaryonic(c1->antiColParticle(), c2->particle((i+1)%3), c2->particle((i+2)%3)); + } + if (tmp_displ < min_displ_sum) { + // if minimal displacement select the MBtoMB CR option + min_displ_sum = tmp_displ; + swap = i; + kind_of_reco = "MBtoMB"; + } + } + return; +} + +void ColourReconnector::_3MtoXreconnectionfinder(std::vector & cv, int & swap0, int & swap1, + int & swap2, double min_displ_sum, string & kind_of_reco) const { + // case of 3M->BbarB CR + double _tmp_displ; + _tmp_displ = _displacementBaryonic((*cv[0])->colParticle(), (*cv[1])->colParticle(), (*cv[2])->colParticle()); + _tmp_displ += _displacementBaryonic((*cv[0])->antiColParticle(), (*cv[1])->antiColParticle(), (*cv[2])->antiColParticle()); + if (_tmp_displ < min_displ_sum) { + // if minimal displacement select the 3Mto2B CR option + kind_of_reco = "3Mto2B"; + min_displ_sum = _tmp_displ; + } + // case for 3M->3M CR + /** + * if 3Mto3M reco probability (_preco3M_3M) is 0 we skip this loop + * since no 3Mto3M CR shall be performed + */ + int i,j; + int i1,i2,i3; + for (i = 0; _preco3M_3M && i<3; i++) { + // veto mesonic octets + if (_isColour8((*cv[0])->colParticle(), (*cv[i])->antiColParticle())) continue; + + // factor _mesonToBaryonFactor to compare baryonic an mesonic cluster + _tmp_displ = _mesonToBaryonFactor * _displacement((*cv[0])->colParticle(), (*cv[i])->antiColParticle()); + for (j=1; j<3; j++) { + // i1, i2, i3 are pairwise distinct + i1=i; + i2=((j+i)%3); + if (i1==0 && i2==1) continue; + i3=(3-i-((j+i)%3)); + + // veto mesonic octets + if (_isColour8((*cv[1])->colParticle(), (*cv[i2])->antiColParticle())) continue; + if (_isColour8((*cv[2])->colParticle(), (*cv[i3])->antiColParticle())) continue; + + _tmp_displ += _mesonToBaryonFactor * _displacement((*cv[1])->colParticle(), (*cv[i2])->antiColParticle()); + _tmp_displ += _mesonToBaryonFactor * _displacement((*cv[2])->colParticle(), (*cv[i3])->antiColParticle()); + + if (_tmp_displ < min_displ_sum) { + // if minimal displacement select the 3Mto3M CR option + kind_of_reco = "3Mto3M"; + min_displ_sum = _tmp_displ; + swap0 = i1; + swap1 = i2; + swap2 = i3; + } + } + } +} + pair ColourReconnector::_shuffle - (const PVector & q, const PVector & aq, unsigned maxtries) const { +(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); + 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 ); + 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 { + } 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 + // 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; + os + << _clreco + << _algorithm + << _annealingFactor + << _annealingSteps + << _triesPerStepFactor + << _initTemp + << _preco + << _precoBaryonic + << _preco3M_3M + << _preco3M_BBbar + << _precoBbarB_3M + << _preco2B_2B + << _precoMB_MB + << _stepFactor + << _mesonToBaryonFactor + << ounit(_maxDistance, femtometer) + << _octetOption + << _cr2BeamClusters + << _prePlainCR + << _localCR + << _causalCR + << _debug + << _junctionMBCR + ; } void ColourReconnector::persistentInput(PersistentIStream & is, int) { - is >> _clreco >> _preco >> _precoBaryonic >> _algorithm >> _initTemp >> _annealingFactor - >> _annealingSteps >> _triesPerStepFactor >> iunit(_maxDistance,femtometer) - >> _octetOption; + is + >> _clreco + >> _algorithm + >> _annealingFactor + >> _annealingSteps + >> _triesPerStepFactor + >> _initTemp + >> _preco + >> _precoBaryonic + >> _preco3M_3M + >> _preco3M_BBbar + >> _precoBbarB_3M + >> _preco2B_2B + >> _precoMB_MB + >> _stepFactor + >> _mesonToBaryonFactor + >> iunit(_maxDistance, femtometer) + >> _octetOption + >> _cr2BeamClusters + >> _prePlainCR + >> _localCR + >> _causalCR + >> _debug + >> _junctionMBCR + ; } void ColourReconnector::Init() { - static ClassDocumentation documentation - ("This class is responsible of the colour reconnection."); + ("This class is responsible of the colour reconnection."); static Switch interfaceColourReconnection - ("ColourReconnection", - "Colour reconnections", - &ColourReconnector::_clreco, 0, true, false); + ("ColourReconnection", + "Colour reconnections", + &ColourReconnector::_clreco, 0, true, false); static SwitchOption interfaceColourReconnectionNo - (interfaceColourReconnection, - "No", - "Colour reconnections off", - 0); + (interfaceColourReconnection, + "No", + "Colour reconnections off", + 0); static SwitchOption interfaceColourReconnectionYes - (interfaceColourReconnection, - "Yes", - "Colour reconnections on", - 1); - + (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); + // Algorithm interface + 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 SwitchOption interfaceAlgorithmBaryonicMesonic + (interfaceAlgorithm, + "BaryonicMesonic", + "Baryonic cluster reconnection with reconnections to and from Mesonic Clusters", + 3); + + + + // Statistical CR Parameters: + 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); + ("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); + ("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); + ("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); + + + // Plain and Baryonic CR Paramters + static Parameter interfaceRecoProb + ("ReconnectionProbability", + "Probability that a found two meson to two meson reconnection possibility is actually accepted (used in Plain & Baryonic)", + &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); + ("ReconnectionProbabilityBaryonic", + "Probability that a found reconnection possibility is actually accepted (used in Baryonic)", + &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); + // BaryonicMesonic CR Paramters + static Parameter interfaceReconnectionProbability3Mto3M + ("ReconnectionProbability3Mto3M", + "Probability that a reconnection candidate is accepted for reconnecting 3M -> 3M\'", + &ColourReconnector::_preco3M_3M, 0.5, 0.0, 1.0, + false, false, Interface::limited); + static Parameter interfaceReconnectionProbability3MtoBBbar + ("ReconnectionProbability3MtoBBbar", + "Probability that a reconnection candidate is accepted for reconnecting 3M -> B,Bbar", + &ColourReconnector::_preco3M_BBbar, 0.5, 0.0, 1.0, + false, false, Interface::limited); + static Parameter interfaceReconnectionProbabilityBbarBto3M + ("ReconnectionProbabilityBbarBto3M", + "Probability that a reconnection candidate is accepted for reconnecting B,Bbar -> 3M", + &ColourReconnector::_precoBbarB_3M, 0.5, 0.0, 1.0, + false, false, Interface::limited); + static Parameter interfaceReconnectionProbability2Bto2B + ("ReconnectionProbability2Bto2B", + "Probability that a reconnection candidate is accepted for reconnecting 2B -> 2B\' or 2Bbar -> 2Bbar\'", + &ColourReconnector::_preco2B_2B, 0.5, 0.0, 1.0, + false, false, Interface::limited); + static Parameter interfaceReconnectionProbabilityMBtoMB + ("ReconnectionProbabilityMBtoMB", + "Probability that a reconnection candidate is accepted for reconnecting M,B -> M\',B\' or M,Bbar -> M\',Bbar\'", + &ColourReconnector::_precoMB_MB, 0.5, 0.0, 1.0, + false, false, Interface::limited); - 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 Parameter interfaceFactorforStep + ("StepFactor", + "Factor for how many reconnection-tries are made in the BaryonicMesonic algorithm", + &ColourReconnector::_stepFactor, 1.0, 0.11111, 10., + false, false, Interface::limited);// at least 3 Clusters -> _stepFactorMin=1/9 + static Parameter interfaceMesonToBaryonFactor + ("MesonToBaryonFactor", + "Factor for comparing mesonic clusters to baryonic clusters in the displacement if BaryonicMesonic CR model is chosen", + &ColourReconnector::_mesonToBaryonFactor, 2.0, 0.5, 3.0, + false, false, Interface::limited); - static Switch interfaceOctetTreatment - ("OctetTreatment", - "Which octets are not allowed to be reconnected", - &ColourReconnector::_octetOption, 0, false, false); + + + // General Parameters and switches + static Parameter interfaceMaxDistance + ("MaxDistance", + "Maximum distance between the clusters at which to consider rearrangement" + " to avoid colour reconneections of displaced vertices (used in all Algorithms). No unit means femtometer", + &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 (used in all Algorithms)", + &ColourReconnector::_octetOption, 0, false, false); static SwitchOption interfaceOctetTreatmentFinal - (interfaceOctetTreatment, - "Final", - "Only prevent for the final (usuaslly non-perturbative) g -> q qbar splitting", - 0); + (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); + (interfaceOctetTreatment, + "All", + "Prevent for all octets", + 1); + static Switch interfaceCR2BeamClusters + ("CR2BeamClusters", + "Option for colour reconnecting 2 beam remnant clusters if the number of clusters is 2.", + &ColourReconnector::_cr2BeamClusters, 0, true, false); + static SwitchOption interfaceCR2BeamClustersYes + (interfaceCR2BeamClusters, + "Yes", + "If possible CR 2 beam clusters", + 1); + static SwitchOption interfaceCR2BeamClustersNo + (interfaceCR2BeamClusters, + "No", + "If possible do not CR 2 beam clusters", + 0); + static Switch interfacePrePlainCR + ("PrePlainCR", + "Perform a Plain CR before another algorithm is performed", + &ColourReconnector::_prePlainCR, 0, true, false); + static SwitchOption interfacePrePlainCRYes + (interfacePrePlainCR, + "Yes", + "enable pre-Plain CR", + 1); + static SwitchOption interfacePrePlainCRNo + (interfacePrePlainCR, + "No", + "disable pre-Plain CR", + 0); + static Switch interfaceLocalCR + ("LocalCR", + "Option for colour reconnecting only if clusters are less distant than MaxDistance", + &ColourReconnector::_localCR, 0, true, false); + static SwitchOption interfaceLocalCRYes + (interfaceLocalCR, + "Yes", + "activate spatial veto", + 1); + static SwitchOption interfaceLocalCRNo + (interfaceLocalCR, + "No", + "deactivate spatial veto", + 0); + static Switch interfaceCausalCR + ("CausalCR", + "Option for colour reconnecting only if clusters their vertices " + "have a positive spacetime difference", + &ColourReconnector::_causalCR, 0, true, false); + static SwitchOption interfaceCausalCRYes + (interfaceCausalCR, + "Yes", + "enable causal veto", + 1); + static SwitchOption interfaceCausalCRNo + (interfaceCausalCR, + "No", + "disable causal veto", + 0); + static Switch interfaceJunction + ("Junction", + "Option for using Junction-like displacement in rapidity-phi plane to compare baryonic cluster " + "instead of pairwise distance (for BaryonicMesonic model)", + &ColourReconnector::_junctionMBCR, 1, true, false); + static SwitchOption interfaceJunctionYes + (interfaceJunction, + "Yes", + "Using junction-like model instead of pairwise distance model", + 1); + static SwitchOption interfaceJunctionNo + (interfaceJunction, + "No", + "Using pairwise distance model instead of junction-like model", + 0); + + // Debug + static Switch interfaceDebug + ("Debug", + "Make a file with some Information of the BaryonicMesonic Algorithm", + &ColourReconnector::_debug, 0, true, false); + static SwitchOption interfaceDebugNo + (interfaceDebug, + "No", + "Debug Information for ColourReconnector Off", + 0); + static SwitchOption interfaceDebugYes + (interfaceDebug, + "Yes", + "Debug Information for ColourReconnector On", + 1); } diff --git a/Hadronization/ColourReconnector.h b/Hadronization/ColourReconnector.h --- a/Hadronization/ColourReconnector.h +++ b/Hadronization/ColourReconnector.h @@ -1,271 +1,537 @@ // -*- C++ -*- // // ColourReconnector.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_ColourReconnector_H #define HERWIG_ColourReconnector_H #include #include "CluHadConfig.h" #include "ColourReconnector.fh" namespace Herwig { using namespace ThePEG; /** \ingroup Hadronization * \class ColourReconnector * \brief Class for changing colour reconnections of partons. * \author Alberto Ribon, Christian Roehr * * This class does the nonperturbative colour rearrangement, after the * nonperturbative gluon splitting and the "normal" cluster formation. * It uses the list of particles in the event record, and the collections of * "usual" clusters which is passed to the main method. If the colour * reconnection is actually accepted, then the previous collections of "usual" * clusters is first deleted and then the new one is created. * * * @see \ref ColourReconnectorInterfaces "The interfaces" * defined for ColourReconnector. */ class ColourReconnector: public Interfaced { public: /** * Does the colour rearrangement, starting out from the list of particles in * the event record and the collection of "usual" clusters passed as * arguments. If the actual rearrangement is accepted, the initial collection of * clusters is overridden by the old ones. */ void rearrange(ClusterVector & clusters); using CluVecIt = ClusterVector::iterator; private: /** PRIVATE MEMBER FUNCTIONS */ /** * @brief Calculates the sum of the squared cluster masses. * @arguments q, aq vectors containing the quarks and antiquarks respectively * @return Sum of cluster squared masses M^2_{q[i],aq[i]}. */ Energy2 _clusterMassSum(const PVector & q, const PVector & aq) const; + + + + /** + * @brief calculates the "euclidean" distance of two quarks in the + * rapdidity-phi plane + * @arguments p, q the two quarks + * @return the dimensionless distance: + * \deltaR_12=sqrt(\deltaY_12^2 + \delta\phi_12^2) + */ + double _displacement(tcPPtr p, tcPPtr q) const; + + + /** + * @brief calculates the "euclidean" distance of a the 3 (anti)quarks + * (anti)baryonic cluster in the rapdidity-phi plane + * @arguments p, q the two quarks + * @return the dimensionless distance: + * if Junction is enabled the difference of all 3 quarks + * with respect to their mean point is calculated: + * = sum_i Y_i/3 + * <\phi> = sum_i \phi_i/3 + * \deltaR = sum_i sqrt( (Y_i - )^2 + (\phi_i - )^2) + * + * if Junction is disabled the difference of all distinct + * pairing of the 3 quarks is computed: + * \deltaR_ij = sqrt(\deltaY_ij^2 + \delta\phi_ij^2) + * \deltaR_tot = sum_{i,j & P) const; /** * @brief A Metropolis-type algorithm which finds a local minimum in the * total sum of cluster masses * @arguments cv cluster vector */ void _doRecoStatistical(ClusterVector & cv) const; /** * @brief Plain colour reconnection as used in Herwig 2.5.0 * @arguments cv cluster vector */ void _doRecoPlain(ClusterVector & cv) const; /** * Baryonic Colour Reconnection model */ void _doRecoBaryonic(ClusterVector & cv) const; + + + /** + * @brief BaryonicMesonic colour reconnection model + * @arguments cv cluster vector + * BaryonicMesonic Colour Reconnection model with reconnections from mesonic clusters + * to baryonic cluster and the contrary. Allows also reconnection between mesonic + * and baryonic Clusters + */ + void _doRecoBaryonicMesonic(ClusterVector & cv) const; + + void _makeBaryonicClusters(ClusterPtr &c1, ClusterPtr &c2, ClusterPtr &c3, ClusterPtr &newcluster1, ClusterPtr &newcluster2) const; /** * @brief Finds the cluster in cv which, if reconnected with the given * cluster cl, would result in the smallest sum of cluster masses. * If no reconnection partner can be found, a pointer to the * original Cluster cl is returned. * @arguments cv cluster vector * cl cluster iterator (must be from cv) which wants to have a reconnection partner * @return iterator to the found cluster, or the original cluster pointer if * no mass-reducing combination can be found */ CluVecIt _findRecoPartner(CluVecIt cl, ClusterVector & cv) const; CluVecIt _findPartnerRapidity(CluVecIt cl, ClusterVector & cv) const; CluVecIt _findPartnerBaryonic(CluVecIt cl, ClusterVector & cv, bool & tetraCand, const ClusterVector& a, CluVecIt &baryonic1, CluVecIt &baryonic2 ) const; - - - + + /** + * @brief Finds best CR option for the BaryonicMesonic CR model + * @arguments cls vector of selected clusters, baryonic is the number of baryonic + * clusters in selection, kind_of_reco is output string denoting best + * CR option and reco_info is output storage for information on which + * (anti-)quarks to exchange and in which way. + * BaryonicMesonic Colour Reconnection model with reconnections from mesonic clusters + * to baryonic cluster and the contrary. Allows also reconnection between mesonic + * and baryonic Clusters + */ + void _findbestreconnectionoption(std::vector &cls, + const unsigned &baryonic, + string &kind_of_reco, + int (&reco_info)[3]) const; /** * @brief Reconnects the constituents of the given clusters to the (only) * other possible cluster combination. * @return pair of pointers to the two new clusters + * Used for Plain and Baryonic Colour Reconnection models */ pair _reconnect(ClusterPtr &c1, ClusterPtr &c2) const; /** - * Reconnection method for baryonic reconenction model + * @brief Reconnects (2B->2B) the constituents of the given clusters to + another possible cluster combination whose information is given + in s1 and s2. + * @arguments c1 and c2 are baryonic clusters and s1 and s2 are the respective + indices of the constituents of c1 and c2 respectively + * @return pair of pointers to the two new clusters + * Used only in BaryonicMesonic algorithm and will exchange constituent s1 of + * c1 with constituent s2 of c2 + */ + pair _reconnect2Bto2B(ClusterPtr &c1, ClusterPtr &c2, const int s1, const int s2) const; + + /** + * @brief Reconnects (B,Bbar->3M) the constituents of the given clusters to + another possible cluster combination whose information is given in + s0, s1 and s2. + * @arguments c1 and c2 are baryonic clusters. s0, s1 and s2 are the respective + indices which determine the CR + * @return tuple of pointers to the 3 new mesonic clusters + * Used only in BaryonicMesonic algorithm and will form 3 mesonic clusters according + * to the indices s0, s1 and s2. The i-th constituent of c1 is connected to the si-th + * constituent of c2 + */ + std::tuple _reconnectBbarBto3M(ClusterPtr &c1, ClusterPtr &c2, const int s0, const int s1, const int s2 ) const; + + /** + * @brief Reconnects (3M->3M) the constituents of the given clusters to + another possible cluster combination whose information is given in + sinfos + * @arguments c1, c2 and c3 are mesonic clusters. infos[3] are the respective + indices which determine the CR + * @return tuple of pointers to the 3 CR'd mesonic clusters + * Used only in BaryonicMesonic algorithm and will reconnect 3 mesonic clusters according + * to the infos, which determine the CR. The coloured quark of the i-th cluster forms + * a new cluster with the anticoloured quark of the info[i]-th cluster + */ + std::tuple _reconnect3Mto3M(ClusterPtr &c1, ClusterPtr &c2, ClusterPtr &c3, const int infos[3] ) const; + + + /** + * Reconnection method for a Baryonic and a Mesonic Cluster to a Baryonic and a Mesonic Cluster + * s0 is the Number of the (Anti)Patrticle of the Baryonic Cluster , which should be swapped with the Anti(Particle) of the Mesonic Cluster */ - pair _reconnectBaryonic(ClusterPtr &c1, ClusterPtr &c2) const; + /** + * @brief Reconnects the constituents of the given clusters to + another possible cluster combination whose information + is given in s0. + * @arguments c1 and c2 are one baryonic and one mesonic cluster respectively + and s0 is the respective index of the constituents of the baryonic + cluster which is to be exchangeed with the mesonic cluster. + * @return pair of pointers to the two new clusters + * Used only in BaryonicMesonic algorithm and will exchange constituent s0 of + * the baryonic cluster with the (anti)-quark of the mesonic cluster + */ + pair _reconnectMBtoMB(ClusterPtr &c1, ClusterPtr &c2, const int s0) const; + + + + + /** + * Methods for the BaryonicMesonic CR algorithm + * Find the best reconnection option for the respective cluster-combination + * + */ + /** + * @brief veto for far apart clusters + * @arguments expects at most 3 CluVecIt in clu vector + * @return returns true if clusters are more distant than _maxDistance + * in space + * TODO: problematic maybe add option to turn off + */ + bool _clustersFarApart( const std::vector & clu ) const; + + /** + * @brief Does reconnect 2 beam clusters for BaryonicMesonic CR model + * if option CR2BeamClusters is enabled + * @arguments cv cluster vector + */ + void _doReco2BeamClusters(ClusterVector & cv) const; + + + /** + * @brief finds the best reconnection option and stores it in bswap1 + * and bswap2 (2B->2B colour reconnection) + * @arguments c1 and c2 cluster pointers and kind_of_reco will output + * the best found reconnection option for c1 and c2 + */ + void _2Bto2BreconnectionFinder(ClusterPtr &c1, ClusterPtr &c2, int &bswap1, int &bswap2, double mindisplsum, string &kind_of_reco ) const; + + /** + * @brief finds the best reconnection option and stores it in mswap0 + * mswap1 and mswap2 (BbarB->3M colour reconnection) + * @arguments c1 and c2 cluster pointers and kind_of_reco will output + * the best found reconnection option for c1 and c2 + */ + void _BbarBto3MreconnectionFinder(ClusterPtr &c1, ClusterPtr &c2, int &mswap0, int &mswap1, int &mswap2, + double min_displ_sum, string & kind_of_reco) const; + + /** + * @brief finds the best reconnection option and stores it in swap + * (BM->BM colour reconnection) + * @arguments c1 and c2 cluster pointers and kind_of_reco will output + * the best found reconnection option for c1 and c2 + */ + void _BMtoBMreconnectionfinder(ClusterPtr &c1, ClusterPtr &c2, int &swap, + double min_displ_sum, string &kind_of_reco) const; + + /** + * @brief finds the best reconnection option and stores it in swap0, + * swap1 and swap2 (3M->{3M,BbarB} colour reconnection) + * @arguments c1 and c2 cluster pointers and kind_of_reco will output + * the best found reconnection option for c1 and c2 + */ + void _3MtoXreconnectionfinder(std::vector &cv, int &swap0, int &swap1, + int &swap2, double min_displ_sum, string &kind_of_reco) const; /** * @brief At random, swap two antiquarks, if not excluded by the * constraint that there must not be any colour-octet clusters. * @arguments q, aq vectors containing the quarks and antiquarks respectively * maxtries maximal number of tries to find a non-colour-octet * reconfiguration * @return Pair of ints indicating the indices of the antiquarks to be * swapped. Returns (-1,-1) if no valid reconfiguration could be * found after maxtries trials */ pair _shuffle(const PVector & q, const PVector & aq, unsigned maxtries = 10) const; /** DATA MEMBERS */ /** * Specifies the colour reconnection algorithm to be used. */ int _algorithm = 0; /** * Do we do colour reconnections? */ int _clreco = 0; + /** + * Do we want debug informations? + */ + int _debug = 0; + + /** + * @brief Junction-like model for BaryonicMesonic model + * Do we want to use the junction-like model for + * computing the displacements of BaryonicMesonic model? + * otherwise pairwise distances are used. + * If _junctionMBCR is activated the displacements are computed in the + * rapidity-Phi plane by difference to the average rapidity and phi: + * DeltaR_i^2 = (rapidity_i - meanRap)^2 + (phi_i - meanPhi)^2 + * DeltaR = sum_i DeltaR_i + * if _junctionMBCR=0 the displacements are computed: + * DeltaR_ij^2 = (rapidity_i - rapidity_j)^2 + (phi_i - phi_j)^2 + * DeltaR = sum_i,j 3M' + * used in BaryonicMesonic + * NOTE: if 0 this type of reconnection is not even tried */ + double _preco3M_3M = 0.5; + + /** + * Probability that a found reconnection possibility is actually accepted. + * For reconnecting 3M -> B,Bbar + * used in BaryonicMesonic + */ + double _preco3M_BBbar = 0.5; + + /** + * Probability that a found reconnection possibility is actually accepted. + * For reconnecting Bbar,B -> 3M + * used in BaryonicMesonic + */ + double _precoBbarB_3M = 0.5; + + /** + * Probability that a found reconnection possibility is actually accepted. + * For reconnecting 2B -> 2B' or 2Bbar -> 2Bbar' + * used in BaryonicMesonic + * NOTE: if 0 this type of reconnection is not even tried + */ + double _preco2B_2B = 0.5; + + /** + * Probability that a found reconnection possibility is actually accepted. + * For reconnecting M,B -> M',B' or M,Bbar -> M',Bbar' + * used in BaryonicMesonic + * NOTE: if 0 this type of reconnection is not even tried + */ + double _precoMB_MB = 0.5; + + /** + * For the BaryonicMesonic algorithm + * How many times do suggest cluster for reconnection? + * n(reconnectionstries) = _stepFactor * n(clusters)*n(clusters); + */ + double _stepFactor = 5.0; + + /** + * Factor for comparing mesonic clusters to baryonic clusters + */ + double _mesonToBaryonFactor = 2.0; /** * Maximium distance for reconnections + * TODO: remove if issues with anticausality are discussed and resolved */ - Length _maxDistance = picometer; + Length _maxDistance = femtometer; /** * @return true, if the two partons are splitting products of the same * gluon */ bool _isColour8(tcPPtr p, tcPPtr q) const; /** * Option for handling octets */ unsigned int _octetOption = 0; + /** + * Option for colour reconnecting 2 Beam Clusters if no others are present + */ + int _cr2BeamClusters = 0; + + /** + * Option for performing Plain colour reconnection before the Statistical, + * Baryonic or BaryonicMesonic algorithm is performed + */ + int _prePlainCR = 0; + + /** + * Option for colour reconnecting Clusters only if their vertex 3-distance + * is less than _maxDistance + */ + int _localCR = 0; + + /** + * Option for colour reconnecting Clusters only if their spacetime difference + * is bigger than 0 + */ + int _causalCR = 0; + public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * Standard Init function used to initialize the interfaces. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} private: /** * Private and non-existent assignment operator. */ ColourReconnector & operator=(const ColourReconnector &) = delete; }; } #endif /* HERWIG_ColourReconnector_H */ diff --git a/src/defaults/Hadronization.in b/src/defaults/Hadronization.in --- a/src/defaults/Hadronization.in +++ b/src/defaults/Hadronization.in @@ -1,98 +1,127 @@ # -*- ThePEG-repository -*- ############################################################ # Setup of default hadronization # # There are no user servicable parts inside. # # Anything that follows below should only be touched if you # know what you're doing. ############################################################# cd /Herwig/Particles create ThePEG::ParticleData Cluster setup Cluster 81 Cluster 0.00990 0.0 0.0 0.0 0 0 0 1 create ThePEG::ParticleData Remnant setup Remnant 82 Remnant 0.00990 0.0 0.0 0.0 0 0 0 1 mkdir /Herwig/Hadronization cd /Herwig/Hadronization create Herwig::ClusterHadronizationHandler ClusterHadHandler create Herwig::PartonSplitter PartonSplitter create Herwig::ClusterFinder ClusterFinder create Herwig::ColourReconnector ColourReconnector create Herwig::ClusterFissioner ClusterFissioner create Herwig::LightClusterDecayer LightClusterDecayer create Herwig::ClusterDecayer ClusterDecayer create Herwig::HwppSelector HadronSelector newdef ClusterHadHandler:PartonSplitter PartonSplitter newdef ClusterHadHandler:ClusterFinder ClusterFinder newdef ClusterHadHandler:ColourReconnector ColourReconnector newdef ClusterHadHandler:ClusterFissioner ClusterFissioner newdef ClusterHadHandler:LightClusterDecayer LightClusterDecayer newdef ClusterHadHandler:ClusterDecayer ClusterDecayer newdef ClusterHadHandler:MinVirtuality2 0.1*GeV2 newdef ClusterHadHandler:MaxDisplacement 1.0e-10*millimeter newdef ClusterHadHandler:UnderlyingEventHandler NULL newdef ClusterFissioner:HadronSelector HadronSelector newdef LightClusterDecayer:HadronSelector HadronSelector newdef ClusterDecayer:HadronSelector HadronSelector +# ColourReconnector Default Parameters newdef ColourReconnector:ColourReconnection Yes +newdef ColourReconnector:Algorithm Baryonic + +# Statistical CR Parameters: +newdef ColourReconnector:AnnealingFactor 0.9 +newdef ColourReconnector:AnnealingSteps 50 +newdef ColourReconnector:TriesPerStepFactor 5.0 +newdef ColourReconnector:InitialTemperature 0.1 + +# Plain and Baryonic CR Paramters +newdef ColourReconnector:ReconnectionProbability 0.95 newdef ColourReconnector:ReconnectionProbabilityBaryonic 0.7 -newdef ColourReconnector:ReconnectionProbability 0.95 -newdef ColourReconnector:Algorithm Baryonic + +# BaryonicMesonic and BaryonicMesonic CR Paramters +newdef ColourReconnector:ReconnectionProbability3Mto3M 0.5 +newdef ColourReconnector:ReconnectionProbability3MtoBBbar 0.5 +newdef ColourReconnector:ReconnectionProbabilityBbarBto3M 0.5 +newdef ColourReconnector:ReconnectionProbability2Bto2B 0.05 +newdef ColourReconnector:ReconnectionProbabilityMBtoMB 0.5 +newdef ColourReconnector:StepFactor 1.0 +newdef ColourReconnector:MesonToBaryonFactor 1.333 + +# General Parameters and switches +newdef ColourReconnector:MaxDistance 1.0e50 newdef ColourReconnector:OctetTreatment All +newdef ColourReconnector:CR2BeamClusters No +newdef ColourReconnector:Junction Yes +newdef ColourReconnector:PrePlainCR No +newdef ColourReconnector:LocalCR No +newdef ColourReconnector:CausalCR No +# Debugging +newdef ColourReconnector:Debug No + # Clustering parameters for light quarks newdef ClusterFissioner:ClMaxLight 3.649 newdef ClusterFissioner:ClPowLight 2.780 newdef ClusterFissioner:PSplitLight 0.899 newdef ClusterDecayer:ClDirLight 1 newdef ClusterDecayer:ClSmrLight 0.78 # Clustering parameters for b-quarks newdef ClusterFissioner:ClMaxBottom 3.757 newdef ClusterFissioner:ClPowBottom 0.547 newdef ClusterFissioner:PSplitBottom 0.625 newdef ClusterDecayer:ClDirBottom 1 newdef ClusterDecayer:ClSmrBottom 0.078 newdef HadronSelector:SingleHadronLimitBottom 0.000 # Clustering parameters for c-quarks newdef ClusterFissioner:ClMaxCharm 3.950 newdef ClusterFissioner:ClPowCharm 2.559 newdef ClusterFissioner:PSplitCharm 0.994 newdef ClusterDecayer:ClDirCharm 1 newdef ClusterDecayer:ClSmrCharm 0.163 newdef HadronSelector:SingleHadronLimitCharm 0.000 # Clustering parameters for exotic quarks # (e.g. hadronizing Susy particles) newdef ClusterFissioner:ClMaxExotic 2.7*GeV newdef ClusterFissioner:ClPowExotic 1.46 newdef ClusterFissioner:PSplitExotic 1.00 newdef ClusterDecayer:ClDirExotic 1 newdef ClusterDecayer:ClSmrExotic 0. newdef HadronSelector:SingleHadronLimitExotic 0. # newdef PartonSplitter:SplitPwtSquark 0.824135 newdef PartonSplitter:Split uds # newdef HadronSelector:PwtDquark 1.0 newdef HadronSelector:PwtUquark 1.0 newdef HadronSelector:PwtSquark 0.291717 newdef HadronSelector:PwtCquark 0.0 newdef HadronSelector:PwtBquark 0.0 newdef HadronSelector:PwtDIquark 0.298 newdef HadronSelector:SngWt 0.74 newdef HadronSelector:DecWt 0.62 newdef HadronSelector:Mode 1 newdef HadronSelector:BelowThreshold All create Herwig::SpinHadronizer SpinHadronizer