Page MenuHomeHEPForge
Files F8308485

No OneTemporary
Actions

View Options

diff --git a/Hadronization/ColourReconnector.cc b/Hadronization/ColourReconnector.cc
--- a/Hadronization/ColourReconnector.cc
+++ b/Hadronization/ColourReconnector.cc
@@ -1,2515 +1,2515 @@
// -*- 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/Utilities/Maths.h"
#include "Herwig/Utilities/expm-1.h"
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
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;
// do the colour reconnection
switch (_algorithm) {
case 0:
_doRecoPlain(clusters);
break;
case 1:
_doRecoStatistical(clusters);
break;
case 2:
_doRecoBaryonic(clusters);
break;
case 3:
_doRecoBaryonicMesonic(clusters);
break;
case 4:
_doRecoBaryonicDiquarkCluster(clusters);
break;
}
}
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);
}
// 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;
return -1;
}
unsigned int scalarProducts(int i, int j);
unsigned int scalarProducts(int i, int j)
{
if (i>j) return scalarProducts(j,i);
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;
}
}
return Nc;
}
double ColourReconnector::_kinematicRecoProbabilityFixedScale(const ClusterPtr c1, const ClusterPtr c2, bool diquarkCR) const{
// Verified
Lorentz5Momentum p1_c = c1->colParticle()->momentum();
Lorentz5Momentum p1_a = c1->antiColParticle()->momentum();
Lorentz5Momentum p2_c = c2->colParticle()->momentum();
Lorentz5Momentum p2_a = c2->antiColParticle()->momentum();
double M12 = (p1_c + p1_a).m2()/sqr(_dynamicCRscale);
double M24 = (p1_a + p2_a).m2()/sqr(_dynamicCRscale);
double M13 = (p1_c + p2_c).m2()/sqr(_dynamicCRscale);
double M23 = (p1_a + p2_c).m2()/sqr(_dynamicCRscale);
double M14 = (p1_c + p2_a).m2()/sqr(_dynamicCRscale);
double M34 = (p2_c + p2_a).m2()/sqr(_dynamicCRscale);
double alphaQCD=_dynamicCRalphaS;
double logSqrOmega12=alphaQCD*pow(log(M12),2)/twopi;
double logSqrOmega24=alphaQCD*pow(log(M24),2)/twopi;
double logSqrOmega13=alphaQCD*pow(log(M13),2)/twopi;
double logSqrOmega23=alphaQCD*pow(log(M23),2)/twopi;
double logSqrOmega14=alphaQCD*pow(log(M14),2)/twopi;
double logSqrOmega34=alphaQCD*pow(log(M34),2)/twopi;
double a=pow(log(sqr(logSqrOmega23))+log(sqr(logSqrOmega14)),2);
double b=pow(log(sqr(logSqrOmega13))+log(sqr(logSqrOmega24)),2);
double c=pow(log(sqr(logSqrOmega12))+log(sqr(logSqrOmega34)),2);
double sqrtDelta=sqrt(9*a*a-4*c*(a+b)-14*a*b+9*b*b+4*c*c);
double U11=sqrtDelta/tanh(sqrtDelta/2.0)+3*(b-a);
double U21=2*(c-b);
double N=3.0;
double denominator=sqr(U11+N*U21)+sqr(N*U11+U21)+(N+1.0)/(2*(N-1.0))*sqr(U11-U21);
double pNormalCR=sqr(U11+N*U21)/denominator;
double pDiquarkCR=(N+1.0)/(2*(N-1.0))*sqr(U11-U21)/denominator;
// if (p<0.1 || p>0.9) std::cout << "p = "<< p << std::endl;
assert( pNormalCR<=1.0 && pNormalCR>=0.0);
assert( pDiquarkCR<=1.0 && pDiquarkCR>=0.0);
if (_debug)
{
ofstream out("WriteOut/kinematicRecoProbability.dat", std::ios::app);
out << pNormalCR << "\t" << pDiquarkCR << "\n";
out.close();
}
if (diquarkCR) return pDiquarkCR;
return pNormalCR;
}
std::vector<double> ColourReconnector::_precoKinematicSGE(const ClusterVector & clusters) const {
std::vector<double> preco;
int size=clusters.size();
assert(clusters.size()<4);
switch(size){
case 0:
break;
case 1:
break;
case 2:
{
// std::cout << "BETTER NEED TODO That" << std::endl;
preco.push_back(_kinematicRecoProbabilityFixedScale(clusters[0],clusters[1]));
break;
}
case 3:
{
// test
if (clusters[0]->numComponents()!=2 ||
clusters[1]->numComponents()!=2 ||
clusters[2]->numComponents()!=2 )
return preco;
const int N=6; // 2*Nclu;
double omIJ[N][N],scalarProds[N][N];
double alphaQCD=_dynamicCRalphaS;
Lorentz5Momentum mom_i, mom_j;
for (int i = 0; i < N; i++)
{
if ((i+1)%2==0) 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==0) mom_j=clusters[j/2]->colParticle()->momentum();
else mom_j=clusters[(j-1)/2]->antiColParticle()->momentum();
// omIJ[i][j]=alphaQCD*pow(log(sqr(clusters[i]->momentum()+clusters[j]->momentum())/sqr(_dynamicCRscale)),2)/twopi;
omIJ[i][j]=alphaQCD*pow(log(sqr(mom_i+mom_j)/sqr(_dynamicCRscale)),2)/twopi;
scalarProds[i][j]=scalarProducts(i,j);
}
}
boost::numeric::ublas::matrix<double> * Uevolve = new boost::numeric::ublas::matrix<double>(N,N);
boost::numeric::ublas::matrix<double> Omega(N,N);
int Nc=3;
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[4][1]);
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[0][1]+omIJ[2][3]+omIJ[4][5]);
+ Omega(5,5) = - 0.5 * Nc * (omIJ[2][3]+omIJ[1][4]+omIJ[0][5]);
*Uevolve=expm_pad(Omega);
// std::cout << *Uevolve << std::endl;
std::vector<double> amp1toJ; // <J|1>
double sumOfAll=0;
double sumBaryon=0;
for (int J = 0; J < N; J++)
{
double sum=0;
for (int i = 0; i < N; i++)
{
//TODO
sum+=pow((*Uevolve)(i,0)*scalarProds[i][J],2);
if (!J)
{
double NB=2.0/sqrt(3.0);
for (int k = 0; k < N; k++)
{
sumBaryon+=(*Uevolve)(i,0)*signPermutationState(k)*scalarProds[i][k]/NB;
}
}
}
sumOfAll+=sum;
amp1toJ.push_back(sum);
}
sumBaryon=sumBaryon*sumBaryon; //square
sumOfAll+=sumBaryon;
amp1toJ.push_back(sumBaryon);
double onetest=0.0;
for (int i = 0; i < N; i++)
{
amp1toJ[i]/=sumOfAll;
onetest+=amp1toJ[i];
// TODO needs decent review
// std::cout << "CR from 1 to "<<i+1<< "state probability: "<< std::scientific<< amp1toJ[i] << std::endl;
}
amp1toJ[N]/=sumOfAll; //Baryon
onetest+=amp1toJ[N];
if (fabs(onetest-1)>1e-14)
{
std::cout << "Not sum to 1 but 1-sum P_i="<<std::scientific<< 1-onetest << std::endl;
for (int i = 0; i < N+1; i++)
{
std::cout << "\t\t\tp["<<i+1<<"]"<<amp1toJ[i]<<"\n";
}
}
// std::cout << "# CR sanity 1="<< std::scientific<< onetest << "state probability: " << std::endl;
delete Uevolve;
return amp1toJ;
break;
}
default:
{
std::cerr << "ERROR : in CR _precoKinematicSGE" << std::endl;
break;
}
}
return preco;
}
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;
}
void ColourReconnector::_doRecoPlain(ClusterVector & cv) const {
ClusterVector newcv = cv;
// 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
// NB this method returns *cit if no reconnection partner can be found
CluVecIt candidate = _findRecoPartner(cit, newcv);
// skip this cluster if no possible reshuffling partner can be found
if (candidate == cit) continue;
// accept the reconnection with probability PrecoProb.
double PrecoProb = _dynamicCR ? _kinematicRecoProbabilityFixedScale(*cit,*candidate):_preco;
if (UseRandom::rnd() < PrecoProb) {
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;
}
}
swap(cv,newcv);
return;
}
namespace {
inline bool hasDiquark(CluVecIt cit) {
for (unsigned int i = 0; i<(*cit)->numComponents(); i++) {
if (DiquarkMatcher::Check(*((*cit)->particle(i)->dataPtr())))
return true;
}
return false;
}
}
void ColourReconnector::_doRecoBaryonicDiquarkCluster(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;
double ProbabilityDiquark = _precoDiquark;
// 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;
// Find a candidate suitable for reconnection
CluVecIt candidate1, candidate2;
unsigned typeOfReconnection=0;
_findPartnerBaryonicDiquarkCluster(cit, newcv,
typeOfReconnection,
deleted,
candidate1,
candidate2);
if (_dynamicCR) {
switch (typeOfReconnection)
{
case 1:
{
ProbabilityMesonic = _kinematicRecoProbabilityFixedScale(*cit, *candidate1);
break;
}
case 2:
{
ClusterVector cluvec={*cit,*candidate1,*candidate2};
std::vector<double> precos=_precoKinematicSGE(cluvec);
ProbabilityBaryonic = precos[6];
break;
}
case 3:
{
ProbabilityDiquark = _kinematicRecoProbabilityFixedScale(*cit, *candidate1, true);
break;
}
}
}
switch (typeOfReconnection)
{
// Mesonic CR
case 1:
if (UseRandom::rnd() < ProbabilityMesonic) {
auto reconnected = _reconnect(*cit, *candidate1);
*cit = reconnected.first;
*candidate1 = reconnected.second;
}
break;
// Baryonic CR
case 2:
if (UseRandom::rnd() < ProbabilityBaryonic) {
deleted.push_back(*candidate2);
// Function that does the reconnection from 3 -> 2 clusters
ClusterPtr b1, b2;
_makeBaryonicClusters(*cit,*candidate1,*candidate2, b1, b2);
*cit = b1;
*candidate1 = b2;
}
break;
// Diquark CR
case 3:
if (UseRandom::rnd() < ProbabilityDiquark){
// We will delete the candidate1 mesonic clusters
// to form a diquark cluster
ClusterPtr DiqCluster;
if (_makeDiquarkCluster(*cit, *candidate1, DiqCluster)){
deleted.push_back(*candidate1);
*cit = DiqCluster;
}
}
break;
// No CR found
case 0:
continue;
default:
assert(false);
}
}
// 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);
}
// 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;
// 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) {
if (isBaryonicCandidate) {
ClusterVector cluvec={*cit,*baryonic1,*baryonic2};
std::vector<double> precos=_precoKinematicSGE(cluvec);
ProbabilityBaryonic = precos[6];
}
else {
ProbabilityMesonic = _kinematicRecoProbabilityFixedScale(*cit, *candidate);
}
}
// 3 aligned meson case
if ( isBaryonicCandidate
&& UseRandom::rnd() < ProbabilityBaryonic ) {
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) {
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
// TODO: after discussion maybe to be omitted
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
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
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) {
// 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
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
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);
}
namespace {
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 Energy pt = sqrt(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() + (pt*pt));
// Correct formula
const double y2 = log((p2.t() + abs(pz))/mtrans);
return ( pz < ZERO ) ? -y2 : y2;
}
}
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(
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();
// TODO Boost not covariant!! ERROR
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();
// TODO Boost not covariant!! ERROR
// 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;
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 normal 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) {
if (candIsOctet1)
typeOfReconnection = 0;
else
typeOfReconnection = 1;
}
else if (candIsQQ1)
{
if (candIsQQ2)
typeOfReconnection = 2;
else
typeOfReconnection = 3;
}
if (!foundCR) typeOfReconnection = 0;
// veto reconnection if cannot make a Diquark Cluster
bool failDCR=false;
if (typeOfReconnection == 3) {
if (!_canMakeDiquarkCluster(*cl, *candidate1)) {
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 (_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);
}
}
}
CluVecIt ColourReconnector::_findPartnerBaryonic(
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::_findRecoPartner(CluVecIt cl,
ClusterVector & cv) const {
CluVecIt candidate = cl;
Energy minMass = 1*TeV;
for (CluVecIt cit=cv.begin(); cit != cv.end(); ++cit) {
// 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());
}
// forms a four-quark cluster
bool ColourReconnector::_canMakeDiquarkCluster(
const ClusterPtr &c1, const ClusterPtr &c2) const {
//make sure they all have 2 components
assert(c1->numComponents()==2);
assert(c2->numComponents()==2);
// Stop Heavy quarks from entering
if (abs(c1->colParticle()->dataPtr()->id()) > 3 ||
abs(c1->antiColParticle()->dataPtr()->id()) > 3 ||
abs(c2->colParticle()->dataPtr()->id()) > 3 ||
abs(c2->antiColParticle()->dataPtr()->id()) > 3 ){
// std::cout << "heavy reject" << std::endl;
return false;
}
// ClusterPtr newClusterCheck=new_ptr(Cluster(c1->colParticle(),c2->colParticle(), c1->antiColParticle(), c2->antiColParticle()));
// TODO replace this with simply getting the diquark data Ptr
// ClusterPtr DiquarkCluster=_handleDiquarkCluster(newClusterCheck);
tcPDPtr dataDiquark = _hadronSpectrum->makeDiquark(c1->colParticle()->dataPtr(), c2->colParticle()->dataPtr());
tcPDPtr dataDiquarkBar = _hadronSpectrum->makeDiquark(c1->antiColParticle()->dataPtr(), c2->antiColParticle()->dataPtr());
if (!dataDiquark){
throw Exception() << "Could not make a diquark from"
<< c1->colParticle()->dataPtr()->PDGName() << " and "
<< c2->colParticle()->dataPtr()->PDGName()
<< " in ClusterFinder::handleDiquarkCluster()"
<< Exception::eventerror;
}
if (!dataDiquarkBar){
throw Exception() << "Could not make an anti-diquark from"
<< c1->antiColParticle()->dataPtr()->PDGName() << " and "
<< c2->antiColParticle()->dataPtr()->PDGName()
<< " in ClusterFinder::handleDiquarkCluster()"
<< Exception::eventerror;
}
Energy minMass=_hadronSpectrum->massLightestHadronPair(dataDiquark,dataDiquarkBar);
Lorentz5Momentum Ptot=c1->colParticle()->momentum()+c2->colParticle()->momentum()
+ c1->antiColParticle()->momentum()+c2->antiColParticle()->momentum();
// if (Ptot.m()!=sqrt(sqr(Ptot)))
// std::cout << "Ptot.m " << ounit(Ptot.m(),GeV) <<"\t sqrt(sqr Pttot) "<< ounit(sqrt(sqr(Ptot)),GeV) << std::endl;;
// Energy Mass=Ptot.m();
Energy Mass=Ptot.mass();
Energy Mass2=sqrt(Ptot*Ptot);
Energy Mass3=sqrt(sqr(Ptot));
if ( fabs((Mass-Mass2)/GeV) > 1e-6) std::cout << "DMass12 = "<< std::scientific << ounit(Mass-Mass2,GeV) <<" Mass = " << ounit(Mass,GeV) << " GeV Mass2 " << ounit(Mass2,GeV) << "\n";
if ( fabs((Mass2-Mass3)/GeV) > 1e-6) std::cout << "DMass23 = "<< std::scientific << ounit(Mass2-Mass3,GeV) <<" Mass2 = " << ounit(Mass2,GeV) << " GeV Mass3 " << ounit(Mass3,GeV) << "\n";
if ( fabs((Mass-Mass3)/GeV) > 1e-6) std::cout << "DMass13 = "<< std::scientific << ounit(Mass-Mass3,GeV) <<" Mass = " << ounit(Mass,GeV) << " GeV Mass3 " << ounit(Mass3,GeV) << "\n";
// if (fabs(Ptot.massError() )>1e-14) std::cout << "Mass inconsistency : " << std::scientific << (Ptot.massError()) <<"\n";
if ( Mass<minMass ) {
// std::cout << "reject Diquark" << std::endl;
return false;
}
return true;
}
bool ColourReconnector::_makeDiquarkCluster(
ClusterPtr &c1, ClusterPtr &c2,
ClusterPtr &newcluster) const{
if (!_canMakeDiquarkCluster(c1,c2)) return false;
//abandon children
c1->colParticle()->abandonChild(c1);
c1->antiColParticle()->abandonChild(c1);
c2->colParticle()->abandonChild(c2);
c2->antiColParticle()->abandonChild(c2);
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::_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 {
// 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, 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++) {
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::_isColour8(tcPPtr p, tcPPtr q) const {
bool octet = false;
// 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
<< _clreco
<< _algorithm
<< _annealingFactor
<< _annealingSteps
<< _triesPerStepFactor
<< _initTemp
<< _preco
<< _precoBaryonic
<< _preco3M_3M
<< _preco3M_BBbar
<< _precoBbarB_3M
<< _preco2B_2B
<< _precoMB_MB
<< _stepFactor
<< _mesonToBaryonFactor
<< ounit(_maxDistance, femtometer)
<< _octetOption
<< _cr2BeamClusters
<< _localCR
<< _causalCR
<< _debug
<< _junctionMBCR
<< _precoDiquark
<< _dynamicCR
<< ounit(_dynamicCRscale, GeV)
<< _dynamicCRalphaS
;
}
void ColourReconnector::persistentInput(PersistentIStream & is, int) {
is
>> _hadronSpectrum
>> _clreco
>> _algorithm
>> _annealingFactor
>> _annealingSteps
>> _triesPerStepFactor
>> _initTemp
>> _preco
>> _precoBaryonic
>> _preco3M_3M
>> _preco3M_BBbar
>> _precoBbarB_3M
>> _preco2B_2B
>> _precoMB_MB
>> _stepFactor
>> _mesonToBaryonFactor
>> iunit(_maxDistance, femtometer)
>> _octetOption
>> _cr2BeamClusters
>> _localCR
>> _causalCR
>> _debug
>> _junctionMBCR
>> _precoDiquark
>> _dynamicCR
>> iunit(_dynamicCRscale, GeV)
>> _dynamicCRalphaS
;
}
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",
&Herwig::ColourReconnector::_hadronSpectrum,
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);
// 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> interfaceColourDynamicCR
("DynamicCR",
"Use dynamic weight for Colour reconnections defined by soft gluon evolution"
"\nNOTE: Only for Mesonic CR so far",
&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);
// General Parameters and switches
static Parameter<ColourReconnector, Energy> interfaceDynamicScale
("DynamicScale",
"Choose dynamic scale of soft gluon evolution for DynamicCR",
&ColourReconnector::_dynamicCRscale, GeV, 1.*GeV, 0.001*GeV, 1e4*GeV,
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, 1.0, 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, unsigned 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 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);
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 12:13 PM (1 d, 17 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4022712
Default Alt Text
(85 KB)

Event Timeline