Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F8310387
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
235 KB
Subscribers
None
View Options
diff --git a/Hadronization/ClusterFinder.cc b/Hadronization/ClusterFinder.cc
--- a/Hadronization/ClusterFinder.cc
+++ b/Hadronization/ClusterFinder.cc
@@ -1,710 +1,742 @@
// -*- 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);
- continue;
+ // Handle diquark clusters
+ ClusterPtr DiquarkClusterReduced = handleDiquarkCluster(*cluIter);
+ if (DiquarkClusterReduced) {
+ assert(DiquarkClusterReduced->numComponents()==2);
+ redefinedClusters.push_back(DiquarkClusterReduced);
+ } else {
+ throw Exception()
+ << "Could not reduce Diquark Cluster M = " << DiquarkClusterReduced->momentum().m()/GeV
+ << "\n m1 = " << DiquarkClusterReduced->particle(0)->dataPtr()->constituentMass()/GeV
+ << "\n m2 = " << DiquarkClusterReduced->particle(1)->dataPtr()->constituentMass()/GeV
+ << "\n m3 = " << DiquarkClusterReduced->particle(2)->dataPtr()->constituentMass()/GeV
+ << "\n m4 = " << DiquarkClusterReduced->particle(3)->dataPtr()->constituentMass()/GeV
+ << Exception::runerror;
+ }
+ 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);
+ // Should now only have baryonic clusters
+ // Handle baryonic clusters
+ ClusterPtr BaryonicClusterReduced = handleBaryonicCluster(*cluIter);
+ if (BaryonicClusterReduced) {
+ assert(BaryonicClusterReduced->numComponents()==2);
+ redefinedClusters.push_back(BaryonicClusterReduced);
+ } else {
+ throw Exception()
+ << "Could not reduce baryonic cluster M = " << BaryonicClusterReduced->momentum().m()/GeV
+ << "\n m1 = " << BaryonicClusterReduced->particle(0)->dataPtr()->constituentMass()/GeV
+ << "\n m2 = " << BaryonicClusterReduced->particle(1)->dataPtr()->constituentMass()/GeV
+ << "\n m3 = " << BaryonicClusterReduced->particle(2)->dataPtr()->constituentMass()/GeV
+ << Exception::runerror;
}
-
- // 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/ClusterFinder.h b/Hadronization/ClusterFinder.h
--- a/Hadronization/ClusterFinder.h
+++ b/Hadronization/ClusterFinder.h
@@ -1,169 +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 */
diff --git a/Hadronization/ColourReconnector.cc b/Hadronization/ColourReconnector.cc
--- a/Hadronization/ColourReconnector.cc
+++ b/Hadronization/ColourReconnector.cc
@@ -1,4124 +1,4220 @@
// -*- C++ -*-
//
// ColourReconnector.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 ColourReconnector class.
//
#include "ColourReconnector.h"
#include "Cluster.h"
#include <ThePEG/Utilities/DescribeClass.h>
#include <ThePEG/Repository/UseRandom.h>
#include <ThePEG/PDT/StandardMatchers.h>
#include <ThePEG/Persistency/PersistentOStream.h>
#include <ThePEG/Persistency/PersistentIStream.h>
#include <ThePEG/Interface/Switch.h>
#include <ThePEG/Interface/Reference.h>
#include <ThePEG/Interface/Parameter.h>
#include "Herwig/MatrixElement/Matchbox/CVolver/ColourFlows.h"
#include "Herwig/Utilities/Maths.h"
#include "Herwig/Utilities/expm-1.h"
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <cassert>
using namespace Herwig;
using CluVecIt = ColourReconnector::CluVecIt;
using Constants::pi;
using Constants::twopi;
DescribeClass<ColourReconnector,Interfaced>
describeColourReconnector("Herwig::ColourReconnector","Herwig.so");
IBPtr ColourReconnector::clone() const {
return new_ptr(*this);
}
IBPtr ColourReconnector::fullclone() const {
return new_ptr(*this);
}
void ColourReconnector::rearrange(ClusterVector & clusters) {
if (_clreco == 0) return;
// need at least two clusters
if (clusters.size() < 2) return;
if (0) {
CVolver::ColourFlow c1({0,1,2});
CVolver::ColourFlow c2({1,0,2});
for (const auto & p : c1.permutation())
{
std::cout << "c1 = " << p<< std::endl;
}
for (const auto & p : c2.permutation())
{
std::cout << "c2 = " << p<< std::endl;
}
std::cout << "<c2|c1> = " << c2.scalarProduct(c1) << std::endl;
std::cout << "<c1|c2> = " << c1.scalarProduct(c2) << std::endl;
std::set<CVolver::ColourFlow> setCF = CVolver::ColourFlow::allFlows(4);
for (const auto & perm : setCF) {
for (const auto & p : perm.permutation())
std::cout << p << " ";
std::cout << std::endl;
}
for (const auto & perm1 : setCF) {
for (const auto & perm2 : setCF) {
std::cout << perm1.scalarProduct(perm2) << " ";
}
std::cout << std::endl;
}
std::cout << "Correct matrix?:" << std::endl;
for (const auto & perm1 : setCF) {
for (auto perm2 : setCF) {
std::cout << perm1.scalarProduct(perm2.conjugate()) << " ";
}
std::cout << std::endl;
}
CVolver::ColourFlow c4({1,2,0});
std::cout << "<c4|c4> = " << c4.scalarProduct(c4) << std::endl;
// std::cout << "s(c1,c2) = Nc^" << c1.scalarProduct(c2)<< std::endl;
}
// std::cout << "size before = "<< clusters.size() << std::endl;
for (unsigned int i = 0; i < _crIterations; i++)
{
// do the colour reconnection
switch (_algorithm) {
case 0:
_doRecoPlain(clusters);
break;
case 1:
// TODO This algorithm has no dynamic CR
_doRecoStatistical(clusters);
break;
case 2:
_doRecoBaryonic(clusters);
break;
case 3:
// TODO This algorithm has no dynamic CR
_doRecoBaryonicMesonic(clusters);
break;
case 4:
_doRecoBaryonicDiquarkCluster(clusters);
break;
default:
assert(false);
}
}
// std::cout << "size after = "<< clusters.size() << std::endl;
}
namespace {
inline int hasDiquark(const ClusterPtr & cit) {
int res = 0;
for (unsigned int i = 0; i<(cit)->numComponents(); i++) {
if (DiquarkMatcher::Check(*((cit)->particle(i)->dataPtr()))) {
res++;
}
}
return res;
}
double calculateRapidityRF(const Lorentz5Momentum & q1,
const Lorentz5Momentum & p2) {
//calculate rapidity wrt the direction of q1
//angle between the particles in the RF of cluster of q1
// calculate the z component of p2 w.r.t the direction of q1
if(q1.rho2()==ZERO) return 0.;
const Energy pz = p2.vect() * q1.vect().unit();
if ( pz == ZERO ) return 0.;
// Transverse momentum of p2 w.r.t the direction of q1
const Energy2 pt2 = p2.vect().mag2() - sqr(pz);
// Transverse mass pf p2 w.r.t to the direction of q1
const Energy mtrans = sqrt(p2.mass()*p2.mass() + pt2);
// Correct formula
const double y2 = log((p2.t() + abs(pz))/mtrans);
return ( pz < ZERO ) ? -y2 : y2;
}
}
Energy2 ColourReconnector::_clusterMassSum(const PVector & q,
const PVector & aq) const {
const size_t nclusters = q.size();
assert (aq.size() == nclusters);
Energy2 sum = ZERO;
for (size_t i = 0; i < nclusters; i++)
sum += ( q[i]->momentum() + aq[i]->momentum() ).m2();
return sum;
}
double ColourReconnector::_displacement(tcPPtr p, tcPPtr q) const {
double deltaRap = (p->rapidity() - q->rapidity());
double deltaPhi = fabs(p->momentum().phi() - q->momentum().phi());
// keep deltaPhi's below Pi due to periodicity
if (deltaPhi > M_PI) deltaPhi-=M_PI;
return sqrt(deltaRap * deltaRap + deltaPhi * deltaPhi);
}
/**
* Computes circular Mean of three angles alpha, beta, gamma
* */
static double circularMean(double alpha, double beta, double gamma) {
double xMean=(cos(alpha)+cos(beta)+cos(gamma))/3.0;
double yMean=(sin(alpha)+sin(beta)+sin(gamma))/3.0;
// to make the function fail-save
if (xMean==0 && yMean==0) return M_PI;
return atan2(yMean,xMean);
}
namespace {
// ColourFlow for 3 colour flows for baryon state
// ordered permutations by one transposition
// sign of permutations # of difference
// |1> = |123> + 0
// |2> = |132> - 1
// |3> = |213> - 1
// |4> = |231> + 2
// |5> = |312> + 2
// |6> = |321> - 1
int signPermutationState(int i);
int signPermutationState(int i)
{
if (i==0 || i==3 || i==4) return 1;
else if (i==1 || i==2 || i==5) return -1;
else assert(false);
}
// ColourFlow scalar product matrix for 3
// colour flows in the following basis
// |1> = |123> + 0
// |2> = |132> - 1
// |3> = |213> - 1
// |4> = |231> + 2
// |5> = |312> + 2
// |6> = |321> - 1
unsigned int scalarProducts(int i, int j);
unsigned int scalarProducts(int i, int j)
{
// Verified for i,j < 3! and 3 colour flows
if (i>j) return scalarProducts(j,i);
// TODO need to restore Nc dependence
unsigned int Nc=3;
if (i==j) return Nc*Nc*Nc;
switch(i)
{
case 0:
{
if (j==1 || j==2 || j==5) return Nc*Nc;
else if (j==3 || j==4) return Nc;
else return Nc*Nc*Nc;
break;
}
case 1:
{
if (j==0 || j==3 || j==4) return Nc*Nc;
else if (j==2 || j==5) return Nc;
else return Nc*Nc*Nc;
break;
}
case 2:
{
if (j==0 || j==3 || j==4) return Nc*Nc;
else if (j==1 || j==5) return Nc;
else return Nc*Nc*Nc;
break;
}
case 3:
{
if (j==1 || j==2 || j==5) return Nc*Nc;
else if (j==0 || j==4) return Nc;
else return Nc*Nc*Nc;
break;
}
case 4:
{
if (j==1 || j==2 || j==5) return Nc*Nc;
else if (j==0 || j==3) return Nc;
else return Nc*Nc*Nc;
break;
}
case 5:
{
if (j==0 || j==3 || j==4) return Nc*Nc;
else if (j==1 || j==2) return Nc;
else return Nc*Nc*Nc;
break;
}
default:
assert(false);
}
return Nc;
}
}
// TODO add nonTrivialInitialState
std::unordered_map<int,double> ColourReconnector::_reconnectionAmplitudesCF2 (const ClusterPtr & c1, const ClusterPtr & c2, const int) const{
// Verified according to convention of analytics/matrices2_BCR.nb and not
// according to https://arxiv.org/abs/1808.06770
// The same convention as in https://arxiv.org/abs/1808.06770 can be obtained by
// the mapping 1->1;2->4;3->2;4->3 (here->paper)
Lorentz5Momentum p1 = c1->colParticle()->momentum();
Lorentz5Momentum p2 = c1->antiColParticle()->momentum();
Lorentz5Momentum p3 = c2->colParticle()->momentum();
Lorentz5Momentum p4 = c2->antiColParticle()->momentum();
Energy mLightestConstMass = 1000*GeV;
for (const auto & id : _hadronSpectrum->lightHadronizingQuarks() ) {
if (getParticleData(id)->constituentMass()<mLightestConstMass)
mLightestConstMass = getParticleData(id)->constituentMass();
}
// TODO REFACTOR THIS TO ALLOW FOR DIFFERENT SCALE CHOICES
double M12 = (p1 + p2).m2()/sqr(2*mLightestConstMass*_dynamicCRscale);
double M24 = (p2 + p4).m2()/sqr(2*mLightestConstMass*_dynamicCRscale);
double M13 = (p1 + p3).m2()/sqr(2*mLightestConstMass*_dynamicCRscale);
double M23 = (p2 + p3).m2()/sqr(2*mLightestConstMass*_dynamicCRscale);
double M14 = (p1 + p4).m2()/sqr(2*mLightestConstMass*_dynamicCRscale);
double M34 = (p3 + p4).m2()/sqr(2*mLightestConstMass*_dynamicCRscale);
if (
M12 < 1.0 ||
M24 < 1.0 ||
M13 < 1.0 ||
M23 < 1.0 ||
M14 < 1.0 ||
M34 < 1.0
) {
throw Exception()
// TODO REFACTOR THIS TO ALLOW FOR DIFFERENT SCALE CHOICES
<< "DynamicCR scale "<< _dynamicCRscale << " too high in ColourReconnector"
<< " must be less than 1"
<< " invariant masses:\n"
<< "M12 = " << sqrt(M12*sqr(2*mLightestConstMass*_dynamicCRscale))/GeV << "\n"
<< "M24 = " << sqrt(M24*sqr(2*mLightestConstMass*_dynamicCRscale))/GeV << "\n"
<< "M13 = " << sqrt(M13*sqr(2*mLightestConstMass*_dynamicCRscale))/GeV << "\n"
<< "M23 = " << sqrt(M23*sqr(2*mLightestConstMass*_dynamicCRscale))/GeV << "\n"
<< "M14 = " << sqrt(M14*sqr(2*mLightestConstMass*_dynamicCRscale))/GeV << "\n"
<< "M34 = " << sqrt(M34*sqr(2*mLightestConstMass*_dynamicCRscale))/GeV << "\n"
<< Exception::eventerror;
}
double alphaQCD=_dynamicCRalphaS;
// TODO missing factor of two due to sum_i!=j
// can be fixed by twopi->pi
double logSqrOmega12=alphaQCD*pow(log(M12),2)/(2.0*twopi);
double logSqrOmega24=alphaQCD*pow(log(M24),2)/(2.0*twopi);
double logSqrOmega13=alphaQCD*pow(log(M13),2)/(2.0*twopi);
double logSqrOmega23=alphaQCD*pow(log(M23),2)/(2.0*twopi);
double logSqrOmega14=alphaQCD*pow(log(M14),2)/(2.0*twopi);
double logSqrOmega34=alphaQCD*pow(log(M34),2)/(2.0*twopi);
// TODO need to restore Nc dependence
double Nc=3.0;
double U11,U21; // relevant matrix elements
switch (_dynamicCR)
{
case 1:
{
double a = (logSqrOmega34 + logSqrOmega12)/2.0;
double b = (logSqrOmega14 + logSqrOmega23)/2.0;
double c = (logSqrOmega13 + logSqrOmega24)/2.0;
double sqrtDelta=sqrt(Nc*Nc*a*a-4*c*(a+b)-(2*Nc*Nc-4)*a*b+Nc*Nc*b*b+4*c*c);
U11=sqrtDelta/tanh(sqrtDelta/2.0)+3*(b-a);
U21=2*(c-b);
break;
}
case 2:
{
// not Exponentiated soft anomalous dimension
// i.e. eq (5.2) Omega matrix in https://arxiv.org/pdf/1808.06770
U11=-Nc*0.5*(logSqrOmega34+logSqrOmega12);
U21= 0.5*(logSqrOmega13+logSqrOmega24-(logSqrOmega14+logSqrOmega23));
break;
}
default:
assert(false);
}
// TODO need to restore Nc dependence
double TransAmpNoCR = Nc*(Nc * U11 + U21);
double TransAmpMesonicCR = Nc*(U11 + Nc * U21);
std::unordered_map<int,double> amplitudes;
// No Colour Reconnection amplitude
amplitudes[123] = TransAmpNoCR;
// Mesonic Colour Reconnection amplitude
amplitudes[213] = TransAmpMesonicCR;
return amplitudes;
}
std::tuple<double,double,double> ColourReconnector::_dynamicRecoProbabilitiesCF2(const ClusterPtr & c1, const ClusterPtr & c2, bool diquarkCR) const{
std::unordered_map<int,double> amplitudes = _reconnectionAmplitudesCF2 (c1, c2);
double pNoCR;
double pMesonicCR;
double pDiquarkCR = 0.0;
double TransAmpNoCR = sqr(amplitudes[123]);
double TransAmpMesonicCR = _becomesColour8Cluster(c1,c2) ? 0.0:sqr(amplitudes[213]);
double sum = TransAmpNoCR + TransAmpMesonicCR;
double PhaseSpace=0.0;
if (diquarkCR && _canMakeDiquarkCluster(c1,c2,PhaseSpace) && !hasDiquark(c1) && !hasDiquark(c2)) {
// TODO need to restore Nc dependence
double ND = 2.0/sqrt(3.0); // sqrt(2*(Nc-1.0)/Nc);
double TransAmpDiquarkCR = sqr((amplitudes[123]-amplitudes[213])/ND);
sum += _phaseSpaceDiquarkFission ? PhaseSpace*TransAmpDiquarkCR:TransAmpDiquarkCR;
pDiquarkCR = TransAmpDiquarkCR/sum;
if (_phaseSpaceDiquarkFission) pDiquarkCR*=PhaseSpace;
assert( pDiquarkCR<=1.0 && pDiquarkCR>=0.0);
}
pNoCR = TransAmpNoCR/sum;
pMesonicCR = TransAmpMesonicCR/sum;
assert( pNoCR<=1.0 && pNoCR>=0.0);
assert( pMesonicCR<=1.0 && pMesonicCR>=0.0);
if (_debug)
{
Lorentz5Momentum p1col = (c1)->colParticle()->momentum();
Lorentz5Momentum p1acol = (c1)->antiColParticle()->momentum();
Lorentz5Momentum p2col = (c2)->colParticle()->momentum();
Lorentz5Momentum p2acol = (c2)->antiColParticle()->momentum();
const Boost boostv1=-c1->momentum().boostVector();
const Boost boostv2=-c2->momentum().boostVector();
p1col .boost(boostv1);
p1acol.boost(boostv1);
p2col .boost(boostv1);
p2acol.boost(boostv1);
double rap1c= calculateRapidityRF(p1acol,p2col);
double rap1a= calculateRapidityRF(p1acol,p2acol);
double pT1c = p2col.vect() .perp(p1acol.vect())/GeV;
double pT1a = p2acol.vect().perp(p1acol.vect())/GeV;
p1col .boost(-boostv1);
p1acol.boost(-boostv1);
p2col .boost(-boostv1);
p2acol.boost(-boostv1);
p1col .boost(boostv2);
p1acol.boost(boostv2);
p2col .boost(boostv2);
p2acol.boost(boostv2);
double rap2c= calculateRapidityRF(p2acol,p1col);
double rap2a= calculateRapidityRF(p2acol,p1acol);
double pT2c = p1col.vect() .perp(p1acol.vect())/GeV;
double pT2a = p1acol.vect().perp(p1acol.vect())/GeV;
// calculate the rapidity of the other constituents of the clusters
// w.r.t axis of p1anticol.vect.unit
// const double rapq = calculateRapidityRF(p1anticol,p2col);
ofstream out("WriteOut/kinematicRecoProbability.dat", std::ios::app);
out << pNoCR << "\t" << pMesonicCR << "\t" << pDiquarkCR << "\t";
out << rap1c << "\t" << rap1a << "\t" << pT1c << "\t" << pT1a << "\t";
out << rap2c << "\t" << rap2a << "\t" << pT2c << "\t" << pT2a << "\n";
out.close();
}
return {pNoCR, pMesonicCR, pDiquarkCR};
}
int ColourReconnector::_stateToPermutation(const int i) const {
switch (i)
{
case 0:
// Baryonic CR
return 0;
// Mesonic CR's
case 1:
return 123;
case 2:
return 132;
case 3:
return 213;
case 4:
return 231;
case 5:
return 312;
case 6:
return 321;
default:
assert(false);
}
}
Selector<int> ColourReconnector::_selector(const ClusterVector & clusters, bool diquarkCR) const {
switch (clusters.size())
{
case 2:
return getProbabilities2CF(clusters[0], clusters[1], diquarkCR);
break;
case 3:
{
if (_dynamicCR) {
return _selectorCF3(clusters, diquarkCR);
}
else
{
Selector<int> CRoptions;
- double sum=_preco+_precoBaryonic;
+ double sum = 0;
const int NpossibilitiesMCR = 6;
int state;
for (int i = 0; i < NpossibilitiesMCR; i++)
{
state = _stateToPermutation(i+1);
if (_isColour8Forbidden(state,clusters))
continue;
CRoptions.insert(_preco, state);
+ sum += _preco;
}
- CRoptions.insert(_precoBaryonic, 0);
+ if (_canMakeBaryonicCluster(clusters[0], clusters[1], clusters[2])) {
+ CRoptions.insert(_precoBaryonic, 0);
+ sum += _precoBaryonic;
+ }
if (diquarkCR) {
const int N=3;
bool first=false;
// Here i is the index of the quark and j of the antiquark
// which will be connected to each other (other partons will be
// forming the diquark cluster if kinematically viable)
double PhaseSpace;
for (int i = 1; i <= N; i++) {
for (int j = 1; j <= N; j++) {
state = -(10*i+j);
if (_isColour8Forbidden(state,clusters))
continue;
if (_canMakeDiquarkCluster(
clusters[i%3]->colParticle(), clusters[(i+1)%3]->colParticle(),
clusters[j%3]->antiColParticle(), clusters[(j+1)%3]->antiColParticle(),PhaseSpace)){
CRoptions.insert(_precoDiquark*PhaseSpace, state);
if (first){
sum+=_precoDiquark*PhaseSpace;
first=false;
}
}
}
}
}
// no CR probability
CRoptions.insert((1-sum) > 0 ? (1-sum):0, _stateToPermutation(1));
return CRoptions;
}
}
default:
assert(false);
}
}
/*
namespace {
int orientation(const ClusterPtr & c1,const ClusterPtr & c2){
Lorentz5Momentum cl1 = c1->momentum();
const Boost boostv(-cl1.boostVector());
// boost constituents of cl into RF of cl
Lorentz5Momentum p1anticol = c1->antiColParticle()->momentum();
p1anticol.boost(boostv);
// boost constituents of cit into RF of cl
Lorentz5Momentum p2col = c2->colParticle()->momentum();
Lorentz5Momentum p2anticol = c2->antiColParticle()->momentum();
p2col.boost(boostv);
p2anticol.boost(boostv);
// calculate the rapidity of the other constituents of the clusters
// w.r.t axis of p1anticol.vect.unit
const double rapq = calculateRapidityRF(p1anticol,p2col);
const double rapqbar = calculateRapidityRF(p1anticol,p2anticol);
if (rapq>0.0 && rapqbar<0.0 ) {
// Mesonic Type
return 1;
} else if (rapq<0.0 && rapqbar>0.0 ) {
// diquark Type
return -1;
}
else {
return 0;
}
}
}
*/
bool ColourReconnector::_isColour8Forbidden(int state, const ClusterVector & clusters) const{
if (state>0){
int c1 = ((state/100) % 4) -1;
int c2 = ((state-(c1+1)*100)/10 % 4) -1;
int c3 = ((state-(c1+1)*100-(c2+1)*10) % 4) -1;
assert(c1>=0 && c1<3);
assert(c2>=0 && c2<3);
assert(c3>=0 && c3<3);
assert(c1!=c2 && c1!=c3 && c2 !=c3);
assert(state==(c1+1)*100+(c2+1)*10+(c3+1));
if (c1!=0 && _isColour8(clusters[0]->colParticle(),clusters[c1]->antiColParticle()))
return true;
if (c2!=1 && _isColour8(clusters[1]->colParticle(),clusters[c2]->antiColParticle()))
return true;
if (c3!=2 && _isColour8(clusters[2]->colParticle(),clusters[c3]->antiColParticle()))
return true;
// int config1 = (0==c1) ? 1:orientation(clusters[0],clusters[c1]);
// int config2 = (1==c2) ? 1:orientation(clusters[1],clusters[c2]);
// int config3 = (2==c3) ? 1:orientation(clusters[2],clusters[c3]);
// if (!(config1>0 && config2>0 && config3>0))
// return true;
// if (!(config1<0 || config2<0 || config3<0))
// return true;
// if (!(config1>=0 && config2>=0 && config3>=0))
// return true;
}
else if (state<0){
int i = -state/10 -1;
int j = -state - 10*(i+1) - 1;
assert(-state==((i+1)*10+(j+1)));
if (_isColour8(clusters[i]->colParticle(),clusters[j]->antiColParticle()))
return true;
// TODO NOT GOOD orientation
// int config=orientation(clusters[i],clusters[j]);
// if (!(config>0))
// return true;
}
return false;
}
Selector<int> ColourReconnector::_selectorCF3(const ClusterVector & clusters, bool diquarkCR) const {
std::unordered_map<int, double> amplitudes = _reconnectionAmplitudesSGE(clusters);
double sum2 = 0;
double sumBaryon = 0;
double NB=2.0/sqrt(3.0);
Selector<int> CRoptions;
double amp2;
const int NpossibilitiesMCR = 6;
int state;
for (int i = 0; i < NpossibilitiesMCR; i++)
{
state=_stateToPermutation(i+1);
sumBaryon += amplitudes[state]*signPermutationState(i)/NB;
if (_isColour8Forbidden(state,clusters))
continue;
amp2 = sqr(amplitudes[state]);
sum2 += amp2;
CRoptions.insert(amp2, state);
}
sum2+=sumBaryon*sumBaryon;
if (diquarkCR) {
// TODO need to restore Nc dependence
double ND=2.0/sqrt(3.0); // sqrt(2*(Nc-1.0)/Nc)
const int N=3;
// TODO need to restore Nc dependence
double NDiqFACT=1.0; // sqrt(2*(Nc-1.0)/Nc)
double PhaseSpace;
int perm1;
int perm2;
// Here i is the index of the quark and j of the antiquark
// which will be connected to each other (other partons will be
// forming the diquark cluster if kinematically viable)
for (int i = 1; i <= N; i++) {
for (int j = 1; j <= N; j++) {
// TODO add here a check if we can make Diquarkclsuter
// std::cout << "i = "<<i-1 <<" j = "<<j-1 << std::endl;
// std::cout << "c1= "<<i%3 <<" a1= "<<j%3 << std::endl;
// std::cout << "c2= "<<(i+1)%3 <<" a2= "<<(j+1)%3 << std::endl;
assert(i%3!=(i-1));
assert((i+1)%3!=(i-1));
assert(j%3!=(j-1));
assert((j+1)%3!=(j-1));
if (_canMakeDiquarkCluster(
clusters[i%3]->colParticle(), clusters[(i+1)%3]->colParticle(),
clusters[j%3]->antiColParticle(), clusters[(j+1)%3]->antiColParticle() ,PhaseSpace))
{
state = -(10*i+j);
if (_isColour8Forbidden(state,clusters))
continue;
switch (i)
{
case 1:
{
perm1 = 100*j + 10*(((j+1)%3)+1) + (((j)%3)+1);
perm2 = 100*j + 10*(((j)%3)+1) + (((j+1)%3)+1);
amp2 = sqr((amplitudes[perm1] - amplitudes[perm2])/ND);
amp2*=NDiqFACT;
sum2 += amp2;
if (amp2==0 || amplitudes[perm2]==0 || amplitudes[perm1]==0 ){
std::cout << "amp23 = "<< amp2 << std::endl;
std::cout << "perm13 = "<< perm1 << std::endl;
std::cout << "perm23 = "<< perm2 << std::endl;
}
assert((((j)%3)+1)!=j);
assert((((j+1)%3)+1)!=j);
CRoptions.insert(PhaseSpace*amp2, state);
break;
}
case 2:
{
perm1 = 100*(((j)%3)+1) + 10*j + (((j+1)%3)+1);
perm2 = 100*(((j+1)%3)+1) + 10*j + (((j)%3)+1);
amp2 = sqr((amplitudes[perm1] - amplitudes[perm2])/ND);
amp2*=NDiqFACT;
sum2 += amp2;
if (amp2==0 || amplitudes[perm2]==0 || amplitudes[perm1]==0 ){
std::cout << "amp22 = "<< amp2 << std::endl;
std::cout << "perm12 = "<< perm1 << std::endl;
std::cout << "perm22 = "<< perm2 << std::endl;
}
assert((((j)%3)+1)!=j);
assert((((j+1)%3)+1)!=j);
CRoptions.insert(PhaseSpace*amp2, state);
break;
}
case 3:
{
perm1 = 100*(((j+1)%3)+1) + 10*(((j)%3)+1) + j;
perm2 = 100*(((j)%3)+1) + 10*(((j+1)%3)+1) + j;
amp2 = sqr((amplitudes[perm1] - amplitudes[perm2])/ND);
amp2*=NDiqFACT;
sum2 += amp2;
if (amp2==0 || amplitudes[perm2]==0 || amplitudes[perm1]==0 ){
std::cout << "amp21 = "<< amp2 << std::endl;
std::cout << "perm11 = "<< perm1 << std::endl;
std::cout << "perm21 = "<< perm2 << std::endl;
}
assert((((j)%3)+1)!=j);
assert((((j+1)%3)+1)!=j);
CRoptions.insert(PhaseSpace*amp2, state);
break;
}
default:
assert(false);
}
// std::cout << CRoptions << std::endl;
}
}
}
}
// double BaryonRate = sumBaryon*sumBaryon/sum2;
// std::cout << "BaryonRate "<<BaryonRate <<"\n";
// std::cout << "NoRecoRate "<<sqr(amplitudes[_stateToPermutation(1)])/sum2 <<"\n";
// std::cout << "MeRecoRate "<<1-BaryonRate-sqr(amplitudes[_stateToPermutation(1)])/sum2 <<"\n";
CRoptions.insert(sumBaryon*sumBaryon, _stateToPermutation(0));
//TODO insert diquarks
return CRoptions;
}
// TODO HERE add context to the int nonTrivialInitialState -> |D> or similar
std::unordered_map<int, double> ColourReconnector::_reconnectionAmplitudesSGE(const ClusterVector & clusters, int) const {
int size=clusters.size();
assert(clusters.size()<4);
switch(size){
case 2:
{
const std::unordered_map<int, double> & amplitudes = _reconnectionAmplitudesCF2(clusters[0],clusters[1]);
return amplitudes;
}
case 3:
{
std::unordered_map<int, double> amplitudes;
if (clusters[0]->numComponents()!=2 ||
clusters[1]->numComponents()!=2 ||
clusters[2]->numComponents()!=2 ||
hasDiquark(clusters[0]) ||
hasDiquark(clusters[1]) ||
hasDiquark(clusters[2]) ) {
std::cout << "reject config. should reject before somehow\n";
return amplitudes;
}
const int N=6; // 2*Nclu;
double omIJ[N][N];
double alphaQCD=_dynamicCRalphaS;
Lorentz5Momentum mom_i, mom_j;
// Assumption in notebook analytics/matrices2_BCR.nb
// Ordering of omIJ and scalarProds is according to {clu1_col, clu1_anti, clu2_col, clu2_anti,...}
Energy mLightestConstMass = 1000*GeV;
for (auto id : _hadronSpectrum->lightHadronizingQuarks() ) {
if (getParticleData(id)->constituentMass()<mLightestConstMass)
mLightestConstMass = getParticleData(id)->constituentMass();
}
for (int i = 0; i < N; i++)
{
if ((i+1)%2==1) mom_i=clusters[i/2]->colParticle()->momentum();
else mom_i=clusters[(i-1)/2]->antiColParticle()->momentum();
for (int j = i+1; j < N; j++)
{
if ((j+1)%2==1) mom_j=clusters[j/2]->colParticle()->momentum();
else mom_j=clusters[(j-1)/2]->antiColParticle()->momentum();
// TODO REFACTOR THIS TO ALLOW FOR DIFFERENT SCALE CHOICES
double rho = (mom_i+mom_j).m2()/sqr(2*mLightestConstMass*_dynamicCRscale);
if (rho < 1.0) {
throw Exception()
<< "DynamicCR scale "<< _dynamicCRscale << " too high in ColourReconnector"
<< " must be less than 1"
<< " Found parton invariant mass combinations with invariant mass ratio "<< (mom_i+mom_j).m2()/sqr((mom_i.m()+mom_j.m()))
<< Exception::eventerror;
}
// TODO missing factor of two due to sum_i!=j
// can be fixed by twopi->pi
omIJ[i][j]=alphaQCD*pow(log(rho),2)/(2.0*twopi);
}
}
boost::numeric::ublas::matrix<double> * Uevolve = new boost::numeric::ublas::matrix<double>(N,N);
boost::numeric::ublas::matrix<double> Omega(N,N);
// TODO need to restore Nc dependence
int Nc = 3;
// Verified Omega Matrix with analytics/matrices2_BCR.nb
Omega(0,0) = - 0.5 * Nc * (omIJ[0][1]+omIJ[2][3]+omIJ[4][5]);
Omega(1,0) = 0.5 * (omIJ[2][4]-omIJ[3][4]-omIJ[2][5]+omIJ[3][5]);
Omega(2,0) = 0.5 * (omIJ[0][2]-omIJ[1][2]-omIJ[0][3]+omIJ[1][3]);
Omega(3,0) = 0.0;
Omega(4,0) = 0.0;
Omega(5,0) = 0.5 * (omIJ[0][4]-omIJ[1][4]-omIJ[0][5]+omIJ[1][5]);
Omega(0,1) = - 0.5 * (omIJ[2][3]-omIJ[2][4]-omIJ[3][5]+omIJ[4][5]);
Omega(1,1) = - 0.5 * Nc * (omIJ[0][1]+omIJ[3][4]+omIJ[2][5]);
Omega(2,1) = 0.0;
Omega(3,1) = - 0.5 * (omIJ[0][3]-omIJ[1][3]-omIJ[0][4]+omIJ[1][4]);
Omega(4,1) = 0.5 * (omIJ[0][2]-omIJ[1][2]-omIJ[0][5]+omIJ[1][5]);
Omega(5,1) = 0.0;
Omega(0,2) = - 0.5 * (omIJ[0][1]-omIJ[0][2]-omIJ[1][3]+omIJ[2][3]);
Omega(1,2) = 0.0;
Omega(2,2) = - 0.5 * Nc * (omIJ[1][2]+omIJ[0][3]+omIJ[4][5]);
Omega(3,2) = - 0.5 * (omIJ[1][4]-omIJ[2][4]-omIJ[1][5]+omIJ[2][5]);
Omega(4,2) = 0.5 * (omIJ[0][4]-omIJ[3][4]-omIJ[0][5]+omIJ[3][5]);
Omega(5,2) = 0.0;
Omega(0,3) = 0.0;
Omega(1,3) = - 0.5 * (omIJ[0][1]-omIJ[1][3]-omIJ[0][4]+omIJ[3][4]);
Omega(2,3) = - 0.5 * (omIJ[1][2]-omIJ[2][4]-omIJ[1][5]+omIJ[4][5]);
Omega(3,3) = - 0.5 * Nc * (omIJ[0][3]+omIJ[1][4]+omIJ[2][5]);
Omega(4,3) = 0.0;
Omega(5,3) = 0.5 * (omIJ[0][2]-omIJ[2][3]-omIJ[0][5]+omIJ[3][5]);
Omega(0,4) = 0.0;
Omega(1,4) = - 0.5 * (omIJ[0][1]-omIJ[0][2]-omIJ[1][5]+omIJ[2][5]);
Omega(2,4) = - 0.5 * (omIJ[0][3]-omIJ[0][4]-omIJ[3][5]+omIJ[4][5]);
Omega(3,4) = 0.0;
Omega(4,4) = - 0.5 * Nc * (omIJ[1][2]+omIJ[3][4]+omIJ[0][5]);
Omega(5,4) = 0.5 * (omIJ[1][3]-omIJ[2][3]-omIJ[1][4]+omIJ[2][4]);
Omega(0,5) = - 0.5 * (omIJ[0][1]-omIJ[0][4]-omIJ[1][5]+omIJ[4][5]);
Omega(1,5) = 0.0;
Omega(2,5) = 0.0;
Omega(3,5) = 0.5 * (omIJ[0][2]-omIJ[0][3]-omIJ[2][5]+omIJ[3][5]);
Omega(4,5) = - 0.5 * (omIJ[1][2]-omIJ[1][3]-omIJ[2][4]+omIJ[3][4]);
Omega(5,5) = - 0.5 * Nc * (omIJ[2][3]+omIJ[1][4]+omIJ[0][5]);
switch (_dynamicCR)
{
case 1:
// Exponentiated
*Uevolve=expm_pad(Omega);
break;
case 2:
// NotExponentiated
*Uevolve=Omega;
break;
default:
assert(false);
}
// std::cout << *Uevolve << std::endl;
// std::vector<double> Transition1toJ; // |<J|U|1>|^2
double amp1toJ;
for (int J = 0; J < N; J++)
{
// amp1toJ is here for each J transition amplitude 1->J
amp1toJ=0;
for (int i = 0; i < N; i++)
{
// amp1toJ corresponds to the Uevolve operator applied to the |1> state
// and projected by fixed <J|
// The resulting amp1toJ is <J|U|1>
// Should be correct this way
// TODO : Generalize here to general starting states i.e. |1> -> |initial>
amp1toJ+=(*Uevolve)(i,0)*scalarProducts(i,J);
}
if (std::isnan(amp1toJ) || std::isinf(amp1toJ)){
throw Exception() << "nan or inf transition probability in ColourReconnector::_reconnectionAmplitudesSGE"
<< Exception::runerror;
}
amplitudes[_stateToPermutation(J+1)] = amp1toJ;
}
delete Uevolve;
return amplitudes;
break;
}
default:
{
throw Exception() << "Found cluster set of "<<size<<" in ColourReconnector::_reconnectionAmplitudesSGE (can only handle 2 or 3 colour Flows)"
<< Exception::runerror;
}
}
std::unordered_map<int, double> amplitudes;
return amplitudes;
}
double ColourReconnector::_displacementBaryonic(tcPPtr q1, tcPPtr q2, tcPPtr q3) const {
if (_junctionMBCR) {
/**
* Junction-like option i.e. displacement
* from "junction centre" (mean rapidity/phi)
*/
double rap1=q1->rapidity();
double rap2=q2->rapidity();
double rap3=q3->rapidity();
double phi1=q1->momentum().phi();
double phi2=q2->momentum().phi();
double phi3=q3->momentum().phi();
double meanRap=(rap1 + rap2 + rap3)/3.0;
// Use circularMean for defining a sensible mean of an angle
double meanPhi=circularMean(phi1,phi2,phi3);
double deltaPhi1=fabs(phi1-meanPhi);
double deltaPhi2=fabs(phi2-meanPhi);
double deltaPhi3=fabs(phi3-meanPhi);
// keep deltaPhi's below Pi due to periodicity
if (deltaPhi1>M_PI) deltaPhi1-=M_PI;
if (deltaPhi2>M_PI) deltaPhi2-=M_PI;
if (deltaPhi3>M_PI) deltaPhi3-=M_PI;
double delR;
delR = sqrt( (rap1-meanRap)*(rap1-meanRap) + deltaPhi1*deltaPhi1 );
delR += sqrt( (rap2-meanRap)*(rap2-meanRap) + deltaPhi2*deltaPhi2 );
delR += sqrt( (rap3-meanRap)*(rap3-meanRap) + deltaPhi3*deltaPhi3 );
return delR;
} else {
/* just summing up all possible 2 quark displacements */
return _displacement(q1, q2) + _displacement(q1, q3) + _displacement(q2, q3);
}
}
bool ColourReconnector::_containsColour8(const ClusterVector & cv,
const vector<size_t> & P) const {
assert (P.size() == cv.size());
for (size_t i = 0; i < cv.size(); i++) {
tcPPtr p = cv[i]->colParticle();
tcPPtr q = cv[P[i]]->antiColParticle();
if (_isColour8(p, q)) return true;
}
return false;
}
void ColourReconnector::_doRecoStatistical(ClusterVector & cv) const {
const size_t nclusters = cv.size();
// initially, enumerate (anti)quarks as given in the cluster vector
ParticleVector q, aq;
for (size_t i = 0; i < nclusters; i++) {
q.push_back( cv[i]->colParticle() );
aq.push_back( cv[i]->antiColParticle() );
}
// annealing scheme
Energy2 t, delta;
Energy2 lambda = _clusterMassSum(q,aq);
const unsigned _ntries = _triesPerStepFactor * nclusters;
// find appropriate starting temperature by measuring the largest lambda
// difference in some dry-run random rearrangements
{
vector<Energy2> typical;
for (int i = 0; i < 10; i++) {
const pair <int,int> toswap = _shuffle(q,aq,5);
ParticleVector newaq = aq;
swap (newaq[toswap.first], newaq[toswap.second]);
Energy2 newlambda = _clusterMassSum(q,newaq);
typical.push_back( abs(newlambda - lambda) );
}
t = _initTemp * Math::median(typical);
}
// anneal in up to _annealingSteps temperature steps
for (unsigned step = 0; step < _annealingSteps; step++) {
// For this temperature step, try to reconnect _ntries times. Stop the
// algorithm if no successful reconnection happens.
unsigned nSuccess = 0;
for (unsigned it = 0; it < _ntries; it++) {
// make a random rearrangement
const unsigned maxtries = 10;
const pair <int,int> toswap = _shuffle(q,aq,maxtries);
const int i = toswap.first;
const int j = toswap.second;
// stop here if we cannot find any allowed reconfiguration
if (i == -1) break;
// create a new antiquark vector with the two partons swapped
ParticleVector newaq = aq;
swap (newaq[i], newaq[j]);
// Check if lambda would decrease. If yes, accept the reconnection. If no,
// accept it only with a probability given by the current Boltzmann
// factor. In the latter case we set p = 0 if the temperature is close to
// 0, to avoid division by 0.
Energy2 newlambda = _clusterMassSum(q,newaq);
delta = newlambda - lambda;
double prob = 1.0;
if (delta > ZERO) prob = ( abs(t) < 1e-8*MeV2 ) ? 0.0 : exp(-delta/t);
if (UseRandom::rnd() < prob) {
lambda = newlambda;
swap (newaq, aq);
nSuccess++;
}
}
if (nSuccess == 0) break;
// reduce temperature
t *= _annealingFactor;
}
// construct the new cluster vector
ClusterVector newclusters;
for (size_t i = 0; i < nclusters; i++) {
ClusterPtr cl = new_ptr( Cluster( q[i], aq[i] ) );
newclusters.push_back(cl);
}
swap(newclusters,cv);
return;
}
Selector <int> ColourReconnector::getProbabilities2CF(const ClusterPtr & c1, const ClusterPtr & c2, bool diquarkCR) const{
Selector <int> res;
if (_dynamicCR) {
double pNoCR;
double pMCR;
double pDCR;
std::tie(pNoCR, pMCR, pDCR) = _dynamicRecoProbabilitiesCF2(c1, c2, diquarkCR);
// TODO NOT GOOD orientation
// int orient= orientation(c1,c2);
// std::cout <<"orient = "<<orient<<"\t" << pMCR << "\t"<<pDCR << "\t" << pNoCR<< "\n";
// Mesonic CR
// if ( orient>0 && !( _isColour8(c1->colParticle(), c2->antiColParticle())
if ( !( _isColour8(c1->colParticle(), c2->antiColParticle())
|| _isColour8(c2->colParticle(), c1->antiColParticle()) ) ) {
res.insert(pMCR, 213);
}
// if ( orient<0 && diquarkCR && pDCR>0.0 ) {
if (diquarkCR && pDCR>0.0 ) {
// Diquark CR
res.insert(pDCR, -213);
}
// No CR
res.insert(pNoCR, 123);
}
else {
double wMCR = _preco;
double sum = 0.0;
// Mesonic CR
if ( !( _isColour8(c1->colParticle(), c2->antiColParticle())
|| _isColour8(c2->colParticle(), c1->antiColParticle()) ) ) {
res.insert(wMCR, 213);
sum+=wMCR;
}
double PhaseSpace;
if (diquarkCR && _canMakeDiquarkCluster(c1,c2,PhaseSpace)) {
double wDCR = _precoDiquark*PhaseSpace;
// Diquark CR
res.insert(wDCR, -213);
sum+=wDCR;
}
// No CR
res.insert((1.0-sum)>0 ? (1.0-sum):0.0, 123);
}
return res;
}
void ColourReconnector::_doRecoPlain(ClusterVector & cv) const {
ClusterVector newcv = cv;
ClusterVector deleted; deleted.reserve(cv.size());
// try to avoid systematic errors by randomising the reconnection order
long (*p_irnd)(long) = UseRandom::irnd;
random_shuffle( newcv.begin(), newcv.end(), p_irnd );
// iterate over all clusters
for (CluVecIt cit = newcv.begin(); cit != newcv.end(); cit++) {
// find the cluster which, if reconnected with *cit, would result in the
// smallest sum of cluster masses
// skip the diquark clusters to be deleted later 2->1 cluster
if (find(deleted.begin(), deleted.end(), *cit) != deleted.end())
continue;
// skip diquark clusters
if ((*cit)->numComponents()>2 || (_diquarkCR && hasDiquark(*cit))) continue;
// NB this method returns *cit if no reconnection partner can be found
CluVecIt candidate = _dynamicCR ? _findRecoPartnerPlainDynamic(cit, newcv, deleted, _diquarkCR>0):_findRecoPartnerPlain(cit, newcv, deleted);
// skip this cluster if no possible reshuffling partner can be found
if (candidate == cit) continue;
// accept the reconnection with probability PrecoProb.
ClusterVector cluvec = {*cit, *candidate};
// const Selector <int> & sel = getProbabilities2CF(*cit, *candidate, _diquarkCR>0);
const Selector <int> & sel = _selector(cluvec, _diquarkCR>0);
enum Selection2ColourFlows {
NoReconnection = 123,
MesonicReconnection = 213,
DiquarkReconnection = -213,
};
switch (sel.select(UseRandom::rnd()))
{
case MesonicReconnection:
{
// Mesonic Colour Reconnection
pair <ClusterPtr,ClusterPtr> reconnected = _reconnect(*cit, *candidate);
// Replace the clusters in the ClusterVector. The order of the
// colour-triplet partons in the cluster vector is retained here.
// replace *cit by reconnected.first
*cit = reconnected.first;
// replace candidate by reconnected.second
*candidate = reconnected.second;
break;
}
case DiquarkReconnection:
{
ClusterPtr DiqCluster;
if (_makeDiquarkCluster(*cit, *candidate, DiqCluster)){
deleted.push_back(*candidate);
// Note that these must be the cit,candidate and not cluvec
*cit = DiqCluster;
}
break;
}
case NoReconnection:
// No colour Reconnection
break;
default:
assert(false);
}
}
if (deleted.size()==0) {
swap(cv, newcv);
}
else {
// create a new vector of clusters except for the ones which are "deleted" during
// Diquark reconnection
ClusterVector clustervector;
for (const auto & cluster : newcv)
if (find(deleted.begin(),
deleted.end(), cluster) == deleted.end())
clustervector.push_back(cluster);
swap(cv, clustervector);
}
}
void ColourReconnector::_doRecoBaryonicDiquarkCluster(ClusterVector & cv) const {
/*START REVIEW*/
ClusterVector newcv = cv;
ClusterVector deleted; deleted.reserve(cv.size());
// try to avoid systematic errors by randomising the reconnection order
long (*p_irnd)(long) = UseRandom::irnd;
random_shuffle(newcv.begin(), newcv.end(), p_irnd);
Selector <int> sel;
ClusterVector cluvec;
int selection;
// iterate over all clusters
for (CluVecIt cit = newcv.begin(); cit != newcv.end(); ++cit) {
//avoid clusters already containing diuarks
if (hasDiquark(*cit)) continue;
//skip the cluster to be deleted later 3->2 cluster
if (find(deleted.begin(), deleted.end(), *cit) != deleted.end())
continue;
// Skip all found baryonic and Tetra clusters, this biases the
// algorithm but implementing something like re-reconnection
// is ongoing work
if ((*cit)->numComponents()>=3) continue;
// switch ((*cit)->numComponents())
// {
// case 2:
// break;
// case 3:
// // leave full baryonic clusters in peace
// continue;
// case 4:
// // Allow diquark clusters to split again (for multiple iterations)
// {
// // TODO: replace magic numbers with their dynamic version
// double probDiquarkSplit = 0.8;
// double probabilityCRafterSplit = 0.5;
// if (UseRandom::rnd()>probDiquarkSplit){
// const auto & res=_splitDiquarkCluster(*cit, UseRandom::rnd()>probabilityCRafterSplit);
// *cit = res.first;
// newcv.push_back(res.second);
// }
// continue;
// }
// default:
// assert(false);
// }
// Find a candidate suitable for reconnection
CluVecIt candidate1, candidate2;
unsigned typeOfReconnection = 0;
// TODO fix the function below yields Cluster in final state
switch (_selectionChoice)
{
case 1:
// Default: rapidity based similar to Baryonic
_findPartnerBaryonicDiquarkCluster(cit, newcv,
typeOfReconnection,
deleted,
candidate1,
candidate2);
break;
case 2:
// Experimental2: Sorting in rapidity sum
_findPartnerBaryonicDiquarkClusterTEST2(cit, newcv,
typeOfReconnection,
deleted,
candidate1,
candidate2);
break;
case 3:
case 4:
// Experimental3: Sorting in mass sum
_findPartnerBaryonicDiquarkClusterTEST3(cit, newcv,
typeOfReconnection,
deleted,
candidate1,
candidate2);
break;
default:
assert(false);
}
switch (typeOfReconnection)
{
case 0:
// No CR found
continue;
case 3:
case 1:
// Mesonic or Diquark CR with 2 Colour Flows
cluvec={*cit,*candidate1};
break;
case 2:
// Baryonic CR with 3 Colour Flows
cluvec={*cit,*candidate1,*candidate2};
break;
default:
assert(false);
}
if (candidate2!=cit) {
cluvec={*cit,*candidate1,*candidate2};
}
else if (candidate1!=cit) {
cluvec={*cit,*candidate1};
}
else {
std::cout << "no CR" << std::endl;
continue;
}
if (_dynamicCR) {
sel = _selector(cluvec, _diquarkCR>0);
if (typeOfReconnection!=0 && sel.empty()){
throw Exception()
<< "No Selection availible"
<< "ColourReconnector::_doRecoBaryonicDiquarkCluster()"
<< Exception::runerror;
}
selection = sel.select(UseRandom::rnd());
sel.clear();
assert(sel.empty());
}
else {
switch (typeOfReconnection)
{
case 0:
// No CR
std::cout << "Should never execute!!!" << std::endl;
selection = 123;
break;
case 1:
// Mesonic CR with 2 Colour Flows
if (UseRandom::rnd() < _preco)
selection = 213;
else
selection = 123;
break;
case 2:
// Baryonic CR with 3 Colour Flows
- if (UseRandom::rnd() < _precoBaryonic)
+ if (UseRandom::rnd() < _precoBaryonic
+ && _canMakeBaryonicCluster(*cit,*candidate1,*candidate2))
selection = 0;
else
selection = 123;
break;
case 3:
// Diquark CR with 2 Colour Flows
if (UseRandom::rnd() < _precoDiquark)
selection = -1000;
else
selection = 123;
break;
default:
assert(false);
}
}
if (selection == 0){
// Baryonic CR (only one option for 3 colourflows)
deleted.push_back(*candidate2);
// Function that does the reconnection from 3 -> 2 clusters
ClusterPtr b1, b2;
_makeBaryonicClusters(*cit,*candidate1,*candidate2, b1, b2);
// Note that these must be the cit,candidate1 and not cluvec
*cit = b1;
*candidate1 = b2;
}
else if (selection > 0) {
// Mesonic CR
int c1 = ((selection/100) % 4) -1;
int c2 = ((selection-(c1+1)*100)/10 % 4) -1;
int c3 = ((selection-(c1+1)*100-(c2+1)*10) % 4) -1;
if ( c1==0
&& c2==1
&& c3==2){
// noCR
continue;
}
assert(c1>=0 && c1<3);
assert(c2>=0 && c2<3);
assert(c3>=0 && c3<3);
assert(c1!=c2 && c1!=c3 && c2 !=c3);
// if last cluster is untouched we do only reconnect the first two clusters
if (c3==2) {
const auto & reconnected = _reconnect(*cit,*candidate1);
// Note that these must be the cit, candidate1 and not cluvec
*cit = reconnected.first;
*candidate1 = reconnected.second;
}
else {
// Form clusters (0,c1) (1,c2) (2,c3)
const int infoMCR[3] = {c1, c2, c3};
const auto & reconnected = _reconnect3Mto3M(*cit,*candidate1,*candidate2, infoMCR);
// Note that these must be the cit,candidateX and not cluvec
*cit = std::get<0>(reconnected);
*candidate1 = std::get<1>(reconnected);
*candidate2 = std::get<2>(reconnected);
}
}
else if (selection < 0) {
//TODO DiquarkCR
// We will delete the candidate1 mesonic clusters
// need to check with 2CF solution
// to form a diquark cluster
if (cluvec.size()==3 && -selection<999) {
const auto & reconnected = _reconnect3MtoMD(cluvec, -selection);
*cit = std::get<0>(reconnected);
*candidate1 = std::get<1>(reconnected);
deleted.push_back(*candidate2);
// *candidate2 = std::get<2>(reconnected);
}
else if (cluvec.size()==2 || -selection>999){
ClusterPtr DiqCluster;
if (_makeDiquarkCluster(*cit, *candidate1, DiqCluster)){
deleted.push_back(*candidate1);
// Note that these must be the cit,candidate and not cluvec
*cit = DiqCluster;
}
}
else {
assert(false);
}
// auto pq1 = (*cit)->particle(0)->momentum();
// auto pq2 = (*cit)->particle(1)->momentum();
// auto pqbar1 = (*cit)->particle(2)->momentum();
// auto pqbar2 = (*cit)->particle(3)->momentum();
// // auto pDi = pq1 + pq2;
// // auto pantiDi = pqbar1 + pqbar2;
// double mDiff1 = sqrt(pq1*pq2/(pq1.m()*pq2.m()) - 1.0);
// double mDiff2 = sqrt(pqbar1*pqbar2/(pqbar1.m()*pqbar2.m()) - 1.0);
// Boost clboost = (*cit)->momentum().boostVector();
// pq1.boost(-clboost);
// pq2.boost(-clboost);
// pqbar1.boost(-clboost);
// pqbar2.boost(-clboost);
// double phi12 = acos(pq1.vect().cosTheta(pq2.vect()));
// double phi34 = acos(pqbar1.vect().cosTheta(pqbar2.vect()));
// if ( phi12>M_PI/2.0 || phi34>M_PI/2.0){
// std::cout << "WARNING: mDiff1 = " << mDiff1 << std::endl;
// std::cout << "WARNING: mDiff2 = " << mDiff2 << std::endl;
// std::cout << "WARNING: phi12 = " << phi12 << std::endl;
// std::cout << "WARNING: phi34 = " << phi34 << std::endl;
// std::cout << "WARNING: cosTheta1 = " << pq1.vect().cosTheta(pq2.vect()) << std::endl;
// std::cout << "WARNING: cosTheta2 = " << pqbar1.vect().cosTheta(pqbar2.vect()) << std::endl;
// }
}
else {
std::cout << "\nError in selection = "<<selection<<"\n" << std::endl;
assert(false);
}
}
ClusterVector addedcv;
// bool splitSomeDiquarkCluster=true;
// bool splitSomeDiquarkCluster=false;
// if (splitSomeDiquarkCluster) {
// for (CluVecIt cit = newcv.begin(); cit != newcv.end(); ++cit){
// if (!*cit || hasDiquark(*cit)) continue;
// if (find(deleted.begin(), deleted.end(), *cit) != deleted.end())
// continue;
// if ((*cit)->numComponents()==4 && UseRandom::rnd()>0.33){
// const auto & res=_splitDiquarkCluster(*cit, UseRandom::rnd()>0.5);
// *cit = res.first;
// addedcv.push_back(res.second);
// }
// }
// }
// create a new vector of clusters except for the ones which are "deleted" during
// baryonic reconnection
ClusterVector clustervector;
// add new clusters
for (const auto & c : addedcv)
clustervector.push_back(c);
// delete deleted clusters
for (const auto & cluster : newcv)
if (find(deleted.begin(),
deleted.end(), cluster) == deleted.end())
clustervector.push_back(cluster);
swap(cv, clustervector);
}
// Implementation of the baryonic reconnection algorithm
void ColourReconnector::_doRecoBaryonic(ClusterVector & cv) const {
ClusterVector newcv = cv;
ClusterVector deleted; deleted.reserve(cv.size());
// try to avoid systematic errors by randomising the reconnection order
long (*p_irnd)(long) = UseRandom::irnd;
random_shuffle(newcv.begin(), newcv.end(), p_irnd);
double ProbabilityMesonic = _preco;
double ProbabilityBaryonic = _precoBaryonic;
ClusterVector cluvec;
cluvec.reserve(3);
// iterate over all clusters
for (CluVecIt cit = newcv.begin(); cit != newcv.end(); ++cit) {
//avoid clusters already containing diuarks
if (hasDiquark(*cit)) continue;
//skip the cluster to be deleted later 3->2 cluster
if (find(deleted.begin(), deleted.end(), *cit) != deleted.end())
continue;
// Skip all found baryonic clusters, this biases the algorithm but implementing
// something like re-reconnection is ongoing work
if ((*cit)->numComponents()>=3) continue;
// Find a candidate suitable for reconnection
CluVecIt baryonic1, baryonic2;
bool isBaryonicCandidate = false;
CluVecIt candidate = _findPartnerBaryonic(cit, newcv,
isBaryonicCandidate,
deleted,
baryonic1, baryonic2);
// skip this cluster if no possible reconnection partner can be found
if ( !isBaryonicCandidate && *candidate==*cit)
continue;
if (_dynamicCR) {
cluvec.clear();
if (isBaryonicCandidate) {
cluvec={*cit,*baryonic1,*baryonic2};
}
else{
cluvec={*cit,*candidate};
}
const Selector <int> & sel = _selector(cluvec, _diquarkCR>0);
int selection = sel.select(UseRandom::rnd());
// std::cout << "selector print = " << sel <<"\n";
// sel.clear();
// assert(sel.empty());
// 3 aligned meson case
// Normal 2->2 Colour reconnection
if (selection == 0){
// Baryonic CR only one option
deleted.push_back(*baryonic2);
// Function that does the reconnection from 3 -> 2 clusters
ClusterPtr b1, b2;
_makeBaryonicClusters(*cit,*baryonic1,*baryonic2, b1, b2);
*cit = b1;
*baryonic1 = b2;
}
else if (selection > 0) {
//TODO Mesonic CR
if (isBaryonicCandidate) {
int c1 = ((selection/100) % 4) -1;
int c2 = ((selection-(c1+1)*100)/10 % 4) -1;
int c3 = ((selection-(c1+1)*100-(c2+1)*10) % 4) -1;
if ( c1==0
&& c2==1
&& c3==2){
// noCR
continue;
}
assert(c1>=0 && c1<3);
assert(c2>=0 && c2<3);
assert(c3>=0 && c3<3);
assert(c1!=c2 && c1!=c3 && c2 !=c3);
if (c3==2 ) {
const auto & reconnected = _reconnect(*cit,*baryonic1);
*cit = reconnected.first;
*baryonic1 = reconnected.second;
}
else {
// Form clusters (0,c1) (1,c2) (2,c3)
const int infoMCR[3] = {c1, c2, c3};
const auto & reconnected = _reconnect3Mto3M(*cit,*baryonic1,*baryonic2, infoMCR);
*cit = std::get<0>(reconnected);
*baryonic1 = std::get<1>(reconnected);
*baryonic2 = std::get<2>(reconnected);
}
}
else {
const auto & reconnected = _reconnect(*cit,*candidate);
*cit = reconnected.first;
*candidate = reconnected.second;
}
}
else if (selection < 0) {
//TODO DiquarkCR
// We will delete the candidate1 mesonic clusters
// need to check with 2CF solution
// to form a diquark cluster
if (cluvec.size()==3) {
const auto & reconnected = _reconnect3MtoMD(cluvec, -selection);
*cit = std::get<0>(reconnected);
*baryonic1 = std::get<1>(reconnected);
deleted.push_back(*baryonic2);
}
else if (cluvec.size()==2) {
ClusterPtr DiqCluster;
if (_makeDiquarkCluster(*cit, *candidate, DiqCluster)){
deleted.push_back(*candidate);
// Note that these must be the cit,candidate and not cluvec
*cit = DiqCluster;
}
}
}
else {
std::cout << "\nError in selection = "<<selection<<"\n" << std::endl;
assert(false);
}
}
else {
// 3 aligned meson case
if ( isBaryonicCandidate
- && UseRandom::rnd() < ProbabilityBaryonic ) {
+ && UseRandom::rnd() < ProbabilityBaryonic
+ && _canMakeBaryonicCluster(*cit, *baryonic1, *baryonic2)) {
deleted.push_back(*baryonic2);
// Function that does the reconnection from 3 -> 2 clusters
ClusterPtr b1, b2;
_makeBaryonicClusters(*cit, *baryonic1, *baryonic2, b1, b2);
*cit = b1;
*baryonic1 = b2;
// Baryonic2 is easily skipped in the next loop
}
// Normal 2->2 Colour reconnection
if (!isBaryonicCandidate
&& UseRandom::rnd() < ProbabilityMesonic) {
const auto & reconnected = _reconnect(*cit, *candidate);
*cit = reconnected.first;
*candidate = reconnected.second;
}
}
}
// create a new vector of clusters except for the ones which are "deleted" during
// baryonic reconnection
ClusterVector clustervector;
for (const auto & cluster : newcv)
if (find(deleted.begin(),
deleted.end(), cluster) == deleted.end())
clustervector.push_back(cluster);
swap(cv, clustervector);
}
bool ColourReconnector::_clustersFarApart( const std::vector<CluVecIt> & clu ) const {
int Ncl=clu.size();
assert(Ncl<=3);
if (Ncl==1) {
return false;
} else if (Ncl==2) {
// veto if Clusters further apart than _maxDistance
if (_localCR && ((*clu[0])->vertex().vect()-(*clu[1])->vertex().vect()).mag() > _maxDistance) return true;
// veto if Clusters have negative spacetime difference
if (_causalCR && ((*clu[0])->vertex()-(*clu[1])->vertex()).m() < ZERO) return true;
} else if (Ncl==3) {
// veto if Clusters further apart than _maxDistance
if (_localCR && ((*clu[0])->vertex().vect()-(*clu[1])->vertex().vect()).mag() > _maxDistance) return true;
if (_localCR && ((*clu[1])->vertex().vect()-(*clu[2])->vertex().vect()).mag() > _maxDistance) return true;
if (_localCR && ((*clu[0])->vertex().vect()-(*clu[2])->vertex().vect()).mag() > _maxDistance) return true;
// veto if Clusters have negative spacetime difference
if (_causalCR && ((*clu[0])->vertex()-(*clu[1])->vertex()).m() < ZERO) return true;
if (_causalCR && ((*clu[1])->vertex()-(*clu[2])->vertex()).m() < ZERO) return true;
if (_causalCR && ((*clu[0])->vertex()-(*clu[2])->vertex()).m() < ZERO) return true;
}
return false;
}
void ColourReconnector::_doReco2BeamClusters(ClusterVector & cv) const {
// try other option
tPPtr p1Di=(cv[0])->colParticle();
tPPtr p2Di=(cv[1])->colParticle();
tPPtr p1Q=(cv[0])->antiColParticle();
tPPtr p2Q=(cv[1])->antiColParticle();
double min_dist=_displacement(p1Di,p1Q)+_displacement(p2Di,p2Q);
if ((_displacement(p1Di,p2Q)+_displacement(p1Di,p2Q))<min_dist) {
_reconnect(cv[0],cv[1]);
}
return;
}
void ColourReconnector::_doRecoBaryonicMesonic(ClusterVector & cv) const {
if (cv.size() < 3) {
/*
* if the option _cr2BeamClusters!=0 is chosen then we try to
* colour reconnect the special case of 2 beam clusters with
* probability 1.0 if there is a better _displacement
* */
if( _cr2BeamClusters && cv.size()==2 ) _doReco2BeamClusters(cv);
return;
}
ClusterVector newcv = cv;
newcv.reserve(2*cv.size());
ClusterVector deleted;
deleted.reserve(cv.size());
// counters for numbers of mesons and baryons selected
unsigned num_meson = 0;
unsigned num_baryon = 0;
// vector of selected clusters
std::vector<CluVecIt> sel;
unsigned number_of_tries = _stepFactor*cv.size()*cv.size();
if (number_of_tries<1) number_of_tries=1;
long (*p_irnd)(long) = UseRandom::irnd;
for (unsigned reconnections_tries = 0; reconnections_tries < number_of_tries; reconnections_tries++) {
num_meson = 0;
num_baryon = 0;
// flag if we are able to find a suitable combinations of clusters
bool _found = false;
// Shuffle list of clusters to avoid systematic bias in cluster selection
random_shuffle(newcv.begin(), newcv.end(), p_irnd);
// loop over clustervector to find CR candidates
for (CluVecIt cit = newcv.begin(); cit != newcv.end(); ++cit) {
// skip the clusters to be deleted later from 3->2 cluster CR
if (find(deleted.begin(), deleted.end(), *cit) != deleted.end()) continue;
// avoid clusters already containing diuarks
if (hasDiquark(*cit)) continue;
// add to selection
sel.push_back(cit);
if (_clustersFarApart(sel)) {
// reject far appart CR
sel.pop_back();
continue;
}
bool isMeson=((*cit)->numComponents() == 2);
if ( isMeson && (num_meson ==0|| num_meson==1) && num_baryon ==0) {
num_meson++;
/**
* now we habe either 1 or 2 mesonic clusters and have to continue
*/
continue;
} else if ( isMeson && (num_baryon == 1 || num_meson ==2)) {
num_meson++;
_found = true;
/**
* we have either 3 mesonic or 1 mesonic and 1 baryonic cluster
* and try to colour reconnect
*/
break;
} else if (num_baryon ==0 && num_meson==0) {
num_baryon++;
/**
* now we have 1 baryonic cluster and have to continue
*/
continue;
} else if (num_meson == 2) {
/**
* we already have 2 mesonic clusters and dont want a baryonic one
* since this is an invalid selection
*/
// remove previously added cluster
sel.pop_back();
continue;
} else {
num_baryon++;
_found = true;
/**
* now we have either 2 baryonic clusters or 1 mesonic and 1 baryonic cluster
* and try to colour reconnect
*/
break;
}
}
// added for more efficent rejection if some reco probabilities are 0
if ( _found ) {
// reject MBtoMB candidates if _precoMB_MB=0
if ( _precoMB_MB == 0 && (num_baryon == 1 && num_meson == 1) ) {
_found=false;
}
// reject BbarBto3M candidates if _precoBbarB_3M=0
if ( _precoBbarB_3M== 0 && num_baryon == 2 ) {
bool isBbarBto3Mcandidate=(
(*sel[0])->particle(0)->hasColour() && (*sel[1])->particle(0)->hasColour(true) )
|| ( (*sel[0])->particle(0)->hasColour(true) && (*sel[1])->particle(0)->hasColour() );
if ( isBbarBto3Mcandidate) _found=false;
}
// reject 2Bto2B candidates if _preco2B_2B=0
if ( _preco2B_2B == 0 && num_baryon == 2 ) {
bool is2Bto2Bcandidate=(
(*sel[0])->particle(0)->hasColour() && (*sel[1])->particle(0)->hasColour() )
|| ( (*sel[0])->particle(0)->hasColour(true) && (*sel[1])->particle(0)->hasColour(true) );
if ( is2Bto2Bcandidate ) _found=false;
}
}
// were we able to find a combination?
if (_found==false) {
// clear the selection if we did not find a valid set of clusters
sel.erase(sel.begin(), sel.end());
continue;
}
assert(sel.size()<4);
assert(sel.size()>1);
string kind_of_reco = "";
int reco_info[3];
// find best CR option for the selection
_findbestreconnectionoption(sel, num_baryon, kind_of_reco, reco_info);
if (kind_of_reco == "") {
// no reconnection was found
sel.erase(sel.begin(), sel.end());
continue;
} else if (kind_of_reco == "3Mto3M" && UseRandom::rnd() < _preco3M_3M) {
// 3Mto3M colour reconnection
const auto & reconnected = _reconnect3Mto3M(*sel[0], *sel[1], *sel[2],
reco_info);
(*sel[0]) = std::get<0>(reconnected);
(*sel[1]) = std::get<1>(reconnected);
(*sel[2]) = std::get<2>(reconnected);
} else if (kind_of_reco=="2Bto3M" && UseRandom::rnd() < _precoBbarB_3M) {
// antibaryonic and baryonic to 3 mesonic reconnecion
const auto & reconnected = _reconnectBbarBto3M(*sel[0], *sel[1],
reco_info[0], reco_info[1], reco_info[2]);
(*sel[0]) = std::get<0>(reconnected);
(*sel[1]) = std::get<1>(reconnected);
newcv.push_back(std::get<2>(reconnected));
- } else if (kind_of_reco=="3Mto2B" && UseRandom::rnd() < _preco3M_BBbar) {
+ } else if (kind_of_reco=="3Mto2B"
+ && _canMakeBaryonicCluster(*sel[0], *sel[1], *sel[2])
+ && UseRandom::rnd() < _preco3M_BBbar) {
// 3 mesonic to antibaryonic and baryonic reconnection
ClusterPtr b1, b2;
_makeBaryonicClusters(*sel[0], *sel[1], *sel[2], b1, b2);
(*sel[0]) = b1;
(*sel[1]) = b2;
deleted.push_back(*sel[2]);
} else if (kind_of_reco=="2Bto2B" && UseRandom::rnd() < _preco2B_2B) {
// 2 (anti)baryonic to 2 (anti)baryonic reconnection
const auto & reconnected = _reconnect2Bto2B(*sel[0], *sel[1],
reco_info[0], reco_info[1]);
(*sel[0]) = reconnected.first;
(*sel[1]) = reconnected.second;
} else if (kind_of_reco=="MBtoMB" && UseRandom::rnd() < _precoMB_MB) {
// (anti)baryonic and mesonic to (anti)baryonic and mesonic reconnection
const auto & reconnected = _reconnectMBtoMB(*sel[0], *sel[1],
reco_info[0]);
(*sel[0]) = reconnected.first;
(*sel[1]) = reconnected.second;
}
// erase the sel-vector
sel.erase(sel.begin(), sel.end());
}
// write to clustervector new CR'd clusters and deleting
// all deleted clusters
ClusterVector clustervector;
for (const auto & cluster : newcv)
if (find(deleted.begin(), deleted.end(), cluster) == deleted.end())
clustervector.push_back(cluster);
swap(cv, clustervector);
}
void ColourReconnector::_findbestreconnectionoption(std::vector<CluVecIt> & cls, const unsigned & baryonic,
string & kind_of_reco, int (&reco_info)[3]) const {
double min_displacement;
if (baryonic==0) {
// case with 3 mesonic clusters
assert(cls.size()==3);
// calculate the initial displacement sum
min_displacement = _mesonToBaryonFactor * _displacement((*cls[0])->particle(0), (*cls[0])->particle(1));
min_displacement += _mesonToBaryonFactor * _displacement((*cls[1])->particle(0), (*cls[1])->particle(1));
min_displacement += _mesonToBaryonFactor * _displacement((*cls[2])->particle(0), (*cls[2])->particle(1));
// find best CR reco_info and kind_of_reco
_3MtoXreconnectionfinder(cls,
reco_info[0], reco_info[1], reco_info[2], min_displacement, kind_of_reco);
/**
* kind_of_reco either "3Mto3M" or "3Mto2B" (or "" if no better configuration is found)
* case 3Mto3M: the coloured particle of the i-th cluster forms a new cluster with the
* antiparticle of the reco_info[i]-th cluster
* case 3MtoBbarB: all 3 (anti)coloured particle form a new (anti)baryonic cluster
*/
} else if (baryonic == 1) {
// case 1 baryonic and 1 mesonic cluster
assert(cls.size()==2);
// make mesonic cluster always the cls[0]
if ((*cls[0])->numComponents() == 3) {
ClusterPtr zw = *cls[0];
*cls[0] = *cls[1];
*cls[1] = zw;
}
// calculate the initial displacement sum
min_displacement = _mesonToBaryonFactor *_displacement((*cls[0])->particle(0), (*cls[0])->particle(1));
min_displacement += _displacementBaryonic((*cls[1])->particle(0), (*cls[1])->particle(1), (*cls[1])->particle(2));
// find best CR reco_info and kind_of_reco
_BMtoBMreconnectionfinder(*cls[0], *cls[1],
reco_info[0], min_displacement, kind_of_reco);
/**
* reco_info[0] is the index of the (anti)quarks of the baryonic cluster cls[1], which should
* be swapped with the (anti)quarks of the mesonic cluster cls[0]
*/
} else {
assert(baryonic==2);
assert(cls.size()==2);
// calculate the initial displacement sum
min_displacement = _displacementBaryonic((*cls[0])->particle(0), (*cls[0])->particle(1), (*cls[0])->particle(2));
min_displacement += _displacementBaryonic((*cls[1])->particle(0), (*cls[1])->particle(1), (*cls[1])->particle(2));
// case 2 (anti)baryonic clusters to 2 other (anti)baryonic clusters
if ( ( (*cls[0])->particle(0)->hasColour() && (*cls[1])->particle(0)->hasColour() )
|| ( (*cls[0])->particle(0)->hasColour(true) && (*cls[1])->particle(0)->hasColour(true) ) ) {
// find best CR reco_info and kind_of_reco
_2Bto2BreconnectionFinder(*cls[0], *cls[1],
reco_info[0], reco_info[1], min_displacement, kind_of_reco);
/**
* swap the reco_info[0]-th particle of the first cluster in the vector with the
* reco_info[1]-th particle of the second cluster
*/
} else {
// case 1 baryonic and 1 antibaryonic cluster to 3 mesonic clusters
// find best CR reco_info and kind_of_reco
_BbarBto3MreconnectionFinder(*cls[0], *cls[1],
reco_info[0], reco_info[1], reco_info[2], min_displacement, kind_of_reco);
/**
* the i-th particle of the first cluster form a new mesonic cluster with the
* reco_info[i]-th particle of the second cluster
*/
}
}
return;
}
void ColourReconnector::_findPartnerBaryonicDiquarkCluster(
const CluVecIt & cl, ClusterVector & cv,
unsigned & typeOfReconnection,
const ClusterVector& deleted,
CluVecIt &candidate1,
CluVecIt &candidate2 ) const {
typeOfReconnection=0; // no Reconnection found
using Constants::pi;
using Constants::twopi;
candidate1=cl;
candidate2=cl;
bool candIsOctet1 = false;
bool candIsOctet2 = false;
bool candIsQQ1 = false;
bool candIsQQ2 = false;
bool foundCR = false;
double maxrap1 = 0.0;
double maxrap2 = 0.0;
double minrap1 = 0.0;
double minrap2 = 0.0;
double maxsum1 = 0.0;
double maxsum2 = 0.0;
double NegativeRapidtyThreshold = 0.0;
double PositiveRapidtyThreshold = 0.0;
// boost into RF of cl
Lorentz5Momentum cl1 = (*cl)->momentum();
const Boost boostv(-cl1.boostVector());
// boost constituent of cl into RF of cl
Lorentz5Momentum p1anticol = (*cl)->antiColParticle()->momentum();
p1anticol.boost(boostv);
for (CluVecIt cit=cv.begin(); cit != cv.end(); ++cit) {
//avoid looping over clusters containing diquarks
if ( (*cit)->numComponents()>=3 ) continue;
//skip the cluster to be deleted later 3->2 cluster
if ( find(deleted.begin(), deleted.end(), *cit) != deleted.end() )
continue;
if ( hasDiquark(*cit) ) continue;
if ( (*cl)->isBeamCluster() && (*cit)->isBeamCluster() )
continue;
if ( cit==cl ) continue;
// veto if Clusters further apart than _maxDistance
if (_localCR && ((**cl).vertex().vect()-(**cit).vertex().vect()).mag() > _maxDistance) continue;
// veto if Clusters have negative spacetime difference
if (_causalCR && ((**cl).vertex()-(**cit).vertex()).m() < ZERO) continue;
bool octetNormalCR =
(_isColour8( (*cl)->colParticle(), (*cit)->antiColParticle() )
||
_isColour8( (*cit)->colParticle(), (*cl)->antiColParticle() ) );
// boost constituents of cit into RF of cl
Lorentz5Momentum p2col = (*cit)->colParticle()->momentum();
Lorentz5Momentum p2anticol = (*cit)->antiColParticle()->momentum();
p2col.boost(boostv);
p2anticol.boost(boostv);
// calculate the rapidity of the other constituents of the clusters
// w.r.t axis of p1anticol.vect.unit
const double rapq = calculateRapidityRF(p1anticol,p2col);
const double rapqbar = calculateRapidityRF(p1anticol,p2anticol);
// std::cout << "\nPRE\n"
// << "\t octet = " << octetNormalCR << "\n"
// << "\t rapq = " << rapq << "\n"
// << "\t rapqbar = " << rapqbar << "\n";
// configuration for Mesonic CR
if ( rapq > 0.0 && rapqbar < 0.0
&& rapq > PositiveRapidtyThreshold
&& rapqbar < NegativeRapidtyThreshold) {
//sum of rapidities of quarks
const double sumQQbar = abs(rapq) + abs(rapqbar);
if ( sumQQbar > maxsum2 ) {
if ( sumQQbar > maxsum1 ) {
double factor = candIsQQ1 ? _mesonToBaryonFactor:1.0;
maxsum2 = (factor*maxsum1) > sumQQbar ? sumQQbar:(factor*maxsum1);
candidate2 = candidate1;
candIsQQ2 = candIsQQ1;
candIsOctet2 = candIsOctet1;
maxrap2 = maxrap1;
minrap2 = minrap1;
maxsum1 = sumQQbar;
candidate1 = cit;
candIsQQ1 = false;
candIsOctet1 = octetNormalCR;
maxrap1 = rapq;
minrap1 = rapqbar;
} else {
maxsum2 = sumQQbar;
candidate2 = cit;
candIsQQ2 = false;
candIsOctet2 = octetNormalCR;
maxrap2 = rapq;
minrap2 = rapqbar;
}
// choose the less stringent threshold for further iterations
PositiveRapidtyThreshold = maxrap1 > maxrap2 ? maxrap2:maxrap1;
NegativeRapidtyThreshold = minrap1 < minrap2 ? minrap2:minrap1;
foundCR=true;
}
}
assert(PositiveRapidtyThreshold<=maxrap1);
assert(PositiveRapidtyThreshold<=maxrap2);
assert(NegativeRapidtyThreshold>=minrap1);
assert(NegativeRapidtyThreshold>=minrap2);
assert(maxsum1>=maxsum2);
if ( rapq < 0.0 && rapqbar > 0.0
&& rapqbar > PositiveRapidtyThreshold/_mesonToBaryonFactor
&& rapq < NegativeRapidtyThreshold/_mesonToBaryonFactor ) {
//sum of rapidities of quarks
const double sumQQ = abs(rapq) + abs(rapqbar);
if ( sumQQ > maxsum2/_mesonToBaryonFactor ) {
if ( sumQQ > maxsum1 ) {
double factor = candIsQQ1 ? _mesonToBaryonFactor:1.0;
maxsum2 = (factor*maxsum1) > sumQQ ? sumQQ:(factor*maxsum1);
candidate2 = candidate1;
candIsQQ2 = candIsQQ1;
candIsOctet2 = candIsOctet1;
maxrap2 = maxrap1;
minrap2 = minrap1;
maxsum1 = sumQQ;
candidate1 = cit;
candIsQQ1 = true;
candIsOctet1 = octetNormalCR;
maxrap1 = rapqbar;
minrap1 = rapq;
} else {
maxsum2 = (_mesonToBaryonFactor*maxsum1) > sumQQ ? sumQQ:(_mesonToBaryonFactor*maxsum1);
candidate2 = cit;
candIsQQ2 = true;
candIsOctet2 = octetNormalCR;
maxrap2 = rapqbar;
minrap2 = rapq;
}
// choose the less stringent threshold for further iterations
PositiveRapidtyThreshold = maxrap1 > maxrap2 ? maxrap2:maxrap1;
NegativeRapidtyThreshold = minrap1 < minrap2 ? minrap2:minrap1;
foundCR=true;
}
}
assert(PositiveRapidtyThreshold<=maxrap1);
assert(PositiveRapidtyThreshold<=maxrap2);
assert(NegativeRapidtyThreshold>=minrap1);
assert(NegativeRapidtyThreshold>=minrap2);
assert(maxsum1>=maxsum2);
}
// determine the type
if (!candIsQQ1) {
// Mesonic CR
if (candIsOctet1) {
if (!candIsQQ2 && !candIsOctet2) {
// TODO Here is the problem
if (candidate2!=cl){
swap(candidate2,candidate1);
typeOfReconnection = 1;
}
else typeOfReconnection = 0;
// typeOfReconnection = 0;
}
else
typeOfReconnection = 0;
}
else
typeOfReconnection = 1;
}
else if (candIsQQ1)
{
if (candIsQQ2)
// Baryonic CR
typeOfReconnection = 2;
else
// Diquark CR
typeOfReconnection = 3;
}
if (!foundCR) typeOfReconnection = 0;
// veto reconnection if cannot make a Diquark Cluster
bool failDCR=false;
if (typeOfReconnection == 3) {
if (_diquarkCR && !_canMakeDiquarkCluster(*cl, *candidate1)) {
/* TODO Here is the problem
if (!candIsQQ2 && !candIsOctet2 && candidate2!=cl) {
// if second nearest is candidate for Mesonic CR
// allow MCR
candidate1=candidate2;
typeOfReconnection=1;
}
else if ( _canMakeDiquarkCluster(*cl, *candidate2) && candidate2!=cl) {
// if second nearest is allowed for DCR
// allow DCR
candidate1=candidate2;
typeOfReconnection=3;
}
else
{
*/
// No CR
typeOfReconnection = 0;
failDCR=true;
// }
}
}
+ if (typeOfReconnection == 2 && !_canMakeBaryonicCluster(*cl,*candidate1,*candidate2)) {
+ if (_canMakeDiquarkCluster(*cl,*candidate1)){
+ // if we cannot make baryonic CR try diquark CR
+ typeOfReconnection = 3;
+ }
+ else {
+ // reject CR if we cannot make neither baryonic nor diquark CR
+ typeOfReconnection = 0;
+ }
+ }
if (_debug) {
std::ofstream outTypes("WriteOut/TypesOfDCR.dat", std::ios::app);
outTypes << (failDCR ? 4:typeOfReconnection) << "\n";
outTypes.close();
switch (typeOfReconnection)
{
// Mesonic CR
case 1:
{
std::ofstream outMCR("WriteOut/MCR.dat", std::ios::app);
outMCR << minrap1 << "\t"
<< maxrap1 << "\t"
<< minrap2 << "\t"
<< maxrap2 << "\t"
<<"\n";
outMCR.close();
break;
}
// Baryonic CR
case 2:
{
std::ofstream outBCR("WriteOut/BCR.dat", std::ios::app);
outBCR << minrap1 << "\t"
<< maxrap1 << "\t"
<< minrap2 << "\t"
<< maxrap2 << "\t"
<<"\n";
outBCR.close();
break;
}
// Diquark CR
case 3:
{
std::ofstream outDCR("WriteOut/DCR.dat", std::ios::app);
outDCR << minrap1 << "\t"
<< maxrap1 << "\t"
<< minrap2 << "\t"
<< maxrap2 << "\t"
<<"\n";
outDCR.close();
break;
}
// No CR found
case 0:
{
std::ofstream outNoCR("WriteOut/NoCR.dat", std::ios::app);
outNoCR<< minrap1 << "\t"
<< maxrap1 << "\t"
<< minrap2 << "\t"
<< maxrap2 << "\t"
<<"\n";
outNoCR.close();
break;
}
default:
assert(false);
}
}
}
namespace {
struct ClusterInfo{
double closenessMeasure;
bool isQQ;
bool isOctetWithOriginal;
CluVecIt cluIt;
} ;
bool compare_Clusters(ClusterInfo c1, ClusterInfo c2){
return c1.closenessMeasure>c2.closenessMeasure;
}
}
void ColourReconnector::_findPartnerBaryonicDiquarkClusterTEST3(
const CluVecIt & cl, ClusterVector & cv,
unsigned & typeOfReconnection,
const ClusterVector& deleted,
CluVecIt &candidate1,
CluVecIt &candidate2 ) const {
typeOfReconnection=0; // no Reconnection found
using Constants::pi;
using Constants::twopi;
candidate1=cl;
candidate2=cl;
// boost into RF of cl
// Lorentz5Momentum cl1 = (*cl)->momentum();
Lorentz5Momentum p1anticol = (*cl)->antiColParticle()->momentum();
Lorentz5Momentum p1col = (*cl)->colParticle()->momentum();
std::vector<ClusterInfo> cluVec;
cluVec.reserve(cv.size());
for (CluVecIt cit=cv.begin(); cit != cv.end(); ++cit) {
//avoid looping over clusters containing diquarks
if ( (*cit)->numComponents()>=3 ) continue;
//skip the cluster to be deleted later 3->2 cluster
if ( find(deleted.begin(), deleted.end(), *cit) != deleted.end() )
continue;
if ( hasDiquark(*cit) ) continue;
if ( (*cl)->isBeamCluster() && (*cit)->isBeamCluster() )
continue;
if ( cit==cl ) continue;
// veto if Clusters further apart than _maxDistance
if (_localCR && ((**cl).vertex().vect()-(**cit).vertex().vect()).mag() > _maxDistance) continue;
// veto if Clusters have negative spacetime difference
if (_causalCR && ((**cl).vertex()-(**cit).vertex()).m() < ZERO) continue;
bool octetNormalCR =
(_isColour8( (*cl)->colParticle(), (*cit)->antiColParticle() )
||
_isColour8( (*cit)->colParticle(), (*cl)->antiColParticle() ) );
Lorentz5Momentum p2col = (*cit)->colParticle()->momentum();
Lorentz5Momentum p2anticol = (*cit)->antiColParticle()->momentum();
double Morig = _selectionChoice==3 ? double(((p1col+p1anticol).m() + (p2col+p2anticol).m())/GeV)
:(pow(log(p1col*p1anticol/(p1col.m()*p1anticol.m())),2) + pow(log(p2col*p2anticol/(p2col.m()*p2anticol.m())),2));
double MmesonicCR = _selectionChoice==3 ? double(((p1col+p2anticol).m() + (p2col+p1anticol).m())/GeV)
:(pow(log(p1col*p2anticol/(p1col.m()*p2anticol.m())),2) + pow(log(p2col*p1anticol/(p2col.m()*p1anticol.m())),2));
double MdiquarkCR = _selectionChoice==3 ? double(((p1col+p2col).m() + (p1anticol+p2anticol).m())/GeV)
:(pow(log(p1col*p2col/(p1col.m()*p2col.m())),2) + pow(log(p1anticol*p2anticol/(p1anticol.m()*p2anticol.m())),2));
if ( Morig < MmesonicCR && Morig < MdiquarkCR )
continue;
// configuration for Mesonic CR
if ( MdiquarkCR > MmesonicCR ) {
//sum of rapidities of quarks
ClusterInfo cluMesInfo;
cluMesInfo.closenessMeasure=1.0/MmesonicCR;
cluMesInfo.isQQ=false;
cluMesInfo.isOctetWithOriginal=octetNormalCR;
cluMesInfo.cluIt=cit;
cluVec.push_back(cluMesInfo);
}
else if ( MdiquarkCR <= MmesonicCR ) {
//sum of rapidities of quarks
ClusterInfo cluMesInfo;
cluMesInfo.closenessMeasure=1.0/MdiquarkCR;
cluMesInfo.isQQ=true;
cluMesInfo.isOctetWithOriginal=false; // never important for DCR
cluMesInfo.cluIt=cit;
cluVec.push_back(cluMesInfo);
}
else {
assert(false);
}
}
std::sort(cluVec.begin(), cluVec.end(), compare_Clusters);
switch (cluVec.size())
{
case 0:
typeOfReconnection=0;
return;
case 1:
{
if (cluVec[0].isQQ) {
// diquark CR if possible
if (_canMakeDiquarkCluster(*cl, *(cluVec[0].cluIt))) {
typeOfReconnection=3;
candidate1=cluVec[0].cluIt;
return;
}
else {
typeOfReconnection=0;
return;
}
}
else {
// Mesonic CR if not Octet
if (cluVec[0].isOctetWithOriginal){
typeOfReconnection=0;
return;
}
else {
candidate1=cluVec[0].cluIt;
typeOfReconnection=1;
return;
}
}
break;
}
}
bool candidate1isQQ=false;
// bool candidate2isQQ=false;
typeOfReconnection=0;
for (const auto & c : cluVec) {
if (candidate1!=cl && candidate2!=cl)
break;
if (c.isQQ) {
// diquark CR if possible
if (_canMakeDiquarkCluster(*cl, *(c.cluIt))) {
if (candidate1==cl){
candidate1=c.cluIt;
candidate1isQQ=true;
typeOfReconnection=3;
}
else {
candidate2=c.cluIt;
// candidate2isQQ=true;
typeOfReconnection=2;
}
}
}
else {
// Mesonic CR if not Octet
if (c.isOctetWithOriginal) {
continue;
}
else {
if (candidate1==cl){
candidate1=c.cluIt;
candidate1isQQ=false;
typeOfReconnection=1;
}
else {
candidate2=c.cluIt;
// candidate2isQQ=false;
if (candidate1isQQ)
typeOfReconnection=3;
else
typeOfReconnection=1;
break;
}
}
}
}
return;
}
void ColourReconnector::_findPartnerBaryonicDiquarkClusterTEST2(
const CluVecIt & cl, ClusterVector & cv,
unsigned & typeOfReconnection,
const ClusterVector& deleted,
CluVecIt &candidate1,
CluVecIt &candidate2 ) const {
typeOfReconnection=0; // no Reconnection found
using Constants::pi;
using Constants::twopi;
candidate1=cl;
candidate2=cl;
// boost into RF of cl
Lorentz5Momentum cl1 = (*cl)->momentum();
const Boost boostv(-cl1.boostVector());
Lorentz5Momentum p1anticol = (*cl)->antiColParticle()->momentum();
// direction
p1anticol.boost(boostv);
std::vector<ClusterInfo> cluVec;
cluVec.reserve(cv.size());
double closeness;
for (CluVecIt cit=cv.begin(); cit != cv.end(); ++cit) {
//avoid looping over clusters containing diquarks
if ( (*cit)->numComponents()>=3 ) continue;
//skip the cluster to be deleted later 3->2 cluster
if ( find(deleted.begin(), deleted.end(), *cit) != deleted.end() )
continue;
if ( hasDiquark(*cit) ) continue;
if ( (*cl)->isBeamCluster() && (*cit)->isBeamCluster() )
continue;
if ( cit==cl ) continue;
// veto if Clusters further apart than _maxDistance
if (_localCR && ((**cl).vertex().vect()-(**cit).vertex().vect()).mag() > _maxDistance) continue;
// veto if Clusters have negative spacetime difference
if (_causalCR && ((**cl).vertex()-(**cit).vertex()).m() < ZERO) continue;
bool octetNormalCR =
(_isColour8( (*cl)->colParticle(), (*cit)->antiColParticle() )
||
_isColour8( (*cit)->colParticle(), (*cl)->antiColParticle() ) );
// boost constituents of cit into RF of cl
Lorentz5Momentum p2col = (*cit)->colParticle()->momentum();
Lorentz5Momentum p2anticol = (*cit)->antiColParticle()->momentum();
p2col.boost(boostv);
p2anticol.boost(boostv);
// calculate the rapidity of the other constituents of the clusters
// w.r.t axis of p1anticol.vect.unit
const double rapq = calculateRapidityRF(p1anticol,p2col);
const double rapqbar = calculateRapidityRF(p1anticol,p2anticol);
// std::cout << "\nPRE\n"
// << "\t octet = " << octetNormalCR << "\n"
// << "\t rapq = " << rapq << "\n"
// << "\t rapqbar = " << rapqbar << "\n";
closeness = abs(rapq) + abs(rapqbar);
// configuration for Mesonic CR
if ( rapq > 0.0 && rapqbar < 0.0 ) {
//sum of rapidities of quarks
ClusterInfo cluMesInfo;
cluMesInfo.closenessMeasure=closeness;
cluMesInfo.isQQ=false;
cluMesInfo.isOctetWithOriginal=octetNormalCR;
cluMesInfo.cluIt=cit;
cluVec.push_back(cluMesInfo);
}
else if ( rapq < 0.0 && rapqbar > 0.0 ) {
//sum of rapidities of quarks
ClusterInfo cluMesInfo;
cluMesInfo.closenessMeasure=closeness;
cluMesInfo.isQQ=true;
cluMesInfo.isOctetWithOriginal=false; // never important for DCR
cluMesInfo.cluIt=cit;
cluVec.push_back(cluMesInfo);
}
}
std::sort(cluVec.begin(), cluVec.end(), compare_Clusters);
switch (cluVec.size())
{
case 0:
typeOfReconnection=0;
return;
case 1:
{
if (cluVec[0].isQQ) {
// diquark CR if possible
if (_canMakeDiquarkCluster(*cl, *(cluVec[0].cluIt))) {
typeOfReconnection=3;
candidate1=cluVec[0].cluIt;
return;
}
else {
typeOfReconnection=0;
return;
}
}
else {
// Mesonic CR if not Octet
if (cluVec[0].isOctetWithOriginal){
typeOfReconnection=0;
return;
}
else {
candidate1=cluVec[0].cluIt;
typeOfReconnection=1;
return;
}
}
break;
}
}
bool candidate1isQQ=false;
// bool candidate2isQQ=false;
typeOfReconnection=0;
for (const auto & c : cluVec) {
if (candidate1!=cl && candidate2!=cl)
break;
if (c.isQQ) {
// diquark CR if possible
if (_canMakeDiquarkCluster(*cl, *(c.cluIt))) {
if (candidate1==cl){
candidate1=c.cluIt;
candidate1isQQ=true;
typeOfReconnection=3;
}
else {
candidate2=c.cluIt;
// candidate2isQQ=true;
typeOfReconnection=2;
}
}
}
else {
// Mesonic CR if not Octet
if (c.isOctetWithOriginal) {
continue;
}
else {
if (candidate1==cl){
candidate1=c.cluIt;
candidate1isQQ=false;
typeOfReconnection=1;
}
else {
candidate2=c.cluIt;
// candidate2isQQ=false;
if (candidate1isQQ)
typeOfReconnection=3;
else
typeOfReconnection=1;
break;
}
}
}
}
return;
}
CluVecIt ColourReconnector::_findPartnerBaryonic(
const CluVecIt & cl, ClusterVector & cv,
bool & baryonicCand,
const ClusterVector& deleted,
CluVecIt &baryonic1,
CluVecIt &baryonic2 ) const {
using Constants::pi;
using Constants::twopi;
// Returns a candidate for possible reconnection
CluVecIt candidate = cl;
bool bcand = false;
double maxrap = 0.0;
double minrap = 0.0;
double maxrapNormal = 0.0;
double minrapNormal = 0.0;
double maxsumnormal = 0.0;
double maxsum = 0.0;
double secondsum = 0.0;
// boost into RF of cl
Lorentz5Momentum cl1 = (*cl)->momentum();
const Boost boostv(-cl1.boostVector());
// cl1.boost(boostv);
// boost constituents of cl into RF of cl
// Lorentz5Momentum p1col = (*cl)->colParticle()->momentum();
Lorentz5Momentum p1anticol = (*cl)->antiColParticle()->momentum();
// p1col.boost(boostv);
p1anticol.boost(boostv);
for (CluVecIt cit=cv.begin(); cit != cv.end(); ++cit) {
//avoid looping over clusters containing diquarks
if ( hasDiquark(*cit) ) continue;
if ( (*cit)->numComponents()>=3 ) continue;
if ( cit==cl ) continue;
//skip the cluster to be deleted later 3->2 cluster
if ( find(deleted.begin(), deleted.end(), *cit) != deleted.end() )
continue;
if ( (*cl)->isBeamCluster() && (*cit)->isBeamCluster() )
continue;
// veto if Clusters further apart than _maxDistance
if (_localCR && ((**cl).vertex().vect()-(**cit).vertex().vect()).mag() > _maxDistance) continue;
// veto if Clusters have negative spacetime difference
if (_causalCR && ((**cl).vertex()-(**cit).vertex()).m() < ZERO) continue;
const bool Colour8 =
_isColour8( (*cl)->colParticle(), (*cit)->antiColParticle() )
||
_isColour8( (*cit)->colParticle(), (*cl)->antiColParticle() ) ;
// boost constituents of cit into RF of cl
Lorentz5Momentum p2col = (*cit)->colParticle()->momentum();
Lorentz5Momentum p2anticol = (*cit)->antiColParticle()->momentum();
p2col.boost(boostv);
p2anticol.boost(boostv);
// calculate the rapidity of the other constituents of the clusters
// w.r.t axis of p1anticol.vect.unit
const double rapq = calculateRapidityRF(p1anticol,p2col);
const double rapqbar = calculateRapidityRF(p1anticol,p2anticol);
// configuration for normal CR
if ( !Colour8
&& rapq > 0.0 && rapqbar < 0.0
&& rapq > maxrap
&& rapqbar < minrap ) {
maxrap = rapq;
minrap = rapqbar;
//sum of rapidities of quarks
const double normalsum = abs(rapq) + abs(rapqbar);
if ( normalsum > maxsumnormal ) {
maxsumnormal = normalsum;
maxrapNormal = rapq;
minrapNormal = rapqbar;
bcand = false;
candidate = cit;
}
}
if ( rapq < 0.0 && rapqbar >0.0
&& rapqbar > maxrapNormal
&& rapq < minrapNormal ) {
maxrap = rapqbar;
minrap = rapq;
const double sumrap = abs(rapqbar) + abs(rapq);
// first candidate gets here. If second baryonic candidate has higher Ysum than the first
// one, the second candidate becomes the first one and the first the second.
if (sumrap > maxsum) {
if (maxsum != 0) {
baryonic2 = baryonic1;
baryonic1 = cit;
bcand = true;
} else {
baryonic1 = cit;
}
maxsum = sumrap;
} else {
if (sumrap > secondsum && sumrap != maxsum) {
secondsum = sumrap;
bcand = true;
baryonic2 = cit;
}
}
}
}
if (bcand == true) {
baryonicCand = true;
}
if (_debug) {
std::ofstream outTypes("WriteOut/TypesOfBCR.dat", std::ios::app);
outTypes << (baryonicCand ? 2:(candidate==cl ? 0:1)) << "\n";
outTypes.close();
}
return candidate;
}
CluVecIt ColourReconnector::_findRecoPartnerPlainDynamic(const CluVecIt & cl,
ClusterVector & cv, const ClusterVector & deleted, bool diquarkCR) const {
CluVecIt candidate = cl;
double pNoRecoMin=1.0;
double pNoReco,preco,precoDiquark;
// boost into RF of cl
Lorentz5Momentum cl1 = (*cl)->momentum();
const Boost boostv(-cl1.boostVector());
// boost constituents of cl into RF of cl
Lorentz5Momentum p1anticol = (*cl)->antiColParticle()->momentum();
p1anticol.boost(boostv);
for (CluVecIt cit=cv.begin(); cit != cv.end(); ++cit) {
// Skip deleted clusters
if (find(deleted.begin(), deleted.end(), *cit) != deleted.end())
continue;
// skip diquark clusters
if ((*cit)->numComponents()>2 || (_diquarkCR && hasDiquark(*cit))) continue;
// don't even look at original cluster
if (cit==cl) continue;
// don't allow colour octet clusters
bool Colour8 = ( _isColour8( (*cl)->colParticle(),
(*cit)->antiColParticle() ) ||
_isColour8( (*cit)->colParticle(),
(*cl)->antiColParticle() ) );
// stop it putting beam remnants together
if ((*cl)->isBeamCluster() && (*cit)->isBeamCluster()) continue;
// veto if Clusters further apart than _maxDistance
if (_localCR && ((**cl).vertex().vect()-(**cit).vertex().vect()).mag() > _maxDistance) continue;
// veto if Clusters have negative spacetime difference
if (_causalCR && ((**cl).vertex()-(**cit).vertex()).m() < ZERO) continue;
// boost constituents of cit into RF of cl
Lorentz5Momentum p2col = (*cit)->colParticle()->momentum();
Lorentz5Momentum p2anticol = (*cit)->antiColParticle()->momentum();
p2col.boost(boostv);
p2anticol.boost(boostv);
// calculate the rapidity of the other constituents of the clusters
// w.r.t axis of p1anticol.vect.unit
const double rapq = calculateRapidityRF(p1anticol,p2col);
const double rapqbar = calculateRapidityRF(p1anticol,p2anticol);
// configuration for normal CR
if ( !Colour8
&& rapq > 0.0 && rapqbar < 0.0) {
}
// configuration for diquark CR
else if ( _diquarkCR && rapq < 0.0 && rapqbar > 0.0) {
}
else
continue;
// Get dynamic CR probabilities and try to minimize the non-reco probability
std::tie( pNoReco, preco, precoDiquark) = _dynamicRecoProbabilitiesCF2(*cl,*cit,diquarkCR);
if (pNoReco<pNoRecoMin) {
candidate = cit;
pNoRecoMin=pNoReco;
}
}
if (_debug) {
ofstream out("WriteOut/pNoReco.dat", std::ios::app);
out << pNoRecoMin << "\t" << cv.size() << "\n";
out.close();
}
return candidate;
}
CluVecIt ColourReconnector::_findRecoPartnerPlain(const CluVecIt & cl,
ClusterVector & cv, const ClusterVector & deleted) const {
CluVecIt candidate = cl;
Energy minMass = 1*TeV;
for (CluVecIt cit=cv.begin(); cit != cv.end(); ++cit) {
// Skip deleted clusters
if (find(deleted.begin(), deleted.end(), *cit) != deleted.end())
continue;
// skip diquark clusters
if ((*cit)->numComponents()>2 || (_diquarkCR && hasDiquark(*cit))) continue;
// don't even look at original cluster
if (cit==cl) continue;
// don't allow colour octet clusters
if ( _isColour8( (*cl)->colParticle(),
(*cit)->antiColParticle() ) ||
_isColour8( (*cit)->colParticle(),
(*cl)->antiColParticle() ) ) {
continue;
}
// stop it putting beam remnants together
if ((*cl)->isBeamCluster() && (*cit)->isBeamCluster()) continue;
// veto if Clusters further apart than _maxDistance
if (_localCR && ((**cl).vertex().vect()-(**cit).vertex().vect()).mag() > _maxDistance) continue;
// veto if Clusters have negative spacetime difference
if (_causalCR && ((**cl).vertex()-(**cit).vertex()).m() < ZERO) continue;
// momenta of the old clusters
Lorentz5Momentum p1 = (*cl)->colParticle()->momentum() +
(*cl)->antiColParticle()->momentum();
Lorentz5Momentum p2 = (*cit)->colParticle()->momentum() +
(*cit)->antiColParticle()->momentum();
// momenta of the new clusters
Lorentz5Momentum p3 = (*cl)->colParticle()->momentum() +
(*cit)->antiColParticle()->momentum();
Lorentz5Momentum p4 = (*cit)->colParticle()->momentum() +
(*cl)->antiColParticle()->momentum();
Energy oldMass = abs( p1.m() ) + abs( p2.m() );
Energy newMass = abs( p3.m() ) + abs( p4.m() );
if ( newMass < oldMass && newMass < minMass ) {
minMass = newMass;
candidate = cit;
}
}
return candidate;
}
// forms two baryonic clusters from three clusters
void ColourReconnector::_makeBaryonicClusters(
ClusterPtr &c1, ClusterPtr &c2,
ClusterPtr &c3,
ClusterPtr &newcluster1,
ClusterPtr &newcluster2) const {
//make sure they all have 2 components
assert(c1->numComponents()==2);
assert(c2->numComponents()==2);
assert(c3->numComponents()==2);
//abandon children
c1->colParticle()->abandonChild(c1);
c1->antiColParticle()->abandonChild(c1);
c2->colParticle()->abandonChild(c2);
c2->antiColParticle()->abandonChild(c2);
c3->colParticle()->abandonChild(c3);
c3->antiColParticle()->abandonChild(c3);
newcluster1 = new_ptr(Cluster(c1->colParticle(),c2->colParticle(), c3->colParticle()));
c1->colParticle()->addChild(newcluster1);
c2->colParticle()->addChild(newcluster1);
c3->colParticle()->addChild(newcluster1);
newcluster1->setVertex(LorentzPoint());
newcluster2 = new_ptr(Cluster(c1->antiColParticle(), c2->antiColParticle(),
c3->antiColParticle()));
c1->antiColParticle()->addChild(newcluster2);
c2->antiColParticle()->addChild(newcluster2);
c3->antiColParticle()->addChild(newcluster2);
newcluster2->setVertex(LorentzPoint());
// auto q1 = c1->colParticle()->momentum();
// auto q2 = c2->colParticle()->momentum();
// auto q3 = c3->colParticle()->momentum();
// auto qbar1 = c1->antiColParticle()->momentum();
// auto qbar2 = c2->antiColParticle()->momentum();
// auto qbar3 = c3->antiColParticle()->momentum();
// double mB12 = q1*q2/(q1.m()*q2.m())-1;
// double mB23 = q2*q3/(q2.m()*q3.m())-1;
// double mB13 = q1*q3/(q1.m()*q3.m())-1;
// double mAB12 = qbar1*qbar2/(qbar1.m()*qbar2.m())-1;
// double mAB23 = qbar2*qbar3/(qbar2.m()*qbar3.m())-1;
// double mAB13 = qbar1*qbar3/(qbar1.m()*qbar3.m())-1;
// double mBmin = mB12;
// double mABmin = mAB12;
// int choiceB = 2;
// int choiceAB = 2;
// if (mBmin > mB23){
// choiceB = 0;
// mBmin = mB23;
// }
// if (mBmin > mB13) {
// choiceB = 1;
// mBmin = mB13;
// }
// if (mABmin > mAB23) {
// choiceAB = 0;
// mABmin = mAB23;
// }
// if (mABmin > mAB13) {
// choiceAB = 1;
// mABmin = mAB13;
// }
// if (mBmin>_cutDiquarkClusterFormation){
// std::cout << "mBmin = " << mBmin << std::endl;
// }
// if (mABmin>_cutDiquarkClusterFormation){
// std::cout << "mABmin = " << mABmin << std::endl;
// }
}
bool ColourReconnector::_canMakeDiquarkCluster(tcPPtr pCol1, tcPPtr pCol2,tcPPtr pAntiCol1, tcPPtr pAntiCol2) const{
double a;
return _canMakeDiquarkCluster( pCol1, pCol2, pAntiCol1, pAntiCol2,a);
}
bool ColourReconnector::_canMakeDiquarkCluster(tcPPtr pCol1, tcPPtr pCol2,tcPPtr pAntiCol1, tcPPtr pAntiCol2, double & PhaseSpace) const{
tcPDPtr dataDiquark = _hadronSpectrum->makeDiquark(pCol1->dataPtr(), pCol2->dataPtr());
tcPDPtr dataDiquarkBar = _hadronSpectrum->makeDiquark(pAntiCol1->dataPtr(), pAntiCol2->dataPtr());
if (!dataDiquark){
throw Exception() << "Could not make a diquark from"
<< pCol1->dataPtr()->PDGName() << " and "
<< pCol2->dataPtr()->PDGName()
<< " in ColourReconnector::_canMakeDiquarkCluster()"
<< Exception::eventerror;
}
if (!dataDiquarkBar){
throw Exception() << "Could not make an anti-diquark from"
<< pAntiCol1->dataPtr()->PDGName() << " and "
<< pAntiCol2->dataPtr()->PDGName()
<< " in ColourReconnector::_canMakeDiquarkCluster()"
<< Exception::eventerror;
}
-
- Energy minMassCD = _hadronSpectrum->massLightestHadronPair(dataDiquark,dataDiquarkBar);
-
+ int diqTreatment = _clusterFinder->diquarkOnShell();
Lorentz5Momentum Ptot=pCol1->momentum() + pCol2->momentum() + pAntiCol1->momentum() + pAntiCol2->momentum();
Energy Mass=Ptot.m();
+ Energy minMassCD = _hadronSpectrum->massLightestHadronPair(dataDiquark,dataDiquarkBar);
+ // We need to guarantee that a diquark cluster can decay to the two lightest hadrons
if ( Mass<=minMassCD ) {
+ // DiQuarkOnShell = Mixed (Default) -> diqTreatment==-1
+ // DiQuarkOnShell = No -> diqTreatment==0
return false;
}
- Energy minMassOnShell = dataDiquark->constituentMass() + dataDiquarkBar->constituentMass();
- if ( Mass<=minMassOnShell ) {
- return false;
- }
+
+ if (diqTreatment == 1){
+ // DiQuarkOnShell = Yes
+ Energy minMassOnShell = dataDiquark->constituentMass() + dataDiquarkBar->constituentMass();
+
+ // We need to guarantee that a diquark cluster can have its two future diquarks on shell
+ if ( Mass<=minMassOnShell ) {
+ return false;
+ }
+ }
+
if (_cutDiquarkClusterFormation>0) {
double cut = _cutDiquarkClusterFormation;
double valDiq = pCol1->momentum()*pCol2->momentum()/(pCol1->momentum().m()*pCol2->momentum().m())-1.0;
double valAntiDiq = pAntiCol1->momentum()*pAntiCol2->momentum()/(pAntiCol1->momentum().m()*pAntiCol2->momentum().m())-1.0;
if (valDiq>cut)
return false;
if (valAntiDiq>cut)
return false;
}
if (_phaseSpaceDiquarkFission){
// Tried to add a factor suppressing to low mass Diquark Clusters
double factor;
switch (_phaseSpaceDiquarkFission)
{
case 1:
{
const auto & BaryonPair=_hadronSpectrum->lightestHadronPair(dataDiquark,dataDiquarkBar);
if (Mass-(BaryonPair.first->mass()+BaryonPair.second->mass())<=ZERO)
return false;
factor = 2.0*Kinematics::pstarTwoBodyDecay(Mass,BaryonPair.first->mass(),BaryonPair.second->mass())/Mass;
break;
}
case 2:
{
if (Mass-(dataDiquark->constituentMass()+dataDiquarkBar->constituentMass())<ZERO)
return false;
factor = 2.0*Kinematics::pstarTwoBodyDecay(Mass,dataDiquark->constituentMass(),dataDiquarkBar->constituentMass())/Mass;
break;
}
default:
assert(false);
}
assert(factor>=0.0);
assert(factor<=1.0);
PhaseSpace = factor;
}
else {
PhaseSpace = 1.0;
}
return true;
}
+bool ColourReconnector::_canMakeBaryonicCluster(const ClusterPtr &c1, const ClusterPtr &c2, const ClusterPtr &c3) const {
+ //make sure they all have 2 components
+ assert(c1->numComponents()==2);
+ assert(c2->numComponents()==2);
+ assert(c3->numComponents()==2);
+ vector<tcPPtr> Baryon = {c1->colParticle(), c2->colParticle(), c3->colParticle()};
+ vector<tcPPtr> AntiBaryon = {c1->antiColParticle(), c2->antiColParticle(), c3->antiColParticle()};
+ return (_canMakeBaryonicCluster(Baryon)
+ && _canMakeBaryonicCluster(AntiBaryon));
+}
+
+bool ColourReconnector::_canMakeBaryonicCluster(vector<tcPPtr> pCol) const {
+ tcPDPtr dataDiquark12 = _hadronSpectrum->makeDiquark(pCol[0]->dataPtr(), pCol[1]->dataPtr());
+ tcPDPtr dataDiquark13 = _hadronSpectrum->makeDiquark(pCol[0]->dataPtr(), pCol[2]->dataPtr());
+ tcPDPtr dataDiquark23 = _hadronSpectrum->makeDiquark(pCol[1]->dataPtr(), pCol[2]->dataPtr());
+
+ if (!dataDiquark12){
+ throw Exception() << "Could not make a diquark from"
+ << pCol[0]->dataPtr()->PDGName() << " and "
+ << pCol[1]->dataPtr()->PDGName()
+ << " in ColourReconnector::_canMakeBaryonicCluster()"
+ << Exception::eventerror;
+ }
+ if (!dataDiquark13){
+ throw Exception() << "Could not make a diquark from"
+ << pCol[0]->dataPtr()->PDGName() << " and "
+ << pCol[2]->dataPtr()->PDGName()
+ << " in ColourReconnector::_canMakeBaryonicCluster()"
+ << Exception::eventerror;
+ }
+ if (!dataDiquark23){
+ throw Exception() << "Could not make a diquark from"
+ << pCol[1]->dataPtr()->PDGName() << " and "
+ << pCol[2]->dataPtr()->PDGName()
+ << " in ColourReconnector::_canMakeBaryonicCluster()"
+ << Exception::eventerror;
+ }
+ int diqTreatment = _clusterFinder->diquarkOnShell();
+ Lorentz5Momentum Ptot=pCol[0]->momentum() + pCol[1]->momentum() + pCol[2]->momentum();
+ Energy Mass=Ptot.m();
+ Energy minMassCD = _hadronSpectrum->massLightestHadron(dataDiquark12, pCol[2]->dataPtr());
+ if ( Mass<=minMassCD ) {
+ // DiQuarkOnShell = Mixed (Default) -> diqTreatment==-1
+ // DiQuarkOnShell = No -> diqTreatment==0
+ return false;
+ }
+
+ if (diqTreatment == 1){
+ // DiQuarkOnShell = Yes
+ Energy minMassOnShell12 = dataDiquark12->constituentMass() + pCol[2]->dataPtr()->constituentMass();
+ Energy minMassOnShell13 = dataDiquark13->constituentMass() + pCol[1]->dataPtr()->constituentMass();
+ Energy minMassOnShell23 = dataDiquark23->constituentMass() + pCol[0]->dataPtr()->constituentMass();
+
+ if ( Mass<=minMassOnShell12 || Mass<=minMassOnShell13 || Mass<=minMassOnShell23 ) {
+ return false;
+ }
+ }
+ return true;
+}
+
+
+
bool ColourReconnector::_canMakeDiquarkCluster(const ClusterPtr &c1, const ClusterPtr &c2) const {
double a;
return _canMakeDiquarkCluster(c1,c2,a);
}
// forms a four-quark cluster
bool ColourReconnector::_canMakeDiquarkCluster(const ClusterPtr &c1, const ClusterPtr &c2, double & PhaseSpace) const {
//make sure they all have 2 components
assert(c1->numComponents()==2);
assert(c2->numComponents()==2);
-
-
-
return _canMakeDiquarkCluster(c1->colParticle(),c2->colParticle(),c1->antiColParticle(),c2->antiColParticle(),PhaseSpace);
}
bool ColourReconnector::_makeDiquarkCluster(
ClusterPtr &c1, ClusterPtr &c2,
ClusterPtr &newcluster) const{
// std::cout << "MakeDiq" << std::endl;
if (!_canMakeDiquarkCluster(c1,c2)){
std::cout << "Should never execute! check this earlier" << std::endl;
return false;
}
//abandon children
c1->colParticle()->abandonChild(c1);
c1->antiColParticle()->abandonChild(c1);
c2->colParticle()->abandonChild(c2);
c2->antiColParticle()->abandonChild(c2);
// Note: the convention that c1->newcluster.particle(0,2) and c2->newcluster.particle(1,3)
// TODO probably need to use particleB to access it
newcluster = new_ptr(Cluster(c1->colParticle(),c2->colParticle(),
c1->antiColParticle(), c2->antiColParticle()));
c1->colParticle()->addChild(newcluster);
c2->colParticle()->addChild(newcluster);
c1->antiColParticle()->addChild(newcluster);
c2->antiColParticle()->addChild(newcluster);
newcluster->setVertex(LorentzPoint());
return true;
}
pair <ClusterPtr,ClusterPtr> ColourReconnector::_splitDiquarkCluster(
ClusterPtr &diquarkCluster, bool colourReconnect) const{
// Works just fine!
// std::cout << "MakeDiq" << std::endl;
// if (!_canMakeDiquarkCluster(c1,c2)){
// std::cout << "Should never execute! check this earlier" << std::endl;
// return false;
// }
ClusterPtr newcluster1,newcluster2;
//abandon children
diquarkCluster->particleB(0)->abandonChild(diquarkCluster);
diquarkCluster->particleB(1)->abandonChild(diquarkCluster);
diquarkCluster->particleB(2)->abandonChild(diquarkCluster);
diquarkCluster->particleB(3)->abandonChild(diquarkCluster);
// Note: the convention that c1->newcluster.particle(0,2) and c2->newcluster.particle(1,3)
// TODO probably need to use particleB to access it
// Here we decide if we want to colour reconnect the original clusters (0,2) and (1,3) to
// the clusters (0,3) and (1,2)
if (colourReconnect) {
// colour reconnect
newcluster1 = new_ptr(Cluster(diquarkCluster->particleB(0), diquarkCluster->particleB(3)));
newcluster2 = new_ptr(Cluster(diquarkCluster->particleB(1), diquarkCluster->particleB(2)));
diquarkCluster->particleB(0)->addChild(newcluster1);
diquarkCluster->particleB(1)->addChild(newcluster2);
diquarkCluster->particleB(2)->addChild(newcluster2);
diquarkCluster->particleB(3)->addChild(newcluster1);
}
else {
// keep original colour connection
newcluster1 = new_ptr(Cluster(diquarkCluster->particleB(0), diquarkCluster->particleB(2)));
newcluster2 = new_ptr(Cluster(diquarkCluster->particleB(1), diquarkCluster->particleB(3)));
diquarkCluster->particleB(0)->addChild(newcluster1);
diquarkCluster->particleB(1)->addChild(newcluster2);
diquarkCluster->particleB(2)->addChild(newcluster1);
diquarkCluster->particleB(3)->addChild(newcluster2);
}
newcluster1->setVertex(LorentzPoint());
newcluster2->setVertex(LorentzPoint());
return pair <ClusterPtr, ClusterPtr> (newcluster1, newcluster2);
}
pair <ClusterPtr,ClusterPtr>
ColourReconnector::_reconnect2Bto2B(ClusterPtr &c1, ClusterPtr &c2, const int s1, const int s2) const {
// form the first new cluster
// separate the quarks from their original cluster
c1->particleB((s1+1)%3)->abandonChild(c1);
c1->particleB((s1+2)%3)->abandonChild(c1);
c2->particleB(s2)->abandonChild(c2);
// now the new cluster
ClusterPtr newCluster1 = new_ptr(Cluster(c1->particleB((s1+1)%3), c1->particleB((s1+2)%3), c2->particleB(s2)));
c1->particleB((s1+1)%3)->addChild(newCluster1);
c1->particleB((s1+2)%3)->addChild(newCluster1);
c2->particleB(s2)->addChild(newCluster1);
// set new vertex
newCluster1->setVertex(LorentzPoint());
// set beam remnants for new cluster
if (c1->isBeamRemnant((s1+1)%3)) newCluster1->setBeamRemnant(0, true);
if (c1->isBeamRemnant((s1+2)%3)) newCluster1->setBeamRemnant(1, true);
if (c2->isBeamRemnant(s2)) newCluster1->setBeamRemnant(2, true);
// for the second cluster same procedure
c2->particleB((s2+1)%3)->abandonChild(c2);
c2->particleB((s2+2)%3)->abandonChild(c2);
c1->particleB(s1)->abandonChild(c1);
ClusterPtr newCluster2 = new_ptr(Cluster(c2->particleB((s2+1)%3), c2->particleB((s2+2)%3), c1->particleB(s1)));
c2->particleB((s2+1)%3)->addChild(newCluster2);
c2->particleB((s2+2)%3)->addChild(newCluster2);
c1->particleB(s1)->addChild(newCluster2);
newCluster2->setVertex(LorentzPoint());
if (c2->isBeamRemnant((s2+1)%3)) newCluster2->setBeamRemnant(0, true);
if (c2->isBeamRemnant((s2+2)%3)) newCluster2->setBeamRemnant(1, true);
if (c1->isBeamRemnant(s1)) newCluster2->setBeamRemnant(2, true);
return pair <ClusterPtr, ClusterPtr> (newCluster1, newCluster2);
}
std::tuple <ClusterPtr, ClusterPtr, ClusterPtr>
ColourReconnector::_reconnectBbarBto3M(ClusterPtr & c1, ClusterPtr & c2, const int s0, const int s1, const int s2) const {
// make sure they all have 3 components
assert(c1->numComponents()==3);
assert(c2->numComponents()==3);
// first Cluster
c1->particleB(0)->abandonChild(c1);
c2->particleB(s0)->abandonChild(c2);
ClusterPtr newCluster1 = new_ptr(Cluster(c1->particleB(0), c2->particleB(s0)));
c1->particleB(0)->addChild(newCluster1);
c2->particleB(s0)->addChild(newCluster1);
// set new vertex
newCluster1->setVertex(0.5*(c1->particleB(0)->vertex() + c2->particleB(s0)->vertex()));
// set beam remnants for new cluster
if (c1->isBeamRemnant(0)) newCluster1->setBeamRemnant(0, true);
if (c2->isBeamRemnant(s0)) newCluster1->setBeamRemnant(1, true);
// same for second cluster
c1->particleB(1)->abandonChild(c1);
c2->particleB(s1)->abandonChild(c2);
ClusterPtr newCluster2 = new_ptr(Cluster(c1->particleB(1), c2->particleB(s1)));
c1->particleB(1)->addChild(newCluster2);
c2->particleB(s1)->addChild(newCluster2);
newCluster2->setVertex(0.5*(c1->particleB(1)->vertex() + c2->particleB(s1)->vertex()));
if (c1->isBeamRemnant(1)) newCluster2->setBeamRemnant(0, true);
if (c2->isBeamRemnant(s1)) newCluster2->setBeamRemnant(1, true);
// same for third cluster
c1->particleB(2)->abandonChild(c1);
c2->particleB(s2)->abandonChild(c2);
ClusterPtr newCluster3 = new_ptr(Cluster(c1->particleB(2), c2->particleB(s2)));
c1->particleB(2)->addChild(newCluster3);
c2->particleB(s2)->addChild(newCluster3);
newCluster3->setVertex(0.5*(c1->particleB(2)->vertex() + c2->particleB(s2)->vertex()));
if (c1->isBeamRemnant(2)) newCluster3->setBeamRemnant(0, true);
if (c2->isBeamRemnant(s2)) newCluster3->setBeamRemnant(1, true);
return std::tuple <ClusterPtr, ClusterPtr, ClusterPtr> (newCluster1, newCluster2, newCluster3);
}
pair <ClusterPtr,ClusterPtr>
ColourReconnector::_reconnect(ClusterPtr &c1, ClusterPtr &c2) const {
// if (_isColour8(c1->colParticle(), c2->antiColParticle())
// || _isColour8(c2->colParticle(), c1->antiColParticle())) {
// std::cout << "reconnect _isColour8" << std::endl;
// }
if (_becomesColour8Cluster(c1,c2)){
std::cout << "Should never execute! _isColour8 in reconnect check this earlier" << std::endl;
}
// choose the other possibility to form two clusters from the given
// constituents
assert(c1->numComponents()==2);
assert(c2->numComponents()==2);
int c1_col(-1),c1_anti(-1),c2_col(-1),c2_anti(-1);
for(unsigned int ix=0; ix<2; ++ix) {
if (c1->particle(ix)->hasColour(false)) c1_col = ix;
else if(c1->particle(ix)->hasColour(true )) c1_anti = ix;
if (c2->particle(ix)->hasColour(false)) c2_col = ix;
else if(c2->particle(ix)->hasColour(true )) c2_anti = ix;
}
assert(c1_col>=0&&c2_col>=0&&c1_anti>=0&&c2_anti>=0);
c1->colParticle()->abandonChild(c1);
c2->antiColParticle()->abandonChild(c2);
ClusterPtr newCluster1
= new_ptr( Cluster( c1->colParticle(), c2->antiColParticle() ) );
c1->colParticle()->addChild(newCluster1);
c2->antiColParticle()->addChild(newCluster1);
/*
* TODO: Questionable setting of the vertex
* */
newCluster1->setVertex(0.5*(c1->colParticle()->vertex() +
c2->antiColParticle()->vertex()));
if(c1->isBeamRemnant(c1_col )) newCluster1->setBeamRemnant(0,true);
if(c2->isBeamRemnant(c2_anti)) newCluster1->setBeamRemnant(1,true);
c1->antiColParticle()->abandonChild(c1);
c2->colParticle()->abandonChild(c2);
ClusterPtr newCluster2
= new_ptr( Cluster( c2->colParticle(), c1->antiColParticle() ) );
c1->antiColParticle()->addChild(newCluster2);
c2->colParticle()->addChild(newCluster2);
/*
* TODO: Questionable setting of the vertex
* */
newCluster2->setVertex(0.5*(c2->colParticle()->vertex() +
c1->antiColParticle()->vertex()));
if(c2->isBeamRemnant(c2_col )) newCluster2->setBeamRemnant(0,true);
if(c1->isBeamRemnant(c1_anti)) newCluster2->setBeamRemnant(1,true);
return pair <ClusterPtr,ClusterPtr> (newCluster1, newCluster2);
}
std::tuple <ClusterPtr, ClusterPtr> ColourReconnector::_reconnect3MtoMD(ClusterVector & cluvec, const int topology) const {
// std::cout << "MakeDiq3MtoMD" << std::endl;
assert(cluvec.size()==3);
assert(topology>0 && topology<=33);
int colIndexMCR = (topology/10)-1;
int antiColIndexMCR = (topology-(colIndexMCR+1)*10)-1;
assert(colIndexMCR<3 && colIndexMCR>=0);
assert(antiColIndexMCR<3 && antiColIndexMCR>=0);
// std::cout << "topo = "<< topology <<"\t==\t" << colIndexMCR+1<<antiColIndexMCR+1 << std::endl;
/*
if (colIndexMCR==antiColIndexMCR) {
ClusterPtr DiquarkClu;
assert((colIndexMCR+1)%3!=(colIndexMCR+2)%3);
assert((colIndexMCR+1)%3!=colIndexMCR);
assert((colIndexMCR+2)%3!=colIndexMCR);
if (!_makeDiquarkCluster(cluvec[(colIndexMCR+1)%3], cluvec[(colIndexMCR+2)%3], DiquarkClu)){
std::cout << "could not make diquark cluster of index " << (colIndexMCR+1)%3 <<", " <<(colIndexMCR+2)%3<<"\nColIdx "<<colIndexMCR << std::endl;
assert(false);
// return std::tuple <ClusterPtr, ClusterPtr, ClusterPtr> (cluvec[0],cluvec[1],cluvec[2]);
}
// return (original untouched mesonic cluster, Diquark cluster, to be deleted cluster)
std::cout << "simple DCR" << std::endl;
// TODO check that the ordering is still maintained
return std::tuple <ClusterPtr, ClusterPtr> (cluvec[colIndexMCR], DiquarkClu);
}
*/
ClusterPtr cMCR1 = cluvec[colIndexMCR];
ClusterPtr cMCR2 = cluvec[antiColIndexMCR];
if (_isColour8(cMCR1->colParticle(),cMCR2->antiColParticle())){
std::cout << "Should never execute! _isColour8 in _reconnect3MtoMD check this earlier" << std::endl;
}
// std::cout << "Never " << std::endl;
// ClusterVector cluvec = cluvec;
ClusterVector newclusters;
int colIdx[3]={-1,-1,-1};
int anticolIdx[3]={-1,-1,-1};
for (int i = 0; i < 3; i++) {
assert(cluvec[i]->numComponents()==2);
for (unsigned ip= 0; ip < 2; ip++){
if (cluvec[i]->particle(ip)->hasColour(false)) colIdx[i] = ip;
if (cluvec[i]->particle(ip)->hasColour(true)) anticolIdx[i] = ip;
}
assert(colIdx[i]>=0);
assert(anticolIdx[i]>=0);
// abbandon all children
cluvec[i]->colParticle()->abandonChild(cluvec[i]);
cluvec[i]->antiColParticle()->abandonChild(cluvec[i]);
}
// make the MCR cluster with indices colIndexMCR,antiColIndexMCR
// form new mesonic cluster
ClusterPtr newClusterMCR = new_ptr(Cluster(cluvec[colIndexMCR]->colParticle(), cluvec[antiColIndexMCR]->antiColParticle()));
newClusterMCR->colParticle()->addChild(newClusterMCR);
newClusterMCR->antiColParticle()->addChild(newClusterMCR);
// set new vertex
newClusterMCR->setVertex(0.5*(newClusterMCR->colParticle()->vertex() +
newClusterMCR->antiColParticle()->vertex()));
// set beam remnants for new cluster
if (cluvec[colIndexMCR]->isBeamRemnant(colIdx[colIndexMCR])) newClusterMCR->setBeamRemnant(0, true);
if (cluvec[antiColIndexMCR]->isBeamRemnant(anticolIdx[antiColIndexMCR])) newClusterMCR->setBeamRemnant(1, true);
// make the DCR cluster with remaining indices
int indexDCRcol1 = (colIndexMCR+1)%3;
int indexDCRcol2 = (colIndexMCR+2)%3;
int indexDCRacol1 = (antiColIndexMCR+1)%3;
int indexDCRacol2 = (antiColIndexMCR+2)%3;
assert(indexDCRcol1!=indexDCRcol2);
assert(indexDCRcol1!=colIndexMCR);
assert(indexDCRcol2!=colIndexMCR);
assert(indexDCRacol1!=indexDCRacol2);
assert(indexDCRacol1!=antiColIndexMCR);
assert(indexDCRacol2!=antiColIndexMCR);
ClusterPtr newClusterDCR = new_ptr(Cluster(
cluvec[indexDCRcol1]->colParticle(),
cluvec[indexDCRcol2]->colParticle(),
cluvec[indexDCRacol1]->antiColParticle(),
cluvec[indexDCRacol2]->antiColParticle()
));
cluvec[indexDCRcol1]->colParticle()->addChild(newClusterDCR);
cluvec[indexDCRcol2]->colParticle()->addChild(newClusterDCR);
cluvec[indexDCRacol1]->antiColParticle()->addChild(newClusterDCR);
cluvec[indexDCRacol2]->antiColParticle()->addChild(newClusterDCR);
// set new vertex
newClusterDCR->setVertex(0.25*(
cluvec[indexDCRcol1]->colParticle()->vertex() +
cluvec[indexDCRcol2]->colParticle()->vertex() +
cluvec[indexDCRacol1]->antiColParticle()->vertex() +
cluvec[indexDCRacol2]->antiColParticle()->vertex()
));
// set beam remnants for new cluster
if (cluvec[indexDCRcol1]->isBeamRemnant(colIdx[indexDCRcol1])) newClusterDCR->setBeamRemnant(0, true);
if (cluvec[indexDCRcol2]->isBeamRemnant(colIdx[indexDCRcol2])) newClusterDCR->setBeamRemnant(1, true);
if (cluvec[indexDCRacol1]->isBeamRemnant(anticolIdx[indexDCRacol1])) newClusterDCR->setBeamRemnant(2, true);
if (cluvec[indexDCRacol2]->isBeamRemnant(anticolIdx[indexDCRacol2])) newClusterDCR->setBeamRemnant(3, true);
return std::tuple <ClusterPtr, ClusterPtr> (newClusterMCR,newClusterDCR);
}
std::tuple <ClusterPtr, ClusterPtr, ClusterPtr>
ColourReconnector::_reconnect3Mto3M(ClusterPtr & c1, ClusterPtr & c2, ClusterPtr & c3, const int infos [3]) const {
// check if mesonic clusters
assert(c1->numComponents()==2);
assert(c2->numComponents()==2);
assert(c3->numComponents()==2);
ClusterVector oldclusters = {c1, c2, c3};
ClusterVector newclusters;
for (int i=0; i<3; i++) {
if ( i!=infos[i] && _isColour8(oldclusters[i]->colParticle(),oldclusters[infos[i]]->antiColParticle())){
std::cout << "Should never execute! _isColour8 in _reconnect3Mto3M "<< i <<" "<< infos[i]<<" check this earlier" << std::endl;
}
int c1_col=-1;
int c2_anticol=-1;
// get which index is coloured and which anticolour
for (unsigned int ix=0; ix<2; ++ix) {
if (oldclusters[i]->particle(ix)->hasColour(false)) c1_col = ix;
if (oldclusters[infos[i]]->particle(ix)->hasColour(true)) c2_anticol = ix;
}
assert(c1_col>=0);
assert(c2_anticol>=0);
oldclusters[i]->colParticle()->abandonChild(oldclusters[i]);
oldclusters[infos[i]]->antiColParticle()->abandonChild(oldclusters[infos[i]]);
// form new cluster
ClusterPtr newCluster = new_ptr(Cluster(oldclusters[i]->colParticle(), oldclusters[infos[i]]->antiColParticle()));
oldclusters[i]->colParticle()->addChild(newCluster);
oldclusters[infos[i]]->antiColParticle()->addChild(newCluster);
// set new vertex
newCluster->setVertex(0.5*(oldclusters[i]->colParticle()->vertex() +
oldclusters[infos[i]]->antiColParticle()->vertex()));
// set beam remnants for new cluster
if (oldclusters[i]->isBeamRemnant(c1_col)) newCluster->setBeamRemnant(0, true);
if (oldclusters[infos[i]]->isBeamRemnant(c2_anticol)) newCluster->setBeamRemnant(1, true);
newclusters.push_back(newCluster);
}
return std::tuple <ClusterPtr, ClusterPtr, ClusterPtr> (newclusters[0], newclusters[1], newclusters[2]);
}
pair <ClusterPtr, ClusterPtr>
ColourReconnector::_reconnectMBtoMB(ClusterPtr & c1, ClusterPtr & c2, const int s0) const {
// make c1 the mesonic cluster
if (c1->numComponents()==2) {
assert(c2->numComponents()==3);
} else {
return _reconnectMBtoMB(c2,c1,s0);
}
int c1_col=-1;
int c1_anti=-1;
// get which index is coloured and which anticolour
for (unsigned int ix=0; ix<2; ++ix) {
if (c1->particle(ix)->hasColour(false)) c1_col = ix;
else if (c1->particle(ix)->hasColour(true)) c1_anti = ix;
}
assert(c1_col>=0);
assert(c1_anti>=0);
// pointers for the new clusters
ClusterPtr newCluster1;
ClusterPtr newCluster2;
if (c2->particle(0)->hasColour()==true) {
// first case: we have a baryonic clusters
// first make the new mesonic cluster
c1->antiColParticle()->abandonChild(c1);
c2->particleB(s0)->abandonChild(c2);
newCluster1 = new_ptr(Cluster(c1->antiColParticle(), c2->particleB(s0)));
c1->antiColParticle()->addChild(newCluster1);
c2->particleB(s0)->addChild(newCluster1);
// set new vertex
newCluster1->setVertex(0.5*(c1->antiColParticle()->vertex() +
c2->particleB(s0)->vertex()));
// set beam remnants for new cluster
if (c1->isBeamRemnant(c1_anti)) newCluster1->setBeamRemnant(0, true);
if (c2->isBeamRemnant(s0)) newCluster1->setBeamRemnant(1, true);
// then the baryonic one
c1->colParticle()->abandonChild(c1);
c2->particleB((s0+1)%3)->abandonChild(c2);
c2->particleB((s0+2)%3)->abandonChild(c2);
newCluster2 = new_ptr(Cluster(c1->colParticle(), c2->particleB((s0+1)%3), c2->particleB((s0+2)%3)));
c1->colParticle()->addChild(newCluster2);
c2->particleB((s0+1)%3)->addChild(newCluster2);
c2->particleB((s0+2)%3)->addChild(newCluster2);
// set new vertex
newCluster2->setVertex(LorentzPoint());
} else {
// second case we have an antibaryonic cluster
// first make the new mesonic cluster
c1->colParticle()->abandonChild(c1);
c2->particleB(s0)->abandonChild(c2);
newCluster1 = new_ptr(Cluster(c1->colParticle(), c2->particleB(s0)));
c1->colParticle()->addChild(newCluster1);
c2->particleB(s0)->addChild(newCluster1);
// set new vertex
newCluster1->setVertex(0.5*(c1->colParticle()->vertex() +
c2->particleB(s0)->vertex()));
// set beam remnants for new cluster
if (c1->isBeamRemnant(c1_col)) newCluster1->setBeamRemnant(0, true);
if (c2->isBeamRemnant(s0)) newCluster1->setBeamRemnant(1, true);
// then the baryonic one
c1->antiColParticle()->abandonChild(c1);
c2->particleB((s0+1)%3)->abandonChild(c2);
c2->particleB((s0+2)%3)->abandonChild(c2);
newCluster2 = new_ptr(Cluster(c1->antiColParticle(), c2->particleB((s0+1)%3), c2->particleB((s0+2)%3)));
c1->antiColParticle()->addChild(newCluster2);
c2->particleB((s0+1)%3)->addChild(newCluster2);
c2->particleB((s0+2)%3)->addChild(newCluster2);
// set new vertex
newCluster2->setVertex(LorentzPoint());
}
return pair <ClusterPtr, ClusterPtr> (newCluster1, newCluster2);
}
void ColourReconnector::_2Bto2BreconnectionFinder(ClusterPtr & c1, ClusterPtr & c2,
int & bswap1, int & bswap2, double min_displ_sum, string & kind_of_reco) const {
double tmp_delta;
for (int i=0; i<3; i++) {
for (int j=0; j<3; j++) {
// try swapping particle i of c1 with particle j of c2
tmp_delta = _displacementBaryonic(c2->particle(j), c1->particle((i+1)%3), c1->particle((i+2)%3));
tmp_delta += _displacementBaryonic(c1->particle(i), c2->particle((j+1)%3), c2->particle((j+2)%3));
if (tmp_delta < min_displ_sum) {
// if minimal displacement select the 2Bto2B CR option
min_displ_sum = tmp_delta;
bswap1 = i;
bswap2 = j;
kind_of_reco = "2Bto2B";
}
}
}
}
void ColourReconnector::_BbarBto3MreconnectionFinder(ClusterPtr & c1, ClusterPtr & c2, int & mswap0, int & mswap1, int & mswap2,
double min_displ_sum, string & kind_of_reco) const {
double pre_tmp_delta;
double tmp_delta;
for (int p1=0; p1 <3; p1++) {
// make sure not to form a mesonic octet
if (_isColour8(c1->particle(0), c2->particle(p1))) continue;
pre_tmp_delta = _displacement(c1->particle(0), c2->particle(p1));
for (int p2=1; p2<3; p2++) {
// make sure not to form a mesonic octet
if (_isColour8(c1->particle(1), c2->particle((p1+p2)%3))) continue;
if (_isColour8(c1->particle(2), c2->particle(3-p1-((p1+p2)%3)))) continue;
tmp_delta = pre_tmp_delta + _displacement(c1->particle(1), c2->particle((p1+p2)%3));
tmp_delta += _displacement(c1->particle(2), c2->particle(3-p1-((p1+p2)%3)));
// factor _mesonToBaryonFactor to compare Baryonic an mesonic cluster
tmp_delta *=_mesonToBaryonFactor;
if (tmp_delta < min_displ_sum) {
// if minimal displacement select the 2Bto3M CR option
min_displ_sum = tmp_delta;
mswap0 = p1;
mswap1 = (p1+p2)%3;
mswap2 = 3-p1-((p1+p2)%3);
kind_of_reco = "2Bto3M";
}
}
}
}
void ColourReconnector::_BMtoBMreconnectionfinder(ClusterPtr & c1, ClusterPtr & c2, int & swap, double min_displ_sum,
string & kind_of_reco) const {
assert(c1->numComponents()==2);
assert(c2->numComponents()==3);
double tmp_displ = 0;
for (int i=0; i<3; i++) {
// Differ if the second cluster is baryonic or antibaryonic
if (c2->particle(0)->hasColour()) {
// c2 is baryonic
// veto mesonic octets
if (_isColour8(c2->particle(i), c1->antiColParticle())) continue;
// factor _mesonToBaryonFactor to compare Baryonic an mesonic cluster
tmp_displ = _mesonToBaryonFactor * _displacement(c2->particle(i), c1->antiColParticle());
tmp_displ += _displacementBaryonic(c1->colParticle(), c2->particle((i+1)%3), c2->particle((i+2)%3));
} else {
// c2 is antibaryonic
// veto mesonic octets
if (_isColour8(c2->particle(i), c1->colParticle())) continue;
// factor _mesonToBaryonFactor to compare Baryonic an mesonic cluster
tmp_displ = _mesonToBaryonFactor * _displacement(c2->particle(i), c1->colParticle());
tmp_displ *= _displacementBaryonic(c1->antiColParticle(), c2->particle((i+1)%3), c2->particle((i+2)%3));
}
if (tmp_displ < min_displ_sum) {
// if minimal displacement select the MBtoMB CR option
min_displ_sum = tmp_displ;
swap = i;
kind_of_reco = "MBtoMB";
}
}
return;
}
void ColourReconnector::_3MtoXreconnectionfinder(std::vector<CluVecIt> & cv, int & swap0, int & swap1,
int & swap2, double min_displ_sum, string & kind_of_reco) const {
// case of 3M->BbarB CR
double _tmp_displ;
_tmp_displ = _displacementBaryonic((*cv[0])->colParticle(), (*cv[1])->colParticle(), (*cv[2])->colParticle());
_tmp_displ += _displacementBaryonic((*cv[0])->antiColParticle(), (*cv[1])->antiColParticle(), (*cv[2])->antiColParticle());
if (_tmp_displ < min_displ_sum) {
// if minimal displacement select the 3Mto2B CR option
kind_of_reco = "3Mto2B";
min_displ_sum = _tmp_displ;
}
// case for 3M->3M CR
/**
* if 3Mto3M reco probability (_preco3M_3M) is 0 we skip this loop
* since no 3Mto3M CR shall be performed
*/
int i,j;
int i1,i2,i3;
for (i = 0; _preco3M_3M && i<3; i++) {
// veto mesonic octets
if (_isColour8((*cv[0])->colParticle(), (*cv[i])->antiColParticle())) continue;
// factor _mesonToBaryonFactor to compare baryonic an mesonic cluster
_tmp_displ = _mesonToBaryonFactor * _displacement((*cv[0])->colParticle(), (*cv[i])->antiColParticle());
for (j=1; j<3; j++) {
// i1, i2, i3 are pairwise distinct
i1=i;
i2=((j+i)%3);
if (i1==0 && i2==1) continue;
i3=(3-i-((j+i)%3));
// veto mesonic octets
if (_isColour8((*cv[1])->colParticle(), (*cv[i2])->antiColParticle())) continue;
if (_isColour8((*cv[2])->colParticle(), (*cv[i3])->antiColParticle())) continue;
_tmp_displ += _mesonToBaryonFactor * _displacement((*cv[1])->colParticle(), (*cv[i2])->antiColParticle());
_tmp_displ += _mesonToBaryonFactor * _displacement((*cv[2])->colParticle(), (*cv[i3])->antiColParticle());
if (_tmp_displ < min_displ_sum) {
// if minimal displacement select the 3Mto3M CR option
kind_of_reco = "3Mto3M";
min_displ_sum = _tmp_displ;
swap0 = i1;
swap1 = i2;
swap2 = i3;
}
}
}
}
pair <int,int> ColourReconnector::_shuffle
(const PVector & q, const PVector & aq, unsigned maxtries) const {
const size_t nclusters = q.size();
assert (nclusters > 1);
assert (aq.size() == nclusters);
int i, j;
unsigned tries = 0;
bool octet=false;
do {
// find two different random integers in the range [0, nclusters)
i = UseRandom::irnd( nclusters );
do {
j = UseRandom::irnd( nclusters );
} while (i == j);
// check if one of the two potential clusters would be a colour octet state
octet = _isColour8( q[i], aq[j] ) || _isColour8( q[j], aq[i] ) ;
tries++;
} while (octet && tries < maxtries);
if (octet) i = j = -1;
return make_pair(i,j);
}
bool ColourReconnector::_becomesColour8Cluster(const ClusterPtr & c1, const ClusterPtr & c2) const {
return _isColour8(c1->colParticle(),c2->antiColParticle()) || _isColour8(c2->colParticle(),c1->antiColParticle());
}
bool ColourReconnector::_isColour8(tcPPtr p, tcPPtr q) const {
bool octet = false;
if(_octetOption<0) return octet;
// make sure we have a triplet and an anti-triplet
if ( ( p->hasColour() && q->hasAntiColour() ) ||
( p->hasAntiColour() && q->hasColour() ) ) {
// true if p and q are originated from a colour octet
if ( !p->parents().empty() && !q->parents().empty() ) {
octet = ( p->parents()[0] == q->parents()[0] ) &&
( p->parents()[0]->data().iColour() == PDT::Colour8 );
}
// (Final) option: check if same colour8 parent
// or already found an octet.
if(_octetOption==0||octet) return octet;
// (All) option handling more octets
// by browsing particle history/colour lines.
tColinePtr cline,aline;
// Get colourlines form final states.
if(p->hasColour() && q->hasAntiColour()) {
cline = p-> colourLine();
aline = q->antiColourLine();
} else {
cline = q-> colourLine();
aline = p->antiColourLine();
}
// Follow the colourline of p.
if ( !p->parents().empty() ) {
tPPtr parent = p->parents()[0];
while (parent) {
if(parent->data().iColour() == PDT::Colour8) {
// Coulour8 particles should have a colour
// and an anticolour line. Currently the
// remnant has none of those. Since the children
// of the remnant are not allowed to emit currently,
// the colour octet remnant is handled by the return
// statement above. The assert also catches other
// colour octets without clines. If the children of
// a remnant should be allowed to emit, the remnant
// should get appropriate colour lines and
// colour states.
// See Ticket: #407
// assert(parent->colourLine()&&parent->antiColourLine());
octet = (parent-> colourLine()==cline &&
parent->antiColourLine()==aline);
}
if(octet||parent->parents().empty()) break;
parent = parent->parents()[0];
}
}
}
return octet;
}
void ColourReconnector::persistentOutput(PersistentOStream & os) const {
os
<< _hadronSpectrum
+ << _clusterFinder
<< _clreco
<< _crIterations
<< _algorithm
<< _annealingFactor
<< _annealingSteps
<< _triesPerStepFactor
<< _initTemp
<< _preco
<< _precoBaryonic
<< _precoDiquark
<< _preco3M_3M
<< _preco3M_BBbar
<< _precoBbarB_3M
<< _preco2B_2B
<< _precoMB_MB
<< _stepFactor
<< _mesonToBaryonFactor
<< ounit(_maxDistance, femtometer)
<< _octetOption
<< _cr2BeamClusters
<< _localCR
<< _causalCR
<< _debug
<< _junctionMBCR
<< _dynamicCR
<< _diquarkCR
<< _dynamicCRscale
<< _dynamicCRalphaS
<< _phaseSpaceDiquarkFission
<< _cutDiquarkClusterFormation
<< _selectionChoice
;
}
void ColourReconnector::persistentInput(PersistentIStream & is, int) {
is
>> _hadronSpectrum
+ >> _clusterFinder
>> _clreco
>> _crIterations
>> _algorithm
>> _annealingFactor
>> _annealingSteps
>> _triesPerStepFactor
>> _initTemp
>> _preco
>> _precoBaryonic
>> _precoDiquark
>> _preco3M_3M
>> _preco3M_BBbar
>> _precoBbarB_3M
>> _preco2B_2B
>> _precoMB_MB
>> _stepFactor
>> _mesonToBaryonFactor
>> iunit(_maxDistance, femtometer)
>> _octetOption
>> _cr2BeamClusters
>> _localCR
>> _causalCR
>> _debug
>> _junctionMBCR
>> _dynamicCR
>> _diquarkCR
>> _dynamicCRscale
>> _dynamicCRalphaS
>> _phaseSpaceDiquarkFission
>> _cutDiquarkClusterFormation
>> _selectionChoice
;
}
void ColourReconnector::Init() {
static ClassDocumentation<ColourReconnector> documentation
("This class is responsible of the colour reconnection.");
static Reference<ColourReconnector,HadronSpectrum>
interfaceHadronSpectrum("HadronSpectrum",
- "A reference to the object HadronSpectrum",
+ "A reference to the object HadronSpectrum"
+ " Needed for thresholds",
&Herwig::ColourReconnector::_hadronSpectrum,
false, false, true, false);
+ static Reference<ColourReconnector,ClusterFinder>
+ interfaceClusterFinder("ClusterFinder",
+ "A reference to the object ClusterFinder "
+ "Needed for Diquark On shell treatment",
+ &Herwig::ColourReconnector::_clusterFinder,
+ false, false, true, false);
+
static Switch<ColourReconnector,int> interfaceColourReconnection
("ColourReconnection",
"Colour reconnections",
&ColourReconnector::_clreco, 0, true, false);
static SwitchOption interfaceColourReconnectionNo
(interfaceColourReconnection,
"No",
"Colour reconnections off",
0);
static SwitchOption interfaceColourReconnectionYes
(interfaceColourReconnection,
"Yes",
"Colour reconnections on",
1);
static Parameter<ColourReconnector, unsigned int> interfaceColourReconnectionIterations
("ColourReconnectionIterations",
"Choose the number of iterations the chosen CR algorithm is performed",
&ColourReconnector::_crIterations, 1, 1, 1, 100,
false, false, Interface::limited);
// Algorithm interface
static Switch<ColourReconnector, int> interfaceAlgorithm
("Algorithm",
"Specifies the colour reconnection algorithm",
&ColourReconnector::_algorithm, 0, true, false);
static SwitchOption interfaceAlgorithmPlain
(interfaceAlgorithm,
"Plain",
"Plain colour reconnection as in Herwig 2.5.0",
0);
static SwitchOption interfaceAlgorithmStatistical
(interfaceAlgorithm,
"Statistical",
"Statistical colour reconnection using simulated annealing",
1);
static SwitchOption interfaceAlgorithmBaryonic
(interfaceAlgorithm,
"Baryonic",
"Baryonic cluster reconnection",
2);
static SwitchOption interfaceAlgorithmBaryonicMesonic
(interfaceAlgorithm,
"BaryonicMesonic",
"Baryonic cluster reconnection with reconnections to and from Mesonic Clusters",
3);
static SwitchOption interfaceAlgorithmBaryonicDiquarkCluster
(interfaceAlgorithm,
"BaryonicDiquarkCluster",
"Baryonic colour reconnection which allows for the formation of DiquarkCluster-like CR",
4);
static Switch<ColourReconnector,int> interfaceColourDiquarkCR
("DiquarkCR",
"Allow diquark type colour Reconnection."
"NOTE: Necessary to be Yes if BaryonicDiquarkCluster algorithm is chosen.",
&ColourReconnector::_diquarkCR, 0, true, false);
static SwitchOption interfaceDiquarkCRNo
(interfaceColourDiquarkCR,
"No",
"Forbid diquark type colour reconnections",
0);
static SwitchOption interfaceDiquarkCRYes
(interfaceColourDiquarkCR,
"Yes",
"Allow diquark type colour reconnections",
1);
static Switch<ColourReconnector,int> interfaceColourDynamicCR
("DynamicCR",
"Use dynamic weight for Colour reconnections defined by soft gluon evolution"
"\nNOTE: Only availible for Plain, Baryonic, BaryonicDiquarkCluster algorithm.",
&ColourReconnector::_dynamicCR, 0, true, false);
static SwitchOption interfaceDynamicCRNo
(interfaceColourDynamicCR,
"No",
"Use regular CR with fixed probabilities",
0);
static SwitchOption interfaceDynamicCRYes
(interfaceColourDynamicCR,
"Yes",
"Use dynamic CR with kinematic dependent probabilities",
1);
static SwitchOption interfaceDynamicCRNotExponentiated
(interfaceColourDynamicCR,
"NotExponentiated",
"Use dynamic CR with kinematic dependent probabilities"
", but without exponentiated soft anomalous dimension."
"NOTE: Only for testing.",
2);
// General Parameters and switches
static Parameter<ColourReconnector, double> interfaceDynamicScale
("DynamicScale",
"Choose dynamic scale of soft gluon evolution for DynamicCR where"
" mu = DynamicScale*(mConstU+mConstU)",
&ColourReconnector::_dynamicCRscale, 1.0, 1e-14, 1.0,
false, false, Interface::limited);
static Parameter<ColourReconnector, double> interfaceDynamicAlphaS
("DynamicAlphaS",
"Choose dynamic alphaS of soft gluon evolution for DynamicCR",
&ColourReconnector::_dynamicCRalphaS, 0.8, 0.001, 10.0,
false, false, Interface::limited);
// Statistical CR Parameters:
static Parameter<ColourReconnector, double> interfaceMtrpAnnealingFactor
("AnnealingFactor",
"The annealing factor is the ratio of the temperatures in two successive "
"temperature steps.",
&ColourReconnector::_annealingFactor, 0.9, 0.0, 1.0,
false, false, Interface::limited);
static Parameter<ColourReconnector,unsigned> interfaceMtrpAnnealingSteps
("AnnealingSteps",
"Number of temperature steps in the statistical annealing algorithm",
&ColourReconnector::_annealingSteps, 50, 1, 10000,
false, false, Interface::limited);
static Parameter<ColourReconnector,double> interfaceMtrpTriesPerStepFactor
("TriesPerStepFactor",
"The number of reconnection tries per temperature steps is the number of "
"clusters times this factor.",
&ColourReconnector::_triesPerStepFactor, 5.0, 0.0, 100.0,
false, false, Interface::limited);
static Parameter<ColourReconnector,double> interfaceMtrpInitialTemp
("InitialTemperature",
"Factor used to determine the initial temperature from the median of the "
"energy change in a few random rearrangements.",
&ColourReconnector::_initTemp, 0.1, 0.00001, 100.0,
false, false, Interface::limited);
// Plain and Baryonic CR Paramters
static Parameter<ColourReconnector, double> interfaceRecoProb
("ReconnectionProbability",
"Probability that a found two meson to two meson reconnection possibility is actually accepted (used in Plain & Baryonic)",
&ColourReconnector::_preco, 0.5, 0.0, 1.0,
false, false, Interface::limited);
static Parameter<ColourReconnector,double> interfaceRecoProbBaryonic
("ReconnectionProbabilityBaryonic",
"Probability that a found reconnection possibility is actually accepted (used in Baryonic)",
&ColourReconnector::_precoBaryonic, 0.5, 0.0, 1.0,
false, false, Interface::limited);
static Parameter<ColourReconnector,double> interfaceRecoProbDiquark
("ReconnectionProbabilityDiquark",
"Probability for forming a tetra-quark cluster",
&ColourReconnector::_precoDiquark, 0.5, 0.0, 1.0,
false, false, Interface::limited);
// BaryonicMesonic CR Paramters
static Parameter<ColourReconnector, double> interfaceReconnectionProbability3Mto3M
("ReconnectionProbability3Mto3M",
"Probability that a reconnection candidate is accepted for reconnecting 3M -> 3M\'",
&ColourReconnector::_preco3M_3M, 0.5, 0.0, 1.0,
false, false, Interface::limited);
static Parameter<ColourReconnector, double> interfaceReconnectionProbability3MtoBBbar
("ReconnectionProbability3MtoBBbar",
"Probability that a reconnection candidate is accepted for reconnecting 3M -> B,Bbar",
&ColourReconnector::_preco3M_BBbar, 0.5, 0.0, 1.0,
false, false, Interface::limited);
static Parameter<ColourReconnector, double> interfaceReconnectionProbabilityBbarBto3M
("ReconnectionProbabilityBbarBto3M",
"Probability that a reconnection candidate is accepted for reconnecting B,Bbar -> 3M",
&ColourReconnector::_precoBbarB_3M, 0.5, 0.0, 1.0,
false, false, Interface::limited);
static Parameter<ColourReconnector, double> interfaceReconnectionProbability2Bto2B
("ReconnectionProbability2Bto2B",
"Probability that a reconnection candidate is accepted for reconnecting 2B -> 2B\' or 2Bbar -> 2Bbar\'",
&ColourReconnector::_preco2B_2B, 0.5, 0.0, 1.0,
false, false, Interface::limited);
static Parameter<ColourReconnector, double> interfaceReconnectionProbabilityMBtoMB
("ReconnectionProbabilityMBtoMB",
"Probability that a reconnection candidate is accepted for reconnecting M,B -> M\',B\' or M,Bbar -> M\',Bbar\'",
&ColourReconnector::_precoMB_MB, 0.5, 0.0, 1.0,
false, false, Interface::limited);
static Parameter<ColourReconnector, double> interfaceFactorforStep
("StepFactor",
"Factor for how many reconnection-tries are made in the BaryonicMesonic algorithm",
&ColourReconnector::_stepFactor, 1.0, 0.11111, 10.,
false, false, Interface::limited);// at least 3 Clusters -> _stepFactorMin=1/9
static Parameter<ColourReconnector, double> interfaceMesonToBaryonFactor
("MesonToBaryonFactor",
"Factor for comparing mesonic clusters to baryonic clusters in the displacement if BaryonicMesonic CR model is chosen",
&ColourReconnector::_mesonToBaryonFactor, 2.0, 0.01, 100.0,
false, false, Interface::limited);
// General Parameters and switches
static Parameter<ColourReconnector, Length> interfaceMaxDistance
("MaxDistance",
"Maximum distance between the clusters at which to consider rearrangement"
" to avoid colour reconneections of displaced vertices (used in all Algorithms). No unit means femtometer",
&ColourReconnector::_maxDistance, femtometer, 1000.*femtometer, 0.0*femtometer, 1e100*femtometer,
false, false, Interface::limited);
static Switch<ColourReconnector, int> interfaceOctetTreatment
("OctetTreatment",
"Which octets are not allowed to be reconnected (used in all Algorithms)",
&ColourReconnector::_octetOption, 0, false, false);
static SwitchOption interfaceOctetTreatmentFinal
(interfaceOctetTreatment,
"Final",
"Only prevent for the final (usuaslly non-perturbative) g -> q qbar splitting",
0);
static SwitchOption interfaceOctetTreatmentAll
(interfaceOctetTreatment,
"All",
"Prevent for all octets",
1);
static SwitchOption interfaceOctetTreatmentNone
(interfaceOctetTreatment,
"None",
"Accept all octets. "
"NOTE: If a static gluon constituent mass is chosen this option is unphysical. "
"It will lead to mConstGluon fixed mass clusters!",
-1);
static Switch<ColourReconnector, int> interfaceCR2BeamClusters
("CR2BeamClusters",
"Option for colour reconnecting 2 beam remnant clusters if the number of clusters is 2.",
&ColourReconnector::_cr2BeamClusters, 0, true, false);
static SwitchOption interfaceCR2BeamClustersYes
(interfaceCR2BeamClusters,
"Yes",
"If possible CR 2 beam clusters",
1);
static SwitchOption interfaceCR2BeamClustersNo
(interfaceCR2BeamClusters,
"No",
"If possible do not CR 2 beam clusters",
0);
static Switch<ColourReconnector, int> interfaceLocalCR
("LocalCR",
"Option for colour reconnecting only if clusters are less distant than MaxDistance",
&ColourReconnector::_localCR, 0, true, false);
static SwitchOption interfaceLocalCRYes
(interfaceLocalCR,
"Yes",
"activate spatial veto",
1);
static SwitchOption interfaceLocalCRNo
(interfaceLocalCR,
"No",
"deactivate spatial veto",
0);
static Switch<ColourReconnector, int> interfaceCausalCR
("CausalCR",
"Option for colour reconnecting only if clusters their vertices "
"have a positive spacetime difference",
&ColourReconnector::_causalCR, 0, true, false);
static SwitchOption interfaceCausalCRYes
(interfaceCausalCR,
"Yes",
"enable causal veto",
1);
static SwitchOption interfaceCausalCRNo
(interfaceCausalCR,
"No",
"disable causal veto",
0);
static Switch<ColourReconnector, int> interfaceJunction
("Junction",
"Option for using Junction-like displacement in rapidity-phi plane to compare baryonic cluster "
"instead of pairwise distance (for BaryonicMesonic model)",
&ColourReconnector::_junctionMBCR, 1, true, false);
static SwitchOption interfaceJunctionYes
(interfaceJunction,
"Yes",
"Using junction-like model instead of pairwise distance model",
1);
static SwitchOption interfaceJunctionNo
(interfaceJunction,
"No",
"Using pairwise distance model instead of junction-like model",
0);
// Debug
static Switch<ColourReconnector, int> interfaceDebug
("Debug",
"Make a file with some Information of the BaryonicMesonic Algorithm",
&ColourReconnector::_debug, 0, true, false);
static SwitchOption interfaceDebugNo
(interfaceDebug,
"No",
"Debug Information for ColourReconnector Off",
0);
static SwitchOption interfaceDebugYes
(interfaceDebug,
"Yes",
"Debug Information for ColourReconnector On",
1);
static Switch<ColourReconnector,int> interfacePhaseSpaceDiquarkFission
("PhaseSpaceDiquarkFission",
"Only for dynamic colour reconnection choose if capturing cluster decay"
" phase space for formation of a diquark cluster in the transition probabilities",
&ColourReconnector::_phaseSpaceDiquarkFission, 0, true, false);
static SwitchOption interfacePhaseSpaceDiquarkFissionNo
(interfacePhaseSpaceDiquarkFission,
"No",
"Not adding the decay phasespace to Diquark Colour Reconnection",
0);
static SwitchOption interfacePhaseSpaceDiquarkFissionYes
(interfacePhaseSpaceDiquarkFission,
"Yes",
"Adding the decay phasespace to Diquark Colour Reconnection",
1);
static SwitchOption interfacePhaseSpaceDiquarkFissionConstituentMasses
(interfacePhaseSpaceDiquarkFission,
"ConstituentMasses",
"Adding the decay phasespace to Diquark Colour Reconnection",
2);
static Parameter<ColourReconnector, double> interfaceCutDiquarkClusterFormation
("CutDiquarkClusterFormation",
"Cut on diquark cluster formation such that clusters are accepted only if "
"(p1*p2/(m1*m2)-1)<cut for diquark and antidiquark. Note that setting this to 0"
" applies no cut!",
&ColourReconnector::_cutDiquarkClusterFormation, 0.5, 0.0, 100.0,
false, false, Interface::limited);
static Switch<ColourReconnector,int> interfaceSelectionAlgorithm
("SelectionAlgorithm",
"Selection algorithm for choices of subsystems of clusters",
&ColourReconnector::_selectionChoice, 1, true, false);
static SwitchOption interfaceSelectionAlgorithmDefault
(interfaceSelectionAlgorithm,
"Default",
"Default rapidity based similar to Baryonic",
1);
static SwitchOption interfaceSelectionAlgorithmExperimental2
(interfaceSelectionAlgorithm,
"Experimental2",
"Experimental rapidity based on sorting for best rapidity sum",
2);
static SwitchOption interfaceSelectionAlgorithmExperimental3
(interfaceSelectionAlgorithm,
"Experimental3",
"Experimental rapidity based on sorting for best invariant mass sum",
3);
static SwitchOption interfaceSelectionAlgorithmExperimental4
(interfaceSelectionAlgorithm,
"Experimental4",
"Experimental rapidity based on sorting for best rescaled log^2(p1*p2/(m1*m2))",
4);
}
diff --git a/Hadronization/ColourReconnector.h b/Hadronization/ColourReconnector.h
--- a/Hadronization/ColourReconnector.h
+++ b/Hadronization/ColourReconnector.h
@@ -1,718 +1,732 @@
// -*- C++ -*-
//
// ColourReconnector.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_ColourReconnector_H
#define HERWIG_ColourReconnector_H
#include <ThePEG/Interface/Interfaced.h>
#include <unordered_map>
#include "CluHadConfig.h"
#include "ColourReconnector.fh"
#include "HadronSpectrum.h"
+#include "ClusterFinder.h"
#include "Herwig/Utilities/Kinematics.h"
#include <bits/stdc++.h>
namespace Herwig {
using namespace ThePEG;
/** \ingroup Hadronization
* \class ColourReconnector
* \brief Class for changing colour reconnections of partons.
* \author Alberto Ribon, Christian Roehr
*
* This class does the nonperturbative colour rearrangement, after the
* nonperturbative gluon splitting and the "normal" cluster formation.
* It uses the list of particles in the event record, and the collections of
* "usual" clusters which is passed to the main method. If the colour
* reconnection is actually accepted, then the previous collections of "usual"
* clusters is first deleted and then the new one is created.
*
* * @see \ref ColourReconnectorInterfaces "The interfaces"
* defined for ColourReconnector.
*/
class ColourReconnector: public Interfaced {
public:
/**
* Does the colour rearrangement, starting out from the list of particles in
* the event record and the collection of "usual" clusters passed as
* arguments. If the actual rearrangement is accepted, the initial collection of
* clusters is overridden by the old ones.
*/
void rearrange(ClusterVector & clusters);
using CluVecIt = ClusterVector::iterator;
private:
/*
class ColourFlowState2 {
};
class ColourFlowBasis2 {
// represents the colour flow state for 2 colour flows:
// Note: i,j are either (1,2)
// |ij> === |(1,2)->(i,j)>
// where 1,2 are coloured partons and i,j are their
// anticoloured partons
public:
ColourFlowBasis2();
ColourFlowBasis2(const ColourFlowBasis2 & cf2){
setIndex(cf2.getIndex(0),cf2.getIndex(1));
};
bool operator=(const ColourFlowBasis2& cf2){
return setIndex(cf2.getIndex(0),cf2.getIndex(1));
}
bool setIndex(int ind[2]) {
return setIndex(ind[0],ind[1]);
};
bool setIndex(int i, int j) {
if (i==j || i<0 || j<0 || i>1 || j>1){
std::cout << "Error cannot create ColourFlowBasis2 = |"<<i+1<<j+1<<">\n";
return true;
}
_index[0]=i+1;
_index[1]=j+1;
return false;
};
int getIndex(int i) const {
if (i<0 || i>1 ){
std::cout << "Error i = "<<i << std::endl;
assert(false);
}
return _index[i];
};
private:
int _index[2] = {1,2};
};
class ColourFlowBasis3 {
// represents the colour flow state for 2 colour flows:
// Note: i,j,k are either (1,2,3)
// |ijk> === |(1,2,3)->(i,j,k)>
// where 1,2,3 are coloured partons and i,j,k are their
// anticoloured partons
public:
ColourFlowBasis3();
ColourFlowBasis3(const ColourFlowBasis3 & cf3){
setIndex(cf3.getIndex(0),cf3.getIndex(1),cf3.getIndex(2));
};
bool operator=(const ColourFlowBasis3& cf3){
return setIndex(cf3.getIndex(0),cf3.getIndex(1),cf3.getIndex(2));
}
bool setIndex(int ind[3]) {
return setIndex(ind[0],ind[1],ind[2]);
};
bool setIndex(int i, int j, int k) {
if (i==j || i==k || j==k || i<0 || j<0 || k<0 || i>2 || j>2 || k>2){
std::cout << "Error cannot create ColourFlowBasis3 = |"<<i+1<<j+1<<k+1<<">\n";
return true;
}
_index[0]=i+1;
_index[1]=j+1;
_index[2]=k+1;
return false;
};
int getIndex(int i) const {
if (i<0 || i>2 ){
std::cout << "Error i = "<<i << std::endl;
assert(false);
}
return _index[i];
};
private:
int _index[3] = {1,2,3};
};*/
- HadronSpectrumPtr _hadronSpectrum;
+ /*
+ * Reference to HadronSpectrum for getting
+ * threshold settings
+ * */
+ HadronSpectrumPtr _hadronSpectrum;
+
+ /*
+ * Reference to ClusterFinder for getting
+ * Diquark treatment
+ * */
+ ClusterFinderPtr _clusterFinder;
/** PRIVATE MEMBER FUNCTIONS */
/**
* @brief Calculates the sum of the squared cluster masses.
* @arguments q, aq vectors containing the quarks and antiquarks respectively
* @return Sum of cluster squared masses M^2_{q[i],aq[i]}.
*/
Energy2 _clusterMassSum(const PVector & q, const PVector & aq) const;
Selector <int> getProbabilities2CF(const ClusterPtr & c1, const ClusterPtr & c2, bool diquarkCR) const;
Selector <int> _selector(const ClusterVector & clusters, bool diquarkCR) const;
Selector <int> _selectorCF3(const ClusterVector & clusters, bool diquarkCR) const;
std::unordered_map<int,double> _reconnectionAmplitudesSGE(const ClusterVector & clusters, int nonTrivialInitialState=123) const;
int _stateToPermutation(const int i) const;
/**
* @brief Calculates the Mesonic reconnection Probability from soft gluon evolution
* @arguments c1, c2 cluster pointer refs, whose kinematics determine the probability
* @return probability in [0:1]
*/
std::tuple<double,double,double> _dynamicRecoProbabilitiesCF2(const ClusterPtr & c1, const ClusterPtr & c2, bool diquarkCR) const;
std::unordered_map<int,double> _reconnectionAmplitudesCF2 (const ClusterPtr & c1, const ClusterPtr & c2, const int nonTrivialInitialState = 0) const;
/**
* @brief calculates the "euclidean" distance of two quarks in the
* rapdidity-phi plane
* @arguments p, q the two quarks
* @return the dimensionless distance:
* \deltaR_12=sqrt(\deltaY_12^2 + \delta\phi_12^2)
*/
double _displacement(tcPPtr p, tcPPtr q) const;
/**
* @brief calculates the "euclidean" distance of a the 3 (anti)quarks
* (anti)baryonic cluster in the rapdidity-phi plane
* @arguments p, q the two quarks
* @return the dimensionless distance:
* if Junction is enabled the difference of all 3 quarks
* with respect to their mean point is calculated:
* <Y> = sum_i Y_i/3
* <\phi> = sum_i \phi_i/3
* \deltaR = sum_i sqrt( (Y_i - <Y>)^2 + (\phi_i - <phi>)^2)
*
* if Junction is disabled the difference of all distinct
* pairing of the 3 quarks is computed:
* \deltaR_ij = sqrt(\deltaY_ij^2 + \delta\phi_ij^2)
* \deltaR_tot = sum_{i,j<i} \deltaR_ij
*
* NOTE: switch the Junction option will necessarily
* need to be combined with a change (or re-tune)
* of _mesonToBaryonFactor otherwise no well
* description is to be expected
* TODO: maybe add a different p-norm option to get
* more phenomenology
*/
double _displacementBaryonic(tcPPtr p1, tcPPtr p2, tcPPtr p3) const;
/**
* @brief Examines whether the cluster vector (under the given permutation of
* the antiquarks) contains colour-octet clusters
* @param cv Cluster vector
* @param P Permutation, a vector of permutated indices from 0 to
* cv.size()-1
*/
bool _containsColour8(const ClusterVector & cv, const vector<size_t> & P) const;
/**
* @brief A Metropolis-type algorithm which finds a local minimum in the
* total sum of cluster masses
* @arguments cv cluster vector
*/
void _doRecoStatistical(ClusterVector & cv) const;
/**
* @brief Plain colour reconnection as used in Herwig 2.5.0
* @arguments cv cluster vector
*/
void _doRecoPlain(ClusterVector & cv) const;
/**
* Baryonic Colour Reconnection model
*/
void _doRecoBaryonic(ClusterVector & cv) const;
void _doRecoBaryonicDiquarkCluster(ClusterVector & cv) const;
/**
* @brief BaryonicMesonic colour reconnection model
* @arguments cv cluster vector
* BaryonicMesonic Colour Reconnection model with reconnections from mesonic clusters
* to baryonic cluster and the contrary. Allows also reconnection between mesonic
* and baryonic Clusters
*/
void _doRecoBaryonicMesonic(ClusterVector & cv) const;
void _makeBaryonicClusters(ClusterPtr &c1, ClusterPtr &c2, ClusterPtr &c3,
ClusterPtr &newcluster1, ClusterPtr &newcluster2) const;
/**
* Method to allow one to make four (light) quark clusters
*/
bool _makeDiquarkCluster(ClusterPtr &c1, ClusterPtr &c2, ClusterPtr &newcluster) const;
pair <ClusterPtr,ClusterPtr> _splitDiquarkCluster(ClusterPtr &diquarkCluster, bool colourReconnect) const;
+
bool _canMakeDiquarkCluster(tcPPtr pCol1, tcPPtr pCol2,tcPPtr pAntiCol1, tcPPtr pAntiCol2) const;
bool _canMakeDiquarkCluster(const ClusterPtr &c1, const ClusterPtr &c2) const;
bool _canMakeDiquarkCluster(tcPPtr pCol1, tcPPtr pCol2,tcPPtr pAntiCol1, tcPPtr pAntiCol2, double &PhaseSpace) const;
bool _canMakeDiquarkCluster(const ClusterPtr &c1, const ClusterPtr &c2, double &PhaseSpace) const;
+ bool _canMakeBaryonicCluster(vector<tcPPtr> pCol) const;
+ bool _canMakeBaryonicCluster(const ClusterPtr &c1, const ClusterPtr &c2, const ClusterPtr &c3) const;
/**
* @brief Finds the cluster in cv which, if reconnected with the given
* cluster cl, would result in the smallest sum of cluster masses.
* If no reconnection partner can be found, a pointer to the
* original Cluster cl is returned.
* @arguments cv cluster vector
* cl cluster iterator (must be from cv) which wants to have a reconnection partner
* @return iterator to the found cluster, or the original cluster pointer if
* no mass-reducing combination can be found
*/
CluVecIt _findRecoPartnerPlain(const CluVecIt & cl, ClusterVector & cv, const ClusterVector & deleted) const;
CluVecIt _findRecoPartnerPlainDynamic(const CluVecIt & cl, ClusterVector & cv, const ClusterVector & deleted, bool diquarkCR) const;
CluVecIt _findPartnerRapidity(const CluVecIt & cl, ClusterVector & cv) const;
CluVecIt _findPartnerBaryonic(const CluVecIt & cl, ClusterVector & cv,
bool & tetraCand,
const ClusterVector& a,
CluVecIt &baryonic1,
CluVecIt &baryonic2 ) const;
void _findPartnerBaryonicDiquarkClusterTEST2( const CluVecIt & cl, ClusterVector & cv,
unsigned & typeOfReconnection,
const ClusterVector& deleted,
CluVecIt &candidate1,
CluVecIt &candidate2 ) const;
void _findPartnerBaryonicDiquarkClusterTEST3( const CluVecIt & cl, ClusterVector & cv,
unsigned & typeOfReconnection,
const ClusterVector& deleted,
CluVecIt &candidate1,
CluVecIt &candidate2 ) const;
void _findPartnerBaryonicDiquarkCluster( const CluVecIt & cl, ClusterVector & cv,
unsigned & typeOfReconnection,
const ClusterVector& deleted,
CluVecIt &candidate1,
CluVecIt &candidate2 ) const;
/**
* @brief Finds best CR option for the BaryonicMesonic CR model
* @arguments cls vector of selected clusters, baryonic is the number of baryonic
* clusters in selection, kind_of_reco is output string denoting best
* CR option and reco_info is output storage for information on which
* (anti-)quarks to exchange and in which way.
* BaryonicMesonic Colour Reconnection model with reconnections from mesonic clusters
* to baryonic cluster and the contrary. Allows also reconnection between mesonic
* and baryonic Clusters
*/
void _findbestreconnectionoption(std::vector<CluVecIt> &cls,
const unsigned &baryonic,
string &kind_of_reco,
int (&reco_info)[3]) const;
/**
* @brief Reconnects the constituents of the given clusters to the (only)
* other possible cluster combination.
* @return pair of pointers to the two new clusters
* Used for Plain and Baryonic Colour Reconnection models
*/
pair <ClusterPtr,ClusterPtr> _reconnect(ClusterPtr &c1, ClusterPtr &c2) const;
std::tuple <ClusterPtr, ClusterPtr> _reconnect3MtoMD(ClusterVector & cluvec, const int topology) const;
/**
* @brief Reconnects (2B->2B) the constituents of the given clusters to
another possible cluster combination whose information is given
in s1 and s2.
* @arguments c1 and c2 are baryonic clusters and s1 and s2 are the respective
indices of the constituents of c1 and c2 respectively
* @return pair of pointers to the two new clusters
* Used only in BaryonicMesonic algorithm and will exchange constituent s1 of
* c1 with constituent s2 of c2
*/
pair <ClusterPtr,ClusterPtr> _reconnect2Bto2B(ClusterPtr &c1, ClusterPtr &c2, const int s1, const int s2) const;
/**
* @brief Reconnects (B,Bbar->3M) the constituents of the given clusters to
another possible cluster combination whose information is given in
s0, s1 and s2.
* @arguments c1 and c2 are baryonic clusters. s0, s1 and s2 are the respective
indices which determine the CR
* @return tuple of pointers to the 3 new mesonic clusters
* Used only in BaryonicMesonic algorithm and will form 3 mesonic clusters according
* to the indices s0, s1 and s2. The i-th constituent of c1 is connected to the si-th
* constituent of c2
*/
std::tuple <ClusterPtr, ClusterPtr, ClusterPtr> _reconnectBbarBto3M(ClusterPtr &c1, ClusterPtr &c2, const int s0, const int s1, const int s2 ) const;
/**
* @brief Reconnects (3M->3M) the constituents of the given clusters to
another possible cluster combination whose information is given in
sinfos
* @arguments c1, c2 and c3 are mesonic clusters. infos[3] are the respective
indices which determine the CR
* @return tuple of pointers to the 3 CR'd mesonic clusters
* Used only in BaryonicMesonic algorithm and will reconnect 3 mesonic clusters according
* to the infos, which determine the CR. The coloured quark of the i-th cluster forms
* a new cluster with the anticoloured quark of the info[i]-th cluster
*/
std::tuple <ClusterPtr, ClusterPtr, ClusterPtr> _reconnect3Mto3M(ClusterPtr &c1, ClusterPtr &c2, ClusterPtr &c3, const int infos[3] ) const;
/**
* Reconnection method for a Baryonic and a Mesonic Cluster to a Baryonic and a Mesonic Cluster
* s0 is the Number of the (Anti)Patrticle of the Baryonic Cluster , which should be swapped with the Anti(Particle) of the Mesonic Cluster
*/
/**
* @brief Reconnects the constituents of the given clusters to
another possible cluster combination whose information
is given in s0.
* @arguments c1 and c2 are one baryonic and one mesonic cluster respectively
and s0 is the respective index of the constituents of the baryonic
cluster which is to be exchangeed with the mesonic cluster.
* @return pair of pointers to the two new clusters
* Used only in BaryonicMesonic algorithm and will exchange constituent s0 of
* the baryonic cluster with the (anti)-quark of the mesonic cluster
*/
pair <ClusterPtr, ClusterPtr> _reconnectMBtoMB(ClusterPtr &c1, ClusterPtr &c2, const int s0) const;
/**
* Methods for the BaryonicMesonic CR algorithm
* Find the best reconnection option for the respective cluster-combination
*
*/
/**
* @brief veto for far apart clusters
* @arguments expects at most 3 CluVecIt in clu vector
* @return returns true if clusters are more distant than _maxDistance
* in space
* TODO: problematic maybe add option to turn off
*/
bool _clustersFarApart( const std::vector<CluVecIt> & clu ) const;
/**
* @brief Does reconnect 2 beam clusters for BaryonicMesonic CR model
* if option CR2BeamClusters is enabled
* @arguments cv cluster vector
*/
void _doReco2BeamClusters(ClusterVector & cv) const;
/**
* @brief finds the best reconnection option and stores it in bswap1
* and bswap2 (2B->2B colour reconnection)
* @arguments c1 and c2 cluster pointers and kind_of_reco will output
* the best found reconnection option for c1 and c2
*/
void _2Bto2BreconnectionFinder(ClusterPtr &c1, ClusterPtr &c2, int &bswap1, int &bswap2, double mindisplsum, string &kind_of_reco ) const;
/**
* @brief finds the best reconnection option and stores it in mswap0
* mswap1 and mswap2 (BbarB->3M colour reconnection)
* @arguments c1 and c2 cluster pointers and kind_of_reco will output
* the best found reconnection option for c1 and c2
*/
void _BbarBto3MreconnectionFinder(ClusterPtr &c1, ClusterPtr &c2, int &mswap0, int &mswap1, int &mswap2,
double min_displ_sum, string & kind_of_reco) const;
/**
* @brief finds the best reconnection option and stores it in swap
* (BM->BM colour reconnection)
* @arguments c1 and c2 cluster pointers and kind_of_reco will output
* the best found reconnection option for c1 and c2
*/
void _BMtoBMreconnectionfinder(ClusterPtr &c1, ClusterPtr &c2, int &swap,
double min_displ_sum, string &kind_of_reco) const;
/**
* @brief finds the best reconnection option and stores it in swap0,
* swap1 and swap2 (3M->{3M,BbarB} colour reconnection)
* @arguments c1 and c2 cluster pointers and kind_of_reco will output
* the best found reconnection option for c1 and c2
*/
void _3MtoXreconnectionfinder(std::vector<CluVecIt> &cv, int &swap0, int &swap1,
int &swap2, double min_displ_sum, string &kind_of_reco) const;
/**
* @brief At random, swap two antiquarks, if not excluded by the
* constraint that there must not be any colour-octet clusters.
* @arguments q, aq vectors containing the quarks and antiquarks respectively
* maxtries maximal number of tries to find a non-colour-octet
* reconfiguration
* @return Pair of ints indicating the indices of the antiquarks to be
* swapped. Returns (-1,-1) if no valid reconfiguration could be
* found after maxtries trials
*/
pair <int,int>
_shuffle(const PVector & q, const PVector & aq, unsigned maxtries = 10) const;
/**
* @return true, if we would reconnect two mesonic clusters to a single gluon
* */
bool _isColour8Forbidden(int state, const ClusterVector & clusters) const;
/**
* @return true, if we would reconnect two mesonic clusters to at least
* one gluon
*/
bool _becomesColour8Cluster(const ClusterPtr & c1, const ClusterPtr & c2) const;
/**
* @return true, if the two partons are splitting products of the same
* gluon
*/
bool _isColour8(tcPPtr p, tcPPtr q) const;
/** DATA MEMBERS */
/**
* Specifies the colour reconnection algorithm to be used.
*/
int _algorithm = 0;
/**
* Do we do colour reconnections?
*/
int _clreco = 0;
/**
* Number of iterations of the Colour reconnecion algorithm
*/
unsigned int _crIterations = 1;
/**
* Do we want debug informations?
*/
int _debug = 0;
/**
* @brief Junction-like model for BaryonicMesonic model
* Do we want to use the junction-like model for
* computing the displacements of BaryonicMesonic model?
* otherwise pairwise distances are used.
* If _junctionMBCR is activated the displacements are computed in the
* rapidity-Phi plane by difference to the average rapidity and phi:
* DeltaR_i^2 = (rapidity_i - meanRap)^2 + (phi_i - meanPhi)^2
* DeltaR = sum_i DeltaR_i
* if _junctionMBCR=0 the displacements are computed:
* DeltaR_ij^2 = (rapidity_i - rapidity_j)^2 + (phi_i - phi_j)^2
* DeltaR = sum_i,j<i DeltaR_ij
*/
int _junctionMBCR = 1;
/**
* Statistical Reco:
* Factor used to determine the initial temperature according to
* InitialTemperature = _initTemp * median {energy changes in a few random
* rearrangements}
*/
double _initTemp = 0.01;
/**
* Statistical Reco:
* The annealing factor is the ratio of two successive temperature steps:
* T_n = _annealingFactor * T_(n-1)
*/
double _annealingFactor = 0.21;
/**
* Statistical Reco:
* Number of temperature steps in the statistical annealing algorithm
*/
unsigned int _annealingSteps = 10;
/**
* Statistical Reco:
* The number of tries per temperature steps is the number of clusters times
* this factor.
*/
double _triesPerStepFactor = 0.66;
/**
* Probability that a found reconnection possibility is actually accepted.
* used in Plain & Baryonic CR
*/
double _preco = 0.5;
/**
* Probability that a found reconnection possibility is actually accepted.
* used in Baryonic CR
*/
double _precoBaryonic = 0.5;
/**
* Probability that a found reconnection possibility is actually accepted.
* used in Baryonic CR
*/
double _precoDiquark = 0.5;
/**
* Probability that a found reconnection possibility is actually accepted.
* For reconnecting 3M -> 3M'
* used in BaryonicMesonic
* NOTE: if 0 this type of reconnection is not even tried
*/
double _preco3M_3M = 0.5;
/**
* Probability that a found reconnection possibility is actually accepted.
* For reconnecting 3M -> B,Bbar
* used in BaryonicMesonic
*/
double _preco3M_BBbar = 0.5;
/**
* Probability that a found reconnection possibility is actually accepted.
* For reconnecting Bbar,B -> 3M
* used in BaryonicMesonic
*/
double _precoBbarB_3M = 0.5;
/**
* Probability that a found reconnection possibility is actually accepted.
* For reconnecting 2B -> 2B' or 2Bbar -> 2Bbar'
* used in BaryonicMesonic
* NOTE: if 0 this type of reconnection is not even tried
*/
double _preco2B_2B = 0.5;
/**
* Probability that a found reconnection possibility is actually accepted.
* For reconnecting M,B -> M',B' or M,Bbar -> M',Bbar'
* used in BaryonicMesonic
* NOTE: if 0 this type of reconnection is not even tried
*/
double _precoMB_MB = 0.5;
/**
* For the BaryonicMesonic algorithm
* How many times do suggest cluster for reconnection?
* n(reconnectionstries) = _stepFactor * n(clusters)*n(clusters);
*/
double _stepFactor = 5.0;
/**
* Factor for comparing mesonic clusters to baryonic clusters
*/
double _mesonToBaryonFactor = 2.0;
/**
* Maximium distance for reconnections
* TODO: remove if issues with anticausality are discussed and resolved
*/
Length _maxDistance = femtometer;
/**
* Option for handling octets
*/
int _octetOption = 0;
/**
* Option for colour reconnecting 2 Beam Clusters if no others are present
*/
int _cr2BeamClusters = 0;
/**
* Option for colour reconnecting Clusters only if their vertex 3-distance
* is less than _maxDistance
*/
int _localCR = 0;
/**
* Option for colour reconnecting Clusters only if their spacetime difference
* is bigger than 0
*/
int _causalCR = 0;
/**
* Option for doing the Dynamic CR probability depending on soft
* gluon evolution
*/
int _dynamicCR = 0;
/**
* Option for doing the diquark colour Reconnection
*/
int _diquarkCR = 0;
/**
* Choose Dynamic CR probability scale depending on soft
* gluon evolution
*/
double _dynamicCRscale = 1.0;
/**
* Choose Dynamic CR probability scale depending on soft
* gluon evolution
*/
double _dynamicCRalphaS = 0.8;
/**
* switch for including PhaseSpace in
* Diquark transition probability
* */
int _phaseSpaceDiquarkFission = 1;
/**
* Cut on diquark formation such that clusters are accepted
* only if (p1*p2/(m1*m2)-1)<cut for diquark and antidiquark.
* Note that setting this to 0 applies no cut!
* */
double _cutDiquarkClusterFormation = 0.0;
/**
* Selecting the Selection algorithm for Diquark CR
* */
int _selectionChoice = 1;
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);
//@}
/**
* Standard Init function used to initialize the interfaces.
*/
static void Init();
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.
*/
ColourReconnector & operator=(const ColourReconnector &) = delete;
};
}
#endif /* HERWIG_ColourReconnector_H */
diff --git a/Hadronization/LightClusterDecayer.cc b/Hadronization/LightClusterDecayer.cc
--- a/Hadronization/LightClusterDecayer.cc
+++ b/Hadronization/LightClusterDecayer.cc
@@ -1,408 +1,422 @@
// -*- C++ -*-
//
// LightClusterDecayer.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 LightClusterDecayer class.
//
#include "LightClusterDecayer.h"
#include <ThePEG/Interface/ClassDocumentation.h>
#include <ThePEG/Interface/Parameter.h>
#include <ThePEG/Interface/Reference.h>
#include <ThePEG/Interface/Switch.h>
#include <ThePEG/Persistency/PersistentOStream.h>
#include <ThePEG/Persistency/PersistentIStream.h>
#include <ThePEG/PDT/EnumParticles.h>
#include <ThePEG/Repository/EventGenerator.h>
#include "Cluster.h"
#include "Herwig/Utilities/Kinematics.h"
#include <ThePEG/Utilities/DescribeClass.h>
using namespace Herwig;
DescribeClass<LightClusterDecayer,Interfaced>
describeLightClusterDecayer("Herwig::LightClusterDecayer","Herwig.so");
IBPtr LightClusterDecayer::clone() const {
return new_ptr(*this);
}
IBPtr LightClusterDecayer::fullclone() const {
return new_ptr(*this);
}
void LightClusterDecayer::persistentOutput(PersistentOStream & os) const {
os << _hadronSpectrum;
}
void LightClusterDecayer::persistentInput(PersistentIStream & is, int) {
is >> _hadronSpectrum;
}
void LightClusterDecayer::Init() {
static ClassDocumentation<LightClusterDecayer> documentation
("There is the class responsible for the one-hadron decay of light clusters");
static Reference<LightClusterDecayer,HadronSpectrum>
interfaceHadronSpectrum("HadronSpectrum",
"A reference to the HadronSpectrum object",
&Herwig::LightClusterDecayer::_hadronSpectrum,
false, false, true, false);
}
bool LightClusterDecayer::decay(ClusterVector & clusters, tPVector & finalhadrons) {
// Loop over all clusters, and for those that were not heavy enough
// to undergo to fission, check if they are below the threshold
// for normal two-hadron decays. If this is the case, then the cluster
// should be decayed into a single hadron: this can happen only if
// it is possible to reshuffle momenta between the cluster and
// another one; in the rare occasions in which such exchange of momenta
// is not possible (because all of the clusters are too light) then
// the event is skipped.
// Notice that, differently from what happens in Fortran Herwig,
// light (that is below the threshold for the production of the lightest
// pair of hadrons with the proper flavours) fission products, produced
// by the fission of heavy clusters in class ClusterFissioner
// have been already "decayed" into single hadron (the lightest one
// with proper flavour) by the same latter class, without requiring
// any reshuffling. Therefore the light clusters that are treated in
// this LightClusterDecayer class are produced directly
// (originally) by the class ClusterFinder.
// To preserve all of the information, the cluster partner with which
// the light cluster (that decays into a single hadron) exchanges
// momentum in the reshuffling procedure is redefined and inserted
// in the vector vecNewRedefinedCluPtr. Only at the end, when all
// light clusters have been examined, the elements this vector will be
// copied in collecCluPtr (the reason is that it is not allowed to
// modify a STL container while iterating over it. At the same time,
// this ensures that a cluster can be redefined only once, which seems
// sensible although not strictly necessary).
// Notice that the cluster reshuffling partner is normally redefined
// and inserted in the vector vecNewRedefinedCluPtr, but not always:
// in the case it is also light, then it is also decayed immediately
// into a single hadron, without redefining it (the reason being that,
// otherwise, the would-be redefined cluster could have undefined
// components).
vector<tClusterPtr> redefinedClusters;
for (ClusterVector::const_iterator it = clusters.begin();
it != clusters.end(); ++it) {
// Skip the clusters that are not available or that are
// heavy, intermediate, clusters that have undergone to fission,
if ( ! (*it)->isAvailable() || ! (*it)->isReadyToDecay() ){
continue;
}
// We need to require (at least at the moment, maybe in the future we
// could change it) that the cluster has exactly two components,
// because otherwise we don't know how to deal with the kinematics.
// If this is not the case, then send a warning because it is not suppose
// to happen, and then do nothing with (ignore) such cluster.
if ( (*it)->numComponents() != 2 ) {
generator()->logWarning( Exception("LightClusterDecayer::decay "
"***Still cluster with not exactly"
" 2 components*** ",
Exception::warning) );
continue;
}
if ( DiquarkMatcher::Check((*it)->particle(0)->dataPtr()->id()) && DiquarkMatcher::Check((*it)->particle(1)->dataPtr()->id())) {
// We ignore the Diquark Clusters in the LightClusterDecayer as potential reshuffling partners
// throw Exception() << "LightClusterDecayer::decay\n"
// "*** Diquark Cluster in LightClusterDecayer ***\n"
// "Cluster = ( "<< (*it)->particle(0)->dataPtr()->id()<<", " << (*it)->particle(1)->dataPtr()->id()<<" )\nMC = " << (*it)->mass()/GeV << " GeV\nMLHP = "
// << _hadronSpectrum->massLightestHadronPair((*it)->particle(0)->dataPtr(),(*it)->particle(1)->dataPtr())/GeV <<" GeV"
// << Exception::runerror;
continue;
}
// select the hadron for single hadron decay
tcPDPtr hadron = _hadronSpectrum->chooseSingleHadron((*it)->particle(0)->dataPtr(),
(*it)->particle(1)->dataPtr(),
(**it).mass());
// if not single decay continue
if(!hadron){
continue;
}
// We assume that the candidate reshuffling cluster partner,
// with whom the light cluster can exchange momenta,
// is chosen as the closest in space-time between the available
// clusters. Notice that an alternative, sensible approach
// could be to consider instead the "closeness" in the colour
// structure...
// Notice that nor a light cluster (which decays into a single hadron)
// neither its cluster reshuffling partner (which either has a
// redefined cluster or also decays into a single hadron) can be
// a reshuffling partner of another light cluster.
// This because we are requiring that the considered candidate cluster
// reshuffling partner has the status "isAvailable && isReadyToDecay" true;
// furthermore, the new redefined clusters are not added to the collection
// of cluster before the end of the entire reshuffling procedure, avoiding
// in this way that the redefined cluster of a cluster reshuffling partner
// is used again later. Needless to say, this is just an assumption,
// although reasonable, but nothing more than that!
// Build a multimap of available reshuffling cluster partners,
// with key given by the module of the invariant space-time distance
// w.r.t. the light cluster, so that this new collection is automatically
// ordered in increasing distance values.
// We use a multimap, rather than a map, just for precaution against not properly
// defined cluster positions which could produce all identical (null) distances.
multimap<Length,tClusterPtr> candidates;
for ( ClusterVector::iterator jt = clusters.begin();
jt != clusters.end(); ++jt ) {
- if ( (*jt)->numComponents() != 2 )
+ if ( (*jt)->numComponents() != 2 ){
continue;
+ }
+ if ( DiquarkMatcher::Check((*jt)->particle(0)->dataPtr()->id())
+ || DiquarkMatcher::Check((*jt)->particle(1)->dataPtr()->id())) {
+ // Avoid reshuffling with diquark or baryonic clusters.
+ // This could lead to momentum violation in rare cases.
+ // TODO: implement this consistently
+ continue;
+ }
if ((*jt)->isAvailable() && (*jt)->isReadyToDecay() && jt != it) {
Length distance = abs (((*it)->vertex() - (*jt)->vertex()).m());
candidates.insert(pair<Length,tClusterPtr>(distance,*jt));
}
}
// Loop sequentially the multimap.
multimap<Length,tClusterPtr>::const_iterator mmapIt = candidates.begin();
bool found = false;
while (!found && mmapIt != candidates.end()) {
found = reshuffling(hadron, *it, (*mmapIt).second, redefinedClusters, finalhadrons);
if (!found) ++mmapIt;
}
- if (!found) return partonicReshuffle(hadron,*it,finalhadrons);
+ if (!found) {
+ bool success = partonicReshuffle(hadron,*it,finalhadrons);
+ if (!success){
+ return false;
+ }
+ // return partonicReshuffle(hadron,*it,finalhadrons);
+ }
} // end loop over collecCluPtr
// Add to collecCluPtr all of the redefined new clusters (indeed the
// pointers to them are added) contained in vecNewRedefinedCluPtr.
for (tClusterVector::const_iterator it = redefinedClusters.begin();
it != redefinedClusters.end(); ++it) {
clusters.push_back(*it);
}
return true;
}
bool LightClusterDecayer::reshuffling(const tcPDPtr pdata1,
tClusterPtr cluPtr1,
tClusterPtr cluPtr2,
tClusterVector & redefinedClusters,
tPVector & finalhadrons)
{
// don't reshuffle with beam clusters
if(cluPtr2->isBeamCluster()) return false;
// This method does the reshuffling of momenta between the cluster "1",
// that must decay into a single hadron (with id equal to idhad1), and
// the candidate cluster "2". It returns true if the reshuffling succeed,
// false otherwise.
PPtr ptrhad1 = pdata1->produceParticle();
if ( ! ptrhad1 ) {
generator()->logWarning( Exception("LightClusterDecayer::reshuffling"
"***Cannot create a particle with specified id***",
Exception::warning) );
return false;
}
Energy mhad1 = ptrhad1->mass();
// Let's call "3" and "4" the two constituents of the second cluster
tPPtr part3 = cluPtr2->particle(0);
tPPtr part4 = cluPtr2->particle(1);
// Check if the system of the two clusters can kinematically be replaced by
// an hadron of mass mhad1 (which is the lightest single hadron with the
// same flavour numbers as the first cluster) and the second cluster.
// If not, then try to replace the second cluster with the lightest hadron
// with the same flavour numbers; if it still fails, then give up!
Lorentz5Momentum pSystem = cluPtr1->momentum() + cluPtr2->momentum();
pSystem.rescaleMass(); // set the mass as the invariant of the quadri-vector
Energy mSystem = pSystem.mass();
Energy mclu2 = cluPtr2->mass();
bool singleHadron = false;
bool isDiquarkCluster = DiquarkMatcher::Check(part3->dataPtr()->id()) && DiquarkMatcher::Check(part4->dataPtr()->id());
Energy mLHP2 = _hadronSpectrum->massLightestHadronPair(part3->dataPtr(),part4->dataPtr());
// avoid calling massLightestHadron for Diquark clusters and only allow kinematic reshuffling
// for diquark clusters (no singleHadron)
Energy mLH2 = isDiquarkCluster ? mSystem:_hadronSpectrum->massLightestHadron(part3->dataPtr(),part4->dataPtr());
if(mSystem > mhad1 + mclu2 && mclu2 > mLHP2) { singleHadron = false; }
else if(mSystem > mhad1 + mLH2) { singleHadron = true; mclu2 = mLH2; }
else return false;
// Let's call from now on "Sys" the system of the two clusters, and
// had1 (of mass mhad1) the lightest hadron in which the first
// cluster decays, and clu2 (of mass mclu2) either the second
// cluster or the lightest hadron in which it decays (depending
// which one is kinematically allowed, see above).
// The idea behind the reshuffling is to replace the system of the
// two clusters by the system of the hadron had1 and (cluster or hadron) clu2,
// but leaving the overall system unchanged. Furthermore, the motion
// of had1 and clu2 in the Sys frame is assumed to be parallel to, respectively,
// those of the original cluster1 and cluster2 in the same Sys frame.
// Calculate the unit three-vector, in the frame "Sys" along which the
// two initial clusters move.
Lorentz5Momentum u( cluPtr1->momentum() );
u.boost( - pSystem.boostVector() ); // boost from LAB to Sys
// Calculate the momenta of had1 and clu2 in the Sys frame first,
// and then boost back in the LAB frame.
Lorentz5Momentum phad1, pclu2;
if (pSystem.m() < mhad1 + mclu2 ) {
throw Exception() << "Impossible Kinematics in LightClusterDecayer::reshuffling()"
<< Exception::eventerror;
}
if (!(u.vect().mag2() > ZERO) ) {
generator()->logWarning( Exception("Impossible Kinematics in LightClusterDecayer::reshuffling()\n"
"Cluster System in Rest frame has |p3| = 0 GeV",
Exception::warning) );
return false;
}
Kinematics::twoBodyDecay(pSystem, mhad1, mclu2, u.vect().unit(), phad1, pclu2);
ptrhad1->set5Momentum( phad1 ); // set momentum of first hadron.
ptrhad1->setVertex(cluPtr1->vertex()); // set hadron vertex position to the
// parent cluster position.
cluPtr1->addChild(ptrhad1);
finalhadrons.push_back(ptrhad1);
cluPtr1->flagAsReshuffled();
cluPtr2->flagAsReshuffled();
if(singleHadron) {
// In the case that also the cluster reshuffling partner is light
// it is decayed into a single hadron, *without* creating the
// redefined cluster (this choice is justified in order to avoid
// clusters that could have undefined components).
PPtr ptrhad2 = _hadronSpectrum->lightestHadron(part3->dataPtr(),part4->dataPtr())
->produceParticle();
ptrhad2->set5Momentum( pclu2 );
ptrhad2->setVertex( cluPtr2->vertex() ); // set hadron vertex position to the
// parent cluster position.
cluPtr2->addChild(ptrhad2);
finalhadrons.push_back(ptrhad2);
} else {
// Create the new cluster which is the redefinitions of the cluster
// partner (cluster "2") used in the reshuffling procedure of the
// light cluster (cluster "1").
// The rationale of this is to preserve completely all of the information.
ClusterPtr cluPtr2new = ClusterPtr();
if(part3 && part4) cluPtr2new = new_ptr(Cluster(part3,part4));
cluPtr2new->set5Momentum( pclu2 );
cluPtr2new->setVertex( cluPtr2->vertex() );
cluPtr2->addChild( cluPtr2new );
redefinedClusters.push_back( cluPtr2new );
// Set consistently the momenta of the two components of the second cluster
// after the reshuffling. To do that we first calculate the momenta of the
// constituents in the initial cluster rest frame; then we boost them back
// in the lab but using this time the new cluster rest frame. Finally we store
// these information in the new cluster. Notice that we do *not* set
// consistently also the momenta of the (eventual) particles pointed by the
// two components: that's because we do not need to do so, being the momentum
// an explicit private member of the class Component (which is set equal
// to the momentum of the eventual particle pointed only in the constructor,
// but then later should not necessary be the same), and furthermore it allows
// us not to loose any information, in the sense that we can always, later on,
// to find the original momenta of the two components before the reshuffling.
Lorentz5Momentum p3 = part3->momentum(); //p3new->momentum();
p3.boost( - (cluPtr2->momentum()).boostVector() ); // from LAB to clu2 (old) frame
p3.boost( pclu2.boostVector() ); // from clu2 (new) to LAB frame
Lorentz5Momentum p4 = part4->momentum(); //p4new->momentum();
p4.boost( - (cluPtr2->momentum()).boostVector() ); // from LAB to clu2 (old) frame
p4.boost( pclu2.boostVector() ); // from clu2 (new) to LAB frame
cluPtr2new->particle(0)->set5Momentum(p3);
cluPtr2new->particle(1)->set5Momentum(p4);
} // end of if (singleHadron)
return true;
}
bool LightClusterDecayer::partonicReshuffle(const tcPDPtr had,
const PPtr cluster,
tPVector & finalhadrons) {
tPPtr meson(cluster);
if(!meson->parents().empty()) meson=meson->parents()[0];
if(!meson->parents().empty()) meson=meson->parents()[0];
// check b/c hadron decay
int ptype(abs(meson->id())%10000);
bool heavy = (ptype/1000 == 5 || ptype/1000 ==4 );
heavy |= (ptype/100 == 5 || ptype/100 ==4 );
heavy |= (ptype/10 == 5 || ptype/10 ==4 );
if(!heavy) return false;
// find the leptons
tPVector leptons;
for(unsigned int ix=0;ix<meson->children().size();++ix) {
if(!(meson->children()[ix]->dataPtr()->coloured())) {
leptons.push_back(meson->children()[ix]);
}
}
if(leptons.size()==1) {
tPPtr w=leptons[0];
leptons.pop_back();
for(unsigned int ix=0;ix<w->children().size();++ix) {
if(!w->children()[ix]->dataPtr()->coloured()) {
leptons.push_back(w->children()[ix]);
}
}
}
if(leptons.size()!=2) return false;
// get momentum of leptonic system and the its minimum possible mass
Energy mmin(ZERO);
Lorentz5Momentum pw;
for(unsigned int ix=0;ix<leptons.size();++ix) {
pw+=leptons[ix]->momentum();
mmin+=leptons[ix]->mass();
}
pw.rescaleMass();
// check we can do the reshuffling
PPtr ptrhad = had->produceParticle();
// total momentum fo the system
Lorentz5Momentum pSystem = pw + cluster->momentum();
pSystem.rescaleMass();
// normal case get additional energy by rescaling momentum in rest frame of
// system
if(pSystem.mass()>ptrhad->mass()+pw.mass()&&pw.mass()>mmin) {
// Calculate the unit three-vector, in the frame "Sys" along which the
// two initial clusters move.
Lorentz5Momentum u(cluster->momentum());
u.boost( - pSystem.boostVector() );
// Calculate the momenta of had1 and clu2 in the Sys frame first,
// and then boost back in the LAB frame.
Lorentz5Momentum phad1, pclu2;
Kinematics::twoBodyDecay(pSystem, ptrhad->mass(), pw.mass(),
u.vect().unit(), phad1, pclu2);
// set momentum of first hadron.
ptrhad->set5Momentum( phad1 );
// set hadron vertex position to the parent cluster position.
ptrhad->setLabVertex(cluster->vertex());
// add hadron
cluster->addChild(ptrhad);
finalhadrons.push_back(ptrhad);
// reshuffle the leptons
// boost the leptons to the rest frame of the system
Boost boost1(-pw.boostVector());
Boost boost2( pclu2.boostVector());
for(unsigned int ix=0;ix<leptons.size();++ix) {
leptons[ix]->deepBoost(boost1);
leptons[ix]->deepBoost(boost2);
}
return true;
}
else {
return false;
}
}
diff --git a/src/defaults/Hadronization.in b/src/defaults/Hadronization.in
--- a/src/defaults/Hadronization.in
+++ b/src/defaults/Hadronization.in
@@ -1,185 +1,186 @@
# -*- ThePEG-repository -*-
############################################################
# Setup of default hadronization
#
# There are no user servicable parts inside.
#
# Anything that follows below should only be touched if you
# know what you're doing.
#############################################################
cd /Herwig/Particles
create ThePEG::ParticleData Cluster
setup Cluster 81 Cluster 0.00990 0.0 0.0 0.0 0 0 0 1
create ThePEG::ParticleData Remnant
setup Remnant 82 Remnant 0.00990 0.0 0.0 0.0 0 0 0 1
mkdir /Herwig/Hadronization
cd /Herwig/Hadronization
create Herwig::ClusterHadronizationHandler ClusterHadHandler
create Herwig::PartonSplitter PartonSplitter
create Herwig::ClusterFinder ClusterFinder
create Herwig::ColourReconnector ColourReconnector
create Herwig::ClusterFissioner ClusterFissioner
create Herwig::MatrixElementClusterFissioner MEClusterFissioner
create Herwig::LightClusterDecayer LightClusterDecayer
create Herwig::ClusterDecayer ClusterDecayer
create Herwig::HwppSelector SMHadronSpectrum
newdef ClusterHadHandler:PartonSplitter PartonSplitter
newdef ClusterHadHandler:ClusterFinder ClusterFinder
newdef ClusterHadHandler:ColourReconnector ColourReconnector
newdef ClusterHadHandler:ClusterFissioner ClusterFissioner
newdef ClusterHadHandler:LightClusterDecayer LightClusterDecayer
newdef ClusterHadHandler:ClusterDecayer ClusterDecayer
do ClusterHadHandler:UseHandlersForInteraction QCD
newdef ClusterHadHandler:MinVirtuality2 0.1*GeV2
newdef ClusterHadHandler:MaxDisplacement 1.0e-10*millimeter
newdef ClusterHadHandler:UnderlyingEventHandler NULL
newdef PartonSplitter:HadronSpectrum SMHadronSpectrum
newdef ClusterFinder:HadronSpectrum SMHadronSpectrum
newdef ClusterFissioner:HadronSpectrum SMHadronSpectrum
newdef ClusterDecayer:HadronSpectrum SMHadronSpectrum
newdef LightClusterDecayer:HadronSpectrum SMHadronSpectrum
# ColourReconnector Default Parameters
newdef ColourReconnector:HadronSpectrum SMHadronSpectrum
+newdef ColourReconnector:ClusterFinder ClusterFinder
newdef ColourReconnector:ColourReconnection Yes
newdef ColourReconnector:Algorithm Baryonic
# Statistical CR Parameters:
newdef ColourReconnector:AnnealingFactor 0.9
newdef ColourReconnector:AnnealingSteps 50
newdef ColourReconnector:TriesPerStepFactor 5.0
newdef ColourReconnector:InitialTemperature 0.1
# Plain and Baryonic CR Paramters
newdef ColourReconnector:ReconnectionProbability 0.95
newdef ColourReconnector:ReconnectionProbabilityBaryonic 0.7
# BaryonicMesonic and BaryonicMesonic CR Paramters
newdef ColourReconnector:ReconnectionProbability3Mto3M 0.5
newdef ColourReconnector:ReconnectionProbability3MtoBBbar 0.5
newdef ColourReconnector:ReconnectionProbabilityBbarBto3M 0.5
newdef ColourReconnector:ReconnectionProbability2Bto2B 0.05
newdef ColourReconnector:ReconnectionProbabilityMBtoMB 0.5
newdef ColourReconnector:StepFactor 1.0
newdef ColourReconnector:MesonToBaryonFactor 1.333
# General Parameters and switches
newdef ColourReconnector:MaxDistance 1.0e50
newdef ColourReconnector:OctetTreatment Final
newdef ColourReconnector:CR2BeamClusters No
newdef ColourReconnector:Junction Yes
newdef ColourReconnector:LocalCR No
newdef ColourReconnector:CausalCR No
# Debugging
newdef ColourReconnector:Debug No
# set ClusterFissioner parameters
newdef ClusterFissioner:KinematicThreshold Dynamic
newdef ClusterFissioner:KineticThresholdShift 0.08844
newdef ClusterFissioner:ProbabilityPowerFactor 6.486
newdef ClusterFissioner:ProbabilityShift -0.87875
# Clustering parameters for light quarks
newdef ClusterFissioner:ClMaxLight 3.528693*GeV
newdef ClusterFissioner:ClPowLight 1.849375
newdef ClusterFissioner:PSplitLight 0.914156
insert ClusterFissioner:FissionPwt 1 1.0
insert ClusterFissioner:FissionPwt 2 1.0
insert ClusterFissioner:FissionPwt 3 0.374094
newdef ClusterDecayer:ClDirLight 1
newdef ClusterDecayer:ClSmrLight 0.78
#
# Cluster Paramters for light Diquark Cluster
# currently set according to Light quark defaults
newdef ClusterFissioner:ClMaxDiquark 3.528693*GeV
newdef ClusterFissioner:ClPowDiquark 1.849375
# Clustering parameters for b-quarks
insert ClusterFissioner:ClMaxHeavy 5 3.757*GeV
insert ClusterFissioner:ClPowHeavy 5 0.547
insert ClusterFissioner:PSplitHeavy 5 0.625
insert ClusterDecayer:ClDirHeavy 5 1
insert ClusterDecayer:ClSmrHeavy 5 0.078
newdef SMHadronSpectrum:SingleHadronLimitBottom 0.000
# Clustering parameters for c-quarks
insert ClusterFissioner:ClMaxHeavy 4 3.950*GeV
insert ClusterFissioner:ClPowHeavy 4 2.559
insert ClusterFissioner:PSplitHeavy 4 0.994
insert ClusterDecayer:ClDirHeavy 4 1
insert ClusterDecayer:ClSmrHeavy 4 0.163
newdef SMHadronSpectrum:SingleHadronLimitCharm 0.000
# Clustering parameters for exotic quarks
# (e.g. hadronizing Susy particles)
newdef ClusterFissioner:ClMaxExotic 2.7*GeV
newdef ClusterFissioner:ClPowExotic 1.46
newdef ClusterFissioner:PSplitExotic 1.00
newdef ClusterDecayer:ClDirExotic 1
newdef ClusterDecayer:ClSmrExotic 0.
newdef SMHadronSpectrum:SingleHadronLimitExotic 0.
################################################
# BEG MEClusterFissioner initialization: #
################################################
# Note that we copy the tuned values from the default
# ClusterFissioner but we need to set them here explicitly
newdef MEClusterFissioner:HadronSpectrum SMHadronSpectrum
newdef MEClusterFissioner:KinematicThreshold Dynamic
newdef MEClusterFissioner:KineticThresholdShift 0.08844
newdef MEClusterFissioner:ProbabilityPowerFactor 6.486
newdef MEClusterFissioner:ProbabilityShift -0.87875
# Clustering parameters for light quarks
newdef MEClusterFissioner:ClMaxLight 3.528693*GeV
newdef MEClusterFissioner:ClPowLight 1.849375
newdef MEClusterFissioner:PSplitLight 0.914156
insert MEClusterFissioner:FissionPwt 1 1.0
insert MEClusterFissioner:FissionPwt 2 1.0
insert MEClusterFissioner:FissionPwt 3 0.374094
# Clustering parameters for b-quarks
insert MEClusterFissioner:ClMaxHeavy 5 3.757*GeV
insert MEClusterFissioner:ClPowHeavy 5 0.547
insert MEClusterFissioner:PSplitHeavy 5 0.625
# Clustering parameters for c-quarks
insert MEClusterFissioner:ClMaxHeavy 4 3.950*GeV
insert MEClusterFissioner:ClPowHeavy 4 2.559
insert MEClusterFissioner:PSplitHeavy 4 0.994
# Cluster Paramters for light Diquark Cluster
# currently set according to Light quark defaults
newdef MEClusterFissioner:ClMaxDiquark 3.528693*GeV
newdef MEClusterFissioner:ClPowDiquark 1.849375
# Clustering parameters for exotic quarks
# (e.g. hadronizing Susy particles)
newdef MEClusterFissioner:ClMaxExotic 2.7*GeV
newdef MEClusterFissioner:ClPowExotic 1.46
newdef MEClusterFissioner:PSplitExotic 1.00
################################################
# END MEClusterFissioner initialization: #
################################################
#
insert PartonSplitter:SplitPwt 1 1.0
insert PartonSplitter:SplitPwt 2 1.0
insert PartonSplitter:SplitPwt 3 0.824135
newdef PartonSplitter:Split Light
#
newdef SMHadronSpectrum:PwtDquark 1.0
newdef SMHadronSpectrum:PwtUquark 1.0
newdef SMHadronSpectrum:PwtSquark 0.374094
newdef SMHadronSpectrum:PwtCquark 0.0
newdef SMHadronSpectrum:PwtBquark 0.0
newdef SMHadronSpectrum:PwtDIquark 0.33107
newdef SMHadronSpectrum:SngWt 0.89050
newdef SMHadronSpectrum:DecWt 0.41628
newdef SMHadronSpectrum:Mode 1
newdef SMHadronSpectrum:BelowThreshold All
create Herwig::SpinHadronizer SpinHadronizer
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Sat, Dec 21, 6:20 PM (8 h, 41 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023772
Default Alt Text
(235 KB)
Attached To
rHERWIGHG herwighg
Event Timeline
Log In to Comment