diff --git a/Hadronization/CheckId.cc b/Hadronization/CheckId.cc
--- a/Hadronization/CheckId.cc
+++ b/Hadronization/CheckId.cc
@@ -1,164 +1,164 @@
// -*- C++ -*-
//
// CheckId.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the CheckId class.
//
#include "CheckId.h"
#include <cassert>
using namespace Herwig;
namespace {
/**
* 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 ) );
}
}
long CheckId::makeDiquarkID(long id1, long id2, long pspin) {
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";
+ //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;
}
PDPtr CheckId::makeDiquark(tcPDPtr par1, tcPDPtr par2, long pspin) {
long id1 = par1->id();
long id2 = par2->id();
long idnew = makeDiquarkID(id1,id2, pspin);
assert(!CurrentGenerator::isVoid());
return CurrentGenerator::current().getParticleData(idnew);
}
bool CheckId::canBeMeson(tcPDPtr par1,tcPDPtr par2) {
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 ( abs(int(par1->iColour()))== 3 &&
abs(int(par2->iColour())) == 3 &&
id1*id2 < 0);
}
bool CheckId::canBeBaryon(tcPDPtr par1, tcPDPtr par2 , tcPDPtr par3) {
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);
}
}
bool CheckId::hasBottom(tcPDPtr par1, tcPDPtr par2, tcPDPtr par3) {
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 CheckId::hasCharm(tcPDPtr par1, tcPDPtr par2, tcPDPtr par3) {
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 CheckId::isExotic(tcPDPtr par1, tcPDPtr par2, tcPDPtr par3) {
/// \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;
}
diff --git a/Hadronization/ClusterFinder.cc b/Hadronization/ClusterFinder.cc
--- a/Hadronization/ClusterFinder.cc
+++ b/Hadronization/ClusterFinder.cc
@@ -1,486 +1,488 @@
// -*- C++ -*-
//
// ClusterFinder.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the ClusterFinder class.
//
#include "ClusterFinder.h"
#include <ThePEG/Interface/ClassDocumentation.h>
#include <ThePEG/Interface/Switch.h>
#include <ThePEG/Persistency/PersistentOStream.h>
#include <ThePEG/Persistency/PersistentIStream.h>
#include <ThePEG/PDT/StandardMatchers.h>
#include <ThePEG/PDT/EnumParticles.h>
#include <ThePEG/Repository/EventGenerator.h>
#include <ThePEG/EventRecord/Collision.h>
#include "CheckId.h"
#include "Herwig/Utilities/EnumParticles.h"
#include "Herwig/Utilities/Kinematics.h"
#include "Cluster.h"
#include <ThePEG/Utilities/DescribeClass.h>
using namespace Herwig;
DescribeClass<ClusterFinder,Interfaced>
describeClusterFinder("Herwig::ClusterFinder","Herwig.so");
IBPtr ClusterFinder::clone() const {
return new_ptr(*this);
}
IBPtr ClusterFinder::fullclone() const {
return new_ptr(*this);
}
void ClusterFinder::persistentOutput(PersistentOStream & os) const {
os << heavyDiquarks_ << diQuarkSelection_ << diQuarkOnShell_;
}
void ClusterFinder::persistentInput(PersistentIStream & is, int) {
is >> heavyDiquarks_ >> diQuarkSelection_ >> diQuarkOnShell_;
}
void ClusterFinder::Init() {
static ClassDocumentation<ClusterFinder> documentation
("This class is responsible of finding clusters.");
static Switch<ClusterFinder,unsigned int> interfaceHeavyDiquarks
("HeavyDiquarks",
"How to treat heavy quarks in baryon number violating clusters",
&ClusterFinder::heavyDiquarks_, 2, false, false);
static SwitchOption interfaceHeavyDiquarksDefault
(interfaceHeavyDiquarks,
"Allow",
"No special treatment, allow both heavy and doubly heavy diquarks",
0);
static SwitchOption interfaceHeavyDiquarksNoDoublyHeavy
(interfaceHeavyDiquarks,
"NoDoublyHeavy",
"Avoid having diquarks with twoo heavy quarks",
1);
static SwitchOption interfaceHeavyDiquarksNoHeavy
(interfaceHeavyDiquarks,
"NoHeavy",
"Try and avoid diquarks contain c and b altogether",
2);
static Switch<ClusterFinder,unsigned int> interfaceDiQuarkSelection
("DiQuarkSelection",
"Option controlling the selection of quarks to merge into a diquark in baryon-number violating clusters",
&ClusterFinder::diQuarkSelection_, 1, false, false);
static SwitchOption interfaceDiQuarkSelectionRandom
(interfaceDiQuarkSelection,
"Random",
"Randomly pick a pair to combine",
0);
static SwitchOption interfaceDiQuarkSelectionLowestMass
(interfaceDiQuarkSelection,
"LowestMass",
"Combine the lowest mass pair",
1);
static Switch<ClusterFinder,bool> interfaceDiQuarkOnShell
("DiQuarkOnShell",
"Force the diquark produced in baryon-number violating clusters to be on-shell",
&ClusterFinder::diQuarkOnShell_, false, false, false);
static SwitchOption interfaceDiQuarkOnShellYes
(interfaceDiQuarkOnShell,
"Yes",
"Force to be on-shell",
true);
static SwitchOption interfaceDiQuarkOnShellNo
(interfaceDiQuarkOnShell,
"No",
"Leave off-shell",
false);
}
ClusterVector ClusterFinder::formClusters(const PVector & partons) {
set<tPPtr> examinedSet; // colour particles already included in a cluster
map<tColinePtr, pair<tPPtr,tPPtr> > quarkQuark; // quark quark
map<tColinePtr, pair<tPPtr,tPPtr> > aQuarkQuark; // anti quark anti quark
ParticleSet inputParticles(partons.begin(),partons.end());
ClusterVector clusters;
// Loop over all current particles.
for(PVector::const_iterator pit=partons.begin();pit!=partons.end();++pit){
// Skip to the next particle if it is not coloured or already examined.
assert(*pit);
assert((*pit)->dataPtr());
if(!(**pit).data().coloured()
|| examinedSet.find(*pit) != examinedSet.end()) {
continue;
}
// We assume that a cluster is made of, at most, 3 constituents although
// in most cases the number will be 2; however, for baryon violating decays
// (for example in Susy model without R parity conservation) we can have 3
// constituents. In the latter case, a quark (antiquark) do not have an
// anticolour (colour) partner as usual, but its colour line either stems
// from a colour source, or ends in a colour sink. In the case of double
// baryon violating decays, but with overall baryon conservation
// ( for instance:
// tilde_u_R -> dbar_1 + dbar_2
// tilde_u_R_star -> d1 + d2
// where tilde_u_R and tilde_u_R_star are colour connected )
// a special treatment is needed, because first we have to process all
// partons in the current step, and then for each left pair of quarks which
// stem from a colour source we have to find the corresponding pair of
// anti-quarks which ends in a colour sink and is connected with the
// above colour source. These special pairs are kept into the maps:
// spec/CluHadConfig.hialQuarkQuarkMap and specialAntiQuarkAntiQuarkMap.
tParticleVector connected(3);
int iElement = 0;
connected[iElement++] = *pit;
bool specialCase = false;
if((*pit)->hasColour()) {
tPPtr partner =
(*pit)->colourLine()->getColouredParticle(partons.begin(),
partons.end(),
true);
if(partner) {
connected[iElement++]= partner;
}
// colour source : baryon-violating process
else {
if((*pit)->colourLine()->sourceNeighbours() != tColinePair()) {
tColinePair sourcePair = (*pit)->colourLine()->sourceNeighbours();
tColinePtr intCL = tColinePtr();
for(int i = 0; i < 2; ++i) {
tColinePtr pLine = i==0 ? sourcePair.first : sourcePair.second;
int saveNumElements = iElement;
for(tPVector::const_iterator cit = pLine->coloured().begin();
cit != pLine->coloured().end(); ++cit ) {
ParticleSet::const_iterator cjt = inputParticles.find(*cit);
if(cjt!=inputParticles.end()) connected[iElement++]= (*cit);
}
if(iElement == saveNumElements) intCL = pLine;
}
if(intCL && iElement == 2) {
specialCase = true;
pair<tPPtr,tPPtr> qp=pair<tPPtr,tPPtr>(connected[0],connected[1]);
quarkQuark.insert(pair<tColinePtr,pair<tPPtr,tPPtr> >(intCL,qp));
}
else if(iElement != 3) {
throw Exception() << "Colour connections fail in the hadronization for "
<< **pit << "in ClusterFinder::formClusters"
<< " for a coloured particle."
<< " Failed to find particles from a source"
<< Exception::runerror;
}
}
else {
throw Exception() << "Colour connections fail in the hadronization for "
<< **pit << "in ClusterFinder::formClusters for"
<< " a coloured particle"
<< Exception::runerror;
}
}
}
if((*pit)->hasAntiColour()) {
tPPtr partner =
(*pit)->antiColourLine()->getColouredParticle(partons.begin(),
partons.end(),
false);
if(partner) {
connected[iElement++]=partner;
}
// colour sink : baryon-violating process
else {
if((*pit)->antiColourLine()->sinkNeighbours() != tColinePair()) {
tColinePair sinkPair = (*pit)->antiColourLine()->sinkNeighbours();
tColinePtr intCL = tColinePtr();
for(int i = 0; i < 2; ++i) {
tColinePtr pLine = i==0 ? sinkPair.first : sinkPair.second;
int saveNumElements = iElement;
for(tPVector::const_iterator cit = pLine->antiColoured().begin();
cit != pLine->antiColoured().end(); ++cit ) {
ParticleSet::const_iterator cjt = inputParticles.find(*cit);
if(cjt!=inputParticles.end()) connected[iElement++]= (*cit);
}
if(iElement == saveNumElements) intCL = pLine;
}
if(intCL && iElement == 2) {
specialCase = true;
pair<tPPtr,tPPtr> aqp=pair<tPPtr,tPPtr>(connected[0],connected[1]);
aQuarkQuark.insert(pair<tColinePtr,pair<tPPtr,tPPtr> >(intCL,aqp));
}
else if( iElement !=3) {
throw Exception() << "Colour connections fail in the hadronization for "
<< **pit << "in ClusterFinder::formClusters for"
<< " an anti-coloured particle."
<< " Failed to find particles from a sink"
<< Exception::runerror;
}
}
else {
throw Exception() << "Colour connections fail in the hadronization for "
<< **pit << "in ClusterFinder::formClusters for"
<< " an anti-coloured particle"
<< Exception::runerror;
}
}
}
if(!specialCase) {
// Tag the components of the found cluster as already examined.
for (int i=0; i<iElement; ++i) examinedSet.insert(connected[i]);
// Create the cluster object with the colour connected particles
ClusterPtr cluPtr = new_ptr(Cluster(connected[0],connected[1],
connected[2]));
// add to the step
connected[0]->addChild(cluPtr);
connected[1]->addChild(cluPtr);
if(connected[2]) connected[2]->addChild(cluPtr);
clusters.push_back(cluPtr);
// Check if any of the components is a beam remnant, and if this
// is the case then inform the cluster.
// this will only work for baryon collisions
for (int i=0; i<iElement; ++i) {
bool fromRemnant = false;
tPPtr parent=connected[i];
while(parent) {
if(parent->id()==ParticleID::Remnant) {
fromRemnant = true;
break;
}
parent = parent->parents().empty() ? tPPtr() : parent->parents()[0];
}
if(fromRemnant&&DiquarkMatcher::Check(connected[i]->id()))
cluPtr->isBeamCluster(connected[i]);
}
}
}
// Treat now the special cases, if any. The idea is to find for each pair
// of quarks coming from a common colour source the corresponding pair of
// antiquarks coming from a common colour sink, connected to the above
// colour source via the same colour line. Then, randomly couple one of
// the two quarks with one of the two antiquarks, and do the same with the
// quark and antiquark left.
for(map<tColinePtr, pair<tPPtr,tPPtr> >::const_iterator
cit = quarkQuark.begin(); cit != quarkQuark.end(); ++cit ) {
tColinePtr coline = cit->first;
pair<tPPtr,tPPtr> quarkPair = cit->second;
if(aQuarkQuark.find( coline ) != aQuarkQuark.end()) {
pair<tPPtr,tPPtr> antiQuarkPair = aQuarkQuark.find(coline)->second;
ClusterPtr cluPtr1, cluPtr2;
if ( UseRandom::rndbool() ) {
cluPtr1 = new_ptr(Cluster(quarkPair.first , antiQuarkPair.first));
cluPtr2 = new_ptr(Cluster(quarkPair.second , antiQuarkPair.second));
quarkPair.first->addChild(cluPtr1);
antiQuarkPair.first->addChild(cluPtr1);
quarkPair.second->addChild(cluPtr2);
antiQuarkPair.second->addChild(cluPtr2);
} else {
cluPtr1 = new_ptr(Cluster(quarkPair.first , antiQuarkPair.second));
cluPtr2 = new_ptr(Cluster(quarkPair.second , antiQuarkPair.first));
quarkPair.second->addChild(cluPtr2);
antiQuarkPair.first->addChild(cluPtr2);
quarkPair.first->addChild(cluPtr1);
antiQuarkPair.second->addChild(cluPtr1);
}
clusters.push_back(cluPtr1);
clusters.push_back(cluPtr2);
}
else {
throw Exception() << "ClusterFinder::formClusters : "
<< "***Skip event: unable to match pairs in "
<< "Baryon-violating processes***"
<< Exception::eventerror;
}
}
return clusters;
}
namespace {
bool PartOrdering(tPPtr p1,tPPtr p2) {
return abs(p1->id())<abs(p2->id());
}
}
void ClusterFinder::reduceToTwoComponents(ClusterVector & clusters) {
// In order to preserve all of the information, we do not modify the
// directly the 3-component clusters, but instead we define new clusters,
// which are related to the original ones by a child-parent relationship,
// by considering two randomly chosen components as a diquark (or anti-diquark).
// These new clusters are first added to the vector vecNewRedefinedCluPtr,
// and at the end, when all input clusters have been examined, the elements of
// this vector will be copied in collecCluPtr (the reason is that it is not
// allowed to modify a STL container while iterating over it).
vector<tClusterPtr> redefinedClusters;
for(ClusterVector::iterator cluIter = clusters.begin() ;
cluIter != clusters.end() ; ++cluIter) {
tParticleVector vec;
if ( (*cluIter)->numComponents() != 3 ||
! (*cluIter)->isAvailable() ) continue;
tPPtr other;
for(unsigned int i = 0; i<(*cluIter)->numComponents(); i++) {
tPPtr part = (*cluIter)->particle(i);
if(!DiquarkMatcher::Check(*(part->dataPtr())))
vec.push_back(part);
else
other = part;
}
if(vec.size()<2) {
throw Exception() << "Could not make a diquark for a baryonic cluster decay from "
<< (*cluIter)->particle(0)->PDGName() << " "
<< (*cluIter)->particle(1)->PDGName() << " "
<< (*cluIter)->particle(2)->PDGName() << " "
<< " in ClusterFinder::reduceToTwoComponents()."
<< Exception::eventerror;
}
// order the vector so heaviest at the end
std::sort(vec.begin(),vec.end(),PartOrdering);
// Special treatment of heavy quarks
// avoid doubly heavy diquarks
if(heavyDiquarks_>=1 && vec.size()>2 &&
abs(vec[1]->id())>3 && abs(vec[0]->id())<=3) {
if(UseRandom::rndbool()) swap(vec[1],vec[2]);
other = vec[2];
vec.pop_back();
}
// avoid singly heavy diquarks
if(heavyDiquarks_==2 && vec.size()>2 &&
abs(vec[2]->id())>3 && abs(vec[1]->id())<=3) {
other = vec[2];
vec.pop_back();
}
// if there's a choice pick the pair to make a diquark from
if(vec.size()>2) {
unsigned int ichoice(0);
// random choice
if(diQuarkSelection_==0) {
ichoice = UseRandom::rnd3(1.0, 1.0, 1.0);
}
// pick the lightest quark pair
else if(diQuarkSelection_==1) {
Energy m12 = (vec[0]->momentum()+vec[1]->momentum()).m();
Energy m13 = (vec[0]->momentum()+vec[2]->momentum()).m();
Energy m23 = (vec[1]->momentum()+vec[2]->momentum()).m();
if (m13<=m12&&m13<=m23) ichoice = 2;
else if(m23<=m12&&m23<=m13) ichoice = 1;
}
else
assert(false);
// make the swaps so select pair first
switch (ichoice) {
case 0:
break;
case 1:
swap(vec[2],vec[0]);
break;
case 2:
swap(vec[2],vec[1]);
break;
}
}
// set up
tcPDPtr temp1 = vec[0]->dataPtr();
tcPDPtr temp2 = vec[1]->dataPtr();
if(!other) other = vec[2];
- long rndSpin = UseRandom::rnd() > 0.5 ? 1 : 3;
+ long rndSpin;
+ if(temp1 == temp2) rndSpin = 3;
+ else rndSpin = UseRandom::rnd() > 0.5 ? 1 : 3;
tcPDPtr dataDiquark = CheckId::makeDiquark(temp1,temp2,rndSpin);
if(!dataDiquark)
throw Exception() << "Could not make a diquark from"
<< temp1->PDGName() << " and "
<< temp2->PDGName()
<< " in ClusterFinder::reduceToTwoComponents()"
<< Exception::eventerror;
// Create the new cluster (with two components) and assign to it the same
// momentum and position of the original (with three components) one.
// Furthermore, assign to the diquark component a momentum given by the
// sum of the two original components from which has been formed; for the
// position, we are assuming, very simply, that the diquark position is
// the average positions of the two original components.
// Notice that the mass (5-th component of the 5-momentum) of the diquark
// is set by hand to the constituent mass of the diquark (which is equal
// to the sum of the constituent masses of the two quarks which form the
// diquark) because the sum of 5-component vectors do add only the "normal"
// 4-components, not the 5-th one. After that, the 5-momentum of the diquark
// is in an inconsistent state, because the mass (5-th component) is not
// equal to the invariant mass obtained from the 4-momemtum. This is not
// unique to this kind of component (all perturbative components are in
// a similar situation), but it is not harmful.
// construct the diquark
PPtr diquark = dataDiquark->produceParticle();
vec[0]->addChild(diquark);
vec[1]->addChild(diquark);
diquark->set5Momentum(Lorentz5Momentum(vec[0]->momentum() + vec[1]->momentum(),
dataDiquark->constituentMass()));
// use the same method as for cluster to determine the diquark position
diquark->setVertex(Cluster::calculateX(vec[0],vec[1]));
// put on-shell if required
if(diQuarkOnShell_) {
Lorentz5Momentum psum = diquark->momentum()+other->momentum();
psum.rescaleMass();
Boost boost = psum.boostVector();
Lorentz5Momentum pother = other->momentum();
Lorentz5Momentum pdiquark = diquark->momentum();
pother.boost(-boost);
pdiquark.boost(-boost);
Energy pcm = Kinematics::pstarTwoBodyDecay(psum.mass(),
other->dataPtr()->constituentMass(),
diquark->dataPtr()->constituentMass());
if(pcm>ZERO) {
double fact = pcm/pother.vect().mag();
pother *= fact;
pdiquark *= fact;
pother .setMass(other->dataPtr()->constituentMass());
pdiquark.setMass(dataDiquark ->constituentMass());
pother .rescaleEnergy();
pdiquark.rescaleEnergy();
pother .boost(boost);
pdiquark.boost(boost);
other->set5Momentum(pother);
diquark->set5Momentum(pdiquark);
}
}
// make the new cluster
ClusterPtr nclus = new_ptr(Cluster(other,diquark));
//vec[0]->addChild(nclus);
//diquark->addChild(nclus);
// Set the parent/children relationship between the original cluster
// (the one with three components) with the new one (the one with two components)
// and add the latter to the vector of new redefined clusters.
(*cluIter)->addChild(nclus);
redefinedClusters.push_back(nclus);
}
// Add to collecCluPtr all of the redefined new clusters (indeed the
// pointers to them are added) contained in vecNewRedefinedCluPtr.
/// \todo why do we keep the original of the redefined clusters?
for (tClusterVector::const_iterator it = redefinedClusters.begin();
it != redefinedClusters.end(); ++it) {
clusters.push_back(*it);
}
}
diff --git a/Hadronization/ClusterFissioner.cc b/Hadronization/ClusterFissioner.cc
--- a/Hadronization/ClusterFissioner.cc
+++ b/Hadronization/ClusterFissioner.cc
@@ -1,1122 +1,1124 @@
// -*- C++ -*-
//
// ClusterFissioner.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// Thisk is the implementation of the non-inlined, non-templated member
// functions of the ClusterFissioner class.
//
#include "ClusterFissioner.h"
#include <ThePEG/Interface/ClassDocumentation.h>
#include <ThePEG/Interface/Reference.h>
#include <ThePEG/Interface/Parameter.h>
#include <ThePEG/Interface/Switch.h>
#include <ThePEG/Persistency/PersistentOStream.h>
#include <ThePEG/Persistency/PersistentIStream.h>
#include <ThePEG/PDT/EnumParticles.h>
#include "Herwig/Utilities/Kinematics.h"
#include "CheckId.h"
#include "Cluster.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include <ThePEG/Utilities/DescribeClass.h>
using namespace Herwig;
DescribeClass<ClusterFissioner,Interfaced>
describeClusterFissioner("Herwig::ClusterFissioner","Herwig.so");
ClusterFissioner::ClusterFissioner() :
_clMaxLight(3.35*GeV),
_clMaxBottom(3.35*GeV),
_clMaxCharm(3.35*GeV),
_clMaxExotic(3.35*GeV),
_clPowLight(2.0),
_clPowBottom(2.0),
_clPowCharm(2.0),
_clPowExotic(2.0),
_pSplitLight(1.0),
_pSplitBottom(1.0),
_pSplitCharm(1.0),
_pSplitExotic(1.0),
_fissionPwtUquark(1),
_fissionPwtDquark(1),
_fissionPwtSquark(0.5),
_fissionCluster(0),
_btClM(1.0*GeV),
_iopRem(1),
_kappa(1.0e15*GeV/meter),
_enhanceSProb(0),
_m0Fission(2.*GeV),
_massMeasure(0)
{}
IBPtr ClusterFissioner::clone() const {
return new_ptr(*this);
}
IBPtr ClusterFissioner::fullclone() const {
return new_ptr(*this);
}
void ClusterFissioner::persistentOutput(PersistentOStream & os) const {
os << _hadronsSelector << ounit(_clMaxLight,GeV)
<< ounit(_clMaxBottom,GeV) << ounit(_clMaxCharm,GeV)
<< ounit(_clMaxExotic,GeV) << _clPowLight << _clPowBottom
<< _clPowCharm << _clPowExotic << _pSplitLight
<< _pSplitBottom << _pSplitCharm << _pSplitExotic
<< _fissionCluster << _fissionPwtUquark << _fissionPwtDquark << _fissionPwtSquark
<< ounit(_btClM,GeV)
<< _iopRem << ounit(_kappa, GeV/meter)
<< _enhanceSProb << ounit(_m0Fission,GeV) << _massMeasure;
}
void ClusterFissioner::persistentInput(PersistentIStream & is, int) {
is >> _hadronsSelector >> iunit(_clMaxLight,GeV)
>> iunit(_clMaxBottom,GeV) >> iunit(_clMaxCharm,GeV)
>> iunit(_clMaxExotic,GeV) >> _clPowLight >> _clPowBottom
>> _clPowCharm >> _clPowExotic >> _pSplitLight
>> _pSplitBottom >> _pSplitCharm >> _pSplitExotic
>> _fissionCluster >> _fissionPwtUquark >> _fissionPwtDquark >> _fissionPwtSquark
>> iunit(_btClM,GeV) >> _iopRem
>> iunit(_kappa, GeV/meter)
>> _enhanceSProb >> iunit(_m0Fission,GeV) >> _massMeasure;
}
void ClusterFissioner::Init() {
static ClassDocumentation<ClusterFissioner> documentation
("Class responsibles for chopping up the clusters");
static Reference<ClusterFissioner,HadronSelector>
interfaceHadronSelector("HadronSelector",
"A reference to the HadronSelector object",
&Herwig::ClusterFissioner::_hadronsSelector,
false, false, true, false);
// ClMax for light, Bottom, Charm and exotic (e.g. Susy) quarks
static Parameter<ClusterFissioner,Energy>
interfaceClMaxLight ("ClMaxLight","cluster max mass for light quarks (unit [GeV])",
&ClusterFissioner::_clMaxLight, GeV, 3.35*GeV, ZERO, 10.0*GeV,
false,false,false);
static Parameter<ClusterFissioner,Energy>
interfaceClMaxBottom ("ClMaxBottom","cluster max mass for b quarks (unit [GeV])",
&ClusterFissioner::_clMaxBottom, GeV, 3.35*GeV, ZERO, 10.0*GeV,
false,false,false);
static Parameter<ClusterFissioner,Energy>
interfaceClMaxCharm ("ClMaxCharm","cluster max mass for c quarks (unit [GeV])",
&ClusterFissioner::_clMaxCharm, GeV, 3.35*GeV, ZERO, 10.0*GeV,
false,false,false);
static Parameter<ClusterFissioner,Energy>
interfaceClMaxExotic ("ClMaxExotic","cluster max mass for exotic quarks (unit [GeV])",
&ClusterFissioner::_clMaxExotic, GeV, 3.35*GeV, ZERO, 10.0*GeV,
false,false,false);
// ClPow for light, Bottom, Charm and exotic (e.g. Susy) quarks
static Parameter<ClusterFissioner,double>
interfaceClPowLight ("ClPowLight","cluster mass exponent for light quarks",
&ClusterFissioner::_clPowLight, 0, 2.0, 0.0, 10.0,false,false,false);
static Parameter<ClusterFissioner,double>
interfaceClPowBottom ("ClPowBottom","cluster mass exponent for b quarks",
&ClusterFissioner::_clPowBottom, 0, 2.0, 0.0, 10.0,false,false,false);
static Parameter<ClusterFissioner,double>
interfaceClPowCharm ("ClPowCharm","cluster mass exponent for c quarks",
&ClusterFissioner::_clPowCharm, 0, 2.0, 0.0, 10.0,false,false,false);
static Parameter<ClusterFissioner,double>
interfaceClPowExotic ("ClPowExotic","cluster mass exponent for exotic quarks",
&ClusterFissioner::_clPowExotic, 0, 2.0, 0.0, 10.0,false,false,false);
// PSplit for light, Bottom, Charm and exotic (e.g. Susy) quarks
static Parameter<ClusterFissioner,double>
interfacePSplitLight ("PSplitLight","cluster mass splitting param for light quarks",
&ClusterFissioner::_pSplitLight, 0, 1.0, 0.0, 10.0,false,false,false);
static Parameter<ClusterFissioner,double>
interfacePSplitBottom ("PSplitBottom","cluster mass splitting param for b quarks",
&ClusterFissioner::_pSplitBottom, 0, 1.0, 0.0, 10.0,false,false,false);
static Parameter<ClusterFissioner,double>
interfacePSplitCharm ("PSplitCharm","cluster mass splitting param for c quarks",
&ClusterFissioner::_pSplitCharm, 0, 1.0, 0.0, 10.0,false,false,false);
static Parameter<ClusterFissioner,double>
interfacePSplitExotic ("PSplitExotic","cluster mass splitting param for exotic quarks",
&ClusterFissioner::_pSplitExotic, 0, 1.0, 0.0, 10.0,false,false,false);
static Switch<ClusterFissioner,int> interfaceFission
("Fission",
"Option for different Fission options",
&ClusterFissioner::_fissionCluster, 1, false, false);
static SwitchOption interfaceFissionDefault
(interfaceFission,
"default",
"Normal cluster fission which depends on the hadron selector class.",
0);
static SwitchOption interfaceFissionNew
(interfaceFission,
"new",
"Alternative cluster fission which does not depend on the hadron selector class",
1);
static Parameter<ClusterFissioner,double> interfaceFissionPwtUquark
("FissionPwtUquark",
"Weight for fission in U quarks",
&ClusterFissioner::_fissionPwtUquark, 1, 0.0, 1.0,
false, false, Interface::limited);
static Parameter<ClusterFissioner,double> interfaceFissionPwtDquark
("FissionPwtDquark",
"Weight for fission in D quarks",
&ClusterFissioner::_fissionPwtDquark, 1, 0.0, 1.0,
false, false, Interface::limited);
static Parameter<ClusterFissioner,double> interfaceFissionPwtSquark
("FissionPwtSquark",
"Weight for fission in S quarks",
&ClusterFissioner::_fissionPwtSquark, 0.5, 0.0, 1.0,
false, false, Interface::limited);
static Switch<ClusterFissioner,int> interfaceRemnantOption
("RemnantOption",
"Option for the treatment of remnant clusters",
&ClusterFissioner::_iopRem, 1, false, false);
static SwitchOption interfaceRemnantOptionSoft
(interfaceRemnantOption,
"Soft",
"Both clusters produced in the fission of the beam cluster"
" are treated as soft clusters.",
0);
static SwitchOption interfaceRemnantOptionHard
(interfaceRemnantOption,
"Hard",
"Only the cluster containing the remnant is treated as a soft cluster.",
1);
static SwitchOption interfaceRemnantOptionVeryHard
(interfaceRemnantOption,
"VeryHard",
"Even remnant clusters are treated as hard, i.e. all clusters the same",
2);
static Parameter<ClusterFissioner,Energy> interfaceBTCLM
("SoftClusterFactor",
"Parameter for the mass spectrum of remnant clusters",
&ClusterFissioner::_btClM, GeV, 1.*GeV, 0.1*GeV, 10.0*GeV,
false, false, Interface::limited);
static Parameter<ClusterFissioner,Tension> interfaceStringTension
("StringTension",
"String tension used in vertex displacement calculation",
&ClusterFissioner::_kappa, GeV/meter,
1.0e15*GeV/meter, ZERO, ZERO,
false, false, Interface::lowerlim);
static Switch<ClusterFissioner,int> interfaceEnhanceSProb
("EnhanceSProb",
"Option for enhancing strangeness",
&ClusterFissioner::_enhanceSProb, 0, false, false);
static SwitchOption interfaceEnhanceSProbNo
(interfaceEnhanceSProb,
"No",
"No strangeness enhancement.",
0);
static SwitchOption interfaceEnhanceSProbScaled
(interfaceEnhanceSProb,
"Scaled",
"Scaled strangeness enhancement",
1);
static SwitchOption interfaceEnhanceSProbExponential
(interfaceEnhanceSProb,
"Exponential",
"Exponential strangeness enhancement",
2);
static Switch<ClusterFissioner,int> interfaceMassMeasure
("MassMeasure",
"Option to use different mass measures",
&ClusterFissioner::_massMeasure,0,false,false);
static SwitchOption interfaceMassMeasureMass
(interfaceMassMeasure,
"Mass",
"Mass Measure",
0);
static SwitchOption interfaceMassMeasureLambda
(interfaceMassMeasure,
"Lambda",
"Lambda Measure",
1);
static Parameter<ClusterFissioner,Energy> interfaceFissionMassScale
("FissionMassScale",
"Cluster fission mass scale",
&ClusterFissioner::_m0Fission, GeV, 2.0*GeV, 0.1*GeV, 50.*GeV,
false, false, Interface::limited);
}
tPVector ClusterFissioner::fission(ClusterVector & clusters, bool softUEisOn) {
// return if no clusters
if (clusters.empty()) return tPVector();
/*****************
* Loop over the (input) collection of cluster pointers, and store in
* the vector splitClusters all the clusters that need to be split
* (these are beam clusters, if soft underlying event is off, and
* heavy non-beam clusters).
********************/
stack<ClusterPtr> splitClusters;
for(ClusterVector::iterator it = clusters.begin() ;
it != clusters.end() ; ++it) {
/**************
* Skip 3-component clusters that have been redefined (as 2-component
* clusters) or not available clusters. The latter check is indeed
* redundant now, but it is used for possible future extensions in which,
* for some reasons, some of the clusters found by ClusterFinder are tagged
* straight away as not available.
**************/
if((*it)->isRedefined() || !(*it)->isAvailable()) continue;
// if the cluster is a beam cluster add it to the vector of clusters
// to be split or if it is heavy
if((*it)->isBeamCluster() || isHeavy(*it)) splitClusters.push(*it);
}
tPVector finalhadrons;
cut(splitClusters, clusters, finalhadrons, softUEisOn);
return finalhadrons;
}
void ClusterFissioner::cut(stack<ClusterPtr> & clusterStack,
ClusterVector &clusters, tPVector & finalhadrons,
bool softUEisOn) {
/**************************************************
* This method does the splitting of the cluster pointed by cluPtr
* and "recursively" by all of its cluster children, if heavy. All of these
* new children clusters are added (indeed the pointers to them) to the
* collection of cluster pointers collecCluPtr. The method works as follows.
* Initially the vector vecCluPtr contains just the input pointer to the
* cluster to be split. Then it will be filled "recursively" by all
* of the cluster's children that are heavy enough to require, in their turn,
* to be split. In each loop, the last element of the vector vecCluPtr is
* considered (only once because it is then removed from the vector).
* This approach is conceptually recursive, but avoid the overhead of
* a concrete recursive function. Furthermore it requires minimal changes
* in the case that the fission of an heavy cluster could produce more
* than two cluster children as assumed now.
*
* Draw the masses: for normal, non-beam clusters a power-like mass dist
* is used, whereas for beam clusters a fast-decreasing exponential mass
* dist is used instead (to avoid many iterative splitting which could
* produce an unphysical large transverse energy from a supposed soft beam
* remnant process).
****************************************/
// Here we recursively loop over clusters in the stack and cut them
while (!clusterStack.empty()) {
// take the last element of the vector
ClusterPtr iCluster = clusterStack.top(); clusterStack.pop();
// split it
cutType ct = iCluster->numComponents() == 2 ?
cutTwo(iCluster, finalhadrons, softUEisOn) :
cutThree(iCluster, finalhadrons, softUEisOn);
// There are cases when we don't want to split, even if it fails mass test
if(!ct.first.first || !ct.second.first) {
// if an unsplit beam cluster leave if for the underlying event
if(iCluster->isBeamCluster() && softUEisOn)
iCluster->isAvailable(false);
continue;
}
// check if clusters
ClusterPtr one = dynamic_ptr_cast<ClusterPtr>(ct.first.first);
ClusterPtr two = dynamic_ptr_cast<ClusterPtr>(ct.second.first);
// is a beam cluster must be split into two clusters
if(iCluster->isBeamCluster() && (!one||!two) && softUEisOn) {
iCluster->isAvailable(false);
continue;
}
// There should always be a intermediate quark(s) from the splitting
assert(ct.first.second && ct.second.second);
/// \todo sort out motherless quark pairs here. Watch out for 'quark in final state' errors
iCluster->addChild(ct.first.first);
// iCluster->addChild(ct.first.second);
// ct.first.second->addChild(ct.first.first);
iCluster->addChild(ct.second.first);
// iCluster->addChild(ct.second.second);
// ct.second.second->addChild(ct.second.first);
// Sometimes the clusters decay C -> H + C' rather then C -> C' + C''
if(one) {
clusters.push_back(one);
if(one->isBeamCluster() && softUEisOn)
one->isAvailable(false);
if(isHeavy(one) && one->isAvailable())
clusterStack.push(one);
}
if(two) {
clusters.push_back(two);
if(two->isBeamCluster() && softUEisOn)
two->isAvailable(false);
if(isHeavy(two) && two->isAvailable())
clusterStack.push(two);
}
}
}
namespace {
/**
* Check if can't make a hadron from the partons
*/
bool cantMakeHadron(tcPPtr p1, tcPPtr p2) {
return ! CheckId::canBeHadron(p1->dataPtr(), p2->dataPtr());
}
/**
* Check if can't make a diquark from the partons
*/
bool cantMakeDiQuark(tcPPtr p1, tcPPtr p2) {
long id1 = p1->id(), id2 = p2->id();
return ! (QuarkMatcher::Check(id1) && QuarkMatcher::Check(id2) && id1*id2>0);
}
}
ClusterFissioner::cutType
ClusterFissioner::cutTwo(ClusterPtr & cluster, tPVector & finalhadrons,
bool softUEisOn) {
// need to make sure only 2-cpt clusters get here
assert(cluster->numComponents() == 2);
tPPtr ptrQ1 = cluster->particle(0);
tPPtr ptrQ2 = cluster->particle(1);
Energy Mc = cluster->mass();
assert(ptrQ1);
assert(ptrQ2);
// And check if those particles are from a beam remnant
bool rem1 = cluster->isBeamRemnant(0);
bool rem2 = cluster->isBeamRemnant(1);
// workout which distribution to use
bool soft1(false),soft2(false);
switch (_iopRem) {
case 0:
soft1 = rem1 || rem2;
soft2 = rem2 || rem1;
break;
case 1:
soft1 = rem1;
soft2 = rem2;
break;
}
// Initialization for the exponential ("soft") mass distribution.
static const int max_loop = 1000;
int counter = 0;
Energy Mc1 = ZERO, Mc2 = ZERO,m1=ZERO,m2=ZERO,m=ZERO;
tcPDPtr toHadron1, toHadron2;
PPtr newPtr1 = PPtr ();
PPtr newPtr2 = PPtr ();
bool succeeded = false;
do
{
succeeded = false;
++counter;
if (_enhanceSProb == 0){
drawNewFlavour(newPtr1,newPtr2);
}
else {
Energy2 mass2 = clustermass(cluster);
drawNewFlavourEnhanced(newPtr1,newPtr2,mass2);
}
// check for right ordering
assert (ptrQ2);
assert (newPtr2);
assert (ptrQ2->dataPtr());
assert (newPtr2->dataPtr());
if(cantMakeHadron(ptrQ1, newPtr1) || cantMakeHadron(ptrQ2, newPtr2)) {
swap(newPtr1, newPtr2);
// check again
if(cantMakeHadron(ptrQ1, newPtr1) || cantMakeHadron(ptrQ2, newPtr2)) {
throw Exception()
<< "ClusterFissioner cannot split the cluster ("
<< ptrQ1->PDGName() << ' ' << ptrQ2->PDGName()
<< ") into hadrons.\n" << Exception::runerror;
}
}
// Check that new clusters can produce particles and there is enough
// phase space to choose the drawn flavour
m1 = ptrQ1->data().constituentMass();
m2 = ptrQ2->data().constituentMass();
m = newPtr1->data().constituentMass();
// Do not split in the case there is no phase space available
if(Mc < m1+m + m2+m) continue;
// power for splitting
double exp1=_pSplitLight;
double exp2=_pSplitLight;
if (CheckId::isExotic(ptrQ1->dataPtr())) exp1 = _pSplitExotic;
else if(CheckId::hasBottom(ptrQ1->dataPtr()))exp1 = _pSplitBottom;
else if(CheckId::hasCharm(ptrQ1->dataPtr())) exp1 = _pSplitCharm;
if (CheckId::isExotic(ptrQ2->dataPtr())) exp2 = _pSplitExotic;
else if(CheckId::hasBottom(ptrQ2->dataPtr())) exp2 = _pSplitBottom;
else if(CheckId::hasCharm(ptrQ2->dataPtr())) exp2 = _pSplitCharm;
// If, during the drawing of candidate masses, too many attempts fail
// (because the phase space available is tiny)
/// \todo run separate loop here?
Mc1 = drawChildMass(Mc,m1,m2,m,exp1,soft1);
Mc2 = drawChildMass(Mc,m2,m1,m,exp2,soft2);
if(Mc1 < m1+m || Mc2 < m+m2 || Mc1+Mc2 > Mc) continue;
/**************************
* New (not present in Fortran Herwig):
* check whether the fragment masses Mc1 and Mc2 are above the
* threshold for the production of the lightest pair of hadrons with the
* right flavours. If not, then set by hand the mass to the lightest
* single hadron with the right flavours, in order to solve correctly
* the kinematics, and (later in this method) create directly such hadron
* and add it to the children hadrons of the cluster that undergoes the
* fission (i.e. the one pointed by iCluPtr). Notice that in this special
* case, the heavy cluster that undergoes the fission has one single
* cluster child and one single hadron child. We prefer this approach,
* rather than to create a light cluster, with the mass set equal to
* the lightest hadron, and let then the class LightClusterDecayer to do
* the job to decay it to that single hadron, for two reasons:
* First, because the sum of the masses of the two constituents can be,
* in this case, greater than the mass of that hadron, hence it would
* be impossible to solve the kinematics for such two components, and
* therefore we would have a cluster whose components are undefined.
* Second, the algorithm is faster, because it avoids the reshuffling
* procedure that would be necessary if we used LightClusterDecayer
* to decay the light cluster to the lightest hadron.
****************************/
toHadron1 = _hadronsSelector->chooseSingleHadron(ptrQ1->dataPtr(), newPtr1->dataPtr(),Mc1);
if(toHadron1) Mc1 = toHadron1->mass();
toHadron2 = _hadronsSelector->chooseSingleHadron(ptrQ2->dataPtr(), newPtr2->dataPtr(),Mc2);
if(toHadron2) Mc2 = toHadron2->mass();
// if a beam cluster not allowed to decay to hadrons
if(cluster->isBeamCluster() && (toHadron1||toHadron2) && softUEisOn)
continue;
// Check if the decay kinematics is still possible: if not then
// force the one-hadron decay for the other cluster as well.
if(Mc1 + Mc2 > Mc) {
if(!toHadron1) {
toHadron1 = _hadronsSelector->chooseSingleHadron(ptrQ1->dataPtr(), newPtr1->dataPtr(),Mc-Mc2);
if(toHadron1) Mc1 = toHadron1->mass();
}
else if(!toHadron2) {
toHadron2 = _hadronsSelector->chooseSingleHadron(ptrQ2->dataPtr(), newPtr2->dataPtr(),Mc-Mc1);
if(toHadron2) Mc2 = toHadron2->mass();
}
}
succeeded = (Mc >= Mc1+Mc2);
}
while (!succeeded && counter < max_loop);
if(counter >= max_loop) {
static const PPtr null = PPtr();
return cutType(PPair(null,null),PPair(null,null));
}
// Determined the (5-components) momenta (all in the LAB frame)
Lorentz5Momentum pClu = cluster->momentum(); // known
Lorentz5Momentum p0Q1 = ptrQ1->momentum(); // known (mom Q1 before fission)
Lorentz5Momentum pClu1, pClu2, pQ1, pQone, pQtwo, pQ2; //unknown
pClu1.setMass(Mc1);
pClu2.setMass(Mc2);
pQ1.setMass(m1);
pQ2.setMass(m2);
pQone.setMass(m);
pQtwo.setMass(m);
calculateKinematics(pClu,p0Q1,toHadron1,toHadron2,
pClu1,pClu2,pQ1,pQone,pQtwo,pQ2); // out
/******************
* The previous methods have determined the kinematics and positions
* of C -> C1 + C2.
* In the case that one of the two product is light, that means either
* decayOneHadronClu1 or decayOneHadronClu2 is true, then the momenta
* of the components of that light product have not been determined,
* and a (light) cluster will not be created: the heavy father cluster
* decays, in this case, into a single (not-light) cluster and a
* single hadron. In the other, "normal", cases the father cluster
* decays into two clusters, each of which has well defined components.
* Notice that, in the case of components which point to particles, the
* momenta of the components is properly set to the new values, whereas
* we do not change the momenta of the pointed particles, because we
* want to keep all of the information (that is the new momentum of a
* component after the splitting, which is contained in the _momentum
* member of the Component class, and the (old) momentum of that component
* before the splitting, which is contained in the momentum of the
* pointed particle). Please not make confusion of this only apparent
* inconsistency!
********************/
LorentzPoint posC,pos1,pos2;
posC = cluster->vertex();
calculatePositions(pClu, posC, pClu1, pClu2, pos1, pos2);
cutType rval;
if(toHadron1) {
rval.first = produceHadron(toHadron1, newPtr1, pClu1, pos1);
finalhadrons.push_back(rval.first.first);
}
else {
rval.first = produceCluster(ptrQ1, newPtr1, pClu1, pos1, pQ1, pQone, rem1);
}
if(toHadron2) {
rval.second = produceHadron(toHadron2, newPtr2, pClu2, pos2);
finalhadrons.push_back(rval.second.first);
}
else {
rval.second = produceCluster(ptrQ2, newPtr2, pClu2, pos2, pQ2, pQtwo, rem2);
}
return rval;
}
ClusterFissioner::cutType
ClusterFissioner::cutThree(ClusterPtr & cluster, tPVector & finalhadrons,
bool softUEisOn) {
// need to make sure only 3-cpt clusters get here
assert(cluster->numComponents() == 3);
// extract quarks
tPPtr ptrQ[3] = {cluster->particle(0),cluster->particle(1),cluster->particle(2)};
assert( ptrQ[0] && ptrQ[1] && ptrQ[2] );
// find maximum mass pair
Energy mmax(ZERO);
Lorentz5Momentum pDiQuark;
int iq1(-1),iq2(-1);
Lorentz5Momentum psum;
for(int q1=0;q1<3;++q1) {
psum+= ptrQ[q1]->momentum();
for(int q2=q1+1;q2<3;++q2) {
Lorentz5Momentum ptest = ptrQ[q1]->momentum()+ptrQ[q2]->momentum();
ptest.rescaleMass();
Energy mass = ptest.m();
if(mass>mmax) {
mmax = mass;
pDiQuark = ptest;
iq1 = q1;
iq2 = q2;
}
}
}
// and the spectators
int iother(-1);
for(int ix=0;ix<3;++ix) if(ix!=iq1&&ix!=iq2) iother=ix;
assert(iq1>=0&&iq2>=0&&iother>=0);
// And check if those particles are from a beam remnant
bool rem1 = cluster->isBeamRemnant(iq1);
bool rem2 = cluster->isBeamRemnant(iq2);
// workout which distribution to use
bool soft1(false),soft2(false);
switch (_iopRem) {
case 0:
soft1 = rem1 || rem2;
soft2 = rem2 || rem1;
break;
case 1:
soft1 = rem1;
soft2 = rem2;
break;
}
// Initialization for the exponential ("soft") mass distribution.
static const int max_loop = 1000;
int counter = 0;
Energy Mc1 = ZERO, Mc2 = ZERO, m1=ZERO, m2=ZERO, m=ZERO;
tcPDPtr toHadron;
bool toDiQuark(false);
PPtr newPtr1 = PPtr(),newPtr2 = PPtr();
PDPtr diquark;
bool succeeded = false;
do {
succeeded = false;
++counter;
if (_enhanceSProb == 0) {
drawNewFlavour(newPtr1,newPtr2);
}
else {
Energy2 mass2 = clustermass(cluster);
drawNewFlavourEnhanced(newPtr1,newPtr2, mass2);
}
// randomly pick which will be (anti)diquark and which a mesonic cluster
if(UseRandom::rndbool()) {
swap(iq1,iq2);
swap(rem1,rem2);
}
// check first order
if(cantMakeHadron(ptrQ[iq1], newPtr1) || cantMakeDiQuark(ptrQ[iq2], newPtr2)) {
swap(newPtr1,newPtr2);
}
// check again
if(cantMakeHadron(ptrQ[iq1], newPtr1) || cantMakeDiQuark(ptrQ[iq2], newPtr2)) {
throw Exception()
<< "ClusterFissioner cannot split the cluster ("
<< ptrQ[iq1]->PDGName() << ' ' << ptrQ[iq2]->PDGName()
<< ") into a hadron and diquark.\n" << Exception::runerror;
}
// Check that new clusters can produce particles and there is enough
// phase space to choose the drawn flavour
m1 = ptrQ[iq1]->data().constituentMass();
m2 = ptrQ[iq2]->data().constituentMass();
m = newPtr1->data().constituentMass();
// Do not split in the case there is no phase space available
if(mmax < m1+m + m2+m) continue;
// power for splitting
double exp1(_pSplitLight),exp2(_pSplitLight);
if (CheckId::isExotic (ptrQ[iq1]->dataPtr())) exp1 = _pSplitExotic;
else if(CheckId::hasBottom(ptrQ[iq1]->dataPtr())) exp1 = _pSplitBottom;
else if(CheckId::hasCharm (ptrQ[iq1]->dataPtr())) exp1 = _pSplitCharm;
if (CheckId::isExotic (ptrQ[iq2]->dataPtr())) exp2 = _pSplitExotic;
else if(CheckId::hasBottom(ptrQ[iq2]->dataPtr())) exp2 = _pSplitBottom;
else if(CheckId::hasCharm (ptrQ[iq2]->dataPtr())) exp2 = _pSplitCharm;
// If, during the drawing of candidate masses, too many attempts fail
// (because the phase space available is tiny)
/// \todo run separate loop here?
Mc1 = drawChildMass(mmax,m1,m2,m,exp1,soft1);
Mc2 = drawChildMass(mmax,m2,m1,m,exp2,soft2);
if(Mc1 < m1+m || Mc2 < m+m2 || Mc1+Mc2 > mmax) continue;
// check if need to force meson clster to hadron
toHadron = _hadronsSelector->chooseSingleHadron(ptrQ[iq1]->dataPtr(), newPtr1->dataPtr(),Mc1);
if(toHadron) Mc1 = toHadron->mass();
// check if need to force diquark cluster to be on-shell
toDiQuark = false;
- long rndSpin = UseRandom::rnd() > 0.5 ? 1 : 3;
+ long rndSpin;
+ if(ptrQ[iq2]->dataPtr()->id() == newPtr2->dataPtr()->id()) rndSpin = 3;
+ else rndSpin = UseRandom::rnd() > 0.5 ? 1 : 3;
diquark = CheckId::makeDiquark(ptrQ[iq2]->dataPtr(),newPtr2->dataPtr(),rndSpin);
if(Mc2 < diquark->constituentMass()) {
Mc2 = diquark->constituentMass();
toDiQuark = true;
}
// if a beam cluster not allowed to decay to hadrons
if(cluster->isBeamCluster() && toHadron && softUEisOn)
continue;
// Check if the decay kinematics is still possible: if not then
// force the one-hadron decay for the other cluster as well.
if(Mc1 + Mc2 > mmax) {
if(!toHadron) {
toHadron = _hadronsSelector->chooseSingleHadron(ptrQ[iq1]->dataPtr(), newPtr1->dataPtr(),mmax-Mc2);
if(toHadron) Mc1 = toHadron->mass();
}
else if(!toDiQuark) {
Mc2 = _hadronsSelector->massLightestHadron(ptrQ[iq2]->dataPtr(), newPtr2->dataPtr());
toDiQuark = true;
}
}
succeeded = (mmax >= Mc1+Mc2);
}
while (!succeeded && counter < max_loop);
// check no of tries
if(counter >= max_loop) return cutType();
// Determine the (5-components) momenta (all in the LAB frame)
Lorentz5Momentum p0Q1 = ptrQ[iq1]->momentum();
// to be determined
Lorentz5Momentum pClu1(Mc1), pClu2(Mc2), pQ1(m1), pQone(m), pQtwo(m), pQ2(m2);
calculateKinematics(pDiQuark,p0Q1,toHadron,toDiQuark,
pClu1,pClu2,pQ1,pQone,pQtwo,pQ2);
// positions of the new clusters
LorentzPoint pos1,pos2;
Lorentz5Momentum pBaryon = pClu2+ptrQ[iother]->momentum();
calculatePositions(cluster->momentum(), cluster->vertex(), pClu1, pBaryon, pos1, pos2);
// first the mesonic cluster/meson
cutType rval;
if(toHadron) {
rval.first = produceHadron(toHadron, newPtr1, pClu1, pos1);
finalhadrons.push_back(rval.first.first);
}
else {
rval.first = produceCluster(ptrQ[iq1], newPtr1, pClu1, pos1, pQ1, pQone, rem1);
}
if(toDiQuark) {
rem2 |= cluster->isBeamRemnant(iother);
PPtr newDiQuark = diquark->produceParticle(pClu2);
rval.second = produceCluster(newDiQuark, ptrQ[iother], pBaryon, pos2, pClu2,
ptrQ[iother]->momentum(), rem2);
}
else {
rval.second = produceCluster(ptrQ[iq2], newPtr2, pBaryon, pos2, pQ2, pQtwo, rem2,
ptrQ[iother],cluster->isBeamRemnant(iother));
}
cluster->isAvailable(false);
return rval;
}
ClusterFissioner::PPair
ClusterFissioner::produceHadron(tcPDPtr hadron, tPPtr newPtr, const Lorentz5Momentum &a,
const LorentzPoint &b) const {
PPair rval;
if(hadron->coloured()) {
rval.first = (_hadronsSelector->lightestHadron(hadron,newPtr->dataPtr()))->produceParticle();
}
else
rval.first = hadron->produceParticle();
rval.second = newPtr;
rval.first->set5Momentum(a);
rval.first->setVertex(b);
return rval;
}
ClusterFissioner::PPair ClusterFissioner::produceCluster(tPPtr ptrQ, tPPtr newPtr,
const Lorentz5Momentum & a,
const LorentzPoint & b,
const Lorentz5Momentum & c,
const Lorentz5Momentum & d,
bool isRem,
tPPtr spect, bool remSpect) const {
PPair rval;
rval.second = newPtr;
ClusterPtr cluster = !spect ? new_ptr(Cluster(ptrQ,rval.second)) : new_ptr(Cluster(ptrQ,rval.second,spect));
rval.first = cluster;
cluster->set5Momentum(a);
cluster->setVertex(b);
assert(cluster->particle(0)->id() == ptrQ->id());
cluster->particle(0)->set5Momentum(c);
cluster->particle(1)->set5Momentum(d);
cluster->setBeamRemnant(0,isRem);
if(remSpect) cluster->setBeamRemnant(2,remSpect);
return rval;
}
void ClusterFissioner::drawNewFlavour(PPtr& newPtrPos,PPtr& newPtrNeg) const {
// Flavour is assumed to be only u, d, s, with weights
// (which are not normalized probabilities) given
// by the same weights as used in HadronsSelector for
// the decay of clusters into two hadrons.
double prob_d;
double prob_u;
double prob_s;
switch(_fissionCluster){
case 0:
prob_d = _hadronsSelector->pwtDquark();
prob_u = _hadronsSelector->pwtUquark();
prob_s = _hadronsSelector->pwtSquark();
break;
case 1:
prob_d = _fissionPwtDquark;
prob_u = _fissionPwtUquark;
prob_s = _fissionPwtSquark;
break;
default :
assert(false);
}
int choice = UseRandom::rnd3(prob_u, prob_d, prob_s);
long idNew = 0;
switch (choice) {
case 0: idNew = ThePEG::ParticleID::u; break;
case 1: idNew = ThePEG::ParticleID::d; break;
case 2: idNew = ThePEG::ParticleID::s; break;
}
newPtrPos = getParticle(idNew);
newPtrNeg = getParticle(-idNew);
assert (newPtrPos);
assert(newPtrNeg);
assert (newPtrPos->dataPtr());
assert(newPtrNeg->dataPtr());
}
void ClusterFissioner::drawNewFlavourEnhanced(PPtr& newPtrPos,PPtr& newPtrNeg,
Energy2 mass2) const {
// Flavour is assumed to be only u, d, s, with weights
// (which are not normalized probabilities) given
// by the same weights as used in HadronsSelector for
// the decay of clusters into two hadrons.
double prob_d;
double prob_u;
double prob_s = 0.;
double scale = abs(double(sqr(_m0Fission)/mass2));
// Choose which splitting weights you wish to use
switch(_fissionCluster){
// 0: ClusterFissioner and ClusterDecayer use the same weights
case 0:
prob_d = _hadronsSelector->pwtDquark();
prob_u = _hadronsSelector->pwtUquark();
/* Strangeness enhancement:
Case 1: probability scaling
Case 2: Exponential scaling
*/
if (_enhanceSProb == 1)
prob_s = (_maxScale < scale) ? 0. : pow(_hadronsSelector->pwtSquark(),scale);
else if (_enhanceSProb == 2)
prob_s = (_maxScale < scale) ? 0. : exp(-scale);
break;
/* 1: ClusterFissioner uses its own unique set of weights,
i.e. decoupled from ClusterDecayer */
case 1:
prob_d = _fissionPwtDquark;
prob_u = _fissionPwtUquark;
if (_enhanceSProb == 1)
prob_s = (_maxScale < scale) ? 0. : pow(_fissionPwtSquark,scale);
else if (_enhanceSProb == 2)
prob_s = (_maxScale < scale) ? 0. : exp(-scale);
break;
}
int choice = UseRandom::rnd3(prob_u, prob_d, prob_s);
long idNew = 0;
switch (choice) {
case 0: idNew = ThePEG::ParticleID::u; break;
case 1: idNew = ThePEG::ParticleID::d; break;
case 2: idNew = ThePEG::ParticleID::s; break;
}
newPtrPos = getParticle(idNew);
newPtrNeg = getParticle(-idNew);
assert (newPtrPos);
assert(newPtrNeg);
assert (newPtrPos->dataPtr());
assert(newPtrNeg->dataPtr());
}
Energy2 ClusterFissioner::clustermass(const ClusterPtr & cluster){
Lorentz5Momentum pIn = cluster->momentum();
Energy2 endpointmass2 = sqr(cluster->particle(0)->mass() +
cluster->particle(1)->mass());
Energy2 singletm2 = pIn.m2();
// Return either the cluster mass, or the lambda measure
return (_massMeasure == 0) ? singletm2 : singletm2 - endpointmass2;
}
Energy ClusterFissioner::drawChildMass(const Energy M, const Energy m1,
const Energy m2, const Energy m,
const double expt, const bool soft) const {
/***************************
* This method, given in input the cluster mass Mclu of an heavy cluster C,
* made of consituents of masses m1 and m2, draws the masses Mclu1 and Mclu2
* of, respectively, the children cluster C1, made of constituent masses m1
* and m, and cluster C2, of mass Mclu2 and made of constituent masses m2
* and m. The mass is extracted from one of the two following mass
* distributions:
* --- power-like ("normal" distribution)
* d(Prob) / d(M^exponent) = const
* where the exponent can be different from the two children C1 (exp1)
* and C2 (exponent2).
* --- exponential ("soft" distribution)
* d(Prob) / d(M^2) = exp(-b*M)
* where b = 2.0 / average.
* Such distributions are limited below by the masses of
* the constituents quarks, and above from the mass of decaying cluster C.
* The choice of which of the two mass distributions to use for each of the
* two cluster children is dictated by iRemnant (see below).
* If the number of attempts to extract a pair of mass values that are
* kinematically acceptable is above some fixed number (max_loop, see below)
* the method gives up and returns false; otherwise, when it succeeds, it
* returns true.
*
* These distributions have been modified from HERWIG:
* Before these were:
* Mclu1 = m1 + (Mclu - m1 - m2)*pow( rnd(), 1.0/exponent1 );
* The new one coded here is a more efficient version, same density
* but taking into account 'in phase space from' beforehand
***************************/
// hard cluster
if(!soft) {
return pow(UseRandom::rnd(pow((M-m1-m2-m)*UnitRemoval::InvE, expt),
pow(m*UnitRemoval::InvE, expt)), 1./expt
)*UnitRemoval::E + m1;
}
// Otherwise it uses a soft mass distribution
else {
static const InvEnergy b = 2.0 / _btClM;
Energy max = M-m1-m2-2.0*m;
double rmin = b*max;
rmin = ( rmin < 50 ) ? exp(-rmin) : 0.;
double r1;
do {
r1 = UseRandom::rnd(rmin, 1.0) * UseRandom::rnd(rmin, 1.0);
}
while (r1 < rmin);
return m1 + m - log(r1)/b;
}
}
void ClusterFissioner::calculateKinematics(const Lorentz5Momentum & pClu,
const Lorentz5Momentum & p0Q1,
const bool toHadron1,
const bool toHadron2,
Lorentz5Momentum & pClu1,
Lorentz5Momentum & pClu2,
Lorentz5Momentum & pQ1,
Lorentz5Momentum & pQbar,
Lorentz5Momentum & pQ,
Lorentz5Momentum & pQ2bar) const {
/******************
* This method solves the kinematics of the two body cluster decay:
* C (Q1 Q2bar) ---> C1 (Q1 Qbar) + C2 (Q Q2bar)
* In input we receive the momentum of C, pClu, and the momentum
* of the quark Q1 (constituent of C), p0Q1, both in the LAB frame.
* Furthermore, two boolean variables inform whether the two fission
* products (C1, C2) decay immediately into a single hadron (in which
* case the cluster itself is identify with that hadron) and we do
* not have to solve the kinematics of the components (Q1,Qbar) for
* C1 and (Q,Q2bar) for C2.
* The output is given by the following momenta (all 5-components,
* and all in the LAB frame):
* pClu1 , pClu2 respectively of C1 , C2
* pQ1 , pQbar respectively of Q1 , Qbar in C1
* pQ , pQ2bar respectively of Q , Q2 in C2
* The assumption, suggested from the string model, is that, in C frame,
* C1 and its constituents Q1 and Qbar are collinear, and collinear to
* the direction of Q1 in C (that is before cluster decay); similarly,
* (always in the C frame) C2 and its constituents Q and Q2bar are
* collinear (and therefore anti-collinear with C1,Q1,Qbar).
* The solution is then obtained by using Lorentz boosts, as follows.
* The kinematics of C1 and C2 is solved in their parent C frame,
* and then boosted back in the LAB. The kinematics of Q1 and Qbar
* is solved in their parent C1 frame and then boosted back in the LAB;
* similarly, the kinematics of Q and Q2bar is solved in their parent
* C2 frame and then boosted back in the LAB. In each of the three
* "two-body decay"-like cases, we use the fact that the direction
* of the motion of the decay products is known in the rest frame of
* their parent. This is obvious for the first case in which the
* parent rest frame is C; but it is also true in the other two cases
* where the rest frames are C1 and C2. This is because C1 and C2
* are boosted w.r.t. C in the same direction where their components,
* respectively (Q1,Qbar) and (Q,Q2bar) move in C1 and C2 rest frame
* respectively.
* Of course, although the notation used assumed that C = (Q1 Q2bar)
* where Q1 is a quark and Q2bar an antiquark, indeed everything remain
* unchanged also in all following cases:
* Q1 quark, Q2bar antiquark; --> Q quark;
* Q1 antiquark , Q2bar quark; --> Q antiquark;
* Q1 quark, Q2bar diquark; --> Q quark
* Q1 antiquark, Q2bar anti-diquark; --> Q antiquark
* Q1 diquark, Q2bar quark --> Q antiquark
* Q1 anti-diquark, Q2bar antiquark; --> Q quark
**************************/
// Calculate the unit three-vector, in the C frame, along which
// all of the constituents and children clusters move.
Lorentz5Momentum u(p0Q1);
u.boost( -pClu.boostVector() ); // boost from LAB to C
// the unit three-vector is then u.vect().unit()
// Calculate the momenta of C1 and C2 in the (parent) C frame first,
// where the direction of C1 is u.vect().unit(), and then boost back in the
// LAB frame.
if (pClu.m() < pClu1.mass() + pClu2.mass() ) {
throw Exception() << "Impossible Kinematics in ClusterFissioner::calculateKinematics() (A)"
<< Exception::eventerror;
}
Kinematics::twoBodyDecay(pClu, pClu1.mass(), pClu2.mass(),
u.vect().unit(), pClu1, pClu2);
// In the case that cluster1 does not decay immediately into a single hadron,
// calculate the momenta of Q1 (as constituent of C1) and Qbar in the
// (parent) C1 frame first, where the direction of Q1 is u.vect().unit(),
// and then boost back in the LAB frame.
if(!toHadron1) {
if (pClu1.m() < pQ1.mass() + pQbar.mass() ) {
throw Exception() << "Impossible Kinematics in ClusterFissioner::calculateKinematics() (B)"
<< Exception::eventerror;
}
Kinematics::twoBodyDecay(pClu1, pQ1.mass(), pQbar.mass(),
u.vect().unit(), pQ1, pQbar);
}
// In the case that cluster2 does not decay immediately into a single hadron,
// Calculate the momenta of Q and Q2bar (as constituent of C2) in the
// (parent) C2 frame first, where the direction of Q is u.vect().unit(),
// and then boost back in the LAB frame.
if(!toHadron2) {
if (pClu2.m() < pQ.mass() + pQ2bar.mass() ) {
throw Exception() << "Impossible Kinematics in ClusterFissioner::calculateKinematics() (C)"
<< Exception::eventerror;
}
Kinematics::twoBodyDecay(pClu2, pQ.mass(), pQ2bar.mass(),
u.vect().unit(), pQ, pQ2bar);
}
}
void ClusterFissioner::calculatePositions(const Lorentz5Momentum & pClu,
const LorentzPoint & positionClu,
const Lorentz5Momentum & pClu1,
const Lorentz5Momentum & pClu2,
LorentzPoint & positionClu1,
LorentzPoint & positionClu2) const {
// Determine positions of cluster children.
// See Marc Smith's thesis, page 127, formulas (4.122) and (4.123).
Energy Mclu = pClu.m();
Energy Mclu1 = pClu1.m();
Energy Mclu2 = pClu2.m();
// Calculate the unit three-vector, in the C frame, along which
// children clusters move.
Lorentz5Momentum u(pClu1);
u.boost( -pClu.boostVector() ); // boost from LAB to C frame
// the unit three-vector is then u.vect().unit()
Energy pstarChild = Kinematics::pstarTwoBodyDecay(Mclu,Mclu1,Mclu2);
// First, determine the relative positions of the children clusters
// in the parent cluster reference frame.
Length x1 = ( 0.25*Mclu + 0.5*( pstarChild + (sqr(Mclu2) - sqr(Mclu1))/(2.0*Mclu)))/_kappa;
Length t1 = Mclu/_kappa - x1;
LorentzDistance distanceClu1( x1 * u.vect().unit(), t1 );
Length x2 = (-0.25*Mclu + 0.5*(-pstarChild + (sqr(Mclu2) - sqr(Mclu1))/(2.0*Mclu)))/_kappa;
Length t2 = Mclu/_kappa + x2;
LorentzDistance distanceClu2( x2 * u.vect().unit(), t2 );
// Then, transform such relative positions from the parent cluster
// reference frame to the Lab frame.
distanceClu1.boost( pClu.boostVector() );
distanceClu2.boost( pClu.boostVector() );
// Finally, determine the absolute positions in the Lab frame.
positionClu1 = positionClu + distanceClu1;
positionClu2 = positionClu + distanceClu2;
}
bool ClusterFissioner::isHeavy(tcClusterPtr clu) {
// default
double clpow = _clPowLight;
Energy clmax = _clMaxLight;
// particle data for constituents
tcPDPtr cptr[3]={tcPDPtr(),tcPDPtr(),tcPDPtr()};
for(int ix=0;ix<min(clu->numComponents(),3);++ix) {
cptr[ix]=clu->particle(ix)->dataPtr();
}
// different parameters for exotic, bottom and charm clusters
if(CheckId::isExotic(cptr[0],cptr[1],cptr[1])) {
clpow = _clPowExotic;
clmax = _clMaxExotic;
}
else if(CheckId::hasBottom(cptr[0],cptr[1],cptr[1])) {
clpow = _clPowBottom;
clmax = _clMaxBottom;
}
else if(CheckId::hasCharm(cptr[0],cptr[1],cptr[1])) {
clpow = _clPowCharm;
clmax = _clMaxCharm;
}
bool aboveCutoff = (
pow(clu->mass()*UnitRemoval::InvE , clpow)
>
pow(clmax*UnitRemoval::InvE, clpow)
+ pow(clu->sumConstituentMasses()*UnitRemoval::InvE, clpow)
);
// required test for SUSY clusters, since aboveCutoff alone
// cannot guarantee (Mc > m1 + m2 + 2*m) in cut()
static const Energy minmass
= getParticleData(ParticleID::d)->constituentMass();
bool canSplitMinimally
= clu->mass() > clu->sumConstituentMasses() + 2.0 * minmass;
return aboveCutoff && canSplitMinimally;
}
diff --git a/Hadronization/HadronSelector.cc b/Hadronization/HadronSelector.cc
--- a/Hadronization/HadronSelector.cc
+++ b/Hadronization/HadronSelector.cc
@@ -1,1022 +1,1078 @@
// -*- C++ -*-
//
// HadronSelector.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the HadronSelector class.
//
#include "HadronSelector.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include <ThePEG/PDT/EnumParticles.h>
#include <ThePEG/Repository/EventGenerator.h>
#include <ThePEG/Repository/CurrentGenerator.h>
#include <ThePEG/Repository/Repository.h>
#include "CheckId.h"
#include <ThePEG/Utilities/DescribeClass.h>
using namespace Herwig;
DescribeAbstractClass<HadronSelector,Interfaced>
describeHadronSelector("Herwig::HadronSelector","Herwig.so");
namespace {
// // debug helper
// void dumpTable(const HadronSelector::HadronTable & tbl) {
// typedef HadronSelector::HadronTable::const_iterator TableIter;
// for (TableIter it = tbl.begin(); it != tbl.end(); ++it) {
// cerr << it->first.first << ' '
// << it->first.second << '\n';
// for (HadronSelector::KupcoData::const_iterator jt = it->second.begin();
// jt != it->second.end(); ++jt) {
// cerr << '\t' << *jt << '\n';
// }
// }
// }
bool weightIsLess (pair<long,double> a, pair<long,double> b) {
return a.second < b.second;
}
}
ostream & operator<< (ostream & os,
const HadronSelector::HadronInfo & hi ) {
os << std::scientific << std::showpoint
<< std::setprecision(4)
<< setw(2)
<< hi.id << '\t'
// << hi.ptrData << ' '
<< hi.swtef << '\t'
<< hi.wt << '\t'
<< hi.overallWeight << '\t'
<< ounit(hi.mass,GeV);
return os;
}
HadronSelector::HadronSelector(unsigned int opt)
: _pwtDquark( 1.0 ),_pwtUquark( 1.0 ),_pwtSquark( 1.0 ),_pwtCquark( 0.0 ),
_pwtBquark( 0.0 ),_pwtDIquarkS0( 1.0 ),_pwtDIquarkS1( 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.),
_repwt(Lmax,vector<vector<double> >(Jmax,vector<double>(Nmax))),
_sngWt( 1.0 ),_decWt( 1.0 ),
_topt(opt),_trial(0),
_limBottom(), _limCharm(), _limExotic(), belowThreshold_(0)
{
// The mixing angles
// the ideal mixing angle
const double idealAngleMix = atan( sqrt(0.5) ) * 180.0 / Constants::pi;
// \eta-\eta' mixing angle
_etamix = -23.0;
// phi-omega mixing angle
_phimix = +36.0;
// h_1'-h_1 mixing angle
_h1mix = idealAngleMix;
// f_0(1710)-f_0(1370) mixing angle
_f0mix = idealAngleMix;
// f_1(1420)-f_1(1285)\f$ mixing angle
_f1mix = idealAngleMix;
// f'_2-f_2\f$ mixing angle
_f2mix = +26.0;
// eta_2(1870)-eta_2(1645) mixing angle
_eta2mix = idealAngleMix;
// phi(???)-omega(1650) mixing angle
_omhmix = idealAngleMix;
// phi_3-omega_3 mixing angle
_ph3mix = +28.0;
// eta(1475)-eta(1295) mixing angle
_eta2Smix = idealAngleMix;
// phi(1680)-omega(1420) mixing angle
_phi2Smix = idealAngleMix;
}
void HadronSelector::persistentOutput(PersistentOStream & os) const {
os << _partons << _pwtDquark << _pwtUquark << _pwtSquark
<< _pwtCquark << _pwtBquark << _pwtDIquarkS0 << _pwtDIquarkS1
<< _etamix << _phimix << _h1mix << _f0mix << _f1mix << _f2mix
<< _eta2mix << _omhmix << _ph3mix << _eta2Smix << _phi2Smix
<< _weight1S0 << _weight3S1 << _weight1P1 << _weight3P0 << _weight3P1
<< _weight3P2 << _weight1D2 << _weight3D1 << _weight3D2 << _weight3D3
<< _forbidden << _sngWt << _decWt << _repwt << _pwt
<< _limBottom << _limCharm << _limExotic << belowThreshold_
<< _table;
}
void HadronSelector::persistentInput(PersistentIStream & is, int) {
is >> _partons >> _pwtDquark >> _pwtUquark >> _pwtSquark
>> _pwtCquark >> _pwtBquark >> _pwtDIquarkS0 >> _pwtDIquarkS1
>> _etamix >> _phimix >> _h1mix >> _f0mix >> _f1mix >> _f2mix
>> _eta2mix >> _omhmix >> _ph3mix >> _eta2Smix >> _phi2Smix
>> _weight1S0 >> _weight3S1 >> _weight1P1 >> _weight3P0 >> _weight3P1
>> _weight3P2 >> _weight1D2 >> _weight3D1 >> _weight3D2 >> _weight3D3
>> _forbidden >> _sngWt >> _decWt >> _repwt >> _pwt
>> _limBottom >> _limCharm >> _limExotic >> belowThreshold_
>> _table;
}
void HadronSelector::Init() {
static ClassDocumentation<HadronSelector> documentation
("There is no documentation for the HadronSelector class");
static Parameter<HadronSelector,double>
interfacePwtDquark("PwtDquark","Weight for choosing a quark D",
&HadronSelector::_pwtDquark, 0, 1.0, 0.0, 10.0,
false,false,false);
static Parameter<HadronSelector,double>
interfacePwtUquark("PwtUquark","Weight for choosing a quark U",
&HadronSelector::_pwtUquark, 0, 1.0, 0.0, 10.0,
false,false,false);
static Parameter<HadronSelector,double>
interfacePwtSquark("PwtSquark","Weight for choosing a quark S",
&HadronSelector::_pwtSquark, 0, 1.0, 0.0, 10.0,
false,false,false);
static Parameter<HadronSelector,double>
interfacePwtCquark("PwtCquark","Weight for choosing a quark C",
&HadronSelector::_pwtCquark, 0, 0.0, 0.0, 10.0,
false,false,false);
static Parameter<HadronSelector,double>
interfacePwtBquark("PwtBquark","Weight for choosing a quark B",
&HadronSelector::_pwtBquark, 0, 0.0, 0.0, 10.0,
false,false,false);
static Parameter<HadronSelector,double>
interfacePwtDIquarkS0("PwtDIquarkS0","Weight for choosing a spin-0 DIquark",
&HadronSelector::_pwtDIquarkS0, 0, 1.0, 0.0, 100.0,
false,false,false);
static Parameter<HadronSelector,double>
interfacePwtDIquarkS1("PwtDIquarkS1","Weight for choosing a spin-1 DIquark",
&HadronSelector::_pwtDIquarkS1, 0, 1.0, 0.0, 100.0,
false,false,false);
static Parameter<HadronSelector,double>
interfaceSngWt("SngWt","Weight for singlet baryons",
&HadronSelector::_sngWt, 0, 1.0, 0.0, 10.0,
false,false,false);
static Parameter<HadronSelector,double>
interfaceDecWt("DecWt","Weight for decuplet baryons",
&HadronSelector::_decWt, 0, 1.0, 0.0, 10.0,
false,false,false);
static RefVector<HadronSelector,ParticleData> interfacePartons
("Partons",
"The partons which are to be considered as the consistuents of the hadrons.",
&HadronSelector::_partons, -1, false, false, true, false, false);
static RefVector<HadronSelector,ParticleData> interfaceForbidden
("Forbidden",
"The PDG codes of the particles which cannot be produced in the hadronization.",
&HadronSelector::_forbidden, -1, false, false, true, false, false);
//
// mixing angles
//
// the ideal mixing angle
const double idealAngleMix = atan( sqrt(0.5) ) * 180.0 / Constants::pi;
static Parameter<HadronSelector,double> interface11S0Mixing
("11S0Mixing",
"The mixing angle for the I=0 mesons from the 1 1S0 multiplet,"
" i.e. eta and etaprime.",
&HadronSelector::_etamix, -23., -180., 180.,
false, false, Interface::limited);
static Parameter<HadronSelector,double> interface13S1Mixing
("13S1Mixing",
"The mixing angle for the I=0 mesons from the 1 3S1 multiplet,"
" i.e. phi and omega.",
&HadronSelector::_phimix, +36., -180., 180.,
false, false, Interface::limited);
static Parameter<HadronSelector,double> interface11P1Mixing
("11P1Mixing",
"The mixing angle for the I=0 mesons from the 1 1P1 multiplet,"
" i.e. h_1' and h_1.",
&HadronSelector::_h1mix, idealAngleMix, -180., 180.,
false, false, Interface::limited);
static Parameter<HadronSelector,double> interface13P0Mixing
("13P0Mixing",
"The mixing angle for the I=0 mesons from the 1 3P0 multiplet,"
" i.e. f_0(1710) and f_0(1370).",
&HadronSelector::_f0mix, idealAngleMix, -180., 180.,
false, false, Interface::limited);
static Parameter<HadronSelector,double> interface13P1Mixing
("13P1Mixing",
"The mixing angle for the I=0 mesons from the 1 3P1 multiplet,"
" i.e. f_1(1420) and f_1(1285).",
&HadronSelector::_f1mix, idealAngleMix, -180., 180.,
false, false, Interface::limited);
static Parameter<HadronSelector,double> interface13P2Mixing
("13P2Mixing",
"The mixing angle for the I=0 mesons from the 1 3P2 multiplet,"
" i.e. f'_2 and f_2.",
&HadronSelector::_f2mix, 26.0, -180., 180.,
false, false, Interface::limited);
static Parameter<HadronSelector,double> interface11D2Mixing
("11D2Mixing",
"The mixing angle for the I=0 mesons from the 1 1D2 multiplet,"
" i.e. eta_2(1870) and eta_2(1645).",
&HadronSelector::_eta2mix, idealAngleMix, -180., 180.,
false, false, Interface::limited);
static Parameter<HadronSelector,double> interface13D0Mixing
("13D0Mixing",
"The mixing angle for the I=0 mesons from the 1 3D0 multiplet,"
" i.e. eta_2(1870) phi(?) and omega(1650).",
&HadronSelector::_omhmix, idealAngleMix, -180., 180.,
false, false, Interface::limited);
static Parameter<HadronSelector,double> interface13D1Mixing
("13D1Mixing",
"The mixing angle for the I=0 mesons from the 1 3D1 multiplet,"
" i.e. phi_3 and omega_3.",
&HadronSelector::_ph3mix, 28.0, -180., 180.,
false, false, Interface::limited);
static Parameter<HadronSelector,double> interface21S0Mixing
("21S0Mixing",
"The mixing angle for the I=0 mesons from the 2 1S0 multiplet,"
" i.e. eta(1475) and eta(1295).",
&HadronSelector::_eta2Smix, idealAngleMix, -180., 180.,
false, false, Interface::limited);
static Parameter<HadronSelector,double> interface23S1Mixing
("23S1Mixing",
"The mixing angle for the I=0 mesons from the 1 3S1 multiplet,"
" i.e. phi(1680) and omega(1420).",
&HadronSelector::_phi2Smix, idealAngleMix, -180., 180.,
false, false, Interface::limited);
//
// the meson weights
//
static ParVector<HadronSelector,double> interface1S0Weights
("1S0Weights",
"The weights for the 1S0 multiplets start with n=1.",
&HadronSelector::_weight1S0, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<HadronSelector,double> interface3S1Weights
("3S1Weights",
"The weights for the 3S1 multiplets start with n=1.",
&HadronSelector::_weight3S1, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<HadronSelector,double> interface1P1Weights
("1P1Weights",
"The weights for the 1P1 multiplets start with n=1.",
&HadronSelector::_weight1P1, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<HadronSelector,double> interface3P0Weights
("3P0Weights",
"The weights for the 3P0 multiplets start with n=1.",
&HadronSelector::_weight3P0, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<HadronSelector,double> interface3P1Weights
("3P1Weights",
"The weights for the 3P1 multiplets start with n=1.",
&HadronSelector::_weight3P1, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<HadronSelector,double> interface3P2Weights
("3P2Weights",
"The weights for the 3P2 multiplets start with n=1.",
&HadronSelector::_weight3P2, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<HadronSelector,double> interface1D2Weights
("1D2Weights",
"The weights for the 1D2 multiplets start with n=1.",
&HadronSelector::_weight1D2, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<HadronSelector,double> interface3D1Weights
("3D1Weights",
"The weights for the 3D1 multiplets start with n=1.",
&HadronSelector::_weight3D1, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<HadronSelector,double> interface3D2Weights
("3D2Weights",
"The weights for the 3D2 multiplets start with n=1.",
&HadronSelector::_weight3D2, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static ParVector<HadronSelector,double> interface3D3Weights
("3D3Weights",
"The weights for the 3D3 multiplets start with n=1.",
&HadronSelector::_weight3D3, Nmax, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static Switch<HadronSelector,unsigned int> interfaceTrial
("Trial",
"A Debugging option to only produce certain types of hadrons",
&HadronSelector::_trial, 0, false, false);
static SwitchOption interfaceTrialAll
(interfaceTrial,
"All",
"Produce all the hadrons",
0);
static SwitchOption interfaceTrialPions
(interfaceTrial,
"Pions",
"Only produce pions",
1);
static SwitchOption interfaceTrialSpin2
(interfaceTrial,
"Spin2",
"Only mesons with spin less than or equal to two are produced",
2);
static SwitchOption interfaceTrialSpin3
(interfaceTrial,
"Spin3",
"Only hadrons with spin less than or equal to three are produced",
3);
static Parameter<HadronSelector,double>
interfaceSingleHadronLimitBottom ("SingleHadronLimitBottom",
"Threshold for one-hadron decay of b-cluster",
&HadronSelector::_limBottom,
0, 0.0, 0.0, 100.0,false,false,false);
static Parameter<HadronSelector,double>
interfaceSingleHadronLimitCharm ("SingleHadronLimitCharm",
"threshold for one-hadron decay of c-cluster",
&HadronSelector::_limCharm,
0, 0.0, 0.0, 100.0,false,false,false);
static Parameter<HadronSelector,double>
interfaceSingleHadronLimitExotic ("SingleHadronLimitExotic",
"threshold for one-hadron decay of exotic cluster",
&HadronSelector::_limExotic,
0, 0.0, 0.0, 100.0,false,false,false);
static Switch<HadronSelector,unsigned int> interfaceBelowThreshold
("BelowThreshold",
"Option fo the selection of the hadrons if the cluster is below the pair threshold",
&HadronSelector::belowThreshold_, 0, false, false);
static SwitchOption interfaceBelowThresholdLightest
(interfaceBelowThreshold,
"Lightest",
"Force cluster to decay to the lightest hadron with the appropriate flavours",
0);
static SwitchOption interfaceBelowThresholdAll
(interfaceBelowThreshold,
"All",
"Select from all the hadrons below the two hadron threshold according to their spin weights",
1);
}
double HadronSelector::mixingStateWeight(long id) const {
switch(id) {
case ParticleID::eta: return 0.5*probabilityMixing(_etamix ,1);
case ParticleID::etaprime: return 0.5*probabilityMixing(_etamix ,2);
case ParticleID::phi: return 0.5*probabilityMixing(_phimix ,1);
case ParticleID::omega: return 0.5*probabilityMixing(_phimix ,2);
case ParticleID::hprime_1: return 0.5*probabilityMixing(_h1mix ,1);
case ParticleID::h_1: return 0.5*probabilityMixing(_h1mix ,2);
case 10331: return 0.5*probabilityMixing(_f0mix ,1);
case 10221: return 0.5*probabilityMixing(_f0mix ,2);
case ParticleID::fprime_1: return 0.5*probabilityMixing(_f1mix ,1);
case ParticleID::f_1: return 0.5*probabilityMixing(_f1mix ,2);
case ParticleID::fprime_2: return 0.5*probabilityMixing(_f2mix ,1);
case ParticleID::f_2: return 0.5*probabilityMixing(_f2mix ,2);
case 10335: return 0.5*probabilityMixing(_eta2mix ,1);
case 10225: return 0.5*probabilityMixing(_eta2mix ,2);
// missing phi member of 13D1 should be here
case 30223: return 0.5*probabilityMixing(_omhmix ,2);
case 337: return 0.5*probabilityMixing(_ph3mix ,1);
case 227: return 0.5*probabilityMixing(_ph3mix ,2);
case 100331: return 0.5*probabilityMixing(_eta2mix ,1);
case 100221: return 0.5*probabilityMixing(_eta2mix ,2);
case 100333: return 0.5*probabilityMixing(_phi2Smix,1);
case 100223: return 0.5*probabilityMixing(_phi2Smix,2);
default: return 1./3.;
}
}
void HadronSelector::doinit() {
Interfaced::doinit();
// the default partons allowed
// the quarks
for ( int ix=1; ix<=5; ++ix ) {
_partons.push_back(getParticleData(ix));
}
// the diquarks
for(unsigned int ix=1;ix<=5;++ix) {
for(unsigned int iy=1; iy<=ix;++iy) {
_partons.push_back(getParticleData(CheckId::makeDiquarkID(ix,iy,long(3))));
if(ix!=iy)
_partons.push_back(getParticleData(CheckId::makeDiquarkID(ix,iy,long(1))));
}
}
// 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];
// weights for the different quarks etc
for(unsigned int ix=0; ix<_partons.size(); ++ix) {
_pwt[_partons[ix]->id()]=1.;
}
_pwt[1] = _pwtDquark;
_pwt[2] = _pwtUquark;
_pwt[3] = _pwtSquark;
_pwt[4] = _pwtCquark;
_pwt[5] = _pwtBquark;
_pwt[1103] = _pwtDIquarkS1 * _pwtDquark * _pwtDquark;
_pwt[2101] = 0.5 * _pwtDIquarkS0 * _pwtUquark * _pwtDquark;
_pwt[2103] = 0.5 * _pwtDIquarkS1 * _pwtUquark * _pwtDquark;
_pwt[2203] = _pwtDIquarkS1 * _pwtUquark * _pwtUquark;
_pwt[3101] = 0.5 * _pwtDIquarkS0 * _pwtSquark * _pwtDquark;
_pwt[3103] = 0.5 * _pwtDIquarkS1 * _pwtSquark * _pwtDquark;
_pwt[3201] = 0.5 * _pwtDIquarkS0 * _pwtSquark * _pwtUquark;
_pwt[3203] = 0.5 * _pwtDIquarkS1 * _pwtSquark * _pwtUquark;
_pwt[3303] = _pwtDIquarkS1 * _pwtSquark * _pwtSquark;
// Commenting out heavy di-quark weights
_pwt[4101] = 0.0;
_pwt[4201] = 0.0;
_pwt[4301] = 0.0;
_pwt[4403] = 0.0;
_pwt[5101] = 0.0;
_pwt[5201] = 0.0;
_pwt[5301] = 0.0;
_pwt[5401] = 0.0;
_pwt[5503] = 0.0;
// find the maximum
map<long,double>::iterator pit =
max_element(_pwt.begin(),_pwt.end(),weightIsLess);
const double pmax = pit->second;
for(pit=_pwt.begin(); pit!=_pwt.end(); ++pit) {
pit->second/=pmax;
}
// construct the hadron tables
constructHadronTable();
// for debugging
// dumpTable(table());
}
void HadronSelector::constructHadronTable() {
// initialise the table
_table.clear();
for(unsigned int ix=0; ix<_partons.size(); ++ix) {
for(unsigned int iy=0; iy<_partons.size(); ++iy) {
if (!(DiquarkMatcher::Check(_partons[ix]->id())
&& DiquarkMatcher::Check(_partons[iy]->id())))
_table[make_pair(_partons[ix]->id(),_partons[iy]->id())] = KupcoData();
}
}
// get the particles from the event generator
ParticleMap particles = generator()->particles();
// loop over the particles
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;
int flav1;
int flav2;
// Skip non-hadrons (susy particles, etc...)
if(x3 == 0 || x2 == 0) continue;
else if(x4 == 0) { // meson
flav1 = x2;
flav2 = x3;
}
else { // baryon
- long rndSpin = UseRandom::rnd() > 0.5 ? 1 : 3;
+ long rndSpin;
+ if(x2 == x3) rndSpin = 3;
+ else rndSpin = UseRandom::rnd() > 0.5 ? 1 : 3;
flav1 = CheckId::makeDiquarkID(x2,x3,rndSpin);
flav2 = x4;
}
if (wantSusy) flav2 += 1000000 * x7;
HadronInfo a(pid,
particle,
- specialWeight(pid),
+ specialWeight(pid,flav1),
particle->mass());
// set the weight to the number of spin states
a.overallWeight = pspin;
// identical light flavours
if(flav1 == flav2 && flav1<=3) {
// ddbar> uubar> admixture states
if(flav1==1) {
if(_topt != 0) a.overallWeight *= 0.5*a.swtef;
_table[make_pair(1,1)].insert(a);
_table[make_pair(2,2)].insert(a);
if(_topt == 0 && a.overallWeight > maxdd) maxdd = a.overallWeight;
}
// load up ssbar> uubar> ddbar> admixture states
else {
a.wt = mixingStateWeight(pid);
a.overallWeight *= a.wt;
if(_topt != 0) a.overallWeight *= a.swtef;
_table[make_pair(1,1)].insert(a);
_table[make_pair(2,2)].insert(a);
if(_topt == 0 && a.overallWeight > maxdd) maxdd = a.overallWeight;
a.wt = (_topt != 0) ? 1.- 2.*a.wt : 1 - a.wt;
if(a.wt > 0) {
a.overallWeight = a.wt * a.swtef * pspin;
_table[make_pair(3,3)].insert(a);
if(_topt == 0 && a.overallWeight > maxss) maxss = a.overallWeight;
}
}
}
// light baryons with all quarks identical
else if((flav1 == 1 && flav2 == 1103) || (flav1 == 1103 && flav2 == 1) ||
(flav1 == 2 && flav2 == 2203) || (flav1 == 2203 && flav2 == 2) ||
(flav1 == 3 && flav2 == 3303) || (flav1 == 3303 && flav2 == 3)) {
if(_topt != 0) a.overallWeight *= 1.5*a.swtef;
_table[make_pair(flav1,flav2)].insert(a);
_table[make_pair(flav2,flav1)].insert(a);
if(_topt == 0 && a.overallWeight > maxrest) maxrest = a.overallWeight;
}
// all other cases
else {
if(_topt != 0) a.overallWeight *=a.swtef;
_table[make_pair(flav1,flav2)].insert(a);
if(flav1 != flav2) _table[make_pair(flav2,flav1)].insert(a);
if(_topt == 0 && a.overallWeight > maxrest) maxrest = a.overallWeight;
}
}
// Account for identical combos of diquark/quarks and symmetrical elements
// e.g. U UD = D UU
HadronTable::iterator tit;
for(tit=_table.begin();tit!=_table.end();++tit) {
if(tit->first.first>ParticleID::c) continue;
if(!DiquarkMatcher::Check(tit->first.second)) continue;
long k, l, sub;
if(tit->first.second>=ParticleID::bd_0) {
k = ParticleID::b;
sub = ParticleID::bd_0/100;
}
else if(tit->first.second>=ParticleID::cd_0) {
k = ParticleID::c;
sub = ParticleID::cd_0/100;
}
else if(tit->first.second>=ParticleID::sd_0) {
k = ParticleID::s;
sub = ParticleID::sd_0/100;
}
else if(tit->first.second>=ParticleID::ud_0) {
k = ParticleID::u;
sub = ParticleID::ud_0/100;
}
else if(tit->first.second==ParticleID::dd_1) {
k = ParticleID::d;
sub = ParticleID::dd_1/100;
}
else continue;
sub=tit->first.second/100-sub+1;
if(sub > tit->first.first) {
l = 1000*sub+100*tit->first.first+1;
}
else if(sub==tit->first.first) {
l = 1000*sub+ 100*tit->first.first+3;
}
else {
l = 100*sub +1000*tit->first.first+1;
}
if(tit->second.empty()) {
pair<long,long> newpair(k,l);
tit->second=_table[newpair];
newpair=make_pair(tit->first.second,tit->first.first);
_table[newpair]=tit->second;
};
}
// normalise weights to one for first option
if(_topt == 0) {
HadronTable::const_iterator tit;
KupcoData::iterator it;
for(tit=_table.begin();tit!=_table.end();++tit) {
double weight;
if(tit->first.first==tit->first.second) {
if(tit->first.first==1||tit->first.first==2) weight=1./maxdd;
else if (tit->first.first==3) weight=1./maxss;
else weight=1./maxrest;
}
else weight=1./maxrest;
for(it = tit->second.begin(); it!=tit->second.end(); ++it) {
it->rescale(weight);
}
}
}
}
-double HadronSelector::specialWeight(long id) const {
+double HadronSelector::specialWeight(long id, int part1) const {
const int pspin = id % 10;
+ // Get the flavours
+ const int x4 = (id/1000)%10;
+ const int x3 = (id/100 )%10;
+ const int x2 = (id/10 )%10;
+ // Clebsch-Gordan coefficients of baryons
+ double CGcoefficient(1.);
// Only K0L and K0S have pspin == 0, should
// not get them until Decay step
assert( pspin != 0 );
// Baryon : J = 1/2 or 3/2
if(pspin == 2) {
// Singlet (Lambda-like) baryon
if( (id/100)%10 < (id/10 )%10 ) return sqr(_sngWt);
// octet
- else return 1.;
+ else {
+ // return 1 if part1 is not a diquark
+ if(part1 < 1100) return 1.;
+ // if part1 is a diquark
+ int diquarkSpin = (part1%10 == 1) ? 1 : 3;
+ // ud0+u
+ if(x4 != x3 && (x2 == x3 || x2 == x4) && diquarkSpin==1)
+ CGcoefficient = 3./4.;
+ // ud0+s
+ else if(x4 != x3 && (x2 != x3 || x2 != x4) && diquarkSpin==1)
+ CGcoefficient = 1./2.;
+ // uu1+u
+ else if(x4 == x3 && x4 == x2 && diquarkSpin==3)
+ CGcoefficient = 0.;
+ // uu1+d
+ else if(x4 == x3 && x4 != x2 && diquarkSpin==3)
+ CGcoefficient = 1./6.;
+ // ud1+u
+ else if(x4 != x3 && (x2 == x3 || x2 == x4) && diquarkSpin==3)
+ CGcoefficient = 1./12.;
+ // ud1+s
+ else if(x4 != x3 && (x2 != x3 || x2 != x4) && diquarkSpin==3)
+ CGcoefficient = 1./6.;
+
+ return CGcoefficient;
+ }
}
// Decuplet baryon
else if (pspin == 4) {
- return sqr(_decWt);
+ // return 1 if part1 is not a diquark
+ if(part1 < 1100) return sqr(_decWt);
+ // if part1 is a diquark
+ int diquarkSpin = (part1%10 == 1) ? 1 : 3;
+ // ud0+u
+ if(x4 != x3 && (x2 == x3 || x2 == x4) && diquarkSpin==1)
+ CGcoefficient = 0.;
+ // ud0+s
+ else if(x4 != x3 && (x2 != x3 || x2 != x4) && diquarkSpin==1)
+ CGcoefficient = 0.;
+ // uu1+u
+ else if(x4 == x3 && x4 == x2 && diquarkSpin==3)
+ CGcoefficient = 1.;
+ // uu1+d
+ else if(x4 == x3 && x4 != x2 && diquarkSpin==3)
+ CGcoefficient = 1./3.;
+ // ud1+u
+ else if(x4 != x3 && (x2 == x3 || x2 == x4) && diquarkSpin==3)
+ CGcoefficient = 2./3.;
+ // ud1+s
+ else if(x4 != x3 && (x2 != x3 || x2 != x4) && diquarkSpin==3)
+ CGcoefficient = 1./3.;
+
+ return sqr(_decWt)*CGcoefficient;
}
// Meson
else if(pspin % 2 == 1) {
// Total angular momentum
int j = (pspin - 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
return 1.0;
}
int HadronSelector::signHadron(tcPDPtr idQ1, tcPDPtr idQ2,
tcPDPtr hadron) const {
// This method receives in input three PDG ids, whose the
// first two have proper signs (corresponding to particles, id > 0,
// or antiparticles, id < 0 ), whereas the third one must
// be always positive (particle not antiparticle),
// corresponding to:
// --- quark-antiquark, or antiquark-quark, or
// quark-diquark, or diquark-quark, or
// antiquark-antidiquark, or antidiquark-antiquark
// for the first two input (idQ1, idQ2);
// --- meson or baryon for the third input (idHad):
// The method returns:
// --- + 1 if the two partons (idQ1, idQ2) are exactly
// the constituents for the hadron idHad;
// --- - 1 if the two partons (idQ1, idQ2) are exactly
// the constituents for the anti-hadron -idHad;
// --- + 0 otherwise.
// The method it is therefore useful to decide the
// sign of the id of the produced hadron as appeared
// in the vector _vecHad (where only hadron idHad > 0 are present)
// given the two constituent partons.
int sign = 0;
long idHad = hadron->id();
assert(idHad > 0);
int chargeIn = idQ1->iCharge() + idQ2->iCharge();
int chargeOut = hadron->iCharge();
// same charge
if( chargeIn == chargeOut && chargeIn !=0 ) sign = +1;
else if(chargeIn == -chargeOut && chargeIn !=0 ) sign = -1;
else if(chargeIn == 0 && chargeOut == 0 ) {
// In the case of same null charge, there are four cases:
// i) K0-like mesons, B0-like mesons, Bs-like mesons
// the PDG convention is to consider them "antiparticle" (idHad < 0)
// if the "dominant" (heavier) flavour (respectively, s, b)
// is a quark (idQ > 0): for instance, B0s = (b, sbar) has id < 0
// Remember that there is an important exception for K0L (id=130) and
// K0S (id=310): they don't have antiparticles, therefore idHad > 0
// always. We use below the fact that K0L and K0S are the unique
// hadrons having 0 the first (less significant) digit of their id.
// 2) D0-like mesons: the PDG convention is to consider them "particle"
// (idHad > 0) if the charm flavour is carried by a c: (c,ubar) has id>0
// 3) the remaining mesons should not have antiparticle, therefore their
// sign is always positive.
// 4) for baryons, that is when one of idQ1 and idQ2 is a (anti-) quark and
// the other one is a (anti-) diquark the sign is negative when both
// constituents are "anti", that is both with id < 0; positive otherwise.
// meson
if(abs(int(idQ1->iColour()))== 3 && abs(int(idQ2->iColour())) == 3 &&
!DiquarkMatcher::Check(idQ1->id()) && !DiquarkMatcher::Check(idQ2->id()))
{
int idQa = abs(idQ1->id());
int idQb = abs(idQ2->id());
int dominant = idQ2->id();
if(idQa > idQb) {
swap(idQa,idQb);
dominant = idQ1->id();
}
if((idQa==ParticleID::d && idQb==ParticleID::s) ||
(idQa==ParticleID::d && idQb==ParticleID::b) ||
(idQa==ParticleID::s && idQb==ParticleID::b)) {
// idHad%10 is zero for K0L,K0S
if (dominant < 0 || idHad%10 == 0) sign = +1;
else if(dominant > 0) sign = -1;
}
else if((idQa==ParticleID::u && idQb==ParticleID::c) ||
(idQa==ParticleID::u && idQb==ParticleID::t) ||
(idQa==ParticleID::c && idQb==ParticleID::t)) {
if (dominant > 0) sign = +1;
else if(dominant < 0) sign = -1;
}
else if(idQa==idQb) sign = +1;
// sets sign for Susy particles
else sign = (dominant > 0) ? +1 : -1;
}
// baryon
else if(DiquarkMatcher::Check(idQ1->id()) || DiquarkMatcher::Check(idQ2->id())) {
if (idQ1->id() > 0 && idQ2->id() > 0) sign = +1;
else if(idQ1->id() < 0 && idQ2->id() < 0) sign = -1;
}
}
if (sign == 0) {
cerr << "Could not work out sign for "
<< idQ1->PDGName() << ' '
<< idQ2->PDGName() << " => "
<< hadron->PDGName() << '\n';
assert(false);
}
return sign;
}
pair<tcPDPtr,tcPDPtr> HadronSelector::lightestHadronPair(tcPDPtr ptr1, tcPDPtr ptr2,
tcPDPtr ptr3) const {
// throw exception of id3!=0 as doesn't work
if ( ptr3 ) throw Exception()
<< "ptr3!=0 not yet implemented in HadronSelector::lightestHadronPair"
<< Exception::abortnow;
// charge
int totalcharge = ptr1->iCharge() + ptr2->iCharge();
if ( ptr3 ) totalcharge += ptr3->iCharge();
tcPDPtr vIdHad1[2]={tcPDPtr(),tcPDPtr()},vIdHad2[2]={tcPDPtr(),tcPDPtr()};
bool vOk[2] = {false, false};
Energy vMassPair[2] = { ZERO, ZERO };
for (int i = 0; i < 2; i++) {
tcPDPtr idPartner = i==0 ? getParticleData(ParticleID::d) : getParticleData(ParticleID::u);
// Change sign to idPartner (transform it into a anti-quark) if it is not
// possible to form a meson or a baryon.
assert (ptr1 && idPartner);
if (!CheckId::canBeHadron(ptr1, idPartner)) idPartner = idPartner->CC();
vIdHad1[i] = lightestHadron(ptr1, idPartner);
vIdHad2[i] = lightestHadron(ptr2, idPartner->CC());
if ( vIdHad1[i] && vIdHad2[i] &&
vIdHad1[i]->iCharge() + vIdHad2[i]->iCharge() == totalcharge ) {
vOk[i] = true;
vMassPair[i] = vIdHad1[i]->mass() + vIdHad2[i]->mass();
}
}
// Take the lightest pair compatible with charge conservation.
if ( vOk[0] && ( ! vOk[1] || vMassPair[0] <= vMassPair[1] ) ) {
return make_pair(vIdHad1[0],vIdHad2[0]);
}
else if ( vOk[1] && ( ! vOk[0] || vMassPair[1] < vMassPair[0] ) ) {
return make_pair(vIdHad1[1],vIdHad2[1]);
}
else {
return make_pair(tcPDPtr(),tcPDPtr());
}
}
Energy HadronSelector::massLightestBaryonPair(tcPDPtr ptr1, tcPDPtr ptr2) const {
// Make sure that we don't have any diquarks as input, return arbitrarily
// large value if we do
Energy currentSum = Constants::MaxEnergy;
for(unsigned int ix=0; ix<_partons.size(); ++ix) {
if(!DiquarkMatcher::Check(_partons[ix]->id())) continue;
HadronTable::const_iterator
tit1=_table.find(make_pair(abs(ptr1->id()),_partons[ix]->id())),
tit2=_table.find(make_pair(_partons[ix]->id(),abs(ptr2->id())));
if( tit1==_table.end() || tit2==_table.end()) continue;
if(tit1->second.empty()||tit2->second.empty()) continue;
Energy s = tit1->second.begin()->mass + tit2->second.begin()->mass;
if(currentSum > s) currentSum = s;
}
return currentSum;
}
tcPDPtr HadronSelector::lightestHadron(tcPDPtr ptr1, tcPDPtr ptr2,
#ifndef NDEBUG
tcPDPtr ptr3) const {
#else
tcPDPtr ) const {
#endif
// The method assumes ptr3 == 0 rest not implemented
assert(ptr1 && ptr2 && !ptr3);
// find entry in the table
pair<long,long> ids = make_pair(abs(ptr1->id()),abs(ptr2->id()));
HadronTable::const_iterator tit=_table.find(ids);
// throw exception if flavours wrong
if (tit==_table.end())
throw Exception() << "Could not find "
<< ids.first << ' ' << ids.second
<< " in _table. "
<< "In HadronSelector::lightestHadron()"
<< Exception::eventerror;
if(tit->second.empty())
throw Exception() << "HadronSelector::lightestHadron "
<< "could not find any hadrons containing "
<< ptr1->id() << ' ' << ptr2->id() << '\n'
<< tit->first.first << ' '
<< tit->first.second << Exception::eventerror;
// find the lightest hadron
int sign = signHadron(ptr1,ptr2,tit->second.begin()->ptrData);
tcPDPtr candidate = sign > 0 ?
tit->second.begin()->ptrData : tit->second.begin()->ptrData->CC();
// \todo 20 GeV limit is temporary fudge to let SM particles go through.
// \todo Use isExotic instead?
if (candidate->mass() > 20*GeV
&& candidate->mass() < ptr1->constituentMass() + ptr2->constituentMass()) {
generator()->log() << "HadronSelector::lightestHadron: "
<< "chosen candidate " << candidate->PDGName()
<< " is lighter than its constituents "
<< ptr1->PDGName() << ", " << ptr2->PDGName() << '\n'
<< candidate->mass()/GeV << " < " << ptr1->constituentMass()/GeV
<< " + " << ptr2->constituentMass()/GeV << '\n'
<< "Check your particle data tables.\n";
assert(false);
}
return candidate;
}
vector<pair<tcPDPtr,double> >
HadronSelector::hadronsBelowThreshold(Energy threshold, tcPDPtr ptr1,
tcPDPtr ptr2,
#ifndef NDEBUG
tcPDPtr ptr3) const {
#else
tcPDPtr ) const {
#endif
// The method assumes ptr3 == 0 rest not implemented
assert(ptr1 && ptr2 && !ptr3);
// find entry in the table
pair<long,long> ids = make_pair(abs(ptr1->id()),abs(ptr2->id()));
HadronTable::const_iterator tit=_table.find(ids);
// throw exception if flavours wrong
if (tit==_table.end())
throw Exception() << "Could not find "
<< ids.first << ' ' << ids.second
<< " in _table. "
<< "In HadronSelector::hadronsBelowThreshold()"
<< Exception::eventerror;
if(tit->second.empty())
throw Exception() << "HadronSelector::hadronsBelowThreshold() "
<< "could not find any hadrons containing "
<< ptr1->id() << ' ' << ptr2->id() << '\n'
<< tit->first.first << ' '
<< tit->first.second << Exception::eventerror;
vector<pair<tcPDPtr,double> > candidates;
KupcoData::const_iterator hit = tit->second.begin();
// find the hadrons
while(hit!=tit->second.end()&&hit->mass<threshold) {
// find the hadron
int sign = signHadron(ptr1,ptr2,hit->ptrData);
tcPDPtr candidate = sign > 0 ? hit->ptrData : hit->ptrData->CC();
// \todo 20 GeV limit is temporary fudge to let SM particles go through.
// \todo Use isExotic instead?
if (candidate->mass() > 20*GeV
&& candidate->mass() < ptr1->constituentMass() + ptr2->constituentMass()) {
generator()->log() << "HadronSelector::hadronsBelowTheshold: "
<< "chosen candidate " << candidate->PDGName()
<< " is lighter than its constituents "
<< ptr1->PDGName() << ", " << ptr2->PDGName() << '\n'
<< candidate->mass()/GeV << " < " << ptr1->constituentMass()/GeV
<< " + " << ptr2->constituentMass()/GeV << '\n'
<< "Check your particle data tables.\n";
assert(false);
}
candidates.push_back(make_pair(candidate,hit->overallWeight));
++hit;
}
return candidates;
}
tcPDPtr HadronSelector::chooseSingleHadron(tcPDPtr par1, tcPDPtr par2,
Energy mass) const {
// Determine the sum of the nominal masses of the two lightest hadrons
// with the right flavour numbers as the cluster under consideration.
// Notice that we don't need real masses (drawn by a Breit-Wigner
// distribution) because the lightest pair of hadrons does not involve
// any broad resonance.
Energy threshold = massLightestHadronPair(par1,par2);
// Special: it allows one-hadron decays also above threshold.
if (CheckId::isExotic(par1,par2))
threshold *= (1.0 + UseRandom::rnd()*_limExotic);
else if (CheckId::hasBottom(par1,par2))
threshold *= (1.0 + UseRandom::rnd()*_limBottom);
else if (CheckId::hasCharm(par1,par2))
threshold *= (1.0 + UseRandom::rnd()*_limCharm);
// only do one hadron decay is mass less than the threshold
if(mass>=threshold) return tcPDPtr();
// select the hadron
tcPDPtr hadron;
// old option pick the lightest hadron
if(belowThreshold_ == 0) {
hadron= lightestHadron(par1,par2);
}
// new option select from those available
else if(belowThreshold_ == 1) {
vector<pair<tcPDPtr,double> > hadrons =
hadronsBelowThreshold(threshold,par1,par2);
if(hadrons.size()==1) {
hadron = hadrons[0].first;
}
else if(hadrons.empty()) {
hadron= lightestHadron(par1,par2);
}
else {
double totalWeight=0.;
for(unsigned int ix=0;ix<hadrons.size();++ix) {
totalWeight += hadrons[ix].second;
}
totalWeight *= UseRandom::rnd();
for(unsigned int ix=0;ix<hadrons.size();++ix) {
if(totalWeight<=hadrons[ix].second) {
hadron = hadrons[ix].first;
break;
}
else
totalWeight -= hadrons[ix].second;
}
assert(hadron);
}
}
else
assert(false);
return hadron;
}
diff --git a/Hadronization/HadronSelector.h b/Hadronization/HadronSelector.h
--- a/Hadronization/HadronSelector.h
+++ b/Hadronization/HadronSelector.h
@@ -1,826 +1,827 @@
// -*- C++ -*-
//
// HadronSelector.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_HadronSelector_H
#define HERWIG_HadronSelector_H
//
// This is the declaration of the HadronSelector class.
//
#include "ThePEG/Interface/Interfaced.h"
#include <ThePEG/Persistency/PersistentOStream.h>
#include <ThePEG/Persistency/PersistentIStream.h>
#include <ThePEG/PDT/ParticleData.h>
#include <ThePEG/PDT/StandardMatchers.h>
#include <ThePEG/Repository/EventGenerator.h>
#include "HadronSelector.fh"
namespace Herwig {
using namespace ThePEG;
/**\ingroup Hadronization
* \class HadronSelector
* \brief This class selects the hadron flavours of a cluster decay.
* \author Philip Stephens
* \author Alberto Ribon
* \author Peter Richardson
*
* This is the base class for the selection of either a pair of hadrons, or
* in some cases a single hadron. The different approaches which were
* previously implemented in this class are now implemented in the
* HwppSelector and Hw64Selector which inherit from this class.
*
* This class implements a number of methods which are needed by all models
* and in addition contains the weights for the different meson multiplets and
* mixing of the light \f$I=0\f$ mesons.
*
* @see \ref HadronSelectorInterfaces "The interfaces"
* defined for HadronSelector.
* @see HwppSelector
* @see Hw64Selector
*/
class HadronSelector: public Interfaced {
public:
/** \ingroup Hadronization
* \class HadronInfo
* \brief Class used to store all the hadron information for easy access.
* \author Philip Stephens
*
* Note that:
* - the hadrons in _table can be filled in any ordered
* w.r.t. the mass value, and flavours for different
* groups (for instance, (u,s) hadrons don't need to
* be placed after (d,s) or any other flavour), but
* all hadrons with the same flavours must be consecutive
* ( for instance you cannot alternate hadrons of type
* (d,s) with those of flavour (u,s) ).
* Furthermore, it is assumed that particle and antiparticle
* have the same weights, and therefore only one of them
* must be entered in the table: we have chosen to refer
* to the particle, defined as PDG id > 0, although if
* an anti-particle is provided in input it is automatically
* transform to its particle, simply by taking the modulus
* of its id.
*/
class HadronInfo {
public:
/**
* Constructor
* @param idin The PDG code of the hadron
* @param datain The pointer to the ParticleData object
* @param swtin The singlet/decuplet/orbital factor
* @param massin The mass of the hadron
*/
HadronInfo(long idin=0, tPDPtr datain=tPDPtr(),
double swtin=1., Energy massin=ZERO)
: id(idin), ptrData(datain), swtef(swtin), wt(1.0), overallWeight(0.0),
mass(massin)
{}
/**
* Comparision operator on mass
*/
bool operator<(const HadronInfo &x) const {
if(mass!=x.mass) return mass < x.mass;
else return id < x.id;
}
/**
* The hadrons id.
*/
long id;
/**
* pointer to ParticleData, to get the spin, etc...
*/
tPDPtr ptrData;
/**
* singlet/decuplet/orbital factor
*/
double swtef;
/**
* mixing factor
*/
double wt;
/**
* (2*J+1)*wt*swtef
*/
double overallWeight;
/**
* The hadrons mass
*/
Energy mass;
/**
* Rescale the weight for a given hadron
*/
void rescale(double x) const {
const_cast<HadronInfo*>(this)->overallWeight *= x;
}
/**
* Friend method used to print the value of a table element.
*/
friend PersistentOStream & operator<< (PersistentOStream & os,
const HadronInfo & hi ) {
os << hi.id << hi.ptrData << hi.swtef << hi.wt
<< hi.overallWeight << ounit(hi.mass,GeV);
return os;
}
/**
* debug output
*/
friend ostream & operator<< (ostream & os, const HadronInfo & hi );
/**
* Friend method used to read in the value of a table element.
*/
friend PersistentIStream & operator>> (PersistentIStream & is,
HadronInfo & hi ) {
is >> hi.id >> hi.ptrData >> hi.swtef >> hi.wt
>> hi.overallWeight >> iunit(hi.mass,GeV);
return is;
}
};
/** \ingroup Hadronization
* \class Kupco
* \brief Class designed to make STL routines easy to use.
* \author Philip Stephens
*
* This class is used to generate a list of the hadron pairs which can
* be produced that allows easy traversal and quick access.
*/
class Kupco {
public:
/**
* Constructor
* @param inidQ PDG code of the quark drawn from the vacuum.
* @param inhad1 ParticleData for the first hadron produced.
* @param inhad2 ParticleData for the second hadron produced.
* @param inwgt The weight for the hadron pair
*/
Kupco(tcPDPtr inidQ,tcPDPtr inhad1,tcPDPtr inhad2, Energy inwgt)
: idQ(inidQ),hadron1(inhad1),hadron2(inhad2),weight(inwgt)
{}
/**
* id of the quark drawn from the vacuum.
*/
tcPDPtr idQ;
/**
* The ParticleData object for the first hadron produced.
*/
tcPDPtr hadron1;
/**
* The ParticleData object for the second hadron produced.
*/
tcPDPtr hadron2;
/**
* Weight factor of this componation.
*/
Energy weight;
};
public:
/**
* The helper classes
*/
//@{
/**
* The type is used to contain all the hadrons info of a given flavour.
*/
typedef set<HadronInfo> KupcoData;
//@}
/**
* The hadron table type.
*/
typedef map<pair<long,long>,KupcoData> HadronTable;
public:
/**
* The default constructor.
*/
HadronSelector(unsigned int);
/**
* Method to return a pair of hadrons given the PDG codes of
* two or three constituents
* @param cluMass The mass of the cluster
* @param par1 The first constituent
* @param par2 The second constituent
* @param par3 The third constituent
*/
virtual pair<tcPDPtr,tcPDPtr> chooseHadronPair(const Energy cluMass, tcPDPtr par1,
tcPDPtr par2,tcPDPtr par3 = PDPtr()) const
= 0;
/**
* Select the single hadron for a cluster decay
* return null pointer if not a single hadron decay
* @param par1 1st constituent
* @param par2 2nd constituent
* @param mass Mass of the cluster
*/
tcPDPtr chooseSingleHadron(tcPDPtr par1, tcPDPtr par2, Energy mass) const;
/**
* This returns the lightest pair of hadrons given by the flavours.
*
* Given the two (or three) constituents of a cluster, it returns
* the two lightest hadrons with proper flavour numbers.
* Furthermore, the first of the two hadrons must have the constituent with
* par1, and the second must have the constituent with par2.
* \todo At the moment it does *nothing* in the case that also par3 is present.
*
* The method is implemented by calling twice lightestHadron,
* once with (par1,quarktopick->CC()) ,and once with (par2,quarktopick)
* where quarktopick is either the pointer to
* d or u quarks . In fact, the idea is that whatever the flavour of par1
* and par2, no matter if (anti-)quark or (anti-)diquark, the lightest
* pair of hadrons containing flavour par1 and par2 will have either
* flavour d or u, being the lightest quarks.
* The method returns the pair (PDPtr(),PDPtr()) if anything goes wrong.
*
* \todo The method assumes par3 == PDPtr() (otherwise we don't know how to proceed: a
* possible, trivial way would be to randomly select two of the three
* (anti-)quarks and treat them as a (anti-)diquark, reducing the problem
* to two components as treated below.
* In the normal (two components) situation, the strategy is the following:
* treat in the same way the two possibilities: (d dbar) (i=0) and
* (u ubar) (i=1) as the pair quark-antiquark necessary to form a
* pair of hadrons containing the input flavour par1 and par2; finally,
* select the one that produces the lightest pair of hadrons, compatible
* with the charge conservation constraint.
*/
pair<tcPDPtr,tcPDPtr> lightestHadronPair(tcPDPtr ptr1, tcPDPtr ptr2,
tcPDPtr ptr3 = PDPtr ()) const;
/**
* Returns the mass of the lightest pair of hadrons with the given particles
* @param ptr1 is the first constituent
* @param ptr2 is the second constituent
* @param ptr3 is the third constituent
*/
Energy massLightestHadronPair(tcPDPtr ptr1, tcPDPtr ptr2,
tcPDPtr ptr3 = PDPtr ()) const {
pair<tcPDPtr,tcPDPtr> pairData = lightestHadronPair(ptr1, ptr2, ptr3);
if ( ! pairData.first || ! pairData.second ) return ZERO;
return ( pairData.first->mass() + pairData.second->mass() );
}
/**
* Returns the lightest hadron formed by the given particles.
*
* Given the id of two (or three) constituents of a cluster, it returns
* the lightest hadron with proper flavour numbers.
* At the moment it does *nothing* in the case that also 'ptr3' present.
* @param ptr1 is the first constituent
* @param ptr2 is the second constituent
* @param ptr3 is the third constituent
*/
tcPDPtr lightestHadron(tcPDPtr ptr1, tcPDPtr ptr2,
tcPDPtr ptr3 = PDPtr ()) const;
/**
* Returns the hadrons below the constituent mass threshold formed by the given particles,
* together with their total weight
*
* Given the id of two (or three) constituents of a cluster, it returns
* the lightest hadron with proper flavour numbers.
* At the moment it does *nothing* in the case that also 'ptr3' present.
* @param threshold The theshold
* @param ptr1 is the first constituent
* @param ptr2 is the second constituent
* @param ptr3 is the third constituent
*/
vector<pair<tcPDPtr,double> >
hadronsBelowThreshold(Energy threshold,
tcPDPtr ptr1, tcPDPtr ptr2,
tcPDPtr ptr3 = PDPtr ()) const;
/**
* Return the nominal mass of the hadron returned by lightestHadron()
* @param ptr1 is the first constituent
* @param ptr2 is the second constituent
* @param ptr3 is the third constituent
*/
Energy massLightestHadron(tcPDPtr ptr1, tcPDPtr ptr2,
#ifndef NDEBUG
tcPDPtr ptr3 = PDPtr ()) const {
#else
tcPDPtr = PDPtr ()) const {
#endif
// The method assumes ptr3 == empty
assert(!(ptr3));
// find entry in the table
pair<long,long> ids(abs(ptr1->id()),abs(ptr2->id()));
HadronTable::const_iterator tit=_table.find(ids);
// throw exception if flavours wrong
if(tit==_table.end()||tit->second.empty())
throw Exception() << "HadronSelector::massLightestHadron "
<< "failed for particle" << ptr1->id() << " "
<< ptr2->id()
<< Exception::eventerror;
// return the mass
return tit->second.begin()->mass;
}
/**
* Returns the mass of the lightest pair of baryons.
* @param ptr1 is the first constituent
* @param ptr2 is the second constituent
*/
Energy massLightestBaryonPair(tcPDPtr ptr1, tcPDPtr ptr2) const;
/**
* Return the weights for the different quarks and diquarks
*/
//@{
/**
* 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 for spin-0 diquarks.
*/
double pwtDIquarkS0() const {
return _pwtDIquarkS0;
}
/**
* The diquark weight for spin-1 diquarks.
*/
double pwtDIquarkS1() const {
return _pwtDIquarkS1;
}
//@}
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
*
* The array _repwt is initialized using the interfaces to set different
* weights for different meson multiplets and the constructHadronTable()
* method called to complete the construction of the hadron tables.
*
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
//@}
protected:
/**
* Construct the table of hadron data
* This is the main method to initialize the hadron data (mainly the
* weights associated to each hadron, taking into account its spin,
* eventual isoscalar-octect mixing, singlet-decuplet factor). This is
* the method that one should update when new or updated hadron data is
* available.
*
* This class implements the construction of the basic table but can be
* overridden if needed in inheriting classes.
*
* The rationale for factors used for diquarks involving different quarks can
* be can be explained by taking a prototype example that in the exact SU(2) limit,
* in which:
* \f[m_u=m_d\f]
* \f[M_p=M_n=M_\Delta\f]
* and we will have equal numbers of u and d quarks produced.
* Suppose that we weight 1 the diquarks made of the same
* quark and 1/2 those made of different quarks, the fractions
* of u and d baryons (p, n, Delta) we get are the following:
* - \f$\Delta^{++}\f$: 1 possibility only u uu with weight 1
* - \f$\Delta^- \f$: 1 possibility only d dd with weight 1
* - \f$p,\Delta^+ \f$: 2 possibilities u ud with weight 1/2
* d uu with weight 1
* - \f$n,\Delta^0 \f$: 2 possibilities d ud with weight 1/2
* u dd with weight 1
* In the latter two cases, we have to take into account the
* fact that p and n have spin 1/2 whereas Delta+ and Delta0
* have spin 3/2 therefore from phase space we get a double weight
* for Delta+ and Delta0 relative to p and n respectively.
* Therefore the relative amount of these baryons that is
* produced is the following:
* # p = # n = ( 1/2 + 1 ) * 1/3 = 1/2
* # Delta++ = # Delta- = 1 = ( 1/2 + 1) * 2/3 # Delta+ = # Delta0
* which is correct, and therefore the weight 1/2 for the
* diquarks of different types of quarks is justified (at least
* in this limit of exact SU(2) ).
*/
virtual void constructHadronTable();
/**
* Access to the table of hadrons
*/
const HadronTable & table() const {
return _table;
}
/**
* Access to the list of partons
*/
const vector<PDPtr> & partons() const {
return _partons;
}
/**
* Access the parton weights
*/
double pwt(long pid) const {
map<long,double>::const_iterator it = _pwt.find(abs(pid));
assert( it != _pwt.end() );
return it->second;
}
/**
* Methods for the mixing of \f$I=0\f$ mesons
*/
//@{
/**
* Return the probability of mixing for Octet-Singlet isoscalar mixing,
* the probability of the
* \f$\frac1{\sqrt{2}}(|u\bar{u}\rangle + |d\bar{d}\rangle)\f$ component
* is returned.
* @param angleMix The mixing angle in degrees (not radians)
* @param order is 0 for no mixing, 1 for the first resonance of a pair,
* 2 for the second one.
* The mixing is defined so that for example with \f$\eta-\eta'\f$ mixing where
* the mixing angle is \f$\theta=-23^0$ with $\eta\f$ as the first particle
* and \f$\eta'\f$ the second one.
* The convention used is
* \f[\eta = \cos\theta|\eta_{\rm octet }\rangle
* -\sin\theta|\eta_{\rm singlet}\rangle\f]
* \f[\eta' = \sin\theta|\eta_{\rm octet }\rangle
* -\cos\theta|\eta_{\rm singlet}\rangle\f]
* with
* \f[|\eta_{\rm singlet}\rangle = \frac1{\sqrt{3}}
* \left[|u\bar{u}\rangle + |d\bar{d}\rangle + |s\bar{s}\rangle\right]\f]
* \f[|\eta_{\rm octet }\rangle = \frac1{\sqrt{6}}
* \left[|u\bar{u}\rangle + |d\bar{d}\rangle - 2|s\bar{s}\rangle\right]\f]
*/
double probabilityMixing(const double angleMix,
const int order) const {
static double convert=Constants::pi/180.0;
if (order == 1)
return sqr( cos( angleMix*convert + atan( sqrt(2.0) ) ) );
else if (order == 2)
return sqr( sin( angleMix*convert + atan( sqrt(2.0) ) ) );
else
return 1.;
}
/**
* Returns the weight of given mixing state.
* @param id The PDG code of the meson
*/
virtual double mixingStateWeight(long id) const;
//@}
/**
* Calculates a special weight specific to a given hadron.
* @param id The PDG code of the hadron
+ * @param part1 The PDG code of the first particle componant of the hadron
*/
- double specialWeight(long id) const;
+ double specialWeight(long id, int part1) const;
/**
* This method returns the proper sign ( > 0 hadron; < 0 anti-hadron )
* for the input PDG id idHad > 0, suppose to be made by the
* two constituent particle pointers: par1 and par2 (both with proper sign).
*/
int signHadron(tcPDPtr ptr1, tcPDPtr ptr2, tcPDPtr hadron) const;
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
HadronSelector & operator=(const HadronSelector &) = delete;
private:
/**
* The PDG codes of the constituent particles allowed
*/
vector<PDPtr> _partons;
/**
* The PDG codes of the hadrons which cannot be produced in the hadronization
*/
vector<PDPtr> _forbidden;
/**
* The weights for the different quarks and diquarks
*/
//@{
/**
* The probability of producting a down quark.
*/
double _pwtDquark;
/**
* The probability of producting an up quark.
*/
double _pwtUquark;
/**
* The probability of producting a strange quark.
*/
double _pwtSquark;
/**
* The probability of producting a charm quark.
*/
double _pwtCquark;
/**
* The probability of producting a bottom quark.
*/
double _pwtBquark;
/**
* The probability of producting a spin-0 diquark.
*/
double _pwtDIquarkS0;
/**
* The probability of producting a spin-1 diquark.
*/
double _pwtDIquarkS1;
/**
* Weights for quarks and diquarks.
*/
map<long,double> _pwt;
//@}
/**
* The mixing angles for the \f$I=0\f$ mesons containing light quarks
*/
//@{
/**
* The \f$\eta-\eta'\f$ mixing angle
*/
double _etamix;
/**
* The \f$\phi-\omega\f$ mixing angle
*/
double _phimix;
/**
* The \f$h_1'-h_1\f$ mixing angle
*/
double _h1mix;
/**
* The \f$f_0(1710)-f_0(1370)\f$ mixing angle
*/
double _f0mix;
/**
* The \f$f_1(1420)-f_1(1285)\f$ mixing angle
*/
double _f1mix;
/**
* The \f$f'_2-f_2\f$ mixing angle
*/
double _f2mix;
/**
* The \f$\eta_2(1870)-\eta_2(1645)\f$ mixing angle
*/
double _eta2mix;
/**
* The \f$\phi(???)-\omega(1650)\f$ mixing angle
*/
double _omhmix;
/**
* The \f$\phi_3-\omega_3\f$ mixing angle
*/
double _ph3mix;
/**
* The \f$\eta(1475)-\eta(1295)\f$ mixing angle
*/
double _eta2Smix;
/**
* The \f$\phi(1680)-\omega(1420)\f$ mixing angle
*/
double _phi2Smix;
//@}
/**
* The weights for the various meson multiplets to be used to supress the
* production of particular states
*/
//@{
/**
* The weights for the \f$\phantom{1}^1S_0\f$ multiplets
*/
vector<double> _weight1S0;
/**
* The weights for the \f$\phantom{1}^3S_1\f$ multiplets
*/
vector<double> _weight3S1;
/**
* The weights for the \f$\phantom{1}^1P_1\f$ multiplets
*/
vector<double> _weight1P1;
/**
* The weights for the \f$\phantom{1}^3P_0\f$ multiplets
*/
vector<double> _weight3P0;
/**
* The weights for the \f$\phantom{1}^3P_1\f$ multiplets
*/
vector<double> _weight3P1;
/**
* The weights for the \f$\phantom{1}^3P_2\f$ multiplets
*/
vector<double> _weight3P2;
/**
* The weights for the \f$\phantom{1}^1D_2\f$ multiplets
*/
vector<double> _weight1D2;
/**
* The weights for the \f$\phantom{1}^3D_1\f$ multiplets
*/
vector<double> _weight3D1;
/**
* The weights for the \f$\phantom{1}^3D_2\f$ multiplets
*/
vector<double> _weight3D2;
/**
* The weights for the \f$\phantom{1}^3D_3\f$ multiplets
*/
vector<double> _weight3D3;
//@}
/**
* The weights for the excited meson multiplets
*/
vector<vector<vector<double> > > _repwt;
/**
* Singlet and Decuplet weights
*/
//@{
/**
* The singlet weight
*/
double _sngWt;
/**
* The decuplet weight
*/
double _decWt;
//@}
/**
* The table of hadron data
*/
HadronTable _table;
/**
* Enums so arrays can be statically allocated
*/
//@{
/**
* Defines values for array sizes. L,J,N max values for excited mesons.
*/
enum MesonMultiplets { Lmax = 3, Jmax = 4, Nmax = 4};
//@}
/**
* Option for the construction of the tables
*/
unsigned int _topt;
/**
* Which particles to produce for debugging purposes
*/
unsigned int _trial;
/**
* @name A parameter used for determining when clusters are too light.
*
* This parameter is used for setting the lower threshold, \f$ t \f$ as
* \f[ t' = t(1 + r B^1_{\rm lim}) \f]
* where \f$ r \f$ is a random number [0,1].
*/
//@{
double _limBottom;
double _limCharm;
double _limExotic;
//@}
/**
* Option for the selection of hadrons below the pair threshold
*/
unsigned int belowThreshold_;
};
}
#endif /* HERWIG_HadronSelector_H */