Page MenuHomeHEPForge
Files F10881699

No OneTemporary
Actions

View Options

diff --git a/Hadronization/OverlapPythiaHandler.cc b/Hadronization/OverlapPythiaHandler.cc
--- a/Hadronization/OverlapPythiaHandler.cc
+++ b/Hadronization/OverlapPythiaHandler.cc
@@ -1,203 +1,213 @@
// -*- 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(){
}
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 << "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->reuse==0&&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();
}
\ No newline at end of file
diff --git a/Hadronization/RandomHandler.h b/Hadronization/RandomHandler.h
--- a/Hadronization/RandomHandler.h
+++ b/Hadronization/RandomHandler.h
@@ -1,114 +1,122 @@
// -*- C++ -*-
#include <iostream>
#include <vector>
#include <math.h>
#include "ThePEG/Repository/RandomGenerator.h"
using namespace std;
#ifndef THEP8I_Plet
#define THEP8I_Plet
class Plet{
public:
Plet(int p, int q) : _p(p), _q(q) {
-
+ // Construct SU(3) multiplet given quantum numbers p and q
}
Plet add(int p, int q) {
return(Plet(_p + p,_q + q));
}
Plet add(Plet t) {
return add(t.p(),t.q());
}
double multiplicity(){
+ // The multiplicity of an SU(3) multiplet
if(_p<0 || _q<0) return 0;
return 0.5*(_p + 1)*(_q + 1)*(_p + _q + 2);
}
Plet addRandomPlet(vector<Plet> plets, double ran){
+ // Take a Plet at random from the plets vector, and add it to *this plet
+ // using multiplicities from above to calculate probability
vector<Plet> ret;
for(vector<Plet>::iterator it = plets.begin(); it!=plets.end(); ++it) ret.push_back(add(*it));
-
+
+ // Calculate normalization constant
double norm = 0;
for(vector<Plet>::iterator it = ret.begin();it!=ret.end();++it) norm+=it->multiplicity();
+ // Sanity check
if(ret.size() == 0){
cout << "Could not walk!" << endl;
return *this;
}
+ // Calculate cumulated probability to walk to a given state, and check whether the random number is less than that.
+ // If yes, return said state
double cumprob = 0;
for(vector<Plet>::iterator it = ret.begin();it!=ret.end();++it){
cumprob += it->multiplicity()/norm;
if(ran < cumprob) return *it;
}
+ cout << "We should never reach this point!" << endl;
return ret.back();
}
int p(){
return _p;
}
int q(){
return _q;
}
int sum(){
return _p+_q;
}
private:
int _p,_q;
};
#endif
#ifndef THEP8I_RandomHandler_H
#define THEP8I_RandomHandler_H
#include "OverlapEnhancementHandler.h"
#include "TheP8EventShapes.h"
#include "StringPipe.h"
namespace TheP8I {
using namespace ThePEG;
class RandomHandler : public OverlapEnhancementHandler {
public:
RandomHandler();
~RandomHandler();
RandomHandler(Length stringR0, TheP8EventShapes * es);
double KappaEnhancement(StringPipe& pipe);
void SetEvent(vector<StringPipe> event);
void AddPipe(StringPipe& pipe);
void RecalculateOverlaps();
bool RemovePipe(StringPipe& pipe);
void clear();
private:
vector<StringPipe> _pipes;
Length _stringR0;
TheP8EventShapes * _es;
vector<Plet> addAntiTriplet;
vector<Plet> addTriplet;
};
}
#endif
\ No newline at end of file
diff --git a/Hadronization/StringPipe.cc b/Hadronization/StringPipe.cc
--- a/Hadronization/StringPipe.cc
+++ b/Hadronization/StringPipe.cc
@@ -1,145 +1,154 @@
#include "StringPipe.h"
using namespace TheP8I;
StringPipe::StringPipe() {
}
StringPipe::~StringPipe(){
}
StringPipe::StringPipe(ColourSinglet* singlet, Length r0, TheP8EventShapes * es) : _externalOverlap(make_pair(0,0)), _es(es), theSinglet(singlet) {
tcPPtr pmin = singlet->partons()[0];
tcPPtr pmax = singlet->partons()[0];
for (tcPVector::iterator pItr = singlet->partons().begin(); pItr!=singlet->partons().end();++pItr) {
// Find the lowest and the highest rapidities to span the full string between
double y = _es ? _es->yT((*pItr)->momentum()) : (*pItr)->rapidity();
if(y<(_es ? _es->yT(pmin->momentum()) : pmin->rapidity())) pmin = *pItr;
else if(y>(_es ? _es->yT(pmax->momentum()) : pmax->rapidity())) pmax = *pItr;
}
_originb.first = (pmin->vertex().x() + pmax->vertex().x()) / 2;
_originb.second = (pmin->vertex().y() + pmax->vertex().y()) / 2;
_ymin = _es ? _es->yT(pmin->momentum()) : pmin->rapidity();
_ymax = _es ? _es->yT(pmax->momentum()) : pmax->rapidity();
pair<Area,Area> Vcs = make_pair(0*femtometer*femtometer,0*femtometer*femtometer);
tcPVector::iterator pminus = singlet->partons().begin();
+ // Calculate pipe radius using distance from parton to mean line
+ // between partons with max and min probability
for (tcPVector::iterator pItr = singlet->partons().begin(); pItr!=singlet->partons().end();++pItr) {
- // Grow the cylinder radius
+ // Grow the pipe radius
Length dx = (*pItr)->vertex().x() - _originb.first;
Length dy = (*pItr)->vertex().y() - _originb.second;
_radius2 += dx*dx + dy*dy;
// Grow the volume of the colour singlet
if(pItr!=singlet->partons().begin()){
// Length dbx = (*pItr)->vertex().x() - (*pminus)->vertex().x();
// Length dby = (*pItr)->vertex().y() - (*pminus)->vertex().y();
double dy = _es ? _es->yT((*pItr)->momentum()) - _es->yT((*pminus)->momentum()) :
(*pItr)->rapidity() - (*pminus)->rapidity();
- // ORIGINAL BELOW -- THIS IS CHANGE PER 20.05.2014
- if(dy>0) Vcs.first += r0 * r0 * dy;
- if(dy<0) Vcs.second += r0 * r0 * dy;
+ // ORIGINAL BELOW -- THIS IS CHANGED PER 20.05.2014, pipe ends should now be parallel
+ // to impact parameter plane. Below they are allowed not to be.
+ if(dy>0) Vcs.first += r0 * r0 * abs(dy);
+ if(dy<0) Vcs.second += r0 * r0 * abs(dy);
// Area db2 = dbx * dbx + dby * dby;
// if(dy>0) Vcs.first += (r0 * r0 * dy * dy) / (sqrt(db2/(r0*r0) +dy*dy));
// if(dy<0) Vcs.second += (r0 * r0 * dy * dy) / (sqrt(db2/(r0*r0) +dy*dy));
++pminus;
}
}
+ // Use radius squared instead of volume to save us from a factor of pi.
_radius2 /= (singlet->partons().end() - singlet->partons().begin());
_radius2 += r0*r0;
+ // Overlap is given in both directions
_internalOverlap.first = Vcs.first / (_radius2 * (_ymax - _ymin));
_internalOverlap.second = Vcs.second / (_radius2 * (_ymax - _ymin));
}
pair<double,double> StringPipe::ExternalOverlap(StringPipe& other){
- if(this==&other){
+ if(this==&other || other.GetVolume() == 0*femtometer*femtometer){
return make_pair(0,0);
}
return make_pair(other._internalOverlap.first * OverlapY(other) * OverlapArea(other) / other.GetVolume(),
other._internalOverlap.second * OverlapY(other) * OverlapArea(other) / other.GetVolume());
}
Area StringPipe::GetVolume(){
return _radius2*M_PI*(_ymax - _ymin);
}
Area StringPipe::OverlapArea(StringPipe& other){
+ // This is the overlap between two cylinders. Purely geometrical.
Area d2 = (_originb.first - other._originb.first) * (_originb.first - other._originb.first) +
(_originb.second - other._originb.second) * (_originb.second - other._originb.second);
Length d = sqrt(d2);
Length r = sqrt(_radius2);
Length R = sqrt(other._radius2);
Qty<4,0,0,1,1,1> argu = (-d + r + R) * (d + r - R) * (d - r + R) * (d + r + R);
-
+ // Save us from dividing by zero
if(argu <= 0*femtometer*femtometer*femtometer*femtometer) return 0*femtometer*femtometer;
return _radius2 * acos((d2 + _radius2 - other._radius2)/(2*d*r)) + other._radius2 * acos((d2 + other._radius2 - _radius2)/(2*d*R)) -
0.5 * sqrt(argu);
}
double StringPipe::IntevalOverlap(double min1, double min2, double max1, double max2){
+ // Helper function to calculate overlap between two 1D intervals.
+ // This could very well exist somewhere in ThePEG already.
if(max1<min2||max2<min1||(min1==min2&&max1==max2))
return 0.;
double ret = 0;
if(min1<min2){
if(max1<max2) ret = abs(max1 - min2);
else {
ret = abs(max2 - min2);
}
}
else {
if(max2<max1) ret = abs(max2 - min1);
else {
ret = abs(max1 - min1);
}
}
return ret;
}
double StringPipe::MaxpTRapidity(){
+ // Gives the rapidity of the parton with the maximal pT
Energy pT = 0*GeV;
double y = 0;
for (tcPVector::iterator pItr = theSinglet->partons().begin(); pItr!=theSinglet->partons().end();++pItr) {
if(pT < (*pItr)->momentum().perp()){
pT = (*pItr)->momentum().perp();
y = (*pItr)->rapidity();
}
}
return y;
}
Energy StringPipe::MeanpT(){
Energy pT = 0*GeV;
for (tcPVector::iterator pItr = theSinglet->partons().begin(); pItr!=theSinglet->partons().end();++pItr) {
pT += (*pItr)->momentum().perp();
}
return pT/(distance(theSinglet->partons().begin(),theSinglet->partons().end()));
}
Energy StringPipe::MaxpT(){
Energy pT = 0*GeV;
for (tcPVector::iterator pItr = theSinglet->partons().begin(); pItr!=theSinglet->partons().end();++pItr) {
if(pT < (*pItr)->momentum().perp()) pT = (*pItr)->momentum().perp();
}
return pT;
}
double StringPipe::OverlapY(StringPipe& other){
return IntevalOverlap(_ymin, other._ymin, _ymax, other._ymax);
}
ColourSinglet * StringPipe::GetSingletPtr(){
return theSinglet;
}
\ No newline at end of file

File Metadata

Mime Type
text/x-diff
Expires
Sat, May 3, 6:42 AM (1 d, 7 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4977175
Default Alt Text
(16 KB)

Event Timeline