Page MenuHomeHEPForge

No OneTemporary

This file is larger than 256 KB, so syntax highlighting was skipped.
diff --git a/Hadronization/ClusterFissioner.cc b/Hadronization/ClusterFissioner.cc
--- a/Hadronization/ClusterFissioner.cc
+++ b/Hadronization/ClusterFissioner.cc
@@ -1,3125 +1,3149 @@
// -*- 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 "Herwig/Utilities/AlphaS.h"
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <vector>
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(0),
_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(_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(_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
;
}
namespace{
void printV(Lorentz5Momentum p) {
std::cout << "("<<p.e()/GeV<<"|"<<p.vect().x()/GeV<<","<<p.vect().y()/GeV<<","<<p.vect().z()/GeV<<") Mass = "<<p.mass()/GeV<<" m = "<<p.m()/GeV<<"\n";
}
}
void ClusterFissioner::doinit() {
- Interfaced::doinit();
+ // Interfaced::doinit();
+ // Energy m1=1.2*GeV;
+ // Energy m2=2.3*GeV;
+ // std::cout <<std::setprecision(18)<< sqr(m1+m2)/GeV2 << std::endl;
+ // std::cout <<std::setprecision(18)<< ((m1+m2)*(m1+m2))/GeV2 << std::endl;
+ // std::cout <<std::setprecision(18)<< ((sqr(m1+m2))-((m1+m2)*(m1+m2)))/GeV2 << std::endl;
/* TESTs for Two particle boost
* TODO
* */
/*
Energy M,M1,M2;
Energy m,m1,m2;
M=20.1*GeV;
M1=10.1*GeV;
M2=5.1*GeV;
m2=0.325*GeV;
m1=0.625*GeV;
m=1.1*GeV;
// double cosTheta1=0.3;
// double phi1=0.7*3.14;
// double cosTheta2=-0.4;
// double phi2=-0.2*3.14;
// Axis LABz(0,0,1);
// Energy Plab=80.1*GeV;
// Lorentz5Momentum P(M,Plab*LABz);
// Energy PCOM=P.boost(-P.boostVector()).vect().mag();
double cosTheta=0.0;
double phi=0.7*3.14;
Energy P1=1.2*GeV;
Energy P2=10.2*GeV;
Energy Q=1.2*GeV;
Lorentz5Momentum p1(m1,P1*Axis(cosTheta*cos(phi),cosTheta*sin(phi),sin(acos(cosTheta))));
Lorentz5Momentum p2(m1*2,Momentum3(-p1.vect().x(),-p1.vect().y(),P2));
Lorentz5Momentum qOld(m,Q*sampleDirectionIsotropic());
Lorentz5Momentum qNew=qOld;
Lorentz5Momentum pC=p1+p2;
printV(p1);
printV(p2);
printV(p1+p2);
printV(pC);
printV(qOld);
printV(qNew);
// p1.boost(-pC.boostVector());
// p2.boost(-pC.boostVector());
// q.boost(-pC.boostVector());
// printV(q);
std::vector<Lorentz5Momentum *> v;
v.push_back(&qNew);
Kinematics::BoostIntoTwoParticleFrame(pC.m(),p1,p2,v);
// Kinematics::twoBodDecay(pC,);
printV(qOld);
printV(qNew);
// std::cout << p1 << std::endl;
// std::cout << p2 << std::endl;
// std::cout << p1+p2 << std::endl;
// std::cout << Q << std::endl;
*/
// 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 SwitchOption interfaceFissionNewDiquarkSuppression
(interfaceFission,
"NewDiquarkSuppression",
"Alternative cluster fission which does not depend on the hadron spectrum class"
" and includes a suppression of AlphaS^2(Mc) for Diquark Production during "
"Cluster Fission",
-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);
static SwitchOption interfaceFissionApproachPreservePert
(interfaceFissionApproach,
"PreservePert",
"New cluster fission which allows to choose MassSampler"
", PhaseSpaceSampler and MatrixElement, but uses Default"
" Approach for BeamClusters",
3);
static SwitchOption interfaceFissionApproachPreserveNonPert
(interfaceFissionApproach,
"PreserveNonPert",
"New cluster fission which allows to choose MassSampler"
", PhaseSpaceSampler and MatrixElement, but uses Default"
" Approach for BeamClusters",
4);
static SwitchOption interfaceFissionApproachPreserveFirstPert
(interfaceFissionApproach,
"PreserveFirstPert",
"New cluster fission which allows to choose MassSampler"
", PhaseSpaceSampler and MatrixElement, but uses Default"
" Approach for BeamClusters",
5);
// 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);
+ static SwitchOption interfaceMatrixElementTest
+ (interfaceMatrixElement,
+ "Test",
+ "Choose Soft p1,p2->q1,q2,g*->q1,q2,q,qbar matrix element "
+ "NOTE: Here the matrix element depends on pi.q(bar)",
+ 4);
// 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,int> interfacePhaseSpaceWeights
("PhaseSpaceWeights",
"Include phase space weights.",
&ClusterFissioner::_phaseSpaceWeights, 0, false, false);
static SwitchOption interfacePhaseSpaceWeightsNo
(interfacePhaseSpaceWeights,
"No",
"Do not include the effect of cluster phase space",
0);
static SwitchOption interfacePhaseSpaceWeightsYes
(interfacePhaseSpaceWeights,
"Yes",
"Do include the effect of cluster fission phase space "
"related to constituent masses."
"Note: Need static Threshold choice",
1);
static SwitchOption interfacePhaseSpaceWeightsUseHadronMasses
(interfacePhaseSpaceWeights,
"UseHadronMasses",
"Do include the effect of cluster fission phase space "
"related to hadron masses."
"Note: Need static Threshold choice",
2);
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;
case 3:
{
bool isPert=false;
for (int i = 0; i < cluster->numComponents(); i++)
{
if (cluster->isPerturbative(i))
isPert=true;
}
if (isPert)
return cutTwoDefault(cluster, finalhadrons, softUEisOn);
else
return cutTwoNew(cluster, finalhadrons, softUEisOn);
break;
}
case 4:
{
bool isPert=false;
for (int i = 0; i < cluster->numComponents(); i++)
{
if (cluster->isPerturbative(i))
isPert=true;
}
if (!isPert)
return cutTwoDefault(cluster, finalhadrons, softUEisOn);
else
return cutTwoNew(cluster, finalhadrons, softUEisOn);
break;
}
case 5:
{
bool isPert=true;
for (int i = 0; i < cluster->numComponents(); i++)
{
if (!cluster->isPerturbative(i)) {
isPert=false;
break;
}
}
if (isPert)
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;
if (_phaseSpaceWeights==2 &&
( Mc1 < spectrum()->massLightestHadronPair(ptrQ1->dataPtr(),newPtr1->dataPtr())
|| Mc2 < spectrum()->massLightestHadronPair(ptrQ2->dataPtr(),newPtr2->dataPtr()) ))
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, ptrQ1, ptrQ2, newPtr1) ) {
// reduce counter as it regards only the mass sampling
counter--;
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;
}
namespace {
int areNotSame(const Lorentz5Momentum & p1,const Lorentz5Momentum & p2){
double tol=1e-5;
if (abs(p1.vect().x()-p2.vect().x())>tol*GeV){
std::cout << "disagreeing x " <<std::setprecision(-log10(tol)+2) << abs(p1.vect().x()-p2.vect().x())/GeV<< std::endl;
return 1;
}
if (abs(p1.vect().y()-p2.vect().y())>tol*GeV){
std::cout << "disagreeing y " <<std::setprecision(-log10(tol)+2)<< abs(p1.vect().y()-p2.vect().y())/GeV<< std::endl;
return 2;
}
if (abs(p1.vect().z()-p2.vect().z())>tol*GeV){
std::cout << "disagreeing z " <<std::setprecision(-log10(tol)+2)<< abs(p1.vect().z()-p2.vect().z())/GeV<< std::endl;
return 3;
}
if (abs(p1.e()-p2.e())>tol*GeV){
std::cout << "disagreeing e " <<std::setprecision(-log10(tol)+2)<< abs(p1.e()-p2.e())/GeV<< std::endl;
return 4;
}
if (abs(p1.m()-p2.m())>tol*GeV){
std::cout << "disagreeing m " <<std::setprecision(-log10(tol)+2)<< abs(p1.m()-p2.m())/GeV<< std::endl;
return 4;
}
if (abs(p1.m()-p1.mass())>tol*GeV){
std::cout << "disagreeing p1 m - mass " <<std::setprecision(-log10(tol)+2)<< abs(p1.m()-p1.mass())/GeV<< std::endl;
return 4;
}
if (abs(p2.m()-p2.mass())>tol*GeV){
std::cout << "disagreeing p2 m - mass " <<std::setprecision(-log10(tol)+2)<< abs(p2.m()-p2.mass())/GeV<< std::endl;
return 4;
}
return 0;
}
}
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)
+ const Lorentz5Momentum p0Q1 = ptrQ1->momentum(); // known (mom Q1 before fission)
+ const 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();
+ // 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);
{ //sanity
if (areNotSame(pClu1+pClu2,pClu)) {
std::cout << "PCLU" << std::endl;
}
if (areNotSame(pClu1,pQ1+pQone)) {
std::cout << "PCLU1" << std::endl;
}
if (areNotSame(pClu2,pQ2+pQtwo)) {
std::cout << "PCLU2" << std::endl;
}
if (areNotSame(pClu,pQ1+pQone+pQ2+pQtwo)) {
std::cout << "PTOT" << std::endl;
}
}
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 = _matrixElement==0 ? 1.0: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);
}
assert(abs(Mc1-pClu1.m())<1e-2*GeV);
assert(abs(Mc2-pClu2.m())<1e-2*GeV);
assert(abs(Mc1-pClu1.mass())<1e-2*GeV);
assert(abs(Mc2-pClu2.mass())<1e-2*GeV);
{
- Lorentz5Momentum pp1(ptrQ1->momentum());
+ Lorentz5Momentum pp1(p0Q1);
Lorentz5Momentum pclu1(pClu1);
pclu1.boost(-pClu.boostVector());
pp1.boost(-pClu.boostVector());
- // std::cout << "angle between initial Clu and final "<< pclu1.vect().cosTheta(pp1.vect())<< std::endl;
+ // std::ofstream out("testCF.dat", std::ios::app);
+ // out << pclu1.vect().cosTheta(pp1.vect())<<"\n";
+ // out.close();
}
// ==> 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) ) {
// reduce counter as it regards only the mass sampling
counter--;
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, tcPPtr pQ1, tcPPtr pQ2, tcPPtr pQ) const {
switch (_phaseSpaceWeights)
{
case 1:
return phaseSpaceVetoConstituentMasses(Mc, Mc1, Mc2, m, m1, m2);
case 2:
return phaseSpaceVetoHadronPairs(Mc, Mc1, Mc2, pQ1, pQ2, pQ);
default:
assert(false);
}
}
/**
* Veto for the phase space weight
* returns true if proposed Masses are rejected
* else returns false
*/
bool ClusterFissioner::phaseSpaceVetoConstituentMasses(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));
}
bool ClusterFissioner::phaseSpaceVetoHadronPairs(const Energy Mc, const Energy Mc1, const Energy Mc2, tcPPtr pQ1, tcPPtr pQ2, tcPPtr pQ, const Energy , const Energy , const Energy ) const {
auto LHP1 = spectrum()->lightestHadronPair(pQ1->dataPtr(),pQ->dataPtr());
auto LHP2 = spectrum()->lightestHadronPair(pQ2->dataPtr(),pQ->dataPtr());
// double weigthHadrons = Kinematics::pstarTwoBodyDecay(Mc1,m1,m)*Kinematics::pstarTwoBodyDecay(Mc2,m2,m)/(Kinematics::pstarTwoBodyDecay(Mc1,LHP1.first->mass(),LHP1.second->mass())*Kinematics::pstarTwoBodyDecay(Mc2,LHP2.first->mass(),LHP2.second->mass()));
// double tot= weightFlatPhaseSpace(Mc, Mc1, Mc2, m, m1, m2)/weigthHadrons;
// if (std::isinf(tot) || std::isnan(tot) || tot<0.0 || tot>1.0)
// std::cout << "tot = " <<
double lam1 = sqrt(Kinematics::kaellen(Mc1/GeV, LHP1.first->mass()/GeV, LHP1.second->mass()/GeV));
double lam2 = sqrt(Kinematics::kaellen(Mc2/GeV, LHP2.first->mass()/GeV, LHP2.second->mass()/GeV));
double lam3 = sqrt(Kinematics::kaellen(Mc/GeV, Mc1/GeV, Mc2/GeV));
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(Mc1*Mc2/GeV2,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(Mc/GeV, 4.*_dim-12.);
ratio = PSweight/overEstimate;
// if (!(ratio>=0)) std::cout << "M "<<Mc/GeV<<" M1 "<<Mc1/GeV<<" M2 "<<Mc2/GeV<<" m1 "<<m1_temp<<" m2 "<<m2_temp<<" m "<<m_temp<<"\n\n";
// if (ratio > 0.9) std::cout << "ratio = " << ratio <<"\n";
if (std::isinf(ratio) || std::isnan(ratio) || ratio<0.0 || ratio>1.0)
std::cout << "ratio = " << ratio<<std::endl;
assert (ratio >= 0);
assert (ratio <= 1);
return (UseRandom::rnd()>ratio);
}
/**
* Calculate the masses and possibly kinematics of the cluster
* fission at hand; if claculateKineamtics is perfomring non-trivial
* steps kinematics calulcated here will be overriden. Currently resorts to the default
*/
bool ClusterFissioner::drawNewMasses(const Energy Mc, const bool soft1, const bool soft2,
Lorentz5Momentum& pClu1, Lorentz5Momentum& pClu2,
tcPPtr ptrQ1, const Lorentz5Momentum& pQ1,
tcPPtr, const Lorentz5Momentum& pQone,
tcPPtr, const Lorentz5Momentum& pQtwo,
tcPPtr ptrQ2, const Lorentz5Momentum& pQ2) const {
// TODO 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;
double weight;
// double factorPS;
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++;
minMass=spectrum()->massLightestHadronPair(pD1,pD2);
if (_fissionCluster==0) choice.insert(_hadronSpectrum->pwtQuark(id),id);
else if (abs(_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);
Energy mH1=spectrum()->massLightestHadron(pD2,cand);
Energy mH2=spectrum()->massLightestHadron(cand,pD1);
// factorPS = Kinematics::pstarTwoBodyDecay(Mc,mH1,mH2)/(Mc/2.0);
// factorPS = 1.0;
minMass = mH1 + mH2;
}
else {
minMass = spectrum()->massLightestBaryonPair(pD1,pD2);
// factorPS = 1.0;
}
if (Mc < minMass) continue;
// countDiQ++;
if (_fissionCluster==0) weight = _hadronSpectrum->pwtQuark(id);
else if (abs(_fissionCluster)==1) weight = _fissionPwt.find(id)->second;
else assert(false);
if (_fissionCluster==-1)
weight*=sqr(Herwig::Math::alphaS(Mc, 0.25*GeV,3, 2));
// weight/=factorPS;
choice.insert(weight,id);
}
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) weight = _hadronSpectrum->pwtQuark(id);
else if (abs(_fissionCluster)==1) weight = _fissionPwt.find(id)->second;
else assert(false);
if (_fissionCluster==-1)
weight*=sqr(Herwig::Math::alphaS(Mc, 0.25*GeV,3, 2));
choice.insert(weight,id);
}
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) weight = _hadronSpectrum->pwtQuark(id);
else if (abs(_fissionCluster)==1) weight = _fissionPwt.find(id)->second;
else assert(false);
if (_fissionCluster==-1)
weight*=sqr(Herwig::Math::alphaS(Mc, 0.25*GeV,3, 2));
choice.insert(weight,id);
}
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(abs(_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(abs(_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::sampleDirectionAlignedSmeared(const Lorentz5Momentum & pQ, const Lorentz5Momentum & pClu) const {
Lorentz5Momentum pQinCOM(pQ);
pQinCOM.setMass(pQ.m());
pQinCOM.boost( -pClu.boostVector() ); // boost from LAB to C
if (pQinCOM.vect().mag2()<=ZERO){
return sampleDirectionAlignedSmeared(pClu-pQ,pClu);
}
const Axis dir = pQinCOM.vect().unit();
double cluSmear = 0.5;
double cosTheta;
do {
cosTheta = 1.0 + cluSmear*log( UseRandom::rnd() );
}
while (fabs(cosTheta)>1.0 || std::isnan(cosTheta) || std::isinf(cosTheta));
if (!(dir.mag2()>ZERO)){
std::cout << "\nDRI = 0"<< dir.mag2() << std::endl;
}
Axis dirSmeared = dir;
double phi;
if (fabs(cosTheta-1.0)>1e-10) {
// rotate to sampled angles
dirSmeared.rotate(acos(cosTheta),dir.orthogonal());
phi = UseRandom::rnd(-M_PI,M_PI);
dirSmeared.rotate(phi,dir);
}
if (!(dirSmeared.mag2()>ZERO)){
std::cout << "\nDRI SMR = 0" <<cosTheta<< "\t" <<acos(cosTheta)<< "\t"<< phi<< "\t"<< dirSmeared.mag2() << std::endl;
std::cout << dir.orthogonal() << std::endl;
}
return dirSmeared;
}
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)*sinTheta, sin(Phi)*sinTheta, cosTheta);
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;
}
+ case 4:
+ {
+ SQME = sqr(GeV2)/((p1*(q1+q))*(p2*(q2+qbar)));
+ 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& m1,
+ const Energy& m2,
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;
}
+ case 4:
+ {
+ // Fit factor for guess of best overestimate
+ SQME_OverEstimate = _safetyFactorMatrixElement*sqr(GeV2)/(m1*m2*(m1+mq)*(m2+mq));
+ 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 {
Lorentz5Momentum pClutmp(pClu);
pClutmp.boost(-pClu.boostVector());
Kinematics::twoBodyDecay(pClutmp, 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;
Lorentz5Momentum pClu1tmp(pClu1);
// pClu1tmp.boost(-pClu.boostVector());
switch (_phaseSpaceSampler)
{
case 0:
// Default aligned
DirClu1 = _covariantBoost ? pClu1tmp.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);
Kinematics::twoBodyDecay(pClu1tmp, 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;
Lorentz5Momentum pClu2tmp(pClu2);
// pClu2tmp.boost(-pClu.boostVector());
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(pClu2tmp, 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);
}
else {
// Kinematics::twoBodyDecay(pClu, pClu1.mass(), pClu2.mass(),pClu1.vect().unit(), pClu1, pClu2);
// Kinematics::twoBodyDecay(pClu, pQ1.mass(), pQbar.mass(), pQ1.vect().unit(), pQ1, pQbar);
// Kinematics::twoBodyDecay(pClu, pQ.mass(), pQ2bar.mass(), pQ.vect().unit(), pQ, pQ2bar);
pClu1.boost(pClu.boostVector());
pClu2.boost(pClu.boostVector());
pQ1.boost(pClu.boostVector());
pQ2bar.boost(pClu.boostVector());
pQ.boost(pClu.boostVector());
pQbar.boost(pClu.boostVector());
}
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()};
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];
}
for ( const long& id : spectrum()->heavyHadronizingQuarks() ) {
if ( spectrum()->hasHeavy(id,cptr[0],cptr[1],cptr[1])) {
clpow = _clPowHeavy[id];
clmax = _clMaxHeavy[id];
}
}
// required test for SUSY clusters, since aboveCutoff alone
// cannot guarantee (Mc > m1 + m2 + 2*m) in cut()
static const Energy minmass
= getParticleData(ParticleID::d)->constituentMass();
bool aboveCutoff = false, canSplitMinimally = false;
// static kinematic threshold
if(_kinematicThresholdChoice == 0) {
aboveCutoff = (
pow(clu->mass()*UnitRemoval::InvE , clpow)
>
pow(clmax*UnitRemoval::InvE, clpow)
+ pow(clu->sumConstituentMasses()*UnitRemoval::InvE, clpow)
);
canSplitMinimally = clu->mass() > clu->sumConstituentMasses() + 2.0 * minmass;
}
// dynamic kinematic threshold
else if(_kinematicThresholdChoice == 1) {
//some smooth probablity function to create a dynamic thershold
double scale = pow(clu->mass()/GeV , clpow);
double threshold = pow(clmax/GeV, clpow)
+ pow(clu->sumConstituentMasses()/GeV, clpow);
aboveCutoff = ProbablityFunction(scale,threshold);
scale = clu->mass()/GeV;
threshold = clu->sumConstituentMasses()/GeV + 2.0 * minmass/GeV;
canSplitMinimally = ProbablityFunction(scale,threshold);
}
return aboveCutoff && canSplitMinimally;
}
diff --git a/Hadronization/ClusterHadronizationHandler.cc b/Hadronization/ClusterHadronizationHandler.cc
--- a/Hadronization/ClusterHadronizationHandler.cc
+++ b/Hadronization/ClusterHadronizationHandler.cc
@@ -1,592 +1,592 @@
// -*- C++ -*-
//
// ClusterHadronizationHandler.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the ClusterHadronizationHandler class.
//
#include "ClusterHadronizationHandler.h"
#include <ThePEG/Interface/ClassDocumentation.h>
#include <ThePEG/Persistency/PersistentOStream.h>
#include <ThePEG/Persistency/PersistentIStream.h>
#include <ThePEG/Interface/Switch.h>
#include <ThePEG/Interface/Parameter.h>
#include <ThePEG/Interface/Reference.h>
#include <ThePEG/Handlers/EventHandler.h>
#include <ThePEG/Handlers/Hint.h>
#include <ThePEG/PDT/ParticleData.h>
#include <ThePEG/EventRecord/Particle.h>
#include <ThePEG/EventRecord/Step.h>
#include <ThePEG/PDT/PDT.h>
#include <ThePEG/PDT/EnumParticles.h>
#include <ThePEG/Utilities/Throw.h>
#include "Herwig/Utilities/EnumParticles.h"
#include "CluHadConfig.h"
#include "Cluster.h"
#include <ThePEG/Utilities/DescribeClass.h>
#include <ThePEG/Interface/Command.h>
using namespace Herwig;
ClusterHadronizationHandler * ClusterHadronizationHandler::currentHandler_ = 0;
DescribeClass<ClusterHadronizationHandler,HadronizationHandler>
describeClusterHadronizationHandler("Herwig::ClusterHadronizationHandler","Herwig.so");
IBPtr ClusterHadronizationHandler::clone() const {
return new_ptr(*this);
}
IBPtr ClusterHadronizationHandler::fullclone() const {
return new_ptr(*this);
}
void ClusterHadronizationHandler::persistentOutput(PersistentOStream & os)
const {
os << _partonSplitters << _clusterFinders << _colourReconnectors
<< _clusterFissioners << _lightClusterDecayers << _clusterDecayers
<< _gluonMassGenerators
<< _partonSplitter << _clusterFinder << _colourReconnector
<< _clusterFissioner << _lightClusterDecayer << _clusterDecayer
<< reshuffle_ << _gluonMassGenerator
<< ounit(_minVirtuality2,GeV2) << ounit(_maxDisplacement,mm)
<< _underlyingEventHandler << _reduceToTwoComponents;
}
void ClusterHadronizationHandler::persistentInput(PersistentIStream & is, int) {
is >> _partonSplitters >> _clusterFinders >> _colourReconnectors
>> _clusterFissioners >> _lightClusterDecayers >> _clusterDecayers
>> _gluonMassGenerators
>> _partonSplitter >> _clusterFinder >> _colourReconnector
>> _clusterFissioner >> _lightClusterDecayer >> _clusterDecayer
>> reshuffle_ >> _gluonMassGenerator
>> iunit(_minVirtuality2,GeV2) >> iunit(_maxDisplacement,mm)
>> _underlyingEventHandler >> _reduceToTwoComponents;
}
void ClusterHadronizationHandler::Init() {
static ClassDocumentation<ClusterHadronizationHandler> documentation
("This is the main handler class for the Cluster Hadronization",
"The hadronization was performed using the cluster model of \\cite{Webber:1983if}.",
"%\\cite{Webber:1983if}\n"
"\\bibitem{Webber:1983if}\n"
" B.~R.~Webber,\n"
" ``A QCD Model For Jet Fragmentation Including Soft Gluon Interference,''\n"
" Nucl.\\ Phys.\\ B {\\bf 238}, 492 (1984).\n"
" %%CITATION = NUPHA,B238,492;%%\n"
// main manual
);
static Reference<ClusterHadronizationHandler,PartonSplitter>
interfacePartonSplitter("PartonSplitter",
"A reference to the PartonSplitter object",
&Herwig::ClusterHadronizationHandler::_partonSplitter,
false, false, true, false);
static Reference<ClusterHadronizationHandler,GluonMassGenerator> interfaceGluonMassGenerator
("GluonMassGenerator",
"Set a reference to a gluon mass generator.",
&ClusterHadronizationHandler::_gluonMassGenerator, false, false, true, true, false);
static Switch<ClusterHadronizationHandler,int> interfaceReshuffle
("Reshuffle",
"Perform reshuffling if constituent masses have not yet been included by the shower",
&ClusterHadronizationHandler::reshuffle_, 0, false, false);
static SwitchOption interfaceReshuffleColourConnected
(interfaceReshuffle,
"ColourConnected",
"Separate reshuffling for colour connected partons.",
2);
static SwitchOption interfaceReshuffleGlobal
(interfaceReshuffle,
"Global",
"Reshuffle globally.",
1);
static SwitchOption interfaceReshuffleNo
(interfaceReshuffle,
"No",
"Do not reshuffle.",
0);
static Reference<ClusterHadronizationHandler,ClusterFinder>
interfaceClusterFinder("ClusterFinder",
"A reference to the ClusterFinder object",
&Herwig::ClusterHadronizationHandler::_clusterFinder,
false, false, true, false);
static Reference<ClusterHadronizationHandler,ColourReconnector>
interfaceColourReconnector("ColourReconnector",
"A reference to the ColourReconnector object",
&Herwig::ClusterHadronizationHandler::_colourReconnector,
false, false, true, false);
static Reference<ClusterHadronizationHandler,ClusterFissioner>
interfaceClusterFissioner("ClusterFissioner",
"A reference to the ClusterFissioner object",
&Herwig::ClusterHadronizationHandler::_clusterFissioner,
false, false, true, false);
static Reference<ClusterHadronizationHandler,LightClusterDecayer>
interfaceLightClusterDecayer("LightClusterDecayer",
"A reference to the LightClusterDecayer object",
&Herwig::ClusterHadronizationHandler::_lightClusterDecayer,
false, false, true, false);
static Reference<ClusterHadronizationHandler,ClusterDecayer>
interfaceClusterDecayer("ClusterDecayer",
"A reference to the ClusterDecayer object",
&Herwig::ClusterHadronizationHandler::_clusterDecayer,
false, false, true, false);
// functions to move the handlers so they are set
static Parameter<ClusterHadronizationHandler,Energy2> interfaceMinVirtuality2
("MinVirtuality2",
"Minimum virtuality^2 of partons to use in calculating distances (unit [GeV2]).",
&ClusterHadronizationHandler::_minVirtuality2, GeV2, 0.1*GeV2, ZERO, 10.0*GeV2,false,false,false);
static Parameter<ClusterHadronizationHandler,Length> interfaceMaxDisplacement
("MaxDisplacement",
"Maximum displacement that is allowed for a particle (unit [millimeter]).",
&ClusterHadronizationHandler::_maxDisplacement, mm, 1.0e-10*mm,
0.0*mm, 1.0e-9*mm,false,false,false);
static Reference<ClusterHadronizationHandler,StepHandler> interfaceUnderlyingEventHandler
("UnderlyingEventHandler",
"Pointer to the handler for the Underlying Event. "
"Set to NULL to disable.",
&ClusterHadronizationHandler::_underlyingEventHandler, false, false, true, true, false);
static Switch<ClusterHadronizationHandler,bool> interfaceReduceToTwoComponents
("ReduceToTwoComponents",
"Whether or not to reduce three component baryon-number violating clusters to two components before cluster splitting or leave"
" this till after the cluster splitting",
&ClusterHadronizationHandler::_reduceToTwoComponents, true, false, false);
static SwitchOption interfaceReduceToTwoComponentsYes
(interfaceReduceToTwoComponents,
"BeforeSplitting",
"Reduce to two components",
true);
static SwitchOption interfaceReduceToTwoComponentsNo
(interfaceReduceToTwoComponents,
"AfterSplitting",
"Treat as three components",
false);
static Command<ClusterHadronizationHandler> interfaceUseHandlersForInteraction
("UseHandlersForInteraction",
"Use the current set of handlers for the given interaction.",
&ClusterHadronizationHandler::_setHandlersForInteraction, false);
}
void ClusterHadronizationHandler::rebind(const TranslationMap & trans) {
for (auto iter : _partonSplitters) {
_partonSplitters[iter.first] = trans.translate(iter.second);
}
for (auto iter : _clusterFinders) {
_clusterFinders[iter.first] = trans.translate(iter.second);
}
for (auto iter : _colourReconnectors) {
_colourReconnectors[iter.first] = trans.translate(iter.second);
}
for (auto iter : _clusterFissioners) {
_clusterFissioners[iter.first] = trans.translate(iter.second);
}
for (auto iter : _lightClusterDecayers) {
_lightClusterDecayers[iter.first] = trans.translate(iter.second);
}
for (auto iter : _clusterDecayers) {
_clusterDecayers[iter.first] = trans.translate(iter.second);
}
for (auto iter : _gluonMassGenerators) {
_gluonMassGenerators[iter.first] = trans.translate(iter.second);
}
HandlerBase::rebind(trans);
}
IVector ClusterHadronizationHandler::getReferences() {
IVector ret = HandlerBase::getReferences();
for (auto iter : _partonSplitters) {
ret.push_back(iter.second);
}
for (auto iter : _clusterFinders) {
ret.push_back(iter.second);
}
for (auto iter : _colourReconnectors) {
ret.push_back(iter.second);
}
for (auto iter : _clusterFissioners) {
ret.push_back(iter.second);
}
for (auto iter : _lightClusterDecayers) {
ret.push_back(iter.second);
}
for (auto iter : _clusterDecayers) {
ret.push_back(iter.second);
}
for (auto iter : _gluonMassGenerators) {
ret.push_back(iter.second);
}
return ret;
}
void ClusterHadronizationHandler::doinit() {
HadronizationHandler::doinit();
// put current handlers into action for QCD to guarantee default behaviour
if ( _partonSplitters.empty() ) {
_partonSplitters[PDT::ColouredQCD] = _partonSplitter;
_clusterFinders[PDT::ColouredQCD] = _clusterFinder;
_colourReconnectors[PDT::ColouredQCD] = _colourReconnector;
_clusterFissioners[PDT::ColouredQCD] = _clusterFissioner;
_lightClusterDecayers[PDT::ColouredQCD] = _lightClusterDecayer;
_clusterDecayers[PDT::ColouredQCD] = _clusterDecayer;
_gluonMassGenerators[PDT::ColouredQCD] = _gluonMassGenerator;
}
}
string ClusterHadronizationHandler::_setHandlersForInteraction(string interaction) {
int interactionId = -1;
if ( interaction == " QCD" ) {
interactionId = PDT::ColouredQCD;
} else if ( interaction == " Dark" ) {
interactionId = PDT::ColouredDark;
} else {
return "unknown interaction " + interaction;
}
_partonSplitters[interactionId] = _partonSplitter;
_clusterFinders[interactionId] = _clusterFinder;
_colourReconnectors[interactionId] = _colourReconnector;
_clusterFissioners[interactionId] = _clusterFissioner;
_lightClusterDecayers[interactionId] = _lightClusterDecayer;
_clusterDecayers[interactionId] = _clusterDecayer;
_gluonMassGenerators[interactionId] = _gluonMassGenerator;
return "";
}
namespace {
void extractChildren(tPPtr p, set<PPtr> & all) {
if (p->children().empty()) return;
for (PVector::const_iterator child = p->children().begin();
child != p->children().end(); ++child) {
all.insert(*child);
extractChildren(*child, all);
}
}
}
void ClusterHadronizationHandler::
handle(EventHandler & ch, const tPVector & tagged,
const Hint &) {
useMe();
currentHandler_ = this;
PVector theList(tagged.begin(),tagged.end());
// set the scale for coloured particles to just above the gluon mass squared
// if less than this so they are classed as perturbative
// TODO: Should this be hard-coded here, or defined in inputs?
map<int,int> interactions = {{PDT::ColouredQCD, ParticleID::g}, {PDT::ColouredDark, ParticleID::darkg}};
// PVector currentlist(tagged.begin(),tagged.end());
map<int,PVector> currentlists;
for ( const auto & p : theList ) {
for ( auto i : interactions ) {
if ( p->data().colouredInteraction() == i.first ) {
currentlists[i.first].push_back(p);
Energy2 Q02 = 1.01*sqr(getParticleData(i.second)->constituentMass());
if(currentlists[i.first].back()->scale()<Q02) currentlists[i.first].back()->scale(Q02);
}
}
}
map<int,ClusterVector> clusters;
// ATTENTION this needs to have changes to include the new hadron spectrum functionality (gluon mass generator)
for ( auto i : interactions ) {
// reshuffle to the constituents
if ( reshuffle_ ) {
vector<PVector> reshufflelists;
if ( reshuffle_ == 1 ) { // global reshuffling
reshufflelists.push_back(currentlists[i.first]);
}
else if ( reshuffle_ == 2 ) {// colour connected reshuffling
splitIntoColourSinglets(currentlists[i.first], reshufflelists, i.first);
}
Ptr<GluonMassGenerator>::ptr gMassGenerator;
auto git = _gluonMassGenerators.find(i.first);
if ( git != _gluonMassGenerators.end() )
gMassGenerator = git->second;
for (auto currentlist : reshufflelists){
// get available energy and energy needed for constituent mass shells
LorentzMomentum totalQ;
Energy needQ = ZERO;
size_t nGluons = 0; // number of gluons for which a mass need be generated
for ( auto p : currentlist ) {
totalQ += p->momentum();
if ( p->id() == i.second && gMassGenerator ) {
++nGluons;
needQ += gMassGenerator->minGluonMass();
continue;
}
needQ += p->dataPtr()->constituentMass();
}
Energy Q = totalQ.m();
if ( needQ > Q )
throw Exception() << "cannot reshuffle to constituent mass shells" << Exception::eventerror;
// generate gluon masses if needed
list<Energy> gluonMasses;
if ( nGluons && gMassGenerator )
gluonMasses = gMassGenerator->generateMany(nGluons,Q-needQ);
// set masses for inidividual particles
vector<Energy> masses;
for ( auto p : currentlist ) {
if ( p->id() == i.second && gMassGenerator ) {
list<Energy>::const_iterator it = gluonMasses.begin();
advance(it,UseRandom::irnd(gluonMasses.size()));
masses.push_back(*it);
gluonMasses.erase(it);
}
else {
masses.push_back(p->dataPtr()->constituentMass());
}
}
// reshuffle to new masses
reshuffle(currentlist,masses);
}
}
// split the gluons
if ( !currentlists[i.first].empty() ) {
assert(_partonSplitters.find(i.first) != _partonSplitters.end());
_partonSplitters[i.first]->split(currentlists[i.first]);
}
// form the clusters
if ( !currentlists[i.first].empty() ) {
assert(_clusterFinders.find(i.first) != _clusterFinders.end());
clusters[i.first] = _clusterFinders[i.first]->formClusters(currentlists[i.first]);
// reduce BV clusters to two components now if needed
if(_reduceToTwoComponents)
_clusterFinders[i.first]->reduceToTwoComponents(clusters[i.first]);
}
}
// perform colour reconnection if needed and then
// decay the clusters into one hadron
for ( auto i : interactions ) {
if ( clusters[i.first].empty() )
continue;
bool lightOK = false;
short tried = 0;
const ClusterVector savedclusters = clusters[i.first];
tPVector finalHadrons; // only needed for partonic decayer
while (!lightOK && tried++ < 10) {
// no colour reconnection with baryon-number-violating (BV) clusters
ClusterVector CRclusters, BVclusters;
CRclusters.reserve( clusters[i.first].size() );
BVclusters.reserve( clusters[i.first].size() );
for (size_t ic = 0; ic < clusters[i.first].size(); ++ic) {
ClusterPtr cl = clusters[i.first].at(ic);
bool hasClusterParent = false;
for (unsigned int ix=0; ix < cl->parents().size(); ++ix) {
if (cl->parents()[ix]->id() == ParticleID::Cluster) {
hasClusterParent = true;
break;
}
}
if (cl->numComponents() > 2 || hasClusterParent) BVclusters.push_back(cl);
else CRclusters.push_back(cl);
}
// colour reconnection
assert(_colourReconnectors.find(i.first) != _colourReconnectors.end());
_colourReconnectors[i.first]->rearrange(CRclusters);
// tag new clusters as children of the partons to hadronize
_setChildren(CRclusters);
// forms diquarks
// we should already have used _clusterFinder[i.first]
_clusterFinders[i.first]->reduceToTwoComponents(CRclusters);
// recombine vectors of (possibly) reconnected and BV clusters
clusters[i.first].clear();
clusters[i.first].insert( clusters[i.first].end(), CRclusters.begin(), CRclusters.end() );
clusters[i.first].insert( clusters[i.first].end(), BVclusters.begin(), BVclusters.end() );
// fission of heavy clusters
// NB: during cluster fission, light hadrons might be produced straight away
assert(_clusterFissioners.find(i.first) != _clusterFissioners.end());
finalHadrons = _clusterFissioners[i.first]->fission(clusters[i.first],i.first == PDT::ColouredQCD ? isSoftUnderlyingEventON() : false);
// if clusters not previously reduced to two components do it now
if(!_reduceToTwoComponents)
_clusterFinders[i.first]->reduceToTwoComponents(clusters[i.first]);
assert(_lightClusterDecayers.find(i.first) != _lightClusterDecayers.end());
lightOK = _lightClusterDecayers[i.first]->decay(clusters[i.first],finalHadrons);
// if the decay of the light clusters was not successful, undo the cluster
// fission and decay steps and revert to the original state of the event
// record
if (!lightOK) {
- // std::cout << "failed LightClusterDecayer " << std::endl;
- clusters[i.first] = savedclusters;
- for_each(clusters[i.first].begin(),
- clusters[i.first].end(),
- std::mem_fn(&Particle::undecay));
+ std::cout << "failed LightClusterDecayer " << std::endl;
+ clusters[i.first] = savedclusters;
+ for_each(clusters[i.first].begin(),
+ clusters[i.first].end(),
+ std::mem_fn(&Particle::undecay));
}
}
if (!lightOK) {
throw Exception( "CluHad::handle(): tried LightClusterDecayer 10 times!",
Exception::eventerror);
}
// decay the remaining clusters
assert(_clusterDecayers[i.first]);
_clusterDecayers[i.first]->decay(clusters[i.first],finalHadrons);
}
// *****************************************
// *****************************************
// *****************************************
bool finalStateCluster=false;
StepPtr pstep = newStep();
set<PPtr> allDecendants;
for (tPVector::const_iterator it = tagged.begin();
it != tagged.end(); ++it) {
extractChildren(*it, allDecendants);
}
for(set<PPtr>::const_iterator it = allDecendants.begin();
it != allDecendants.end(); ++it) {
// this is a workaround because the set sometimes
// re-orders parents after their children
if ((*it)->children().empty()){
// If there is a cluster in the final state throw an event error
if((*it)->id()==81) {
finalStateCluster=true;
}
pstep->addDecayProduct(*it);
}
else {
pstep->addDecayProduct(*it);
pstep->addIntermediate(*it);
}
}
// For very small center of mass energies it might happen that baryonic clusters cannot decay into hadrons
if (finalStateCluster){
throw Exception( "CluHad::Handle(): Cluster in the final state",
Exception::eventerror);
}
// *****************************************
// *****************************************
// *****************************************
// soft underlying event if needed
if (isSoftUnderlyingEventON()) {
assert(_underlyingEventHandler);
ch.performStep(_underlyingEventHandler,Hint::Default());
}
}
// Sets parent child relationship of all clusters with two components
// Relationships for clusters with more than two components are set elsewhere in the Colour Reconnector
void ClusterHadronizationHandler::_setChildren(const ClusterVector & clusters) const {
// erase existing information about the partons' children
tPVector partons;
for ( const auto & cl : clusters ) {
if ( cl->numComponents() > 2 ) continue;
partons.push_back( cl->colParticle() );
partons.push_back( cl->antiColParticle() );
}
// erase all previous information about parent child relationship
for_each(partons.begin(), partons.end(), std::mem_fn(&Particle::undecay));
// give new parents to the clusters: their constituents
for ( const auto & cl : clusters ) {
if ( cl->numComponents() > 2 ) continue;
cl->colParticle()->addChild(cl);
cl->antiColParticle()->addChild(cl);
}
}
void ClusterHadronizationHandler::splitIntoColourSinglets(PVector copylist,
vector<PVector>& reshufflelists,
int){
PVector currentlist;
bool gluonloop;
PPtr firstparticle, temp;
reshufflelists.clear();
while (copylist.size()>0){
gluonloop=false;
currentlist.clear();
firstparticle=copylist.back();
copylist.pop_back();
if (!firstparticle->coloured()){
continue; //non-coloured particles are not included
}
currentlist.push_back(firstparticle);
//go up the anitColourLine and check if we are in a gluon loop
temp=firstparticle;
while( temp->hasAntiColour()){
temp = temp->antiColourLine()->endParticle();
if(temp==firstparticle){
gluonloop=true;
break;
}
else{
currentlist.push_back(temp);
copylist.erase(remove(copylist.begin(),copylist.end(), temp), copylist.end());
}
}
//if not a gluon loop, go up the ColourLine
if(!gluonloop){
temp=firstparticle;
while( temp->hasColour()){
temp=temp->colourLine()->startParticle();
currentlist.push_back(temp);
copylist.erase(remove(copylist.begin(),copylist.end(), temp), copylist.end());
}
}
reshufflelists.push_back(currentlist);
}
}
diff --git a/Hadronization/ColourReconnector.cc b/Hadronization/ColourReconnector.cc
--- a/Hadronization/ColourReconnector.cc
+++ b/Hadronization/ColourReconnector.cc
@@ -1,3382 +1,3392 @@
// -*- C++ -*-
//
// ColourReconnector.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the ColourReconnector class.
//
#include "ColourReconnector.h"
#include "Cluster.h"
#include <ThePEG/Utilities/DescribeClass.h>
#include <ThePEG/Repository/UseRandom.h>
#include <ThePEG/PDT/StandardMatchers.h>
#include <ThePEG/Persistency/PersistentOStream.h>
#include <ThePEG/Persistency/PersistentIStream.h>
#include <ThePEG/Interface/Switch.h>
#include <ThePEG/Interface/Reference.h>
#include <ThePEG/Interface/Parameter.h>
#include "Herwig/Utilities/Maths.h"
#include "Herwig/Utilities/expm-1.h"
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
using namespace Herwig;
using CluVecIt = ColourReconnector::CluVecIt;
using Constants::pi;
using Constants::twopi;
DescribeClass<ColourReconnector,Interfaced>
describeColourReconnector("Herwig::ColourReconnector","Herwig.so");
IBPtr ColourReconnector::clone() const {
return new_ptr(*this);
}
IBPtr ColourReconnector::fullclone() const {
return new_ptr(*this);
}
void ColourReconnector::rearrange(ClusterVector & clusters) {
if (_clreco == 0) return;
// need at least two clusters
if (clusters.size() < 2) return;
// std::cout << "size before = "<< clusters.size() << std::endl;
// do the colour reconnection
switch (_algorithm) {
case 0:
_doRecoPlain(clusters);
break;
case 1:
_doRecoStatistical(clusters);
break;
case 2:
_doRecoBaryonic(clusters);
break;
case 3:
_doRecoBaryonicMesonic(clusters);
break;
case 4:
_doRecoBaryonicDiquarkCluster(clusters);
break;
}
// std::cout << "size after = "<< clusters.size() << std::endl;
}
namespace {
double calculateRapidityRF(const Lorentz5Momentum & q1,
const Lorentz5Momentum & p2) {
//calculate rapidity wrt the direction of q1
//angle between the particles in the RF of cluster of q1
// calculate the z component of p2 w.r.t the direction of q1
if(q1.rho2()==ZERO) return 0.;
const Energy pz = p2.vect() * q1.vect().unit();
if ( pz == ZERO ) return 0.;
// Transverse momentum of p2 w.r.t the direction of q1
const Energy pt = sqrt(p2.vect().mag2() - sqr(pz));
// Transverse mass pf p2 w.r.t to the direction of q1
const Energy mtrans = sqrt(p2.mass()*p2.mass() + (pt*pt));
// Correct formula
const double y2 = log((p2.t() + abs(pz))/mtrans);
return ( pz < ZERO ) ? -y2 : y2;
}
}
Energy2 ColourReconnector::_clusterMassSum(const PVector & q,
const PVector & aq) const {
const size_t nclusters = q.size();
assert (aq.size() == nclusters);
Energy2 sum = ZERO;
for (size_t i = 0; i < nclusters; i++)
sum += ( q[i]->momentum() + aq[i]->momentum() ).m2();
return sum;
}
double ColourReconnector::_displacement(tcPPtr p, tcPPtr q) const {
double deltaRap = (p->rapidity() - q->rapidity());
double deltaPhi = fabs(p->momentum().phi() - q->momentum().phi());
// keep deltaPhi's below Pi due to periodicity
if (deltaPhi > M_PI) deltaPhi-=M_PI;
return sqrt(deltaRap * deltaRap + deltaPhi * deltaPhi);
}
/**
* Computes circular Mean of three angles alpha, beta, gamma
* */
static double circularMean(double alpha, double beta, double gamma) {
double xMean=(cos(alpha)+cos(beta)+cos(gamma))/3.0;
double yMean=(sin(alpha)+sin(beta)+sin(gamma))/3.0;
// to make the function fail-save
if (xMean==0 && yMean==0) return M_PI;
return atan2(yMean,xMean);
}
namespace {
// ColourFlow for 3 colour flows for baryon state
// ordered permutations by one transposition
// sign of permutations # of difference
// |1> = |123> + 0
// |2> = |132> - 1
// |3> = |213> - 1
// |4> = |231> + 2
// |5> = |312> + 2
// |6> = |321> - 1
int signPermutationState(int i);
int signPermutationState(int i)
{
if (i==0 || i==3 || i==4) return 1;
else if (i==1 || i==2 || i==5) return -1;
else assert(false);
}
// ColourFlow scalar product matrix for 3
// colour flows in the following basis
// |1> = |123> + 0
// |2> = |132> - 1
// |3> = |213> - 1
// |4> = |231> + 2
// |5> = |312> + 2
// |6> = |321> - 1
unsigned int scalarProducts(int i, int j);
unsigned int scalarProducts(int i, int j)
{
// Verified for i,j < 3! and 3 colour flows
if (i>j) return scalarProducts(j,i);
unsigned int Nc=3;
if (i==j) return Nc*Nc*Nc;
switch(i)
{
case 0:
{
if (j==1 || j==2 || j==5) return Nc*Nc;
else if (j==3 || j==4) return Nc;
else return Nc*Nc*Nc;
break;
}
case 1:
{
if (j==0 || j==3 || j==4) return Nc*Nc;
else if (j==2 || j==5) return Nc;
else return Nc*Nc*Nc;
break;
}
case 2:
{
if (j==0 || j==3 || j==4) return Nc*Nc;
else if (j==1 || j==5) return Nc;
else return Nc*Nc*Nc;
break;
}
case 3:
{
if (j==1 || j==2 || j==5) return Nc*Nc;
else if (j==0 || j==4) return Nc;
else return Nc*Nc*Nc;
break;
}
case 4:
{
if (j==1 || j==2 || j==5) return Nc*Nc;
else if (j==0 || j==3) return Nc;
else return Nc*Nc*Nc;
break;
}
case 5:
{
if (j==0 || j==3 || j==4) return Nc*Nc;
else if (j==1 || j==2) return Nc;
else return Nc*Nc*Nc;
break;
}
}
return Nc;
}
}
std::unordered_map<int,double> ColourReconnector::_reconnectionAmplitudesCF2 (const ClusterPtr & c1, const ClusterPtr & c2) const{
// Verified according to convention of analytics/matrices2_BCR.nb and not
// according to https://arxiv.org/abs/1808.06770
// The same convention as in https://arxiv.org/abs/1808.06770 can be obtained by
// the mapping 1->1;2->4;3->2;4->3
Lorentz5Momentum p1 = c1->colParticle()->momentum();
Lorentz5Momentum p2 = c1->antiColParticle()->momentum();
Lorentz5Momentum p3 = c2->colParticle()->momentum();
Lorentz5Momentum p4 = c2->antiColParticle()->momentum();
double M12 = (p1 + p2).m2()/sqr(_dynamicCRscale);
double M24 = (p2 + p4).m2()/sqr(_dynamicCRscale);
double M13 = (p1 + p3).m2()/sqr(_dynamicCRscale);
double M23 = (p2 + p3).m2()/sqr(_dynamicCRscale);
double M14 = (p1 + p4).m2()/sqr(_dynamicCRscale);
double M34 = (p3 + p4).m2()/sqr(_dynamicCRscale);
double alphaQCD=_dynamicCRalphaS;
double logSqrOmega12=alphaQCD*pow(log(M12),2)/(2.0*twopi);
double logSqrOmega24=alphaQCD*pow(log(M24),2)/(2.0*twopi);
double logSqrOmega13=alphaQCD*pow(log(M13),2)/(2.0*twopi);
double logSqrOmega23=alphaQCD*pow(log(M23),2)/(2.0*twopi);
double logSqrOmega14=alphaQCD*pow(log(M14),2)/(2.0*twopi);
double logSqrOmega34=alphaQCD*pow(log(M34),2)/(2.0*twopi);
double a = (logSqrOmega34 + logSqrOmega12)/2.0;
double b = (logSqrOmega14 + logSqrOmega23)/2.0;
double c = (logSqrOmega13 + logSqrOmega24)/2.0;
double sqrtDelta=sqrt(9*a*a-4*c*(a+b)-14*a*b+9*b*b+4*c*c);
double U11=sqrtDelta/tanh(sqrtDelta/2.0)+3*(b-a);
double U21=2*(c-b);
double Nc=3.0;
double TransAmpNoCR = Nc*(Nc * U11 + U21);
double TransAmpMesonicCR = Nc*(U11 + Nc * U21);
std::unordered_map<int,double> amplitudes;
// No Colour Reconnection amplitude
amplitudes[123] = TransAmpNoCR;
// Mesonic Colour Reconnection amplitude
amplitudes[213] = TransAmpMesonicCR;
return amplitudes;
}
std::tuple<double,double,double> ColourReconnector::_dynamicRecoProbabilitiesCF2(const ClusterPtr & c1, const ClusterPtr & c2, bool diquarkCR) const{
std::unordered_map<int,double> amplitudes = _reconnectionAmplitudesCF2 (c1, c2);
double pNoCR;
double pMesonicCR;
double pDiquarkCR = 0.0;
double TransAmpNoCR = sqr(amplitudes[123]);
double TransAmpMesonicCR = sqr(amplitudes[213]);
double sum = TransAmpNoCR + TransAmpMesonicCR;
if (diquarkCR && _canMakeDiquarkCluster(c1,c2)) {
double ND = 2.0/sqrt(3.0); // sqrt(2*(Nc-1.0)/Nc);
double TransAmpDiquarkCR = sqr((amplitudes[123]-amplitudes[213])/ND);
sum += TransAmpDiquarkCR;
pDiquarkCR = TransAmpDiquarkCR/sum;
assert( pDiquarkCR<=1.0 && pDiquarkCR>=0.0);
}
pNoCR = TransAmpNoCR/sum;
pMesonicCR = TransAmpMesonicCR/sum;
assert( pNoCR<=1.0 && pNoCR>=0.0);
assert( pMesonicCR<=1.0 && pMesonicCR>=0.0);
if (_debug)
{
Lorentz5Momentum p1col = (c1)->colParticle()->momentum();
Lorentz5Momentum p1acol = (c1)->antiColParticle()->momentum();
Lorentz5Momentum p2col = (c2)->colParticle()->momentum();
Lorentz5Momentum p2acol = (c2)->antiColParticle()->momentum();
const Boost boostv1=-c1->momentum().boostVector();
const Boost boostv2=-c2->momentum().boostVector();
p1col .boost(boostv1);
p1acol.boost(boostv1);
p2col .boost(boostv1);
p2acol.boost(boostv1);
double rap1c= calculateRapidityRF(p1acol,p2col);
double rap1a= calculateRapidityRF(p1acol,p2acol);
double pT1c = p2col.vect() .perp(p1acol.vect())/GeV;
double pT1a = p2acol.vect().perp(p1acol.vect())/GeV;
p1col .boost(-boostv1);
p1acol.boost(-boostv1);
p2col .boost(-boostv1);
p2acol.boost(-boostv1);
p1col .boost(boostv2);
p1acol.boost(boostv2);
p2col .boost(boostv2);
p2acol.boost(boostv2);
double rap2c= calculateRapidityRF(p2acol,p1col);
double rap2a= calculateRapidityRF(p2acol,p1acol);
double pT2c = p1col.vect() .perp(p1acol.vect())/GeV;
double pT2a = p1acol.vect().perp(p1acol.vect())/GeV;
// calculate the rapidity of the other constituents of the clusters
// w.r.t axis of p1anticol.vect.unit
// const double rapq = calculateRapidityRF(p1anticol,p2col);
ofstream out("WriteOut/kinematicRecoProbability.dat", std::ios::app);
out << pNoCR << "\t" << pMesonicCR << "\t" << pDiquarkCR << "\t";
out << rap1c << "\t" << rap1a << "\t" << pT1c << "\t" << pT1a << "\t";
out << rap2c << "\t" << rap2a << "\t" << pT2c << "\t" << pT2a << "\n";
out.close();
}
return {pNoCR, pMesonicCR, pDiquarkCR};
}
int ColourReconnector::_stateToPermutation(const int i) const {
switch (i)
{
case 0:
// Baryonic CR
return 0;
// Mesonic CR's
case 1:
return 123;
case 2:
return 132;
case 3:
return 213;
case 4:
return 231;
case 5:
return 312;
case 6:
return 321;
default:
assert(false);
}
}
Selector<int> ColourReconnector::_selector(const ClusterVector & clusters, bool diquarkCR) const {
switch (clusters.size())
{
case 2:
return getProbabilities2CF(clusters[0], clusters[1], diquarkCR);
break;
case 3:
{
if (_dynamicCR) {
return _selectorCF3(clusters, diquarkCR);
// return _selectorCF3(clusters, false);
}
else
{
Selector<int> CRoptions;
double sum=_preco+_precoBaryonic;
const int NpossibilitiesMCR = 6;
for (int i = 0; i < NpossibilitiesMCR; i++)
{
CRoptions.insert(_preco, _stateToPermutation(i+1));
}
CRoptions.insert(_precoBaryonic, 0);
if (diquarkCR) {
const int N=3;
bool first=false;
// Here i is the index of the quark and j of the antiquark
// which will be connected to each other (other partons will be
// forming the diquark cluster if kinematically viable)
double PhaseSpace;
for (int i = 1; i <= N; i++) {
for (int j = 1; j <= N; j++) {
if (_canMakeDiquarkCluster(
clusters[i%3]->colParticle(), clusters[(i+1)%3]->colParticle(),
clusters[j%3]->antiColParticle(), clusters[(j+1)%3]->antiColParticle(),PhaseSpace)){
CRoptions.insert(_precoDiquark*PhaseSpace, -(10*i+j));
if (first){
sum+=_precoDiquark*PhaseSpace;
first=false;
}
}
}
}
}
// no CR probability
CRoptions.insert((1-sum) > 0 ? (1-sum):0, _stateToPermutation(1));
return CRoptions;
}
}
default:
assert(false);
}
}
Selector<int> ColourReconnector::_selectorCF3(const ClusterVector & clusters, bool diquarkCR) const {
std::unordered_map<int, double> amplitudes = _reconnectionAmplitudesSGE(clusters);
double sum2 = 0;
double sumBaryon = 0;
double NB=2.0/sqrt(3.0);
Selector<int> CRoptions;
double amp2;
const int NpossibilitiesMCR = 6;
for (int i = 0; i < NpossibilitiesMCR; i++)
{
sumBaryon += amplitudes[_stateToPermutation(i+1)]*signPermutationState(i)/NB;
amp2 = sqr(amplitudes[_stateToPermutation(i+1)]);
sum2 += amp2;
CRoptions.insert(amp2, _stateToPermutation(i+1));
}
sum2+=sumBaryon*sumBaryon;
if (diquarkCR) {
double ND=2.0/sqrt(3.0); // sqrt(2*(Nc-1.0)/Nc)
const int N=3;
double NDiqFACT=1.0; // sqrt(2*(Nc-1.0)/Nc)
double PhaseSpace;
int perm1;
int perm2;
// Here i is the index of the quark and j of the antiquark
// which will be connected to each other (other partons will be
// forming the diquark cluster if kinematically viable)
for (int i = 1; i <= N; i++) {
for (int j = 1; j <= N; j++) {
// TODO add here a check if we can make Diquarkclsuter
// std::cout << "i = "<<i-1 <<" j = "<<j-1 << std::endl;
// std::cout << "c1= "<<i%3 <<" a1= "<<j%3 << std::endl;
// std::cout << "c2= "<<(i+1)%3 <<" a2= "<<(j+1)%3 << std::endl;
assert(i%3!=(i-1));
assert((i+1)%3!=(i-1));
assert(j%3!=(j-1));
assert((j+1)%3!=(j-1));
if (_canMakeDiquarkCluster(
clusters[i%3]->colParticle(), clusters[(i+1)%3]->colParticle(),
clusters[j%3]->antiColParticle(), clusters[(j+1)%3]->antiColParticle() ,PhaseSpace))
{
switch (i)
{
case 1:
{
perm1 = 100*j + 10*(((j+1)%3)+1) + (((j)%3)+1);
perm2 = 100*j + 10*(((j)%3)+1) + (((j+1)%3)+1);
amp2 = sqr((amplitudes[perm1] - amplitudes[perm2])/ND);
amp2*=NDiqFACT;
sum2 += amp2;
if (amp2==0 || amplitudes[perm2]==0 || amplitudes[perm1]==0 ){
std::cout << "amp23 = "<< amp2 << std::endl;
std::cout << "perm13 = "<< perm1 << std::endl;
std::cout << "perm23 = "<< perm2 << std::endl;
}
// std::cout << "Mesonic CR = "<< i<<j << std::endl;
// std::cout << "Diquark CR = "<< i%3+1<<(i+1)%3 +1<<j%3+1<<(j+1)%3 +1<< std::endl;
// std::cout << "perm1 = "<< perm1 << std::endl;
// std::cout << "perm2 = "<< perm2 << std::endl;
assert((((j)%3)+1)!=j);
assert((((j+1)%3)+1)!=j);
CRoptions.insert(PhaseSpace*amp2, -(10*i+j));
break;
}
case 2:
{
perm1 = 100*(((j)%3)+1) + 10*j + (((j+1)%3)+1);
perm2 = 100*(((j+1)%3)+1) + 10*j + (((j)%3)+1);
amp2 = sqr((amplitudes[perm1] - amplitudes[perm2])/ND);
amp2*=NDiqFACT;
sum2 += amp2;
if (amp2==0 || amplitudes[perm2]==0 || amplitudes[perm1]==0 ){
std::cout << "amp22 = "<< amp2 << std::endl;
std::cout << "perm12 = "<< perm1 << std::endl;
std::cout << "perm22 = "<< perm2 << std::endl;
}
// std::cout << "Mesonic CR = "<< i<<j << std::endl;
// std::cout << "Diquark CR = "<< i%3+1<<(i+1)%3 +1<<j%3+1<<(j+1)%3 +1<< std::endl;
// std::cout << "perm1 = "<< perm1 << std::endl;
// std::cout << "perm2 = "<< perm2 << std::endl;
assert((((j)%3)+1)!=j);
assert((((j+1)%3)+1)!=j);
CRoptions.insert(PhaseSpace*amp2, -(10*i+j));
break;
}
case 3:
{
perm1 = 100*(((j+1)%3)+1) + 10*(((j)%3)+1) + j;
perm2 = 100*(((j)%3)+1) + 10*(((j+1)%3)+1) + j;
amp2 = sqr((amplitudes[perm1] - amplitudes[perm2])/ND);
amp2*=NDiqFACT;
sum2 += amp2;
if (amp2==0 || amplitudes[perm2]==0 || amplitudes[perm1]==0 ){
std::cout << "amp21 = "<< amp2 << std::endl;
std::cout << "perm11 = "<< perm1 << std::endl;
std::cout << "perm21 = "<< perm2 << std::endl;
}
// std::cout << "Mesonic CR = "<< i<<j << std::endl;
// std::cout << "Diquark CR = "<< i%3+1<<(i+1)%3 +1<<j%3+1<<(j+1)%3 +1<< std::endl;
// std::cout << "perm1 = "<< perm1 << std::endl;
// std::cout << "perm2 = "<< perm2 << std::endl;
assert((((j)%3)+1)!=j);
assert((((j+1)%3)+1)!=j);
CRoptions.insert(PhaseSpace*amp2, -(10*i+j));
break;
}
default:
assert(false);
}
// std::cout << CRoptions << std::endl;
}
}
}
}
// double BaryonRate = sumBaryon*sumBaryon/sum2;
// std::cout << "BaryonRate "<<BaryonRate <<"\n";
// std::cout << "NoRecoRate "<<sqr(amplitudes[_stateToPermutation(1)])/sum2 <<"\n";
// std::cout << "MeRecoRate "<<1-BaryonRate-sqr(amplitudes[_stateToPermutation(1)])/sum2 <<"\n";
CRoptions.insert(sumBaryon*sumBaryon, _stateToPermutation(0));
//TODO insert diquarks
return CRoptions;
}
namespace {
inline int hasDiquark(const ClusterPtr & cit) {
int res = 0;
for (unsigned int i = 0; i<(cit)->numComponents(); i++) {
if (DiquarkMatcher::Check(*((cit)->particle(i)->dataPtr()))) {
res++;
}
}
return res;
}
}
std::unordered_map<int, double> ColourReconnector::_reconnectionAmplitudesSGE(const ClusterVector & clusters) const {
std::unordered_map<int, double> amplitudes;
int size=clusters.size();
assert(clusters.size()<4);
switch(size){
case 2:
{
amplitudes = _reconnectionAmplitudesCF2(clusters[0],clusters[1]);
break;
}
case 3:
{
if (clusters[0]->numComponents()!=2 ||
clusters[1]->numComponents()!=2 ||
clusters[2]->numComponents()!=2 ||
hasDiquark(clusters[0]) ||
hasDiquark(clusters[1]) ||
hasDiquark(clusters[2]) ) {
std::cout << "reject config. should reject before somehow\n";
return amplitudes;
}
const int N=6; // 2*Nclu;
double omIJ[N][N];
double alphaQCD=_dynamicCRalphaS;
Lorentz5Momentum mom_i, mom_j;
// Assumption in notebook analytics/matrices2_BCR.nb
// Ordering of omIJ and scalarProds is according to {clu1_col, clu1_anti, clu2_col, clu2_anti,...}
for (int i = 0; i < N; i++)
{
if ((i+1)%2==1) mom_i=clusters[i/2]->colParticle()->momentum();
else mom_i=clusters[(i-1)/2]->antiColParticle()->momentum();
for (int j = i+1; j < N; j++)
{
if ((j+1)%2==1) mom_j=clusters[j/2]->colParticle()->momentum();
else mom_j=clusters[(j-1)/2]->antiColParticle()->momentum();
omIJ[i][j]=alphaQCD*pow(log((mom_i+mom_j).m2()/sqr(_dynamicCRscale)),2)/(2.0*twopi);
}
}
boost::numeric::ublas::matrix<double> * Uevolve = new boost::numeric::ublas::matrix<double>(N,N);
boost::numeric::ublas::matrix<double> Omega(N,N);
int Nc = 3;
// Verified Omega Matrix with analytics/matrices2_BCR.nb
Omega(0,0) = - 0.5 * Nc * (omIJ[0][1]+omIJ[2][3]+omIJ[4][5]);
Omega(1,0) = 0.5 * (omIJ[2][4]-omIJ[3][4]-omIJ[2][5]+omIJ[3][5]);
Omega(2,0) = 0.5 * (omIJ[0][2]-omIJ[1][2]-omIJ[0][3]+omIJ[1][3]);
Omega(3,0) = 0.0;
Omega(4,0) = 0.0;
Omega(5,0) = 0.5 * (omIJ[0][4]-omIJ[1][4]-omIJ[0][5]+omIJ[1][5]);
Omega(0,1) = - 0.5 * (omIJ[2][3]-omIJ[2][4]-omIJ[3][5]+omIJ[4][5]);
Omega(1,1) = - 0.5 * Nc * (omIJ[0][1]+omIJ[3][4]+omIJ[2][5]);
Omega(2,1) = 0.0;
Omega(3,1) = - 0.5 * (omIJ[0][3]-omIJ[1][3]-omIJ[0][4]+omIJ[1][4]);
Omega(4,1) = 0.5 * (omIJ[0][2]-omIJ[1][2]-omIJ[0][5]+omIJ[1][5]);
Omega(5,1) = 0.0;
Omega(0,2) = - 0.5 * (omIJ[0][1]-omIJ[0][2]-omIJ[1][3]+omIJ[2][3]);
Omega(1,2) = 0.0;
Omega(2,2) = - 0.5 * Nc * (omIJ[1][2]+omIJ[0][3]+omIJ[4][5]);
Omega(3,2) = - 0.5 * (omIJ[1][4]-omIJ[2][4]-omIJ[1][5]+omIJ[2][5]);
Omega(4,2) = 0.5 * (omIJ[0][4]-omIJ[3][4]-omIJ[0][5]+omIJ[3][5]);
Omega(5,2) = 0.0;
Omega(0,3) = 0.0;
Omega(1,3) = - 0.5 * (omIJ[0][1]-omIJ[1][3]-omIJ[0][4]+omIJ[3][4]);
Omega(2,3) = - 0.5 * (omIJ[1][2]-omIJ[2][4]-omIJ[1][5]+omIJ[4][5]);
Omega(3,3) = - 0.5 * Nc * (omIJ[0][3]+omIJ[1][4]+omIJ[2][5]);
Omega(4,3) = 0.0;
Omega(5,3) = 0.5 * (omIJ[0][2]-omIJ[2][3]-omIJ[0][5]+omIJ[3][5]);
Omega(0,4) = 0.0;
Omega(1,4) = - 0.5 * (omIJ[0][1]-omIJ[0][2]-omIJ[1][5]+omIJ[2][5]);
Omega(2,4) = - 0.5 * (omIJ[0][3]-omIJ[0][4]-omIJ[3][5]+omIJ[4][5]);
Omega(3,4) = 0.0;
Omega(4,4) = - 0.5 * Nc * (omIJ[1][2]+omIJ[3][4]+omIJ[0][5]);
Omega(5,4) = 0.5 * (omIJ[1][3]-omIJ[2][3]-omIJ[1][4]+omIJ[2][4]);
Omega(0,5) = - 0.5 * (omIJ[0][1]-omIJ[0][4]-omIJ[1][5]+omIJ[4][5]);
Omega(1,5) = 0.0;
Omega(2,5) = 0.0;
Omega(3,5) = 0.5 * (omIJ[0][2]-omIJ[0][3]-omIJ[2][5]+omIJ[3][5]);
Omega(4,5) = - 0.5 * (omIJ[1][2]-omIJ[1][3]-omIJ[2][4]+omIJ[3][4]);
Omega(5,5) = - 0.5 * Nc * (omIJ[2][3]+omIJ[1][4]+omIJ[0][5]);
*Uevolve=expm_pad(Omega);
// std::cout << *Uevolve << std::endl;
// std::vector<double> Transition1toJ; // |<J|U|1>|^2
double amp1toJ;
for (int J = 0; J < N; J++)
{
// amp1toJ is here for each J transition amplitude 1->J
amp1toJ=0;
for (int i = 0; i < N; i++)
{
// amp1toJ corresponds to the Uevolve operator applied to the |1> state
// and projected by fixed <J|
// The resulting amp1toJ is <J|U|1>
// Should be correct this way
amp1toJ+=(*Uevolve)(i,0)*scalarProducts(i,J);
}
if (std::isnan(amp1toJ) || std::isinf(amp1toJ)){
throw Exception() << "nan or inf transition probability in ColourReconnector::_reconnectionAmplitudesSGE"
<< Exception::runerror;
}
amplitudes[_stateToPermutation(J+1)] = amp1toJ;
}
delete Uevolve;
break;
}
default:
{
throw Exception() << "Found cluster set of "<<size<<" in ColourReconnector::_reconnectionAmplitudesSGE (can only handle 2 or 3 colour Flows)"
<< Exception::runerror;
}
}
return amplitudes;
}
double ColourReconnector::_displacementBaryonic(tcPPtr q1, tcPPtr q2, tcPPtr q3) const {
if (_junctionMBCR) {
/**
* Junction-like option i.e. displacement
* from "junction centre" (mean rapidity/phi)
*/
double rap1=q1->rapidity();
double rap2=q2->rapidity();
double rap3=q3->rapidity();
double phi1=q1->momentum().phi();
double phi2=q2->momentum().phi();
double phi3=q3->momentum().phi();
double meanRap=(rap1 + rap2 + rap3)/3.0;
// Use circularMean for defining a sensible mean of an angle
double meanPhi=circularMean(phi1,phi2,phi3);
double deltaPhi1=fabs(phi1-meanPhi);
double deltaPhi2=fabs(phi2-meanPhi);
double deltaPhi3=fabs(phi3-meanPhi);
// keep deltaPhi's below Pi due to periodicity
if (deltaPhi1>M_PI) deltaPhi1-=M_PI;
if (deltaPhi2>M_PI) deltaPhi2-=M_PI;
if (deltaPhi3>M_PI) deltaPhi3-=M_PI;
double delR;
delR = sqrt( (rap1-meanRap)*(rap1-meanRap) + deltaPhi1*deltaPhi1 );
delR += sqrt( (rap2-meanRap)*(rap2-meanRap) + deltaPhi2*deltaPhi2 );
delR += sqrt( (rap3-meanRap)*(rap3-meanRap) + deltaPhi3*deltaPhi3 );
return delR;
} else {
/* just summing up all possible 2 quark displacements */
return _displacement(q1, q2) + _displacement(q1, q3) + _displacement(q2, q3);
}
}
bool ColourReconnector::_containsColour8(const ClusterVector & cv,
const vector<size_t> & P) const {
assert (P.size() == cv.size());
for (size_t i = 0; i < cv.size(); i++) {
tcPPtr p = cv[i]->colParticle();
tcPPtr q = cv[P[i]]->antiColParticle();
if (_isColour8(p, q)) return true;
}
return false;
}
void ColourReconnector::_doRecoStatistical(ClusterVector & cv) const {
const size_t nclusters = cv.size();
// initially, enumerate (anti)quarks as given in the cluster vector
ParticleVector q, aq;
for (size_t i = 0; i < nclusters; i++) {
q.push_back( cv[i]->colParticle() );
aq.push_back( cv[i]->antiColParticle() );
}
// annealing scheme
Energy2 t, delta;
Energy2 lambda = _clusterMassSum(q,aq);
const unsigned _ntries = _triesPerStepFactor * nclusters;
// find appropriate starting temperature by measuring the largest lambda
// difference in some dry-run random rearrangements
{
vector<Energy2> typical;
for (int i = 0; i < 10; i++) {
const pair <int,int> toswap = _shuffle(q,aq,5);
ParticleVector newaq = aq;
swap (newaq[toswap.first], newaq[toswap.second]);
Energy2 newlambda = _clusterMassSum(q,newaq);
typical.push_back( abs(newlambda - lambda) );
}
t = _initTemp * Math::median(typical);
}
// anneal in up to _annealingSteps temperature steps
for (unsigned step = 0; step < _annealingSteps; step++) {
// For this temperature step, try to reconnect _ntries times. Stop the
// algorithm if no successful reconnection happens.
unsigned nSuccess = 0;
for (unsigned it = 0; it < _ntries; it++) {
// make a random rearrangement
const unsigned maxtries = 10;
const pair <int,int> toswap = _shuffle(q,aq,maxtries);
const int i = toswap.first;
const int j = toswap.second;
// stop here if we cannot find any allowed reconfiguration
if (i == -1) break;
// create a new antiquark vector with the two partons swapped
ParticleVector newaq = aq;
swap (newaq[i], newaq[j]);
// Check if lambda would decrease. If yes, accept the reconnection. If no,
// accept it only with a probability given by the current Boltzmann
// factor. In the latter case we set p = 0 if the temperature is close to
// 0, to avoid division by 0.
Energy2 newlambda = _clusterMassSum(q,newaq);
delta = newlambda - lambda;
double prob = 1.0;
if (delta > ZERO) prob = ( abs(t) < 1e-8*MeV2 ) ? 0.0 : exp(-delta/t);
if (UseRandom::rnd() < prob) {
lambda = newlambda;
swap (newaq, aq);
nSuccess++;
}
}
if (nSuccess == 0) break;
// reduce temperature
t *= _annealingFactor;
}
// construct the new cluster vector
ClusterVector newclusters;
for (size_t i = 0; i < nclusters; i++) {
ClusterPtr cl = new_ptr( Cluster( q[i], aq[i] ) );
newclusters.push_back(cl);
}
swap(newclusters,cv);
return;
}
Selector <int> ColourReconnector::getProbabilities2CF(const ClusterPtr & c1, const ClusterPtr & c2, bool diquarkCR) const{
Selector <int> res;
if (_dynamicCR) {
double pNoCR;
double pMCR;
double pDCR;
std::tie(pNoCR, pMCR, pDCR) = _dynamicRecoProbabilitiesCF2(c1, c2, diquarkCR);
// Mesonic CR
if ( !( _isColour8(c1->colParticle(), c2->antiColParticle())
|| _isColour8(c2->colParticle(), c1->antiColParticle()) ) ) {
res.insert(pMCR, 213);
}
double PhaseSpace;
if (diquarkCR && _canMakeDiquarkCluster(c1,c2,PhaseSpace) ) {
// Diquark CR
res.insert(pDCR*PhaseSpace, -213);
}
// No CR
res.insert(pNoCR, 123);
}
else {
double wMCR = _preco;
// Mesonic CR
res.insert(wMCR, 213);
double sum = wMCR;
double PhaseSpace;
if (diquarkCR && _canMakeDiquarkCluster(c1,c2,PhaseSpace)) {
double wDCR = _precoDiquark*PhaseSpace;
// Diquark CR
res.insert(wDCR, -213);
sum+=wDCR;
}
// No CR
res.insert((1.0-sum)>0 ? (1.0-sum):0.0, 123);
}
return res;
}
void ColourReconnector::_doRecoPlain(ClusterVector & cv) const {
ClusterVector newcv = cv;
ClusterVector deleted; deleted.reserve(cv.size());
// try to avoid systematic errors by randomising the reconnection order
long (*p_irnd)(long) = UseRandom::irnd;
random_shuffle( newcv.begin(), newcv.end(), p_irnd );
// iterate over all clusters
for (CluVecIt cit = newcv.begin(); cit != newcv.end(); cit++) {
// find the cluster which, if reconnected with *cit, would result in the
// smallest sum of cluster masses
// skip the diquark clusters to be deleted later 2->1 cluster
if (find(deleted.begin(), deleted.end(), *cit) != deleted.end())
continue;
// skip diquark clusters
if ((*cit)->numComponents()>2 || (_diquarkCR && hasDiquark(*cit))) continue;
// NB this method returns *cit if no reconnection partner can be found
- CluVecIt candidate = _dynamicCR ? _findRecoPartnerPlainDynamic(cit, newcv, deleted):_findRecoPartnerPlain(cit, newcv, deleted);
+ CluVecIt candidate = _dynamicCR ? _findRecoPartnerPlainDynamic(cit, newcv, deleted, _diquarkCR>0):_findRecoPartnerPlain(cit, newcv, deleted);
// skip this cluster if no possible reshuffling partner can be found
if (candidate == cit) continue;
// accept the reconnection with probability PrecoProb.
ClusterVector cluvec = {*cit, *candidate};
// Selector <int> sel = getProbabilities2CF(*cit, *candidate, _diquarkCR>0);
Selector <int> sel = _selector(cluvec, _diquarkCR>0);
enum Selection2ColourFlows {
NoReconnection = 123,
MesonicReconnection = 213,
DiquarkReconnection = -213,
};
switch (sel.select(UseRandom::rnd()))
{
case MesonicReconnection:
{
// Mesonic Colour Reconnection
pair <ClusterPtr,ClusterPtr> reconnected = _reconnect(*cit, *candidate);
// Replace the clusters in the ClusterVector. The order of the
// colour-triplet partons in the cluster vector is retained here.
// replace *cit by reconnected.first
*cit = reconnected.first;
// replace candidate by reconnected.second
*candidate = reconnected.second;
break;
}
case DiquarkReconnection:
{
ClusterPtr DiqCluster;
if (_makeDiquarkCluster(*cit, *candidate, DiqCluster)){
deleted.push_back(*candidate);
// Note that these must be the cit,candidate and not cluvec
*cit = DiqCluster;
}
break;
}
case NoReconnection:
// No colour Reconnection
break;
default:
assert(false);
}
}
if (deleted.size()==0) {
swap(cv, newcv);
}
else {
// create a new vector of clusters except for the ones which are "deleted" during
// Diquark reconnection
ClusterVector clustervector;
for (const auto & cluster : newcv)
if (find(deleted.begin(),
deleted.end(), cluster) == deleted.end())
clustervector.push_back(cluster);
swap(cv, clustervector);
}
}
void ColourReconnector::_doRecoBaryonicDiquarkCluster(ClusterVector & cv) const {
/*START REVIEW*/
ClusterVector newcv = cv;
ClusterVector deleted; deleted.reserve(cv.size());
// try to avoid systematic errors by randomising the reconnection order
long (*p_irnd)(long) = UseRandom::irnd;
random_shuffle(newcv.begin(), newcv.end(), p_irnd);
Selector <int> sel;
ClusterVector cluvec;
// iterate over all clusters
for (CluVecIt cit = newcv.begin(); cit != newcv.end(); ++cit) {
//avoid clusters already containing diuarks
if (hasDiquark(*cit)) continue;
//skip the cluster to be deleted later 3->2 cluster
if (find(deleted.begin(), deleted.end(), *cit) != deleted.end())
continue;
// Skip all found baryonic and Tetra clusters, this biases the
// algorithm but implementing something like re-reconnection
// is ongoing work
if ((*cit)->numComponents()>=3) continue;
// Find a candidate suitable for reconnection
CluVecIt candidate1, candidate2;
unsigned typeOfReconnection = 0;
// _findPartnerBaryonicDiquarkCluster(cit, newcv,
_findPartnerBaryonicDiquarkClusterTEST(cit, newcv,
typeOfReconnection,
deleted,
candidate1,
candidate2);
switch (typeOfReconnection)
{
case 0:
// No CR found
continue;
case 3:
case 1:
// Mesonic or Diquark CR with 2 Colour Flows
cluvec={*cit,*candidate1};
break;
case 2:
// Baryonic CR with 3 Colour Flows
cluvec={*cit,*candidate1,*candidate2};
break;
default:
assert(false);
}
if (candidate2!=cit) {
cluvec={*cit,*candidate1,*candidate2};
}
else if (candidate1!=cit) {
cluvec={*cit,*candidate1};
}
else {
std::cout << "no CR" << std::endl;
continue;
}
sel = _selector(cluvec, _diquarkCR>0);
if (typeOfReconnection!=0 && sel.empty()){
throw Exception()
<< "No Selection availible"
<< "ColourReconnector::_doRecoBaryonicDiquarkCluster()" << Exception::runerror;
}
int selection = sel.select(UseRandom::rnd());
sel.clear();
assert(sel.empty());
if (selection == 0){
// Baryonic CR (only one option for 3 colourflows)
deleted.push_back(*candidate2);
// Function that does the reconnection from 3 -> 2 clusters
ClusterPtr b1, b2;
_makeBaryonicClusters(*cit,*candidate1,*candidate2, b1, b2);
// Note that these must be the cit,candidate1 and not cluvec
*cit = b1;
*candidate1 = b2;
}
else if (selection > 0) {
// Mesonic CR
int c1 = ((selection/100) % 4) -1;
int c2 = ((selection-(c1+1)*100)/10 % 4) -1;
int c3 = ((selection-(c1+1)*100-(c2+1)*10) % 4) -1;
if ( c1==0
&& c2==1
&& c3==2){
// noCR
continue;
}
assert(c1>=0 && c1<3);
assert(c2>=0 && c2<3);
assert(c3>=0 && c3<3);
assert(c1!=c2 && c1!=c3 && c2 !=c3);
// if last cluster is untouched we do only reconnect the first two clusters
if (c3==2) {
auto reconnected = _reconnect(*cit,*candidate1);
// Note that these must be the cit, candidate1 and not cluvec
*cit = reconnected.first;
*candidate1 = reconnected.second;
}
else {
// Form clusters (0,c1) (1,c2) (2,c3)
const int infoMCR[3] = {c1, c2, c3};
auto reconnected = _reconnect3Mto3M(*cit,*candidate1,*candidate2, infoMCR);
// Note that these must be the cit,candidateX and not cluvec
*cit = std::get<0>(reconnected);
*candidate1 = std::get<1>(reconnected);
*candidate2 = std::get<2>(reconnected);
}
}
else if (selection < 0) {
//TODO DiquarkCR
// We will delete the candidate1 mesonic clusters
// need to check with 2CF solution
// to form a diquark cluster
if (cluvec.size()==3) {
auto reconnected = _reconnect3MtoMD(cluvec, -selection);
*cit = std::get<0>(reconnected);
*candidate1 = std::get<1>(reconnected);
deleted.push_back(*candidate2);
// *candidate2 = std::get<2>(reconnected);
}
else if (cluvec.size()==2){
ClusterPtr DiqCluster;
if (_makeDiquarkCluster(*cit, *candidate1, DiqCluster)){
deleted.push_back(*candidate1);
// Note that these must be the cit,candidate and not cluvec
*cit = DiqCluster;
}
}
else {
assert(false);
}
}
else {
std::cout << "\nError in selection = "<<selection<<"\n" << std::endl;
assert(false);
}
}
// create a new vector of clusters except for the ones which are "deleted" during
// baryonic reconnection
ClusterVector clustervector;
for (const auto & cluster : newcv)
if (find(deleted.begin(),
deleted.end(), cluster) == deleted.end())
clustervector.push_back(cluster);
swap(cv, clustervector);
}
// Implementation of the baryonic reconnection algorithm
void ColourReconnector::_doRecoBaryonic(ClusterVector & cv) const {
ClusterVector newcv = cv;
ClusterVector deleted; deleted.reserve(cv.size());
// try to avoid systematic errors by randomising the reconnection order
long (*p_irnd)(long) = UseRandom::irnd;
random_shuffle(newcv.begin(), newcv.end(), p_irnd);
double ProbabilityMesonic = _preco;
double ProbabilityBaryonic = _precoBaryonic;
// iterate over all clusters
for (CluVecIt cit = newcv.begin(); cit != newcv.end(); ++cit) {
//avoid clusters already containing diuarks
if (hasDiquark(*cit)) continue;
//skip the cluster to be deleted later 3->2 cluster
if (find(deleted.begin(), deleted.end(), *cit) != deleted.end())
continue;
// Skip all found baryonic clusters, this biases the algorithm but implementing
// something like re-reconnection is ongoing work
if ((*cit)->numComponents()>=3) continue;
// Find a candidate suitable for reconnection
CluVecIt baryonic1, baryonic2;
bool isBaryonicCandidate = false;
CluVecIt candidate = _findPartnerBaryonic(cit, newcv,
isBaryonicCandidate,
deleted,
baryonic1, baryonic2);
// skip this cluster if no possible reconnection partner can be found
if ( !isBaryonicCandidate && *candidate==*cit)
continue;
if (_dynamicCR) {
ClusterVector cluvec;
if (isBaryonicCandidate) {
cluvec={*cit,*baryonic1,*baryonic2};
}
else{
cluvec={*cit,*candidate};
}
Selector <int> sel = _selector(cluvec, _diquarkCR>0);
int selection = sel.select(UseRandom::rnd());
// std::cout << "selector print = " << sel <<"\n";
// sel.clear();
// assert(sel.empty());
// 3 aligned meson case
// Normal 2->2 Colour reconnection
if (selection == 0){
// Baryonic CR only one option
deleted.push_back(*baryonic2);
// Function that does the reconnection from 3 -> 2 clusters
ClusterPtr b1, b2;
_makeBaryonicClusters(*cit,*baryonic1,*baryonic2, b1, b2);
*cit = b1;
*baryonic1 = b2;
}
else if (selection > 0) {
//TODO Mesonic CR
if (isBaryonicCandidate) {
int c1 = ((selection/100) % 4) -1;
int c2 = ((selection-(c1+1)*100)/10 % 4) -1;
int c3 = ((selection-(c1+1)*100-(c2+1)*10) % 4) -1;
if ( c1==0
&& c2==1
&& c3==2){
// noCR
continue;
}
assert(c1>=0 && c1<3);
assert(c2>=0 && c2<3);
assert(c3>=0 && c3<3);
assert(c1!=c2 && c1!=c3 && c2 !=c3);
if (c3==2 ) {
auto reconnected = _reconnect(*cit,*baryonic1);
*cit = reconnected.first;
*baryonic1 = reconnected.second;
}
else {
// Form clusters (0,c1) (1,c2) (2,c3)
const int infoMCR[3] = {c1, c2, c3};
auto reconnected = _reconnect3Mto3M(*cit,*baryonic1,*baryonic2, infoMCR);
*cit = std::get<0>(reconnected);
*baryonic1 = std::get<1>(reconnected);
*baryonic2 = std::get<2>(reconnected);
}
}
else {
auto reconnected = _reconnect(*cit,*candidate);
*cit = reconnected.first;
*candidate = reconnected.second;
}
}
else if (selection < 0) {
//TODO DiquarkCR
// We will delete the candidate1 mesonic clusters
// need to check with 2CF solution
// to form a diquark cluster
if (cluvec.size()==3) {
auto reconnected = _reconnect3MtoMD(cluvec, -selection);
*cit = std::get<0>(reconnected);
*baryonic1 = std::get<1>(reconnected);
deleted.push_back(*baryonic2);
}
else if (cluvec.size()==2) {
ClusterPtr DiqCluster;
if (_makeDiquarkCluster(*cit, *candidate, DiqCluster)){
deleted.push_back(*candidate);
// Note that these must be the cit,candidate and not cluvec
*cit = DiqCluster;
}
}
}
else {
std::cout << "\nError in selection = "<<selection<<"\n" << std::endl;
assert(false);
}
}
else {
// 3 aligned meson case
if ( isBaryonicCandidate
&& UseRandom::rnd() < ProbabilityBaryonic ) {
deleted.push_back(*baryonic2);
// Function that does the reconnection from 3 -> 2 clusters
ClusterPtr b1, b2;
_makeBaryonicClusters(*cit, *baryonic1, *baryonic2, b1, b2);
*cit = b1;
*baryonic1 = b2;
// Baryonic2 is easily skipped in the next loop
}
// Normal 2->2 Colour reconnection
if (!isBaryonicCandidate
&& UseRandom::rnd() < ProbabilityMesonic) {
auto reconnected = _reconnect(*cit, *candidate);
*cit = reconnected.first;
*candidate = reconnected.second;
}
}
}
// create a new vector of clusters except for the ones which are "deleted" during
// baryonic reconnection
ClusterVector clustervector;
for (const auto & cluster : newcv)
if (find(deleted.begin(),
deleted.end(), cluster) == deleted.end())
clustervector.push_back(cluster);
swap(cv, clustervector);
}
bool ColourReconnector::_clustersFarApart( const std::vector<CluVecIt> & clu ) const {
int Ncl=clu.size();
assert(Ncl<=3);
if (Ncl==1) {
return false;
} else if (Ncl==2) {
// veto if Clusters further apart than _maxDistance
if (_localCR && ((*clu[0])->vertex().vect()-(*clu[1])->vertex().vect()).mag() > _maxDistance) return true;
// veto if Clusters have negative spacetime difference
if (_causalCR && ((*clu[0])->vertex()-(*clu[1])->vertex()).m() < ZERO) return true;
} else if (Ncl==3) {
// veto if Clusters further apart than _maxDistance
if (_localCR && ((*clu[0])->vertex().vect()-(*clu[1])->vertex().vect()).mag() > _maxDistance) return true;
if (_localCR && ((*clu[1])->vertex().vect()-(*clu[2])->vertex().vect()).mag() > _maxDistance) return true;
if (_localCR && ((*clu[0])->vertex().vect()-(*clu[2])->vertex().vect()).mag() > _maxDistance) return true;
// veto if Clusters have negative spacetime difference
if (_causalCR && ((*clu[0])->vertex()-(*clu[1])->vertex()).m() < ZERO) return true;
if (_causalCR && ((*clu[1])->vertex()-(*clu[2])->vertex()).m() < ZERO) return true;
if (_causalCR && ((*clu[0])->vertex()-(*clu[2])->vertex()).m() < ZERO) return true;
}
return false;
}
void ColourReconnector::_doReco2BeamClusters(ClusterVector & cv) const {
// try other option
tPPtr p1Di=(cv[0])->colParticle();
tPPtr p2Di=(cv[1])->colParticle();
tPPtr p1Q=(cv[0])->antiColParticle();
tPPtr p2Q=(cv[1])->antiColParticle();
double min_dist=_displacement(p1Di,p1Q)+_displacement(p2Di,p2Q);
if ((_displacement(p1Di,p2Q)+_displacement(p1Di,p2Q))<min_dist) {
_reconnect(cv[0],cv[1]);
}
return;
}
void ColourReconnector::_doRecoBaryonicMesonic(ClusterVector & cv) const {
if (cv.size() < 3) {
/*
* if the option _cr2BeamClusters!=0 is chosen then we try to
* colour reconnect the special case of 2 beam clusters with
* probability 1.0 if there is a better _displacement
* */
if( _cr2BeamClusters && cv.size()==2 ) _doReco2BeamClusters(cv);
return;
}
ClusterVector newcv = cv;
newcv.reserve(2*cv.size());
ClusterVector deleted;
deleted.reserve(cv.size());
// counters for numbers of mesons and baryons selected
unsigned num_meson = 0;
unsigned num_baryon = 0;
// vector of selected clusters
std::vector<CluVecIt> sel;
unsigned number_of_tries = _stepFactor*cv.size()*cv.size();
if (number_of_tries<1) number_of_tries=1;
long (*p_irnd)(long) = UseRandom::irnd;
for (unsigned reconnections_tries = 0; reconnections_tries < number_of_tries; reconnections_tries++) {
num_meson = 0;
num_baryon = 0;
// flag if we are able to find a suitable combinations of clusters
bool _found = false;
// Shuffle list of clusters to avoid systematic bias in cluster selection
random_shuffle(newcv.begin(), newcv.end(), p_irnd);
// loop over clustervector to find CR candidates
for (CluVecIt cit = newcv.begin(); cit != newcv.end(); ++cit) {
// skip the clusters to be deleted later from 3->2 cluster CR
if (find(deleted.begin(), deleted.end(), *cit) != deleted.end()) continue;
// avoid clusters already containing diuarks
if (hasDiquark(*cit)) continue;
// add to selection
sel.push_back(cit);
if (_clustersFarApart(sel)) {
// reject far appart CR
// TODO: after discussion maybe to be omitted
sel.pop_back();
continue;
}
bool isMeson=((*cit)->numComponents() == 2);
if ( isMeson && (num_meson ==0|| num_meson==1) && num_baryon ==0) {
num_meson++;
/**
* now we habe either 1 or 2 mesonic clusters and have to continue
*/
continue;
} else if ( isMeson && (num_baryon == 1 || num_meson ==2)) {
num_meson++;
_found = true;
/**
* we have either 3 mesonic or 1 mesonic and 1 baryonic cluster
* and try to colour reconnect
*/
break;
} else if (num_baryon ==0 && num_meson==0) {
num_baryon++;
/**
* now we have 1 baryonic cluster and have to continue
*/
continue;
} else if (num_meson == 2) {
/**
* we already have 2 mesonic clusters and dont want a baryonic one
* since this is an invalid selection
*/
// remove previously added cluster
sel.pop_back();
continue;
} else {
num_baryon++;
_found = true;
/**
* now we have either 2 baryonic clusters or 1 mesonic and 1 baryonic cluster
* and try to colour reconnect
*/
break;
}
}
// added for more efficent rejection if some reco probabilities are 0
if ( _found ) {
// reject MBtoMB candidates if _precoMB_MB=0
if ( _precoMB_MB == 0 && (num_baryon == 1 && num_meson == 1) ) {
_found=false;
}
// reject BbarBto3M candidates if _precoBbarB_3M=0
if ( _precoBbarB_3M== 0 && num_baryon == 2 ) {
bool isBbarBto3Mcandidate=(
(*sel[0])->particle(0)->hasColour() && (*sel[1])->particle(0)->hasColour(true) )
|| ( (*sel[0])->particle(0)->hasColour(true) && (*sel[1])->particle(0)->hasColour() );
if ( isBbarBto3Mcandidate) _found=false;
}
// reject 2Bto2B candidates if _preco2B_2B=0
if ( _preco2B_2B == 0 && num_baryon == 2 ) {
bool is2Bto2Bcandidate=(
(*sel[0])->particle(0)->hasColour() && (*sel[1])->particle(0)->hasColour() )
|| ( (*sel[0])->particle(0)->hasColour(true) && (*sel[1])->particle(0)->hasColour(true) );
if ( is2Bto2Bcandidate ) _found=false;
}
}
// were we able to find a combination?
if (_found==false) {
// clear the selection if we did not find a valid set of clusters
sel.erase(sel.begin(), sel.end());
continue;
}
assert(sel.size()<4);
assert(sel.size()>1);
string kind_of_reco = "";
int reco_info[3];
// find best CR option for the selection
_findbestreconnectionoption(sel, num_baryon, kind_of_reco, reco_info);
if (kind_of_reco == "") {
// no reconnection was found
sel.erase(sel.begin(), sel.end());
continue;
} else if (kind_of_reco == "3Mto3M" && UseRandom::rnd() < _preco3M_3M) {
// 3Mto3M colour reconnection
auto reconnected = _reconnect3Mto3M(*sel[0], *sel[1], *sel[2],
reco_info);
(*sel[0]) = std::get<0>(reconnected);
(*sel[1]) = std::get<1>(reconnected);
(*sel[2]) = std::get<2>(reconnected);
} else if (kind_of_reco=="2Bto3M" && UseRandom::rnd() < _precoBbarB_3M) {
// antibaryonic and baryonic to 3 mesonic reconnecion
auto reconnected = _reconnectBbarBto3M(*sel[0], *sel[1],
reco_info[0], reco_info[1], reco_info[2]);
(*sel[0]) = std::get<0>(reconnected);
(*sel[1]) = std::get<1>(reconnected);
newcv.push_back(std::get<2>(reconnected));
} else if (kind_of_reco=="3Mto2B" && UseRandom::rnd() < _preco3M_BBbar) {
// 3 mesonic to antibaryonic and baryonic reconnection
ClusterPtr b1, b2;
_makeBaryonicClusters(*sel[0], *sel[1], *sel[2], b1, b2);
(*sel[0]) = b1;
(*sel[1]) = b2;
deleted.push_back(*sel[2]);
} else if (kind_of_reco=="2Bto2B" && UseRandom::rnd() < _preco2B_2B) {
// 2 (anti)baryonic to 2 (anti)baryonic reconnection
auto reconnected = _reconnect2Bto2B(*sel[0], *sel[1],
reco_info[0], reco_info[1]);
(*sel[0]) = reconnected.first;
(*sel[1]) = reconnected.second;
} else if (kind_of_reco=="MBtoMB" && UseRandom::rnd() < _precoMB_MB) {
// (anti)baryonic and mesonic to (anti)baryonic and mesonic reconnection
auto reconnected = _reconnectMBtoMB(*sel[0], *sel[1],
reco_info[0]);
(*sel[0]) = reconnected.first;
(*sel[1]) = reconnected.second;
}
// erase the sel-vector
sel.erase(sel.begin(), sel.end());
}
// write to clustervector new CR'd clusters and deleting
// all deleted clusters
ClusterVector clustervector;
for (const auto & cluster : newcv)
if (find(deleted.begin(), deleted.end(), cluster) == deleted.end())
clustervector.push_back(cluster);
swap(cv, clustervector);
}
void ColourReconnector::_findbestreconnectionoption(std::vector<CluVecIt> & cls, const unsigned & baryonic,
string & kind_of_reco, int (&reco_info)[3]) const {
double min_displacement;
if (baryonic==0) {
// case with 3 mesonic clusters
assert(cls.size()==3);
// calculate the initial displacement sum
min_displacement = _mesonToBaryonFactor * _displacement((*cls[0])->particle(0), (*cls[0])->particle(1));
min_displacement += _mesonToBaryonFactor * _displacement((*cls[1])->particle(0), (*cls[1])->particle(1));
min_displacement += _mesonToBaryonFactor * _displacement((*cls[2])->particle(0), (*cls[2])->particle(1));
// find best CR reco_info and kind_of_reco
_3MtoXreconnectionfinder(cls,
reco_info[0], reco_info[1], reco_info[2], min_displacement, kind_of_reco);
/**
* kind_of_reco either "3Mto3M" or "3Mto2B" (or "" if no better configuration is found)
* case 3Mto3M: the coloured particle of the i-th cluster forms a new cluster with the
* antiparticle of the reco_info[i]-th cluster
* case 3MtoBbarB: all 3 (anti)coloured particle form a new (anti)baryonic cluster
*/
} else if (baryonic == 1) {
// case 1 baryonic and 1 mesonic cluster
assert(cls.size()==2);
// make mesonic cluster always the cls[0]
if ((*cls[0])->numComponents() == 3) {
ClusterPtr zw = *cls[0];
*cls[0] = *cls[1];
*cls[1] = zw;
}
// calculate the initial displacement sum
min_displacement = _mesonToBaryonFactor *_displacement((*cls[0])->particle(0), (*cls[0])->particle(1));
min_displacement += _displacementBaryonic((*cls[1])->particle(0), (*cls[1])->particle(1), (*cls[1])->particle(2));
// find best CR reco_info and kind_of_reco
_BMtoBMreconnectionfinder(*cls[0], *cls[1],
reco_info[0], min_displacement, kind_of_reco);
/**
* reco_info[0] is the index of the (anti)quarks of the baryonic cluster cls[1], which should
* be swapped with the (anti)quarks of the mesonic cluster cls[0]
*/
} else {
assert(baryonic==2);
assert(cls.size()==2);
// calculate the initial displacement sum
min_displacement = _displacementBaryonic((*cls[0])->particle(0), (*cls[0])->particle(1), (*cls[0])->particle(2));
min_displacement += _displacementBaryonic((*cls[1])->particle(0), (*cls[1])->particle(1), (*cls[1])->particle(2));
// case 2 (anti)baryonic clusters to 2 other (anti)baryonic clusters
if ( ( (*cls[0])->particle(0)->hasColour() && (*cls[1])->particle(0)->hasColour() )
|| ( (*cls[0])->particle(0)->hasColour(true) && (*cls[1])->particle(0)->hasColour(true) ) ) {
// find best CR reco_info and kind_of_reco
_2Bto2BreconnectionFinder(*cls[0], *cls[1],
reco_info[0], reco_info[1], min_displacement, kind_of_reco);
/**
* swap the reco_info[0]-th particle of the first cluster in the vector with the
* reco_info[1]-th particle of the second cluster
*/
} else {
// case 1 baryonic and 1 antibaryonic cluster to 3 mesonic clusters
// find best CR reco_info and kind_of_reco
_BbarBto3MreconnectionFinder(*cls[0], *cls[1],
reco_info[0], reco_info[1], reco_info[2], min_displacement, kind_of_reco);
/**
* the i-th particle of the first cluster form a new mesonic cluster with the
* reco_info[i]-th particle of the second cluster
*/
}
}
return;
}
void ColourReconnector::_findPartnerBaryonicDiquarkCluster(
CluVecIt cl, ClusterVector & cv,
unsigned & typeOfReconnection,
const ClusterVector& deleted,
CluVecIt &candidate1,
CluVecIt &candidate2 ) const {
typeOfReconnection=0; // no Reconnection found
using Constants::pi;
using Constants::twopi;
candidate1=cl;
candidate2=cl;
bool candIsOctet1 = false;
bool candIsOctet2 = false;
bool candIsQQ1 = false;
bool candIsQQ2 = false;
bool foundCR = false;
double maxrap1 = 0.0;
double maxrap2 = 0.0;
double minrap1 = 0.0;
double minrap2 = 0.0;
double maxsum1 = 0.0;
double maxsum2 = 0.0;
double NegativeRapidtyThreshold = 0.0;
double PositiveRapidtyThreshold = 0.0;
// boost into RF of cl
Lorentz5Momentum cl1 = (*cl)->momentum();
// TODO Boost not covariant!! ERROR
const Boost boostv(-cl1.boostVector());
// cl1.boost(boostv);
// boost constituents of cl into RF of cl
// Lorentz5Momentum p1col = (*cl)->colParticle()->momentum();
Lorentz5Momentum p1anticol = (*cl)->antiColParticle()->momentum();
// TODO Boost not covariant!! ERROR
// p1col.boost(boostv);
p1anticol.boost(boostv);
for (CluVecIt cit=cv.begin(); cit != cv.end(); ++cit) {
//avoid looping over clusters containing diquarks
if ( (*cit)->numComponents()>=3 ) continue;
//skip the cluster to be deleted later 3->2 cluster
if ( find(deleted.begin(), deleted.end(), *cit) != deleted.end() )
continue;
if ( hasDiquark(*cit) ) continue;
if ( (*cl)->isBeamCluster() && (*cit)->isBeamCluster() )
continue;
if ( cit==cl ) continue;
// veto if Clusters further apart than _maxDistance
if (_localCR && ((**cl).vertex().vect()-(**cit).vertex().vect()).mag() > _maxDistance) continue;
// veto if Clusters have negative spacetime difference
if (_causalCR && ((**cl).vertex()-(**cit).vertex()).m() < ZERO) continue;
bool octetNormalCR =
(_isColour8( (*cl)->colParticle(), (*cit)->antiColParticle() )
||
_isColour8( (*cit)->colParticle(), (*cl)->antiColParticle() ) );
// boost constituents of cit into RF of cl
Lorentz5Momentum p2col = (*cit)->colParticle()->momentum();
Lorentz5Momentum p2anticol = (*cit)->antiColParticle()->momentum();
p2col.boost(boostv);
p2anticol.boost(boostv);
// calculate the rapidity of the other constituents of the clusters
// w.r.t axis of p1anticol.vect.unit
const double rapq = calculateRapidityRF(p1anticol,p2col);
const double rapqbar = calculateRapidityRF(p1anticol,p2anticol);
// std::cout << "\nPRE\n"
// << "\t octet = " << octetNormalCR << "\n"
// << "\t rapq = " << rapq << "\n"
// << "\t rapqbar = " << rapqbar << "\n";
// configuration for Mesonic CR
if ( rapq > 0.0 && rapqbar < 0.0
&& rapq > PositiveRapidtyThreshold
&& rapqbar < NegativeRapidtyThreshold) {
//sum of rapidities of quarks
const double sumQQbar = abs(rapq) + abs(rapqbar);
if ( sumQQbar > maxsum2 ) {
if ( sumQQbar > maxsum1 ) {
double factor = candIsQQ1 ? _mesonToBaryonFactor:1.0;
maxsum2 = (factor*maxsum1) > sumQQbar ? sumQQbar:(factor*maxsum1);
candidate2 = candidate1;
candIsQQ2 = candIsQQ1;
candIsOctet2 = candIsOctet1;
maxrap2 = maxrap1;
minrap2 = minrap1;
maxsum1 = sumQQbar;
candidate1 = cit;
candIsQQ1 = false;
candIsOctet1 = octetNormalCR;
maxrap1 = rapq;
minrap1 = rapqbar;
} else {
maxsum2 = sumQQbar;
candidate2 = cit;
candIsQQ2 = false;
candIsOctet2 = octetNormalCR;
maxrap2 = rapq;
minrap2 = rapqbar;
}
// choose the less stringent threshold for further iterations
PositiveRapidtyThreshold = maxrap1 > maxrap2 ? maxrap2:maxrap1;
NegativeRapidtyThreshold = minrap1 < minrap2 ? minrap2:minrap1;
foundCR=true;
}
}
assert(PositiveRapidtyThreshold<=maxrap1);
assert(PositiveRapidtyThreshold<=maxrap2);
assert(NegativeRapidtyThreshold>=minrap1);
assert(NegativeRapidtyThreshold>=minrap2);
assert(maxsum1>=maxsum2);
if ( rapq < 0.0 && rapqbar > 0.0
&& rapqbar > PositiveRapidtyThreshold/_mesonToBaryonFactor
&& rapq < NegativeRapidtyThreshold/_mesonToBaryonFactor ) {
//sum of rapidities of quarks
const double sumQQ = abs(rapq) + abs(rapqbar);
if ( sumQQ > maxsum2/_mesonToBaryonFactor ) {
if ( sumQQ > maxsum1 ) {
double factor = candIsQQ1 ? _mesonToBaryonFactor:1.0;
maxsum2 = (factor*maxsum1) > sumQQ ? sumQQ:(factor*maxsum1);
candidate2 = candidate1;
candIsQQ2 = candIsQQ1;
candIsOctet2 = candIsOctet1;
maxrap2 = maxrap1;
minrap2 = minrap1;
maxsum1 = sumQQ;
candidate1 = cit;
candIsQQ1 = true;
candIsOctet1 = octetNormalCR;
maxrap1 = rapqbar;
minrap1 = rapq;
} else {
maxsum2 = (_mesonToBaryonFactor*maxsum1) > sumQQ ? sumQQ:(_mesonToBaryonFactor*maxsum1);
candidate2 = cit;
candIsQQ2 = true;
candIsOctet2 = octetNormalCR;
maxrap2 = rapqbar;
minrap2 = rapq;
}
// choose the less stringent threshold for further iterations
PositiveRapidtyThreshold = maxrap1 > maxrap2 ? maxrap2:maxrap1;
NegativeRapidtyThreshold = minrap1 < minrap2 ? minrap2:minrap1;
foundCR=true;
}
}
assert(PositiveRapidtyThreshold<=maxrap1);
assert(PositiveRapidtyThreshold<=maxrap2);
assert(NegativeRapidtyThreshold>=minrap1);
assert(NegativeRapidtyThreshold>=minrap2);
assert(maxsum1>=maxsum2);
}
// determine the type
if (!candIsQQ1) {
if (candIsOctet1) {
if (!candIsQQ2 && !candIsOctet2) {
swap(candidate2,candidate1);
typeOfReconnection = 1;
}
else
typeOfReconnection = 0;
}
else
typeOfReconnection = 1;
}
else if (candIsQQ1)
{
if (candIsQQ2)
typeOfReconnection = 2;
else
typeOfReconnection = 3;
}
if (!foundCR) typeOfReconnection = 0;
// veto reconnection if cannot make a Diquark Cluster
bool failDCR=false;
if (typeOfReconnection == 3) {
if (_diquarkCR && !_canMakeDiquarkCluster(*cl, *candidate1)) {
if (!candIsQQ2 && !candIsOctet2 && candidate2!=cl) {
// if second nearest is candidate for Mesonic CR
// allow MCR
candidate1=candidate2;
typeOfReconnection=1;
}
else if ( _canMakeDiquarkCluster(*cl, *candidate2) && candidate2!=cl) {
// if second nearest is allowed for DCR
// allow DCR
candidate1=candidate2;
typeOfReconnection=3;
}
else
{
// No CR
typeOfReconnection = 0;
failDCR=true;
}
}
}
if (_debug) {
std::ofstream outTypes("WriteOut/TypesOfDCR.dat", std::ios::app);
outTypes << (failDCR ? 4:typeOfReconnection) << "\n";
outTypes.close();
switch (typeOfReconnection)
{
// Mesonic CR
case 1:
{
std::ofstream outMCR("WriteOut/MCR.dat", std::ios::app);
outMCR << minrap1 << "\t"
<< maxrap1 << "\t"
<< minrap2 << "\t"
<< maxrap2 << "\t"
<<"\n";
outMCR.close();
break;
}
// Baryonic CR
case 2:
{
std::ofstream outBCR("WriteOut/BCR.dat", std::ios::app);
outBCR << minrap1 << "\t"
<< maxrap1 << "\t"
<< minrap2 << "\t"
<< maxrap2 << "\t"
<<"\n";
outBCR.close();
break;
}
// Diquark CR
case 3:
{
std::ofstream outDCR("WriteOut/DCR.dat", std::ios::app);
outDCR << minrap1 << "\t"
<< maxrap1 << "\t"
<< minrap2 << "\t"
<< maxrap2 << "\t"
<<"\n";
outDCR.close();
break;
}
// No CR found
case 0:
{
std::ofstream outNoCR("WriteOut/NoCR.dat", std::ios::app);
outNoCR<< minrap1 << "\t"
<< maxrap1 << "\t"
<< minrap2 << "\t"
<< maxrap2 << "\t"
<<"\n";
outNoCR.close();
break;
}
default:
assert(false);
}
}
}
// namespace {
struct ClusterInfo{
double rapSumMin;
bool isQQ;
bool isOctetWithOriginal;
CluVecIt cluIt;
} ;
bool compare_Clusters(ClusterInfo c1, ClusterInfo c2){
return c1.rapSumMin>c2.rapSumMin;
}
// }
void ColourReconnector::_findPartnerBaryonicDiquarkClusterTEST(
CluVecIt cl, ClusterVector & cv,
unsigned & typeOfReconnection,
const ClusterVector& deleted,
CluVecIt &candidate1,
CluVecIt &candidate2 ) const {
typeOfReconnection=0; // no Reconnection found
using Constants::pi;
using Constants::twopi;
candidate1=cl;
candidate2=cl;
// boost into RF of cl
Lorentz5Momentum cl1 = (*cl)->momentum();
const Boost boostv(-cl1.boostVector());
Lorentz5Momentum p1anticol = (*cl)->antiColParticle()->momentum();
// direction
p1anticol.boost(boostv);
std::vector<ClusterInfo> cluVec;
for (CluVecIt cit=cv.begin(); cit != cv.end(); ++cit) {
//avoid looping over clusters containing diquarks
if ( (*cit)->numComponents()>=3 ) continue;
//skip the cluster to be deleted later 3->2 cluster
if ( find(deleted.begin(), deleted.end(), *cit) != deleted.end() )
continue;
if ( hasDiquark(*cit) ) continue;
if ( (*cl)->isBeamCluster() && (*cit)->isBeamCluster() )
continue;
if ( cit==cl ) continue;
// veto if Clusters further apart than _maxDistance
if (_localCR && ((**cl).vertex().vect()-(**cit).vertex().vect()).mag() > _maxDistance) continue;
// veto if Clusters have negative spacetime difference
if (_causalCR && ((**cl).vertex()-(**cit).vertex()).m() < ZERO) continue;
bool octetNormalCR =
(_isColour8( (*cl)->colParticle(), (*cit)->antiColParticle() )
||
_isColour8( (*cit)->colParticle(), (*cl)->antiColParticle() ) );
// boost constituents of cit into RF of cl
Lorentz5Momentum p2col = (*cit)->colParticle()->momentum();
Lorentz5Momentum p2anticol = (*cit)->antiColParticle()->momentum();
p2col.boost(boostv);
p2anticol.boost(boostv);
// calculate the rapidity of the other constituents of the clusters
// w.r.t axis of p1anticol.vect.unit
const double rapq = calculateRapidityRF(p1anticol,p2col);
const double rapqbar = calculateRapidityRF(p1anticol,p2anticol);
// std::cout << "\nPRE\n"
// << "\t octet = " << octetNormalCR << "\n"
// << "\t rapq = " << rapq << "\n"
// << "\t rapqbar = " << rapqbar << "\n";
const double sumAbsRap = abs(rapq) + abs(rapqbar);
// configuration for Mesonic CR
if ( rapq > 0.0 && rapqbar < 0.0 ) {
//sum of rapidities of quarks
ClusterInfo cluMesInfo;
cluMesInfo.rapSumMin=sumAbsRap;
cluMesInfo.isQQ=false;
cluMesInfo.isOctetWithOriginal=octetNormalCR;
cluMesInfo.cluIt=cit;
cluVec.push_back(cluMesInfo);
}
else if ( rapq < 0.0 && rapqbar > 0.0 ) {
//sum of rapidities of quarks
ClusterInfo cluMesInfo;
cluMesInfo.rapSumMin=sumAbsRap;
cluMesInfo.isQQ=false;
cluMesInfo.isOctetWithOriginal=octetNormalCR;
cluMesInfo.cluIt=cit;
cluVec.push_back(cluMesInfo);
}
}
std::sort(cluVec.begin(), cluVec.end(), compare_Clusters);
switch (cluVec.size())
{
case 0:
typeOfReconnection=0;
return;
case 1:
{
if (cluVec[0].isQQ) {
// diquark CR if possible
if (_canMakeDiquarkCluster(*cl, *(cluVec[0].cluIt))) {
typeOfReconnection=3;
candidate1=cluVec[0].cluIt;
return;
}
else {
typeOfReconnection=0;
return;
}
}
else {
// Mesonic CR if not Octet
if (cluVec[0].isOctetWithOriginal){
typeOfReconnection=0;
return;
}
else {
candidate1=cluVec[0].cluIt;
typeOfReconnection=1;
return;
}
}
break;
}
}
bool candidate1isQQ=false;
bool candidate2isQQ=false;
typeOfReconnection=0;
for (auto c : cluVec) {
if (candidate1!=cl && candidate2!=cl)
break;
if (c.isQQ) {
// diquark CR if possible
if (_canMakeDiquarkCluster(*cl, *(c.cluIt))) {
if (candidate1==cl){
candidate1=c.cluIt;
candidate1isQQ=true;
typeOfReconnection=3;
}
else {
candidate2=c.cluIt;
candidate2isQQ=true;
typeOfReconnection=2;
}
}
}
else {
// Mesonic CR if not Octet
if (c.isOctetWithOriginal) {
continue;
}
else {
if (candidate1==cl){
candidate1=c.cluIt;
candidate1isQQ=false;
typeOfReconnection=1;
}
else {
candidate2=c.cluIt;
candidate2isQQ=false;
if (candidate1isQQ)
typeOfReconnection=3;
else
typeOfReconnection=1;
break;
}
}
}
}
return;
}
CluVecIt ColourReconnector::_findPartnerBaryonic(
CluVecIt cl, ClusterVector & cv,
bool & baryonicCand,
const ClusterVector& deleted,
CluVecIt &baryonic1,
CluVecIt &baryonic2 ) const {
using Constants::pi;
using Constants::twopi;
// Returns a candidate for possible reconnection
CluVecIt candidate = cl;
bool bcand = false;
double maxrap = 0.0;
double minrap = 0.0;
double maxrapNormal = 0.0;
double minrapNormal = 0.0;
double maxsumnormal = 0.0;
double maxsum = 0.0;
double secondsum = 0.0;
// boost into RF of cl
Lorentz5Momentum cl1 = (*cl)->momentum();
const Boost boostv(-cl1.boostVector());
// cl1.boost(boostv);
// boost constituents of cl into RF of cl
// Lorentz5Momentum p1col = (*cl)->colParticle()->momentum();
Lorentz5Momentum p1anticol = (*cl)->antiColParticle()->momentum();
// p1col.boost(boostv);
p1anticol.boost(boostv);
for (CluVecIt cit=cv.begin(); cit != cv.end(); ++cit) {
//avoid looping over clusters containing diquarks
if ( hasDiquark(*cit) ) continue;
if ( (*cit)->numComponents()>=3 ) continue;
if ( cit==cl ) continue;
//skip the cluster to be deleted later 3->2 cluster
if ( find(deleted.begin(), deleted.end(), *cit) != deleted.end() )
continue;
if ( (*cl)->isBeamCluster() && (*cit)->isBeamCluster() )
continue;
// veto if Clusters further apart than _maxDistance
if (_localCR && ((**cl).vertex().vect()-(**cit).vertex().vect()).mag() > _maxDistance) continue;
// veto if Clusters have negative spacetime difference
if (_causalCR && ((**cl).vertex()-(**cit).vertex()).m() < ZERO) continue;
const bool Colour8 =
_isColour8( (*cl)->colParticle(), (*cit)->antiColParticle() )
||
_isColour8( (*cit)->colParticle(), (*cl)->antiColParticle() ) ;
// boost constituents of cit into RF of cl
Lorentz5Momentum p2col = (*cit)->colParticle()->momentum();
Lorentz5Momentum p2anticol = (*cit)->antiColParticle()->momentum();
p2col.boost(boostv);
p2anticol.boost(boostv);
// calculate the rapidity of the other constituents of the clusters
// w.r.t axis of p1anticol.vect.unit
const double rapq = calculateRapidityRF(p1anticol,p2col);
const double rapqbar = calculateRapidityRF(p1anticol,p2anticol);
// configuration for normal CR
if ( !Colour8
&& rapq > 0.0 && rapqbar < 0.0
&& rapq > maxrap
&& rapqbar < minrap ) {
maxrap = rapq;
minrap = rapqbar;
//sum of rapidities of quarks
const double normalsum = abs(rapq) + abs(rapqbar);
if ( normalsum > maxsumnormal ) {
maxsumnormal = normalsum;
maxrapNormal = rapq;
minrapNormal = rapqbar;
bcand = false;
candidate = cit;
}
}
if ( rapq < 0.0 && rapqbar >0.0
&& rapqbar > maxrapNormal
&& rapq < minrapNormal ) {
maxrap = rapqbar;
minrap = rapq;
const double sumrap = abs(rapqbar) + abs(rapq);
// first candidate gets here. If second baryonic candidate has higher Ysum than the first
// one, the second candidate becomes the first one and the first the second.
if (sumrap > maxsum) {
if (maxsum != 0) {
baryonic2 = baryonic1;
baryonic1 = cit;
bcand = true;
} else {
baryonic1 = cit;
}
maxsum = sumrap;
} else {
if (sumrap > secondsum && sumrap != maxsum) {
secondsum = sumrap;
bcand = true;
baryonic2 = cit;
}
}
}
}
if (bcand == true) {
baryonicCand = true;
}
if (_debug) {
std::ofstream outTypes("WriteOut/TypesOfBCR.dat", std::ios::app);
outTypes << (baryonicCand ? 2:(candidate==cl ? 0:1)) << "\n";
outTypes.close();
}
return candidate;
}
CluVecIt ColourReconnector::_findRecoPartnerPlainDynamic(CluVecIt cl,
- ClusterVector & cv, const ClusterVector & deleted) const {
+ ClusterVector & cv, const ClusterVector & deleted, bool diquarkCR) const {
CluVecIt candidate = cl;
double pNoRecoMin=1.0;
double pNoReco,preco,precoDiquark;
for (CluVecIt cit=cv.begin(); cit != cv.end(); ++cit) {
// Skip deleted clusters
if (find(deleted.begin(), deleted.end(), *cit) != deleted.end())
continue;
// skip diquark clusters
if ((*cit)->numComponents()>2 || (_diquarkCR && hasDiquark(*cit))) continue;
// don't even look at original cluster
if (cit==cl) continue;
// don't allow colour octet clusters
if ( _isColour8( (*cl)->colParticle(),
(*cit)->antiColParticle() ) ||
_isColour8( (*cit)->colParticle(),
(*cl)->antiColParticle() ) ) {
continue;
}
// stop it putting beam remnants together
if ((*cl)->isBeamCluster() && (*cit)->isBeamCluster()) continue;
// veto if Clusters further apart than _maxDistance
if (_localCR && ((**cl).vertex().vect()-(**cit).vertex().vect()).mag() > _maxDistance) continue;
// veto if Clusters have negative spacetime difference
if (_causalCR && ((**cl).vertex()-(**cit).vertex()).m() < ZERO) continue;
// Get dynamic CR probabilities and try to minimize the non-reco probability
- std::tie( pNoReco, preco, precoDiquark) = _dynamicRecoProbabilitiesCF2(*cl,*cit,false);
+ std::tie( pNoReco, preco, precoDiquark) = _dynamicRecoProbabilitiesCF2(*cl,*cit,diquarkCR);
+ double PS;
+ if (!_canMakeDiquarkCluster(*cl,*cit, PS))
+ PS=0;
+ if (precoDiquark){ //Reweight
+ double newSum=pNoReco+preco+PS*precoDiquark;
+ pNoReco=pNoReco/newSum;
+ preco=preco/newSum;
+ precoDiquark=PS*precoDiquark/newSum;
+ assert(fabs(pNoReco+preco+precoDiquark-1)<1e-6);
+ }
if (pNoReco<pNoRecoMin) {
candidate = cit;
pNoRecoMin=pNoReco;
}
}
if (_debug) {
ofstream out("WriteOut/pNoReco.dat", std::ios::app);
out << pNoRecoMin << "\t" << cv.size() << "\n";
out.close();
}
return candidate;
}
CluVecIt ColourReconnector::_findRecoPartnerPlain(CluVecIt cl,
ClusterVector & cv, const ClusterVector & deleted) const {
CluVecIt candidate = cl;
Energy minMass = 1*TeV;
for (CluVecIt cit=cv.begin(); cit != cv.end(); ++cit) {
// Skip deleted clusters
if (find(deleted.begin(), deleted.end(), *cit) != deleted.end())
continue;
// skip diquark clusters
if ((*cit)->numComponents()>2 || (_diquarkCR && hasDiquark(*cit))) continue;
// don't even look at original cluster
if (cit==cl) continue;
// don't allow colour octet clusters
if ( _isColour8( (*cl)->colParticle(),
(*cit)->antiColParticle() ) ||
_isColour8( (*cit)->colParticle(),
(*cl)->antiColParticle() ) ) {
continue;
}
// stop it putting beam remnants together
if ((*cl)->isBeamCluster() && (*cit)->isBeamCluster()) continue;
// veto if Clusters further apart than _maxDistance
if (_localCR && ((**cl).vertex().vect()-(**cit).vertex().vect()).mag() > _maxDistance) continue;
// veto if Clusters have negative spacetime difference
if (_causalCR && ((**cl).vertex()-(**cit).vertex()).m() < ZERO) continue;
// momenta of the old clusters
Lorentz5Momentum p1 = (*cl)->colParticle()->momentum() +
(*cl)->antiColParticle()->momentum();
Lorentz5Momentum p2 = (*cit)->colParticle()->momentum() +
(*cit)->antiColParticle()->momentum();
// momenta of the new clusters
Lorentz5Momentum p3 = (*cl)->colParticle()->momentum() +
(*cit)->antiColParticle()->momentum();
Lorentz5Momentum p4 = (*cit)->colParticle()->momentum() +
(*cl)->antiColParticle()->momentum();
Energy oldMass = abs( p1.m() ) + abs( p2.m() );
Energy newMass = abs( p3.m() ) + abs( p4.m() );
if ( newMass < oldMass && newMass < minMass ) {
minMass = newMass;
candidate = cit;
}
}
return candidate;
}
// forms two baryonic clusters from three clusters
void ColourReconnector::_makeBaryonicClusters(
ClusterPtr &c1, ClusterPtr &c2,
ClusterPtr &c3,
ClusterPtr &newcluster1,
ClusterPtr &newcluster2) const {
//make sure they all have 2 components
assert(c1->numComponents()==2);
assert(c2->numComponents()==2);
assert(c3->numComponents()==2);
//abandon children
c1->colParticle()->abandonChild(c1);
c1->antiColParticle()->abandonChild(c1);
c2->colParticle()->abandonChild(c2);
c2->antiColParticle()->abandonChild(c2);
c3->colParticle()->abandonChild(c3);
c3->antiColParticle()->abandonChild(c3);
newcluster1 = new_ptr(Cluster(c1->colParticle(),c2->colParticle(), c3->colParticle()));
c1->colParticle()->addChild(newcluster1);
c2->colParticle()->addChild(newcluster1);
c3->colParticle()->addChild(newcluster1);
newcluster1->setVertex(LorentzPoint());
newcluster2 = new_ptr(Cluster(c1->antiColParticle(), c2->antiColParticle(),
c3->antiColParticle()));
c1->antiColParticle()->addChild(newcluster2);
c2->antiColParticle()->addChild(newcluster2);
c3->antiColParticle()->addChild(newcluster2);
newcluster2->setVertex(LorentzPoint());
}
bool ColourReconnector::_canMakeDiquarkCluster(tcPPtr pCol1, tcPPtr pCol2,tcPPtr pAntiCol1, tcPPtr pAntiCol2) const{
double a;
return _canMakeDiquarkCluster( pCol1, pCol2, pAntiCol1, pAntiCol2,a);
}
bool ColourReconnector::_canMakeDiquarkCluster(tcPPtr pCol1, tcPPtr pCol2,tcPPtr pAntiCol1, tcPPtr pAntiCol2, double & PhaseSpace) const{
tcPDPtr dataDiquark = _hadronSpectrum->makeDiquark(pCol1->dataPtr(), pCol2->dataPtr());
tcPDPtr dataDiquarkBar = _hadronSpectrum->makeDiquark(pAntiCol1->dataPtr(), pAntiCol2->dataPtr());
if (!dataDiquark){
throw Exception() << "Could not make a diquark from"
<< pCol1->dataPtr()->PDGName() << " and "
<< pCol2->dataPtr()->PDGName()
<< " in ColourReconnector::_canMakeDiquarkCluster()"
<< Exception::eventerror;
}
if (!dataDiquarkBar){
throw Exception() << "Could not make an anti-diquark from"
<< pAntiCol1->dataPtr()->PDGName() << " and "
<< pAntiCol2->dataPtr()->PDGName()
<< " in ColourReconnector::_canMakeDiquarkCluster()"
<< Exception::eventerror;
}
Energy minMass=_hadronSpectrum->massLightestHadronPair(dataDiquark,dataDiquarkBar);
Lorentz5Momentum Ptot=pCol1->momentum() + pCol2->momentum() + pAntiCol1->momentum() + pAntiCol2->momentum();
Energy Mass=Ptot.m();
if ( Mass<=minMass ) {
return false;
}
if (_phaseSpaceDiquarkFission){
// Tried to add a factor suppressing to low mass Diquark Clusters
double factor;
switch (_phaseSpaceDiquarkFission)
{
case 1:
{
auto BaryonPair=_hadronSpectrum->lightestHadronPair(dataDiquark,dataDiquarkBar);
factor = 2.0*Kinematics::pstarTwoBodyDecay(Mass,BaryonPair.first->mass(),BaryonPair.second->mass())/Mass;
break;
}
case 2:
{
factor = 2.0*Kinematics::pstarTwoBodyDecay(Mass,dataDiquark->constituentMass(),dataDiquarkBar->constituentMass())/Mass;
break;
}
default:
assert(false);
}
assert(factor>=0.0);
assert(factor<=1.0);
PhaseSpace = factor;
}
else {
PhaseSpace = 1.0;
}
return true;
}
bool ColourReconnector::_canMakeDiquarkCluster(const ClusterPtr &c1, const ClusterPtr &c2) const {
double a;
return _canMakeDiquarkCluster(c1,c2,a);
}
// forms a four-quark cluster
bool ColourReconnector::_canMakeDiquarkCluster(const ClusterPtr &c1, const ClusterPtr &c2, double & PhaseSpace) const {
//make sure they all have 2 components
assert(c1->numComponents()==2);
assert(c2->numComponents()==2);
return _canMakeDiquarkCluster(c1->colParticle(),c2->colParticle(),c1->antiColParticle(),c2->antiColParticle(),PhaseSpace);
}
bool ColourReconnector::_makeDiquarkCluster(
ClusterPtr &c1, ClusterPtr &c2,
ClusterPtr &newcluster) const{
// std::cout << "MakeDiq" << std::endl;
if (!_canMakeDiquarkCluster(c1,c2)){
std::cout << "Should never execute! check this earlier" << std::endl;
return false;
}
//abandon children
c1->colParticle()->abandonChild(c1);
c1->antiColParticle()->abandonChild(c1);
c2->colParticle()->abandonChild(c2);
c2->antiColParticle()->abandonChild(c2);
newcluster = new_ptr(Cluster(c1->colParticle(),c2->colParticle(),
c1->antiColParticle(), c2->antiColParticle()));
c1->colParticle()->addChild(newcluster);
c2->colParticle()->addChild(newcluster);
c1->antiColParticle()->addChild(newcluster);
c2->antiColParticle()->addChild(newcluster);
newcluster->setVertex(LorentzPoint());
return true;
}
pair <ClusterPtr,ClusterPtr>
ColourReconnector::_reconnect2Bto2B(ClusterPtr &c1, ClusterPtr &c2, const int s1, const int s2) const {
// form the first new cluster
// separate the quarks from their original cluster
c1->particleB((s1+1)%3)->abandonChild(c1);
c1->particleB((s1+2)%3)->abandonChild(c1);
c2->particleB(s2)->abandonChild(c2);
// now the new cluster
ClusterPtr newCluster1 = new_ptr(Cluster(c1->particleB((s1+1)%3), c1->particleB((s1+2)%3), c2->particleB(s2)));
c1->particleB((s1+1)%3)->addChild(newCluster1);
c1->particleB((s1+2)%3)->addChild(newCluster1);
c2->particleB(s2)->addChild(newCluster1);
// set new vertex
newCluster1->setVertex(LorentzPoint());
// set beam remnants for new cluster
if (c1->isBeamRemnant((s1+1)%3)) newCluster1->setBeamRemnant(0, true);
if (c1->isBeamRemnant((s1+2)%3)) newCluster1->setBeamRemnant(1, true);
if (c2->isBeamRemnant(s2)) newCluster1->setBeamRemnant(2, true);
// for the second cluster same procedure
c2->particleB((s2+1)%3)->abandonChild(c2);
c2->particleB((s2+2)%3)->abandonChild(c2);
c1->particleB(s1)->abandonChild(c1);
ClusterPtr newCluster2 = new_ptr(Cluster(c2->particleB((s2+1)%3), c2->particleB((s2+2)%3), c1->particleB(s1)));
c2->particleB((s2+1)%3)->addChild(newCluster2);
c2->particleB((s2+2)%3)->addChild(newCluster2);
c1->particleB(s1)->addChild(newCluster2);
newCluster2->setVertex(LorentzPoint());
if (c2->isBeamRemnant((s2+1)%3)) newCluster2->setBeamRemnant(0, true);
if (c2->isBeamRemnant((s2+2)%3)) newCluster2->setBeamRemnant(1, true);
if (c1->isBeamRemnant(s1)) newCluster2->setBeamRemnant(2, true);
return pair <ClusterPtr, ClusterPtr> (newCluster1, newCluster2);
}
std::tuple <ClusterPtr, ClusterPtr, ClusterPtr>
ColourReconnector::_reconnectBbarBto3M(ClusterPtr & c1, ClusterPtr & c2, const int s0, const int s1, const int s2) const {
// make sure they all have 3 components
assert(c1->numComponents()==3);
assert(c2->numComponents()==3);
// first Cluster
c1->particleB(0)->abandonChild(c1);
c2->particleB(s0)->abandonChild(c2);
ClusterPtr newCluster1 = new_ptr(Cluster(c1->particleB(0), c2->particleB(s0)));
c1->particleB(0)->addChild(newCluster1);
c2->particleB(s0)->addChild(newCluster1);
// set new vertex
newCluster1->setVertex(0.5*(c1->particleB(0)->vertex() + c2->particleB(s0)->vertex()));
// set beam remnants for new cluster
if (c1->isBeamRemnant(0)) newCluster1->setBeamRemnant(0, true);
if (c2->isBeamRemnant(s0)) newCluster1->setBeamRemnant(1, true);
// same for second cluster
c1->particleB(1)->abandonChild(c1);
c2->particleB(s1)->abandonChild(c2);
ClusterPtr newCluster2 = new_ptr(Cluster(c1->particleB(1), c2->particleB(s1)));
c1->particleB(1)->addChild(newCluster2);
c2->particleB(s1)->addChild(newCluster2);
newCluster2->setVertex(0.5*(c1->particleB(1)->vertex() + c2->particleB(s1)->vertex()));
if (c1->isBeamRemnant(1)) newCluster2->setBeamRemnant(0, true);
if (c2->isBeamRemnant(s1)) newCluster2->setBeamRemnant(1, true);
// same for third cluster
c1->particleB(2)->abandonChild(c1);
c2->particleB(s2)->abandonChild(c2);
ClusterPtr newCluster3 = new_ptr(Cluster(c1->particleB(2), c2->particleB(s2)));
c1->particleB(2)->addChild(newCluster3);
c2->particleB(s2)->addChild(newCluster3);
newCluster3->setVertex(0.5*(c1->particleB(2)->vertex() + c2->particleB(s2)->vertex()));
if (c1->isBeamRemnant(2)) newCluster3->setBeamRemnant(0, true);
if (c2->isBeamRemnant(s2)) newCluster3->setBeamRemnant(1, true);
return std::tuple <ClusterPtr, ClusterPtr, ClusterPtr> (newCluster1, newCluster2, newCluster3);
}
pair <ClusterPtr,ClusterPtr>
ColourReconnector::_reconnect(ClusterPtr &c1, ClusterPtr &c2) const {
// choose the other possibility to form two clusters from the given
// constituents
assert(c1->numComponents()==2);
assert(c2->numComponents()==2);
int c1_col(-1),c1_anti(-1),c2_col(-1),c2_anti(-1);
for(unsigned int ix=0; ix<2; ++ix) {
if (c1->particle(ix)->hasColour(false)) c1_col = ix;
else if(c1->particle(ix)->hasColour(true )) c1_anti = ix;
if (c2->particle(ix)->hasColour(false)) c2_col = ix;
else if(c2->particle(ix)->hasColour(true )) c2_anti = ix;
}
assert(c1_col>=0&&c2_col>=0&&c1_anti>=0&&c2_anti>=0);
c1->colParticle()->abandonChild(c1);
c2->antiColParticle()->abandonChild(c2);
ClusterPtr newCluster1
= new_ptr( Cluster( c1->colParticle(), c2->antiColParticle() ) );
c1->colParticle()->addChild(newCluster1);
c2->antiColParticle()->addChild(newCluster1);
/*
* TODO: Questionable setting of the vertex
* */
newCluster1->setVertex(0.5*(c1->colParticle()->vertex() +
c2->antiColParticle()->vertex()));
if(c1->isBeamRemnant(c1_col )) newCluster1->setBeamRemnant(0,true);
if(c2->isBeamRemnant(c2_anti)) newCluster1->setBeamRemnant(1,true);
c1->antiColParticle()->abandonChild(c1);
c2->colParticle()->abandonChild(c2);
ClusterPtr newCluster2
= new_ptr( Cluster( c2->colParticle(), c1->antiColParticle() ) );
c1->antiColParticle()->addChild(newCluster2);
c2->colParticle()->addChild(newCluster2);
/*
* TODO: Questionable setting of the vertex
* */
newCluster2->setVertex(0.5*(c2->colParticle()->vertex() +
c1->antiColParticle()->vertex()));
if(c2->isBeamRemnant(c2_col )) newCluster2->setBeamRemnant(0,true);
if(c1->isBeamRemnant(c1_anti)) newCluster2->setBeamRemnant(1,true);
return pair <ClusterPtr,ClusterPtr> (newCluster1, newCluster2);
}
std::tuple <ClusterPtr, ClusterPtr> ColourReconnector::_reconnect3MtoMD(ClusterVector & cluvec, const int topology) const {
// std::cout << "MakeDiq3MtoMD" << std::endl;
assert(cluvec.size()==3);
assert(topology>0 && topology<=33);
int colIndexMCR = (topology/10)-1;
int antiColIndexMCR = (topology-(colIndexMCR+1)*10)-1;
assert(colIndexMCR<3 && colIndexMCR>=0);
assert(antiColIndexMCR<3 && antiColIndexMCR>=0);
// std::cout << "topo = "<< topology <<"\t==\t" << colIndexMCR+1<<antiColIndexMCR+1 << std::endl;
/*
if (colIndexMCR==antiColIndexMCR) {
ClusterPtr DiquarkClu;
assert((colIndexMCR+1)%3!=(colIndexMCR+2)%3);
assert((colIndexMCR+1)%3!=colIndexMCR);
assert((colIndexMCR+2)%3!=colIndexMCR);
if (!_makeDiquarkCluster(cluvec[(colIndexMCR+1)%3], cluvec[(colIndexMCR+2)%3], DiquarkClu)){
std::cout << "could not make diquark cluster of index " << (colIndexMCR+1)%3 <<", " <<(colIndexMCR+2)%3<<"\nColIdx "<<colIndexMCR << std::endl;
assert(false);
// return std::tuple <ClusterPtr, ClusterPtr, ClusterPtr> (cluvec[0],cluvec[1],cluvec[2]);
}
// return (original untouched mesonic cluster, Diquark cluster, to be deleted cluster)
std::cout << "simple DCR" << std::endl;
// TODO check that the ordering is still maintained
return std::tuple <ClusterPtr, ClusterPtr> (cluvec[colIndexMCR], DiquarkClu);
}
*/
ClusterPtr cMCR1 = cluvec[colIndexMCR];
ClusterPtr cMCR2 = cluvec[antiColIndexMCR];
// std::cout << "Never " << std::endl;
// ClusterVector cluvec = cluvec;
ClusterVector newclusters;
int colIdx[3]={-1,-1,-1};
int anticolIdx[3]={-1,-1,-1};
for (int i = 0; i < 3; i++) {
assert(cluvec[i]->numComponents()==2);
for (unsigned ip= 0; ip < 2; ip++){
if (cluvec[i]->particle(ip)->hasColour(false)) colIdx[i] = ip;
if (cluvec[i]->particle(ip)->hasColour(true)) anticolIdx[i] = ip;
}
assert(colIdx[i]>=0);
assert(anticolIdx[i]>=0);
// abbandon all children
cluvec[i]->colParticle()->abandonChild(cluvec[i]);
cluvec[i]->antiColParticle()->abandonChild(cluvec[i]);
}
// make the MCR cluster with indices colIndexMCR,antiColIndexMCR
// form new mesonic cluster
ClusterPtr newClusterMCR = new_ptr(Cluster(cluvec[colIndexMCR]->colParticle(), cluvec[antiColIndexMCR]->antiColParticle()));
newClusterMCR->colParticle()->addChild(newClusterMCR);
newClusterMCR->antiColParticle()->addChild(newClusterMCR);
// set new vertex
newClusterMCR->setVertex(0.5*(newClusterMCR->colParticle()->vertex() +
newClusterMCR->antiColParticle()->vertex()));
// set beam remnants for new cluster
if (cluvec[colIndexMCR]->isBeamRemnant(colIdx[colIndexMCR])) newClusterMCR->setBeamRemnant(0, true);
if (cluvec[antiColIndexMCR]->isBeamRemnant(anticolIdx[antiColIndexMCR])) newClusterMCR->setBeamRemnant(1, true);
// make the DCR cluster with remaining indices
int indexDCRcol1 = (colIndexMCR+1)%3;
int indexDCRcol2 = (colIndexMCR+2)%3;
int indexDCRacol1 = (antiColIndexMCR+1)%3;
int indexDCRacol2 = (antiColIndexMCR+2)%3;
assert(indexDCRcol1!=indexDCRcol2);
assert(indexDCRcol1!=colIndexMCR);
assert(indexDCRcol2!=colIndexMCR);
assert(indexDCRacol1!=indexDCRacol2);
assert(indexDCRacol1!=antiColIndexMCR);
assert(indexDCRacol2!=antiColIndexMCR);
ClusterPtr newClusterDCR = new_ptr(Cluster(
cluvec[indexDCRcol1]->colParticle(),
cluvec[indexDCRcol2]->colParticle(),
cluvec[indexDCRacol1]->antiColParticle(),
cluvec[indexDCRacol2]->antiColParticle()
));
cluvec[indexDCRcol1]->colParticle()->addChild(newClusterDCR);
cluvec[indexDCRcol2]->colParticle()->addChild(newClusterDCR);
cluvec[indexDCRacol1]->antiColParticle()->addChild(newClusterDCR);
cluvec[indexDCRacol2]->antiColParticle()->addChild(newClusterDCR);
// set new vertex
newClusterDCR->setVertex(0.25*(
cluvec[indexDCRcol1]->colParticle()->vertex() +
cluvec[indexDCRcol2]->colParticle()->vertex() +
cluvec[indexDCRacol1]->antiColParticle()->vertex() +
cluvec[indexDCRacol2]->antiColParticle()->vertex()
));
// set beam remnants for new cluster
if (cluvec[indexDCRcol1]->isBeamRemnant(colIdx[indexDCRcol1])) newClusterDCR->setBeamRemnant(0, true);
if (cluvec[indexDCRcol2]->isBeamRemnant(colIdx[indexDCRcol2])) newClusterDCR->setBeamRemnant(1, true);
if (cluvec[indexDCRacol1]->isBeamRemnant(anticolIdx[indexDCRacol1])) newClusterDCR->setBeamRemnant(2, true);
if (cluvec[indexDCRacol2]->isBeamRemnant(anticolIdx[indexDCRacol2])) newClusterDCR->setBeamRemnant(3, true);
return std::tuple <ClusterPtr, ClusterPtr> (newClusterMCR,newClusterDCR);
}
std::tuple <ClusterPtr, ClusterPtr, ClusterPtr>
ColourReconnector::_reconnect3Mto3M(ClusterPtr & c1, ClusterPtr & c2, ClusterPtr & c3, const int infos [3]) const {
// check if mesonic clusters
assert(c1->numComponents()==2);
assert(c2->numComponents()==2);
assert(c3->numComponents()==2);
ClusterVector oldclusters = {c1, c2, c3};
ClusterVector newclusters;
for (int i=0; i<3; i++) {
int c1_col=-1;
int c2_anticol=-1;
// get which index is coloured and which anticolour
for (unsigned int ix=0; ix<2; ++ix) {
if (oldclusters[i]->particle(ix)->hasColour(false)) c1_col = ix;
if (oldclusters[infos[i]]->particle(ix)->hasColour(true)) c2_anticol = ix;
}
assert(c1_col>=0);
assert(c2_anticol>=0);
oldclusters[i]->colParticle()->abandonChild(oldclusters[i]);
oldclusters[infos[i]]->antiColParticle()->abandonChild(oldclusters[infos[i]]);
// form new cluster
ClusterPtr newCluster = new_ptr(Cluster(oldclusters[i]->colParticle(), oldclusters[infos[i]]->antiColParticle()));
oldclusters[i]->colParticle()->addChild(newCluster);
oldclusters[infos[i]]->antiColParticle()->addChild(newCluster);
// set new vertex
newCluster->setVertex(0.5*(oldclusters[i]->colParticle()->vertex() +
oldclusters[infos[i]]->antiColParticle()->vertex()));
// set beam remnants for new cluster
if (oldclusters[i]->isBeamRemnant(c1_col)) newCluster->setBeamRemnant(0, true);
if (oldclusters[infos[i]]->isBeamRemnant(c2_anticol)) newCluster->setBeamRemnant(1, true);
newclusters.push_back(newCluster);
}
return std::tuple <ClusterPtr, ClusterPtr, ClusterPtr> (newclusters[0], newclusters[1], newclusters[2]);
}
pair <ClusterPtr, ClusterPtr>
ColourReconnector::_reconnectMBtoMB(ClusterPtr & c1, ClusterPtr & c2, const int s0) const {
// make c1 the mesonic cluster
if (c1->numComponents()==2) {
assert(c2->numComponents()==3);
} else {
return _reconnectMBtoMB(c2,c1,s0);
}
int c1_col=-1;
int c1_anti=-1;
// get which index is coloured and which anticolour
for (unsigned int ix=0; ix<2; ++ix) {
if (c1->particle(ix)->hasColour(false)) c1_col = ix;
else if (c1->particle(ix)->hasColour(true)) c1_anti = ix;
}
assert(c1_col>=0);
assert(c1_anti>=0);
// pointers for the new clusters
ClusterPtr newCluster1;
ClusterPtr newCluster2;
if (c2->particle(0)->hasColour()==true) {
// first case: we have a baryonic clusters
// first make the new mesonic cluster
c1->antiColParticle()->abandonChild(c1);
c2->particleB(s0)->abandonChild(c2);
newCluster1 = new_ptr(Cluster(c1->antiColParticle(), c2->particleB(s0)));
c1->antiColParticle()->addChild(newCluster1);
c2->particleB(s0)->addChild(newCluster1);
// set new vertex
newCluster1->setVertex(0.5*(c1->antiColParticle()->vertex() +
c2->particleB(s0)->vertex()));
// set beam remnants for new cluster
if (c1->isBeamRemnant(c1_anti)) newCluster1->setBeamRemnant(0, true);
if (c2->isBeamRemnant(s0)) newCluster1->setBeamRemnant(1, true);
// then the baryonic one
c1->colParticle()->abandonChild(c1);
c2->particleB((s0+1)%3)->abandonChild(c2);
c2->particleB((s0+2)%3)->abandonChild(c2);
newCluster2 = new_ptr(Cluster(c1->colParticle(), c2->particleB((s0+1)%3), c2->particleB((s0+2)%3)));
c1->colParticle()->addChild(newCluster2);
c2->particleB((s0+1)%3)->addChild(newCluster2);
c2->particleB((s0+2)%3)->addChild(newCluster2);
// set new vertex
newCluster2->setVertex(LorentzPoint());
} else {
// second case we have an antibaryonic cluster
// first make the new mesonic cluster
c1->colParticle()->abandonChild(c1);
c2->particleB(s0)->abandonChild(c2);
newCluster1 = new_ptr(Cluster(c1->colParticle(), c2->particleB(s0)));
c1->colParticle()->addChild(newCluster1);
c2->particleB(s0)->addChild(newCluster1);
// set new vertex
newCluster1->setVertex(0.5*(c1->colParticle()->vertex() +
c2->particleB(s0)->vertex()));
// set beam remnants for new cluster
if (c1->isBeamRemnant(c1_col)) newCluster1->setBeamRemnant(0, true);
if (c2->isBeamRemnant(s0)) newCluster1->setBeamRemnant(1, true);
// then the baryonic one
c1->antiColParticle()->abandonChild(c1);
c2->particleB((s0+1)%3)->abandonChild(c2);
c2->particleB((s0+2)%3)->abandonChild(c2);
newCluster2 = new_ptr(Cluster(c1->antiColParticle(), c2->particleB((s0+1)%3), c2->particleB((s0+2)%3)));
c1->antiColParticle()->addChild(newCluster2);
c2->particleB((s0+1)%3)->addChild(newCluster2);
c2->particleB((s0+2)%3)->addChild(newCluster2);
// set new vertex
newCluster2->setVertex(LorentzPoint());
}
return pair <ClusterPtr, ClusterPtr> (newCluster1, newCluster2);
}
void ColourReconnector::_2Bto2BreconnectionFinder(ClusterPtr & c1, ClusterPtr & c2,
int & bswap1, int & bswap2, double min_displ_sum, string & kind_of_reco) const {
double tmp_delta;
for (int i=0; i<3; i++) {
for (int j=0; j<3; j++) {
// try swapping particle i of c1 with particle j of c2
tmp_delta = _displacementBaryonic(c2->particle(j), c1->particle((i+1)%3), c1->particle((i+2)%3));
tmp_delta += _displacementBaryonic(c1->particle(i), c2->particle((j+1)%3), c2->particle((j+2)%3));
if (tmp_delta < min_displ_sum) {
// if minimal displacement select the 2Bto2B CR option
min_displ_sum = tmp_delta;
bswap1 = i;
bswap2 = j;
kind_of_reco = "2Bto2B";
}
}
}
}
void ColourReconnector::_BbarBto3MreconnectionFinder(ClusterPtr & c1, ClusterPtr & c2, int & mswap0, int & mswap1, int & mswap2,
double min_displ_sum, string & kind_of_reco) const {
double pre_tmp_delta;
double tmp_delta;
for (int p1=0; p1 <3; p1++) {
// make sure not to form a mesonic octet
if (_isColour8(c1->particle(0), c2->particle(p1))) continue;
pre_tmp_delta = _displacement(c1->particle(0), c2->particle(p1));
for (int p2=1; p2<3; p2++) {
// make sure not to form a mesonic octet
if (_isColour8(c1->particle(1), c2->particle((p1+p2)%3))) continue;
if (_isColour8(c1->particle(2), c2->particle(3-p1-((p1+p2)%3)))) continue;
tmp_delta = pre_tmp_delta + _displacement(c1->particle(1), c2->particle((p1+p2)%3));
tmp_delta += _displacement(c1->particle(2), c2->particle(3-p1-((p1+p2)%3)));
// factor _mesonToBaryonFactor to compare Baryonic an mesonic cluster
tmp_delta *=_mesonToBaryonFactor;
if (tmp_delta < min_displ_sum) {
// if minimal displacement select the 2Bto3M CR option
min_displ_sum = tmp_delta;
mswap0 = p1;
mswap1 = (p1+p2)%3;
mswap2 = 3-p1-((p1+p2)%3);
kind_of_reco = "2Bto3M";
}
}
}
}
void ColourReconnector::_BMtoBMreconnectionfinder(ClusterPtr & c1, ClusterPtr & c2, int & swap, double min_displ_sum,
string & kind_of_reco) const {
assert(c1->numComponents()==2);
assert(c2->numComponents()==3);
double tmp_displ = 0;
for (int i=0; i<3; i++) {
// Differ if the second cluster is baryonic or antibaryonic
if (c2->particle(0)->hasColour()) {
// c2 is baryonic
// veto mesonic octets
if (_isColour8(c2->particle(i), c1->antiColParticle())) continue;
// factor _mesonToBaryonFactor to compare Baryonic an mesonic cluster
tmp_displ = _mesonToBaryonFactor * _displacement(c2->particle(i), c1->antiColParticle());
tmp_displ += _displacementBaryonic(c1->colParticle(), c2->particle((i+1)%3), c2->particle((i+2)%3));
} else {
// c2 is antibaryonic
// veto mesonic octets
if (_isColour8(c2->particle(i), c1->colParticle())) continue;
// factor _mesonToBaryonFactor to compare Baryonic an mesonic cluster
tmp_displ = _mesonToBaryonFactor * _displacement(c2->particle(i), c1->colParticle());
tmp_displ *= _displacementBaryonic(c1->antiColParticle(), c2->particle((i+1)%3), c2->particle((i+2)%3));
}
if (tmp_displ < min_displ_sum) {
// if minimal displacement select the MBtoMB CR option
min_displ_sum = tmp_displ;
swap = i;
kind_of_reco = "MBtoMB";
}
}
return;
}
void ColourReconnector::_3MtoXreconnectionfinder(std::vector<CluVecIt> & cv, int & swap0, int & swap1,
int & swap2, double min_displ_sum, string & kind_of_reco) const {
// case of 3M->BbarB CR
double _tmp_displ;
_tmp_displ = _displacementBaryonic((*cv[0])->colParticle(), (*cv[1])->colParticle(), (*cv[2])->colParticle());
_tmp_displ += _displacementBaryonic((*cv[0])->antiColParticle(), (*cv[1])->antiColParticle(), (*cv[2])->antiColParticle());
if (_tmp_displ < min_displ_sum) {
// if minimal displacement select the 3Mto2B CR option
kind_of_reco = "3Mto2B";
min_displ_sum = _tmp_displ;
}
// case for 3M->3M CR
/**
* if 3Mto3M reco probability (_preco3M_3M) is 0 we skip this loop
* since no 3Mto3M CR shall be performed
*/
int i,j;
int i1,i2,i3;
for (i = 0; _preco3M_3M && i<3; i++) {
// veto mesonic octets
if (_isColour8((*cv[0])->colParticle(), (*cv[i])->antiColParticle())) continue;
// factor _mesonToBaryonFactor to compare baryonic an mesonic cluster
_tmp_displ = _mesonToBaryonFactor * _displacement((*cv[0])->colParticle(), (*cv[i])->antiColParticle());
for (j=1; j<3; j++) {
// i1, i2, i3 are pairwise distinct
i1=i;
i2=((j+i)%3);
if (i1==0 && i2==1) continue;
i3=(3-i-((j+i)%3));
// veto mesonic octets
if (_isColour8((*cv[1])->colParticle(), (*cv[i2])->antiColParticle())) continue;
if (_isColour8((*cv[2])->colParticle(), (*cv[i3])->antiColParticle())) continue;
_tmp_displ += _mesonToBaryonFactor * _displacement((*cv[1])->colParticle(), (*cv[i2])->antiColParticle());
_tmp_displ += _mesonToBaryonFactor * _displacement((*cv[2])->colParticle(), (*cv[i3])->antiColParticle());
if (_tmp_displ < min_displ_sum) {
// if minimal displacement select the 3Mto3M CR option
kind_of_reco = "3Mto3M";
min_displ_sum = _tmp_displ;
swap0 = i1;
swap1 = i2;
swap2 = i3;
}
}
}
}
pair <int,int> ColourReconnector::_shuffle
(const PVector & q, const PVector & aq, unsigned maxtries) const {
const size_t nclusters = q.size();
assert (nclusters > 1);
assert (aq.size() == nclusters);
int i, j;
unsigned tries = 0;
bool octet=false;
do {
// find two different random integers in the range [0, nclusters)
i = UseRandom::irnd( nclusters );
do {
j = UseRandom::irnd( nclusters );
} while (i == j);
// check if one of the two potential clusters would be a colour octet state
octet = _isColour8( q[i], aq[j] ) || _isColour8( q[j], aq[i] ) ;
tries++;
} while (octet && tries < maxtries);
if (octet) i = j = -1;
return make_pair(i,j);
}
bool ColourReconnector::_isColour8(tcPPtr p, tcPPtr q) const {
bool octet = false;
// make sure we have a triplet and an anti-triplet
if ( ( p->hasColour() && q->hasAntiColour() ) ||
( p->hasAntiColour() && q->hasColour() ) ) {
// true if p and q are originated from a colour octet
if ( !p->parents().empty() && !q->parents().empty() ) {
octet = ( p->parents()[0] == q->parents()[0] ) &&
( p->parents()[0]->data().iColour() == PDT::Colour8 );
}
// (Final) option: check if same colour8 parent
// or already found an octet.
if(_octetOption==0||octet) return octet;
// (All) option handling more octets
// by browsing particle history/colour lines.
tColinePtr cline,aline;
// Get colourlines form final states.
if(p->hasColour() && q->hasAntiColour()) {
cline = p-> colourLine();
aline = q->antiColourLine();
} else {
cline = q-> colourLine();
aline = p->antiColourLine();
}
// Follow the colourline of p.
if ( !p->parents().empty() ) {
tPPtr parent = p->parents()[0];
while (parent) {
if(parent->data().iColour() == PDT::Colour8) {
// Coulour8 particles should have a colour
// and an anticolour line. Currently the
// remnant has none of those. Since the children
// of the remnant are not allowed to emit currently,
// the colour octet remnant is handled by the return
// statement above. The assert also catches other
// colour octets without clines. If the children of
// a remnant should be allowed to emit, the remnant
// should get appropriate colour lines and
// colour states.
// See Ticket: #407
// assert(parent->colourLine()&&parent->antiColourLine());
octet = (parent-> colourLine()==cline &&
parent->antiColourLine()==aline);
}
if(octet||parent->parents().empty()) break;
parent = parent->parents()[0];
}
}
}
return octet;
}
void ColourReconnector::persistentOutput(PersistentOStream & os) const {
os
<< _hadronSpectrum
<< _clreco
<< _algorithm
<< _annealingFactor
<< _annealingSteps
<< _triesPerStepFactor
<< _initTemp
<< _preco
<< _precoBaryonic
<< _precoDiquark
<< _preco3M_3M
<< _preco3M_BBbar
<< _precoBbarB_3M
<< _preco2B_2B
<< _precoMB_MB
<< _stepFactor
<< _mesonToBaryonFactor
<< ounit(_maxDistance, femtometer)
<< _octetOption
<< _cr2BeamClusters
<< _localCR
<< _causalCR
<< _debug
<< _junctionMBCR
<< _dynamicCR
<< _diquarkCR
<< ounit(_dynamicCRscale, GeV)
<< _dynamicCRalphaS
<< _phaseSpaceDiquarkFission
;
}
void ColourReconnector::persistentInput(PersistentIStream & is, int) {
is
>> _hadronSpectrum
>> _clreco
>> _algorithm
>> _annealingFactor
>> _annealingSteps
>> _triesPerStepFactor
>> _initTemp
>> _preco
>> _precoBaryonic
>> _precoDiquark
>> _preco3M_3M
>> _preco3M_BBbar
>> _precoBbarB_3M
>> _preco2B_2B
>> _precoMB_MB
>> _stepFactor
>> _mesonToBaryonFactor
>> iunit(_maxDistance, femtometer)
>> _octetOption
>> _cr2BeamClusters
>> _localCR
>> _causalCR
>> _debug
>> _junctionMBCR
>> _dynamicCR
>> _diquarkCR
>> iunit(_dynamicCRscale, GeV)
>> _dynamicCRalphaS
>> _phaseSpaceDiquarkFission
;
}
void ColourReconnector::Init() {
static ClassDocumentation<ColourReconnector> documentation
("This class is responsible of the colour reconnection.");
static Reference<ColourReconnector,HadronSpectrum>
interfaceHadronSpectrum("HadronSpectrum",
"A reference to the object HadronSpectrum",
&Herwig::ColourReconnector::_hadronSpectrum,
false, false, true, false);
static Switch<ColourReconnector,int> interfaceColourReconnection
("ColourReconnection",
"Colour reconnections",
&ColourReconnector::_clreco, 0, true, false);
static SwitchOption interfaceColourReconnectionNo
(interfaceColourReconnection,
"No",
"Colour reconnections off",
0);
static SwitchOption interfaceColourReconnectionYes
(interfaceColourReconnection,
"Yes",
"Colour reconnections on",
1);
// Algorithm interface
static Switch<ColourReconnector, int> interfaceAlgorithm
("Algorithm",
"Specifies the colour reconnection algorithm",
&ColourReconnector::_algorithm, 0, true, false);
static SwitchOption interfaceAlgorithmPlain
(interfaceAlgorithm,
"Plain",
"Plain colour reconnection as in Herwig 2.5.0",
0);
static SwitchOption interfaceAlgorithmStatistical
(interfaceAlgorithm,
"Statistical",
"Statistical colour reconnection using simulated annealing",
1);
static SwitchOption interfaceAlgorithmBaryonic
(interfaceAlgorithm,
"Baryonic",
"Baryonic cluster reconnection",
2);
static SwitchOption interfaceAlgorithmBaryonicMesonic
(interfaceAlgorithm,
"BaryonicMesonic",
"Baryonic cluster reconnection with reconnections to and from Mesonic Clusters",
3);
static SwitchOption interfaceAlgorithmBaryonicDiquarkCluster
(interfaceAlgorithm,
"BaryonicDiquarkCluster",
"Baryonic colour reconnection which allows for the formation of DiquarkCluster-like CR",
4);
static Switch<ColourReconnector,int> interfaceColourDiquarkCR
("DiquarkCR",
"Allow diquark type colour Reconnection",
&ColourReconnector::_diquarkCR, 0, true, false);
static SwitchOption interfaceDiquarkCRNo
(interfaceColourDiquarkCR,
"No",
"Use regular CR without diquark CR",
0);
static SwitchOption interfaceDiquarkCRYes
(interfaceColourDiquarkCR,
"Yes",
"Use CR with diquark CR",
1);
static Switch<ColourReconnector,int> interfaceColourDynamicCR
("DynamicCR",
"Use dynamic weight for Colour reconnections defined by soft gluon evolution"
"\nNOTE: Only for Mesonic CR so far",
&ColourReconnector::_dynamicCR, 0, true, false);
static SwitchOption interfaceDynamicCRNo
(interfaceColourDynamicCR,
"No",
"Use regular CR with fixed probabilities",
0);
static SwitchOption interfaceDynamicCRYes
(interfaceColourDynamicCR,
"Yes",
"Use dynamic CR with kinematic dependent probabilities",
1);
// General Parameters and switches
static Parameter<ColourReconnector, Energy> interfaceDynamicScale
("DynamicScale",
"Choose dynamic scale of soft gluon evolution for DynamicCR",
&ColourReconnector::_dynamicCRscale, GeV, 1.*GeV, 0.001*GeV, 1e4*GeV,
false, false, Interface::limited);
static Parameter<ColourReconnector, double> interfaceDynamicAlphaS
("DynamicAlphaS",
"Choose dynamic alphaS of soft gluon evolution for DynamicCR",
&ColourReconnector::_dynamicCRalphaS, 0.8, 0.001, 10.0,
false, false, Interface::limited);
// Statistical CR Parameters:
static Parameter<ColourReconnector, double> interfaceMtrpAnnealingFactor
("AnnealingFactor",
"The annealing factor is the ratio of the temperatures in two successive "
"temperature steps.",
&ColourReconnector::_annealingFactor, 0.9, 0.0, 1.0,
false, false, Interface::limited);
static Parameter<ColourReconnector,unsigned> interfaceMtrpAnnealingSteps
("AnnealingSteps",
"Number of temperature steps in the statistical annealing algorithm",
&ColourReconnector::_annealingSteps, 50, 1, 10000,
false, false, Interface::limited);
static Parameter<ColourReconnector,double> interfaceMtrpTriesPerStepFactor
("TriesPerStepFactor",
"The number of reconnection tries per temperature steps is the number of "
"clusters times this factor.",
&ColourReconnector::_triesPerStepFactor, 5.0, 0.0, 100.0,
false, false, Interface::limited);
static Parameter<ColourReconnector,double> interfaceMtrpInitialTemp
("InitialTemperature",
"Factor used to determine the initial temperature from the median of the "
"energy change in a few random rearrangements.",
&ColourReconnector::_initTemp, 0.1, 0.00001, 100.0,
false, false, Interface::limited);
// Plain and Baryonic CR Paramters
static Parameter<ColourReconnector, double> interfaceRecoProb
("ReconnectionProbability",
"Probability that a found two meson to two meson reconnection possibility is actually accepted (used in Plain & Baryonic)",
&ColourReconnector::_preco, 0.5, 0.0, 1.0,
false, false, Interface::limited);
static Parameter<ColourReconnector,double> interfaceRecoProbBaryonic
("ReconnectionProbabilityBaryonic",
"Probability that a found reconnection possibility is actually accepted (used in Baryonic)",
&ColourReconnector::_precoBaryonic, 0.5, 0.0, 1.0,
false, false, Interface::limited);
static Parameter<ColourReconnector,double> interfaceRecoProbDiquark
("ReconnectionProbabilityDiquark",
"Probability for forming a tetra-quark cluster",
&ColourReconnector::_precoDiquark, 0.5, 0.0, 1.0,
false, false, Interface::limited);
// BaryonicMesonic CR Paramters
static Parameter<ColourReconnector, double> interfaceReconnectionProbability3Mto3M
("ReconnectionProbability3Mto3M",
"Probability that a reconnection candidate is accepted for reconnecting 3M -> 3M\'",
&ColourReconnector::_preco3M_3M, 0.5, 0.0, 1.0,
false, false, Interface::limited);
static Parameter<ColourReconnector, double> interfaceReconnectionProbability3MtoBBbar
("ReconnectionProbability3MtoBBbar",
"Probability that a reconnection candidate is accepted for reconnecting 3M -> B,Bbar",
&ColourReconnector::_preco3M_BBbar, 0.5, 0.0, 1.0,
false, false, Interface::limited);
static Parameter<ColourReconnector, double> interfaceReconnectionProbabilityBbarBto3M
("ReconnectionProbabilityBbarBto3M",
"Probability that a reconnection candidate is accepted for reconnecting B,Bbar -> 3M",
&ColourReconnector::_precoBbarB_3M, 0.5, 0.0, 1.0,
false, false, Interface::limited);
static Parameter<ColourReconnector, double> interfaceReconnectionProbability2Bto2B
("ReconnectionProbability2Bto2B",
"Probability that a reconnection candidate is accepted for reconnecting 2B -> 2B\' or 2Bbar -> 2Bbar\'",
&ColourReconnector::_preco2B_2B, 0.5, 0.0, 1.0,
false, false, Interface::limited);
static Parameter<ColourReconnector, double> interfaceReconnectionProbabilityMBtoMB
("ReconnectionProbabilityMBtoMB",
"Probability that a reconnection candidate is accepted for reconnecting M,B -> M\',B\' or M,Bbar -> M\',Bbar\'",
&ColourReconnector::_precoMB_MB, 0.5, 0.0, 1.0,
false, false, Interface::limited);
static Parameter<ColourReconnector, double> interfaceFactorforStep
("StepFactor",
"Factor for how many reconnection-tries are made in the BaryonicMesonic algorithm",
&ColourReconnector::_stepFactor, 1.0, 0.11111, 10.,
false, false, Interface::limited);// at least 3 Clusters -> _stepFactorMin=1/9
static Parameter<ColourReconnector, double> interfaceMesonToBaryonFactor
("MesonToBaryonFactor",
"Factor for comparing mesonic clusters to baryonic clusters in the displacement if BaryonicMesonic CR model is chosen",
&ColourReconnector::_mesonToBaryonFactor, 2.0, 0.01, 100.0,
false, false, Interface::limited);
// General Parameters and switches
static Parameter<ColourReconnector, Length> interfaceMaxDistance
("MaxDistance",
"Maximum distance between the clusters at which to consider rearrangement"
" to avoid colour reconneections of displaced vertices (used in all Algorithms). No unit means femtometer",
&ColourReconnector::_maxDistance, femtometer, 1000.*femtometer, 0.0*femtometer, 1e100*femtometer,
false, false, Interface::limited);
static Switch<ColourReconnector, unsigned int> interfaceOctetTreatment
("OctetTreatment",
"Which octets are not allowed to be reconnected (used in all Algorithms)",
&ColourReconnector::_octetOption, 0, false, false);
static SwitchOption interfaceOctetTreatmentFinal
(interfaceOctetTreatment,
"Final",
"Only prevent for the final (usuaslly non-perturbative) g -> q qbar splitting",
0);
static SwitchOption interfaceOctetTreatmentAll
(interfaceOctetTreatment,
"All",
"Prevent for all octets",
1);
static Switch<ColourReconnector, int> interfaceCR2BeamClusters
("CR2BeamClusters",
"Option for colour reconnecting 2 beam remnant clusters if the number of clusters is 2.",
&ColourReconnector::_cr2BeamClusters, 0, true, false);
static SwitchOption interfaceCR2BeamClustersYes
(interfaceCR2BeamClusters,
"Yes",
"If possible CR 2 beam clusters",
1);
static SwitchOption interfaceCR2BeamClustersNo
(interfaceCR2BeamClusters,
"No",
"If possible do not CR 2 beam clusters",
0);
static Switch<ColourReconnector, int> interfaceLocalCR
("LocalCR",
"Option for colour reconnecting only if clusters are less distant than MaxDistance",
&ColourReconnector::_localCR, 0, true, false);
static SwitchOption interfaceLocalCRYes
(interfaceLocalCR,
"Yes",
"activate spatial veto",
1);
static SwitchOption interfaceLocalCRNo
(interfaceLocalCR,
"No",
"deactivate spatial veto",
0);
static Switch<ColourReconnector, int> interfaceCausalCR
("CausalCR",
"Option for colour reconnecting only if clusters their vertices "
"have a positive spacetime difference",
&ColourReconnector::_causalCR, 0, true, false);
static SwitchOption interfaceCausalCRYes
(interfaceCausalCR,
"Yes",
"enable causal veto",
1);
static SwitchOption interfaceCausalCRNo
(interfaceCausalCR,
"No",
"disable causal veto",
0);
static Switch<ColourReconnector, int> interfaceJunction
("Junction",
"Option for using Junction-like displacement in rapidity-phi plane to compare baryonic cluster "
"instead of pairwise distance (for BaryonicMesonic model)",
&ColourReconnector::_junctionMBCR, 1, true, false);
static SwitchOption interfaceJunctionYes
(interfaceJunction,
"Yes",
"Using junction-like model instead of pairwise distance model",
1);
static SwitchOption interfaceJunctionNo
(interfaceJunction,
"No",
"Using pairwise distance model instead of junction-like model",
0);
// Debug
static Switch<ColourReconnector, int> interfaceDebug
("Debug",
"Make a file with some Information of the BaryonicMesonic Algorithm",
&ColourReconnector::_debug, 0, true, false);
static SwitchOption interfaceDebugNo
(interfaceDebug,
"No",
"Debug Information for ColourReconnector Off",
0);
static SwitchOption interfaceDebugYes
(interfaceDebug,
"Yes",
"Debug Information for ColourReconnector On",
1);
static Switch<ColourReconnector,int> interfacePhaseSpaceDiquarkFission
("PhaseSpaceDiquarkFission",
"Colour reconnections",
&ColourReconnector::_phaseSpaceDiquarkFission, 0, true, false);
static SwitchOption interfacePhaseSpaceDiquarkFissionNo
(interfacePhaseSpaceDiquarkFission,
"No",
"Not adding the decay phasespace to Diquark Colour Reconnection",
0);
static SwitchOption interfacePhaseSpaceDiquarkFissionYes
(interfacePhaseSpaceDiquarkFission,
"Yes",
"Adding the decay phasespace to Diquark Colour Reconnection",
1);
static SwitchOption interfacePhaseSpaceDiquarkFissionConstituentMasses
(interfacePhaseSpaceDiquarkFission,
"ConstituentMasses",
"Adding the decay phasespace to Diquark Colour Reconnection",
2);
}
diff --git a/Hadronization/ColourReconnector.h b/Hadronization/ColourReconnector.h
--- a/Hadronization/ColourReconnector.h
+++ b/Hadronization/ColourReconnector.h
@@ -1,682 +1,682 @@
// -*- C++ -*-
//
// ColourReconnector.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_ColourReconnector_H
#define HERWIG_ColourReconnector_H
#include <ThePEG/Interface/Interfaced.h>
#include <unordered_map>
#include "CluHadConfig.h"
#include "ColourReconnector.fh"
#include "HadronSpectrum.h"
#include "Herwig/Utilities/Kinematics.h"
#include <bits/stdc++.h>
namespace Herwig {
using namespace ThePEG;
/** \ingroup Hadronization
* \class ColourReconnector
* \brief Class for changing colour reconnections of partons.
* \author Alberto Ribon, Christian Roehr
*
* This class does the nonperturbative colour rearrangement, after the
* nonperturbative gluon splitting and the "normal" cluster formation.
* It uses the list of particles in the event record, and the collections of
* "usual" clusters which is passed to the main method. If the colour
* reconnection is actually accepted, then the previous collections of "usual"
* clusters is first deleted and then the new one is created.
*
* * @see \ref ColourReconnectorInterfaces "The interfaces"
* defined for ColourReconnector.
*/
class ColourReconnector: public Interfaced {
public:
/**
* Does the colour rearrangement, starting out from the list of particles in
* the event record and the collection of "usual" clusters passed as
* arguments. If the actual rearrangement is accepted, the initial collection of
* clusters is overridden by the old ones.
*/
void rearrange(ClusterVector & clusters);
using CluVecIt = ClusterVector::iterator;
private:
/*
class ColourFlowState2 {
};
class ColourFlowBasis2 {
// represents the colour flow state for 2 colour flows:
// Note: i,j are either (1,2)
// |ij> === |(1,2)->(i,j)>
// where 1,2 are coloured partons and i,j are their
// anticoloured partons
public:
ColourFlowBasis2();
ColourFlowBasis2(const ColourFlowBasis2 & cf2){
setIndex(cf2.getIndex(0),cf2.getIndex(1));
};
bool operator=(const ColourFlowBasis2& cf2){
return setIndex(cf2.getIndex(0),cf2.getIndex(1));
}
bool setIndex(int ind[2]) {
return setIndex(ind[0],ind[1]);
};
bool setIndex(int i, int j) {
if (i==j || i<0 || j<0 || i>1 || j>1){
std::cout << "Error cannot create ColourFlowBasis2 = |"<<i+1<<j+1<<">\n";
return true;
}
_index[0]=i+1;
_index[1]=j+1;
return false;
};
int getIndex(int i) const {
if (i<0 || i>1 ){
std::cout << "Error i = "<<i << std::endl;
assert(false);
}
return _index[i];
};
private:
int _index[2] = {1,2};
};
class ColourFlowBasis3 {
// represents the colour flow state for 2 colour flows:
// Note: i,j,k are either (1,2,3)
// |ijk> === |(1,2,3)->(i,j,k)>
// where 1,2,3 are coloured partons and i,j,k are their
// anticoloured partons
public:
ColourFlowBasis3();
ColourFlowBasis3(const ColourFlowBasis3 & cf3){
setIndex(cf3.getIndex(0),cf3.getIndex(1),cf3.getIndex(2));
};
bool operator=(const ColourFlowBasis3& cf3){
return setIndex(cf3.getIndex(0),cf3.getIndex(1),cf3.getIndex(2));
}
bool setIndex(int ind[3]) {
return setIndex(ind[0],ind[1],ind[2]);
};
bool setIndex(int i, int j, int k) {
if (i==j || i==k || j==k || i<0 || j<0 || k<0 || i>2 || j>2 || k>2){
std::cout << "Error cannot create ColourFlowBasis3 = |"<<i+1<<j+1<<k+1<<">\n";
return true;
}
_index[0]=i+1;
_index[1]=j+1;
_index[2]=k+1;
return false;
};
int getIndex(int i) const {
if (i<0 || i>2 ){
std::cout << "Error i = "<<i << std::endl;
assert(false);
}
return _index[i];
};
private:
int _index[3] = {1,2,3};
};*/
HadronSpectrumPtr _hadronSpectrum;
/** PRIVATE MEMBER FUNCTIONS */
/**
* @brief Calculates the sum of the squared cluster masses.
* @arguments q, aq vectors containing the quarks and antiquarks respectively
* @return Sum of cluster squared masses M^2_{q[i],aq[i]}.
*/
Energy2 _clusterMassSum(const PVector & q, const PVector & aq) const;
Selector <int> getProbabilities2CF(const ClusterPtr & c1, const ClusterPtr & c2, bool diquarkCR) const;
Selector <int> _selector(const ClusterVector & clusters, bool diquarkCR) const;
Selector <int> _selectorCF3(const ClusterVector & clusters, bool diquarkCR) const;
std::unordered_map<int,double> _reconnectionAmplitudesSGE(const ClusterVector & clusters) const;
int _stateToPermutation(const int i) const;
/**
* @brief Calculates the Mesonic reconnection Probability from soft gluon evolution
* @arguments c1, c2 cluster pointer refs, whose kinematics determine the probability
* @return probability in [0:1]
*/
std::tuple<double,double,double> _dynamicRecoProbabilitiesCF2(const ClusterPtr & c1, const ClusterPtr & c2, bool diquarkCR) const;
std::unordered_map<int,double> _reconnectionAmplitudesCF2 (const ClusterPtr & c1, const ClusterPtr & c2) const;
/**
* @brief calculates the "euclidean" distance of two quarks in the
* rapdidity-phi plane
* @arguments p, q the two quarks
* @return the dimensionless distance:
* \deltaR_12=sqrt(\deltaY_12^2 + \delta\phi_12^2)
*/
double _displacement(tcPPtr p, tcPPtr q) const;
/**
* @brief calculates the "euclidean" distance of a the 3 (anti)quarks
* (anti)baryonic cluster in the rapdidity-phi plane
* @arguments p, q the two quarks
* @return the dimensionless distance:
* if Junction is enabled the difference of all 3 quarks
* with respect to their mean point is calculated:
* <Y> = sum_i Y_i/3
* <\phi> = sum_i \phi_i/3
* \deltaR = sum_i sqrt( (Y_i - <Y>)^2 + (\phi_i - <phi>)^2)
*
* if Junction is disabled the difference of all distinct
* pairing of the 3 quarks is computed:
* \deltaR_ij = sqrt(\deltaY_ij^2 + \delta\phi_ij^2)
* \deltaR_tot = sum_{i,j<i} \deltaR_ij
*
* NOTE: switch the Junction option will necessarily
* need to be combined with a change (or re-tune)
* of _mesonToBaryonFactor otherwise no well
* description is to be expected
* TODO: maybe add a different p-norm option to get
* more phenomenology
*/
double _displacementBaryonic(tcPPtr p1, tcPPtr p2, tcPPtr p3) const;
/**
* @brief Examines whether the cluster vector (under the given permutation of
* the antiquarks) contains colour-octet clusters
* @param cv Cluster vector
* @param P Permutation, a vector of permutated indices from 0 to
* cv.size()-1
*/
bool _containsColour8(const ClusterVector & cv, const vector<size_t> & P) const;
/**
* @brief A Metropolis-type algorithm which finds a local minimum in the
* total sum of cluster masses
* @arguments cv cluster vector
*/
void _doRecoStatistical(ClusterVector & cv) const;
/**
* @brief Plain colour reconnection as used in Herwig 2.5.0
* @arguments cv cluster vector
*/
void _doRecoPlain(ClusterVector & cv) const;
/**
* Baryonic Colour Reconnection model
*/
void _doRecoBaryonic(ClusterVector & cv) const;
void _doRecoBaryonicDiquarkCluster(ClusterVector & cv) const;
/**
* @brief BaryonicMesonic colour reconnection model
* @arguments cv cluster vector
* BaryonicMesonic Colour Reconnection model with reconnections from mesonic clusters
* to baryonic cluster and the contrary. Allows also reconnection between mesonic
* and baryonic Clusters
*/
void _doRecoBaryonicMesonic(ClusterVector & cv) const;
void _makeBaryonicClusters(ClusterPtr &c1, ClusterPtr &c2, ClusterPtr &c3,
ClusterPtr &newcluster1, ClusterPtr &newcluster2) const;
/**
* Method to allow one to make four (light) quark clusters
*/
bool _makeDiquarkCluster(ClusterPtr &c1, ClusterPtr &c2, ClusterPtr &newcluster) const;
bool _canMakeDiquarkCluster(tcPPtr pCol1, tcPPtr pCol2,tcPPtr pAntiCol1, tcPPtr pAntiCol2) const;
bool _canMakeDiquarkCluster(const ClusterPtr &c1, const ClusterPtr &c2) const;
bool _canMakeDiquarkCluster(tcPPtr pCol1, tcPPtr pCol2,tcPPtr pAntiCol1, tcPPtr pAntiCol2, double &PhaseSpace) const;
bool _canMakeDiquarkCluster(const ClusterPtr &c1, const ClusterPtr &c2, double &PhaseSpace) const;
/**
* @brief Finds the cluster in cv which, if reconnected with the given
* cluster cl, would result in the smallest sum of cluster masses.
* If no reconnection partner can be found, a pointer to the
* original Cluster cl is returned.
* @arguments cv cluster vector
* cl cluster iterator (must be from cv) which wants to have a reconnection partner
* @return iterator to the found cluster, or the original cluster pointer if
* no mass-reducing combination can be found
*/
CluVecIt _findRecoPartnerPlain(CluVecIt cl, ClusterVector & cv, const ClusterVector & deleted) const;
- CluVecIt _findRecoPartnerPlainDynamic(CluVecIt cl, ClusterVector & cv, const ClusterVector & deleted) const;
+ CluVecIt _findRecoPartnerPlainDynamic(CluVecIt cl, ClusterVector & cv, const ClusterVector & deleted, bool diquarkCR) const;
CluVecIt _findPartnerRapidity(CluVecIt cl, ClusterVector & cv) const;
CluVecIt _findPartnerBaryonic(CluVecIt cl, ClusterVector & cv,
bool & tetraCand,
const ClusterVector& a,
CluVecIt &baryonic1,
CluVecIt &baryonic2 ) const;
void _findPartnerBaryonicDiquarkClusterTEST( CluVecIt cl, ClusterVector & cv,
unsigned & typeOfReconnection,
const ClusterVector& deleted,
CluVecIt &candidate1,
CluVecIt &candidate2 ) const;
void _findPartnerBaryonicDiquarkCluster( CluVecIt cl, ClusterVector & cv,
unsigned & typeOfReconnection,
const ClusterVector& deleted,
CluVecIt &candidate1,
CluVecIt &candidate2 ) const;
/**
* @brief Finds best CR option for the BaryonicMesonic CR model
* @arguments cls vector of selected clusters, baryonic is the number of baryonic
* clusters in selection, kind_of_reco is output string denoting best
* CR option and reco_info is output storage for information on which
* (anti-)quarks to exchange and in which way.
* BaryonicMesonic Colour Reconnection model with reconnections from mesonic clusters
* to baryonic cluster and the contrary. Allows also reconnection between mesonic
* and baryonic Clusters
*/
void _findbestreconnectionoption(std::vector<CluVecIt> &cls,
const unsigned &baryonic,
string &kind_of_reco,
int (&reco_info)[3]) const;
/**
* @brief Reconnects the constituents of the given clusters to the (only)
* other possible cluster combination.
* @return pair of pointers to the two new clusters
* Used for Plain and Baryonic Colour Reconnection models
*/
pair <ClusterPtr,ClusterPtr> _reconnect(ClusterPtr &c1, ClusterPtr &c2) const;
std::tuple <ClusterPtr, ClusterPtr> _reconnect3MtoMD(ClusterVector & cluvec, const int topology) const;
/**
* @brief Reconnects (2B->2B) the constituents of the given clusters to
another possible cluster combination whose information is given
in s1 and s2.
* @arguments c1 and c2 are baryonic clusters and s1 and s2 are the respective
indices of the constituents of c1 and c2 respectively
* @return pair of pointers to the two new clusters
* Used only in BaryonicMesonic algorithm and will exchange constituent s1 of
* c1 with constituent s2 of c2
*/
pair <ClusterPtr,ClusterPtr> _reconnect2Bto2B(ClusterPtr &c1, ClusterPtr &c2, const int s1, const int s2) const;
/**
* @brief Reconnects (B,Bbar->3M) the constituents of the given clusters to
another possible cluster combination whose information is given in
s0, s1 and s2.
* @arguments c1 and c2 are baryonic clusters. s0, s1 and s2 are the respective
indices which determine the CR
* @return tuple of pointers to the 3 new mesonic clusters
* Used only in BaryonicMesonic algorithm and will form 3 mesonic clusters according
* to the indices s0, s1 and s2. The i-th constituent of c1 is connected to the si-th
* constituent of c2
*/
std::tuple <ClusterPtr, ClusterPtr, ClusterPtr> _reconnectBbarBto3M(ClusterPtr &c1, ClusterPtr &c2, const int s0, const int s1, const int s2 ) const;
/**
* @brief Reconnects (3M->3M) the constituents of the given clusters to
another possible cluster combination whose information is given in
sinfos
* @arguments c1, c2 and c3 are mesonic clusters. infos[3] are the respective
indices which determine the CR
* @return tuple of pointers to the 3 CR'd mesonic clusters
* Used only in BaryonicMesonic algorithm and will reconnect 3 mesonic clusters according
* to the infos, which determine the CR. The coloured quark of the i-th cluster forms
* a new cluster with the anticoloured quark of the info[i]-th cluster
*/
std::tuple <ClusterPtr, ClusterPtr, ClusterPtr> _reconnect3Mto3M(ClusterPtr &c1, ClusterPtr &c2, ClusterPtr &c3, const int infos[3] ) const;
/**
* Reconnection method for a Baryonic and a Mesonic Cluster to a Baryonic and a Mesonic Cluster
* s0 is the Number of the (Anti)Patrticle of the Baryonic Cluster , which should be swapped with the Anti(Particle) of the Mesonic Cluster
*/
/**
* @brief Reconnects the constituents of the given clusters to
another possible cluster combination whose information
is given in s0.
* @arguments c1 and c2 are one baryonic and one mesonic cluster respectively
and s0 is the respective index of the constituents of the baryonic
cluster which is to be exchangeed with the mesonic cluster.
* @return pair of pointers to the two new clusters
* Used only in BaryonicMesonic algorithm and will exchange constituent s0 of
* the baryonic cluster with the (anti)-quark of the mesonic cluster
*/
pair <ClusterPtr, ClusterPtr> _reconnectMBtoMB(ClusterPtr &c1, ClusterPtr &c2, const int s0) const;
/**
* Methods for the BaryonicMesonic CR algorithm
* Find the best reconnection option for the respective cluster-combination
*
*/
/**
* @brief veto for far apart clusters
* @arguments expects at most 3 CluVecIt in clu vector
* @return returns true if clusters are more distant than _maxDistance
* in space
* TODO: problematic maybe add option to turn off
*/
bool _clustersFarApart( const std::vector<CluVecIt> & clu ) const;
/**
* @brief Does reconnect 2 beam clusters for BaryonicMesonic CR model
* if option CR2BeamClusters is enabled
* @arguments cv cluster vector
*/
void _doReco2BeamClusters(ClusterVector & cv) const;
/**
* @brief finds the best reconnection option and stores it in bswap1
* and bswap2 (2B->2B colour reconnection)
* @arguments c1 and c2 cluster pointers and kind_of_reco will output
* the best found reconnection option for c1 and c2
*/
void _2Bto2BreconnectionFinder(ClusterPtr &c1, ClusterPtr &c2, int &bswap1, int &bswap2, double mindisplsum, string &kind_of_reco ) const;
/**
* @brief finds the best reconnection option and stores it in mswap0
* mswap1 and mswap2 (BbarB->3M colour reconnection)
* @arguments c1 and c2 cluster pointers and kind_of_reco will output
* the best found reconnection option for c1 and c2
*/
void _BbarBto3MreconnectionFinder(ClusterPtr &c1, ClusterPtr &c2, int &mswap0, int &mswap1, int &mswap2,
double min_displ_sum, string & kind_of_reco) const;
/**
* @brief finds the best reconnection option and stores it in swap
* (BM->BM colour reconnection)
* @arguments c1 and c2 cluster pointers and kind_of_reco will output
* the best found reconnection option for c1 and c2
*/
void _BMtoBMreconnectionfinder(ClusterPtr &c1, ClusterPtr &c2, int &swap,
double min_displ_sum, string &kind_of_reco) const;
/**
* @brief finds the best reconnection option and stores it in swap0,
* swap1 and swap2 (3M->{3M,BbarB} colour reconnection)
* @arguments c1 and c2 cluster pointers and kind_of_reco will output
* the best found reconnection option for c1 and c2
*/
void _3MtoXreconnectionfinder(std::vector<CluVecIt> &cv, int &swap0, int &swap1,
int &swap2, double min_displ_sum, string &kind_of_reco) const;
/**
* @brief At random, swap two antiquarks, if not excluded by the
* constraint that there must not be any colour-octet clusters.
* @arguments q, aq vectors containing the quarks and antiquarks respectively
* maxtries maximal number of tries to find a non-colour-octet
* reconfiguration
* @return Pair of ints indicating the indices of the antiquarks to be
* swapped. Returns (-1,-1) if no valid reconfiguration could be
* found after maxtries trials
*/
pair <int,int>
_shuffle(const PVector & q, const PVector & aq, unsigned maxtries = 10) const;
/** DATA MEMBERS */
/**
* Specifies the colour reconnection algorithm to be used.
*/
int _algorithm = 0;
/**
* Do we do colour reconnections?
*/
int _clreco = 0;
/**
* Do we want debug informations?
*/
int _debug = 0;
/**
* @brief Junction-like model for BaryonicMesonic model
* Do we want to use the junction-like model for
* computing the displacements of BaryonicMesonic model?
* otherwise pairwise distances are used.
* If _junctionMBCR is activated the displacements are computed in the
* rapidity-Phi plane by difference to the average rapidity and phi:
* DeltaR_i^2 = (rapidity_i - meanRap)^2 + (phi_i - meanPhi)^2
* DeltaR = sum_i DeltaR_i
* if _junctionMBCR=0 the displacements are computed:
* DeltaR_ij^2 = (rapidity_i - rapidity_j)^2 + (phi_i - phi_j)^2
* DeltaR = sum_i,j<i DeltaR_ij
*/
int _junctionMBCR = 1;
/**
* Statistical Reco:
* Factor used to determine the initial temperature according to
* InitialTemperature = _initTemp * median {energy changes in a few random
* rearrangements}
*/
double _initTemp = 0.01;
/**
* Statistical Reco:
* The annealing factor is the ratio of two successive temperature steps:
* T_n = _annealingFactor * T_(n-1)
*/
double _annealingFactor = 0.21;
/**
* Statistical Reco:
* Number of temperature steps in the statistical annealing algorithm
*/
unsigned int _annealingSteps = 10;
/**
* Statistical Reco:
* The number of tries per temperature steps is the number of clusters times
* this factor.
*/
double _triesPerStepFactor = 0.66;
/**
* Probability that a found reconnection possibility is actually accepted.
* used in Plain & Baryonic CR
*/
double _preco = 0.5;
/**
* Probability that a found reconnection possibility is actually accepted.
* used in Baryonic CR
*/
double _precoBaryonic = 0.5;
/**
* Probability that a found reconnection possibility is actually accepted.
* used in Baryonic CR
*/
double _precoDiquark = 0.5;
/**
* Probability that a found reconnection possibility is actually accepted.
* For reconnecting 3M -> 3M'
* used in BaryonicMesonic
* NOTE: if 0 this type of reconnection is not even tried
*/
double _preco3M_3M = 0.5;
/**
* Probability that a found reconnection possibility is actually accepted.
* For reconnecting 3M -> B,Bbar
* used in BaryonicMesonic
*/
double _preco3M_BBbar = 0.5;
/**
* Probability that a found reconnection possibility is actually accepted.
* For reconnecting Bbar,B -> 3M
* used in BaryonicMesonic
*/
double _precoBbarB_3M = 0.5;
/**
* Probability that a found reconnection possibility is actually accepted.
* For reconnecting 2B -> 2B' or 2Bbar -> 2Bbar'
* used in BaryonicMesonic
* NOTE: if 0 this type of reconnection is not even tried
*/
double _preco2B_2B = 0.5;
/**
* Probability that a found reconnection possibility is actually accepted.
* For reconnecting M,B -> M',B' or M,Bbar -> M',Bbar'
* used in BaryonicMesonic
* NOTE: if 0 this type of reconnection is not even tried
*/
double _precoMB_MB = 0.5;
/**
* For the BaryonicMesonic algorithm
* How many times do suggest cluster for reconnection?
* n(reconnectionstries) = _stepFactor * n(clusters)*n(clusters);
*/
double _stepFactor = 5.0;
/**
* Factor for comparing mesonic clusters to baryonic clusters
*/
double _mesonToBaryonFactor = 2.0;
/**
* Maximium distance for reconnections
* TODO: remove if issues with anticausality are discussed and resolved
*/
Length _maxDistance = femtometer;
/**
* @return true, if the two partons are splitting products of the same
* gluon
*/
bool _isColour8(tcPPtr p, tcPPtr q) const;
/**
* Option for handling octets
*/
unsigned int _octetOption = 0;
/**
* Option for colour reconnecting 2 Beam Clusters if no others are present
*/
int _cr2BeamClusters = 0;
/**
* Option for colour reconnecting Clusters only if their vertex 3-distance
* is less than _maxDistance
*/
int _localCR = 0;
/**
* Option for colour reconnecting Clusters only if their spacetime difference
* is bigger than 0
*/
int _causalCR = 0;
/**
* Option for doing the Dynamic CR probability depending on soft
* gluon evolution
*/
int _dynamicCR = 0;
/**
* Option for doing the diquark colour Reconnection
*/
int _diquarkCR = 0;
/**
* Choose Dynamic CR probability scale depending on soft
* gluon evolution
*/
Energy _dynamicCRscale = 1.0*GeV;
/**
* Choose Dynamic CR probability scale depending on soft
* gluon evolution
*/
double _dynamicCRalphaS = 0.8;
/**
* switch for including PhaseSpace in
* Diquark transition probability
* */
int _phaseSpaceDiquarkFission = 1;
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* Standard Init function used to initialize the interfaces.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
private:
/**
* Private and non-existent assignment operator.
*/
ColourReconnector & operator=(const ColourReconnector &) = delete;
};
}
#endif /* HERWIG_ColourReconnector_H */

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 5:31 PM (1 d, 14 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805403
Default Alt Text
(277 KB)

Event Timeline