diff --git a/Hadronization/HadronSpectrum.cc b/Hadronization/HadronSpectrum.cc
--- a/Hadronization/HadronSpectrum.cc
+++ b/Hadronization/HadronSpectrum.cc
@@ -1,609 +1,609 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the HadronSpectrum class.
//
#include "HadronSpectrum.h"
#include "ClusterHadronizationHandler.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include <ThePEG/Repository/CurrentGenerator.h>
#include "Herwig/Utilities/Kinematics.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
namespace {
// debug helper
void dumpTable(const HadronSpectrum::HadronTable & tbl) {
typedef HadronSpectrum::HadronTable::const_iterator TableIter;
for (TableIter it = tbl.begin(); it != tbl.end(); ++it) {
cerr << it->first.first << ' '
<< it->first.second << '\n';
for (HadronSpectrum::KupcoData::const_iterator jt = it->second.begin();
jt != it->second.end(); ++jt) {
cerr << '\t' << *jt << '\n';
}
}
}
}
HadronSpectrum::HadronSpectrum()
: Interfaced(),
belowThreshold_(0),
_repwt(Lmax,vector<vector<double> >(Jmax,vector<double>(Nmax))) {}
HadronSpectrum::~HadronSpectrum() {}
void HadronSpectrum::doinit() {
Interfaced::doinit();
// construct the hadron tables
constructHadronTable();
// lightest members (hadrons)
for(const PDPtr & p1 : partons()) {
for(const PDPtr & p2 : partons()) {
tcPDPair lp = lightestHadronPair(p1,p2);
if(lp.first && lp.second)
lightestHadrons_[make_pair(p1->id(),p2->id())] = lp;
}
}
// for debugging
- if(Debug::level >= 10 )
+ if (Debug::level >= 10)
dumpTable(table());
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void HadronSpectrum::persistentOutput(PersistentOStream & os) const {
os << _table << _partons << _forbidden
<< belowThreshold_ << _repwt << _pwt << lightestHadrons_;
}
void HadronSpectrum::persistentInput(PersistentIStream & is, int) {
is >> _table >> _partons >> _forbidden
>> belowThreshold_ >> _repwt >> _pwt >> lightestHadrons_;
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeAbstractClass<HadronSpectrum,Interfaced>
describeHerwigHadronSpectrum("Herwig::HadronSpectrum", "Herwig.so");
void HadronSpectrum::Init() {
static ClassDocumentation<HadronSpectrum> documentation
("There is no documentation for the HadronSpectrum class");
static RefVector<HadronSpectrum,ParticleData> interfacePartons
("Partons",
"The partons which are to be considered as the consistuents of the hadrons.",
&HadronSpectrum::_partons, -1, false, false, true, false, false);
static RefVector<HadronSpectrum,ParticleData> interfaceForbidden
("Forbidden",
"The PDG codes of the particles which cannot be produced in the hadronization.",
&HadronSpectrum::_forbidden, -1, false, false, true, false, false);
}
void HadronSpectrum::insertToHadronTable(tPDPtr &particle, int flav1, int flav2) {
// inserting a new Hadron in the hadron table.
long pid = particle->id();
int pspin = particle->iSpin();
HadronInfo a(pid, particle,specialWeight(pid),particle->mass());
// set the weight to the number of spin states
a.overallWeight = pspin*a.swtef;
// mesons
if(pspin%2==1) insertMeson(a,flav1,flav2);
// spin-1/2 baryons
else if(pspin==2) insertOneHalf(a,flav1,flav2);
// spin -3/2 baryons
else if(pspin==4) insertThreeHalf(a,flav1,flav2);
// all other cases
else {
assert(false);
}
}
void HadronSpectrum::insertOneHalf(HadronInfo a, int flav1, int flav2) {
assert(DiquarkMatcher::Check(flav1));
long iq1 = flav1/1000;
long iq2 = (flav1/100)%10;
if(iq1!=iq2 && flav1%10==3) flav1-=2;
if(iq1==iq2) {
if(iq1==flav2) {
a.overallWeight *= 1.5;
_table[make_pair(flav1,flav2)].insert(a);
_table[make_pair(flav2,flav1)].insert(a);
}
else {
_table[make_pair(flav1,flav2)].insert(a);
_table[make_pair(flav2,flav1)].insert(a);
long f3 = makeDiquarkID(iq1,flav2,1);
_table[make_pair(iq1,f3 )].insert(a);
_table[make_pair(f3 ,iq1)].insert(a);
}
}
else if(iq1==flav2) {
// ud1 u type
_table[make_pair(flav1,flav2)].insert(a);
_table[make_pair(flav2,flav1)].insert(a);
// and uu1 d type
long f3 = makeDiquarkID(iq1,iq1,3);
a.overallWeight *= a.wt;
_table[make_pair(f3 ,iq2)].insert(a);
_table[make_pair(iq2, f3)].insert(a);
}
else if(iq2==flav2) assert(false);
else {
_table[make_pair(flav1,flav2)].insert(a);
_table[make_pair(flav2,flav1)].insert(a);
long f3 = makeDiquarkID(iq1,flav2,1);
_table[make_pair(iq2,f3)].insert(a);
_table[make_pair(f3,iq2)].insert(a);
// 3rd perm
f3 = makeDiquarkID(iq2,flav2,1);
_table[make_pair(iq1,f3)].insert(a);
_table[make_pair(f3,iq1)].insert(a);
}
}
void HadronSpectrum::insertThreeHalf(HadronInfo a, int flav1, int flav2) {
assert(DiquarkMatcher::Check(flav1));
long iq1 = flav1/1000;
long iq2 = (flav1/100)%10;
if(iq1!=iq2 && flav1%10==3) flav1-=2;
if(iq1==iq2) {
if(iq1==flav2) {
a.overallWeight *= 1.5;
_table[make_pair(flav1,flav2)].insert(a);
_table[make_pair(flav2,flav1)].insert(a);
}
else {
_table[make_pair(flav1,flav2)].insert(a);
_table[make_pair(flav2,flav1)].insert(a);
long f3 = makeDiquarkID(iq1,flav2,1);
_table[make_pair(iq1,f3 )].insert(a);
_table[make_pair(f3 ,iq1)].insert(a);
}
}
else if(iq1==flav2) {
// ud1 u type
_table[make_pair(flav1,flav2)].insert(a);
_table[make_pair(flav2,flav1)].insert(a);
// and uu1 d type
long f3 = makeDiquarkID(iq1,iq1,3);
a.overallWeight *= a.wt;
_table[make_pair(f3 ,iq2)].insert(a);
_table[make_pair(iq2, f3)].insert(a);
}
else {
_table[make_pair(flav1,flav2)].insert(a);
_table[make_pair(flav2,flav1)].insert(a);
long f3 = makeDiquarkID(iq1,flav2,1);
_table[make_pair(iq2,f3)].insert(a);
_table[make_pair(f3,iq2)].insert(a);
// 3rd perm
f3 = makeDiquarkID(iq2,flav2,1);
_table[make_pair(iq1,f3)].insert(a);
_table[make_pair(f3,iq1)].insert(a);
}
}
tcPDPtr HadronSpectrum::chooseSingleHadron(tcPDPtr par1, tcPDPtr par2,
Energy mass) const {
Energy threshold = hadronPairThreshold(par1,par2);
// only do one hadron decay if mass less than the threshold
if(mass>=threshold) return tcPDPtr();
// select the hadron
tcPDPtr hadron;
// old option pick the lightest hadron
if(belowThreshold_ == 0) {
hadron= lightestHadron(par1,par2);
}
// new option select from those available
else if(belowThreshold_ == 1) {
vector<pair<tcPDPtr,double> > hadrons =
hadronsBelowThreshold(threshold,par1,par2);
if(hadrons.size()==1) {
hadron = hadrons[0].first;
}
else if(hadrons.empty()) {
hadron= lightestHadron(par1,par2);
}
else {
double totalWeight=0.;
for(unsigned int ix=0;ix<hadrons.size();++ix) {
totalWeight += hadrons[ix].second;
}
totalWeight *= UseRandom::rnd();
for(unsigned int ix=0;ix<hadrons.size();++ix) {
if(totalWeight<=hadrons[ix].second) {
hadron = hadrons[ix].first;
break;
}
else
totalWeight -= hadrons[ix].second;
}
assert(hadron);
}
}
else
assert(false);
return hadron;
}
tcPDPair HadronSpectrum::chooseHadronPair(const Energy cluMass,
tcPDPtr par1, tcPDPtr par2) const {
useMe();
// if either of the input partons is a diquark don't allow diquarks to be
// produced
bool diquark0 = !(DiquarkMatcher::Check(par1->id()) || DiquarkMatcher::Check(par2->id()));
bool diquark1 = diquark0;
bool quark = true;
// decide is baryon or meson production
if(diquark0) std::tie(quark,diquark0,diquark1) = selectBaryon(cluMass,par1,par2);
// weights for the different possibilities
Energy weight, wgtsum(ZERO);
// loop over all hadron pairs with the allowed flavours
static vector<Kupco> hadrons;
hadrons.clear();
for(unsigned int ix=0;ix<partons().size();++ix) {
tcPDPtr quarktopick = partons()[ix];
if(!quark && std::find(hadronizingQuarks().begin(), hadronizingQuarks().end(),
abs(quarktopick->id())) != hadronizingQuarks().end()) continue;
if(DiquarkMatcher::Check(quarktopick->id()) &&
((!diquark0 && quarktopick->iSpin()==1) ||
(!diquark1 && quarktopick->iSpin()==3))) continue;
HadronTable::const_iterator
tit1 = table().find(make_pair(abs(par1->id()),quarktopick->id()));
HadronTable::const_iterator
tit2 = table().find(make_pair(quarktopick->id(),abs(par2->id())));
// If not in table skip
if(tit1 == table().end()||tit2==table().end()) continue;
// tables empty skip
const KupcoData & T1 = tit1->second;
const KupcoData & T2 = tit2->second;
if(T1.empty()||T2.empty()) continue;
// if too massive skip
if(cluMass <= T1.begin()->mass +
T2.begin()->mass) continue;
// quark weight
double quarkWeight = pwt(quarktopick->id());
quarkWeight = specialQuarkWeight(quarkWeight,quarktopick->id(),
cluMass,par1,par2);
// loop over the hadrons
KupcoData::const_iterator H1,H2;
for(H1 = T1.begin();H1 != T1.end(); ++H1) {
for(H2 = T2.begin();H2 != T2.end(); ++H2) {
// break if cluster too light
if(cluMass < H1->mass + H2->mass) break;
weight = quarkWeight * H1->overallWeight * H2->overallWeight *
Kinematics::pstarTwoBodyDecay(cluMass, H1->mass, H2->mass);
int signQ = 0;
assert (par1 && quarktopick);
assert (par2);
assert(quarktopick->CC());
if(canBeHadron(par1, quarktopick->CC())
&& canBeHadron(quarktopick, par2))
signQ = +1;
else if(canBeHadron(par1, quarktopick)
&& canBeHadron(quarktopick->CC(), par2))
signQ = -1;
else {
cerr << "Could not make sign for" << par1->id()<< " " << quarktopick->id()
<< " " << par2->id() << "\n";
assert(false);
}
if (signQ == -1)
quarktopick = quarktopick->CC();
// construct the object with the info
Kupco a(quarktopick, H1->ptrData, H2->ptrData, weight);
hadrons.push_back(a);
wgtsum += weight;
}
}
}
if (hadrons.empty())
return make_pair(tcPDPtr(),tcPDPtr());
// select the hadron
wgtsum *= UseRandom::rnd();
unsigned int ix=0;
do {
wgtsum-= hadrons[ix].weight;
++ix;
}
while(wgtsum > ZERO && ix < hadrons.size());
if(ix == hadrons.size() && wgtsum > ZERO)
return make_pair(tcPDPtr(),tcPDPtr());
--ix;
assert(hadrons[ix].idQ);
int signHad1 = signHadron(par1, hadrons[ix].idQ->CC(), hadrons[ix].hadron1);
int signHad2 = signHadron(par2, hadrons[ix].idQ, hadrons[ix].hadron2);
assert( signHad1 != 0 && signHad2 != 0 );
return make_pair
( signHad1 > 0 ? hadrons[ix].hadron1 : tcPDPtr(hadrons[ix].hadron1->CC()),
signHad2 > 0 ? hadrons[ix].hadron2 : tcPDPtr(hadrons[ix].hadron2->CC()));
}
std::tuple<bool,bool,bool> HadronSpectrum::selectBaryon(const Energy, tcPDPtr, tcPDPtr ) const {
assert(false);
}
tcPDPair HadronSpectrum::lightestHadronPair(tcPDPtr ptr1, tcPDPtr ptr2) const {
Energy currentSum = Constants::MaxEnergy;
tcPDPair output;
for(unsigned int ix=0; ix<partons().size(); ++ix) {
HadronTable::const_iterator
tit1=table().find(make_pair(abs(ptr1->id()),partons()[ix]->id())),
tit2=table().find(make_pair(partons()[ix]->id(),abs(ptr2->id())));
if( tit1==table().end() || tit2==table().end()) continue;
if(tit1->second.empty()||tit2->second.empty()) continue;
Energy s = tit1->second.begin()->mass + tit2->second.begin()->mass;
if(currentSum > s) {
currentSum = s;
output.first = tit1->second.begin()->ptrData;
output.second = tit2->second.begin()->ptrData;
}
}
return output;
}
tcPDPtr HadronSpectrum::lightestHadron(tcPDPtr ptr1, tcPDPtr ptr2) const {
assert(ptr1 && ptr2);
// find entry in the table
pair<long,long> ids = make_pair(abs(ptr1->id()),abs(ptr2->id()));
HadronTable::const_iterator tit=_table.find(ids);
// throw exception if flavours wrong
if (tit==_table.end())
throw Exception() << "Could not find "
<< ids.first << ' ' << ids.second
<< " in _table. "
<< "In HadronSpectrum::lightestHadron()"
<< Exception::eventerror;
if(tit->second.empty())
throw Exception() << "HadronSpectrum::lightestHadron "
<< "could not find any hadrons containing "
<< ptr1->id() << ' ' << ptr2->id() << '\n'
<< tit->first.first << ' '
<< tit->first.second << Exception::eventerror;
// find the lightest hadron
int sign = signHadron(ptr1,ptr2,tit->second.begin()->ptrData);
tcPDPtr candidate = sign > 0 ?
tit->second.begin()->ptrData : tit->second.begin()->ptrData->CC();
// \todo 20 GeV limit is temporary fudge to let SM particles go through.
// \todo Use isExotic instead?
if (candidate->mass() > 20*GeV
&& candidate->mass() < ptr1->constituentMass() + ptr2->constituentMass()) {
generator()->log() << "HadronSpectrum::lightestHadron: "
<< "chosen candidate " << candidate->PDGName()
<< " is lighter than its constituents "
<< ptr1->PDGName() << ", " << ptr2->PDGName() << '\n'
<< candidate->mass()/GeV << " < " << ptr1->constituentMass()/GeV
<< " + " << ptr2->constituentMass()/GeV << '\n'
<< "Check your particle data tables.\n";
assert(false);
}
return candidate;
}
vector<pair<tcPDPtr,double> >
HadronSpectrum::hadronsBelowThreshold(Energy threshold, tcPDPtr ptr1,
tcPDPtr ptr2) const {
assert(ptr1 && ptr2);
// find entry in the table
pair<long,long> ids = make_pair(abs(ptr1->id()),abs(ptr2->id()));
HadronTable::const_iterator tit=_table.find(ids);
// throw exception if flavours wrong
if (tit==_table.end())
throw Exception() << "Could not find "
<< ids.first << ' ' << ids.second
<< " in _table. "
<< "In HadronSpectrum::hadronsBelowThreshold()"
<< Exception::eventerror;
if(tit->second.empty())
throw Exception() << "HadronSpectrum::hadronsBelowThreshold() "
<< "could not find any hadrons containing "
<< ptr1->id() << ' ' << ptr2->id() << '\n'
<< tit->first.first << ' '
<< tit->first.second << Exception::eventerror;
vector<pair<tcPDPtr,double> > candidates;
KupcoData::const_iterator hit = tit->second.begin();
// find the hadrons
while(hit!=tit->second.end()&&hit->mass<threshold) {
// find the hadron
int sign = signHadron(ptr1,ptr2,hit->ptrData);
tcPDPtr candidate = sign > 0 ? hit->ptrData : hit->ptrData->CC();
// \todo 20 GeV limit is temporary fudge to let SM particles go through.
// \todo Use isExotic instead?
if (candidate->mass() > 20*GeV
&& candidate->mass() < ptr1->constituentMass() + ptr2->constituentMass()) {
generator()->log() << "HadronSpectrum::hadronsBelowTheshold: "
<< "chosen candidate " << candidate->PDGName()
<< " is lighter than its constituents "
<< ptr1->PDGName() << ", " << ptr2->PDGName() << '\n'
<< candidate->mass()/GeV << " < " << ptr1->constituentMass()/GeV
<< " + " << ptr2->constituentMass()/GeV << '\n'
<< "Check your particle data tables.\n";
assert(false);
}
candidates.push_back(make_pair(candidate,hit->overallWeight));
++hit;
}
return candidates;
}
Energy HadronSpectrum::massLightestBaryonPair(tcPDPtr ptr1, tcPDPtr ptr2) const {
// Make sure that we don't have any diquarks as input, return arbitrarily
// large value if we do
Energy currentSum = Constants::MaxEnergy;
for(unsigned int ix=0; ix<_partons.size(); ++ix) {
if(!DiquarkMatcher::Check(_partons[ix]->id())) continue;
HadronTable::const_iterator
tit1=_table.find(make_pair(abs(ptr1->id()),_partons[ix]->id())),
tit2=_table.find(make_pair(_partons[ix]->id(),abs(ptr2->id())));
if( tit1==_table.end() || tit2==_table.end()) continue;
if(tit1->second.empty()||tit2->second.empty()) continue;
Energy s = tit1->second.begin()->mass + tit2->second.begin()->mass;
if(currentSum > s) currentSum = s;
}
return currentSum;
}
double HadronSpectrum::mesonWeight(long id) const {
// Total angular momentum
int j = ((id % 10) - 1) / 2;
// related to Orbital angular momentum l
int nl = (id/10000 )%10;
int l = -999;
int n = (id/100000)%10; // Radial excitation
if(j == 0) l = nl;
else if(nl == 0) l = j - 1;
else if(nl == 1 || nl == 2) l = j;
else if(nl == 3) l = j + 1;
// Angular or Radial excited meson
if((l||j||n) && l>=0 && l<Lmax && j<Jmax && n<Nmax) {
return sqr(_repwt[l][j][n]);
}
// rest is not excited or
// has spin >= 5/2 (ispin >= 6), haven't got those
else
return 1.0;
}
int HadronSpectrum::signHadron(tcPDPtr idQ1, tcPDPtr idQ2,
tcPDPtr hadron) const {
// This method receives in input three PDG ids, whose the
// first two have proper signs (corresponding to particles, id > 0,
// or antiparticles, id < 0 ), whereas the third one must
// be always positive (particle not antiparticle),
// corresponding to:
// --- quark-antiquark, or antiquark-quark, or
// quark-diquark, or diquark-quark, or
// antiquark-antidiquark, or antidiquark-antiquark
// for the first two input (idQ1, idQ2);
// --- meson or baryon for the third input (idHad):
// The method returns:
// --- + 1 if the two partons (idQ1, idQ2) are exactly
// the constituents for the hadron idHad;
// --- - 1 if the two partons (idQ1, idQ2) are exactly
// the constituents for the anti-hadron -idHad;
// --- + 0 otherwise.
// The method it is therefore useful to decide the
// sign of the id of the produced hadron as appeared
// in the vector _vecHad (where only hadron idHad > 0 are present)
// given the two constituent partons.
int sign = 0;
long idHad = hadron->id();
assert(idHad > 0);
int chargeIn = idQ1->iCharge() + idQ2->iCharge();
int chargeOut = hadron->iCharge();
// same charge
if( chargeIn == chargeOut && chargeIn !=0 ) sign = +1;
else if(chargeIn == -chargeOut && chargeIn !=0 ) sign = -1;
else if(chargeIn == 0 && chargeOut == 0 ) {
// In the case of same null charge, there are four cases:
// i) K0-like mesons, B0-like mesons, Bs-like mesons
// the PDG convention is to consider them "antiparticle" (idHad < 0)
// if the "dominant" (heavier) flavour (respectively, s, b)
// is a quark (idQ > 0): for instance, B0s = (b, sbar) has id < 0
// Remember that there is an important exception for K0L (id=130) and
// K0S (id=310): they don't have antiparticles, therefore idHad > 0
// always. We use below the fact that K0L and K0S are the unique
// hadrons having 0 the first (less significant) digit of their id.
// 2) D0-like mesons: the PDG convention is to consider them "particle"
// (idHad > 0) if the charm flavour is carried by a c: (c,ubar) has id>0
// 3) the remaining mesons should not have antiparticle, therefore their
// sign is always positive.
// 4) for baryons, that is when one of idQ1 and idQ2 is a (anti-) quark and
// the other one is a (anti-) diquark the sign is negative when both
// constituents are "anti", that is both with id < 0; positive otherwise.
// meson
if(std::find(hadronizingQuarks().begin(), hadronizingQuarks().end(),
abs(idQ1->id())) != hadronizingQuarks().end() &&
std::find(hadronizingQuarks().begin(), hadronizingQuarks().end(),
abs(idQ2->id())) != hadronizingQuarks().end())
{
int idQa = abs(idQ1->id());
int idQb = abs(idQ2->id());
int dominant = idQ2->id();
if(idQa > idQb) {
swap(idQa,idQb);
dominant = idQ1->id();
}
if((idQa==ParticleID::d && idQb==ParticleID::s) ||
(idQa==ParticleID::d && idQb==ParticleID::b) ||
(idQa==ParticleID::s && idQb==ParticleID::b)) {
// idHad%10 is zero for K0L,K0S
if (dominant < 0 || idHad%10 == 0) sign = +1;
else if(dominant > 0) sign = -1;
}
else if((idQa==ParticleID::u && idQb==ParticleID::c) ||
(idQa==ParticleID::u && idQb==ParticleID::t) ||
(idQa==ParticleID::c && idQb==ParticleID::t)) {
if (dominant > 0) sign = +1;
else if(dominant < 0) sign = -1;
}
else if(idQa==idQb) sign = +1;
// sets sign for Susy particles
else sign = (dominant > 0) ? +1 : -1;
}
// baryon
else if(DiquarkMatcher::Check(idQ1->id()) || DiquarkMatcher::Check(idQ2->id())) {
if (idQ1->id() > 0 && idQ2->id() > 0) sign = +1;
else if(idQ1->id() < 0 && idQ2->id() < 0) sign = -1;
}
}
if (sign == 0) {
cerr << "Could not work out sign for "
<< idQ1->PDGName() << ' '
<< idQ2->PDGName() << " => "
<< hadron->PDGName() << '\n';
assert(false);
}
return sign;
}
/*
PDPtr HadronSpectrum::makeDiquark(tcPDPtr par1, tcPDPtr par2) const {
long id1 = par1->id();
long id2 = par2->id();
long pspin = id1==id2 ? 3 : 1;
long idnew = makeDiquarkID(id1,id2, pspin);
return getParticleData(idnew);
}
*/
bool HadronSpectrum::canBeMeson(tcPDPtr par1,tcPDPtr par2) const {
assert(par1 && par2);
long id1 = par1->id();
long id2 = par2->id();
// a Meson must not have any diquarks
if(DiquarkMatcher::Check(id1) || DiquarkMatcher::Check(id2)) return false;
return (std::find(hadronizingQuarks().begin(), hadronizingQuarks().end(),
abs(id1)) != hadronizingQuarks().end() &&
std::find(hadronizingQuarks().begin(), hadronizingQuarks().end(),
abs(id2)) != hadronizingQuarks().end() &&
id1*id2 < 0);
}
diff --git a/Hadronization/StandardModelHadronSpectrum.cc b/Hadronization/StandardModelHadronSpectrum.cc
--- a/Hadronization/StandardModelHadronSpectrum.cc
+++ b/Hadronization/StandardModelHadronSpectrum.cc
@@ -1,696 +1,695 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the StandardModelHadronSpectrum class.
//
#include "StandardModelHadronSpectrum.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include <ThePEG/PDT/EnumParticles.h>
#include <ThePEG/Repository/EventGenerator.h>
#include <ThePEG/Repository/Repository.h>
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
namespace {
bool weightIsLess (pair<long,double> a, pair<long,double> b) {
return a.second < b.second;
}
/**
* Return true if the particle pointer corresponds to a diquark
* or anti-diquark carrying b flavour; false otherwise.
*/
inline bool isDiquarkWithB(tcPDPtr par1) {
if (!par1) return false;
long id1 = par1->id();
return DiquarkMatcher::Check(id1) && (abs(id1)/1000)%10 == ParticleID::b;
}
/**
* Return true if the particle pointer corresponds to a diquark
* or anti-diquark carrying c flavour; false otherwise.
*/
inline bool isDiquarkWithC(tcPDPtr par1) {
if (!par1) return false;
long id1 = par1->id();
return ( DiquarkMatcher::Check(id1) &&
( (abs(id1)/1000)%10 == ParticleID::c
|| (abs(id1)/100)%10 == ParticleID::c ) );
}
}
StandardModelHadronSpectrum::StandardModelHadronSpectrum(unsigned int opt)
: HadronSpectrum(),
_pwtDquark( 1.0 ),_pwtUquark( 1.0 ),_pwtSquark( 1.0 ),_pwtCquark( 0.0 ),
_pwtBquark( 0.0 ),
_sngWt( 1.0 ),_decWt( 1.0 ),
_weight1S0(Nmax,1.),_weight3S1(Nmax,1.),_weight1P1(Nmax,1.),_weight3P0(Nmax,1.),
_weight3P1(Nmax,1.),_weight3P2(Nmax,1.),_weight1D2(Nmax,1.),_weight3D1(Nmax,1.),
_weight3D2(Nmax,1.),_weight3D3(Nmax,1.),
_topt(opt),_trial(0),
_limBottom(), _limCharm(), _limExotic()
{
// The mixing angles
// the ideal mixing angle
const double idealAngleMix = atan( sqrt(0.5) ) * 180.0 / Constants::pi;
// \eta-\eta' mixing angle
_etamix = -23.0;
// phi-omega mixing angle
_phimix = +36.0;
// h_1'-h_1 mixing angle
_h1mix = idealAngleMix;
// f_0(1710)-f_0(1370) mixing angle
_f0mix = idealAngleMix;
// f_1(1420)-f_1(1285)\f$ mixing angle
_f1mix = idealAngleMix;
// f'_2-f_2\f$ mixing angle
_f2mix = +26.0;
// eta_2(1870)-eta_2(1645) mixing angle
_eta2mix = idealAngleMix;
// phi(???)-omega(1650) mixing angle
_omhmix = idealAngleMix;
// phi_3-omega_3 mixing angle
_ph3mix = +28.0;
// eta(1475)-eta(1295) mixing angle
_eta2Smix = idealAngleMix;
// phi(1680)-omega(1420) mixing angle
_phi2Smix = idealAngleMix;
}
StandardModelHadronSpectrum::~StandardModelHadronSpectrum() {}
void StandardModelHadronSpectrum::persistentOutput(PersistentOStream & os) const {
os << _pwtDquark << _pwtUquark << _pwtSquark
<< _pwtCquark << _pwtBquark
<< _etamix << _phimix << _h1mix << _f0mix << _f1mix << _f2mix
<< _eta2mix << _omhmix << _ph3mix << _eta2Smix << _phi2Smix
<< _weight1S0 << _weight3S1 << _weight1P1 << _weight3P0 << _weight3P1
<< _weight3P2 << _weight1D2 << _weight3D1 << _weight3D2 << _weight3D3
<< _sngWt << _decWt << _repwt
<< _limBottom << _limCharm << _limExotic;
}
void StandardModelHadronSpectrum::persistentInput(PersistentIStream & is, int) {
is >> _pwtDquark >> _pwtUquark >> _pwtSquark
>> _pwtCquark >> _pwtBquark
>> _etamix >> _phimix >> _h1mix >> _f0mix >> _f1mix >> _f2mix
>> _eta2mix >> _omhmix >> _ph3mix >> _eta2Smix >> _phi2Smix
>> _weight1S0 >> _weight3S1 >> _weight1P1 >> _weight3P0 >> _weight3P1
>> _weight3P2 >> _weight1D2 >> _weight3D1 >> _weight3D2 >> _weight3D3
>> _sngWt >> _decWt >> _repwt
>> _limBottom >> _limCharm >> _limExotic;
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeAbstractClass<StandardModelHadronSpectrum,HadronSpectrum>
describeHerwigStandardModelHadronSpectrum("Herwig::StandardModelHadronSpectrum", "Herwig.so");
void StandardModelHadronSpectrum::Init() {
static ClassDocumentation<StandardModelHadronSpectrum> documentation
("There is no documentation for the StandardModelHadronSpectrum class");
static Parameter<StandardModelHadronSpectrum,double>
interfacePwtDquark("PwtDquark","Weight for choosing a quark D",
&StandardModelHadronSpectrum::_pwtDquark, 0, 1.0, 0.0, 10.0,
false,false,false);
static Parameter<StandardModelHadronSpectrum,double>
interfacePwtUquark("PwtUquark","Weight for choosing a quark U",
&StandardModelHadronSpectrum::_pwtUquark, 0, 1.0, 0.0, 10.0,
false,false,false);
static Parameter<StandardModelHadronSpectrum,double>
interfacePwtSquark("PwtSquark","Weight for choosing a quark S",
&StandardModelHadronSpectrum::_pwtSquark, 0, 1.0, 0.0, 10.0,
false,false,false);
static Parameter<StandardModelHadronSpectrum,double>
interfacePwtCquark("PwtCquark","Weight for choosing a quark C",
&StandardModelHadronSpectrum::_pwtCquark, 0, 0.0, 0.0, 10.0,
false,false,false);
static Parameter<StandardModelHadronSpectrum,double>
interfacePwtBquark("PwtBquark","Weight for choosing a quark B",
&StandardModelHadronSpectrum::_pwtBquark, 0, 0.0, 0.0, 10.0,
false,false,false);
static Parameter<StandardModelHadronSpectrum,double>
interfaceSngWt("SngWt","Weight for singlet baryons",
&StandardModelHadronSpectrum::_sngWt, 0, 1.0, 0.0, 10.0,
false,false,false);
static Parameter<StandardModelHadronSpectrum,double>
interfaceDecWt("DecWt","Weight for decuplet baryons",
&StandardModelHadronSpectrum::_decWt, 0, 1.0, 0.0, 10.0,
false,false,false);
//
// mixing angles
//
// the ideal mixing angle
const double idealAngleMix = atan( sqrt(0.5) ) * 180.0 / Constants::pi;
static Parameter<StandardModelHadronSpectrum,double> interface11S0Mixing
("11S0Mixing",
"The mixing angle for the I=0 mesons from the 1 1S0 multiplet,"
" i.e. eta and etaprime.",
&StandardModelHadronSpectrum::_etamix, -23., -180., 180.,
false, false, Interface::limited);
static Parameter<StandardModelHadronSpectrum,double> interface13S1Mixing
("13S1Mixing",
"The mixing angle for the I=0 mesons from the 1 3S1 multiplet,"
" i.e. phi and omega.",
&StandardModelHadronSpectrum::_phimix, +36., -180., 180.,
false, false, Interface::limited);
static Parameter<StandardModelHadronSpectrum,double> interface11P1Mixing
("11P1Mixing",
"The mixing angle for the I=0 mesons from the 1 1P1 multiplet,"
" i.e. h_1' and h_1.",
&StandardModelHadronSpectrum::_h1mix, idealAngleMix, -180., 180.,
false, false, Interface::limited);
static Parameter<StandardModelHadronSpectrum,double> interface13P0Mixing
("13P0Mixing",
"The mixing angle for the I=0 mesons from the 1 3P0 multiplet,"
" i.e. f_0(1710) and f_0(1370).",
&StandardModelHadronSpectrum::_f0mix, idealAngleMix, -180., 180.,
false, false, Interface::limited);
static Parameter<StandardModelHadronSpectrum,double> interface13P1Mixing
("13P1Mixing",
"The mixing angle for the I=0 mesons from the 1 3P1 multiplet,"
" i.e. f_1(1420) and f_1(1285).",
&StandardModelHadronSpectrum::_f1mix, idealAngleMix, -180., 180.,
false, false, Interface::limited);
static Parameter<StandardModelHadronSpectrum,double> interface13P2Mixing
("13P2Mixing",
"The mixing angle for the I=0 mesons from the 1 3P2 multiplet,"
" i.e. f'_2 and f_2.",
&StandardModelHadronSpectrum::_f2mix, 26.0, -180., 180.,
false, false, Interface::limited);
static Parameter<StandardModelHadronSpectrum,double> interface11D2Mixing
("11D2Mixing",
"The mixing angle for the I=0 mesons from the 1 1D2 multiplet,"
" i.e. eta_2(1870) and eta_2(1645).",
&StandardModelHadronSpectrum::_eta2mix, idealAngleMix, -180., 180.,
false, false, Interface::limited);
static Parameter<StandardModelHadronSpectrum,double> interface13D0Mixing
("13D0Mixing",
"The mixing angle for the I=0 mesons from the 1 3D0 multiplet,"
" i.e. eta_2(1870) phi(?) and omega(1650).",
&StandardModelHadronSpectrum::_omhmix, idealAngleMix, -180., 180.,
false, false, Interface::limited);
static Parameter<StandardModelHadronSpectrum,double> interface13D1Mixing
("13D1Mixing",
"The mixing angle for the I=0 mesons from the 1 3D1 multiplet,"
" i.e. phi_3 and omega_3.",
&StandardModelHadronSpectrum::_ph3mix, 28.0, -180., 180.,
false, false, Interface::limited);
static Parameter<StandardModelHadronSpectrum,double> interface21S0Mixing
("21S0Mixing",
"The mixing angle for the I=0 mesons from the 2 1S0 multiplet,"
" i.e. eta(1475) and eta(1295).",
&StandardModelHadronSpectrum::_eta2Smix, idealAngleMix, -180., 180.,
false, false, Interface::limited);
static Parameter<StandardModelHadronSpectrum,double> interface23S1Mixing
("23S1Mixing",
"The mixing angle for the I=0 mesons from the 1 3S1 multiplet,"
" i.e. phi(1680) and omega(1420).",
&StandardModelHadronSpectrum::_phi2Smix, idealAngleMix, -180., 180.,
false, false, Interface::limited);
//
// the meson weights
//
static ParVector<StandardModelHadronSpectrum,double> interface1S0Weights
("1S0Weights",
"The weights for the 1S0 multiplets start with n=1.",
&StandardModelHadronSpectrum::_weight1S0, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<StandardModelHadronSpectrum,double> interface3S1Weights
("3S1Weights",
"The weights for the 3S1 multiplets start with n=1.",
&StandardModelHadronSpectrum::_weight3S1, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<StandardModelHadronSpectrum,double> interface1P1Weights
("1P1Weights",
"The weights for the 1P1 multiplets start with n=1.",
&StandardModelHadronSpectrum::_weight1P1, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<StandardModelHadronSpectrum,double> interface3P0Weights
("3P0Weights",
"The weights for the 3P0 multiplets start with n=1.",
&StandardModelHadronSpectrum::_weight3P0, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<StandardModelHadronSpectrum,double> interface3P1Weights
("3P1Weights",
"The weights for the 3P1 multiplets start with n=1.",
&StandardModelHadronSpectrum::_weight3P1, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<StandardModelHadronSpectrum,double> interface3P2Weights
("3P2Weights",
"The weights for the 3P2 multiplets start with n=1.",
&StandardModelHadronSpectrum::_weight3P2, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<StandardModelHadronSpectrum,double> interface1D2Weights
("1D2Weights",
"The weights for the 1D2 multiplets start with n=1.",
&StandardModelHadronSpectrum::_weight1D2, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<StandardModelHadronSpectrum,double> interface3D1Weights
("3D1Weights",
"The weights for the 3D1 multiplets start with n=1.",
&StandardModelHadronSpectrum::_weight3D1, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<StandardModelHadronSpectrum,double> interface3D2Weights
("3D2Weights",
"The weights for the 3D2 multiplets start with n=1.",
&StandardModelHadronSpectrum::_weight3D2, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<StandardModelHadronSpectrum,double> interface3D3Weights
("3D3Weights",
"The weights for the 3D3 multiplets start with n=1.",
&StandardModelHadronSpectrum::_weight3D3, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static Switch<StandardModelHadronSpectrum,unsigned int> interfaceTrial
("Trial",
"A Debugging option to only produce certain types of hadrons",
&StandardModelHadronSpectrum::_trial, 0, false, false);
static SwitchOption interfaceTrialAll
(interfaceTrial,
"All",
"Produce all the hadrons",
0);
static SwitchOption interfaceTrialPions
(interfaceTrial,
"Pions",
"Only produce pions",
1);
static SwitchOption interfaceTrialSpin2
(interfaceTrial,
"Spin2",
"Only mesons with spin less than or equal to two are produced",
2);
static SwitchOption interfaceTrialSpin3
(interfaceTrial,
"Spin3",
"Only hadrons with spin less than or equal to three are produced",
3);
static Parameter<StandardModelHadronSpectrum,double>
interfaceSingleHadronLimitBottom ("SingleHadronLimitBottom",
"Threshold for one-hadron decay of b-cluster",
&StandardModelHadronSpectrum::_limBottom,
0, 0.0, 0.0, 100.0,false,false,false);
static Parameter<StandardModelHadronSpectrum,double>
interfaceSingleHadronLimitCharm ("SingleHadronLimitCharm",
"threshold for one-hadron decay of c-cluster",
&StandardModelHadronSpectrum::_limCharm,
0, 0.0, 0.0, 100.0,false,false,false);
static Parameter<StandardModelHadronSpectrum,double>
interfaceSingleHadronLimitExotic ("SingleHadronLimitExotic",
"threshold for one-hadron decay of exotic cluster",
&StandardModelHadronSpectrum::_limExotic,
0, 0.0, 0.0, 100.0,false,false,false);
static Switch<StandardModelHadronSpectrum,unsigned int> interfaceBelowThreshold
("BelowThreshold",
"Option fo the selection of the hadrons if the cluster is below the pair threshold",
&StandardModelHadronSpectrum::belowThreshold_, 0, false, false);
static SwitchOption interfaceBelowThresholdLightest
(interfaceBelowThreshold,
"Lightest",
"Force cluster to decay to the lightest hadron with the appropriate flavours",
0);
static SwitchOption interfaceBelowThresholdAll
(interfaceBelowThreshold,
"All",
"Select from all the hadrons below the two hadron threshold according to their spin weights",
1);
}
PDPtr StandardModelHadronSpectrum::makeDiquark(tcPDPtr par1, tcPDPtr par2) const {
long id1 = par1->id();
long id2 = par2->id();
long pspin = id1==id2 ? 3 : 1;
long idnew = makeDiquarkID(id1,id2, pspin);
return getParticleData(idnew);
}
Energy StandardModelHadronSpectrum::hadronPairThreshold(tcPDPtr par1, tcPDPtr par2) const {
// Determine the sum of the nominal masses of the two lightest hadrons
// with the right flavour numbers as the cluster under consideration.
// Notice that we don't need real masses (drawn by a Breit-Wigner
// distribution) because the lightest pair of hadrons does not involve
// any broad resonance.
Energy threshold = massLightestHadronPair(par1,par2);
// Special: it allows one-hadron decays also above threshold.
if (isExotic(par1,par2))
threshold *= (1.0 + UseRandom::rnd()*_limExotic);
else if (hasBottom(par1,par2))
threshold *= (1.0 + UseRandom::rnd()*_limBottom);
else if (hasCharm(par1,par2))
threshold *= (1.0 + UseRandom::rnd()*_limCharm);
return threshold;
}
double StandardModelHadronSpectrum::mixingStateWeight(long id) const {
switch(id) {
case ParticleID::eta: return 0.5*probabilityMixing(_etamix ,1);
case ParticleID::etaprime: return 0.5*probabilityMixing(_etamix ,2);
case ParticleID::phi: return 0.5*probabilityMixing(_phimix ,1);
case ParticleID::omega: return 0.5*probabilityMixing(_phimix ,2);
case ParticleID::hprime_1: return 0.5*probabilityMixing(_h1mix ,1);
case ParticleID::h_1: return 0.5*probabilityMixing(_h1mix ,2);
case 10331: return 0.5*probabilityMixing(_f0mix ,1);
case 10221: return 0.5*probabilityMixing(_f0mix ,2);
case ParticleID::fprime_1: return 0.5*probabilityMixing(_f1mix ,1);
case ParticleID::f_1: return 0.5*probabilityMixing(_f1mix ,2);
case ParticleID::fprime_2: return 0.5*probabilityMixing(_f2mix ,1);
case ParticleID::f_2: return 0.5*probabilityMixing(_f2mix ,2);
case 10335: return 0.5*probabilityMixing(_eta2mix ,1);
case 10225: return 0.5*probabilityMixing(_eta2mix ,2);
// 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 StandardModelHadronSpectrum::doinit() {
- HadronSpectrum::doinit();
-
// set the weights for the various excited mesons
// set all to one to start with
for (int l = 0; l < Lmax; ++l ) {
for (int j = 0; j < Jmax; ++j) {
for (int n = 0; n < Nmax; ++n) {
_repwt[l][j][n] = 1.0;
}
}
}
// set the others from the relevant vectors
for( int ix=0;ix<max(int(_weight1S0.size()),int(Nmax));++ix)
_repwt[0][0][ix]=_weight1S0[ix];
for( int ix=0;ix<max(int(_weight3S1.size()),int(Nmax));++ix)
_repwt[0][1][ix]=_weight3S1[ix];
for( int ix=0;ix<max(int(_weight1P1.size()),int(Nmax));++ix)
_repwt[1][1][ix]=_weight1P1[ix];
for( int ix=0;ix<max(int(_weight3P0.size()),int(Nmax));++ix)
_repwt[1][0][ix]=_weight3P0[ix];
for( int ix=0;ix<max(int(_weight3P1.size()),int(Nmax));++ix)
_repwt[1][1][ix]=_weight3P1[ix];
for( int ix=0;ix<max(int(_weight3P2.size()),int(Nmax));++ix)
_repwt[1][2][ix]=_weight3P2[ix];
for( int ix=0;ix<max(int(_weight1D2.size()),int(Nmax));++ix)
_repwt[2][2][ix]=_weight1D2[ix];
for( int ix=0;ix<max(int(_weight3D1.size()),int(Nmax));++ix)
_repwt[2][1][ix]=_weight3D1[ix];
for( int ix=0;ix<max(int(_weight3D2.size()),int(Nmax));++ix)
_repwt[2][2][ix]=_weight3D2[ix];
for( int ix=0;ix<max(int(_weight3D3.size()),int(Nmax));++ix)
_repwt[2][3][ix]=_weight3D3[ix];
// find the maximum
map<long,double>::iterator pit =
max_element(_pwt.begin(),_pwt.end(),weightIsLess);
const double pmax = pit->second;
for(pit=_pwt.begin(); pit!=_pwt.end(); ++pit) {
pit->second/=pmax;
}
+ HadronSpectrum::doinit();
}
void StandardModelHadronSpectrum::constructHadronTable() {
// initialise the table
_table.clear();
for(unsigned int ix=0; ix<_partons.size(); ++ix) {
for(unsigned int iy=0; iy<_partons.size(); ++iy) {
if (!(DiquarkMatcher::Check(_partons[ix]->id())
&& DiquarkMatcher::Check(_partons[iy]->id())))
_table[make_pair(_partons[ix]->id(),_partons[iy]->id())] = KupcoData();
}
}
// get the particles from the event generator
ParticleMap particles = generator()->particles();
// loop over the particles
//double maxdd(0.),maxss(0.),maxrest(0.);
for(ParticleMap::iterator it=particles.begin();
it!=particles.end(); ++it) {
long pid = it->first;
tPDPtr particle = it->second;
int pspin = particle->iSpin();
// Don't include hadrons which are explicitly forbidden
if(find(_forbidden.begin(),_forbidden.end(),particle)!=_forbidden.end())
continue;
// Don't include non-hadrons or antiparticles
if(pid < 100) continue;
// remove diffractive particles
if(pspin == 0) continue;
// K_0S and K_0L not made make K0 and Kbar0
if(pid==ParticleID::K_S0||pid==ParticleID::K_L0) continue;
// Debugging options
// Only include those with 2J+1 less than...5
if(_trial==2 && pspin >= 5) continue;
// Only include those with 2J+1 less than...7
if(_trial==3 && pspin >= 7) continue;
// Only include pions
if(_trial==1 && pid!=111 && pid!=211) continue;
// shouldn't be coloured
if(particle->coloured()) continue;
// Get the flavours
const int x4 = (pid/1000)%10;
const int x3 = (pid/100 )%10;
const int x2 = (pid/10 )%10;
const int x7 = (pid/1000000)%10;
const bool wantSusy = x7 == 1 || x7 == 2;
// Skip non-hadrons (susy particles, etc...)
if(x3 == 0 || x2 == 0) continue;
// Skip particles which are neither SM nor SUSY
- if(x7 >= 3) continue;
+ if(x7 >= 3 && x7 != 9) continue;
int flav1,flav2;
// meson
if(x4 == 0) {
flav1 = x2;
flav2 = x3;
}
// baryon
else {
flav2 = x4;
// insert the spin 1 diquark, sort out the rest later
flav1 = makeDiquarkID(x2,x3,3);
}
if (wantSusy) flav2 += 1000000 * x7;
insertToHadronTable(particle,flav1,flav2);
}
// normalise the weights
if(_topt == 0) {
HadronTable::const_iterator tit;
KupcoData::iterator it;
for(tit=_table.begin();tit!=_table.end();++tit) {
double weight=0;
for(it = tit->second.begin(); it!=tit->second.end(); ++it)
weight=max(weight,it->overallWeight);
weight = 1./weight;
}
// double weight;
// if(tit->first.first==tit->first.second) {
// if(tit->first.first==1||tit->first.first==2) weight=1./maxdd;
// else if (tit->first.first==3) weight=1./maxss;
// else weight=1./maxrest;
// }
// else weight=1./maxrest;
// for(it = tit->second.begin(); it!=tit->second.end(); ++it) {
// it->rescale(weight);
// }
// }
}
}
double StandardModelHadronSpectrum::strangeWeight(const Energy, tcPDPtr, tcPDPtr) const {
assert(false);
}
void StandardModelHadronSpectrum::insertMeson(HadronInfo a, int flav1, int flav2) {
// identical light flavours
if(flav1 == flav2 && flav1<=3) {
// ddbar> uubar> admixture states
if(flav1==1) {
a.overallWeight *= 0.5;
_table[make_pair(1,1)].insert(a);
_table[make_pair(2,2)].insert(a);
}
// load up ssbar> uubar> ddbar> admixture states
else {
// uubar ddbar pieces
a.wt = mixingStateWeight(a.id);
a.overallWeight *= a.wt;
_table[make_pair(1,1)].insert(a);
_table[make_pair(2,2)].insert(a);
a.overallWeight /=a.wt;
// ssbar piece
a.wt = 1.- 2.*a.wt;
if(a.wt > 0) {
a.overallWeight *= a.wt;
_table[make_pair(3,3)].insert(a);
}
}
}
else {
_table[make_pair(flav1,flav2)].insert(a);
if(flav1 != flav2) _table[make_pair(flav2,flav1)].insert(a);
}
}
long StandardModelHadronSpectrum::makeDiquarkID(long id1, long id2, long pspin) const {
assert( id1 * id2 > 0
&& QuarkMatcher::Check(id1)
&& QuarkMatcher::Check(id2)) ;
long ida = abs(id1);
long idb = abs(id2);
if (ida < idb) swap(ida,idb);
if (pspin != 1 && pspin != 3) assert(false);
long idnew = ida*1000 + idb*100 + pspin;
// Diquarks made of quarks of the same type: uu, dd, ss, cc, bb,
// have spin 1, and therefore the less significant digit (which
// corresponds to 2*J+1) is 3 rather than 1 as all other Diquarks.
if (id1 == id2 && pspin == 1) {
//cerr<<"WARNING: spin-0 diquiark of the same type cannot exist."
// <<" Switching to spin-1 diquark.\n";
idnew = ida*1000 + idb*100 + 3;
}
return id1 > 0 ? idnew : -idnew;
}
bool StandardModelHadronSpectrum::hasBottom(tcPDPtr par1, tcPDPtr par2, tcPDPtr par3) const {
long id1 = par1 ? par1->id() : 0;
if ( !par2 && !par3 ) {
return
abs(id1) == ThePEG::ParticleID::b ||
isDiquarkWithB(par1) ||
( MesonMatcher::Check(id1)
&& (abs(id1)/100)%10 == ThePEG::ParticleID::b ) ||
( BaryonMatcher::Check(id1)
&& (abs(id1)/1000)%10 == ThePEG::ParticleID::b );
}
else {
long id2 = par2 ? par2->id() : 0;
long id3 = par3 ? par3->id() : 0;
return
abs(id1) == ThePEG::ParticleID::b || isDiquarkWithB(par1) ||
abs(id2) == ThePEG::ParticleID::b || isDiquarkWithB(par2) ||
abs(id3) == ThePEG::ParticleID::b || isDiquarkWithB(par3);
}
}
bool StandardModelHadronSpectrum::hasCharm(tcPDPtr par1, tcPDPtr par2, tcPDPtr par3) const {
long id1 = par1 ? par1->id(): 0;
if (!par2 && !par3) {
return
abs(id1) == ThePEG::ParticleID::c ||
isDiquarkWithC(par1) ||
( MesonMatcher::Check(id1) &&
((abs(id1)/100)%10 == ThePEG::ParticleID::c ||
(abs(id1)/10)%10 == ThePEG::ParticleID::c) ) ||
( BaryonMatcher::Check(id1) &&
((abs(id1)/1000)%10 == ThePEG::ParticleID::c ||
(abs(id1)/100)%10 == ThePEG::ParticleID::c ||
(abs(id1)/10)%10 == ThePEG::ParticleID::c) );
}
else {
long id2 = par2 ? par1->id(): 0;
long id3 = par3 ? par1->id(): 0;
return
abs(id1) == ThePEG::ParticleID::c || isDiquarkWithC(par1) ||
abs(id2) == ThePEG::ParticleID::c || isDiquarkWithC(par2) ||
abs(id3) == ThePEG::ParticleID::c || isDiquarkWithC(par3);
}
}
bool StandardModelHadronSpectrum::isExotic(tcPDPtr par1, tcPDPtr par2, tcPDPtr par3) const {
/// \todo make this more general
long id1 = par1 ? par1->id(): 0;
long id2 = par2 ? par2->id(): 0;
long id3 = par3 ? par3->id(): 0;
return
( (id1/1000000)% 10 != 0 && (id1/1000000)% 10 != 9 ) ||
( (id2/1000000)% 10 != 0 && (id2/1000000)% 10 != 9 ) ||
( (id3/1000000)% 10 != 0 && (id3/1000000)% 10 != 9 ) ||
abs(id1)==6||abs(id2)==6;
}
bool StandardModelHadronSpectrum::canBeBaryon(tcPDPtr par1, tcPDPtr par2 , tcPDPtr par3) const {
assert(par1 && par2);
long id1 = par1->id(), id2 = par2->id();
if (!par3) {
if( id1*id2 < 0) return false;
if(DiquarkMatcher::Check(id1))
return abs(int(par2->iColour())) == 3 && !DiquarkMatcher::Check(id2);
if(DiquarkMatcher::Check(id2))
return abs(int(par1->iColour())) == 3;
return false;
}
else {
// In this case, to be a baryon, all three components must be (anti-)quarks
// and with the same sign.
return (par1->iColour() == 3 && par2->iColour() == 3 && par3->iColour() == 3) ||
(par1->iColour() == -3 && par2->iColour() == -3 && par3->iColour() == -3);
}
}
diff --git a/Hadronization/StandardModelHadronSpectrum.h b/Hadronization/StandardModelHadronSpectrum.h
--- a/Hadronization/StandardModelHadronSpectrum.h
+++ b/Hadronization/StandardModelHadronSpectrum.h
@@ -1,581 +1,581 @@
// -*- C++ -*-
#ifndef Herwig_StandardModelHadronSpectrum_H
#define Herwig_StandardModelHadronSpectrum_H
//
// This is the declaration of the StandardModelHadronSpectrum class.
//
#include "Herwig/Hadronization/HadronSpectrum.h"
#include <ThePEG/PDT/ParticleData.h>
#include <ThePEG/PDT/StandardMatchers.h>
#include <ThePEG/Repository/EventGenerator.h>
#include <ThePEG/PDT/EnumParticles.h>
#include "ThePEG/Repository/CurrentGenerator.h"
#include <ThePEG/Persistency/PersistentOStream.h>
#include <ThePEG/Persistency/PersistentIStream.h>
namespace Herwig {
using namespace ThePEG;
/**
* Here is the documentation of the StandardModelHadronSpectrum class.
*
* @see \ref StandardModelHadronSpectrumInterfaces "The interfaces"
* defined for StandardModelHadronSpectrum.
*/
class StandardModelHadronSpectrum: public HadronSpectrum {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
StandardModelHadronSpectrum(unsigned int opt);
/**
* The destructor.
*/
virtual ~StandardModelHadronSpectrum();
//@}
public:
/** @name Partonic content */
//@{
/**
* Return the id of the gluon
*/
virtual long gluonId() const { return ParticleID::g; }
/**
* Return the ids of all hadronizing quarks
*/
virtual const vector<long>& hadronizingQuarks() const {
static vector<long> hadronizing =
{ ParticleID::d, ParticleID::u, ParticleID::s, ParticleID::c, ParticleID::b };
return hadronizing;
}
/**
* The light hadronizing quarks
*/
virtual const vector<long>& lightHadronizingQuarks() const {
static vector<long> light =
{ ParticleID::d, ParticleID::u, ParticleID::s };
return light;
}
/**
* The light hadronizing diquarks
*/
virtual const vector<long>& lightHadronizingDiquarks() const {
/**
* Diquarks q==q_0 are not allowed as they need to have antisymmetric
* spin wave-function, which forces the spin to 1
* Diquarks q!=q'_1 are not allowed as they need to have antisymmetric
* spin wave-function, which forces the spin to 1
* */
static vector<long> light = {
// TODO: strange diquarks are turned off for the moment
// since in combination with the current ClusterFission
// they fail (overshoot) to reproduce the Xi and Lambda
// pT spectra.
// One may enable these after the ClusterFission
// kinematics are settled
// ParticleID::sd_0,
// ParticleID::sd_1, // Not allowed by exceptions
// ParticleID::su_0,
// ParticleID::su_1, // Not allowed by exceptions
// ParticleID::ss_1,
ParticleID::uu_1,
ParticleID::dd_1,
ParticleID::ud_0
// TODO why ud_1 not allowed?
// exceptions: Could not find 2103 1 in _table
// but no problem for 2103 2 ???
// ParticleID::ud_1
};
return light;
}
/**
* The heavy hadronizing quarks
*/
virtual const vector<long>& heavyHadronizingQuarks() const {
static vector<long> heavy =
{ ParticleID::c, ParticleID::b };
return heavy;
}
/**
* Return true if any of the possible three input particles contains
* the indicated heavy quark. false otherwise. In the case that
* only the first particle is specified, it can be: an (anti-)quark,
* an (anti-)diquark an (anti-)meson, an (anti-)baryon; in the other
* cases, each pointer is assumed to be either (anti-)quark or
* (anti-)diquark.
*/
virtual bool hasHeavy(long id, tcPDPtr par1, tcPDPtr par2 = PDPtr(), tcPDPtr par3 = PDPtr()) const {
if ( abs(id) == ParticleID::c )
return hasCharm(par1,par2,par3);
if ( abs(id) == ParticleID::b )
return hasBottom(par1,par2,par3);
return false;
}
//@}
/**
* Return the threshold for a cluster to split into a pair of hadrons.
* This is normally the mass of the lightest hadron Pair, but can be
* higher for heavy and exotic clusters
*/
virtual Energy hadronPairThreshold(tcPDPtr par1, tcPDPtr par2) const;
/**
* Return the weight for the given flavour
*/
virtual double pwtQuark(const long& id) const {
switch(id) {
case ParticleID::d: return pwtDquark(); break;
case ParticleID::u: return pwtUquark(); break;
case ParticleID::s: return pwtSquark(); break;
case ParticleID::c: return pwtCquark(); break;
case ParticleID::b: return pwtBquark(); break;
}
return 0.;
}
/**
* The down quark weight.
*/
double pwtDquark() const {
return _pwtDquark;
}
/**
* The up quark weight.
*/
double pwtUquark() const {
return _pwtUquark;
}
/**
* The strange quark weight.
*/
double pwtSquark() const {
return _pwtSquark;
}
/**
* The charm quark weight.
*/
double pwtCquark() const {
return _pwtCquark;
}
/**
* The bottom quark weight.
*/
double pwtBquark() const {
return _pwtBquark;
}
/**
* The diquark weight.
*/
double pwtDIquark() const {
return _pwtDIquark;
}
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
/**
* Return the particle data of the diquark (anti-diquark) made by the two
* quarks (antiquarks) par1, par2.
* @param par1 (anti-)quark data pointer
* @param par2 (anti-)quark data pointer
*/
PDPtr makeDiquark(tcPDPtr par1, tcPDPtr par2) const;
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
*
* The array _repwt is initialized using the interfaces to set different
* weights for different meson multiplets and the constructHadronTable()
* method called to complete the construction of the hadron tables.
*
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
//@}
/**
* Return the id of the diquark (anti-diquark) made by the two
* quarks (antiquarks) of id specified in input (id1, id2).
* Caller must ensure that id1 and id2 are quarks.
*/
long makeDiquarkID(long id1, long id2, long pspin) const;
/**
* Return true if any of the possible three input particles has
* b-flavour;
* false otherwise. In the case that only the first particle is specified,
* it can be: an (anti-)quark, an (anti-)diquark
* an (anti-)meson, an (anti-)baryon; in the other cases, each pointer
* is assumed to be either (anti-)quark or (anti-)diquark.
*/
bool hasBottom(tcPDPtr par1, tcPDPtr par2 = PDPtr(), tcPDPtr par3 = PDPtr()) const;
/**
* Return true if any of the possible three input particles has
* c-flavour;
* false otherwise.In the case that only the first pointer is specified,
* it can be: a (anti-)quark, a (anti-)diquark
* a (anti-)meson, a (anti-)baryon; in the other cases, each pointer
* is assumed to be either (anti-)quark or (anti-)diquark.
*/
bool hasCharm(tcPDPtr par1, tcPDPtr par2 = PDPtr(), tcPDPtr par3 = PDPtr()) const;
/**
* Return true, if any of the possible input particle pointer is an exotic quark, e.g. Susy quark;
* false otherwise.
*/
bool isExotic(tcPDPtr par1, tcPDPtr par2 = PDPtr(), tcPDPtr par3 = PDPtr()) const;
/**
* Return true if the two or three particles in input can be the components
* of a baryon; false otherwise.
*/
virtual bool canBeBaryon(tcPDPtr par1, tcPDPtr par2 , tcPDPtr par3 = PDPtr()) const;
protected:
/**
* Construct the table of hadron data
* This is the main method to initialize the hadron data (mainly the
* weights associated to each hadron, taking into account its spin,
* eventual isoscalar-octect mixing, singlet-decuplet factor). This is
* the method that one should update when new or updated hadron data is
* available.
*
* This class implements the construction of the basic table but can be
* overridden if needed in inheriting classes.
*
* The rationale for factors used for diquarks involving different quarks can
* be can be explained by taking a prototype example that in the exact SU(2) limit,
* in which:
* \f[m_u=m_d\f]
* \f[M_p=M_n=M_\Delta\f]
* and we will have equal numbers of u and d quarks produced.
* Suppose that we weight 1 the diquarks made of the same
* quark and 1/2 those made of different quarks, the fractions
* of u and d baryons (p, n, Delta) we get are the following:
* - \f$\Delta^{++}\f$: 1 possibility only u uu with weight 1
* - \f$\Delta^- \f$: 1 possibility only d dd with weight 1
* - \f$p,\Delta^+ \f$: 2 possibilities u ud with weight 1/2
* d uu with weight 1
* - \f$n,\Delta^0 \f$: 2 possibilities d ud with weight 1/2
* u dd with weight 1
* In the latter two cases, we have to take into account the
* fact that p and n have spin 1/2 whereas Delta+ and Delta0
* have spin 3/2 therefore from phase space we get a double weight
* for Delta+ and Delta0 relative to p and n respectively.
* Therefore the relative amount of these baryons that is
* produced is the following:
* # p = # n = ( 1/2 + 1 ) * 1/3 = 1/2
* # Delta++ = # Delta- = 1 = ( 1/2 + 1) * 2/3 # Delta+ = # Delta0
* which is correct, and therefore the weight 1/2 for the
* diquarks of different types of quarks is justified (at least
* in this limit of exact SU(2) ).
*/
virtual void constructHadronTable();
/**
* Insert a meson in the table
*/
virtual void insertMeson(HadronInfo a, int flav1, int flav2);
/**
* Methods for the mixing of \f$I=0\f$ mesons
*/
//@{
/**
* Return the probability of mixing for Octet-Singlet isoscalar mixing,
* the probability of the
* \f$\frac1{\sqrt{2}}(|u\bar{u}\rangle + |d\bar{d}\rangle)\f$ component
* is returned.
* @param angleMix The mixing angle in degrees (not radians)
* @param order is 0 for no mixing, 1 for the first resonance of a pair,
* 2 for the second one.
* The mixing is defined so that for example with \f$\eta-\eta'\f$ mixing where
* the mixing angle is \f$\theta=-23^0$ with $\eta\f$ as the first particle
* and \f$\eta'\f$ the second one.
* The convention used is
* \f[\eta = \cos\theta|\eta_{\rm octet }\rangle
* -\sin\theta|\eta_{\rm singlet}\rangle\f]
* \f[\eta' = \sin\theta|\eta_{\rm octet }\rangle
* -\cos\theta|\eta_{\rm singlet}\rangle\f]
* with
* \f[|\eta_{\rm singlet}\rangle = \frac1{\sqrt{3}}
* \left[|u\bar{u}\rangle + |d\bar{d}\rangle + |s\bar{s}\rangle\right]\f]
* \f[|\eta_{\rm octet }\rangle = \frac1{\sqrt{6}}
* \left[|u\bar{u}\rangle + |d\bar{d}\rangle - 2|s\bar{s}\rangle\right]\f]
*/
double probabilityMixing(const double angleMix,
const int order) const {
static double convert=Constants::pi/180.0;
if (order == 1)
return sqr( cos( angleMix*convert + atan( sqrt(2.0) ) ) );
else if (order == 2)
return sqr( sin( angleMix*convert + atan( sqrt(2.0) ) ) );
else
return 1.;
}
/**
* Returns the weight of given mixing state.
* @param id The PDG code of the meson
*/
virtual double mixingStateWeight(long id) const;
//@}
virtual double specialQuarkWeight(double quarkWeight, long id,
const Energy cluMass, tcPDPtr par1, tcPDPtr par2) const {
// special for strange
- if(id == 3)
+ if(abs(id) == 3)
return strangeWeight(cluMass,par1,par2);
else
return quarkWeight;
}
/**
* Strange quark weight
*/
virtual double strangeWeight(const Energy cluMass, tcPDPtr par1, tcPDPtr par2) const;
/**
* The weights for the different quarks and diquarks
*/
//@{
/**
* The probability of producting a down quark.
*/
double _pwtDquark;
/**
* The probability of producting an up quark.
*/
double _pwtUquark;
/**
* The probability of producting a strange quark.
*/
double _pwtSquark;
/**
* The probability of producting a charm quark.
*/
double _pwtCquark;
/**
* The probability of producting a bottom quark.
*/
double _pwtBquark;
//@}
/**
* The probability of producting a diquark.
*/
double _pwtDIquark;
/**
* Singlet and Decuplet weights
*/
//@{
/**
* The singlet weight
*/
double _sngWt;
/**
* The decuplet weight
*/
double _decWt;
//@}
/**
* The mixing angles for the \f$I=0\f$ mesons containing light quarks
*/
//@{
/**
* The \f$\eta-\eta'\f$ mixing angle
*/
double _etamix;
/**
* The \f$\phi-\omega\f$ mixing angle
*/
double _phimix;
/**
* The \f$h_1'-h_1\f$ mixing angle
*/
double _h1mix;
/**
* The \f$f_0(1710)-f_0(1370)\f$ mixing angle
*/
double _f0mix;
/**
* The \f$f_1(1420)-f_1(1285)\f$ mixing angle
*/
double _f1mix;
/**
* The \f$f'_2-f_2\f$ mixing angle
*/
double _f2mix;
/**
* The \f$\eta_2(1870)-\eta_2(1645)\f$ mixing angle
*/
double _eta2mix;
/**
* The \f$\phi(???)-\omega(1650)\f$ mixing angle
*/
double _omhmix;
/**
* The \f$\phi_3-\omega_3\f$ mixing angle
*/
double _ph3mix;
/**
* The \f$\eta(1475)-\eta(1295)\f$ mixing angle
*/
double _eta2Smix;
/**
* The \f$\phi(1680)-\omega(1420)\f$ mixing angle
*/
double _phi2Smix;
//@}
/**
* The weights for the various meson multiplets to be used to supress the
* production of particular states
*/
//@{
/**
* The weights for the \f$\phantom{1}^1S_0\f$ multiplets
*/
vector<double> _weight1S0;
/**
* The weights for the \f$\phantom{1}^3S_1\f$ multiplets
*/
vector<double> _weight3S1;
/**
* The weights for the \f$\phantom{1}^1P_1\f$ multiplets
*/
vector<double> _weight1P1;
/**
* The weights for the \f$\phantom{1}^3P_0\f$ multiplets
*/
vector<double> _weight3P0;
/**
* The weights for the \f$\phantom{1}^3P_1\f$ multiplets
*/
vector<double> _weight3P1;
/**
* The weights for the \f$\phantom{1}^3P_2\f$ multiplets
*/
vector<double> _weight3P2;
/**
* The weights for the \f$\phantom{1}^1D_2\f$ multiplets
*/
vector<double> _weight1D2;
/**
* The weights for the \f$\phantom{1}^3D_1\f$ multiplets
*/
vector<double> _weight3D1;
/**
* The weights for the \f$\phantom{1}^3D_2\f$ multiplets
*/
vector<double> _weight3D2;
/**
* The weights for the \f$\phantom{1}^3D_3\f$ multiplets
*/
vector<double> _weight3D3;
//@}
/**
* Option for the construction of the tables
*/
unsigned int _topt;
/**
* Which particles to produce for debugging purposes
*/
unsigned int _trial;
/**
* @name A parameter used for determining when clusters are too light.
*
* This parameter is used for setting the lower threshold, \f$ t \f$ as
* \f[ t' = t(1 + r B^1_{\rm lim}) \f]
* where \f$ r \f$ is a random number [0,1].
*/
//@{
double _limBottom;
double _limCharm;
double _limExotic;
//@}
};
}
#endif /* Herwig_StandardModelHadronSpectrum_H */