diff --git a/Hadronization/ClusterFissioner.cc b/Hadronization/ClusterFissioner.cc
--- a/Hadronization/ClusterFissioner.cc
+++ b/Hadronization/ClusterFissioner.cc
@@ -1,2644 +1,2642 @@
// -*- C++ -*-
//
// ClusterFissioner.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// 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 "Cluster.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include <ThePEG/Utilities/DescribeClass.h>
#include "ThePEG/Interface/ParMap.h"
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/lu.hpp>
using namespace Herwig;
DescribeClass<ClusterFissioner,Interfaced>
describeClusterFissioner("Herwig::ClusterFissioner","Herwig.so");
ClusterFissioner::ClusterFissioner() :
_clMaxLight(3.35*GeV),
_clMaxExotic(3.35*GeV),
_clPowLight(2.0),
_clPowExotic(2.0),
_pSplitLight(1.0),
_pSplitExotic(1.0),
_phaseSpaceWeights(false),
_dim(4),
_fissionCluster(0),
_kinematicThresholdChoice(0),
_pwtDIquark(0.0),
_diquarkClusterFission(0),
_btClM(1.0*GeV),
_iopRem(1),
_kappa(1.0e15*GeV/meter),
_enhanceSProb(0),
_m0Fission(2.*GeV),
_massMeasure(0),
_probPowFactor(4.0),
_probShift(0.0),
_kinThresholdShift(1.0*sqr(GeV)),
_strictDiquarkKinematics(0),
_covariantBoost(false),
_allowHadronFinalStates(2),
_massSampler(0),
_phaseSpaceSampler(0),
_matrixElement(0),
_fissionApproach(0)
{}
IBPtr ClusterFissioner::clone() const {
return new_ptr(*this);
}
IBPtr ClusterFissioner::fullclone() const {
return new_ptr(*this);
}
void ClusterFissioner::persistentOutput(PersistentOStream & os) const {
os << ounit(_clMaxLight,GeV)
<< ounit(_clMaxHeavy,GeV) << ounit(_clMaxExotic,GeV) << _clPowLight << _clPowHeavy
<< _clPowExotic << _pSplitLight
<< _pSplitHeavy << _pSplitExotic
<< _fissionCluster << _fissionPwt
<< _pwtDIquark
<< _diquarkClusterFission
<< ounit(_btClM,GeV)
<< _iopRem << ounit(_kappa, GeV/meter)
<< _enhanceSProb << ounit(_m0Fission,GeV) << _massMeasure
<< _dim << _phaseSpaceWeights
<< _hadronSpectrum << _kinematicThresholdChoice
<< _probPowFactor << _probShift << ounit(_kinThresholdShift,sqr(GeV))
<< _strictDiquarkKinematics
<< _covariantBoost
<< _allowHadronFinalStates
<< _massSampler
<< _phaseSpaceSampler
<< _matrixElement
<< _fissionApproach;
}
void ClusterFissioner::persistentInput(PersistentIStream & is, int) {
is >> iunit(_clMaxLight,GeV)
>> iunit(_clMaxHeavy,GeV) >> iunit(_clMaxExotic,GeV) >> _clPowLight >> _clPowHeavy
>> _clPowExotic >> _pSplitLight
>> _pSplitHeavy >> _pSplitExotic
>> _fissionCluster >> _fissionPwt
>> _pwtDIquark
>> _diquarkClusterFission
>> iunit(_btClM,GeV)
>> _iopRem >> iunit(_kappa, GeV/meter)
>> _enhanceSProb >> iunit(_m0Fission,GeV) >> _massMeasure
>> _dim >> _phaseSpaceWeights
>> _hadronSpectrum >> _kinematicThresholdChoice
>> _probPowFactor >> _probShift >> iunit(_kinThresholdShift,sqr(GeV))
>> _strictDiquarkKinematics
>> _covariantBoost
>> _allowHadronFinalStates
>> _massSampler
>> _phaseSpaceSampler
>> _matrixElement
>> _fissionApproach;
}
void ClusterFissioner::doinit() {
Interfaced::doinit();
for ( const long& id : spectrum()->heavyHadronizingQuarks() ) {
if ( _pSplitHeavy.find(id) == _pSplitHeavy.end() ||
_clPowHeavy.find(id) == _clPowHeavy.end() ||
_clMaxHeavy.find(id) == _clMaxHeavy.end() )
throw InitException() << "not all parameters have been set for heavy quark cluster fission";
}
// for default Pwts not needed to initialize
if (_fissionCluster==0) return;
for ( const long& id : spectrum()->lightHadronizingQuarks() ) {
if ( _fissionPwt.find(id) == _fissionPwt.end() )
// check that all relevant weights are set
throw InitException() << "fission weights for light quarks have not been set";
}
/*
// Not needed since we set Diquark weights from quark weights
for ( const long& id : spectrum()->lightHadronizingDiquarks() ) {
if ( _fissionPwt.find(id) == _fissionPwt.end() )
throw InitException() << "fission weights for light diquarks have not been set";
}*/
double pwtDquark=_fissionPwt.find(ParticleID::d)->second;
double pwtUquark=_fissionPwt.find(ParticleID::u)->second;
double pwtSquark=_fissionPwt.find(ParticleID::s)->second;
// ERROR: TODO makeDiquarkID is protected function ?
// _fissionPwt[spectrum()->makeDiquarkID(ParticleID::d,ParticleID::d,3)] = _pwtDIquark * pwtDquark * pwtDquark;
// _fissionPwt[spectrum()->makeDiquarkID(ParticleID::u,ParticleID::d,1)] = 0.5 * _pwtDIquark * pwtUquark * pwtDquark;
// _fissionPwt[spectrum()->makeDiquarkID(ParticleID::u,ParticleID::u,3)] = _pwtDIquark * pwtUquark * pwtUquark;
// _fissionPwt[spectrum()->makeDiquarkID(ParticleID::s,ParticleID::d,1)] = 0.5 * _pwtDIquark * pwtSquark * pwtDquark;
// _fissionPwt[spectrum()->makeDiquarkID(ParticleID::s,ParticleID::u,1)] = 0.5 * _pwtDIquark * pwtSquark * pwtUquark;
// _fissionPwt[spectrum()->makeDiquarkID(ParticleID::s,ParticleID::s,3)] = _pwtDIquark * pwtSquark * pwtSquark;
// TODO better solution for this magic number alternative
_fissionPwt[1103] = _pwtDIquark * pwtDquark * pwtDquark;
_fissionPwt[2101] = 0.5 * _pwtDIquark * pwtUquark * pwtDquark;
_fissionPwt[2203] = _pwtDIquark * pwtUquark * pwtUquark;
_fissionPwt[3101] = 0.5 * _pwtDIquark * pwtSquark * pwtDquark;
_fissionPwt[3201] = 0.5 * _pwtDIquark * pwtSquark * pwtUquark;
_fissionPwt[3303] = _pwtDIquark * pwtSquark * pwtSquark;
}
void ClusterFissioner::Init() {
static ClassDocumentation<ClusterFissioner> documentation
("Class responsibles for chopping up the clusters");
static Reference<ClusterFissioner,HadronSpectrum> interfaceHadronSpectrum
("HadronSpectrum",
"Set the Hadron spectrum for this cluster fissioner.",
&ClusterFissioner::_hadronSpectrum, 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 ParMap<ClusterFissioner,Energy> interfaceClMaxHeavy
("ClMaxHeavy",
"ClMax for heavy quarkls",
&ClusterFissioner::_clMaxHeavy, GeV, -1, 3.35*GeV, ZERO, 10.0*GeV,
false, false, Interface::upperlim);
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 ParMap<ClusterFissioner,double> interfaceClPowHeavy
("ClPowHeavy",
"ClPow for heavy quarks",
&ClusterFissioner::_clPowHeavy, -1, 1.0, 0.0, 10.0,
false, false, Interface::upperlim);
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 ParMap<ClusterFissioner,double> interfacePSplitHeavy
("PSplitHeavy",
"PSplit for heavy quarks",
&ClusterFissioner::_pSplitHeavy, -1, 1.0, 0.0, 10.0,
false, false, Interface::upperlim);
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, 0, false, false);
+ &ClusterFissioner::_fissionCluster, 1, false, false);
static SwitchOption interfaceFissionDefault
(interfaceFission,
"Default",
"Normal cluster fission which depends on the hadron spectrum class.",
0);
static SwitchOption interfaceFissionNew
(interfaceFission,
"New",
"Alternative cluster fission which does not depend on the hadron spectrum class",
1);
static Switch<ClusterFissioner,int> interfaceFissionApproach
("FissionApproach",
"Option for different Cluster Fission approaches",
&ClusterFissioner::_fissionApproach, 0, false, false);
static SwitchOption interfaceFissionApproachDefault
(interfaceFissionApproach,
"Default",
"Default Herwig-7.3.0 cluster fission without restructuring",
0);
static SwitchOption interfaceFissionApproachNew
(interfaceFissionApproach,
"New",
"New cluster fission which allows to choose MassSampler"
", PhaseSpaceSampler and MatrixElement",
1);
// Switch C->H1,C2 C->H1,H2 on or off
static Switch<ClusterFissioner,int> interfaceAllowHadronFinalStates
("AllowHadronFinalStates",
"Option for allowing hadron final states of cluster fission",
&ClusterFissioner::_allowHadronFinalStates, 0, false, false);
static SwitchOption interfaceAllowHadronFinalStatesNone
(interfaceAllowHadronFinalStates,
"None",
"Option for disabling hadron final states of cluster fission",
0);
static SwitchOption interfaceAllowHadronFinalStatesSemiHadronicOnly
(interfaceAllowHadronFinalStates,
"SemiHadronicOnly",
"Option for allowing hadron final states of cluster fission of type C->H1,C2 "
"(NOT YET USABLE WITH MatrixElement only use option None)",
1);
static SwitchOption interfaceAllowHadronFinalStatesAll
(interfaceAllowHadronFinalStates,
"All",
"Option for allowing hadron final states of cluster fission "
"of type C->H1,C2 or C->H1,H2 "
"(NOT YET USABLE WITH MatrixElement only use option None)",
2);
// Mass Sampler Switch
static Switch<ClusterFissioner,int> interfaceMassSampler
("MassSampler",
"Option for different Mass sampling options",
&ClusterFissioner::_massSampler, 0, false, false);
static SwitchOption interfaceMassSamplerDefault
(interfaceMassSampler,
"Default",
"Choose H7.2.3 default mass sampling using PSplitX",
0);
static SwitchOption interfaceMassSamplerUniform
(interfaceMassSampler,
"Uniform",
"Choose Uniform Mass sampling in M1,M2 space",
1);
static SwitchOption interfaceMassSamplerFlatPhaseSpace
(interfaceMassSampler,
"FlatPhaseSpace",
"Choose Flat Phase Space sampling of Mass in M1,M2 space",
2);
static SwitchOption interfaceMassSampleSoftMEPheno
(interfaceMassSampler,
"SoftMEPheno",
"Choose skewed Phase Space sampling of Masses in M1,M2 space "
"for greater efficiency in soft matrix element sampling",
3);
// Phase Space Sampler Switch
static Switch<ClusterFissioner,int> interfacePhaseSpaceSampler
("PhaseSpaceSampler",
"Option for different phase space sampling options",
&ClusterFissioner::_phaseSpaceSampler, 0, false, false);
static SwitchOption interfacePhaseSpaceSamplerDefault
(interfacePhaseSpaceSampler,
"FullyAligned",
"Herwig H7.2.3 default Cluster fission of all partons "
"aligned to the relative momentum of the mother cluster",
0);
static SwitchOption interfacePhaseSpaceSamplerAlignedIsotropic
(interfacePhaseSpaceSampler,
"AlignedIsotropic",
"Aligned Clusters but isotropic partons in their respective rest frame",
1);
static SwitchOption interfacePhaseSpaceSamplerFullyIsotropic
(interfacePhaseSpaceSampler,
"FullyIsotropic",
"Isotropic Clusters and isotropic partons in their respective rest frame "
"NOTE: Testing only!!",
2);
// Matrix Element Choice Switch
static Switch<ClusterFissioner,int> interfaceMatrixElement
("MatrixElement",
"Option for different Matrix Element options",
&ClusterFissioner::_matrixElement, 0, false, false);
static SwitchOption interfaceMatrixElementDefault
(interfaceMatrixElement,
"Default",
"Choose trivial matrix element i.e. whatever comes from the mass and "
"phase space sampling",
0);
static SwitchOption interfaceMatrixElementSoftQQbar
(interfaceMatrixElement,
"SoftQQbar",
"Choose Soft q1,q2->q1,q2,g*->q1,q2,q,qbar matrix element",
1);
static Switch<ClusterFissioner,int> interfaceDiquarkClusterFission
("DiquarkClusterFission",
"Allow clusters to fission to 1 or 2 diquark Clusters",
&ClusterFissioner::_diquarkClusterFission, 0, false, false);
static SwitchOption interfaceDiquarkClusterFissionAll
(interfaceDiquarkClusterFission,
"All",
"Allow diquark clusters and baryon clusters to fission to new diquark Clusters",
2);
static SwitchOption interfaceDiquarkClusterFissionOnlyBaryonClusters
(interfaceDiquarkClusterFission,
"OnlyBaryonClusters",
"Allow only baryon clusters to fission to new diquark Clusters",
1);
static SwitchOption interfaceDiquarkClusterFissionNo
(interfaceDiquarkClusterFission,
"No",
"Don't allow clusters to fission to new diquark Clusters",
0);
static ParMap<ClusterFissioner,double> interfaceFissionPwt
("FissionPwt",
"The weights for quarks in the fission process.",
&ClusterFissioner::_fissionPwt, -1, 1.0, 0.0, 10.0,
false, false, Interface::upperlim);
-
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,
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);
static Parameter<ClusterFissioner,double> interfaceProbPowFactor
("ProbablityPowerFactor",
"Power factor in ClausterFissioner bell probablity function",
&ClusterFissioner::_probPowFactor, 2.0, 1.0, 20.0,
false, false, Interface::limited);
static Parameter<ClusterFissioner,double> interfaceProbShift
("ProbablityShift",
"Shifts from the center in ClausterFissioner bell probablity function",
&ClusterFissioner::_probShift, 0.0, -10.0, 10.0,
false, false, Interface::limited);
static Parameter<ClusterFissioner,Energy2> interfaceKineticThresholdShift
("KineticThresholdShift",
"Shifts from the kinetic threshold in ClausterFissioner",
&ClusterFissioner::_kinThresholdShift, sqr(GeV), 0.*sqr(GeV), -10.0*sqr(GeV), 10.0*sqr(GeV),
false, false, Interface::limited);
static Switch<ClusterFissioner,int> interfaceKinematicThreshold
("KinematicThreshold",
"Option for using static or dynamic kinematic thresholds in cluster splittings",
&ClusterFissioner::_kinematicThresholdChoice, 0, false, false);
static SwitchOption interfaceKinematicThresholdStatic
(interfaceKinematicThreshold,
"Static",
"Set static kinematic thresholds for cluster splittings.",
0);
static SwitchOption interfaceKinematicThresholdDynamic
(interfaceKinematicThreshold,
"Dynamic",
"Set dynamic kinematic thresholds for cluster splittings.",
1);
static Switch<ClusterFissioner,bool> interfaceCovariantBoost
("CovariantBoost",
"Use single Covariant Boost for Cluster Fission",
&ClusterFissioner::_covariantBoost, false, false, false);
static SwitchOption interfaceCovariantBoostYes
(interfaceCovariantBoost,
"Yes",
"Use Covariant boost",
true);
static SwitchOption interfaceCovariantBoostNo
(interfaceCovariantBoost,
"No",
"Do NOT use Covariant boost",
false);
static Switch<ClusterFissioner,int> interfaceStrictDiquarkKinematics
("StrictDiquarkKinematics",
"Option for selecting different selection criterions of diquarks for ClusterFission",
&ClusterFissioner::_strictDiquarkKinematics, 0, false, false);
static SwitchOption interfaceStrictDiquarkKinematicsLoose
(interfaceStrictDiquarkKinematics,
"Loose",
"No kinematic threshold for diquark selection except for Mass bigger than 2 baryons",
0);
static SwitchOption interfaceStrictDiquarkKinematicsStrict
(interfaceStrictDiquarkKinematics,
"Strict",
"Resulting clusters are at least as heavy as 2 lightest baryons",
1);
static Parameter<ClusterFissioner,double> interfacePwtDIquark
("PwtDIquark",
"specific probability for choosing a d diquark",
&ClusterFissioner::_pwtDIquark, 0.0, 0.0, 10.0,
false, false, Interface::limited);
static Switch<ClusterFissioner,bool> interfacePhaseSpaceWeights
("PhaseSpaceWeights",
"Include phase space weights.",
&ClusterFissioner::_phaseSpaceWeights, false, false, false);
static SwitchOption interfacePhaseSpaceWeightsYes
(interfacePhaseSpaceWeights,
"Yes",
"Do include the effect of cluster fission phase space",
true);
static SwitchOption interfacePhaseSpaceWeightsNo
(interfacePhaseSpaceWeights,
"No",
"Do not include the effect of cluster phase space",
false);
static Parameter<ClusterFissioner,double>
interfaceDim ("Dimension","Dimension in which phase space weights are calculated",
&ClusterFissioner::_dim, 0, 4.0, 0.0, 10.0,false,false,false);
}
tPVector ClusterFissioner::fission(ClusterVector & clusters, bool softUEisOn) {
// return if no clusters
if (clusters.empty()) return tPVector();
/*****************
* Loop over the (input) collection of cluster pointers, and store in
* the vector splitClusters all the clusters that need to be split
* (these are beam clusters, if soft underlying event is off, and
* heavy non-beam clusters).
********************/
stack<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,
* 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
- // TODO maybe BeamClusters must not necessarily be split since can be very light
- // i.e. also lighter than the lightest constituents one can make of those
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,
bool softUEisOn) {
/**************************************************
* This method does the splitting of the cluster pointed by cluPtr
* and "recursively" by all of its cluster children, if heavy. All of these
* new children clusters are added (indeed the pointers to them) to the
* collection of cluster pointers collecCluPtr. The method works as follows.
* Initially the vector vecCluPtr contains just the input pointer to the
* cluster to be split. Then it will be filled "recursively" by all
* of the cluster's children that are heavy enough to require, in their turn,
* to be split. In each loop, the last element of the vector vecCluPtr is
* considered (only once because it is then removed from the vector).
* This approach is conceptually recursive, but avoid the overhead of
* a concrete recursive function. Furthermore it requires minimal changes
* in the case that the fission of an heavy cluster could produce more
* than two cluster children as assumed now.
*
* Draw the masses: for normal, non-beam clusters a power-like mass dist
* is used, whereas for beam clusters a fast-decreasing exponential mass
* dist is used instead (to avoid many iterative splitting which could
* produce an unphysical large transverse energy from a supposed soft beam
* remnant process).
****************************************/
// Here we recursively loop over clusters in the stack and cut them
while (!clusterStack.empty()) {
// take the last element of the vector
ClusterPtr iCluster = clusterStack.top(); clusterStack.pop();
// split it
cutType ct = iCluster->numComponents() == 2 ?
cutTwo(iCluster, finalhadrons, softUEisOn) :
cutThree(iCluster, finalhadrons, softUEisOn);
// There are cases when we don't want to split, even if it fails mass test
if(!ct.first.first || !ct.second.first) {
// if an unsplit beam cluster leave if for the underlying event
if(iCluster->isBeamCluster() && softUEisOn)
iCluster->isAvailable(false);
continue;
}
// check if clusters
ClusterPtr one = dynamic_ptr_cast<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' or C -> H + H' 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);
}
}
}
ClusterFissioner::cutType
ClusterFissioner::cutTwo(ClusterPtr & cluster, tPVector & finalhadrons,
bool softUEisOn) {
switch (_fissionApproach)
{
case 0:
return cutTwoDefault(cluster, finalhadrons, softUEisOn);
break;
case 1:
return cutTwoNew(cluster, finalhadrons, softUEisOn);
break;
default:
assert(false);
}
}
ClusterFissioner::cutType
ClusterFissioner::cutTwoDefault(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;
Lorentz5Momentum pClu1, pClu2, pQ1, pQone, pQtwo, pQ2;
do
{
succeeded = false;
++counter;
// get a flavour for the qqbar pair
drawNewFlavour(newPtr1,newPtr2,cluster);
// check for right ordering
assert (ptrQ2);
assert (newPtr2);
assert (ptrQ2->dataPtr());
assert (newPtr2->dataPtr());
if(cantMakeHadron(ptrQ1, newPtr1) || cantMakeHadron(ptrQ2, newPtr2)) {
swap(newPtr1, newPtr2);
// check again
if(cantMakeHadron(ptrQ1, newPtr1) || cantMakeHadron(ptrQ2, newPtr2)) {
throw Exception()
<< "ClusterFissioner cannot split the cluster ("
<< ptrQ1->PDGName() << ' ' << ptrQ2->PDGName()
<< ") into hadrons.\n" << Exception::runerror;
}
}
// Check that new clusters can produce particles and there is enough
// phase space to choose the drawn flavour
m1 = ptrQ1->data().constituentMass();
m2 = ptrQ2->data().constituentMass();
m = newPtr1->data().constituentMass();
// Do not split in the case there is no phase space available
if(Mc < m1+m + m2+m) continue;
pQ1.setMass(m1);
pQone.setMass(m);
pQtwo.setMass(m);
pQ2.setMass(m2);
- // pair<Energy,Energy> res = drawNewMasses(Mc, soft1, soft2, pClu1, pClu2,
- // ptrQ1, pQ1, newPtr1, pQone,
- // newPtr2, pQtwo, ptrQ2, pQ2);
- drawNewMasses(Mc, soft1, soft2, pClu1, pClu2,
- ptrQ1, pQ1, newPtr1, pQone,
- newPtr2, pQtwo, ptrQ2, pQ2);
+ // pair<Energy,Energy> res = drawNewMasses(Mc, soft1, soft2, pClu1, pClu2,
+ // ptrQ1, pQ1, newPtr1, pQone,
+ // newPtr2, pQtwo, ptrQ2, pQ2);
+ drawNewMasses(Mc, soft1, soft2, pClu1, pClu2, ptrQ1, pQ1, newPtr1, pQone, newPtr2, pQtwo, ptrQ2, pQ2);
// derive the masses of the children
Mc1 = pClu1.mass();
Mc2 = pClu2.mass();
// static kinematic threshold
if(_kinematicThresholdChoice == 0) {
if(Mc1 < m1+m || Mc2 < m+m2 || Mc1+Mc2 > Mc) continue;
// dynamic kinematic threshold
} else if(_kinematicThresholdChoice == 1) {
bool C1 = ( sqr(Mc1) )/( sqr(m1) + sqr(m) + _kinThresholdShift ) < 1.0 ? true : false;
bool C2 = ( sqr(Mc2) )/( sqr(m2) + sqr(m) + _kinThresholdShift ) < 1.0 ? true : false;
bool C3 = ( sqr(Mc1) + sqr(Mc2) )/( sqr(Mc) ) > 1.0 ? true : false;
if( C1 || C2 || C3 ) continue;
}
if ( _phaseSpaceWeights ) {
if ( phaseSpaceVeto(Mc,Mc1,Mc2,m,m1,m2) )
continue;
}
/**************************
* New (not present in Fortran Herwig):
* check whether the fragment masses Mc1 and Mc2 are above the
* threshold for the production of the lightest pair of hadrons with the
* right flavours. If not, then set by hand the mass to the lightest
* single hadron with the right flavours, in order to solve correctly
* the kinematics, and (later in this method) create directly such hadron
* and add it to the children hadrons of the cluster that undergoes the
* fission (i.e. the one pointed by iCluPtr). Notice that in this special
* case, the heavy cluster that undergoes the fission has one single
* cluster child and one single hadron child. We prefer this approach,
* rather than to create a light cluster, with the mass set equal to
* the lightest hadron, and let then the class LightClusterDecayer to do
* the job to decay it to that single hadron, for two reasons:
* First, because the sum of the masses of the two constituents can be,
* in this case, greater than the mass of that hadron, hence it would
* be impossible to solve the kinematics for such two components, and
* therefore we would have a cluster whose components are undefined.
* Second, the algorithm is faster, because it avoids the reshuffling
* procedure that would be necessary if we used LightClusterDecayer
* to decay the light cluster to the lightest hadron.
****************************/
// override chosen masses if needed
toHadron1 = _hadronSpectrum->chooseSingleHadron(ptrQ1->dataPtr(), newPtr1->dataPtr(),Mc1);
if(toHadron1) { Mc1 = toHadron1->mass(); pClu1.setMass(Mc1); }
toHadron2 = _hadronSpectrum->chooseSingleHadron(ptrQ2->dataPtr(), newPtr2->dataPtr(),Mc2);
if(toHadron2) { Mc2 = toHadron2->mass(); pClu2.setMass(Mc2); }
// if a beam cluster not allowed to decay to hadrons
if(cluster->isBeamCluster() && (toHadron1||toHadron2) && softUEisOn)
continue;
// Check if the decay kinematics is still possible: if not then
// force the one-hadron decay for the other cluster as well.
if(Mc1 + Mc2 > Mc) {
if(!toHadron1) {
toHadron1 = _hadronSpectrum->chooseSingleHadron(ptrQ1->dataPtr(), newPtr1->dataPtr(),Mc-Mc2);
if(toHadron1) { Mc1 = toHadron1->mass(); pClu1.setMass(Mc1); }
}
else if(!toHadron2) {
toHadron2 = _hadronSpectrum->chooseSingleHadron(ptrQ2->dataPtr(), newPtr2->dataPtr(),Mc-Mc1);
if(toHadron2) { Mc2 = toHadron2->mass(); pClu2.setMass(Mc2); }
}
}
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)
calculateKinematics(pClu,p0Q1,toHadron1,toHadron2,
pClu1,pClu2,pQ1,pQone,pQtwo,pQ2);
/******************
* The previous methods have determined the kinematics and positions
* of C -> C1 + C2.
* In the case that one of the two product is light, that means either
* decayOneHadronClu1 or decayOneHadronClu2 is true, then the momenta
* of the components of that light product have not been determined,
* and a (light) cluster will not be created: the heavy father cluster
* decays, in this case, into a single (not-light) cluster and a
* single hadron. In the other, "normal", cases the father cluster
* decays into two clusters, each of which has well defined components.
* Notice that, in the case of components which point to particles, the
* momenta of the components is properly set to the new values, whereas
* we do not change the momenta of the pointed particles, because we
* want to keep all of the information (that is the new momentum of a
* component after the splitting, which is contained in the _momentum
* member of the Component class, and the (old) momentum of that component
* before the splitting, which is contained in the momentum of the
* pointed particle). Please not make confusion of this only apparent
* inconsistency!
********************/
LorentzPoint posC,pos1,pos2;
posC = cluster->vertex();
calculatePositions(pClu, posC, pClu1, pClu2, pos1, pos2);
cutType rval;
if(toHadron1) {
rval.first = produceHadron(toHadron1, newPtr1, pClu1, pos1);
finalhadrons.push_back(rval.first.first);
}
else {
rval.first = produceCluster(ptrQ1, newPtr1, pClu1, pos1, pQ1, pQone, rem1);
}
if(toHadron2) {
rval.second = produceHadron(toHadron2, newPtr2, pClu2, pos2);
finalhadrons.push_back(rval.second.first);
}
else {
rval.second = produceCluster(ptrQ2, newPtr2, pClu2, pos2, pQ2, pQtwo, rem2);
}
return rval;
}
ClusterFissioner::cutType
ClusterFissioner::cutTwoNew(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();
if ( Mc < spectrum()->massLightestHadronPair(ptrQ1->dataPtr(),ptrQ2->dataPtr())) {
static const PPtr null = PPtr();
return cutType(PPair(null,null),PPair(null,null));
}
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_ME= 5000000;
int counter_MEtry = 0;
Energy Mc1 = ZERO;
Energy Mc2 = ZERO;
Energy m=ZERO;
Energy m1 = ptrQ1->data().constituentMass();
Energy m2 = ptrQ2->data().constituentMass();
Energy mMin = getParticleData(ParticleID::d)->constituentMass();
// TODO get rid of magic number
double eps = 0.1;
// Minimal threshold for non-zero Mass PhaseSpace
if ( Mc < (1.0+eps)*(m1 + m2 + 2*mMin )) {
static const PPtr null = PPtr();
return cutType(PPair(null,null),PPair(null,null));
}
tcPDPtr toHadron1, toHadron2;
PPtr newPtr1 = PPtr ();
PPtr newPtr2 = PPtr ();
Lorentz5Momentum pClu1, pClu2, pQ1, pQone, pQtwo, pQ2;
Lorentz5Momentum pClu = cluster->momentum(); // known
Lorentz5Momentum p0Q1 = ptrQ1->momentum(); // known (mom Q1 before fission)
Lorentz5Momentum p0Q2 = ptrQ2->momentum(); // known (mom Q2 before fission)
// where to return to in case of rejected sample
enum returnTo {
FlavourSampling,
MassSampling,
PhaseSpaceSampling,
MatrixElementSampling,
Done
};
// start with FlavourSampling
returnTo retTo=FlavourSampling;
int letFissionToXHadrons = _allowHadronFinalStates;
// if a beam cluster not allowed to decay to hadrons
if (cluster->isBeamCluster() && softUEisOn) letFissionToXHadrons = 0;
// ### Flavour, Mass, PhaseSpace and MatrixElement Sampling loop until accepted: ###
do
{
switch (retTo)
{
case FlavourSampling:
{
// ## Flavour sampling and kinematic constraints ##
drawNewFlavour(newPtr1,newPtr2,cluster);
// get a flavour for the qqbar pair
// check for right ordering
assert (ptrQ2);
assert (newPtr2);
assert (ptrQ2->dataPtr());
assert (newPtr2->dataPtr());
// careful if DiquarkClusters can exist
bool Q1diq = DiquarkMatcher::Check(ptrQ1->id());
bool Q2diq = DiquarkMatcher::Check(ptrQ2->id());
bool newQ1diq = DiquarkMatcher::Check(newPtr1->id());
bool newQ2diq = DiquarkMatcher::Check(newPtr2->id());
bool diqClu1 = Q1diq && newQ1diq;
bool diqClu2 = Q2diq && newQ2diq;
// DEBUG only:
// std::cout << "Make Clusters: ( " << ptrQ1->PDGName() << " " << newPtr1->PDGName() << " ), ( "
// << ptrQ2->PDGName() << " " << newPtr2->PDGName() << " )\n";
// check if Hadron formation is possible
if (!( diqClu1 || diqClu2 )
&& (cantMakeHadron(ptrQ1, newPtr1) || cantMakeHadron(ptrQ2, newPtr2))) {
swap(newPtr1, newPtr2);
// check again
if(cantMakeHadron(ptrQ1, newPtr1) || cantMakeHadron(ptrQ2, newPtr2)) {
throw Exception()
<< "ClusterFissioner cannot split the cluster ("
<< ptrQ1->PDGName() << ' ' << ptrQ2->PDGName()
<< ") into hadrons.\n"
<< "drawn Flavour: "<< newPtr1->PDGName()<<"\n"<< Exception::runerror;
}
}
else if ( diqClu1 || diqClu2 ){
bool swapped=false;
if ( !diqClu1 && cantMakeHadron(ptrQ1,newPtr1) ) {
swap(newPtr1, newPtr2);
swapped=true;
}
if ( !diqClu2 && cantMakeHadron(ptrQ2,newPtr2) ) {
assert(!swapped);
swap(newPtr1, newPtr2);
}
if ( diqClu2 && diqClu1 && ptrQ1->id()*newPtr1->id()>0 ) {
assert(!swapped);
swap(newPtr1, newPtr2);
}
if (!diqClu1) assert(!cantMakeHadron(ptrQ1,newPtr1));
if (!diqClu2) assert(!cantMakeHadron(ptrQ2,newPtr2));
}
// Check that new clusters can produce particles and there is enough
// phase space to choose the drawn flavour
m = newPtr1->data().constituentMass();
// Do not split in the case there is no phase space available + permille security
if(Mc < (m1 + m + m2 + m)) {
retTo = FlavourSampling;
// escape if no flavour phase space possibile without fission
if (fabs((m - getParticleData(ThePEG::ParticleID::d)->constituentMass())/GeV) < 1e-3) {
counter_MEtry = max_loop_ME;
retTo = Done;
}
continue;
}
pQ1.setMass(m1);
pQone.setMass(m);
pQtwo.setMass(m);
pQ2.setMass(m2);
Energy MLHP1 = spectrum()->hadronPairThreshold(ptrQ1->dataPtr(),newPtr1->dataPtr());
Energy MLHP2 = spectrum()->hadronPairThreshold(ptrQ2->dataPtr(),newPtr2->dataPtr());
Energy MLH1 = _hadronSpectrum->lightestHadron(ptrQ1->dataPtr(),newPtr1->dataPtr())->mass();
Energy MLH2 = _hadronSpectrum->lightestHadron(ptrQ2->dataPtr(),newPtr2->dataPtr())->mass();
bool canBeSingleHadron1 = (m1 + m) < MLHP1;
bool canBeSingleHadron2 = (m2 + m) < MLHP2;
Energy mThresh1 = (m1 + m);
Energy mThresh2 = (m2 + m);
if (canBeSingleHadron1) mThresh1 = MLHP1;
if (canBeSingleHadron2) mThresh2 = MLHP2;
/*
// TODO LEAVE this since in case of threshold problems
std::cout << "######################\n";
std::cout << "Mc = "<< Mc/GeV <<"\n";
std::cout << "m1 = "<< m1/GeV <<"\n";
std::cout << "m2 = "<< m2/GeV <<"\n";
std::cout << "m = "<< m/GeV <<"\n";
std::cout << "MLHP1 = "<< MLHP1/GeV <<"\n";
std::cout << "MLHP2 = "<< MLHP2/GeV <<"\n";
std::cout << "MLH1 = "<< MLH1/GeV <<"\n";
std::cout << "MLH2 = "<< MLH2/GeV <<"\n";
std::cout << "######################\n";
std::cout << "Mc-M1m= "<< (Mc-(m1+m))/GeV <<"\n";
std::cout << "Mc-M2m= "<< (Mc-(m2+m))/GeV <<"\n";
std::cout << "BEAM = "<< cluster->isBeamCluster() <<"\n";
*/
switch (letFissionToXHadrons)
{
case 0:
{
// Option None: only C->C1,C2 allowed
// check if mass phase space is non-zero
// resample or escape if only allowed mass phase space is for C->H1,H2 or C->H1,C2
if ( Mc < (mThresh1 + mThresh2)) {
// escape if not even the lightest flavour phase space is possibile
// TODO make this independent of explicit u/d quark
if (
fabs((m - getParticleData(ThePEG::ParticleID::u)->constituentMass())/GeV) < 1e-3
||
fabs((m - getParticleData(ThePEG::ParticleID::d)->constituentMass())/GeV) < 1e-3
) {
counter_MEtry = max_loop_ME;
retTo = Done;
continue;
}
else {
retTo = FlavourSampling;
continue;
}
}
break;
}
case 1:
{
// Option SemiHadronicOnly: C->H,C allowed
// NOTE: TODO implement matrix element for this
// resample or escape if only allowed mass phase space is for C->H1,H2
// First case is for ensuring the enough mass to be available and second one rejects disjoint mass regions
if ( ( (canBeSingleHadron1 && canBeSingleHadron2)
&& Mc < (mThresh1 + mThresh2) )
||
( (canBeSingleHadron1 || canBeSingleHadron2)
&& (canBeSingleHadron1 ? Mc-(m2+m) < MLH1:false || canBeSingleHadron2 ? Mc-(m1+m) < MLH2:false) )
){
// escape if not even the lightest flavour phase space is possibile
// TODO make this independent of explicit u/d quark
if (
fabs((m - getParticleData(ThePEG::ParticleID::u)->constituentMass())/GeV) < 1e-3
||
fabs((m - getParticleData(ThePEG::ParticleID::d)->constituentMass())/GeV) < 1e-3
) {
counter_MEtry = max_loop_ME;
retTo = Done;
continue;
}
else {
retTo = FlavourSampling;
continue;
}
}
break;
}
case 2:
{
// Option All: C->H,C and C->H,H allowed
// NOTE: TODO implement matrix element for this
// Mass Phase space for all option can always be found if cluster massive enough to go
// to the lightest 2 hadrons
if ( ( (canBeSingleHadron1 || canBeSingleHadron2)
&& (canBeSingleHadron1 ? Mc-(m2+m) < MLH1:false || canBeSingleHadron2 ? Mc-(m1+m) < MLH2:false) )
){
// escape if not even the lightest flavour phase space is possibile
// TODO make this independent of explicit u/d quark
if (
fabs((m - getParticleData(ThePEG::ParticleID::u)->constituentMass())/GeV) < 1e-3
||
fabs((m - getParticleData(ThePEG::ParticleID::d)->constituentMass())/GeV) < 1e-3
) {
counter_MEtry = max_loop_ME;
retTo = Done;
continue;
}
else {
retTo = FlavourSampling;
continue;
}
}
break;
}
default:
assert(false);
}
// Note: want to fallthrough (in C++17 one could uncomment
// the line below to show that this is intentional)
[[fallthrough]];
}
/*
* MassSampling choices:
* - Default (default)
* - Uniform
* - FlatPhaseSpace
* - SoftMEPheno
* */
case MassSampling:
{
bool failure = drawNewMasses(Mc, soft1, soft2, pClu1, pClu2,
ptrQ1, pQ1, newPtr1, pQone,
newPtr2, pQtwo, ptrQ2, pQ2);
// TODO IF C->C1,C2 (and not C->C,H or H1,H2) masses sampled and is in PhaseSpace must push through
// because otherwise no matrix element
if(failure) {
// TODO check which option is better
// retTo = FlavourSampling;
retTo = MassSampling;
continue;
}
// derive the masses of the children
Mc1 = pClu1.mass();
Mc2 = pClu2.mass();
// static kinematic threshold
if(_kinematicThresholdChoice == 0) {
if(Mc1 < m1+m || Mc2 < m+m2 || Mc1+Mc2 > Mc){
// TODO check which option is better
// retTo = FlavourSampling;
retTo = MassSampling;
continue;
}
// dynamic kinematic threshold
} else if(_kinematicThresholdChoice == 1) {
bool C1 = ( sqr(Mc1) )/( sqr(m1) + sqr(m) + _kinThresholdShift ) < 1.0 ? true : false;
bool C2 = ( sqr(Mc2) )/( sqr(m2) + sqr(m) + _kinThresholdShift ) < 1.0 ? true : false;
bool C3 = ( sqr(Mc1) + sqr(Mc2) )/( sqr(Mc) ) > 1.0 ? true : false;
if( C1 || C2 || C3 ) {
// TODO check which option is better
// retTo = FlavourSampling;
retTo = MassSampling;
continue;
}
}
/**************************
* New (not present in Fortran Herwig):
* check whether the fragment masses Mc1 and Mc2 are above the
* threshold for the production of the lightest pair of hadrons with the
* right flavours. If not, then set by hand the mass to the lightest
* single hadron with the right flavours, in order to solve correctly
* the kinematics, and (later in this method) create directly such hadron
* and add it to the children hadrons of the cluster that undergoes the
* fission (i.e. the one pointed by iCluPtr). Notice that in this special
* case, the heavy cluster that undergoes the fission has one single
* cluster child and one single hadron child. We prefer this approach,
* rather than to create a light cluster, with the mass set equal to
* the lightest hadron, and let then the class LightClusterDecayer to do
* the job to decay it to that single hadron, for two reasons:
* First, because the sum of the masses of the two constituents can be,
* in this case, greater than the mass of that hadron, hence it would
* be impossible to solve the kinematics for such two components, and
* therefore we would have a cluster whose components are undefined.
* Second, the algorithm is faster, because it avoids the reshuffling
* procedure that would be necessary if we used LightClusterDecayer
* to decay the light cluster to the lightest hadron.
****************************/
// override chosen masses if needed
toHadron1 = _hadronSpectrum->chooseSingleHadron(ptrQ1->dataPtr(), newPtr1->dataPtr(),Mc1);
if ( letFissionToXHadrons == 0 && toHadron1 ) {
// reject mass samples which would force C->C,H or C->H1,H2 fission if desired
// TODO check which option is better
// retTo = FlavourSampling;
retTo = MassSampling;
continue;
}
if(toHadron1) {
Mc1 = toHadron1->mass();
pClu1.setMass(Mc1);
}
toHadron2 = _hadronSpectrum->chooseSingleHadron(ptrQ2->dataPtr(), newPtr2->dataPtr(),Mc2);
if ( letFissionToXHadrons == 0 && toHadron2 ) {
// reject mass samples which would force C->C,H or C->H1,H2 fission if desired
// TODO check which option is better
// retTo = FlavourSampling;
retTo = MassSampling;
continue;
}
if(toHadron2) {
Mc2 = toHadron2->mass();
pClu2.setMass(Mc2);
}
if (letFissionToXHadrons == 1 && (toHadron1 && toHadron2) ) {
// reject mass samples which would force C->H1,H2 fission if desired
// TODO check which option is better
// retTo = FlavourSampling;
retTo = MassSampling;
continue;
}
// Check if the decay kinematics is still possible: if not then
// force the one-hadron decay for the other cluster as well.
if(Mc1 + Mc2 > Mc) {
// reject if we would need to create a C->H,H process if letFissionToXHadrons<2
if (letFissionToXHadrons < 2) {
// TODO check which option is better
// retTo = FlavourSampling;
retTo = MassSampling;
continue;
}
// TODO forbid other cluster!!!! to be also at hadron
if(!toHadron1) {
toHadron1 = _hadronSpectrum->chooseSingleHadron(ptrQ1->dataPtr(), newPtr1->dataPtr(),Mc-Mc2);
// toHadron1 = _hadronSpectrum->chooseSingleHadron(ptrQ1->dataPtr(), newPtr1->dataPtr(),ZERO);
if(toHadron1) {
Mc1 = toHadron1->mass();
pClu1.setMass(Mc1);
}
}
else if(!toHadron2) {
toHadron2 = _hadronSpectrum->chooseSingleHadron(ptrQ2->dataPtr(), newPtr2->dataPtr(),Mc-Mc1);
// toHadron2 = _hadronSpectrum->chooseSingleHadron(ptrQ2->dataPtr(), newPtr2->dataPtr(),ZERO);
if(toHadron2) {
Mc2 = toHadron2->mass();
pClu2.setMass(Mc2);
}
}
}
if (Mc <= Mc1+Mc2){
// escape if already lightest quark drawn
if (
fabs((m - getParticleData(ThePEG::ParticleID::u)->constituentMass())/GeV) < 1e-3
||
fabs((m - getParticleData(ThePEG::ParticleID::d)->constituentMass())/GeV) < 1e-3
) {
counter_MEtry = max_loop_ME;
retTo = Done;
}
else {
// Try again with lighter quark
retTo = FlavourSampling;
}
continue;
}
// Determined the (5-components) momenta (all in the LAB frame)
p0Q1.setMass(ptrQ1->mass()); // known (mom Q1 before fission)
p0Q1.rescaleEnergy(); // TODO check if needed
p0Q2.setMass(ptrQ2->mass()); // known (mom Q2 before fission)
p0Q2.rescaleEnergy();// TODO check if needed
pClu.rescaleMass();
// Note: want to fallthrough (in C++17 one could uncomment
// the line below to show that this is intentional)
[[fallthrough]];
}
/*
* PhaseSpaceSampling choices:
* - FullyAligned (default)
* - AlignedIsotropic
* - FullyIsotropic
* */
case PhaseSpaceSampling:
{
// ### Sample the Phase Space with respect to Matrix Element: ###
// TODO insert here PhaseSpace sampler
bool failure = drawKinematics(pClu,p0Q1,p0Q2,toHadron1,toHadron2,
pClu1,pClu2,pQ1,pQone,pQtwo,pQ2);
if(failure) {
// TODO check which option is better
retTo = MassSampling;
continue;
}
// Should be precise i.e. no rejection expected
// Note: want to fallthrough (in C++17 one could uncomment
// the line below to show that this is intentional)
[[fallthrough]];
}
/*
* MatrixElementSampling choices:
* - Default (default)
* - SoftQQbar
* */
case MatrixElementSampling:
{
counter_MEtry++;
// TODO maybe bridge this to work more neatly
// Ignore matrix element for C->C,H or C->H1,H2 fission
if (toHadron1 || toHadron2) {
retTo = Done;
break;
}
// TODO insert here partonic Matrix Element rejection
double SQME = calculateSQME(p0Q1,p0Q2,pQ1,pQone,pQ2,pQtwo);
double SQMEoverEstimate = calculateSQME_OverEstimate(p0Q1,p0Q2,pQ1,pQone,pQ2,pQtwo);
double Pacc = SQME/SQMEoverEstimate;
if (Pacc > 1.0 || std::isnan(Pacc) || std::isinf(Pacc)) throw Exception() << "Pacc = "<< Pacc << " > 1 in ClusterFissioner::cutTwo" << Exception::runerror;
if (UseRandom::rnd()<Pacc) {
if (0 && _matrixElement!=0) {
// std::cout << "\nAccept Pacc = "<<Pacc<<"\n";
// std::ofstream out("data_CluFis.dat", std::ios::app | std::ios::out);
// out << Pacc << "\t"
// << Mc/GeV << "\t"
// << pClu1.mass()/GeV << "\t"
// << pClu2.mass()/GeV << "\t"
// << pQone*pQtwo/GeV2 << "\t"
// << pQ1*pQ2/GeV2 << "\t"
// << pQ1*pQtwo/GeV2 << "\t"
// << pQ2*pQone/GeV2 << "\t"
// << m/GeV
// << "\n";
// out.close();
Energy MbinStart=2.0*GeV;
Energy MbinEnd=91.0*GeV;
Energy dMbin=1.0*GeV;
Energy MbinLow;
Energy MbinHig;
if ( fabs((m-0.325*GeV)/GeV)<1e-5
&& fabs((m1-0.325*GeV)/GeV)<1e-5
&& fabs((m2-0.325*GeV)/GeV)<1e-5) {
int cnt = 0;
for (Energy MbinIt=MbinStart; MbinIt < MbinEnd; MbinIt+=dMbin)
{
if (0){
std::cout << "\nFirst\n" << std::endl;
int ctr=0;
for (Energy MbinIt=MbinStart; MbinIt < MbinEnd; MbinIt+=dMbin)
{
std::string name = "WriteOut/data_CluFis_BinM_"+std::to_string(ctr)+".dat";
std::ofstream out2(name,std::ios::out);
out2 << "# Binned from "<<MbinIt/GeV << "\tto\t" << (MbinIt+dMbin)/GeV<<"\n";
out2 << "# m=m1=m2=0.325\n";
out2 << "# (1) Pacc\t"
<< "(2) Mc/GeV\t"
<< "(3) M1/GeV\t"
<< "(4) M2/GeV\t"
<< "(5) q.qbar/GeV2\t"
<< "(6) q1.q2/GeV2\t"
<< "(7) q1.qbar/GeV2\t"
<< "(8) q2.q/GeV2\t"
<< "(9) m/GeV\n";
out2.close();
ctr++;
}
// first=0;
}
MbinLow = MbinIt;
MbinHig = MbinLow+dMbin;
if (Mc>MbinLow && Mc<MbinHig) {
std::string name = "WriteOut/data_CluFis_BinM_"+std::to_string(cnt)+".dat";
std::ofstream out2(name, std::ios::app );
out2 << Pacc << "\t"
<< Mc/GeV << "\t"
<< pClu1.mass()/GeV << "\t"
<< pClu2.mass()/GeV << "\t"
<< pQone*pQtwo/GeV2 << "\t"
<< pQ1*pQ2/GeV2 << "\t"
<< pQ1*pQtwo/GeV2 << "\t"
<< pQ2*pQone/GeV2 << "\t"
<< m/GeV
<< "\n";
out2.close();
break;
}
// if (Mc<MbinIt) break;
cnt++;
}
}
}
retTo=Done;
break;
}
retTo = MassSampling;
continue;
}
default:
{
assert(false);
}
}
}
while (retTo!=Done && counter_MEtry < max_loop_ME);
// while (retTo!=Done);
if(counter_MEtry >= max_loop_ME) {
// happens if we get at too light cluster to begin with
// or if we get too massive clusters where the matrix element
// is very inefficiently sampled
// TODO exclude this by ensuring that there is always enough phase space for M>m1+m2+2*m and maybe other conditions
// std::cout << "ERROR: Max Looped\n";
static const PPtr null = PPtr();
return cutType(PPair(null,null),PPair(null,null));
}
// ==> full sample generated
/******************
* The previous methods have determined the kinematics and positions
* of C -> C1 + C2.
* In the case that one of the two product is light, that means either
* decayOneHadronClu1 or decayOneHadronClu2 is true, then the momenta
* of the components of that light product have not been determined,
* and a (light) cluster will not be created: the heavy father cluster
* decays, in this case, into a single (not-light) cluster and a
* single hadron. In the other, "normal", cases the father cluster
* decays into two clusters, each of which has well defined components.
* Notice that, in the case of components which point to particles, the
* momenta of the components is properly set to the new values, whereas
* we do not change the momenta of the pointed particles, because we
* want to keep all of the information (that is the new momentum of a
* component after the splitting, which is contained in the _momentum
* member of the Component class, and the (old) momentum of that component
* before the splitting, which is contained in the momentum of the
* pointed particle). Please not make confusion of this only apparent
* inconsistency!
********************/
LorentzPoint posC,pos1,pos2;
posC = cluster->vertex();
calculatePositions(pClu, posC, pClu1, pClu2, pos1, pos2);
cutType rval;
if(toHadron1) {
rval.first = produceHadron(toHadron1, newPtr1, pClu1, pos1);
finalhadrons.push_back(rval.first.first);
}
else {
rval.first = produceCluster(ptrQ1, newPtr1, pClu1, pos1, pQ1, pQone, rem1);
}
if(toHadron2) {
rval.second = produceHadron(toHadron2, newPtr2, pClu2, pos2);
finalhadrons.push_back(rval.second.first);
}
else {
rval.second = produceCluster(ptrQ2, newPtr2, pClu2, pos2, pQ2, pQtwo, rem2);
}
return rval;
}
ClusterFissioner::cutType
ClusterFissioner::cutThree(ClusterPtr & cluster, tPVector & finalhadrons,
bool softUEisOn) {
// need to make sure only 3-cpt clusters get here
assert(cluster->numComponents() == 3);
// extract quarks
tPPtr ptrQ[3] = {cluster->particle(0),cluster->particle(1),cluster->particle(2)};
assert( ptrQ[0] && ptrQ[1] && ptrQ[2] );
// find maximum mass pair
Energy mmax(ZERO);
Lorentz5Momentum pDiQuark;
int iq1(-1),iq2(-1);
Lorentz5Momentum psum;
for(int q1=0;q1<3;++q1) {
psum+= ptrQ[q1]->momentum();
for(int q2=q1+1;q2<3;++q2) {
Lorentz5Momentum ptest = ptrQ[q1]->momentum()+ptrQ[q2]->momentum();
ptest.rescaleMass();
Energy mass = ptest.m();
if(mass>mmax) {
mmax = mass;
pDiQuark = ptest;
iq1 = q1;
iq2 = q2;
}
}
}
// and the spectators
int iother(-1);
for(int ix=0;ix<3;++ix) if(ix!=iq1&&ix!=iq2) iother=ix;
assert(iq1>=0&&iq2>=0&&iother>=0);
// And check if those particles are from a beam remnant
bool rem1 = cluster->isBeamRemnant(iq1);
bool rem2 = cluster->isBeamRemnant(iq2);
// workout which distribution to use
bool soft1(false),soft2(false);
switch (_iopRem) {
case 0:
soft1 = rem1 || rem2;
soft2 = rem2 || rem1;
break;
case 1:
soft1 = rem1;
soft2 = rem2;
break;
}
// Initialization for the exponential ("soft") mass distribution.
static const int max_loop = 1000;
int counter = 0;
Energy Mc1 = ZERO, Mc2 = ZERO, m1=ZERO, m2=ZERO, m=ZERO;
tcPDPtr toHadron;
bool toDiQuark(false);
PPtr newPtr1 = PPtr(),newPtr2 = PPtr();
PDPtr diquark;
bool succeeded = false;
Lorentz5Momentum pClu1, pClu2, pQ1, pQone, pQtwo, pQ2;
do {
succeeded = false;
++counter;
// get a flavour for the qqbar pair
drawNewFlavour(newPtr1,newPtr2,cluster);
// 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) || !spectrum()->canMakeDiQuark(ptrQ[iq2], newPtr2)) {
swap(newPtr1,newPtr2);
}
// check again
if(cantMakeHadron(ptrQ[iq1], newPtr1) || !spectrum()->canMakeDiQuark(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;
pQ1.setMass(m1);
pQone.setMass(m);
pQtwo.setMass(m);
pQ2.setMass(m2);
bool failure = drawNewMasses(mmax, soft1, soft2, pClu1, pClu2,
ptrQ[iq1], pQ1, newPtr1, pQone,
newPtr2, pQtwo, ptrQ[iq1], pQ2);
if (failure) continue;
Mc1 = pClu1.mass();
Mc2 = pClu2.mass();
if(Mc1 < m1+m || Mc2 < m+m2 || Mc1+Mc2 > mmax) continue;
+ if ( _phaseSpaceWeights ) {
+ if ( phaseSpaceVeto(mmax,Mc1,Mc2,m,m1,m2) )
+ continue;
+ }
+
// check if need to force meson clster to hadron
toHadron = _hadronSpectrum->chooseSingleHadron(ptrQ[iq1]->dataPtr(), newPtr1->dataPtr(),Mc1);
if(toHadron) { Mc1 = toHadron->mass(); pClu1.setMass(Mc1); }
// check if need to force diquark cluster to be on-shell
toDiQuark = false;
diquark = spectrum()->makeDiquark(ptrQ[iq2]->dataPtr(), newPtr2->dataPtr());
if(Mc2 < diquark->constituentMass()) {
Mc2 = diquark->constituentMass(); pClu2.setMass(Mc2);
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 = _hadronSpectrum->chooseSingleHadron(ptrQ[iq1]->dataPtr(), newPtr1->dataPtr(),mmax-Mc2);
if(toHadron) { Mc1 = toHadron->mass(); pClu1.setMass(Mc1); }
}
else if(!toDiQuark) {
Mc2 = _hadronSpectrum->massLightestHadron(ptrQ[iq2]->dataPtr(), newPtr2->dataPtr()); pClu2.setMass(Mc2);
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();
- Lorentz5Momentum p0Q2 = ptrQ[iq2]->momentum(); //dummy
- p0Q1.setMass(ptrQ[iq1]->momentum().mass());
- p0Q2.setMass(ptrQ[iq2]->momentum().mass());
- drawKinematics(pDiQuark,p0Q1,p0Q2,toHadron,toDiQuark,
+ calculateKinematics(pDiQuark,p0Q1,toHadron,toDiQuark,
pClu1,pClu2,pQ1,pQone,pQtwo,pQ2);
// positions of the new clusters
LorentzPoint pos1,pos2;
Lorentz5Momentum pBaryon = pClu2+ptrQ[iother]->momentum();
calculatePositions(cluster->momentum(), cluster->vertex(), pClu1, pBaryon, pos1, pos2);
// first the mesonic cluster/meson
cutType rval;
if(toHadron) {
rval.first = produceHadron(toHadron, newPtr1, pClu1, pos1);
finalhadrons.push_back(rval.first.first);
}
else {
rval.first = produceCluster(ptrQ[iq1], newPtr1, pClu1, pos1, pQ1, pQone, rem1);
}
if(toDiQuark) {
rem2 |= cluster->isBeamRemnant(iother);
PPtr newDiQuark = diquark->produceParticle(pClu2);
rval.second = produceCluster(newDiQuark, ptrQ[iother], pBaryon, pos2, pClu2,
ptrQ[iother]->momentum(), rem2);
}
else {
rval.second = produceCluster(ptrQ[iq2], newPtr2, pBaryon, pos2, pQ2, pQtwo, rem2,
ptrQ[iother],cluster->isBeamRemnant(iother));
}
cluster->isAvailable(false);
return rval;
}
ClusterFissioner::PPair
ClusterFissioner::produceHadron(tcPDPtr hadron, tPPtr newPtr, const Lorentz5Momentum &a,
const LorentzPoint &b) const {
PPair rval;
if(hadron->coloured()) {
rval.first = (_hadronSpectrum->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;
}
/**
* Calculate a veto for the phase space weight
*/
bool ClusterFissioner::phaseSpaceVeto(const Energy Mc, const Energy Mc1, const Energy Mc2,
const Energy m, const Energy m1, const Energy m2) const {
double M_temp = Mc/GeV;
double M1_temp = Mc1/GeV;
double M2_temp = Mc2/GeV;
double m_temp = m/GeV;
double m1_temp = m1/GeV;
double m2_temp = m2/GeV;
double lam1 = sqrt(Kinematics::kaellen(M1_temp, m1_temp, m_temp));
double lam2 = sqrt(Kinematics::kaellen(M2_temp, m2_temp, m_temp));
double lam3 = sqrt(Kinematics::kaellen(M_temp, M1_temp, M2_temp));
double ratio;
double PSweight = pow(lam1*lam2*lam3,_dim-3.)*pow(M1_temp*M2_temp,2.-_dim);
// overestimate only possible for dim>=3.0
assert(_dim>=3.0);
double overEstimate = _dim>=4.0 ? pow(M_temp,4.*_dim-14.):pow(M_temp,2*(_dim-3.0))/pow((m1_temp+m_temp)*(m2_temp+m_temp),4.0-_dim);
ratio = PSweight/overEstimate;
// if (!(ratio>=0)) std::cout << "M "<<M_temp<<" M1 "<<M1_temp<<" M2 "<<M2_temp<<" m1 "<<m1_temp<<" m2 "<<m2_temp<<" m "<<m_temp<<"\n\n";
assert (ratio >= 0);
assert (ratio <= 1);
return (UseRandom::rnd()>ratio);
}
/**
* Calculate the masses and possibly kinematics of the cluster
* fission at hand; if claculateKineamtics is perfomring non-trivial
* steps kinematics calulcated here will be overriden. Currently resorts to the default
*/
bool ClusterFissioner::drawNewMasses(const Energy Mc, const bool soft1, const bool soft2,
Lorentz5Momentum& pClu1, Lorentz5Momentum& pClu2,
tcPPtr ptrQ1, const Lorentz5Momentum& pQ1,
tcPPtr, const Lorentz5Momentum& pQone,
tcPPtr, const Lorentz5Momentum& pQtwo,
tcPPtr ptrQ2, const Lorentz5Momentum& pQ2) const {
// TODO in all for better efficiency treat
// letFissionToXHadrons!=0 as seperate cases for reduced phasespace
// when mi+m < MassLightestHadronPair
switch (_massSampler)
{
case 0:
return drawNewMassesDefault(Mc, soft1, soft2, pClu1, pClu2, ptrQ1, pQ1, pQone, pQtwo, ptrQ2, pQ2);
break;
case 1:
return drawNewMassesUniform(Mc, soft1, soft2, pClu1, pClu2, pQ1, pQone, pQ2);
break;
case 2:
return drawNewMassesPhaseSpace(Mc, soft1, soft2, pClu1, pClu2, pQ1, pQone, pQ2);
break;
case 3:
return drawNewMassesPhaseSpaceExtended(Mc, soft1, soft2, pClu1, pClu2, pQ1, pQone, pQ2);
break;
default:
assert(false);
}
return true;// failure
}
/**
* Calculate the masses and possibly kinematics of the cluster
* fission at hand; if claculateKineamtics is perfomring non-trivial
* steps kinematics claulcated here will be overriden. Currentl;y resorts to the default
*/
bool ClusterFissioner::drawNewMassesDefault(const Energy Mc, const bool soft1, const bool soft2,
Lorentz5Momentum& pClu1, Lorentz5Momentum& pClu2,
tcPPtr ptrQ1, const Lorentz5Momentum& pQ1,
const Lorentz5Momentum& pQone,
const Lorentz5Momentum& pQtwo,
tcPPtr ptrQ2, const Lorentz5Momentum& pQ2) const {
// power for splitting
double exp1 = !spectrum()->isExotic(ptrQ1->dataPtr()) ? _pSplitLight : _pSplitExotic;
double exp2 = !spectrum()->isExotic(ptrQ2->dataPtr()) ? _pSplitLight : _pSplitExotic;
for ( const long& id : spectrum()->heavyHadronizingQuarks() ) {
assert(_pSplitHeavy.find(id) != _pSplitHeavy.end());
if ( spectrum()->hasHeavy(id,ptrQ1->dataPtr()) ) exp1 = _pSplitHeavy.find(id)->second;
if ( spectrum()->hasHeavy(id,ptrQ2->dataPtr()) ) exp2 = _pSplitHeavy.find(id)->second;
}
Energy M1 = drawChildMass(Mc,pQ1.mass(),pQ2.mass(),pQone.mass(),exp1,soft1);
Energy M2 = drawChildMass(Mc,pQ2.mass(),pQ1.mass(),pQtwo.mass(),exp2,soft2);
pClu1.setMass(M1);
pClu2.setMass(M2);
return false; // succeeds
}
/**
* Sample the masses for flat phase space
* */
bool ClusterFissioner::drawNewMassesUniform(const Energy Mc, const bool soft1, const bool soft2,
Lorentz5Momentum& pClu1, Lorentz5Momentum& pClu2,
const Lorentz5Momentum& pQ1,
const Lorentz5Momentum& pQ,
const Lorentz5Momentum& pQ2) const {
Energy M1,M2;
const Energy m1 = pQ1.mass();
const Energy m2 = pQ2.mass();
const Energy m = pQ.mass();
const Energy M1min = m1 + m;
const Energy M2min = m2 + m;
const Energy M1max = Mc - M2min;
const Energy M2max = Mc - M1min;
assert(M1max-M1min>ZERO);
assert(M2max-M2min>ZERO);
double r1;
double r2;
int counter = 0;
const int max_counter = 100;
while (counter < max_counter) {
r1 = UseRandom::rnd();
r2 = UseRandom::rnd();
M1 = (M1max-M1min)*r1 + M1min;
M2 = (M2max-M2min)*r2 + M2min;
counter++;
if ( Mc > M1 + M2) break;
}
if (counter==max_counter
|| Mc < M1 + M2
|| M1 <= M1min
|| M2 <= M2min ) return true; // failure
pClu1.setMass(M1);
pClu2.setMass(M2);
return false; // succeeds
}
/**
* Sample the masses for flat phase space
* */
bool ClusterFissioner::drawNewMassesPhaseSpace(const Energy Mc, const bool soft1, const bool soft2,
Lorentz5Momentum& pClu1, Lorentz5Momentum& pClu2,
const Lorentz5Momentum& pQ1,
const Lorentz5Momentum& pQ,
const Lorentz5Momentum& pQ2) const {
Energy M1,M2,MuS;
const Energy m1 = pQ1.mass();
const Energy m2 = pQ2.mass();
const Energy m = pQ.mass();
const Energy M1min = m1 + m;
const Energy M2min = m2 + m;
// const Energy M1max = Mc - M2min;
// const Energy M2max = Mc - M1min;
// assert(M1max-M1min>ZERO);
// assert(M2max-M2min>ZERO);
double r1;
double r2;
int counter = 0;
const int max_counter = 200;
// const Energy MuMminS = M1min + M2min;
// const Energy MuMminD = M1min - M2min;
const Energy MuMax = Mc - (M1min+M2min);
while (counter < max_counter) {
r1 = UseRandom::rnd();
r2 = UseRandom::rnd();
// TODO make this more efficient
// M1 = (M1max-M1min)*r1 + M1min;
// M2 = (M2max-M2min)*r2 + M2min;
MuS = MuMax * sqrt(r1);
// MD = 2*MS*r2 - MS + MuMminD;
M1 = M1min + MuS * r2;
M2 = M2min + MuS * (1.0 - r2);
counter++;
// Automatically satisfied
// if ( Mc <= M1 + M2) std::cout << "Mc "<< Mc/GeV << " M1 "<< M1/GeV <<" M2 " <<M2/GeV << std::endl;;
// if ( M1 <= M1min ) std::cout << "M1 "<< M1/GeV <<" M1min " <<M1min/GeV << std::endl;;
// if ( M2 <= M2min ) std::cout << "M2 "<< M2/GeV <<" M2min " <<M2min/GeV << std::endl;;
// assert( Mc > M1 + M2) ;
// assert( M1 > M1min ) ;
// assert( M2 > M2min ) ;
// if ( Mc <= M1 + M2) continue;
// if ( M1 <= M1min ) continue;
// if ( M2 <= M2min ) continue;
if (!phaseSpaceVeto(Mc,M1,M2,m,m1,m2) ) break; // For FlatPhaseSpace sampling vetoing
}
if (counter==max_counter) return true; // failure
pClu1.setMass(M1);
pClu2.setMass(M2);
return false; // succeeds
}
/**
* Sample the masses for flat phase space with modulation e.g. here
* FlatPhaseSpace*(M1*M2)**alpha
* */
bool ClusterFissioner::drawNewMassesPhaseSpaceExtended(const Energy Mc,
const bool soft1, const bool soft2,
Lorentz5Momentum& pClu1, Lorentz5Momentum& pClu2,
const Lorentz5Momentum& pQ1,
const Lorentz5Momentum& pQ,
const Lorentz5Momentum& pQ2) const {
Energy M1,M2;
const Energy m1 = pQ1.mass();
const Energy m2 = pQ2.mass();
const Energy m = pQ.mass();
const Energy M1min = m1 + m;
const Energy M2min = m2 + m;
const Energy M1max = Mc - M2min;
const Energy M2max = Mc - M1min;
assert(M1max-M1min>ZERO);
assert(M2max-M2min>ZERO);
double r1;
double r2;
int counter = 0;
const int max_counter = 200;
const int alpha = 2;
while (counter < max_counter) {
r1 = UseRandom::rnd();
r2 = UseRandom::rnd();
// TODO make this more efficient
// Uniform sampling
// M1 = (M1max-M1min)*r1 + M1min;
// M2 = (M2max-M2min)*r2 + M2min;
// Power sampling giving (M1*M2)**alpha
M1 = M1min * pow( 1.0 - r1 + r1 * pow(M1max/M1min,alpha+1.0) , 1.0/(alpha+1.0));
M2 = M2min * pow( 1.0 - r2 + r2 * pow(M2max/M2min,alpha+1.0) , 1.0/(alpha+1.0));
counter++;
if ( Mc <= M1 + M2) continue;
if ( M1 <= M1min ) continue;
if ( M2 <= M2min ) continue;
if (!phaseSpaceVeto(Mc,M1,M2,m,m1,m2) ) break; // For FlatPhaseSpace sampling vetoing
}
if (counter==max_counter
|| Mc < M1 + M2
|| M1 <= M1min
|| M2 <= M2min ) return true; // failure
pClu1.setMass(M1);
pClu2.setMass(M2);
return false; // succeeds
}
void ClusterFissioner::drawNewFlavourDiquarks(PPtr& newPtrPos,PPtr& newPtrNeg,
const ClusterPtr & clu) 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.
unsigned hasDiquarks=0;
assert(clu->numComponents()==2);
tcPDPtr pD1=clu->particle(0)->dataPtr();
tcPDPtr pD2=clu->particle(1)->dataPtr();
bool isDiq1=DiquarkMatcher::Check(pD1->id());
if (isDiq1) hasDiquarks++;
bool isDiq2=DiquarkMatcher::Check(pD2->id());
if (isDiq2) hasDiquarks++;
assert(hasDiquarks<=2);
Energy Mc=(clu->momentum().mass());
// if (fabs(clu->momentum().massError() )>1e-14) std::cout << "Mass inconsistency CF : " << std::scientific << clu->momentum().massError() <<"\n";
// Not allow yet Diquark Clusters
// if ( hasDiquarks>=1 || Mc < spectrum()->massLightestBaryonPair(pD1,pD2) )
// return drawNewFlavour(newPtrPos,newPtrNeg);
Energy minMass;
Selector<long> choice;
// int countQ=0;
// int countDiQ=0;
// adding quark-antiquark pairs to the selection list
for ( const long& id : spectrum()->lightHadronizingQuarks() ) {
// TODO uncommenting below gives sometimes 0 selection possibility,
// maybe need to be checked in the LightClusterDecayer and ColourReconnector
// if (Mc < spectrum()->massLightestHadronPair(pD1,pD2)) continue;
// countQ++;
if (_fissionCluster==0) choice.insert(_hadronSpectrum->pwtQuark(id),id);
else if (_fissionCluster==1) choice.insert(_fissionPwt.find(id)->second,id);
else assert(false);
}
// adding diquark-antidiquark pairs to the selection list
switch (hasDiquarks)
{
case 0:
for ( const long& id : spectrum()->lightHadronizingDiquarks() ) {
if (_strictDiquarkKinematics) {
tPDPtr cand = getParticleData(id);
minMass = spectrum()->massLightestHadron(pD2,cand)
+ spectrum()->massLightestHadron(cand,pD1);
}
else minMass = spectrum()->massLightestBaryonPair(pD1,pD2);
if (Mc < minMass) continue;
// countDiQ++;
if (_fissionCluster==0) choice.insert(_hadronSpectrum->pwtQuark(id),id);
else if (_fissionCluster==1) choice.insert(_fissionPwt.find(id)->second,id);
else assert(false);
}
break;
case 1:
if (_diquarkClusterFission<1) break;
for ( const long& id : spectrum()->lightHadronizingDiquarks() ) {
tPDPtr diq = getParticleData(id);
if (isDiq1)
minMass = spectrum()->massLightestHadron(pD2,diq)
+ spectrum()->massLightestBaryonPair(diq,pD1);
else
minMass = spectrum()->massLightestHadron(pD1,diq)
+ spectrum()->massLightestBaryonPair(diq,pD2);
if (Mc < minMass) continue;
// countDiQ++;
if (_fissionCluster==0) choice.insert(_hadronSpectrum->pwtQuark(id),id);
else if (_fissionCluster==1) choice.insert(_fissionPwt.find(id)->second,id);
else assert(false);
}
break;
case 2:
if (_diquarkClusterFission<2) break;
for ( const long& id : spectrum()->lightHadronizingDiquarks() ) {
tPDPtr diq = getParticleData(id);
if (Mc < spectrum()->massLightestBaryonPair(pD1,pD2)) {
throw Exception() << "Found Diquark Cluster:\n" << *clu << "\nwith MassCluster = "
<< ounit(Mc,GeV) <<" GeV MassLightestBaryonPair = "
<< ounit(spectrum()->massLightestBaryonPair(pD1,pD2) ,GeV)
<< " GeV cannot decay" << Exception::eventerror;
}
minMass = spectrum()->massLightestBaryonPair(pD1,diq)
+ spectrum()->massLightestBaryonPair(diq,pD2);
if (Mc < minMass) continue;
// countDiQ++;
if (_fissionCluster==0) choice.insert(_hadronSpectrum->pwtQuark(id),id);
else if (_fissionCluster==1) choice.insert(_fissionPwt.find(id)->second,id);
else assert(false);
}
break;
default:
assert(false);
}
assert(choice.size()>0);
long idNew = choice.select(UseRandom::rnd());
newPtrPos = getParticle(idNew);
newPtrNeg = getParticle(-idNew);
assert(newPtrPos);
assert(newPtrNeg);
assert(newPtrPos->dataPtr());
assert(newPtrNeg->dataPtr());
}
void ClusterFissioner::drawNewFlavourQuarks(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.
Selector<long> choice;
switch(_fissionCluster){
case 0:
for ( const long& id : spectrum()->lightHadronizingQuarks() )
choice.insert(_hadronSpectrum->pwtQuark(id),id);
break;
case 1:
for ( const long& id : spectrum()->lightHadronizingQuarks() )
choice.insert(_fissionPwt.find(id)->second,id);
break;
default :
assert(false);
}
long idNew = choice.select(UseRandom::rnd());
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 {
if ( spectrum()->gluonId() != ParticleID::g )
throw Exception() << "strange enhancement only working with Standard Model hadronization"
<< Exception::runerror;
// 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 = 0.;
double prob_u = 0.;
double prob_s = 0.;
double scale = abs(double(sqr(_m0Fission)/mass2));
// Choose which splitting weights you wish to use
switch(_fissionCluster){
// 0: ClusterFissioner and ClusterDecayer use the same weights
case 0:
prob_d = _hadronSpectrum->pwtQuark(ParticleID::d);
prob_u = _hadronSpectrum->pwtQuark(ParticleID::u);
/* Strangeness enhancement:
Case 1: probability scaling
Case 2: Exponential scaling
*/
if (_enhanceSProb == 1)
prob_s = (_maxScale < scale) ? 0. : pow(_hadronSpectrum->pwtQuark(ParticleID::s),scale);
else if (_enhanceSProb == 2)
prob_s = (_maxScale < scale) ? 0. : exp(-scale);
break;
/* 1: ClusterFissioner uses its own unique set of weights,
i.e. decoupled from ClusterDecayer */
case 1:
prob_d = _fissionPwt.find(ParticleID::d)->second;
prob_u = _fissionPwt.find(ParticleID::u)->second;
if (_enhanceSProb == 1)
prob_s = (_maxScale < scale) ? 0. : pow(_fissionPwt.find(ParticleID::s)->second,scale);
else if (_enhanceSProb == 2)
prob_s = (_maxScale < scale) ? 0. : exp(-scale);
break;
default:
assert(false);
}
int choice = UseRandom::rnd3(prob_u, prob_d, prob_s);
long idNew = 0;
switch (choice) {
case 0: idNew = ThePEG::ParticleID::u; break;
case 1: idNew = ThePEG::ParticleID::d; break;
case 2: idNew = ThePEG::ParticleID::s; break;
}
newPtrPos = getParticle(idNew);
newPtrNeg = getParticle(-idNew);
assert (newPtrPos);
assert(newPtrNeg);
assert (newPtrPos->dataPtr());
assert(newPtrNeg->dataPtr());
}
Energy2 ClusterFissioner::clustermass(const ClusterPtr & cluster) const {
Lorentz5Momentum pIn = cluster->momentum();
Energy2 endpointmass2 = sqr(cluster->particle(0)->mass() +
cluster->particle(1)->mass());
Energy2 singletm2 = pIn.m2();
// Return either the cluster mass, or the lambda measure
return (_massMeasure == 0) ? singletm2 : singletm2 - endpointmass2;
}
Energy ClusterFissioner::drawChildMass(const Energy M, const Energy m1,
const Energy m2, const Energy m,
const double expt, const bool soft) const {
/***************************
* This method, given in input the cluster mass Mclu of an heavy cluster C,
* made of consituents of masses m1 and m2, draws the masses Mclu1 and Mclu2
* of, respectively, the children cluster C1, made of constituent masses m1
* and m, and cluster C2, of mass Mclu2 and made of constituent masses m2
* and m. The mass is extracted from one of the two following mass
* distributions:
* --- power-like ("normal" distribution)
* d(Prob) / d(M^exponent) = const
* where the exponent can be different from the two children C1 (exp1)
* and C2 (exponent2).
* --- exponential ("soft" distribution)
* d(Prob) / d(M^2) = exp(-b*M)
* where b = 2.0 / average.
* Such distributions are limited below by the masses of
* the constituents quarks, and above from the mass of decaying cluster C.
* The choice of which of the two mass distributions to use for each of the
* two cluster children is dictated by iRemnant (see below).
* If the number of attempts to extract a pair of mass values that are
* kinematically acceptable is above some fixed number (max_loop, see below)
* the method gives up and returns false; otherwise, when it succeeds, it
* returns true.
*
* These distributions have been modified from HERWIG:
* Before these were:
* Mclu1 = m1 + (Mclu - m1 - m2)*pow( rnd(), 1.0/exponent1 );
* The new one coded here is a more efficient version, same density
* but taking into account 'in phase space from' beforehand
***************************/
// hard cluster
if(!soft) {
return pow(UseRandom::rnd(pow((M-m1-m2-m)*UnitRemoval::InvE, expt),
pow(m*UnitRemoval::InvE, expt)), 1./expt
)*UnitRemoval::E + m1;
}
// Otherwise it uses a soft mass distribution
else {
static const InvEnergy b = 2.0 / _btClM;
Energy max = M-m1-m2-2.0*m;
double rmin = b*max;
rmin = ( rmin < 50 ) ? exp(-rmin) : 0.;
double r1;
do {
r1 = UseRandom::rnd(rmin, 1.0) * UseRandom::rnd(rmin, 1.0);
}
while (r1 < rmin);
return m1 + m - log(r1)/b;
}
}
Axis ClusterFissioner::sampleDirectionAligned(const Lorentz5Momentum & pQ, const Lorentz5Momentum pClu) const {
Lorentz5Momentum pQinCOM(pQ);
pQinCOM.setMass(pQ.m());
pQinCOM.boost( -pClu.boostVector() ); // boost from LAB to C
return pQinCOM.vect().unit();
}
Axis ClusterFissioner::sampleDirectionIsotropic() const {
double cosTheta = -1 + 2.0 * UseRandom::rnd();
double sinTheta = sqrt(1.0-cosTheta*cosTheta);
double Phi = 2.0 * Constants::pi * UseRandom::rnd();
Axis uClusterUniform(cos(Phi)*cosTheta, sin(Phi)*cosTheta, sinTheta);
return uClusterUniform.unit();
}
Axis ClusterFissioner::sampleDirectionSemiUniform(const Lorentz5Momentum & pQ, const Lorentz5Momentum pClu) const {
Axis dir = sampleDirectionAligned(pQ,pClu);
Axis res;
do {
res=sampleDirectionIsotropic();
}
while (dir*res<0);
return res;
}
/* SQME for p1,p2->C1(q1,q),C2(q2,qbar)
* Note that:
* p0Q1 -> p1
* p0Q2 -> p2
* pQ1 -> q1
* pQone -> q
* pQ2 -> q2
* pQone -> qbar
* */
double ClusterFissioner::calculateSQME(
const Lorentz5Momentum & p1,
const Lorentz5Momentum & p2,
const Lorentz5Momentum & q1,
const Lorentz5Momentum & q,
const Lorentz5Momentum & q2,
const Lorentz5Momentum & qbar) const {
double SQME;
switch (_matrixElement)
{
case 0:
SQME = 1.0;
break;
case 1:
{
// Energy2 p1p2 = p1*p2;
Energy2 q1q2 = q1 * q2;
Energy2 q1q = q1 * q ;
Energy2 q2qbar = q2 * qbar;
Energy2 q2q = q2 * q ;
Energy2 q1qbar = q1 * qbar;
Energy2 qqbar = q * qbar;
Energy2 mq2 = q.mass2();
double Numerator = q1q2 * (qqbar + mq2)/sqr(GeV2);
Numerator += 0.5 * (q1q - q1qbar)*(q2q - q2qbar)/sqr(GeV2);
double Denominator = sqr(qqbar + mq2)*(q1q + q1qbar)*(q2q + q2qbar)/sqr(sqr(GeV2));
SQME = Numerator/Denominator;
break;
}
default:
assert(false);
}
if (SQME < 0) throw Exception()
<< "Squared Matrix Element = "<< SQME <<" < 0 in ClusterFissioner::calculateSQME() "
<< Exception::runerror;
return SQME;
}
/* Overestimate for SQME for p1,p2->C1(q1,q),C2(q2,qbar)
* Note that:
* p0Q1 -> p1
* p0Q2 -> p2
* pQ1 -> q1
* pQone -> q
* pQ2 -> q2
* pQone -> qbar
* */
double ClusterFissioner::calculateSQME_OverEstimate(
const Lorentz5Momentum & p1,
const Lorentz5Momentum & p2,
const Lorentz5Momentum & q1,
const Lorentz5Momentum & q,
const Lorentz5Momentum & q2,
const Lorentz5Momentum & qbar) const {
double SQME_OverEstimate;
switch (_matrixElement)
{
case 0:
SQME_OverEstimate = 1.0;
break;
case 1:
{
// Energy2 p1p2 = p1*p2;
Energy2 mq2 = q .mass2();
Energy2 m12 = q1.mass2();
Energy2 m22 = q2.mass2();
// Energy2 M2 = sqr(p1 + p2);
Energy2 M2 = sqr(q1 + q + q2 + qbar);
Energy2 M12 = sqr(q1 + q);
Energy2 M22 = sqr(q2 + qbar);
Energy2 q1q = (M12 - m12 - mq2)/2.0; // fixed by masses
Energy2 q2qbar = (M22 - m22 - mq2)/2.0; // fixed by masses
Energy2 q1q2Max = (M2 - M12 - M22)/2.0;
// Energy2 qqbarMax = (M2 - M12 - M22)/2.0;
Energy2 q2qMax = (M2 - M12 - M22)/2.0;
Energy2 q1qbarMax = (M2 - M12 - M22)/2.0;
// Energy2 q1q2Min = sqrt(m12*m22);
// Energy2 qqbarMin = mq2;
Energy2 q2qMin = sqrt(m22*mq2);
Energy2 q1qbarMin = sqrt(m12*mq2);
double Numerator = q1q2Max * (2.0*mq2)/sqr(GeV2);
Numerator += 0.5 * (q1q*q2qMax + q1qbarMax*q2qbar)/sqr(GeV2);
Numerator -= 0.5 * (q1q*q2qbar + q1qbarMin*q2qMin)/sqr(GeV2);
double Denominator = sqr(2.0*mq2)*(q1q + q1qbarMin)*(q2qMin + q2qbar)/sqr(sqr(GeV2));
SQME_OverEstimate = Numerator/Denominator;
if (SQME_OverEstimate < 0) throw Exception()
<< "Squared Matrix Element_Overestimate = "<< SQME_OverEstimate <<" < 0 in ClusterFissioner::calculateSQME_OverEstimate()\n"
<< "\tM = " << sqrt(M2)/GeV
<< "\tM1 = " << sqrt(M12)/GeV
<< "\tM2 = " << sqrt(M22)/GeV
<< "\n\tm1 = " << sqrt(m12)/GeV
<< "\tm2 = " << sqrt(m22)/GeV
<< "\tm = " << sqrt(mq2)/GeV
<< Exception::runerror;
break;
}
default:
assert(false);
}
return SQME_OverEstimate;
}
bool ClusterFissioner::drawKinematics(const Lorentz5Momentum & pClu,
const Lorentz5Momentum & p0Q1,
const Lorentz5Momentum & p0Q2,
const bool toHadron1,
const bool toHadron2,
Lorentz5Momentum & pClu1,
Lorentz5Momentum & pClu2,
Lorentz5Momentum & pQ1,
Lorentz5Momentum & pQbar,
Lorentz5Momentum & pQ,
Lorentz5Momentum & pQ2bar) const {
if (pClu.mass() < pClu1.mass() + pClu2.mass()
|| pClu1.mass()<ZERO
|| pClu2.mass()<ZERO ) {
throw Exception() << "Impossible Kinematics in ClusterFissioner::drawKinematics() (A)"
<< Exception::eventerror;
}
// Sample direction of the daughter clusters
Axis DirToClu = sampleDirectionCluster(p0Q1, pClu);
if (_covariantBoost) {
const Energy M = pClu.mass();
const Energy M1 = pClu1.mass();
const Energy M2 = pClu2.mass();
const Energy PcomClu=Kinematics::pstarTwoBodyDecay(M,M1,M2);
Momentum3 pClu1sampVect( PcomClu*DirToClu);
Momentum3 pClu2sampVect(-PcomClu*DirToClu);
pClu1.setMass(M1);
pClu1.setVect(pClu1sampVect);
pClu1.rescaleEnergy();
pClu2.setMass(M2);
pClu2.setVect(pClu2sampVect);
pClu2.rescaleEnergy();
}
else {
Kinematics::twoBodyDecay(pClu, pClu1.mass(), pClu2.mass(),DirToClu, pClu1, pClu2);
}
// In the case that cluster1 does not decay immediately into a single hadron,
// calculate the momenta of Q1 (as constituent of C1) and Qbar in the
// (parent) C1 frame first, where the direction of Q1 is u.vect().unit(),
// and then boost back in the LAB frame.
if(!toHadron1) {
if (pClu1.m() < pQ1.mass() + pQbar.mass() ) {
throw Exception() << "Impossible Kinematics in ClusterFissioner::drawKinematics() (B)"
<< Exception::eventerror;
}
Axis DirClu1;
switch (_phaseSpaceSampler)
{
case 0:
// Default aligned
DirClu1 = _covariantBoost ? pClu1.vect().unit():sampleDirectionAligned(pClu1, pClu);
break;
case 1:
// Isotropic
DirClu1 = sampleDirectionIsotropic();
break;
case 2:
// Isotropic
DirClu1 = sampleDirectionIsotropic();
break;
default:
assert(false);
}
// Need to boost constituents first into the pClu rest frame
// boost from Cluster1 rest frame to Cluster COM Frame
Kinematics::twoBodyDecay(pClu1, pQ1.mass(), pQbar.mass(), DirClu1, pQ1, pQbar);
}
// In the case that cluster2 does not decay immediately into a single hadron,
// Calculate the momenta of Q and Q2bar (as constituent of C2) in the
// (parent) C2 frame first, where the direction of Q is u.vect().unit(),
// and then boost back in the LAB frame.
if(!toHadron2) {
if (pClu2.m() < pQ.mass() + pQ2bar.mass() ) {
throw Exception() << "Impossible Kinematics in ClusterFissioner::drawKinematics() (C)"
<< Exception::eventerror;
}
Axis DirClu2;
switch (_phaseSpaceSampler)
{
case 0:
// Default aligned
DirClu2 = _covariantBoost ? pClu2.vect().unit():sampleDirectionAligned(pClu2, pClu);
break;
case 1:
// Isotropic
DirClu2 = sampleDirectionIsotropic();
break;
case 2:
// Isotropic
DirClu2 = sampleDirectionIsotropic();
break;
default:
assert(false);
}
// boost from Cluster2 rest frame to Cluster COM Frame
Kinematics::twoBodyDecay(pClu2, pQ.mass(), pQ2bar.mass(), DirClu2, pQ, pQ2bar);
}
// Boost all momenta from the Cluster COM frame to the Lab frame
if (_covariantBoost) {
std::vector<Lorentz5Momentum *> momenta;
momenta.push_back(&pClu1);
momenta.push_back(&pClu2);
if (!toHadron1) {
momenta.push_back(&pQ1);
momenta.push_back(&pQbar);
}
if (!toHadron2) {
momenta.push_back(&pQ);
momenta.push_back(&pQ2bar);
}
Kinematics::BoostIntoTwoParticleFrame(pClu.mass(),p0Q1, p0Q2, momenta);
}
return false; // success
}
void ClusterFissioner::calculateKinematics(const Lorentz5Momentum & pClu,
const Lorentz5Momentum & p0Q1,
const bool toHadron1,
const bool toHadron2,
Lorentz5Momentum & pClu1,
Lorentz5Momentum & pClu2,
Lorentz5Momentum & pQ1,
Lorentz5Momentum & pQbar,
Lorentz5Momentum & pQ,
Lorentz5Momentum & pQ2bar) const {
/******************
* This method solves the kinematics of the two body cluster decay:
* C (Q1 Q2bar) ---> C1 (Q1 Qbar) + C2 (Q Q2bar)
* In input we receive the momentum of C, pClu, and the momentum
* of the quark Q1 (constituent of C), p0Q1, both in the LAB frame.
* Furthermore, two boolean variables inform whether the two fission
* products (C1, C2) decay immediately into a single hadron (in which
* case the cluster itself is identify with that hadron) and we do
* not have to solve the kinematics of the components (Q1,Qbar) for
* C1 and (Q,Q2bar) for C2.
* The output is given by the following momenta (all 5-components,
* and all in the LAB frame):
* pClu1 , pClu2 respectively of C1 , C2
* pQ1 , pQbar respectively of Q1 , Qbar in C1
* pQ , pQ2bar respectively of Q , Q2 in C2
* The assumption, suggested from the string model, is that, in C frame,
* C1 and its constituents Q1 and Qbar are collinear, and collinear to
* the direction of Q1 in C (that is before cluster decay); similarly,
* (always in the C frame) C2 and its constituents Q and Q2bar are
* collinear (and therefore anti-collinear with C1,Q1,Qbar).
* The solution is then obtained by using Lorentz boosts, as follows.
* The kinematics of C1 and C2 is solved in their parent C frame,
* and then boosted back in the LAB. The kinematics of Q1 and Qbar
* is solved in their parent C1 frame and then boosted back in the LAB;
* similarly, the kinematics of Q and Q2bar is solved in their parent
* C2 frame and then boosted back in the LAB. In each of the three
* "two-body decay"-like cases, we use the fact that the direction
* of the motion of the decay products is known in the rest frame of
* their parent. This is obvious for the first case in which the
* parent rest frame is C; but it is also true in the other two cases
* where the rest frames are C1 and C2. This is because C1 and C2
* are boosted w.r.t. C in the same direction where their components,
* respectively (Q1,Qbar) and (Q,Q2bar) move in C1 and C2 rest frame
* respectively.
* Of course, although the notation used assumed that C = (Q1 Q2bar)
* where Q1 is a quark and Q2bar an antiquark, indeed everything remain
* unchanged also in all following cases:
* Q1 quark, Q2bar antiquark; --> Q quark;
* Q1 antiquark , Q2bar quark; --> Q antiquark;
* Q1 quark, Q2bar diquark; --> Q quark
* Q1 antiquark, Q2bar anti-diquark; --> Q antiquark
* Q1 diquark, Q2bar quark --> Q antiquark
* Q1 anti-diquark, Q2bar antiquark; --> Q quark
**************************/
// Calculate the unit three-vector, in the C frame, along which
// all of the constituents and children clusters move.
Lorentz5Momentum u(p0Q1);
u.boost( -pClu.boostVector() ); // boost from LAB to C
// the unit three-vector is then u.vect().unit()
// Calculate the momenta of C1 and C2 in the (parent) C frame first,
// where the direction of C1 is u.vect().unit(), and then boost back in the
// LAB frame.
if (pClu.m() < pClu1.mass() + pClu2.mass() ) {
throw Exception() << "Impossible Kinematics in ClusterFissioner::calculateKinematics() (A)"
<< Exception::eventerror;
}
Kinematics::twoBodyDecay(pClu, pClu1.mass(), pClu2.mass(),
u.vect().unit(), pClu1, pClu2);
// In the case that cluster1 does not decay immediately into a single hadron,
// calculate the momenta of Q1 (as constituent of C1) and Qbar in the
// (parent) C1 frame first, where the direction of Q1 is u.vect().unit(),
// and then boost back in the LAB frame.
if(!toHadron1) {
if (pClu1.m() < pQ1.mass() + pQbar.mass() ) {
throw Exception() << "Impossible Kinematics in ClusterFissioner::calculateKinematics() (B)"
<< Exception::eventerror;
}
Kinematics::twoBodyDecay(pClu1, pQ1.mass(), pQbar.mass(),
u.vect().unit(), pQ1, pQbar);
}
// In the case that cluster2 does not decay immediately into a single hadron,
// Calculate the momenta of Q and Q2bar (as constituent of C2) in the
// (parent) C2 frame first, where the direction of Q is u.vect().unit(),
// and then boost back in the LAB frame.
if(!toHadron2) {
if (pClu2.m() < pQ.mass() + pQ2bar.mass() ) {
throw Exception() << "Impossible Kinematics in ClusterFissioner::calculateKinematics() (C)"
<< Exception::eventerror;
}
Kinematics::twoBodyDecay(pClu2, pQ.mass(), pQ2bar.mass(),
u.vect().unit(), pQ, pQ2bar);
}
}
void ClusterFissioner::calculatePositions(const Lorentz5Momentum & pClu,
const LorentzPoint & positionClu,
const Lorentz5Momentum & pClu1,
const Lorentz5Momentum & pClu2,
LorentzPoint & positionClu1,
LorentzPoint & positionClu2) const {
// Determine positions of cluster children.
// See Marc Smith's thesis, page 127, formulas (4.122) and (4.123).
Energy Mclu = pClu.m();
Energy Mclu1 = pClu1.m();
Energy Mclu2 = pClu2.m();
// Calculate the unit three-vector, in the C frame, along which
// children clusters move.
Lorentz5Momentum u(pClu1);
u.boost( -pClu.boostVector() ); // boost from LAB to C frame
// the unit three-vector is then u.vect().unit()
Energy pstarChild = Kinematics::pstarTwoBodyDecay(Mclu,Mclu1,Mclu2);
// First, determine the relative positions of the children clusters
// in the parent cluster reference frame.
Energy2 mag2 = u.vect().mag2();
InvEnergy fact = mag2>ZERO ? 1./sqrt(mag2) : 1./GeV;
Length x1 = ( 0.25*Mclu + 0.5*( pstarChild + (sqr(Mclu2) - sqr(Mclu1))/(2.0*Mclu)))/_kappa;
Length t1 = Mclu/_kappa - x1;
LorentzDistance distanceClu1( x1 * fact * u.vect(), 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 * fact * u.vect(), 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::ProbablityFunction(double scale, double threshold) {
double cut = UseRandom::rnd(0.0,1.0);
return 1./(1.+pow(abs((threshold-_probShift)/scale),_probPowFactor)) > cut ? true : false;
}
bool ClusterFissioner::isHeavy(tcClusterPtr clu) {
// particle data for constituents
tcPDPtr cptr[3]={tcPDPtr(),tcPDPtr(),tcPDPtr()};
for(size_t ix=0;ix<min(clu->numComponents(),3);++ix) {
cptr[ix]=clu->particle(ix)->dataPtr();
}
// different parameters for exotic, bottom and charm clusters
double clpow = !spectrum()->isExotic(cptr[0],cptr[1],cptr[1]) ? _clPowLight : _clPowExotic;
Energy clmax = !spectrum()->isExotic(cptr[0],cptr[1],cptr[1]) ? _clMaxLight : _clMaxExotic;
for ( const long& id : spectrum()->heavyHadronizingQuarks() ) {
if ( spectrum()->hasHeavy(id,cptr[0],cptr[1],cptr[1]) ) {
clpow = _clPowHeavy[id];
clmax = _clMaxHeavy[id];
}
}
// required test for SUSY clusters, since aboveCutoff alone
// cannot guarantee (Mc > m1 + m2 + 2*m) in cut()
static const Energy minmass
= getParticleData(ParticleID::d)->constituentMass();
bool aboveCutoff = false, canSplitMinimally = false;
// static kinematic threshold
if(_kinematicThresholdChoice == 0) {
aboveCutoff = (
pow(clu->mass()*UnitRemoval::InvE , clpow)
>
pow(clmax*UnitRemoval::InvE, clpow)
+ pow(clu->sumConstituentMasses()*UnitRemoval::InvE, clpow)
);
canSplitMinimally = clu->mass() > clu->sumConstituentMasses() + 2.0 * minmass;
}
// dynamic kinematic threshold
else if(_kinematicThresholdChoice == 1) {
//some smooth probablity function to create a dynamic thershold
double scale = pow(clu->mass()/GeV , clpow);
double threshold = pow(clmax/GeV, clpow)
+ pow(clu->sumConstituentMasses()/GeV, clpow);
aboveCutoff = ProbablityFunction(scale,threshold);
scale = clu->mass()/GeV;
threshold = clu->sumConstituentMasses()/GeV + 2.0 * minmass/GeV;
canSplitMinimally = ProbablityFunction(scale,threshold);
}
return aboveCutoff && canSplitMinimally;
}
diff --git a/Hadronization/ClusterHadronizationHandler.cc b/Hadronization/ClusterHadronizationHandler.cc
--- a/Hadronization/ClusterHadronizationHandler.cc
+++ b/Hadronization/ClusterHadronizationHandler.cc
@@ -1,593 +1,590 @@
// -*- C++ -*-
//
// ClusterHadronizationHandler.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the 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/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 <ThePEG/Utilities/DescribeClass.h>
#include <ThePEG/Interface/Command.h>
using namespace Herwig;
ClusterHadronizationHandler * ClusterHadronizationHandler::currentHandler_ = 0;
DescribeClass<ClusterHadronizationHandler,HadronizationHandler>
describeClusterHadronizationHandler("Herwig::ClusterHadronizationHandler","Herwig.so");
IBPtr ClusterHadronizationHandler::clone() const {
return new_ptr(*this);
}
IBPtr ClusterHadronizationHandler::fullclone() const {
return new_ptr(*this);
}
void ClusterHadronizationHandler::persistentOutput(PersistentOStream & os)
const {
os << _partonSplitters << _clusterFinders << _colourReconnectors
<< _clusterFissioners << _lightClusterDecayers << _clusterDecayers
<< _gluonMassGenerators
<< _partonSplitter << _clusterFinder << _colourReconnector
<< _clusterFissioner << _lightClusterDecayer << _clusterDecayer
<< reshuffle_ << _gluonMassGenerator
<< ounit(_minVirtuality2,GeV2) << ounit(_maxDisplacement,mm)
<< _underlyingEventHandler << _reduceToTwoComponents;
}
void ClusterHadronizationHandler::persistentInput(PersistentIStream & is, int) {
is >> _partonSplitters >> _clusterFinders >> _colourReconnectors
>> _clusterFissioners >> _lightClusterDecayers >> _clusterDecayers
>> _gluonMassGenerators
>> _partonSplitter >> _clusterFinder >> _colourReconnector
>> _clusterFissioner >> _lightClusterDecayer >> _clusterDecayer
>> reshuffle_ >> _gluonMassGenerator
>> 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",
&Herwig::ClusterHadronizationHandler::_partonSplitter,
false, false, true, false);
static Reference<ClusterHadronizationHandler,GluonMassGenerator> interfaceGluonMassGenerator
("GluonMassGenerator",
"Set a reference to a gluon mass generator.",
&ClusterHadronizationHandler::_gluonMassGenerator, false, false, true, true, false);
static Switch<ClusterHadronizationHandler,int> interfaceReshuffle
("Reshuffle",
"Perform reshuffling if constituent masses have not yet been included by the shower",
&ClusterHadronizationHandler::reshuffle_, 0, false, false);
static SwitchOption interfaceReshuffleColourConnected
(interfaceReshuffle,
"ColourConnected",
"Separate reshuffling for colour connected partons.",
2);
static SwitchOption interfaceReshuffleGlobal
(interfaceReshuffle,
"Global",
"Reshuffle globally.",
1);
static SwitchOption interfaceReshuffleNo
(interfaceReshuffle,
"No",
"Do not reshuffle.",
0);
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",
&Herwig::ClusterHadronizationHandler::_colourReconnector,
false, false, true, false);
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",
&Herwig::ClusterHadronizationHandler::_lightClusterDecayer,
false, false, true, false);
static Reference<ClusterHadronizationHandler,ClusterDecayer>
interfaceClusterDecayer("ClusterDecayer",
"A reference to the ClusterDecayer object",
&Herwig::ClusterHadronizationHandler::_clusterDecayer,
false, false, true, false);
// functions to move the handlers so they are set
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
("MaxDisplacement",
"Maximum displacement that is allowed for a particle (unit [millimeter]).",
&ClusterHadronizationHandler::_maxDisplacement, mm, 1.0e-10*mm,
0.0*mm, 1.0e-9*mm,false,false,false);
static Reference<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);
static Command<ClusterHadronizationHandler> interfaceUseHandlersForInteraction
("UseHandlersForInteraction",
"Use the current set of handlers for the given interaction.",
&ClusterHadronizationHandler::_setHandlersForInteraction, false);
}
void ClusterHadronizationHandler::rebind(const TranslationMap & trans) {
for (auto iter : _partonSplitters) {
_partonSplitters[iter.first] = trans.translate(iter.second);
}
for (auto iter : _clusterFinders) {
_clusterFinders[iter.first] = trans.translate(iter.second);
}
for (auto iter : _colourReconnectors) {
_colourReconnectors[iter.first] = trans.translate(iter.second);
}
for (auto iter : _clusterFissioners) {
_clusterFissioners[iter.first] = trans.translate(iter.second);
}
for (auto iter : _lightClusterDecayers) {
_lightClusterDecayers[iter.first] = trans.translate(iter.second);
}
for (auto iter : _clusterDecayers) {
_clusterDecayers[iter.first] = trans.translate(iter.second);
}
for (auto iter : _gluonMassGenerators) {
_gluonMassGenerators[iter.first] = trans.translate(iter.second);
}
HandlerBase::rebind(trans);
}
IVector ClusterHadronizationHandler::getReferences() {
IVector ret = HandlerBase::getReferences();
for (auto iter : _partonSplitters) {
ret.push_back(iter.second);
}
for (auto iter : _clusterFinders) {
ret.push_back(iter.second);
}
for (auto iter : _colourReconnectors) {
ret.push_back(iter.second);
}
for (auto iter : _clusterFissioners) {
ret.push_back(iter.second);
}
for (auto iter : _lightClusterDecayers) {
ret.push_back(iter.second);
}
for (auto iter : _clusterDecayers) {
ret.push_back(iter.second);
}
for (auto iter : _gluonMassGenerators) {
ret.push_back(iter.second);
}
return ret;
}
void ClusterHadronizationHandler::doinit() {
HadronizationHandler::doinit();
// put current handlers into action for QCD to guarantee default behaviour
if ( _partonSplitters.empty() ) {
_partonSplitters[PDT::ColouredQCD] = _partonSplitter;
_clusterFinders[PDT::ColouredQCD] = _clusterFinder;
_colourReconnectors[PDT::ColouredQCD] = _colourReconnector;
_clusterFissioners[PDT::ColouredQCD] = _clusterFissioner;
_lightClusterDecayers[PDT::ColouredQCD] = _lightClusterDecayer;
_clusterDecayers[PDT::ColouredQCD] = _clusterDecayer;
_gluonMassGenerators[PDT::ColouredQCD] = _gluonMassGenerator;
}
}
string ClusterHadronizationHandler::_setHandlersForInteraction(string interaction) {
int interactionId = -1;
if ( interaction == " QCD" ) {
interactionId = PDT::ColouredQCD;
} else {
return "unknown interaction " + interaction;
}
_partonSplitters[interactionId] = _partonSplitter;
_clusterFinders[interactionId] = _clusterFinder;
_colourReconnectors[interactionId] = _colourReconnector;
_clusterFissioners[interactionId] = _clusterFissioner;
_lightClusterDecayers[interactionId] = _lightClusterDecayer;
_clusterDecayers[interactionId] = _clusterDecayer;
_gluonMassGenerators[interactionId] = _gluonMassGenerator;
return "";
}
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;
- for (auto iter : _clusterFinders) {
- iter.second->spectrum()->setGenerator(currentHandler_->generator());
- }
PVector theList(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
// TODO: Should this be hard-coded here, or defined in inputs?
map<int,int> interactions = {{PDT::ColouredQCD, ParticleID::g}};
// PVector currentlist(tagged.begin(),tagged.end());
map<int,PVector> currentlists;
for ( const auto & p : theList ) {
for ( auto i : interactions ) {
if ( p->data().colouredInteraction() == i.first ) {
currentlists[i.first].push_back(p);
Energy2 Q02 = 1.01*sqr(getParticleData(i.second)->constituentMass());
if(currentlists[i.first].back()->scale()<Q02) currentlists[i.first].back()->scale(Q02);
}
}
}
map<int,ClusterVector> clusters;
// ATTENTION this needs to have changes to include the new hadron spectrum functionality (gluon mass generator)
for ( auto i : interactions ) {
// reshuffle to the constituents
if ( reshuffle_ ) {
vector<PVector> reshufflelists;
if ( reshuffle_ == 1 ) { // global reshuffling
reshufflelists.push_back(currentlists[i.first]);
}
else if ( reshuffle_ == 2 ) {// colour connected reshuffling
splitIntoColourSinglets(currentlists[i.first], reshufflelists, i.first);
}
Ptr<GluonMassGenerator>::ptr gMassGenerator;
auto git = _gluonMassGenerators.find(i.first);
if ( git != _gluonMassGenerators.end() )
gMassGenerator = git->second;
for (auto currentlist : reshufflelists){
// get available energy and energy needed for constituent mass shells
LorentzMomentum totalQ;
Energy needQ = ZERO;
size_t nGluons = 0; // number of gluons for which a mass need be generated
for ( auto p : currentlist ) {
totalQ += p->momentum();
if ( p->id() == i.second && gMassGenerator ) {
++nGluons;
needQ += gMassGenerator->minGluonMass();
continue;
}
needQ += p->dataPtr()->constituentMass();
}
Energy Q = totalQ.m();
if ( needQ > Q )
throw Exception() << "cannot reshuffle to constituent mass shells" << Exception::eventerror;
// generate gluon masses if needed
list<Energy> gluonMasses;
if ( nGluons && gMassGenerator )
gluonMasses = gMassGenerator->generateMany(nGluons,Q-needQ);
// set masses for inidividual particles
vector<Energy> masses;
for ( auto p : currentlist ) {
if ( p->id() == i.second && gMassGenerator ) {
list<Energy>::const_iterator it = gluonMasses.begin();
advance(it,UseRandom::irnd(gluonMasses.size()));
masses.push_back(*it);
gluonMasses.erase(it);
}
else {
masses.push_back(p->dataPtr()->constituentMass());
}
}
// reshuffle to new masses
reshuffle(currentlist,masses);
}
}
// split the gluons
if ( !currentlists[i.first].empty() ) {
assert(_partonSplitters.find(i.first) != _partonSplitters.end());
_partonSplitters[i.first]->split(currentlists[i.first]);
}
// form the clusters
if ( !currentlists[i.first].empty() ) {
assert(_clusterFinders.find(i.first) != _clusterFinders.end());
clusters[i.first] = _clusterFinders[i.first]->formClusters(currentlists[i.first]);
// reduce BV clusters to two components now if needed
if(_reduceToTwoComponents)
_clusterFinders[i.first]->reduceToTwoComponents(clusters[i.first]);
}
}
// perform colour reconnection if needed and then
// decay the clusters into one hadron
for ( auto i : interactions ) {
if ( clusters[i.first].empty() )
continue;
bool lightOK = false;
short tried = 0;
const ClusterVector savedclusters = clusters[i.first];
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[i.first].size() );
BVclusters.reserve( clusters[i.first].size() );
for (size_t ic = 0; ic < clusters[i.first].size(); ++ic) {
ClusterPtr cl = clusters[i.first].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
assert(_colourReconnectors.find(i.first) != _colourReconnectors.end());
_colourReconnectors[i.first]->rearrange(CRclusters);
// tag new clusters as children of the partons to hadronize
_setChildren(CRclusters);
// forms diquarks
// we should already have used _clusterFinder[i.first]
_clusterFinders[i.first]->reduceToTwoComponents(CRclusters);
// recombine vectors of (possibly) reconnected and BV clusters
clusters[i.first].clear();
clusters[i.first].insert( clusters[i.first].end(), CRclusters.begin(), CRclusters.end() );
clusters[i.first].insert( clusters[i.first].end(), BVclusters.begin(), BVclusters.end() );
// fission of heavy clusters
// NB: during cluster fission, light hadrons might be produced straight away
assert(_clusterFissioners.find(i.first) != _clusterFissioners.end());
finalHadrons = _clusterFissioners[i.first]->fission(clusters[i.first],i.first == PDT::ColouredQCD ? isSoftUnderlyingEventON() : false);
// if clusters not previously reduced to two components do it now
if(!_reduceToTwoComponents)
_clusterFinders[i.first]->reduceToTwoComponents(clusters[i.first]);
assert(_lightClusterDecayers.find(i.first) != _lightClusterDecayers.end());
lightOK = _lightClusterDecayers[i.first]->decay(clusters[i.first],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) {
// std::cout << "failed LightClusterDecayer " << std::endl;
clusters[i.first] = savedclusters;
for_each(clusters[i.first].begin(),
clusters[i.first].end(),
std::mem_fn(&Particle::undecay));
}
}
if (!lightOK) {
throw Exception( "CluHad::handle(): tried LightClusterDecayer 10 times!",
Exception::eventerror);
}
// decay the remaining clusters
assert(_clusterDecayers[i.first]);
_clusterDecayers[i.first]->decay(clusters[i.first],finalHadrons);
}
// *****************************************
// *****************************************
// *****************************************
bool finalStateCluster=false;
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
// re-orders parents after their children
if ((*it)->children().empty()){
// If there is a cluster in the final state throw an event error
if((*it)->id()==81) {
finalStateCluster=true;
}
pstep->addDecayProduct(*it);
}
else {
pstep->addDecayProduct(*it);
pstep->addIntermediate(*it);
}
}
// For very small center of mass energies it might happen that baryonic clusters cannot decay into hadrons
if (finalStateCluster){
throw Exception( "CluHad::Handle(): Cluster in the final state",
Exception::eventerror);
}
// *****************************************
// *****************************************
// *****************************************
// soft underlying event if needed
if (isSoftUnderlyingEventON()) {
assert(_underlyingEventHandler);
ch.performStep(_underlyingEventHandler,Hint::Default());
}
}
// Sets parent child relationship of all clusters with two components
// Relationships for clusters with more than two components are set elsewhere in the Colour Reconnector
void ClusterHadronizationHandler::_setChildren(const ClusterVector & clusters) const {
// erase existing information about the partons' children
tPVector partons;
for ( const auto & cl : clusters ) {
if ( cl->numComponents() > 2 ) continue;
partons.push_back( cl->colParticle() );
partons.push_back( cl->antiColParticle() );
}
// erase all previous information about parent child relationship
for_each(partons.begin(), partons.end(), std::mem_fn(&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);
}
}
void ClusterHadronizationHandler::splitIntoColourSinglets(PVector copylist,
vector<PVector>& reshufflelists,
int){
PVector currentlist;
bool gluonloop;
PPtr firstparticle, temp;
reshufflelists.clear();
while (copylist.size()>0){
gluonloop=false;
currentlist.clear();
firstparticle=copylist.back();
copylist.pop_back();
if (!firstparticle->coloured()){
continue; //non-coloured particles are not included
}
currentlist.push_back(firstparticle);
//go up the anitColourLine and check if we are in a gluon loop
temp=firstparticle;
while( temp->hasAntiColour()){
temp = temp->antiColourLine()->endParticle();
if(temp==firstparticle){
gluonloop=true;
break;
}
else{
currentlist.push_back(temp);
copylist.erase(remove(copylist.begin(),copylist.end(), temp), copylist.end());
}
}
//if not a gluon loop, go up the ColourLine
if(!gluonloop){
temp=firstparticle;
while( temp->hasColour()){
temp=temp->colourLine()->startParticle();
currentlist.push_back(temp);
copylist.erase(remove(copylist.begin(),copylist.end(), temp), copylist.end());
}
}
reshufflelists.push_back(currentlist);
}
}
diff --git a/Hadronization/HadronSpectrum.cc b/Hadronization/HadronSpectrum.cc
--- a/Hadronization/HadronSpectrum.cc
+++ b/Hadronization/HadronSpectrum.cc
@@ -1,612 +1,613 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the HadronSpectrum class.
//
#include "HadronSpectrum.h"
#include "ClusterHadronizationHandler.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include <ThePEG/Repository/CurrentGenerator.h>
#include "Herwig/Utilities/Kinematics.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
namespace {
// debug helper
void dumpTable(const HadronSpectrum::HadronTable & tbl) {
typedef HadronSpectrum::HadronTable::const_iterator TableIter;
for (TableIter it = tbl.begin(); it != tbl.end(); ++it) {
cerr << it->first.first << ' '
<< it->first.second << '\n';
for (HadronSpectrum::KupcoData::const_iterator jt = it->second.begin();
jt != it->second.end(); ++jt) {
cerr << '\t' << *jt << '\n';
}
}
}
}
HadronSpectrum::HadronSpectrum()
: Interfaced(),
belowThreshold_(0),
_repwt(Lmax,vector<vector<double> >(Jmax,vector<double>(Nmax))) {}
HadronSpectrum::~HadronSpectrum() {}
void HadronSpectrum::doinit() {
Interfaced::doinit();
// construct the hadron tables
constructHadronTable();
// lightest members (hadrons)
for(const PDPtr & p1 : partons()) {
for(const PDPtr & p2 : partons()) {
tcPDPair lp = lightestHadronPair(p1,p2);
if(lp.first && lp.second)
lightestHadrons_[make_pair(p1->id(),p2->id())] = lp;
}
}
// for debugging
if (Debug::level >= 10)
dumpTable(table());
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void HadronSpectrum::persistentOutput(PersistentOStream & os) const {
os << _table << _partons << _forbidden
<< belowThreshold_ << _repwt << _pwt << lightestHadrons_;
}
void HadronSpectrum::persistentInput(PersistentIStream & is, int) {
is >> _table >> _partons >> _forbidden
>> belowThreshold_ >> _repwt >> _pwt >> lightestHadrons_;
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeAbstractClass<HadronSpectrum,Interfaced>
describeHerwigHadronSpectrum("Herwig::HadronSpectrum", "Herwig.so");
void HadronSpectrum::Init() {
static ClassDocumentation<HadronSpectrum> documentation
("There is no documentation for the HadronSpectrum class");
static RefVector<HadronSpectrum,ParticleData> interfacePartons
("Partons",
"The partons which are to be considered as the consistuents of the hadrons.",
&HadronSpectrum::_partons, -1, false, false, true, false, false);
static RefVector<HadronSpectrum,ParticleData> interfaceForbidden
("Forbidden",
"The PDG codes of the particles which cannot be produced in the hadronization.",
&HadronSpectrum::_forbidden, -1, false, false, true, false, false);
}
void HadronSpectrum::insertToHadronTable(tPDPtr &particle, int flav1, int flav2) {
// inserting a new Hadron in the hadron table.
long pid = particle->id();
int pspin = particle->iSpin();
HadronInfo a(pid, particle,specialWeight(pid),particle->mass());
// set the weight to the number of spin states
a.overallWeight = pspin*a.swtef;
// mesons
if(pspin%2==1) insertMeson(a,flav1,flav2);
// spin-1/2 baryons
else if(pspin==2) insertOneHalf(a,flav1,flav2);
// spin -3/2 baryons
else if(pspin==4) insertThreeHalf(a,flav1,flav2);
// all other cases
else {
assert(false);
}
}
void HadronSpectrum::insertOneHalf(HadronInfo a, int flav1, int flav2) {
assert(DiquarkMatcher::Check(flav1));
long iq1 = flav1/1000;
long iq2 = (flav1/100)%10;
if(iq1!=iq2 && flav1%10==3) flav1-=2;
if(iq1==iq2) {
if(iq1==flav2) {
a.overallWeight *= 1.5;
_table[make_pair(flav1,flav2)].insert(a);
_table[make_pair(flav2,flav1)].insert(a);
}
else {
_table[make_pair(flav1,flav2)].insert(a);
_table[make_pair(flav2,flav1)].insert(a);
long f3 = makeDiquarkID(iq1,flav2,1);
_table[make_pair(iq1,f3 )].insert(a);
_table[make_pair(f3 ,iq1)].insert(a);
}
}
else if(iq1==flav2) {
// ud1 u type
_table[make_pair(flav1,flav2)].insert(a);
_table[make_pair(flav2,flav1)].insert(a);
// and uu1 d type
long f3 = makeDiquarkID(iq1,iq1,3);
a.overallWeight *= a.wt;
_table[make_pair(f3 ,iq2)].insert(a);
_table[make_pair(iq2, f3)].insert(a);
}
else if(iq2==flav2) assert(false);
else {
_table[make_pair(flav1,flav2)].insert(a);
_table[make_pair(flav2,flav1)].insert(a);
long f3 = makeDiquarkID(iq1,flav2,1);
_table[make_pair(iq2,f3)].insert(a);
_table[make_pair(f3,iq2)].insert(a);
// 3rd perm
f3 = makeDiquarkID(iq2,flav2,1);
_table[make_pair(iq1,f3)].insert(a);
_table[make_pair(f3,iq1)].insert(a);
}
}
void HadronSpectrum::insertThreeHalf(HadronInfo a, int flav1, int flav2) {
assert(DiquarkMatcher::Check(flav1));
long iq1 = flav1/1000;
long iq2 = (flav1/100)%10;
if(iq1!=iq2 && flav1%10==3) flav1-=2;
if(iq1==iq2) {
if(iq1==flav2) {
a.overallWeight *= 1.5;
_table[make_pair(flav1,flav2)].insert(a);
_table[make_pair(flav2,flav1)].insert(a);
}
else {
_table[make_pair(flav1,flav2)].insert(a);
_table[make_pair(flav2,flav1)].insert(a);
long f3 = makeDiquarkID(iq1,flav2,1);
_table[make_pair(iq1,f3 )].insert(a);
_table[make_pair(f3 ,iq1)].insert(a);
}
}
else if(iq1==flav2) {
// ud1 u type
_table[make_pair(flav1,flav2)].insert(a);
_table[make_pair(flav2,flav1)].insert(a);
// and uu1 d type
long f3 = makeDiquarkID(iq1,iq1,3);
a.overallWeight *= a.wt;
_table[make_pair(f3 ,iq2)].insert(a);
_table[make_pair(iq2, f3)].insert(a);
}
else {
_table[make_pair(flav1,flav2)].insert(a);
_table[make_pair(flav2,flav1)].insert(a);
long f3 = makeDiquarkID(iq1,flav2,1);
_table[make_pair(iq2,f3)].insert(a);
_table[make_pair(f3,iq2)].insert(a);
// 3rd perm
f3 = makeDiquarkID(iq2,flav2,1);
_table[make_pair(iq1,f3)].insert(a);
_table[make_pair(f3,iq1)].insert(a);
}
}
tcPDPtr HadronSpectrum::chooseSingleHadron(tcPDPtr par1, tcPDPtr par2,
Energy mass) const {
Energy threshold = hadronPairThreshold(par1,par2);
// only do one hadron decay if mass less than the threshold
if(mass>=threshold) return tcPDPtr();
// select the hadron
tcPDPtr hadron;
// old option pick the lightest hadron
if(belowThreshold_ == 0) {
hadron= lightestHadron(par1,par2);
}
// new option select from those available
else if(belowThreshold_ == 1) {
vector<pair<tcPDPtr,double> > hadrons =
hadronsBelowThreshold(threshold,par1,par2);
if(hadrons.size()==1) {
hadron = hadrons[0].first;
}
else if(hadrons.empty()) {
hadron= lightestHadron(par1,par2);
}
else {
double totalWeight=0.;
for(unsigned int ix=0;ix<hadrons.size();++ix) {
totalWeight += hadrons[ix].second;
}
totalWeight *= UseRandom::rnd();
for(unsigned int ix=0;ix<hadrons.size();++ix) {
if(totalWeight<=hadrons[ix].second) {
hadron = hadrons[ix].first;
break;
}
else
totalWeight -= hadrons[ix].second;
}
assert(hadron);
}
}
else
assert(false);
return hadron;
}
tcPDPair HadronSpectrum::chooseHadronPair(const Energy cluMass,
tcPDPtr par1, tcPDPtr par2) const {
useMe();
// if either of the input partons is a diquark don't allow diquarks to be
// produced
bool isDiquark1 = DiquarkMatcher::Check(par1->id());
bool isDiquark2 = DiquarkMatcher::Check(par2->id());
bool noDiquarkInCluster = !(isDiquark1 || isDiquark2);
bool oneDiquarkInCluster = (isDiquark1 != isDiquark2);
bool quark = true;
// decide is baryon or meson production
if(noDiquarkInCluster) std::tie(quark,noDiquarkInCluster,oneDiquarkInCluster)
= selectBaryon(cluMass,par1,par2);
// 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 && std::find(hadronizingQuarks().begin(), hadronizingQuarks().end(),
abs(quarktopick->id())) != hadronizingQuarks().end()) continue;
if(DiquarkMatcher::Check(quarktopick->id()) &&
((!noDiquarkInCluster && quarktopick->iSpin()==1) ||
(!oneDiquarkInCluster && quarktopick->iSpin()==3))) continue;
HadronTable::const_iterator
tit1 = table().find(make_pair(abs(par1->id()),quarktopick->id()));
HadronTable::const_iterator
tit2 = table().find(make_pair(quarktopick->id(),abs(par2->id())));
// If not in table skip
if(tit1 == table().end()||tit2==table().end()) continue;
// tables empty skip
const KupcoData & T1 = tit1->second;
const KupcoData & T2 = tit2->second;
if(T1.empty()||T2.empty()) continue;
// if too massive skip
if(cluMass <= T1.begin()->mass +
T2.begin()->mass) continue;
// quark weight
double quarkWeight = pwt(quarktopick->id());
quarkWeight = specialQuarkWeight(quarkWeight,quarktopick->id(),
cluMass,par1,par2);
// loop over the hadrons
KupcoData::const_iterator H1,H2;
for(H1 = T1.begin();H1 != T1.end(); ++H1) {
for(H2 = T2.begin();H2 != T2.end(); ++H2) {
// break if cluster too light
if(cluMass < H1->mass + H2->mass) break;
weight = quarkWeight * H1->overallWeight * H2->overallWeight *
Kinematics::pstarTwoBodyDecay(cluMass, H1->mass, H2->mass);
int signQ = 0;
assert (par1 && quarktopick);
assert (par2);
assert(quarktopick->CC());
if(canBeHadron(par1, quarktopick->CC())
&& canBeHadron(quarktopick, par2))
signQ = +1;
else if(canBeHadron(par1, quarktopick)
&& canBeHadron(quarktopick->CC(), par2))
signQ = -1;
else {
cerr << "Could not make sign for" << par1->id()<< " " << quarktopick->id()
<< " " << par2->id() << "\n";
assert(false);
}
if (signQ == -1)
quarktopick = quarktopick->CC();
// construct the object with the info
Kupco a(quarktopick, H1->ptrData, H2->ptrData, weight);
hadrons.push_back(a);
wgtsum += weight;
}
}
}
if (hadrons.empty())
return make_pair(tcPDPtr(),tcPDPtr());
// select the hadron
wgtsum *= UseRandom::rnd();
unsigned int ix=0;
do {
wgtsum-= hadrons[ix].weight;
++ix;
}
while(wgtsum > ZERO && ix < hadrons.size());
if(ix == hadrons.size() && wgtsum > ZERO)
return make_pair(tcPDPtr(),tcPDPtr());
--ix;
assert(hadrons[ix].idQ);
int signHad1 = signHadron(par1, hadrons[ix].idQ->CC(), hadrons[ix].hadron1);
int signHad2 = signHadron(par2, hadrons[ix].idQ, hadrons[ix].hadron2);
assert( signHad1 != 0 && signHad2 != 0 );
return make_pair
( signHad1 > 0 ? hadrons[ix].hadron1 : tcPDPtr(hadrons[ix].hadron1->CC()),
signHad2 > 0 ? hadrons[ix].hadron2 : tcPDPtr(hadrons[ix].hadron2->CC()));
}
std::tuple<bool,bool,bool> HadronSpectrum::selectBaryon(const Energy, tcPDPtr, tcPDPtr ) const {
assert(false);
}
tcPDPair HadronSpectrum::lightestHadronPair(tcPDPtr ptr1, tcPDPtr ptr2) const {
Energy currentSum = Constants::MaxEnergy;
tcPDPair output;
for(unsigned int ix=0; ix<partons().size(); ++ix) {
HadronTable::const_iterator
tit1=table().find(make_pair(abs(ptr1->id()),partons()[ix]->id())),
tit2=table().find(make_pair(partons()[ix]->id(),abs(ptr2->id())));
if( tit1==table().end() || tit2==table().end()) continue;
if(tit1->second.empty()||tit2->second.empty()) continue;
Energy s = tit1->second.begin()->mass + tit2->second.begin()->mass;
if(currentSum > s) {
currentSum = s;
output.first = tit1->second.begin()->ptrData;
output.second = tit2->second.begin()->ptrData;
}
}
return output;
}
tcPDPtr HadronSpectrum::lightestHadron(tcPDPtr ptr1, tcPDPtr ptr2) const {
assert(ptr1 && ptr2);
// find entry in the table
pair<long,long> ids = make_pair(abs(ptr1->id()),abs(ptr2->id()));
HadronTable::const_iterator tit=_table.find(ids);
// throw exception if flavours wrong
if (tit==_table.end())
throw Exception() << "Could not find "
<< ids.first << ' ' << ids.second
<< " in _table. "
<< "In HadronSpectrum::lightestHadron()"
<< Exception::eventerror;
if(tit->second.empty())
throw Exception() << "HadronSpectrum::lightestHadron "
<< "could not find any hadrons containing "
<< ptr1->id() << ' ' << ptr2->id() << '\n'
<< tit->first.first << ' '
<< tit->first.second << Exception::eventerror;
// find the lightest hadron
int sign = signHadron(ptr1,ptr2,tit->second.begin()->ptrData);
tcPDPtr candidate = sign > 0 ?
tit->second.begin()->ptrData : tit->second.begin()->ptrData->CC();
// \todo 20 GeV limit is temporary fudge to let SM particles go through.
// \todo Use isExotic instead?
if (candidate->mass() > 20*GeV
&& candidate->mass() < ptr1->constituentMass() + ptr2->constituentMass()) {
generator()->log() << "HadronSpectrum::lightestHadron: "
<< "chosen candidate " << candidate->PDGName()
<< " is lighter than its constituents "
<< ptr1->PDGName() << ", " << ptr2->PDGName() << '\n'
<< candidate->mass()/GeV << " < " << ptr1->constituentMass()/GeV
<< " + " << ptr2->constituentMass()/GeV << '\n'
<< "Check your particle data tables.\n";
assert(false);
}
return candidate;
}
vector<pair<tcPDPtr,double> >
HadronSpectrum::hadronsBelowThreshold(Energy threshold, tcPDPtr ptr1,
tcPDPtr ptr2) const {
assert(ptr1 && ptr2);
// find entry in the table
pair<long,long> ids = make_pair(abs(ptr1->id()),abs(ptr2->id()));
HadronTable::const_iterator tit=_table.find(ids);
// throw exception if flavours wrong
if (tit==_table.end())
throw Exception() << "Could not find "
<< ids.first << ' ' << ids.second
<< " in _table. "
<< "In HadronSpectrum::hadronsBelowThreshold()"
<< Exception::eventerror;
if(tit->second.empty())
throw Exception() << "HadronSpectrum::hadronsBelowThreshold() "
<< "could not find any hadrons containing "
<< ptr1->id() << ' ' << ptr2->id() << '\n'
<< tit->first.first << ' '
<< tit->first.second << Exception::eventerror;
vector<pair<tcPDPtr,double> > candidates;
KupcoData::const_iterator hit = tit->second.begin();
// find the hadrons
while(hit!=tit->second.end()&&hit->mass<threshold) {
// find the hadron
int sign = signHadron(ptr1,ptr2,hit->ptrData);
tcPDPtr candidate = sign > 0 ? hit->ptrData : hit->ptrData->CC();
// \todo 20 GeV limit is temporary fudge to let SM particles go through.
// \todo Use isExotic instead?
if (candidate->mass() > 20*GeV
&& candidate->mass() < ptr1->constituentMass() + ptr2->constituentMass()) {
generator()->log() << "HadronSpectrum::hadronsBelowTheshold: "
<< "chosen candidate " << candidate->PDGName()
<< " is lighter than its constituents "
<< ptr1->PDGName() << ", " << ptr2->PDGName() << '\n'
<< candidate->mass()/GeV << " < " << ptr1->constituentMass()/GeV
<< " + " << ptr2->constituentMass()/GeV << '\n'
<< "Check your particle data tables.\n";
assert(false);
}
candidates.push_back(make_pair(candidate,hit->overallWeight));
++hit;
}
return candidates;
}
Energy HadronSpectrum::massLightestBaryonPair(tcPDPtr ptr1, tcPDPtr ptr2) const {
// Make sure that we don't have any diquarks as input, return arbitrarily
// large value if we do
Energy currentSum = Constants::MaxEnergy;
for(unsigned int ix=0; ix<_partons.size(); ++ix) {
if(!DiquarkMatcher::Check(_partons[ix]->id())) continue;
HadronTable::const_iterator
tit1=_table.find(make_pair(abs(ptr1->id()),_partons[ix]->id())),
tit2=_table.find(make_pair(_partons[ix]->id(),abs(ptr2->id())));
if( tit1==_table.end() || tit2==_table.end()) continue;
if(tit1->second.empty()||tit2->second.empty()) continue;
Energy s = tit1->second.begin()->mass + tit2->second.begin()->mass;
if(currentSum > s) currentSum = s;
}
return currentSum;
}
double HadronSpectrum::mesonWeight(long id) const {
// Total angular momentum
int j = ((id % 10) - 1) / 2;
// related to Orbital angular momentum l
int nl = (id/10000 )%10;
int l = -999;
int n = (id/100000)%10; // Radial excitation
if(j == 0) l = nl;
else if(nl == 0) l = j - 1;
else if(nl == 1 || nl == 2) l = j;
else if(nl == 3) l = j + 1;
// Angular or Radial excited meson
if((l||j||n) && l>=0 && l<Lmax && j<Jmax && n<Nmax) {
return sqr(_repwt[l][j][n]);
}
// rest is not excited or
// has spin >= 5/2 (ispin >= 6), haven't got those
else
return 1.0;
}
int HadronSpectrum::signHadron(tcPDPtr idQ1, tcPDPtr idQ2,
tcPDPtr hadron) const {
// This method receives in input three PDG ids, whose the
// first two have proper signs (corresponding to particles, id > 0,
// or antiparticles, id < 0 ), whereas the third one must
// be always positive (particle not antiparticle),
// corresponding to:
// --- quark-antiquark, or antiquark-quark, or
// quark-diquark, or diquark-quark, or
// antiquark-antidiquark, or antidiquark-antiquark
// for the first two input (idQ1, idQ2);
// --- meson or baryon for the third input (idHad):
// The method returns:
// --- + 1 if the two partons (idQ1, idQ2) are exactly
// the constituents for the hadron idHad;
// --- - 1 if the two partons (idQ1, idQ2) are exactly
// the constituents for the anti-hadron -idHad;
// --- + 0 otherwise.
// The method it is therefore useful to decide the
// sign of the id of the produced hadron as appeared
// in the vector _vecHad (where only hadron idHad > 0 are present)
// given the two constituent partons.
int sign = 0;
long idHad = hadron->id();
assert(idHad > 0);
int chargeIn = idQ1->iCharge() + idQ2->iCharge();
int chargeOut = hadron->iCharge();
// same charge
if( chargeIn == chargeOut && chargeIn !=0 ) sign = +1;
else if(chargeIn == -chargeOut && chargeIn !=0 ) sign = -1;
else if(chargeIn == 0 && chargeOut == 0 ) {
// In the case of same null charge, there are four cases:
// i) K0-like mesons, B0-like mesons, Bs-like mesons
// the PDG convention is to consider them "antiparticle" (idHad < 0)
// if the "dominant" (heavier) flavour (respectively, s, b)
// is a quark (idQ > 0): for instance, B0s = (b, sbar) has id < 0
// Remember that there is an important exception for K0L (id=130) and
// K0S (id=310): they don't have antiparticles, therefore idHad > 0
// always. We use below the fact that K0L and K0S are the unique
// hadrons having 0 the first (less significant) digit of their id.
// 2) D0-like mesons: the PDG convention is to consider them "particle"
// (idHad > 0) if the charm flavour is carried by a c: (c,ubar) has id>0
// 3) the remaining mesons should not have antiparticle, therefore their
// sign is always positive.
// 4) for baryons, that is when one of idQ1 and idQ2 is a (anti-) quark and
// the other one is a (anti-) diquark the sign is negative when both
// constituents are "anti", that is both with id < 0; positive otherwise.
// meson
if(std::find(hadronizingQuarks().begin(), hadronizingQuarks().end(),
abs(idQ1->id())) != hadronizingQuarks().end() &&
std::find(hadronizingQuarks().begin(), hadronizingQuarks().end(),
abs(idQ2->id())) != hadronizingQuarks().end())
{
int idQa = abs(idQ1->id());
int idQb = abs(idQ2->id());
int dominant = idQ2->id();
if(idQa > idQb) {
swap(idQa,idQb);
dominant = idQ1->id();
}
if((idQa==ParticleID::d && idQb==ParticleID::s) ||
(idQa==ParticleID::d && idQb==ParticleID::b) ||
(idQa==ParticleID::s && idQb==ParticleID::b)) {
// idHad%10 is zero for K0L,K0S
if (dominant < 0 || idHad%10 == 0) sign = +1;
else if(dominant > 0) sign = -1;
}
else if((idQa==ParticleID::u && idQb==ParticleID::c) ||
(idQa==ParticleID::u && idQb==ParticleID::t) ||
(idQa==ParticleID::c && idQb==ParticleID::t)) {
if (dominant > 0) sign = +1;
else if(dominant < 0) sign = -1;
}
else if(idQa==idQb) sign = +1;
// sets sign for Susy particles
else sign = (dominant > 0) ? +1 : -1;
}
// baryon
else if(DiquarkMatcher::Check(idQ1->id()) || DiquarkMatcher::Check(idQ2->id())) {
if (idQ1->id() > 0 && idQ2->id() > 0) sign = +1;
else if(idQ1->id() < 0 && idQ2->id() < 0) sign = -1;
}
}
if (sign == 0) {
cerr << "Could not work out sign for "
<< idQ1->PDGName() << ' '
<< idQ2->PDGName() << " => "
<< hadron->PDGName() << '\n';
assert(false);
}
return sign;
}
-/*
+
PDPtr HadronSpectrum::makeDiquark(tcPDPtr par1, tcPDPtr par2) const {
long id1 = par1->id();
long id2 = par2->id();
long pspin = id1==id2 ? 3 : 1;
long idnew = makeDiquarkID(id1,id2, pspin);
- return getParticleData(idnew);
+ assert(!CurrentGenerator::isVoid());
+ return CurrentGenerator::current().getParticleData(idnew);
}
-*/
+
bool HadronSpectrum::canBeMeson(tcPDPtr par1,tcPDPtr par2) const {
assert(par1 && par2);
long id1 = par1->id();
long id2 = par2->id();
// a Meson must not have any diquarks
if(DiquarkMatcher::Check(id1) || DiquarkMatcher::Check(id2)) return false;
return (std::find(hadronizingQuarks().begin(), hadronizingQuarks().end(),
abs(id1)) != hadronizingQuarks().end() &&
std::find(hadronizingQuarks().begin(), hadronizingQuarks().end(),
abs(id2)) != hadronizingQuarks().end() &&
id1*id2 < 0);
}
diff --git a/Hadronization/HadronSpectrum.h b/Hadronization/HadronSpectrum.h
--- a/Hadronization/HadronSpectrum.h
+++ b/Hadronization/HadronSpectrum.h
@@ -1,649 +1,649 @@
// -*- C++ -*-
#ifndef Herwig_HadronSpectrum_H
#define Herwig_HadronSpectrum_H
//
// This is the declaration of the HadronSpectrum class.
//
#include "ThePEG/Interface/Interfaced.h"
#include "HadronSpectrum.fh"
#include <ThePEG/Persistency/PersistentOStream.h>
#include <ThePEG/Persistency/PersistentIStream.h>
#include <ThePEG/PDT/ParticleData.h>
#include "Kupco.h"
/* These last two imports don't seem to be used here, but are needed for other
classes which import this. Should tidy up at some point*/
#include <ThePEG/PDT/StandardMatchers.h>
#include <ThePEG/Repository/EventGenerator.h>
namespace Herwig {
using namespace ThePEG;
/**
* Here is the documentation of the HadronSpectrum class.
*
* @see \ref HadronSpectrumInterfaces "The interfaces"
* defined for HadronSpectrum.
*/
class HadronSpectrum: public Interfaced {
public:
/** \ingroup Hadronization
* \class HadronInfo
* \brief Class used to store all the hadron information for easy access.
* \author Philip Stephens
*
* Note that:
* - the hadrons in _table can be filled in any ordered
* w.r.t. the mass value, and flavours for different
* groups (for instance, (u,s) hadrons don't need to
* be placed after (d,s) or any other flavour), but
* all hadrons with the same flavours must be consecutive
* ( for instance you cannot alternate hadrons of type
* (d,s) with those of flavour (u,s) ).
* Furthermore, it is assumed that particle and antiparticle
* have the same weights, and therefore only one of them
* must be entered in the table: we have chosen to refer
* to the particle, defined as PDG id > 0, although if
* an anti-particle is provided in input it is automatically
* transform to its particle, simply by taking the modulus
* of its id.
*/
class HadronInfo {
public:
/**
* Constructor
* @param idin The PDG code of the hadron
* @param datain The pointer to the ParticleData object
* @param swtin The singlet/decuplet/orbital factor
* @param massin The mass of the hadron
*/
HadronInfo(long idin=0, tPDPtr datain=tPDPtr(),
double swtin=1., Energy massin=ZERO)
: id(idin), ptrData(datain), swtef(swtin), wt(1.0), overallWeight(0.0),
mass(massin)
{}
/**
* Comparision operator on mass
*/
bool operator<(const HadronInfo &x) const {
if(mass!=x.mass) return mass < x.mass;
else return id < x.id;
}
/**
* The hadrons id.
*/
long id;
/**
* pointer to ParticleData, to get the spin, etc...
*/
tPDPtr ptrData;
/**
* singlet/decuplet/orbital factor
*/
double swtef;
/**
* mixing factor
*/
double wt;
/**
* (2*J+1)*wt*swtef
*/
double overallWeight;
/**
* The hadrons mass
*/
Energy mass;
/**
* Rescale the weight for a given hadron
*/
void rescale(double x) const {
const_cast<HadronInfo*>(this)->overallWeight *= x;
}
/**
* Friend method used to print the value of a table element.
*/
friend PersistentOStream & operator<< (PersistentOStream & os,
const HadronInfo & hi ) {
os << hi.id << hi.ptrData << hi.swtef << hi.wt
<< hi.overallWeight << ounit(hi.mass,GeV);
return os;
}
/**
* debug output
*/
friend ostream & operator<< (ostream & os, const HadronInfo & hi ) {
os << std::scientific << std::showpoint
<< std::setprecision(4)
<< setw(2)
<< hi.id << '\t'
<< hi.swtef << '\t'
<< hi.wt << '\t'
<< hi.overallWeight << '\t'
<< ounit(hi.mass,GeV);
return os;
}
/**
* Friend method used to read in the value of a table element.
*/
friend PersistentIStream & operator>> (PersistentIStream & is,
HadronInfo & hi ) {
is >> hi.id >> hi.ptrData >> hi.swtef >> hi.wt
>> hi.overallWeight >> iunit(hi.mass,GeV);
return is;
}
};
public:
/**
* The helper classes
*/
//@{
/**
* The type is used to contain all the hadrons info of a given flavour.
*/
typedef set<HadronInfo> KupcoData;
//@}
/**
* The hadron table type.
*/
typedef map<pair<long,long>,KupcoData> HadronTable;
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
HadronSpectrum();
/**
* The destructor.
*/
virtual ~HadronSpectrum();
//@}
public:
/** @name Partonic content */
//@{
/**
* Return the id of the gluon
*/
virtual long gluonId() const = 0;
/**
* Return the ids of all hadronizing quarks
*/
virtual const vector<long>& hadronizingQuarks() const = 0;
/**
* The light hadronizing quarks
*/
virtual const vector<long>& lightHadronizingQuarks() const = 0;
/**
* The light hadronizing diquarks
*/
virtual const vector<long>& lightHadronizingDiquarks() const = 0;
/**
* The heavy hadronizing quarks
*/
virtual const vector<long>& heavyHadronizingQuarks() const = 0;
/**
* Return true if any of the possible three input particles contains
* the indicated heavy quark. false otherwise. In the case that
* only the first particle is specified, it can be: an (anti-)quark,
* an (anti-)diquark an (anti-)meson, an (anti-)baryon; in the other
* cases, each pointer is assumed to be either (anti-)quark or
* (anti-)diquark.
*/
virtual bool hasHeavy(long id, tcPDPtr par1, tcPDPtr par2 = PDPtr(), tcPDPtr par3 = PDPtr()) const = 0;
/**
* Return true, if any of the possible input particle pointer is an
* exotic quark, e.g. Susy quark; false otherwise.
*/
virtual bool isExotic(tcPDPtr par1, tcPDPtr par2 = PDPtr(), tcPDPtr par3 = PDPtr()) const = 0;
//@}
/**
* Access the parton weights
*/
double pwt(long pid) const {
map<long,double>::const_iterator it = _pwt.find(abs(pid));
assert( it != _pwt.end() );
return it->second;
}
/**
* Return true if the two or three particles in input can be the components
* of a hadron; false otherwise.
*/
inline bool canBeHadron(tcPDPtr par1, tcPDPtr par2 , tcPDPtr par3 = PDPtr()) const {
return (!par3 && canBeMeson(par1,par2)) || canBeBaryon(par1,par2,par3);
}
/**
* Check if can't make a diquark from the partons
*/
bool canMakeDiQuark(tcPPtr p1, tcPPtr p2) const {
long id1 = p1->id(), id2 = p2->id();
return QuarkMatcher::Check(id1) && QuarkMatcher::Check(id2) && id1*id2>0;
}
/**
* Return the particle data of the diquark (anti-diquark) made by the two
* quarks (antiquarks) par1, par2.
* @param par1 (anti-)quark data pointer
* @param par2 (anti-)quark data pointer
*/
- virtual PDPtr makeDiquark(tcPDPtr par1, tcPDPtr par2) const = 0;
+ PDPtr makeDiquark(tcPDPtr par1, tcPDPtr par2) const;
/**
* Method to return a pair of hadrons given the PDG codes of
* two or three constituents
* @param cluMass The mass of the cluster
* @param par1 The first constituent
* @param par2 The second constituent
* @param par3 The third constituent
*/
virtual pair<tcPDPtr,tcPDPtr> chooseHadronPair(const Energy cluMass, tcPDPtr par1,
tcPDPtr par2) const;
/**
* Select the single hadron for a cluster decay
* return null pointer if not a single hadron decay
* @param par1 1st constituent
* @param par2 2nd constituent
* @param mass Mass of the cluster
*/
virtual tcPDPtr chooseSingleHadron(tcPDPtr par1, tcPDPtr par2, Energy mass) const;
/**
* This returns the lightest pair of hadrons given by the flavours.
*
* Given the two (or three) constituents of a cluster, it returns
* the two lightest hadrons with proper flavour numbers.
* Furthermore, the first of the two hadrons must have the constituent with
* par1, and the second must have the constituent with par2.
* \todo At the moment it does *nothing* in the case that also par3 is present.
*
* The method is implemented by calling twice lightestHadron,
* once with (par1,quarktopick->CC()) ,and once with (par2,quarktopick)
* where quarktopick is either the pointer to
* d or u quarks . In fact, the idea is that whatever the flavour of par1
* and par2, no matter if (anti-)quark or (anti-)diquark, the lightest
* pair of hadrons containing flavour par1 and par2 will have either
* flavour d or u, being the lightest quarks.
* The method returns the pair (PDPtr(),PDPtr()) if anything goes wrong.
*
* \todo The method assumes par3 == PDPtr() (otherwise we don't know how to proceed: a
* possible, trivial way would be to randomly select two of the three
* (anti-)quarks and treat them as a (anti-)diquark, reducing the problem
* to two components as treated below.
* In the normal (two components) situation, the strategy is the following:
* treat in the same way the two possibilities: (d dbar) (i=0) and
* (u ubar) (i=1) as the pair quark-antiquark necessary to form a
* pair of hadrons containing the input flavour par1 and par2; finally,
* select the one that produces the lightest pair of hadrons, compatible
* with the charge conservation constraint.
*/
tcPDPair lightestHadronPair(tcPDPtr ptr1, tcPDPtr ptr2) const;
/**
* Returns the mass of the lightest pair of hadrons with the given particles
* @param ptr1 is the first constituent
* @param ptr2 is the second constituent
*/
Energy massLightestHadronPair(tcPDPtr ptr1, tcPDPtr ptr2) const {
map<pair<long,long>,tcPDPair>::const_iterator lightest =
lightestHadrons_.find(make_pair(abs(ptr1->id()),abs(ptr2->id())));
if(lightest!=lightestHadrons_.end())
return lightest->second.first->mass()+lightest->second.second->mass();
else
return ZERO;
}
/**
* Returns the lightest hadron formed by the given particles.
*
* Given the id of two (or three) constituents of a cluster, it returns
* the lightest hadron with proper flavour numbers.
* @param ptr1 is the first constituent
* @param ptr2 is the second constituent
*/
tcPDPtr lightestHadron(tcPDPtr ptr1, tcPDPtr ptr2) const;
/**
* Return the threshold for a cluster to split into a pair of hadrons.
* This is normally the mass of the lightest hadron Pair, but can be
* higher for heavy and exotic clusters
*/
virtual Energy hadronPairThreshold(tcPDPtr par1, tcPDPtr par2) const=0;
/**
* Returns the hadrons below the constituent mass threshold formed by the given particles,
* together with their total weight
*
* Given the id of two (or three) constituents of a cluster, it returns
* the lightest hadron with proper flavour numbers.
* At the moment it does *nothing* in the case that also 'ptr3' present.
* @param threshold The theshold
* @param ptr1 is the first constituent
* @param ptr2 is the second constituent
* @param ptr3 is the third constituent
*/
vector<pair<tcPDPtr,double> > hadronsBelowThreshold(Energy threshold,
tcPDPtr ptr1, tcPDPtr ptr2) const;
/**
* Return the nominal mass of the hadron returned by lightestHadron()
* @param ptr1 is the first constituent
* @param ptr2 is the second constituent
* @param ptr3 is the third constituent
*/
Energy massLightestHadron(tcPDPtr ptr1, tcPDPtr ptr2) const {
// find entry in the table
pair<long,long> ids(abs(ptr1->id()),abs(ptr2->id()));
HadronTable::const_iterator tit=_table.find(ids);
// throw exception if flavours wrong
if(tit==_table.end()||tit->second.empty())
throw Exception() << "HadronSpectrum::massLightestHadron "
<< "failed for particle" << ptr1->id() << " "
<< ptr2->id()
<< Exception::eventerror;
// return the mass
return tit->second.begin()->mass;
}
/**
* Force baryon/meson selection
*/
virtual std::tuple<bool,bool,bool> selectBaryon(const Energy cluMass, tcPDPtr par1, tcPDPtr par2) const;
/**
* Returns the mass of the lightest pair of baryons.
* @param ptr1 is the first constituent
* @param ptr2 is the second constituent
*/
Energy massLightestBaryonPair(tcPDPtr ptr1, tcPDPtr ptr2) const;
/**
* Return the weight for the given flavour
*/
virtual double pwtQuark(const long& id) const = 0;
virtual double specialQuarkWeight(double quarkWeight, long,
const Energy, tcPDPtr, tcPDPtr) const {
return quarkWeight;
}
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();
void setGenerator(tEGPtr generator) {
Interfaced::setGenerator(generator);
}
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
*
* The array _repwt is initialized using the interfaces to set different
* weights for different meson multiplets and the constructHadronTable()
* method called to complete the construction of the hadron tables.
*
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
//@}
/**
* Return the id of the diquark (anti-diquark) made by the two
* quarks (antiquarks) of id specified in input (id1, id2).
* Caller must ensure that id1 and id2 are quarks.
*/
virtual long makeDiquarkID(long id1, long id2, long pspin) const = 0;
protected:
/**
* Return true if the two particles in input can be the components of a meson;
*false otherwise.
*/
bool canBeMeson(tcPDPtr par1,tcPDPtr par2) const;
/**
* Return true if the two or three particles in input can be the components
* of a baryon; false otherwise.
*/
virtual bool canBeBaryon(tcPDPtr par1, tcPDPtr par2 , tcPDPtr par3 = PDPtr()) const = 0;
/**
* A sub-function of HadronSpectrum::constructHadronTable().
* It receives the information of a prospective Hadron and inserts it
* into the hadron table construct.
* @param particle is a particle data pointer to the hadron
* @param flav1 is the first constituent of the hadron
* @param flav2 is the second constituent of the hadron
*/
void insertToHadronTable(tPDPtr &particle, int flav1, int flav2);
/**
* Construct the table of hadron data
* This is the main method to initialize the hadron data (mainly the
* weights associated to each hadron, taking into account its spin,
* eventual isoscalar-octect mixing, singlet-decuplet factor). This is
* the method that one should update when new or updated hadron data is
* available.
*
* This class implements the construction of the basic table but can be
* overridden if needed in inheriting classes.
*
* The rationale for factors used for diquarks involving different quarks can
* be can be explained by taking a prototype example that in the exact SU(2) limit,
* in which:
* \f[m_u=m_d\f]
* \f[M_p=M_n=M_\Delta\f]
* and we will have equal numbers of u and d quarks produced.
* Suppose that we weight 1 the diquarks made of the same
* quark and 1/2 those made of different quarks, the fractions
* of u and d baryons (p, n, Delta) we get are the following:
* - \f$\Delta^{++}\f$: 1 possibility only u uu with weight 1
* - \f$\Delta^- \f$: 1 possibility only d dd with weight 1
* - \f$p,\Delta^+ \f$: 2 possibilities u ud with weight 1/2
* d uu with weight 1
* - \f$n,\Delta^0 \f$: 2 possibilities d ud with weight 1/2
* u dd with weight 1
* In the latter two cases, we have to take into account the
* fact that p and n have spin 1/2 whereas Delta+ and Delta0
* have spin 3/2 therefore from phase space we get a double weight
* for Delta+ and Delta0 relative to p and n respectively.
* Therefore the relative amount of these baryons that is
* produced is the following:
* # p = # n = ( 1/2 + 1 ) * 1/3 = 1/2
* # Delta++ = # Delta- = 1 = ( 1/2 + 1) * 2/3 # Delta+ = # Delta0
* which is correct, and therefore the weight 1/2 for the
* diquarks of different types of quarks is justified (at least
* in this limit of exact SU(2) ).
*/
virtual void constructHadronTable() = 0;
/**
* The table of hadron data
*/
HadronTable _table;
/**
* The PDG codes of the constituent particles allowed
*/
vector<PDPtr> _partons;
/**
* The PDG codes of the hadrons which cannot be produced in the hadronization
*/
vector<PDPtr> _forbidden;
/**
* Access to the table of hadrons
*/
const HadronTable & table() const {
return _table;
}
/**
* Access to the list of partons
*/
const vector<PDPtr> & partons() const {
return _partons;
}
/**
* Calculates a special weight specific to a given hadron.
* @param id The PDG code of the hadron
*/
double specialWeight(long id) const {
const int pspin = id % 10;
// Only K0L and K0S have pspin == 0, should
// not get them until Decay step
assert( pspin != 0 );
// Baryon : J = 1/2 or 3/2
if(pspin%2==0)
return baryonWeight(id);
// Meson
else
return mesonWeight(id);
}
/**
* Weights for mesons
*/
virtual double mesonWeight(long id) const;
/**
* Weights for baryons
*/
virtual double baryonWeight(long id) const = 0;
/**
* This method returns the proper sign ( > 0 hadron; < 0 anti-hadron )
* for the input PDG id idHad > 0, suppose to be made by the
* two constituent particle pointers: par1 and par2 (both with proper sign).
*/
int signHadron(tcPDPtr ptr1, tcPDPtr ptr2, tcPDPtr hadron) const;
/**
* Insert a meson in the table
*/
virtual void insertMeson(HadronInfo a, int flav1, int flav2) = 0;
/**
* Insert a spin\f$\frac12\f$ baryon in the table
*/
virtual void insertOneHalf(HadronInfo a, int flav1, int flav2);
/**
* Insert a spin\f$\frac32\f$ baryon in the table
*/
virtual void insertThreeHalf(HadronInfo a, int flav1, int flav2);
/**
* Option for the selection of hadrons below the pair threshold
*/
unsigned int belowThreshold_;
/**
* The weights for the excited meson multiplets
*/
vector<vector<vector<double> > > _repwt;
/**
* Weights for quarks and diquarks.
*/
map<long,double> _pwt;
/**
* Enums so arrays can be statically allocated
*/
//@{
/**
* Defines values for array sizes. L,J,N max values for excited mesons.
*/
enum MesonMultiplets { Lmax = 3, Jmax = 4, Nmax = 4};
//@}
/**
* Caches of lightest pairs for speed
*/
//@{
/**
* Masses of lightest hadron pair
*/
map<pair<long,long>,tcPDPair> lightestHadrons_;
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
HadronSpectrum & operator=(const HadronSpectrum &) = delete;
};
}
#endif /* Herwig_HadronSpectrum_H */
diff --git a/Hadronization/StandardModelHadronSpectrum.cc b/Hadronization/StandardModelHadronSpectrum.cc
--- a/Hadronization/StandardModelHadronSpectrum.cc
+++ b/Hadronization/StandardModelHadronSpectrum.cc
@@ -1,718 +1,709 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the StandardModelHadronSpectrum class.
//
#include "StandardModelHadronSpectrum.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include <ThePEG/PDT/EnumParticles.h>
#include <ThePEG/Repository/EventGenerator.h>
#include <ThePEG/Repository/Repository.h>
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
namespace {
bool weightIsLess (pair<long,double> a, pair<long,double> b) {
return a.second < b.second;
}
/**
* Return true if the particle pointer corresponds to a diquark
* or anti-diquark carrying b flavour; false otherwise.
*/
inline bool isDiquarkWithB(tcPDPtr par1) {
if (!par1) return false;
long id1 = par1->id();
return DiquarkMatcher::Check(id1) && (abs(id1)/1000)%10 == ParticleID::b;
}
/**
* Return true if the particle pointer corresponds to a diquark
* or anti-diquark carrying c flavour; false otherwise.
*/
inline bool isDiquarkWithC(tcPDPtr par1) {
if (!par1) return false;
long id1 = par1->id();
return ( DiquarkMatcher::Check(id1) &&
( (abs(id1)/1000)%10 == ParticleID::c
|| (abs(id1)/100)%10 == ParticleID::c ) );
}
}
StandardModelHadronSpectrum::StandardModelHadronSpectrum(unsigned int opt)
: HadronSpectrum(),
_hadronizingStrangeDiquarks(1),
_pwtDquark( 1.0 ),_pwtUquark( 1.0 ),_pwtSquark( 1.0 ),_pwtCquark( 0.0 ),
_pwtBquark( 0.0 ),
_sngWt( 1.0 ),_decWt( 1.0 ),
_weight1S0(Nmax,1.),_weight3S1(Nmax,1.),_weight1P1(Nmax,1.),_weight3P0(Nmax,1.),
_weight3P1(Nmax,1.),_weight3P2(Nmax,1.),_weight1D2(Nmax,1.),_weight3D1(Nmax,1.),
_weight3D2(Nmax,1.),_weight3D3(Nmax,1.),
_topt(opt),_trial(0),
_limBottom(0.0), _limCharm(0.0), _limExotic(0.0)
{
// The mixing angles
// the ideal mixing angle
const double idealAngleMix = atan( sqrt(0.5) ) * 180.0 / Constants::pi;
// \eta-\eta' mixing angle
_etamix = -23.0;
// phi-omega mixing angle
_phimix = +36.0;
// h_1'-h_1 mixing angle
_h1mix = idealAngleMix;
// f_0(1710)-f_0(1370) mixing angle
_f0mix = idealAngleMix;
// f_1(1420)-f_1(1285)\f$ mixing angle
_f1mix = idealAngleMix;
// f'_2-f_2\f$ mixing angle
_f2mix = +26.0;
// eta_2(1870)-eta_2(1645) mixing angle
_eta2mix = idealAngleMix;
// phi(???)-omega(1650) mixing angle
_omhmix = idealAngleMix;
// phi_3-omega_3 mixing angle
_ph3mix = +28.0;
// eta(1475)-eta(1295) mixing angle
_eta2Smix = idealAngleMix;
// phi(1680)-omega(1420) mixing angle
_phi2Smix = idealAngleMix;
}
StandardModelHadronSpectrum::~StandardModelHadronSpectrum() {}
void StandardModelHadronSpectrum::persistentOutput(PersistentOStream & os) const {
os << _hadronizingStrangeDiquarks
<< _pwtDquark << _pwtUquark << _pwtSquark
<< _pwtCquark << _pwtBquark
<< _etamix << _phimix << _h1mix << _f0mix << _f1mix << _f2mix
<< _eta2mix << _omhmix << _ph3mix << _eta2Smix << _phi2Smix
<< _weight1S0 << _weight3S1 << _weight1P1 << _weight3P0 << _weight3P1
<< _weight3P2 << _weight1D2 << _weight3D1 << _weight3D2 << _weight3D3
<< _sngWt << _decWt << _repwt
<< _limBottom << _limCharm << _limExotic;
}
void StandardModelHadronSpectrum::persistentInput(PersistentIStream & is, int) {
is >> _hadronizingStrangeDiquarks
>> _pwtDquark >> _pwtUquark >> _pwtSquark
>> _pwtCquark >> _pwtBquark
>> _etamix >> _phimix >> _h1mix >> _f0mix >> _f1mix >> _f2mix
>> _eta2mix >> _omhmix >> _ph3mix >> _eta2Smix >> _phi2Smix
>> _weight1S0 >> _weight3S1 >> _weight1P1 >> _weight3P0 >> _weight3P1
>> _weight3P2 >> _weight1D2 >> _weight3D1 >> _weight3D2 >> _weight3D3
>> _sngWt >> _decWt >> _repwt
>> _limBottom >> _limCharm >> _limExotic;
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeAbstractClass<StandardModelHadronSpectrum,HadronSpectrum>
describeHerwigStandardModelHadronSpectrum("Herwig::StandardModelHadronSpectrum", "Herwig.so");
void StandardModelHadronSpectrum::Init() {
static ClassDocumentation<StandardModelHadronSpectrum> documentation
("There is no documentation for the StandardModelHadronSpectrum class");
static Parameter<StandardModelHadronSpectrum,double>
interfacePwtDquark("PwtDquark","Weight for choosing a quark D",
&StandardModelHadronSpectrum::_pwtDquark, 0, 1.0, 0.0, 10.0,
false,false,false);
static Parameter<StandardModelHadronSpectrum,double>
interfacePwtUquark("PwtUquark","Weight for choosing a quark U",
&StandardModelHadronSpectrum::_pwtUquark, 0, 1.0, 0.0, 10.0,
false,false,false);
static Parameter<StandardModelHadronSpectrum,double>
interfacePwtSquark("PwtSquark","Weight for choosing a quark S",
&StandardModelHadronSpectrum::_pwtSquark, 0, 1.0, 0.0, 10.0,
false,false,false);
static Parameter<StandardModelHadronSpectrum,double>
interfacePwtCquark("PwtCquark","Weight for choosing a quark C",
&StandardModelHadronSpectrum::_pwtCquark, 0, 0.0, 0.0, 10.0,
false,false,false);
static Parameter<StandardModelHadronSpectrum,double>
interfacePwtBquark("PwtBquark","Weight for choosing a quark B",
&StandardModelHadronSpectrum::_pwtBquark, 0, 0.0, 0.0, 10.0,
false,false,false);
static Parameter<StandardModelHadronSpectrum,double>
interfaceSngWt("SngWt","Weight for singlet baryons",
&StandardModelHadronSpectrum::_sngWt, 0, 1.0, 0.0, 10.0,
false,false,false);
static Parameter<StandardModelHadronSpectrum,double>
interfaceDecWt("DecWt","Weight for decuplet baryons",
&StandardModelHadronSpectrum::_decWt, 0, 1.0, 0.0, 10.0,
false,false,false);
//
// mixing angles
//
// the ideal mixing angle
const double idealAngleMix = atan( sqrt(0.5) ) * 180.0 / Constants::pi;
static Parameter<StandardModelHadronSpectrum,double> interface11S0Mixing
("11S0Mixing",
"The mixing angle for the I=0 mesons from the 1 1S0 multiplet,"
" i.e. eta and etaprime.",
&StandardModelHadronSpectrum::_etamix, -23., -180., 180.,
false, false, Interface::limited);
static Parameter<StandardModelHadronSpectrum,double> interface13S1Mixing
("13S1Mixing",
"The mixing angle for the I=0 mesons from the 1 3S1 multiplet,"
" i.e. phi and omega.",
&StandardModelHadronSpectrum::_phimix, +36., -180., 180.,
false, false, Interface::limited);
static Parameter<StandardModelHadronSpectrum,double> interface11P1Mixing
("11P1Mixing",
"The mixing angle for the I=0 mesons from the 1 1P1 multiplet,"
" i.e. h_1' and h_1.",
&StandardModelHadronSpectrum::_h1mix, idealAngleMix, -180., 180.,
false, false, Interface::limited);
static Parameter<StandardModelHadronSpectrum,double> interface13P0Mixing
("13P0Mixing",
"The mixing angle for the I=0 mesons from the 1 3P0 multiplet,"
" i.e. f_0(1710) and f_0(1370).",
&StandardModelHadronSpectrum::_f0mix, idealAngleMix, -180., 180.,
false, false, Interface::limited);
static Parameter<StandardModelHadronSpectrum,double> interface13P1Mixing
("13P1Mixing",
"The mixing angle for the I=0 mesons from the 1 3P1 multiplet,"
" i.e. f_1(1420) and f_1(1285).",
&StandardModelHadronSpectrum::_f1mix, idealAngleMix, -180., 180.,
false, false, Interface::limited);
static Parameter<StandardModelHadronSpectrum,double> interface13P2Mixing
("13P2Mixing",
"The mixing angle for the I=0 mesons from the 1 3P2 multiplet,"
" i.e. f'_2 and f_2.",
&StandardModelHadronSpectrum::_f2mix, 26.0, -180., 180.,
false, false, Interface::limited);
static Parameter<StandardModelHadronSpectrum,double> interface11D2Mixing
("11D2Mixing",
"The mixing angle for the I=0 mesons from the 1 1D2 multiplet,"
" i.e. eta_2(1870) and eta_2(1645).",
&StandardModelHadronSpectrum::_eta2mix, idealAngleMix, -180., 180.,
false, false, Interface::limited);
static Parameter<StandardModelHadronSpectrum,double> interface13D0Mixing
("13D0Mixing",
"The mixing angle for the I=0 mesons from the 1 3D0 multiplet,"
" i.e. eta_2(1870) phi(?) and omega(1650).",
&StandardModelHadronSpectrum::_omhmix, idealAngleMix, -180., 180.,
false, false, Interface::limited);
static Parameter<StandardModelHadronSpectrum,double> interface13D1Mixing
("13D1Mixing",
"The mixing angle for the I=0 mesons from the 1 3D1 multiplet,"
" i.e. phi_3 and omega_3.",
&StandardModelHadronSpectrum::_ph3mix, 28.0, -180., 180.,
false, false, Interface::limited);
static Parameter<StandardModelHadronSpectrum,double> interface21S0Mixing
("21S0Mixing",
"The mixing angle for the I=0 mesons from the 2 1S0 multiplet,"
" i.e. eta(1475) and eta(1295).",
&StandardModelHadronSpectrum::_eta2Smix, idealAngleMix, -180., 180.,
false, false, Interface::limited);
static Parameter<StandardModelHadronSpectrum,double> interface23S1Mixing
("23S1Mixing",
"The mixing angle for the I=0 mesons from the 1 3S1 multiplet,"
" i.e. phi(1680) and omega(1420).",
&StandardModelHadronSpectrum::_phi2Smix, idealAngleMix, -180., 180.,
false, false, Interface::limited);
//
// the meson weights
//
static ParVector<StandardModelHadronSpectrum,double> interface1S0Weights
("1S0Weights",
"The weights for the 1S0 multiplets start with n=1.",
&StandardModelHadronSpectrum::_weight1S0, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<StandardModelHadronSpectrum,double> interface3S1Weights
("3S1Weights",
"The weights for the 3S1 multiplets start with n=1.",
&StandardModelHadronSpectrum::_weight3S1, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<StandardModelHadronSpectrum,double> interface1P1Weights
("1P1Weights",
"The weights for the 1P1 multiplets start with n=1.",
&StandardModelHadronSpectrum::_weight1P1, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<StandardModelHadronSpectrum,double> interface3P0Weights
("3P0Weights",
"The weights for the 3P0 multiplets start with n=1.",
&StandardModelHadronSpectrum::_weight3P0, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<StandardModelHadronSpectrum,double> interface3P1Weights
("3P1Weights",
"The weights for the 3P1 multiplets start with n=1.",
&StandardModelHadronSpectrum::_weight3P1, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<StandardModelHadronSpectrum,double> interface3P2Weights
("3P2Weights",
"The weights for the 3P2 multiplets start with n=1.",
&StandardModelHadronSpectrum::_weight3P2, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<StandardModelHadronSpectrum,double> interface1D2Weights
("1D2Weights",
"The weights for the 1D2 multiplets start with n=1.",
&StandardModelHadronSpectrum::_weight1D2, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<StandardModelHadronSpectrum,double> interface3D1Weights
("3D1Weights",
"The weights for the 3D1 multiplets start with n=1.",
&StandardModelHadronSpectrum::_weight3D1, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<StandardModelHadronSpectrum,double> interface3D2Weights
("3D2Weights",
"The weights for the 3D2 multiplets start with n=1.",
&StandardModelHadronSpectrum::_weight3D2, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<StandardModelHadronSpectrum,double> interface3D3Weights
("3D3Weights",
"The weights for the 3D3 multiplets start with n=1.",
&StandardModelHadronSpectrum::_weight3D3, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static Switch<StandardModelHadronSpectrum,unsigned int> interfaceTrial
("Trial",
"A Debugging option to only produce certain types of hadrons",
&StandardModelHadronSpectrum::_trial, 0, false, false);
static SwitchOption interfaceTrialAll
(interfaceTrial,
"All",
"Produce all the hadrons",
0);
static SwitchOption interfaceTrialPions
(interfaceTrial,
"Pions",
"Only produce pions",
1);
static SwitchOption interfaceTrialSpin2
(interfaceTrial,
"Spin2",
"Only mesons with spin less than or equal to two are produced",
2);
static SwitchOption interfaceTrialSpin3
(interfaceTrial,
"Spin3",
"Only hadrons with spin less than or equal to three are produced",
3);
static Parameter<StandardModelHadronSpectrum,double>
interfaceSingleHadronLimitBottom ("SingleHadronLimitBottom",
"Threshold for one-hadron decay of b-cluster",
&StandardModelHadronSpectrum::_limBottom,
0, 0.0, 0.0, 100.0,false,false,false);
static Parameter<StandardModelHadronSpectrum,double>
interfaceSingleHadronLimitCharm ("SingleHadronLimitCharm",
"threshold for one-hadron decay of c-cluster",
&StandardModelHadronSpectrum::_limCharm,
0, 0.0, 0.0, 100.0,false,false,false);
static Parameter<StandardModelHadronSpectrum,double>
interfaceSingleHadronLimitExotic ("SingleHadronLimitExotic",
"threshold for one-hadron decay of exotic cluster",
&StandardModelHadronSpectrum::_limExotic,
0, 0.0, 0.0, 100.0,false,false,false);
static Switch<StandardModelHadronSpectrum,unsigned int> interfaceBelowThreshold
("BelowThreshold",
"Option fo the selection of the hadrons if the cluster is below the pair threshold",
&StandardModelHadronSpectrum::belowThreshold_, 0, false, false);
static SwitchOption interfaceBelowThresholdLightest
(interfaceBelowThreshold,
"Lightest",
"Force cluster to decay to the lightest hadron with the appropriate flavours",
0);
static SwitchOption interfaceBelowThresholdAll
(interfaceBelowThreshold,
"All",
"Select from all the hadrons below the two hadron threshold according to their spin weights",
1);
static Switch<StandardModelHadronSpectrum,unsigned int> interfaceHadronizingStrangeDiquarks
("HadronizingStrangeDiquarks",
"Option for adding strange diquarks to Hadronization",
&StandardModelHadronSpectrum::_hadronizingStrangeDiquarks, 0, false, false);
static SwitchOption interfaceHadronizingStrangeDiquarksNo
(interfaceHadronizingStrangeDiquarks,
"No",
"No strangeness containing diquarks in Hadronization",
0);
static SwitchOption interfaceHadronizingStrangeDiquarksOnlySingleStrange
(interfaceHadronizingStrangeDiquarks,
"OnlySingleStrange",
"Only one strangeness containing diquarks in Hadronization i.e. su,sd",
1);
static SwitchOption interfaceHadronizingStrangeDiquarksAll
(interfaceHadronizingStrangeDiquarks,
"All",
"All strangeness containing diquarks in Hadronization i.e. su,sd,ss",
2);
}
-
-PDPtr StandardModelHadronSpectrum::makeDiquark(tcPDPtr par1, tcPDPtr par2) const {
- long id1 = par1->id();
- long id2 = par2->id();
- long pspin = id1==id2 ? 3 : 1;
- long idnew = makeDiquarkID(id1,id2, pspin);
- return getParticleData(idnew);
-}
-
Energy StandardModelHadronSpectrum::hadronPairThreshold(tcPDPtr par1, tcPDPtr par2) const {
// Determine the sum of the nominal masses of the two lightest hadrons
// with the right flavour numbers as the cluster under consideration.
// Notice that we don't need real masses (drawn by a Breit-Wigner
// distribution) because the lightest pair of hadrons does not involve
// any broad resonance.
Energy threshold = massLightestHadronPair(par1,par2);
// Special: it allows one-hadron decays also above threshold.
if (isExotic(par1,par2))
threshold *= (1.0 + UseRandom::rnd()*_limExotic);
else if (hasBottom(par1,par2))
threshold *= (1.0 + UseRandom::rnd()*_limBottom);
else if (hasCharm(par1,par2))
threshold *= (1.0 + UseRandom::rnd()*_limCharm);
return threshold;
}
double StandardModelHadronSpectrum::mixingStateWeight(long id) const {
switch(id) {
case ParticleID::eta: return 0.5*probabilityMixing(_etamix ,1);
case ParticleID::etaprime: return 0.5*probabilityMixing(_etamix ,2);
case ParticleID::phi: return 0.5*probabilityMixing(_phimix ,1);
case ParticleID::omega: return 0.5*probabilityMixing(_phimix ,2);
case ParticleID::hprime_1: return 0.5*probabilityMixing(_h1mix ,1);
case ParticleID::h_1: return 0.5*probabilityMixing(_h1mix ,2);
case 10331: return 0.5*probabilityMixing(_f0mix ,1);
case 10221: return 0.5*probabilityMixing(_f0mix ,2);
case ParticleID::fprime_1: return 0.5*probabilityMixing(_f1mix ,1);
case ParticleID::f_1: return 0.5*probabilityMixing(_f1mix ,2);
case ParticleID::fprime_2: return 0.5*probabilityMixing(_f2mix ,1);
case ParticleID::f_2: return 0.5*probabilityMixing(_f2mix ,2);
case 10335: return 0.5*probabilityMixing(_eta2mix ,1);
case 10225: return 0.5*probabilityMixing(_eta2mix ,2);
case 30333: return 0.5*probabilityMixing(_omhmix ,1);
case 30223: return 0.5*probabilityMixing(_omhmix ,2);
case 337: return 0.5*probabilityMixing(_ph3mix ,1);
case 227: return 0.5*probabilityMixing(_ph3mix ,2);
case 100331: return 0.5*probabilityMixing(_eta2mix ,1);
case 100221: return 0.5*probabilityMixing(_eta2mix ,2);
case 100333: return 0.5*probabilityMixing(_phi2Smix,1);
case 100223: return 0.5*probabilityMixing(_phi2Smix,2);
default: return 1./3.;
}
}
void StandardModelHadronSpectrum::doinit() {
// set the weights for the various excited mesons
// set all to one to start with
for (int l = 0; l < Lmax; ++l ) {
for (int j = 0; j < Jmax; ++j) {
for (int n = 0; n < Nmax; ++n) {
_repwt[l][j][n] = 1.0;
}
}
}
// set the others from the relevant vectors
for( int ix=0;ix<max(int(_weight1S0.size()),int(Nmax));++ix)
_repwt[0][0][ix]=_weight1S0[ix];
for( int ix=0;ix<max(int(_weight3S1.size()),int(Nmax));++ix)
_repwt[0][1][ix]=_weight3S1[ix];
for( int ix=0;ix<max(int(_weight1P1.size()),int(Nmax));++ix)
_repwt[1][1][ix]=_weight1P1[ix];
for( int ix=0;ix<max(int(_weight3P0.size()),int(Nmax));++ix)
_repwt[1][0][ix]=_weight3P0[ix];
for( int ix=0;ix<max(int(_weight3P1.size()),int(Nmax));++ix)
_repwt[1][1][ix]=_weight3P1[ix];
for( int ix=0;ix<max(int(_weight3P2.size()),int(Nmax));++ix)
_repwt[1][2][ix]=_weight3P2[ix];
for( int ix=0;ix<max(int(_weight1D2.size()),int(Nmax));++ix)
_repwt[2][2][ix]=_weight1D2[ix];
for( int ix=0;ix<max(int(_weight3D1.size()),int(Nmax));++ix)
_repwt[2][1][ix]=_weight3D1[ix];
for( int ix=0;ix<max(int(_weight3D2.size()),int(Nmax));++ix)
_repwt[2][2][ix]=_weight3D2[ix];
for( int ix=0;ix<max(int(_weight3D3.size()),int(Nmax));++ix)
_repwt[2][3][ix]=_weight3D3[ix];
// find the maximum
map<long,double>::iterator pit =
max_element(_pwt.begin(),_pwt.end(),weightIsLess);
const double pmax = pit->second;
for(pit=_pwt.begin(); pit!=_pwt.end(); ++pit) {
pit->second/=pmax;
}
HadronSpectrum::doinit();
}
void StandardModelHadronSpectrum::constructHadronTable() {
// initialise the table
_table.clear();
for(unsigned int ix=0; ix<_partons.size(); ++ix) {
for(unsigned int iy=0; iy<_partons.size(); ++iy) {
if (!(DiquarkMatcher::Check(_partons[ix]->id())
&& DiquarkMatcher::Check(_partons[iy]->id())))
_table[make_pair(_partons[ix]->id(),_partons[iy]->id())] = KupcoData();
}
}
// get the particles from the event generator
ParticleMap particles = generator()->particles();
// loop over the particles
//double maxdd(0.),maxss(0.),maxrest(0.);
for(ParticleMap::iterator it=particles.begin();
it!=particles.end(); ++it) {
long pid = it->first;
tPDPtr particle = it->second;
int pspin = particle->iSpin();
// Don't include hadrons which are explicitly forbidden
if(find(_forbidden.begin(),_forbidden.end(),particle)!=_forbidden.end())
continue;
// Don't include non-hadrons or antiparticles
if(pid < 100) continue;
// remove diffractive particles
if(pspin == 0) continue;
// K_0S and K_0L not made make K0 and Kbar0
if(pid==ParticleID::K_S0||pid==ParticleID::K_L0) continue;
// Debugging options
// Only include those with 2J+1 less than...5
if(_trial==2 && pspin >= 5) continue;
// Only include those with 2J+1 less than...7
if(_trial==3 && pspin >= 7) continue;
// Only include pions
if(_trial==1 && pid!=111 && pid!=211) continue;
// shouldn't be coloured
if(particle->coloured()) continue;
// Get the flavours
const int x4 = (pid/1000)%10;
const int x3 = (pid/100 )%10;
const int x2 = (pid/10 )%10;
const int x7 = (pid/1000000)%10;
const bool wantSusy = x7 == 1 || x7 == 2;
// Skip non-hadrons (susy particles, etc...)
if(x3 == 0 || x2 == 0) continue;
// Skip particles which are neither SM nor SUSY
if(x7 >= 3 && x7 != 9) continue;
int flav1,flav2;
// meson
if(x4 == 0) {
flav1 = x2;
flav2 = x3;
}
// baryon
else {
flav2 = x4;
// insert the spin 1 diquark, sort out the rest later
flav1 = makeDiquarkID(x2,x3,3);
}
if (wantSusy) flav2 += 1000000 * x7;
insertToHadronTable(particle,flav1,flav2);
}
// normalise the weights
if(_topt == 0) {
HadronTable::const_iterator tit;
KupcoData::iterator it;
for(tit=_table.begin();tit!=_table.end();++tit) {
double weight=0;
for(it = tit->second.begin(); it!=tit->second.end(); ++it)
weight=max(weight,it->overallWeight);
weight = 1./weight;
}
// double weight;
// if(tit->first.first==tit->first.second) {
// if(tit->first.first==1||tit->first.first==2) weight=1./maxdd;
// else if (tit->first.first==3) weight=1./maxss;
// else weight=1./maxrest;
// }
// else weight=1./maxrest;
// for(it = tit->second.begin(); it!=tit->second.end(); ++it) {
// it->rescale(weight);
// }
// }
}
}
double StandardModelHadronSpectrum::strangeWeight(const Energy, tcPDPtr, tcPDPtr) const {
assert(false);
}
void StandardModelHadronSpectrum::insertMeson(HadronInfo a, int flav1, int flav2) {
// identical light flavours
if(flav1 == flav2 && flav1<=3) {
// ddbar> uubar> admixture states
if(flav1==1) {
a.overallWeight *= 0.5;
_table[make_pair(1,1)].insert(a);
_table[make_pair(2,2)].insert(a);
}
// load up ssbar> uubar> ddbar> admixture states
else {
// uubar ddbar pieces
a.wt = mixingStateWeight(a.id);
a.overallWeight *= a.wt;
_table[make_pair(1,1)].insert(a);
_table[make_pair(2,2)].insert(a);
a.overallWeight /=a.wt;
// ssbar piece
a.wt = 1.- 2.*a.wt;
if(a.wt > 0) {
a.overallWeight *= a.wt;
_table[make_pair(3,3)].insert(a);
}
}
}
else {
_table[make_pair(flav1,flav2)].insert(a);
if(flav1 != flav2) _table[make_pair(flav2,flav1)].insert(a);
}
}
long StandardModelHadronSpectrum::makeDiquarkID(long id1, long id2, long pspin) const {
assert( id1 * id2 > 0
&& QuarkMatcher::Check(id1)
&& QuarkMatcher::Check(id2)) ;
long ida = abs(id1);
long idb = abs(id2);
if (ida < idb) swap(ida,idb);
if (pspin != 1 && pspin != 3) assert(false);
long idnew = ida*1000 + idb*100 + pspin;
// Diquarks made of quarks of the same type: uu, dd, ss, cc, bb,
// have spin 1, and therefore the less significant digit (which
// corresponds to 2*J+1) is 3 rather than 1 as all other Diquarks.
if (id1 == id2 && pspin == 1) {
//cerr<<"WARNING: spin-0 diquiark of the same type cannot exist."
// <<" Switching to spin-1 diquark.\n";
idnew = ida*1000 + idb*100 + 3;
}
return id1 > 0 ? idnew : -idnew;
}
bool StandardModelHadronSpectrum::hasBottom(tcPDPtr par1, tcPDPtr par2, tcPDPtr par3) const {
long id1 = par1 ? par1->id() : 0;
if ( !par2 && !par3 ) {
return
abs(id1) == ThePEG::ParticleID::b ||
isDiquarkWithB(par1) ||
( MesonMatcher::Check(id1)
&& (abs(id1)/100)%10 == ThePEG::ParticleID::b ) ||
( BaryonMatcher::Check(id1)
&& (abs(id1)/1000)%10 == ThePEG::ParticleID::b );
}
else {
long id2 = par2 ? par2->id() : 0;
long id3 = par3 ? par3->id() : 0;
return
abs(id1) == ThePEG::ParticleID::b || isDiquarkWithB(par1) ||
abs(id2) == ThePEG::ParticleID::b || isDiquarkWithB(par2) ||
abs(id3) == ThePEG::ParticleID::b || isDiquarkWithB(par3);
}
}
bool StandardModelHadronSpectrum::hasCharm(tcPDPtr par1, tcPDPtr par2, tcPDPtr par3) const {
long id1 = par1 ? par1->id(): 0;
if (!par2 && !par3) {
return
abs(id1) == ThePEG::ParticleID::c ||
isDiquarkWithC(par1) ||
( MesonMatcher::Check(id1) &&
((abs(id1)/100)%10 == ThePEG::ParticleID::c ||
(abs(id1)/10)%10 == ThePEG::ParticleID::c) ) ||
( BaryonMatcher::Check(id1) &&
((abs(id1)/1000)%10 == ThePEG::ParticleID::c ||
(abs(id1)/100)%10 == ThePEG::ParticleID::c ||
(abs(id1)/10)%10 == ThePEG::ParticleID::c) );
}
else {
long id2 = par2 ? par1->id(): 0;
long id3 = par3 ? par1->id(): 0;
return
abs(id1) == ThePEG::ParticleID::c || isDiquarkWithC(par1) ||
abs(id2) == ThePEG::ParticleID::c || isDiquarkWithC(par2) ||
abs(id3) == ThePEG::ParticleID::c || isDiquarkWithC(par3);
}
}
bool StandardModelHadronSpectrum::isExotic(tcPDPtr par1, tcPDPtr par2, tcPDPtr par3) const {
/// \todo make this more general
long id1 = par1 ? par1->id(): 0;
long id2 = par2 ? par2->id(): 0;
long id3 = par3 ? par3->id(): 0;
return
( (id1/1000000)% 10 != 0 && (id1/1000000)% 10 != 9 ) ||
( (id2/1000000)% 10 != 0 && (id2/1000000)% 10 != 9 ) ||
( (id3/1000000)% 10 != 0 && (id3/1000000)% 10 != 9 ) ||
abs(id1)==6||abs(id2)==6;
}
bool StandardModelHadronSpectrum::canBeBaryon(tcPDPtr par1, tcPDPtr par2 , tcPDPtr par3) const {
assert(par1 && par2);
long id1 = par1->id(), id2 = par2->id();
if (!par3) {
if( id1*id2 < 0) return false;
if(DiquarkMatcher::Check(id1))
return abs(int(par2->iColour())) == 3 && !DiquarkMatcher::Check(id2);
if(DiquarkMatcher::Check(id2))
return abs(int(par1->iColour())) == 3;
return false;
}
else {
// In this case, to be a baryon, all three components must be (anti-)quarks
// and with the same sign.
return (par1->iColour() == 3 && par2->iColour() == 3 && par3->iColour() == 3) ||
(par1->iColour() == -3 && par2->iColour() == -3 && par3->iColour() == -3);
}
}
diff --git a/Hadronization/StandardModelHadronSpectrum.h b/Hadronization/StandardModelHadronSpectrum.h
--- a/Hadronization/StandardModelHadronSpectrum.h
+++ b/Hadronization/StandardModelHadronSpectrum.h
@@ -1,622 +1,614 @@
// -*- C++ -*-
#ifndef Herwig_StandardModelHadronSpectrum_H
#define Herwig_StandardModelHadronSpectrum_H
//
// This is the declaration of the StandardModelHadronSpectrum class.
//
#include "Herwig/Hadronization/HadronSpectrum.h"
#include <ThePEG/PDT/ParticleData.h>
#include <ThePEG/PDT/StandardMatchers.h>
#include <ThePEG/Repository/EventGenerator.h>
#include <ThePEG/PDT/EnumParticles.h>
#include "ThePEG/Repository/CurrentGenerator.h"
#include <ThePEG/Persistency/PersistentOStream.h>
#include <ThePEG/Persistency/PersistentIStream.h>
namespace Herwig {
using namespace ThePEG;
/**
* Here is the documentation of the StandardModelHadronSpectrum class.
*
* @see \ref StandardModelHadronSpectrumInterfaces "The interfaces"
* defined for StandardModelHadronSpectrum.
*/
class StandardModelHadronSpectrum: public HadronSpectrum {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
StandardModelHadronSpectrum(unsigned int opt);
/**
* The destructor.
*/
virtual ~StandardModelHadronSpectrum();
//@}
public:
/** @name Partonic content */
//@{
/**
* Return the id of the gluon
*/
virtual long gluonId() const { return ParticleID::g; }
/**
* Return the ids of all hadronizing quarks
*/
virtual const vector<long>& hadronizingQuarks() const {
static vector<long> hadronizing =
{ ParticleID::d, ParticleID::u, ParticleID::s, ParticleID::c, ParticleID::b };
return hadronizing;
}
/**
* The light hadronizing quarks
*/
virtual const vector<long>& lightHadronizingQuarks() const {
static vector<long> light =
{ ParticleID::d, ParticleID::u, ParticleID::s };
return light;
}
/**
* The light hadronizing diquarks
*/
virtual const vector<long>& lightHadronizingDiquarks() const {
/**
* Diquarks q==q_0 are not allowed as they need to have antisymmetric
* spin wave-function, which forces the spin to 1
* Diquarks q!=q'_1 are not allowed as they need to have antisymmetric
* spin wave-function, which forces the spin to 1
* */
// TODO: strange diquarks are turned off for the moment
// since in combination with the current ClusterFission
// they fail (overshoot) to reproduce the Xi and Lambda
// pT spectra.
// One may enable these after the ClusterFission
// kinematics are settled
// TODO why ud_1 not allowed?
// exceptions: Could not find 2103 1 in _table
// but no problem for 2103 2 ???
// ParticleID::ud_1
switch(_hadronizingStrangeDiquarks) {
case 0:
{
static vector<long> light = {
ParticleID::uu_1,
ParticleID::dd_1,
ParticleID::ud_0
// ParticleID::ud_1,
};
return light;
// break;
}
case 1:
{
static vector<long> light = {
ParticleID::uu_1,
ParticleID::dd_1,
ParticleID::ud_0,
// ParticleID::ud_1,
ParticleID::su_0,
// ParticleID::su_1,
ParticleID::sd_0
// ParticleID::sd_1
};
return light;
}
case 2:
{
static vector<long> light = {
ParticleID::uu_1,
ParticleID::dd_1,
ParticleID::ud_0,
// ParticleID::ud_1,
ParticleID::su_0,
// ParticleID::su_1,
ParticleID::sd_0,
// ParticleID::sd_1,
ParticleID::ss_1
};
return light;
// break;
}
default:
assert(false);
}
static vector<long> light;
return light;
}
/**
* The heavy hadronizing quarks
*/
virtual const vector<long>& heavyHadronizingQuarks() const {
static vector<long> heavy =
{ ParticleID::c, ParticleID::b };
return heavy;
}
/**
* Return true if any of the possible three input particles contains
* the indicated heavy quark. false otherwise. In the case that
* only the first particle is specified, it can be: an (anti-)quark,
* an (anti-)diquark an (anti-)meson, an (anti-)baryon; in the other
* cases, each pointer is assumed to be either (anti-)quark or
* (anti-)diquark.
*/
virtual bool hasHeavy(long id, tcPDPtr par1, tcPDPtr par2 = PDPtr(), tcPDPtr par3 = PDPtr()) const {
if ( abs(id) == ParticleID::c )
return hasCharm(par1,par2,par3);
if ( abs(id) == ParticleID::b )
return hasBottom(par1,par2,par3);
return false;
}
//@}
/**
* Return the threshold for a cluster to split into a pair of hadrons.
* This is normally the mass of the lightest hadron Pair, but can be
* higher for heavy and exotic clusters
*/
virtual Energy hadronPairThreshold(tcPDPtr par1, tcPDPtr par2) const;
/**
* Return the weight for the given flavour
*/
virtual double pwtQuark(const long& id) const {
switch(id) {
case ParticleID::d: return pwtDquark(); break;
case ParticleID::u: return pwtUquark(); break;
case ParticleID::s: return pwtSquark(); break;
case ParticleID::c: return pwtCquark(); break;
case ParticleID::b: return pwtBquark(); break;
}
return 0.;
}
/**
* The down quark weight.
*/
double pwtDquark() const {
return _pwtDquark;
}
/**
* The up quark weight.
*/
double pwtUquark() const {
return _pwtUquark;
}
/**
* The strange quark weight.
*/
double pwtSquark() const {
return _pwtSquark;
}
/**
* The charm quark weight.
*/
double pwtCquark() const {
return _pwtCquark;
}
/**
* The bottom quark weight.
*/
double pwtBquark() const {
return _pwtBquark;
}
/**
* The diquark weight.
*/
double pwtDIquark() const {
return _pwtDIquark;
}
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
- /**
- * Return the particle data of the diquark (anti-diquark) made by the two
- * quarks (antiquarks) par1, par2.
- * @param par1 (anti-)quark data pointer
- * @param par2 (anti-)quark data pointer
- */
- PDPtr makeDiquark(tcPDPtr par1, tcPDPtr par2) const;
-
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
*
* The array _repwt is initialized using the interfaces to set different
* weights for different meson multiplets and the constructHadronTable()
* method called to complete the construction of the hadron tables.
*
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
//@}
/**
* Return the id of the diquark (anti-diquark) made by the two
* quarks (antiquarks) of id specified in input (id1, id2).
* Caller must ensure that id1 and id2 are quarks.
*/
long makeDiquarkID(long id1, long id2, long pspin) const;
/**
* Return true if any of the possible three input particles has
* b-flavour;
* false otherwise. In the case that only the first particle is specified,
* it can be: an (anti-)quark, an (anti-)diquark
* an (anti-)meson, an (anti-)baryon; in the other cases, each pointer
* is assumed to be either (anti-)quark or (anti-)diquark.
*/
bool hasBottom(tcPDPtr par1, tcPDPtr par2 = PDPtr(), tcPDPtr par3 = PDPtr()) const;
/**
* Return true if any of the possible three input particles has
* c-flavour;
* false otherwise.In the case that only the first pointer is specified,
* it can be: a (anti-)quark, a (anti-)diquark
* a (anti-)meson, a (anti-)baryon; in the other cases, each pointer
* is assumed to be either (anti-)quark or (anti-)diquark.
*/
bool hasCharm(tcPDPtr par1, tcPDPtr par2 = PDPtr(), tcPDPtr par3 = PDPtr()) const;
/**
* Return true, if any of the possible input particle pointer is an exotic quark, e.g. Susy quark;
* false otherwise.
*/
bool isExotic(tcPDPtr par1, tcPDPtr par2 = PDPtr(), tcPDPtr par3 = PDPtr()) const;
/**
* Return true if the two or three particles in input can be the components
* of a baryon; false otherwise.
*/
virtual bool canBeBaryon(tcPDPtr par1, tcPDPtr par2 , tcPDPtr par3 = PDPtr()) const;
protected:
/**
* Construct the table of hadron data
* This is the main method to initialize the hadron data (mainly the
* weights associated to each hadron, taking into account its spin,
* eventual isoscalar-octect mixing, singlet-decuplet factor). This is
* the method that one should update when new or updated hadron data is
* available.
*
* This class implements the construction of the basic table but can be
* overridden if needed in inheriting classes.
*
* The rationale for factors used for diquarks involving different quarks can
* be can be explained by taking a prototype example that in the exact SU(2) limit,
* in which:
* \f[m_u=m_d\f]
* \f[M_p=M_n=M_\Delta\f]
* and we will have equal numbers of u and d quarks produced.
* Suppose that we weight 1 the diquarks made of the same
* quark and 1/2 those made of different quarks, the fractions
* of u and d baryons (p, n, Delta) we get are the following:
* - \f$\Delta^{++}\f$: 1 possibility only u uu with weight 1
* - \f$\Delta^- \f$: 1 possibility only d dd with weight 1
* - \f$p,\Delta^+ \f$: 2 possibilities u ud with weight 1/2
* d uu with weight 1
* - \f$n,\Delta^0 \f$: 2 possibilities d ud with weight 1/2
* u dd with weight 1
* In the latter two cases, we have to take into account the
* fact that p and n have spin 1/2 whereas Delta+ and Delta0
* have spin 3/2 therefore from phase space we get a double weight
* for Delta+ and Delta0 relative to p and n respectively.
* Therefore the relative amount of these baryons that is
* produced is the following:
* # p = # n = ( 1/2 + 1 ) * 1/3 = 1/2
* # Delta++ = # Delta- = 1 = ( 1/2 + 1) * 2/3 # Delta+ = # Delta0
* which is correct, and therefore the weight 1/2 for the
* diquarks of different types of quarks is justified (at least
* in this limit of exact SU(2) ).
*/
virtual void constructHadronTable();
/**
* Insert a meson in the table
*/
virtual void insertMeson(HadronInfo a, int flav1, int flav2);
/**
* Methods for the mixing of \f$I=0\f$ mesons
*/
//@{
/**
* Return the probability of mixing for Octet-Singlet isoscalar mixing,
* the probability of the
* \f$\frac1{\sqrt{2}}(|u\bar{u}\rangle + |d\bar{d}\rangle)\f$ component
* is returned.
* @param angleMix The mixing angle in degrees (not radians)
* @param order is 0 for no mixing, 1 for the first resonance of a pair,
* 2 for the second one.
* The mixing is defined so that for example with \f$\eta-\eta'\f$ mixing where
* the mixing angle is \f$\theta=-23^0$ with $\eta\f$ as the first particle
* and \f$\eta'\f$ the second one.
* The convention used is
* \f[\eta = \cos\theta|\eta_{\rm octet }\rangle
* -\sin\theta|\eta_{\rm singlet}\rangle\f]
* \f[\eta' = \sin\theta|\eta_{\rm octet }\rangle
* -\cos\theta|\eta_{\rm singlet}\rangle\f]
* with
* \f[|\eta_{\rm singlet}\rangle = \frac1{\sqrt{3}}
* \left[|u\bar{u}\rangle + |d\bar{d}\rangle + |s\bar{s}\rangle\right]\f]
* \f[|\eta_{\rm octet }\rangle = \frac1{\sqrt{6}}
* \left[|u\bar{u}\rangle + |d\bar{d}\rangle - 2|s\bar{s}\rangle\right]\f]
*/
double probabilityMixing(const double angleMix,
const int order) const {
static double convert=Constants::pi/180.0;
if (order == 1)
return sqr( cos( angleMix*convert + atan( sqrt(2.0) ) ) );
else if (order == 2)
return sqr( sin( angleMix*convert + atan( sqrt(2.0) ) ) );
else
return 1.;
}
/**
* Returns the weight of given mixing state.
* @param id The PDG code of the meson
*/
virtual double mixingStateWeight(long id) const;
//@}
virtual double specialQuarkWeight(double quarkWeight, long id,
const Energy cluMass, tcPDPtr par1, tcPDPtr par2) const {
// special for strange
if(abs(id) == 3)
return strangeWeight(cluMass,par1,par2);
else
return quarkWeight;
}
/**
* Strange quark weight
*/
virtual double strangeWeight(const Energy cluMass, tcPDPtr par1, tcPDPtr par2) const;
/**
* Strange diquark option for Hadronization
*/
unsigned int _hadronizingStrangeDiquarks;
/**
* The weights for the different quarks and diquarks
*/
//@{
/**
* The probability of producting a down quark.
*/
double _pwtDquark;
/**
* The probability of producting an up quark.
*/
double _pwtUquark;
/**
* The probability of producting a strange quark.
*/
double _pwtSquark;
/**
* The probability of producting a charm quark.
*/
double _pwtCquark;
/**
* The probability of producting a bottom quark.
*/
double _pwtBquark;
//@}
/**
* The probability of producting a diquark.
*/
double _pwtDIquark;
/**
* Singlet and Decuplet weights
*/
//@{
/**
* The singlet weight
*/
double _sngWt;
/**
* The decuplet weight
*/
double _decWt;
//@}
/**
* The mixing angles for the \f$I=0\f$ mesons containing light quarks
*/
//@{
/**
* The \f$\eta-\eta'\f$ mixing angle
*/
double _etamix;
/**
* The \f$\phi-\omega\f$ mixing angle
*/
double _phimix;
/**
* The \f$h_1'-h_1\f$ mixing angle
*/
double _h1mix;
/**
* The \f$f_0(1710)-f_0(1370)\f$ mixing angle
*/
double _f0mix;
/**
* The \f$f_1(1420)-f_1(1285)\f$ mixing angle
*/
double _f1mix;
/**
* The \f$f'_2-f_2\f$ mixing angle
*/
double _f2mix;
/**
* The \f$\eta_2(1870)-\eta_2(1645)\f$ mixing angle
*/
double _eta2mix;
/**
* The \f$\phi(???)-\omega(1650)\f$ mixing angle
*/
double _omhmix;
/**
* The \f$\phi_3-\omega_3\f$ mixing angle
*/
double _ph3mix;
/**
* The \f$\eta(1475)-\eta(1295)\f$ mixing angle
*/
double _eta2Smix;
/**
* The \f$\phi(1680)-\omega(1420)\f$ mixing angle
*/
double _phi2Smix;
//@}
/**
* The weights for the various meson multiplets to be used to supress the
* production of particular states
*/
//@{
/**
* The weights for the \f$\phantom{1}^1S_0\f$ multiplets
*/
vector<double> _weight1S0;
/**
* The weights for the \f$\phantom{1}^3S_1\f$ multiplets
*/
vector<double> _weight3S1;
/**
* The weights for the \f$\phantom{1}^1P_1\f$ multiplets
*/
vector<double> _weight1P1;
/**
* The weights for the \f$\phantom{1}^3P_0\f$ multiplets
*/
vector<double> _weight3P0;
/**
* The weights for the \f$\phantom{1}^3P_1\f$ multiplets
*/
vector<double> _weight3P1;
/**
* The weights for the \f$\phantom{1}^3P_2\f$ multiplets
*/
vector<double> _weight3P2;
/**
* The weights for the \f$\phantom{1}^1D_2\f$ multiplets
*/
vector<double> _weight1D2;
/**
* The weights for the \f$\phantom{1}^3D_1\f$ multiplets
*/
vector<double> _weight3D1;
/**
* The weights for the \f$\phantom{1}^3D_2\f$ multiplets
*/
vector<double> _weight3D2;
/**
* The weights for the \f$\phantom{1}^3D_3\f$ multiplets
*/
vector<double> _weight3D3;
//@}
/**
* Option for the construction of the tables
*/
unsigned int _topt;
/**
* Which particles to produce for debugging purposes
*/
unsigned int _trial;
/**
* @name A parameter used for determining when clusters are too light.
*
* This parameter is used for setting the lower threshold, \f$ t \f$ as
* \f[ t' = t(1 + r B^1_{\rm lim}) \f]
* where \f$ r \f$ is a random number [0,1].
*/
//@{
double _limBottom;
double _limCharm;
double _limExotic;
//@}
};
}
#endif /* Herwig_StandardModelHadronSpectrum_H */
diff --git a/src/defaults/Hadronization.in b/src/defaults/Hadronization.in
--- a/src/defaults/Hadronization.in
+++ b/src/defaults/Hadronization.in
@@ -1,143 +1,143 @@
# -*- ThePEG-repository -*-
############################################################
# Setup of default hadronization
#
# There are no user servicable parts inside.
#
# Anything that follows below should only be touched if you
# know what you're doing.
#############################################################
cd /Herwig/Particles
create ThePEG::ParticleData Cluster
setup Cluster 81 Cluster 0.00990 0.0 0.0 0.0 0 0 0 1
create ThePEG::ParticleData Remnant
setup Remnant 82 Remnant 0.00990 0.0 0.0 0.0 0 0 0 1
mkdir /Herwig/Hadronization
cd /Herwig/Hadronization
create Herwig::ClusterHadronizationHandler ClusterHadHandler
create Herwig::PartonSplitter PartonSplitter
create Herwig::ClusterFinder ClusterFinder
create Herwig::ColourReconnector ColourReconnector
create Herwig::ClusterFissioner ClusterFissioner
create Herwig::LightClusterDecayer LightClusterDecayer
create Herwig::ClusterDecayer ClusterDecayer
create Herwig::HwppSelector SMHadronSpectrum
newdef ClusterHadHandler:PartonSplitter PartonSplitter
newdef ClusterHadHandler:ClusterFinder ClusterFinder
newdef ClusterHadHandler:ColourReconnector ColourReconnector
newdef ClusterHadHandler:ClusterFissioner ClusterFissioner
newdef ClusterHadHandler:LightClusterDecayer LightClusterDecayer
newdef ClusterHadHandler:ClusterDecayer ClusterDecayer
do ClusterHadHandler:UseHandlersForInteraction QCD
newdef ClusterHadHandler:MinVirtuality2 0.1*GeV2
newdef ClusterHadHandler:MaxDisplacement 1.0e-10*millimeter
newdef ClusterHadHandler:UnderlyingEventHandler NULL
newdef PartonSplitter:HadronSpectrum SMHadronSpectrum
newdef ClusterFinder:HadronSpectrum SMHadronSpectrum
newdef ClusterFissioner:HadronSpectrum SMHadronSpectrum
newdef ClusterDecayer:HadronSpectrum SMHadronSpectrum
newdef LightClusterDecayer:HadronSpectrum SMHadronSpectrum
# ColourReconnector Default Parameters
newdef ColourReconnector:HadronSpectrum SMHadronSpectrum
newdef ColourReconnector:ColourReconnection Yes
newdef ColourReconnector:Algorithm Baryonic
# Statistical CR Parameters:
newdef ColourReconnector:AnnealingFactor 0.9
newdef ColourReconnector:AnnealingSteps 50
newdef ColourReconnector:TriesPerStepFactor 5.0
newdef ColourReconnector:InitialTemperature 0.1
# Plain and Baryonic CR Paramters
newdef ColourReconnector:ReconnectionProbability 0.95
newdef ColourReconnector:ReconnectionProbabilityBaryonic 0.7
# BaryonicMesonic and BaryonicMesonic CR Paramters
newdef ColourReconnector:ReconnectionProbability3Mto3M 0.5
newdef ColourReconnector:ReconnectionProbability3MtoBBbar 0.5
newdef ColourReconnector:ReconnectionProbabilityBbarBto3M 0.5
newdef ColourReconnector:ReconnectionProbability2Bto2B 0.05
newdef ColourReconnector:ReconnectionProbabilityMBtoMB 0.5
newdef ColourReconnector:StepFactor 1.0
newdef ColourReconnector:MesonToBaryonFactor 1.333
# General Parameters and switches
newdef ColourReconnector:MaxDistance 1.0e50
newdef ColourReconnector:OctetTreatment All
newdef ColourReconnector:CR2BeamClusters No
newdef ColourReconnector:Junction Yes
newdef ColourReconnector:LocalCR No
newdef ColourReconnector:CausalCR No
# Debugging
newdef ColourReconnector:Debug No
# set ClusterFissioner parameters
-set /Herwig/Hadronization/ClusterFissioner:KinematicThreshold Dynamic
-set /Herwig/Hadronization/ClusterFissioner:ProbablityPowerFactor 3.3462
-set /Herwig/Hadronization/ClusterFissioner:ProbablityShift 6.0442
-set /Herwig/Hadronization/ClusterFissioner:KineticThresholdShift -1.5321
+newdef ClusterFissioner:KinematicThreshold Dynamic
+newdef ClusterFissioner:ProbablityPowerFactor 3.3462
+newdef ClusterFissioner:ProbablityShift 6.0442
+newdef ClusterFissioner:KineticThresholdShift -1.5321
# Clustering parameters for light quarks
newdef ClusterFissioner:ClMaxLight 3.2344
newdef ClusterFissioner:ClPowLight 2.6464
newdef ClusterFissioner:PSplitLight 0.7235
newdef ClusterFissioner:Fission Default
insert ClusterFissioner:FissionPwt 1 1.0
insert ClusterFissioner:FissionPwt 2 1.0
insert ClusterFissioner:FissionPwt 3 0.35743
newdef ClusterDecayer:ClDirLight 1
newdef ClusterDecayer:ClSmrLight 0.78
# Clustering parameters for b-quarks
insert ClusterFissioner:ClMaxHeavy 5 3.757*GeV
insert ClusterFissioner:ClPowHeavy 5 0.547*GeV
insert ClusterFissioner:PSplitHeavy 5 0.625*GeV
insert ClusterDecayer:ClDirHeavy 5 1
insert ClusterDecayer:ClSmrHeavy 5 0.078
newdef SMHadronSpectrum:SingleHadronLimitBottom 0.000
# Clustering parameters for c-quarks
insert ClusterFissioner:ClMaxHeavy 4 3.950
insert ClusterFissioner:ClPowHeavy 4 2.559
insert ClusterFissioner:PSplitHeavy 4 0.994
insert ClusterDecayer:ClDirHeavy 4 1
insert ClusterDecayer:ClSmrHeavy 4 0.163
newdef SMHadronSpectrum:SingleHadronLimitCharm 0.000
# Clustering parameters for exotic quarks
# (e.g. hadronizing Susy particles)
newdef ClusterFissioner:ClMaxExotic 2.7*GeV
newdef ClusterFissioner:ClPowExotic 1.46
newdef ClusterFissioner:PSplitExotic 1.00
newdef ClusterDecayer:ClDirExotic 1
newdef ClusterDecayer:ClSmrExotic 0.
newdef SMHadronSpectrum:SingleHadronLimitExotic 0.
#
insert PartonSplitter:SplitPwt 1 1.0
insert PartonSplitter:SplitPwt 2 1.0
insert PartonSplitter:SplitPwt 3 0.824135
newdef PartonSplitter:Split Light
#
newdef SMHadronSpectrum:PwtDquark 1.0
newdef SMHadronSpectrum:PwtUquark 1.0
newdef SMHadronSpectrum:PwtSquark 0.35743
newdef SMHadronSpectrum:PwtCquark 0.0
newdef SMHadronSpectrum:PwtBquark 0.0
newdef SMHadronSpectrum:PwtDIquark 0.36452
newdef SMHadronSpectrum:SngWt 0.88022
newdef SMHadronSpectrum:DecWt 0.34622
newdef SMHadronSpectrum:Mode 1
newdef SMHadronSpectrum:BelowThreshold All
create Herwig::SpinHadronizer SpinHadronizer