Page MenuHomeHEPForge
Files F8725617

No OneTemporary
Actions

View Options

diff --git a/Hadronization/Ropewalk.cc b/Hadronization/Ropewalk.cc
--- a/Hadronization/Ropewalk.cc
+++ b/Hadronization/Ropewalk.cc
@@ -1,1207 +1,1207 @@
// -*- 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,
Energy gaIn, double geIn, double ycIn, double dyIn, Length tsIn, Length dtIn,
bool lorentzIn, bool cylinderIn, 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), 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 ){
stringdipoles[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();
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++;
+ ++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,true) && dip1->recoil(deltap1,true)){
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());
}
dip2->recoil(deltap2);
dip1->recoil(deltap1);
pp1->setMomentum(p1new);
pp2->setMomentum(p2new);
}
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){
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*deltat*gAmplitude/R0/gExponent/sqrt(2*M_PI)*
exp(-dist*dist/2/gExponent/gExponent/R0/R0);
}
Energy dpx = gain*direction.x()/dist;
Energy dpy = gain*direction.y()/dist;
// 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;
}
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;
}
}
diff --git a/Hadronization/Ropewalk.h b/Hadronization/Ropewalk.h
--- a/Hadronization/Ropewalk.h
+++ b/Hadronization/Ropewalk.h
@@ -1,727 +1,733 @@
// -*- C++ -*-
#ifndef TheP8I_Ropewalk_H
#define TheP8I_Ropewalk_H
//
// This is the declaration of the Ropewalk class.
//
#include "TheP8I/Config/TheP8I.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/EventRecord/ColourSinglet.h"
#include "ThePEG/Utilities/UtilityBase.h"
+#include "ThePEG/Utilities/Throw.h"
#include "HistKeeper.h"
namespace std {
template <>
struct less<ThePEG::ColourSinglet *> :
public binary_function<ThePEG::ColourSinglet *,
ThePEG::ColourSinglet *, bool>
{
/**
* This is the function called when comparing two pointers to
* type_info.
*/
bool operator()(const ThePEG::ColourSinglet * x,
const ThePEG::ColourSinglet * y) const {
if ( x && y ) return x->partons() < y->partons();
return !y;
}
};
}
namespace TheP8I {
using namespace ThePEG;
/**
* Here is the documentation of the Ropewalk class.
*/
class Ropewalk {
public:
double lambdaBefore;
vector<double> mspec;
/**
* Forward declaration of a Dipole.
*/
struct Dipole;
/**
* Container class storing information of a dipole overlapping with another;
*/
struct OverlappingDipole {
/**
* Standard constructor.
*/
OverlappingDipole(const Dipole & d,
const LorentzRotation R, const Ropewalk * rw);
/**
* Check if ovelapping at specific rapidity.
*/
bool overlap(double y, LorentzPoint b1, Length R0) {
if ( y < min(ya, yc) || y > max(ya, yc) ) return false;
LorentzPoint b2 = ba + (bc - ba)*(y - ya)/(yc - ya);
return (b1 - b2).perp() <= R0;
}
LorentzPoint pointAtRap(double y){
return ba + (bc - ba)*(y - ya)/(yc - ya);
}
/**
* Pointer to the dipole.
*/
const Dipole * dipole;
/**
* The boosted rapidities.
*/
double yc, ya;
/**
* The propagated and boosted positions.
*/
LorentzPoint bc, ba;
/**
* Relative direction.
*/
int dir;
/**
* Initially estimated fraction of overlap.
*/
double yfrac;
};
/**
* Container class containing information about overlaps for a dipole.
* *** ATTENTION *** the meaning of pa and pc is reversed: pa is
* always the COLOURED and pc is always the ANTI-COLOURED
*/
struct Dipole {
struct DipoleException : public Exception {};
/**
* Constructor.
*/
Dipole(tcPPtr p1, tcPPtr p2, ColourSinglet & str, Energy m0In)
: pc(p2), pa(p1), string(&str), n(0), m(0), p(1), q(0), nb(0), broken(false), hadr(false), m0(m0In) {
if(pc->colourLine() != pa->antiColourLine())
swap(pc,pa);
LorentzMomentum pcm = pc->momentum();
LorentzMomentum pam = pa->momentum();
//if(pcm.x()/GeV == 0 && pcm.y()/GeV == 0 && pcm.z()/GeV == 0 && pcm.e()/GeV == 0)
//cout << pc->id() << " " << pa->id() << endl;
//cout << pam.x()/GeV << " " << pam.y()/GeV << " " << pam.z()/GeV << " " << pam.e()/GeV << endl;
+ if ( (pcm + pam).m2() < 0.000001*GeV2 )
+ Throw<DipoleException>()
+ << "Ropewalk found a minisculedipole. This should not happen, "
+ << "but hoping it is just a fluke and throwing away the event."
+ << Exception::eventerror;
R = Utilities::boostToCM(make_pair(&pcm,&pam) );
}
/**
* Return true if this dipole is breakable. Only dipoles between
* gluons (not belonging to other dipoles which have broken) and
* which has a mass above 2 GeV can break.
*/
bool breakable() const {
return !broken && pc->id() == ParticleID::g && pa->id() == ParticleID::g &&
!pc->decayed() && !pa->decayed() && s() > 4.0*GeV2;
}
/**
* Calculate the probability that this dipole should break. Taking
* into account the number of negative steps in the colour space
* and the fraction of overlapping dipoles which are able to
* break.
*/
double breakupProbability() const;
void clear(){
overlaps.clear();
n = 0, m = 0, p = 1, q = 0, nb = 0, broken = false, hadr = false;
}
// Propagate the particles and the excitations
void propagate(Length deltat, double deltay);
/**
* The minimal rapidity in the lab system of the excited dipole ends
*/
double minRapidity() const {
//double y1 = 0.5*log(epc->Pplus()/epc->Pminus());
//double y2 = 0.5*log(epa->Pplus()/epa->Pminus());
double y1 = epc->momentum().rapidity();
double y2 = epa->momentum().rapidity();
return min(y1,y2);
}
/**
* The maximal rapidity in the lab system
*/
double maxRapidity() const {
//double y1 = 0.5*log(epc->Pplus()/epc->Pminus());
//double y2 = 0.5*log(epa->Pplus()/epa->Pminus());
double y1 = epc->momentum().rapidity();
double y2 = epa->momentum().rapidity();
return max(y1,y2);
}
/*
* The b-coordinate (in lab frame) of a point
* between the propagated particles, at lab-rapidity yIn
*/
LorentzPoint pointAtRap(double yIn) const {
double ymin = minRapidity();
double ymax = maxRapidity();
LorentzPoint bMin, bMax;
if (pc->momentum().rapidity() == ymin )
bMin = pc->vertex(), bMax = pa->vertex();
else
bMin = pa->vertex(), bMax = pc->vertex();
// We assume that the rapidity fraction stays the same,
// as we cannot calculate yIn in dipole rest frame
return bMin + (bMax - bMin)*(yIn - ymin)/(ymax - ymin);
}
/**
* Recalculate overlaps assuming a position a fraction ry from
* the coloured parton.
*/
double reinit(double ry, Length R0, Energy m0);
/**
* Set and return the effective multiplet.
*/
void initMultiplet();
void initNewMultiplet();
void initWeightMultiplet();
/**
* Calculate the multiplicity given pp and qq;
*/
static double multiplicity(int pp, int qq) {
return ( pp < 0 || qq < 0 || pp + qq == 0 )? 0.0:
0.5*(pp + 1)*(qq + 1)*(pp + qq + 2);
}
/**
* Calculate the average kappa-enhancement.
*/
double kappaEnhancement() const {
// return 1.0;
return 0.25*(p + q + 3 - p*q/double(p + q));
}
/**
* Calculate the kappa-enhancement in the first break-up.
*/
double firstKappa() const {
return 0.25*(2 + 2*p + q);
}
/**
* Calculate boosted kappa-enhancement.
*/
double firstKappa(double alpha) const {
if (alpha > 0)
return alpha*firstKappa();
return firstKappa();
}
/*
LorentzMomentum shoveOverlappers(double dipFrac, double m2Had, Length R0, Lorentz5Momentum pH){
// We first find the b-space point at the breaking point
LorentzPoint thisPoint = ba + (bc - ba)*dipFrac;
double thisy = pa->momentum().rapidity() +
(pc->momentum().rapidity() - pa->momentum().rapidity())*dipFrac;
Energy mp = sqrt(pH.mass2());
Energy mtd = sqrt(s());
Energy mSumH = sqrt(m2Had)*GeV;
double deltay = mp/(mtd - (mSumH - dipFrac*mtd)/(1 - dipFrac) );
// We then go through all dipoles
// overlapping in rapidity
Energy ex = 0*GeV, ey = 0*GeV;
for(int i = 0, N = overlaps.size(); i < N; ++i){
OverlappingDipole& od = overlaps[i];
// Must overlap at the breaking point
if (thisy < min(od.ya, od.yc) ||
thisy > max(od.ya, od.yc))
continue;
LorentzPoint otherPoint = od.pointAtRap(thisy);
LorentzPoint direction = thisPoint - otherPoint;
// The distance between the two dipoles
Length dist = direction.perp();
Energy gain = 2*fabs(deltay)*exp(-dist*dist/4/R0/R0)*GeV;
ex+=gain*direction.x()/dist;
ey+=gain*direction.y()/dist;
od.dipole->recoil(dipFrac,gain,direction);
}
Lorentz5Momentum ret(ex, ey, 0*GeV, 0*GeV);
return ret+pH;
}
*/
/**
* Push a dipole
*/
/*void recoil(double dipFrac, Energy gain, LorentzPoint direction) const {
LorentzMomentum pPush(-gain*direction.x()/direction.perp(), -gain*direction.y()/direction.perp(), 0*GeV, 0*GeV);
(void) pPush;
// We go to the dipole rest frame
// and calculate the push on the ends
//Lorentz5Momentum pcPush = R*pPush + pc->momentum();
//Lorentz5Momentum paPush = R*pPush + pa->momentum();
// const_cast<ThePEG::Particle& >((*pc)).setMomentum(pcPush);
// const_cast<ThePEG::Particle& >((*pa)).setMomentum(paPush);
}
*/
/**
* return the momentum of this dipole
*/
LorentzMomentum momentum() const {
return pc->momentum() + pa->momentum();
}
/**
* Return the squared invariant mass of this dipole.
*/
Energy2 s() const {
return (pc->momentum() + pa->momentum()).m2();
}
/**
* Return the effective rapidityspan of this dipole.
*/
double yspan(Energy m0) const {
return s() > sqr(m0)? 0.5*log(s()/sqr(m0)): 0.0;
}
/**
* The coloured particle.
*/
mutable tcPPtr pc;
/**
* The anti-coloured particle.
*/
mutable tcPPtr pa;
/**
* The coloured, excited particle
*/
PPtr epc;
/**
* The anti-coloured, excited particle
*/
PPtr epa;
/*
* All excitations belonging to this dipole ordered in
* rapidity in lab system
*/
typedef map<double, PPtr> DipEx;
DipEx excitations;
void printExcitations() {
for(DipEx::iterator itr = excitations.begin(); itr != excitations.end(); ++itr){
cout << itr->first << " " << itr->second->rapidity() << " " << itr->second->momentum().perp()/GeV << endl;
}
}
/*
* Recoil the dipole from adding a gluon
* this function will not actually add the gluon, only the
* recoil, if possible, and return true or false if it was
*/
bool recoil(LorentzMomentum& pg, bool dummy = false){
// Keep track of direction
int sign = 1;
if(epc->rapidity() < epa->rapidity())
sign = -1;
// lc momenta after inserting the gluon
Energy pplus = epc->Pplus() + epa->Pplus() - pg.plus();
Energy pminus = epc->Pminus() + epa->Pminus() - pg.minus();
// The new lightcone momenta of the dipole ends
Energy ppa = 0.0*GeV;
Energy ppc = 0.0*GeV;
Energy pma = 0.0*GeV;
Energy pmc = 0.0*GeV;
Energy4 sqarg = sqr(pplus*pminus + epa->mt2() - epc->mt2()) -
4.0*epa->mt2()*pplus*pminus;
// Calculate the new momenta
if(sqarg >= 0.0*GeV*GeV*GeV*GeV){
if(sign > 0){
ppa = (pplus*pminus + epa->mt2() - epc->mt2() +
sqrt(sqarg))*0.5/pminus;
pma = epa->mt2()/ppa;
pmc = pminus - pma;
ppc = epc->mt2()/pmc;
}
else{
pma = (pplus*pminus + epa->mt2() - epc->mt2() +
sqrt(sqarg))*0.5/pplus;
ppa = epa->mt2()/pma;
ppc = pplus - ppa;
pmc = epc->mt2()/ppc;
}
LorentzMomentum shifta = LorentzMomentum(epa->momentum().x(),epa->momentum().y(),
0.5*(ppa - pma),0.5*(ppa + pma));
LorentzMomentum shiftc = LorentzMomentum(epc->momentum().x(), epc->momentum().y(),
0.5*(ppc - pmc), 0.5*(ppc + pmc));
if(shifta.t() / GeV < 0 || shiftc.t() / GeV < 0)
return false;
// Assign the new momenta
if(!dummy){
epa->setMomentum(shifta);
epc->setMomentum(shiftc);
}
return true;
}
// We could not shift the momenta
return false;
}
/*
* Insert an excitation provided that it is not already there
*/
void addExcitation(double ylab, PPtr ex){
auto ret = excitations.equal_range(ylab);
for(auto itr = ret.first; itr != ret.second; ++itr )
if(ex == itr->second){
return;
}
excitations.insert(make_pair(ylab,ex));
}
/*
* Add all excitations to the dipole
*/
void excitationsToString(bool gluonMass, double yCutOff, HistKeeper* hkPtr = NULL);
/*
* Extract the distance between all adjecant excitation pairs on a dipole
*/
multimap< double, pair<double, double> > getDistPairs();
/**
* Return the current step.
*/
Step & step() const;
/**
* The propagated and boosted positions.
*/
LorentzPoint bc, ba;
/**
* The Lorentz rotation taking us to this dipole's rest system
*/
LorentzRotation R;
/**
* The overlapping dipoles with the amount of overlap (negative
* means anti overlap).
*/
vector<OverlappingDipole> overlaps;
/**
* The string to which this Dipole belongs.
*/
ColourSinglet * string;
/**
* The summed parallel (n) and anti-parallel (m) overlaps from
* other dipoles.
*/
int n, m;
/**
* The multiplet.
*/
int p, q;
/**
* The number of baryons from junctions
*/
int nb;
/**
* Indicate that this dipole has been broken.
*/
mutable bool broken;
mutable bool hadr;
/*
* The m0 parameter
*/
Energy m0;
};
/**
* Convenient typedef
*/
typedef map<ColourSinglet *, vector<Dipole *> > DipoleMap;
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
Ropewalk(const vector<ColourSinglet> & singlets, Length width,
Energy scale, double jDC, bool throwaway = true, bool shoving = true, double rcIn = 5.0, Energy gaIn = 1.0*GeV, double geIn = 1.0, double ycIn = 4.0, double dyIn = 0.1, Length tsIn = 1.0*femtometer, Length dtIn = 1.0*femtometer, bool lorentzIn = true, bool cylinderIn = false, HistKeeper* hk = NULL, bool deb = false);
/**
* The destructor.
*/
virtual ~Ropewalk();
//@}
/**
* Return all ColourSinglet objects with their kappa enhancement.
*/
vector< pair<ColourSinglet,double> > getSinglets(double DeltaY = 0.0) const;
/**
* Calculate the average kappa-enhancement for the given colour
* singlet, optionally only taking into account dipoles the
* rapidity interval |y| < deltaY. Note that \a cs must point to a
* object in the stringdipoles map. @return -1 if cs is not found.
*/
double averageKappaEnhancement(ColourSinglet * cs, double deltaY = 0.0) const {
DipoleMap::const_iterator it = stringdipoles.find(cs);
if ( it == stringdipoles.end() ) return -1.0;
return averageKappaEnhancement(it, deltaY);
}
/**
* Calculate the average kappa-enhancement for the given iterator
* pointing into the stringdipoles map, optionally only taking
* into account dipoles the rapidity interval |y| < deltaY. Note
* that \a cs must point to a object in the stringdipoles
* map. @return -1 if something went wrong.
*/
double averageKappaEnhancement(DipoleMap::const_iterator it,
double deltaY = 0.0) const {
return averageKappaEnhancement(it->second, deltaY);
}
/**
* Calculate the average kappa-enhancement the given vector of
* dipoles, optionally only taking into account dipoles the
* rapidity interval |y| < deltaY. @return -1 if something went
* wrong.
*/
double averageKappaEnhancement(const vector<Dipole *> & dipoles,
double deltaY = 0.0) const;
/**
* Propagate a parton a finite time inthe lab-system.
*/
LorentzPoint propagate(const LorentzPoint & b,
const LorentzMomentum & p) const;
/**
/ Get number of juctions
*/
double getNb();
double getkW();
/**
/ Get m,n for all strings
*/
pair<double,double> getMN();
/**
/ Get lambda measure for event
*/
double lambdaSum(double cutoff = 0.1);
/**
* Calculate the rapidity of a Lorentz vector but limit the
* transverse mass to be above m0.
*/
double limitrap(const LorentzMomentum & p) const;
/**
* Sort the given vector of dipoles so that the particle \a p is
* first and the others follows in sequence. If \a p is a
* (anti-)triplet, it will be the pc(pa) of the first dipole, and
* the pc(pa) of the second will be the pa(pc) of the first and so
* on. If \a p is an octet, the direction will be as if it was a
* triplet.
*/
static void sort(vector<Dipole *> & dips, tcPPtr p);
protected:
/**
* Extract all dipoles.
*/
void getDipoles();
/**
* Calculate all overlaps and initialize multiplets in all dipoles.
*/
void calculateOverlaps(bool doPropagate);
/**
* Break dipoles (and strings)
*/
void doBreakups();
/**
* Shove the dipoles away from each other
*/
void shoveTheDipoles();
/**
* Return the current step.
*/
Step & step() const;
/**
* Make a copy of a ColourSinglet making sure all partons are final.
*/
ColourSinglet * cloneToFinal(const ColourSinglet & cs);
public:
DipoleMap::iterator begin() {
return stringdipoles.begin();
}
DipoleMap::iterator end() {
return stringdipoles.end();
}
/**
* Exception class.
*/
struct ColourException: public Exception {};
private:
/**
* The assumed radius of a string.
*/
Length R0;
/**
* The assumed mass for calculating rapidities
*/
Energy m0;
/*
* The junction colour fluctuation parameter
*/
double junctionDiquarkProb;
/*
* Cutoff distance where we no longer consider
* shoving
*/
Length rCutOff;
/*
* Amplitude af shoving function
*/
Energy gAmplitude;
/*
* Shoving function (inverse) multiplier
*/
double gExponent;
/*
* Cutoff in effective rapidity span when clustering
*/
double yCutOff;
/*
* Rapidity slicing for shoving
*/
double deltay;
/*
* Shoving time
*/
Length tshove;
/*
* Shoving time slicing
*/
Length deltat;
/**
* The vector of all Dipoles
*/
vector<Dipole> dipoles;
/**
* All Dipoles in all string
*/
DipoleMap stringdipoles;
/*
* Lorentz contraction of fields
*/
bool lorentz;
bool cylinderShove;
// For analysis
HistKeeper* hkPtr;
/**
* Debug flag.
*/
bool debug;
public:
mutable double strl0;
mutable double strl;
mutable double avh;
mutable double avp;
mutable double avq;
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
Ropewalk & operator=(const Ropewalk &);
double nb;
};
}
#endif /* TheP8I_Ropewalk_H */

File Metadata

Mime Type
text/x-diff
Expires
Tue, Jan 21, 2:09 AM (1 d, 14 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4243591
Default Alt Text
(60 KB)

Event Timeline