Page MenuHomeHEPForge
Files F8309438

No OneTemporary
Actions

View Options

diff --git a/Hadronization/StringFragmentation.cc b/Hadronization/StringFragmentation.cc
--- a/Hadronization/StringFragmentation.cc
+++ b/Hadronization/StringFragmentation.cc
@@ -1,1089 +1,1100 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the StringFragmentation class.
//
#include "StringFragmentation.h"
#include "Ropewalk.h"
#include "RandomAverageHandler.h"
#include "RandomHandler.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/Utilities/EnumIO.h"
#include "ThePEG/Utilities/MaxCmp.h"
#include "ThePEG/Utilities/HoldFlag.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include <algorithm>
using namespace TheP8I;
StringFragmentation::StringFragmentation()
: pythia(), fScheme(0), stringR0(0.5*femtometer), stringm0(1.0*GeV), junctionDiquark(1.0), alpha(1.0), average(true),
stringyCut(2.0), stringpTCut(6*GeV), fragmentationMass(0.134),
- baryonSuppression(0.5), window(false), minStringy(8.0), longSoft(false), throwaway(false), shoving(false), deltay(5.0), tshove(1.0*femtometer),
+ baryonSuppression(0.5), window(false), minStringy(8.0), longSoft(false), longStringsOnly(false), throwaway(false), shoving(false), deltay(5.0), tshove(1.0*femtometer),
deltat(0.1*femtometer), rCutOff(4.0), gAmplitude(1.0*GeV/femtometer), gExponent(1.0), yCutOff(4.0), lorentz(true), cylinder(false), _eventShapes(0),
useThrustAxis(0), opHandler(0), maxTries(2), doAnalysis(0), analysisPath(""),
#include "StringFragmentation-init.h"
{
}
StringFragmentation::~StringFragmentation() {
if ( opHandler ) delete opHandler;
}
void StringFragmentation::handle(EventHandler & eh, const tPVector & tagged,
const Hint & h) {
++nev;
static unsigned long lastEventId = 0;
bool secondary = false;
// Get the event
tPVector all = RemnantDecayer::decayRemnants(tagged, *newStep());
HoldFlag<int> oldFS(fScheme, fScheme);
// Prevent funny behavior if hadronizing decays of eg. Upsilon.
if ( newStep()->collision()->uniqueId == lastEventId ){
secondary = true;
if(fScheme == 4 || fScheme == 5 || fScheme == 1) fScheme = 99;
else fScheme = 0;
}
else
lastEventId = newStep()->collision()->uniqueId;
if(_eventShapes) _eventShapes->reset(all);
vector<ColourSinglet> singlets;
singlets.clear();
if ( theCollapser && !secondary )
singlets = theCollapser->collapse(all, newStep());
else
singlets = ColourSinglet::getSinglets(all.begin(), all.end());
// Goto correct hadronization scheme
switch(fScheme){
case 0: // Pythia
{
hadronizeSystems(*opHandler->GetPythiaPtr(1.0,nev%10000==0), singlets, all);
break;
}
case 1: // For playing around and printing stuff
{
// Cartoon of string shoving
//
if(singlets.size() > 10){
Ropewalk ropewalkInit(singlets, stringR0, stringm0,
baryonSuppression, throwaway, false);
vector< pair<ColourSinglet,double> > allStrings =
ropewalkInit.getSinglets(stringyCut);
tPVector allParticles;
for ( vector< pair<ColourSinglet,double> >::iterator sItr = allStrings.begin();
sItr != allStrings.end(); ++sItr )
for( tcPVector::iterator pItr = sItr->first.partons().begin();
pItr != sItr->first.partons().end(); ++pItr ){
allParticles.push_back(const_ptr_cast<tPPtr>(*pItr));
}
singlets = theCollapser->collapse(allParticles, newStep());
for(double time = 0.0; time < 1.0; time += 0.1){
Ropewalk ropewalkOne(singlets, stringR0, stringm0,
baryonSuppression, false, shoving, rCutOff, gAmplitude, gExponent,
yCutOff, deltay, time*femtometer, 0.1*femtometer, lorentz, cylinder, stringpTCut, minStringy, longSoft, NULL, false);
auto css = ropewalkOne.getSinglets(1.0);
vector<ColourSinglet> cs;
for(auto itr = css.begin(); itr != css.end(); ++itr)
cs.push_back(itr->first);
ofstream out((analysisPath + "../lhcevent.m").c_str(),ios::app);
out << PrintStringsToMatlab(cs) << endl;
out << "figure;" << endl;
out << "tubeplot(a,radius);" << endl;
out << "hold on; tubeplot(b,radius);" << endl;
out << "hold on; tubeplot(c,radius);" << endl;;
out << "hold on; tubeplot(d,radius);" << endl;;
out << "hold on; tubeplot(e,radius);" << endl;;
out << "hold on; tubeplot(f,radius);" << endl;;
out << "hold on; tubeplot(g,radius);" << endl;;
out << "hold on; tubeplot(h,radius);" << endl;;
out << "hold on; tubeplot(i,radius);" << endl;;
out << "hold on; tubeplot(j,radius);" << endl;;
out << "hold on; tubeplot(k,radius);" << endl;;
out << "hold on; tubeplot(l,radius);" << endl;;
out << "hold on; tubeplot(m,radius);" << endl;;
}
exit(1);
}
/*if(throwaway){
Ropewalk ropewalk_init(singlets, stringR0, stringm0, baryonSuppression, throwaway, false);
vector< pair<ColourSinglet,double> > allStrings = ropewalk_init.getSinglets(stringyCut);
tPVector allParticles;
for(vector< pair<ColourSinglet,double> >::iterator sItr = allStrings.begin(); sItr != allStrings.end(); ++sItr)
for(tcPVector::iterator pItr = sItr->first.partons().begin(); pItr != sItr->first.partons().end(); ++pItr)
allParticles.push_back(const_ptr_cast<tPPtr>(*pItr));
singlets = theCollapser->collapse(allParticles, newStep());
}
Ropewalk ropewalk(singlets, stringR0, stringm0, baryonSuppression, false, false);
vector< pair<ColourSinglet,double> > strings = ropewalk.getSinglets(stringyCut);
static ofstream out(( generator()->runName() + ".txt").c_str());
if(!std::isnan(ropewalk.lambdaSum()+ropewalk.getkW()+ropewalk.getNb()))
out << generator()->currentEvent()->weight() << " " << ropewalk.lambdaSum() << " " << ropewalk.getkW() << " " << ropewalk.getNb() << endl;
//vector<double> mspec = ropewalk.mspec;
//for(int i = 0, N = mspec.size(); i < N; ++i) out << generator()->currentEvent()->weight() << " " << mspec[i] << endl;
*/
break;
}
case 2: // Do shoving and construct parton level histograms.
{
double weight = generator()->currentEvent()->weight();
if(doAnalysis)
hk.setWeight(weight);
if (throwaway) {
Ropewalk ropewalkInit(singlets, stringR0, stringm0, baryonSuppression, throwaway, false);
if(doAnalysis){
hk.hist("stringLengthBefore")->fill(ropewalkInit.lambdaSum(),weight);
hk.hist("stringLengthNoCutBefore")->fill(ropewalkInit.lambdaSum(0.0),weight);
}
vector< pair<ColourSinglet,double> > allStrings = ropewalkInit.getSinglets(stringyCut);
tPVector allParticles;
for ( vector< pair<ColourSinglet,double> >::iterator sItr = allStrings.begin(); sItr != allStrings.end(); ++sItr )
for( tcPVector::iterator pItr = sItr->first.partons().begin(); pItr != sItr->first.partons().end(); ++pItr ){
allParticles.push_back(const_ptr_cast<tPPtr>(*pItr));
}
singlets = theCollapser->collapse(allParticles, newStep());
}
if(doAnalysis){
// string length before shoving
Ropewalk rr(singlets,stringR0,stringm0,baryonSuppression,false,false);
hk.hist("stringLengthCollapse")->fill(rr.lambdaSum(),weight);
hk.hist("stringLengthNoCutCollapse")->fill(rr.lambdaSum(0.0),weight);
}
typedef multimap<Energy,Ropewalk::DipoleMap::const_iterator> OrderedMap;
OrderedMap ordered;
Ropewalk ropewalk(singlets, stringR0, stringm0,
baryonSuppression, false, shoving, rCutOff, gAmplitude, gExponent,
yCutOff, deltay, tshove, deltat, lorentz, cylinder, stringpTCut, minStringy, longSoft, (doAnalysis ? &hk : NULL), false);
if(doAnalysis){
hk.hist("stringLengthShove")->fill(ropewalk.lambdaSum(),weight);
hk.hist("stringLengthNoCutShove")->fill(ropewalk.lambdaSum(0.0),weight);
}
for(Ropewalk::DipoleMap::iterator itr = ropewalk.begin(); itr != ropewalk.end(); ++itr)
ordered.insert(make_pair(maxPT(*itr->first, stringyCut), itr));
for(OrderedMap::iterator it = ordered.begin(); it != ordered.end(); ++it ) {
if(doAnalysis){
hk.hist("averageKappa")->fill(ropewalk.averageKappaEnhancement(it->second, stringyCut));
}
if(!pythia.getRopeUserHooksPtr()->setDipoles(&(*it->second), stringm0, stringR0, !average, alpha)) {
pythia.getRopeUserHooksPtr()->setEnhancement(ropewalk.averageKappaEnhancement(it->second, stringyCut));
for (int i = 0, N = it->second->second.size(); i < N; ++i)
it->second->second[i]->hadr = true;
}
vector<ColourSinglet> toHadronization(1, *(it->second->first));
forcerun = false;
if(!hadronizeSystems(pythia, toHadronization, all))
hadronizeSystems(pythia, toHadronization, all);
}
break;
}
case 3: // Pipes with average
{
RandomAverageHandler avg_enhancer(throwaway);
// TODO: Implement throwaway functionality in this
avg_enhancer.clear();
vector<StringPipe> pipes;
// Make all the pipes
for(vector<ColourSinglet>::iterator sItr = singlets.begin(); sItr!=singlets.end(); ++sItr)
pipes.push_back(StringPipe(&(*sItr),stringR0,_eventShapes));
avg_enhancer.SetEvent(pipes);
for(vector<StringPipe>::iterator it = pipes.begin(); it!=pipes.end(); ++it){
vector<ColourSinglet> toHadronization;
double h = avg_enhancer.KappaEnhancement(*it);
if(h > 0){
toHadronization.push_back(*(*it).GetSingletPtr());
forcerun = false;
if(!hadronizeSystems(*(opHandler->GetPythiaPtr(h,nev%10000==0)),toHadronization,all))
hadronizeSystems(*(opHandler->GetPythiaPtr(1.0,false)),toHadronization,all);
}
}
break;
}
case 4: // Average kappa over whole strings
{
Ropewalk ropewalk(singlets, stringR0, stringm0, baryonSuppression,throwaway, false);
vector< pair<ColourSinglet,double> > strings = ropewalk.getSinglets(stringyCut);
vector<ColourSinglet> toHadronization;
double avge = 0;
for(int i = 0, N = strings.size(); i < N; ++i ){
avge += strings[i].second;
pythia.getRopeUserHooksPtr()->setEnhancement(strings[i].second);
toHadronization.clear();
toHadronization.push_back(strings[i].first);
hadronizeSystems(pythia,toHadronization,all);
}
// ofstream out(( "/scratch/galette/bierlich/work/dipsy/pprange/" + generator()->runName() + "/stringplot.txt").c_str(),ios::app);
// out << generator()->currentEvent()->weight() << " " << strings.size() << " " << avge << endl;
break;
}
case 5: // Altering kappa per dipole
{
double weight = generator()->currentEvent()->weight();
if(doAnalysis)
hk.setWeight(weight);
if (throwaway) {
Ropewalk ropewalkInit(singlets, stringR0, stringm0, baryonSuppression, throwaway, false);
if(doAnalysis){
hk.hist("stringLengthBefore")->fill(ropewalkInit.lambdaSum(),weight);
hk.hist("stringLengthNoCutBefore")->fill(ropewalkInit.lambdaSum(0.0),weight);
}
vector< pair<ColourSinglet,double> > allStrings = ropewalkInit.getSinglets(stringyCut);
tPVector allParticles;
for ( vector< pair<ColourSinglet,double> >::iterator sItr = allStrings.begin(); sItr != allStrings.end(); ++sItr )
for( tcPVector::iterator pItr = sItr->first.partons().begin(); pItr != sItr->first.partons().end(); ++pItr ){
allParticles.push_back(const_ptr_cast<tPPtr>(*pItr));
}
singlets = theCollapser->collapse(allParticles, newStep());
}
if(doAnalysis){
// string length before shoving
Ropewalk rr(singlets,stringR0,stringm0,baryonSuppression,false,false);
hk.hist("stringLengthCollapse")->fill(rr.lambdaSum(),weight);
hk.hist("stringLengthNoCutCollapse")->fill(rr.lambdaSum(0.0),weight);
}
typedef multimap<Energy,Ropewalk::DipoleMap::const_iterator> OrderedMap;
OrderedMap ordered;
Ropewalk ropewalk(singlets, stringR0, stringm0,
baryonSuppression, false, shoving, rCutOff, gAmplitude, gExponent,
yCutOff, deltay, tshove, deltat, lorentz, cylinder, stringpTCut, minStringy, longSoft, (doAnalysis ? &hk : NULL), false);
if(doAnalysis){
hk.hist("stringLengthShove")->fill(ropewalk.lambdaSum(),weight);
hk.hist("stringLengthNoCutShove")->fill(ropewalk.lambdaSum(0.0),weight);
}
for(Ropewalk::DipoleMap::iterator itr = ropewalk.begin(); itr != ropewalk.end(); ++itr)
ordered.insert(make_pair(maxPT(*itr->first, stringyCut), itr));
for(OrderedMap::iterator it = ordered.begin(); it != ordered.end(); ++it ) {
if(doAnalysis){
hk.hist("averageKappa")->fill(ropewalk.averageKappaEnhancement(it->second, stringyCut));
}
else if(!pythia.getRopeUserHooksPtr()->setDipoles(&(*it->second), stringm0, stringR0, !average, alpha)) {
pythia.getRopeUserHooksPtr()->setEnhancement(ropewalk.averageKappaEnhancement(it->second, stringyCut));
for (int i = 0, N = it->second->second.size(); i < N; ++i)
it->second->second[i]->hadr = true;
}
vector<ColourSinglet> toHadronization(1, *(it->second->first));
forcerun = false;
if(!hadronizeSystems(pythia, toHadronization, all))
hadronizeSystems(pythia, toHadronization, all);
}
- ordered.clear();
- for(Ropewalk::DipoleMap::iterator itr = ropewalk.beginSpectators(); itr != ropewalk.endSpectators(); ++itr)
- ordered.insert(make_pair(maxPT(*itr->first, stringyCut), itr));
- for(OrderedMap::iterator it = ordered.begin(); it != ordered.end(); ++it ) {
- if(doAnalysis){
- hk.hist("averageKappa")->fill(ropewalk.averageKappaEnhancement(it->second, stringyCut));
+ if(!longStringsOnly){
+ ordered.clear();
+ for(Ropewalk::DipoleMap::iterator itr = ropewalk.beginSpectators(); itr != ropewalk.endSpectators(); ++itr)
+ ordered.insert(make_pair(maxPT(*itr->first, stringyCut), itr));
+ for(OrderedMap::iterator it = ordered.begin(); it != ordered.end(); ++it ) {
+ if(doAnalysis){
+ hk.hist("averageKappa")->fill(ropewalk.averageKappaEnhancement(it->second, stringyCut));
+ }
+ pythia.getRopeUserHooksPtr()->setEnhancement(1.0);
+ vector<ColourSinglet> toHadronization(1, *(it->second->first));
+ forcerun = false;
+ if(!hadronizeSystems(pythia, toHadronization, all))
+ hadronizeSystems(pythia, toHadronization, all);
}
- pythia.getRopeUserHooksPtr()->setEnhancement(1.0);
- vector<ColourSinglet> toHadronization(1, *(it->second->first));
- forcerun = false;
- if(!hadronizeSystems(pythia, toHadronization, all))
- hadronizeSystems(pythia, toHadronization, all);
}
-
- break;
+ break;
}
case 99: // New Pythia
{
pythia.getRopeUserHooksPtr()->setEnhancement(1.0);
hadronizeSystems(pythia, singlets, all);
break;
}
default:
cout << "We should really not be here. This is bad..." << endl;
}
// fScheme = oldFS;
// Do short analysis here:
if(doAnalysis){
ofstream out((analysisPath + "analysis.txt").c_str(),ios::app);
// out << "<event>" << endl;
// out << PrintStringsToMatlab(singlets) << endl;
// out << "</event>" << endl;
out.close();
}
}
void StringFragmentation::dofinish(){
if(doAnalysis){
ofstream of((analysisPath + generator()->runName()) + ".histo");
of << hk;
of.close();
}
//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());
//}
}
Energy StringFragmentation::maxPT(const ColourSinglet & cs, double deltaY) {
MaxCmp<Energy2> maxpt2;
for ( int i = 0, N = cs.partons().size(); i < N; ++i ) {
if ( deltaY > 0.0 && abs(cs.partons()[i]->eta()) > deltaY ) continue;
maxpt2(cs.partons()[i]->momentum().perp2());
}
if ( !maxpt2 ) return ZERO;
return sqrt(maxpt2.value());
}
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;
//event.list(cout);
//event.list(cerr);
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 ){
/* if(abs(p->id()) == 310 || abs(p->id()) == 3122){
static ofstream fout(( generator()->runName() + ".txt").c_str(),ios::app);
tcPVector pVector = singlets[0].partons();
fout << p->id() << " " << generator()->currentEvent()->weight() << " ";
for(size_t j=0;j<pVector.size();++j){
//if(pVector[j]->id()<20){
fout << pVector[j]->id() << " " << pVector[j]->momentum().perp()/GeV;
fout.close();
//}
}
*/
forcerun = true;
return false;
}
}
}
}
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 ) {
event.list(cout);
cout << i << endl;
Throw<StringFragError>()
<< "Failed to reconstruct hadronized state from Pythia8:\n"
<< stdout.str() << "This event will be discarded!\n" << Exception::warning;
throw Veto();
}
if ( std::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();
if(fScheme == 2){
// book histograms
hk.bookHisto("impactParameter",10,0,2);
// hk.bookSurfHisto("twopPeriph",100,0,10)
}
if(doAnalysis && fScheme!=2){
// Book some histograms
hk.bookHisto("stringlength",100,0,20);
hk.bookHisto("stringptMax",100,0,15);
hk.bookHisto("stringpTMaxCut",100,0,15);
hk.bookHisto("averageKappa",10,1,2);
hk.bookHisto("averageKappa",10,1,2);
hk.bookHisto("averageKappa",10,1,2);
hk.bookHisto("stringLengthBefore",100,0,1000);
hk.bookHisto("stringLengthCollapse",100,0,1000);
hk.bookHisto("stringLengthShove",100,0,1000);
hk.bookHisto("stringLengthNoCutBefore",100,0,1000);
hk.bookHisto("stringLengthNoCutCollapse",100,0,1000);
hk.bookHisto("stringLengthNoCutShove",100,0,1000);
hk.bookHisto("yDistBefore",100,0,10);
hk.bookHisto("yDistAfter",100,0,10);
hk.bookHisto("excitationPt",100,0,10);
hk.bookHisto("excitationRestFramePt",100,0,10);
hk.bookHisto("normalGluonPt",100,0,10);
hk.bookHisto("normalQuarkPt",100,0,10);
hk.bookHisto("normalDiQuarkPt",100,0,10);
hk.bookHisto("normalGluonRestFramePt",100,0,10);
hk.bookHisto("normalQuarkRestFramePt",100,0,10);
hk.bookHisto("normalDiQuarkRestFramePt",100,0,10);
hk.bookHisto("excitedEndGluonPt",100,0,10);
hk.bookHisto("excitedQuarkPt",100,0,10);
hk.bookHisto("excitedEndGluonRestFramePt",100,0,10);
hk.bookHisto("excitedQuarkRestFramePt",100,0,10);
hk.bookHisto("excitedDiQuarkPt",100,0,10);
hk.bookHisto("excitedDiQuarkRestFramePt",100,0,10);
hk.bookHisto("nExcitations",100,0,200);
hk.bookHisto("nPartons",100,0,200);
hk.bookHisto("shove1",100,0,1);
hk.bookHisto("shove2",100,0,1);
hk.bookHisto("shove3",100,0,1);
hk.bookHisto("shove4",100,0,1);
hk.bookHisto("shove5",100,0,1);
hk.bookHisto("dist1",100,0,1);
hk.bookHisto("dist2",100,0,1);
hk.bookHisto("dist3",100,0,1);
hk.bookHisto("dist4",100,0,1);
hk.bookHisto("dist5",100,0,1);
hk.bookHisto("shoveInRestFrame",1000,0,10);
hk.bookHisto("rejectedShoveInRestFrame",1000,0,10);
}
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");
if(fScheme == 4 || fScheme == 5 || fScheme == 1){ // We don't need multiple Pythia objects anymore!
pythia.enableHooks();
pythia.init(*this,theAdditionalP8Settings);
//if(pythia.version() > 0 && pythia.version() - 1.234 > 1.0){
// cout << "Pythia version: " << pythia.version() << endl;
// cout << "The chosen fragmentation scheeme requires a tweaked "
// << "Pythia v. 1.234 for option. I will default you to "
// "option 'pythia'." << endl;
// fScheme = 0;
// }
// else{
PytPars p;
p.insert(make_pair("StringPT:sigma",theStringPT_sigma));
p.insert(make_pair("StringZ:aLund",theStringZ_aLund));
p.insert(make_pair("StringZ:bLund",theStringZ_bLund));
p.insert(make_pair("StringFlav:probStoUD",theStringFlav_probStoUD));
p.insert(make_pair("StringFlav:probSQtoQQ",theStringFlav_probSQtoQQ));
p.insert(make_pair("StringFlav:probQQ1toQQ0",theStringFlav_probQQ1toQQ0));
p.insert(make_pair("StringFlav:probQQtoQ",theStringFlav_probQQtoQ));
phandler.init(fragmentationMass*fragmentationMass,baryonSuppression,p);
pythia.getRopeUserHooksPtr()->setParameterHandler(&phandler);
//pythia.getRopeUserHooksPtr()->setWindow(window,stringpTCut);
// }
}
if( !( fScheme == 4 || fScheme == 5 || fScheme == 1) ){ // Not 'else' as it can change above...
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));
// Initialize the OverlapPythia Handler
if ( opHandler ) delete opHandler;
opHandler = new OverlapPythiaHandler(this,moresettings);
}
// Should we do event shapes?
if(_eventShapes) delete _eventShapes;
_eventShapes = NULL;
if(useThrustAxis==1){
_eventShapes = new TheP8EventShapes();
}
//if( doAnalysis ) {
// _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"));
//}
nev = 0;
}
void StringFragmentation::persistentOutput(PersistentOStream & os) const {
os
#include "StringFragmentation-output.h"
<< fScheme << ounit(stringR0,femtometer) << ounit(stringm0,GeV) << junctionDiquark << alpha << average << stringyCut
<< ounit(stringpTCut,GeV) << fragmentationMass << baryonSuppression
- << oenum(window) << minStringy << oenum(longSoft) << oenum(throwaway) << oenum(shoving) << deltay << ounit(tshove,femtometer) << ounit(deltat,femtometer)
+ << oenum(window) << minStringy << oenum(longSoft) << oenum(longStringsOnly) << oenum(throwaway) << oenum(shoving) << deltay << ounit(tshove,femtometer) << ounit(deltat,femtometer)
<< rCutOff << ounit(gAmplitude,GeV/femtometer) << gExponent << yCutOff << oenum(lorentz) << oenum(cylinder) << useThrustAxis << maxTries
<< theCollapser << doAnalysis << analysisPath;
}
void StringFragmentation::persistentInput(PersistentIStream & is, int) {
is
#include "StringFragmentation-input.h"
>> fScheme >> iunit(stringR0,femtometer) >> iunit(stringm0,GeV) >> junctionDiquark >> alpha >> average >> stringyCut
>> iunit(stringpTCut,GeV) >> fragmentationMass >> baryonSuppression
- >> ienum(window) >> minStringy >> ienum(longSoft) >> ienum(throwaway) >> ienum(shoving) >> deltay >> iunit(tshove,femtometer) >> iunit(deltat,femtometer)
+ >> ienum(window) >> minStringy >> ienum(longSoft) >> ienum(longStringsOnly) >> ienum(throwaway) >> ienum(shoving) >> deltay >> iunit(tshove,femtometer) >> iunit(deltat,femtometer)
>> rCutOff >> iunit(gAmplitude,GeV/femtometer) >> gExponent >> yCutOff >> ienum(lorentz) >> ienum(cylinder) >> useThrustAxis >> 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 Switch<StringFragmentation,int> interfaceFragmentationScheme
("FragmentationScheme",
"Different options for how to handle overlapping strings.",
&StringFragmentation::fScheme, 0, true, false);
static SwitchOption interfaceFragmentationSchemepythia
(interfaceFragmentationScheme,
"pythia",
"Plain old Pythia fragmentation",
0);
static SwitchOption interfaceFragmentationSchemedep1
(interfaceFragmentationScheme,
"dep1",
"Not sure about this.",
1);
static SwitchOption interfaceFragmentationSchemedep2
(interfaceFragmentationScheme,
"dep2",
"Not sure about this.",
2);
static SwitchOption interfaceFragmentationSchemenone
(interfaceFragmentationScheme,
"none",
"Plain old Pythia fragmentation",
0);
static SwitchOption interfaceFragmentationSchemepipe
(interfaceFragmentationScheme,
"pipe",
"Each string is enclosed by a cylinder. Effective parameters are calculated from the overlap of cylinders.",
3);
static SwitchOption interfaceFragmentationSchemeaverage
(interfaceFragmentationScheme,
"average",
"The overlap is calculated for each dipole in all strings, and for each string the effective parameters are obtained from the average overlap.",
4);
static SwitchOption interfaceFragmentationSchemedipole
(interfaceFragmentationScheme,
"dipole",
"Effective parameters are calculated for each breakup by determining the overlap of the corresponding dipole with other dipoles.",
5);
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> interfaceThrowAway
("ThrowAway",
"Throw away strings with weight:"
"Use 1 - (p + q)/(m + n) for (0,-1) configuration." ,
&StringFragmentation::throwaway, false, true, false);
static SwitchOption interfaceThrowAwayTrue
(interfaceThrowAway,"True","enabled.",true);
static SwitchOption interfaceThrowAwayFalse
(interfaceThrowAway,"False","disabled.",false);
static Switch<StringFragmentation,bool> interfaceShoving
("Shoving",
"Strings will shove each other around.",
&StringFragmentation::shoving, false, true, false);
static SwitchOption interfaceShovingTrue
(interfaceShoving,"True","enabled.",true);
static SwitchOption interfaceShovingFalse
(interfaceShoving,"False","disabled.",false);
static Parameter<StringFragmentation,double> interfaceDeltaY
("DeltaY",
"The size of rapidity slicing intervals",
&StringFragmentation::deltay, 0, 0.1,
0.0, 0,
true, false, Interface::lowerlim);
static Parameter<StringFragmentation,Length> interfaceTShove
("TShove",
"The total shoving -and propagation time",
&StringFragmentation::tshove, femtometer, 1.0*femtometer,
0.0*femtometer, 0*femtometer,
true, false, Interface::lowerlim);
static Parameter<StringFragmentation,Length> interfaceDeltaT
("DeltaT",
"Shoving time slicing -- size of a slice",
&StringFragmentation::deltat, femtometer, 0.1*femtometer,
0.0*femtometer, 0*femtometer,
true, false, Interface::lowerlim);
static Parameter<StringFragmentation,double> interfaceRCutOff
("RCutOff",
"Cutoff radius at which we dont let strings shove each"
"other around anymore -- given as a multiple of StringR0",
&StringFragmentation::rCutOff, 0, 4.0,
0.0, 0,
true, false, Interface::lowerlim);
static Parameter<StringFragmentation,StringTension> interfaceGAmplitude
("GAmplitude",
"Amplitude of shoving energy density function, should be approx."
"the average effective string tension",
&StringFragmentation::gAmplitude, GeV/femtometer, 1.0*GeV/femtometer,
0.0*GeV/femtometer, 0*GeV/femtometer,
true, false, Interface::lowerlim);
static Parameter<StringFragmentation,double> interfaceGExponent
("GExponent",
"Parameter to divide the exponent of energy density function for"
"shoving.",
&StringFragmentation::gExponent, 0, 10.0,
0.0, 0,
true, false, Interface::lowerlim);
static Parameter<StringFragmentation,double> interfaceYCutOff
("YCutOff",
"Cutoff rapidity for excitation clustering. This is the minimal"
"effective rapidity span for a dipole after clustering",
&StringFragmentation::yCutOff, 0, 5.0,
0.0, 0,
true, false, Interface::lowerlim);
static Switch<StringFragmentation,bool> interfaceLorentz
("Lorentz",
"Do Lorentz contraction of string fields moving with a transverse velocity",
&StringFragmentation::lorentz, true, true, false);
static SwitchOption interfaceLorentzTrue
(interfaceLorentz,"True","enabled.",true);
static SwitchOption interfaceLorentzFalse
(interfaceLorentz,"False","disabled.",false);
static Switch<StringFragmentation,bool> interfaceCylinder
("Cylinder",
"Do Cylinder Area Shoving",
&StringFragmentation::cylinder, false, true, false);
static SwitchOption interfaceCylinderTrue
(interfaceCylinder,"True","enabled.",true);
static SwitchOption interfaceCylinderFalse
(interfaceCylinder,"False","disabled.",false);
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",
"Strings with a constituent pT higher than StringpTCut (GeV), will not participate "
" in the ropewalk",
&StringFragmentation::stringpTCut, GeV, 3.0*GeV,
0.0*GeV, 0*GeV,
true, false, Interface::lowerlim);
static Parameter<StringFragmentation,double> interfaceMinStringy
("MinStringy",
"String smaller than MinStringy (in rapidity) will not participate "
" in the ropewalk",
&StringFragmentation::minStringy, 0, 8.0,
0.0, 0,
true, false, Interface::lowerlim);
static Switch<StringFragmentation,bool> interfaceLongSoft
("LongSoft",
"Only do ropewalk with long and soft strings, set cuts accordingly" ,
&StringFragmentation::longSoft, false, true, false);
static SwitchOption interfaceLongSoftTrue
(interfaceLongSoft,"True","enabled.",true);
static SwitchOption interfaceLongSoftFalse
(interfaceLongSoft,"False","disabled.",false);
- static Parameter<StringFragmentation,Energy> interfaceStringm0
+ static Switch<StringFragmentation,bool> interfaceLongStringsOnly
+ ("LongStringsOnly",
+ "Hadronize only long and soft strings, no spectators, set cuts accordingly" ,
+ &StringFragmentation::longStringsOnly, false, true, false);
+ static SwitchOption interfaceLongStringsOnlyTrue
+ (interfaceLongStringsOnly,"True","enabled.",true);
+ static SwitchOption interfaceLongStringsOnlyFalse
+ (interfaceLongStringsOnly,"False","disabled.",false);
+
+
+ static Parameter<StringFragmentation,Energy> interfaceStringm0
("Stringm0",
".",
&StringFragmentation::stringm0, GeV, 1.0*GeV,
0.0*GeV, 0*GeV,
true, false, Interface::lowerlim);
static Parameter<StringFragmentation,double> interfaceJunctionDiquark
("JunctionDiquark",
"Suppress diquark production in breaking of junctions with this amount",
&StringFragmentation::junctionDiquark, 0, 1.0,
0.0, 0,
true, false, Interface::lowerlim);
static Parameter<StringFragmentation,double> interfaceAlpha
("Alpha",
"Boost the enhanced string tension with additional factor",
&StringFragmentation::alpha, 0, 1.0,
1.0, 0,
true, false, Interface::lowerlim);
static Switch<StringFragmentation,bool> interfaceAverage
("Average",
"Disable to try and hadronise strings with individual tensions"
"instead of average value." ,
&StringFragmentation::average, false, true, false);
static SwitchOption interfaceAverageTrue
(interfaceAverage,"True","enabled.",true);
static SwitchOption interfaceAverageFalse
(interfaceAverage,"False","disabled.",false);
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, 0.5,
0.0, 1.0,
true, false, Interface::limited);
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();
}
diff --git a/Hadronization/StringFragmentation.h b/Hadronization/StringFragmentation.h
--- a/Hadronization/StringFragmentation.h
+++ b/Hadronization/StringFragmentation.h
@@ -1,389 +1,390 @@
// -*- C++ -*-
#ifndef THEP8I_StringFragmentation_H
#define THEP8I_StringFragmentation_H
//
// This is the declaration of the StringFragmentation class.
//
#include "ThePEG/Handlers/HadronizationHandler.h"
#include "ThePEG/Handlers/ClusterCollapser.h"
#include "TheP8I/Config/Pythia8Interface.h"
#include "OverlapPythiaHandler.h"
#include "TheP8EventShapes.h"
#include <sstream>
#include "histtools/HistKeeper.h"
//#include "YODA/Histo1D.h"
//#include "YODA/Histo2D.h"
//#include "YODA/Writer.h"
//#include "YODA/WriterYODA.h"
namespace TheP8I {
bool sorting_principle(pair<ThePEG::ColourSinglet*, vector<TheP8I::Ropewalk::Dipole*> > c1, pair<ThePEG::ColourSinglet*, vector<TheP8I::Ropewalk::Dipole*> > c2){
return (c1.first->momentum().perp() < c2.first->momentum().perp());
}
using namespace ThePEG;
/**
* Here is the documentation of the StringFragmentation class.
*
* @see \ref StringFragmentationInterfaces "The interfaces"
* defined for StringFragmentation.
*/
class StringFragmentation: public HadronizationHandler {
typedef map<string, double> PytPars;
friend class Ropewalk::Dipole;
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
StringFragmentation();
/**
* The destructor.
*/
virtual ~StringFragmentation();
//@}
public:
/**
* The main function called by the EventHandler class to
* perform the Hadronization step.
* @param eh the EventHandler in charge of the Event generation.
* @param tagged if not empty these are the only particles which should
* be considered by the StepHandler.
* @param hint a Hint object with possible information from previously
* performed steps.
* @throws Veto if the StepHandler requires the current step to be
* discarded.
* @throws Stop if the generation of the current Event should be stopped
* after this call.
* @throws Exception if something goes wrong.
*/
virtual void handle(EventHandler & eh, const tPVector & tagged,
const Hint & hint);
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Let the given Pythia8Interface hadronize the ColourSinglet
* systems provided.
*/
bool hadronizeSystems(Pythia8Interface & pyt, const vector<ColourSinglet> & singlets,
const tPVector & all);
void hadronizeTweak(Pythia8Interface & pyt, const ColourSinglet & singlet,
const tPVector & all);
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Helper function to get the maximum transverse momenta of any
* parton in a string. Optionally only consider partons within a
* rapidity interval |eta| < deltaY.
*/
static Energy maxPT(const ColourSinglet & cs, double deltaY = 0.0);
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
//@}
private:
HistKeeper hk;
/**
* The interface to the Pythia 8 object,
*/
Pythia8Interface pythia;
/**
* Internal switch for the fragmentation scheme
*/
int fScheme;
/**
* Object that handles calculation of new Pythia parameters
*/
ParameterHandler phandler;
/**
* The intrinsic string radius
*/
Length stringR0;
/**
* The assumed mass for calculating rapidities in the reeperbahn scheme
*/
Energy stringm0;
/**
* Supression of diquarks emerging from breaking junctions
*/
double junctionDiquark;
double alpha;
bool average;
/**
*????????? If **window** is set, this determines the abs(y) window for which no enhancement will be made
*/
double stringyCut;
/**
* Cut on maxpT for a string to go into the ropewalk. REquires **longsoft** to be set
*/
Energy stringpTCut;
/**
* Assumed m0 value in calculation of normalization of f(z) the Lund splitting function.
*/
double fragmentationMass;
/**
* Parameter surpressing baryonic content further in calculation of effective parameters
*/
double baryonSuppression;
// Force hadronize systems to run to the end, regardless of windows
bool forcerun;
/**
* Can provide a window in y and pT where no enhancement is made, in order to not include
* very hard gluons in central rapidity region in a rope
*/
bool window;
/*
* Minimal rapidity span for a string to go into the ropewalk. REquires **longsoft** to be set
*/
double minStringy;
/*
* Only allow long and soft strings to go into the ropewalk
*/
bool longSoft;
+ bool longStringsOnly;
/**
* Throw away some string entirely while making ropes
*/
bool throwaway;
/*
* Let the strings shove each other around
*/
bool shoving;
/*
* The interval with which we slice in rapidity
*/
double deltay;
/*
* Total shoving -and propagation time
*/
Length tshove;
/*
* Time slicing
*/
Length deltat;
/*
* Cutoff radius at which we don't let dipoles affect each other any more
* given as a multiple of R0
*/
double rCutOff;
typedef Qty<-1,1,0> StringTension;
/*
* Amplitude of energy density function for string shoving
*/
StringTension gAmplitude;
/*
* Exponent multiple of energy density function for string shoving
*/
double gExponent;
/*
* Cutoff rapidity for merging of string excitations
*/
double yCutOff;
/*Lorentz contraction of fields
*
*/
bool lorentz;
/*
* Cylinder shoving
*/
bool cylinder;
TheP8EventShapes * _eventShapes;
/**
* Use thrust axis to calculate rapidities; useful for LEP validation
*/
int useThrustAxis;
/**
* Additional interfaces to Pythia8 objects in case the string
* overlap model is to be used.
*/
OverlapPythiaHandler * opHandler;
/**
* Sometimes Pythia gives up on an event too easily. We allow it to
* re-try a couple of times.
*/
int maxTries;
// Keep track of how many events we already hadronized using this object, in order to clean up
int nev;
// Can do a small analysis
int doAnalysis;
//vector<YODA::Histo1D*> _histograms;
//vector<YODA::Histo2D*> _histograms2D;
string analysisPath;
void dofinish();
// Can print the strings of the event in a matlab friendly format
string PrintStringsToMatlab(vector<ColourSinglet>& singlets);
string convert(double d);
#include "StringFragmentation-var.h"
/**
* The object used to avoid too small strings in the hadronization.
*/
ClusterCollapserPtr theCollapser;
public:
/**
* Exception class declaration.
*/
struct StringFragError: public Exception {};
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<StringFragmentation> initStringFragmentation;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
StringFragmentation & operator=(const StringFragmentation &);
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of StringFragmentation. */
template <>
struct BaseClassTrait<TheP8I::StringFragmentation,1> {
/** Typedef of the first base class of StringFragmentation. */
typedef HadronizationHandler NthBase;
};
/** This template specialization informs ThePEG about the name of
* the StringFragmentation class and the shared object where it is defined. */
template <>
struct ClassTraits<TheP8I::StringFragmentation>
: public ClassTraitsBase<TheP8I::StringFragmentation> {
/** Return a platform-independent class name */
static string className() { return "TheP8I::StringFragmentation"; }
/**
* The name of a file containing the dynamic library where the class
* StringFragmentation is implemented. It may also include several, space-separated,
* libraries if the class StringFragmentation depends on other classes (base classes
* excepted). In this case the listed libraries will be dynamically
* linked in the order they are specified.
*/
static string library() { return "libTheP8I.so"; }
};
/** @endcond */
}
#endif /* THEP8I_StringFragmentation_H */

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 3:22 PM (1 d, 50 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023272
Default Alt Text
(57 KB)

Event Timeline