Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F7877614
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
65 KB
Subscribers
None
View Options
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 */
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Tue, Nov 19, 3:54 PM (1 d, 21 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805047
Default Alt Text
(65 KB)
Attached To
rHERWIGHG herwighg
Event Timeline
Log In to Comment