// from AliRoot
#include "AliFastGlauber.h"
// from root
-#include <Riostream.h>
#include <TCanvas.h>
#include <TF1.h>
#include <TF2.h>
#include <TH2F.h>
#include <TLegend.h>
#include <TMath.h>
-#include <TROOT.h>
#include <TRandom.h>
#include <TStyle.h>
Float_t AliFastGlauber::fgBMax = 0.;
TF1* AliFastGlauber::fgWSb = NULL;
+TF1* AliFastGlauber::fgRWSb = NULL;
TF2* AliFastGlauber::fgWSbz = NULL;
TF1* AliFastGlauber::fgWSz = NULL;
TF1* AliFastGlauber::fgWSta = NULL;
TF2* AliFastGlauber::fgWAlmondCurrent = NULL;
TF2* AliFastGlauber::fgWAlmondFixedB[40];
const Int_t AliFastGlauber::fgkMCInts = 100000;
-Int_t AliFastGlauber::fgCounter = 0;
+AliFastGlauber* AliFastGlauber::fgGlauber = NULL;
+
AliFastGlauber::AliFastGlauber():
fWSr0(0.),
fName()
{
// Default Constructor
- fgCounter++;
- if(fgCounter>1)
- Error("AliFastGlauber","More than one instance (%d) is not supported, check your code!",fgCounter);
-
// Defaults for Pb
SetMaxImpact();
SetLengthDefinition();
SetPbPbLHC();
+ fXY[0] = fXY[1] = 0;
+ fI0I1[0] = fI0I1[1] = 0;
}
AliFastGlauber::AliFastGlauber(const AliFastGlauber & gl)
{
// Copy constructor
gl.Copy(*this);
+ fXY[0] = fXY[1] = 0;
+ fI0I1[0] = fI0I1[1] = 0;
+}
+
+AliFastGlauber* AliFastGlauber::Instance()
+{
+// Set random number generator
+ if (fgGlauber) {
+ return fgGlauber;
+ } else {
+ fgGlauber = new AliFastGlauber();
+ return fgGlauber;
+ }
}
AliFastGlauber::~AliFastGlauber()
{
// Destructor
- fgCounter--;
for(Int_t k=0; k<40; k++) delete fgWAlmondFixedB[k];
}
fgWSb->SetParameter(2, fWSw);
fgWSb->SetParameter(3, fWSn);
+ fgRWSb = new TF1("RWSb", RWSb, 0, fgBMax, 4);
+ fgRWSb->SetParameter(0, fWSr0);
+ fgRWSb->SetParameter(1, fWSd);
+ fgRWSb->SetParameter(2, fWSw);
+ fgRWSb->SetParameter(3, fWSn);
+
fgWSbz = new TF2("WSbz", WSbz, 0, fgBMax, 0, fgBMax, 4);
fgWSbz->SetParameter(0, fWSr0);
fgWSbz->SetParameter(1, fWSd);
//
// Hard collisions per event
//
- fgWSN = new TF1("WSN", WSN, 0., fgBMax, 1);
+ fgWSN = new TF1("WSN", WSN, 0.01, fgBMax, 1);
fgWSN->SetNpx(100);
//
Char_t almondName[100];
TFile* ff = new TFile(fName.Data());
for(Int_t k=0; k<40; k++) {
- sprintf(almondName,"WAlmondFixedB%d",k);
+ snprintf(almondName,100, "WAlmondFixedB%d",k);
fgWAlmondCurrent = (TF2*)ff->Get(almondName);
fgWAlmondFixedB[k] = fgWAlmondCurrent;
}
// in case init is called twice
if(fgWSb) delete fgWSb;
+ if(fgRWSb) delete fgRWSb;
if(fgWSbz) delete fgWSbz;
if(fgWSz) delete fgWSz;
if(fgWSta) delete fgWSta;
//
TCanvas *c1 = new TCanvas("c1","Wood Saxon",400,10,600,700);
c1->cd();
- Double_t max=fgWSb->GetMaximum(0,fgBMax)*1.01;
+ Double_t max=fgWSb->GetMaximum(0.,fgBMax)*1.01;
TH2F *h2f=new TH2F("h2fwsb","Wood Saxon: #rho(r) = n (1-#omega(r/r_{0})^2)/(1+exp((r-r_{0})/d)) [fm^{-3}]",2,0,fgBMax,2,0,max);
h2f->SetStats(0);
h2f->GetXaxis()->SetTitle("r [fm]");
l1a->SetFillStyle(0);
l1a->SetBorderSize(0);
Char_t label[100];
- sprintf(label,"r_{0} = %.2f fm",fWSr0);
+ snprintf(label,100, "r_{0} = %.2f fm",fWSr0);
l1a->AddEntry(fgWSb,label,"");
- sprintf(label,"d = %.2f fm",fWSd);
+ snprintf(label,100, "d = %.2f fm",fWSd);
l1a->AddEntry(fgWSb,label,"");
- sprintf(label,"n = %.2e fm^{-3}",fWSn);
+ snprintf(label,100, "n = %.2e fm^{-3}",fWSn);
l1a->AddEntry(fgWSb,label,"");
- sprintf(label,"#omega = %.2f",fWSw);
+ snprintf(label,100, "#omega = %.2f",fWSw);
l1a->AddEntry(fgWSb,label,"");
l1a->Draw();
c1->Update();
l1a->SetFillStyle(0);
l1a->SetBorderSize(0);
Char_t label[100];
- sprintf(label,"#sigma^{inel.}_{NN} = %.1f mbarn",fSigmaNN);
+ snprintf(label,100, "#sigma^{inel.}_{NN} = %.1f mbarn",fSigmaNN);
l1a->AddEntry(fgWParticipants,label,"");
l1a->Draw();
c3->Update();
l1a->SetFillStyle(0);
l1a->SetBorderSize(0);
Char_t label[100];
- sprintf(label,"#sigma_{NN}^{inel.} = %.1f mbarn",fSigmaNN);
+ snprintf(label,100, "#sigma_{NN}^{inel.} = %.1f mbarn",fSigmaNN);
l1a->AddEntry(fgWSgeo,label,"");
l1a->Draw();
c5->Update();
l1a->SetFillStyle(0);
l1a->SetBorderSize(0);
Char_t label[100];
- sprintf(label,"#sigma_{NN}^{hard} = %.1f mbarn",fSigmaHard);
+ snprintf(label,100, "#sigma_{NN}^{hard} = %.1f mbarn",fSigmaHard);
l1a->AddEntry(fgWSb,label,"");
l1a->Draw();
c6->Update();
//
TCanvas *c7 = new TCanvas("c7","Binaries per event",400,10,600,700);
c7->cd();
- Double_t max=fgWSN->GetMaximum(0,fgBMax)*1.01;
+ Double_t max=fgWSN->GetMaximum(0.01,fgBMax)*1.01;
TH2F *h2f=new TH2F("h2fhardcols","Number of hard collisions: T_{AB} #sigma^{hard}_{NN}/#sigma_{AB}^{geo}",2,0,fgBMax,2,0,max);
h2f->SetStats(0);
h2f->GetXaxis()->SetTitle("b [fm]");
l1a->SetFillStyle(0);
l1a->SetBorderSize(0);
Char_t label[100];
- sprintf(label,"#sigma^{hard}_{NN} = %.1f mbarn",fSigmaHard);
+ snprintf(label,100, "#sigma^{hard}_{NN} = %.1f mbarn",fSigmaHard);
l1a->AddEntry(fgWSN,label,"");
- sprintf(label,"#sigma^{inel.}_{NN} = %.1f mbarn",fSigmaNN);
+ snprintf(label,100, "#sigma^{inel.}_{NN} = %.1f mbarn",fSigmaNN);
l1a->AddEntry(fgWSN,label,"");
l1a->Draw();
c7->Update();
l1a->SetFillStyle(0);
l1a->SetBorderSize(0);
Char_t label[100];
- sprintf(label,"b = %.1f fm",b);
+ snprintf(label, 100, "b = %.1f fm",b);
l1a->AddEntry(fgWStarfi,label,"");
l1a->Draw();
c8->Update();
TCanvas *c9 = new TCanvas("c9","Almond",400,10,600,700);
c9->cd();
fgWAlmond->SetParameter(0, b);
- TH2F *h2f=new TH2F("h2falmond","Interaction Almond [fm^{-4}]",2,0,fgBMax,2,0,fgBMax);
+ TH2F *h2f=new TH2F("h2falmond","Interaction Almond [fm^{-4}]",2,-fgBMax, fgBMax, 2, -fgBMax, fgBMax);
h2f->SetStats(0);
h2f->GetXaxis()->SetTitle("x [fm]");
h2f->GetYaxis()->SetTitle("y [fm]");
- h2f->Draw();
- fgWAlmond->Draw("same");
+ h2f->Draw("");
+ gStyle->SetPalette(1);
+ fgWAlmond->Draw("colzsame");
TLegend *l1a = new TLegend(0.65,0.8,.90,0.9);
l1a->SetFillStyle(0);
l1a->SetBorderSize(0);
Char_t label[100];
- sprintf(label,"b = %.1f fm",b);
+ snprintf(label, 100, "b = %.1f fm",b);
l1a->AddEntry(fgWAlmond,label,"");
l1a->Draw();
c9->Update();
fgWIntRadius->Draw("same");
}
-Double_t AliFastGlauber::WSb(Double_t* x, Double_t* par)
+Double_t AliFastGlauber::WSb(const Double_t* x, const Double_t* par)
{
//
// Woods-Saxon Parameterisation
return y; //fm^-3
}
-Double_t AliFastGlauber::WSbz(Double_t* x, Double_t* par)
+Double_t AliFastGlauber::RWSb(const Double_t* x, const Double_t* par)
+{
+ //
+ // Woods-Saxon Parameterisation
+ // as a function of radius (xx)
+ // times r**2
+ const Double_t kxx = x[0]; //fm
+ const Double_t kr0 = par[0]; //fm
+ const Double_t kd = par[1]; //fm
+ const Double_t kw = par[2]; //no units
+ const Double_t kn = par[3]; //fm^-3 (used to normalize integral to one)
+ Double_t y = kxx * kxx * kn * (1.+kw*(kxx/kr0)*(kxx/kr0))/(1.+TMath::Exp((kxx-kr0)/kd));
+
+ return y; //fm^-1
+}
+
+Double_t AliFastGlauber::WSbz(const Double_t* x, const Double_t* par)
{
//
// Wood Saxon Parameterisation
return y; //fm^-3
}
-Double_t AliFastGlauber::WSz(Double_t* x, Double_t* par)
+Double_t AliFastGlauber::WSz(const Double_t* x, const Double_t* par)
{
//
// Wood Saxon Parameterisation
return y; //fm^-3
}
-Double_t AliFastGlauber::WSta(Double_t* x, Double_t* /*par*/)
+Double_t AliFastGlauber::WSta(const Double_t* x, const Double_t* /*par*/)
{
//
// Thickness function T_A
return y; //fm^-2
}
-Double_t AliFastGlauber::WStarfi(Double_t* x, Double_t* par)
+Double_t AliFastGlauber::WStarfi(const Double_t* x, const Double_t* par)
{
//
// Kernel for overlap function: T_A(s)*T_A(s-b)
return y; //fm^-3
}
-Double_t AliFastGlauber::WStaa(Double_t* x, Double_t* par)
+Double_t AliFastGlauber::WStaa(const Double_t* x, const Double_t* par)
{
//
// Overlap function
return y;
}
-Double_t AliFastGlauber::WKParticipants(Double_t* x, Double_t* par)
+Double_t AliFastGlauber::WKParticipants(const Double_t* x, const Double_t* par)
{
//
// Kernel for number of participants
return y; //fm^-1
}
-Double_t AliFastGlauber::WParticipants(Double_t* x, Double_t* par)
+Double_t AliFastGlauber::WParticipants(const Double_t* x, const Double_t* par)
{
//
// Number of Participants as
return y; //no units
}
-Double_t AliFastGlauber::WSgeo(Double_t* x, Double_t* par)
+Double_t AliFastGlauber::WSgeo(const Double_t* x, const Double_t* par)
{
//
// Geometrical Cross-Section
return y; //fm
}
-Double_t AliFastGlauber::WSbinary(Double_t* x, Double_t* par)
+Double_t AliFastGlauber::WSbinary(const Double_t* x, const Double_t* par)
{
//
// Number of binary hard collisions
return y; //fm
}
-Double_t AliFastGlauber::WSN(Double_t* x, Double_t* /*par*/)
+Double_t AliFastGlauber::WSN(const Double_t* x, const Double_t* /*par*/)
{
//
// Number of hard processes per event
return y; //no units
}
-Double_t AliFastGlauber::WEnergyDensity(Double_t* x, Double_t* par)
+Double_t AliFastGlauber::WEnergyDensity(const Double_t* x, const Double_t* par)
{
//
// Initial energy density
return y; //fm^-4
}
-Double_t AliFastGlauber::WAlmond(Double_t* x, Double_t* par)
+Double_t AliFastGlauber::WAlmond(const Double_t* x, const Double_t* par)
{
//
// Almond shaped interaction region
return y; //fm^-4
}
-Double_t AliFastGlauber::WIntRadius(Double_t* x, Double_t* par)
+Double_t AliFastGlauber::WIntRadius(const Double_t* x, const Double_t* par)
{
//
// Average interaction density over radius
return y; //fm^-3
}
-Double_t AliFastGlauber::WPathLength0(Double_t* x, Double_t* par)
+Double_t AliFastGlauber::WPathLength0(const Double_t* x, const Double_t* par)
{
//
// Path Length as a function of phi
return y; //fm
}
-Double_t AliFastGlauber::WPathLength(Double_t* x, Double_t* par)
+Double_t AliFastGlauber::WPathLength(const Double_t* x, const Double_t* par)
{
//
// Path Length as a function of phi
return fgWSbinary->Integral(b1, b2)/fgWSbinary->Integral(0., 100.);
}
-Double_t AliFastGlauber::NHard(Double_t b1, Double_t b2) const
+Double_t AliFastGlauber::NHard(const Double_t b1, const Double_t b2) const
{
//
// Number of binary hard collisions
//
// Return number of binary hard collisions normalized to 1 at b=0
//
- if(b==0) b=1e-4;
+ if(b < 1.e-4) b = 1e-4;
return fgWSN->Eval(b)/fgWSN->Eval(1e-4);
}
//
// Return number of binary hard collisions at b
//
- if(b==0) b=1e-4;
+ if(b<1.e-4) b=1e-4;
return fgWSN->Eval(b);
}
//
// Return the number of participants normalized to 1 at b=0
//
- if(b==0) b=1e-4;
+ if(b<1.e-4) b=1e-4;
return (fgWParticipants->Eval(b)/fgWParticipants->Eval(1e-4));
}
//
// Return the number of participants for impact parameter b
//
- if(b==0) b=1e-4;
+ if(b<1.e-4) b=1e-4;
return (fgWParticipants->Eval(b));
}
//
// Return the number of collisions for impact parameter b
//
- if(b==0) b=1e-4;
+ if(b<1.e-4) b=1e-4;
return (fgWStaa->Eval(b)*fSigmaNN);
}
Char_t almondName[100];
TFile* ff = new TFile(fName.Data(),"update");
for(Int_t k=0; k<40; k++) {
- sprintf(almondName,"WAlmondFixedB%d",k);
+ snprintf(almondName, 100, "WAlmondFixedB%d",k);
Double_t b = 0.25+k*0.5;
Info("StoreAlmonds"," b = %f\n",b);
fgWAlmond->SetParameter(0,b);
return;
}
-void AliFastGlauber::GetLengthsForPythia(Int_t n,Double_t* phi,Double_t* ell, Double_t b)
+void AliFastGlauber::GetLengthsForPythia(Int_t n,Double_t* const phi,Double_t* ell, Double_t b)
{
//
// Returns lenghts for n partons with azimuthal angles phi[n]