diff --git a/Hadronization/Ropewalk.cc b/Hadronization/Ropewalk.cc
--- a/Hadronization/Ropewalk.cc
+++ b/Hadronization/Ropewalk.cc
@@ -1,1259 +1,1261 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the Ropewalk class.
//
#include "Ropewalk.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/EventRecord/Step.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/CurrentGenerator.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "ThePEG/Utilities/Selector.h"
#include "ThePEG/Utilities/SimplePhaseSpace.h"
#include "ThePEG/Utilities/UtilityBase.h"
#include "ThePEG/Utilities/Throw.h"
#include "ThePEG/Handlers/StepHandler.h"
using namespace TheP8I;
Ropewalk::Ropewalk(const vector<ColourSinglet> & singlets, Length width,
Energy scale, double jDP, bool throwaway, bool shoving, double rcIn,
StringTension gaIn, double geIn, double ycIn, double dyIn, Length tsIn, Length dtIn,
bool lorentzIn, bool cylinderIn, Energy pTcutIn, double mstringyIn, bool longSoftIn,
HistKeeper* hk, bool deb)
: R0(width), m0(scale), junctionDiquarkProb(jDP), rCutOff(rcIn*width), gAmplitude(gaIn),
gExponent(geIn), yCutOff(ycIn), deltay(dyIn), tshove(tsIn), deltat(dtIn), lorentz(lorentzIn),
cylinderShove(cylinderIn), stringpTCut(pTcutIn), minStringy(mstringyIn), longSoft(longSoftIn),
hkPtr(hk), debug(deb), strl0(0.0), strl(0.0), avh(0.0), avp(0.0), avq(0.0) {
for ( int is = 0, Ns = singlets.size(); is < Ns; ++is ){
if(!longSoft)
stringdipoles[cloneToFinal(singlets[is])];
else{
tcPVector pVec = singlets[is].partons();
double maxrap = 0;
double minrap = 0;
Energy pTmax = 0*GeV;
for(int id = 0, Nd = int(pVec.size()); id < Nd; ++id){
double y = pVec[id]->rapidity();
Energy pT = pVec[id]->momentum().perp();
if(id == 0){
maxrap = y;
minrap = y;
pTmax = pT;
}
else{
if(pTmax < pT)
pTmax = pT;
if(maxrap < y)
maxrap = y;
if(minrap > y)
minrap = y;
}
}
if(hkPtr){
hkPtr->hist("stringlength")->fill(maxrap-minrap,hkPtr->currentWeight());
hkPtr->hist("stringptMax")->fill(pTmax/GeV,hkPtr->currentWeight());
if(maxrap - minrap > minStringy)
hkPtr->hist("stringpTMaxCut")->fill(pTmax/GeV,hkPtr->currentWeight());
}
if(maxrap-minrap > minStringy && pTmax < stringpTCut)
stringdipoles[cloneToFinal(singlets[is])];
else
spectatorDipoles[cloneToFinal(singlets[is])];
}
}
getDipoles();
calculateOverlaps(false);
if(throwaway)
doBreakups();
lambdaBefore = lambdaSum();
if ( debug ) CurrentGenerator::log() << singlets.size() << "s "
<< dipoles.size() << "d ";
if(hkPtr){
for(DipoleMap::iterator itr = stringdipoles.begin(); itr!=stringdipoles.end();++itr){
vector<Dipole*> dips = itr->second;
for(int id = 0, Nd = dips.size(); id < Nd; ++id){
Dipole* d = dips[id];
LorentzMomentum mpa = d->pa->momentum();
LorentzMomentum mpc = d->pc->momentum();
LorentzMomentum mpaRot = (d->R)*mpa;
LorentzMomentum mpcRot = (d->R)*mpc;
if(d->pa->id() == 21){
hk->hist("normalGluonPt")->fill(mpa.perp()/GeV,
hk->currentWeight());
hk->hist("normalGluonRestFramePt")->fill(
mpaRot.perp()/GeV,hk->currentWeight());
}
if(abs(d->pa->id()) < 10 ){
hk->hist("normalQuarkPt")->fill(mpa.perp()/GeV,
hk->currentWeight());
hk->hist("normalQuarkRestFramePt")->fill(
mpaRot.perp()/GeV,hk->currentWeight());
}
if(abs(d->pc->id()) < 10 ){
hk->hist("normalQuarkPt")->fill(mpc.perp()/GeV,
hk->currentWeight());
hk->hist("normalQuarkRestFramePt")->fill(
mpcRot.perp()/GeV,hk->currentWeight());
}
if(abs(d->pc->id()) > 100 ){
hk->hist("normalDiQuarkPt")->fill(mpc.perp()/GeV,
hk->currentWeight());
hk->hist("normalDiQuarkRestFramePt")->fill(
mpcRot.perp()/GeV,hk->currentWeight());
}
if(abs(d->pa->id()) > 100 ){
hk->hist("normalDiQuarkPt")->fill(mpa.perp()/GeV,
hk->currentWeight());
hk->hist("normalDiQuarkRestFramePt")->fill(
mpaRot.perp()/GeV,hk->currentWeight());
}
}
}
}
if(shoving){
shoveTheDipoles();
// update color singlets
DipoleMap newStringDipoles;
for(DipoleMap::iterator itr = stringdipoles.begin(); itr != stringdipoles.end(); ++itr){
delete itr->first;
// The old dipoles
vector<Dipole*> oldDipoles = itr->second;
tcParticleSet updatedParticles;
for(int id = 0, Nd = oldDipoles.size(); id < Nd; ++id){
Dipole* d = oldDipoles[id];
const ParticleVector children = d->pa->children();
if(children.empty())
updatedParticles.insert(d->pa);
for(int ic = 0, Nc = children.size(); ic < Nc; ++ic)
updatedParticles.insert(children[ic]);
if(d->pc->id() != 21){
const ParticleVector children = d->pc->children();
if(children.empty())
updatedParticles.insert(d->pc);
for(int ic = 0, Nc = children.size(); ic < Nc; ++ic)
updatedParticles.insert(children[ic]);
}
}
tcPPtr p = *updatedParticles.begin();
if(p->colourLine())
newStringDipoles[new ColourSinglet(p->colourLine(), updatedParticles)];
else if(p->antiColourLine())
newStringDipoles[new ColourSinglet(p->antiColourLine(), updatedParticles)];
else
Throw<ColourException>()
<< "Updating colour singlets after shoving failed."
<< "This is a serious error, please contact the authors."
<< Exception::abortnow;
}
stringdipoles = newStringDipoles;
getDipoles();
calculateOverlaps(false);
if(hkPtr){
for(DipoleMap::iterator itr = stringdipoles.begin(); itr!=stringdipoles.end();++itr){
vector<Dipole*> dips = itr->second;
for(int id = 0, Nd = dips.size(); id < Nd; ++id){
Dipole* d = dips[id];
LorentzMomentum mpa = d->pa->momentum();
LorentzMomentum mpc = d->pc->momentum();
LorentzMomentum mpaRot = (d->R)*mpa;
LorentzMomentum mpcRot = (d->R)*mpc;
if(d->pa->id() == 21){
hk->hist("excitedEndGluonPt")->fill(mpa.perp()/GeV,
hk->currentWeight());
hk->hist("excitedEndGluonRestFramePt")->fill(
mpaRot.perp()/GeV,hk->currentWeight());
}
if(abs(d->pa->id()) < 10 ){
hk->hist("excitedQuarkPt")->fill(mpa.perp()/GeV,
hk->currentWeight());
hk->hist("excitedQuarkRestFramePt")->fill(
mpaRot.perp()/GeV,hk->currentWeight());
}
if(abs(d->pc->id()) < 10 ){
hk->hist("excitedQuarkPt")->fill(mpc.perp()/GeV,
hk->currentWeight());
hk->hist("excitedQuarkRestFramePt")->fill(
mpcRot.perp()/GeV,hk->currentWeight());
}
if(abs(d->pc->id()) > 100 ){
hk->hist("excitedDiQuarkPt")->fill(mpc.perp()/GeV,
hk->currentWeight());
hk->hist("excitedDiQuarkRestFramePt")->fill(
mpcRot.perp()/GeV,hk->currentWeight());
}
if(abs(d->pa->id()) > 100 ){
hk->hist("excitedDiQuarkPt")->fill(mpa.perp()/GeV,
hk->currentWeight());
hk->hist("excitedDiQuarkRestFramePt")->fill(
mpaRot.perp()/GeV,hk->currentWeight());
}
}
}
}
}
}
Ropewalk::~Ropewalk() {
for ( DipoleMap::iterator it = stringdipoles.begin();
it != stringdipoles.end(); ++it ) delete it->first;
}
ColourSinglet * Ropewalk::cloneToFinal(const ColourSinglet & cs) {
tcParticleSet final;
for ( int i = 0, N = cs.partons().size(); i < N; ++i )
final.insert(cs.partons()[i]->final());
tcPPtr p = *final.begin();
if(hkPtr){
hkPtr->hist("nPartons")->fill(final.size(), hkPtr->currentWeight());
}
if ( p->colourLine() ) return new ColourSinglet(p->colourLine(), final);
if ( p->antiColourLine() ) return new ColourSinglet(p->antiColourLine(), final);
Throw<ColourException>()
<< "Cloning ColourSinglets failed in Ropewalk. "
<< "This is a serious error - please contact the authors."
<< Exception::abortnow;
return 0;
}
LorentzPoint Ropewalk::
propagate(const LorentzPoint & b, const LorentzMomentum & p) const {
return b + p/(p.e()*m0/hbarc);
// return b;
}
void Ropewalk::getDipoles() {
dipoles.clear();
// First extract all dipoles from all strings (pieces).
for ( DipoleMap::iterator it = stringdipoles.begin();
it != stringdipoles.end(); ++it ) {
ColourSinglet & string = *it->first;
for ( int isp = 1, Nsp = string.nPieces(); isp <= Nsp; ++isp ) {
const ColourSinglet::StringPiece & sp = string.piece(isp);
if ( sp.size() < 2 ) continue;
bool forward = sp[0]->colourLine() &&
sp[0]->colourLine() == sp[1]->antiColourLine();
for ( int i = 0, N = sp.size(); i < (N - 1); ++i ) {
if ( forward )
dipoles.push_back(Dipole(sp[i], sp[i + 1], string, m0));
else
dipoles.push_back(Dipole(sp[i + 1], sp[i], string, m0));
}
if ( !forward && sp.front()->colourLine() &&
sp.front()->colourLine() == sp.back()->antiColourLine() )
dipoles.push_back(Dipole(sp.front(), sp.back(), string, m0));
else if ( forward && sp.front()->antiColourLine() &&
sp.front()->antiColourLine() == sp.back()->colourLine() )
dipoles.push_back(Dipole(sp.back(), sp.front(), string, m0));
}
}
for ( int i = 0, N = dipoles.size(); i < N; ++i ) {
stringdipoles[dipoles[i].string].push_back(&dipoles[i]);
strl0 += dipoles[i].yspan(1.0*GeV);
}
}
double Ropewalk::limitrap(const LorentzMomentum & p) const {
if ( p.z() == ZERO ) return 0.0;
Energy mt = sqrt(max(p.perp2() + p.m2(), sqr(m0)));
double rap = log((p.t() + abs(p.z()))/mt);
return p.z() > ZERO? rap: -rap;
}
Ropewalk::OverlappingDipole::
OverlappingDipole(const Dipole & d, const LorentzRotation R, const Ropewalk * rw)
: dipole(&d), dir(1) {
// Get the boost to other dipole's rest frame and calculate
// rapidities.
bc = R*rw->propagate(d.pc->vertex(), d.pc->momentum());
ba = R*rw->propagate(d.pa->vertex(), d.pa->momentum());
yc = rw->limitrap(R*d.pc->momentum());
ya = rw->limitrap(R*d.pa->momentum());
if ( yc < ya ) dir = -1;
}
void Ropewalk::calculateOverlaps(bool doPropagate) {
// Then calculate all overlaps between pairs of dipoles.
for ( int i1 = 0, Nd = dipoles.size(); i1 < Nd; ++i1 ) {
Dipole & d1 = dipoles[i1];
d1.clear();
if ( d1.s() <= sqr(m0) ) continue;
// Get the boost to dipoles rest frame and calculate rapidities.
LorentzMomentum pc1 = d1.pc->momentum();
LorentzMomentum pa1 = d1.pa->momentum();
//LorentzRotation R = Utilities::boostToCM(make_pair(&pc1, &pa1));
if(doPropagate){
d1.bc = d1.R*propagate(d1.pc->vertex(), d1.pc->momentum());
d1.ba = d1.R*propagate(d1.pa->vertex(), d1.pa->momentum());
}
else{
d1.bc = d1.pc->vertex();
d1.ba = d1.pa->vertex();
}
double yc1 = limitrap(pc1);
double ya1 = limitrap(pa1);
double n = 0.0;
double m = 0.0;
if ( yc1 <= ya1 ) continue; // don't care about too small strings.
for ( int i2 = 0; i2 < Nd; ++i2 ) {
if ( i1 == i2 ) continue;
// Boost to the rest frame of dipole 1.
if ( dipoles[i2].s() <= sqr(m0) ) continue;
OverlappingDipole od(dipoles[i2], d1.R, this);
// Ignore if not overlapping in rapidity.
if ( min(od.yc, od.ya) > yc1 || max(od.yc, od.ya) < ya1 || od.yc == od.ya )
continue;
// Calculate rapidity overlap.
double yomax = min(max(od.yc, od.ya), yc1);
double yomin = max(min(od.yc, od.ya), ya1);
// Sample a random point in the rapidity overlap region and get
// positions in space and check overlap.
double yfrac = (yomax - yomin)/(yc1 - ya1);
double y = UseRandom::rnd(yomin, yomax);
if ( !od.overlap(y, d1.ba + (d1.bc - d1.ba)*(y - ya1)/(yc1 - ya1), R0) )
yfrac = 0.0;
// Sum overlaps.
if ( od.dir > 0 ) {
n += yfrac;
} else {
m += yfrac;
}
od.yfrac = yfrac*od.dir;
d1.overlaps.push_back(od);
}
d1.n = int(n + UseRandom::rnd());
d1.m = int(m + UseRandom::rnd());
d1.initWeightMultiplet();
avp += d1.p*d1.yspan(1.0*GeV);
avq += d1.q*d1.yspan(1.0*GeV);
}
}
double Ropewalk::Dipole::reinit(double ry, Length R0, Energy m0) {
// First calculate the rapidity in this restframe.
double y = yspan(m0)*(ry - 0.5);
// Then get the position in space of the string at this rapidity.
LorentzPoint b = ba + (bc - ba)*ry;
p = 1;
q = 0;
// Now go through all overlapping dipoles.
for ( int i = 0, N = overlaps.size(); i < N; ++i ) {
if ( overlaps[i].dipole->broken ) continue;
if( overlaps[i].dipole->hadr ) continue;
if ( overlaps[i].overlap(y, b, R0) ) {
if ( overlaps[i].dir > 0 ) ++p;
else ++q;
}
}
return kappaEnhancement();
}
void Ropewalk::Dipole::initMultiplet() {
typedef pair<int,int> Plet;
int ns = 0;
int ms = 0;
while ( ns < n || ms < m ) {
Plet diff;
double octetveto = double(ns + ms + 1 - p - q)/double(ns + ms + 1);
if ( double(n - ns) > UseRandom::rnd()*double(n + m - ns - ms) ) {
Selector<Plet> sel;
sel.insert(multiplicity(p + 1, q), Plet(1, 0));
sel.insert(multiplicity(p, q - 1)*octetveto, Plet(0, -1));
sel.insert(multiplicity(p - 1, q + 1), Plet(-1, 1));
diff = sel[UseRandom::rnd()];
++ns;
} else {
Selector<Plet> sel;
sel.insert(multiplicity(p, q + 1), Plet(0, 1));
sel.insert(multiplicity(p - 1, q)*octetveto, Plet(-1, 0));
sel.insert(multiplicity(p + 1, q - 1), Plet(1, -1));
diff = sel[UseRandom::rnd()];
++ms;
}
if ( diff.first == -diff.second ) nb++;
p += diff.first;
q += diff.second;
}
/* *** ATTENTION *** maybe better if Christian implements this */
}
void Ropewalk::Dipole::initNewMultiplet() {
typedef pair<int,int> Plet;
int ns = 0;
int ms = 0;
for ( int i = 0, N = overlaps.size(); i < N; ++i ) {
if ( abs(overlaps[i].yfrac) < UseRandom::rnd() ) continue;
Plet diff(0,0);
double octetveto = double(ns + ms + 1 - p - q)/double(ns + ms + 1);
if ( overlaps[i].dir > 0 ) {
Selector<Plet> sel;
sel.insert(multiplicity(p + 1, q), Plet(1, 0));
sel.insert(multiplicity(p, q - 1)*octetveto, Plet(0, -1));
sel.insert(multiplicity(p - 1, q + 1), Plet(-1, 1));
diff = sel[UseRandom::rnd()];
++ns;
} else {
Selector<Plet> sel;
sel.insert(multiplicity(p, q + 1), Plet(0, 1));
sel.insert(multiplicity(p - 1, q)*octetveto, Plet(-1, 0));
sel.insert(multiplicity(p + 1, q - 1), Plet(1, -1));
diff = sel[UseRandom::rnd()];
++ms;
}
if ( diff.first == -diff.second ) nb++;
p += diff.first;
q += diff.second;
}
n = ns;
m = ms;
/* *** ATTENTION *** maybe better if Christian implements this */
}
void Ropewalk::Dipole::initWeightMultiplet() {
typedef pair<int,int> Plet;
int ns = 0;
int ms = 0;
double mab = sqrt(s())/GeV;
for ( int i = 0, N = overlaps.size(); i < N; ++i ) {
if ( abs(overlaps[i].yfrac) < UseRandom::rnd() ) continue;
const Dipole * od = overlaps[i].dipole;
double mcd = sqrt(od->s())/GeV;
double w8 = 1.0/sqr(mab*mcd); //The weight for 3 + 3bar -> 8
double w6 = w8; // The weight for 3 + 3 -> 6
Plet diff(0,0);
double octetveto = double(ns + ms + 1 - p - q)/double(ns + ms + 1);
if ( overlaps[i].dir > 0 ) {
double mac = (pc->momentum() + od->pc->momentum()).m()/GeV;
double mbd = (pa->momentum() + od->pa->momentum()).m()/GeV;
double w1 = octetveto/sqr(mac*mbd); // The weight for 3 + 3bar -> 1
double w3 = 1.0/(mab*mcd*mac*mbd); //The weight for 3 + 3 -> 3bar
Selector<Plet> sel;
sel.insert(multiplicity(p + 1, q) +
multiplicity(p, q - 1)*w8/(w8 + w1) +
multiplicity(p - 1, q + 1)*w6/(w6 + w3), Plet(1, 0));
sel.insert(multiplicity(p, q - 1)*w1/(w1 + w8), Plet(0, -1));
sel.insert(multiplicity(p - 1, q + 1)*w3/(w3 + w6), Plet(-1, 1));
diff = sel[UseRandom::rnd()];
++ns;
} else {
double mac = (pa->momentum() + od->pc->momentum()).m()/GeV;
double mbd = (pc->momentum() + od->pa->momentum()).m()/GeV;
if ( mac*mbd <= 0.0 ) {
++q; ++ms;
continue;
}
double w1 = octetveto/sqr(mac*mbd); // The weight for 3 + 3bar -> 1
double w3 = 1.0/(mab*mcd*mac*mbd); //The weight for 3 + 3 -> 3bar
Selector<Plet> sel;
sel.insert(multiplicity(p, q + 1) +
multiplicity(p - 1, q)*w8/(w8 + w1) +
multiplicity(p + 1, q - 1)*w6/(w6 + w3), Plet(0, 1));
sel.insert(multiplicity(p - 1, q)*w1/(w1 + w8), Plet(-1, 0));
sel.insert(multiplicity(p + 1, q - 1)*w3/(w3 + w6), Plet(1, -1));
diff = sel[UseRandom::rnd()];
++ms;
}
if ( diff.first == -diff.second ) nb++;
p += diff.first;
q += diff.second;
}
n = ns;
m = ms;
/* *** ATTENTION *** maybe better if Christian implements this */
}
double Ropewalk::getNb(){
double nba = 0;
for ( int i = 0, N = dipoles.size(); i < N; ++i ){
nba += double(dipoles[i].nb);
//cout << dipoles[i].n << " " << dipoles[i].m << " " << dipoles[i].p << " " << dipoles[i].q << endl;
//cout << double(dipoles[i].n + dipoles[i].m + 1 - dipoles[i].p - dipoles[i].q) << endl;
//cout << " --- " << endl;
//nba += double(dipoles[i].n + dipoles[i].m + 1 - dipoles[i].p - dipoles[i].q);
}
// cout << nba << endl;
return nba;
}
double Ropewalk::getkW(){
double kw = 0;
for ( int i = 0, N = dipoles.size(); i < N; ++i )
kw += dipoles[i].kappaEnhancement()*log(dipoles[i].s()/sqr(m0));
return kw;
}
pair<double,double> Ropewalk::getMN(){
pair<double,double> ret = make_pair<double,double>(0,0);
for (int i = 0, N = dipoles.size(); i < N; ++i){
ret.first += dipoles[i].m;
ret.second += dipoles[i].n;
}
return ret;
}
double Ropewalk::lambdaSum(double cutoff){
double ret = 0;
for (int i = 0, N = dipoles.size(); i < N; ++i)
ret += dipoles[i].s() > cutoff*GeV2 ? log(dipoles[i].s()/sqr(m0)) : 0;
return ret;
}
double Ropewalk::Dipole::breakupProbability() const {
if ( !breakable() || n + m <= 0.0 || n + m + 1 == p + q ) return 0.0;
double sum = 0.0;
for (int j = 0, N = overlaps.size(); j < N; ++j )
if ( overlaps[j].dipole->breakable() )
sum += abs(overlaps[j].yfrac);
if ( sum <= 0.0 ) return 1.0;
return double(n + m + 1 - p - q)/(sum + 1.0);
}
Step & Ropewalk::step() const {
return *(CurrentGenerator()->currentStepHandler()->newStep());
}
Step & Ropewalk::Dipole::step() const {
return *(CurrentGenerator()->currentStepHandler()->newStep());
}
// Helper class to describe a pair of excitations
struct Exc {
// The constructor
Exc(double yIn, Energy gmassIn, PPtr p1In, PPtr p2In, Ropewalk::Dipole* dip1In, Ropewalk::Dipole* dip2In) : y(yIn), gluonMass(gmassIn), pp1(p1In), pp2(p2In), dip1(dip1In), dip2(dip2In) {
// Set a sensible starting vertex
pp1->setVertex(dip1->pointAtRap(y));
pp2->setVertex(dip2->pointAtRap(y));
// Set a pointer to the excitation in the dipole
dip1->addExcitation(y, pp1);
dip2->addExcitation(y, pp2);
}
// Give the excitation a kick in the x and y direction,
void shove(Energy dpx, Energy dpy, HistKeeper* hkPtr){
// The old momenta
LorentzMomentum p2 = pp2->momentum();
LorentzMomentum p1 = pp1->momentum();
// The new momenta, after the shove
Energy mt2new = sqrt(sqr(p2.x() - dpx) + sqr(p2.y() - dpy) + sqr(gluonMass));
Energy e2new = mt2new*cosh(y);
Energy p2znew = mt2new*sinh(y);
LorentzMomentum p2new(p2.x() - dpx, p2.y() - dpy, p2znew, e2new);
Energy mt1new = sqrt(sqr(p1.x() + dpx) + sqr(p1.y() + dpy) + sqr(gluonMass));
Energy e1new = mt1new*cosh(y);
Energy p1znew = mt1new*sinh(y);
LorentzMomentum p1new(p1.x() + dpx, p1.y() + dpy, p1znew, e1new);
// The differences between the two
LorentzMomentum deltap1 = p1new - p1;
LorentzMomentum deltap2 = p2new - p2;
// We now want to check if we can add these
// two excitations to the dipoles
if( dip2->recoil(deltap2) ) {
if ( dip1->recoil(deltap1) ) {
pp1->setMomentum(p1new);
pp2->setMomentum(p2new);
if(hkPtr){
LorentzMomentum ex2Rot = (dip2->R)*deltap2;
LorentzMomentum ex1Rot = (dip1->R)*deltap1;
hkPtr->hist("shoveInRestFrame")->fill(ex2Rot.perp()/GeV,hkPtr->currentWeight());
hkPtr->hist("shoveInRestFrame")->fill(ex1Rot.perp()/GeV,hkPtr->currentWeight());
}
} else {
dip2->recoil(-deltap2);
if(hkPtr){
LorentzMomentum ex2Rot = (dip2->R)*deltap2;
LorentzMomentum ex1Rot = (dip1->R)*deltap1;
hkPtr->hist("rejectedShoveInRestFrame")->fill(ex2Rot.perp()/GeV,hkPtr->currentWeight());
hkPtr->hist("rejectedShoveInRestFrame")->fill(ex1Rot.perp()/GeV,hkPtr->currentWeight());
}
}
}
else if(hkPtr){
LorentzMomentum ex2Rot = (dip2->R)*deltap2;
LorentzMomentum ex1Rot = (dip1->R)*deltap1;
hkPtr->hist("rejectedShoveInRestFrame")->fill(ex2Rot.perp()/GeV,hkPtr->currentWeight());
hkPtr->hist("rejectedShoveInRestFrame")->fill(ex1Rot.perp()/GeV,hkPtr->currentWeight());
}
}
LorentzPoint direction(){
// This gives the push direction
return (pp2->vertex()-pp1->vertex());
}
double y;
Energy gluonMass;
PPtr pp1, pp2;
Ropewalk::Dipole* dip1;
Ropewalk::Dipole* dip2;
};
void Ropewalk::shoveTheDipoles() {
typedef multimap<double, Ropewalk::Dipole *> RMap;
RMap rapDipoles;
// Order the dipoles in maxrapidity and insert new, excited dipole ends
double ymin = 0;
double ymax = 0;
for(int i = 0, N = dipoles.size(); i < N ;++i){
if(dipoles[i].pc->decayed()){
// Set epc to the decay product of pc here
dipoles[i].epc = dipoles[i].pc->children()[0];
}
else{
// Make a new particle
dipoles[i].epc = CurrentGenerator()->getParticle(dipoles[i].pc->id());
dipoles[i].epc->setMomentum(dipoles[i].pc->momentum());
dipoles[i].epc->setVertex(dipoles[i].pc->vertex());
step().addDecayProduct(dipoles[i].pc,dipoles[i].epc, false);
}
if(dipoles[i].pa->decayed()){
// Set epa to the decay product of pa here
dipoles[i].epa = dipoles[i].pa->children()[0];
}
else{
dipoles[i].epa = CurrentGenerator()->getParticle(dipoles[i].pa->id());
dipoles[i].epa->setMomentum(dipoles[i].pa->momentum());
dipoles[i].epa->setVertex(dipoles[i].pa->vertex());
step().addDecayProduct(dipoles[i].pa,dipoles[i].epa, false);
}
// order the dipoles in max rapidity
rapDipoles.insert(make_pair(dipoles[i].maxRapidity(), &dipoles[i]));
// find maximal and minimal rapidity to sample
if(dipoles[i].minRapidity() < ymin)
ymin = dipoles[i].minRapidity();
if(dipoles[i].maxRapidity() > ymax)
ymax = dipoles[i].maxRapidity();
}
vector<double> rapidities;
// cout << rCutOff/femtometer << " " << gAmplitude/GeV << " " << gExponent << " " << yCutOff << " " << deltay << " " << tshove/femtometer << " " << deltat/femtometer << endl;
// We can let the gluon go off-shell between the pushes,
// and on-shell again towards the end
// Not fully implemented yet! Do not toggle this on!
bool gluonMass = false;
// Maybe do better sampling
for(double y = ymin; y < ymax; y+=deltay)
rapidities.push_back(y);
// For each value of ySample, we should have
// an excitation pair per dipole we want to treat
typedef map<double, vector<Exc> > ExMap;
ExMap exPairs;
Energy gMass = 0.001*keV;
for(int i = 0, N = rapidities.size(); i < N; ++i){
// find dipoles that are sampled here
// and store them temporarily
double ySample = rapidities[i];
vector<Ropewalk::Dipole*> tmp;
for(RMap::iterator rItr = rapDipoles.lower_bound(ySample);
rItr != rapDipoles.end(); ++rItr){
if(rItr->second->minRapidity() < ySample)
tmp.push_back(rItr->second);
}
// Construct the excitation particles, one for each sampled dipole
vector<PPtr> eParticles;
// dipoles to erase
vector<int> eraseDipoles;
for(int j = 0, M = tmp.size(); j < M; ++j){
// Test if it is possible to add the small excitation to
// the dipole
Energy pz = gMass*sinh(ySample);
Energy e = gMass*cosh(ySample);
LorentzMomentum ex = LorentzMomentum(0.0*GeV,
0.0*GeV,pz,e);
if(tmp[j]->recoil(ex,true)){
// Add the excitation
tmp[j]->recoil(ex);
eParticles.push_back(CurrentGenerator()->getParticle(21));
eParticles.back()->setMomentum(ex);
}
else{
// if it was not possible, remove the dipole from the list
// but NOT inside the loop
eraseDipoles.push_back(j);
}
}
// Erase dipoles which could not bear an excitation
for(int j = 0, M = eraseDipoles.size(); j < M; ++j){
tmp.erase(tmp.begin() + (eraseDipoles[j]-j) );
}
// Construct all pairs of possible excitations in this slice
exPairs[ySample] = vector<Exc>();
for(int j = 0, M = tmp.size(); j < M; ++j)
for(int k = j + 1; k < M; ++k){
// Both dipoles should span the full rapidity interval
//if(tmp[j]->maxRapidity() >= ySample+deltay/2.0 && tmp[k]->maxRapidity() >= ySample+deltay/2.0
// && tmp[j]->minRapidity() <= ySample-deltay/2.0
// && tmp[k]->minRapidity() <= ySample-deltay/2.0){
exPairs[ySample].push_back(Exc(ySample,gMass,eParticles[j],eParticles[k],tmp[j],tmp[k]));
//}
}
}
// For all times
int ii = 0;
for(Length t = 0*femtometer; t < tshove; t+=deltat){
++ii;
// For all slices
for(ExMap::iterator slItr = exPairs.begin(); slItr != exPairs.end(); ++slItr)
// For all excitation pairs
for(int i = 0, N = slItr->second.size(); i < N; ++i){
Exc& ep = slItr->second[i];
LorentzPoint direction = ep.direction();
Length dist = direction.perp();
if(hkPtr && ii < 6){
Histogram* hd = hkPtr->hist("dist"+to_string(ii));
hd->fill(dist/femtometer,hkPtr->currentWeight());
}
if(dist < rCutOff){
double gamma = 1.0;
if(lorentz){
LorentzMomentum pp = ep.dip1->momentum() + ep.dip2->momentum();
LorentzMomentum p1 = (ep.dip1->R)*pp;
LorentzMomentum p2 = (ep.dip2->R)*pp;
double gamma1 = sqrt(1 - sqr(p1.perp()/p1.e()));
double gamma2 = sqrt(1 - sqr(p2.perp()/p2.e()));
gamma = min(gamma1,gamma2);
}
Energy gain = 0*GeV;
if(cylinderShove){
+ Length R = R0*gExponent;
+ gain = gamma*deltay*deltat*t*gAmplitude*dist/R/R*exp(-dist*dist/2/R/R);
/*if(dist > R0*gExponent){
Length R = R0*gExponent;
Area area = 2*R*R*acos(dist/2/R) - 0.5*dist*sqrt(4*R*R - dist*dist);
Energy inner = gAmplitude*area/M_PI/R/R;
gain = gamma*(sqrt((1*GeV + inner)*(1*GeV + inner)) - 1*GeV);
}
else{
continue;
}
*/
}
else{
gain = gamma*deltay*(t+deltat)*deltat*gAmplitude/R0/gExponent/sqrt(2*M_PI)*
exp(-dist*dist/2/gExponent/gExponent/R0/R0);
}
Energy dpx = dist > ZERO ? gain*direction.x()/dist: 0.0*GeV;
Energy dpy = dist > ZERO ? gain*direction.y()/dist: 0.0*GeV;
// Print stuff
if(hkPtr && ii < 6){
Histogram* h = hkPtr->hist("shove"+to_string(ii));
h->fill(gain/MeV,hkPtr->currentWeight());
}
ep.shove(dpx,dpy,hkPtr);
}
}
// Propagate the dipoles
for(int i = 0, Nd = dipoles.size(); i < Nd; ++i)
dipoles[i].propagate(deltat, deltay);
}
// Now add the excitations to the dipoles
for(int i = 0, Nd = dipoles.size(); i < Nd; ++i){
if(dipoles[i].excitations.size() > 0)
dipoles[i].excitationsToString(gluonMass, yCutOff, hkPtr);
else{
dipoles[i].epc->antiColourConnect(dipoles[i].epa);
}
}
}
multimap< double, pair<double, double> > Ropewalk::Dipole::getDistPairs(){
// Fill a container with the particle pair rapidities in lab system,
// ordered by their spans
multimap< double, pair<double, double> > distPairs;
DipEx::iterator iplusone = excitations.begin();
++iplusone;
for(DipEx::iterator itr = excitations.begin(); itr != excitations.end(); ++itr,++iplusone){
if(itr == excitations.begin()){
// First -- special
LorentzMomentum p1 = (pa->momentum().rapidity() == minRapidity() ?
pa->momentum() : pc->momentum());
Energy2 s = (p1 + itr->second->momentum()).m2();
double ysep = (s > sqr(m0) ? 0.5*log(s/sqr(m0)) : 0.0);
distPairs.insert(make_pair(ysep, make_pair(minRapidity(),itr->first)) );
}
if(iplusone == excitations.end()){
// Last -- special
LorentzMomentum p1 = (pa->momentum().rapidity() == maxRapidity() ?
pa->momentum() : pc->momentum());
Energy2 s = (p1 + itr->second->momentum()).m2();
double ysep = (s > sqr(m0) ? 0.5*log(s/sqr(m0)) : 0.0);
distPairs.insert(make_pair(ysep, make_pair(itr->first,maxRapidity())) );
}
else{
Energy2 s = (itr->second->momentum() + iplusone->second->momentum()).m2();
double ysep = (s > sqr(m0) ? 0.5*log(s/sqr(m0)) : 0.0);
distPairs.insert(make_pair(ysep, make_pair(itr->first,iplusone->first)));
}
}
return distPairs;
}
void Ropewalk::Dipole::excitationsToString(bool gluonMass, double yCutOff, HistKeeper* hkPtr){
if(gluonMass){
cout << "Oops! Not implemented yet!" << endl;
return;
}
// Clustering should always be turned on
// only turn it off for debugging purposes.
// Change the parameters DeltaY or YCutOff instead
bool clustering = true;
if(clustering){
// Clustering
//
// For all excitation neighbors in the lab-system,
// we calculate their effective rapidity span
// and keep only those above yCutOff
multimap<double, pair<double, double> > distPairs = getDistPairs();
if(hkPtr){
for(auto itr = distPairs.begin(); itr != distPairs.end(); ++itr)
hkPtr->hist("yDistBefore")->fill(itr->first,hkPtr->currentWeight());
}
if(distPairs.begin() != distPairs.end())
while(distPairs.begin()->first < yCutOff){
// There is only one excitation on the string
// The ends will share the momentum provided
// that both distances are small
multimap<double, pair<double, double> >::iterator itr = distPairs.begin();
if(distPairs.size() == 2 && distPairs.rbegin()->first < yCutOff){
DipEx::iterator exP = excitations.find(itr->second.second);
if(exP == excitations.end())
exP = excitations.find(itr->second.first);
LorentzMomentum panew = epa->momentum() +
0.5*exP->second->momentum();
LorentzMomentum pcnew = epc->momentum() +
0.5*exP->second->momentum();
epa->setMomentum(panew);
epc->setMomentum(pcnew);
excitations.erase(exP);
}
// If it is the first piece, put the momentum on the dipole end
else if(itr->second.first == minRapidity()){
LorentzMomentum pnew = (epa->momentum().rapidity() == minRapidity() ?
epa->momentum() : epc->momentum());
pnew+=excitations[itr->second.second]->momentum();
if(epa->momentum().rapidity() == minRapidity())
epa->setMomentum(pnew);
else
epc->setMomentum(pnew);
// Erase the excitation
excitations.erase(itr->second.second);
}
// If it is the last piece ditto...
else if(itr->second.second == maxRapidity()){
LorentzMomentum pnew = (epa->momentum().rapidity() == maxRapidity() ?
epa->momentum() : epc->momentum());
pnew+=excitations[itr->second.first]->momentum();
if(epa->momentum().rapidity() == maxRapidity())
epa->setMomentum(pnew);
else
epc->setMomentum(pnew);
// Erase the excitation
excitations.erase(itr->second.first);
}
else{
LorentzMomentum pnew = excitations[itr->second.first]->momentum() +
excitations[itr->second.second]->momentum();
// update the first
excitations[itr->second.first]->setMomentum(pnew);
// erase the second
excitations.erase(itr->second.second);
}
// We can surely do better than recalculating everything
// each time - lets see if it kills us
if(excitations.empty())
break;
distPairs = getDistPairs();
}
if(hkPtr){
hkPtr->hist("nExcitations")->fill(excitations.size(),hkPtr->currentWeight());
for(auto itr = distPairs.begin(); itr != distPairs.end(); ++itr)
hkPtr->hist("yDistAfter")->fill(itr->first,hkPtr->currentWeight());
}
//
}
// Go through the remaining excitations and sanity check them
for(DipEx::iterator itr = excitations.begin(); itr != excitations.end(); ++itr){
// In the dipole rest frame, we should have pz == 0
//LorentzMomentum ppmom = R*itr->second->momentum();
//LorentzMomentum ppn(ppmom.x(),ppmom.y(),0*GeV,ppmom.t());
//itr->second->setMomentum(R.inverse()*ppn);
LorentzMomentum ppm = itr->second->momentum();
if(fabs(ppm.x()/GeV) < 1e-5 && fabs(ppm.y()/GeV) < 1e-5 && fabs(ppm.z()/GeV) < 1e-5 )
excitations.erase(itr);
else if(hkPtr){
hkPtr->hist("excitationPt")->fill(ppm.perp()/GeV,hkPtr->currentWeight());
LorentzMomentum ppr = R*ppm;
hkPtr->hist("excitationRestFramePt")->fill(ppr.perp()/GeV,hkPtr->currentWeight());
}
}
// If we dont have any excitations, we are done
if(excitations.empty()){
epa->colourConnect(epc);
return;
}
// Flag if we start from the colored parton
bool color = false;
// Link the first excitation to the dipole end with the smallest rapidity
if(minRapidity() == epc->momentum().rapidity()){
epa->colourConnect(excitations.rbegin()->second);
epc->antiColourConnect(excitations.begin()->second);
//excitations.begin()->second->colourConnect(epc);
//excitations.rbegin()->second->antiColourConnect(epa);
step().addDecayProduct(pc,excitations.begin()->second,false);
color = true;
}
else{
epa->colourConnect(excitations.begin()->second);
epc->antiColourConnect(excitations.rbegin()->second);
//excitations.begin()->second->antiColourConnect(epa);
//excitations.rbegin()->second->colourConnect(epc);
step().addDecayProduct(pa,excitations.begin()->second,false);
}
// If we have only one excitation, we are quickly done
if(excitations.size() == 1)
return;
// We have more...
DipEx::iterator iplusone = excitations.begin();
++iplusone;
for(DipEx::iterator itr = excitations.begin(); iplusone != excitations.end(); ++itr,++iplusone){
// Are we halfway through?
bool halfway = double(distance(excitations.begin(),iplusone))/double(excitations.size()) > 0.5;
if(color){
itr->second->antiColourConnect(iplusone->second);
//iplusone->second->colourConnect(itr->second);
step().addDecayProduct((halfway ? pa : pc), iplusone->second,false);
}
else{
itr->second->colourConnect(iplusone->second);
//iplusone->second->antiColourConnect(itr->second);
step().addDecayProduct((halfway ? pc : pa), iplusone->second,false);
}
}
}
void Ropewalk::Dipole::propagate(Length deltat, double deltay){
const Lorentz5Momentum& pcm = epc->momentum();
const Lorentz5Momentum& pam = epa->momentum();
if(pcm.e()/GeV == 0 || pam.e()/GeV == 0)
Throw<ColourException>()
<< "Tried to propagate a dipole end with energy 0; this should never happen."
<< "This is a serious error - please contact the authors."
<< Exception::abortnow;
// Propagate in dipole rest frame
bc = R*(bc + deltat * pcm/(pcm.e()));
ba = R*(ba + deltat * pam/(pam.e()));
// Propagate in the lab frame
LorentzPoint bclab = epc->vertex() + deltat*pcm/pcm.e();
LorentzPoint balab = epa->vertex() + deltat*pam/pam.e();
epc->setVertex(bclab);
epa->setVertex(balab);
//const_cast<ThePEG::Particle& >((*pc)).setVertex(bclab);
//const_cast<ThePEG::Particle& >((*pa)).setVertex(balab);
for(DipEx::iterator eItr = excitations.begin(); eItr != excitations.end(); ++eItr){
const Lorentz5Momentum& em = eItr->second->momentum();
// Propagate excitations only in lab frame
// only if it has been pushed
if(em.e()/GeV > 0){
eItr->second->setVertex(eItr->second->vertex() + deltat*em/(em.perp() + m0*deltay*deltay));
}
else{
eItr->second->setVertex(pointAtRap(eItr->first));
}
}
}
void Ropewalk::doBreakups() {
typedef multimap<double, const Dipole *> BMap;
BMap breakups;
// First order alldipoles in increasing probability of breaking
for ( int i = 0, N = dipoles.size(); i < N; ++i )
breakups.insert(make_pair(dipoles[i].breakupProbability(), &dipoles[i]));
// Then start breaking in reverse order
mspec.clear();
for ( BMap::reverse_iterator it = breakups.rbegin();
it != breakups.rend(); ++it ) {
const Dipole & d = *it->second;
if ( d.breakupProbability() < UseRandom::rnd() ) continue;
mspec.push_back(d.s()/GeV2);
double hinv = 1.0 / d.kappaEnhancement();
// The default parameters should be set
// properly from the outside. Now just Pythia8 defaults
double rho = pow(0.215, hinv);
double x = pow(1.0,hinv);
double y = pow(0.027,hinv);
Selector<int> qsel;
qsel.clear();
if(UseRandom::rnd() < junctionDiquarkProb){
// Take just ordinary quarks
qsel.insert(1.0, ParticleID::ubar);
qsel.insert(1.0, ParticleID::dbar);
qsel.insert(rho, ParticleID::sbar);
}
else {
// Take diquarks
// TODO (Pythia 8.2) Set status code of these diquarks to 74
qsel.insert(1.0, ParticleID::ud_0);
qsel.insert(3.0*y, ParticleID::dd_1);
qsel.insert(3.0*y, ParticleID::ud_1);
qsel.insert(3.0*y, ParticleID::uu_1);
qsel.insert(rho*x, ParticleID::su_0);
qsel.insert(3*x*rho*y, ParticleID::su_1);
qsel.insert(3*y*x*x*rho*rho, ParticleID::ss_1);
qsel.insert(rho*x, ParticleID::sd_0);
qsel.insert(3*x*rho*y, ParticleID::sd_1);
}
int flav = qsel[UseRandom::rnd()];
PPtr dic = CurrentGenerator()->getParticle(flav);
PPtr dia = CurrentGenerator()->getParticle(-flav);
//cout << (dic->momentum() + dia->momentum()).m2() / GeV2 << " " << d.s() / GeV2 << endl;
if((dic->momentum() + dia->momentum()).m2() > d.s()) continue;
// Get boost to dipole rest frame and place the di-quarks there.
LorentzRotation R = Utilities::getBoostFromCM(make_pair(d.pc, d.pa));
SimplePhaseSpace::CMS(dic, dia, d.s(), 1.0, 0.0);
dia->transform(R);
dic->transform(R);
// Add the di-quarks as children to the decayed gluons and split
// the resulting string.
if ( !( step().addDecayProduct(d.pc, dic, true) &&
step().addDecayProduct(d.pa, dia, true) &&
step().addDecayProduct(d.pc, dia, false) &&
step().addDecayProduct(d.pa, dic, false)) )
Throw<ColourException>()
<< "Adding decay products failed in Ropewalk. "
<< "This is a serious error - please contact the authors."
<< Exception::abortnow;
tcParticleSet pset(d.string->partons().begin(), d.string->partons().end());
pset.erase(d.pc);
pset.erase(d.pa);
pset.insert(dic);
pset.insert(dia);
// Reset the neighboring dipoles to point to the new diquarks.
d.broken = true;
ColourSinglet * olds = d.string;
const vector<Dipole *> & oldd = stringdipoles[olds];
for ( int id = 0, Nd = oldd.size(); id < Nd; ++id ) {
if ( oldd[id]->pc == d.pa ) oldd[id]->pc = dia;
if ( oldd[id]->pa == d.pc ) oldd[id]->pa = dic;
}
// remove the old string and insert the new ones fixing the
// pointers from the dipoles.
vector<ColourSinglet> newstrings =
ColourSinglet::getSinglets(pset.begin(), pset.end());
for ( int is = 0, Ns = newstrings.size(); is < Ns; ++is ) {
ColourSinglet * news = new ColourSinglet(newstrings[is]);
tcParticleSet pset(news->partons().begin(), news->partons().end());
vector<Dipole *> & newd = stringdipoles[news];
for ( int id = 0, Nd = oldd.size(); id < Nd; ++id ) {
if ( !oldd[id]->broken && pset.find(oldd[id]->pc) != pset.end() ) {
newd.push_back(oldd[id]);
oldd[id]->string = news;
}
}
sort(newd, news->partons()[0]);
}
stringdipoles.erase(olds);
delete olds;
}
}
vector< pair<ColourSinglet,double> > Ropewalk::getSinglets(double DeltaY) const {
vector< pair<ColourSinglet,double> > ret;
ret.reserve(stringdipoles.size());
double toty = 0.0;
double totyh = 0.0;
double toth = 0.0;
double toth2 = 0.0;
double minh = 10.0;
double maxh = 0.0;
// Calculate the kappa-enhancement of the remaining strings.
for ( DipoleMap::const_iterator it = stringdipoles.begin();
it != stringdipoles.end(); ++it ) {
double sumyh = 0.0;
double sumy = 0.0;
for ( int id = 0, Nd = it->second.size(); id < Nd; ++id ) {
double y = it->second[id]->yspan(m0);
double h = it->second[id]->kappaEnhancement();
avh += h*it->second[id]->yspan(1.0*GeV);
strl += it->second[id]->yspan(1.0*GeV);
if ( DeltaY > 0.0 &&
abs(it->second[id]->pc->eta()) > DeltaY &&
abs(it->second[id]->pa->eta()) > DeltaY ) continue;
sumy += y;
sumyh += y*h;
}
toty += sumy;
totyh += sumyh;
if ( sumy > 0.0 ) {
double h = 0.1*int(10.0*sumyh/sumy + UseRandom::rnd());
ret.push_back(make_pair(*it->first, h));
toth += sumyh/sumy;
toth2 += sqr(sumyh/sumy);
minh = min(minh, sumyh/sumy);
maxh = max(maxh, sumyh/sumy);
} else {
ret.push_back(make_pair(*it->first, 1.0));
toth += 1.0;
toth2 += 1.0;
minh = min(minh, 1.0);
maxh = max(maxh, 1.0);
}
}
if ( debug )
CurrentGenerator::log() << ret.size() << "ns "
<< totyh/max(toty, 0.00001) << "hwa "
<< minh << "<h "
<< toth/ret.size() << "ha "
<< maxh << ">h "
<< sqrt(toth2/ret.size() - sqr(toth/ret.size())) << "hsd " << endl;
if ( avh > 0.0 ) avh /= strl;
if ( avp > 0.0 ) avp /= strl0;
if ( avq > 0.0 ) avq /= strl0;
return ret;
}
double Ropewalk::averageKappaEnhancement(const vector<Dipole *> & dipoles,
double DeltaY) const {
double sumyh = 0.0;
double sumy = 0.0;
for ( int id = 0, Nd = dipoles.size(); id < Nd; ++id ) {
dipoles[id]->reinit(UseRandom::rnd(), R0, m0);
double y = dipoles[id]->yspan(m0);
double h = dipoles[id]->firstKappa();
if ( DeltaY > 0.0 &&
abs(dipoles[id]->pc->eta()) > DeltaY &&
abs(dipoles[id]->pa->eta()) > DeltaY ) continue;
sumy += y;
sumyh += y*h;
}
if ( sumy <= 0.0 ) return 1.0;
return sumyh/sumy;
}
void Ropewalk::sort(vector<Dipole *> & dips, tcPPtr p) {
if ( p->antiColourLine() ) {
map<tcPPtr, Dipole *> dmap;
Dipole * current = 0;
for ( int i = 0, N = dips.size(); i < N; ++i ) {
if ( dips[i]->pa == p ) current = dips[i];
else dmap[dips[i]->pa] = dips[i];
}
dips.clear();
while ( current ) {
dips.push_back(current);
map<tcPPtr, Dipole *>::iterator it = dmap.find(current->pc);
if ( it == dmap.end() ) current = 0;
else {
current = it->second;
dmap.erase(it);
}
}
if ( !dmap.empty() )
Throw<ColourException>()
<< "Failed to sort dipoles in Ropewalk. "
<< "This is a serious error - please contact the authors."
<< Exception::abortnow;
} else {
map<tcPPtr, Dipole *> dmap;
Dipole * current = 0;
for ( int i = 0, N = dips.size(); i < N; ++i ) {
if ( dips[i]->pc == p ) current = dips[i];
else dmap[dips[i]->pc] = dips[i];
}
dips.clear();
while ( current ) {
dips.push_back(current);
map<tcPPtr, Dipole *>::iterator it = dmap.find(current->pa);
if ( it == dmap.end() ) current = 0;
else {
current = it->second;
dmap.erase(it);
}
}
if ( !dmap.empty() )
Throw<ColourException>()
<< "Failed to sort dipoles in Ropewalk. "
<< "This is a serious error - please contact the authors."
<< Exception::abortnow;
}
}