diff --git a/Hadronization/ClusterDecayer.cc b/Hadronization/ClusterDecayer.cc --- a/Hadronization/ClusterDecayer.cc +++ b/Hadronization/ClusterDecayer.cc @@ -1,470 +1,480 @@ // -*- C++ -*- // // ClusterDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the ClusterDecayer class. // #include "ClusterDecayer.h" #include #include #include #include #include #include #include #include #include "Herwig/Utilities/Kinematics.h" #include "CheckId.h" #include "Cluster.h" #include #include +#include using namespace Herwig; DescribeClass describeClusterDecayer("Herwig::ClusterDecayer",""); ClusterDecayer::ClusterDecayer() : - _clDirLight(1), + _clDirLight(1), _clDirBottom(1), _clDirCharm(1), - _clDirExotic(1), - _clSmrLight(0.0), + _clDirExotic(1), + _clSmrLight(0.0), _clSmrBottom(0.0), _clSmrCharm(0.0), _clSmrExotic(0.0), _onshell(false), _masstry(20) {} IBPtr ClusterDecayer::clone() const { return new_ptr(*this); } IBPtr ClusterDecayer::fullclone() const { return new_ptr(*this); } -void ClusterDecayer::persistentOutput(PersistentOStream & os) const +void ClusterDecayer::persistentOutput(PersistentOStream & os) const { os << _hadronsSelector << _clDirLight << _clDirBottom - << _clDirCharm << _clDirExotic << _clSmrLight << _clSmrBottom + << _clDirCharm << _clDirExotic << _clSmrLight << _clSmrBottom << _clSmrCharm << _clSmrExotic << _onshell << _masstry; } void ClusterDecayer::persistentInput(PersistentIStream & is, int) { is >> _hadronsSelector >> _clDirLight >> _clDirBottom - >> _clDirCharm >> _clDirExotic >> _clSmrLight >> _clSmrBottom + >> _clDirCharm >> _clDirExotic >> _clSmrLight >> _clSmrBottom >> _clSmrCharm >> _clSmrExotic >> _onshell >> _masstry; } void ClusterDecayer::Init() { static ClassDocumentation documentation ("This class is responsible for the two-body decays of normal clusters"); - static Reference - interfaceHadronSelector("HadronSelector", - "A reference to the HadronSelector object", + static Reference + interfaceHadronSelector("HadronSelector", + "A reference to the HadronSelector object", &Herwig::ClusterDecayer::_hadronsSelector, false, false, true, false); //ClDir for light, Bottom, Charm and exotic (e.g Susy) quarks static Switch interfaceClDirLight ("ClDirLight", "Cluster direction for light quarks", &ClusterDecayer::_clDirLight, true, false, false); static SwitchOption interfaceClDirLightPreserve (interfaceClDirLight, "Preserve", "Preserve the direction of the quark from the perturbative stage" " as the direction of the meson produced from it", true); static SwitchOption interfaceClDirLightIsotropic (interfaceClDirLight, "Isotropic", "Assign the direction of the meson randomly", false); static Switch interfaceClDirBottom ("ClDirBottom", "Cluster direction for bottom quarks", &ClusterDecayer::_clDirBottom, true, false, false); static SwitchOption interfaceClDirBottomPreserve (interfaceClDirBottom, "Preserve", "Preserve the direction of the quark from the perturbative stage" " as the direction of the meson produced from it", true); static SwitchOption interfaceClDirBottomIsotropic (interfaceClDirBottom, "Isotropic", "Assign the direction of the meson randomly", false); static Switch interfaceClDirCharm ("ClDirCharm", "Cluster direction for charm quarks", &ClusterDecayer::_clDirCharm, true, false, false); static SwitchOption interfaceClDirCharmPreserve (interfaceClDirCharm, "Preserve", "Preserve the direction of the quark from the perturbative stage" " as the direction of the meson produced from it", true); static SwitchOption interfaceClDirCharmIsotropic (interfaceClDirCharm, "Isotropic", "Assign the direction of the meson randomly", false); static Switch interfaceClDirExotic ("ClDirExotic", "Cluster direction for exotic quarks", &ClusterDecayer::_clDirExotic, true, false, false); static SwitchOption interfaceClDirExoticPreserve (interfaceClDirExotic, "Preserve", "Preserve the direction of the quark from the perturbative stage" " as the direction of the meson produced from it", true); static SwitchOption interfaceClDirExoticIsotropic (interfaceClDirExotic, "Isotropic", "Assign the direction of the meson randomly", false); // ClSmr for ligth, Bottom, Charm and exotic (e.g. Susy) quarks - static Parameter + static Parameter interfaceClSmrLight ("ClSmrLight", "cluster direction Gaussian smearing for light quark", &ClusterDecayer::_clSmrLight, 0, 0.0 , 0.0 , 2.0,false,false,false); - static Parameter + static Parameter interfaceClSmrBottom ("ClSmrBottom", "cluster direction Gaussian smearing for b quark", - &ClusterDecayer::_clSmrBottom, 0, 0.0 , 0.0 , 2.0,false,false,false); -static Parameter + &ClusterDecayer::_clSmrBottom, 0, 0.0 , 0.0 , 2.0,false,false,false); +static Parameter interfaceClSmrCharm ("ClSmrCharm", "cluster direction Gaussian smearing for c quark", - &ClusterDecayer::_clSmrCharm, 0, 0.0 , 0.0 , 2.0,false,false,false); -static Parameter + &ClusterDecayer::_clSmrCharm, 0, 0.0 , 0.0 , 2.0,false,false,false); +static Parameter interfaceClSmrExotic ("ClSmrExotic", "cluster direction Gaussian smearing for exotic quark", - &ClusterDecayer::_clSmrExotic, 0, 0.0 , 0.0 , 2.0,false,false,false); - + &ClusterDecayer::_clSmrExotic, 0, 0.0 , 0.0 , 2.0,false,false,false); + static Switch interfaceOnShell ("OnShell", "Whether or not the hadrons produced should by on shell or generated using the" " mass generator.", &ClusterDecayer::_onshell, false, false, false); static SwitchOption interfaceOnShellOnShell (interfaceOnShell, "Yes", "Produce the hadrons on shell", true); static SwitchOption interfaceOnShellOffShell (interfaceOnShell, "No", "Generate the masses using the mass generator.", false); static Parameter interfaceMassTry ("MassTry", "The number attempts to generate the masses of the hadrons produced" " in the cluster decay.", &ClusterDecayer::_masstry, 20, 1, 50, false, false, Interface::limited); } -void ClusterDecayer::decay(const ClusterVector & clusters, tPVector & finalhadrons) +void ClusterDecayer::decay(const ClusterVector & clusters, tPVector & finalhadrons) { // Loop over all clusters, and if they are not too heavy (that is - // intermediate clusters that have undergone to fission) or not - // too light (that is final clusters that have been already decayed + // intermediate clusters that have undergone to fission) or not + // too light (that is final clusters that have been already decayed // into single hadron) then decay them into two hadrons. + for (ClusterVector::const_iterator it = clusters.begin(); it != clusters.end(); ++it) { - if ((*it)->isAvailable() && !(*it)->isStatusFinal() - && (*it)->isReadyToDecay()) { + if ((*it)->isAvailable() && !(*it)->isStatusFinal() + && (*it)->isReadyToDecay()) { pair prod = decayIntoTwoHadrons(*it); (*it)->addChild(prod.first); (*it)->addChild(prod.second); finalhadrons.push_back(prod.first); finalhadrons.push_back(prod.second); } } } pair ClusterDecayer::decayIntoTwoHadrons(tClusterPtr ptr) { using Constants::pi; using Constants::twopi; // To decay the cluster into two hadrons one must distinguish between // constituent quarks (or diquarks) that originate from perturbative // processes (hard process or parton shower) from those that are generated // by the non-perturbative gluon splitting or from fission of heavy clusters. // In the latter case the two body decay is assumed to be *isotropic*. - // In the former case instead, if proper flag are activated, the two body - // decay is assumed to "remember" the direction of the constituents inside - // the cluster, in the cluster frame. The actual smearing of the hadron - // directions around the direction of the constituents, in the cluster + // In the former case instead, if proper flag are activated, the two body + // decay is assumed to "remember" the direction of the constituents inside + // the cluster, in the cluster frame. The actual smearing of the hadron + // directions around the direction of the constituents, in the cluster // frame, can be different between non-b hadrons and b-hadrons, but is given // by the same functional form: // cosThetaSmearing = 1 + smearFactor * log( rnd() ) // (repeated until cosThetaSmearing > -1) // where the factor smearFactor is different between b- and non-b hadrons. // - // We need to require (at least at the moment, maybe in the future we - // could change it) that the cluster has exactly two components. - // If this is not the case, then send a warning because it is not suppose + // We need to require (at least at the moment, maybe in the future we + // could change it) that the cluster has exactly two components. + // If this is not the case, then send a warning because it is not suppose // to happen, and then return. if ( ptr->numComponents() != 2 ) { generator()->logWarning( Exception("ClusterDecayer::decayIntoTwoHadrons " - "***Still cluster with not exactly 2 components*** ", + "***Still cluster with not exactly 2 components*** ", Exception::warning) ); return pair(); } - + // Extract the id and particle pointer of the two components of the cluster. tPPtr ptr1 = ptr->particle(0); tPPtr ptr2 = ptr->particle(1); tcPDPtr ptr1data = ptr1->dataPtr(); tcPDPtr ptr2data = ptr2->dataPtr(); - + bool isHad1FlavSpecial = false; bool cluDirHad1 = _clDirLight; double cluSmearHad1 = _clSmrLight; bool isHad2FlavSpecial = false; bool cluDirHad2 = _clDirLight; double cluSmearHad2 = _clSmrLight; if (CheckId::isExotic(ptr1data)) { isHad1FlavSpecial = true; cluDirHad1 = _clDirExotic; cluSmearHad1 = _clSmrExotic; - } + } else if (CheckId::hasBottom(ptr1data)) { isHad1FlavSpecial = true; cluDirHad1 = _clDirBottom; cluSmearHad1 = _clSmrBottom; - } + } else if (CheckId::hasCharm(ptr1data)) { isHad1FlavSpecial = true; cluDirHad1 = _clDirCharm; cluSmearHad1 = _clSmrCharm; - } + } if (CheckId::isExotic(ptr2data)) { isHad2FlavSpecial = true; cluDirHad2 = _clDirExotic; cluSmearHad2 = _clSmrExotic; - } + } else if (CheckId::hasBottom(ptr2data)) { isHad2FlavSpecial = true; cluDirHad2 = _clDirBottom; cluSmearHad2 = _clSmrBottom; - } + } else if (CheckId::hasCharm(ptr2data)) { isHad2FlavSpecial = true; cluDirHad2 = _clDirCharm; cluSmearHad2 = _clSmrCharm; - } + } bool isOrigin1Perturbative = ptr->isPerturbative(0); bool isOrigin2Perturbative = ptr->isPerturbative(1); - // We have to decide which, if any, of the two hadrons will have + // We have to decide which, if any, of the two hadrons will have // the momentum, in the cluster parent frame, smeared around the // direction of its constituent (for Had1 is the one pointed by // ptr1, and for Had2 is the one pointed by ptr2). // This happens only if the flag _ClDirX is 1 and the constituent is // perturbative (that is not coming from nonperturbative gluon splitting // or cluster fission). In the case that both the hadrons satisfy this // two requirements (of course only one must be treated, because the other - // one will have the momentum automatically fixed by the momentum + // one will have the momentum automatically fixed by the momentum // conservation) then more priority is given in the case of a b-hadron. // Finally, in the case that the two hadrons have same priority, then - // we choose randomly, with equal probability, one of the two. + // we choose randomly, with equal probability, one of the two. int priorityHad1 = 0; if ( cluDirHad1 && isOrigin1Perturbative ) { priorityHad1 = isHad1FlavSpecial ? 2 : 1; } int priorityHad2 = 0; if ( cluDirHad2 && isOrigin2Perturbative ) { priorityHad2 = isHad2FlavSpecial ? 2 : 1; } if ( priorityHad2 && priorityHad1 == priorityHad2 && UseRandom::rndbool() ) { priorityHad2 = 0; } Lorentz5Momentum pClu = ptr->momentum(); bool secondHad = false; Axis uSmear_v3; - if ( priorityHad1 || priorityHad2 ) { + if ( priorityHad1 || priorityHad2 ) { double cluSmear; Lorentz5Momentum pQ; if ( priorityHad1 > priorityHad2 ) { pQ = ptr1->momentum(); cluSmear = cluSmearHad1; - } else { + } else { pQ = ptr2->momentum(); cluSmear = cluSmearHad2; secondHad = true; } - + // To find the momenta of the two hadrons in the parent cluster frame // we proceed as follows. First, we determine the unit vector parallel // to the direction of the constituent in the cluster frame. Then we // have to smear such direction using the following prescription: // --- in theta angle w.r.t. such direction (not the z-axis), // the drawing of the cosine of such angle is done via: // 1.0 + cluSmear*log( rnd() ) // (repeated until it gives a value above -1.0) // --- in phi angle w.r.t. such direction, the drawing is simply flat. // Then, given the direction in the parent cluster frame of one of the // two hadrons, it is straighforward to get the momenta of both hadrons // (always in the same parent cluster frame). - pQ.boost( -pClu.boostVector() ); // boost from Lab to Cluster frame + pQ.boost( -pClu.boostVector() ); // boost from Lab to Cluster frame uSmear_v3 = pQ.vect().unit(); // skip if cluSmear is too small if ( cluSmear > 0.001 ) { - // generate the smearing angle + // generate the smearing angle double cosSmear; do cosSmear = 1.0 + cluSmear*log( UseRandom::rnd() ); while ( cosSmear < -1.0 ); double sinSmear = sqrt( 1.0 - sqr(cosSmear) ); // calculate rotation to z axis LorentzRotation rot; double sinth(sqrt(1.-sqr(uSmear_v3.z()))); if(abs(uSmear_v3.perp2()/uSmear_v3.z())>1e-10) rot.setRotate(-acos(uSmear_v3.z()), Axis(-uSmear_v3.y()/sinth,uSmear_v3.x()/sinth,0.)); // + random azimuthal rotation rot.rotateZ(UseRandom::rnd()*twopi); // set direction in rotated frame Lorentz5Vector ltemp(0.,sinSmear,cosSmear,0.); // rotate back rot.invert(); ltemp *= rot; uSmear_v3 = ltemp.vect(); } - } + } else { - - // Isotropic decay: flat in cosTheta and phi. + + // Isotropic decay: flat in cosTheta and phi. uSmear_v3 = Axis(1.0, 0.0, 0.0); // just to set the rho to 1 uSmear_v3.setTheta( acos( UseRandom::rnd( -1.0 , 1.0 ) ) ); - uSmear_v3.setPhi( UseRandom::rnd( -pi , pi ) ); - + uSmear_v3.setPhi( UseRandom::rnd( -pi , pi ) ); + } - pair dataPair + pair dataPair = _hadronsSelector->chooseHadronPair(ptr->mass(), ptr1data, ptr2data); - if(dataPair.first == tcPDPtr() || + if(dataPair.first == tcPDPtr() || dataPair.second == tcPDPtr()) return pair(); + /*ofstream afile("mCDLEP.txt", std::ios_base::app | std::ios_base::out); + afile<mass()/GeV<<"\n"; + afile.close(); + Lorentz5Momentum mom = ptr1->momentum() + ptr2->momentum(); + Energy pT = abs(sqrt( sqr(mom.x()) + sqr(mom.y()) ) ); + ofstream bfile("pTCDLEP.txt", std::ios_base::app | std::ios_base::out); + bfile<produceParticle(dataPair.first ->mass()); ptrHad2 = dataPair.second->produceParticle(dataPair.second->mass()); } // produce the hadrons with mass given by the mass generator else { unsigned int ntry(0); do { ptrHad1 = dataPair.first ->produceParticle(); ptrHad2 = dataPair.second->produceParticle(); ++ntry; } while(ntry<_masstry&&ptrHad1->mass()+ptrHad2->mass()>ptr->mass()); // if fails produce on shell and issue warning (should never happen??) if( ptrHad1->mass() + ptrHad2->mass() > ptr->mass() ) { generator()->log() << "Failed to produce off-shell hadrons in " << "ClusterDecayer::decayIntoTwoHadrons producing hadrons " << "on-shell" << endl; ptrHad1 = dataPair.first ->produceParticle(dataPair.first ->mass()); ptrHad2 = dataPair.second->produceParticle(dataPair.second->mass()); } } - + if (!ptrHad1 || !ptrHad2) { ostringstream s; s << "ClusterDecayer::decayIntoTwoHadrons ***Cannot create the two hadrons***" - << dataPair.first ->PDGName() << " and " + << dataPair.first ->PDGName() << " and " << dataPair.second->PDGName(); cerr << s.str() << endl; generator()->logWarning( Exception(s.str(), Exception::warning) ); } else { Lorentz5Momentum pHad1, pHad2; // 5-momentum vectors that we have to determine if ( secondHad ) uSmear_v3 *= -1.0; if (pClu.m() < ptrHad1->mass()+ptrHad2->mass() ) { - throw Exception() << "Impossible Kinematics in ClusterDecayer::decayIntoTwoHadrons()" + throw Exception() << "Impossible Kinematics in ClusterDecayer::decayIntoTwoHadrons()" << Exception::eventerror; } Kinematics::twoBodyDecay(pClu,ptrHad1->mass(),ptrHad2->mass(),uSmear_v3, pHad1,pHad2); ptrHad1->set5Momentum(pHad1); ptrHad2->set5Momentum(pHad2); // Determine the positions of the two children clusters. LorentzPoint positionHad1 = LorentzPoint(); LorentzPoint positionHad2 = LorentzPoint(); calculatePositions(pClu, ptr->vertex(), pHad1, pHad2, positionHad1, positionHad2); ptrHad1->setVertex(positionHad1); ptrHad2->setVertex(positionHad2); } return pair(ptrHad1,ptrHad2); } void ClusterDecayer:: -calculatePositions(const Lorentz5Momentum &pClu, - const LorentzPoint &positionClu, - const Lorentz5Momentum &, - const Lorentz5Momentum &, - LorentzPoint &positionHad1, +calculatePositions(const Lorentz5Momentum &pClu, + const LorentzPoint &positionClu, + const Lorentz5Momentum &, + const Lorentz5Momentum &, + LorentzPoint &positionHad1, LorentzPoint &positionHad2 ) const { // First, determine the relative positions of the children hadrons // with respect to their parent cluster, in the cluster reference frame, - // assuming gaussian smearing with width inversely proportional to the + // assuming gaussian smearing with width inversely proportional to the // parent cluster mass. Length smearingWidth = hbarc / pClu.m(); LorentzDistance distanceHad[2]; for ( int i = 0; i < 2; i++ ) { // i==0 is the first hadron; i==1 is the second one Length delta[4]={ZERO,ZERO,ZERO,ZERO}; // smearing of the four components of the LorentzDistance, two at the same time to improve speed for ( int j = 0; j < 3; j += 2 ) { delta[j] = UseRandom::rndGauss(smearingWidth, Length(ZERO)); delta[j+1] = UseRandom::rndGauss(smearingWidth, Length(ZERO)); } // set the distance delta[0] = abs(delta[0]) +sqrt(sqr(delta[1])+sqr(delta[2])+sqr(delta[3])); distanceHad[i] = LorentzDistance(delta[1],delta[2],delta[3],delta[0]); // Boost such relative positions of the children hadrons, // with respect to their parent cluster, // from the cluster reference frame to the Lab frame. distanceHad[i].boost(pClu.boostVector()); } // Finally, determine the absolute positions of the children hadrons // in the Lab frame. positionHad1 = distanceHad[0] + positionClu; positionHad2 = distanceHad[1] + positionClu; } - diff --git a/Hadronization/ClusterDecayer.h b/Hadronization/ClusterDecayer.h --- a/Hadronization/ClusterDecayer.h +++ b/Hadronization/ClusterDecayer.h @@ -1,166 +1,167 @@ // -*- C++ -*- // // ClusterDecayer.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_ClusterDecayer_H #define HERWIG_ClusterDecayer_H #include #include #include "CluHadConfig.h" #include "HadronSelector.h" #include "ClusterDecayer.fh" namespace Herwig { using namespace ThePEG; /** \ingroup Hadronization * \class ClusterDecayer * \brief This class decays the "normal" clusters * \author Philip Stephens * \author Alberto Ribon * * This class decays the "normal" clusters, e.g. ones that are not heavy - * enough for fission, and not too light to decay into one hadron. + * enough for fission, and not too light to decay into one hadron. * - * This class is directs the production of hadrons via 2-body cluster decays. + * This class is directs the production of hadrons via 2-body cluster decays. * The selection of the hadron flavours is given by Herwig::HadronSelector. * - * @see HadronSelector + * @see HadronSelector * @see \ref ClusterDecayerInterfaces "The interfaces" * defined for ClusterDecayer. */ class ClusterDecayer: public Interfaced { public: /** @name Standard constructors and destructors. */ //@{ /** * Default constructor. */ - ClusterDecayer(); + ClusterDecayer(); //@} - /** Decays all remaining clusters into hadrons. + /** Decays all remaining clusters into hadrons. * This routine decays the clusters that are left after * Herwig::ClusterFissioner::fission and * Herwig::LightClusterDecayer::decay have been called. These are all * the "normal" clusters which are not forced into hadrons by * the other functions. */ - void decay(const ClusterVector & clusters, tPVector & finalhadrons) + void decay(const ClusterVector & clusters, tPVector & finalhadrons) ; public: /** * Standard ThePEG function for writing a persistent stream. */ void persistentOutput(PersistentOStream &) const; /** * Standard ThePEG function for reading from a persistent stream. */ void persistentInput(PersistentIStream &, int); /** * Standard Init function used to initialize the interfaces. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} private: /** * Private and non-existent assignment operator. */ ClusterDecayer & operator=(const ClusterDecayer &) = delete; public: - /** Decays the cluster into two hadrons. + /** Decays the cluster into two hadrons. * * This routine is used to take a given cluster and decay it into * two hadrons which are returned. If one of the constituents is from * the perturbative regime then the direction of the perturbative parton - * is remembered and the decay is preferentially in that direction. The + * is remembered and the decay is preferentially in that direction. The * direction of the decay is given by * \f[ \cos \theta = 1 + S \log r_1 \f] * where \f$ S \f$ is a parameter of the model and \f$ r_1 \f$ is a random * number [0,1]. */ pair decayIntoTwoHadrons(tClusterPtr ptr); private: /** Compute the positions of the new hadrons based on the clusters position. * * This method calculates the positions of the children hadrons by a - * call to ThePEG::RandomGenerator::rndGaussTwoNumbers with width inversely + * call to ThePEG::RandomGenerator::rndGaussTwoNumbers with width inversely * proportional to the cluster mass, around the parent cluster position. */ - void calculatePositions( const Lorentz5Momentum &, const LorentzPoint &, + void calculatePositions( const Lorentz5Momentum &, const LorentzPoint &, const Lorentz5Momentum &, const Lorentz5Momentum &, LorentzPoint &, LorentzPoint &) const; /** * Pointer to a Herwig::HadronSelector for choosing decay types */ Ptr::pointer _hadronsSelector; - //@{ + //@{ /** * Whether a cluster decays along the perturbative parton direction. */ bool _clDirLight; bool _clDirBottom; bool _clDirCharm; bool _clDirExotic; /** * The S parameter from decayIntoTwoHadrons */ double _clSmrLight; double _clSmrBottom; double _clSmrCharm; double _clSmrExotic; //@} /** * Whether or not the hadrons produced should be on-shell * or generated used the MassGenerator */ bool _onshell; /** * Number of tries to generate the masses of the decay products */ unsigned int _masstry; + }; } #endif /* HERWIG_ClusterDecayer_H */ diff --git a/Hadronization/ClusterFissioner.cc b/Hadronization/ClusterFissioner.cc --- a/Hadronization/ClusterFissioner.cc +++ b/Hadronization/ClusterFissioner.cc @@ -1,943 +1,1110 @@ // -*- C++ -*- // // ClusterFissioner.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // Thisk is the implementation of the non-inlined, non-templated member // functions of the ClusterFissioner class. // #include "ClusterFissioner.h" #include #include #include #include #include #include #include #include "Herwig/Utilities/Kinematics.h" #include "CheckId.h" #include "Cluster.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include using namespace Herwig; DescribeClass describeClusterFissioner("Herwig::ClusterFissioner",""); ClusterFissioner::ClusterFissioner() : _clMaxLight(3.35*GeV), _clMaxBottom(3.35*GeV), _clMaxCharm(3.35*GeV), _clMaxExotic(3.35*GeV), _clPowLight(2.0), _clPowBottom(2.0), _clPowCharm(2.0), _clPowExotic(2.0), _pSplitLight(1.0), _pSplitBottom(1.0), _pSplitCharm(1.0), _pSplitExotic(1.0), + _fissionCluster(0), + _fissionPwtUquark(1), + _fissionPwtDquark(1), + _fissionPwtSquark(0.5), _btClM(1.0*GeV), _iopRem(1), - _kappa(1.0e15*GeV/meter) + _kappa(1.0e15*GeV/meter), + _enhanceSProb(0), + _m0Fission(2.*GeV), + _massMeasure(0) {} IBPtr ClusterFissioner::clone() const { return new_ptr(*this); } IBPtr ClusterFissioner::fullclone() const { return new_ptr(*this); } void ClusterFissioner::persistentOutput(PersistentOStream & os) const { os << _hadronsSelector << ounit(_clMaxLight,GeV) << ounit(_clMaxBottom,GeV) << ounit(_clMaxCharm,GeV) << ounit(_clMaxExotic,GeV) << _clPowLight << _clPowBottom - << _clPowCharm << _clPowExotic << _pSplitLight - << _pSplitBottom << _pSplitCharm << _pSplitExotic << ounit(_btClM,GeV) - << _iopRem << ounit(_kappa, GeV/meter); + << _clPowCharm << _clPowExotic << _pSplitLight + << _pSplitBottom << _pSplitCharm << _pSplitExotic + << _fissionCluster << _fissionPwtUquark << _fissionPwtDquark << _fissionPwtSquark + << ounit(_btClM,GeV) + << _iopRem << ounit(_kappa, GeV/meter) + << _enhanceSProb << ounit(_m0Fission,GeV) << _massMeasure; } void ClusterFissioner::persistentInput(PersistentIStream & is, int) { is >> _hadronsSelector >> iunit(_clMaxLight,GeV) >> iunit(_clMaxBottom,GeV) >> iunit(_clMaxCharm,GeV) >> iunit(_clMaxExotic,GeV) >> _clPowLight >> _clPowBottom - >> _clPowCharm >> _clPowExotic >> _pSplitLight + >> _clPowCharm >> _clPowExotic >> _pSplitLight >> _pSplitBottom >> _pSplitCharm >> _pSplitExotic + >> _fissionCluster >> _fissionPwtUquark >> _fissionPwtDquark >> _fissionPwtSquark >> iunit(_btClM,GeV) >> _iopRem - >> iunit(_kappa, GeV/meter); + >> iunit(_kappa, GeV/meter) + >> _enhanceSProb >> iunit(_m0Fission,GeV) >> _massMeasure; } void ClusterFissioner::Init() { static ClassDocumentation documentation ("Class responsibles for chopping up the clusters"); - static Reference - interfaceHadronSelector("HadronSelector", - "A reference to the HadronSelector object", + static Reference + interfaceHadronSelector("HadronSelector", + "A reference to the HadronSelector object", &Herwig::ClusterFissioner::_hadronsSelector, false, false, true, false); - + // ClMax for light, Bottom, Charm and exotic (e.g. Susy) quarks static Parameter interfaceClMaxLight ("ClMaxLight","cluster max mass for light quarks (unit [GeV])", &ClusterFissioner::_clMaxLight, GeV, 3.35*GeV, ZERO, 10.0*GeV, false,false,false); static Parameter interfaceClMaxBottom ("ClMaxBottom","cluster max mass for b quarks (unit [GeV])", &ClusterFissioner::_clMaxBottom, GeV, 3.35*GeV, ZERO, 10.0*GeV, false,false,false); static Parameter interfaceClMaxCharm ("ClMaxCharm","cluster max mass for c quarks (unit [GeV])", &ClusterFissioner::_clMaxCharm, GeV, 3.35*GeV, ZERO, 10.0*GeV, false,false,false); static Parameter interfaceClMaxExotic ("ClMaxExotic","cluster max mass for exotic quarks (unit [GeV])", &ClusterFissioner::_clMaxExotic, GeV, 3.35*GeV, ZERO, 10.0*GeV, false,false,false); - + // ClPow for light, Bottom, Charm and exotic (e.g. Susy) quarks static Parameter interfaceClPowLight ("ClPowLight","cluster mass exponent for light quarks", &ClusterFissioner::_clPowLight, 0, 2.0, 0.0, 10.0,false,false,false); static Parameter interfaceClPowBottom ("ClPowBottom","cluster mass exponent for b quarks", &ClusterFissioner::_clPowBottom, 0, 2.0, 0.0, 10.0,false,false,false); static Parameter interfaceClPowCharm ("ClPowCharm","cluster mass exponent for c quarks", &ClusterFissioner::_clPowCharm, 0, 2.0, 0.0, 10.0,false,false,false); static Parameter interfaceClPowExotic ("ClPowExotic","cluster mass exponent for exotic quarks", &ClusterFissioner::_clPowExotic, 0, 2.0, 0.0, 10.0,false,false,false); // PSplit for light, Bottom, Charm and exotic (e.g. Susy) quarks static Parameter interfacePSplitLight ("PSplitLight","cluster mass splitting param for light quarks", &ClusterFissioner::_pSplitLight, 0, 1.0, 0.0, 10.0,false,false,false); static Parameter interfacePSplitBottom ("PSplitBottom","cluster mass splitting param for b quarks", &ClusterFissioner::_pSplitBottom, 0, 1.0, 0.0, 10.0,false,false,false); static Parameter interfacePSplitCharm ("PSplitCharm","cluster mass splitting param for c quarks", &ClusterFissioner::_pSplitCharm, 0, 1.0, 0.0, 10.0,false,false,false); static Parameter interfacePSplitExotic ("PSplitExotic","cluster mass splitting param for exotic quarks", &ClusterFissioner::_pSplitExotic, 0, 1.0, 0.0, 10.0,false,false,false); + static Switch interfaceFission + ("Fission", + "Option for different Fission options", + &ClusterFissioner::_fissionCluster, 1, false, false); + static SwitchOption interfaceFissionDefault + (interfaceFission, + "default", + "Normal cluster fission which depends on the hadron selector class.", + 0); + static SwitchOption interfaceFissionNew + (interfaceFission, + "new", + "Alternative cluster fission which does not depend on the hadron selector class", + 1); + + + static Parameter interfaceFissionPwtUquark + ("FissionPwtUquark", + "Weight for fission in U quarks", + &ClusterFissioner::_fissionPwtUquark, 1, 0.0, 1.0, + false, false, Interface::limited); + static Parameter interfaceFissionPwtDquark + ("FissionPwtDquark", + "Weight for fission in D quarks", + &ClusterFissioner::_fissionPwtDquark, 1, 0.0, 1.0, + false, false, Interface::limited); + static Parameter interfaceFissionPwtSquark + ("FissionPwtSquark", + "Weight for fission in S quarks", + &ClusterFissioner::_fissionPwtSquark, 0.5, 0.0, 1.0, + false, false, Interface::limited); + + static Switch interfaceRemnantOption ("RemnantOption", "Option for the treatment of remnant clusters", &ClusterFissioner::_iopRem, 1, false, false); static SwitchOption interfaceRemnantOptionSoft (interfaceRemnantOption, "Soft", "Both clusters produced in the fission of the beam cluster" " are treated as soft clusters.", 0); static SwitchOption interfaceRemnantOptionHard (interfaceRemnantOption, "Hard", "Only the cluster containing the remnant is treated as a soft cluster.", 1); static SwitchOption interfaceRemnantOptionVeryHard (interfaceRemnantOption, "VeryHard", "Even remnant clusters are treated as hard, i.e. all clusters the same", 2); static Parameter interfaceBTCLM ("SoftClusterFactor", "Parameter for the mass spectrum of remnant clusters", &ClusterFissioner::_btClM, GeV, 1.*GeV, 0.1*GeV, 10.0*GeV, false, false, Interface::limited); - + static Parameter interfaceStringTension ("StringTension", "String tension used in vertex displacement calculation", - &ClusterFissioner::_kappa, GeV/meter, + &ClusterFissioner::_kappa, GeV/meter, 1.0e15*GeV/meter, ZERO, ZERO, false, false, Interface::lowerlim); + static Switch interfaceEnhanceSProb + ("EnhanceSProb", + "Option for enhancing strangeness", + &ClusterFissioner::_enhanceSProb, 0, false, false); + static SwitchOption interfaceEnhanceSProbNo + (interfaceEnhanceSProb, + "No", + "No strangeness enhancement.", + 0); + static SwitchOption interfaceEnhanceSProbScaled + (interfaceEnhanceSProb, + "Scaled", + "Scaled strangeness enhancement", + 1); + static SwitchOption interfaceEnhanceSProbExponential + (interfaceEnhanceSProb, + "Exponential", + "Exponential strangeness enhancement", + 2); + + static Switch interfaceMassMeasure + ("MassMeasure", + "Option to use different mass measures", + &ClusterFissioner::_massMeasure,0,false,false); + static SwitchOption interfaceMassMeasureMass + (interfaceMassMeasure, + "Mass", + "Mass Measure", + 0); + static SwitchOption interfaceMassMeasureLambda + (interfaceMassMeasure, + "Lambda", + "Lambda Measure", + 1); + + static Parameter interfaceFissionMassScale + ("FissionMassScale", + "Cluster fission mass scale", + &ClusterFissioner::_m0Fission, GeV, 2.0*GeV, 0.1*GeV, 50.*GeV, + false, false, Interface::limited); + } tPVector ClusterFissioner::fission(ClusterVector & clusters, bool softUEisOn) { // return if no clusters if (clusters.empty()) return tPVector(); /***************** - * Loop over the (input) collection of cluster pointers, and store in + * Loop over the (input) collection of cluster pointers, and store in * the vector splitClusters all the clusters that need to be split - * (these are beam clusters, if soft underlying event is off, and + * (these are beam clusters, if soft underlying event is off, and * heavy non-beam clusters). ********************/ - - stack splitClusters; - for(ClusterVector::iterator it = clusters.begin() ; + + stack splitClusters; + for(ClusterVector::iterator it = clusters.begin() ; it != clusters.end() ; ++it) { - /************** - * Skip 3-component clusters that have been redefined (as 2-component - * clusters) or not available clusters. The latter check is indeed - * redundant now, but it is used for possible future extensions in which, + * Skip 3-component clusters that have been redefined (as 2-component + * clusters) or not available clusters. The latter check is indeed + * redundant now, but it is used for possible future extensions in which, * for some reasons, some of the clusters found by ClusterFinder are tagged * straight away as not available. **************/ if((*it)->isRedefined() || !(*it)->isAvailable()) continue; // if the cluster is a beam cluster add it to the vector of clusters // to be split or if it is heavy if((*it)->isBeamCluster() || isHeavy(*it)) splitClusters.push(*it); - } + } tPVector finalhadrons; cut(splitClusters, clusters, finalhadrons, softUEisOn); return finalhadrons; } -void ClusterFissioner::cut(stack & clusterStack, - ClusterVector &clusters, tPVector & finalhadrons, +void ClusterFissioner::cut(stack & clusterStack, + ClusterVector &clusters, tPVector & finalhadrons, bool softUEisOn) { /************************************************** * This method does the splitting of the cluster pointed by cluPtr * and "recursively" by all of its cluster children, if heavy. All of these * new children clusters are added (indeed the pointers to them) to the * collection of cluster pointers collecCluPtr. The method works as follows. * Initially the vector vecCluPtr contains just the input pointer to the * cluster to be split. Then it will be filled "recursively" by all * of the cluster's children that are heavy enough to require, in their turn, - * to be split. In each loop, the last element of the vector vecCluPtr is + * to be split. In each loop, the last element of the vector vecCluPtr is * considered (only once because it is then removed from the vector). * This approach is conceptually recursive, but avoid the overhead of * a concrete recursive function. Furthermore it requires minimal changes * in the case that the fission of an heavy cluster could produce more - * than two cluster children as assumed now. + * than two cluster children as assumed now. * * Draw the masses: for normal, non-beam clusters a power-like mass dist - * is used, whereas for beam clusters a fast-decreasing exponential mass - * dist is used instead (to avoid many iterative splitting which could + * is used, whereas for beam clusters a fast-decreasing exponential mass + * dist is used instead (to avoid many iterative splitting which could * produce an unphysical large transverse energy from a supposed soft beam * remnant process). ****************************************/ // Here we recursively loop over clusters in the stack and cut them while (!clusterStack.empty()) { // take the last element of the vector - ClusterPtr iCluster = clusterStack.top(); clusterStack.pop(); + ClusterPtr iCluster = clusterStack.top(); clusterStack.pop(); // split it cutType ct = iCluster->numComponents() == 2 ? - cutTwo(iCluster, finalhadrons, softUEisOn) : + cutTwo(iCluster, finalhadrons, softUEisOn) : cutThree(iCluster, finalhadrons, softUEisOn); // There are cases when we don't want to split, even if it fails mass test if(!ct.first.first || !ct.second.first) { // if an unsplit beam cluster leave if for the underlying event if(iCluster->isBeamCluster() && softUEisOn) iCluster->isAvailable(false); continue; } // check if clusters ClusterPtr one = dynamic_ptr_cast(ct.first.first); ClusterPtr two = dynamic_ptr_cast(ct.second.first); // is a beam cluster must be split into two clusters if(iCluster->isBeamCluster() && (!one||!two) && softUEisOn) { iCluster->isAvailable(false); continue; } // There should always be a intermediate quark(s) from the splitting assert(ct.first.second && ct.second.second); /// \todo sort out motherless quark pairs here. Watch out for 'quark in final state' errors iCluster->addChild(ct.first.first); // iCluster->addChild(ct.first.second); // ct.first.second->addChild(ct.first.first); - + iCluster->addChild(ct.second.first); // iCluster->addChild(ct.second.second); // ct.second.second->addChild(ct.second.first); - + // Sometimes the clusters decay C -> H + C' rather then C -> C' + C'' if(one) { clusters.push_back(one); if(one->isBeamCluster() && softUEisOn) one->isAvailable(false); if(isHeavy(one) && one->isAvailable()) clusterStack.push(one); } if(two) { clusters.push_back(two); if(two->isBeamCluster() && softUEisOn) two->isAvailable(false); if(isHeavy(two) && two->isAvailable()) clusterStack.push(two); } } } namespace { /** * Check if can't make a hadron from the partons */ bool cantMakeHadron(tcPPtr p1, tcPPtr p2) { return ! CheckId::canBeHadron(p1->dataPtr(), p2->dataPtr()); } - + /** * Check if can't make a diquark from the partons */ bool cantMakeDiQuark(tcPPtr p1, tcPPtr p2) { long id1 = p1->id(), id2 = p2->id(); return ! (QuarkMatcher::Check(id1) && QuarkMatcher::Check(id2) && id1*id2>0); } } ClusterFissioner::cutType ClusterFissioner::cutTwo(ClusterPtr & cluster, tPVector & finalhadrons, bool softUEisOn) { // need to make sure only 2-cpt clusters get here assert(cluster->numComponents() == 2); tPPtr ptrQ1 = cluster->particle(0); tPPtr ptrQ2 = cluster->particle(1); Energy Mc = cluster->mass(); assert(ptrQ1); assert(ptrQ2); - + // And check if those particles are from a beam remnant bool rem1 = cluster->isBeamRemnant(0); bool rem2 = cluster->isBeamRemnant(1); // workout which distribution to use bool soft1(false),soft2(false); switch (_iopRem) { case 0: soft1 = rem1 || rem2; soft2 = rem2 || rem1; break; case 1: soft1 = rem1; soft2 = rem2; break; } // Initialization for the exponential ("soft") mass distribution. static const int max_loop = 1000; int counter = 0; Energy Mc1 = ZERO, Mc2 = ZERO,m1=ZERO,m2=ZERO,m=ZERO; tcPDPtr toHadron1, toHadron2; PPtr newPtr1 = PPtr (); PPtr newPtr2 = PPtr (); bool succeeded = false; do { succeeded = false; ++counter; - - drawNewFlavour(newPtr1,newPtr2); + if (_enhanceSProb == 0){ + drawNewFlavour(newPtr1,newPtr2); + } + else { + Energy2 mass2 = clustermass(cluster); + drawNewFlavourEnhanced(newPtr1,newPtr2,mass2); + } // check for right ordering assert (ptrQ2); assert (newPtr2); assert (ptrQ2->dataPtr()); assert (newPtr2->dataPtr()); if(cantMakeHadron(ptrQ1, newPtr1) || cantMakeHadron(ptrQ2, newPtr2)) { swap(newPtr1, newPtr2); // check again if(cantMakeHadron(ptrQ1, newPtr1) || cantMakeHadron(ptrQ2, newPtr2)) { - throw Exception() + throw Exception() << "ClusterFissioner cannot split the cluster (" << ptrQ1->PDGName() << ' ' << ptrQ2->PDGName() << ") into hadrons.\n" << Exception::runerror; } } // Check that new clusters can produce particles and there is enough // phase space to choose the drawn flavour m1 = ptrQ1->data().constituentMass(); m2 = ptrQ2->data().constituentMass(); m = newPtr1->data().constituentMass(); // Do not split in the case there is no phase space available if(Mc < m1+m + m2+m) continue; // power for splitting double exp1=_pSplitLight; double exp2=_pSplitLight; if (CheckId::isExotic(ptrQ1->dataPtr())) exp1 = _pSplitExotic; else if(CheckId::hasBottom(ptrQ1->dataPtr()))exp1 = _pSplitBottom; else if(CheckId::hasCharm(ptrQ1->dataPtr())) exp1 = _pSplitCharm; if (CheckId::isExotic(ptrQ2->dataPtr())) exp2 = _pSplitExotic; else if(CheckId::hasBottom(ptrQ2->dataPtr())) exp2 = _pSplitBottom; else if(CheckId::hasCharm(ptrQ2->dataPtr())) exp2 = _pSplitCharm; - // If, during the drawing of candidate masses, too many attempts fail + // If, during the drawing of candidate masses, too many attempts fail // (because the phase space available is tiny) /// \todo run separate loop here? Mc1 = drawChildMass(Mc,m1,m2,m,exp1,soft1); Mc2 = drawChildMass(Mc,m2,m1,m,exp2,soft2); if(Mc1 < m1+m || Mc2 < m+m2 || Mc1+Mc2 > Mc) continue; /************************** * New (not present in Fortran Herwig): - * check whether the fragment masses Mc1 and Mc2 are above the - * threshold for the production of the lightest pair of hadrons with the - * right flavours. If not, then set by hand the mass to the lightest + * check whether the fragment masses Mc1 and Mc2 are above the + * threshold for the production of the lightest pair of hadrons with the + * right flavours. If not, then set by hand the mass to the lightest * single hadron with the right flavours, in order to solve correctly * the kinematics, and (later in this method) create directly such hadron * and add it to the children hadrons of the cluster that undergoes the * fission (i.e. the one pointed by iCluPtr). Notice that in this special - * case, the heavy cluster that undergoes the fission has one single + * case, the heavy cluster that undergoes the fission has one single * cluster child and one single hadron child. We prefer this approach, * rather than to create a light cluster, with the mass set equal to - * the lightest hadron, and let then the class LightClusterDecayer to do - * the job to decay it to that single hadron, for two reasons: - * First, because the sum of the masses of the two constituents can be, + * the lightest hadron, and let then the class LightClusterDecayer to do + * the job to decay it to that single hadron, for two reasons: + * First, because the sum of the masses of the two constituents can be, * in this case, greater than the mass of that hadron, hence it would * be impossible to solve the kinematics for such two components, and * therefore we would have a cluster whose components are undefined. * Second, the algorithm is faster, because it avoids the reshuffling * procedure that would be necessary if we used LightClusterDecayer - * to decay the light cluster to the lightest hadron. + * to decay the light cluster to the lightest hadron. ****************************/ toHadron1 = _hadronsSelector->chooseSingleHadron(ptrQ1->dataPtr(), newPtr1->dataPtr(),Mc1); if(toHadron1) Mc1 = toHadron1->mass(); toHadron2 = _hadronsSelector->chooseSingleHadron(ptrQ2->dataPtr(), newPtr2->dataPtr(),Mc2); if(toHadron2) Mc2 = toHadron2->mass(); // if a beam cluster not allowed to decay to hadrons if(cluster->isBeamCluster() && (toHadron1||toHadron2) && softUEisOn) continue; - // Check if the decay kinematics is still possible: if not then + // Check if the decay kinematics is still possible: if not then // force the one-hadron decay for the other cluster as well. if(Mc1 + Mc2 > Mc) { if(!toHadron1) { toHadron1 = _hadronsSelector->chooseSingleHadron(ptrQ1->dataPtr(), newPtr1->dataPtr(),Mc-Mc2); if(toHadron1) Mc1 = toHadron1->mass(); - } + } else if(!toHadron2) { toHadron2 = _hadronsSelector->chooseSingleHadron(ptrQ2->dataPtr(), newPtr2->dataPtr(),Mc-Mc1); if(toHadron2) Mc2 = toHadron2->mass(); } } succeeded = (Mc >= Mc1+Mc2); } while (!succeeded && counter < max_loop); if(counter >= max_loop) { static const PPtr null = PPtr(); return cutType(PPair(null,null),PPair(null,null)); } // Determined the (5-components) momenta (all in the LAB frame) Lorentz5Momentum pClu = cluster->momentum(); // known Lorentz5Momentum p0Q1 = ptrQ1->momentum(); // known (mom Q1 before fission) Lorentz5Momentum pClu1, pClu2, pQ1, pQone, pQtwo, pQ2; //unknown pClu1.setMass(Mc1); pClu2.setMass(Mc2); pQ1.setMass(m1); pQ2.setMass(m2); - pQone.setMass(m); + pQone.setMass(m); pQtwo.setMass(m); calculateKinematics(pClu,p0Q1,toHadron1,toHadron2, pClu1,pClu2,pQ1,pQone,pQtwo,pQ2); // out /****************** * The previous methods have determined the kinematics and positions - * of C -> C1 + C2. + * of C -> C1 + C2. * In the case that one of the two product is light, that means either * decayOneHadronClu1 or decayOneHadronClu2 is true, then the momenta * of the components of that light product have not been determined, * and a (light) cluster will not be created: the heavy father cluster * decays, in this case, into a single (not-light) cluster and a * single hadron. In the other, "normal", cases the father cluster * decays into two clusters, each of which has well defined components. * Notice that, in the case of components which point to particles, the * momenta of the components is properly set to the new values, whereas - * we do not change the momenta of the pointed particles, because we + * we do not change the momenta of the pointed particles, because we * want to keep all of the information (that is the new momentum of a * component after the splitting, which is contained in the _momentum * member of the Component class, and the (old) momentum of that component * before the splitting, which is contained in the momentum of the * pointed particle). Please not make confusion of this only apparent * inconsistency! ********************/ LorentzPoint posC,pos1,pos2; posC = cluster->vertex(); calculatePositions(pClu, posC, pClu1, pClu2, pos1, pos2); cutType rval; if(toHadron1) { rval.first = produceHadron(toHadron1, newPtr1, pClu1, pos1); finalhadrons.push_back(rval.first.first); } else { rval.first = produceCluster(ptrQ1, newPtr1, pClu1, pos1, pQ1, pQone, rem1); } if(toHadron2) { rval.second = produceHadron(toHadron2, newPtr2, pClu2, pos2); finalhadrons.push_back(rval.second.first); } else { rval.second = produceCluster(ptrQ2, newPtr2, pClu2, pos2, pQ2, pQtwo, rem2); } return rval; } ClusterFissioner::cutType ClusterFissioner::cutThree(ClusterPtr & cluster, tPVector & finalhadrons, bool softUEisOn) { // need to make sure only 3-cpt clusters get here assert(cluster->numComponents() == 3); // extract quarks tPPtr ptrQ[3] = {cluster->particle(0),cluster->particle(1),cluster->particle(2)}; assert( ptrQ[0] && ptrQ[1] && ptrQ[2] ); // find maximum mass pair Energy mmax(ZERO); Lorentz5Momentum pDiQuark; int iq1(-1),iq2(-1); Lorentz5Momentum psum; for(int q1=0;q1<3;++q1) { psum+= ptrQ[q1]->momentum(); for(int q2=q1+1;q2<3;++q2) { Lorentz5Momentum ptest = ptrQ[q1]->momentum()+ptrQ[q2]->momentum(); ptest.rescaleMass(); Energy mass = ptest.m(); if(mass>mmax) { mmax = mass; pDiQuark = ptest; iq1 = q1; iq2 = q2; } } } // and the spectators int iother(-1); for(int ix=0;ix<3;++ix) if(ix!=iq1&&ix!=iq2) iother=ix; assert(iq1>=0&&iq2>=0&&iother>=0); // And check if those particles are from a beam remnant bool rem1 = cluster->isBeamRemnant(iq1); bool rem2 = cluster->isBeamRemnant(iq2); // workout which distribution to use bool soft1(false),soft2(false); switch (_iopRem) { case 0: soft1 = rem1 || rem2; soft2 = rem2 || rem1; break; case 1: soft1 = rem1; soft2 = rem2; break; } // Initialization for the exponential ("soft") mass distribution. static const int max_loop = 1000; int counter = 0; Energy Mc1 = ZERO, Mc2 = ZERO, m1=ZERO, m2=ZERO, m=ZERO; tcPDPtr toHadron; bool toDiQuark(false); PPtr newPtr1 = PPtr(),newPtr2 = PPtr(); PDPtr diquark; bool succeeded = false; do { succeeded = false; ++counter; - drawNewFlavour(newPtr1,newPtr2); + if (_enhanceSProb == 0) { + drawNewFlavour(newPtr1,newPtr2); + } + else { + Energy2 mass2 = clustermass(cluster); + drawNewFlavourEnhanced(newPtr1,newPtr2, mass2); + } // randomly pick which will be (anti)diquark and which a mesonic cluster if(UseRandom::rndbool()) { swap(iq1,iq2); swap(rem1,rem2); } // check first order if(cantMakeHadron(ptrQ[iq1], newPtr1) || cantMakeDiQuark(ptrQ[iq2], newPtr2)) { swap(newPtr1,newPtr2); } // check again if(cantMakeHadron(ptrQ[iq1], newPtr1) || cantMakeDiQuark(ptrQ[iq2], newPtr2)) { throw Exception() << "ClusterFissioner cannot split the cluster (" << ptrQ[iq1]->PDGName() << ' ' << ptrQ[iq2]->PDGName() << ") into a hadron and diquark.\n" << Exception::runerror; } // Check that new clusters can produce particles and there is enough // phase space to choose the drawn flavour m1 = ptrQ[iq1]->data().constituentMass(); m2 = ptrQ[iq2]->data().constituentMass(); m = newPtr1->data().constituentMass(); // Do not split in the case there is no phase space available if(mmax < m1+m + m2+m) continue; // power for splitting double exp1(_pSplitLight),exp2(_pSplitLight); if (CheckId::isExotic (ptrQ[iq1]->dataPtr())) exp1 = _pSplitExotic; else if(CheckId::hasBottom(ptrQ[iq1]->dataPtr())) exp1 = _pSplitBottom; else if(CheckId::hasCharm (ptrQ[iq1]->dataPtr())) exp1 = _pSplitCharm; if (CheckId::isExotic (ptrQ[iq2]->dataPtr())) exp2 = _pSplitExotic; else if(CheckId::hasBottom(ptrQ[iq2]->dataPtr())) exp2 = _pSplitBottom; else if(CheckId::hasCharm (ptrQ[iq2]->dataPtr())) exp2 = _pSplitCharm; // If, during the drawing of candidate masses, too many attempts fail // (because the phase space available is tiny) /// \todo run separate loop here? Mc1 = drawChildMass(mmax,m1,m2,m,exp1,soft1); Mc2 = drawChildMass(mmax,m2,m1,m,exp2,soft2); if(Mc1 < m1+m || Mc2 < m+m2 || Mc1+Mc2 > mmax) continue; // check if need to force meson clster to hadron toHadron = _hadronsSelector->chooseSingleHadron(ptrQ[iq1]->dataPtr(), newPtr1->dataPtr(),Mc1); if(toHadron) Mc1 = toHadron->mass(); // check if need to force diquark cluster to be on-shell toDiQuark = false; diquark = CheckId::makeDiquark(ptrQ[iq2]->dataPtr(), newPtr2->dataPtr()); if(Mc2 < diquark->constituentMass()) { Mc2 = diquark->constituentMass(); toDiQuark = true; } // if a beam cluster not allowed to decay to hadrons if(cluster->isBeamCluster() && toHadron && softUEisOn) continue; // Check if the decay kinematics is still possible: if not then // force the one-hadron decay for the other cluster as well. if(Mc1 + Mc2 > mmax) { if(!toHadron) { toHadron = _hadronsSelector->chooseSingleHadron(ptrQ[iq1]->dataPtr(), newPtr1->dataPtr(),mmax-Mc2); if(toHadron) Mc1 = toHadron->mass(); } else if(!toDiQuark) { Mc2 = _hadronsSelector->massLightestHadron(ptrQ[iq2]->dataPtr(), newPtr2->dataPtr()); toDiQuark = true; } } succeeded = (mmax >= Mc1+Mc2); } while (!succeeded && counter < max_loop); // check no of tries if(counter >= max_loop) return cutType(); // Determine the (5-components) momenta (all in the LAB frame) Lorentz5Momentum p0Q1 = ptrQ[iq1]->momentum(); // to be determined Lorentz5Momentum pClu1(Mc1), pClu2(Mc2), pQ1(m1), pQone(m), pQtwo(m), pQ2(m2); calculateKinematics(pDiQuark,p0Q1,toHadron,toDiQuark, pClu1,pClu2,pQ1,pQone,pQtwo,pQ2); // positions of the new clusters LorentzPoint pos1,pos2; Lorentz5Momentum pBaryon = pClu2+ptrQ[iother]->momentum(); calculatePositions(cluster->momentum(), cluster->vertex(), pClu1, pBaryon, pos1, pos2); // first the mesonic cluster/meson cutType rval; if(toHadron) { rval.first = produceHadron(toHadron, newPtr1, pClu1, pos1); finalhadrons.push_back(rval.first.first); } else { rval.first = produceCluster(ptrQ[iq1], newPtr1, pClu1, pos1, pQ1, pQone, rem1); } if(toDiQuark) { rem2 |= cluster->isBeamRemnant(iother); PPtr newDiQuark = diquark->produceParticle(pClu2); rval.second = produceCluster(newDiQuark, ptrQ[iother], pBaryon, pos2, pClu2, ptrQ[iother]->momentum(), rem2); } else { rval.second = produceCluster(ptrQ[iq2], newPtr2, pBaryon, pos2, pQ2, pQtwo, rem2, ptrQ[iother],cluster->isBeamRemnant(iother)); } cluster->isAvailable(false); return rval; } -ClusterFissioner::PPair +ClusterFissioner::PPair ClusterFissioner::produceHadron(tcPDPtr hadron, tPPtr newPtr, const Lorentz5Momentum &a, const LorentzPoint &b) const { PPair rval; if(hadron->coloured()) { rval.first = (_hadronsSelector->lightestHadron(hadron,newPtr->dataPtr()))->produceParticle(); } else rval.first = hadron->produceParticle(); rval.second = newPtr; rval.first->set5Momentum(a); rval.first->setVertex(b); return rval; } ClusterFissioner::PPair ClusterFissioner::produceCluster(tPPtr ptrQ, tPPtr newPtr, const Lorentz5Momentum & a, const LorentzPoint & b, const Lorentz5Momentum & c, const Lorentz5Momentum & d, bool isRem, tPPtr spect, bool remSpect) const { PPair rval; rval.second = newPtr; ClusterPtr cluster = !spect ? new_ptr(Cluster(ptrQ,rval.second)) : new_ptr(Cluster(ptrQ,rval.second,spect)); rval.first = cluster; cluster->set5Momentum(a); cluster->setVertex(b); assert(cluster->particle(0)->id() == ptrQ->id()); cluster->particle(0)->set5Momentum(c); cluster->particle(1)->set5Momentum(d); cluster->setBeamRemnant(0,isRem); if(remSpect) cluster->setBeamRemnant(2,remSpect); return rval; } void ClusterFissioner::drawNewFlavour(PPtr& newPtrPos,PPtr& newPtrNeg) const { // Flavour is assumed to be only u, d, s, with weights // (which are not normalized probabilities) given // by the same weights as used in HadronsSelector for - // the decay of clusters into two hadrons. - double prob_d = _hadronsSelector->pwtDquark(); - double prob_u = _hadronsSelector->pwtUquark(); - double prob_s = _hadronsSelector->pwtSquark(); + // the decay of clusters into two hadrons. + + + double prob_d; + double prob_u; + double prob_s; +switch(_fissionCluster){ + case 0: + prob_d = _hadronsSelector->pwtDquark(); + prob_u = _hadronsSelector->pwtUquark(); + prob_s = _hadronsSelector->pwtSquark(); + break; + case 1: + prob_d = _fissionPwtDquark; + prob_u = _fissionPwtUquark; + prob_s = _fissionPwtSquark; + break; +} int choice = UseRandom::rnd3(prob_u, prob_d, prob_s); long idNew = 0; switch (choice) { - case 0: idNew = ThePEG::ParticleID::u; break; + case 0: idNew = ThePEG::ParticleID::u; break; case 1: idNew = ThePEG::ParticleID::d; break; case 2: idNew = ThePEG::ParticleID::s; break; } newPtrPos = getParticle(idNew); newPtrNeg = getParticle(-idNew); assert (newPtrPos); assert(newPtrNeg); assert (newPtrPos->dataPtr()); assert(newPtrNeg->dataPtr()); } +void ClusterFissioner::drawNewFlavourEnhanced(PPtr& newPtrPos,PPtr& newPtrNeg, + Energy2 mass2) const { + + // Flavour is assumed to be only u, d, s, with weights + // (which are not normalized probabilities) given + // by the same weights as used in HadronsSelector for + // the decay of clusters into two hadrons. + + double prob_d; + double prob_u; + double prob_s = 0.; + double scale = abs(double(sqr(_m0Fission)/mass2)); +switch(_fissionCluster){ + case 0: + prob_d = _hadronsSelector->pwtDquark(); + prob_u = _hadronsSelector->pwtUquark(); + if (_enhanceSProb == 1) + prob_s = (20. < scale || scale < 0.) ? 0. : pow(_hadronsSelector->pwtSquark(),scale); + else if (_enhanceSProb == 2) + prob_s = (20. < scale || scale < 0.) ? 0. : exp(-scale); + break; + case 1: + prob_d = _fissionPwtDquark; + prob_u = _fissionPwtUquark; + if (_enhanceSProb == 1) + prob_s = (20. < scale || scale < 0.) ? 0. : pow(_fissionPwtSquark,scale); + else if (_enhanceSProb == 2) + prob_s = (20. < scale || scale < 0.) ? 0. : exp(-scale); + break; +} + + int choice = UseRandom::rnd3(prob_u, prob_d, prob_s); + long idNew = 0; + switch (choice) { + case 0: idNew = ThePEG::ParticleID::u; break; + case 1: idNew = ThePEG::ParticleID::d; break; + case 2: idNew = ThePEG::ParticleID::s; break; + } + newPtrPos = getParticle(idNew); + newPtrNeg = getParticle(-idNew); + assert (newPtrPos); + assert(newPtrNeg); + assert (newPtrPos->dataPtr()); + assert(newPtrNeg->dataPtr()); + +} + + +Energy2 ClusterFissioner::clustermass(const ClusterPtr & cluster){ + Lorentz5Momentum pIn = cluster->momentum(); + Energy2 endpointmass2 = sqr(cluster->particle(0)->mass() + + cluster->particle(1)->mass()); + Energy2 singletm2 = pIn.m2(); + return (_massMeasure == 0) ? singletm2 : singletm2 - endpointmass2; +} + + Energy ClusterFissioner::drawChildMass(const Energy M, const Energy m1, const Energy m2, const Energy m, const double expt, const bool soft) const { /*************************** * This method, given in input the cluster mass Mclu of an heavy cluster C, * made of consituents of masses m1 and m2, draws the masses Mclu1 and Mclu2 - * of, respectively, the children cluster C1, made of constituent masses m1 - * and m, and cluster C2, of mass Mclu2 and made of constituent masses m2 - * and m. The mass is extracted from one of the two following mass + * of, respectively, the children cluster C1, made of constituent masses m1 + * and m, and cluster C2, of mass Mclu2 and made of constituent masses m2 + * and m. The mass is extracted from one of the two following mass * distributions: * --- power-like ("normal" distribution) * d(Prob) / d(M^exponent) = const * where the exponent can be different from the two children C1 (exp1) * and C2 (exponent2). - * --- exponential ("soft" distribution) - * d(Prob) / d(M^2) = exp(-b*M) + * --- exponential ("soft" distribution) + * d(Prob) / d(M^2) = exp(-b*M) * where b = 2.0 / average. * Such distributions are limited below by the masses of - * the constituents quarks, and above from the mass of decaying cluster C. + * the constituents quarks, and above from the mass of decaying cluster C. * The choice of which of the two mass distributions to use for each of the * two cluster children is dictated by iRemnant (see below). - * If the number of attempts to extract a pair of mass values that are + * If the number of attempts to extract a pair of mass values that are * kinematically acceptable is above some fixed number (max_loop, see below) * the method gives up and returns false; otherwise, when it succeeds, it - * returns true. + * returns true. * - * These distributions have been modified from HERWIG: + * These distributions have been modified from HERWIG: * Before these were: * Mclu1 = m1 + (Mclu - m1 - m2)*pow( rnd(), 1.0/exponent1 ); - * The new one coded here is a more efficient version, same density - * but taking into account 'in phase space from' beforehand + * The new one coded here is a more efficient version, same density + * but taking into account 'in phase space from' beforehand ***************************/ // hard cluster if(!soft) { - return pow(UseRandom::rnd(pow((M-m1-m2-m)*UnitRemoval::InvE, expt), + return pow(UseRandom::rnd(pow((M-m1-m2-m)*UnitRemoval::InvE, expt), pow(m*UnitRemoval::InvE, expt)), 1./expt )*UnitRemoval::E + m1; } // Otherwise it uses a soft mass distribution - else { + else { static const InvEnergy b = 2.0 / _btClM; Energy max = M-m1-m2-2.0*m; double rmin = b*max; rmin = ( rmin < 50 ) ? exp(-rmin) : 0.; double r1; do { r1 = UseRandom::rnd(rmin, 1.0) * UseRandom::rnd(rmin, 1.0); } while (r1 < rmin); return m1 + m - log(r1)/b; } } void ClusterFissioner::calculateKinematics(const Lorentz5Momentum & pClu, const Lorentz5Momentum & p0Q1, const bool toHadron1, const bool toHadron2, - Lorentz5Momentum & pClu1, - Lorentz5Momentum & pClu2, + Lorentz5Momentum & pClu1, + Lorentz5Momentum & pClu2, Lorentz5Momentum & pQ1, Lorentz5Momentum & pQbar, - Lorentz5Momentum & pQ, + Lorentz5Momentum & pQ, Lorentz5Momentum & pQ2bar) const { /****************** * This method solves the kinematics of the two body cluster decay: * C (Q1 Q2bar) ---> C1 (Q1 Qbar) + C2 (Q Q2bar) - * In input we receive the momentum of C, pClu, and the momentum + * In input we receive the momentum of C, pClu, and the momentum * of the quark Q1 (constituent of C), p0Q1, both in the LAB frame. * Furthermore, two boolean variables inform whether the two fission * products (C1, C2) decay immediately into a single hadron (in which - * case the cluster itself is identify with that hadron) and we do + * case the cluster itself is identify with that hadron) and we do * not have to solve the kinematics of the components (Q1,Qbar) for * C1 and (Q,Q2bar) for C2. - * The output is given by the following momenta (all 5-components, + * The output is given by the following momenta (all 5-components, * and all in the LAB frame): - * pClu1 , pClu2 respectively of C1 , C2 + * pClu1 , pClu2 respectively of C1 , C2 * pQ1 , pQbar respectively of Q1 , Qbar in C1 * pQ , pQ2bar respectively of Q , Q2 in C2 * The assumption, suggested from the string model, is that, in C frame, * C1 and its constituents Q1 and Qbar are collinear, and collinear to * the direction of Q1 in C (that is before cluster decay); similarly, * (always in the C frame) C2 and its constituents Q and Q2bar are * collinear (and therefore anti-collinear with C1,Q1,Qbar). * The solution is then obtained by using Lorentz boosts, as follows. * The kinematics of C1 and C2 is solved in their parent C frame, * and then boosted back in the LAB. The kinematics of Q1 and Qbar - * is solved in their parent C1 frame and then boosted back in the LAB; - * similarly, the kinematics of Q and Q2bar is solved in their parent + * is solved in their parent C1 frame and then boosted back in the LAB; + * similarly, the kinematics of Q and Q2bar is solved in their parent * C2 frame and then boosted back in the LAB. In each of the three - * "two-body decay"-like cases, we use the fact that the direction - * of the motion of the decay products is known in the rest frame of - * their parent. This is obvious for the first case in which the + * "two-body decay"-like cases, we use the fact that the direction + * of the motion of the decay products is known in the rest frame of + * their parent. This is obvious for the first case in which the * parent rest frame is C; but it is also true in the other two cases - * where the rest frames are C1 and C2. This is because C1 and C2 - * are boosted w.r.t. C in the same direction where their components, + * where the rest frames are C1 and C2. This is because C1 and C2 + * are boosted w.r.t. C in the same direction where their components, * respectively (Q1,Qbar) and (Q,Q2bar) move in C1 and C2 rest frame - * respectively. + * respectively. * Of course, although the notation used assumed that C = (Q1 Q2bar) * where Q1 is a quark and Q2bar an antiquark, indeed everything remain * unchanged also in all following cases: * Q1 quark, Q2bar antiquark; --> Q quark; - * Q1 antiquark , Q2bar quark; --> Q antiquark; + * Q1 antiquark , Q2bar quark; --> Q antiquark; * Q1 quark, Q2bar diquark; --> Q quark * Q1 antiquark, Q2bar anti-diquark; --> Q antiquark * Q1 diquark, Q2bar quark --> Q antiquark * Q1 anti-diquark, Q2bar antiquark; --> Q quark **************************/ // Calculate the unit three-vector, in the C frame, along which // all of the constituents and children clusters move. Lorentz5Momentum u(p0Q1); u.boost( -pClu.boostVector() ); // boost from LAB to C // the unit three-vector is then u.vect().unit() - - // Calculate the momenta of C1 and C2 in the (parent) C frame first, - // where the direction of C1 is u.vect().unit(), and then boost back in the + + // Calculate the momenta of C1 and C2 in the (parent) C frame first, + // where the direction of C1 is u.vect().unit(), and then boost back in the // LAB frame. if (pClu.m() < pClu1.mass() + pClu2.mass() ) { - throw Exception() << "Impossible Kinematics in ClusterFissioner::calculateKinematics() (A)" + throw Exception() << "Impossible Kinematics in ClusterFissioner::calculateKinematics() (A)" << Exception::eventerror; } - Kinematics::twoBodyDecay(pClu, pClu1.mass(), pClu2.mass(), + Kinematics::twoBodyDecay(pClu, pClu1.mass(), pClu2.mass(), u.vect().unit(), pClu1, pClu2); // In the case that cluster1 does not decay immediately into a single hadron, // calculate the momenta of Q1 (as constituent of C1) and Qbar in the - // (parent) C1 frame first, where the direction of Q1 is u.vect().unit(), - // and then boost back in the LAB frame. + // (parent) C1 frame first, where the direction of Q1 is u.vect().unit(), + // and then boost back in the LAB frame. if(!toHadron1) { if (pClu1.m() < pQ1.mass() + pQbar.mass() ) { - throw Exception() << "Impossible Kinematics in ClusterFissioner::calculateKinematics() (B)" + throw Exception() << "Impossible Kinematics in ClusterFissioner::calculateKinematics() (B)" << Exception::eventerror; } - Kinematics::twoBodyDecay(pClu1, pQ1.mass(), pQbar.mass(), + Kinematics::twoBodyDecay(pClu1, pQ1.mass(), pQbar.mass(), u.vect().unit(), pQ1, pQbar); } // In the case that cluster2 does not decay immediately into a single hadron, // Calculate the momenta of Q and Q2bar (as constituent of C2) in the - // (parent) C2 frame first, where the direction of Q is u.vect().unit(), - // and then boost back in the LAB frame. + // (parent) C2 frame first, where the direction of Q is u.vect().unit(), + // and then boost back in the LAB frame. if(!toHadron2) { if (pClu2.m() < pQ.mass() + pQ2bar.mass() ) { - throw Exception() << "Impossible Kinematics in ClusterFissioner::calculateKinematics() (C)" + throw Exception() << "Impossible Kinematics in ClusterFissioner::calculateKinematics() (C)" << Exception::eventerror; } - Kinematics::twoBodyDecay(pClu2, pQ.mass(), pQ2bar.mass(), + Kinematics::twoBodyDecay(pClu2, pQ.mass(), pQ2bar.mass(), u.vect().unit(), pQ, pQ2bar); } } void ClusterFissioner::calculatePositions(const Lorentz5Momentum & pClu, const LorentzPoint & positionClu, const Lorentz5Momentum & pClu1, const Lorentz5Momentum & pClu2, LorentzPoint & positionClu1, LorentzPoint & positionClu2) const { // Determine positions of cluster children. // See Marc Smith's thesis, page 127, formulas (4.122) and (4.123). Energy Mclu = pClu.m(); Energy Mclu1 = pClu1.m(); Energy Mclu2 = pClu2.m(); - // Calculate the unit three-vector, in the C frame, along which + // Calculate the unit three-vector, in the C frame, along which // children clusters move. Lorentz5Momentum u(pClu1); u.boost( -pClu.boostVector() ); // boost from LAB to C frame // the unit three-vector is then u.vect().unit() Energy pstarChild = Kinematics::pstarTwoBodyDecay(Mclu,Mclu1,Mclu2); // First, determine the relative positions of the children clusters // in the parent cluster reference frame. Length x1 = ( 0.25*Mclu + 0.5*( pstarChild + (sqr(Mclu2) - sqr(Mclu1))/(2.0*Mclu)))/_kappa; - Length t1 = Mclu/_kappa - x1; + Length t1 = Mclu/_kappa - x1; LorentzDistance distanceClu1( x1 * u.vect().unit(), t1 ); Length x2 = (-0.25*Mclu + 0.5*(-pstarChild + (sqr(Mclu2) - sqr(Mclu1))/(2.0*Mclu)))/_kappa; Length t2 = Mclu/_kappa + x2; LorentzDistance distanceClu2( x2 * u.vect().unit(), t2 ); // Then, transform such relative positions from the parent cluster // reference frame to the Lab frame. distanceClu1.boost( pClu.boostVector() ); distanceClu2.boost( pClu.boostVector() ); // Finally, determine the absolute positions in the Lab frame. positionClu1 = positionClu + distanceClu1; positionClu2 = positionClu + distanceClu2; } bool ClusterFissioner::isHeavy(tcClusterPtr clu) { // default double clpow = _clPowLight; Energy clmax = _clMaxLight; // particle data for constituents tcPDPtr cptr[3]={tcPDPtr(),tcPDPtr(),tcPDPtr()}; for(int ix=0;ixnumComponents(),3);++ix) { cptr[ix]=clu->particle(ix)->dataPtr(); } // different parameters for exotic, bottom and charm clusters if(CheckId::isExotic(cptr[0],cptr[1],cptr[1])) { clpow = _clPowExotic; clmax = _clMaxExotic; } else if(CheckId::hasBottom(cptr[0],cptr[1],cptr[1])) { clpow = _clPowBottom; clmax = _clMaxBottom; } else if(CheckId::hasCharm(cptr[0],cptr[1],cptr[1])) { clpow = _clPowCharm; clmax = _clMaxCharm; } bool aboveCutoff = ( - pow(clu->mass()*UnitRemoval::InvE , clpow) - > - pow(clmax*UnitRemoval::InvE, clpow) + pow(clu->mass()*UnitRemoval::InvE , clpow) + > + pow(clmax*UnitRemoval::InvE, clpow) + pow(clu->sumConstituentMasses()*UnitRemoval::InvE, clpow) ); // required test for SUSY clusters, since aboveCutoff alone // cannot guarantee (Mc > m1 + m2 + 2*m) in cut() - static const Energy minmass + static const Energy minmass = getParticleData(ParticleID::d)->constituentMass(); - bool canSplitMinimally + bool canSplitMinimally = clu->mass() > clu->sumConstituentMasses() + 2.0 * minmass; return aboveCutoff && canSplitMinimally; } diff --git a/Hadronization/ClusterFissioner.h b/Hadronization/ClusterFissioner.h --- a/Hadronization/ClusterFissioner.h +++ b/Hadronization/ClusterFissioner.h @@ -1,386 +1,412 @@ // -*- C++ -*- // // ClusterFissioner.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_ClusterFissioner_H #define HERWIG_ClusterFissioner_H #include #include "CluHadConfig.h" #include "HadronSelector.h" #include "ClusterFissioner.fh" namespace Herwig { using namespace ThePEG; //class Cluster; // forward declaration /** \ingroup Hadronization * \class ClusterFissioner * \brief This class handles clusters which are too heavy. * \author Philip Stephens * \author Alberto Ribon * \author Stefan Gieseke * - * This class does the job of chopping up either heavy clusters or beam - * clusters in two lighter ones. The procedure is repeated recursively until + * This class does the job of chopping up either heavy clusters or beam + * clusters in two lighter ones. The procedure is repeated recursively until * all of the cluster children have masses below some threshold values. * * For the beam remnant clusters, at the moment what is done is the following. - * In the case that the soft underlying event is switched on, the + * In the case that the soft underlying event is switched on, the * beam remnant clusters are tagged as not available, - * therefore they will not be treated at all during the hadronization. + * therefore they will not be treated at all during the hadronization. * In the case instead that the soft underlying event is switched off, * then the beam remnant clusters are treated exactly as "normal" clusters, * with the only exception of the mass spectrum used to generate the * cluster children masses. For non-beam clusters, the masses of the cluster * children are draw from a power-like mass distribution; for beam clusters, - * according to the value of the flag _IOpRem, either both - * children masses are draw from a fast-decreasing exponential mass - * distribution (case _IOpRem == 0, or, indendently by - * _IOpRem, in the special case that the beam cluster contains two + * according to the value of the flag _IOpRem, either both + * children masses are draw from a fast-decreasing exponential mass + * distribution (case _IOpRem == 0, or, indendently by + * _IOpRem, in the special case that the beam cluster contains two * beam remnants), or one mass from the exponential distribution (corresponding - * of the cluster child with the beam remnant) and the other with the usual - * power-like distribution (case _IOpRem == 1, which is the - * default one, as in Herwig 6.3). + * of the cluster child with the beam remnant) and the other with the usual + * power-like distribution (case _IOpRem == 1, which is the + * default one, as in Herwig 6.3). * - * The reason behind the use of a fast-decreasing exponential distribution + * The reason behind the use of a fast-decreasing exponential distribution * is that to avoid a large transverse energy from the many sequential - * fissions that would otherwise occur due to the typical large cluster - * mass of beam clusters. Using instead an exponential distribution - * the masses of the two cluster children will be very small (order of + * fissions that would otherwise occur due to the typical large cluster + * mass of beam clusters. Using instead an exponential distribution + * the masses of the two cluster children will be very small (order of * GeV). * * The rationale behind the implementation of the splitting of clusters - * has been to preserve *all* of the information about such splitting + * has been to preserve *all* of the information about such splitting * process. More explicitly a ThePEG::Step class is passed in and the * new clusters are added to the step as the decay products of the - * heavy cluster. This approach has the twofold - * advantage to provide all of the information that could be needed - * (expecially in future developments), without any information loss, - * and furthermore it allows a better debugging. + * heavy cluster. This approach has the twofold + * advantage to provide all of the information that could be needed + * (expecially in future developments), without any information loss, + * and furthermore it allows a better debugging. * - * @see HadronSelector + * @see HadronSelector * @see \ref ClusterFissionerInterfaces "The interfaces" * defined for ClusterFissioner. - */ + */ class ClusterFissioner: public Interfaced { public: /** @name Standard constructors and destructors. */ //@{ /** * Default constructor. */ ClusterFissioner(); + //@} /** Splits the clusters which are too heavy. * - * Split either heavy clusters or beam clusters recursively until all + * Split either heavy clusters or beam clusters recursively until all * children have mass below some threshold. Heavy clusters are those that - * satisfy the condition + * satisfy the condition * \f[ M^P > C^P + S^P \f] * where \f$ M \f$ is the clusters mass, \f$ P \f$ is the parameter - * ClPow, \f$ C \f$ is the parameter ClMax and \f$ S \f$ is the + * ClPow, \f$ C \f$ is the parameter ClMax and \f$ S \f$ is the * sum of the clusters constituent partons. * For beam clusters, they are split only if the soft underlying event * is switched off, otherwise these clusters will be tagged as unavailable - * and they will not be treated by the hadronization altogether. + * and they will not be treated by the hadronization altogether. * In the case beam clusters will be split, the procedure is exactly * the same as for normal non-beam clusters, with the only exception - * of the mass spectrum from which to draw the masses of the two + * of the mass spectrum from which to draw the masses of the two * cluster children (see method drawChildrenMasses for details). */ tPVector fission(ClusterVector & clusters, bool softUEisOn); 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. */ ClusterFissioner & operator=(const ClusterFissioner &) = delete; - /** + /** * This method directs the splitting of the heavy clusters * - * This method does the splitting of the clusters and all of its cluster + * This method does the splitting of the clusters and all of its cluster * children, if heavy. All of these new children clusters are added to the * collection of clusters. The method works as follows. * Initially the vector contains just the stack of input pointers to the * clusters to be split. Then it will be filled recursively by all * of the cluster's children that are heavy enough to require - * to be split. In each loop, the last element of the vector is + * to be split. In each loop, the last element of the vector is * considered (only once because it is then removed from the vector). * * \todo is the following still true? * For normal, non-beam clusters, a power-like mass distribution - * is used, whereas for beam clusters a fast-decreasing exponential mass - * distribution is used instead. This avoids many iterative splitting which - * could produce an unphysical large transverse energy from a supposed + * is used, whereas for beam clusters a fast-decreasing exponential mass + * distribution is used instead. This avoids many iterative splitting which + * could produce an unphysical large transverse energy from a supposed * soft beam remnant process. */ void cut(stack &, ClusterVector&, tPVector & finalhadrons, bool softUEisOn); public: /** * Definition for easy passing of two particles. */ typedef pair PPair; /** * Definition for use in the cut function. */ typedef pair cutType; - /** + /** * Splits the input cluster. * * Split the input cluster (which can be either an heavy non-beam * cluster or a beam cluster). The result is two pairs of particles. The * first element of each pair is new cluster/hadron, while the second * element of each pair is the particle drawn from the vacuum to create * the new cluster/hadron. * Notice that this method treats also beam clusters by using a different * mass spectrum used to generate the cluster child masses (see method * drawChildMass). */ //@{ /** * Split two-component cluster */ virtual cutType cutTwo(ClusterPtr &, tPVector & finalhadrons, bool softUEisOn); /** * Split three-component cluster */ virtual cutType cutThree(ClusterPtr &, tPVector & finalhadrons, bool softUEisOn); //@} public: /** * Produces a hadron and returns the flavour drawn from the vacuum. * * This routine produces a new hadron. It * also sets the momentum and vertex to the values given. */ PPair produceHadron(tcPDPtr hadron, tPPtr newPtr, const Lorentz5Momentum &a, const LorentzPoint &b) const; protected: /** * Produces a cluster from the flavours passed in. * * This routine produces a new cluster with the flavours given by ptrQ and newPtr. * The new 5 momentum is a and the parent momentum are c and d. C is for the * ptrQ and d is for the new particle newPtr. rem specifies whether the existing * particle is a beam remnant or not. */ PPair produceCluster(tPPtr ptrQ, tPPtr newPtr, const Lorentz5Momentum &a, const LorentzPoint &b, const Lorentz5Momentum &c, const Lorentz5Momentum &d, const bool rem, tPPtr spect=tPPtr(), bool remSpect=false) const; /** * Returns the new quark-antiquark pair * needed for fission of a heavy cluster. Equal probabilities - * are assumed for producing u, d, or s pairs. + * are assumed for producing u, d, or s pairs. */ void drawNewFlavour(PPtr& newPtrPos,PPtr& newPtrNeg) const; - + + /** * Produces the mass of a child cluster. * - * Draw the masses \f$M'\f$ of the the cluster child produced + * Draw the masses \f$M'\f$ of the the cluster child produced * by the fission of an heavy cluster (of mass M). m1, m2 are the masses - * of the constituents of the cluster; m is the mass of the quark extract + * of the constituents of the cluster; m is the mass of the quark extract * from the vacuum (together with its antiparticle). The algorithm produces * the mass of the cluster formed with consituent m1. * Two mass distributions can be used for the child cluster mass: - * -# power-like mass distribution ("normal" mass) with power exp + * -# power-like mass distribution ("normal" mass) with power exp * \f[ M' = {\rm rnd}((M-m_1-m_2-m)^P, m^p)^{1/P} + m_1 \f] * where \f$ P \f$ is a parameter of the model and \f$ \rm{rnd} \f$ is * the function: * \f[ \rm{rnd}(a,b) = (1-r)a + r b \f] * and here \f$ r \f$ is a random number [0,1]. - * -# fast-decreasing exponential mass distribution ("soft" mass) with - * rmin. rmin is given by + * -# fast-decreasing exponential mass distribution ("soft" mass) with + * rmin. rmin is given by * \f[ r_{\rm min} = \exp(-b (M - m_1 - m_2 - 2 m)) \f] * where \f$ b \f$ is a parameter of the model. The generated mass is * given by * \f[ M' = m_1 + m - \frac{\log\left( * {\rm rnd}(r_{\rm min}, 1-r_{\rm min})\right)}{b} \f]. * * The choice of which mass distribution should be used for each of the two * cluster children is dictated by the parameter soft. */ - Energy drawChildMass(const Energy M, const Energy m1, const Energy m2, + Energy drawChildMass(const Energy M, const Energy m1, const Energy m2, const Energy m, const double exp, const bool soft) const; /** * Determines the kinematics of a heavy cluster decay C->C1 + C2 */ - void calculateKinematics(const Lorentz5Momentum &pClu, - const Lorentz5Momentum &p0Q1, + void calculateKinematics(const Lorentz5Momentum &pClu, + const Lorentz5Momentum &p0Q1, const bool toHadron1, const bool toHadron2, - Lorentz5Momentum &pClu1, Lorentz5Momentum &pClu2, - Lorentz5Momentum &pQ1, Lorentz5Momentum &pQb, + Lorentz5Momentum &pClu1, Lorentz5Momentum &pClu2, + Lorentz5Momentum &pQ1, Lorentz5Momentum &pQb, Lorentz5Momentum &pQ2, Lorentz5Momentum &pQ2b) const; /** * Determine the positions of the two children clusters. * * This routine generates the momentum of the decay products. It also * generates the momentum in the lab frame of the partons drawn out of * the vacuum. */ - void calculatePositions(const Lorentz5Momentum &pClu, + void calculatePositions(const Lorentz5Momentum &pClu, const LorentzPoint & positionClu, - const Lorentz5Momentum & pClu1, - const Lorentz5Momentum & pClu2, - LorentzPoint & positionClu1, + const Lorentz5Momentum & pClu1, + const Lorentz5Momentum & pClu2, + LorentzPoint & positionClu1, LorentzPoint & positionClu2 ) const; protected: /** @name Access members for child classes. */ //@{ /** * Access to the hadron selector */ HadronSelectorPtr hadronsSelector() const {return _hadronsSelector;} /** * Access to soft-cluster parameter */ Energy btClM() const {return _btClM;} /** * Cluster splitting paramater for light quarks */ - double pSplitLight() const {return _pSplitLight;} + double pSplitLight() const {return _pSplitLight;} /** * Cluster splitting paramater for bottom quarks */ double pSplitBottom() const {return _pSplitBottom;} /** * Cluster splitting paramater for charm quarks */ double pSplitCharm() const {return _pSplitCharm;} /** * Cluster splitting paramater for exotic particles */ double pSplitExotic() const {return _pSplitExotic;} //@} - + private: /** * Check if a cluster is heavy enough to split again */ bool isHeavy(tcClusterPtr ); /** * A pointer to a Herwig::HadronSelector object for generating hadrons. */ HadronSelectorPtr _hadronsSelector; /** - * @name The Cluster max mass,dependant on which quarks are involved, used to determine when + * @name The Cluster max mass,dependant on which quarks are involved, used to determine when * fission will occur. */ //@{ Energy _clMaxLight; Energy _clMaxBottom; Energy _clMaxCharm; Energy _clMaxExotic; //@} /** * @name The power used to determine when cluster fission will occur. */ //@{ double _clPowLight; double _clPowBottom; double _clPowCharm; double _clPowExotic; //@} /** * @name The power, dependant on whic quarks are involved, used in the cluster mass generation. */ //@{ double _pSplitLight; double _pSplitBottom; double _pSplitCharm; double _pSplitExotic; + + + // weights for alternaive cluster fission + double _fissionPwtUquark; + double _fissionPwtDquark; + double _fissionPwtSquark; + + /** + * Flag used to determine between normal cluster fission and alternative cluster fission + */ + int _fissionCluster; + //@} /** * Parameter used (2/b) for the beam cluster mass generation. * Currently hard coded value. */ - Energy _btClM; + Energy _btClM; /** * Flag used to determine what distributions to use for the cluster masses. */ int _iopRem; /** * The string constant */ Tension _kappa; + int _enhanceSProb; + + Energy _m0Fission; + + Energy2 clustermass(const ClusterPtr & cluster); + + int _massMeasure; + +protected: + void drawNewFlavourEnhanced(PPtr& newPtrPos,PPtr& newPtrNeg, Energy2 mass2) const; + + }; } #endif /* HERWIG_ClusterFissioner_H */ diff --git a/Hadronization/ClusterHadronizationHandler.cc b/Hadronization/ClusterHadronizationHandler.cc --- a/Hadronization/ClusterHadronizationHandler.cc +++ b/Hadronization/ClusterHadronizationHandler.cc @@ -1,307 +1,307 @@ // -*- C++ -*- // // ClusterHadronizationHandler.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the ClusterHadronizationHandler class. // #include "ClusterHadronizationHandler.h" #include #include #include -#include -#include +#include +#include #include #include #include #include #include #include #include #include #include #include "Herwig/Utilities/EnumParticles.h" #include "CluHadConfig.h" -#include "Cluster.h" +#include "Cluster.h" #include using namespace Herwig; ClusterHadronizationHandler * ClusterHadronizationHandler::currentHandler_ = 0; DescribeClass describeClusterHadronizationHandler("Herwig::ClusterHadronizationHandler",""); IBPtr ClusterHadronizationHandler::clone() const { return new_ptr(*this); } IBPtr ClusterHadronizationHandler::fullclone() const { return new_ptr(*this); } -void ClusterHadronizationHandler::persistentOutput(PersistentOStream & os) +void ClusterHadronizationHandler::persistentOutput(PersistentOStream & os) const { os << _partonSplitter << _clusterFinder << _colourReconnector << _clusterFissioner << _lightClusterDecayer << _clusterDecayer << ounit(_minVirtuality2,GeV2) << ounit(_maxDisplacement,mm) << _underlyingEventHandler << _reduceToTwoComponents; } void ClusterHadronizationHandler::persistentInput(PersistentIStream & is, int) { is >> _partonSplitter >> _clusterFinder >> _colourReconnector >> _clusterFissioner >> _lightClusterDecayer >> _clusterDecayer >> iunit(_minVirtuality2,GeV2) >> iunit(_maxDisplacement,mm) >> _underlyingEventHandler >> _reduceToTwoComponents; } void ClusterHadronizationHandler::Init() { static ClassDocumentation documentation ("This is the main handler class for the Cluster Hadronization", "The hadronization was performed using the cluster model of \\cite{Webber:1983if}.", "%\\cite{Webber:1983if}\n" "\\bibitem{Webber:1983if}\n" " B.~R.~Webber,\n" " ``A QCD Model For Jet Fragmentation Including Soft Gluon Interference,''\n" " Nucl.\\ Phys.\\ B {\\bf 238}, 492 (1984).\n" " %%CITATION = NUPHA,B238,492;%%\n" // main manual ); - static Reference - interfacePartonSplitter("PartonSplitter", - "A reference to the PartonSplitter object", + static Reference + interfacePartonSplitter("PartonSplitter", + "A reference to the PartonSplitter object", &Herwig::ClusterHadronizationHandler::_partonSplitter, false, false, true, false); - static Reference - interfaceClusterFinder("ClusterFinder", - "A reference to the ClusterFinder object", + static Reference + interfaceClusterFinder("ClusterFinder", + "A reference to the ClusterFinder object", &Herwig::ClusterHadronizationHandler::_clusterFinder, false, false, true, false); - static Reference - interfaceColourReconnector("ColourReconnector", - "A reference to the ColourReconnector object", + static Reference + interfaceColourReconnector("ColourReconnector", + "A reference to the ColourReconnector object", &Herwig::ClusterHadronizationHandler::_colourReconnector, false, false, true, false); - static Reference - interfaceClusterFissioner("ClusterFissioner", - "A reference to the ClusterFissioner object", + static Reference + interfaceClusterFissioner("ClusterFissioner", + "A reference to the ClusterFissioner object", &Herwig::ClusterHadronizationHandler::_clusterFissioner, false, false, true, false); - static Reference - interfaceLightClusterDecayer("LightClusterDecayer", - "A reference to the LightClusterDecayer object", + static Reference + interfaceLightClusterDecayer("LightClusterDecayer", + "A reference to the LightClusterDecayer object", &Herwig::ClusterHadronizationHandler::_lightClusterDecayer, false, false, true, false); - static Reference - interfaceClusterDecayer("ClusterDecayer", - "A reference to the ClusterDecayer object", + static Reference + interfaceClusterDecayer("ClusterDecayer", + "A reference to the ClusterDecayer object", &Herwig::ClusterHadronizationHandler::_clusterDecayer, false, false, true, false); - static Parameter interfaceMinVirtuality2 + static Parameter interfaceMinVirtuality2 ("MinVirtuality2", "Minimum virtuality^2 of partons to use in calculating distances (unit [GeV2]).", &ClusterHadronizationHandler::_minVirtuality2, GeV2, 0.1*GeV2, ZERO, 10.0*GeV2,false,false,false); - - static Parameter interfaceMaxDisplacement + + static Parameter interfaceMaxDisplacement ("MaxDisplacement", "Maximum displacement that is allowed for a particle (unit [millimeter]).", - &ClusterHadronizationHandler::_maxDisplacement, mm, 1.0e-10*mm, + &ClusterHadronizationHandler::_maxDisplacement, mm, 1.0e-10*mm, 0.0*mm, 1.0e-9*mm,false,false,false); static Reference interfaceUnderlyingEventHandler ("UnderlyingEventHandler", "Pointer to the handler for the Underlying Event. " "Set to NULL to disable.", &ClusterHadronizationHandler::_underlyingEventHandler, false, false, true, true, false); static Switch interfaceReduceToTwoComponents ("ReduceToTwoComponents", "Whether or not to reduce three component baryon-number violating clusters to two components before cluster splitting or leave" " this till after the cluster splitting", &ClusterHadronizationHandler::_reduceToTwoComponents, true, false, false); static SwitchOption interfaceReduceToTwoComponentsYes (interfaceReduceToTwoComponents, "BeforeSplitting", "Reduce to two components", true); static SwitchOption interfaceReduceToTwoComponentsNo (interfaceReduceToTwoComponents, "AfterSplitting", "Treat as three components", false); } namespace { void extractChildren(tPPtr p, set & all) { if (p->children().empty()) return; for (PVector::const_iterator child = p->children().begin(); child != p->children().end(); ++child) { all.insert(*child); extractChildren(*child, all); } } } void ClusterHadronizationHandler:: handle(EventHandler & ch, const tPVector & tagged, const Hint &) { useMe(); currentHandler_ = this; PVector currentlist(tagged.begin(),tagged.end()); // set the scale for coloured particles to just above the gluon mass squared // if less than this so they are classed as perturbative Energy2 Q02 = 1.01*sqr(getParticleData(ParticleID::g)->constituentMass()); for(unsigned int ix=0;ixscale()scale(Q02); } // split the gluons _partonSplitter->split(currentlist); + // form the clusters ClusterVector clusters = _clusterFinder->formClusters(currentlist); // reduce BV clusters to two components now if needed if(_reduceToTwoComponents) _clusterFinder->reduceToTwoComponents(clusters); // perform colour reconnection if needed and then // decay the clusters into one hadron bool lightOK = false; short tried = 0; const ClusterVector savedclusters = clusters; tPVector finalHadrons; // only needed for partonic decayer while (!lightOK && tried++ < 10) { // no colour reconnection with baryon-number-violating (BV) clusters ClusterVector CRclusters, BVclusters; CRclusters.reserve( clusters.size() ); BVclusters.reserve( clusters.size() ); for (size_t ic = 0; ic < clusters.size(); ++ic) { ClusterPtr cl = clusters.at(ic); bool hasClusterParent = false; for (unsigned int ix=0; ix < cl->parents().size(); ++ix) { if (cl->parents()[ix]->id() == ParticleID::Cluster) { hasClusterParent = true; break; } } if (cl->numComponents() > 2 || hasClusterParent) BVclusters.push_back(cl); else CRclusters.push_back(cl); } // colour reconnection _colourReconnector->rearrange(CRclusters); // tag new clusters as children of the partons to hadronize _setChildren(CRclusters); - // forms diquarks + // forms diquarks _clusterFinder->reduceToTwoComponents(CRclusters); // recombine vectors of (possibly) reconnected and BV clusters clusters.clear(); clusters.insert( clusters.end(), CRclusters.begin(), CRclusters.end() ); clusters.insert( clusters.end(), BVclusters.begin(), BVclusters.end() ); - + // fission of heavy clusters // NB: during cluster fission, light hadrons might be produced straight away finalHadrons = _clusterFissioner->fission(clusters,isSoftUnderlyingEventON()); // if clusters not previously reduced to two components do it now if(!_reduceToTwoComponents) _clusterFinder->reduceToTwoComponents(clusters); lightOK = _lightClusterDecayer->decay(clusters,finalHadrons); // if the decay of the light clusters was not successful, undo the cluster // fission and decay steps and revert to the original state of the event // record if (!lightOK) { clusters = savedclusters; - for_each(clusters.begin(), - clusters.end(), + for_each(clusters.begin(), + clusters.end(), mem_fun(&Particle::undecay)); } - } + } if (!lightOK) { - throw Exception("CluHad::handle(): tried LightClusterDecayer 10 times!", + throw Exception("CluHad::handle(): tried LightClusterDecayer 10 times!", Exception::eventerror); } // decay the remaining clusters _clusterDecayer->decay(clusters,finalHadrons); // ***************************************** // ***************************************** // ***************************************** StepPtr pstep = newStep(); set allDecendants; for (tPVector::const_iterator it = tagged.begin(); it != tagged.end(); ++it) { extractChildren(*it, allDecendants); } for(set::const_iterator it = allDecendants.begin(); it != allDecendants.end(); ++it) { - // this is a workaround because the set sometimes + // this is a workaround because the set sometimes // re-orders parents after their children if ((*it)->children().empty()) pstep->addDecayProduct(*it); else { pstep->addDecayProduct(*it); pstep->addIntermediate(*it); } } // ***************************************** // ***************************************** // ***************************************** // soft underlying event if needed if (isSoftUnderlyingEventON()) { assert(_underlyingEventHandler); ch.performStep(_underlyingEventHandler,Hint::Default()); } } // Sets parent child relationship of all clusters with two components -// Relationships for clusters with more than two components are set elsewhere in the Colour Reconnector +// Relationships for clusters with more than two components are set elsewhere in the Colour Reconnector void ClusterHadronizationHandler::_setChildren(const ClusterVector & clusters) const { // erase existing information about the partons' children tPVector partons; for ( const auto & cl : clusters ) { if ( cl->numComponents() > 2 ) continue; partons.push_back( cl->colParticle() ); partons.push_back( cl->antiColParticle() ); } // erase all previous information about parent child relationship for_each(partons.begin(), partons.end(), mem_fun(&Particle::undecay)); // give new parents to the clusters: their constituents for ( const auto & cl : clusters ) { if ( cl->numComponents() > 2 ) continue; cl->colParticle()->addChild(cl); cl->antiColParticle()->addChild(cl); } } - diff --git a/Hadronization/ColourReconnector.cc b/Hadronization/ColourReconnector.cc --- a/Hadronization/ColourReconnector.cc +++ b/Hadronization/ColourReconnector.cc @@ -1,821 +1,822 @@ // -*- C++ -*- // // ColourReconnector.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the ColourReconnector class. // #include "ColourReconnector.h" #include "Cluster.h" #include #include #include #include #include #include #include #include "Herwig/Utilities/Maths.h" using namespace Herwig; using CluVecIt = ColourReconnector::CluVecIt; using Constants::pi; using Constants::twopi; DescribeClass describeColourReconnector("Herwig::ColourReconnector",""); IBPtr ColourReconnector::clone() const { return new_ptr(*this); } IBPtr ColourReconnector::fullclone() const { return new_ptr(*this); } void ColourReconnector::rearrange(ClusterVector & clusters) { if (_clreco == 0) return; // need at least two clusters if (clusters.size() < 2) return; // do the colour reconnection switch (_algorithm) { case 0: _doRecoPlain(clusters); break; case 1: _doRecoStatistical(clusters); break; case 2: _doRecoBaryonic(clusters); break; } } Energy2 ColourReconnector::_clusterMassSum(const PVector & q, const PVector & aq) const { const size_t nclusters = q.size(); assert (aq.size() == nclusters); Energy2 sum = ZERO; for (size_t i = 0; i < nclusters; i++) sum += ( q[i]->momentum() + aq[i]->momentum() ).m2(); return sum; } bool ColourReconnector::_containsColour8(const ClusterVector & cv, const vector & P) const { assert (P.size() == cv.size()); for (size_t i = 0; i < cv.size(); i++) { tcPPtr p = cv[i]->colParticle(); tcPPtr q = cv[P[i]]->antiColParticle(); if (_isColour8(p, q)) return true; } return false; } void ColourReconnector::_doRecoStatistical(ClusterVector & cv) const { const size_t nclusters = cv.size(); // initially, enumerate (anti)quarks as given in the cluster vector ParticleVector q, aq; for (size_t i = 0; i < nclusters; i++) { q.push_back( cv[i]->colParticle() ); aq.push_back( cv[i]->antiColParticle() ); } // annealing scheme Energy2 t, delta; Energy2 lambda = _clusterMassSum(q,aq); const unsigned _ntries = _triesPerStepFactor * nclusters; // find appropriate starting temperature by measuring the largest lambda // difference in some dry-run random rearrangements { vector typical; for (int i = 0; i < 10; i++) { const pair toswap = _shuffle(q,aq,5); ParticleVector newaq = aq; swap (newaq[toswap.first], newaq[toswap.second]); Energy2 newlambda = _clusterMassSum(q,newaq); typical.push_back( abs(newlambda - lambda) ); } t = _initTemp * Math::median(typical); } // anneal in up to _annealingSteps temperature steps for (unsigned step = 0; step < _annealingSteps; step++) { // For this temperature step, try to reconnect _ntries times. Stop the // algorithm if no successful reconnection happens. unsigned nSuccess = 0; for (unsigned it = 0; it < _ntries; it++) { // make a random rearrangement const unsigned maxtries = 10; const pair toswap = _shuffle(q,aq,maxtries); const int i = toswap.first; const int j = toswap.second; // stop here if we cannot find any allowed reconfiguration if (i == -1) break; // create a new antiquark vector with the two partons swapped ParticleVector newaq = aq; swap (newaq[i], newaq[j]); // Check if lambda would decrease. If yes, accept the reconnection. If no, // accept it only with a probability given by the current Boltzmann // factor. In the latter case we set p = 0 if the temperature is close to // 0, to avoid division by 0. Energy2 newlambda = _clusterMassSum(q,newaq); delta = newlambda - lambda; double prob = 1.0; if (delta > ZERO) prob = ( abs(t) < 1e-8*MeV2 ) ? 0.0 : exp(-delta/t); if (UseRandom::rnd() < prob) { lambda = newlambda; swap (newaq, aq); nSuccess++; } } if (nSuccess == 0) break; // reduce temperature t *= _annealingFactor; } // construct the new cluster vector ClusterVector newclusters; for (size_t i = 0; i < nclusters; i++) { ClusterPtr cl = new_ptr( Cluster( q[i], aq[i] ) ); newclusters.push_back(cl); } swap(newclusters,cv); return; } void ColourReconnector::_doRecoPlain(ClusterVector & cv) const { ClusterVector newcv = cv; // try to avoid systematic errors by randomising the reconnection order long (*p_irnd)(long) = UseRandom::irnd; random_shuffle( newcv.begin(), newcv.end(), p_irnd ); // iterate over all clusters for (CluVecIt cit = newcv.begin(); cit != newcv.end(); cit++) { // find the cluster which, if reconnected with *cit, would result in the // smallest sum of cluster masses // NB this method returns *cit if no reconnection partner can be found CluVecIt candidate = _findRecoPartner(cit, newcv); // skip this cluster if no possible reshuffling partner can be found if (candidate == cit) continue; // accept the reconnection with probability _preco. if (UseRandom::rnd() < _preco) { pair reconnected = _reconnect(*cit, *candidate); // Replace the clusters in the ClusterVector. The order of the // colour-triplet partons in the cluster vector is retained here. // replace *cit by reconnected.first *cit = reconnected.first; // replace candidate by reconnected.second *candidate = reconnected.second; } } swap(cv,newcv); return; } namespace { inline bool hasDiquark(CluVecIt cit) { for(int i = 0; i<(*cit)->numComponents(); i++) { if (DiquarkMatcher::Check(*((*cit)->particle(i)->dataPtr()))) return true; } return false; } } // Implementation of the baryonic reconnection algorithm void ColourReconnector::_doRecoBaryonic(ClusterVector & cv) const { ClusterVector newcv = cv; ClusterVector deleted; deleted.reserve(cv.size()); // try to avoid systematic errors by randomising the reconnection order long (*p_irnd)(long) = UseRandom::irnd; random_shuffle( newcv.begin(), newcv.end(), p_irnd ); // iterate over all clusters for (CluVecIt cit = newcv.begin(); cit != newcv.end(); ++cit) { //avoid clusters already containing diuarks if (hasDiquark(cit)) continue; //skip the cluster to be deleted later 3->2 cluster if (find(deleted.begin(), deleted.end(), *cit) != deleted.end()) continue; // Skip all found baryonic clusters, this biases the algorithm but implementing // something like re-reconnection is ongoing work if ((*cit)->numComponents()==3) continue; // Find a candidate suitable for reconnection CluVecIt baryonic1, baryonic2; bool isBaryonicCandidate = false; CluVecIt candidate = _findPartnerBaryonic(cit, newcv, isBaryonicCandidate, deleted, baryonic1, baryonic2); // skip this cluster if no possible reconnection partner can be found if ( !isBaryonicCandidate && candidate==cit ) continue; if ( isBaryonicCandidate && UseRandom::rnd() < _precoBaryonic ) { deleted.push_back(*baryonic2); // Function that does the reconnection from 3 -> 2 clusters ClusterPtr b1, b2; _makeBaryonicClusters(*cit,*baryonic1,*baryonic2, b1, b2); *cit = b1; *baryonic1 = b2; // Baryonic2 is easily skipped in the next loop } // Normal 2->2 Colour reconnection if ( !isBaryonicCandidate && UseRandom::rnd() < _preco ) { auto reconnected = _reconnectBaryonic(*cit, *candidate); *cit = reconnected.first; *candidate = reconnected.second; } } // create a new vector of clusters except for the ones which are "deleted" during // baryonic reconnection ClusterVector clustervector; for ( const auto & cluster : newcv ) if ( find(deleted.begin(), deleted.end(), cluster) == deleted.end() ) clustervector.push_back(cluster); swap(cv,clustervector); } namespace { double calculateRapidityRF(const Lorentz5Momentum & q1, const Lorentz5Momentum & p2) { //calculate rapidity wrt the direction of q1 //angle between the particles in the RF of cluster of q1 // calculate the z component of p2 w.r.t the direction of q1 const Energy pz = p2.vect() * q1.vect().unit(); if ( pz == ZERO ) return 0.; // Transverse momentum of p2 w.r.t the direction of q1 const Energy pt = sqrt(p2.vect().mag2() - sqr(pz)); // Transverse mass pf p2 w.r.t to the direction of q1 const Energy mtrans = sqrt(p2.mass()*p2.mass() + (pt*pt)); // Correct formula const double y2 = log((p2.t() + abs(pz))/mtrans); return ( pz < ZERO ) ? -y2 : y2; } } CluVecIt ColourReconnector::_findPartnerBaryonic( CluVecIt cl, ClusterVector & cv, bool & baryonicCand, const ClusterVector& deleted, CluVecIt &baryonic1, CluVecIt &baryonic2 ) const { using Constants::pi; using Constants::twopi; // Returns a candidate for possible reconnection CluVecIt candidate = cl; bool bcand = false; double maxrap = 0.0; double minrap = 0.0; double maxrapNormal = 0.0; double minrapNormal = 0.0; double maxsumnormal = 0.0; double maxsum = 0.0; double secondsum = 0.0; // boost into RF of cl Lorentz5Momentum cl1 = (*cl)->momentum(); const Boost boostv(-cl1.boostVector()); cl1.boost(boostv); // boost constituents of cl into RF of cl Lorentz5Momentum p1col = (*cl)->colParticle()->momentum(); Lorentz5Momentum p1anticol = (*cl)->antiColParticle()->momentum(); p1col.boost(boostv); p1anticol.boost(boostv); for (CluVecIt cit=cv.begin(); cit != cv.end(); ++cit) { //avoid looping over clusters containing diquarks if ( hasDiquark(cit) ) continue; if ( (*cit)->numComponents()==3 ) continue; if ( cit==cl ) continue; //skip the cluster to be deleted later 3->2 cluster if ( find(deleted.begin(), deleted.end(), *cit) != deleted.end() ) continue; if ( (*cl)->isBeamCluster() && (*cit)->isBeamCluster() ) continue; // stop it putting far apart clusters together if ( ( (**cl).vertex()-(**cit).vertex() ).m() >_maxDistance ) continue; const bool Colour8 = _isColour8( (*cl)->colParticle(), (*cit)->antiColParticle() ) || _isColour8( (*cit)->colParticle(), (*cl)->antiColParticle() ) ; if ( Colour8 ) continue; // boost constituents of cit into RF of cl Lorentz5Momentum p2col = (*cit)->colParticle()->momentum(); Lorentz5Momentum p2anticol = (*cit)->antiColParticle()->momentum(); p2col.boost(boostv); p2anticol.boost(boostv); // calculate the rapidity of the other constituents of the clusters // w.r.t axis of p1anticol.vect.unit const double rapq = calculateRapidityRF(p1anticol,p2col); const double rapqbar = calculateRapidityRF(p1anticol,p2anticol); // configuration for normal CR if ( rapq > 0.0 && rapqbar < 0.0 && rapq > maxrap && rapqbar < minrap ) { maxrap = rapq; minrap = rapqbar; //sum of rapidities of quarks const double normalsum = abs(rapq) + abs(rapqbar); if ( normalsum > maxsumnormal ) { maxsumnormal = normalsum; maxrapNormal = rapq; minrapNormal = rapqbar; bcand = false; candidate = cit; } } if ( rapq < 0.0 && rapqbar >0.0 && rapqbar > maxrapNormal && rapq < minrapNormal ) { maxrap = rapqbar; minrap = rapq; const double sumrap = abs(rapqbar) + abs(rapq); // first candidate gets here. If second baryonic candidate has higher Ysum than the first // one, the second candidate becomes the first one and the first the second. if (sumrap > maxsum) { if(maxsum != 0){ baryonic2 = baryonic1; baryonic1 = cit; bcand = true; } else { baryonic1 = cit; } maxsum = sumrap; } else { if (sumrap > secondsum && sumrap != maxsum) { secondsum = sumrap; bcand = true; baryonic2 = cit; } } } } if(bcand == true){ baryonicCand = true; } return candidate; } CluVecIt ColourReconnector::_findRecoPartner(CluVecIt cl, ClusterVector & cv) const { CluVecIt candidate = cl; Energy minMass = 1*TeV; for (CluVecIt cit=cv.begin(); cit != cv.end(); ++cit) { // don't even look at original cluster if(cit==cl) continue; // don't allow colour octet clusters if ( _isColour8( (*cl)->colParticle(), (*cit)->antiColParticle() ) || _isColour8( (*cit)->colParticle(), (*cl)->antiColParticle() ) ) { continue; } // stop it putting beam remnants together if((*cl)->isBeamCluster() && (*cit)->isBeamCluster()) continue; // stop it putting far apart clusters together if(((**cl).vertex()-(**cit).vertex()).m()>_maxDistance) continue; // momenta of the old clusters Lorentz5Momentum p1 = (*cl)->colParticle()->momentum() + (*cl)->antiColParticle()->momentum(); Lorentz5Momentum p2 = (*cit)->colParticle()->momentum() + (*cit)->antiColParticle()->momentum(); // momenta of the new clusters Lorentz5Momentum p3 = (*cl)->colParticle()->momentum() + (*cit)->antiColParticle()->momentum(); Lorentz5Momentum p4 = (*cit)->colParticle()->momentum() + (*cl)->antiColParticle()->momentum(); Energy oldMass = abs( p1.m() ) + abs( p2.m() ); Energy newMass = abs( p3.m() ) + abs( p4.m() ); if ( newMass < oldMass && newMass < minMass ) { minMass = newMass; candidate = cit; } } return candidate; } // forms two baryonic clusters from three clusters void ColourReconnector::_makeBaryonicClusters( ClusterPtr &c1, ClusterPtr &c2, ClusterPtr &c3, ClusterPtr &newcluster1, ClusterPtr &newcluster2) const{ //make sure they all have 2 components assert(c1->numComponents()==2); assert(c2->numComponents()==2); assert(c3->numComponents()==2); - //abandon childs + //abandon children c1->colParticle()->abandonChild(c1); c1->antiColParticle()->abandonChild(c1); c2->colParticle()->abandonChild(c2); c2->antiColParticle()->abandonChild(c2); c3->colParticle()->abandonChild(c3); c3->antiColParticle()->abandonChild(c3); newcluster1 = new_ptr(Cluster(c1->colParticle(),c2->colParticle(), c3->colParticle())); c1->colParticle()->addChild(newcluster1); c2->colParticle()->addChild(newcluster1); c3->colParticle()->addChild(newcluster1); newcluster1->setVertex(LorentzPoint()); + newcluster2 = new_ptr(Cluster(c1->antiColParticle(), c2->antiColParticle(), - c3->antiColParticle())); + c3->antiColParticle())); c1->antiColParticle()->addChild(newcluster2); c2->antiColParticle()->addChild(newcluster2); c3->antiColParticle()->addChild(newcluster2); newcluster2->setVertex(LorentzPoint()); } pair ColourReconnector::_reconnect(ClusterPtr &c1, ClusterPtr &c2) const { // choose the other possibility to form two clusters from the given // constituents assert(c1->numComponents()==2); assert(c2->numComponents()==2); int c1_col(-1),c1_anti(-1),c2_col(-1),c2_anti(-1); for(unsigned int ix=0;ix<2;++ix) { if (c1->particle(ix)->hasColour(false)) c1_col = ix; else if(c1->particle(ix)->hasColour(true )) c1_anti = ix; if (c2->particle(ix)->hasColour(false)) c2_col = ix; else if(c2->particle(ix)->hasColour(true )) c2_anti = ix; } assert(c1_col>=0&&c2_col>=0&&c1_anti>=0&&c2_anti>=0); ClusterPtr newCluster1 = new_ptr( Cluster( c1->colParticle(), c2->antiColParticle() ) ); newCluster1->setVertex(0.5*( c1->colParticle()->vertex() + c2->antiColParticle()->vertex() )); if(c1->isBeamRemnant(c1_col )) newCluster1->setBeamRemnant(0,true); if(c2->isBeamRemnant(c2_anti)) newCluster1->setBeamRemnant(1,true); ClusterPtr newCluster2 = new_ptr( Cluster( c2->colParticle(), c1->antiColParticle() ) ); newCluster2->setVertex(0.5*( c2->colParticle()->vertex() + c1->antiColParticle()->vertex() )); if(c2->isBeamRemnant(c2_col )) newCluster2->setBeamRemnant(0,true); if(c1->isBeamRemnant(c1_anti)) newCluster2->setBeamRemnant(1,true); return pair (newCluster1, newCluster2); } pair ColourReconnector::_reconnectBaryonic(ClusterPtr &c1, ClusterPtr &c2) const { // choose the other possibility to form two clusters from the given // constituents assert(c1->numComponents()==2); assert(c2->numComponents()==2); int c1_col(-1),c1_anti(-1),c2_col(-1),c2_anti(-1); for(unsigned int ix=0;ix<2;++ix) { if (c1->particle(ix)->hasColour(false)) c1_col = ix; else if(c1->particle(ix)->hasColour(true )) c1_anti = ix; if (c2->particle(ix)->hasColour(false)) c2_col = ix; else if(c2->particle(ix)->hasColour(true )) c2_anti = ix; } assert(c1_col>=0&&c2_col>=0&&c1_anti>=0&&c2_anti>=0); c1->colParticle()->abandonChild(c1); c2->antiColParticle()->abandonChild(c2); ClusterPtr newCluster1 = new_ptr( Cluster( c1->colParticle(), c2->antiColParticle() ) ); c1->colParticle()->addChild(newCluster1); c2->antiColParticle()->addChild(newCluster1); newCluster1->setVertex(0.5*( c1->colParticle()->vertex() + c2->antiColParticle()->vertex() )); if(c1->isBeamRemnant(c1_col )) newCluster1->setBeamRemnant(0,true); if(c2->isBeamRemnant(c2_anti)) newCluster1->setBeamRemnant(1,true); c1->antiColParticle()->abandonChild(c1); c2->colParticle()->abandonChild(c2); ClusterPtr newCluster2 = new_ptr( Cluster( c2->colParticle(), c1->antiColParticle() ) ); c1->antiColParticle()->addChild(newCluster2); c2->colParticle()->addChild(newCluster2); newCluster2->setVertex(0.5*( c2->colParticle()->vertex() + c1->antiColParticle()->vertex() )); if(c2->isBeamRemnant(c2_col )) newCluster2->setBeamRemnant(0,true); if(c1->isBeamRemnant(c1_anti)) newCluster2->setBeamRemnant(1,true); return pair (newCluster1, newCluster2); } pair ColourReconnector::_shuffle (const PVector & q, const PVector & aq, unsigned maxtries) const { const size_t nclusters = q.size(); assert (nclusters > 1); assert (aq.size() == nclusters); int i, j; unsigned tries = 0; bool octet; do { // find two different random integers in the range [0, nclusters) i = UseRandom::irnd( nclusters ); do { j = UseRandom::irnd( nclusters ); } while (i == j); // check if one of the two potential clusters would be a colour octet state octet = _isColour8( q[i], aq[j] ) || _isColour8( q[j], aq[i] ) ; tries++; } while (octet && tries < maxtries); if (octet) i = j = -1; return make_pair(i,j); } bool ColourReconnector::_isColour8(tcPPtr p, tcPPtr q) const { bool octet = false; // make sure we have a triplet and an anti-triplet if ( ( p->hasColour() && q->hasAntiColour() ) || ( p->hasAntiColour() && q->hasColour() ) ) { // true if p and q are originated from a colour octet if ( !p->parents().empty() && !q->parents().empty() ) { octet = ( p->parents()[0] == q->parents()[0] ) && ( p->parents()[0]->data().iColour() == PDT::Colour8 ); } // (Final) option: check if same colour8 parent // or already found an octet. if(_octetOption==0||octet) return octet; // (All) option handling more octets // by browsing particle history/colour lines. tColinePtr cline,aline; // Get colourlines form final states. if(p->hasColour() && q->hasAntiColour()) { cline = p-> colourLine(); aline = q->antiColourLine(); } else { cline = q-> colourLine(); aline = p->antiColourLine(); } // Follow the colourline of p. if ( !p->parents().empty() ) { tPPtr parent = p->parents()[0]; while (parent) { if(parent->data().iColour() == PDT::Colour8) { // Coulour8 particles should have a colour // and an anticolour line. Currently the // remnant has none of those. Since the children // of the remnant are not allowed to emit currently, // the colour octet remnant is handled by the return // statement above. The assert also catches other // colour octets without clines. If the children of // a remnant should be allowed to emit, the remnant // should get appropriate colour lines and // colour states. // See Ticket: #407 // assert(parent->colourLine()&&parent->antiColourLine()); octet = (parent-> colourLine()==cline && parent->antiColourLine()==aline); } if(octet||parent->parents().empty()) break; parent = parent->parents()[0]; } } } return octet; } void ColourReconnector::persistentOutput(PersistentOStream & os) const { os << _clreco << _preco << _precoBaryonic << _algorithm << _initTemp << _annealingFactor << _annealingSteps << _triesPerStepFactor << ounit(_maxDistance,femtometer) << _octetOption; } void ColourReconnector::persistentInput(PersistentIStream & is, int) { is >> _clreco >> _preco >> _precoBaryonic >> _algorithm >> _initTemp >> _annealingFactor >> _annealingSteps >> _triesPerStepFactor >> iunit(_maxDistance,femtometer) >> _octetOption; } void ColourReconnector::Init() { static ClassDocumentation documentation ("This class is responsible of the colour reconnection."); static Switch interfaceColourReconnection ("ColourReconnection", "Colour reconnections", &ColourReconnector::_clreco, 0, true, false); static SwitchOption interfaceColourReconnectionNo (interfaceColourReconnection, "No", "Colour reconnections off", 0); static SwitchOption interfaceColourReconnectionYes (interfaceColourReconnection, "Yes", "Colour reconnections on", 1); static Parameter interfaceMtrpAnnealingFactor ("AnnealingFactor", "The annealing factor is the ratio of the temperatures in two successive " "temperature steps.", &ColourReconnector::_annealingFactor, 0.9, 0.0, 1.0, false, false, Interface::limited); static Parameter interfaceMtrpAnnealingSteps ("AnnealingSteps", "Number of temperature steps in the statistical annealing algorithm", &ColourReconnector::_annealingSteps, 50, 1, 10000, false, false, Interface::limited); static Parameter interfaceMtrpTriesPerStepFactor ("TriesPerStepFactor", "The number of reconnection tries per temperature steps is the number of " "clusters times this factor.", &ColourReconnector::_triesPerStepFactor, 5.0, 0.0, 100.0, false, false, Interface::limited); static Parameter interfaceMtrpInitialTemp ("InitialTemperature", "Factor used to determine the initial temperature from the median of the " "energy change in a few random rearrangements.", &ColourReconnector::_initTemp, 0.1, 0.00001, 100.0, false, false, Interface::limited); static Parameter interfaceRecoProb ("ReconnectionProbability", "Probability that a found reconnection possibility is actually accepted", &ColourReconnector::_preco, 0.5, 0.0, 1.0, false, false, Interface::limited); static Parameter interfaceRecoProbBaryonic ("ReconnectionProbabilityBaryonic", "Probability that a found reconnection possibility is actually accepted", &ColourReconnector::_precoBaryonic, 0.5, 0.0, 1.0, false, false, Interface::limited); static Switch interfaceAlgorithm ("Algorithm", "Specifies the colour reconnection algorithm", &ColourReconnector::_algorithm, 0, true, false); static SwitchOption interfaceAlgorithmPlain (interfaceAlgorithm, "Plain", "Plain colour reconnection as in Herwig 2.5.0", 0); static SwitchOption interfaceAlgorithmStatistical (interfaceAlgorithm, "Statistical", "Statistical colour reconnection using simulated annealing", 1); static SwitchOption interfaceAlgorithmBaryonic (interfaceAlgorithm, "BaryonicReco", "Baryonic cluster reconnection", 2); static Parameter interfaceMaxDistance ("MaxDistance", "Maximum distance between the clusters at which to consider rearrangement" " to avoid colour reconneections of displaced vertices", &ColourReconnector::_maxDistance, femtometer, 1000.*femtometer, 0.0*femtometer, 1e100*femtometer, false, false, Interface::limited); static Switch interfaceOctetTreatment ("OctetTreatment", "Which octets are not allowed to be reconnected", &ColourReconnector::_octetOption, 0, false, false); static SwitchOption interfaceOctetTreatmentFinal (interfaceOctetTreatment, "Final", "Only prevent for the final (usuaslly non-perturbative) g -> q qbar splitting", 0); static SwitchOption interfaceOctetTreatmentAll (interfaceOctetTreatment, "All", "Prevent for all octets", 1); } diff --git a/Hadronization/HwppSelector.cc b/Hadronization/HwppSelector.cc --- a/Hadronization/HwppSelector.cc +++ b/Hadronization/HwppSelector.cc @@ -1,184 +1,248 @@ // -*- C++ -*- // // HwppSelector.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the HwppSelector class. // #include "HwppSelector.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Switch.h" +#include "ThePEG/Interface/Parameter.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "Herwig/Utilities/Kinematics.h" #include "ThePEG/Utilities/Selector.h" #include "ThePEG/Repository/UseRandom.h" #include "CheckId.h" #include #include using namespace Herwig; DescribeClass describeHwppSelector("Herwig::HwppSelector",""); IBPtr HwppSelector::clone() const { return new_ptr(*this); } IBPtr HwppSelector::fullclone() const { return new_ptr(*this); } void HwppSelector::doinit() { HadronSelector::doinit(); } void HwppSelector::persistentOutput(PersistentOStream & os) const { - os << _mode; + os << _mode << _enhanceSProb << ounit(_m0Decay,GeV) << _massMeasure; } void HwppSelector::persistentInput(PersistentIStream & is, int) { - is >> _mode; + is >> _mode >> _enhanceSProb >> iunit(_m0Decay,GeV) >> _massMeasure; } void HwppSelector::Init() { static ClassDocumentation documentation ("The HwppSelector class implements the Herwig algorithm for selecting" " the hadrons", "The hadronization used the selection algorithm described in \\cite{Kupco:1998fx}.", "%\\cite{Kupco:1998fx}\n" "\\bibitem{Kupco:1998fx}\n" " A.~Kupco,\n" " ``Cluster hadronization in HERWIG 5.9,''\n" " arXiv:hep-ph/9906412.\n" " %%CITATION = HEP-PH/9906412;%%\n" ); // put useMe() only in correct place! static Switch interfaceMode ("Mode", "Which algorithm to use", &HwppSelector::_mode, 1, false, false); static SwitchOption interfaceModeKupco (interfaceMode, "Kupco", "Use the Kupco approach", 0); static SwitchOption interfaceModeHwpp (interfaceMode, "Hwpp", "Use the Herwig approach", 1); + static Switch interfaceEnhanceSProb + ("EnhanceSProb", + "Option for enhancing strangeness", + &HwppSelector::_enhanceSProb, 0, false, false); + static SwitchOption interfaceEnhanceSProbNo + (interfaceEnhanceSProb, + "No", + "No strangeness enhancement.", + 0); + static SwitchOption interfaceEnhanceSProbScaled + (interfaceEnhanceSProb, + "Scaled", + "Scaled strangeness enhancement", + 1); + static SwitchOption interfaceEnhanceSProbExponential + (interfaceEnhanceSProb, + "Exponential", + "Exponential strangeness enhancement", + 2); + + static Switch interfaceMassMeasure + ("MassMeasure", + "Option to use different mass measures", + &HwppSelector::_massMeasure,0,false,false); + static SwitchOption interfaceMassMeasureMass + (interfaceMassMeasure, + "Mass", + "Mass Measure", + 0); + static SwitchOption interfaceMassMeasureLambda + (interfaceMassMeasure, + "Lambda", + "Lambda Measure", + 1); + + static Parameter interfaceDecayMassScale + ("DecayMassScale", + "Cluster decay mass scale", + &HwppSelector::_m0Decay, GeV, 1.0*GeV, 0.1*GeV, 50.*GeV, + false, false, Interface::limited); + } -pair HwppSelector::chooseHadronPair(const Energy cluMass,tcPDPtr par1, +pair HwppSelector::chooseHadronPair(const Energy cluMass,tcPDPtr par1, tcPDPtr par2,tcPDPtr ) const { - // if either of the input partons is a diquark don't allow diquarks to be + + // if either of the input partons is a diquark don't allow diquarks to be // produced bool diquark = !(DiquarkMatcher::Check(par1->id()) || DiquarkMatcher::Check(par2->id())); bool quark = true; - // if the Herwig algorithm + // if the Herwig algorithm if(_mode ==1) { if(UseRandom::rnd() > 1./(1.+pwtDIquark()) &&cluMass > massLightestBaryonPair(par1,par2)) { diquark = true; quark = false; } else { useMe(); diquark = false; quark = true; } } // weights for the different possibilities Energy weight, wgtsum(ZERO); // loop over all hadron pairs with the allowed flavours static vector hadrons; hadrons.clear(); for(unsigned int ix=0;ixiColour())) == 3 && !DiquarkMatcher::Check(quarktopick->id())) continue; if(!diquark && abs(int(quarktopick->iColour())) == 3 && DiquarkMatcher::Check(quarktopick->id())) continue; - HadronTable::const_iterator + HadronTable::const_iterator tit1 = table().find(make_pair(abs(par1->id()),quarktopick->id())); - HadronTable::const_iterator + HadronTable::const_iterator tit2 = table().find(make_pair(quarktopick->id(),abs(par2->id()))); // If not in table skip if(tit1 == table().end()||tit2==table().end()) continue; // tables empty skip const KupcoData & T1 = tit1->second; const KupcoData & T2 = tit2->second; if(T1.empty()||T2.empty()) continue; // if too massive skip - if(cluMass <= T1.begin()->mass + - T2.begin()->mass) continue; + if(cluMass <= T1.begin()->mass + + T2.begin()->mass) continue; // loop over the hadrons KupcoData::const_iterator H1,H2; for(H1 = T1.begin();H1 != T1.end(); ++H1) { for(H2 = T2.begin();H2 != T2.end(); ++H2) { // break if cluster too light if(cluMass < H1->mass + H2->mass) break; // calculate the weight - weight = pwt(quarktopick->id()) * H1->overallWeight * H2->overallWeight * - Kinematics::pstarTwoBodyDecay(cluMass, H1->mass, H2->mass ); + double pwtstrange; + if (quarktopick->id() == 3){ + pwtstrange = pwt(3); + if (_enhanceSProb == 1){ + double scale = double(sqr(_m0Decay/cluMass)); + pwtstrange = (20. < scale || scale < 0.) ? 0. : pow(pwtstrange,scale); + } + else if (_enhanceSProb == 2){ + Energy2 mass2; + Energy endpointmass = par1->mass() + par2->mass(); + mass2 = (_massMeasure == 0) ? sqr(cluMass) : + sqr(cluMass) - sqr(endpointmass); + double scale = double(sqr(_m0Decay)/mass2); + pwtstrange = (20. < scale || scale < 0.) ? 0. : exp(-scale); + } + weight = pwtstrange * H1->overallWeight * H2->overallWeight * + Kinematics::pstarTwoBodyDecay(cluMass, H1->mass, H2->mass ); + } + else { + weight = pwt(quarktopick->id()) * H1->overallWeight * H2->overallWeight * + Kinematics::pstarTwoBodyDecay(cluMass, H1->mass, H2->mass ); + } + int signQ = 0; assert (par1 && quarktopick); assert (par2); assert(quarktopick->CC()); - - if(CheckId::canBeHadron(par1, quarktopick->CC()) + + if(CheckId::canBeHadron(par1, quarktopick->CC()) && CheckId::canBeHadron(quarktopick, par2)) signQ = +1; - else if(CheckId::canBeHadron(par1, quarktopick) + else if(CheckId::canBeHadron(par1, quarktopick) && CheckId::canBeHadron(quarktopick->CC(), par2)) signQ = -1; else { - cerr << "Could not make sign for" << par1->id()<< " " << quarktopick->id() + cerr << "Could not make sign for" << par1->id()<< " " << quarktopick->id() << " " << par2->id() << "\n"; assert(false); } if (signQ == -1) quarktopick = quarktopick->CC(); // construct the object with the info Kupco a(quarktopick, H1->ptrData, H2->ptrData, weight); hadrons.push_back(a); wgtsum += weight; } } } - if (hadrons.empty()) + if (hadrons.empty()) return make_pair(tcPDPtr(),tcPDPtr()); // select the hadron wgtsum *= UseRandom::rnd(); unsigned int ix=0; do { wgtsum-= hadrons[ix].weight; ++ix; } while(wgtsum > ZERO && ix < hadrons.size()); - if(ix == hadrons.size() && wgtsum > ZERO) + if(ix == hadrons.size() && wgtsum > ZERO) return make_pair(tcPDPtr(),tcPDPtr()); --ix; assert(hadrons[ix].idQ); int signHad1 = signHadron(par1, hadrons[ix].idQ->CC(), hadrons[ix].hadron1); int signHad2 = signHadron(par2, hadrons[ix].idQ, hadrons[ix].hadron2); assert( signHad1 != 0 && signHad2 != 0 ); return make_pair ( signHad1 > 0 ? hadrons[ix].hadron1 : tcPDPtr(hadrons[ix].hadron1->CC()), signHad2 > 0 ? hadrons[ix].hadron2 : tcPDPtr(hadrons[ix].hadron2->CC())); } diff --git a/Hadronization/HwppSelector.h b/Hadronization/HwppSelector.h --- a/Hadronization/HwppSelector.h +++ b/Hadronization/HwppSelector.h @@ -1,144 +1,151 @@ // -*- C++ -*- // // HwppSelector.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_HwppSelector_H #define HERWIG_HwppSelector_H // // This is the declaration of the HwppSelector class. // #include "HadronSelector.h" #include "HwppSelector.fh" namespace Herwig { using namespace ThePEG; /** \ingroup hadronization * The HwppSelector class selects the hadrons produced in cluster decay using * the Herwig variant of the cluster model. * * @see \ref HwppSelectorInterfaces "The interfaces" * defined for HwppSelector. */ class HwppSelector: public HadronSelector { public: /** * The default constructor. */ - HwppSelector() : HadronSelector(1), _mode(1) + HwppSelector() : HadronSelector(1), _mode(1), _enhanceSProb(0), _m0Decay(1.*GeV) {} /** * * This method is used to choose a pair of hadrons. * - * Given the mass of a cluster and the particle pointers of its + * Given the mass of a cluster and the particle pointers of its * two (or three) constituents, this returns the pair of particle pointers of - * the two hadrons with proper flavour numbers. - * Furthermore, the first of the two hadron must have the + * the two hadrons with proper flavour numbers. + * Furthermore, the first of the two hadron must have the * constituent with par1, and the second must have the constituent with par2. * At the moment it does *nothing* in the case that also par3 is present. * * Kupco's method is used, rather than one used in FORTRAN HERWIG * The idea is to build on the fly a table of all possible pairs * of hadrons (Had1,Had2) (that we can call "cluster decay channels") - * which are kinematically above threshold and have flavour + * which are kinematically above threshold and have flavour * Had1=(par1,quarktopick->CC()), Had2=(quarktopick,par2), where quarktopick - * is the poniter of: - * --- d, u, s, c, b - * if either par1 or par2 is a diquark; - * --- d, u, s, c, b, dd, ud, uu, sd, su, ss, + * is the poniter of: + * --- d, u, s, c, b + * if either par1 or par2 is a diquark; + * --- d, u, s, c, b, dd, ud, uu, sd, su, ss, * cd, cu, cs, cc, bd, bu, bs, bc, bb * if both par1 and par2 are quarks. * The weight associated with each channel is given by the product - * of: the phase space available including the spin factor 2*J+1, - * the constant weight factor for chosen idQ, - * the octet-singlet isoscalar mixing factor, and finally + * of: the phase space available including the spin factor 2*J+1, + * the constant weight factor for chosen idQ, + * the octet-singlet isoscalar mixing factor, and finally * the singlet-decuplet weight factor. */ - pair chooseHadronPair(const Energy cluMass,tcPDPtr par1, + pair chooseHadronPair(const Energy cluMass,tcPDPtr par1, tcPDPtr par2,tcPDPtr par3 = PDPtr()) const ; public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); //@} private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ HwppSelector & operator=(const HwppSelector &) = delete; private: /** * Which algorithm to use */ unsigned int _mode; + + int _enhanceSProb; + + Energy _m0Decay; + + int _massMeasure; + }; } #endif /* HERWIG_HwppSelector_H */ diff --git a/Hadronization/Makefile.am b/Hadronization/Makefile.am --- a/Hadronization/Makefile.am +++ b/Hadronization/Makefile.am @@ -1,16 +1,16 @@ noinst_LTLIBRARIES = libHwHadronization.la libHwHadronization_la_SOURCES = \ CheckId.cc CheckId.h \ CluHadConfig.h \ Cluster.h Cluster.cc Cluster.fh \ ClusterDecayer.cc ClusterDecayer.h ClusterDecayer.fh \ ClusterFinder.cc ClusterFinder.h ClusterFinder.fh \ ClusterFissioner.cc ClusterFissioner.h ClusterFissioner.fh \ ClusterHadronizationHandler.cc ClusterHadronizationHandler.h \ ClusterHadronizationHandler.fh \ ColourReconnector.cc ColourReconnector.h ColourReconnector.fh\ HadronSelector.cc HadronSelector.h HadronSelector.fh\ Hw64Selector.cc Hw64Selector.h Hw64Selector.fh\ HwppSelector.cc HwppSelector.h HwppSelector.fh\ LightClusterDecayer.cc LightClusterDecayer.h LightClusterDecayer.fh \ -PartonSplitter.cc PartonSplitter.h PartonSplitter.fh +PartonSplitter.cc PartonSplitter.h PartonSplitter.fh diff --git a/Hadronization/PartonSplitter.cc b/Hadronization/PartonSplitter.cc --- a/Hadronization/PartonSplitter.cc +++ b/Hadronization/PartonSplitter.cc @@ -1,223 +1,568 @@ // -*- C++ -*- // // PartonSplitter.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the PartonSplitter class. // #include "PartonSplitter.h" #include #include #include #include #include #include #include #include "ThePEG/Interface/Parameter.h" #include #include #include "ThePEG/Repository/UseRandom.h" #include "Herwig/Utilities/Kinematics.h" #include #include "ClusterHadronizationHandler.h" +#include +#include #include "CheckId.h" +#include + using namespace Herwig; +static int nEvent = 1; + IBPtr PartonSplitter::clone() const { return new_ptr(*this); } IBPtr PartonSplitter::fullclone() const { return new_ptr(*this); } void PartonSplitter::persistentOutput(PersistentOStream & os) const { os << _quarkSelector << ounit(_gluonDistance,femtometer) - << _splitGluon << _splitPwtUquark << _splitPwtDquark << _splitPwtSquark; + << _splitGluon << _splitPwtUquark << _splitPwtDquark << _splitPwtSquark + << _enhanceSProb << ounit(_m0,GeV) << _massMeasure; } void PartonSplitter::persistentInput(PersistentIStream & is, int) { is >> _quarkSelector >> iunit(_gluonDistance,femtometer) - >> _splitGluon >> _splitPwtUquark >> _splitPwtDquark >> _splitPwtSquark; + >> _splitGluon >> _splitPwtUquark >> _splitPwtDquark >> _splitPwtSquark + >>_enhanceSProb >> iunit(_m0,GeV) >> _massMeasure; } -DescribeClass +DescribeClass describePartonSplitter("Herwig::PartonSplitter",""); void PartonSplitter::Init() { static ClassDocumentation documentation ("This class is reponsible of the nonperturbative splitting of partons"); static Switch interfaceSplit ("Split", "Option for different splitting options", - &PartonSplitter::_splitGluon, 1, false, false); + &PartonSplitter::_splitGluon, 0, false, false); static SwitchOption interfaceSplitDefault (interfaceSplit, "ud", "Normal cluster splitting where only u and d quarks are drawn is used.", 0); static SwitchOption interfaceSplitAll (interfaceSplit, "uds", "Alternative cluster splitting where all light quark pairs (u, d, s) can be drawn.", - 1); + 1); static Parameter interfaceSplitPwtUquark ("SplitPwtUquark", "Weight for splitting in U quarks", &PartonSplitter::_splitPwtUquark, 1, 0.0, 1.0, false, false, Interface::limited); static Parameter interfaceSplitPwtDquark ("SplitPwtDquark", "Weight for splitting in D quarks", &PartonSplitter::_splitPwtDquark, 1, 0.0, 1.0, - false, false, Interface::limited); + false, false, Interface::limited); static Parameter interfaceSplitPwtSquark ("SplitPwtSquark", "Weight for splitting in S quarks", &PartonSplitter::_splitPwtSquark, 0.5, 0.0, 1.0, false, false, Interface::limited); + static Switch interfaceEnhanceSProb + ("EnhanceSProb", + "Option to enhance the strangeness weight (MassSplit Switch needs to be on)", + &PartonSplitter::_enhanceSProb,0,false,false); + static SwitchOption interfaceEnhanceSProbNo + (interfaceEnhanceSProb, + "No", + "No scaling for strangeness", + 0); + static SwitchOption interfaceEnhanceSProbScale + (interfaceEnhanceSProb, + "Scaled", + "Power-law scaling for strangeness", + 1); + static SwitchOption interfaceEnhanceSProbExp + (interfaceEnhanceSProb, + "Exponential", + "Exponential suppression for strangeness", + 2); + + static Switch interfaceMassMeasure + ("MassMeasure", + "Option to use different mass measures", + &PartonSplitter::_massMeasure,0,false,false); + static SwitchOption interfaceMassMeasureMass + (interfaceMassMeasure, + "Mass", + "Mass Measure", + 0); + static SwitchOption interfaceMassMeasureLambda + (interfaceMassMeasure, + "Lambda", + "Lambda Measure", + 1); + + static Parameter interfaceMassScale + ("MassScale", + "Mass scale for g->qqb strangeness enhancement", + &PartonSplitter::_m0, GeV, 20.*GeV, 1.*GeV, 1000.*GeV, + false, false, Interface::limited); } void PartonSplitter::split(PVector & tagged) { // set the gluon c tau once and for all static bool first = true; + nEvent++; + if(first) { _gluonDistance = hbarc*getParticleData(ParticleID::g)->constituentMass()/ ClusterHadronizationHandler::currentHandler()->minVirtuality2(); first = false; } + // Copy incoming for the (possible sorting and) splitting + PVector particles = tagged; + // Switch to enhance strangeness + if (_enhanceSProb >= 1){ + colourSorted(particles); + } + PVector newtag; Energy2 Q02 = 0.99*sqr(getParticleData(ParticleID::g)->constituentMass()); // Loop over all of the particles in the event. - for(PVector::iterator pit = tagged.begin(); pit!=tagged.end(); ++pit) { + // Loop counter for colourSorted + for(PVector::iterator pit = particles.begin(); pit!=particles.end(); ++pit) { // only considering gluons so add other particles to list of particles if( (**pit).data().id() != ParticleID::g ) { newtag.push_back(*pit); continue; } // should not have been called for massless or space-like gluons if((**pit).momentum().m2() <= 0.0*sqr(MeV) ) { throw Exception() << "Spacelike or massless gluon m2= " << (**pit).momentum().m2()/GeV2 << "GeV2 in PartonSplitter::split()" << Exception::eventerror; } // time like gluon gets split PPtr ptrQ = PPtr(); PPtr ptrQbar = PPtr(); - splitTimeLikeGluon(*pit,ptrQ,ptrQbar); + if (_enhanceSProb == 0){ + splitTimeLikeGluon(*pit,ptrQ,ptrQbar); + } + else { + size_t i = pit - particles.begin(); + massSplitTimeLikeGluon(*pit, ptrQ, ptrQbar, i); + } + ptrQ->scale(Q02); ptrQbar->scale(Q02); (*pit)->colourLine()->addColoured(ptrQ); (*pit)->addChild(ptrQ); newtag.push_back(ptrQ); (*pit)->antiColourLine()->addAntiColoured(ptrQbar); (*pit)->addChild(ptrQbar); newtag.push_back(ptrQbar); - + // set the life length of gluon Length distance = UseRandom::rndExp(_gluonDistance); (**pit).setLifeLength((distance/(**pit).mass())*(**pit).momentum()); // assume quarks same position as gluon ptrQ ->setVertex((**pit).decayVertex()); ptrQ ->setLifeLength(Lorentz5Distance()); ptrQbar->setVertex((**pit).decayVertex()); ptrQbar->setLifeLength(Lorentz5Distance()); } swap(tagged,newtag); } -void PartonSplitter::splitTimeLikeGluon(tcPPtr ptrGluon, - PPtr & ptrQ, +void PartonSplitter::splitTimeLikeGluon(tcPPtr ptrGluon, + PPtr & ptrQ, PPtr & ptrQbar){ // select the quark flavour tPDPtr quark; long idNew=0; - switch(_splitGluon){ case 0: quark = _quarkSelector.select(UseRandom::rnd()); break; case 1: if ( ptrGluon->momentum().m() < 2.0 *getParticle(ThePEG::ParticleID::s)->data().constituentMass() ) { throw Exception() << "Impossible Kinematics in PartonSplitter::splitTimeLikeGluon()" << Exception::runerror; } // Only allow light quarks u,d,s with the probabilities - double prob_d = _splitPwtDquark; + double prob_d = _splitPwtDquark; double prob_u = _splitPwtUquark; double prob_s = _splitPwtSquark; - + int choice = UseRandom::rnd3(prob_u, prob_d, prob_s); switch(choice) { case 0: idNew = ThePEG::ParticleID::u; break; case 1: idNew = ThePEG::ParticleID::d; break; case 2: idNew = ThePEG::ParticleID::s; break; } ptrQ = getParticle(idNew); ptrQbar = getParticle(-idNew); break; } // Solve the kinematics of the two body decay G --> Q + Qbar Lorentz5Momentum momentumQ; Lorentz5Momentum momentumQbar; double cosThetaStar = UseRandom::rnd( -1.0 , 1.0 ); using Constants::pi; double phiStar = UseRandom::rnd( -pi , pi ); Energy constituentQmass; if(_splitGluon==0) { constituentQmass = quark->constituentMass(); }else{ constituentQmass = ptrQ->data().constituentMass(); } if (ptrGluon->momentum().m() < 2.0*constituentQmass) { - throw Exception() << "Impossible Kinematics in PartonSplitter::splitTimeLikeGluon()" + throw Exception() << "Impossible Kinematics in PartonSplitter::splitTimeLikeGluon()" << Exception::eventerror; } - Kinematics::twoBodyDecay(ptrGluon->momentum(), constituentQmass, - constituentQmass, cosThetaStar, phiStar, momentumQ, + Kinematics::twoBodyDecay(ptrGluon->momentum(), constituentQmass, + constituentQmass, cosThetaStar, phiStar, momentumQ, momentumQbar ); - // Create quark and anti-quark particles of the chosen flavour + // Create quark and anti-quark particles of the chosen flavour // and set they 5-momentum (the mass is the constituent one). if(_splitGluon==0) { ptrQ = new_ptr(Particle(quark )); ptrQbar = new_ptr(Particle(quark->CC())); } ptrQ ->set5Momentum( momentumQ ); ptrQbar ->set5Momentum( momentumQbar ); } void PartonSplitter::doinit() { Interfaced::doinit(); // calculate the probabilties for the gluon to branch into each quark type // based on the available phase-space, as in fortran. Energy mg=getParticleData(ParticleID::g)->constituentMass(); for( int ix=1; ix<6; ++ix ) { PDPtr quark = getParticleData(ix); Energy pcm = Kinematics::pstarTwoBodyDecay(mg,quark->constituentMass(), quark->constituentMass()); if(pcm>ZERO) _quarkSelector.insert(pcm/GeV,quark); } - if(_quarkSelector.empty()) + if(_quarkSelector.empty()) throw InitException() << "At least one quark must have constituent mass less " << "then the constituent mass of the gluon in " << "PartonSplitter::doinit()" << Exception::runerror; } + +// Method to colour sort the event and calculate the masses of the +// pre-clusters +// Convention is to have the triplet of the colour singlet first, +// then all gluons, then the antitriplet (and repeat for all the +// colour-singlets in the event) +void PartonSplitter::colourSorted(PVector& tagged) { + // Set up the output + PVector sorted; + // Reset the storage variables for doing the mass-based strangeness + // enhancement + colSingletSize.resize(0); + colSingletm2.resize(0); + // Variable to exit out of gluon loops + bool gluonLoop = false; + // Partons left to consider + PVector left = tagged; + // Loop over all the tagged particles + while (int(left.size()) > 0){ + // Pick the first particle available + PPtr p = left[0]; + // Temporary holding pen for momenta + Lorentz5Momentum tempMom(ZERO, ZERO, ZERO, ZERO); + // Don't do anything if the particle has no colour + // Simply add it back into the list of particles + if ( !p->coloured() ) { + sorted.push_back(p); + // nparts is the index of the particle after adding it to the sorted list + int nparts = 1; + // Add on the last entry of colSingletSize if the vector is not empty + // This is essentially the index the particle will have once it + // Get placed into the new colour sorted event + nparts += (colSingletSize.empty()) ? 0 : colSingletSize.back(); + tempMom += p->momentum(); + Energy2 singletm2 = singletMass2(tempMom); + // Store the number of particles and mass of the colour-singlet + colSingletSize.push_back(nparts); + colSingletm2.push_back(singletm2); + // Remove the particle from the list of particles left + left.erase(remove(left.begin(),left.end(), p), left.end()); + } + // Temporary holding pen for partons + PVector temp; + // Variable to sum end-point masses i.e. triplets and anti-triplets + Energy endPointMass = ZERO; + // If the particle in question is a gluon, search up it's antiColourLine + // until we get to the triplet. + // Note there are situations where we have a gluon loop + if ( p->hasColour() && p->hasAntiColour() ){ + // Search up its anticolour line until we get to the start particle + tPPtr partner = p->antiColourLine()->endParticle(); + // Dummy index used to loop over the anticolour line trace + tPPtr dummy = partner; + // Store the partner particle + temp.push_back(partner); + tempMom += partner->momentum(); + left.erase(remove(left.begin(),left.end(), partner), left.end()); + // While loop continues while we still reach a particle with with + // anti-colour, i.e. a gluon + while ( dummy->hasAntiColour() ){ + dummy = dummy->antiColourLine()->endParticle(); + // Check that we haven't already included it via colour indices + // If we have, it is a gluon loop. + if ( find(left.begin(), left.end(), dummy) == left.end() ) { + gluonLoop = true; + break; + } + // Store the dummy partons in reverse + temp.push_back(dummy); + tempMom += dummy->momentum(); + // Remove counted ones from the remaining list + left.erase(remove(left.begin(),left.end(), dummy), left.end()); + } + // Number of particles in this colour singlets so far + int nparts = int(temp.size()); + // Insert the new particles in the reverse order + sorted.insert(sorted.end(), temp.rbegin(), temp.rend()); + endPointMass += ((temp.back())->mass()); + // If it is a gluon loop we've already looped over the colour-singlet + // in its entirety, so we need to end early + if (gluonLoop){ + // Store the index of the final entry + nparts += (colSingletSize.empty()) ? 0 : colSingletSize.back(); + // Insert the new particles in the correct order + // i.e. triplet first, then the gluons we have seen so far + Energy2 singletm2 = singletMass2(tempMom); + colSingletSize.push_back(nparts); + colSingletm2.push_back(singletm2); + continue; + } + // If it is not a gluon loop, we now need to trace the colourLine + // down, until we reach the triplet. + // Works similarly to the first half + // Reset the temp PVector + temp.resize(0); + // Push back the particle in question + temp.push_back(p); + // tempMom hasn't been reset, add the particle we started with + tempMom += p->momentum(); + left.erase(remove(left.begin(),left.end(), p), left.end()); + // Search down its colour line until we get to the end particle + dummy = p->colourLine()->startParticle(); + temp.push_back(dummy); + tempMom += dummy->momentum(); + left.erase(remove(left.begin(),left.end(), dummy), left.end()); + while ( dummy->hasColour() ){ + dummy = dummy->colourLine()->startParticle(); + temp.push_back(dummy); + tempMom += dummy->momentum(); + left.erase(remove(left.begin(),left.end(), dummy), left.end()); + } + endPointMass += ((temp.back())->mass()); + // Update size of colour singlet + nparts += int(temp.size()); + nparts += (colSingletSize.empty()) ? 0 : colSingletSize.back(); + // Insert the new particles in the correct order + Energy2 singletm2 = singletMass2(tempMom); + sorted.insert(sorted.end(), temp.begin(), temp.end()); + endPointMass += ((sorted.back())->mass()); + colSingletSize.push_back(nparts); + // Chooses whether to use invariant mass of singlet + // or to use the lambda measure i.e. m^2 - (\sum m_i)^2 + Energy2 m2 = (_massMeasure == 0) ? singletm2 : singletm2 - sqr(endPointMass); + colSingletm2.push_back(m2); + } + // Else if it's a quark + else if ( p->hasColour() ) { + // Search up its colour line until we get to the start particle + // Works the same way as the second half of the gluon handling + tPPtr partner = p->colourLine()->startParticle(); + tPPtr dummy = partner; + temp.push_back(p); + endPointMass += ((temp.back())->mass()); + temp.push_back(partner); + tempMom += p->momentum(); + tempMom += partner->momentum(); + left.erase(remove(left.begin(),left.end(), p), left.end()); + left.erase(remove(left.begin(),left.end(), partner), left.end()); + while ( dummy->hasColour() ){ + dummy = dummy->colourLine()->startParticle(); + temp.push_back(dummy); + tempMom += dummy->momentum(); + left.erase(remove(left.begin(),left.end(), dummy), left.end()); + } + // Number of particles in this colour singlets + int nparts = int(temp.size()); + nparts += (colSingletSize.empty()) ? 0 : colSingletSize.back(); + // Insert the new particles in the correct order + Energy2 singletm2 = singletMass2(tempMom); + sorted.insert(sorted.end(), temp.begin(), temp.end()); + endPointMass += ((sorted.back())->mass()); + colSingletSize.push_back(nparts); + Energy2 m2 = (_massMeasure == 0) ? singletm2 : singletm2 - sqr(endPointMass); + colSingletm2.push_back(m2); + } + // Else it's an antiquark + else if ( p->hasAntiColour() ) { + // Search along anti-colour line, storing particles, and reversing the order + // at the end + // Works in the same way as the first half of the gluon handling + tPPtr partner = p->antiColourLine()->endParticle(); + tPPtr dummy = partner; + temp.push_back(p); + endPointMass += ((temp.back())->mass()); + temp.push_back(partner); + tempMom += p->momentum(); + tempMom += partner->momentum(); + left.erase(remove(left.begin(),left.end(), p), left.end()); + left.erase(remove(left.begin(),left.end(), partner), left.end()); + while ( dummy->hasAntiColour() ){ + dummy = dummy->antiColourLine()->endParticle(); + temp.push_back(dummy); + tempMom += dummy->momentum(); + left.erase(remove(left.begin(),left.end(), dummy), left.end()); + } + // Number of particles in this colour singlets + int nparts = int(temp.size()); + nparts += (colSingletSize.empty()) ? 0 : colSingletSize.back(); + // Insert the particles in the reverse order + Energy2 singletm2 = singletMass2(tempMom); + sorted.insert(sorted.end(), temp.rbegin(), temp.rend()); + endPointMass += ((temp.back())->mass()); + colSingletSize.push_back(nparts); + Energy2 m2 = (_massMeasure == 0) ? singletm2 : singletm2 - sqr(endPointMass); + colSingletm2.push_back(m2); + } + } + + // Check that the algorithm hasn't missed any particles. + assert( sorted.size() == tagged.size() ); + + swap(sorted,tagged); + +} + +Energy2 PartonSplitter::singletMass2(const Lorentz5Momentum& pIn){ + Energy2 e2 = sqr(pIn.t()); + Energy2 px2 = sqr(pIn.x()); + Energy2 py2 = sqr(pIn.y()); + Energy2 pz2 = sqr(pIn.z()); + return e2 - (px2 + py2 + pz2); +} + +void PartonSplitter::massSplitTimeLikeGluon(tcPPtr ptrGluon, + PPtr & ptrQ, + PPtr & ptrQbar, size_t i){ + // select the quark flavour + tPDPtr quark; + long idNew=0; + + if ( ptrGluon->momentum().m() < + 2.0 *getParticle(ThePEG::ParticleID::s)->data().constituentMass() ) { +throw Exception() << "Impossible Kinematics in PartonSplitter::massSplitTimeLikeGluon()" + << Exception::runerror; + } + // Only allow light quarks u,d,s with the probabilities + double prob_d = _splitPwtDquark; + double prob_u = _splitPwtUquark; + double prob_s = enhanceStrange(i); + int choice = UseRandom::rnd3(prob_u, prob_d, prob_s); + switch(choice) { + case 0: idNew = ThePEG::ParticleID::u; break; + case 1: idNew = ThePEG::ParticleID::d; break; + case 2: idNew = ThePEG::ParticleID::s; break; + } + ptrQ = getParticle(idNew); + ptrQbar = getParticle(-idNew); + + // Solve the kinematics of the two body decay G --> Q + Qbar + Lorentz5Momentum momentumQ; + Lorentz5Momentum momentumQbar; + double cosThetaStar = UseRandom::rnd( -1.0 , 1.0 ); + using Constants::pi; + double phiStar = UseRandom::rnd( -pi , pi ); + + Energy constituentQmass; + constituentQmass = ptrQ->data().constituentMass(); + + if (ptrGluon->momentum().m() < 2.0*constituentQmass) { + throw Exception() << "Impossible Kinematics in PartonSplitter::massSplitTimeLikeGluon()" + << Exception::eventerror; + } + Kinematics::twoBodyDecay(ptrGluon->momentum(), constituentQmass, + constituentQmass, cosThetaStar, phiStar, momentumQ, + momentumQbar ); + // Create quark and anti-quark particles of the chosen flavour + // and set they 5-momentum (the mass is the constituent one). + ptrQ ->set5Momentum( momentumQ ); + ptrQbar ->set5Momentum( momentumQbar ); + +} + +double PartonSplitter::enhanceStrange(size_t i){ + + // Get the m2 of the relevant colour singlet + // First we need to get the index of the colour-singlet the indexed i + // parton has + auto const it = lower_bound(colSingletSize.begin(), colSingletSize.end(), i); + + // Get the index of the colourSinglet mass + int indx = distance(colSingletSize.begin(), it); + + Energy2 mass2 = colSingletm2[indx]; + Energy2 m2 = _m0*_m0; + + // Scaling strangeness enhancement + if (_enhanceSProb == 1){ + // If the mass is significantly smaller than the characteristic mass, + // just set the prob to 0 + double scale = double(m2/mass2); + double prob = (20. < scale || scale < 0.) ? 0. : pow(_splitPwtSquark,scale); + return prob; + } + // Exponential strangeness enhancement + else if (_enhanceSProb == 2){ + double scale = double(m2/mass2); + double prob = (20. < scale || scale < 0.) ? 0. : exp(-scale); + return prob; + } + else return _splitPwtSquark; +} diff --git a/Hadronization/PartonSplitter.h b/Hadronization/PartonSplitter.h --- a/Hadronization/PartonSplitter.h +++ b/Hadronization/PartonSplitter.h @@ -1,172 +1,232 @@ // -*- C++ -*- // // PartonSplitter.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_PartonSplitter_H #define HERWIG_PartonSplitter_H #include "CluHadConfig.h" #include #include #include "PartonSplitter.fh" namespace Herwig { using namespace ThePEG; /** \ingroup Hadronization * \class PartonSplitter * \brief This class splits the gluons from the end of the shower. * \author Philip Stephens * \author Alberto Ribon - * - * This class does all of the nonperturbative parton splittings needed + * + * This class does all of the nonperturbative parton splittings needed * immediately after the end of the showering (both initial and final), * as very first step of the cluster hadronization. * * the quarks are attributed with different weights for the splitting * by default only the splitting in u and d quarks is allowed * the option "set /Herwig/Hadronization/PartonSplitter:Split 1" * allows for additional splitting into s quarks based on some weight - * in order for that to work the mass of the strange quark has to be changed + * in order for that to work the mass of the strange quark has to be changed * from the default value s.t. m_g > 2m_s - * + * * * * @see \ref PartonSplitterInterfaces "The interfaces" * defined for PartonSplitter. */ class PartonSplitter: public Interfaced { public: /** * Default constructor */ PartonSplitter() : _splitPwtUquark(1), _splitPwtDquark(1), _splitPwtSquark(0.5), _gluonDistance(ZERO), - _splitGluon(0) + _splitGluon(0), + _enhanceSProb(0), + _m0(10.*GeV), + _massMeasure(0) {} /** * This method does the nonperturbative splitting of: * time-like gluons. At the end of the shower the gluons should be * on a "physical" mass shell and should therefore be time-like. * @param tagged The tagged particles to be split * @return The particles which were not split and the products of splitting. */ void split(PVector & tagged); - + public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); //@} private: /** * Private and non-existent assignment operator. */ PartonSplitter & operator=(const PartonSplitter &) = delete; /** * Non-perturbatively split a time-like gluon, * if something goes wrong null pointers are returned. * @param gluon The gluon to be split * @param quark The quark produced in the splitting * @param anti The antiquark produced in the splitting */ void splitTimeLikeGluon(tcPPtr gluon, PPtr & quark, PPtr & anti); + + /** + * Non-perturbatively split a time-like gluon using a mass-based + * strangeness enhancement, + * if something goes wrong null pointers are returned. + * @param gluon The gluon to be split + * @param quark The quark produced in the splitting + * @param anti The antiquark produced in the splitting + * @param gluon's index in the colour-sorted Particle vector + */ + void massSplitTimeLikeGluon(tcPPtr gluon, PPtr & quark, PPtr & anti, + size_t i); + + /** + * Colour-sort the event into grouped colour-singlets. + * Convention: triplet - gluons - antitriplet, repeat for + * all other colour-singlets or colourless particles. + */ + void colourSorted(PVector& tagged); + // probabilities for the different quark types double _splitPwtUquark; double _splitPwtDquark; double _splitPwtSquark; private: /** * The selector to pick the type of quark */ Selector _quarkSelector; /** * A pointer to a Herwig::HadronSelector object for generating hadrons. */ /** * c tau for gluon decays */ Length _gluonDistance; - /** + /** * Flag used to determine between normal gluon splitting and alternative gluon splitting */ int _splitGluon; + /** + * Vector that stores the index in the event record of the anti-triplet + * of the colour-singlets, or of a colourless particle + */ + vector colSingletSize; + + /** + * Vector that stores the masses of the colour-singlets or of a colourless + * particle + */ + vector colSingletm2; + + /** + * Method to calculate the mass/lambda measure of a colour singlet + */ + Energy2 singletMass2(const Lorentz5Momentum& pIn); + + /** + * Method to calculate the probability for producing strange quarks + */ + double enhanceStrange(size_t i); + + /** + * Option to choose which functional form to use for strangeness + * production + */ + int _enhanceSProb; + + /** + * Characteristic mass scale for mass-based strangeness enhancement + */ + Energy _m0; + + /** + * Option to choose which mass measure to use + */ + int _massMeasure; }; } #endif /* HERWIG_PartonSplitter_H */