Page MenuHomeHEPForge
Files F9514739

No OneTemporary
Actions

View Options

diff --git a/Hadronization/Makefile.am b/Hadronization/Makefile.am
--- a/Hadronization/Makefile.am
+++ b/Hadronization/Makefile.am
@@ -1,21 +1,21 @@
if PYTHIA8WORKS
-mySOURCES = StringFragmentation.cc TheP8EventShapes.cc StringPipe.cc OverlapPythiaHandler.cc RandomAverageHandler.cc RandomHandler.cc VanillaPythiaHandler.cc
+mySOURCES = StringFragmentation.cc TheP8EventShapes.cc StringPipe.cc OverlapPythiaHandler.cc RandomAverageHandler.cc RandomHandler.cc VanillaPythiaHandler.cc Ropewalk.cc
-DOCFILES = StringFragmentation.h TheP8EventShapes.h OverlapEnhancementHandler.h StringPipe.h OverlapPythiaHandler.h RandomAverageHandler.h RandomHandler.h VanillaPythiaHandler.h
+DOCFILES = StringFragmentation.h TheP8EventShapes.h OverlapEnhancementHandler.h StringPipe.h OverlapPythiaHandler.h RandomAverageHandler.h RandomHandler.h VanillaPythiaHandler.h Ropewalk.h
AUTOGENERATED = StringFragmentation-var.h StringFragmentation-init.h \
StringFragmentation-input.h StringFragmentation-output.h \
StringFragmentation-interfaces.h
INCLUDEFILES = $(DOCFILES) $(AUTOGENERATED)
AM_CPPFLAGS += -I/usr/local/new/include/
noinst_LTLIBRARIES = libTheP8IString.la
libTheP8IString_la_SOURCES = $(mySOURCES) $(INCLUDEFILES)
endif
include $(top_srcdir)/Config/Makefile.aminclude
diff --git a/Hadronization/Ropewalk.cc b/Hadronization/Ropewalk.cc
new file mode 100644
--- /dev/null
+++ b/Hadronization/Ropewalk.cc
@@ -0,0 +1,257 @@
+// -*- 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)
+ : R0(width), m0(scale) {
+ for ( int is = 0, Ns = singlets.size(); is < Ns; ++is )
+ stringdipoles[new ColourSinglet(singlets[is])];
+ getDipoles();
+ calculateOverlaps();
+ doBreakups();
+}
+
+Ropewalk::~Ropewalk() {
+ for ( DipoleMap::iterator it = stringdipoles.begin();
+ it != stringdipoles.end(); ++it ) delete it->first;
+}
+
+LorentzPoint Ropewalk::
+propagate(const LorentzPoint & b, const LorentzMomentum & p) const {
+ return b + p/(p.e()*max(m0, p.m())/hbarc);
+}
+
+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);
+ for ( int i = 0, N = sp.size(); i < (N - 1); ++i ) {
+ if ( sp[i]->colourLine() &&
+ sp[i]->colourLine() == sp[i+1]->antiColourLine() )
+ dipoles.push_back(Dipole(sp[i], sp[i +1 ], string));
+ else if ( sp[i]->antiColourLine() &&
+ sp[i]->antiColourLine() == sp[i+1]->colourLine() )
+ dipoles.push_back(Dipole(sp[i + 1], sp[i], string));
+ }
+ if ( sp.front()->colourLine() &&
+ sp.front()->colourLine() == sp.back()->antiColourLine() )
+ dipoles.push_back(Dipole(sp.front(), sp.back(), string));
+ else if ( 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 )
+ stringdipoles[dipoles[i].string].push_back(&dipoles[i]);
+
+}
+
+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];
+ // Get the boost to dipoles rest frame and calculate rapidities.
+ LorentzMomentum pc1 = d1.pc->momentum();
+ LorentzMomentum pa1 = d1.pa->momentum();
+ LorentzRotation R = Utilities::getBoostToCM(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 = log((pc1.e() + pc1.z())/max(m0, pc1.m()));
+ double ya1 = log((pa1.e() - pa1.z())/max(m0, pa1.m()));
+ 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];
+ // 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 = log((pc2.e() + pc2.z())/max(m0, pc2.m()));
+ double ya2 = log((pa2.e() - pa2.z())/max(m0, pa2.m()));
+ // Ignore if not overlapping in rapidity.
+ if ( min(yc2, ya2) > yc1 || max(yc2, ya2) < ya1 ) 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();
+ }
+}
+
+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 + q);
+ if ( ns < n && ( ms == m || UseRandom::rndbool() ) ) {
+ 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 ( 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 dia = CurrentGenerator()->getParticle(flav);
+ PPtr dic = 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.
+ step().addDecayProduct(d.pc, dic, true);
+ step().addDecayProduct(d.pa, dia, true);
+ step().addDecayProduct(d.pc, dia, false);
+ step().addDecayProduct(d.pa, dic, false);
+ 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 dipoles.
+ 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(newstrings[is].partons().begin(),
+ newstrings[is].partons().end());
+ for ( int id = 0, Nd = oldd.size(); id < Nd; ++id ) {
+ if ( pset.find(oldd[id]->pc) != pset.end() ) {
+ stringdipoles[news].push_back(oldd[id]);
+ oldd[id]->string = news;
+ }
+ }
+ }
+ stringdipoles.erase(olds);
+ delete olds;
+ }
+}
+
+vector< pair<ColourSinglet,double> > Ropewalk::getSingets() const {
+ vector< pair<ColourSinglet,double> > ret;
+ ret.reserve(stringdipoles.size());
+ // 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();
+ sumy += y;
+ sumyh += y*h;
+ }
+ ret.push_back(make_pair(*it->first, sumy/sumyh));
+ }
+ return ret;
+}
+
+
+
diff --git a/Hadronization/Ropewalk.h b/Hadronization/Ropewalk.h
new file mode 100644
--- /dev/null
+++ b/Hadronization/Ropewalk.h
@@ -0,0 +1,210 @@
+// -*- 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 ablle 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 0.5*log(s()/sqr(m0));
+ }
+
+ /**
+ * 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);
+
+ /**
+ * The destructor.
+ */
+ virtual ~Ropewalk();
+ //@}
+
+
+ /**
+ * Return all ColourSinglet objects with their kappa enhancement.
+ */
+ vector< pair<ColourSinglet,double> > getSingets() 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;
+ /**
+ * Return the current step.
+ */
+ Step & step() const;
+
+
+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;
+
+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
Mon, Feb 24, 6:39 AM (1 d, 14 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4494584
Default Alt Text
(14 KB)

Event Timeline