diff --git a/Hadronization/ColourReconnector.cc b/Hadronization/ColourReconnector.cc
--- a/Hadronization/ColourReconnector.cc
+++ b/Hadronization/ColourReconnector.cc
@@ -1,2266 +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]);
+ *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;
}
-
-double ColourReconnector::_kinematicRecoProbabilityFixedScale(const ClusterPtr& c1, const ClusterPtr& c2) const{
- // Verified
- // _dynamicCRscale
- 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 p=1.0/(1.0+pow((N*U11+U21)/(U11+N*U21),2));
- // if (p<0.1 || p>0.9) std::cout << "p = "<< p << std::endl;
- assert( p<=1.0 && p>=0.0);
- assert( !std::isnan(p) && !std::isinf(p));
- // if (_writeOut)
- // {
- // ofstream out("WriteOut/kinematicRecoProbability.dat", std::ios::app);
- // out << p << endl;
- // out.close();
- // }
- return p;
- // random choice for kinematics as toy test
- // return 0.1+0.8*UseRandom::rnd();
-}
-
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;
- // accept the reconnection with probability PrecoProb.
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 Diquark clusters, this biases the
+ // 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() < _preco) {
+ if (UseRandom::rnd() < ProbabilityMesonic) {
auto reconnected = _reconnect(*cit, *candidate1);
*cit = reconnected.first;
*candidate1 = reconnected.second;
}
break;
// Baryonic CR
case 2:
- if (UseRandom::rnd() < _precoBaryonic) {
+ 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() < _precoDiquark){
+ 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 (isBaryonicCandidate
- && UseRandom::rnd() < _precoBaryonic) {
- deleted.push_back(*baryonic2);
+ 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) {
- double PrecoProb = _dynamicCR ? _kinematicRecoProbabilityFixedScale(*cit,*candidate):_preco;
- if ( UseRandom::rnd() < PrecoProb) {
- auto reconnected = _reconnect(*cit, *candidate);
- *cit = reconnected.first;
- *candidate = reconnected.second;
- }
- }
+ 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();
- // TODO Boost not covariant!! ERROR
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();
- // 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();
- // TODO Boost not covariant!! ERROR
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();
- // TODO Boost not covariant!! ERROR
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);
}
diff --git a/Hadronization/ColourReconnector.h b/Hadronization/ColourReconnector.h
--- a/Hadronization/ColourReconnector.h
+++ b/Hadronization/ColourReconnector.h
@@ -1,575 +1,578 @@
// -*- 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 "CluHadConfig.h"
#include "ColourReconnector.fh"
#include "HadronSpectrum.h"
#include "Herwig/Utilities/Kinematics.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:
HadronSpectrumPtr _hadronSpectrum;
/** 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;
-
+
+
+ std::vector<double> _precoKinematicSGE(const ClusterVector & clusters) 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]
*/
- double _kinematicRecoProbabilityFixedScale(const ClusterPtr& c1, const ClusterPtr& c2) const;
+ double _kinematicRecoProbabilityFixedScale(const ClusterPtr c1, const ClusterPtr c2, bool diquarkCR=false) 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;
bool _canMakeDiquarkCluster(const ClusterPtr &c1, const ClusterPtr &c2) 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 _findRecoPartner(CluVecIt cl, ClusterVector & cv) const;
CluVecIt _findPartnerRapidity(CluVecIt cl, ClusterVector & cv) const;
CluVecIt _findPartnerBaryonic(CluVecIt cl, ClusterVector & cv,
bool & tetraCand,
const ClusterVector& a,
CluVecIt &baryonic1,
CluVecIt &baryonic2 ) const;
void _findPartnerBaryonicDiquarkCluster( 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;
/**
* @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;
/** DATA MEMBERS */
/**
* Specifies the colour reconnection algorithm to be used.
*/
int _algorithm = 0;
/**
* Do we do colour reconnections?
*/
int _clreco = 0;
/**
* 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;
/**
* @return true, if the two partons are splitting products of the same
* gluon
*/
bool _isColour8(tcPPtr p, tcPPtr q) const;
/**
* Option for handling octets
*/
unsigned 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;
/**
* Choose Dynamic CR probability scale depending on soft
* gluon evolution
*/
Energy _dynamicCRscale = 1.0*GeV;
/**
* Choose Dynamic CR probability scale depending on soft
* gluon evolution
*/
double _dynamicCRalphaS = 0.8;
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/Utilities/expm-1.h b/Utilities/expm-1.h
--- a/Utilities/expm-1.h
+++ b/Utilities/expm-1.h
@@ -1,147 +1,148 @@
//
// Copyright (c) 2007
// Tsai, Dung-Bang
// National Taiwan University, Department of Physics
//
// E-Mail : dbtsai (at) gmail.com
// Begine : 2007/11/20
// Last modify : 2007/11/22
// Version : v0.1
//
// EXPGM_PAD computes the matrix exponential exp(H) for general matrixs,
// including complex and real matrixs using the irreducible (p,p) degree
// rational Pade approximation to the exponential
// exp(z) = r(z)=(+/-)( I+2*(Q(z)/P(z))).
//
// Usage :
//
// U = expm_pad(H)
// U = expm_pad(H, p)
//
// where p is internally set to 6 (recommended and gererally satisfactory).
//
// See also MATLAB supplied functions, EXPM and EXPM1.
//
// Reference :
// EXPOKIT, Software Package for Computing Matrix Exponentials.
// ACM - Transactions On Mathematical Software, 24(1):130-156, 1998
//
// Permission to use, copy, modify, distribute and sell this software
// and its documentation for any purpose is hereby granted without fee,
// provided that the above copyright notice appear in all copies and
// that both that copyright notice and this permission notice appear
// in supporting documentation. The authors make no representations
// about the suitability of this software for any purpose.
// It is provided "as is" without express or implied warranty.
//
#ifndef _BOOST_UBLAS_EXPM_
#define _BOOST_UBLAS_EXPM_
#include <complex>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/lu.hpp>
+#include <boost/throw_exception.hpp>
namespace boost { namespace numeric { namespace ublas {
template<typename MATRIX> MATRIX expm_pad(const MATRIX &H, const int p = 6) {
typedef typename MATRIX::value_type value_type;
typedef typename MATRIX::size_type size_type;
typedef double real_value_type; // Correct me. Need to modify.
assert(H.size1() == H.size2());
const size_type n = H.size1();
const identity_matrix<value_type> I(n);
matrix<value_type> U(n,n),H2(n,n),P(n,n),Q(n,n);
real_value_type norm = 0.0;
// Calcuate Pade coefficients (1-based instead of 0-based as in the c vector)
vector<real_value_type> c(p+2);
c(1)=1;
- for(size_type i = 1; i <= p; ++i)
+ for(size_type i = 1; i <= (unsigned) p; ++i)
c(i+1) = c(i) * ((p + 1.0 - i)/(i * (2.0 * p + 1 - i)));
// Calcuate the infinty norm of H, which is defined as the largest row sum of a matrix
for(size_type i=0; i<n; ++i)
{
real_value_type temp = 0.0;
for(size_type j=0;j<n;j++)
temp += std::abs<real_value_type>(H(j,i)); // Correct me, if H is complex, can I use that abs?
norm = std::max<real_value_type>(norm, temp);
}
if (norm == 0.0)
{
boost::throw_exception(boost::numeric::ublas::bad_argument());
std::cerr<<"Error! Null input in the routine EXPM_PAD.\n";
exit(0);
}
// Scaling, seek s such that || H*2^(-s) || < 1/2, and set scale = 2^(-s)
int s = 0;
real_value_type scale = 1.0;
if(norm > 0.5) {
s = std::max<int>(0, static_cast<int>((log(norm) / log(2.0) + 2.0)));
scale /= static_cast<real_value_type>(std::pow(2.0, s));
U.assign(scale * H); // Here U is used as temp value due to that H is const
}
else
U.assign(H);
// Horner evaluation of the irreducible fraction, see the following ref above.
// Initialise P (numerator) and Q (denominator)
H2.assign( prod(U, U) );
Q.assign( c(p+1)*I );
P.assign( c(p)*I );
size_type odd = 1;
for( size_type k = p - 1; k > 0; --k)
{
if( odd == 1)
{
Q = ( prod(Q, H2) + c(k) * I );
}
else
{
P = ( prod(P, H2) + c(k) * I );
}
odd = 1 - odd;
}
if( odd == 1)
{
Q = ( prod(Q, U) );
Q -= P ;
//U.assign( -(I + 2*(Q\P)));
}
else
{
P = (prod(P, U));
Q -= P;
//U.assign( I + 2*(Q\P));
}
// In origine expokit package, they use lapack ZGESV to obtain inverse matrix,
// and in that ZGESV routine, it uses LU decomposition for obtaing inverse matrix.
// Since in ublas, there is no matrix inversion template, I simply use the build-in
// LU decompostion package in ublas, and back substitute by myself.
//
//////////////// Implement Matrix Inversion ///////////////////////
permutation_matrix<size_type> pm(n);
int res = lu_factorize(Q, pm);
if( res != 0)
{
std::cerr << "Error in the matrix inversion in template expm_pad.\n";
exit(0);
}
H2 = I; // H2 is not needed anymore, so it is temporary used as identity matrix for substituting.
lu_substitute(Q, pm, H2);
if( odd == 1)
U.assign( -(I + 2.0 * prod(H2, P)));
else
U.assign( I + 2.0 * prod(H2, P));
// Squaring
- for(size_t i = 0; i < s; ++i)
+ for(size_t i = 0; i < (unsigned) s; ++i)
{
U = (prod(U,U));
}
return U;
}
}}}
#endif