diff --git a/Hadronization/HadronSelector.cc b/Hadronization/HadronSelector.cc
--- a/Hadronization/HadronSelector.cc
+++ b/Hadronization/HadronSelector.cc
@@ -1,1030 +1,1134 @@
// -*- C++ -*-
//
// HadronSelector.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 HadronSelector class.
//
#include "HadronSelector.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include <ThePEG/PDT/EnumParticles.h>
#include <ThePEG/Repository/EventGenerator.h>
#include <ThePEG/Repository/CurrentGenerator.h>
#include <ThePEG/Repository/Repository.h>
#include "CheckId.h"
#include <ThePEG/Utilities/DescribeClass.h>
+#include "Herwig/Utilities/Kinematics.h"
using namespace Herwig;
DescribeAbstractClass<HadronSelector,Interfaced>
describeHadronSelector("Herwig::HadronSelector","Herwig.so");
namespace {
// debug helper
void dumpTable(const HadronTable & tbl) {
typedef HadronTable::const_iterator TableIter;
for (TableIter it = tbl.begin(); it != tbl.end(); ++it) {
cerr << it->first.first << ' '
<< it->first.second << '\n';
for (KupcoData::const_iterator jt = it->second.begin();
jt != it->second.end(); ++jt) {
cerr << '\t' << *jt << '\n';
}
}
}
bool weightIsLess (pair<long,double> a, pair<long,double> b) {
return a.second < b.second;
}
}
HadronSelector::HadronSelector(unsigned int opt)
: _weight1S0(Nmax,1.),_weight3S1(Nmax,1.),_weight1P1(Nmax,1.),_weight3P0(Nmax,1.),
_weight3P1(Nmax,1.),_weight3P2(Nmax,1.),_weight1D2(Nmax,1.),_weight3D1(Nmax,1.),
_weight3D2(Nmax,1.),_weight3D3(Nmax,1.),
_repwt(Lmax,vector<vector<double> >(Jmax,vector<double>(Nmax))),
_topt(opt),_trial(0),
_limBottom(), _limCharm(), _limExotic(), belowThreshold_(0)
{
// The mixing angles
// the ideal mixing angle
const double idealAngleMix = atan( sqrt(0.5) ) * 180.0 / Constants::pi;
// \eta-\eta' mixing angle
_etamix = -23.0;
// phi-omega mixing angle
_phimix = +36.0;
// h_1'-h_1 mixing angle
_h1mix = idealAngleMix;
// f_0(1710)-f_0(1370) mixing angle
_f0mix = idealAngleMix;
// f_1(1420)-f_1(1285)\f$ mixing angle
_f1mix = idealAngleMix;
// f'_2-f_2\f$ mixing angle
_f2mix = +26.0;
// eta_2(1870)-eta_2(1645) mixing angle
_eta2mix = idealAngleMix;
// phi(???)-omega(1650) mixing angle
_omhmix = idealAngleMix;
// phi_3-omega_3 mixing angle
_ph3mix = +28.0;
// eta(1475)-eta(1295) mixing angle
_eta2Smix = idealAngleMix;
// phi(1680)-omega(1420) mixing angle
_phi2Smix = idealAngleMix;
}
void HadronSelector::persistentOutput(PersistentOStream & os) const {
os << _partons
<< _etamix << _phimix << _h1mix << _f0mix << _f1mix << _f2mix
<< _eta2mix << _omhmix << _ph3mix << _eta2Smix << _phi2Smix
<< _weight1S0 << _weight3S1 << _weight1P1 << _weight3P0 << _weight3P1
<< _weight3P2 << _weight1D2 << _weight3D1 << _weight3D2 << _weight3D3
<< _forbidden << _repwt << _pwt
<< _limBottom << _limCharm << _limExotic << belowThreshold_
<< _table;
}
void HadronSelector::persistentInput(PersistentIStream & is, int) {
is >> _partons
>> _etamix >> _phimix >> _h1mix >> _f0mix >> _f1mix >> _f2mix
>> _eta2mix >> _omhmix >> _ph3mix >> _eta2Smix >> _phi2Smix
>> _weight1S0 >> _weight3S1 >> _weight1P1 >> _weight3P0 >> _weight3P1
>> _weight3P2 >> _weight1D2 >> _weight3D1 >> _weight3D2 >> _weight3D3
>> _forbidden >> _repwt >> _pwt
>> _limBottom >> _limCharm >> _limExotic >> belowThreshold_
>> _table;
}
void HadronSelector::Init() {
static ClassDocumentation<HadronSelector> documentation
("There is no documentation for the HadronSelector class");
static RefVector<HadronSelector,ParticleData> interfacePartons
("Partons",
"The partons which are to be considered as the consistuents of the hadrons.",
&HadronSelector::_partons, -1, false, false, true, false, false);
static RefVector<HadronSelector,ParticleData> interfaceForbidden
("Forbidden",
"The PDG codes of the particles which cannot be produced in the hadronization.",
&HadronSelector::_forbidden, -1, false, false, true, false, false);
//
// mixing angles
//
// the ideal mixing angle
const double idealAngleMix = atan( sqrt(0.5) ) * 180.0 / Constants::pi;
static Parameter<HadronSelector,double> interface11S0Mixing
("11S0Mixing",
"The mixing angle for the I=0 mesons from the 1 1S0 multiplet,"
" i.e. eta and etaprime.",
&HadronSelector::_etamix, -23., -180., 180.,
false, false, Interface::limited);
static Parameter<HadronSelector,double> interface13S1Mixing
("13S1Mixing",
"The mixing angle for the I=0 mesons from the 1 3S1 multiplet,"
" i.e. phi and omega.",
&HadronSelector::_phimix, +36., -180., 180.,
false, false, Interface::limited);
static Parameter<HadronSelector,double> interface11P1Mixing
("11P1Mixing",
"The mixing angle for the I=0 mesons from the 1 1P1 multiplet,"
" i.e. h_1' and h_1.",
&HadronSelector::_h1mix, idealAngleMix, -180., 180.,
false, false, Interface::limited);
static Parameter<HadronSelector,double> interface13P0Mixing
("13P0Mixing",
"The mixing angle for the I=0 mesons from the 1 3P0 multiplet,"
" i.e. f_0(1710) and f_0(1370).",
&HadronSelector::_f0mix, idealAngleMix, -180., 180.,
false, false, Interface::limited);
static Parameter<HadronSelector,double> interface13P1Mixing
("13P1Mixing",
"The mixing angle for the I=0 mesons from the 1 3P1 multiplet,"
" i.e. f_1(1420) and f_1(1285).",
&HadronSelector::_f1mix, idealAngleMix, -180., 180.,
false, false, Interface::limited);
static Parameter<HadronSelector,double> interface13P2Mixing
("13P2Mixing",
"The mixing angle for the I=0 mesons from the 1 3P2 multiplet,"
" i.e. f'_2 and f_2.",
&HadronSelector::_f2mix, 26.0, -180., 180.,
false, false, Interface::limited);
static Parameter<HadronSelector,double> interface11D2Mixing
("11D2Mixing",
"The mixing angle for the I=0 mesons from the 1 1D2 multiplet,"
" i.e. eta_2(1870) and eta_2(1645).",
&HadronSelector::_eta2mix, idealAngleMix, -180., 180.,
false, false, Interface::limited);
static Parameter<HadronSelector,double> interface13D0Mixing
("13D0Mixing",
"The mixing angle for the I=0 mesons from the 1 3D0 multiplet,"
" i.e. eta_2(1870) phi(?) and omega(1650).",
&HadronSelector::_omhmix, idealAngleMix, -180., 180.,
false, false, Interface::limited);
static Parameter<HadronSelector,double> interface13D1Mixing
("13D1Mixing",
"The mixing angle for the I=0 mesons from the 1 3D1 multiplet,"
" i.e. phi_3 and omega_3.",
&HadronSelector::_ph3mix, 28.0, -180., 180.,
false, false, Interface::limited);
static Parameter<HadronSelector,double> interface21S0Mixing
("21S0Mixing",
"The mixing angle for the I=0 mesons from the 2 1S0 multiplet,"
" i.e. eta(1475) and eta(1295).",
&HadronSelector::_eta2Smix, idealAngleMix, -180., 180.,
false, false, Interface::limited);
static Parameter<HadronSelector,double> interface23S1Mixing
("23S1Mixing",
"The mixing angle for the I=0 mesons from the 1 3S1 multiplet,"
" i.e. phi(1680) and omega(1420).",
&HadronSelector::_phi2Smix, idealAngleMix, -180., 180.,
false, false, Interface::limited);
//
// the meson weights
//
static ParVector<HadronSelector,double> interface1S0Weights
("1S0Weights",
"The weights for the 1S0 multiplets start with n=1.",
&HadronSelector::_weight1S0, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<HadronSelector,double> interface3S1Weights
("3S1Weights",
"The weights for the 3S1 multiplets start with n=1.",
&HadronSelector::_weight3S1, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<HadronSelector,double> interface1P1Weights
("1P1Weights",
"The weights for the 1P1 multiplets start with n=1.",
&HadronSelector::_weight1P1, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<HadronSelector,double> interface3P0Weights
("3P0Weights",
"The weights for the 3P0 multiplets start with n=1.",
&HadronSelector::_weight3P0, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<HadronSelector,double> interface3P1Weights
("3P1Weights",
"The weights for the 3P1 multiplets start with n=1.",
&HadronSelector::_weight3P1, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<HadronSelector,double> interface3P2Weights
("3P2Weights",
"The weights for the 3P2 multiplets start with n=1.",
&HadronSelector::_weight3P2, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<HadronSelector,double> interface1D2Weights
("1D2Weights",
"The weights for the 1D2 multiplets start with n=1.",
&HadronSelector::_weight1D2, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<HadronSelector,double> interface3D1Weights
("3D1Weights",
"The weights for the 3D1 multiplets start with n=1.",
&HadronSelector::_weight3D1, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<HadronSelector,double> interface3D2Weights
("3D2Weights",
"The weights for the 3D2 multiplets start with n=1.",
&HadronSelector::_weight3D2, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<HadronSelector,double> interface3D3Weights
("3D3Weights",
"The weights for the 3D3 multiplets start with n=1.",
&HadronSelector::_weight3D3, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static Switch<HadronSelector,unsigned int> interfaceTrial
("Trial",
"A Debugging option to only produce certain types of hadrons",
&HadronSelector::_trial, 0, false, false);
static SwitchOption interfaceTrialAll
(interfaceTrial,
"All",
"Produce all the hadrons",
0);
static SwitchOption interfaceTrialPions
(interfaceTrial,
"Pions",
"Only produce pions",
1);
static SwitchOption interfaceTrialSpin2
(interfaceTrial,
"Spin2",
"Only mesons with spin less than or equal to two are produced",
2);
static SwitchOption interfaceTrialSpin3
(interfaceTrial,
"Spin3",
"Only hadrons with spin less than or equal to three are produced",
3);
static Parameter<HadronSelector,double>
interfaceSingleHadronLimitBottom ("SingleHadronLimitBottom",
"Threshold for one-hadron decay of b-cluster",
&HadronSelector::_limBottom,
0, 0.0, 0.0, 100.0,false,false,false);
static Parameter<HadronSelector,double>
interfaceSingleHadronLimitCharm ("SingleHadronLimitCharm",
"threshold for one-hadron decay of c-cluster",
&HadronSelector::_limCharm,
0, 0.0, 0.0, 100.0,false,false,false);
static Parameter<HadronSelector,double>
interfaceSingleHadronLimitExotic ("SingleHadronLimitExotic",
"threshold for one-hadron decay of exotic cluster",
&HadronSelector::_limExotic,
0, 0.0, 0.0, 100.0,false,false,false);
static Switch<HadronSelector,unsigned int> interfaceBelowThreshold
("BelowThreshold",
"Option fo the selection of the hadrons if the cluster is below the pair threshold",
&HadronSelector::belowThreshold_, 0, false, false);
static SwitchOption interfaceBelowThresholdLightest
(interfaceBelowThreshold,
"Lightest",
"Force cluster to decay to the lightest hadron with the appropriate flavours",
0);
static SwitchOption interfaceBelowThresholdAll
(interfaceBelowThreshold,
"All",
"Select from all the hadrons below the two hadron threshold according to their spin weights",
1);
}
double HadronSelector::mixingStateWeight(long id) const {
switch(id) {
case ParticleID::eta: return 0.5*probabilityMixing(_etamix ,1);
case ParticleID::etaprime: return 0.5*probabilityMixing(_etamix ,2);
case ParticleID::phi: return 0.5*probabilityMixing(_phimix ,1);
case ParticleID::omega: return 0.5*probabilityMixing(_phimix ,2);
case ParticleID::hprime_1: return 0.5*probabilityMixing(_h1mix ,1);
case ParticleID::h_1: return 0.5*probabilityMixing(_h1mix ,2);
case 10331: return 0.5*probabilityMixing(_f0mix ,1);
case 10221: return 0.5*probabilityMixing(_f0mix ,2);
case ParticleID::fprime_1: return 0.5*probabilityMixing(_f1mix ,1);
case ParticleID::f_1: return 0.5*probabilityMixing(_f1mix ,2);
case ParticleID::fprime_2: return 0.5*probabilityMixing(_f2mix ,1);
case ParticleID::f_2: return 0.5*probabilityMixing(_f2mix ,2);
case 10335: return 0.5*probabilityMixing(_eta2mix ,1);
case 10225: return 0.5*probabilityMixing(_eta2mix ,2);
// missing phi member of 13D1 should be here
case 30223: return 0.5*probabilityMixing(_omhmix ,2);
case 337: return 0.5*probabilityMixing(_ph3mix ,1);
case 227: return 0.5*probabilityMixing(_ph3mix ,2);
case 100331: return 0.5*probabilityMixing(_eta2mix ,1);
case 100221: return 0.5*probabilityMixing(_eta2mix ,2);
case 100333: return 0.5*probabilityMixing(_phi2Smix,1);
case 100223: return 0.5*probabilityMixing(_phi2Smix,2);
default: return 1./3.;
}
}
void HadronSelector::doinit() {
Interfaced::doinit();
// set the weights for the various excited mesons
// set all to one to start with
for (int l = 0; l < Lmax; ++l ) {
for (int j = 0; j < Jmax; ++j) {
for (int n = 0; n < Nmax; ++n) {
_repwt[l][j][n] = 1.0;
}
}
}
// set the others from the relevant vectors
for( int ix=0;ix<max(int(_weight1S0.size()),int(Nmax));++ix)
_repwt[0][0][ix]=_weight1S0[ix];
for( int ix=0;ix<max(int(_weight3S1.size()),int(Nmax));++ix)
_repwt[0][1][ix]=_weight3S1[ix];
for( int ix=0;ix<max(int(_weight1P1.size()),int(Nmax));++ix)
_repwt[1][1][ix]=_weight1P1[ix];
for( int ix=0;ix<max(int(_weight3P0.size()),int(Nmax));++ix)
_repwt[1][0][ix]=_weight3P0[ix];
for( int ix=0;ix<max(int(_weight3P1.size()),int(Nmax));++ix)
_repwt[1][1][ix]=_weight3P1[ix];
for( int ix=0;ix<max(int(_weight3P2.size()),int(Nmax));++ix)
_repwt[1][2][ix]=_weight3P2[ix];
for( int ix=0;ix<max(int(_weight1D2.size()),int(Nmax));++ix)
_repwt[2][2][ix]=_weight1D2[ix];
for( int ix=0;ix<max(int(_weight3D1.size()),int(Nmax));++ix)
_repwt[2][1][ix]=_weight3D1[ix];
for( int ix=0;ix<max(int(_weight3D2.size()),int(Nmax));++ix)
_repwt[2][2][ix]=_weight3D2[ix];
for( int ix=0;ix<max(int(_weight3D3.size()),int(Nmax));++ix)
_repwt[2][3][ix]=_weight3D3[ix];
// find the maximum
map<long,double>::iterator pit =
max_element(_pwt.begin(),_pwt.end(),weightIsLess);
const double pmax = pit->second;
for(pit=_pwt.begin(); pit!=_pwt.end(); ++pit) {
pit->second/=pmax;
}
// construct the hadron tables
constructHadronTable();
// for debugging
if(Debug::level >= 10 )
dumpTable(table());
}
void HadronSelector::insertToHadronTable(tPDPtr &particle, int flav1, int flav2) {
// inserting a new Hadron in the hadron table.
long pid = particle->id();
int pspin = particle->iSpin();
double maxdd(0.),maxss(0.),maxrest(0.);
HadronInfo a(pid, particle,specialWeight(pid),particle->mass());
// set the weight to the number of spin states
a.overallWeight = pspin;
// identical light flavours
if(flav1 == flav2 && flav1<=3) {
// ddbar> uubar> admixture states
if(flav1==1) {
if(_topt != 0) a.overallWeight *= 0.5*a.swtef;
_table[make_pair(1,1)].insert(a);
_table[make_pair(2,2)].insert(a);
if(_topt == 0 && a.overallWeight > maxdd) maxdd = a.overallWeight;
}
// load up ssbar> uubar> ddbar> admixture states
else {
a.wt = mixingStateWeight(pid);
a.overallWeight *= a.wt;
if(_topt != 0) a.overallWeight *= a.swtef;
_table[make_pair(1,1)].insert(a);
_table[make_pair(2,2)].insert(a);
if(_topt == 0 && a.overallWeight > maxdd) maxdd = a.overallWeight;
a.wt = (_topt != 0) ? 1.- 2.*a.wt : 1 - a.wt;
if(a.wt > 0) {
a.overallWeight = a.wt * a.swtef * pspin;
_table[make_pair(3,3)].insert(a);
if(_topt == 0 && a.overallWeight > maxss) maxss = a.overallWeight;
}
}
}
// spin-1/2 baryons
else if(pspin==2) {
// octet bit
assert(DiquarkMatcher::Check(flav1));
long iq1 = flav1/1000;
long iq2 = (flav1/100)%10;
if(_topt != 0) a.overallWeight *= a.swtef;
if(iq1==iq2) {
assert(flav2!=iq1);
// first the uu1 d type piece
a.wt = 1./3.;
a.overallWeight *= a.wt;
_table[make_pair(flav1,flav2)].insert(a);
_table[make_pair(flav2,flav1)].insert(a);
if(_topt == 0 && a.overallWeight > maxrest) maxrest = a.overallWeight;
// also need ud1 u type
long f3 = CheckId::makeDiquarkID(iq1,flav2,3);
a.overallWeight /= a.wt;
a.wt = 1./6.;
a.overallWeight *= a.wt;
_table[make_pair(iq1,f3 )].insert(a);
_table[make_pair(f3 ,iq1)].insert(a);
if(_topt == 0 && a.overallWeight > maxrest) maxrest = a.overallWeight;
// and ud0 u type
f3 = CheckId::makeDiquarkID(iq1,flav2,1);
a.overallWeight /= a.wt;
a.wt = 0.5;
a.overallWeight *= a.wt;
_table[make_pair(iq1,f3 )].insert(a);
_table[make_pair(f3 ,iq1)].insert(a);
if(_topt == 0 && a.overallWeight > maxrest) maxrest = a.overallWeight;
}
else if(iq1==flav2) {
// ud1 u type
a.wt = 1./6.;
a.overallWeight *= a.wt;
_table[make_pair(flav1,flav2)].insert(a);
_table[make_pair(flav2,flav1)].insert(a);
if(_topt == 0 && a.overallWeight > maxrest) maxrest = a.overallWeight;
// also need ud0 u type
long f3 = CheckId::makeDiquarkID(iq1,iq2,1);
a.overallWeight /= a.wt;
a.wt = 0.5;
a.overallWeight *= a.wt;
_table[make_pair(f3 ,flav2)].insert(a);
_table[make_pair(flav2 ,f3 )].insert(a);
if(_topt == 0 && a.overallWeight > maxrest) maxrest = a.overallWeight;
// and uu1 d type
f3 = CheckId::makeDiquarkID(iq1,iq1,3);
a.overallWeight /= a.wt;
a.wt = 1./3.;
a.overallWeight *= a.wt;
_table[make_pair(f3 ,iq2)].insert(a);
_table[make_pair(iq2, f3)].insert(a);
if(_topt == 0 && a.overallWeight > maxrest) maxrest = a.overallWeight;
}
else if(iq2==flav2) assert(false);
else {
// determine if light quarks in spin 0 or spin 1
long it1 = (pid/100)%10;
long it2 = (pid/10 )%10;
// first perm
double wgt0(1./4.),wgt1(1./12.);
// only spin-0
if(it1<it2) {
long f3 = CheckId::makeDiquarkID(iq1,iq2,1);
a.wt = 1./3.;
a.overallWeight *= a.wt;
_table[make_pair(f3,flav2)].insert(a);
_table[make_pair(flav2,f3)].insert(a);
if(_topt == 0 && a.overallWeight > maxrest) maxrest = a.overallWeight;
swap(wgt0,wgt1);
}
// only spin-1
else {
a.wt = 1./3.;
a.overallWeight *= a.wt;
_table[make_pair(flav1,flav2)].insert(a);
_table[make_pair(flav2,flav1)].insert(a);
if(_topt == 0 && a.overallWeight > maxrest) maxrest = a.overallWeight;
}
// not sure here
// spin 0
a.overallWeight /= a.wt;
a.wt = wgt0;
a.overallWeight *= a.wt;
// second perm
long f3 = CheckId::makeDiquarkID(iq1,flav2,1);
_table[make_pair(iq2,f3)].insert(a);
_table[make_pair(f3,iq2)].insert(a);
if(_topt == 0 && a.overallWeight > maxrest) maxrest = a.overallWeight;
// 3rd perm
f3 = CheckId::makeDiquarkID(iq2,flav2,1);
_table[make_pair(iq1,f3)].insert(a);
_table[make_pair(f3,iq1)].insert(a);
if(_topt == 0 && a.overallWeight > maxrest) maxrest = a.overallWeight;
// spin 1
a.overallWeight /= a.wt;
a.wt = wgt1;
a.overallWeight *= a.wt;
// second perm
f3 = CheckId::makeDiquarkID(iq1,flav2,3);
_table[make_pair(iq2,f3)].insert(a);
_table[make_pair(f3,iq2)].insert(a);
if(_topt == 0 && a.overallWeight > maxrest) maxrest = a.overallWeight;
// 3rd perm
f3 = CheckId::makeDiquarkID(iq2,flav2,3);
_table[make_pair(iq1,f3)].insert(a);
_table[make_pair(f3,iq1)].insert(a);
if(_topt == 0 && a.overallWeight > maxrest) maxrest = a.overallWeight;
}
}
// spin -3/2 baryons
else if(pspin==4) {
long iq1 = flav1/1000;
long iq2 = (flav1/100)%10;
if(_topt != 0) a.overallWeight *= a.swtef;
// all the same
if(iq1==iq2 && iq1==flav2) {
a.wt = 1.;
a.overallWeight *= a.wt;
_table[make_pair(flav1,flav2)].insert(a);
_table[make_pair(flav2,flav1)].insert(a);
if(_topt == 0 && a.overallWeight > maxrest) maxrest = a.overallWeight;
}
else if(iq1==iq2) {
// first option uu1 d
a.wt = 1./3.;
a.overallWeight *= a.wt;
_table[make_pair(flav1,flav2)].insert(a);
_table[make_pair(flav2,flav1)].insert(a);
if(_topt == 0 && a.overallWeight > maxrest) maxrest = a.overallWeight;
// also need ud1 u type
long f3 = CheckId::makeDiquarkID(iq1,flav2,3);
a.overallWeight /= a.wt;
a.wt = 2./3.;
a.overallWeight *= a.wt;
_table[make_pair(iq1,f3 )].insert(a);
_table[make_pair(f3 ,iq1)].insert(a);
if(_topt == 0 && a.overallWeight > maxrest) maxrest = a.overallWeight;
}
else if(iq1==flav2) {
// also need ud1 u type
a.wt = 2./3.;
a.overallWeight *= a.wt;
_table[make_pair(flav1,flav2)].insert(a);
_table[make_pair(flav2,flav1)].insert(a);
if(_topt == 0 && a.overallWeight > maxrest) maxrest = a.overallWeight;
// and uu1 d type
long f3 = CheckId::makeDiquarkID(iq1,iq1,3);
a.overallWeight /= a.wt;
a.wt = 1./3.;
a.overallWeight *= a.wt;
_table[make_pair(iq2,f3)].insert(a);
_table[make_pair(f3,iq2)].insert(a);
if(_topt == 0 && a.overallWeight > maxrest) maxrest = a.overallWeight;
}
else if(iq2==flav2) assert(false);
else {
// just need the three different combinations
// first perm
a.wt = 1./3.;
a.overallWeight *= a.wt;
_table[make_pair(flav1,flav2)].insert(a);
_table[make_pair(flav2,flav1)].insert(a);
if(_topt == 0 && a.overallWeight > maxrest) maxrest = a.overallWeight;
// 2nd perm
long f3 = CheckId::makeDiquarkID(iq1,flav2,3);
_table[make_pair(iq2,f3)].insert(a);
_table[make_pair(f3,iq2)].insert(a);
if(_topt == 0 && a.overallWeight > maxrest) maxrest = a.overallWeight;
// 3rd perm
f3 = CheckId::makeDiquarkID(iq2,flav2,3);
_table[make_pair(iq1,f3)].insert(a);
_table[make_pair(f3,iq1)].insert(a);
if(_topt == 0 && a.overallWeight > maxrest) maxrest = a.overallWeight;
}
}
// all other cases
else {
if(_topt != 0) a.overallWeight *=a.swtef;
_table[make_pair(flav1,flav2)].insert(a);
if(flav1 != flav2) _table[make_pair(flav2,flav1)].insert(a);
if(_topt == 0 && a.overallWeight > maxrest) maxrest = a.overallWeight;
}
// normalise weights to one for first option
if(_topt == 0) {
HadronTable::const_iterator tit;
KupcoData::iterator it;
for(tit=_table.begin();tit!=_table.end();++tit) {
double weight;
if(tit->first.first==tit->first.second) {
if(tit->first.first==1||tit->first.first==2) weight=1./maxdd;
else if (tit->first.first==3) weight=1./maxss;
else weight=1./maxrest;
}
else weight=1./maxrest;
for(it = tit->second.begin(); it!=tit->second.end(); ++it) {
it->rescale(weight);
}
}
}
}
void HadronSelector::constructHadronTable() {
// initialise the table
_table.clear();
for(unsigned int ix=0; ix<_partons.size(); ++ix) {
for(unsigned int iy=0; iy<_partons.size(); ++iy) {
if (!(DiquarkMatcher::Check(_partons[ix]->id())
&& DiquarkMatcher::Check(_partons[iy]->id())))
_table[make_pair(_partons[ix]->id(),_partons[iy]->id())] = KupcoData();
}
}
// get the particles from the event generator
ParticleMap particles = generator()->particles();
// loop over the particles
for(ParticleMap::iterator it=particles.begin();
it!=particles.end(); ++it) {
long pid = it->first;
tPDPtr particle = it->second;
int pspin = particle->iSpin();
// Don't include hadrons which are explicitly forbidden
if(find(_forbidden.begin(),_forbidden.end(),particle)!=_forbidden.end())
continue;
// Don't include non-hadrons or antiparticles
if(pid < 100) continue;
// remove diffractive particles
if(pspin == 0) continue;
// K_0S and K_0L not made make K0 and Kbar0
if(pid==ParticleID::K_S0||pid==ParticleID::K_L0) continue;
// Debugging options
// Only include those with 2J+1 less than...5
if(_trial==2 && pspin >= 5) continue;
// Only include those with 2J+1 less than...7
if(_trial==3 && pspin >= 7) continue;
// Only include pions
if(_trial==1 && pid!=111 && pid!=211) continue;
// shouldn't be coloured
if(particle->coloured()) continue;
// Get the flavours
const int x4 = (pid/1000)%10;
const int x3 = (pid/100 )%10;
const int x2 = (pid/10 )%10;
const int x7 = (pid/1000000)%10;
const bool wantSusy = x7 == 1 || x7 == 2;
// Skip non-hadrons (susy particles, etc...)
if(x3 == 0 || x2 == 0) continue;
int flav1,flav2;
// meson
if(x4 == 0) {
flav1 = x2;
flav2 = x3;
if (wantSusy) flav2 += 1000000 * x7;
}
// baryon
else {
flav2 = x4;
if (wantSusy) flav2 += 1000000 * x7;
// insert the spin 1 diquark, sort out the rest later
flav1 = CheckId::makeDiquarkID(x2,x3,3);
}
insertToHadronTable(particle,flav1,flav2);
}
}
double HadronSelector::mesonWeight(long id) const {
// Total angular momentum
int j = ((id % 10) - 1) / 2;
// related to Orbital angular momentum l
int nl = (id/10000 )%10;
int l = -999;
int n = (id/100000)%10; // Radial excitation
if(j == 0) l = nl;
else if(nl == 0) l = j - 1;
else if(nl == 1 || nl == 2) l = j;
else if(nl == 3) l = j + 1;
// Angular or Radial excited meson
if((l||j||n) && l>=0 && l<Lmax && j<Jmax && n<Nmax) {
return sqr(_repwt[l][j][n]);
}
// rest is not excited or
// has spin >= 5/2 (ispin >= 6), haven't got those
else
return 1.0;
}
int HadronSelector::signHadron(tcPDPtr idQ1, tcPDPtr idQ2,
tcPDPtr hadron) const {
// This method receives in input three PDG ids, whose the
// first two have proper signs (corresponding to particles, id > 0,
// or antiparticles, id < 0 ), whereas the third one must
// be always positive (particle not antiparticle),
// corresponding to:
// --- quark-antiquark, or antiquark-quark, or
// quark-diquark, or diquark-quark, or
// antiquark-antidiquark, or antidiquark-antiquark
// for the first two input (idQ1, idQ2);
// --- meson or baryon for the third input (idHad):
// The method returns:
// --- + 1 if the two partons (idQ1, idQ2) are exactly
// the constituents for the hadron idHad;
// --- - 1 if the two partons (idQ1, idQ2) are exactly
// the constituents for the anti-hadron -idHad;
// --- + 0 otherwise.
// The method it is therefore useful to decide the
// sign of the id of the produced hadron as appeared
// in the vector _vecHad (where only hadron idHad > 0 are present)
// given the two constituent partons.
int sign = 0;
long idHad = hadron->id();
assert(idHad > 0);
int chargeIn = idQ1->iCharge() + idQ2->iCharge();
int chargeOut = hadron->iCharge();
// same charge
if( chargeIn == chargeOut && chargeIn !=0 ) sign = +1;
else if(chargeIn == -chargeOut && chargeIn !=0 ) sign = -1;
else if(chargeIn == 0 && chargeOut == 0 ) {
// In the case of same null charge, there are four cases:
// i) K0-like mesons, B0-like mesons, Bs-like mesons
// the PDG convention is to consider them "antiparticle" (idHad < 0)
// if the "dominant" (heavier) flavour (respectively, s, b)
// is a quark (idQ > 0): for instance, B0s = (b, sbar) has id < 0
// Remember that there is an important exception for K0L (id=130) and
// K0S (id=310): they don't have antiparticles, therefore idHad > 0
// always. We use below the fact that K0L and K0S are the unique
// hadrons having 0 the first (less significant) digit of their id.
// 2) D0-like mesons: the PDG convention is to consider them "particle"
// (idHad > 0) if the charm flavour is carried by a c: (c,ubar) has id>0
// 3) the remaining mesons should not have antiparticle, therefore their
// sign is always positive.
// 4) for baryons, that is when one of idQ1 and idQ2 is a (anti-) quark and
// the other one is a (anti-) diquark the sign is negative when both
// constituents are "anti", that is both with id < 0; positive otherwise.
// meson
if(abs(int(idQ1->iColour()))== 3 && abs(int(idQ2->iColour())) == 3 &&
!DiquarkMatcher::Check(idQ1->id()) && !DiquarkMatcher::Check(idQ2->id()))
{
int idQa = abs(idQ1->id());
int idQb = abs(idQ2->id());
int dominant = idQ2->id();
if(idQa > idQb) {
swap(idQa,idQb);
dominant = idQ1->id();
}
if((idQa==ParticleID::d && idQb==ParticleID::s) ||
(idQa==ParticleID::d && idQb==ParticleID::b) ||
(idQa==ParticleID::s && idQb==ParticleID::b)) {
// idHad%10 is zero for K0L,K0S
if (dominant < 0 || idHad%10 == 0) sign = +1;
else if(dominant > 0) sign = -1;
}
else if((idQa==ParticleID::u && idQb==ParticleID::c) ||
(idQa==ParticleID::u && idQb==ParticleID::t) ||
(idQa==ParticleID::c && idQb==ParticleID::t)) {
if (dominant > 0) sign = +1;
else if(dominant < 0) sign = -1;
}
else if(idQa==idQb) sign = +1;
// sets sign for Susy particles
else sign = (dominant > 0) ? +1 : -1;
}
// baryon
else if(DiquarkMatcher::Check(idQ1->id()) || DiquarkMatcher::Check(idQ2->id())) {
if (idQ1->id() > 0 && idQ2->id() > 0) sign = +1;
else if(idQ1->id() < 0 && idQ2->id() < 0) sign = -1;
}
}
if (sign == 0) {
cerr << "Could not work out sign for "
<< idQ1->PDGName() << ' '
<< idQ2->PDGName() << " => "
<< hadron->PDGName() << '\n';
assert(false);
}
return sign;
}
pair<tcPDPtr,tcPDPtr> HadronSelector::lightestHadronPair(tcPDPtr ptr1, tcPDPtr ptr2,
tcPDPtr ptr3) const {
// throw exception of id3!=0 as doesn't work
if ( ptr3 ) throw Exception()
<< "ptr3!=0 not yet implemented in HadronSelector::lightestHadronPair"
<< Exception::abortnow;
// charge
int totalcharge = ptr1->iCharge() + ptr2->iCharge();
if ( ptr3 ) totalcharge += ptr3->iCharge();
tcPDPtr vIdHad1[2]={tcPDPtr(),tcPDPtr()},vIdHad2[2]={tcPDPtr(),tcPDPtr()};
bool vOk[2] = {false, false};
Energy vMassPair[2] = { ZERO, ZERO };
for (int i = 0; i < 2; i++) {
tcPDPtr idPartner = i==0 ? getParticleData(ParticleID::d) : getParticleData(ParticleID::u);
// Change sign to idPartner (transform it into a anti-quark) if it is not
// possible to form a meson or a baryon.
assert (ptr1 && idPartner);
if (!CheckId::canBeHadron(ptr1, idPartner)) idPartner = idPartner->CC();
vIdHad1[i] = lightestHadron(ptr1, idPartner);
vIdHad2[i] = lightestHadron(ptr2, idPartner->CC());
if ( vIdHad1[i] && vIdHad2[i] &&
vIdHad1[i]->iCharge() + vIdHad2[i]->iCharge() == totalcharge ) {
vOk[i] = true;
vMassPair[i] = vIdHad1[i]->mass() + vIdHad2[i]->mass();
}
}
// Take the lightest pair compatible with charge conservation.
if ( vOk[0] && ( ! vOk[1] || vMassPair[0] <= vMassPair[1] ) ) {
return make_pair(vIdHad1[0],vIdHad2[0]);
}
else if ( vOk[1] && ( ! vOk[0] || vMassPair[1] < vMassPair[0] ) ) {
return make_pair(vIdHad1[1],vIdHad2[1]);
}
else {
return make_pair(tcPDPtr(),tcPDPtr());
}
}
Energy HadronSelector::massLightestBaryonPair(tcPDPtr ptr1, tcPDPtr ptr2) const {
// Make sure that we don't have any diquarks as input, return arbitrarily
// large value if we do
Energy currentSum = Constants::MaxEnergy;
for(unsigned int ix=0; ix<_partons.size(); ++ix) {
if(!DiquarkMatcher::Check(_partons[ix]->id())) continue;
HadronTable::const_iterator
tit1=_table.find(make_pair(abs(ptr1->id()),_partons[ix]->id())),
tit2=_table.find(make_pair(_partons[ix]->id(),abs(ptr2->id())));
if( tit1==_table.end() || tit2==_table.end()) continue;
if(tit1->second.empty()||tit2->second.empty()) continue;
Energy s = tit1->second.begin()->mass + tit2->second.begin()->mass;
if(currentSum > s) currentSum = s;
}
return currentSum;
}
tcPDPtr HadronSelector::lightestHadron(tcPDPtr ptr1, tcPDPtr ptr2,
#ifndef NDEBUG
tcPDPtr ptr3) const {
#else
tcPDPtr ) const {
#endif
// The method assumes ptr3 == 0 rest not implemented
assert(ptr1 && ptr2 && !ptr3);
// find entry in the table
pair<long,long> ids = make_pair(abs(ptr1->id()),abs(ptr2->id()));
HadronTable::const_iterator tit=_table.find(ids);
// throw exception if flavours wrong
if (tit==_table.end())
throw Exception() << "Could not find "
<< ids.first << ' ' << ids.second
<< " in _table. "
<< "In HadronSelector::lightestHadron()"
<< Exception::eventerror;
if(tit->second.empty())
throw Exception() << "HadronSelector::lightestHadron "
<< "could not find any hadrons containing "
<< ptr1->id() << ' ' << ptr2->id() << '\n'
<< tit->first.first << ' '
<< tit->first.second << Exception::eventerror;
// find the lightest hadron
int sign = signHadron(ptr1,ptr2,tit->second.begin()->ptrData);
tcPDPtr candidate = sign > 0 ?
tit->second.begin()->ptrData : tit->second.begin()->ptrData->CC();
// \todo 20 GeV limit is temporary fudge to let SM particles go through.
// \todo Use isExotic instead?
if (candidate->mass() > 20*GeV
&& candidate->mass() < ptr1->constituentMass() + ptr2->constituentMass()) {
generator()->log() << "HadronSelector::lightestHadron: "
<< "chosen candidate " << candidate->PDGName()
<< " is lighter than its constituents "
<< ptr1->PDGName() << ", " << ptr2->PDGName() << '\n'
<< candidate->mass()/GeV << " < " << ptr1->constituentMass()/GeV
<< " + " << ptr2->constituentMass()/GeV << '\n'
<< "Check your particle data tables.\n";
assert(false);
}
return candidate;
}
vector<pair<tcPDPtr,double> >
HadronSelector::hadronsBelowThreshold(Energy threshold, tcPDPtr ptr1,
tcPDPtr ptr2,
#ifndef NDEBUG
tcPDPtr ptr3) const {
#else
tcPDPtr ) const {
#endif
// The method assumes ptr3 == 0 rest not implemented
assert(ptr1 && ptr2 && !ptr3);
// find entry in the table
pair<long,long> ids = make_pair(abs(ptr1->id()),abs(ptr2->id()));
HadronTable::const_iterator tit=_table.find(ids);
// throw exception if flavours wrong
if (tit==_table.end())
throw Exception() << "Could not find "
<< ids.first << ' ' << ids.second
<< " in _table. "
<< "In HadronSelector::hadronsBelowThreshold()"
<< Exception::eventerror;
if(tit->second.empty())
throw Exception() << "HadronSelector::hadronsBelowThreshold() "
<< "could not find any hadrons containing "
<< ptr1->id() << ' ' << ptr2->id() << '\n'
<< tit->first.first << ' '
<< tit->first.second << Exception::eventerror;
vector<pair<tcPDPtr,double> > candidates;
KupcoData::const_iterator hit = tit->second.begin();
// find the hadrons
while(hit!=tit->second.end()&&hit->mass<threshold) {
// find the hadron
int sign = signHadron(ptr1,ptr2,hit->ptrData);
tcPDPtr candidate = sign > 0 ? hit->ptrData : hit->ptrData->CC();
// \todo 20 GeV limit is temporary fudge to let SM particles go through.
// \todo Use isExotic instead?
if (candidate->mass() > 20*GeV
&& candidate->mass() < ptr1->constituentMass() + ptr2->constituentMass()) {
generator()->log() << "HadronSelector::hadronsBelowTheshold: "
<< "chosen candidate " << candidate->PDGName()
<< " is lighter than its constituents "
<< ptr1->PDGName() << ", " << ptr2->PDGName() << '\n'
<< candidate->mass()/GeV << " < " << ptr1->constituentMass()/GeV
<< " + " << ptr2->constituentMass()/GeV << '\n'
<< "Check your particle data tables.\n";
assert(false);
}
candidates.push_back(make_pair(candidate,hit->overallWeight));
++hit;
}
return candidates;
}
tcPDPtr HadronSelector::chooseSingleHadron(tcPDPtr par1, tcPDPtr par2,
Energy mass) const {
// Determine the sum of the nominal masses of the two lightest hadrons
// with the right flavour numbers as the cluster under consideration.
// Notice that we don't need real masses (drawn by a Breit-Wigner
// distribution) because the lightest pair of hadrons does not involve
// any broad resonance.
Energy threshold = massLightestHadronPair(par1,par2);
// Special: it allows one-hadron decays also above threshold.
if (CheckId::isExotic(par1,par2))
threshold *= (1.0 + UseRandom::rnd()*_limExotic);
else if (CheckId::hasBottom(par1,par2))
threshold *= (1.0 + UseRandom::rnd()*_limBottom);
else if (CheckId::hasCharm(par1,par2))
threshold *= (1.0 + UseRandom::rnd()*_limCharm);
// only do one hadron decay is mass less than the threshold
if(mass>=threshold) return tcPDPtr();
// select the hadron
tcPDPtr hadron;
// old option pick the lightest hadron
if(belowThreshold_ == 0) {
hadron= lightestHadron(par1,par2);
}
// new option select from those available
else if(belowThreshold_ == 1) {
vector<pair<tcPDPtr,double> > hadrons =
hadronsBelowThreshold(threshold,par1,par2);
if(hadrons.size()==1) {
hadron = hadrons[0].first;
}
else if(hadrons.empty()) {
hadron= lightestHadron(par1,par2);
}
else {
double totalWeight=0.;
for(unsigned int ix=0;ix<hadrons.size();++ix) {
totalWeight += hadrons[ix].second;
}
totalWeight *= UseRandom::rnd();
for(unsigned int ix=0;ix<hadrons.size();++ix) {
if(totalWeight<=hadrons[ix].second) {
hadron = hadrons[ix].first;
break;
}
else
totalWeight -= hadrons[ix].second;
}
assert(hadron);
}
}
else
assert(false);
return hadron;
}
+
+pair<tcPDPtr,tcPDPtr> HadronSelector::chooseHadronPair(const Energy cluMass,
+ tcPDPtr par1, tcPDPtr par2) const {
+ useMe();
+ // if either of the input partons is a diquark don't allow diquarks to be
+ // produced
+ bool diquark = !(DiquarkMatcher::Check(par1->id()) || DiquarkMatcher::Check(par2->id()));
+ bool quark = true;
+ // decide is baryon or meson production
+ if(diquark) std::tie(quark,diquark) = selectBaryon(cluMass,par1,par2);
+ // weights for the different possibilities
+ Energy weight, wgtsum(ZERO);
+ // loop over all hadron pairs with the allowed flavours
+ static vector<Kupco> hadrons;
+ hadrons.clear();
+ for(unsigned int ix=0;ix<partons().size();++ix) {
+ tcPDPtr quarktopick = partons()[ix];
+ if(!quark && abs(int(quarktopick->iColour())) == 3
+ && !DiquarkMatcher::Check(quarktopick->id())) continue;
+ if(!diquark && abs(int(quarktopick->iColour())) == 3
+ && DiquarkMatcher::Check(quarktopick->id())) continue;
+ HadronTable::const_iterator
+ tit1 = table().find(make_pair(abs(par1->id()),quarktopick->id()));
+ HadronTable::const_iterator
+ tit2 = table().find(make_pair(quarktopick->id(),abs(par2->id())));
+ // If not in table skip
+ if(tit1 == table().end()||tit2==table().end()) continue;
+ // tables empty skip
+ const KupcoData & T1 = tit1->second;
+ const KupcoData & T2 = tit2->second;
+ if(T1.empty()||T2.empty()) continue;
+ // if too massive skip
+ if(cluMass <= T1.begin()->mass +
+ T2.begin()->mass) continue;
+ // quark weight
+ double quarkWeight = pwt(quarktopick->id());
+ // special for strange
+ if(abs(quarktopick->id()) == 3)
+ quarkWeight = strangeWeight(cluMass,par1,par2);
+ // loop over the hadrons
+ KupcoData::const_iterator H1,H2;
+ for(H1 = T1.begin();H1 != T1.end(); ++H1) {
+ for(H2 = T2.begin();H2 != T2.end(); ++H2) {
+ // break if cluster too light
+ if(cluMass < H1->mass + H2->mass) break;
+ weight = quarkWeight * H1->overallWeight * H2->overallWeight *
+ Kinematics::pstarTwoBodyDecay(cluMass, H1->mass, H2->mass);
+ int signQ = 0;
+ assert (par1 && quarktopick);
+ assert (par2);
+
+ assert(quarktopick->CC());
+
+ if(CheckId::canBeHadron(par1, quarktopick->CC())
+ && CheckId::canBeHadron(quarktopick, par2))
+ signQ = +1;
+ else if(CheckId::canBeHadron(par1, quarktopick)
+ && CheckId::canBeHadron(quarktopick->CC(), par2))
+ signQ = -1;
+ else {
+ cerr << "Could not make sign for" << par1->id()<< " " << quarktopick->id()
+ << " " << par2->id() << "\n";
+ assert(false);
+ }
+
+ if (signQ == -1)
+ quarktopick = quarktopick->CC();
+ // construct the object with the info
+ Kupco a(quarktopick, H1->ptrData, H2->ptrData, weight);
+ hadrons.push_back(a);
+ wgtsum += weight;
+ }
+ }
+ }
+ if (hadrons.empty())
+ return make_pair(tcPDPtr(),tcPDPtr());
+ // select the hadron
+ wgtsum *= UseRandom::rnd();
+ unsigned int ix=0;
+ do {
+ wgtsum-= hadrons[ix].weight;
+ ++ix;
+ }
+ while(wgtsum > ZERO && ix < hadrons.size());
+ if(ix == hadrons.size() && wgtsum > ZERO)
+ return make_pair(tcPDPtr(),tcPDPtr());
+ --ix;
+ assert(hadrons[ix].idQ);
+ int signHad1 = signHadron(par1, hadrons[ix].idQ->CC(), hadrons[ix].hadron1);
+ int signHad2 = signHadron(par2, hadrons[ix].idQ, hadrons[ix].hadron2);
+ assert( signHad1 != 0 && signHad2 != 0 );
+ return make_pair
+ ( signHad1 > 0 ? hadrons[ix].hadron1 : tcPDPtr(hadrons[ix].hadron1->CC()),
+ signHad2 > 0 ? hadrons[ix].hadron2 : tcPDPtr(hadrons[ix].hadron2->CC()));
+}
+
+pair<bool,bool> HadronSelector::selectBaryon(const Energy, tcPDPtr, tcPDPtr ) const {
+ assert(false);
+}
+
+double HadronSelector::strangeWeight(const Energy, tcPDPtr, tcPDPtr) const {
+ assert(false);
+}
diff --git a/Hadronization/HadronSelector.h b/Hadronization/HadronSelector.h
--- a/Hadronization/HadronSelector.h
+++ b/Hadronization/HadronSelector.h
@@ -1,592 +1,603 @@
// -*- C++ -*-
//
// HadronSelector.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_HadronSelector_H
#define HERWIG_HadronSelector_H
//
// This is the declaration of the HadronSelector class.
//
#include "ThePEG/Interface/Interfaced.h"
#include <ThePEG/Persistency/PersistentOStream.h>
#include <ThePEG/Persistency/PersistentIStream.h>
#include <ThePEG/PDT/ParticleData.h>
#include <ThePEG/PDT/StandardMatchers.h>
#include <ThePEG/Repository/EventGenerator.h>
#include "HadronSelector.fh"
#include "HadronInfo.h"
#include "Kupco.h"
namespace Herwig {
using namespace ThePEG;
/**\ingroup Hadronization
* \class HadronSelector
* \brief This class selects the hadron flavours of a cluster decay.
* \author Philip Stephens
* \author Alberto Ribon
* \author Peter Richardson
*
* This is the base class for the selection of either a pair of hadrons, or
* in some cases a single hadron. The different approaches which were
* previously implemented in this class are now implemented in the
* HwppSelector and Hw64Selector which inherit from this class.
*
* This class implements a number of methods which are needed by all models
* and in addition contains the weights for the different meson multiplets and
* mixing of the light \f$I=0\f$ mesons.
*
* @see \ref HadronSelectorInterfaces "The interfaces"
* defined for HadronSelector.
* @see HwppSelector
* @see Hw64Selector
*/
class HadronSelector: public Interfaced {
public:
/**
* The default constructor.
*/
HadronSelector(unsigned int);
/**
* Method to return a pair of hadrons given the PDG codes of
* two or three constituents
* @param cluMass The mass of the cluster
* @param par1 The first constituent
* @param par2 The second constituent
* @param par3 The third constituent
*/
- virtual pair<tcPDPtr,tcPDPtr> chooseHadronPair(const Energy cluMass, tcPDPtr par1,
- tcPDPtr par2,tcPDPtr par3 = PDPtr()) const = 0;
+ virtual pair<tcPDPtr,tcPDPtr> chooseHadronPair(const Energy cluMass,
+ tcPDPtr par1, tcPDPtr par2) const;
/**
* Select the single hadron for a cluster decay
* return null pointer if not a single hadron decay
* @param par1 1st constituent
* @param par2 2nd constituent
* @param mass Mass of the cluster
*/
tcPDPtr chooseSingleHadron(tcPDPtr par1, tcPDPtr par2, Energy mass) const;
/**
* This returns the lightest pair of hadrons given by the flavours.
*
* Given the two (or three) constituents of a cluster, it returns
* the two lightest hadrons with proper flavour numbers.
* Furthermore, the first of the two hadrons must have the constituent with
* par1, and the second must have the constituent with par2.
* \todo At the moment it does *nothing* in the case that also par3 is present.
*
* The method is implemented by calling twice lightestHadron,
* once with (par1,quarktopick->CC()) ,and once with (par2,quarktopick)
* where quarktopick is either the pointer to
* d or u quarks . In fact, the idea is that whatever the flavour of par1
* and par2, no matter if (anti-)quark or (anti-)diquark, the lightest
* pair of hadrons containing flavour par1 and par2 will have either
* flavour d or u, being the lightest quarks.
* The method returns the pair (PDPtr(),PDPtr()) if anything goes wrong.
*
* \todo The method assumes par3 == PDPtr() (otherwise we don't know how to proceed: a
* possible, trivial way would be to randomly select two of the three
* (anti-)quarks and treat them as a (anti-)diquark, reducing the problem
* to two components as treated below.
* In the normal (two components) situation, the strategy is the following:
* treat in the same way the two possibilities: (d dbar) (i=0) and
* (u ubar) (i=1) as the pair quark-antiquark necessary to form a
* pair of hadrons containing the input flavour par1 and par2; finally,
* select the one that produces the lightest pair of hadrons, compatible
* with the charge conservation constraint.
*/
pair<tcPDPtr,tcPDPtr> lightestHadronPair(tcPDPtr ptr1, tcPDPtr ptr2,
tcPDPtr ptr3 = PDPtr ()) const;
/**
* Returns the mass of the lightest pair of hadrons with the given particles
* @param ptr1 is the first constituent
* @param ptr2 is the second constituent
* @param ptr3 is the third constituent
*/
Energy massLightestHadronPair(tcPDPtr ptr1, tcPDPtr ptr2,
tcPDPtr ptr3 = PDPtr ()) const {
pair<tcPDPtr,tcPDPtr> pairData = lightestHadronPair(ptr1, ptr2, ptr3);
if ( ! pairData.first || ! pairData.second ) return ZERO;
return ( pairData.first->mass() + pairData.second->mass() );
}
/**
* Returns the lightest hadron formed by the given particles.
*
* Given the id of two (or three) constituents of a cluster, it returns
* the lightest hadron with proper flavour numbers.
* At the moment it does *nothing* in the case that also 'ptr3' present.
* @param ptr1 is the first constituent
* @param ptr2 is the second constituent
* @param ptr3 is the third constituent
*/
tcPDPtr lightestHadron(tcPDPtr ptr1, tcPDPtr ptr2,
tcPDPtr ptr3 = PDPtr ()) const;
/**
* Returns the hadrons below the constituent mass threshold formed by the given particles,
* together with their total weight
*
* Given the id of two (or three) constituents of a cluster, it returns
* the lightest hadron with proper flavour numbers.
* At the moment it does *nothing* in the case that also 'ptr3' present.
* @param threshold The theshold
* @param ptr1 is the first constituent
* @param ptr2 is the second constituent
* @param ptr3 is the third constituent
*/
vector<pair<tcPDPtr,double> >
hadronsBelowThreshold(Energy threshold,
tcPDPtr ptr1, tcPDPtr ptr2,
tcPDPtr ptr3 = PDPtr ()) const;
/**
* Return the nominal mass of the hadron returned by lightestHadron()
* @param ptr1 is the first constituent
* @param ptr2 is the second constituent
* @param ptr3 is the third constituent
*/
Energy massLightestHadron(tcPDPtr ptr1, tcPDPtr ptr2,
#ifndef NDEBUG
tcPDPtr ptr3 = PDPtr ()) const {
#else
tcPDPtr = PDPtr ()) const {
#endif
// The method assumes ptr3 == empty
assert(!(ptr3));
// find entry in the table
pair<long,long> ids(abs(ptr1->id()),abs(ptr2->id()));
HadronTable::const_iterator tit=_table.find(ids);
// throw exception if flavours wrong
if(tit==_table.end()||tit->second.empty())
throw Exception() << "HadronSelector::massLightestHadron "
<< "failed for particle" << ptr1->id() << " "
<< ptr2->id()
<< Exception::eventerror;
// return the mass
return tit->second.begin()->mass;
}
/**
* Returns the mass of the lightest pair of baryons.
* @param ptr1 is the first constituent
* @param ptr2 is the second constituent
*/
Energy massLightestBaryonPair(tcPDPtr ptr1, tcPDPtr ptr2) const;
/**
* Access the parton weights
*/
double pwt(long pid) const {
map<long,double>::const_iterator it = _pwt.find(abs(pid));
assert( it != _pwt.end() );
return it->second;
}
+ /**
+ * Force baryon/meson selection
+ */
+ virtual pair<bool,bool> selectBaryon(const Energy cluMass, tcPDPtr par1, tcPDPtr par2) const;
+
+ /**
+ * Strange quark weight
+ */
+ virtual double strangeWeight(const Energy cluMass, tcPDPtr par1, tcPDPtr par2) const;
+
+
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
*
* The array _repwt is initialized using the interfaces to set different
* weights for different meson multiplets and the constructHadronTable()
* method called to complete the construction of the hadron tables.
*
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
//@}
protected:
/**
* A sub-function of HadronSelector::constructHadronTable().
* It receives the information of a prospective Hadron and inserts it
* into the hadron table construct.
* @param particle is a particle data pointer to the hadron
* @param flav1 is the first constituent of the hadron
* @param flav2 is the second constituent of the hadron
*/
void insertToHadronTable(tPDPtr &particle, int flav1, int flav2);
/**
* Construct the table of hadron data
* This is the main method to initialize the hadron data (mainly the
* weights associated to each hadron, taking into account its spin,
* eventual isoscalar-octect mixing, singlet-decuplet factor). This is
* the method that one should update when new or updated hadron data is
* available.
*
* This class implements the construction of the basic table but can be
* overridden if needed in inheriting classes.
*
* The rationale for factors used for diquarks involving different quarks can
* be can be explained by taking a prototype example that in the exact SU(2) limit,
* in which:
* \f[m_u=m_d\f]
* \f[M_p=M_n=M_\Delta\f]
* and we will have equal numbers of u and d quarks produced.
* Suppose that we weight 1 the diquarks made of the same
* quark and 1/2 those made of different quarks, the fractions
* of u and d baryons (p, n, Delta) we get are the following:
* - \f$\Delta^{++}\f$: 1 possibility only u uu with weight 1
* - \f$\Delta^- \f$: 1 possibility only d dd with weight 1
* - \f$p,\Delta^+ \f$: 2 possibilities u ud with weight 1/2
* d uu with weight 1
* - \f$n,\Delta^0 \f$: 2 possibilities d ud with weight 1/2
* u dd with weight 1
* In the latter two cases, we have to take into account the
* fact that p and n have spin 1/2 whereas Delta+ and Delta0
* have spin 3/2 therefore from phase space we get a double weight
* for Delta+ and Delta0 relative to p and n respectively.
* Therefore the relative amount of these baryons that is
* produced is the following:
* # p = # n = ( 1/2 + 1 ) * 1/3 = 1/2
* # Delta++ = # Delta- = 1 = ( 1/2 + 1) * 2/3 # Delta+ = # Delta0
* which is correct, and therefore the weight 1/2 for the
* diquarks of different types of quarks is justified (at least
* in this limit of exact SU(2) ).
*/
virtual void constructHadronTable();
/**
* Access to the table of hadrons
*/
const HadronTable & table() const {
return _table;
}
/**
* Access to the list of partons
*/
const vector<PDPtr> & partons() const {
return _partons;
}
/**
* Access to the list of partons
*/
vector<PDPtr> & partons() {
return _partons;
}
/**
* Access the parton weights
*/
map<long,double> & pwt() {
return _pwt;
}
/**
* Methods for the mixing of \f$I=0\f$ mesons
*/
//@{
/**
* Return the probability of mixing for Octet-Singlet isoscalar mixing,
* the probability of the
* \f$\frac1{\sqrt{2}}(|u\bar{u}\rangle + |d\bar{d}\rangle)\f$ component
* is returned.
* @param angleMix The mixing angle in degrees (not radians)
* @param order is 0 for no mixing, 1 for the first resonance of a pair,
* 2 for the second one.
* The mixing is defined so that for example with \f$\eta-\eta'\f$ mixing where
* the mixing angle is \f$\theta=-23^0$ with $\eta\f$ as the first particle
* and \f$\eta'\f$ the second one.
* The convention used is
* \f[\eta = \cos\theta|\eta_{\rm octet }\rangle
* -\sin\theta|\eta_{\rm singlet}\rangle\f]
* \f[\eta' = \sin\theta|\eta_{\rm octet }\rangle
* -\cos\theta|\eta_{\rm singlet}\rangle\f]
* with
* \f[|\eta_{\rm singlet}\rangle = \frac1{\sqrt{3}}
* \left[|u\bar{u}\rangle + |d\bar{d}\rangle + |s\bar{s}\rangle\right]\f]
* \f[|\eta_{\rm octet }\rangle = \frac1{\sqrt{6}}
* \left[|u\bar{u}\rangle + |d\bar{d}\rangle - 2|s\bar{s}\rangle\right]\f]
*/
double probabilityMixing(const double angleMix,
const int order) const {
static double convert=Constants::pi/180.0;
if (order == 1)
return sqr( cos( angleMix*convert + atan( sqrt(2.0) ) ) );
else if (order == 2)
return sqr( sin( angleMix*convert + atan( sqrt(2.0) ) ) );
else
return 1.;
}
/**
* Returns the weight of given mixing state.
* @param id The PDG code of the meson
*/
virtual double mixingStateWeight(long id) const;
//@}
/**
* Calculates a special weight specific to a given hadron.
* @param id The PDG code of the hadron
*/
double specialWeight(long id) const {
const int pspin = id % 10;
// Only K0L and K0S have pspin == 0, should
// not get them until Decay step
assert( pspin != 0 );
// Baryon : J = 1/2 or 3/2
if(pspin%2==0)
return baryonWeight(id);
// Meson
else
return mesonWeight(id);
}
/**
* Weights for mesons
*/
virtual double mesonWeight(long id) const;
/**
* Weights for baryons
*/
virtual double baryonWeight(long id) const = 0;
/**
* This method returns the proper sign ( > 0 hadron; < 0 anti-hadron )
* for the input PDG id idHad > 0, suppose to be made by the
* two constituent particle pointers: par1 and par2 (both with proper sign).
*/
int signHadron(tcPDPtr ptr1, tcPDPtr ptr2, tcPDPtr hadron) const;
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
HadronSelector & operator=(const HadronSelector &) = delete;
private:
/**
* The PDG codes of the constituent particles allowed
*/
vector<PDPtr> _partons;
/**
* The PDG codes of the hadrons which cannot be produced in the hadronization
*/
vector<PDPtr> _forbidden;
/**
* Weights for quarks and diquarks.
*/
map<long,double> _pwt;
/**
* The mixing angles for the \f$I=0\f$ mesons containing light quarks
*/
//@{
/**
* The \f$\eta-\eta'\f$ mixing angle
*/
double _etamix;
/**
* The \f$\phi-\omega\f$ mixing angle
*/
double _phimix;
/**
* The \f$h_1'-h_1\f$ mixing angle
*/
double _h1mix;
/**
* The \f$f_0(1710)-f_0(1370)\f$ mixing angle
*/
double _f0mix;
/**
* The \f$f_1(1420)-f_1(1285)\f$ mixing angle
*/
double _f1mix;
/**
* The \f$f'_2-f_2\f$ mixing angle
*/
double _f2mix;
/**
* The \f$\eta_2(1870)-\eta_2(1645)\f$ mixing angle
*/
double _eta2mix;
/**
* The \f$\phi(???)-\omega(1650)\f$ mixing angle
*/
double _omhmix;
/**
* The \f$\phi_3-\omega_3\f$ mixing angle
*/
double _ph3mix;
/**
* The \f$\eta(1475)-\eta(1295)\f$ mixing angle
*/
double _eta2Smix;
/**
* The \f$\phi(1680)-\omega(1420)\f$ mixing angle
*/
double _phi2Smix;
//@}
/**
* The weights for the various meson multiplets to be used to supress the
* production of particular states
*/
//@{
/**
* The weights for the \f$\phantom{1}^1S_0\f$ multiplets
*/
vector<double> _weight1S0;
/**
* The weights for the \f$\phantom{1}^3S_1\f$ multiplets
*/
vector<double> _weight3S1;
/**
* The weights for the \f$\phantom{1}^1P_1\f$ multiplets
*/
vector<double> _weight1P1;
/**
* The weights for the \f$\phantom{1}^3P_0\f$ multiplets
*/
vector<double> _weight3P0;
/**
* The weights for the \f$\phantom{1}^3P_1\f$ multiplets
*/
vector<double> _weight3P1;
/**
* The weights for the \f$\phantom{1}^3P_2\f$ multiplets
*/
vector<double> _weight3P2;
/**
* The weights for the \f$\phantom{1}^1D_2\f$ multiplets
*/
vector<double> _weight1D2;
/**
* The weights for the \f$\phantom{1}^3D_1\f$ multiplets
*/
vector<double> _weight3D1;
/**
* The weights for the \f$\phantom{1}^3D_2\f$ multiplets
*/
vector<double> _weight3D2;
/**
* The weights for the \f$\phantom{1}^3D_3\f$ multiplets
*/
vector<double> _weight3D3;
//@}
/**
* The weights for the excited meson multiplets
*/
vector<vector<vector<double> > > _repwt;
/**
* The table of hadron data
*/
HadronTable _table;
/**
* Enums so arrays can be statically allocated
*/
//@{
/**
* Defines values for array sizes. L,J,N max values for excited mesons.
*/
enum MesonMultiplets { Lmax = 3, Jmax = 4, Nmax = 4};
//@}
/**
* Option for the construction of the tables
*/
unsigned int _topt;
/**
* Which particles to produce for debugging purposes
*/
unsigned int _trial;
/**
* @name A parameter used for determining when clusters are too light.
*
* This parameter is used for setting the lower threshold, \f$ t \f$ as
* \f[ t' = t(1 + r B^1_{\rm lim}) \f]
* where \f$ r \f$ is a random number [0,1].
*/
//@{
double _limBottom;
double _limCharm;
double _limExotic;
//@}
/**
* Option for the selection of hadrons below the pair threshold
*/
unsigned int belowThreshold_;
};
}
#endif /* HERWIG_HadronSelector_H */
diff --git a/Hadronization/Hw64Selector.cc b/Hadronization/Hw64Selector.cc
--- a/Hadronization/Hw64Selector.cc
+++ b/Hadronization/Hw64Selector.cc
@@ -1,212 +1,212 @@
// -*- C++ -*-
//
// Hw64Selector.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 Hw64Selector class.
//
#include "Hw64Selector.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Repository/UseRandom.h"
#include "CheckId.h"
#include "Herwig/Utilities/Kinematics.h"
#include <ThePEG/Utilities/DescribeClass.h>
using namespace Herwig;
DescribeClass<Hw64Selector,HadronSelector>
describeHw64Selector("Herwig::Hw64Selector","Herwig.so");
IBPtr Hw64Selector::clone() const {
return new_ptr(*this);
}
IBPtr Hw64Selector::fullclone() const {
return new_ptr(*this);
}
void Hw64Selector::persistentOutput(PersistentOStream & os) const {
os << _pwtDquark << _pwtUquark << _pwtSquark
<< _pwtCquark << _pwtBquark << _pwtDIquarkS0 << _pwtDIquarkS1
<< _sngWt << _decWt ;
}
void Hw64Selector::persistentInput(PersistentIStream & is, int) {
is >> _pwtDquark >> _pwtUquark >> _pwtSquark
>> _pwtCquark >> _pwtBquark >> _pwtDIquarkS0 >> _pwtDIquarkS1
>> _sngWt >> _decWt ;
}
void Hw64Selector::doinit() {
// the default partons allowed
// the quarks
for ( int ix=1; ix<=5; ++ix ) {
partons().push_back(getParticleData(ix));
}
// the diquarks
for(unsigned int ix=1;ix<=5;++ix) {
for(unsigned int iy=1; iy<=ix;++iy) {
partons().push_back(getParticleData(CheckId::makeDiquarkID(ix,iy,long(3))));
if(ix!=iy)
partons().push_back(getParticleData(CheckId::makeDiquarkID(ix,iy,long(1))));
}
}
// weights for the different quarks etc
for(unsigned int ix=0; ix<partons().size(); ++ix) {
pwt()[partons()[ix]->id()]=0.;
}
pwt()[1] = _pwtDquark;
pwt()[2] = _pwtUquark;
pwt()[3] = _pwtSquark;
pwt()[4] = _pwtCquark;
pwt()[5] = _pwtBquark;
pwt()[1103] = _pwtDIquarkS1 * _pwtDquark * _pwtDquark;
pwt()[2101] = _pwtDIquarkS0 * _pwtUquark * _pwtDquark;
pwt()[2103] = _pwtDIquarkS1 * _pwtUquark * _pwtDquark;
pwt()[2203] = _pwtDIquarkS1 * _pwtUquark * _pwtUquark;
pwt()[3101] = _pwtDIquarkS0 * _pwtSquark * _pwtDquark;
pwt()[3103] = _pwtDIquarkS1 * _pwtSquark * _pwtDquark;
pwt()[3201] = _pwtDIquarkS0 * _pwtSquark * _pwtUquark;
pwt()[3203] = _pwtDIquarkS1 * _pwtSquark * _pwtUquark;
pwt()[3303] = _pwtDIquarkS1 * _pwtSquark * _pwtSquark;
HadronSelector::doinit();
}
void Hw64Selector::Init() {
static ClassDocumentation<Hw64Selector> documentation
("The Hw64Selector class implements the hadron selection algorithm of Hw6");
static Parameter<Hw64Selector,double>
interfacePwtDquark("PwtDquark","Weight for choosing a quark D",
&Hw64Selector::_pwtDquark, 0, 1.0, 0.0, 10.0,
false,false,false);
static Parameter<Hw64Selector,double>
interfacePwtUquark("PwtUquark","Weight for choosing a quark U",
&Hw64Selector::_pwtUquark, 0, 1.0, 0.0, 10.0,
false,false,false);
static Parameter<Hw64Selector,double>
interfacePwtSquark("PwtSquark","Weight for choosing a quark S",
&Hw64Selector::_pwtSquark, 0, 1.0, 0.0, 10.0,
false,false,false);
static Parameter<Hw64Selector,double>
interfacePwtCquark("PwtCquark","Weight for choosing a quark C",
&Hw64Selector::_pwtCquark, 0, 0.0, 0.0, 10.0,
false,false,false);
static Parameter<Hw64Selector,double>
interfacePwtBquark("PwtBquark","Weight for choosing a quark B",
&Hw64Selector::_pwtBquark, 0, 0.0, 0.0, 10.0,
false,false,false);
static Parameter<Hw64Selector,double>
interfacePwtDIquarkS0("PwtDIquarkS0","Weight for choosing a spin-0 DIquark",
&Hw64Selector::_pwtDIquarkS0, 0, 1.0, 0.0, 100.0,
false,false,false);
static Parameter<Hw64Selector,double>
interfacePwtDIquarkS1("PwtDIquarkS1","Weight for choosing a spin-1 DIquark",
&Hw64Selector::_pwtDIquarkS1, 0, 1.0, 0.0, 100.0,
false,false,false);
static Parameter<Hw64Selector,double>
interfaceSngWt("SngWt","Weight for singlet baryons",
&Hw64Selector::_sngWt, 0, 1.0, 0.0, 10.0,
false,false,false);
static Parameter<Hw64Selector,double>
interfaceDecWt("DecWt","Weight for decuplet baryons",
&Hw64Selector::_decWt, 0, 1.0, 0.0, 10.0,
false,false,false);
}
-pair<tcPDPtr,tcPDPtr> Hw64Selector::chooseHadronPair(const Energy cluMass,tcPDPtr par1,
- tcPDPtr par2,tcPDPtr) const {
+pair<tcPDPtr,tcPDPtr> Hw64Selector::chooseHadronPair(const Energy cluMass,
+ tcPDPtr par1, tcPDPtr par2) const {
bool diquark = !(DiquarkMatcher::Check(par1->id()) || DiquarkMatcher::Check(par2->id()));
pair<tcPDPtr,tcPDPtr> lighthad = lightestHadronPair(par1, par2);
if(!lighthad.first || !lighthad.second)
throw Exception() << "Hw64Selector::chooseHadronPair "
<< "We have 0's! First id = " << par1->id() << " second = "
<< par2->id() << ". This is probably a problem with either"
<< " undecayed heavy particles or colour connections"
<< Exception::eventerror;
// calculate maximum momentum
Energy PCMax = Kinematics::pstarTwoBodyDecay(cluMass,lighthad.first->mass(),
lighthad.second->mass());
tcPDPtr had1 = tcPDPtr();
tcPDPtr had2 = tcPDPtr();
int ntry = 0;
tcPDPtr quark = tcPDPtr();
const int nmax = 5000;
Energy p;
do {
quark = partons()[UseRandom::irnd(partons().size())];
if(diquark && DiquarkMatcher::Check(quark->id())) continue;
if(pwt(quark->id()) <= UseRandom::rnd()) continue;
pair<long,long> pid(abs(par1->id()),quark->id());
KupcoData::const_iterator it1,it2;
const HadronTable::const_iterator tit = table().find(pid);
assert(tit != table().end());
const KupcoData & hdata = tit->second;
do {
it1 = hdata.begin();
advance(it1,int(hdata.size()*UseRandom::rnd()));
}
while(it1 != hdata.end() && it1->overallWeight < UseRandom::rnd());
had1 = it1->ptrData;
pid = make_pair(quark->id(),abs(par2->id()));
do {
it2 = hdata.begin();
advance(it2,int(hdata.size()*UseRandom::rnd()));
}
while(it2 != hdata.end() && it2->overallWeight < UseRandom::rnd());
had2 = it2->ptrData;
if(had1 && had2) {
p = Kinematics::pstarTwoBodyDecay(cluMass, it1->mass, it2->mass);
if(p/PCMax < UseRandom::rnd()) {
had1 = had2 = tcPDPtr();
ntry++;
}
}
}
while((!had1|| !had2) && ntry < nmax);
if(ntry >= nmax) return lighthad;
int signHad1 = 0;
int signHad2 = 0;
if(CheckId::canBeHadron(par1,quark->CC()) && CheckId::canBeHadron(quark,par2)) {
signHad1 = signHadron(par1, quark->CC(), had1);
signHad2 = signHadron(par2, quark, had2);
}
else if(CheckId::canBeHadron(par1,quark) && CheckId::canBeHadron(quark->CC(),par2)) {
signHad1 = signHadron(par1, quark, had1);
signHad2 = signHadron(par2, quark->CC(), had2);
}
else throw Exception() << "Hw64Selector::chooseHadronPair()"
<< Exception::abortnow;
return make_pair( signHad1 > 0 ? had1 : tcPDPtr(had1->CC()),
signHad2 > 0 ? had2 : tcPDPtr(had2->CC()));
}
double Hw64Selector::baryonWeight(long id) const {
const int pspin = id % 10;
if(pspin == 2) {
// Singlet (Lambda-like) baryon
if( (id/100)%10 < (id/10 )%10 ) return sqr(_sngWt);
}
// Decuplet baryon
else if (pspin == 4) return sqr(_decWt);
return 1.;
}
diff --git a/Hadronization/Hw64Selector.h b/Hadronization/Hw64Selector.h
--- a/Hadronization/Hw64Selector.h
+++ b/Hadronization/Hw64Selector.h
@@ -1,184 +1,184 @@
// -*- C++ -*-
//
// Hw64Selector.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_Hw64Selector_H
#define HERWIG_Hw64Selector_H
//
// This is the declaration of the Hw64Selector class.
//
#include "HadronSelector.h"
#include "Hw64Selector.fh"
namespace Herwig {
using namespace ThePEG;
/** \ingroup hadronization
* The Hw64Selector class selects the hadrons produced in cluster decay using
* the FORTRAN HERWIG variant of the cluster model.
*
* @see \ref Hw64SelectorInterfaces "The interfaces"
* defined for Hw64Selector.
*/
class Hw64Selector: public HadronSelector {
public:
/**
* The default constructor.
*/
Hw64Selector() : HadronSelector(0),
_pwtDquark( 1.0 ),_pwtUquark( 1.0 ),_pwtSquark( 1.0 ),_pwtCquark( 0.0 ),
_pwtBquark( 0.0 ),_pwtDIquarkS0( 1.0 ),_pwtDIquarkS1( 1.0 ),
_sngWt( 1.0 ), _decWt( 1.0 )
{}
/**
* Method to return a pair of hadrons given the PDG codes of
* two or three constituents
* @param cluMass The mass of the cluster
* @param par1 The particle pointer of the first constituent
* @param par2 The particle pointer of the second constituent
* @param par3 The particle pointer of the third constituent
*/
- virtual pair<tcPDPtr,tcPDPtr> chooseHadronPair(const Energy cluMass,tcPDPtr par1,
- tcPDPtr par2,tcPDPtr par3 = PDPtr()) const;
+ virtual pair<tcPDPtr,tcPDPtr> chooseHadronPair(const Energy cluMass,
+ tcPDPtr par1, tcPDPtr par2) const;
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected :
/**
* Weights for baryons
*/
virtual double baryonWeight(long id) const;
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
Hw64Selector & operator=(const Hw64Selector &) = delete;
private:
/**
* The weights for the different quarks and diquarks
*/
//@{
/**
* The probability of producting a down quark.
*/
double _pwtDquark;
/**
* The probability of producting an up quark.
*/
double _pwtUquark;
/**
* The probability of producting a strange quark.
*/
double _pwtSquark;
/**
* The probability of producting a charm quark.
*/
double _pwtCquark;
/**
* The probability of producting a bottom quark.
*/
double _pwtBquark;
/**
* The probability of producting a spin-0 diquark.
*/
double _pwtDIquarkS0;
/**
* The probability of producting a spin-1 diquark.
*/
double _pwtDIquarkS1;
//@}
/**
* Singlet and Decuplet weights
*/
//@{
/**
* The singlet weight
*/
double _sngWt;
/**
* The decuplet weight
*/
double _decWt;
//@}
};
}
#endif /* HERWIG_Hw64Selector_H */
diff --git a/Hadronization/HwppSelector.cc b/Hadronization/HwppSelector.cc
--- a/Hadronization/HwppSelector.cc
+++ b/Hadronization/HwppSelector.cc
@@ -1,364 +1,281 @@
// -*- C++ -*-
//
// HwppSelector.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 HwppSelector class.
//
#include "HwppSelector.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/Utilities/Kinematics.h"
#include "ThePEG/Utilities/Selector.h"
#include "ThePEG/Repository/UseRandom.h"
#include "CheckId.h"
#include <cassert>
#include <ThePEG/Utilities/DescribeClass.h>
using namespace Herwig;
DescribeClass<HwppSelector,HadronSelector>
describeHwppSelector("Herwig::HwppSelector","Herwig.so");
IBPtr HwppSelector::clone() const {
return new_ptr(*this);
}
IBPtr HwppSelector::fullclone() const {
return new_ptr(*this);
}
void HwppSelector::doinit() {
// the default partons allowed
// the quarks
for ( int ix=1; ix<=5; ++ix ) {
partons().push_back(getParticleData(ix));
}
// the diquarks
for(unsigned int ix=1;ix<=5;++ix) {
for(unsigned int iy=1; iy<=ix;++iy) {
partons().push_back(getParticleData(CheckId::makeDiquarkID(ix,iy,long(3))));
if(ix!=iy)
partons().push_back(getParticleData(CheckId::makeDiquarkID(ix,iy,long(1))));
}
}
// weights for the different quarks etc
for(unsigned int ix=0; ix<partons().size(); ++ix) {
pwt()[partons()[ix]->id()]=0.;
}
pwt()[1] = _pwtDquark;
pwt()[2] = _pwtUquark;
pwt()[3] = _pwtSquark;
pwt()[4] = _pwtCquark;
pwt()[5] = _pwtBquark;
pwt()[1103] = _pwtDIquarkS1 * _pwtDquark * _pwtDquark;
pwt()[2101] = _pwtDIquarkS0 * _pwtUquark * _pwtDquark;
pwt()[2103] = _pwtDIquarkS1 * _pwtUquark * _pwtDquark;
pwt()[2203] = _pwtDIquarkS1 * _pwtUquark * _pwtUquark;
pwt()[3101] = _pwtDIquarkS0 * _pwtSquark * _pwtDquark;
pwt()[3103] = _pwtDIquarkS1 * _pwtSquark * _pwtDquark;
pwt()[3201] = _pwtDIquarkS0 * _pwtSquark * _pwtUquark;
pwt()[3203] = _pwtDIquarkS1 * _pwtSquark * _pwtUquark;
pwt()[3303] = _pwtDIquarkS1 * _pwtSquark * _pwtSquark;
HadronSelector::doinit();
}
void HwppSelector::persistentOutput(PersistentOStream & os) const {
os << _pwtDquark << _pwtUquark << _pwtSquark
<< _pwtCquark << _pwtBquark << _pwtDIquarkS0 << _pwtDIquarkS1
<< _sngWt << _decWt
<< _mode << _enhanceSProb << ounit(_m0Decay,GeV) << _massMeasure
<< _scHadronWtFactor << _sbHadronWtFactor;
}
void HwppSelector::persistentInput(PersistentIStream & is, int) {
is >> _pwtDquark >> _pwtUquark >> _pwtSquark
>> _pwtCquark >> _pwtBquark >> _pwtDIquarkS0 >> _pwtDIquarkS1
>> _sngWt >> _decWt
>> _mode >> _enhanceSProb >> iunit(_m0Decay,GeV) >> _massMeasure
>> _scHadronWtFactor >> _sbHadronWtFactor;
}
void HwppSelector::Init() {
static ClassDocumentation<HwppSelector> documentation
("The HwppSelector class implements the Herwig algorithm for selecting"
" the hadrons",
"The hadronization used the selection algorithm described in \\cite{Kupco:1998fx}.",
"%\\cite{Kupco:1998fx}\n"
"\\bibitem{Kupco:1998fx}\n"
" A.~Kupco,\n"
" ``Cluster hadronization in HERWIG 5.9,''\n"
" arXiv:hep-ph/9906412.\n"
" %%CITATION = HEP-PH/9906412;%%\n"
);
static Parameter<HwppSelector,double>
interfacePwtDquark("PwtDquark","Weight for choosing a quark D",
&HwppSelector::_pwtDquark, 0, 1.0, 0.0, 10.0,
false,false,false);
static Parameter<HwppSelector,double>
interfacePwtUquark("PwtUquark","Weight for choosing a quark U",
&HwppSelector::_pwtUquark, 0, 1.0, 0.0, 10.0,
false,false,false);
static Parameter<HwppSelector,double>
interfacePwtSquark("PwtSquark","Weight for choosing a quark S",
&HwppSelector::_pwtSquark, 0, 1.0, 0.0, 10.0,
false,false,false);
static Parameter<HwppSelector,double>
interfacePwtCquark("PwtCquark","Weight for choosing a quark C",
&HwppSelector::_pwtCquark, 0, 0.0, 0.0, 10.0,
false,false,false);
static Parameter<HwppSelector,double>
interfacePwtBquark("PwtBquark","Weight for choosing a quark B",
&HwppSelector::_pwtBquark, 0, 0.0, 0.0, 10.0,
false,false,false);
static Parameter<HwppSelector,double>
interfacePwtDIquarkS0("PwtDIquarkS0","Weight for choosing a spin-0 DIquark",
&HwppSelector::_pwtDIquarkS0, 0, 1.0, 0.0, 100.0,
false,false,false);
static Parameter<HwppSelector,double>
interfacePwtDIquarkS1("PwtDIquarkS1","Weight for choosing a spin-1 DIquark",
&HwppSelector::_pwtDIquarkS1, 0, 1.0, 0.0, 100.0,
false,false,false);
static Parameter<HwppSelector,double>
interfaceSngWt("SngWt","Weight for singlet baryons",
&HwppSelector::_sngWt, 0, 1.0, 0.0, 10.0,
false,false,false);
static Parameter<HwppSelector,double>
interfaceDecWt("DecWt","Weight for decuplet baryons",
&HwppSelector::_decWt, 0, 1.0, 0.0, 10.0,
false,false,false);
static Switch<HwppSelector,unsigned int> interfaceMode
("Mode",
"Which algorithm to use",
&HwppSelector::_mode, 1, false, false);
static SwitchOption interfaceModeKupco
(interfaceMode,
"Kupco",
"Use the Kupco approach",
0);
static SwitchOption interfaceModeHwpp
(interfaceMode,
"Hwpp",
"Use the Herwig approach",
1);
static Switch<HwppSelector,int> interfaceEnhanceSProb
("EnhanceSProb",
"Option for enhancing strangeness",
&HwppSelector::_enhanceSProb, 0, false, false);
static SwitchOption interfaceEnhanceSProbNo
(interfaceEnhanceSProb,
"No",
"No strangeness enhancement.",
0);
static SwitchOption interfaceEnhanceSProbScaled
(interfaceEnhanceSProb,
"Scaled",
"Scaled strangeness enhancement",
1);
static SwitchOption interfaceEnhanceSProbExponential
(interfaceEnhanceSProb,
"Exponential",
"Exponential strangeness enhancement",
2);
static Switch<HwppSelector,int> interfaceMassMeasure
("MassMeasure",
"Option to use different mass measures",
&HwppSelector::_massMeasure,0,false,false);
static SwitchOption interfaceMassMeasureMass
(interfaceMassMeasure,
"Mass",
"Mass Measure",
0);
static SwitchOption interfaceMassMeasureLambda
(interfaceMassMeasure,
"Lambda",
"Lambda Measure",
1);
static Parameter<HwppSelector,double> interfacescHadronWtFactor
("scHadronWtFactor",
"Wight factor for strenge-charm heavy hadrns",
&HwppSelector::_scHadronWtFactor, 1., 0., 10.,
false, false, Interface::limited);
static Parameter<HwppSelector,double> interfacesbHadronWtFactor
("sbHadronWtFactor",
"Wight factor for strenge-bottom heavy hadrns",
&HwppSelector::_sbHadronWtFactor, 1., 0., 10.,
false, false, Interface::limited);
static Parameter<HwppSelector,Energy> interfaceDecayMassScale
("DecayMassScale",
"Cluster decay mass scale",
&HwppSelector::_m0Decay, GeV, 1.0*GeV, 0.1*GeV, 50.*GeV,
false, false, Interface::limited);
}
-pair<tcPDPtr,tcPDPtr> HwppSelector::chooseHadronPair(const Energy cluMass,tcPDPtr par1,
- tcPDPtr par2,tcPDPtr ) const {
- // if either of the input partons is a diquark don't allow diquarks to be
- // produced
- bool diquark = !(DiquarkMatcher::Check(par1->id()) || DiquarkMatcher::Check(par2->id()));
- bool quark = true;
- // if the Herwig algorithm
+double HwppSelector::baryonWeight(long id) const {
+ const int pspin = id % 10;
+ if(pspin == 2) {
+ // Singlet (Lambda-like) baryon
+ if( (id/100)%10 < (id/10 )%10 ) return sqr(_sngWt);
+ }
+ // Decuplet baryon
+ else if (pspin == 4) return sqr(_decWt);
+ return 1.;
+}
+
+pair<bool,bool> HwppSelector::selectBaryon(const Energy cluMass, tcPDPtr par1, tcPDPtr par2) const {
+ bool quark=true, diquark=true;
+ useMe();
if(_mode ==1) {
if(UseRandom::rnd() > 1./(1.+_pwtDIquarkS0+_pwtDIquarkS1)
&& cluMass > massLightestBaryonPair(par1,par2)) {
diquark = true;
quark = false;
}
else {
useMe();
diquark = false;
quark = true;
}
}
- // weights for the different possibilities
- Energy weight, wgtsum(ZERO);
- // loop over all hadron pairs with the allowed flavours
- static vector<Kupco> hadrons;
- hadrons.clear();
- for(unsigned int ix=0;ix<partons().size();++ix) {
- tcPDPtr quarktopick = partons()[ix];
- if(!quark && abs(int(quarktopick->iColour())) == 3
- && !DiquarkMatcher::Check(quarktopick->id())) continue;
- if(!diquark && abs(int(quarktopick->iColour())) == 3
- && DiquarkMatcher::Check(quarktopick->id())) continue;
- HadronTable::const_iterator
- tit1 = table().find(make_pair(abs(par1->id()),quarktopick->id()));
- HadronTable::const_iterator
- tit2 = table().find(make_pair(quarktopick->id(),abs(par2->id())));
- // If not in table skip
- if(tit1 == table().end()||tit2==table().end()) continue;
- // tables empty skip
- const KupcoData & T1 = tit1->second;
- const KupcoData & T2 = tit2->second;
- if(T1.empty()||T2.empty()) continue;
- // if too massive skip
- if(cluMass <= T1.begin()->mass +
- T2.begin()->mass) continue;
- // quark weight
- double quarkWeight = pwt(quarktopick->id());
- if(abs(quarktopick->id()) == 3) {
- // Decoupling the weight of heavy strenge hadrons
- if(_enhanceSProb == 0 && abs(par1->id()) == 4) {
- quarkWeight = pwt(quarktopick->id())*_scHadronWtFactor;
- }
- else if(_enhanceSProb == 0 && abs(par1->id()) == 5) {
- quarkWeight = pwt(quarktopick->id())*_sbHadronWtFactor;
- }
- // Scaling strangeness enhancement
- else if(_enhanceSProb == 1) {
- double scale = double(sqr(_m0Decay/cluMass));
- quarkWeight = (_maxScale < scale) ? 0. : pow(quarkWeight,scale);
- }
- // Exponential strangeness enhancement
- else if(_enhanceSProb == 2) {
- Energy2 mass2;
- Energy endpointmass = par1->mass() + par2->mass();
- // Choose to use either the cluster mass
- // or to use the lambda measure
- mass2 = (_massMeasure == 0) ? sqr(cluMass) :
- sqr(cluMass) - sqr(endpointmass);
- double scale = double(sqr(_m0Decay)/mass2);
- quarkWeight = (_maxScale < scale) ? 0. : exp(-scale);
- }
- }
- // loop over the hadrons
- KupcoData::const_iterator H1,H2;
- for(H1 = T1.begin();H1 != T1.end(); ++H1) {
- for(H2 = T2.begin();H2 != T2.end(); ++H2) {
- // break if cluster too light
- if(cluMass < H1->mass + H2->mass) break;
- weight = quarkWeight * H1->overallWeight * H2->overallWeight *
- Kinematics::pstarTwoBodyDecay(cluMass, H1->mass, H2->mass);
- //cerr<<H1->id<<" "<<H2->id<<" "<<weight/GeV<<"\n";
- int signQ = 0;
- assert (par1 && quarktopick);
- assert (par2);
-
- assert(quarktopick->CC());
-
- if(CheckId::canBeHadron(par1, quarktopick->CC())
- && CheckId::canBeHadron(quarktopick, par2))
- signQ = +1;
- else if(CheckId::canBeHadron(par1, quarktopick)
- && CheckId::canBeHadron(quarktopick->CC(), par2))
- signQ = -1;
- else {
- cerr << "Could not make sign for" << par1->id()<< " " << quarktopick->id()
- << " " << par2->id() << "\n";
- assert(false);
- }
-
- if (signQ == -1)
- quarktopick = quarktopick->CC();
- // construct the object with the info
- Kupco a(quarktopick, H1->ptrData, H2->ptrData, weight);
- hadrons.push_back(a);
- wgtsum += weight;
- }
- }
- }
- if (hadrons.empty())
- return make_pair(tcPDPtr(),tcPDPtr());
- // select the hadron
- wgtsum *= UseRandom::rnd();
- unsigned int ix=0;
- do {
- wgtsum-= hadrons[ix].weight;
- ++ix;
- }
- while(wgtsum > ZERO && ix < hadrons.size());
- if(ix == hadrons.size() && wgtsum > ZERO)
- return make_pair(tcPDPtr(),tcPDPtr());
- --ix;
- assert(hadrons[ix].idQ);
- int signHad1 = signHadron(par1, hadrons[ix].idQ->CC(), hadrons[ix].hadron1);
- int signHad2 = signHadron(par2, hadrons[ix].idQ, hadrons[ix].hadron2);
- assert( signHad1 != 0 && signHad2 != 0 );
- return make_pair
- ( signHad1 > 0 ? hadrons[ix].hadron1 : tcPDPtr(hadrons[ix].hadron1->CC()),
- signHad2 > 0 ? hadrons[ix].hadron2 : tcPDPtr(hadrons[ix].hadron2->CC()));
+ return make_pair(quark,diquark);
}
-double HwppSelector::baryonWeight(long id) const {
- const int pspin = id % 10;
- if(pspin == 2) {
- // Singlet (Lambda-like) baryon
- if( (id/100)%10 < (id/10 )%10 ) return sqr(_sngWt);
+double HwppSelector::strangeWeight(const Energy cluMass, tcPDPtr par1, tcPDPtr par2) const {
+ // Decoupling the weight of heavy strenge hadrons
+ if(_enhanceSProb == 0 && abs(par1->id()) == 4) {
+ return pwt(3)*_scHadronWtFactor;
}
- // Decuplet baryon
- else if (pspin == 4) return sqr(_decWt);
- return 1.;
+ else if(_enhanceSProb == 0 && abs(par1->id()) == 5) {
+ return pwt(3)*_sbHadronWtFactor;
+ }
+ // Scaling strangeness enhancement
+ else if(_enhanceSProb == 1) {
+ double scale = double(sqr(_m0Decay/cluMass));
+ return (_maxScale < scale) ? 0. : pow(pwt(3),scale);
+ }
+ // Exponential strangeness enhancement
+ else if(_enhanceSProb == 2) {
+ Energy2 mass2;
+ Energy endpointmass = par1->mass() + par2->mass();
+ // Choose to use either the cluster mass
+ // or to use the lambda measure
+ mass2 = (_massMeasure == 0) ? sqr(cluMass) :
+ sqr(cluMass) - sqr(endpointmass);
+ double scale = double(sqr(_m0Decay)/mass2);
+ return (_maxScale < scale) ? 0. : exp(-scale);
+ }
+ return pwt(3);
}
diff --git a/Hadronization/HwppSelector.h b/Hadronization/HwppSelector.h
--- a/Hadronization/HwppSelector.h
+++ b/Hadronization/HwppSelector.h
@@ -1,247 +1,225 @@
// -*- C++ -*-
//
// HwppSelector.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_HwppSelector_H
#define HERWIG_HwppSelector_H
//
// This is the declaration of the HwppSelector class.
//
#include "HadronSelector.h"
#include "HwppSelector.fh"
namespace Herwig {
using namespace ThePEG;
/** \ingroup hadronization
* The HwppSelector class selects the hadrons produced in cluster decay using
* the Herwig variant of the cluster model.
*
* @see \ref HwppSelectorInterfaces "The interfaces"
* defined for HwppSelector.
*/
class HwppSelector: public HadronSelector {
public:
/**
* The default constructor.
*/
HwppSelector() : HadronSelector(1),
_pwtDquark( 1.0 ),_pwtUquark( 1.0 ),_pwtSquark( 1.0 ),_pwtCquark( 0.0 ),
_pwtBquark( 0.0 ),_pwtDIquarkS0( 1.0 ),_pwtDIquarkS1( 1.0 ),
_sngWt( 1.0 ), _decWt( 1.0 ),
_mode(1), _enhanceSProb(0), _m0Decay(1.*GeV),
_scHadronWtFactor(1.), _sbHadronWtFactor(1.)
{}
- /**
- *
- * This method is used to choose a pair of hadrons.
- *
- * Given the mass of a cluster and the particle pointers of its
- * two (or three) constituents, this returns the pair of particle pointers of
- * the two hadrons with proper flavour numbers.
- * Furthermore, the first of the two hadron must have the
- * constituent with par1, and the second must have the constituent with par2.
- * At the moment it does *nothing* in the case that also par3 is present.
- *
- * Kupco's method is used, rather than one used in FORTRAN HERWIG
- * The idea is to build on the fly a table of all possible pairs
- * of hadrons (Had1,Had2) (that we can call "cluster decay channels")
- * which are kinematically above threshold and have flavour
- * Had1=(par1,quarktopick->CC()), Had2=(quarktopick,par2), where quarktopick
- * is the poniter of:
- * --- d, u, s, c, b
- * if either par1 or par2 is a diquark;
- * --- d, u, s, c, b, dd, ud, uu, sd, su, ss,
- * cd, cu, cs, cc, bd, bu, bs, bc, bb
- * if both par1 and par2 are quarks.
- * The weight associated with each channel is given by the product
- * of: the phase space available including the spin factor 2*J+1,
- * the constant weight factor for chosen idQ,
- * the octet-singlet isoscalar mixing factor, and finally
- * the singlet-decuplet weight factor.
- */
- pair<tcPDPtr,tcPDPtr> chooseHadronPair(const Energy cluMass,tcPDPtr par1,
- tcPDPtr par2,tcPDPtr par3 = PDPtr()) const
- ;
-
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/**
* Weights for baryons
*/
virtual double baryonWeight(long id) const;
+ /**
+ * Whether to select a meson or a baryon
+ */
+ pair<bool,bool> selectBaryon(const Energy cluMass, tcPDPtr par1, tcPDPtr par2) const;
+
+ /**
+ * Strange quark weight
+ */
+ virtual double strangeWeight(const Energy cluMass, tcPDPtr par1, tcPDPtr par2) const;
+
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
HwppSelector & operator=(const HwppSelector &) = delete;
private:
/**
* The weights for the different quarks and diquarks
*/
//@{
/**
* The probability of producting a down quark.
*/
double _pwtDquark;
/**
* The probability of producting an up quark.
*/
double _pwtUquark;
/**
* The probability of producting a strange quark.
*/
double _pwtSquark;
/**
* The probability of producting a charm quark.
*/
double _pwtCquark;
/**
* The probability of producting a bottom quark.
*/
double _pwtBquark;
/**
* The probability of producting a spin-0 diquark.
*/
double _pwtDIquarkS0;
/**
* The probability of producting a spin-1 diquark.
*/
double _pwtDIquarkS1;
//@}
/**
* Singlet and Decuplet weights
*/
//@{
/**
* The singlet weight
*/
double _sngWt;
/**
* The decuplet weight
*/
double _decWt;
//@}
private:
/**
* Which algorithm to use
*/
unsigned int _mode;
/**
* Flag that switches between no strangeness enhancement, scaling enhancement,
* and exponential enhancement (in numerical order)
*/
int _enhanceSProb;
/**
* Parameter that governs the strangeness enhancement scaling
*/
Energy _m0Decay;
/**
* Flag that switches between mass measures used in strangeness enhancement:
* cluster mass, or the lambda measure - ( m_{clu}^2 - (m_q + m_{qbar})^2 )
*/
int _massMeasure;
/**
* Constant variable that stops the scale in strangeness enhancement from
* becoming too large
*/
const double _maxScale = 20.;
/**
* Heavy strange-charm hadron wight coefficient
*/
double _scHadronWtFactor;
/**
* Heavy strange-bottom hadron wight coefficient
*/
double _sbHadronWtFactor;
};
}
#endif /* HERWIG_HwppSelector_H */
diff --git a/src/defaults/Hadronization.in b/src/defaults/Hadronization.in
--- a/src/defaults/Hadronization.in
+++ b/src/defaults/Hadronization.in
@@ -1,102 +1,101 @@
# -*- ThePEG-repository -*-
############################################################
# Setup of default hadronization
#
# There are no user servicable parts inside.
#
# Anything that follows below should only be touched if you
# know what you're doing.
#############################################################
cd /Herwig/Particles
create ThePEG::ParticleData Cluster
setup Cluster 81 Cluster 0.00990 0.0 0.0 0.0 0 0 0 1
create ThePEG::ParticleData Remnant
setup Remnant 82 Remnant 0.00990 0.0 0.0 0.0 0 0 0 1
mkdir /Herwig/Hadronization
cd /Herwig/Hadronization
create Herwig::ClusterHadronizationHandler ClusterHadHandler
create Herwig::PartonSplitter PartonSplitter
create Herwig::ClusterFinder ClusterFinder
create Herwig::ColourReconnector ColourReconnector
create Herwig::ClusterFissioner ClusterFissioner
create Herwig::LightClusterDecayer LightClusterDecayer
create Herwig::ClusterDecayer ClusterDecayer
create Herwig::HwppSelector HadronSelector
newdef ClusterHadHandler:PartonSplitter PartonSplitter
newdef ClusterHadHandler:ClusterFinder ClusterFinder
newdef ClusterHadHandler:ColourReconnector ColourReconnector
newdef ClusterHadHandler:ClusterFissioner ClusterFissioner
newdef ClusterHadHandler:LightClusterDecayer LightClusterDecayer
newdef ClusterHadHandler:ClusterDecayer ClusterDecayer
newdef ClusterHadHandler:MinVirtuality2 0.1*GeV2
newdef ClusterHadHandler:MaxDisplacement 1.0e-10*millimeter
newdef ClusterHadHandler:UnderlyingEventHandler NULL
newdef ClusterFissioner:HadronSelector HadronSelector
newdef LightClusterDecayer:HadronSelector HadronSelector
newdef ClusterDecayer:HadronSelector HadronSelector
newdef ColourReconnector:ColourReconnection Yes
newdef ColourReconnector:ReconnectionProbabilityBaryonic 0.7
newdef ColourReconnector:ReconnectionProbability 0.95
newdef ColourReconnector:Algorithm Baryonic
newdef ColourReconnector:OctetTreatment All
# Clustering parameters for light quarks
newdef ClusterFissioner:ClMaxLight 3.649
newdef ClusterFissioner:ClPowLight 2.780
newdef ClusterFissioner:PSplitLight 0.899
newdef ClusterDecayer:ClDirLight 1
newdef ClusterDecayer:ClSmrLight 0.78
# Clustering parameters for b-quarks
newdef ClusterFissioner:ClMaxBottom 3.757
newdef ClusterFissioner:ClPowBottom 0.547
newdef ClusterFissioner:PSplitBottom 0.625
newdef ClusterDecayer:ClDirBottom 1
newdef ClusterDecayer:ClSmrBottom 0.078
newdef HadronSelector:SingleHadronLimitBottom 0.000
# Clustering parameters for c-quarks
newdef ClusterFissioner:ClMaxCharm 3.950
newdef ClusterFissioner:ClPowCharm 2.559
newdef ClusterFissioner:PSplitCharm 0.994
newdef ClusterDecayer:ClDirCharm 1
newdef ClusterDecayer:ClSmrCharm 0.163
newdef HadronSelector:SingleHadronLimitCharm 0.000
# Clustering parameters for exotic quarks
# (e.g. hadronizing Susy particles)
newdef ClusterFissioner:ClMaxExotic 2.7*GeV
newdef ClusterFissioner:ClPowExotic 1.46
newdef ClusterFissioner:PSplitExotic 1.00
newdef ClusterDecayer:ClDirExotic 1
newdef ClusterDecayer:ClSmrExotic 0.
newdef HadronSelector:SingleHadronLimitExotic 0.
#
newdef PartonSplitter:SplitPwtSquark 0.824135
newdef PartonSplitter:Split uds
#
newdef HadronSelector:PwtDquark 1.0
newdef HadronSelector:PwtUquark 1.0
newdef HadronSelector:PwtSquark 0.291717
newdef HadronSelector:PwtCquark 0.0
newdef HadronSelector:PwtBquark 0.0
newdef HadronSelector:PwtDIquarkS0 0.2
newdef HadronSelector:PwtDIquarkS1 0.1
newdef HadronSelector:SngWt 0.75
-newdef HadronSelector:OctWt 0.60
newdef HadronSelector:DecWt 0.50
newdef HadronSelector:Mode 1
newdef HadronSelector:BelowThreshold All
newdef HadronSelector:scHadronWtFactor 2.5
newdef HadronSelector:sbHadronWtFactor 2.5
create Herwig::SpinHadronizer SpinHadronizer