diff --git a/Hadronization/ClusterFinder.cc b/Hadronization/ClusterFinder.cc
--- a/Hadronization/ClusterFinder.cc
+++ b/Hadronization/ClusterFinder.cc
@@ -1,710 +1,723 @@
// -*- 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 "Herwig/Utilities/EnumParticles.h"
#include "Herwig/Utilities/Kinematics.h"
#include "Cluster.h"
#include <ThePEG/Utilities/DescribeClass.h>
#include "ThePEG/Interface/Reference.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_ << _hadronSpectrum;
}
void ClusterFinder::persistentInput(PersistentIStream & is, int) {
is >> heavyDiquarks_ >> diQuarkSelection_ >> diQuarkOnShell_ >> _hadronSpectrum;
}
void ClusterFinder::Init() {
static ClassDocumentation<ClusterFinder> documentation
("This class is responsible of finding clusters.");
static Reference<ClusterFinder,HadronSpectrum> interfaceHadronSpectrum
("HadronSpectrum",
"Set the hadron spectrum for this parton splitter.",
&ClusterFinder::_hadronSpectrum, false, false, false, false, false);
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 two 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 SwitchOption interfaceDiQuarkSelectionLowestMassRescaled
(interfaceDiQuarkSelection,
"LowestMassRescaled",
"Combine the quark pair with the lowest value of (p1+p2)^2/(m1+m2)^2",
2);
static Switch<ClusterFinder,int> interfaceDiQuarkOnShell
("DiQuarkOnShell",
"Force the diquark produced in baryon-number violating clusters to be on-shell",
&ClusterFinder::diQuarkOnShell_, -1, false, false);
static SwitchOption interfaceDiQuarkOnShellYes
(interfaceDiQuarkOnShell,
"Yes",
"Force to be on-shell",
1);
static SwitchOption interfaceDiQuarkOnShellMixed
(interfaceDiQuarkOnShell,
"Mixed",
"Set Mass to constituent mass, but leave momentum as is.",
-1);
static SwitchOption interfaceDiQuarkOnShellNo
(interfaceDiQuarkOnShell,
"No",
"Leave off-shell",
0);
}
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());
}
}
+ClusterPtr ClusterFinder::handleBaryonicCluster(ClusterPtr clu) const {
+ // A modification of reduceToTwoComponents that is dedicated
+ // to make a quark-diquark or antiquark-antidiquark out of a
+ // baryonic or antibaryonic cluster
+ tParticleVector vec;
+ tPPtr other;
+ for(unsigned int i = 0; i<clu->numComponents(); i++) {
+ tPPtr part = clu->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 "
+ << clu->particle(0)->PDGName() << " "
+ << clu->particle(1)->PDGName() << " "
+ << clu->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) {
+ // TODO: Note this change is in principle not identical to
+ // Default, but IMO is reasonable
+ // 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) {
+ // TODO refactor this
+ 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;
+ }
+ // pick the lightest quark pair rescaled to the mass
+ else if(diQuarkSelection_==2) {
+ double m12 = (vec[0]->momentum()+vec[1]->momentum()).m()/(vec[0]->momentum().m()+vec[1]->momentum().m());
+ double m13 = (vec[0]->momentum()+vec[2]->momentum()).m()/(vec[0]->momentum().m()+vec[2]->momentum().m());
+ double m23 = (vec[1]->momentum()+vec[2]->momentum()).m()/(vec[1]->momentum().m()+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];
+
+ tcPDPtr dataDiquark = _hadronSpectrum->makeDiquark(temp1,temp2);
+
+ 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);
+ // use the same method as for cluster to determine the diquark position
+ diquark->setVertex(Cluster::calculateX(vec[0],vec[1]));
+ // Set momenta of diquarks
+ Lorentz5Momentum pDiq = vec[0]->momentum() + vec[1]->momentum();
+ pDiq.setMass(dataDiquark->constituentMass());
+ diquark->set5Momentum(Lorentz5Momentum(pDiq, dataDiquark->constituentMass()));
+ switch (diQuarkOnShell_)
+ {
+ case 0:
+ {
+ // No: Set the momentum of the diquarks to the original and mass to the Off-shell mass
+ diquark->set5Momentum(Lorentz5Momentum(pDiq, pDiq.m()));
+ break;
+ }
+ case -1:
+ {
+ // Mixed (default): Set the momentum of the diquarks to the original, but mass to constituentMass
+ // as above
+ break;
+ }
+ case 1:
+ {
+ // Yes: Rescale the momentum and mass of the diquarks to their diquark constituentMass
+ Lorentz5Momentum psum = pDiq+other->momentum();
+ psum.rescaleMass();
+ Boost boost = psum.boostVector();
+ Lorentz5Momentum pother = other->momentum();
+ Lorentz5Momentum pdiquark = pDiq;
+ 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);
+ } else {
+ throw Exception()
+ << "Not enough mass for baryonic cluster M = " << psum.mass()/GeV
+ << " m1 = " << other->dataPtr()->constituentMass()/GeV
+ << " m2 = " << dataDiquark ->constituentMass()/GeV
+ << " to be on shell! Pcom = " << pcm/GeV
+ << Exception::runerror;
+ }
+ break;
+ }
+ default:
+ assert(false);
+ }
+
+ // 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.
+ clu->addChild(nclus);
+ return nclus;
+}
ClusterPtr ClusterFinder::handleDiquarkCluster(ClusterPtr clu) const {
// A modification of reduceToTwoComponents that is dedicated
// to make two diquark pairs out of the diquark cluster
tParticleVector vec;
tPPtr other;
for(unsigned int i = 0; i<clu->numComponents(); i++) {
tPPtr part = clu->particle(i);
// Check if the parton is already a diquark
// if it isn't, add it to the 'to-be-diquarked' list
if(!DiquarkMatcher::Check(*(part->dataPtr()))) vec.push_back(part);
else other = part;
}
assert(clu->particle(0)->id()*clu->particle(1)->id()>0);
assert(clu->particle(2)->id()*clu->particle(3)->id()>0);
if(vec.size()<4) {
throw Exception() << "Could not make diquark pairs for a diquark cluster decay from "
<< clu->particle(0)->PDGName() << " "
<< clu->particle(1)->PDGName() << " "
<< clu->particle(2)->PDGName() << " "
<< clu->particle(3)->PDGName() << " "
<< " in ClusterFinder::handleDiquarkCluster()."
<< Exception::eventerror;
}
// order the vector so heaviest at the end
//std::sort(vec.begin(),vec.end(),PartOrdering);
// TODO: Figure out heavy quark stuff for tetraquark situation
// Special treatment of heavy quarks
// avoid doubly heavy diquarks
// if(heavyDiquarks_>=1 && 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 && abs(vec[2]->id())>3 && abs(vec[1]->id())<=3) {
// other = vec[2];
// vec.pop_back();
// }
// For tetraquarks there is no choice available, there are 2 quarks,
// and 2 antiquarks, and they have to be paired up together.
tcPDPtr tempQ1 = vec[0]->dataPtr();
tcPDPtr tempQ2 = vec[1]->dataPtr();
tcPDPtr tempQb1 = vec[2]->dataPtr();
tcPDPtr tempQb2 = vec[3]->dataPtr();
tcPDPtr dataDiquark = _hadronSpectrum->makeDiquark(tempQ1, tempQ2);
tcPDPtr dataDiquarkBar = _hadronSpectrum->makeDiquark(tempQb1, tempQb2);
if (!dataDiquark){
throw Exception() << "Could not make a diquark from"
<< tempQ1->PDGName() << " and "
<< tempQ2->PDGName()
<< " in ClusterFinder::handleDiquarkCluster()"
<< Exception::eventerror;
}
if (!dataDiquarkBar){
throw Exception() << "Could not make an anti-diquark from"
<< tempQb1->PDGName() << " and "
<< tempQb2->PDGName()
<< " in ClusterFinder::handleDiquarkCluster()"
<< Exception::eventerror;
}
// Create the new cluster (with two components) and assign it the same
// momentum and position of the original (with 4 components).
// Construct the diquarks
PPtr diquark = dataDiquark->produceParticle();
PPtr diquarkBar = dataDiquarkBar->produceParticle();
// update the new child cluster for the partons
vec[0]->addChild(diquark);
vec[1]->addChild(diquark);
vec[2]->addChild(diquarkBar);
vec[3]->addChild(diquarkBar);
// use the same method as for cluster to determine the diquark position
diquark->setVertex(Cluster::calculateX(vec[0],vec[1]));
diquarkBar->setVertex(Cluster::calculateX(vec[2],vec[3]));
// Set momenta of diquarks
Lorentz5Momentum pDiq = vec[0]->momentum() + vec[1]->momentum();
Lorentz5Momentum pAntiDiq = vec[2]->momentum() + vec[3]->momentum();
pDiq.setMass(dataDiquark->constituentMass());
pAntiDiq.setMass(dataDiquarkBar->constituentMass());
// Set the momentum of the diquarks
diquark->set5Momentum(Lorentz5Momentum(pDiq, dataDiquark->constituentMass()));
diquarkBar->set5Momentum(Lorentz5Momentum(pAntiDiq, dataDiquarkBar->constituentMass()));
switch (diQuarkOnShell_)
{
case 0:
{
// No: Set the momentum of the diquarks to the original and mass to the Off-shell mass
diquark->set5Momentum(Lorentz5Momentum(pDiq, pDiq.m()));
diquarkBar->set5Momentum(Lorentz5Momentum(pAntiDiq, pAntiDiq.m()));
break;
}
case -1:
{
// Mixed (default): Set the momentum of the diquarks to the original, but mass to constituentMass
// as above
break;
}
case 1:
{
// Yes: Rescale the momentum and mass of the diquarks to their diquark constituentMass
Lorentz5Momentum psum = pDiq + pAntiDiq;
psum.rescaleMass();
Boost boost = psum.boostVector();
Lorentz5Momentum pdiquarkBar = pAntiDiq;
Lorentz5Momentum pdiquark = pDiq;
pdiquarkBar.boost(-boost);
pdiquark.boost(-boost);
Energy pcm = Kinematics::pstarTwoBodyDecay(psum.mass(),
diquarkBar->dataPtr()->constituentMass(),
diquark->dataPtr()->constituentMass());
if(pcm>ZERO) {
double fact = pcm/pdiquarkBar.vect().mag();
pdiquarkBar *= fact;
pdiquark *= fact;
pdiquarkBar.setMass(dataDiquarkBar->constituentMass());
pdiquark .setMass(dataDiquark ->constituentMass());
pdiquarkBar.rescaleEnergy();
pdiquark .rescaleEnergy();
pdiquarkBar.boost(boost);
pdiquark .boost(boost);
diquarkBar->set5Momentum(pdiquarkBar);
diquark ->set5Momentum(pdiquark);
} else {
throw Exception()
<< "Not enough mass for cluster Diquark cluster M = " << psum.mass()/GeV
<< " m1 = " << diquarkBar->dataPtr()->constituentMass()/GeV
<< " m2 = " << diquark->dataPtr()->constituentMass()/GeV
<< " to be on shell! Pcom = " << pcm/GeV
<< Exception::runerror;
}
break;
}
default:
assert(false);
}
// make the new cluster
ClusterPtr nclus = new_ptr(Cluster(diquark, diquarkBar));
// Set the parent/child relationship
clu->addChild(nclus);
return nclus;
}
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() == 2
&& (
DiquarkMatcher::Check(*((*cluIter)->particle(0)->dataPtr()))
|| DiquarkMatcher::Check(*((*cluIter)->particle(1)->dataPtr())))){
if (fabs((*cluIter)->momentum().m()/GeV-((*cluIter)->particle(0)->momentum()+(*cluIter)->particle(1)->momentum()).m()/GeV)>1e-5){
std::cout << "momentum().m() "<< std::setprecision(18) << (*cluIter)->momentum().m()/GeV << std::endl;
std::cout << "(p1+p2).m() "<< std::setprecision(18) << ((*cluIter)->particle(0)->momentum()+(*cluIter)->particle(1)->momentum()).m()/GeV << std::endl;
}
}
if ( (*cluIter)->numComponents() == 4 && (*cluIter)->isAvailable() ){
- ClusterPtr nclus = handleDiquarkCluster(*cluIter);
- assert(nclus->numComponents()==2);
- redefinedClusters.push_back(nclus);
+ // Handle diquark clusters
+ ClusterPtr DiquarkClusterReduced = handleDiquarkCluster(*cluIter);
+ assert(DiquarkClusterReduced->numComponents()==2);
+ redefinedClusters.push_back(DiquarkClusterReduced);
continue;
}
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) {
- // TODO refactor this
- 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;
- }
- // pick the lightest quark pair rescaled to the mass
- else if(diQuarkSelection_==2) {
- double m12 = (vec[0]->momentum()+vec[1]->momentum()).m()/(vec[0]->momentum().m()+vec[1]->momentum().m());
- double m13 = (vec[0]->momentum()+vec[2]->momentum()).m()/(vec[0]->momentum().m()+vec[2]->momentum().m());
- double m23 = (vec[1]->momentum()+vec[2]->momentum()).m()/(vec[1]->momentum().m()+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];
-
- tcPDPtr dataDiquark = _hadronSpectrum->makeDiquark(temp1,temp2);
-
- 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);
- // use the same method as for cluster to determine the diquark position
- diquark->setVertex(Cluster::calculateX(vec[0],vec[1]));
- // Set momenta of diquarks
- Lorentz5Momentum pDiq = vec[0]->momentum() + vec[1]->momentum();
- pDiq.setMass(dataDiquark->constituentMass());
- diquark->set5Momentum(Lorentz5Momentum(pDiq, dataDiquark->constituentMass()));
- switch (diQuarkOnShell_)
- {
- case 0:
- {
- // No: Set the momentum of the diquarks to the original and mass to the Off-shell mass
- diquark->set5Momentum(Lorentz5Momentum(pDiq, pDiq.m()));
- break;
- }
- case -1:
- {
- // Mixed (default): Set the momentum of the diquarks to the original, but mass to constituentMass
- // as above
- break;
- }
- case 1:
- {
- // Yes: Rescale the momentum and mass of the diquarks to their diquark constituentMass
- Lorentz5Momentum psum = pDiq+other->momentum();
- psum.rescaleMass();
- Boost boost = psum.boostVector();
- Lorentz5Momentum pother = other->momentum();
- Lorentz5Momentum pdiquark = pDiq;
- 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);
- } else {
- throw Exception()
- << "Not enough mass for baryonic cluster M = " << psum.mass()/GeV
- << " m1 = " << other->dataPtr()->constituentMass()/GeV
- << " m2 = " << dataDiquark ->constituentMass()/GeV
- << " to be on shell! Pcom = " << pcm/GeV
- << Exception::runerror;
- }
- break;
- }
- default:
- assert(false);
- }
-
- // 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);
+ // Should now only have baryonic clusters
+ // Handle baryonic clusters
+ ClusterPtr BaryonicClusterReduced = handleBaryonicCluster(*cluIter);
+ assert(BaryonicClusterReduced->numComponents()==2);
+ redefinedClusters.push_back(BaryonicClusterReduced);
}
// 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/ClusterFinder.h b/Hadronization/ClusterFinder.h
--- a/Hadronization/ClusterFinder.h
+++ b/Hadronization/ClusterFinder.h
@@ -1,176 +1,186 @@
// -*- C++ -*-
//
// ClusterFinder.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_ClusterFinder_H
#define HERWIG_ClusterFinder_H
#include <ThePEG/Interface/Interfaced.h>
#include "CluHadConfig.h"
#include "ClusterFinder.fh"
#include "HadronSpectrum.h"
namespace Herwig {
using namespace ThePEG;
/*! \ingroup Hadronization
* \class ClusterFinder
* \brief This class forms clusters from the partons produced in the Shower.
* \author Philip Stephens
* \author Alberto Ribon
*
* This class scans through the particles in the event and produces a
* collection of clusters, defined as a colour-singlet combinations of
* colour-connected particles. There are no assumptions about the type
* (i.e. quark or diquark) or number of the component particles of the
* cluster; however, most of the time clusters are formed by quark-antiquark
* pairs. In special situations, such as baryon-violating processes in
* R-nonconserved Susy, three quarks (or three antiquarks) could form a
* cluster. Because at the moment we don't know how to handle 3-component
* clusters (i.e. how to fission heavy ones, or how to decay clusters), we
* provide also a separate method, reduceToTwoComponents, which
* does the job of redefining these 3-component clusters as "normal"
* 2-component ones, simply by randomly considering two (anti-) quarks as a
* (anti-) diquark. Notice that if in the future the method
* reduceToTwoComponents is modified or even eliminated, the
* main method for finding clusters, formClusters, will not need
* any change.
*
* @see \ref ClusterFinderInterfaces "The interfaces"
* defined for ClusterFinder.
*/
class ClusterFinder: public Interfaced {
public :
/** @name Standard constructors and destructors. */
//@{
/**
* Default constructor.
*/
ClusterFinder() : heavyDiquarks_(2), diQuarkSelection_(1), diQuarkOnShell_(-1)
{}
//@}
public:
/**
* This routine forms the clusters of the event.
*
* Form clusters starting from the list of partons given.
* It also checks if the cluster is a beam cluster, that is if
* at least one of its components is a beam remnant.
*/
ClusterVector formClusters(const PVector & partons);
/**
* Reduces three component clusters into two components.
*
* For the eventual clusters that have three components
* (quark, quark, quark) or (antiquark, antiquark, antiquark),
* it redefines them as "normal" clusters with two components:
* (quark,diquark) or (antiquark,antidiquark), by a random drawing.
* This could be eliminated or changed in the future.
*/
void reduceToTwoComponents(ClusterVector&);
/**
* handling Diquark clusters:
* For the eventual clusters that have four components
* (quark, quark, antiquark, antiquark)
* it redefines them as "normal" clusters with two components:
* (antidiquark,diquark) by a random drawing.
*/
ClusterPtr handleDiquarkCluster(ClusterPtr clu) const;
+
+ /**
+ * handling Baryonic clusters:
+ * For the eventual clusters that have three components
+ * (quark, quark, quark) or (antiquark, antiquark, antiquark)
+ * it redefines them as "normal" clusters with two components:
+ * (quark, diquark) or (antiquark, antidiquark) by a prescription.
+ */
+ ClusterPtr handleBaryonicCluster(ClusterPtr clu) const;
+
/**
* Return the hadron spectrum
*/
Ptr<HadronSpectrum>::tptr spectrum() const {
return _hadronSpectrum;
}
/**
* access for diquark treatment
*/
int diquarkOnShell() {
return diQuarkOnShell_;
}
public:
/**
* Standard Init function used to initialize the interfaces.
*/
static void Init();
/** @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);
//@}
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
private:
/**
* Private and non-existent assignment operator.
*/
ClusterFinder & operator=(const ClusterFinder &) = delete;
private:
/**
* Treatment of diquarks contain heavy quarks in baryon-number violating clusters
*/
unsigned int heavyDiquarks_;
/**
* Option for the selection of which quarks to make into a diquark
*/
unsigned int diQuarkSelection_;
/**
* Force diquarks to be on-shell
*/
int diQuarkOnShell_;
/**
* The hadron spectrum to consider
*/
Ptr<HadronSpectrum>::ptr _hadronSpectrum;
};
}
#endif /* HERWIG_ClusterFinder_H */