Page MenuHomeHEPForge
Files F10881212

No OneTemporary
Actions

View Options

diff --git a/Hadronization/Ropewalk.cc b/Hadronization/Ropewalk.cc
--- a/Hadronization/Ropewalk.cc
+++ b/Hadronization/Ropewalk.cc
@@ -1,311 +1,326 @@
// -*- 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/Handlers/StepHandler.h"
using namespace TheP8I;
Ropewalk::Ropewalk(const vector<ColourSinglet> & singlets, Length width,
Energy scale, bool deb)
- : R0(width), m0(scale), debug(deb) {
+ : R0(width), m0(scale), 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();
if ( debug ) CurrentGenerator::log() << singlets.size() << "s "
<< dipoles.size() << "d ";
calculateOverlaps();
doBreakups();
}
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 ( p->colourLine() ) return new ColourSinglet(p->colourLine(), final);
if ( p->antiColourLine() ) return new ColourSinglet(p->antiColourLine(), final);
cerr << "Oooops! cloning ColourSinglets failed." << endl;
return 0;
}
LorentzPoint Ropewalk::
propagate(const LorentzPoint & b, const LorentzMomentum & p) const {
return b + p/(p.e()*m0/hbarc);
+ // return b;
}
void Ropewalk::getDipoles() {
// 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));
else
dipoles.push_back(Dipole(sp[i + 1], sp[i], string));
}
if ( !forward && sp.front()->colourLine() &&
sp.front()->colourLine() == sp.back()->antiColourLine() )
dipoles.push_back(Dipole(sp.front(), sp.back(), string));
else if ( forward && sp.front()->antiColourLine() &&
sp.front()->antiColourLine() == sp.back()->colourLine() )
dipoles.push_back(Dipole(sp.back(), sp.front(), string));
}
}
- for ( int i = 0, N = dipoles.size(); i < N; ++i )
+ 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;
}
void Ropewalk::calculateOverlaps() {
// Then calculate all overlaps between pairs of dipoles.
for ( int i1 = 0, Nd = dipoles.size(); i1 < Nd; ++i1 ) {
Dipole & d1 = dipoles[i1];
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));
LorentzPoint bc1 = R*propagate(d1.pc->vertex(), d1.pc->momentum());
LorentzPoint ba1 = R*propagate(d1.pa->vertex(), d1.pa->momentum());
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.
Dipole & d2 = dipoles[i2];
if ( d2.s() <= sqr(m0) ) continue;
// Get the boost to dipoles rest frame and calculate rapidities.
LorentzMomentum pc2 = R*d2.pc->momentum();
LorentzMomentum pa2 = R*d2.pa->momentum();
LorentzPoint bc2 = R*propagate(d2.pc->vertex(), d2.pc->momentum());
LorentzPoint ba2 = R*propagate(d2.pa->vertex(), d2.pa->momentum());
double yc2 = limitrap(pc2);
double ya2 = limitrap(pa2);
// Ignore if not overlapping in rapidity.
if ( min(yc2, ya2) > yc1 || max(yc2, ya2) < ya1 || yc2 == ya2 ) continue;
// Calculate rapidity overlap.
double yomax = min(max(yc2, ya2), yc1);
double yomin = max(min(yc2, ya2), ya1);
// Sample a random point in the rapidity overlap region and get
// position in space.
double y = UseRandom::rnd(yomin, yomax);
LorentzPoint b1 = ba1 + (bc1 - ba1)*(y - ya1)/(yc1 - ya1);
LorentzPoint b2 = ba2 + (bc2 - ba2)*(y - ya2)/(yc2 - ya2);
// Skip if there is no overlap in transverse position.
if ( (b1 - b2).perp() > R0 ) continue;
// Register the overlap.
double yfrac = (yomax - yomin)/(yc1 - ya1);
if ( yc2 > ya2 ) {
d1.overlaps.push_back(make_pair(&d2, yfrac));
n += yfrac;
} else {
d1.overlaps.push_back(make_pair(&d2, -yfrac));
m += yfrac;
}
}
d1.n = int(n + UseRandom::rnd());
d1.m = int(m + UseRandom::rnd());
d1.initMultiplet();
+ avp += d1.p*d1.yspan(1.0*GeV);
+ avq += d1.q*d1.yspan(1.0*GeV);
}
}
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;
}
p += diff.first;
q += diff.second;
}
/* *** ATTENTION *** maybe better if Christian implements this */
}
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].first->breakable() )
sum += abs(overlaps[j].second);
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());
}
void Ropewalk::doBreakups() {
// We need to be able to select diquarks.
Selector<int> disel;
disel.insert(1.0, ParticleID::ud_0);
disel.insert(3.0, ParticleID::dd_1);
disel.insert(3.0, ParticleID::ud_1);
disel.insert(3.0, ParticleID::uu_1);
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
for ( BMap::reverse_iterator it = breakups.rbegin();
it != breakups.rend(); ++it ) {
const Dipole & d = *it->second;
if ( d.breakupProbability() < UseRandom::rnd() ) continue;
// Select di-quarks
int flav = disel[UseRandom::rnd()];
PPtr dic = CurrentGenerator()->getParticle(flav);
PPtr dia = CurrentGenerator()->getParticle(-flav);
// 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) ) )
cerr << "Oooops! adding decay products failed" << endl;
tcParticleSet pset(d.string->partons().begin(), d.string->partons().end());
pset.erase(d.pc);
pset.erase(d.pa);
pset.insert(dic);
pset.insert(dia);
// remove the old stringand insert the new ones fixing the
// pointers from the dipolmes.
ColourSinglet * olds = d.string;
const vector<Dipole *> & oldd = stringdipoles[olds];
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 ( pset.find(oldd[id]->pc) != pset.end() ) {
newd.push_back(oldd[id]);
oldd[id]->string = news;
}
}
}
stringdipoles.erase(olds);
delete olds;
}
}
-vector< pair<ColourSinglet,double> > Ropewalk::getSingets() const {
+vector< pair<ColourSinglet,double> > Ropewalk::getSingets(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;
}
diff --git a/Hadronization/Ropewalk.h b/Hadronization/Ropewalk.h
--- a/Hadronization/Ropewalk.h
+++ b/Hadronization/Ropewalk.h
@@ -1,226 +1,235 @@
// -*- C++ -*-
#ifndef TheP8I_Ropewalk_H
#define TheP8I_Ropewalk_H
//
// This is the declaration of the Ropewalk class.
//
#include "TheP8I/Config/Pythia8Interface.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/EventRecord/ColourSinglet.h"
namespace TheP8I {
using namespace ThePEG;
/**
* Here is the documentation of the Ropewalk class.
*/
class Ropewalk {
public:
/**
* Container class containing information about individual dipole overlaps.
*/
struct Dipole {
/**
* Constructor.
*/
Dipole(tcPPtr p1, tcPPtr p2, ColourSinglet & str)
: pc(p1), pa(p2), string(&str), n(0), m(0), p(1), q(0) {}
/**
* 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 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;
/**
* Set and return the effective multiplet.
*/
void initMultiplet();
/**
* 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 kappa-enhancement.
*/
double kappaEnhancement() const {
return 0.25*(p + q + 3 - p*q/double(p + q));
}
/**
* 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.
*/
tcPPtr pc;
/**
* The anti-coloured particle.
*/
tcPPtr pa;
/**
* The overlapping dipoles with the amount of overlap (negative
* means anti overlap).
*/
vector< pair<const Dipole *, double> > 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;
};
/**
* 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, bool deb = false);
/**
* The destructor.
*/
virtual ~Ropewalk();
//@}
/**
* Return all ColourSinglet objects with their kappa enhancement.
*/
- vector< pair<ColourSinglet,double> > getSingets() const;
+ vector< pair<ColourSinglet,double> > getSingets(double DeltaY = 0.0) const;
protected:
/**
* Extract all dipoles.
*/
void getDipoles();
/**
* Calculate all overlaps and initialize multiplets in all dipoles.
*/
void calculateOverlaps();
/**
* Break dipoles (and strings)
*/
void doBreakups();
/**
* Propagate a parton a finite time inthe lab-system.
*/
LorentzPoint propagate(const LorentzPoint & b,
const LorentzMomentum & p) const;
/**
* Calculate the rapidity of a Lorentz vector but limit the
* transverse mass to be above m0.
*/
double limitrap(const LorentzMomentum & p) const;
/**
* Return the current step.
*/
Step & step() const;
/**
* Make a copy of a ColourSinglet making sure all partons are final.
*/
static ColourSinglet * cloneToFinal(const ColourSinglet & cs);
private:
/**
* The assumed radius of a string.
*/
Length R0;
/**
* The assumed mass for calculating rapidities
*/
Energy m0;
/**
* The vector of all Dipoles
*/
vector<Dipole> dipoles;
/**
* All Dipoles in all string
*/
DipoleMap stringdipoles;
/**
* 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 &);
};
}
#endif /* TheP8I_Ropewalk_H */

File Metadata

Mime Type
text/x-diff
Expires
Sat, May 3, 6:01 AM (1 d, 3 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4969089
Default Alt Text
(16 KB)

Event Timeline