diff --git a/patch_added_newColourReconnection.diff b/patch_added_newColourReconnection.diff new file mode 100644 --- /dev/null +++ b/patch_added_newColourReconnection.diff @@ -0,0 +1,1680 @@ +diff -r 1781af1e3113 -r 5d17715d8d26 Hadronization/ColourReconnector.cc +--- a/Hadronization/ColourReconnector.cc Mon Jul 18 09:45:26 2022 +0100 ++++ b/Hadronization/ColourReconnector.cc Wed Jul 20 18:37:09 2022 +0200 +@@ -51,9 +51,18 @@ + + // do the colour reconnection + switch (_algorithm) { +- case 0: _doRecoPlain(clusters); break; +- case 1: _doRecoStatistical(clusters); break; +- case 2: _doRecoBaryonic(clusters); break; ++ case 0: ++ _doRecoPlain(clusters); ++ break; ++ case 1: ++ _doRecoStatistical(clusters); ++ break; ++ case 2: ++ _doRecoBaryonic(clusters); ++ break; ++ case 3: ++ _doRecoBaryonicMesonic(clusters); ++ break; + } + } + +@@ -68,7 +77,40 @@ + return sum; + } + ++double ColourReconnector::_displacement(tcPPtr p, tcPPtr q) const { ++ double deltaRap = (p->rapidity() - q->rapidity()); ++ double deltaPhi = (p->momentum().phi() - q->momentum().phi()); + ++ return sqrt(deltaRap * deltaRap + deltaPhi * deltaPhi); ++} ++ ++ ++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; ++ double meanPhi=(phi1 + phi2 + phi3)/3.0; ++ double delR; ++ ++ delR = sqrt( (rap1-meanRap)*(rap1-meanRap) + (phi1-meanPhi)*(phi1-meanPhi) ); ++ delR += sqrt( (rap2-meanRap)*(rap2-meanRap) + (phi2-meanPhi)*(phi2-meanPhi) ); ++ delR += sqrt( (rap3-meanRap)*(rap3-meanRap) + (phi3-meanPhi)*(phi3-meanPhi) ); ++ 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 & P) const { + assert (P.size() == cv.size()); +@@ -202,13 +244,13 @@ + } + + namespace { +- inline bool hasDiquark(CluVecIt cit) { +- for(int i = 0; i<(*cit)->numComponents(); i++) { ++inline bool hasDiquark(CluVecIt cit) { ++ for (int i = 0; i<(*cit)->numComponents(); i++) { + if (DiquarkMatcher::Check(*((*cit)->particle(i)->dataPtr()))) + return true; + } + return false; +- } ++} + } + + +@@ -221,7 +263,7 @@ + + // try to avoid systematic errors by randomising the reconnection order + long (*p_irnd)(long) = UseRandom::irnd; +- random_shuffle( newcv.begin(), newcv.end(), p_irnd ); ++ random_shuffle(newcv.begin(), newcv.end(), p_irnd); + + // iterate over all clusters + for (CluVecIt cit = newcv.begin(); cit != newcv.end(); ++cit) { +@@ -245,16 +287,16 @@ + baryonic1, baryonic2); + + // skip this cluster if no possible reconnection partner can be found +- if ( !isBaryonicCandidate && candidate==cit ) ++ if (!isBaryonicCandidate && candidate==cit) + continue; + +- if ( isBaryonicCandidate +- && UseRandom::rnd() < _precoBaryonic ) { ++ if (isBaryonicCandidate ++ && UseRandom::rnd() < _precoBaryonic) { + deleted.push_back(*baryonic2); + + // Function that does the reconnection from 3 -> 2 clusters + ClusterPtr b1, b2; +- _makeBaryonicClusters(*cit,*baryonic1,*baryonic2, b1, b2); ++ _makeBaryonicClusters(*cit, *baryonic1, *baryonic2, b1, b2); + + *cit = b1; + *baryonic1 = b2; +@@ -263,9 +305,9 @@ + } + + // Normal 2->2 Colour reconnection +- if ( !isBaryonicCandidate +- && UseRandom::rnd() < _preco ) { +- auto reconnected = _reconnectBaryonic(*cit, *candidate); ++ if (!isBaryonicCandidate ++ && UseRandom::rnd() < _preco) { ++ auto reconnected = _reconnect(*cit, *candidate); + *cit = reconnected.first; + *candidate = reconnected.second; + } +@@ -274,18 +316,254 @@ + // 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() ) ++ for (const auto & cluster : newcv) ++ if (find(deleted.begin(), ++ deleted.end(), cluster) == deleted.end()) + clustervector.push_back(cluster); + +- swap(cv,clustervector); ++ swap(cv, clustervector); + + + + } + + ++bool ColourReconnector::_clustersFarApart( const std::vector & clu ) const { ++ int Ncl=clu.size(); ++ assert(Ncl<=3); ++ if (Ncl==1) { ++ return false; ++ } else if (Ncl==2) { ++ // TODO: keep turned off all until things are more clear ++ // BUG: space-time difference compared to _maxDistance ++ // if (((*clu[0])->vertex()-(*clu[1])->vertex()).m() >_maxDistance) return true; ++ // Causal selection if desired ++ // if (((*clu[0])->vertex()-(*clu[1])->vertex()).m() >ZERO) return true; ++ } else if (Ncl==3) { ++ // TODO: keep turned off all until things are more clear ++ // BUG: space-time difference compared to _maxDistance ++ // if (((*clu[0])->vertex()-(*clu[1])->vertex()).m()> _maxDistance) return true; ++ // Causal selection if desired ++ // if (((*clu[0])->vertex()-(*clu[1])->vertex()).m()> ZERO) return true; ++ // if (((*clu[1])->vertex()-(*clu[2])->vertex()).m()> _maxDistance) return true; ++ // if (((*clu[1])->vertex()-(*clu[2])->vertex()).m()> ZERO) return true; ++ // if (((*clu[0])->vertex()-(*clu[2])->vertex()).m()> _maxDistance) return true; ++ // if (((*clu[0])->vertex()-(*clu[2])->vertex()).m()> ZERO) return true; ++ } ++ ++ return false; ++} ++ ++ ++ ++void ColourReconnector::_doReco2BeamClusters(ClusterVector & cv) const { ++ // try other option ++ PPtr p1Di=(cv[0])->colParticle(); ++ PPtr p2Di=(cv[1])->colParticle(); ++ ++ PPtr p1Q=(cv[0])->antiColParticle(); ++ PPtr p2Q=(cv[1])->antiColParticle(); ++ ++ double min_dist=_displacement(p1Di,p1Q)+_displacement(p2Di,p2Q); ++ ++ if ((_displacement(p1Di,p2Q)+_displacement(p1Di,p2Q)) 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 { + +@@ -314,6 +592,84 @@ + } + + ++void ColourReconnector::_findbestreconnectionoption(std::vector & 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; ++} ++ ++ + CluVecIt ColourReconnector::_findPartnerBaryonic( + CluVecIt cl, ClusterVector & cv, + bool & baryonicCand, +@@ -349,7 +705,6 @@ + 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; +@@ -364,8 +719,10 @@ + continue; + + // stop it putting far apart clusters together +- if ( ( (**cl).vertex()-(**cit).vertex() ).m() >_maxDistance ) +- continue; ++ // BUG: space-time difference compared to _maxDistance ++ // if (((**cl).vertex()-(**cit).vertex()).m() >_maxDistance) ++ // Causal selection if desired ++ // if (((**cl).vertex()-(**cit).vertex()).m() >ZERO) continue; + + const bool Colour8 = + _isColour8( (*cl)->colParticle(), (*cit)->antiColParticle() ) +@@ -412,7 +769,7 @@ + // 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){ ++ if (maxsum != 0) { + baryonic2 = baryonic1; + baryonic1 = cit; + bcand = true; +@@ -431,7 +788,7 @@ + + } + +- if(bcand == true){ ++ if (bcand == true) { + baryonicCand = true; + } + +@@ -447,7 +804,7 @@ + for (CluVecIt cit=cv.begin(); cit != cv.end(); ++cit) { + + // don't even look at original cluster +- if(cit==cl) continue; ++ if (cit==cl) continue; + + // don't allow colour octet clusters + if ( _isColour8( (*cl)->colParticle(), +@@ -458,10 +815,13 @@ + } + + // stop it putting beam remnants together +- if((*cl)->isBeamCluster() && (*cit)->isBeamCluster()) continue; ++ if ((*cl)->isBeamCluster() && (*cit)->isBeamCluster()) continue; + + // stop it putting far apart clusters together +- if(((**cl).vertex()-(**cit).vertex()).m()>_maxDistance) continue; ++ // BUG: space-time difference compared to _maxDistance ++ // if (((**cl).vertex()-(**cit).vertex()).m()>_maxDistance) continue; ++ // Causal selection if desired ++ // if (((**cl).vertex()-(**cit).vertex()).m()>ZERO) continue; + + // momenta of the old clusters + Lorentz5Momentum p1 = (*cl)->colParticle()->momentum() + +@@ -493,7 +853,7 @@ + ClusterPtr &c1, ClusterPtr &c2, + ClusterPtr &c3, + ClusterPtr &newcluster1, +- ClusterPtr &newcluster2) const{ ++ ClusterPtr &newcluster2) const { + + //make sure they all have 2 components + assert(c1->numComponents()==2); +@@ -522,6 +882,105 @@ + } + + pair ++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 (newCluster1, newCluster2); ++} ++ ++ ++std::tuple ++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 (newCluster1, newCluster2, newCluster3); ++} ++ ++pair + ColourReconnector::_reconnect(ClusterPtr &c1, ClusterPtr &c2) const { + + // choose the other possibility to form two clusters from the given +@@ -530,7 +989,7 @@ + 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) { ++ 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; +@@ -538,50 +997,8 @@ + } + assert(c1_col>=0&&c2_col>=0&&c1_anti>=0&&c2_anti>=0); + +- ClusterPtr newCluster1 +- = new_ptr( Cluster( c1->colParticle(), c2->antiColParticle() ) ); +- +- 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); +- +- ClusterPtr newCluster2 +- = new_ptr( Cluster( c2->colParticle(), c1->antiColParticle() ) ); +- +- 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 (newCluster1, newCluster2); +-} +- +- +- +- +- +-pair +-ColourReconnector::_reconnectBaryonic(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); ++ c1->colParticle()->abandonChild(c1); ++ c2->antiColParticle()->abandonChild(c2); + + ClusterPtr newCluster1 + = new_ptr( Cluster( c1->colParticle(), c2->antiColParticle() ) ); +@@ -589,8 +1006,11 @@ + c1->colParticle()->addChild(newCluster1); + c2->antiColParticle()->addChild(newCluster1); + +- newCluster1->setVertex(0.5*( c1->colParticle()->vertex() + +- c2->antiColParticle()->vertex() )); ++ /* ++ * 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); +@@ -604,8 +1024,11 @@ + c1->antiColParticle()->addChild(newCluster2); + c2->colParticle()->addChild(newCluster2); + +- newCluster2->setVertex(0.5*( c2->colParticle()->vertex() + +- c1->antiColParticle()->vertex() )); ++ /* ++ * 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); +@@ -613,9 +1036,289 @@ + return pair (newCluster1, newCluster2); + } + ++std::tuple ++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 (newclusters[0], newclusters[1], newclusters[2]); ++} ++ ++ ++pair ++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 (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 & 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 ColourReconnector::_shuffle +- (const PVector & q, const PVector & aq, unsigned maxtries) const { ++(const PVector & q, const PVector & aq, unsigned maxtries) const { + + const size_t nclusters = q.size(); + assert (nclusters > 1); +@@ -628,7 +1331,9 @@ + do { + // find two different random integers in the range [0, nclusters) + i = UseRandom::irnd( nclusters ); +- do { j = UseRandom::irnd( nclusters ); } while (i == j); ++ 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] ) ; +@@ -666,8 +1371,7 @@ + if(p->hasColour() && q->hasAntiColour()) { + cline = p-> colourLine(); + aline = q->antiColourLine(); +- } +- else { ++ } else { + cline = q-> colourLine(); + aline = p->antiColourLine(); + } +@@ -703,20 +1407,57 @@ + + + void ColourReconnector::persistentOutput(PersistentOStream & os) const { +- os << _clreco << _preco << _precoBaryonic << _algorithm << _initTemp << _annealingFactor +- << _annealingSteps << _triesPerStepFactor << ounit(_maxDistance,femtometer) +- << _octetOption; ++ os ++ << _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 ++ << _debug ++ << _junctionMBCR ++ ; + } + + void ColourReconnector::persistentInput(PersistentIStream & is, int) { +- is >> _clreco >> _preco >> _precoBaryonic >> _algorithm >> _initTemp >> _annealingFactor +- >> _annealingSteps >> _triesPerStepFactor >> iunit(_maxDistance,femtometer) +- >> _octetOption; ++ is ++ >> _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 ++ >> _debug ++ >> _junctionMBCR ++ ; + } + + + void ColourReconnector::Init() { +- + static ClassDocumentation documentation + ("This class is responsible of the colour reconnection."); + +@@ -736,8 +1477,36 @@ + "Colour reconnections on", + 1); + ++ // Algorithm interface ++ static Switch 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 Parameter interfaceMtrpAnnealingFactor ++ ++ ++ // Statistical CR Parameters: ++ static Parameter interfaceMtrpAnnealingFactor + ("AnnealingFactor", + "The annealing factor is the ratio of the temperatures in two successive " + "temperature steps.", +@@ -766,50 +1535,75 @@ + false, false, Interface::limited); + + +- static Parameter interfaceRecoProb ++ ++ ++ // Plain and Baryonic CR Paramters ++ static Parameter interfaceRecoProb + ("ReconnectionProbability", +- "Probability that a found reconnection possibility is actually accepted", ++ "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 interfaceRecoProbBaryonic + ("ReconnectionProbabilityBaryonic", +- "Probability that a found reconnection possibility is actually accepted", ++ "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 Switch 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); ++ // BaryonicMesonic CR Paramters ++ static Parameter 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 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 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 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 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 interfaceMaxDistance ++ static Parameter 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 interfaceMesonToBaryonFactor ++ ("MesonToBaryonFactor", ++ "Factor for comparing mesonic clusters to baryonic clusters in the displacement if BaryonicMesonic CR model is chosen", ++ &ColourReconnector::_mesonToBaryonFactor, 2.0, 0.5, 3.0, ++ false, false, Interface::limited); ++ ++ ++ ++ // General Parameters and switches ++ static Parameter interfaceMaxDistance + ("MaxDistance", + "Maximum distance between the clusters at which to consider rearrangement" +- " to avoid colour reconneections of displaced vertices", ++ " 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 interfaceOctetTreatment ++ static Switch interfaceOctetTreatment + ("OctetTreatment", +- "Which octets are not allowed to be reconnected", ++ "Which octets are not allowed to be reconnected (used in all Algorithms)", + &ColourReconnector::_octetOption, 0, false, false); + static SwitchOption interfaceOctetTreatmentFinal + (interfaceOctetTreatment, +@@ -821,6 +1615,51 @@ + "All", + "Prevent for all octets", + 1); ++ static Switch 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 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 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 -r 1781af1e3113 -r 5d17715d8d26 Hadronization/ColourReconnector.h +--- a/Hadronization/ColourReconnector.h Mon Jul 18 09:45:26 2022 +0100 ++++ b/Hadronization/ColourReconnector.h Wed Jul 20 18:37:09 2022 +0200 +@@ -59,6 +59,43 @@ + Energy2 _clusterMassSum(const PVector & q, const PVector & aq) 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: ++ * = sum_i Y_i/3 ++ * <\phi> = sum_i \phi_i/3 ++ * \deltaR = sum_i sqrt( (Y_i - )^2 + (\phi_i - )^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 &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 _reconnect(ClusterPtr &c1, ClusterPtr &c2) const; + + /** +- * Reconnection method for baryonic reconenction model ++ * @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 _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 _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 _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 _reconnectBaryonic(ClusterPtr &c1, ClusterPtr &c2) const; ++ pair _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 & 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 &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 +@@ -156,6 +329,25 @@ + */ + 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 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 = picometer; ++ Length _maxDistance = femtometer; + + /** + * @return true, if the two partons are splitting products of the same +@@ -213,6 +456,11 @@ + */ + unsigned int _octetOption = 0; + ++ /** ++ * Option for colour reconnecting 2 Beam Clusters if no others are present ++ */ ++ int _cr2BeamClusters = 0; ++ + public: + + /** @name Functions used by the persistent I/O system. */ +diff -r 1781af1e3113 -r 5d17715d8d26 src/defaults/Hadronization.in +--- a/src/defaults/Hadronization.in Mon Jul 18 09:45:26 2022 +0100 ++++ b/src/defaults/Hadronization.in Wed Jul 20 18:37:09 2022 +0200 +@@ -41,11 +41,37 @@ + newdef LightClusterDecayer:HadronSelector HadronSelector + newdef ClusterDecayer:HadronSelector HadronSelector + ++# ColourReconnector Default Parameters + newdef ColourReconnector:ColourReconnection Yes +-newdef ColourReconnector:ReconnectionProbabilityBaryonic 0.7 ++newdef ColourReconnector:Algorithm Baryonic ++ ++# Statistical CR Parameters: ++newdef ColourReconnector:AnnealingFactor 0.9 ++newdef ColourReconnector:AnnealingSteps 50 ++newdef ColourReconnector:TriesPerStepFactor 5.0 ++newdef ColourReconnector:InitialTemperature 0.1 ++ ++# Plain and Baryonic CR Paramters + newdef ColourReconnector:ReconnectionProbability 0.95 +-newdef ColourReconnector:Algorithm Baryonic ++newdef ColourReconnector:ReconnectionProbabilityBaryonic 0.7 ++ ++# BaryonicMesonic and BaryonicMesonic CR Paramters ++newdef ColourReconnector:ReconnectionProbability3Mto3M 0.5 ++newdef ColourReconnector:ReconnectionProbability3MtoBBbar 0.5 ++newdef ColourReconnector:ReconnectionProbabilityBbarBto3M 0.5 ++newdef ColourReconnector:ReconnectionProbability2Bto2B 0.05 ++newdef ColourReconnector:ReconnectionProbabilityMBtoMB 0.5 ++newdef ColourReconnector:StepFactor 1.0 ++newdef ColourReconnector:MesonToBaryonFactor 1.333 ++ ++# General Parameters and switches ++newdef ColourReconnector:MaxDistance 1.0e50 + newdef ColourReconnector:OctetTreatment All ++newdef ColourReconnector:CR2BeamClusters No ++newdef ColourReconnector:Junction Yes ++# Debugging ++newdef ColourReconnector:Debug No ++ + + # Clustering parameters for light quarks + newdef ClusterFissioner:ClMaxLight 3.649