Page MenuHomeHEPForge
Files F10881264

No OneTemporary
Actions

View Options

diff --git a/Hadronization/Ropewalk.cc b/Hadronization/Ropewalk.cc
--- a/Hadronization/Ropewalk.cc
+++ b/Hadronization/Ropewalk.cc
@@ -1,265 +1,296 @@
// -*- 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) {
+ Energy scale, bool deb)
+ : R0(width), m0(scale), debug(deb) {
for ( int is = 0, Ns = singlets.size(); is < Ns; ++is )
stringdipoles[new ColourSinglet(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;
}
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);
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 )
stringdipoles[dipoles[i].string].push_back(&dipoles[i]);
}
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];
// 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 = 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];
// 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 ) 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 + 1);
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;
+ 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 dia = CurrentGenerator()->getParticle(flav);
- PPtr dic = CurrentGenerator()->getParticle(-flav);
+ 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.
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(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> > 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();
sumy += y;
sumyh += y*h;
}
- if ( sumyh > 0.0 ) ret.push_back(make_pair(*it->first, sumy/sumyh));
- else ret.push_back(make_pair(*it->first, 1.0));
+ toty += sumy;
+ totyh += sumyh;
+
+ if ( sumy > 0.0 ) {
+ ret.push_back(make_pair(*it->first, sumyh/sumy));
+ 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;
+
return ret;
+
}
diff --git a/Hadronization/Ropewalk.h b/Hadronization/Ropewalk.h
--- a/Hadronization/Ropewalk.h
+++ b/Hadronization/Ropewalk.h
@@ -1,217 +1,222 @@
// -*- 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
+ * 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 0.5*log(s()/sqr(m0));
+ return max(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);
+ Energy scale, bool deb = false);
/**
* 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;
/**
* 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;
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;
+
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 */
diff --git a/Hadronization/StringFragmentation.cc b/Hadronization/StringFragmentation.cc
--- a/Hadronization/StringFragmentation.cc
+++ b/Hadronization/StringFragmentation.cc
@@ -1,639 +1,639 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the StringFragmentation class.
//
#include "StringFragmentation.h"
#include "Ropewalk.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/PDT/RemnantDecayer.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/EventRecord/Step.h"
#include "ThePEG/Handlers/EventHandler.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DebugItem.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include <algorithm> // std::sort
using namespace TheP8I;
StringFragmentation::StringFragmentation()
: pythia(), overlapN(0), testPow(1.0), stringR0(0.5*femtometer), stringyCut(2.0), stringpTCut(6*GeV), fragmentationMass(0.134), baryonSuppression(2.), window(false), useThrustAxis(0), fragmentationScheme("none"), maxTries(2), doAnalysis(0), analysisPath(""),
#include "StringFragmentation-init.h"
{
}
StringFragmentation::~StringFragmentation() {
}
void StringFragmentation::handle(EventHandler & eh, const tPVector & tagged,
const Hint & h) {
++nev;
// Get the event
tPVector all = RemnantDecayer::decayRemnants(tagged, *newStep());
if(_eventShapes) _eventShapes->reset(all);
vector<ColourSinglet> singlets;
//singlets.clear();
//singlets = ColourSinglet::getSinglets(all.begin(), all.end());
// for (vector<ColourSinglet>::iterator sItr= singlets.begin(); sItr!=singlets.end();++sItr){
// cout << sItr->momentum().m2() / (GeV*GeV) << " " << sItr->momentum().perp() / GeV << endl;
// }
double havg = 0;
singlets.clear();
if ( theCollapser ) singlets = theCollapser->collapse(all, newStep());
else singlets = ColourSinglet::getSinglets(all.begin(), all.end());
// /* How to use the dipole-based overlap calculator */
- // Ropewalk ropewalk(singlets, 0.5*femtometer, 1.0*GeV);
+ // Ropewalk ropewalk(singlets, 0.5*femtometer, 1.0*GeV, true);
// vector< pair<ColourSinglet,double> > strings = ropewalk.getSingets();
// singlets.clear();
// singlets.reserve(strings.size());
// for ( int i = 0, N = strings.size(); i < N; ++i )
// singlets.push_back(strings[i].first);
enhance->clear();
vector<StringPipe> pipes;
vector<double> maxpT;
vector<size_t> position;
size_t ind = 0;
for(vector<ColourSinglet>::iterator sItr = singlets.begin(); sItr!=singlets.end(); ++sItr){
pipes.push_back(StringPipe(&(*sItr),stringR0,_eventShapes));
// maxpT.push_back(pipes.back().MaxpT()/GeV * pipes.back().MaxpT()/GeV + pipes.back().MaxpTRapidity() * pipes.back().MaxpTRapidity()); // Using pT*pT + y*y to sort
maxpT.push_back(pipes.back().MaxpT()/GeV); // using pT to sort
position.push_back(ind++);
}
// std::sort(position.begin(),position.end(),CompareDoubles(maxpT));
// Loop over singlets, find their h-factor using the handler, send them to the correct pythia
enhance->SetEvent(pipes);
//int cut = 0;
//int total = 0;
Pythia8Interface * pytone = opHandler->GetPythiaPtr(1.0,nev%10000==0);
forcerun = false;
// cout << "---" << endl;
int alls = 0;
int keep = 0;
for(vector<size_t>::iterator it = position.begin(); it!=position.end(); ++it){
vector<ColourSinglet> toHadronization;
//double h = 1.0;
//if(pipes[*it].MaxpT() > stringpTCut && abs(pipes[*it].MaxpTRapidity()) < stringyCut){
// cout << "pT = " << pipes[*it].MaxpT()/GeV << " y = " << pipes[*it].MaxpTRapidity() << endl;
//cout << generator()->currentEvent()->weight() << " " << pipes[*it]._internalOverlap.first << " " << pipes[*it]._internalOverlap.second << endl;
// cut++;
//}
//total++;
alls++;
double h = enhance->KappaEnhancement(pipes[*it]);
if(h > -998.0){
keep++;
Pythia8Interface * pytp = opHandler->GetPythiaPtr(h,nev%10000==0);
toHadronization.push_back(*(pipes[*it].GetSingletPtr()));
if(!hadronizeSystems(*pytp,toHadronization,all)) hadronizeSystems(*pytone,toHadronization,all);
}
/*{
int npartons = toHadronization.back().partons().size();
if(toHadronization.back().nPieces() > 1){
// Dont do any tricks if we have a junction
//cout << "Junction!" << endl;
hadronizeSystems(*pytp,toHadronization,all);
}
else if(npartons>10){
hadronizeSystems(*pytp,toHadronization,all);
}
else{
hadronizeSystems(*pytone,toHadronization,all);
}
}
*/
toHadronization.clear();
// Do short analysis here:
if(doAnalysis){
if(singlets.size() > 10){
ofstream out((analysisPath + "matlabstrings.txt").c_str(),ios::app);
out << "<event>" << endl;
out << PrintStringsToMatlab(singlets) << endl;
out << "</event>" << endl;
out.close();
}
}
}
ofstream out3("ratio.txt",ios::app);
// out3 << generator()->currentEvent()->weight() << " " << alls << " " << keep << endl;
if(doAnalysis){
ofstream out2((analysisPath + "/stringplot.txt").c_str(),ios::app);
out2 << generator()->maximumCMEnergy()/GeV << " " << generator()->currentEvent()->weight() << " " << singlets.size() << " " << havg/singlets.size() << endl;
}
}
void StringFragmentation::dofinish(){
if(doAnalysis!=0){
YODA::mkWriter("yoda");
YODA::WriterYODA::write(analysisPath + generator()->runName() + "1D.yoda",_histograms.begin(),_histograms.end());
YODA::WriterYODA::write(analysisPath + generator()->runName() + "2D.yoda",_histograms2D.begin(),_histograms2D.end());
}
}
bool StringFragmentation::
hadronizeSystems(Pythia8Interface & pyt, const vector<ColourSinglet> & singlets, const tPVector & all) {
TheP8EventShapes * _es = NULL;
Pythia8::Event & event = pyt.event();
pyt.clearEvent();
for ( int i = 0, N = singlets.size(); i < N; ++i ) {
if ( singlets[i].nPieces() == 3 ) {
// Simple junction.
// Save place where we will store dummy particle.
int nsave = event.size();
event.append(22, -21, 0, 0, 0, 0, 0, 0, 0.0, 0.0, 0.0, 0.0);
int npart = 0;
for ( int ip = 1; ip <= 3; ++ip )
for ( int j = 0, M = singlets[i].piece(ip).size(); j < M; ++j ) {
pyt.addParticle(singlets[i].piece(ip)[j], 23, nsave, 0);
++npart;
}
event[nsave].daughters(nsave + 1, nsave + npart);
}
else if ( singlets[i].nPieces() == 5 ) {
// Double, connected junction
// Save place where we will store dummy beam particles.
int nb1 = event.size();
int nb2 = nb1 + 1;
event.append(2212, -21, 0, 0, nb1 + 2, nb1 + 4, 0, 0,
0.0, 0.0, 0.0, 0.0);
event.append(-2212, -21, 0, 0, nb2 + 4, nb2 + 6, 0, 0,
0.0, 0.0, 0.0, 0.0);
// Find the string piece connecting the junctions, and the
// other loose string pieces.
int connector = 0;
int q1 = 0;
int q2 = 0;
int aq1 = 0;
int aq2 = 0;
for ( int ip = 1; ip <= 5; ++ip ) {
if ( singlets[i].sink(ip).first && singlets[i].source(ip).first )
connector = ip;
else if ( singlets[i].source(ip).first ) {
if ( q1 ) q2 = ip;
else q1 = ip;
}
else if ( singlets[i].sink(ip).first ) {
if ( aq1 ) aq2 = ip;
else aq1 = ip;
}
}
if ( !connector || !q1 || !q2 || !aq1 || ! aq2 )
Throw<StringFragError>()
<< name() << " found complicated junction string. Although Pythia8 can "
<< "hadronize junction strings, this one was too complicated."
<< Exception::runerror;
// Insert the partons of the loose triplet ends.
int start = event.size();
for ( int j = 0, M = singlets[i].piece(q1).size(); j < M; ++j )
pyt.addParticle(singlets[i].piece(q1)[j], 23, nb1, 0);
for ( int j = 0, M = singlets[i].piece(q2).size(); j < M; ++j )
pyt.addParticle(singlets[i].piece(q2)[j], 23, nb1, 0);
// Insert dummy triplet incoming parton with correct colour code.
int col1 = singlets[i].piece(connector).empty()? event.nextColTag():
pyt.addColourLine(singlets[i].piece(connector).front()->colourLine());
int dum1 = event.size();
event.append(2, -21, nb1, 0, 0, 0, col1, 0, 0.0, 0.0, 0.0, 0.0, 0.0);
event[nb1].daughters(start, start + singlets[i].piece(q1).size() +
singlets[i].piece(q2).size());
// Insert the partons of the loose anti-triplet ends.
start = event.size();
for ( int j = 0, M = singlets[i].piece(aq1).size(); j < M; ++j )
pyt.addParticle(singlets[i].piece(aq1)[j], 23, nb2, 0);
for ( int j = 0, M = singlets[i].piece(aq2).size(); j < M; ++j )
pyt.addParticle(singlets[i].piece(aq2)[j], 23, nb2, 0);
// Insert dummy anti-triplet incoming parton with correct colour code.
int col2 = singlets[i].piece(connector).empty()? col1:
pyt.addColourLine(singlets[i].piece(connector).back()->antiColourLine());
int dum2 = event.size();
event.append(-2, -21, nb2, 0, 0, 0, 0, col2, 0.0, 0.0, 0.0, 0.0, 0.0);
event[nb2].daughters(start, start + singlets[i].piece(aq1).size() +
singlets[i].piece(aq2).size());
// Add the partons from the connecting string piece.
for ( int j = 0, M = singlets[i].piece(connector).size(); j < M; ++j )
pyt.addParticle(singlets[i].piece(connector)[j], 23, dum1, dum2);
}
else if ( singlets[i].nPieces() > 1 ) {
// We don't know how to handle other junctions yet.
Throw<StringFragError>()
<< name() << " found complicated junction string. Although Pythia8 can "
<< "hadronize junction strings, that interface is not ready yet."
<< Exception::runerror;
} else {
// Normal string
for ( int j = 0, M = singlets[i].partons().size(); j < M; ++j )
pyt.addParticle(singlets[i].partons()[j], 23, 0, 0);
}
}
for ( int i = 0, N = all.size(); i < N; ++i )
if ( !all[i]->coloured() )
pyt.addParticle(all[i], 23, 0, 0);
int oldsize = event.size();
Pythia8::Event saveEvent = event;
int ntry = maxTries;
CurrentGenerator::Redirect stdout(cout, false);
while ( !pyt.go() && --ntry ) event = saveEvent;
if ( !ntry ) {
static DebugItem printfailed("TheP8I::PrintFailed", 10);
if ( printfailed ) {
double ymax = -1000.0;
double ymin = 1000.0;
double sumdy = 0.0;
for ( int i = 0, N = singlets.size(); i < N; ++i )
for ( int j = 0, M = singlets[i].partons().size(); j < M; ++j ) {
const Particle & p = *singlets[i].partons()[j];
ymax = max(ymax, (_es ? _es->yT(p.momentum()) : p.momentum().rapidity()));
//ymax = max(ymax, p.momentum().rapidity());
ymin = min(ymin,(_es ? _es->yT(p.momentum()) : p.momentum().rapidity()));
//ymin = min(ymin, p.momentum().rapidity());
cerr << setw(5) << j << setw(14) << p.momentum().rapidity()
<< setw(14) << p.momentum().perp()/GeV
<< setw(14) << p.momentum().m()/GeV;
if ( j == 0 && p.data().iColour() == PDT::Colour8 ) {
cerr << setw(14) << (p.momentum() + singlets[i].partons().back()->momentum()).m()/GeV;
sumdy += abs((_es ? _es->yT(p.momentum()) : p.momentum().rapidity()) ) -(_es ? _es->yT(singlets[i].partons().back()->momentum())
: singlets[i].partons().back()->momentum().rapidity());
//sumdy += abs(p.momentum().rapidity() - singlets[i].partons().back()->momentum().rapidity());
}
else if ( j > 0 ) {
cerr << setw(14) << (p.momentum() + singlets[i].partons()[j-1]->momentum()).m()/GeV;
sumdy += abs((_es ? _es->yT(p.momentum()) : p.momentum().rapidity()) ) -(_es ? _es->yT(singlets[i].partons()[j-i]->momentum())
: singlets[i].partons()[j-i]->momentum().rapidity());
//sumdy += abs(p.momentum().rapidity() - singlets[i].partons()[j-1]->momentum().rapidity());
if ( j + 1 < M )
cerr << setw(14) << (p.momentum() + singlets[i].partons()[j-1]->momentum()).m()*
(p.momentum() + singlets[i].partons()[j+1]->momentum()).m()/
(p.momentum() + singlets[i].partons()[j+1]->momentum() +
singlets[i].partons()[j-1]->momentum()).m()/GeV;
}
cerr << endl;
}
cerr << setw(14) << ymax - ymin << setw(14) << sumdy/(ymax - ymin) << endl;
}
CurrentGenerator::Redirect stdout2(cout, true);
//event.list();
pyt.errorlist();
cout << "ThePEG event listing:\n" << *(generator()->currentEvent());
Throw<StringFragError>()
<< "Pythia8 failed to hadronize partonic state:\n"
<< stdout2.str() << "This event will be discarded!\n"
<< Exception::warning;
throw Veto();
}
if(window){
// cout << stringyCut << " " << stringpTCut/GeV << endl;
if(forcerun) forcerun = false;
else{
for ( int i = 1, N = event.size(); i < N; ++i ) {
tPPtr p = pyt.getParticle(i);
// cout << p->momentum().perp()/GeV << " " << p->rapidity() << " " << p->id() << endl;
if (p->momentum().perp() > stringpTCut && abs(p->rapidity()) < stringyCut ){
forcerun = true;
return false;
}
}
}
}
// event.list(cerr);
map<tPPtr, set<tPPtr> > children;
map<tPPtr, set<tPPtr> > parents;
for ( int i = 1, N = event.size(); i < N; ++i ) {
tPPtr p = pyt.getParticle(i);
int d1 = event[i].daughter1();
if ( d1 <= 0 ) continue;
children[p].insert(pyt.getParticle(d1));
parents[pyt.getParticle(d1)].insert(p);
int d2 = event[i].daughter2();
if ( d2 > 0 ) {
children[p].insert(pyt.getParticle(d2));
parents[pyt.getParticle(d2)].insert(p);
}
for ( int di = d1 + 1; di < d2; ++di ) {
children[p].insert(pyt.getParticle(di));
parents[pyt.getParticle(di)].insert(p);
}
}
for ( int i = oldsize, N = event.size(); i < N; ++i ) {
PPtr p = pyt.getParticle(i);
set<tPPtr> & pars = parents[p];
if ( !p ) {
Throw<StringFragError>()
<< "Failed to reconstruct hadronized state from Pythia8:\n"
<< stdout.str() << "This event will be discarded!\n" << Exception::warning;
throw Veto();
}
if ( isnan(p->momentum().perp()/GeV) ){
Throw<StringFragError>()
<< "Failed to reconstruct hadronized state from Pythia8:\n"
<< stdout.str() << "This event will be discarded!\n" << Exception::warning;
throw Veto();
}
if ( pars.empty() ) {
Pythia8::Particle & pyp = event[i];
if ( pyp.mother1() > 0 ) pars.insert(pyt.getParticle(pyp.mother1()));
if ( pyp.mother2() > 0 ) pars.insert(pyt.getParticle(pyp.mother2()));
for ( int im = pyp.mother1() + 1; im < pyp.mother2(); ++im )
pars.insert(pyt.getParticle(im));
if ( pars.empty() ) {
Throw<StringFragError>()
<< "Failed to reconstruct hadronized state from Pythia8:\n"
<< stdout.str() << "This event will be discarded!\n" << Exception::warning;
throw Veto();
}
}
if ( pars.size() == 1 ) {
tPPtr par = *pars.begin();
if ( children[par].size() == 1 &&
*children[par].begin() == p && par->id() == p->id() )
newStep()->setCopy(par, p);
else
newStep()->addDecayProduct(par, p);
} else {
newStep()->addDecayProduct(pars.begin(), pars.end(), p);
}
}
return true;
}
string StringFragmentation::PrintStringsToMatlab(vector<ColourSinglet>& singlets) {
stringstream ss;
vector<vector<double> > drawit;
vector<char> names;
for(int i=0;i<26;++i) names.push_back(char(i+97));
vector<char>::iterator nItr = names.begin();
for(vector<ColourSinglet>::iterator sItr = singlets.begin(); sItr!=singlets.end(); ++sItr){
drawit.clear();
for (tcPVector::iterator pItr = sItr->partons().begin(); pItr!=sItr->partons().end();++pItr) {
vector<double> tmp;
tmp.clear();
tmp.push_back((*pItr)->rapidity());
tmp.push_back((*pItr)->vertex().x()/femtometer);
tmp.push_back((*pItr)->vertex().y()/femtometer);
drawit.push_back(tmp);
}
ss << *nItr << " = [";
for(int i=0, N=drawit.size();i<N;++i){
ss << drawit[i].at(0) << ", " << drawit[i].at(1) << ", " << drawit[i].at(2) << ";" << endl;
}
ss << "]';\n\n" << endl;
++nItr;
}
return ss.str();
}
IBPtr StringFragmentation::clone() const {
return new_ptr(*this);
}
IBPtr StringFragmentation::fullclone() const {
return new_ptr(*this);
}
void StringFragmentation::doinitrun() {
HadronizationHandler::doinitrun();
theAdditionalP8Settings.push_back("ProcessLevel:all = off");
theAdditionalP8Settings.push_back("HadronLevel:Decay = off");
theAdditionalP8Settings.push_back("Check:event = off");
theAdditionalP8Settings.push_back("Next:numberCount = 0");
theAdditionalP8Settings.push_back("Next:numberShowLHA = 0");
theAdditionalP8Settings.push_back("Next:numberShowInfo = 0");
theAdditionalP8Settings.push_back("Next:numberShowProcess = 0");
theAdditionalP8Settings.push_back("Next:numberShowEvent = 0");
theAdditionalP8Settings.push_back("Init:showChangedSettings = 0");
theAdditionalP8Settings.push_back("Init:showAllSettings = 0");
theAdditionalP8Settings.push_back("Init:showChangedParticleData = 0");
theAdditionalP8Settings.push_back("Init:showChangedResonanceData = 0");
theAdditionalP8Settings.push_back("Init:showAllParticleData = 0");
theAdditionalP8Settings.push_back("Init:showOneParticleData = 0");
theAdditionalP8Settings.push_back("Init:showProcesses = 0");
vector<string> moresettings = theAdditionalP8Settings;
moresettings.push_back("StringPT:sigma = " + convert(theStringPT_sigma));
moresettings.push_back("StringZ:aLund = " + convert(theStringZ_aLund));
moresettings.push_back("StringZ:bLund = " + convert(theStringZ_bLund));
moresettings.push_back("StringFlav:probStoUD = " +
convert(theStringFlav_probStoUD));
moresettings.push_back("StringFlav:probSQtoQQ = " +
convert(theStringFlav_probSQtoQQ));
moresettings.push_back("StringFlav:probQQ1toQQ0 = " +
convert(theStringFlav_probQQ1toQQ0));
moresettings.push_back("StringFlav:probQQtoQ = " +
convert(theStringFlav_probQQtoQ));
moresettings.push_back("OverlapStrings:fragMass = " +
convert(fragmentationMass));
moresettings.push_back("OverlapStrings:baryonSuppression = " + convert(baryonSuppression));
nev = 0;
// Initialize the OverlapPythia Handler
opHandler = new OverlapPythiaHandler(this,moresettings);
// Should we do event shapes?
_eventShapes = NULL;
if(useThrustAxis==1){
_eventShapes = new TheP8EventShapes();
}
// Select hadronization scheme
// test only to remain backwards compatible.
// TODO: change automated run procedure and adapt
if(fragmentationScheme.find("rope")!=string::npos || fragmentationScheme.find("test")!=string::npos){
enhance = new RandomAverageHandler(stringR0,_eventShapes);
}
else if(fragmentationScheme.find("random")!=string::npos){
enhance = new RandomHandler(stringR0,_eventShapes);
}
else{
enhance = new VanillaPythiaHandler();
}
if( doAnalysis ) {
/*_histograms.push_back(new YODA::Histo1D(20,1,10, "/enhancement" ,"enhancement"));
_histograms.push_back(new YODA::Histo1D(100,0,1., "/parameterA", "parameterA"));
_histograms.push_back(new YODA::Histo1D(100,0,1., "/parameterB", "parameterB"));
_histograms.push_back(new YODA::Histo1D(100,0,1., "/parameterRHO", "parameterRHO"));
_histograms.push_back(new YODA::Histo1D(100,0,1., "/parameterXI", "parameterXI"));
_histograms.push_back(new YODA::Histo1D(100,0,1., "/parameterX", "parameterX"));
_histograms.push_back(new YODA::Histo1D(100,0,1., "/parameterY", "parameterY"));
*/
_histograms.push_back(new YODA::Histo1D(100,1,15,"/h_lowpT","h_lowpT"));
_histograms.push_back(new YODA::Histo1D(100,1,15,"/h_midpT","h_midpT"));
_histograms.push_back(new YODA::Histo1D(100,1,15,"/h_highpT","h_highpT"));
_histograms2D.push_back(new YODA::Histo2D(10,0.,15,10,0,15,"/m_vs_pT","m_vs_pT"));
}
}
void StringFragmentation::persistentOutput(PersistentOStream & os) const {
os
#include "StringFragmentation-output.h"
<< overlapN << testPow << ounit(stringR0,femtometer) << stringyCut << ounit(stringpTCut,GeV) << fragmentationMass << baryonSuppression << ounit(window,true) << useThrustAxis << fragmentationScheme << maxTries << theCollapser << doAnalysis << analysisPath;
}
void StringFragmentation::persistentInput(PersistentIStream & is, int) {
is
#include "StringFragmentation-input.h"
>> overlapN >> testPow >> iunit(stringR0,femtometer) >> stringyCut >> iunit(stringpTCut,GeV) >> fragmentationMass >> baryonSuppression >> iunit(window,true) >> useThrustAxis >> fragmentationScheme >> maxTries >> theCollapser >> doAnalysis >> analysisPath;
}
ClassDescription<StringFragmentation> StringFragmentation::initStringFragmentation;
// Definition of the static class description member.
void StringFragmentation::Init() {
#include "StringFragmentation-interfaces.h"
static Reference<StringFragmentation,ClusterCollapser> interfaceCollapser
("Collapser",
"A ThePEG::ClusterCollapser object used to collapse colour singlet "
"clusters which are too small to fragment. If no object is given the "
"MinistringFragmentetion of Pythia8 is used instead.",
&StringFragmentation::theCollapser, true, false, true, true, false);
static Parameter<StringFragmentation, string> interfaceFragmentationScheme
("FragmentationScheme",
" Choose between hadronization schemes. Choices are: "
" One way strings overlapping in central rapidity bin [OneY], "
" One way strings with cylindrical containers in (b,Y)-space [OnebY], "
" Two way strings with cylindrical containers in (b,Y)-space [TwobY], "
" Ordinary Pythia 8 (no overlap) string hadronization [default]." ,
&StringFragmentation::fragmentationScheme, "");
static Parameter<StringFragmentation,int> interfaceOverlapN
("OverlapN",
"If positive, use the string overlap model with different "
"parameters for up to overlapN strings.",
&StringFragmentation::overlapN, 0, 0, 0,
true, false, Interface::lowerlim);
static Parameter<StringFragmentation,int> interfaceUseThrustAxis
("UseThrustAxis",
"Whether or not to use rapidity wrt. Thrust Axis when counting strings "
"(put 1 or 0).",
&StringFragmentation::useThrustAxis, 0, 0, 0,
true, false, Interface::lowerlim);
static Parameter<StringFragmentation,int> interfaceDoAnalysis
("DoAnalysis",
"Can do a short analysis where histograms of the effective parameters are "
"plotted, and the strings in (bx,by,y)-space are printed for a couple of events. "
"(put 1 or 0).",
&StringFragmentation::doAnalysis, 0, 0, 0,
true, false, Interface::lowerlim);
static Parameter<StringFragmentation,string> interfaceAnalysisPath
("AnalysisPath",
"Set the system path where the (optional) analysis output should appear "
" (directory must exist!)" ,
&StringFragmentation::analysisPath, "");
static Switch<StringFragmentation,bool> interfaceWindow
("StringWindow",
"Enable the 'window'-cut procedure, off by default."
"Parameters will only affect if this is switched on. " ,
&StringFragmentation::window, false, true, false);
static SwitchOption interfaceWindowTrue
(interfaceWindow,"True","enabled.",true);
static SwitchOption interfaceWindowFalse
(interfaceWindow,"False","disabled.",false);
static Parameter<StringFragmentation,Energy> interfaceStringpTCut
("StringpTCut",
"No enhancement of strings with a constituent pT higher than StringpTCut (GeV), and within "
" StringyCut (in GeV).",
&StringFragmentation::stringpTCut, GeV, 6.0*GeV,
0.0*GeV, 0*GeV,
true, false, Interface::lowerlim);
static Parameter<StringFragmentation,double> interfaceStringyCut
("StringyCut",
"No enhancement of strings with a constituent pT higher than StringpTCut (GeV), and within "
" StringyCut.",
&StringFragmentation::stringyCut, 0, 2.0,
0.0, 0,
true, false, Interface::lowerlim);
static Parameter<StringFragmentation,Length> interfaceStringR0
("StringR0",
"In the string overlap model the string R0 dictates the minimum radius "
" of a string (in fm).",
&StringFragmentation::stringR0, femtometer, 0.5*femtometer,
0.0*femtometer, 0*femtometer,
true, false, Interface::lowerlim);
static Parameter<StringFragmentation,double> interfaceFragmentationMass
("FragmentationMass",
"Set the mass used for the f(z) integral in the overlap string model "
" default is pion mass (units GeV).",
&StringFragmentation::fragmentationMass, 0, 0.135,
0., 1.,
true, false, Interface::lowerlim);
static Parameter<StringFragmentation,double> interfaceBaryonSuppression
("BaryonSuppression",
"Set the fudge factor used for baryon suppression in the overlap string model. "
" One day this should be calculated properly. This day is not today. Probably around 2...",
&StringFragmentation::baryonSuppression, 0, 2.0,
1., 10.,
true, false, Interface::lowerlim);
static Parameter<StringFragmentation,double> interfaceOverlapTestPow
("OverlapTestPow",
"Initial value for the exponent used for the string overlap model.",
&StringFragmentation::testPow, 0.5, 0.0, 0,
true, false, Interface::lowerlim);
static Parameter<StringFragmentation,int> interfaceMaxTries
("MaxTries",
"Sometimes Pythia gives up on an event too easily. We therefore allow "
"it to re-try a couple of times.",
&StringFragmentation::maxTries, 2, 1, 0,
true, false, Interface::lowerlim);
}
string StringFragmentation::convert(double d) {
ostringstream os;
os << d;
return os.str();
}

File Metadata

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

Event Timeline