Page MenuHomeHEPForge
Files F7879637

No OneTemporary
Actions

View Options

diff --git a/Hadronization/StringFragmentation.cc b/Hadronization/StringFragmentation.cc
--- a/Hadronization/StringFragmentation.cc
+++ b/Hadronization/StringFragmentation.cc
@@ -1,520 +1,521 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the StringFragmentation class.
//
#include "StringFragmentation.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"
using namespace TheP8I;
StringFragmentation::StringFragmentation()
: pythia(), overlapN(0), testPow(1.0), stringR0(0.5*femtometer), fragmentationMass(0.134), baryonSuppression(2.), useThrustAxis(0), fragmentationScheme("none"), maxTries(2), doAnalysis(0), analysisPath(""),
#include "StringFragmentation-init.h"
{
}
StringFragmentation::~StringFragmentation() {
// delete enhance;
// delete opHandler;
}
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();
if ( theCollapser ) singlets = theCollapser->collapse(all, newStep());
else singlets = ColourSinglet::getSinglets(all.begin(), all.end());
// Set surroundings to this event
//vector<ColourSinglet> * singPtr = singlets;
enhance->SetEvent(&singlets);
// Loop over singlets, find their h-factor using the handler, send them to the correct pythia
for(vector<ColourSinglet>::iterator sItr = singlets.begin(); sItr!=singlets.end(); ++sItr){
double h = enhance->KappaEnhancement(sItr);
Pythia8Interface * pytp = opHandler->GetPythiaPtr(h,(nev%10000==0 ? true : false));
// Is this really correct? Or do we delete the event every time?
vector<ColourSinglet> sig;
// Do short analysis here:
if(doAnalysis!=0){
// Print the first 10 events to matlab-friendly format
if(nev<10){
ofstream out((analysisPath + "matlabstrings.txt").c_str(),ios::app);
out << "<event>" << endl;
out << PrintStringsToMatlab(singlets) << endl;
out << "</event>" << endl;
out.close();
}
double stringWeight = sItr->momentum().m2() / (TeV*TeV);
double weight = generator()->currentEvent()->weight();
weight *= stringWeight;
vector<double> values = opHandler->GetPythiaParameters(h);
for(int i=0, N=_histograms.size();i<N;++i){
_histograms[i]->fill(values[i],weight);
}
}
sig.push_back(*sItr);
hadronizeSystems(*pytp,sig,all);
sig.clear();
}
}
void StringFragmentation::dofinish(){
if(doAnalysis!=0){
YODA::mkWriter("yoda");
YODA::WriterYODA::write(analysisPath + generator()->runName() + ".yoda",_histograms.begin(),_histograms.end());
}
}
void 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();
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();
}
// 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 ( 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);
}
}
}
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
if(fragmentationScheme.find("OnebY")!=string::npos){
enhance = new OnebYHandler(stringR0,_eventShapes);
}
else if(fragmentationScheme.find("OneY")!=string::npos){
enhance = new OneYHandler(-1,1);
}
else if(fragmentationScheme.find("TwobY")!=string::npos){
enhance = new TwobYHandler(stringR0,_eventShapes);
}
else{
enhance = new VanillaPythiaHandler();
}
- _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"));
-
+ 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"));
+ }
}
void StringFragmentation::persistentOutput(PersistentOStream & os) const {
os
#include "StringFragmentation-output.h"
<< overlapN << testPow << ounit(stringR0,femtometer) << fragmentationMass << baryonSuppression << useThrustAxis << fragmentationScheme << maxTries << theCollapser << doAnalysis << analysisPath;
}
void StringFragmentation::persistentInput(PersistentIStream & is, int) {
is
#include "StringFragmentation-input.h"
>> overlapN >> testPow >> iunit(stringR0,femtometer) >> fragmentationMass >>baryonSuppression >> 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 Parameter<StringFragmentation,Length> interfaceStringR0
("StringR0",
"In the string overlap model the string R0 dictates the minimum radius "
" of a string (im 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
Tue, Nov 19, 8:39 PM (1 d, 2 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3806106
Default Alt Text
(22 KB)

Event Timeline