Page MenuHomeHEPForge
Files F10881915

No OneTemporary
Actions

View Options

diff --git a/Hadronization/OverlapEnhancementHandler.h b/Hadronization/OverlapEnhancementHandler.h
--- a/Hadronization/OverlapEnhancementHandler.h
+++ b/Hadronization/OverlapEnhancementHandler.h
@@ -1,64 +1,62 @@
// -*- C++ -*-
#ifndef THEP8I_OverlapEnhancementHandler_H
#define THEP8I_OverlapEnhancementHandler_H
// This is the abstract base class for String Overlap Handlers.
// Implementations of this class are implementations of different ways of calculating
// effective local string constants before the event is send further down stream to Pythia.
#include "TheP8I/Config/Pythia8Interface.h"
#include "StringPipe.h"
namespace TheP8I {
using namespace ThePEG;
class OverlapEnhancementHandler {
/** @name Standard constructors and destructors. */
//@{
- /**
- * The default constructoris protected.
- */
+
public:
OverlapEnhancementHandler(){
}
/**
* The destructor.
*/
virtual ~OverlapEnhancementHandler(){
}
//@}
/* Calculate the effective string tension enhancement factor (h)
for a given string, in the surroundings set by SetEvent()
*/
virtual double KappaEnhancement(StringPipe& pipe) = 0;
/* Set the surroundings of the string, ie. all strings of event
*/
virtual void SetEvent(vector<StringPipe> event) = 0;
/* Add a new pipe to the event
*/
virtual void AddPipe(StringPipe& pipe) = 0;
/* Remove a pipe from the event
*/
virtual bool RemovePipe(StringPipe& pipe) = 0;
virtual void clear() = 0;
protected:
private:
};
}
#endif
diff --git a/Hadronization/OverlapPythiaHandler.cc b/Hadronization/OverlapPythiaHandler.cc
--- a/Hadronization/OverlapPythiaHandler.cc
+++ b/Hadronization/OverlapPythiaHandler.cc
@@ -1,216 +1,216 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the OverlapPythiaHandler class.
//
#include "OverlapPythiaHandler.h"
using namespace TheP8I;
OverlapPythiaHandler::OverlapPythiaHandler(){
}
OverlapPythiaHandler::~OverlapPythiaHandler(){
for ( vector<PythiaPtr>::iterator itr = _overlapPythias.begin();
itr != _overlapPythias.end(); ++itr )
itr->nullify();
_overlapPythias.clear();
}
OverlapPythiaHandler::OverlapPythiaHandler(HadronizationHandler * sf, vector<string>arguments) : _sf(sf){
// Set the arguments
for(vector<string>::iterator itr = arguments.begin(); itr!=arguments.end();++itr){
if(itr->find("StringPT:sigma")!=string::npos){
sigma = PythiaParameter(*itr);
}
else if(itr->find("StringZ:aLund")!=string::npos){
a = PythiaParameter(*itr);
}
else if(itr->find("StringZ:bLund")!=string::npos){
b = PythiaParameter(*itr);
}
else if(itr->find("StringFlav:probStoUD")!=string::npos){
rho = PythiaParameter(*itr);
}
else if(itr->find("StringFlav:probSQtoQQ")!=string::npos){
x = PythiaParameter(*itr);
}
else if(itr->find("StringFlav:probQQ1toQQ0")!=string::npos){
y = PythiaParameter(*itr);
}
else if(itr->find("StringFlav:probQQtoQ")!=string::npos){
xi = PythiaParameter(*itr);
}
else if(itr->find("OverlapStrings:fragMass")!=string::npos){
double m = PythiaParameter(*itr);
m2 = m*m;
}
else if(itr->find("OverlapStrings:baryonSuppression")!=string::npos){
_bsparameter = PythiaParameter(*itr);
}
else{
_pythiaSettings.push_back(*itr);
}
}
// Sanity check of parameters
if (rho < 0 || rho > 1 || x < 0 || x > 1 || y < 0 || y > 1 || xi < 0 ||
xi > 1 || sigma < 0 || sigma > 1 || a < 0.0 || a > 2.0 || b < 0.2 || b > 2.0){
cout << "Did you set up sensible initial Pythia 8 values? Remember:" << endl;
cout << "0 < a < 2; 0.2 < b < 2; 0 < sigma < 1; 0 < xi < 1; 0 < y < 1; 0 < x < 1; 0 < rho < 1" << endl;
}
}
void OverlapPythiaHandler::CreateSinglePythia(double h){
if(!CalculateEffectiveParameters(h)) cout << "Unexpected error setting up parameters for overlap string mode!" << endl;
- /*cout << "Settings and parameters for overlap string hadronization" << endl;
+ /* cout << "Settings and parameters for overlap string hadronization" << endl;
cout << "bs = " << _bsparameter << " m_cut = " << sqrt(m2) << endl;
cout << " h = " << h << " rho_eff = " << rho_eff << " xi_eff = " << xi_eff << " a_eff = " << a_eff << " b_eff = " << b_eff << " x_eff = " << x_eff << " y_eff = " << y_eff << endl;
*/
vector<string> newPythiaSettings = _pythiaSettings;
newPythiaSettings.push_back("StringPT:sigma = " + convert(sigma_eff));
newPythiaSettings.push_back("StringZ:aLund = " + convert(a_eff));
newPythiaSettings.push_back("StringZ:bLund = " + convert(b_eff));
newPythiaSettings.push_back("StringFlav:probStoUD = " + convert(rho_eff));
newPythiaSettings.push_back("StringFlav:probSQtoQQ = " + convert(x_eff));
newPythiaSettings.push_back("StringFlav:probQQ1toQQ0 = " + convert(y_eff));
newPythiaSettings.push_back("StringFlav:probQQtoQ = " + convert(xi_eff));
PythiaPtr tmp;
_overlapPythias.push_back(tmp);
_overlapPythias.back().pPtr = new Pythia8Interface();
_overlapPythias.back().pPtr->init(*_sf,newPythiaSettings);
_overlapPythias.back().h = h;
ComparePythias comparePythias;
sort(_overlapPythias.begin(),_overlapPythias.end(),comparePythias);
}
Pythia8Interface * OverlapPythiaHandler::GetPythiaPtr(double h, bool recycle){
// Broken calculations could give h = nan. This recursive function would thus
// loop infinitely, as nothing == nan. In that case, throw a Veto, and discard the event
// after making sufficient noise.
if(isnan(h)){
cout << "OverlapPythiaHandler::GetPythiaPtr: h is nan. Fix preceeding calculation before running again." << endl;
throw Veto();
}
// Option to delete all unused Pythia ptrs every ef. 10000 events.
// Usefull for Heavy Ion, which contains rare circumstances
if(recycle){
vector<PythiaPtr> newVec;
for(vector<PythiaPtr>::iterator itr = _overlapPythias.begin(); itr!=_overlapPythias.end();++itr){
if ( itr->h !=h ){
itr->nullify( );
}
else{
newVec.push_back(*itr);
}
}
_overlapPythias = newVec;
}
// Check if we already have
for(vector<PythiaPtr>::iterator itr = _overlapPythias.begin(); itr!=_overlapPythias.end();++itr){
if(itr->h==h){
itr->used();
return itr->getpPtr();
}
}
//cout << "Creating new Pythia with h = " << h << endl;
CreateSinglePythia(h);
return GetPythiaPtr(h);
}
vector<double> OverlapPythiaHandler::GetPythiaParameters(double h){
vector<double> ret;
ret.clear();
if(!CalculateEffectiveParameters(h)) cout << "Something went wrong calculating effective Pythia parameters!" << endl;
ret.push_back(h);
ret.push_back(a_eff);
ret.push_back(b_eff);
ret.push_back(rho_eff);
ret.push_back(xi_eff);
ret.push_back(x_eff);
ret.push_back(y_eff);
return ret;
}
bool OverlapPythiaHandler::CalculateEffectiveParameters(double h){
if (h <= 0) return false;
double hinv = 1.0 / h;
rho_eff = pow(rho, hinv);
x_eff = pow(x, hinv);
y_eff = pow(3.0 * y, hinv) / 3.0;
sigma_eff = sigma * sqrt(h);
double C1 = _bsparameter * xi * (2 + rho) / (3 + 2*x*rho + 9*y + 6*x*rho*y + x*x*rho*rho +
3*y*x*x*rho*rho);
xi_eff = (pow(C1, hinv)/_bsparameter) * (3 + 2*x_eff*rho_eff + 9*y_eff +
6*x_eff*rho_eff*y_eff + x_eff*x_eff*rho_eff*rho_eff +
3*y_eff*x_eff*x_eff*rho_eff*rho_eff) / (2 + rho_eff);
if (xi_eff > 1.0) xi_eff = 1.0;
C1 = b / (2 + rho);
b_eff = C1 * (2 + rho_eff);
if (b_eff < 0.2) b_eff = 0.2;
if (b_eff > 2.0) b_eff = 2.0;
double N = IFragmentationF(a, b);
double N_eff = IFragmentationF(a, b_eff);
int s = sign(N - N_eff);
double da = 0.1;
a_eff = a - s * da;
do {
N_eff = IFragmentationF(a_eff, b_eff);
if (sign(N - N_eff) != s) {
s = sign(N - N_eff);
da /= 10.0;
}
a_eff -= s * da;
if (a_eff < 0.0) {a_eff = 0.0; break;}
if (a_eff > 2.0) {a_eff = 2.0; break;}
} while (da > 0.00001);
return true;
}
double OverlapPythiaHandler::IFragmentationF(double a, double b){
const int N = 100000;
const double dz = 1.0 / N; //increment of z variable
double Integral = 0.0;
for (double z = dz; z < 1; z += dz) {
//numerical integral calculation
Integral += pow(1 - z, a) * exp(- b * m2 / z) / z;
}
return Integral * dz;
}
double OverlapPythiaHandler::PythiaParameter(string par){
size_t found = par.find("=");
string number = par.substr(found+2,par.size());
string::iterator end_pos = remove(par.begin(), par.end(), ' ');
return atof(number.c_str());
}
int OverlapPythiaHandler::sign(double num){
return (num < 0) ? -1 : 1;
}
string OverlapPythiaHandler::convert(double d) {
ostringstream os;
os << d;
return os.str();
}
diff --git a/Hadronization/RandomAverageHandler.cc b/Hadronization/RandomAverageHandler.cc
--- a/Hadronization/RandomAverageHandler.cc
+++ b/Hadronization/RandomAverageHandler.cc
@@ -1,96 +1,100 @@
#include "RandomAverageHandler.h"
#include "ThePEG/Repository/UseRandom.h"
using namespace TheP8I;
RandomAverageHandler::RandomAverageHandler(){
}
RandomAverageHandler::~RandomAverageHandler(){
}
RandomAverageHandler::RandomAverageHandler(Length stringR0, TheP8EventShapes * es){
}
double RandomAverageHandler::KappaEnhancement(StringPipe& pipe){
for(vector<StringPipe>::iterator sItr = _pipes.begin(); sItr!=_pipes.end(); ++sItr){
if(pipe._ymin == sItr->_ymin && pipe._ymax == sItr->_ymax && pipe._radius2 == sItr->_radius2){
double m = floor(sItr->_internalOverlap.first + sItr->_externalOverlap.first + 0.5);
double n = floor(sItr->_internalOverlap.second + sItr->_externalOverlap.second + 0.5);
// Find out if the +1 should be left or right.
// if(sItr->_internalOverlap.first > sItr->_internalOverlap.second) m+=1;
// else n+=1;
// This will become nan if one of the containers have a volume of zero.
// In that case no overlap should be counted.
if(isnan(m+n)) return 1;
double p = 0;
double q = 0;
// These two^M^M^M three cases are handled separately
if((m==1 && n==0) || (m==0 && n==1) || (m==0 && n==0)){
p = 1;
q = 0;
}
else{
p = pow(m+n,17./24.)/2;
q = p;
}
+
+
+ // Should we throw the string?
+ if(UseRandom::rnd() > (p + q)/(m + n)) return -999.0;
double N = p + q;
return (0.25*(N + 3 - p*q/N));
}
}
cout << "Could not find pipe..." << endl;
AddPipe(pipe);
return KappaEnhancement(pipe);
}
void RandomAverageHandler::SetEvent(vector<StringPipe> event){
// When we set a new event, we will calculate all overlaps once and for all,
// This will be stored as a std::vector of StringContainers _containers, which is a private
// data member of the present class
_pipes = event;
RecalculateOverlaps();
}
void RandomAverageHandler::AddPipe(StringPipe& pipe){
_pipes.push_back(pipe);
RecalculateOverlaps();
}
void RandomAverageHandler::RecalculateOverlaps(){
// Calculate all external overlaps to get the full number of overlapping strings per container
for(vector<StringPipe>::iterator sItr = _pipes.begin(); sItr!=_pipes.end(); ++sItr){
sItr->_externalOverlap.first = 0;
sItr->_externalOverlap.second = 0;
for(vector<StringPipe>::iterator sItr2 = _pipes.begin(); sItr2!=_pipes.end(); ++sItr2){
sItr->_externalOverlap.first += (sItr->ExternalOverlap(*sItr2)).first;
sItr->_externalOverlap.second += (sItr->ExternalOverlap(*sItr2)).second;
}
}
}
bool RandomAverageHandler::RemovePipe(StringPipe& pipe){
for(vector<StringPipe>::iterator sItr = _pipes.begin(); sItr!=_pipes.end(); ++sItr){
if(pipe._ymin == sItr->_ymin && pipe._ymax == sItr->_ymax && pipe._radius2 == sItr->_radius2){
_pipes.erase(sItr);
RecalculateOverlaps();
return true;
}
}
return false;
}
void RandomAverageHandler::clear(){
_pipes.clear();
}
\ No newline at end of file
diff --git a/Hadronization/Ropewalk.cc b/Hadronization/Ropewalk.cc
--- a/Hadronization/Ropewalk.cc
+++ b/Hadronization/Ropewalk.cc
@@ -1,311 +1,311 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the Ropewalk class.
//
#include "Ropewalk.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/EventRecord/Step.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/CurrentGenerator.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "ThePEG/Utilities/Selector.h"
#include "ThePEG/Utilities/SimplePhaseSpace.h"
#include "ThePEG/Utilities/UtilityBase.h"
#include "ThePEG/Handlers/StepHandler.h"
using namespace TheP8I;
Ropewalk::Ropewalk(const vector<ColourSinglet> & singlets, Length width,
Energy scale, bool deb)
: R0(width), m0(scale), debug(deb) {
for ( int is = 0, Ns = singlets.size(); is < Ns; ++is )
stringdipoles[cloneToFinal(singlets[is])];
getDipoles();
if ( debug ) CurrentGenerator::log() << singlets.size() << "s "
<< dipoles.size() << "d ";
calculateOverlaps();
doBreakups();
}
Ropewalk::~Ropewalk() {
for ( DipoleMap::iterator it = stringdipoles.begin();
it != stringdipoles.end(); ++it ) delete it->first;
}
ColourSinglet * Ropewalk::cloneToFinal(const ColourSinglet & cs) {
tcParticleSet final;
for ( int i = 0, N = cs.partons().size(); i < N; ++i )
final.insert(cs.partons()[i]->final());
tcPPtr p = *final.begin();
if ( p->colourLine() ) return new ColourSinglet(p->colourLine(), final);
if ( p->antiColourLine() ) return new ColourSinglet(p->antiColourLine(), final);
cerr << "Oooops! cloning ColourSinglets failed." << endl;
return 0;
}
LorentzPoint Ropewalk::
propagate(const LorentzPoint & b, const LorentzMomentum & p) const {
return b + p/(p.e()*m0/hbarc);
}
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];
if ( d1.s() <= sqr(m0) ) continue;
// Get the boost to dipoles rest frame and calculate rapidities.
LorentzMomentum pc1 = d1.pc->momentum();
LorentzMomentum pa1 = d1.pa->momentum();
LorentzRotation R = Utilities::boostToCM(make_pair(&pc1, &pa1));
LorentzPoint bc1 = R*propagate(d1.pc->vertex(), d1.pc->momentum());
LorentzPoint ba1 = R*propagate(d1.pa->vertex(), d1.pa->momentum());
double yc1 = limitrap(pc1);
double ya1 = limitrap(pa1);
double n = 0.0;
double m = 0.0;
if ( yc1 <= ya1 ) continue; // don't care about too small strings.
for ( int i2 = 0; i2 < Nd; ++i2 ) {
if ( i1 == i2 ) continue;
// Boost to the rest frame of dipole 1.
Dipole & d2 = dipoles[i2];
if ( d2.s() <= sqr(m0) ) continue;
// Get the boost to dipoles rest frame and calculate rapidities.
LorentzMomentum pc2 = R*d2.pc->momentum();
LorentzMomentum pa2 = R*d2.pa->momentum();
LorentzPoint bc2 = R*propagate(d2.pc->vertex(), d2.pc->momentum());
LorentzPoint ba2 = R*propagate(d2.pa->vertex(), d2.pa->momentum());
double yc2 = limitrap(pc2);
double ya2 = limitrap(pa2);
// Ignore if not overlapping in rapidity.
if ( min(yc2, ya2) > yc1 || max(yc2, ya2) < ya1 || yc2 == ya2 ) continue;
// Calculate rapidity overlap.
double yomax = min(max(yc2, ya2), yc1);
double yomin = max(min(yc2, ya2), ya1);
// Sample a random point in the rapidity overlap region and get
// position in space.
double y = UseRandom::rnd(yomin, yomax);
LorentzPoint b1 = ba1 + (bc1 - ba1)*(y - ya1)/(yc1 - ya1);
LorentzPoint b2 = ba2 + (bc2 - ba2)*(y - ya2)/(yc2 - ya2);
// Skip if there is no overlap in transverse position.
if ( (b1 - b2).perp() > R0 ) continue;
// Register the overlap.
double yfrac = (yomax - yomin)/(yc1 - ya1);
if ( yc2 > ya2 ) {
d1.overlaps.push_back(make_pair(&d2, yfrac));
n += yfrac;
} else {
d1.overlaps.push_back(make_pair(&d2, -yfrac));
m += yfrac;
}
}
d1.n = int(n + UseRandom::rnd());
d1.m = int(m + UseRandom::rnd());
d1.initMultiplet();
}
}
void Ropewalk::Dipole::initMultiplet() {
typedef pair<int,int> Plet;
int ns = 0;
int ms = 0;
while ( ns < n || ms < m ) {
Plet diff;
double octetveto = double(ns + ms + 1 - p - q)/double(ns + ms + 1);
if ( double(n - ns) > UseRandom::rnd()*double(n + m - ns - ms) ) {
Selector<Plet> sel;
sel.insert(multiplicity(p + 1, q), Plet(1, 0));
sel.insert(multiplicity(p, q - 1)*octetveto, Plet(0, -1));
sel.insert(multiplicity(p - 1, q + 1), Plet(-1, 1));
diff = sel[UseRandom::rnd()];
++ns;
} else {
Selector<Plet> sel;
sel.insert(multiplicity(p, q + 1), Plet(0, 1));
sel.insert(multiplicity(p - 1, q)*octetveto, Plet(-1, 0));
sel.insert(multiplicity(p + 1, q - 1), Plet(1, -1));
diff = sel[UseRandom::rnd()];
++ms;
}
p += diff.first;
q += diff.second;
}
/* *** ATTENTION *** maybe better if Christian implements this */
}
double Ropewalk::Dipole::breakupProbability() const {
if ( !breakable() || n + m <= 0.0 || n + m + 1 == p + q ) return 0.0;
double sum = 0.0;
for (int j = 0, N = overlaps.size(); j < N; ++j )
if ( overlaps[j].first->breakable() )
sum += abs(overlaps[j].second);
if ( sum <= 0.0 ) return 1.0;
return double(n + m + 1 - p - q)/(sum + 1.0);
}
Step & Ropewalk::step() const {
return *(CurrentGenerator()->currentStepHandler()->newStep());
}
void Ropewalk::doBreakups() {
// We need to be able to select diquarks.
Selector<int> disel;
disel.insert(1.0, ParticleID::ud_0);
disel.insert(3.0, ParticleID::dd_1);
disel.insert(3.0, ParticleID::ud_1);
disel.insert(3.0, ParticleID::uu_1);
typedef multimap<double, const Dipole *> BMap;
BMap breakups;
// First order alldipoles in increasing probability of breaking
for ( int i = 0, N = dipoles.size(); i < N; ++i )
breakups.insert(make_pair(dipoles[i].breakupProbability(), &dipoles[i]));
// Then start breaking in reverse order
for ( BMap::reverse_iterator it = breakups.rbegin();
it != breakups.rend(); ++it ) {
const Dipole & d = *it->second;
if ( d.breakupProbability() < UseRandom::rnd() ) continue;
// Select di-quarks
int flav = disel[UseRandom::rnd()];
PPtr dic = CurrentGenerator()->getParticle(flav);
PPtr dia = CurrentGenerator()->getParticle(-flav);
// Get boost to dipole rest frame and place the di-quarks there.
LorentzRotation R = Utilities::getBoostFromCM(make_pair(d.pc, d.pa));
SimplePhaseSpace::CMS(dic, dia, d.s(), 1.0, 0.0);
dia->transform(R);
dic->transform(R);
// Add the di-quarks as children to the decayed gluons and split
// the resulting string.
if ( !( step().addDecayProduct(d.pc, dic, true) &&
step().addDecayProduct(d.pa, dia, true) &&
step().addDecayProduct(d.pc, dia, false) &&
step().addDecayProduct(d.pa, dic, false) ) )
cerr << "Oooops! adding decay products failed" << endl;
tcParticleSet pset(d.string->partons().begin(), d.string->partons().end());
pset.erase(d.pc);
pset.erase(d.pa);
pset.insert(dic);
pset.insert(dia);
// remove the old stringand insert the new ones fixing the
// pointers from the dipolmes.
ColourSinglet * olds = d.string;
const vector<Dipole *> & oldd = stringdipoles[olds];
vector<ColourSinglet> newstrings =
ColourSinglet::getSinglets(pset.begin(), pset.end());
for ( int is = 0, Ns = newstrings.size(); is < Ns; ++is ) {
ColourSinglet * news = new ColourSinglet(newstrings[is]);
tcParticleSet pset(news->partons().begin(), news->partons().end());
vector<Dipole *> & newd = stringdipoles[news];
for ( int id = 0, Nd = oldd.size(); id < Nd; ++id ) {
if ( pset.find(oldd[id]->pc) != pset.end() ) {
newd.push_back(oldd[id]);
oldd[id]->string = news;
}
}
}
stringdipoles.erase(olds);
delete olds;
}
}
-vector< pair<ColourSinglet,double> > Ropewalk::getSingets() const {
+vector< pair<ColourSinglet,double> > Ropewalk::getSinglets() 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;
}
toty += sumy;
totyh += sumyh;
if ( sumy > 0.0 ) {
double h = 0.1*int(10.0*sumyh/sumy + UseRandom::rnd());
ret.push_back(make_pair(*it->first, h));
toth += sumyh/sumy;
toth2 += sqr(sumyh/sumy);
minh = min(minh, sumyh/sumy);
maxh = max(maxh, sumyh/sumy);
} else {
ret.push_back(make_pair(*it->first, 1.0));
toth += 1.0;
toth2 += 1.0;
minh = min(minh, 1.0);
maxh = max(maxh, 1.0);
}
}
if ( debug )
CurrentGenerator::log() << ret.size() << "ns "
<< totyh/max(toty, 0.00001) << "hwa "
<< minh << "<h "
<< toth/ret.size() << "ha "
<< maxh << ">h "
<< sqrt(toth2/ret.size() - sqr(toth/ret.size())) << "hsd " << endl;
return ret;
}
diff --git a/Hadronization/Ropewalk.h b/Hadronization/Ropewalk.h
--- a/Hadronization/Ropewalk.h
+++ b/Hadronization/Ropewalk.h
@@ -1,226 +1,226 @@
// -*- C++ -*-
#ifndef TheP8I_Ropewalk_H
#define TheP8I_Ropewalk_H
//
// This is the declaration of the Ropewalk class.
//
#include "TheP8I/Config/Pythia8Interface.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/EventRecord/ColourSinglet.h"
namespace TheP8I {
using namespace ThePEG;
/**
* Here is the documentation of the Ropewalk class.
*/
class Ropewalk {
public:
/**
* Container class containing information about individual dipole overlaps.
*/
struct Dipole {
/**
* Constructor.
*/
Dipole(tcPPtr p1, tcPPtr p2, ColourSinglet & str)
: pc(p1), pa(p2), string(&str), n(0), m(0), p(1), q(0) {}
/**
* Return true if this dipole is breakable. Only dipoles between
* gluons (not belonging to other dipoles which have broken) and
* which has a mass above 2 GeV can break.
*/
bool breakable() const {
return pc->id() == ParticleID::g && pa->id() == ParticleID::g &&
!pc->decayed() && !pa->decayed() && s() > 4.0*GeV2;
}
/**
* Calculate the probability that this dipole should break. Taking
* into account the number of negative steps in the colour space
* and the fraction of overlapping dipoles which are able to
* break.
*/
double breakupProbability() const;
/**
* Set and return the effective multiplet.
*/
void initMultiplet();
/**
* Calculate the multiplicity given pp and qq;
*/
static double multiplicity(int pp, int qq) {
return ( pp < 0 || qq < 0 || pp + qq == 0 )? 0.0:
0.5*(pp + 1)*(qq + 1)*(pp + qq + 2);
}
/**
* Calculate the kappa-enhancement.
*/
double kappaEnhancement() const {
return 0.25*(p + q + 3 - p*q/double(p + q));
}
/**
* Return the squared invariant mass of this dipole.
*/
Energy2 s() const {
return (pc->momentum() + pa->momentum()).m2();
}
/**
* Return the effective rapidityspan of this dipole.
*/
double yspan(Energy m0) const {
return s() > sqr(m0)? 0.5*log(s()/sqr(m0)): 0.0;
}
/**
* The coloured particle.
*/
tcPPtr pc;
/**
* The anti-coloured particle.
*/
tcPPtr pa;
/**
* The overlapping dipoles with the amount of overlap (negative
* means anti overlap).
*/
vector< pair<const Dipole *, double> > overlaps;
/**
* The string to which this Dipole belongs.
*/
ColourSinglet * string;
/**
* The summed parallel (n) and anti-parallel (m) overlaps from
* other dipoles.
*/
int n, m;
/**
* The multiplet.
*/
int p, q;
};
/**
* Convenient typedef
*/
typedef map<ColourSinglet *, vector<Dipole *> > DipoleMap;
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
Ropewalk(const vector<ColourSinglet> & singlets, Length width,
Energy scale, bool deb = false);
/**
* The destructor.
*/
virtual ~Ropewalk();
//@}
/**
* Return all ColourSinglet objects with their kappa enhancement.
*/
- vector< pair<ColourSinglet,double> > getSingets() const;
+ vector< pair<ColourSinglet,double> > getSinglets() const;
protected:
/**
* Extract all dipoles.
*/
void getDipoles();
/**
* Calculate all overlaps and initialize multiplets in all dipoles.
*/
void calculateOverlaps();
/**
* Break dipoles (and strings)
*/
void doBreakups();
/**
* Propagate a parton a finite time inthe lab-system.
*/
LorentzPoint propagate(const LorentzPoint & b,
const LorentzMomentum & p) const;
/**
* Calculate the rapidity of a Lorentz vector but limit the
* transverse mass to be above m0.
*/
double limitrap(const LorentzMomentum & p) const;
/**
* Return the current step.
*/
Step & step() const;
/**
* Make a copy of a ColourSinglet making sure all partons are final.
*/
static ColourSinglet * cloneToFinal(const ColourSinglet & cs);
private:
/**
* The assumed radius of a string.
*/
Length R0;
/**
* The assumed mass for calculating rapidities
*/
Energy m0;
/**
* The vector of all Dipoles
*/
vector<Dipole> dipoles;
/**
* All Dipoles in all string
*/
DipoleMap stringdipoles;
/**
* Debug flag.
*/
bool debug;
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
Ropewalk & operator=(const Ropewalk &);
};
}
#endif /* TheP8I_Ropewalk_H */

File Metadata

Mime Type
text/x-diff
Expires
Sat, May 3, 6:56 AM (6 h, 51 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4980847
Default Alt Text
(28 KB)

Event Timeline