diff --git a/Hadronization/ClusterFissioner.cc b/Hadronization/ClusterFissioner.cc
--- a/Hadronization/ClusterFissioner.cc
+++ b/Hadronization/ClusterFissioner.cc
@@ -1,2752 +1,2780 @@
// -*- 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),
_powerLawPower(-2.0),
_maxLoopFissionMatrixElement(5000000),
_safetyFactorMatrixElement(10.0),
_writeOut(0),
_hadronizingStrangeDiquarks(2)
{
}
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
+ os << ounit(_clMaxLight,GeV) << ounit(_clMaxHeavy,GeV) << ounit(_clMaxDiquark,GeV) << ounit(_clMaxExotic,GeV)
+ << _clPowLight << _clPowHeavy << _clPowDiquark << _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
<< _powerLawPower
<< _maxLoopFissionMatrixElement
<< _safetyFactorMatrixElement
<< _writeOut
;
}
void ClusterFissioner::persistentInput(PersistentIStream & is, int) {
- is >> iunit(_clMaxLight,GeV)
- >> iunit(_clMaxHeavy,GeV) >> iunit(_clMaxExotic,GeV) >> _clPowLight >> _clPowHeavy
- >> _clPowExotic >> _pSplitLight
- >> _pSplitHeavy >> _pSplitExotic
+ is >> iunit(_clMaxLight,GeV) >> iunit(_clMaxHeavy,GeV) >> iunit(_clMaxDiquark,GeV) >> iunit(_clMaxExotic,GeV)
+ >> _clPowLight >> _clPowHeavy >> _clPowDiquark >> _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
>> _powerLawPower
>> _maxLoopFissionMatrixElement
>> _safetyFactorMatrixElement
>> _writeOut
;
}
void ClusterFissioner::doinit() {
Interfaced::doinit();
// TODO: Some User warnings/errors but not complete list
if (_covariantBoost && _phaseSpaceSampler==0) throw Exception() << "Cannot use ClusterFissioner:CovariantBoost = Yes and ClusterFissioner:PhaseSpaceSampler = FullyAligned\nNot implemented yet\n" << Exception::runerror;
if (_matrixElement!=0 && _fissionApproach==0) generator()->logWarning( Exception(
"For non-trivial MatrixElement you need to enable FissionApproach=New or Hybrid\n",
Exception::warning));
if (_matrixElement==1 && !(_phaseSpaceSampler==2 && _massSampler==1 ) ) generator()->logWarning( Exception(
"The chosen ClusterFissioner:MatrixElement is only taken into account properly by using:\nMassSampler = Uniform & PhaseSpaceSampler = FullyIsotropic\n",
Exception::warning));
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 ( const long& id : spectrum()->lightHadronizingDiquarks() ) {
+ if (
+ _clPowDiquark.find(id) == _clPowDiquark.end() ||
+ _clMaxDiquark.find(id) == _clMaxDiquark.end() )
+ throw InitException() << "not all parameters have been set for diquarks 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;
if (_hadronizingStrangeDiquarks>0) {
_fissionPwt[3101] = 0.5 * _pwtDIquark * pwtSquark * pwtDquark;
_fissionPwt[3201] = 0.5 * _pwtDIquark * pwtSquark * pwtUquark;
if (_hadronizingStrangeDiquarks==2) {
_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 quarks",
&ClusterFissioner::_clMaxHeavy, GeV, -1, 3.35*GeV, ZERO, 10.0*GeV,
false, false, Interface::upperlim);
+ static ParMap<ClusterFissioner,Energy> interfaceClMaxDiquark
+ ("ClMaxDiquark",
+ "ClMax for light hadronizing diquarks",
+ &ClusterFissioner::_clMaxDiquark, 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 ParMap<ClusterFissioner,double> interfaceClPowDiquark
+ ("ClPowDiquark",
+ "ClPow for light hadronizing diquarks",
+ &ClusterFissioner::_clPowDiquark, -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, 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);
static SwitchOption interfaceFissionApproachHybrid
(interfaceFissionApproach,
"Hybrid",
"New cluster fission which allows to choose MassSampler"
", PhaseSpaceSampler and MatrixElement, but uses Default"
" Approach for BeamClusters",
2);
// 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 (Recommended)",
2);
// 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 interfaceMatrixElementSoftQQbarFinalFinal
(interfaceMatrixElement,
"SoftQQbarFinalFinal",
"Choose Soft p1,p2->q1,q2,g*->q1,q2,q,qbar matrix element"
"NOTE: Here the matrix element depends on qi.q(bar)",
1);
static SwitchOption interfaceMatrixElementSoftQQbarInitialFinal
(interfaceMatrixElement,
"SoftQQbarInitialFinal",
"Choose Soft p1,p2->q1,q2,g*->q1,q2,q,qbar matrix element "
"NOTE: Here the matrix element depends on pi.q(bar)",
2);
// Technical Max loop parameter for New ClusterFission approach and Matrix Element
// Rejection sampling
static Parameter<ClusterFissioner,int> interfaceMaxLoopMatrixElement
("MaxLoopMatrixElement",
"Technical Parameter for how many tries are allowed to sample the "
"Cluster Fission matrix element before reverting to fissioning "
"using the default Fission Aproach",
&ClusterFissioner::_maxLoopFissionMatrixElement, 5000000, 100, 1e8,
false, false, Interface::limited);
static Switch<ClusterFissioner,int> interfaceDiquarkClusterFission
("DiquarkClusterFission",
"Allow clusters to fission to 1 or 2 diquark Clusters or Turn off diquark fission completely",
&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 SwitchOption interfaceDiquarkClusterFissionOff
(interfaceDiquarkClusterFission,
"Off",
"Don't allow clusters fission to draw diquarks ",
-1);
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);
// Matrix Element Choice Switch
static Parameter<ClusterFissioner,double>
interfaceSafetyFactorOverEst("SafetyFactorOverEst","Safety factor with which the numerical overestimate is calculated",
&ClusterFissioner::_safetyFactorMatrixElement, 0, 10.0, 1.0, 1000.0,false,false,false);
// Matrix Element Choice Switch
static Switch<ClusterFissioner,int> interfaceWriteOut
("WriteOut",
"Option for different Matrix Element options",
&ClusterFissioner::_writeOut, 0, false, false);
static SwitchOption interfaceWriteOutNo
(interfaceWriteOut,
"No",
"Choose trivial matrix element i.e. whatever comes from the mass and "
"phase space sampling",
0);
static SwitchOption interfaceWriteOutYes
(interfaceWriteOut,
"Yes",
"Choose Soft q1,q2->q1,q2,g*->q1,q2,q,qbar matrix element",
1);
// Allowing for strange diquarks in the ClusterFission
static Switch<ClusterFissioner,unsigned int> interfaceHadronizingStrangeDiquarks
("HadronizingStrangeDiquarks",
"Option for adding strange diquarks to Cluster Fission (if Fission = New or Hybrid is enabled)",
&ClusterFissioner::_hadronizingStrangeDiquarks, 0, false, false);
static SwitchOption interfaceHadronizingStrangeDiquarksNo
(interfaceHadronizingStrangeDiquarks,
"No",
"No strangeness containing diquarks during Cluster Fission",
0);
static SwitchOption interfaceHadronizingStrangeDiquarksOnlySingleStrange
(interfaceHadronizingStrangeDiquarks,
"OnlySingleStrange",
"Only one strangeness containing diquarks during Cluster Fission i.e. su,sd",
1);
static SwitchOption interfaceHadronizingStrangeDiquarksAll
(interfaceHadronizingStrangeDiquarks,
"All",
"All strangeness containing diquarks during Cluster Fission i.e. su,sd,ss",
2);
}
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
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;
case 2:
if (cluster->isBeamCluster()) {
return cutTwoDefault(cluster, finalhadrons, softUEisOn);
}
else {
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);
bool failure = drawNewMassesDefault(Mc, soft1, soft2, pClu1, pClu2,
ptrQ1, pQ1, newPtr1, pQone,
newPtr2, pQtwo, ptrQ2, pQ2);
if (failure)
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) 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 && !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();
// TODO BEGIN Changed comp to default
if ( Mc < spectrum()->massLightestHadronPair(ptrQ1->dataPtr(),ptrQ2->dataPtr())) {
static const PPtr null = PPtr();
return cutType(PPair(null,null),PPair(null,null));
}
// TODO END Changed comp to default
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.
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();
// Minimal threshold for non-zero Mass PhaseSpace
if ( Mc < (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: ###
bool escape = false;
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 - mMin)/GeV) < 1e-3) {
escape = true;
retTo = Done;
}
continue;
}
pQ1.setMass(m1);
pQone.setMass(m);
pQtwo.setMass(m);
pQ2.setMass(m2);
// 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();
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;
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
) {
escape = true;
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
) {
escape = true;
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
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
//
// Check if Mc1max < MLHP1, in which case we might need to choose a different flavour
// else resampling the masses should be sufficient
Energy MLHP1 = _hadronSpectrum->massLightestHadronPair(ptrQ1->dataPtr(), newPtr1->dataPtr());
Energy MLHP2 = _hadronSpectrum->massLightestHadronPair(ptrQ2->dataPtr(), newPtr2->dataPtr());
Energy Mc1max = (m2+m) > MLHP2 ? (Mc-(m2+m)):(Mc-MLHP2);
// for avoiding inf loops set min threshold
if ( Mc1max - MLHP1 < 1e-2*GeV ) {
retTo = FlavourSampling;
}
else {
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
//
// Check if Mc2max < MLHP2, in which case we might need to choose a different flavour
// else resampling the masses should be sufficient
Energy MLHP1 = _hadronSpectrum->massLightestHadronPair(ptrQ1->dataPtr(), newPtr1->dataPtr());
Energy MLHP2 = _hadronSpectrum->massLightestHadronPair(ptrQ2->dataPtr(), newPtr2->dataPtr());
Energy Mc2max = (m1+m) > MLHP1 ? (Mc-(m1+m)):(Mc-MLHP1);
// for avoiding inf loops set min threshold
if ( Mc2max - MLHP2 < 1e-2*GeV ) {
retTo = FlavourSampling;
}
else {
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
) {
// escape = true;
retTo = Done;
}
else {
// Try again with lighter quark
retTo = FlavourSampling;
}
continue;
}
// 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;
}
// Actual MatrixElement evaluated at sampled PhaseSpace point
double SQME = calculateSQME(p0Q1,p0Q2,pQ1,pQone,pQ2,pQtwo);
// Total overestimate of MatrixElement independent
// of the PhaseSpace point and independent of M1,M2
double SQMEoverEstimate = calculateSQME_OverEstimate(Mc,m1,m2,m);
// weight for MatrixElement*PhaseSpace must be in [0:1]
double weightSQME = SQME/SQMEoverEstimate;
assert(weightSQME>0.0);
// weight(M1,M2) for M1*M2*(Two body PhaseSpace)**3 should be in [0,1]
double weightFlatPS = weightFlatPhaseSpace(Mc, Mc1, Mc2, m, m1, m2);
if (weightFlatPS<0 || weightFlatPS>1.0){
throw Exception() << "weightFlatPS = "<< weightFlatPS << " > 1 or negative in ClusterFissioner::cutTwo"
<< "Mc = " << Mc/GeV
<< "Mc1 = " << Mc1/GeV
<< "Mc2 = " << Mc2/GeV
<< "m1 = " << m1/GeV
<< "m2 = " << m2/GeV
<< "m = " << m/GeV
<< "SQME = " << SQME
<< "SQME_OE = " << SQMEoverEstimate
<< "PS = " << weightFlatPS
<< Exception::runerror;
}
// current phase space point is distributed according to weightSamp
double Pacc = weightFlatPS * weightSQME;
if (!(Pacc >= 0 ) || std::isnan(Pacc) || std::isinf(Pacc)){
throw Exception() << "Pacc = "<< Pacc << " y < 0 in ClusterFissioner::cutTwo"
<< "Mc = " << Mc/GeV
<< "Mc1 = " << Mc1/GeV
<< "Mc2 = " << Mc2/GeV
<< "m1 = " << m1/GeV
<< "m2 = " << m2/GeV
<< "m = " << m/GeV
<< "SQME = " << SQME
<< "SQME_OE = " << SQMEoverEstimate
<< "PS = " << weightFlatPS
<< Exception::runerror;
}
assert(Pacc >= 0.0);
if (Pacc > 1.0){
throw Exception() << "Pacc = "<< Pacc << " > 1 in ClusterFissioner::cutTwo"
<< "Mc = " << Mc/GeV
<< "Mc1 = " << Mc1/GeV
<< "Mc2 = " << Mc2/GeV
<< "m1 = " << m1/GeV
<< "m2 = " << m2/GeV
<< "m = " << m/GeV
<< Exception::warning;
}
static int first=_writeOut;
if (UseRandom::rnd()<Pacc) {
if (_writeOut && _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 (first){
first=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 && !escape && counter_MEtry < _maxLoopFissionMatrixElement);
if(escape) {
// happens if we get at too light cluster to begin with
static const PPtr null = PPtr();
return cutType(PPair(null,null),PPair(null,null));
}
if(counter_MEtry >= _maxLoopFissionMatrixElement) {
// happens if we get too massive clusters where the matrix element
// is very inefficiently sampled
std::stringstream warning;
warning
<< "Matrix Element rejection sampling tried more than "
<< counter_MEtry
<< " times.\nMcluster = " << Mc/GeV
<< "GeV\nm1 = " << m1/GeV
<< "GeV\nm2 = " << m2/GeV
<< "GeV\nm = " << m/GeV
<< "GeV\n"
<< "isBeamCluster = " << cluster->isBeamCluster()
<<"Using default as ClusterFissioner::cutTwoDefault as a fallback.\n"
<< "***This Exception should not happen too often!*** ";
generator()->logWarning( Exception(warning.str(),Exception::warning));
// throw Exception() << warning.str()
// << Exception::eventerror;
return cutTwoDefault(cluster,finalhadrons,softUEisOn);
}
// ==> 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 = drawNewMassesDefault(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 && !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();
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 the phase space weight for M1*M2*(2 body PhaseSpace)^3
*/
double ClusterFissioner::weightFlatPhaseSpace(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;
// old weight of Jan without the Jacobi factor M1*M2 of the Mass integration
// double PSweight = pow(lam1*lam2*lam3,_dim-3.)*pow(M1_temp*M2_temp,2.-_dim);
// new weight with the Jacobi factor M1*M2 of the Mass integration
double PSweight = pow(lam1*lam2*lam3,_dim-3.)*pow(M1_temp*M2_temp,3.-_dim);
// overestimate only possible for dim>=3.0
assert(_dim>=3.0);
// old overestimate of Jan without the Jacobi factor M1*M2 of the Mass integration
// 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);
// new improved overestimate with the Jacobi factor M1*M2 of the Mass integration
double overEstimate = pow(6.0*sqrt(3.0), 3.0 - _dim)*pow(M_temp, 4.*_dim-12.);
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";
// if (ratio > 0.9) std::cout << "ratio = " << ratio <<"\n";
assert (ratio >= 0);
assert (ratio <= 1);
return ratio;
}
/**
* Veto for the phase space weight
* returns true if proposed Masses are rejected
* else returns false
*/
bool ClusterFissioner::phaseSpaceVeto(const Energy Mc, const Energy Mc1, const Energy Mc2,
const Energy m, const Energy m1, const Energy m2) const {
return (UseRandom::rnd()>weightFlatPhaseSpace(Mc, Mc1, Mc2, m, m1, m2));
}
/**
* 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 add precise weightMS that could be used used for improving the rejection Sampling
switch (_massSampler)
{
case 0:
return drawNewMassesDefault(Mc, soft1, soft2, pClu1, pClu2, ptrQ1, pQ1, tcPPtr(), pQone, tcPPtr(),pQtwo, ptrQ2, pQ2);
break;
case 1:
return drawNewMassesUniform(Mc, pClu1, pClu2, pQ1, pQone, pQ2);
break;
case 2:
return drawNewMassesPhaseSpace(Mc, 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,
tcPPtr, const Lorentz5Momentum& pQone,
tcPPtr, 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, 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,
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
}
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::sampleDirectionCluster(
const Lorentz5Momentum & pQ,
const Lorentz5Momentum & pClu) const {
switch (_phaseSpaceSampler)
{
case 0:
// Default aligned
if (_covariantBoost)
// in Covariant Boost the positive z-Axis is defined as the direction of
// the pQ vector in the Cluster rest frame
return Axis(0,0,1);
else
return sampleDirectionAligned(pQ, pClu);
case 1:
// Default aligned
if (_covariantBoost)
// in Covariant Boost the positive z-Axis is defined as the direction of
// the pQ vector in the Cluster rest frame
return Axis(0,0,1);
else
return sampleDirectionAligned(pQ, pClu);
case 2:
// Isotropic
return sampleDirectionIsotropic();
default:
assert(false);
}
}
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;
}
case 2:
{
Energy2 p1p2 = p1 * p2;
Energy2 p1q = p1 * q ;
Energy2 p2qbar = p2 * qbar;
Energy2 p2q = p2 * q ;
Energy2 p1qbar = p1 * qbar;
Energy2 qqbar = q * qbar;
Energy2 mp2 = q.mass2();
double Numerator = p1p2 * (qqbar + mp2)/sqr(GeV2);
Numerator += 0.5 * (p1q - p1qbar)*(p2q - p2qbar)/sqr(GeV2);
double Denominator = sqr(qqbar + mp2)*(p1q + p1qbar)*(p2q + p2qbar)/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 where Mass -> m1
* p0Q2 -> p2 where Mass -> m2
* pQ1 -> q1 where Mass -> m1
* pQone -> q where Mass -> mq
* pQ2 -> q2 where Mass -> m2
* pQone -> qbar where Mass -> mq
* */
double ClusterFissioner::calculateSQME_OverEstimate(
const Energy&,
const Energy&,
const Energy&,
const Energy& mq
) const {
double SQME_OverEstimate;
switch (_matrixElement)
{
case 0:
SQME_OverEstimate = 1.0;
break;
case 1:
{
// Fit factor for guess of best overestimate
double A = 0.25;
SQME_OverEstimate = _safetyFactorMatrixElement*A*pow(mq/GeV,-4);
break;
}
case 2:
{
// Fit factor for guess of best overestimate
double A = 0.25;
SQME_OverEstimate = _safetyFactorMatrixElement*A*pow(mq/GeV,-4);
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 {
// calculateKinematics(pClu,p0Q1,toHadron1,toHadron2,
// pClu1,pClu2,pQ1,pQbar,pQ,pQ2bar);
// return false;
if (pClu.m() < pClu1.mass() + pClu2.mass()
|| pClu1.mass()<ZERO
|| pClu2.mass()<ZERO ) {
throw Exception() << "Impossible Kinematics in ClusterFissioner::drawKinematics() (A)\n"
<< "Mc = "<< pClu.mass()/GeV <<" GeV\n"
<< "Mc1 = "<< pClu1.mass()/GeV <<" GeV\n"
<< "Mc2 = "<< pClu2.mass()/GeV <<" GeV\n"
<< 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():DirToClu;
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():DirToClu;
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)
- );
+ // particle data for constituents
+ tcPDPtr cptr[3]={tcPDPtr(),tcPDPtr(),tcPDPtr()};
+ long heaviestHadronizingDiquark=0;
+ for(size_t ix=0;ix<min(clu->numComponents(),3);++ix) {
+ cptr[ix]=clu->particle(ix)->dataPtr();
+ // Assuming diquark masses are ordered with larger id corresponding to larger masses
+ if (DiquarkMatcher::Check(*(clu->particle(ix)->dataPtr()))
+ && abs(clu->particle(ix)->dataPtr()->id()) > heaviestHadronizingDiquark
+ && abs(clu->particle(ix)->dataPtr()->id()) < getParticleData(ParticleID::ss_1)->id()) {
+ heaviestHadronizingDiquark=abs(clu->particle(ix)->dataPtr()->id());
+ }
+ }
+ // 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;
+ // if no heavy quark is found in the cluster, but diquarks are present use
+ // different ClMax and ClPow
+ if ( heaviestHadronizingDiquark) {
+ clpow = _clPowDiquark.find(heaviestHadronizingDiquark) == _clPowDiquark.end() ? clpow:_clPowDiquark[heaviestHadronizingDiquark];
+ clmax = _clMaxDiquark.find(heaviestHadronizingDiquark) == _clMaxDiquark.end() ? clmax:_clMaxDiquark[heaviestHadronizingDiquark];
+ }
- 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);
+ 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)
+ );
- scale = clu->mass()/GeV;
- threshold = clu->sumConstituentMasses()/GeV + 2.0 * minmass/GeV;
+ 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);
- canSplitMinimally = ProbablityFunction(scale,threshold);
- }
+ scale = clu->mass()/GeV;
+ threshold = clu->sumConstituentMasses()/GeV + 2.0 * minmass/GeV;
- return aboveCutoff && canSplitMinimally;
+ canSplitMinimally = ProbablityFunction(scale,threshold);
+ }
+
+ return aboveCutoff && canSplitMinimally;
}
diff --git a/Hadronization/ClusterFissioner.h b/Hadronization/ClusterFissioner.h
--- a/Hadronization/ClusterFissioner.h
+++ b/Hadronization/ClusterFissioner.h
@@ -1,753 +1,755 @@
// -*- C++ -*-
//
// ClusterFissioner.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_ClusterFissioner_H
#define HERWIG_ClusterFissioner_H
#include <ThePEG/Interface/Interfaced.h>
#include "CluHadConfig.h"
#include "ClusterFissioner.fh"
#include "HadronSpectrum.h"
namespace Herwig {
using namespace ThePEG;
//class Cluster; // forward declaration
/** \ingroup Hadronization
* \class ClusterFissioner
* \brief This class handles clusters which are too heavy.
* \author Philip Stephens
* \author Alberto Ribon
* \author Stefan Gieseke
*
* This class does the job of chopping up either heavy clusters or beam
* clusters in two lighter ones. The procedure is repeated recursively until
* all of the cluster children have masses below some threshold values.
*
* For the beam remnant clusters, at the moment what is done is the following.
* In the case that the soft underlying event is switched on, the
* beam remnant clusters are tagged as not available,
* therefore they will not be treated at all during the hadronization.
* In the case instead that the soft underlying event is switched off,
* then the beam remnant clusters are treated exactly as "normal" clusters,
* with the only exception of the mass spectrum used to generate the
* cluster children masses. For non-beam clusters, the masses of the cluster
* children are draw from a power-like mass distribution; for beam clusters,
* according to the value of the flag _IOpRem, either both
* children masses are draw from a fast-decreasing exponential mass
* distribution (case _IOpRem == 0, or, indendently by
* _IOpRem, in the special case that the beam cluster contains two
* beam remnants), or one mass from the exponential distribution (corresponding
* of the cluster child with the beam remnant) and the other with the usual
* power-like distribution (case _IOpRem == 1, which is the
* default one, as in Herwig 6.3).
*
* The reason behind the use of a fast-decreasing exponential distribution
* is that to avoid a large transverse energy from the many sequential
* fissions that would otherwise occur due to the typical large cluster
* mass of beam clusters. Using instead an exponential distribution
* the masses of the two cluster children will be very small (order of
* GeV).
*
* The rationale behind the implementation of the splitting of clusters
* has been to preserve *all* of the information about such splitting
* process. More explicitly a ThePEG::Step class is passed in and the
* new clusters are added to the step as the decay products of the
* heavy cluster. This approach has the twofold
* advantage to provide all of the information that could be needed
* (expecially in future developments), without any information loss,
* and furthermore it allows a better debugging.
*
* @see \ref ClusterFissionerInterfaces "The interfaces"
* defined for ClusterFissioner.
*/
class ClusterFissioner: public Interfaced {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* Default constructor.
*/
ClusterFissioner();
//@}
/** Splits the clusters which are too heavy.
*
* Split either heavy clusters or beam clusters recursively until all
* children have mass below some threshold. Heavy clusters are those that
* satisfy the condition
* \f[ M^P > C^P + S^P \f]
* where \f$ M \f$ is the clusters mass, \f$ P \f$ is the parameter
* ClPow, \f$ C \f$ is the parameter ClMax and \f$ S \f$ is the
* sum of the clusters constituent partons.
* For beam clusters, they are split only if the soft underlying event
* is switched off, otherwise these clusters will be tagged as unavailable
* and they will not be treated by the hadronization altogether.
* In the case beam clusters will be split, the procedure is exactly
* the same as for normal non-beam clusters, with the only exception
* of the mass spectrum from which to draw the masses of the two
* cluster children (see method drawChildrenMasses for details).
*/
tPVector fission(ClusterVector & clusters, bool softUEisOn);
/**
* Return the hadron spectrum
*/
Ptr<HadronSpectrum>::tptr spectrum() const {
return _hadronSpectrum;
}
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* Standard Init function used to initialize the interfaces.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
private:
/**
* Private and non-existent assignment operator.
*/
ClusterFissioner & operator=(const ClusterFissioner &) = delete;
/**
* This method directs the splitting of the heavy clusters
*
* This method does the splitting of the clusters and all of its cluster
* children, if heavy. All of these new children clusters are added to the
* collection of clusters. The method works as follows.
* Initially the vector contains just the stack of input pointers to the
* clusters to be split. Then it will be filled recursively by all
* of the cluster's children that are heavy enough to require
* to be split. In each loop, the last element of the vector is
* considered (only once because it is then removed from the vector).
*
* \todo is the following still true?
* For normal, non-beam clusters, a power-like mass distribution
* is used, whereas for beam clusters a fast-decreasing exponential mass
* distribution is used instead. This avoids many iterative splitting which
* could produce an unphysical large transverse energy from a supposed
* soft beam remnant process.
*/
void cut(stack<ClusterPtr> &,
ClusterVector&, tPVector & finalhadrons, bool softUEisOn);
public:
/**
* Definition for easy passing of two particles.
*/
typedef pair<PPtr,PPtr> PPair;
/**
* Definition for use in the cut function.
*/
typedef pair<PPair,PPair> cutType;
/**
* Splits the input cluster.
*
* Split the input cluster (which can be either an heavy non-beam
* cluster or a beam cluster). The result is two pairs of particles. The
* first element of each pair is new cluster/hadron, while the second
* element of each pair is the particle drawn from the vacuum to create
* the new cluster/hadron.
* Notice that this method treats also beam clusters by using a different
* mass spectrum used to generate the cluster child masses (see method
* drawChildMass).
*/
//@{
/**
* Split two-component cluster
*/
virtual cutType cutTwo(ClusterPtr &, tPVector & finalhadrons, bool softUEisOn);
/**
* Split two-component cluster using Default FissionApproach
*/
virtual cutType cutTwoDefault(ClusterPtr &, tPVector & finalhadrons, bool softUEisOn);
/**
* Split two-component cluster using New FissionApproach
*/
virtual cutType cutTwoNew(ClusterPtr &, tPVector & finalhadrons, bool softUEisOn);
/**
* Split three-component cluster
*/
virtual cutType cutThree(ClusterPtr &, tPVector & finalhadrons, bool softUEisOn);
//@}
public:
/**
* Produces a hadron and returns the flavour drawn from the vacuum.
*
* This routine produces a new hadron. It
* also sets the momentum and vertex to the values given.
*/
PPair produceHadron(tcPDPtr hadron, tPPtr newPtr, const Lorentz5Momentum &a,
const LorentzPoint &b) const;
protected:
/**
* Produces a cluster from the flavours passed in.
*
* This routine produces a new cluster with the flavours given by ptrQ and newPtr.
* The new 5 momentum is a and the parent momentum are c and d. C is for the
* ptrQ and d is for the new particle newPtr. rem specifies whether the existing
* particle is a beam remnant or not.
*/
PPair produceCluster(tPPtr ptrQ, tPPtr newPtr, const Lorentz5Momentum &a,
const LorentzPoint &b, const Lorentz5Momentum &c,
const Lorentz5Momentum &d, const bool rem,
tPPtr spect=tPPtr(), bool remSpect=false) const;
/**
* Returns the new quark-antiquark pair
* needed for fission of a heavy cluster. Equal probabilities
* are assumed for producing u, d, or s pairs.
*/
void drawNewFlavourQuarks(PPtr& newPtrPos,PPtr& newPtrNeg) const;
/**
* Returns the new quark-antiquark pair or diquark -
* antidiquark pair needed for fission of a heavy cluster.
*/
void drawNewFlavourDiquarks(PPtr& newPtrPos,PPtr& newPtrNeg,
const ClusterPtr & clu) const;
/**
* Returns the new quark-antiquark pair
* needed for fission of a heavy cluster. Equal probabilities
* are assumed for producing u, d, or s pairs.
* Extra argument is used when performing strangeness enhancement
*/
void drawNewFlavourEnhanced(PPtr& newPtrPos,PPtr& newPtrNeg, Energy2 mass2) const;
/**
* Produces the mass of a child cluster.
*
* Draw the masses \f$M'\f$ of the the cluster child produced
* by the fission of an heavy cluster (of mass M). m1, m2 are the masses
* of the constituents of the cluster; m is the mass of the quark extract
* from the vacuum (together with its antiparticle). The algorithm produces
* the mass of the cluster formed with consituent m1.
* Two mass distributions can be used for the child cluster mass:
* -# power-like mass distribution ("normal" mass) with power exp
* \f[ M' = {\rm rnd}((M-m_1-m_2-m)^P, m^p)^{1/P} + m_1 \f]
* where \f$ P \f$ is a parameter of the model and \f$ \rm{rnd} \f$ is
* the function:
* \f[ \rm{rnd}(a,b) = (1-r)a + r b \f]
* and here \f$ r \f$ is a random number [0,1].
* -# fast-decreasing exponential mass distribution ("soft" mass) with
* rmin. rmin is given by
* \f[ r_{\rm min} = \exp(-b (M - m_1 - m_2 - 2 m)) \f]
* where \f$ b \f$ is a parameter of the model. The generated mass is
* given by
* \f[ M' = m_1 + m - \frac{\log\left(
* {\rm rnd}(r_{\rm min}, 1-r_{\rm min})\right)}{b} \f].
*
* The choice of which mass distribution should be used for each of the two
* cluster children is dictated by the parameter soft.
*/
Energy drawChildMass(const Energy M, const Energy m1, const Energy m2,
const Energy m, const double exp, const bool soft) const;
/**
* Determine the positions of the two children clusters.
*
* This routine generates the momentum of the decay products. It also
* generates the momentum in the lab frame of the partons drawn out of
* the vacuum.
*/
void calculatePositions(const Lorentz5Momentum &pClu,
const LorentzPoint & positionClu,
const Lorentz5Momentum & pClu1,
const Lorentz5Momentum & pClu2,
LorentzPoint & positionClu1,
LorentzPoint & positionClu2 ) const;
protected:
/**
* Dimension used to calculate phase space weights
*/
double dim() const {return _dim;}
/**
* Access to soft-cluster parameter
*/
Energy btClM() const {return _btClM;}
/**
* Function that returns either the cluster mass or the lambda measure
*/
Energy2 clustermass(const ClusterPtr & cluster) const;
/**
* old calculateKinematics function for Default FissionApproach
*/
void 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;
/**
* Draw a new flavour for the given cluster; currently defaults to
* the default model
*/
virtual void drawNewFlavour(PPtr& newPtr1, PPtr& newPtr2, const ClusterPtr & cluster) const {
if (_enhanceSProb == 0){
if (_diquarkClusterFission>=0) drawNewFlavourDiquarks(newPtr1,newPtr2,cluster);
else drawNewFlavourQuarks(newPtr1,newPtr2);
}
else {
drawNewFlavourEnhanced(newPtr1,newPtr2,clustermass(cluster));
}
}
/**
* 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 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;
/**
* 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 drawNewMassesDefault(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;
/**
* Sample the masses for flat phase space
* */
bool drawNewMassesUniform(const Energy Mc,
Lorentz5Momentum & pClu1, Lorentz5Momentum & pClu2,
const Lorentz5Momentum & pQ1,
const Lorentz5Momentum & pQ,
const Lorentz5Momentum & pQ2) const;
/**
* Sample the masses for flat phase space
* */
bool drawNewMassesPhaseSpace(const Energy Mc,
Lorentz5Momentum & pClu1, Lorentz5Momentum & pClu2,
const Lorentz5Momentum & pQ1,
const Lorentz5Momentum & pQ,
const Lorentz5Momentum & pQ2) const;
/**
* Calculate the final kinematics of a heavy cluster decay C->C1 +
* C2, if not already performed by drawNewMasses
* @return returns false if failes
*/
bool drawKinematics(const Lorentz5Momentum &pClu,
const Lorentz5Momentum &p0Q1,
const Lorentz5Momentum &p0Q2,
const bool toHadron1, const bool toHadron2,
Lorentz5Momentum &pClu1, Lorentz5Momentum &pClu2,
Lorentz5Momentum &pQ1, Lorentz5Momentum &pQb,
Lorentz5Momentum &pQ, Lorentz5Momentum &pQ2b) const;
bool drawKinematicsFullyAligned(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;
bool drawKinematicsAlignedIsotropic(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;
bool drawKinematicsFullyIsotropic(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;
/**
* Calculation of the squared matrix element from the drawn
* momenta
* @return value of the squared matrix element in units of
* 1/GeV4
* */
double calculateSQME(
const Lorentz5Momentum & p1,
const Lorentz5Momentum & p2,
const Lorentz5Momentum & q1,
const Lorentz5Momentum & q,
const Lorentz5Momentum & q2,
const Lorentz5Momentum & qbar) const;
/**
* Calculation of the overestimate for the squared matrix
* element independent on M1 and M2
* @return value of the overestimate squared matrix element
* in units of 1/GeV4
* */
double calculateSQME_OverEstimate(
const Energy& Mc,
const Energy& m1,
const Energy& m2,
const Energy& mq) const;
protected:
/** @name Access members for child classes. */
//@{
/**
* Access to the hadron selector
*/
HadronSpectrumPtr hadronSpectrum() const {return _hadronSpectrum;}
//@}
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
//@}
private:
Axis sampleDirectionCluster(const Lorentz5Momentum & pQ, const Lorentz5Momentum & pClu) const;
/**
* Samples the direction of Cluster Fission products uniformly
**/
Axis sampleDirectionIsotropic() const;
/**
* Samples the direction of Cluster Fission products uniformly
* but only accepts those flying in the direction of pRelCOM
**/
Axis sampleDirectionSemiUniform(const Lorentz5Momentum & pQ, const Lorentz5Momentum & pClu) const;
/**
* Samples the direction of Cluster Fission products according to Default
* fully aligned D = 1+1 Fission
* */
Axis sampleDirectionAligned(const Lorentz5Momentum & pQ, const Lorentz5Momentum & pClu) const;
/**
* Smooth probability for dynamic threshold cuts:
* @scale the current scale, e.g. the mass of the cluster,
* @threshold the physical threshold,
*/
bool ProbablityFunction(double scale, double threshold);
/**
* Check if a cluster is heavy enough to split again
*/
bool isHeavy(tcClusterPtr );
/**
* Check if a cluster is heavy enough to be at least kinematically able to split
*/
bool canSplitMinimally(tcClusterPtr, Energy);
/**
* Check if can't make a hadron from the partons
*/
inline bool cantMakeHadron(tcPPtr p1, tcPPtr p2) {
return ! spectrum()->canBeHadron(p1->dataPtr(), p2->dataPtr());
}
/**
* Flat PhaseSpace weight for ClusterFission
*/
double weightFlatPhaseSpace(const Energy Mc, const Energy Mc1, const Energy Mc2,
const Energy m, const Energy m1, const Energy m2) const;
/**
* Calculate a veto for the phase space weight
*/
bool phaseSpaceVeto(const Energy Mc, const Energy Mc1, const Energy Mc2,
const Energy m, const Energy m1, const Energy m2) const;
/**
* A pointer to a Herwig::HadronSpectrum object for generating hadrons.
*/
HadronSpectrumPtr _hadronSpectrum;
/**
* @name The Cluster max mass,dependant on which quarks are involved, used to determine when
* fission will occur.
*/
//@{
Energy _clMaxLight;
map<long,Energy> _clMaxHeavy;
+ map<long,Energy> _clMaxDiquark;
Energy _clMaxExotic;
//@}
/**
* @name The power used to determine when cluster fission will occur.
*/
//@{
double _clPowLight;
map<long,double> _clPowHeavy;
+ map<long,double> _clPowDiquark;
double _clPowExotic;
//@}
/**
* @name The power, dependant on whic quarks are involved, used in the cluster mass generation.
*/
//@{
double _pSplitLight;
map<long,double> _pSplitHeavy;
double _pSplitExotic;
/**
* Weights for alternative cluster fission
*/
map<long,double> _fissionPwt;
/**
* Include phase space weights
*/
bool _phaseSpaceWeights;
/**
* Dimensionality of phase space weight
*/
double _dim;
/**
* Flag used to determine between normal cluster fission and alternative cluster fission
*/
int _fissionCluster;
/**
* Flag to choose static or dynamic kinematic thresholds in cluster splittings
*/
int _kinematicThresholdChoice;
/**
* Pwt weight for drawing diquark
*/
double _pwtDIquark;
/**
* allow clusters to fission to 1 (or 2) diquark clusters or not
*/
int _diquarkClusterFission;
//@}
/**
* Parameter used (2/b) for the beam cluster mass generation.
* Currently hard coded value.
*/
Energy _btClM;
/**
* Flag used to determine what distributions to use for the cluster masses.
*/
int _iopRem;
/**
* The string constant
*/
Tension _kappa;
/**
* Width of the gaussian sampling for the FluxTube Kinematics
*/
double _fluxTubeWidth;
/**
* Flag that switches between no strangeness enhancement, scaling enhancement,
* and exponential enhancement (in numerical order)
*/
int _enhanceSProb;
/**
* Parameter that governs the strangeness enhancement scaling
*/
Energy _m0Fission;
/**
* Flag that switches between mass measures used in strangeness enhancement:
* cluster mass, or the lambda measure - ( m_{clu}^2 - (m_q + m_{qbar})^2 )
*/
int _massMeasure;
/**
* Constant variable which stops the scale from being to large, and not worth
* calculating
*/
const double _maxScale = 20.;
/**
* Power factor in ClausterFissioner bell probablity function
*/
double _probPowFactor;
/**
* Shifts from the center in ClausterFissioner bell probablity function
*/
double _probShift;
/**
* Shifts from the kinetic threshold in ClausterFissioner
*/
Energy2 _kinThresholdShift;
/**
* Flag for strict diquark selection according to kinematics
*/
int _strictDiquarkKinematics;
/**
* Use Covariant boost in ClusterFissioner
*/
bool _covariantBoost;
/**
* Flag for allowing Hadron Final states in Cluster Fission
*/
int _allowHadronFinalStates;
/**
* Choice of Mass sampling for ClusterFissioner
* i.e. rejection sampling starting from MassPresampling
* Chooses the to be sampled Mass distribution
* Note: This ideal distribution would be ideally
* exactly the integral over the angles as a
* function of M1,M2
*/
int _massSampler;
/**
* Choice of Phase Space sampling for ClusterFissioner
*/
int _phaseSpaceSampler;
/**
* Choice of Matrix Element for ClusterFissioner
* Note: The choice of the matrix element requires
* to provide one overestimate as a function
* of M1,M2 and another overestimate of the
* previous overestimate independent of
* M1 and M2 i.e. only dependent on M,m1,m2,m
*/
int _matrixElement;
/**
* Choice of ClusterFissioner Approach
*/
int _fissionApproach;
/**
* Power for MassPreSampler = PowerLaw
*/
double _powerLawPower;
/**
* Choice of ClusterFissioner Approach
* Technical Parameter for how many tries are allowed to sample the
* Cluster Fission matrix element before reverting to fissioning
* using the default Fission Aproach
*/
int _maxLoopFissionMatrixElement;
/*
* Safety factor for a better overestimate of the matrix Element
*
* */
double _safetyFactorMatrixElement;
int _writeOut;
/*
* flag for allowing strange Diquarks to be produced during
* Cluster Fission
* */
unsigned int _hadronizingStrangeDiquarks;
};
}
#endif /* HERWIG_ClusterFissioner_H */
diff --git a/Hadronization/LightClusterDecayer.cc b/Hadronization/LightClusterDecayer.cc
--- a/Hadronization/LightClusterDecayer.cc
+++ b/Hadronization/LightClusterDecayer.cc
@@ -1,405 +1,409 @@
// -*- C++ -*-
//
// LightClusterDecayer.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 LightClusterDecayer class.
//
#include "LightClusterDecayer.h"
#include <ThePEG/Interface/ClassDocumentation.h>
#include <ThePEG/Interface/Parameter.h>
#include <ThePEG/Interface/Reference.h>
#include <ThePEG/Interface/Switch.h>
#include <ThePEG/Persistency/PersistentOStream.h>
#include <ThePEG/Persistency/PersistentIStream.h>
#include <ThePEG/PDT/EnumParticles.h>
#include <ThePEG/Repository/EventGenerator.h>
#include "Cluster.h"
#include "Herwig/Utilities/Kinematics.h"
#include <ThePEG/Utilities/DescribeClass.h>
using namespace Herwig;
DescribeClass<LightClusterDecayer,Interfaced>
describeLightClusterDecayer("Herwig::LightClusterDecayer","Herwig.so");
IBPtr LightClusterDecayer::clone() const {
return new_ptr(*this);
}
IBPtr LightClusterDecayer::fullclone() const {
return new_ptr(*this);
}
void LightClusterDecayer::persistentOutput(PersistentOStream & os) const {
os << _hadronSpectrum;
}
void LightClusterDecayer::persistentInput(PersistentIStream & is, int) {
is >> _hadronSpectrum;
}
void LightClusterDecayer::Init() {
static ClassDocumentation<LightClusterDecayer> documentation
("There is the class responsible for the one-hadron decay of light clusters");
static Reference<LightClusterDecayer,HadronSpectrum>
interfaceHadronSpectrum("HadronSpectrum",
"A reference to the HadronSpectrum object",
&Herwig::LightClusterDecayer::_hadronSpectrum,
false, false, true, false);
}
bool LightClusterDecayer::decay(ClusterVector & clusters, tPVector & finalhadrons) {
// Loop over all clusters, and for those that were not heavy enough
// to undergo to fission, check if they are below the threshold
// for normal two-hadron decays. If this is the case, then the cluster
// should be decayed into a single hadron: this can happen only if
// it is possible to reshuffle momenta between the cluster and
// another one; in the rare occasions in which such exchange of momenta
// is not possible (because all of the clusters are too light) then
// the event is skipped.
// Notice that, differently from what happens in Fortran Herwig,
// light (that is below the threshold for the production of the lightest
// pair of hadrons with the proper flavours) fission products, produced
// by the fission of heavy clusters in class ClusterFissioner
// have been already "decayed" into single hadron (the lightest one
// with proper flavour) by the same latter class, without requiring
// any reshuffling. Therefore the light clusters that are treated in
// this LightClusterDecayer class are produced directly
// (originally) by the class ClusterFinder.
// To preserve all of the information, the cluster partner with which
// the light cluster (that decays into a single hadron) exchanges
// momentum in the reshuffling procedure is redefined and inserted
// in the vector vecNewRedefinedCluPtr. Only at the end, when all
// light clusters have been examined, the elements this vector will be
// copied in collecCluPtr (the reason is that it is not allowed to
// modify a STL container while iterating over it. At the same time,
// this ensures that a cluster can be redefined only once, which seems
// sensible although not strictly necessary).
// Notice that the cluster reshuffling partner is normally redefined
// and inserted in the vector vecNewRedefinedCluPtr, but not always:
// in the case it is also light, then it is also decayed immediately
// into a single hadron, without redefining it (the reason being that,
// otherwise, the would-be redefined cluster could have undefined
// components).
vector<tClusterPtr> redefinedClusters;
for (ClusterVector::const_iterator it = clusters.begin();
it != clusters.end(); ++it) {
// Skip the clusters that are not available or that are
// heavy, intermediate, clusters that have undergone to fission,
if ( ! (*it)->isAvailable() || ! (*it)->isReadyToDecay() ){
continue;
}
// We need to require (at least at the moment, maybe in the future we
// could change it) that the cluster has exactly two components,
// because otherwise we don't know how to deal with the kinematics.
// If this is not the case, then send a warning because it is not suppose
// to happen, and then do nothing with (ignore) such cluster.
if ( (*it)->numComponents() != 2 ) {
generator()->logWarning( Exception("LightClusterDecayer::decay "
"***Still cluster with not exactly"
" 2 components*** ",
Exception::warning) );
continue;
}
if ( DiquarkMatcher::Check((*it)->particle(0)->dataPtr()->id()) && DiquarkMatcher::Check((*it)->particle(1)->dataPtr()->id())) {
- // throw Exception() << "LightClusterDecayer::decay\n*** Diquark Cluster in LightClusterDecayer ***"
- // << Exception::runerror;
- continue;
+ // TODO We should never get Diquark Clusters in the LightClusterDecayer
+ throw Exception() << "LightClusterDecayer::decay\n"
+ "*** Diquark Cluster in LightClusterDecayer ***\n"
+ "Cluster = ( "<< (*it)->particle(0)->dataPtr()->id()<<", " << (*it)->particle(1)->dataPtr()->id()<<" )\nMC = " << (*it)->mass()/GeV << " GeV MLHP = "
+ << _hadronSpectrum->massLightestHadronPair((*it)->particle(0)->dataPtr(),(*it)->particle(1)->dataPtr())/GeV <<" GeV"
+ << Exception::runerror;
+ // continue;
}
// select the hadron for single hadron decay
tcPDPtr hadron = _hadronSpectrum->chooseSingleHadron((*it)->particle(0)->dataPtr(),
(*it)->particle(1)->dataPtr(),
(**it).mass());
// if not single decay continue
if(!hadron){
continue;
}
// We assume that the candidate reshuffling cluster partner,
// with whom the light cluster can exchange momenta,
// is chosen as the closest in space-time between the available
// clusters. Notice that an alternative, sensible approach
// could be to consider instead the "closeness" in the colour
// structure...
// Notice that nor a light cluster (which decays into a single hadron)
// neither its cluster reshuffling partner (which either has a
// redefined cluster or also decays into a single hadron) can be
// a reshuffling partner of another light cluster.
// This because we are requiring that the considered candidate cluster
// reshuffling partner has the status "isAvailable && isReadyToDecay" true;
// furthermore, the new redefined clusters are not added to the collection
// of cluster before the end of the entire reshuffling procedure, avoiding
// in this way that the redefined cluster of a cluster reshuffling partner
// is used again later. Needless to say, this is just an assumption,
// although reasonable, but nothing more than that!
// Build a multimap of available reshuffling cluster partners,
// with key given by the module of the invariant space-time distance
// w.r.t. the light cluster, so that this new collection is automatically
// ordered in increasing distance values.
// We use a multimap, rather than a map, just for precaution against not properly
// defined cluster positions which could produce all identical (null) distances.
multimap<Length,tClusterPtr> candidates;
for ( ClusterVector::iterator jt = clusters.begin();
jt != clusters.end(); ++jt ) {
if ( (*jt)->numComponents() != 2 )
continue;
// if (DiquarkMatcher::Check(*(*jt)->particle(0)->dataPtr())
// && DiquarkMatcher::Check(*(*jt)->particle(1)->dataPtr()))
// continue;
// if ( DiquarkMatcher::Check(*(*jt)->particle(0)->dataPtr())
// && DiquarkMatcher::Check(*(*jt)->particle(1)->dataPtr()))
// continue;
if ((*jt)->isAvailable() && (*jt)->isReadyToDecay() && jt != it) {
Length distance = abs (((*it)->vertex() - (*jt)->vertex()).m());
candidates.insert(pair<Length,tClusterPtr>(distance,*jt));
}
}
// Loop sequentially the multimap.
multimap<Length,tClusterPtr>::const_iterator mmapIt = candidates.begin();
bool found = false;
while (!found && mmapIt != candidates.end()) {
found = reshuffling(hadron, *it, (*mmapIt).second, redefinedClusters, finalhadrons);
if (!found) ++mmapIt;
}
if (!found) return partonicReshuffle(hadron,*it,finalhadrons);
} // end loop over collecCluPtr
// Add to collecCluPtr all of the redefined new clusters (indeed the
// pointers to them are added) contained in vecNewRedefinedCluPtr.
for (tClusterVector::const_iterator it = redefinedClusters.begin();
it != redefinedClusters.end(); ++it) {
clusters.push_back(*it);
}
return true;
}
bool LightClusterDecayer::reshuffling(const tcPDPtr pdata1,
tClusterPtr cluPtr1,
tClusterPtr cluPtr2,
tClusterVector & redefinedClusters,
tPVector & finalhadrons)
{
// don't reshuffle with beam clusters
if(cluPtr2->isBeamCluster()) return false;
// This method does the reshuffling of momenta between the cluster "1",
// that must decay into a single hadron (with id equal to idhad1), and
// the candidate cluster "2". It returns true if the reshuffling succeed,
// false otherwise.
PPtr ptrhad1 = pdata1->produceParticle();
if ( ! ptrhad1 ) {
generator()->logWarning( Exception("LightClusterDecayer::reshuffling"
"***Cannot create a particle with specified id***",
Exception::warning) );
return false;
}
Energy mhad1 = ptrhad1->mass();
// Let's call "3" and "4" the two constituents of the second cluster
tPPtr part3 = cluPtr2->particle(0);
tPPtr part4 = cluPtr2->particle(1);
// Check if the system of the two clusters can kinematically be replaced by
// an hadron of mass mhad1 (which is the lightest single hadron with the
// same flavour numbers as the first cluster) and the second cluster.
// If not, then try to replace the second cluster with the lightest hadron
// with the same flavour numbers; if it still fails, then give up!
Lorentz5Momentum pSystem = cluPtr1->momentum() + cluPtr2->momentum();
pSystem.rescaleMass(); // set the mass as the invariant of the quadri-vector
Energy mSystem = pSystem.mass();
Energy mclu2 = cluPtr2->mass();
bool singleHadron = false;
bool isDiquarkCluster = DiquarkMatcher::Check(part3->dataPtr()->id()) && DiquarkMatcher::Check(part4->dataPtr()->id());
Energy mLHP2 = _hadronSpectrum->massLightestHadronPair(part3->dataPtr(),part4->dataPtr());
// avoid calling massLightestHadron for Diquark clusters and only allow kinematic reshuffling
// for diquark clusters (no singleHadron)
Energy mLH2 = isDiquarkCluster ? mSystem:_hadronSpectrum->massLightestHadron(part3->dataPtr(),part4->dataPtr());
if(mSystem > mhad1 + mclu2 && mclu2 > mLHP2) { singleHadron = false; }
else if(mSystem > mhad1 + mLH2) { singleHadron = true; mclu2 = mLH2; }
else return false;
// Let's call from now on "Sys" the system of the two clusters, and
// had1 (of mass mhad1) the lightest hadron in which the first
// cluster decays, and clu2 (of mass mclu2) either the second
// cluster or the lightest hadron in which it decays (depending
// which one is kinematically allowed, see above).
// The idea behind the reshuffling is to replace the system of the
// two clusters by the system of the hadron had1 and (cluster or hadron) clu2,
// but leaving the overall system unchanged. Furthermore, the motion
// of had1 and clu2 in the Sys frame is assumed to be parallel to, respectively,
// those of the original cluster1 and cluster2 in the same Sys frame.
// Calculate the unit three-vector, in the frame "Sys" along which the
// two initial clusters move.
Lorentz5Momentum u( cluPtr1->momentum() );
u.boost( - pSystem.boostVector() ); // boost from LAB to Sys
// Calculate the momenta of had1 and clu2 in the Sys frame first,
// and then boost back in the LAB frame.
Lorentz5Momentum phad1, pclu2;
if (pSystem.m() < mhad1 + mclu2 ) {
throw Exception() << "Impossible Kinematics in LightClusterDecayer::reshuffling()"
<< Exception::eventerror;
}
Kinematics::twoBodyDecay(pSystem, mhad1, mclu2, u.vect().unit(), phad1, pclu2);
ptrhad1->set5Momentum( phad1 ); // set momentum of first hadron.
ptrhad1->setVertex(cluPtr1->vertex()); // set hadron vertex position to the
// parent cluster position.
cluPtr1->addChild(ptrhad1);
finalhadrons.push_back(ptrhad1);
cluPtr1->flagAsReshuffled();
cluPtr2->flagAsReshuffled();
if(singleHadron) {
// In the case that also the cluster reshuffling partner is light
// it is decayed into a single hadron, *without* creating the
// redefined cluster (this choice is justified in order to avoid
// clusters that could have undefined components).
PPtr ptrhad2 = _hadronSpectrum->lightestHadron(part3->dataPtr(),part4->dataPtr())
->produceParticle();
ptrhad2->set5Momentum( pclu2 );
ptrhad2->setVertex( cluPtr2->vertex() ); // set hadron vertex position to the
// parent cluster position.
cluPtr2->addChild(ptrhad2);
finalhadrons.push_back(ptrhad2);
} else {
// Create the new cluster which is the redefinitions of the cluster
// partner (cluster "2") used in the reshuffling procedure of the
// light cluster (cluster "1").
// The rationale of this is to preserve completely all of the information.
ClusterPtr cluPtr2new = ClusterPtr();
if(part3 && part4) cluPtr2new = new_ptr(Cluster(part3,part4));
cluPtr2new->set5Momentum( pclu2 );
cluPtr2new->setVertex( cluPtr2->vertex() );
cluPtr2->addChild( cluPtr2new );
redefinedClusters.push_back( cluPtr2new );
// Set consistently the momenta of the two components of the second cluster
// after the reshuffling. To do that we first calculate the momenta of the
// constituents in the initial cluster rest frame; then we boost them back
// in the lab but using this time the new cluster rest frame. Finally we store
// these information in the new cluster. Notice that we do *not* set
// consistently also the momenta of the (eventual) particles pointed by the
// two components: that's because we do not need to do so, being the momentum
// an explicit private member of the class Component (which is set equal
// to the momentum of the eventual particle pointed only in the constructor,
// but then later should not necessary be the same), and furthermore it allows
// us not to loose any information, in the sense that we can always, later on,
// to find the original momenta of the two components before the reshuffling.
Lorentz5Momentum p3 = part3->momentum(); //p3new->momentum();
p3.boost( - (cluPtr2->momentum()).boostVector() ); // from LAB to clu2 (old) frame
p3.boost( pclu2.boostVector() ); // from clu2 (new) to LAB frame
Lorentz5Momentum p4 = part4->momentum(); //p4new->momentum();
p4.boost( - (cluPtr2->momentum()).boostVector() ); // from LAB to clu2 (old) frame
p4.boost( pclu2.boostVector() ); // from clu2 (new) to LAB frame
cluPtr2new->particle(0)->set5Momentum(p3);
cluPtr2new->particle(1)->set5Momentum(p4);
} // end of if (singleHadron)
return true;
}
bool LightClusterDecayer::partonicReshuffle(const tcPDPtr had,
const PPtr cluster,
tPVector & finalhadrons) {
tPPtr meson(cluster);
if(!meson->parents().empty()) meson=meson->parents()[0];
if(!meson->parents().empty()) meson=meson->parents()[0];
// check b/c hadron decay
int ptype(abs(meson->id())%10000);
bool heavy = (ptype/1000 == 5 || ptype/1000 ==4 );
heavy |= (ptype/100 == 5 || ptype/100 ==4 );
heavy |= (ptype/10 == 5 || ptype/10 ==4 );
if(!heavy) return false;
// find the leptons
tPVector leptons;
for(unsigned int ix=0;ix<meson->children().size();++ix) {
if(!(meson->children()[ix]->dataPtr()->coloured())) {
leptons.push_back(meson->children()[ix]);
}
}
if(leptons.size()==1) {
tPPtr w=leptons[0];
leptons.pop_back();
for(unsigned int ix=0;ix<w->children().size();++ix) {
if(!w->children()[ix]->dataPtr()->coloured()) {
leptons.push_back(w->children()[ix]);
}
}
}
if(leptons.size()!=2) return false;
// get momentum of leptonic system and the its minimum possible mass
Energy mmin(ZERO);
Lorentz5Momentum pw;
for(unsigned int ix=0;ix<leptons.size();++ix) {
pw+=leptons[ix]->momentum();
mmin+=leptons[ix]->mass();
}
pw.rescaleMass();
// check we can do the reshuffling
PPtr ptrhad = had->produceParticle();
// total momentum fo the system
Lorentz5Momentum pSystem = pw + cluster->momentum();
pSystem.rescaleMass();
// normal case get additional energy by rescaling momentum in rest frame of
// system
if(pSystem.mass()>ptrhad->mass()+pw.mass()&&pw.mass()>mmin) {
// Calculate the unit three-vector, in the frame "Sys" along which the
// two initial clusters move.
Lorentz5Momentum u(cluster->momentum());
u.boost( - pSystem.boostVector() );
// Calculate the momenta of had1 and clu2 in the Sys frame first,
// and then boost back in the LAB frame.
Lorentz5Momentum phad1, pclu2;
Kinematics::twoBodyDecay(pSystem, ptrhad->mass(), pw.mass(),
u.vect().unit(), phad1, pclu2);
// set momentum of first hadron.
ptrhad->set5Momentum( phad1 );
// set hadron vertex position to the parent cluster position.
ptrhad->setLabVertex(cluster->vertex());
// add hadron
cluster->addChild(ptrhad);
finalhadrons.push_back(ptrhad);
// reshuffle the leptons
// boost the leptons to the rest frame of the system
Boost boost1(-pw.boostVector());
Boost boost2( pclu2.boostVector());
for(unsigned int ix=0;ix<leptons.size();++ix) {
leptons[ix]->deepBoost(boost1);
leptons[ix]->deepBoost(boost2);
}
return true;
}
else {
return false;
}
}
diff --git a/src/defaults/Hadronization.in b/src/defaults/Hadronization.in
--- a/src/defaults/Hadronization.in
+++ b/src/defaults/Hadronization.in
@@ -1,142 +1,162 @@
# -*- 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:KineticThresholdShift 0.08844
set /Herwig/Hadronization/ClusterFissioner:ProbablityPowerFactor 6.486
set /Herwig/Hadronization/ClusterFissioner:ProbablityShift -0.87875
# Clustering parameters for light quarks
newdef ClusterFissioner:ClMaxLight 3.528693*GeV
newdef ClusterFissioner:ClPowLight 1.849375
newdef ClusterFissioner:PSplitLight 0.914156
insert ClusterFissioner:FissionPwt 1 1.0
insert ClusterFissioner:FissionPwt 2 1.0
insert ClusterFissioner:FissionPwt 3 0.374094
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
insert ClusterFissioner:PSplitHeavy 5 0.625
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*GeV
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
+# Cluster Paramters for light Diquark Cluster
+# currently set according to Light quark defaults
+insert ClusterFissioner:ClMaxDiquark 1103 3.528693*GeV
+insert ClusterFissioner:ClPowDiquark 1103 1.849375
+insert ClusterFissioner:ClMaxDiquark 2101 3.528693*GeV
+insert ClusterFissioner:ClPowDiquark 2101 1.849375
+insert ClusterFissioner:ClMaxDiquark 2103 3.528693*GeV
+insert ClusterFissioner:ClPowDiquark 2103 1.849375
+insert ClusterFissioner:ClMaxDiquark 2203 3.528693*GeV
+insert ClusterFissioner:ClPowDiquark 2203 1.849375
+insert ClusterFissioner:ClMaxDiquark 3101 3.528693*GeV
+insert ClusterFissioner:ClPowDiquark 3101 1.849375
+insert ClusterFissioner:ClMaxDiquark 3103 3.528693*GeV
+insert ClusterFissioner:ClPowDiquark 3103 1.849375
+insert ClusterFissioner:ClMaxDiquark 3201 3.528693*GeV
+insert ClusterFissioner:ClPowDiquark 3201 1.849375
+insert ClusterFissioner:ClMaxDiquark 3203 3.528693*GeV
+insert ClusterFissioner:ClPowDiquark 3203 1.849375
+insert ClusterFissioner:ClMaxDiquark 3303 3.528693*GeV
+insert ClusterFissioner:ClPowDiquark 3303 1.849375
# 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.374094
newdef SMHadronSpectrum:PwtCquark 0.0
newdef SMHadronSpectrum:PwtBquark 0.0
newdef SMHadronSpectrum:PwtDIquark 0.33107
newdef SMHadronSpectrum:SngWt 0.89050
newdef SMHadronSpectrum:DecWt 0.41628
newdef SMHadronSpectrum:Mode 1
newdef SMHadronSpectrum:BelowThreshold All
create Herwig::SpinHadronizer SpinHadronizer