Page MenuHomeHEPForge
Files F8309798

No OneTemporary
Actions

View Options

diff --git a/Hadronization/ClusterDecayer.cc b/Hadronization/ClusterDecayer.cc
--- a/Hadronization/ClusterDecayer.cc
+++ b/Hadronization/ClusterDecayer.cc
@@ -1,470 +1,470 @@
// -*- 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 <ThePEG/Interface/ClassDocumentation.h>
#include <ThePEG/Interface/Reference.h>
#include <ThePEG/Interface/Parameter.h>
#include <ThePEG/Interface/Switch.h>
#include <ThePEG/Persistency/PersistentOStream.h>
#include <ThePEG/Persistency/PersistentIStream.h>
#include <ThePEG/PDT/EnumParticles.h>
#include <ThePEG/Repository/EventGenerator.h>
#include "Herwig/Utilities/Kinematics.h"
#include "CheckId.h"
#include "Cluster.h"
#include <ThePEG/Utilities/DescribeClass.h>
#include <ThePEG/Repository/UseRandom.h>
using namespace Herwig;
DescribeClass<ClusterDecayer,Interfaced>
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<ClusterDecayer> documentation
("This class is responsible for the two-body decays of normal clusters");
- static Reference<ClusterDecayer,HadronSelector>
- interfaceHadronSelector("HadronSelector",
- "A reference to the HadronSelector object",
+ static Reference<ClusterDecayer,HadronSelector>
+ 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<ClusterDecayer,bool> 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<ClusterDecayer,bool> 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<ClusterDecayer,bool> 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<ClusterDecayer,bool> 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<ClusterDecayer,double>
+ static Parameter<ClusterDecayer,double>
interfaceClSmrLight ("ClSmrLight", "cluster direction Gaussian smearing for light quark",
&ClusterDecayer::_clSmrLight, 0, 0.0 , 0.0 , 2.0,false,false,false);
- static Parameter<ClusterDecayer,double>
+ static Parameter<ClusterDecayer,double>
interfaceClSmrBottom ("ClSmrBottom", "cluster direction Gaussian smearing for b quark",
- &ClusterDecayer::_clSmrBottom, 0, 0.0 , 0.0 , 2.0,false,false,false);
-static Parameter<ClusterDecayer,double>
+ &ClusterDecayer::_clSmrBottom, 0, 0.0 , 0.0 , 2.0,false,false,false);
+static Parameter<ClusterDecayer,double>
interfaceClSmrCharm ("ClSmrCharm", "cluster direction Gaussian smearing for c quark",
- &ClusterDecayer::_clSmrCharm, 0, 0.0 , 0.0 , 2.0,false,false,false);
-static Parameter<ClusterDecayer,double>
+ &ClusterDecayer::_clSmrCharm, 0, 0.0 , 0.0 , 2.0,false,false,false);
+static Parameter<ClusterDecayer,double>
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<ClusterDecayer,bool> 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<ClusterDecayer,unsigned int> 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<PPtr,PPtr> prod = decayIntoTwoHadrons(*it);
(*it)->addChild(prod.first);
(*it)->addChild(prod.second);
finalhadrons.push_back(prod.first);
finalhadrons.push_back(prod.second);
}
}
}
pair<PPtr,PPtr> 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<PPtr,PPtr>();
}
-
+
// 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<double> 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<tcPDPtr,tcPDPtr> dataPair
+ pair<tcPDPtr,tcPDPtr> dataPair
= _hadronsSelector->chooseHadronPair(ptr->mass(),
ptr1data,
ptr2data);
- if(dataPair.first == tcPDPtr() ||
+ if(dataPair.first == tcPDPtr() ||
dataPair.second == tcPDPtr()) return pair<PPtr,PPtr>();
// Create the two hadron particle objects with the specified id.
PPtr ptrHad1,ptrHad2;
// produce the hadrons on mass shell
if(_onshell) {
ptrHad1 = dataPair.first ->produceParticle(dataPair.first ->mass());
ptrHad2 = dataPair.second->produceParticle(dataPair.second->mass());
}
// produce the hadrons with mass given by the mass generator
else {
unsigned int ntry(0);
do {
ptrHad1 = dataPair.first ->produceParticle();
ptrHad2 = dataPair.second->produceParticle();
++ntry;
}
while(ntry<_masstry&&ptrHad1->mass()+ptrHad2->mass()>ptr->mass());
// if fails produce on shell and issue warning (should never happen??)
if( ptrHad1->mass() + ptrHad2->mass() > ptr->mass() ) {
generator()->log() << "Failed to produce off-shell hadrons in "
<< "ClusterDecayer::decayIntoTwoHadrons producing hadrons "
<< "on-shell" << endl;
ptrHad1 = dataPair.first ->produceParticle(dataPair.first ->mass());
ptrHad2 = dataPair.second->produceParticle(dataPair.second->mass());
}
}
-
+
if (!ptrHad1 || !ptrHad2) {
ostringstream s;
s << "ClusterDecayer::decayIntoTwoHadrons ***Cannot create the two hadrons***"
- << dataPair.first ->PDGName() << " and "
+ << dataPair.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<PPtr,PPtr>(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 <ThePEG/Interface/Interfaced.h>
#include <ThePEG/EventRecord/Step.h>
#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<PPtr,PPtr> 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<HadronSelector>::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,1112 @@
// -*- 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 <ThePEG/Interface/ClassDocumentation.h>
#include <ThePEG/Interface/Reference.h>
#include <ThePEG/Interface/Parameter.h>
#include <ThePEG/Interface/Switch.h>
#include <ThePEG/Persistency/PersistentOStream.h>
#include <ThePEG/Persistency/PersistentIStream.h>
#include <ThePEG/PDT/EnumParticles.h>
#include "Herwig/Utilities/Kinematics.h"
#include "CheckId.h"
#include "Cluster.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include <ThePEG/Utilities/DescribeClass.h>
using namespace Herwig;
DescribeClass<ClusterFissioner,Interfaced>
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<ClusterFissioner> documentation
("Class responsibles for chopping up the clusters");
- static Reference<ClusterFissioner,HadronSelector>
- interfaceHadronSelector("HadronSelector",
- "A reference to the HadronSelector object",
+ static Reference<ClusterFissioner,HadronSelector>
+ 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<ClusterFissioner,Energy>
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<ClusterFissioner,Energy>
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<ClusterFissioner,Energy>
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<ClusterFissioner,Energy>
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<ClusterFissioner,double>
interfaceClPowLight ("ClPowLight","cluster mass exponent for light quarks",
&ClusterFissioner::_clPowLight, 0, 2.0, 0.0, 10.0,false,false,false);
static Parameter<ClusterFissioner,double>
interfaceClPowBottom ("ClPowBottom","cluster mass exponent for b quarks",
&ClusterFissioner::_clPowBottom, 0, 2.0, 0.0, 10.0,false,false,false);
static Parameter<ClusterFissioner,double>
interfaceClPowCharm ("ClPowCharm","cluster mass exponent for c quarks",
&ClusterFissioner::_clPowCharm, 0, 2.0, 0.0, 10.0,false,false,false);
static Parameter<ClusterFissioner,double>
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<ClusterFissioner,double>
interfacePSplitLight ("PSplitLight","cluster mass splitting param for light quarks",
&ClusterFissioner::_pSplitLight, 0, 1.0, 0.0, 10.0,false,false,false);
static Parameter<ClusterFissioner,double>
interfacePSplitBottom ("PSplitBottom","cluster mass splitting param for b quarks",
&ClusterFissioner::_pSplitBottom, 0, 1.0, 0.0, 10.0,false,false,false);
static Parameter<ClusterFissioner,double>
interfacePSplitCharm ("PSplitCharm","cluster mass splitting param for c quarks",
&ClusterFissioner::_pSplitCharm, 0, 1.0, 0.0, 10.0,false,false,false);
static Parameter<ClusterFissioner,double>
interfacePSplitExotic ("PSplitExotic","cluster mass splitting param for exotic quarks",
&ClusterFissioner::_pSplitExotic, 0, 1.0, 0.0, 10.0,false,false,false);
+ static Switch<ClusterFissioner,int> 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<ClusterFissioner,double> interfaceFissionPwtUquark
+ ("FissionPwtUquark",
+ "Weight for fission in U quarks",
+ &ClusterFissioner::_fissionPwtUquark, 1, 0.0, 1.0,
+ false, false, Interface::limited);
+ static Parameter<ClusterFissioner,double> interfaceFissionPwtDquark
+ ("FissionPwtDquark",
+ "Weight for fission in D quarks",
+ &ClusterFissioner::_fissionPwtDquark, 1, 0.0, 1.0,
+ false, false, Interface::limited);
+ static Parameter<ClusterFissioner,double> interfaceFissionPwtSquark
+ ("FissionPwtSquark",
+ "Weight for fission in S quarks",
+ &ClusterFissioner::_fissionPwtSquark, 0.5, 0.0, 1.0,
+ false, false, Interface::limited);
+
+
static Switch<ClusterFissioner,int> 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<ClusterFissioner,Energy> 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<ClusterFissioner,Tension> 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<ClusterFissioner,int> 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<ClusterFissioner,int> 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<ClusterFissioner,Energy> 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<ClusterPtr> splitClusters;
- for(ClusterVector::iterator it = clusters.begin() ;
+
+ stack<ClusterPtr> 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<ClusterPtr> & clusterStack,
- ClusterVector &clusters, tPVector & finalhadrons,
+void ClusterFissioner::cut(stack<ClusterPtr> & 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<ClusterPtr>(ct.first.first);
ClusterPtr two = dynamic_ptr_cast<ClusterPtr>(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;
+ default :
+ assert(false);
+ }
int choice = UseRandom::rnd3(prob_u, prob_d, prob_s);
long idNew = 0;
switch (choice) {
- case 0: idNew = ThePEG::ParticleID::u; break;
+ case 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;ix<min(clu->numComponents(),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 <ThePEG/Interface/Interfaced.h>
#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<ClusterPtr> &,
ClusterVector&, tPVector & finalhadrons, bool softUEisOn);
public:
/**
* Definition for easy passing of two particles.
*/
typedef pair<PPtr,PPtr> PPair;
/**
* Definition for use in the cut function.
*/
typedef pair<PPair,PPair> 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 <ThePEG/Interface/ClassDocumentation.h>
#include <ThePEG/Persistency/PersistentOStream.h>
#include <ThePEG/Persistency/PersistentIStream.h>
-#include <ThePEG/Interface/Switch.h>
-#include <ThePEG/Interface/Parameter.h>
+#include <ThePEG/Interface/Switch.h>
+#include <ThePEG/Interface/Parameter.h>
#include <ThePEG/Interface/Reference.h>
#include <ThePEG/Handlers/EventHandler.h>
#include <ThePEG/Handlers/Hint.h>
#include <ThePEG/PDT/ParticleData.h>
#include <ThePEG/EventRecord/Particle.h>
#include <ThePEG/EventRecord/Step.h>
#include <ThePEG/PDT/PDT.h>
#include <ThePEG/PDT/EnumParticles.h>
#include <ThePEG/Utilities/Throw.h>
#include "Herwig/Utilities/EnumParticles.h"
#include "CluHadConfig.h"
-#include "Cluster.h"
+#include "Cluster.h"
#include <ThePEG/Utilities/DescribeClass.h>
using namespace Herwig;
ClusterHadronizationHandler * ClusterHadronizationHandler::currentHandler_ = 0;
DescribeClass<ClusterHadronizationHandler,HadronizationHandler>
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<ClusterHadronizationHandler> 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<ClusterHadronizationHandler,PartonSplitter>
- interfacePartonSplitter("PartonSplitter",
- "A reference to the PartonSplitter object",
+ static Reference<ClusterHadronizationHandler,PartonSplitter>
+ interfacePartonSplitter("PartonSplitter",
+ "A reference to the PartonSplitter object",
&Herwig::ClusterHadronizationHandler::_partonSplitter,
false, false, true, false);
- static Reference<ClusterHadronizationHandler,ClusterFinder>
- interfaceClusterFinder("ClusterFinder",
- "A reference to the ClusterFinder object",
+ static Reference<ClusterHadronizationHandler,ClusterFinder>
+ interfaceClusterFinder("ClusterFinder",
+ "A reference to the ClusterFinder object",
&Herwig::ClusterHadronizationHandler::_clusterFinder,
false, false, true, false);
- static Reference<ClusterHadronizationHandler,ColourReconnector>
- interfaceColourReconnector("ColourReconnector",
- "A reference to the ColourReconnector object",
+ static Reference<ClusterHadronizationHandler,ColourReconnector>
+ interfaceColourReconnector("ColourReconnector",
+ "A reference to the ColourReconnector object",
&Herwig::ClusterHadronizationHandler::_colourReconnector,
false, false, true, false);
- static Reference<ClusterHadronizationHandler,ClusterFissioner>
- interfaceClusterFissioner("ClusterFissioner",
- "A reference to the ClusterFissioner object",
+ static Reference<ClusterHadronizationHandler,ClusterFissioner>
+ interfaceClusterFissioner("ClusterFissioner",
+ "A reference to the ClusterFissioner object",
&Herwig::ClusterHadronizationHandler::_clusterFissioner,
false, false, true, false);
- static Reference<ClusterHadronizationHandler,LightClusterDecayer>
- interfaceLightClusterDecayer("LightClusterDecayer",
- "A reference to the LightClusterDecayer object",
+ static Reference<ClusterHadronizationHandler,LightClusterDecayer>
+ interfaceLightClusterDecayer("LightClusterDecayer",
+ "A reference to the LightClusterDecayer object",
&Herwig::ClusterHadronizationHandler::_lightClusterDecayer,
false, false, true, false);
- static Reference<ClusterHadronizationHandler,ClusterDecayer>
- interfaceClusterDecayer("ClusterDecayer",
- "A reference to the ClusterDecayer object",
+ static Reference<ClusterHadronizationHandler,ClusterDecayer>
+ interfaceClusterDecayer("ClusterDecayer",
+ "A reference to the ClusterDecayer object",
&Herwig::ClusterHadronizationHandler::_clusterDecayer,
false, false, true, false);
- static Parameter<ClusterHadronizationHandler,Energy2> interfaceMinVirtuality2
+ static Parameter<ClusterHadronizationHandler,Energy2> 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<ClusterHadronizationHandler,Length> interfaceMaxDisplacement
+
+ static Parameter<ClusterHadronizationHandler,Length> 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<ClusterHadronizationHandler,StepHandler> interfaceUnderlyingEventHandler
("UnderlyingEventHandler",
"Pointer to the handler for the Underlying Event. "
"Set to NULL to disable.",
&ClusterHadronizationHandler::_underlyingEventHandler, false, false, true, true, false);
static Switch<ClusterHadronizationHandler,bool> 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<PPtr> & 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;ix<currentlist.size();++ix) {
if(currentlist[ix]->scale()<Q02) currentlist[ix]->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<PPtr> allDecendants;
for (tPVector::const_iterator it = tagged.begin();
it != tagged.end(); ++it) {
extractChildren(*it, allDecendants);
}
for(set<PPtr>::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 <ThePEG/Utilities/DescribeClass.h>
#include <ThePEG/Repository/UseRandom.h>
#include <ThePEG/PDT/StandardMatchers.h>
#include <ThePEG/Persistency/PersistentOStream.h>
#include <ThePEG/Persistency/PersistentIStream.h>
#include <ThePEG/Interface/Switch.h>
#include <ThePEG/Interface/Parameter.h>
#include "Herwig/Utilities/Maths.h"
using namespace Herwig;
using CluVecIt = ColourReconnector::CluVecIt;
using Constants::pi;
using Constants::twopi;
DescribeClass<ColourReconnector,Interfaced>
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<size_t> & 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<Energy2> typical;
for (int i = 0; i < 10; i++) {
const pair <int,int> 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 <int,int> 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 <ClusterPtr,ClusterPtr> 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 <ClusterPtr,ClusterPtr>
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 <ClusterPtr,ClusterPtr> (newCluster1, newCluster2);
}
pair <ClusterPtr,ClusterPtr>
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 <ClusterPtr,ClusterPtr> (newCluster1, newCluster2);
}
pair <int,int> 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<ColourReconnector> documentation
("This class is responsible of the colour reconnection.");
static Switch<ColourReconnector,int> 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<ColourReconnector,double> 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<ColourReconnector,unsigned> interfaceMtrpAnnealingSteps
("AnnealingSteps",
"Number of temperature steps in the statistical annealing algorithm",
&ColourReconnector::_annealingSteps, 50, 1, 10000,
false, false, Interface::limited);
static Parameter<ColourReconnector,double> 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<ColourReconnector,double> 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<ColourReconnector,double> 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<ColourReconnector,double> 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<ColourReconnector,int> 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<ColourReconnector,Length> 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<ColourReconnector,unsigned int> 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 <cassert>
#include <ThePEG/Utilities/DescribeClass.h>
using namespace Herwig;
DescribeClass<HwppSelector,HadronSelector>
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<HwppSelector> 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<HwppSelector,unsigned int> 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<HwppSelector,int> 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<HwppSelector,int> 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<HwppSelector,Energy> interfaceDecayMassScale
+ ("DecayMassScale",
+ "Cluster decay mass scale",
+ &HwppSelector::_m0Decay, GeV, 1.0*GeV, 0.1*GeV, 50.*GeV,
+ false, false, Interface::limited);
+
}
-pair<tcPDPtr,tcPDPtr> HwppSelector::chooseHadronPair(const Energy cluMass,tcPDPtr par1,
+pair<tcPDPtr,tcPDPtr> 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<Kupco> hadrons;
hadrons.clear();
for(unsigned int ix=0;ix<partons().size();++ix) {
tcPDPtr quarktopick = partons()[ix];
if(!quark && abs(int(quarktopick->iColour())) == 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<tcPDPtr,tcPDPtr> chooseHadronPair(const Energy cluMass,tcPDPtr par1,
+ pair<tcPDPtr,tcPDPtr> 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,555 @@
// -*- C++ -*-
//
// PartonSplitter.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the PartonSplitter class.
//
#include "PartonSplitter.h"
#include <ThePEG/Interface/ClassDocumentation.h>
#include <ThePEG/Interface/Reference.h>
#include <ThePEG/Interface/Switch.h>
#include <ThePEG/Persistency/PersistentOStream.h>
#include <ThePEG/Persistency/PersistentIStream.h>
#include <ThePEG/PDT/EnumParticles.h>
#include <ThePEG/EventRecord/Step.h>
#include "ThePEG/Interface/Parameter.h"
#include <ThePEG/Repository/EventGenerator.h>
#include <ThePEG/Repository/CurrentGenerator.h>
#include "ThePEG/Repository/UseRandom.h"
#include "Herwig/Utilities/Kinematics.h"
#include <ThePEG/Utilities/DescribeClass.h>
#include "ClusterHadronizationHandler.h"
+#include <ThePEG/EventRecord/Particle.h>
+#include <ThePEG/PDT/PDT.h>
#include "CheckId.h"
+
using namespace Herwig;
IBPtr PartonSplitter::clone() const {
return new_ptr(*this);
}
IBPtr PartonSplitter::fullclone() const {
return new_ptr(*this);
}
void PartonSplitter::persistentOutput(PersistentOStream & os) const {
os << _quarkSelector << ounit(_gluonDistance,femtometer)
- << _splitGluon << _splitPwtUquark << _splitPwtDquark << _splitPwtSquark;
+ << _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<PartonSplitter,Interfaced>
+DescribeClass<PartonSplitter,Interfaced>
describePartonSplitter("Herwig::PartonSplitter","");
void PartonSplitter::Init() {
static ClassDocumentation<PartonSplitter> documentation
("This class is reponsible of the nonperturbative splitting of partons");
static Switch<PartonSplitter,int> 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<PartonSplitter,double> interfaceSplitPwtUquark
("SplitPwtUquark",
"Weight for splitting in U quarks",
&PartonSplitter::_splitPwtUquark, 1, 0.0, 1.0,
false, false, Interface::limited);
static Parameter<PartonSplitter,double> 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<PartonSplitter,double> interfaceSplitPwtSquark
("SplitPwtSquark",
"Weight for splitting in S quarks",
&PartonSplitter::_splitPwtSquark, 0.5, 0.0, 1.0,
false, false, Interface::limited);
+ static Switch<PartonSplitter,int> 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<PartonSplitter,int> 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<PartonSplitter,Energy> interfaceMassScale
+ ("MassScale",
+ "Mass scale for g->qqb strangeness enhancement",
+ &PartonSplitter::_m0, GeV, 20.*GeV, 1.*GeV, 1000.*GeV,
+ false, false, Interface::limited);
}
void PartonSplitter::split(PVector & tagged) {
// set the gluon c tau once and for all
static bool first = true;
+
if(first) {
_gluonDistance = hbarc*getParticleData(ParticleID::g)->constituentMass()/
ClusterHadronizationHandler::currentHandler()->minVirtuality2();
first = false;
}
+ // Copy incoming for the (possible sorting and) splitting
+ PVector particles = tagged;
+ // Switch to enhance strangeness
+ if (_enhanceSProb >= 1){
+ colourSorted(particles);
+ }
+
PVector newtag;
Energy2 Q02 = 0.99*sqr(getParticleData(ParticleID::g)->constituentMass());
// Loop over all of the particles in the event.
- 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 = tempMom.m2();
+ // Store the number of particles and mass of the colour-singlet
+ _colSingletSize.push_back(nparts);
+ _colSingletm2.push_back(singletm2);
+ // Remove the particle from the list of particles left
+ left.erase(remove(left.begin(),left.end(), p), left.end());
+ }
+ // Temporary holding pen for partons
+ PVector temp;
+ // Variable to sum end-point masses i.e. triplets and anti-triplets
+ Energy endPointMass = ZERO;
+ // If the particle in question is a gluon, search up it's antiColourLine
+ // until we get to the triplet.
+ // Note there are situations where we have a gluon loop
+ if ( p->hasColour() && p->hasAntiColour() ){
+ // Search up its anticolour line until we get to the start particle
+ tPPtr partner = p->antiColourLine()->endParticle();
+ // Dummy index used to loop over the anticolour line trace
+ tPPtr dummy = partner;
+ // Store the partner particle
+ temp.push_back(partner);
+ tempMom += partner->momentum();
+ left.erase(remove(left.begin(),left.end(), partner), left.end());
+ // While loop continues while we still reach a particle with with
+ // anti-colour, i.e. a gluon
+ while ( dummy->hasAntiColour() ){
+ dummy = dummy->antiColourLine()->endParticle();
+ // Check that we haven't already included it via colour indices
+ // If we have, it is a gluon loop.
+ if ( find(left.begin(), left.end(), dummy) == left.end() ) {
+ gluonLoop = true;
+ break;
+ }
+ // Store the dummy partons in reverse
+ temp.push_back(dummy);
+ tempMom += dummy->momentum();
+ // Remove counted ones from the remaining list
+ left.erase(remove(left.begin(),left.end(), dummy), left.end());
+ }
+ // Number of particles in this colour singlets so far
+ int nparts = int(temp.size());
+ // Insert the new particles in the reverse order
+ sorted.insert(sorted.end(), temp.rbegin(), temp.rend());
+ endPointMass += ((temp.back())->mass());
+ // If it is a gluon loop we've already looped over the colour-singlet
+ // in its entirety, so we need to end early
+ if (gluonLoop){
+ // Store the index of the final entry
+ nparts += (_colSingletSize.empty()) ? 0 : _colSingletSize.back();
+ // Insert the new particles in the correct order
+ // i.e. triplet first, then the gluons we have seen so far
+ Energy2 singletm2 = tempMom.m2();
+ _colSingletSize.push_back(nparts);
+ _colSingletm2.push_back(singletm2);
+ continue;
+ }
+ // If it is not a gluon loop, we now need to trace the colourLine
+ // down, until we reach the triplet.
+ // Works similarly to the first half
+ // Reset the temp PVector
+ temp.resize(0);
+ // Push back the particle in question
+ temp.push_back(p);
+ // tempMom hasn't been reset, add the particle we started with
+ tempMom += p->momentum();
+ left.erase(remove(left.begin(),left.end(), p), left.end());
+ // Search down its colour line until we get to the end particle
+ dummy = p->colourLine()->startParticle();
+ temp.push_back(dummy);
+ tempMom += dummy->momentum();
+ left.erase(remove(left.begin(),left.end(), dummy), left.end());
+ while ( dummy->hasColour() ){
+ dummy = dummy->colourLine()->startParticle();
+ temp.push_back(dummy);
+ tempMom += dummy->momentum();
+ left.erase(remove(left.begin(),left.end(), dummy), left.end());
+ }
+ endPointMass += ((temp.back())->mass());
+ // Update size of colour singlet
+ nparts += int(temp.size());
+ nparts += (_colSingletSize.empty()) ? 0 : _colSingletSize.back();
+ // Insert the new particles in the correct order
+ Energy2 singletm2 = tempMom.m2();
+ sorted.insert(sorted.end(), temp.begin(), temp.end());
+ endPointMass += ((sorted.back())->mass());
+ _colSingletSize.push_back(nparts);
+ // Chooses whether to use invariant mass of singlet
+ // or to use the lambda measure i.e. m^2 - (\sum m_i)^2
+ Energy2 m2 = (_massMeasure == 0) ? singletm2 : singletm2 - sqr(endPointMass);
+ _colSingletm2.push_back(m2);
+ }
+ // Else if it's a quark
+ else if ( p->hasColour() ) {
+ // Search up its colour line until we get to the start particle
+ // Works the same way as the second half of the gluon handling
+ tPPtr partner = p->colourLine()->startParticle();
+ tPPtr dummy = partner;
+ temp.push_back(p);
+ endPointMass += ((temp.back())->mass());
+ temp.push_back(partner);
+ tempMom += p->momentum();
+ tempMom += partner->momentum();
+ left.erase(remove(left.begin(),left.end(), p), left.end());
+ left.erase(remove(left.begin(),left.end(), partner), left.end());
+ while ( dummy->hasColour() ){
+ dummy = dummy->colourLine()->startParticle();
+ temp.push_back(dummy);
+ tempMom += dummy->momentum();
+ left.erase(remove(left.begin(),left.end(), dummy), left.end());
+ }
+ // Number of particles in this colour singlets
+ int nparts = int(temp.size());
+ nparts += (_colSingletSize.empty()) ? 0 : _colSingletSize.back();
+ // Insert the new particles in the correct order
+ Energy2 singletm2 = tempMom.m2();
+ sorted.insert(sorted.end(), temp.begin(), temp.end());
+ endPointMass += ((sorted.back())->mass());
+ _colSingletSize.push_back(nparts);
+ Energy2 m2 = (_massMeasure == 0) ? singletm2 : singletm2 - sqr(endPointMass);
+ _colSingletm2.push_back(m2);
+ }
+ // Else it's an antiquark
+ else if ( p->hasAntiColour() ) {
+ // Search along anti-colour line, storing particles, and reversing the order
+ // at the end
+ // Works in the same way as the first half of the gluon handling
+ tPPtr partner = p->antiColourLine()->endParticle();
+ tPPtr dummy = partner;
+ temp.push_back(p);
+ endPointMass += ((temp.back())->mass());
+ temp.push_back(partner);
+ tempMom += p->momentum();
+ tempMom += partner->momentum();
+ left.erase(remove(left.begin(),left.end(), p), left.end());
+ left.erase(remove(left.begin(),left.end(), partner), left.end());
+ while ( dummy->hasAntiColour() ){
+ dummy = dummy->antiColourLine()->endParticle();
+ temp.push_back(dummy);
+ tempMom += dummy->momentum();
+ left.erase(remove(left.begin(),left.end(), dummy), left.end());
+ }
+ // Number of particles in this colour singlets
+ int nparts = int(temp.size());
+ nparts += (_colSingletSize.empty()) ? 0 : _colSingletSize.back();
+ // Insert the particles in the reverse order
+ Energy2 singletm2 = tempMom.m2();
+ sorted.insert(sorted.end(), temp.rbegin(), temp.rend());
+ endPointMass += ((temp.back())->mass());
+ _colSingletSize.push_back(nparts);
+ Energy2 m2 = (_massMeasure == 0) ? singletm2 : singletm2 - sqr(endPointMass);
+ _colSingletm2.push_back(m2);
+ }
+ }
+
+ // Check that the algorithm hasn't missed any particles.
+ assert( sorted.size() == tagged.size() );
+
+ swap(sorted,tagged);
+
+}
+
+void PartonSplitter::massSplitTimeLikeGluon(tcPPtr ptrGluon,
+ PPtr & ptrQ,
+ PPtr & ptrQbar, size_t i){
+ // select the quark flavour
+ tPDPtr quark;
+ long idNew=0;
+
+ if ( ptrGluon->momentum().m() <
+ 2.0 *getParticle(ThePEG::ParticleID::s)->data().constituentMass() ) {
+throw Exception() << "Impossible Kinematics in PartonSplitter::massSplitTimeLikeGluon()"
+ << Exception::runerror;
+ }
+ // Only allow light quarks u,d,s with the probabilities
+ double prob_d = _splitPwtDquark;
+ double prob_u = _splitPwtUquark;
+ double prob_s = enhanceStrange(i);
+ int choice = UseRandom::rnd3(prob_u, prob_d, prob_s);
+ switch(choice) {
+ case 0: idNew = ThePEG::ParticleID::u; break;
+ case 1: idNew = ThePEG::ParticleID::d; break;
+ case 2: idNew = ThePEG::ParticleID::s; break;
+ }
+ ptrQ = getParticle(idNew);
+ ptrQbar = getParticle(-idNew);
+
+ // Solve the kinematics of the two body decay G --> Q + Qbar
+ Lorentz5Momentum momentumQ;
+ Lorentz5Momentum momentumQbar;
+ double cosThetaStar = UseRandom::rnd( -1.0 , 1.0 );
+ using Constants::pi;
+ double phiStar = UseRandom::rnd( -pi , pi );
+
+ Energy constituentQmass;
+ constituentQmass = ptrQ->data().constituentMass();
+
+ if (ptrGluon->momentum().m() < 2.0*constituentQmass) {
+ throw Exception() << "Impossible Kinematics in PartonSplitter::massSplitTimeLikeGluon()"
+ << Exception::eventerror;
+ }
+ Kinematics::twoBodyDecay(ptrGluon->momentum(), constituentQmass,
+ constituentQmass, cosThetaStar, phiStar, momentumQ,
+ momentumQbar );
+ // Create quark and anti-quark particles of the chosen flavour
+ // and set they 5-momentum (the mass is the constituent one).
+ ptrQ ->set5Momentum( momentumQ );
+ ptrQbar ->set5Momentum( momentumQbar );
+
+}
+
+double PartonSplitter::enhanceStrange(size_t i){
+
+ // Get the m2 of the relevant colour singlet
+ // First we need to get the index of the colour-singlet the indexed i
+ // parton has
+ auto const it = lower_bound(_colSingletSize.begin(), _colSingletSize.end(), i);
+
+ // Get the index of the colourSinglet mass
+ int indx = distance(_colSingletSize.begin(), it);
+
+ Energy2 mass2 = _colSingletm2[indx];
+ Energy2 m2 = _m0*_m0;
+
+ // Scaling strangeness enhancement
+ if (_enhanceSProb == 1){
+ // If the mass is significantly smaller than the characteristic mass,
+ // just set the prob to 0
+ double scale = double(m2/mass2);
+ return (20. < scale || scale < 0.) ? 0. : pow(_splitPwtSquark,scale);
+ }
+ // Exponential strangeness enhancement
+ else if (_enhanceSProb == 2){
+ double scale = double(m2/mass2);
+ return (20. < scale || scale < 0.) ? 0. : exp(-scale);
+ }
+ else
+ return _splitPwtSquark;
+}
diff --git a/Hadronization/PartonSplitter.h b/Hadronization/PartonSplitter.h
--- a/Hadronization/PartonSplitter.h
+++ b/Hadronization/PartonSplitter.h
@@ -1,172 +1,222 @@
// -*- 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 <ThePEG/Interface/Interfaced.h>
#include <ThePEG/Utilities/Selector.h>
#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);
+
+ /**
+ * Method to calculate the probability for producing strange quarks
+ */
+ double enhanceStrange(size_t i);
+
+private:
+
// probabilities for the different quark types
double _splitPwtUquark;
double _splitPwtDquark;
double _splitPwtSquark;
-
-private:
-
/**
* The selector to pick the type of quark
*/
Selector<PDPtr,double> _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<int> _colSingletSize;
+
+ /**
+ * Vector that stores the masses of the colour-singlets or of a colourless
+ * particle
+ */
+ vector<Energy2> _colSingletm2;
+
+ /**
+ * 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 */

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 4:34 PM (22 h, 1 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023465
Default Alt Text
(177 KB)

Event Timeline