diff --git a/Hadronization/ClusterFissioner.cc b/Hadronization/ClusterFissioner.cc
--- a/Hadronization/ClusterFissioner.cc
+++ b/Hadronization/ClusterFissioner.cc
@@ -1,2135 +1,2139 @@
// -*- C++ -*-
//
// ClusterFissioner.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// Thisk is the implementation of the non-inlined, non-templated member
// functions of the ClusterFissioner class.
//
#include "ClusterFissioner.h"
#include <ThePEG/Interface/ClassDocumentation.h>
#include <ThePEG/Interface/Reference.h>
#include <ThePEG/Interface/Parameter.h>
#include <ThePEG/Interface/Switch.h>
#include <ThePEG/Persistency/PersistentOStream.h>
#include <ThePEG/Persistency/PersistentIStream.h>
#include <ThePEG/PDT/EnumParticles.h>
#include "Herwig/Utilities/Kinematics.h"
#include "Cluster.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include <ThePEG/Utilities/DescribeClass.h>
#include "ThePEG/Interface/ParMap.h"
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/lu.hpp>
using namespace Herwig;
DescribeClass<ClusterFissioner,Interfaced>
describeClusterFissioner("Herwig::ClusterFissioner","Herwig.so");
ClusterFissioner::ClusterFissioner() :
_clMaxLight(3.35*GeV),
_clMaxExotic(3.35*GeV),
_clPowLight(2.0),
_clPowExotic(2.0),
_pSplitLight(1.0),
_pSplitExotic(1.0),
_phaseSpaceWeights(false),
_dim(4),
_fissionCluster(0),
_kinematicThresholdChoice(0),
_pwtDIquark(0.0),
_diquarkClusterFission(0),
_btClM(1.0*GeV),
_iopRem(1),
_kappa(1.0e15*GeV/meter),
_kinematics(0),
_fluxTubeWidth(1.0*GeV),
_enhanceSProb(0),
_m0Fission(2.*GeV),
_massMeasure(0),
_probPowFactor(4.0),
_probShift(0.0),
_kinThresholdShift(1.0*sqr(GeV)),
_strictDiquarkKinematics(0),
_covariantBoost(false)
{}
IBPtr ClusterFissioner::clone() const {
return new_ptr(*this);
}
IBPtr ClusterFissioner::fullclone() const {
return new_ptr(*this);
}
void ClusterFissioner::persistentOutput(PersistentOStream & os) const {
os << ounit(_clMaxLight,GeV)
<< ounit(_clMaxHeavy,GeV) << ounit(_clMaxExotic,GeV) << _clPowLight << _clPowHeavy
<< _clPowExotic << _pSplitLight
<< _pSplitHeavy << _pSplitExotic
<< _fissionCluster << _fissionPwt
<< _pwtDIquark
<< _diquarkClusterFission
<< _kinematics
<< ounit(_fluxTubeWidth,GeV)
<< ounit(_btClM,GeV)
<< _iopRem << ounit(_kappa, GeV/meter)
<< _enhanceSProb << ounit(_m0Fission,GeV) << _massMeasure << _dim << _phaseSpaceWeights
<< _hadronSpectrum
<< _probPowFactor << _probShift << ounit(_kinThresholdShift,sqr(GeV))
<< _strictDiquarkKinematics
<< _covariantBoost;
}
void ClusterFissioner::persistentInput(PersistentIStream & is, int) {
is >> iunit(_clMaxLight,GeV)
>> iunit(_clMaxHeavy,GeV) >> iunit(_clMaxExotic,GeV) >> _clPowLight >> _clPowHeavy
>> _clPowExotic >> _pSplitLight
>> _pSplitHeavy >> _pSplitExotic
>> _fissionCluster >> _fissionPwt
>> _pwtDIquark
>> _diquarkClusterFission
>> _kinematics
>> iunit(_fluxTubeWidth,GeV)
>> iunit(_btClM,GeV)
>> _iopRem >> iunit(_kappa, GeV/meter)
>> _enhanceSProb >> iunit(_m0Fission,GeV) >> _massMeasure >> _dim >> _phaseSpaceWeights
>> _hadronSpectrum
>> _probPowFactor >> _probShift >> iunit(_kinThresholdShift,sqr(GeV))
>> _strictDiquarkKinematics
>> _covariantBoost;
}
void ClusterFissioner::doinit() {
Interfaced::doinit();
for ( const long& id : spectrum()->heavyHadronizingQuarks() ) {
if ( _pSplitHeavy.find(id) == _pSplitHeavy.end() ||
_clPowHeavy.find(id) == _clPowHeavy.end() ||
_clMaxHeavy.find(id) == _clMaxHeavy.end() )
throw InitException() << "not all parameters have been set for heavy quark cluster fission";
}
// for default Pwts not needed to initialize
if (_fissionCluster==0) return;
for ( const long& id : spectrum()->lightHadronizingQuarks() ) {
if ( _fissionPwt.find(id) == _fissionPwt.end() )
// check that all relevant weights are set
throw InitException() << "fission weights for light quarks have not been set";
}
/*
// Not needed since we set Diquark weights from quark weights
for ( const long& id : spectrum()->lightHadronizingDiquarks() ) {
if ( _fissionPwt.find(id) == _fissionPwt.end() )
throw InitException() << "fission weights for light diquarks have not been set";
}*/
double pwtDquark=_fissionPwt.find(ParticleID::d)->second;
double pwtUquark=_fissionPwt.find(ParticleID::u)->second;
double pwtSquark=_fissionPwt.find(ParticleID::s)->second;
// ERROR: TODO makeDiquarkID is protected function ?
// _fissionPwt[spectrum()->makeDiquarkID(ParticleID::d,ParticleID::d,3)] = _pwtDIquark * pwtDquark * pwtDquark;
// _fissionPwt[spectrum()->makeDiquarkID(ParticleID::u,ParticleID::d,1)] = 0.5 * _pwtDIquark * pwtUquark * pwtDquark;
// _fissionPwt[spectrum()->makeDiquarkID(ParticleID::u,ParticleID::u,3)] = _pwtDIquark * pwtUquark * pwtUquark;
// _fissionPwt[spectrum()->makeDiquarkID(ParticleID::s,ParticleID::d,1)] = 0.5 * _pwtDIquark * pwtSquark * pwtDquark;
// _fissionPwt[spectrum()->makeDiquarkID(ParticleID::s,ParticleID::u,1)] = 0.5 * _pwtDIquark * pwtSquark * pwtUquark;
// _fissionPwt[spectrum()->makeDiquarkID(ParticleID::s,ParticleID::s,3)] = _pwtDIquark * pwtSquark * pwtSquark;
// TODO magic number alternative
_fissionPwt[1103] = _pwtDIquark * pwtDquark * pwtDquark;
_fissionPwt[2101] = 0.5 * _pwtDIquark * pwtUquark * pwtDquark;
_fissionPwt[2203] = _pwtDIquark * pwtUquark * pwtUquark;
_fissionPwt[3101] = 0.5 * _pwtDIquark * pwtSquark * pwtDquark;
_fissionPwt[3201] = 0.5 * _pwtDIquark * pwtSquark * pwtUquark;
_fissionPwt[3303] = _pwtDIquark * pwtSquark * pwtSquark;
}
void ClusterFissioner::Init() {
static ClassDocumentation<ClusterFissioner> documentation
("Class responsibles for chopping up the clusters");
static Reference<ClusterFissioner,HadronSpectrum> interfaceHadronSpectrum
("HadronSpectrum",
"Set the Hadron spectrum for this cluster fissioner.",
&ClusterFissioner::_hadronSpectrum, false, false, true, false);
// ClMax for light, Bottom, Charm and exotic (e.g. Susy) quarks
static Parameter<ClusterFissioner,Energy>
interfaceClMaxLight ("ClMaxLight","cluster max mass for light quarks (unit [GeV])",
&ClusterFissioner::_clMaxLight, GeV, 3.35*GeV, ZERO, 10.0*GeV,
false,false,false);
static ParMap<ClusterFissioner,Energy> interfaceClMaxHeavy
("ClMaxHeavy",
"ClMax for heavy quarkls",
&ClusterFissioner::_clMaxHeavy, GeV, -1, 3.35*GeV, ZERO, 10.0*GeV,
false, false, Interface::upperlim);
static Parameter<ClusterFissioner,Energy>
interfaceClMaxExotic ("ClMaxExotic","cluster max mass for exotic quarks (unit [GeV])",
&ClusterFissioner::_clMaxExotic, GeV, 3.35*GeV, ZERO, 10.0*GeV,
false,false,false);
// ClPow for light, Bottom, Charm and exotic (e.g. Susy) quarks
static Parameter<ClusterFissioner,double>
interfaceClPowLight ("ClPowLight","cluster mass exponent for light quarks",
&ClusterFissioner::_clPowLight, 0, 2.0, 0.0, 10.0,false,false,false);
static ParMap<ClusterFissioner,double> interfaceClPowHeavy
("ClPowHeavy",
"ClPow for heavy quarks",
&ClusterFissioner::_clPowHeavy, -1, 1.0, 0.0, 10.0,
false, false, Interface::upperlim);
static Parameter<ClusterFissioner,double>
interfaceClPowExotic ("ClPowExotic","cluster mass exponent for exotic quarks",
&ClusterFissioner::_clPowExotic, 0, 2.0, 0.0, 10.0,false,false,false);
// PSplit for light, Bottom, Charm and exotic (e.g. Susy) quarks
static Parameter<ClusterFissioner,double>
interfacePSplitLight ("PSplitLight","cluster mass splitting param for light quarks",
&ClusterFissioner::_pSplitLight, 0, 1.0, 0.0, 10.0,false,false,false);
static ParMap<ClusterFissioner,double> interfacePSplitHeavy
("PSplitHeavy",
"PSplit for heavy quarks",
&ClusterFissioner::_pSplitHeavy, -1, 1.0, 0.0, 10.0,
false, false, Interface::upperlim);
static Parameter<ClusterFissioner,double>
interfacePSplitExotic ("PSplitExotic","cluster mass splitting param for exotic quarks",
&ClusterFissioner::_pSplitExotic, 0, 1.0, 0.0, 10.0,false,false,false);
static Switch<ClusterFissioner,int> interfaceFission
("Fission",
"Option for different Fission options",
&ClusterFissioner::_fissionCluster, 1, false, false);
static SwitchOption interfaceFissionDefault
(interfaceFission,
"default",
"Normal cluster fission which depends on the hadron selector class.",
0);
static SwitchOption interfaceFissionNew
(interfaceFission,
"new",
"Alternative cluster fission which does not depend on the hadron selector class",
1);
static Switch<ClusterFissioner,int> interfaceDiquarkClusterFission
("DiquarkClusterFission",
"Allow clusters to fission to 1 or 2 diquark Clusters",
&ClusterFissioner::_diquarkClusterFission, 0, false, false);
static SwitchOption interfaceDiquarkClusterFissionAll
(interfaceDiquarkClusterFission,
"All",
"Allow diquark clusters and baryon clusters to fission to new diquark Clusters",
2);
static SwitchOption interfaceDiquarkClusterFissionOnlyBaryonClusters
(interfaceDiquarkClusterFission,
"OnlyBaryonClusters",
"Allow only baryon clusters to fission to new diquark Clusters",
1);
static SwitchOption interfaceDiquarkClusterFissionNo
(interfaceDiquarkClusterFission,
"No",
"Don't allow clusters to fission to new diquark Clusters",
0);
static ParMap<ClusterFissioner,double> interfaceFissionPwt
("FissionPwt",
"The weights for quarks in the fission process.",
&ClusterFissioner::_fissionPwt, -1, 1.0, 0.0, 10.0,
false, false, Interface::upperlim);
static Switch<ClusterFissioner,int> interfaceRemnantOption
("RemnantOption",
"Option for the treatment of remnant clusters",
&ClusterFissioner::_iopRem, 1, false, false);
static SwitchOption interfaceRemnantOptionSoft
(interfaceRemnantOption,
"Soft",
"Both clusters produced in the fission of the beam cluster"
" are treated as soft clusters.",
0);
static SwitchOption interfaceRemnantOptionHard
(interfaceRemnantOption,
"Hard",
"Only the cluster containing the remnant is treated as a soft cluster.",
1);
static SwitchOption interfaceRemnantOptionVeryHard
(interfaceRemnantOption,
"VeryHard",
"Even remnant clusters are treated as hard, i.e. all clusters the same",
2);
static Parameter<ClusterFissioner,Energy> interfaceBTCLM
("SoftClusterFactor",
"Parameter for the mass spectrum of remnant clusters",
&ClusterFissioner::_btClM, GeV, 1.*GeV, 0.1*GeV, 10.0*GeV,
false, false, Interface::limited);
static Parameter<ClusterFissioner,Tension> interfaceStringTension
("StringTension",
"String tension used in vertex displacement calculation",
&ClusterFissioner::_kappa, GeV/meter,
1.0e15*GeV/meter, ZERO, ZERO,
false, false, Interface::lowerlim);
static Switch<ClusterFissioner,int> interfaceEnhanceSProb
("EnhanceSProb",
"Option for enhancing strangeness",
&ClusterFissioner::_enhanceSProb, 0, false, false);
static SwitchOption interfaceEnhanceSProbNo
(interfaceEnhanceSProb,
"No",
"No strangeness enhancement.",
0);
static SwitchOption interfaceEnhanceSProbScaled
(interfaceEnhanceSProb,
"Scaled",
"Scaled strangeness enhancement",
1);
static SwitchOption interfaceEnhanceSProbExponential
(interfaceEnhanceSProb,
"Exponential",
"Exponential strangeness enhancement",
2);
static Switch<ClusterFissioner,int> interfaceKinematics
("Kinematics",
"Option for selecting different Kinematics for ClusterFission",
&ClusterFissioner::_kinematics, 0, false, false);
static SwitchOption interfaceKinematicsDefault
(interfaceKinematics,
"Default",
"Fully aligned Cluster Fission along the Original cluster direction",
0);
static SwitchOption interfaceKinematicsIsotropic
(interfaceKinematics,
"Isotropic",
"Fully isotropic two body decay. Not recommended!",
1);
static SwitchOption interfaceKinematicsFluxTube
(interfaceKinematics,
"FluxTube",
"Aligned decay with gaussian pT kick with sigma=ClusterFissioner::FluxTube",
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,Energy> interfaceFluxTubeWidth
("FluxTubeWidth",
"sigma of gaussian sampling of pT for FluxTube kinematics",
&ClusterFissioner::_fluxTubeWidth, GeV, 2.0*GeV, 0.0*GeV, 5.0*GeV,
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> interfacePhaseSpaceWeights
("PhaseSpaceWeights",
"Include phase space weights.",
&ClusterFissioner::_phaseSpaceWeights, false, false, false);
static SwitchOption interfacePhaseSpaceWeightsYes
(interfacePhaseSpaceWeights,
"Yes",
"Do include the effect of cluster fission phase space",
true);
static SwitchOption interfacePhaseSpaceWeightsNo
(interfacePhaseSpaceWeights,
"No",
"Do not include the effect of cluster phase space",
false);
static 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 Parameter<ClusterFissioner,double>
interfaceDim ("Dimension","Dimension in which phase space weights are calculated",
&ClusterFissioner::_dim, 0, 4.0, 3.0, 10.0,false,false, Interface::limited);
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);
}
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) {
// 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();
Energy Mc = cluster->momentum().m();
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());
// TODO 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;
// std::cout << "Make Clusters: ( " << ptrQ1->PDGName() << " " << newPtr1->PDGName() << " ), ( "
// << ptrQ2->PDGName() << " " << newPtr2->PDGName() << " )\n";
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);
}
// std::cout << "Make Clusters: ( " << ptrQ1->PDGName() << " " << newPtr1->PDGName() << " ), ( "
// << ptrQ2->PDGName() << " " << newPtr2->PDGName() << " )\n";
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
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);
// derive the masses of the children
Mc1 = res.first;
Mc2 = res.second;
// static kinematic threshold
if(_kinematicThresholdChoice == 0) {
if(Mc1 < m1+m || Mc2 < m+m2 || Mc1+Mc2 > Mc) continue;
// dynamic kinematic threshold
} else if(_kinematicThresholdChoice == 1) {
bool C1 = ( sqr(Mc1) )/( sqr(m1) + sqr(m) + _kinThresholdShift ) < 1.0 ? true : false;
bool C2 = ( sqr(Mc2) )/( sqr(m2) + sqr(m) + _kinThresholdShift ) < 1.0 ? true : false;
bool C3 = ( sqr(Mc1) + sqr(Mc2) )/( sqr(Mc) ) > 1.0 ? true : false;
if( C1 || C2 || C3 ) continue;
}
if ( _phaseSpaceWeights ) {
if ( phaseSpaceVeto(Mc,Mc1,Mc2,m,m1,m2) )
continue;
}
if ( Mc1 < spectrum()->massLightestHadronPair(ptrQ1->dataPtr(),newPtr1->dataPtr())) {
// std::cout << "Cluster decays to hadron MC"<< Mc1/GeV << "\t MLHP "<< spectrum()->massLightestHadronPair(ptrQ1->dataPtr(),newPtr1->dataPtr())/GeV<<"\n";
// std::cout << "Flavour " << ptrQ1->PDGName() <<", " << newPtr1->PDGName()<< "\n";
continue;
}
if ( Mc2 < spectrum()->massLightestHadronPair(ptrQ2->dataPtr(),newPtr2->dataPtr())) {
// std::cout << "Cluster decays to hadron MC"<< Mc2/GeV << "\t MLHP "<< spectrum()->massLightestHadronPair(ptrQ2->dataPtr(),newPtr2->dataPtr())/GeV<<"\n";
// std::cout << "Flavour " << ptrQ2->PDGName() <<", " << newPtr2->PDGName()<< "\n";
continue;
}
// if (diqClu1) {
// if (Mc1 < spectrum()->massLightestBaryonPair(ptrQ1->dataPtr(),newPtr1->dataPtr())) {
// continue;
// }
// }
// if (diqClu2) {
// if (Mc2 < spectrum()->massLightestBaryonPair(ptrQ2->dataPtr(),newPtr2->dataPtr())) {
// 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) {std::cout << "\nHadron1\n" << std::endl; Mc1 = toHadron1->mass(); pClu1.setMass(Mc1); }
toHadron2 = _hadronSpectrum->chooseSingleHadron(ptrQ2->dataPtr(), newPtr2->dataPtr(),Mc2);
if(toHadron2) {std::cout << "\nHadron2\n" << std::endl; Mc2 = toHadron2->mass(); pClu2.setMass(Mc2); }
if (toHadron1 || toHadron2) throw Exception()
<< "Cluster wants to go to H,C or H,H'" << Exception::runerror;
// 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) {std::cout << "\nHadron1\n" << std::endl; Mc1 = toHadron1->mass(); pClu1.setMass(Mc1); }
}
else if(!toHadron2) {
toHadron2 = _hadronSpectrum->chooseSingleHadron(ptrQ2->dataPtr(), newPtr2->dataPtr(),Mc-Mc1);
if(toHadron2) {std::cout << "\nHadron2\n" << std::endl; 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)
p0Q1.setMass(ptrQ1->mass()); // known (mom Q1 before fission)
p0Q1.rescaleEnergy();
pClu.rescaleMass();
calculateKinematics(pClu,p0Q1,toHadron1,toHadron2,
pClu1,pClu2,pQ1,pQone,pQtwo,pQ2);
/******************
* The previous methods have determined the kinematics and positions
* of C -> C1 + C2.
* In the case that one of the two product is light, that means either
* decayOneHadronClu1 or decayOneHadronClu2 is true, then the momenta
* of the components of that light product have not been determined,
* and a (light) cluster will not be created: the heavy father cluster
* decays, in this case, into a single (not-light) cluster and a
* single hadron. In the other, "normal", cases the father cluster
* decays into two clusters, each of which has well defined components.
* Notice that, in the case of components which point to particles, the
* momenta of the components is properly set to the new values, whereas
* we do not change the momenta of the pointed particles, because we
* want to keep all of the information (that is the new momentum of a
* component after the splitting, which is contained in the _momentum
* member of the Component class, and the (old) momentum of that component
* before the splitting, which is contained in the momentum of the
* pointed particle). Please not make confusion of this only apparent
* inconsistency!
********************/
LorentzPoint posC,pos1,pos2;
posC = cluster->vertex();
calculatePositions(pClu, posC, pClu1, pClu2, pos1, pos2);
cutType rval;
if(toHadron1) {
rval.first = produceHadron(toHadron1, newPtr1, pClu1, pos1);
finalhadrons.push_back(rval.first.first);
}
else {
rval.first = produceCluster(ptrQ1, newPtr1, pClu1, pos1, pQ1, pQone, rem1);
}
if(toHadron2) {
rval.second = produceHadron(toHadron2, newPtr2, pClu2, pos2);
finalhadrons.push_back(rval.second.first);
}
else {
rval.second = produceCluster(ptrQ2, newPtr2, pClu2, pos2, pQ2, pQtwo, rem2);
}
return rval;
}
ClusterFissioner::cutType
ClusterFissioner::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);
pair<Energy,Energy> res = drawNewMasses(mmax, soft1, soft2, pClu1, pClu2,
ptrQ[iq1], pQ1, newPtr1, pQone,
newPtr2, pQtwo, ptrQ[iq1], pQ2);
Mc1 = res.first; Mc2 = res.second;
if(Mc1 < m1+m || Mc2 < m+m2 || Mc1+Mc2 > mmax) continue;
if ( _phaseSpaceWeights ) {
if ( phaseSpaceVeto(mmax,Mc1,Mc2,m,m1,m2) )
continue;
}
// check if need to force meson clster to hadron
toHadron = _hadronSpectrum->chooseSingleHadron(ptrQ[iq1]->dataPtr(), newPtr1->dataPtr(),Mc1);
if(toHadron) { Mc1 = toHadron->mass(); pClu1.setMass(Mc1); }
// check if need to force diquark cluster to be on-shell
toDiQuark = false;
diquark = spectrum()->makeDiquark(ptrQ[iq2]->dataPtr(), newPtr2->dataPtr());
if(Mc2 < diquark->constituentMass()) {
Mc2 = diquark->constituentMass(); pClu2.setMass(Mc2);
toDiQuark = true;
}
// if a beam cluster not allowed to decay to hadrons
if(cluster->isBeamCluster() && toHadron && softUEisOn)
continue;
// Check if the decay kinematics is still possible: if not then
// force the one-hadron decay for the other cluster as well.
if(Mc1 + Mc2 > mmax) {
if(!toHadron) {
toHadron = _hadronSpectrum->chooseSingleHadron(ptrQ[iq1]->dataPtr(), newPtr1->dataPtr(),mmax-Mc2);
if(toHadron) { Mc1 = toHadron->mass(); pClu1.setMass(Mc1); }
}
else if(!toDiQuark) {
Mc2 = _hadronSpectrum->massLightestHadron(ptrQ[iq2]->dataPtr(), newPtr2->dataPtr()); pClu2.setMass(Mc2);
toDiQuark = true;
}
}
succeeded = (mmax >= Mc1+Mc2);
}
while (!succeeded && counter < max_loop);
// check no of tries
if(counter >= max_loop) return cutType();
// Determine the (5-components) momenta (all in the LAB frame)
Lorentz5Momentum p0Q1 = ptrQ[iq1]->momentum();
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;
}
void ClusterFissioner::drawNewFlavourDiquarks(PPtr& newPtrPos,PPtr& newPtrNeg,
const ClusterPtr & clu) const {
// Flavour is assumed to be only u, d, s, with weights
// (which are not normalized probabilities) given
// by the same weights as used in HadronsSelector for
// the decay of clusters into two hadrons.
unsigned hasDiquarks=0;
assert(clu->numComponents()==2);
tcPDPtr pD1=clu->particle(0)->dataPtr();
tcPDPtr pD2=clu->particle(1)->dataPtr();
bool isDiq1=DiquarkMatcher::Check(pD1->id());
if (isDiq1) hasDiquarks++;
bool isDiq2=DiquarkMatcher::Check(pD2->id());
if (isDiq2) hasDiquarks++;
assert(hasDiquarks<=2);
Energy Mc=(clu->momentum().mass());
// if (fabs(clu->momentum().massError() )>1e-14) std::cout << "Mass inconsistency CF : " << std::scientific << clu->momentum().massError() <<"\n";
// Not allow yet Diquark Clusters
// if ( hasDiquarks>=1 || Mc < spectrum()->massLightestBaryonPair(pD1,pD2) )
// return drawNewFlavour(newPtrPos,newPtrNeg);
Energy minMass;
Selector<long> choice;
// int countQ=0;
// int countDiQ=0;
// adding quark-antiquark pairs to the selection list
for ( const long& id : spectrum()->lightHadronizingQuarks() ) {
// TODO uncommenting below gives sometimes 0 selection possibility,
// maybe need to be checked in the LightClusterDecayer and ColourReconnector
// if (Mc < spectrum()->massLightestHadronPair(pD1,pD2)) continue;
// countQ++;
if (_fissionCluster==0) choice.insert(_hadronSpectrum->pwtQuark(id),id);
else if (_fissionCluster==1) choice.insert(_fissionPwt.find(id)->second,id);
else assert(false);
}
// adding diquark-antidiquark pairs to the selection list
switch (hasDiquarks)
{
case 0:
for ( const long& id : spectrum()->lightHadronizingDiquarks() ) {
if (_strictDiquarkKinematics) {
tPDPtr cand = getParticleData(id);
minMass = spectrum()->massLightestHadron(pD2,cand)
+ spectrum()->massLightestHadron(cand,pD1);
}
else minMass = spectrum()->massLightestBaryonPair(pD1,pD2);
if (Mc < minMass) continue;
// countDiQ++;
if (_fissionCluster==0) choice.insert(_hadronSpectrum->pwtQuark(id),id);
else if (_fissionCluster==1) choice.insert(_fissionPwt.find(id)->second,id);
else assert(false);
}
break;
case 1:
if (_diquarkClusterFission<1) break;
for ( const long& id : spectrum()->lightHadronizingDiquarks() ) {
tPDPtr diq = getParticleData(id);
if (isDiq1)
minMass = spectrum()->massLightestHadron(pD2,diq)
+ spectrum()->massLightestBaryonPair(diq,pD1);
else
minMass = spectrum()->massLightestHadron(pD1,diq)
+ spectrum()->massLightestBaryonPair(diq,pD2);
if (Mc < minMass) continue;
// countDiQ++;
if (_fissionCluster==0) choice.insert(_hadronSpectrum->pwtQuark(id),id);
else if (_fissionCluster==1) choice.insert(_fissionPwt.find(id)->second,id);
else assert(false);
}
break;
case 2:
if (_diquarkClusterFission<2) break;
for ( const long& id : spectrum()->lightHadronizingDiquarks() ) {
tPDPtr diq = getParticleData(id);
if (Mc < spectrum()->massLightestBaryonPair(pD1,pD2)) {
throw Exception() << "Found Diquark Cluster:\n" << *clu << "\nwith MassCluster = "
<< ounit(Mc,GeV) <<" GeV MassLightestBaryonPair = "
<< ounit(spectrum()->massLightestBaryonPair(pD1,pD2) ,GeV)
<< " GeV cannot decay" << Exception::eventerror;
}
minMass = spectrum()->massLightestBaryonPair(pD1,diq)
+ spectrum()->massLightestBaryonPair(diq,pD2);
if (Mc < minMass) continue;
// countDiQ++;
if (_fissionCluster==0) choice.insert(_hadronSpectrum->pwtQuark(id),id);
else if (_fissionCluster==1) choice.insert(_fissionPwt.find(id)->second,id);
else assert(false);
}
break;
default:
assert(false);
}
assert(choice.size()>0);
long idNew = choice.select(UseRandom::rnd());
newPtrPos = getParticle(idNew);
newPtrNeg = getParticle(-idNew);
assert (newPtrPos);
assert(newPtrNeg);
assert (newPtrPos->dataPtr());
assert(newPtrNeg->dataPtr());
}
void ClusterFissioner::drawNewFlavour(PPtr& newPtrPos,PPtr& newPtrNeg) const {
// Flavour is assumed to be only u, d, s, with weights
// (which are not normalized probabilities) given
// by the same weights as used in HadronsSelector for
// the decay of clusters into two hadrons.
Selector<long> choice;
switch(_fissionCluster){
case 0:
for ( const long& id : spectrum()->lightHadronizingQuarks() )
choice.insert(_hadronSpectrum->pwtQuark(id),id);
break;
case 1:
for ( const long& id : spectrum()->lightHadronizingQuarks() )
choice.insert(_fissionPwt.find(id)->second,id);
break;
default :
assert(false);
}
long idNew = choice.select(UseRandom::rnd());
newPtrPos = getParticle(idNew);
newPtrNeg = getParticle(-idNew);
assert (newPtrPos);
assert(newPtrNeg);
assert (newPtrPos->dataPtr());
assert(newPtrNeg->dataPtr());
}
void ClusterFissioner::drawNewFlavourEnhanced(PPtr& newPtrPos,PPtr& newPtrNeg,
Energy2 mass2) const {
if ( spectrum()->gluonId() != ParticleID::g )
throw Exception() << "strange enhancement only working with Standard Model hadronization"
<< Exception::runerror;
// Flavour is assumed to be only u, d, s, with weights
// (which are not normalized probabilities) given
// by the same weights as used in HadronsSelector for
// the decay of clusters into two hadrons.
double prob_d = 0.;
double prob_u = 0.;
double prob_s = 0.;
double scale = abs(double(sqr(_m0Fission)/mass2));
// Choose which splitting weights you wish to use
switch(_fissionCluster){
// 0: ClusterFissioner and ClusterDecayer use the same weights
case 0:
prob_d = _hadronSpectrum->pwtQuark(ParticleID::d);
prob_u = _hadronSpectrum->pwtQuark(ParticleID::u);
/* Strangeness enhancement:
Case 1: probability scaling
Case 2: Exponential scaling
*/
if (_enhanceSProb == 1)
prob_s = (_maxScale < scale) ? 0. : pow(_hadronSpectrum->pwtQuark(ParticleID::s),scale);
else if (_enhanceSProb == 2)
prob_s = (_maxScale < scale) ? 0. : exp(-scale);
break;
/* 1: ClusterFissioner uses its own unique set of weights,
i.e. decoupled from ClusterDecayer */
case 1:
prob_d = _fissionPwt.find(ParticleID::d)->second;
prob_u = _fissionPwt.find(ParticleID::u)->second;
if (_enhanceSProb == 1)
prob_s = (_maxScale < scale) ? 0. : pow(_fissionPwt.find(ParticleID::s)->second,scale);
else if (_enhanceSProb == 2)
prob_s = (_maxScale < scale) ? 0. : exp(-scale);
break;
default:
assert(false);
}
int choice = UseRandom::rnd3(prob_u, prob_d, prob_s);
long idNew = 0;
switch (choice) {
case 0: idNew = ThePEG::ParticleID::u; break;
case 1: idNew = ThePEG::ParticleID::d; break;
case 2: idNew = ThePEG::ParticleID::s; break;
}
newPtrPos = getParticle(idNew);
newPtrNeg = getParticle(-idNew);
assert (newPtrPos);
assert(newPtrNeg);
assert (newPtrPos->dataPtr());
assert(newPtrNeg->dataPtr());
}
Energy2 ClusterFissioner::clustermass(const ClusterPtr & cluster) const {
Lorentz5Momentum pIn = cluster->momentum();
Energy2 endpointmass2 = sqr(cluster->particle(0)->mass() +
cluster->particle(1)->mass());
Energy2 singletm2 = pIn.m2();
// Return either the cluster mass, or the lambda measure
return (_massMeasure == 0) ? singletm2 : singletm2 - endpointmass2;
}
Energy ClusterFissioner::drawChildMass(const Energy M, const Energy m1,
const Energy m2, const Energy m,
const double expt, const bool soft) const {
/***************************
* This method, given in input the cluster mass Mclu of an heavy cluster C,
* made of consituents of masses m1 and m2, draws the masses Mclu1 and Mclu2
* of, respectively, the children cluster C1, made of constituent masses m1
* and m, and cluster C2, of mass Mclu2 and made of constituent masses m2
* and m. The mass is extracted from one of the two following mass
* distributions:
* --- power-like ("normal" distribution)
* d(Prob) / d(M^exponent) = const
* where the exponent can be different from the two children C1 (exp1)
* and C2 (exponent2).
* --- exponential ("soft" distribution)
* d(Prob) / d(M^2) = exp(-b*M)
* where b = 2.0 / average.
* Such distributions are limited below by the masses of
* the constituents quarks, and above from the mass of decaying cluster C.
* The choice of which of the two mass distributions to use for each of the
* two cluster children is dictated by iRemnant (see below).
* If the number of attempts to extract a pair of mass values that are
* kinematically acceptable is above some fixed number (max_loop, see below)
* the method gives up and returns false; otherwise, when it succeeds, it
* returns true.
*
* These distributions have been modified from HERWIG:
* Before these were:
* Mclu1 = m1 + (Mclu - m1 - m2)*pow( rnd(), 1.0/exponent1 );
* The new one coded here is a more efficient version, same density
* but taking into account 'in phase space from' beforehand
***************************/
// hard cluster
if(!soft) {
return pow(UseRandom::rnd(pow((M-m1-m2-m)*UnitRemoval::InvE, expt),
pow(m*UnitRemoval::InvE, expt)), 1./expt
)*UnitRemoval::E + m1;
}
// Otherwise it uses a soft mass distribution
else {
static const InvEnergy b = 2.0 / _btClM;
Energy max = M-m1-m2-2.0*m;
double rmin = b*max;
rmin = ( rmin < 50 ) ? exp(-rmin) : 0.;
double r1;
do {
r1 = UseRandom::rnd(rmin, 1.0) * UseRandom::rnd(rmin, 1.0);
}
while (r1 < rmin);
return m1 + m - log(r1)/b;
}
}
Axis ClusterFissioner::sampleDirectionAligned(const Lorentz5Momentum & pRelCOM) const {
return pRelCOM.vect().unit();
}
Axis ClusterFissioner::sampleDirectionUniform() const {
double cosTheta = -1 + 2.0 * UseRandom::rnd();
double sinTheta = sqrt(1.0-cosTheta*cosTheta);
double Phi = 2.0 * Constants::pi * UseRandom::rnd();
Axis uClusterUniform(cos(Phi)*cosTheta, sin(Phi)*cosTheta, sinTheta);
return uClusterUniform.unit();
}
Axis ClusterFissioner::sampleDirectionSemiUniform(const Lorentz5Momentum & pRelCOM) const {
Axis dir=pRelCOM.vect().unit();
Axis res;
do {
res=sampleDirectionUniform();
}
while (dir*res<0);
return res;
}
Energy ClusterFissioner::sample2DGaussianPT(const Energy & Pcom) const {
Energy sigmaPT=_fluxTubeWidth*Pcom/GeV;
if (_fluxTubeWidth==ZERO) return ZERO;
Energy magnitude;
Energy pTx,pTy;
// Energy2 pT2;
const int maxcount=100;
int count=0;
do {
if (count>=maxcount) {
if (Pcom>0.5*sigmaPT) throw Exception() << "Could not sample direction in ClusterFissioner::sampleDirectionGaussianPT() "
<< Exception::eventerror;
// Fallback uniform sampling
// double phi=UseRandom::rnd()*2.0*Constants::pi;
magnitude=UseRandom::rnd()*Pcom;
// pTx = magnitude*cos(phi);
// pTy = magnitude*sin(phi);
break;
}
pTx = UseRandom::rndGauss(sigmaPT);
pTy = UseRandom::rndGauss(sigmaPT);
// pT2 = sqr(pTx)+sqr(pTy);
magnitude=sqrt(sqr(pTx)+sqr(pTy));
count++;
}
while (magnitude>Pcom);
return magnitude;
}
Axis ClusterFissioner::sampleDirectionGaussianPT(const Lorentz5Momentum & pRelCOM,
const Energy Mass, const Energy mass1, const Energy mass2) const {
Energy pstarCOM = Kinematics::pstarTwoBodyDecay(Mass,mass1,mass2);
Energy pT=sample2DGaussianPT(pstarCOM);
Energy2 pT2=pT*pT;
double phi=UseRandom::rnd()*2.0*Constants::pi;
Energy pTx=pT*cos(phi);
Energy pTy=pT*sin(phi);
// const int maxcount=100;
// int count=0;
// Energy sigmaPT=_fluxTubeWidth;
// do {
// if (count>=maxcount) {
// if (pstarCOM>0.5*sigmaPT) throw Exception() << "Could not sample direction in ClusterFissioner::sampleDirectionGaussianPT() "
// << Exception::eventerror;
// //Fallback uniform sampling
// double phi=UseRandom::rnd()*2.0*Constants::pi;
// Energy magnitude=UseRandom::rnd()*pstarCOM;
// pTx = magnitude*cos(phi);
// pTy = magnitude*sin(phi);
// pT2 = sqr(pTx)+sqr(pTy);
// break;
// }
// pTx = UseRandom::rndGauss(sigmaPT);
// pTy = UseRandom::rndGauss(sigmaPT);
// pT2 = sqr(pTx)+sqr(pTy);
// count++;
// }
// while (pT2>sqr(pstarCOM));
Axis pRelativeDir=pRelCOM.vect().unit();
Axis TrvX = pRelativeDir.orthogonal().unit();
Axis TrvY = TrvX.cross(pRelativeDir).unit();
Axis pTkick = pTx/pstarCOM*TrvX + pTy/pstarCOM*TrvY;
Axis uClusterFluxTube = (pRelativeDir*sqrt(1.0-pT2/sqr(pstarCOM)) + pTkick);
return uClusterFluxTube.unit();
}
static void printVec(const Lorentz5Momentum & p);
static void printVec(const Lorentz5Momentum & p){
std::cout << "( "
<< std::setprecision(8)
<< p.t()/GeV << ",\t"
<< p.vect().x()/GeV << ",\t"
<< p.vect().y()/GeV << ",\t"
<< p.vect().z()/GeV
<< " )\t|p| = "
<< sqrt(sqr(p.vect())/GeV2) <<"\tm = "
<< p.m()/GeV <<"\tmass = "
<< p.mass()/GeV
<< std::setprecision(3)
<<"\n";
}
// static void printVecBoost(const boost::numeric::ublas::vector<Energy> & p);
// static void printVecBoost(const boost::numeric::ublas::vector<Energy> & p){
// std::cout << "( "
// << p[0]/GeV << ",\t"
// << p[1]/GeV << ",\t"
// << p[2]/GeV << ",\t"
// << p[3]/GeV
// << " )\t|p| = "
// << sqrt(sqr(p[1]*p[1]+p[2]*p[2]+p[3]*p[3])/GeV2) <<"\tmass = "
// << sqrt(sqr(p[0]*p[0]-(p[1]*p[1]+p[2]*p[2]+p[3]*p[3]))/GeV2) <<"\n";
// }
static void printVecBoostDBL(const boost::numeric::ublas::vector<double> & p);
static void printVecBoostDBL(const boost::numeric::ublas::vector<double> & p){
std::cout << "( "
<< std::setprecision(8)
<< p[0] << ",\t"
<< p[1] << ",\t"
<< p[2] << ",\t"
<< p[3]
<< " )\t|p| = "
<< sqrt(p[1]*p[1]+p[2]*p[2]+p[3]*p[3]) <<"\tmass = "
<< sqrt(p[0]*p[0]-(p[1]*p[1]+p[2]*p[2]+p[3]*p[3]))
<< std::setprecision(3)
<<"\n";
}
static int printIfNotNullBOOST(const boost::numeric::ublas::vector<double> & p1,const boost::numeric::ublas::vector<double> & p2);
static int printIfNotNullBOOST(const boost::numeric::ublas::vector<double> & p1,const boost::numeric::ublas::vector<double> & p2)
{
double eps=1e-5;
double Eavg=(p1(0)+p2(0))/2.0;
if (fabs(p1(0)-p2(0))>eps*Eavg || std::isnan(p1(0)-p2(0)))
{
std::cout << std::setprecision(int(-log10(eps))) << "DeltaE not 0" << std::endl;
printVecBoostDBL(p1);
printVecBoostDBL(p2);
std::cout << std::setprecision(3);
return 1;
}
double P1=sqrt(pow(p1(1),2)+pow(p1(2),2)+pow(p1(3),2));
double P2=sqrt(pow(p2(1),2)+pow(p2(2),2)+pow(p2(3),2));
double Pavg=(P1+P2)/2.0;
if (fabs(p1(1)-p2(1))>eps*Pavg || std::isnan(p1(1)-p2(1)))
{
std::cout << std::setprecision(int(-log10(eps))) << "DeltaPx not 0" << std::endl;
printVecBoostDBL(p1);
printVecBoostDBL(p2);
std::cout << std::setprecision(3);
return 1;
}
if (fabs(p1(2)-p2(2))>eps*Pavg|| std::isnan(p1(2)-p2(2))
)
{
std::cout << std::setprecision(int(-log10(eps))) << "DeltaPy not 0" << std::endl;
printVecBoostDBL(p1);
printVecBoostDBL(p2);
std::cout << std::setprecision(3);
return 1;
}
if (fabs(p1(3)-p2(3))>eps*Pavg|| std::isnan(p1(3)-p2(3))
)
{
std::cout << std::setprecision(int(-log10(eps))) << "DeltaPz not 0" << std::endl;
printVecBoostDBL(p1);
printVecBoostDBL(p2);
std::cout << std::setprecision(3);
return 1;
}
return 0;
}
static int printIfNotNull(const Lorentz5Momentum & p);
static int printIfNotNull(const Lorentz5Momentum & p)
{
double eps=1e-5;
if (abs(p.t()/GeV)>eps || std::isnan(p.t()/GeV) )
{
std::cout << std::setprecision(int(-log10(eps))) << "DeltaE not 0" << std::endl;
printVec(p);
std::cout << std::setprecision(3);
return 1;
}
if (abs(p.x()/GeV)>eps|| std::isnan(p.x()/GeV) )
{
std::cout << std::setprecision(int(-log10(eps))) << "DeltaPx not 0" << std::endl;
printVec(p);
std::cout << std::setprecision(3);
return 1;
}
if (abs(p.y()/GeV)>eps|| std::isnan(p.y()/GeV) )
{
std::cout << std::setprecision(int(-log10(eps))) << "DeltaPy not 0" << std::endl;
printVec(p);
std::cout << std::setprecision(3);
return 1;
}
if (abs(p.z()/GeV)>eps|| std::isnan(p.z()/GeV) )
{
std::cout << std::setprecision(int(-log10(eps))) << "DeltaPz not 0" << std::endl;
printVec(p);
std::cout << std::setprecision(3);
return 1;
}
return 0;
}
static int printIfNotNullComp(const Lorentz5Momentum & p1,const boost::numeric::ublas::vector<double> & p2);
static int printIfNotNullComp(const Lorentz5Momentum & p1,const boost::numeric::ublas::vector<double> & p2)
{
double eps=1e-5;
double Eavg=(p1.t()/GeV+p2(0))/2.0;
if (fabs(p1.t()/GeV-p2(0))>eps*Eavg|| std::isnan(p1.t()/GeV-p2(0)) )
{
std::cout << std::setprecision(int(-log10(eps))) << "DeltaE not 0" << std::endl;
printVec(p1);
printVecBoostDBL(p2);
std::cout << std::setprecision(3);
return 1;
}
double P1=p1.vect().mag()/GeV;
double P2=sqrt(pow(p2(1),2)+pow(p2(2),2)+pow(p2(3),2));
double Pavg=(P1+P2)/2.0;
if (abs(p1.x()/GeV-p2(1))>eps*Pavg|| std::isnan(p1.x()/GeV-p2(1)) )
{
std::cout << std::setprecision(int(-log10(eps))) << "DeltaPx not 0" << std::endl;
printVec(p1);
printVecBoostDBL(p2);
std::cout << std::setprecision(3);
return 1;
}
if (abs(p1.y()/GeV-p2(2))>eps*Pavg|| std::isnan(p1.y()/GeV-p2(2)) )
{
std::cout << std::setprecision(int(-log10(eps))) << "DeltaPy not 0" << std::endl;
printVec(p1);
printVecBoostDBL(p2);
std::cout << std::setprecision(3);
return 1;
}
if (abs(p1.z()/GeV-p2(3))>eps*Pavg|| std::isnan(p1.z()/GeV-p2(3)) )
{
std::cout << std::setprecision(int(-log10(eps))) << "DeltaPz not 0" << std::endl;
printVec(p1);
printVecBoostDBL(p2);
std::cout << std::setprecision(3);
return 1;
}
return 0;
}
void ClusterFissioner::calculateKinematicsClusterCOM(const Energy M, const Lorentz5Momentum & piLab,
const Lorentz5Momentum & pjLab,
Lorentz5Momentum & pClu1,
Lorentz5Momentum & pClu2) const {
// std::cout << std::setprecision(18)<< "######### Info #########" << std::endl;
// std::cout << std::setprecision(18)<< "M = "<< M/GeV << std::endl;
// std::cout << std::setprecision(18)<< "pi+pj mass = "<< (piLab+pjLab).m()/GeV << std::endl;
// std::cout << std::setprecision(18)<< "piLab.mass = "<< piLab.mass()/GeV << std::endl;
// std::cout << std::setprecision(18)<< "pjLab.mass = "<< pjLab.mass()/GeV << std::endl;
// std::cout << std::setprecision(18)<< "piLab.m = "<< piLab.m()/GeV << std::endl;
// std::cout << std::setprecision(18)<< "pjLab.m = "<< pjLab.m()/GeV << std::endl;
// std::cout << std::setprecision(18)<< "pClu1.mass = "<< pClu1.mass()/GeV << std::endl;
// std::cout << std::setprecision(18)<< "pClu2.mass = "<< pClu2.mass()/GeV << std::endl;
// std::cout << std::setprecision(18)<< "pClu1.m = "<< pClu1.m()/GeV << std::endl;
// std::cout << std::setprecision(18)<< "pClu2.m = "<< pClu2.m()/GeV << std::endl;
// printIfNotNull(pClu1+pClu2);
// std::cout << "######### End #########" << std::endl;
// const Lorentz5Momentum piLab(p0Q1);
// const Lorentz5Momentum pjLab(p0Q2);
// piLab.setMass(p0Q1.mass());
// pjLab.setMass(p0Q2.mass());
// Energy Ei=sqrt(piLab.mass()*piLab.mass()+sqr(piLab.vect()));
// Energy Ej=sqrt(pjLab.mass()*pjLab.mass()+sqr(pjLab.vect()));
Energy Ei=piLab.e();
Energy Ej=pjLab.e();
// Final angle choice :
double cosPhi=piLab.vect().cosTheta(pjLab.vect());
double Phi=acos(cosPhi);
double sinPhi=sin(Phi);
if (std::isnan(Phi) || std::isinf(Phi)) throw Exception() << "NAN in Phi\n"
<< Exception::runerror;
Energy2 Ei2=Ei*Ei;
Energy2 Ej2=Ej*Ej;
Energy mi=piLab.mass();
Energy mj=pjLab.mass();
Energy2 mi2=mi*mi;
Energy2 mj2=mj*mj;
Energy Pi=piLab.vect().mag();
Energy Pj=pjLab.vect().mag();
- if (fabs((Ei-Pi)/GeV)<1e-5) std::cout<< std::setprecision(18) << "WARNING: Ei-Pi " << (Ei-Pi)/GeV<< std::endl;
- if (fabs((Ej-Pj)/GeV)<1e-5) std::cout<< std::setprecision(18) << "WARNING: Ej-Pj " << (Ej-Pj)/GeV<< std::endl;
+ // if (fabs((Ei-Pi)/GeV)<1e-5) std::cout<< std::setprecision(18) << "WARNING: Ei-Pi " << (Ei-Pi)/GeV<< std::endl;
+ // if (fabs((Ej-Pj)/GeV)<1e-5) std::cout<< std::setprecision(18) << "WARNING: Ej-Pj " << (Ej-Pj)/GeV<< std::endl;
boost::numeric::ublas::vector<double> Pclu1Hat(4);
boost::numeric::ublas::vector<double> Pclu2Hat(4);
// const Energy M = (piLab+pjLab).m();
const Energy M1 = pClu1.mass();
const Energy M2 = pClu2.mass();
// const Energy PcomClu=Kinematics::pstarTwoBodyDecay(M,M1,M2);
// Energy pT = sample2DGaussianPT(PcomClu);
// Energy pZ = sqrt(PcomClu*PcomClu - pT*pT);
// Pclu1Hat(0)=sqrt(M1*M1+PcomClu*PcomClu)/GeV;
// Pclu1Hat(1)=pT/GeV*cos(phi);
// Pclu1Hat(2)=pT/GeV*sin(phi);
// Pclu1Hat(3)=pZ/GeV;
// Pclu2Hat(0)=sqrt(M2*M2+PcomClu*PcomClu)/GeV;
// Pclu2Hat(1)=-pT/GeV*cos(phi);
// Pclu2Hat(2)=-pT/GeV*sin(phi);
// Pclu2Hat(3)=-pZ/GeV;
Pclu1Hat(0)=pClu1.e()/GeV;
Pclu1Hat(1)=pClu1.x()/GeV;
Pclu1Hat(2)=pClu1.y()/GeV;
Pclu1Hat(3)=pClu1.z()/GeV;
Pclu2Hat(0)=pClu2.e()/GeV;
Pclu2Hat(1)=pClu2.x()/GeV;
Pclu2Hat(2)=pClu2.y()/GeV;
Pclu2Hat(3)=pClu2.z()/GeV;
// const Energy2 pipj=piLab*pjLab;
// std::cout << "pipjDiff = " << sqrt(fabs((mi2+mj2+2*pipj-(M1*M1+M2*M2+2*pClu1*pClu2))/GeV2))<< std::endl;
// std::cout << "pipjDiff23 = " << sqrt(fabs((mi2+mj2+2*pipj-((piLab+pjLab).m2()))/GeV2))<< std::endl;
// std::cout << "pipjDiff3 = " << sqrt(fabs(((M1*M1+M2*M2+2*pClu1*pClu2)-((piLab+pjLab).m2()))/GeV2))<< std::endl;
// if (std::isnan(pipj/GeV2) || std::isinf(pipj/GeV2)) throw Exception() << "NAN in pipj/GeV2\n"
// << Exception::runerror;
boost::numeric::ublas::matrix<double> Lambda(4,4);
// Energy sqrtS=sqrt(mi2+mj2+2*pipj);
// const Energy sqrtS=(pClu1+pClu2).m();
// const Energy sqrtS=(piLab+pjLab).m();
const Energy sqrtS=M;
// TODO better numerics for pipj:
const Energy2 pipj=(sqrtS*sqrtS-mi2-mj2)/2.0;
// lambda(sqrtS,mi,mj)
// Energy4 A4=(pipj*pipj-mi2*mj2);
Energy pstar=Kinematics::pstarTwoBodyDecay(sqrtS,mi,mj);
Energy2 A2=pstar*sqrtS;
Energy2 B2=Pi*Pj*sinPhi;
Energy B2divPi=Pj*sinPhi;
+ double sinSqrPhidiv2=pow(sin(Phi/2.0),2);
if (std::isnan(A2/GeV2) || std::isinf(A2/GeV2)) throw Exception() << "NAN in A2/GeV2\n"
<< Exception::runerror;
Lambda(0,0) = (Ei+Ej)/sqrtS;
Lambda(0,1) = B2/A2;
Lambda(0,2) = 0;
Lambda(0,3) = (Ei-Ej)/(2.0*pstar)-((mi2-mj2)*(Ei+Ej))/(2.0*sqrtS*A2);
Lambda(1,0) = B2divPi/sqrtS;
- Lambda(1,1) = (Pi*Ej-Ei*Pj*cosPhi)/A2;
+ // Lambda(1,1) = (Pi*Ej-Ei*Pj*cosPhi)/A2;
+ // // TODO
+ Lambda(1,1) = Phi<1e-10? 1:(Pi*Ej-Ei*Pj*cosPhi)/A2;
Lambda(1,2) = 0;
Lambda(1,3) = -(sqrtS*sqrtS-(mj2-mi2))*B2divPi/(2.0*sqrtS*A2);
Lambda(2,0) = 0;
Lambda(2,1) = 0;
Lambda(2,2) = 1;
Lambda(2,3) = 0;
Lambda(3,0) = (Pi+Pj*cosPhi)/sqrtS;
Lambda(3,1) = Ei*B2divPi/A2;
Lambda(3,2) = 0;
// Lambda(3,3) = (mj2*Pi+Pj*cosPhi*(Pi*Pj*cosPhi-Ei2-Ei*Ej)+Pi*Ei*Ej)/(sqrtS*A2);
assert(Pi>ZERO);
// Lambda(3,3) = (A2*A2+0.5*Ei*(sqrtS*sqrtS*(Ei-Ej)-(mi2-mj2)*(Ei+Ej)))/(Pi*sqrtS*A2);
+ // // TODO
Lambda(3,3) = (A2*A2-0.5*Ei*(Ej*(sqrtS*sqrtS-(mj2-mi2))-Ei*(sqrtS*sqrtS-(mi2-mj2))))/(Pi*sqrtS*A2);
double det=0;
det += Lambda(0,0)*Lambda(1,1)*Lambda(3,3);
det += Lambda(1,0)*Lambda(3,1)*Lambda(0,3);
det += Lambda(3,0)*Lambda(0,1)*Lambda(1,3);
det -= Lambda(0,3)*Lambda(1,1)*Lambda(3,0);
det -= Lambda(1,0)*Lambda(0,1)*Lambda(3,3);
det -= Lambda(1,3)*Lambda(3,1)*Lambda(0,0);
- if (fabs(det-1.0)>1e-3 || std::isnan(det)) {
+ if (fabs(det-1.0)>1e-1 || std::isnan(det)) {
std::cout << "A2 = "<<std::setprecision(18)<< A2/(GeV2) << std::endl;
std::cout << "B2 = "<< B2/(GeV2)<< std::endl;
std::cout << "DET-1 = "<< det-1.0 << std::setprecision(3)<< std::endl;
std::cout << "Lambda:" << std::endl;
for (int i = 0; i < 4; i++)
{
for (int j = 0; j < 4; j++)
{
std::cout<< std::setprecision(18) << Lambda(i,j)<< std::setprecision(3) << "\t";
}
std::cout << "\n";
}
std::cout << "pi,pj in Lab" << std::endl;
printVec(piLab);
printVec(pjLab);
std::cout << "Phi = " << std::setprecision(18) << Phi << std::endl;
std::cout << "Phi - pi = " << std::setprecision(18) << Phi - M_PI << std::endl;
}
boost::numeric::ublas::vector<double> Pclu1=boost::numeric::ublas::prod(Lambda,Pclu1Hat);
boost::numeric::ublas::vector<double> Pclu2=boost::numeric::ublas::prod(Lambda,Pclu2Hat);
Axis zAxis(0,0,1);
Axis xAxis(1,0,0);
Lorentz5Momentum piRes(mi,Pi*zAxis);
Lorentz5Momentum pjRes(mj,Pj*(xAxis*sinPhi+zAxis*cosPhi));
Lorentz5Momentum pClu1Res(M1,GeV*Axis(Pclu1[1],Pclu1[2],Pclu1[3]));
Lorentz5Momentum pClu2Res(M2,GeV*Axis(Pclu2[1],Pclu2[2],Pclu2[3]));
Axis omega1=piRes.vect().unit().cross(piLab.vect().unit());
double cosAngle1=piRes.vect().unit()*piLab.vect().unit();
double angle1=acos(cosAngle1);
piRes.rotate(angle1, omega1);
pjRes.rotate(angle1, omega1);
pClu1Res.rotate(angle1, omega1);
pClu2Res.rotate(angle1, omega1);
Axis omega2=piRes.vect().unit();
// if (fabs(sinPhi)<1e-13)
// {
// std::cout << "sinPhi "<<std::setprecision(18)<< sinPhi<<std::setprecision(3)<<"\n";
// std::cout << "piLab\n";
// printVec(piLab);
// std::cout << "pjLab\n";
// printVec(pjLab);
// std::cout << "piRes\n";
// printVec(piRes);
// std::cout << "pjRes\n";
// printVec(pjRes);
// if (fabs(sinPhi)<1e-14) {
// std::cout << "r1 = "<<std::setprecision(18)<< (pjLab.vect()-piRes.vect().unit()*(pjLab.vect()*piRes.vect().unit())).mag()/GeV << "\n";
// std::cout << "r2 = "<< (pjRes.vect()-piRes.vect().unit()*(pjRes.vect()*piRes.vect().unit())).mag()/GeV << std::setprecision(3)<<std::endl;
// }
// }
Momentum3 r1dim=(pjLab.vect()-piRes.vect().unit()*(pjLab.vect()*piRes.vect().unit()));
Momentum3 r2dim=(pjRes.vect()-piRes.vect().unit()*(pjRes.vect()*piRes.vect().unit()));
if (r1dim.mag()==ZERO || r2dim.mag()==ZERO || fabs(sinPhi)<1e-14)
{
// std::cout << "!!!!!!!!!!!RETURNN " << std::endl;
pClu1=pClu1Res;
// pClu1.setMass(M1);
// pClu1.rescaleMass();
// pClu1.rescaleEnergy();
// if (fabs((pClu1.mass()-pClu1.m())/GeV)>1e-6*pClu1.mass()/GeV) {
// std::cout << "M1.mass == " <<std::setprecision(18)<< pClu1.mass()/GeV << std::endl;
// std::cout << "M1.m == " <<std::setprecision(18)<< pClu1.m()/GeV << std::endl;
// std::cout << "DeltaM1 == " <<std::setprecision(18)<< pClu1.mass()/GeV - pClu1.m()/GeV << std::endl;
// std::cout << "Phi = " << std::setprecision(18) << Phi << std::endl;
// std::cout << "Phi - pi = " << std::setprecision(18) << Phi - M_PI << std::endl;
// }
pClu2=pClu2Res;
// pClu2.setMass(M2);
// pClu2.rescaleMass();
// pClu2.rescaleEnergy();
//
// if (fabs((pClu2.mass()-pClu2.m())/GeV)>1e-10*pClu2.mass()/GeV) {
// std::cout << "M2.mass == " <<std::setprecision(18)<< pClu2.mass()/GeV << std::endl;
// std::cout << "M2.m == " <<std::setprecision(18)<< pClu2.m()/GeV << std::endl;
// std::cout << "DeltaM2 == " <<std::setprecision(18)<< pClu2.mass()/GeV - pClu2.m()/GeV << std::endl;
// std::cout << "Phi = " << std::setprecision(18) << Phi << std::endl;
// std::cout << "Phi - pi = " << std::setprecision(18) << Phi - M_PI << std::endl;
// }
return;
}
Axis r1=r1dim.unit();
Axis r2=r2dim.unit();
int signToPi = (piRes.vect()*pjLab.vect())/GeV2 > 0 ? 1:-1;
int signToR1R2 = signToPi*(r2.unit().cross(r1.unit())*piRes.vect())/GeV> 0 ? 1:-1;
double angle2=acos(r1.unit()*r2.unit());
if (signToR1R2<0) angle2=-angle2;
// std::cout << "angle2 = " << angle2 << "\t cosAngle2 = "<< cos(angle2) << std::endl;
// piRes.rotate(-angle2, omega2.unit());
// pjRes.rotate(-angle2, signToPi*omega2.unit());
pjRes.rotate(angle2, signToPi*omega2);
pClu1Res.rotate(angle2, signToPi*omega2);
pClu2Res.rotate(angle2, signToPi*omega2);
if (printIfNotNull(piRes-piLab) ){std::cout << "^- piRes after rot\n";
std::cout << "omega2 = ("
<< omega2.unit().x() << "\t"
<< omega2.unit().y() << "\t"
<< omega2.unit().z() << ")\n";
std::cout << "Phi = "<<std::setprecision(18)<< Phi<<std::setprecision(3)<< std::endl;
std::cout << "Phi -pi = "<<std::setprecision(18)<< Phi-M_PI<<std::setprecision(3)<<std::endl;
printVec(piRes);
printVec(piLab);
}
if (printIfNotNull(pjRes-pjLab) ){std::cout << "^- pjRes after rot\n";
std::cout << "omega2 = ("
<< omega2.unit().x() << "\t"
<< omega2.unit().y() << "\t"
<< omega2.unit().z() << ")\n";
std::cout << "Phi = "<<std::setprecision(18)<< Phi<<std::setprecision(3)<< std::endl;
std::cout << "Phi -pi = "<<std::setprecision(18)<< Phi-M_PI<<std::setprecision(3)<<std::endl;
printVec(pjRes);
printVec(pjLab);
}
if (printIfNotNull(piRes+pjRes-pClu1Res-pClu2Res) ){std::cout << "^- sanity total mom conserv\n";
std::cout << "Phi = "<<std::setprecision(18)<< Phi<<std::setprecision(3)<<std::endl;
std::cout << "Phi - pi= "<<std::setprecision(18)<< Phi-M_PI<<std::setprecision(3)<<std::endl;
printVec(piRes+pjRes);
printVec(pClu1Res+pClu2Res);
std::cout << "M = " <<std::setprecision(18)<< M/GeV << std::endl;
std::cout << "M1 = "<<std::setprecision(18) << M1/GeV << std::endl;
std::cout << "M2 = " <<std::setprecision(18)<< M2/GeV << std::endl;
}
pClu1=pClu1Res;
// pClu1.setMass(M1);
// pClu1.rescaleMass();
// pClu1.rescaleEnergy();
- if (fabs((pClu1.mass()-pClu1.m())/GeV)>1e-10*pClu1.mass()/GeV) {
- std::cout << "M1.mass == " <<std::setprecision(18)<< pClu1.mass()/GeV << std::endl;
- std::cout << "M1.m == " <<std::setprecision(18)<< pClu1.m()/GeV << std::endl;
- std::cout << "DeltaM1 == " <<std::setprecision(18)<< pClu1.mass()/GeV - pClu1.m()/GeV << std::endl;
- }
+ // if (fabs((pClu1.mass()-pClu1.m())/GeV)>1e-10*pClu1.mass()/GeV) {
+ // std::cout << "M1.mass == " <<std::setprecision(18)<< pClu1.mass()/GeV << std::endl;
+ // std::cout << "M1.m == " <<std::setprecision(18)<< pClu1.m()/GeV << std::endl;
+ // std::cout << "DeltaM1 == " <<std::setprecision(18)<< pClu1.mass()/GeV - pClu1.m()/GeV << std::endl;
+ // }
pClu2=pClu2Res;
// pClu2.setMass(M2);
// pClu2.rescaleMass();
// pClu2.rescaleEnergy();
- if (fabs((pClu2.mass()-pClu2.m())/GeV)>1e-10*pClu2.mass()/GeV) {
- std::cout << "M2.mass == " <<std::setprecision(18)<< pClu2.mass()/GeV << std::endl;
- std::cout << "M2.m == " <<std::setprecision(18)<< pClu2.m()/GeV << std::endl;
- std::cout << "DeltaM2 == " <<std::setprecision(18)<< pClu2.mass()/GeV - pClu2.m()/GeV << std::endl;
- }
+ // if (fabs((pClu2.mass()-pClu2.m())/GeV)>1e-10*pClu2.mass()/GeV) {
+ // std::cout << "M2.mass == " <<std::setprecision(18)<< pClu2.mass()/GeV << std::endl;
+ // std::cout << "M2.m == " <<std::setprecision(18)<< pClu2.m()/GeV << std::endl;
+ // std::cout << "DeltaM2 == " <<std::setprecision(18)<< pClu2.mass()/GeV - pClu2.m()/GeV << std::endl;
+ // }
// std::cout << "############################\n";
}
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
**************************/
if (pClu.mass() < pClu1.mass() + pClu2.mass()
|| pClu1.mass()<ZERO
|| pClu2.mass()<ZERO ) {
throw Exception() << "Impossible Kinematics in ClusterFissioner::calculateKinematics[FluxTube]() (A)"
<< Exception::eventerror;
}
// Calculate the unit three-vector, in the Cluster frame,
// using the sampling funtion sampleDirection in the respective
// clusters rest frame
// Calculate the momenta of C1 and C2 in the (parent) C frame first,
// where the direction of C1 is samples according to sampleDirection,
// and then boost back in the LAB frame.
//
// relative momentum in centre of mass of cluster
// OLD
Lorentz5Momentum uRelCluster(p0Q1);
uRelCluster.boost( -pClu.boostVector() ); // boost from LAB to C
// NEW ? :
// Lorentz5Momentum uRelCluster(p0Q1);
// Lorentz5Momentum pCluTmp(pClu);
// uRelCluster.boost(-p0Q1.boostVector());
// pCluTmp.boost(-p0Q1.boostVector());
// uRelCluster.boost(-pCluTmp.boostVector());
// sample direction (options = Default(aligned), Isotropic
// or FluxTube(gaussian pT kick))
Axis DirClu = sampleDirection(uRelCluster,
pClu.mass(), pClu1.mass(), pClu2.mass());
Axis DirClu1;
Axis DirClu2;
// int _covariantBoost=1;
if (_covariantBoost) {
Lorentz5Momentum p0Q2(pQ2bar.mass(),(pClu-p0Q1).vect());
if (fabs((p0Q1+p0Q2).m()/GeV- pClu.mass()/GeV)>1e-5*pClu.mass()/GeV
){
// printVec(p0Q1);
// printVec(p0Q2);
// printVec(p0Q1+p0Q2);
// printVec(pClu);
- Lorentz5Momentum pTottmp(p0Q1+p0Q2,pClu.mass());
- std::cout << "M tot = "<<std::setprecision(18)<< (pTottmp).massError() << std::endl;
- std::cout << "M.m clu = "<<std::setprecision(18)<< pClu.m()/GeV << std::endl;
- std::cout << "M.mass clu = "<< pClu.mass()/GeV << std::endl;
- std::cout << "MassError pTottmp = "<<std::setprecision(18)<< (pTottmp).massError() << std::endl;
- std::cout << "MassError pClu = "<<std::setprecision(18)<< pClu.massError() << std::endl;
- std::cout << "Mclu1 = "<<std::setprecision(18)<< pClu1.mass()/GeV << std::endl;
- std::cout << "Mclu2 = "<<std::setprecision(18)<< pClu2.mass()/GeV << std::endl;
- std::cout << "MassError pClu1 = "<<std::setprecision(18)<< pClu1.massError() << std::endl;
- std::cout << "MassError pClu2 = "<<std::setprecision(18)<< pClu2.massError() << std::endl;
+ // Lorentz5Momentum pTottmp(p0Q1+p0Q2,pClu.mass());
+ // std::cout << "M tot = "<<std::setprecision(18)<< (pTottmp).massError() << std::endl;
+ // std::cout << "M.m clu = "<<std::setprecision(18)<< pClu.m()/GeV << std::endl;
+ // std::cout << "M.mass clu = "<< pClu.mass()/GeV << std::endl;
+ // std::cout << "MassError pTottmp = "<<std::setprecision(18)<< (pTottmp).massError() << std::endl;
+ // std::cout << "MassError pClu = "<<std::setprecision(18)<< pClu.massError() << std::endl;
+ // std::cout << "Mclu1 = "<<std::setprecision(18)<< pClu1.mass()/GeV << std::endl;
+ // std::cout << "Mclu2 = "<<std::setprecision(18)<< pClu2.mass()/GeV << std::endl;
+ // std::cout << "MassError pClu1 = "<<std::setprecision(18)<< pClu1.massError() << std::endl;
+ // std::cout << "MassError pClu2 = "<<std::setprecision(18)<< pClu2.massError() << std::endl;
}
const Energy M = pClu.mass();
const Energy M1 = pClu1.mass();
const Energy M2 = pClu2.mass();
const Energy PcomClu=Kinematics::pstarTwoBodyDecay(M,M1,M2);
Energy pT = sample2DGaussianPT(PcomClu);
Energy pZ = sqrt(PcomClu*PcomClu - pT*pT);
double phi = 2*Constants::pi*UseRandom::rnd();
Momentum3 pClu1sampVect( pT*cos(phi), pT*sin(phi), pZ);
Momentum3 pClu2sampVect(-pT*cos(phi),-pT*sin(phi),-pZ);
// Lorentz5Momentum pClu1Samp(M1,pClu1sampVect);
// Lorentz5Momentum pClu2Samp(M2,pClu2sampVect);
// pClu1=pClu1Samp;
pClu1.setMass(M1);
pClu1.setVect(pClu1sampVect);
pClu1.rescaleEnergy();
// pClu2=pClu2Samp;
pClu2.setMass(M2);
pClu2.setVect(pClu2sampVect);
pClu2.rescaleEnergy();
// std::cout << std::setprecision(18)<< "######### Info2 #########" << std::endl;
// std::cout << std::setprecision(18)<< "M = "<< pClu.mass()/GeV << std::endl;
// std::cout << std::setprecision(18)<< "pi+pj mass = "<< (p0Q1+p0Q2).m()/GeV << std::endl;
// std::cout << std::setprecision(18)<< "p0Q1.mass = "<< p0Q1.mass()/GeV << std::endl;
// std::cout << std::setprecision(18)<< "p0Q2.mass = "<< p0Q2.mass()/GeV << std::endl;
// std::cout << std::setprecision(18)<< "p0Q1.m = "<< p0Q1.m()/GeV << std::endl;
// std::cout << std::setprecision(18)<< "p0Q2.m = "<< p0Q2.m()/GeV << std::endl;
// std::cout << std::setprecision(18)<< "pClu1.mass = "<< pClu1.mass()/GeV << std::endl;
// std::cout << std::setprecision(18)<< "pClu2.mass = "<< pClu2.mass()/GeV << std::endl;
// std::cout << std::setprecision(18)<< "pClu1.m = "<< pClu1.m()/GeV << std::endl;
// std::cout << std::setprecision(18)<< "pClu2.m = "<< pClu2.m()/GeV << std::endl;
// printIfNotNull(pClu1+pClu2);
// std::cout << "######### End #########" << std::endl;
calculateKinematicsClusterCOM(pClu.mass(),p0Q1, p0Q2, pClu1, pClu2);
// Kinematics::BoostIntoTwoParticleFrame(p0Q1, p0Q2, pClu1, pClu2);
// Lorentz5Momentum pClu1_tmp(pClu1);
// pClu1_tmp.rescaleMass();
// pClu1_tmp.boost( -pClu.boostVector() );
// if ((pClu1_tmp.vect().mag2())>ZERO){
// }
// else {
// std::cout << "Danger:!!!!" << std::endl;
// std::cout << "pClu1_tmp.vect().mag2() == " << pClu1_tmp.vect().mag2()/GeV2<< std::endl;
// std::cout << "M.mass clu = "<< pClu.mass()/GeV << std::endl;
// std::cout << "Mclu1 = "<<std::setprecision(18)<< pClu1.mass()/GeV << std::endl;
// std::cout << "Mclu2 = "<<std::setprecision(18)<< pClu2.mass()/GeV << std::endl;
// std::cout << "pClu.m = "<<pClu.m()/GeV << std::endl;
// std::cout << "pClu.mass = "<<pClu.mass()/GeV << std::endl;
// std::cout << "pClu.t = "<<pClu.t()/GeV << std::endl;
// std::cout << "pClu.x = "<<pClu.vect().x()/GeV << std::endl;
// std::cout << "pClu.y = "<<pClu.vect().y()/GeV << std::endl;
// std::cout << "pClu.z = "<<pClu.vect().z()/GeV << std::endl;
// std::cout << "pClu2.m = "<<pClu2.m()/GeV << std::endl;
// std::cout << "pClu2.mass = "<<pClu2.mass()/GeV << std::endl;
// std::cout << "pClu2.t = "<<pClu2.t()/GeV << std::endl;
// std::cout << "pClu2.x = "<<pClu2.vect().x()/GeV << std::endl;
// std::cout << "pClu2.y = "<<pClu2.vect().y()/GeV << std::endl;
// std::cout << "pClu2.z = "<<pClu2.vect().z()/GeV << std::endl;
// std::cout << "pClu1.m = "<<pClu1.m()/GeV << std::endl;
// std::cout << "pClu1.mass = "<<pClu1.mass()/GeV << std::endl;
// std::cout << "pClu1.t = "<<pClu1.t()/GeV << std::endl;
// std::cout << "pClu1.x = "<<pClu1.vect().x()/GeV << std::endl;
// std::cout << "pClu1.y = "<<pClu1.vect().y()/GeV << std::endl;
// std::cout << "pClu1.z = "<<pClu1.vect().z()/GeV << std::endl;
// std::cout << "pClu1_tmp.m = "<<pClu1_tmp.m()/GeV << std::endl;
// std::cout << "pClu1_tmp.mass = "<<pClu1_tmp.mass()/GeV << std::endl;
// std::cout << "pClu1_tmp.t = "<<pClu1_tmp.t()/GeV << std::endl;
// std::cout << "pClu1_tmp.x = "<<pClu1_tmp.vect().x()/GeV << std::endl;
// std::cout << "pClu1_tmp.y = "<<pClu1_tmp.vect().y()/GeV << std::endl;
// std::cout << "pClu1_tmp.z = "<<pClu1_tmp.vect().z()/GeV << std::endl;
// }
// DirClu=pClu1_tmp.vect().mag2()>ZERO ? pClu1_tmp.vect().unit():sampleDirectionUniform();
DirClu1=sampleDirectionUniform();
DirClu2=sampleDirectionUniform();
}
else {
Kinematics::twoBodyDecay(pClu, pClu1.mass(), pClu2.mass(),DirClu, pClu1, pClu2);
DirClu1=sampleDirectionUniform();
DirClu2=sampleDirectionUniform();
}
// 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[FluxTube]() (B)"
<< Exception::eventerror;
}
// sample direction (options = Default(aligned), Isotropic
// or FluxTube(gaussian pT kick))
// Axis DirClu1 = sampleDirection(uRelCluster,
// pClu1.mass(), pQ1.mass(), pQbar.mass());
Kinematics::twoBodyDecay(pClu1, pQ1.mass(), pQbar.mass(),
DirClu1, pQ1, pQbar);
}
// In the case that cluster2 does not decay immediately into a single hadron,
// Calculate the momenta of Q and Q2bar (as constituent of C2) in the
// (parent) C2 frame first, where the direction of Q is u.vect().unit(),
// and then boost back in the LAB frame.
if(!toHadron2) {
if (pClu2.m() < pQ.mass() + pQ2bar.mass() ) {
throw Exception() << "Impossible Kinematics in ClusterFissioner::calculateKinematics[FluxTube]() (C)"
<< Exception::eventerror;
}
// sample direction (options = Default(aligned), Isotropic
// or FluxTube(gaussian pT kick))
// Axis DirClu2 = sampleDirection(uRelCluster,
// pClu2.mass(), pQ2bar.mass(), pQ.mass());
Kinematics::twoBodyDecay(pClu2, pQ.mass(), pQ2bar.mass(),
DirClu2, pQ, pQ2bar);
}
}
void ClusterFissioner::calculatePositions(const Lorentz5Momentum & pClu,
const LorentzPoint & positionClu,
const Lorentz5Momentum & pClu1,
const Lorentz5Momentum & pClu2,
LorentzPoint & positionClu1,
LorentzPoint & positionClu2) const {
// Determine positions of cluster children.
// See Marc Smith's thesis, page 127, formulas (4.122) and (4.123).
Energy Mclu = pClu.m();
Energy Mclu1 = pClu1.m();
Energy Mclu2 = pClu2.m();
// Calculate the unit three-vector, in the C frame, along which
// children clusters move.
Lorentz5Momentum u(pClu1);
u.boost( -pClu.boostVector() ); // boost from LAB to C frame
// the unit three-vector is then u.vect().unit()
Energy pstarChild = Kinematics::pstarTwoBodyDecay(Mclu,Mclu1,Mclu2);
// First, determine the relative positions of the children clusters
// in the parent cluster reference frame.
Energy2 mag2 = u.vect().mag2();
InvEnergy fact = mag2>ZERO ? 1./sqrt(mag2) : 1./GeV;
Length x1 = ( 0.25*Mclu + 0.5*( pstarChild + (sqr(Mclu2) - sqr(Mclu1))/(2.0*Mclu)))/_kappa;
Length t1 = Mclu/_kappa - x1;
LorentzDistance distanceClu1( x1 * fact * u.vect(), t1 );
Length x2 = (-0.25*Mclu + 0.5*(-pstarChild + (sqr(Mclu2) - sqr(Mclu1))/(2.0*Mclu)))/_kappa;
Length t2 = Mclu/_kappa + x2;
LorentzDistance distanceClu2( x2 * fact * u.vect(), t2 );
// Then, transform such relative positions from the parent cluster
// reference frame to the Lab frame.
distanceClu1.boost( pClu.boostVector() );
distanceClu2.boost( pClu.boostVector() );
// Finally, determine the absolute positions in the Lab frame.
positionClu1 = positionClu + distanceClu1;
positionClu2 = positionClu + distanceClu2;
}
bool ClusterFissioner::ProbablityFunction(double scale, double threshold) {
double cut = UseRandom::rnd(0.0,1.0);
return 1./(1.+pow(abs((threshold-_probShift)/scale),_probPowFactor)) > cut ? true : false;
}
bool ClusterFissioner::isHeavy(tcClusterPtr clu) {
// particle data for constituents
tcPDPtr cptr[3]={tcPDPtr(),tcPDPtr(),tcPDPtr()};
for(size_t ix=0;ix<min(clu->numComponents(),3);++ix) {
cptr[ix]=clu->particle(ix)->dataPtr();
}
// different parameters for exotic, bottom and charm clusters
double clpow = !spectrum()->isExotic(cptr[0],cptr[1],cptr[1]) ? _clPowLight : _clPowExotic;
Energy clmax = !spectrum()->isExotic(cptr[0],cptr[1],cptr[1]) ? _clMaxLight : _clMaxExotic;
for ( const long& id : spectrum()->heavyHadronizingQuarks() ) {
if ( spectrum()->hasHeavy(id,cptr[0],cptr[1],cptr[1]) ) {
clpow = _clPowHeavy[id];
clmax = _clMaxHeavy[id];
}
}
// required test for SUSY clusters, since aboveCutoff alone
// cannot guarantee (Mc > m1 + m2 + 2*m) in cut()
static const Energy minmass
= getParticleData(ParticleID::d)->constituentMass();
bool aboveCutoff = false, canSplitMinimally = false;
// static kinematic threshold
if(_kinematicThresholdChoice == 0) {
aboveCutoff = (
pow(clu->mass()*UnitRemoval::InvE , clpow)
>
pow(clmax*UnitRemoval::InvE, clpow)
+ pow(clu->sumConstituentMasses()*UnitRemoval::InvE, clpow)
);
canSplitMinimally = clu->mass() > clu->sumConstituentMasses() + 2.0 * minmass;
}
// dynamic kinematic threshold
else if(_kinematicThresholdChoice == 1) {
//some smooth probablity function to create a dynamic thershold
double scale = pow(clu->mass()/GeV , clpow);
double threshold = pow(clmax/GeV, clpow)
+ pow(clu->sumConstituentMasses()/GeV, clpow);
aboveCutoff = ProbablityFunction(scale,threshold);
scale = clu->mass()/GeV;
threshold = clu->sumConstituentMasses()/GeV + 2.0 * minmass/GeV;
canSplitMinimally = ProbablityFunction(scale,threshold);
}
return aboveCutoff && canSplitMinimally;
}
diff --git a/Hadronization/ClusterHadronizationHandler.cc b/Hadronization/ClusterHadronizationHandler.cc
--- a/Hadronization/ClusterHadronizationHandler.cc
+++ b/Hadronization/ClusterHadronizationHandler.cc
@@ -1,594 +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 {
return "unknown interaction " + interaction;
}
_partonSplitters[interactionId] = _partonSplitter;
_clusterFinders[interactionId] = _clusterFinder;
_colourReconnectors[interactionId] = _colourReconnector;
_clusterFissioners[interactionId] = _clusterFissioner;
_lightClusterDecayers[interactionId] = _lightClusterDecayer;
_clusterDecayers[interactionId] = _clusterDecayer;
_gluonMassGenerators[interactionId] = _gluonMassGenerator;
return "";
}
namespace {
void extractChildren(tPPtr p, set<PPtr> & all) {
if (p->children().empty()) return;
for (PVector::const_iterator child = p->children().begin();
child != p->children().end(); ++child) {
all.insert(*child);
extractChildren(*child, all);
}
}
}
void ClusterHadronizationHandler::
handle(EventHandler & ch, const tPVector & tagged,
const Hint &) {
useMe();
currentHandler_ = this;
for (auto iter : _clusterFinders) {
iter.second->spectrum()->setGenerator(currentHandler_->generator());
}
PVector theList(tagged.begin(),tagged.end());
// set the scale for coloured particles to just above the gluon mass squared
// if less than this so they are classed as perturbative
// TODO: Should this be hard-coded here, or defined in inputs?
map<int,int> interactions = {{PDT::ColouredQCD, ParticleID::g}};
// PVector currentlist(tagged.begin(),tagged.end());
map<int,PVector> currentlists;
for ( const auto & p : theList ) {
for ( auto i : interactions ) {
if ( p->data().colouredInteraction() == i.first ) {
currentlists[i.first].push_back(p);
Energy2 Q02 = 1.01*sqr(getParticleData(i.second)->constituentMass());
if(currentlists[i.first].back()->scale()<Q02) currentlists[i.first].back()->scale(Q02);
}
}
}
map<int,ClusterVector> clusters;
// ATTENTION this needs to have changes to include the new hadron spectrum functionality (gluon mass generator)
for ( auto i : interactions ) {
// reshuffle to the constituents
if ( reshuffle_ ) {
vector<PVector> reshufflelists;
if ( reshuffle_ == 1 ) { // global reshuffling
reshufflelists.push_back(currentlists[i.first]);
}
else if ( reshuffle_ == 2 ) {// colour connected reshuffling
splitIntoColourSinglets(currentlists[i.first], reshufflelists, i.first);
}
Ptr<GluonMassGenerator>::ptr gMassGenerator;
auto git = _gluonMassGenerators.find(i.first);
if ( git != _gluonMassGenerators.end() )
gMassGenerator = git->second;
for (auto currentlist : reshufflelists){
// get available energy and energy needed for constituent mass shells
LorentzMomentum totalQ;
Energy needQ = ZERO;
size_t nGluons = 0; // number of gluons for which a mass need be generated
for ( auto p : currentlist ) {
totalQ += p->momentum();
if ( p->id() == i.second && gMassGenerator ) {
++nGluons;
needQ += gMassGenerator->minGluonMass();
continue;
}
needQ += p->dataPtr()->constituentMass();
}
Energy Q = totalQ.m();
if ( needQ > Q )
throw Exception() << "cannot reshuffle to constituent mass shells" << Exception::eventerror;
// generate gluon masses if needed
list<Energy> gluonMasses;
if ( nGluons && gMassGenerator )
gluonMasses = gMassGenerator->generateMany(nGluons,Q-needQ);
// set masses for inidividual particles
vector<Energy> masses;
for ( auto p : currentlist ) {
if ( p->id() == i.second && gMassGenerator ) {
list<Energy>::const_iterator it = gluonMasses.begin();
advance(it,UseRandom::irnd(gluonMasses.size()));
masses.push_back(*it);
gluonMasses.erase(it);
}
else {
masses.push_back(p->dataPtr()->constituentMass());
}
}
// reshuffle to new masses
reshuffle(currentlist,masses);
}
}
// split the gluons
if ( !currentlists[i.first].empty() ) {
assert(_partonSplitters.find(i.first) != _partonSplitters.end());
_partonSplitters[i.first]->split(currentlists[i.first]);
}
// form the clusters
if ( !currentlists[i.first].empty() ) {
assert(_clusterFinders.find(i.first) != _clusterFinders.end());
clusters[i.first] = _clusterFinders[i.first]->formClusters(currentlists[i.first]);
// reduce BV clusters to two components now if needed
if(_reduceToTwoComponents)
_clusterFinders[i.first]->reduceToTwoComponents(clusters[i.first]);
}
}
// perform colour reconnection if needed and then
// decay the clusters into one hadron
for ( auto i : interactions ) {
if ( clusters[i.first].empty() )
continue;
bool lightOK = false;
short tried = 0;
const ClusterVector savedclusters = clusters[i.first];
tPVector finalHadrons; // only needed for partonic decayer
while (!lightOK && tried++ < 10) {
// no colour reconnection with baryon-number-violating (BV) clusters
ClusterVector CRclusters, BVclusters;
CRclusters.reserve( clusters[i.first].size() );
BVclusters.reserve( clusters[i.first].size() );
for (size_t ic = 0; ic < clusters[i.first].size(); ++ic) {
ClusterPtr cl = clusters[i.first].at(ic);
bool hasClusterParent = false;
for (unsigned int ix=0; ix < cl->parents().size(); ++ix) {
if (cl->parents()[ix]->id() == ParticleID::Cluster) {
hasClusterParent = true;
break;
}
}
if (cl->numComponents() > 2 || hasClusterParent) BVclusters.push_back(cl);
else CRclusters.push_back(cl);
}
// colour reconnection
assert(_colourReconnectors.find(i.first) != _colourReconnectors.end());
_colourReconnectors[i.first]->rearrange(CRclusters);
// tag new clusters as children of the partons to hadronize
_setChildren(CRclusters);
// forms diquarks
// we should already have used _clusterFinder[i.first]
_clusterFinders[i.first]->reduceToTwoComponents(CRclusters);
// recombine vectors of (possibly) reconnected and BV clusters
clusters[i.first].clear();
clusters[i.first].insert( clusters[i.first].end(), CRclusters.begin(), CRclusters.end() );
clusters[i.first].insert( clusters[i.first].end(), BVclusters.begin(), BVclusters.end() );
// fission of heavy clusters
// NB: during cluster fission, light hadrons might be produced straight away
assert(_clusterFissioners.find(i.first) != _clusterFissioners.end());
- if (clusters[i.first].size()<=2) std::cout << "Diffractive event maybe" << std::endl;
finalHadrons = _clusterFissioners[i.first]->fission(clusters[i.first],i.first == PDT::ColouredQCD ? isSoftUnderlyingEventON() : false);
- if (clusters[i.first].size()<=2) std::cout << "Diffractive event finished" << std::endl;
// 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) {
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);
}
}