#include "AliFastGlauber.h"
#include <TCanvas.h>
#include <TF1.h>
#include <TF2.h>
#include <TFile.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TLegend.h>
#include <TMath.h>
#include <TRandom.h>
#include <TStyle.h>
ClassImp(AliFastGlauber)
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::fgWStarfi = NULL;
TF2* AliFastGlauber::fgWAlmond = NULL;
TF1* AliFastGlauber::fgWStaa = NULL;
TF1* AliFastGlauber::fgWSgeo = NULL;
TF1* AliFastGlauber::fgWSbinary = NULL;
TF1* AliFastGlauber::fgWSN = NULL;
TF1* AliFastGlauber::fgWPathLength0 = NULL;
TF1* AliFastGlauber::fgWPathLength = NULL;
TF1* AliFastGlauber::fgWEnergyDensity = NULL;
TF1* AliFastGlauber::fgWIntRadius = NULL;
TF2* AliFastGlauber::fgWKParticipants = NULL;
TF1* AliFastGlauber::fgWParticipants = NULL;
TF2* AliFastGlauber::fgWAlmondCurrent = NULL;
TF2* AliFastGlauber::fgWAlmondFixedB[40];
const Int_t AliFastGlauber::fgkMCInts = 100000;
AliFastGlauber* AliFastGlauber::fgGlauber = NULL;
AliFastGlauber::AliFastGlauber():
fWSr0(0.),
fWSd(0.),
fWSw(0.),
fWSn(0.),
fSigmaHard(0.),
fSigmaNN(0.),
fA(0),
fBmin(0.),
fBmax(0.),
fEllDef(0),
fName()
{
SetMaxImpact();
SetLengthDefinition();
SetPbPbLHC();
fXY[0] = fXY[1] = 0;
fI0I1[0] = fI0I1[1] = 0;
}
AliFastGlauber::AliFastGlauber(const AliFastGlauber & gl)
:TObject(gl),
fWSr0(0.),
fWSd(0.),
fWSw(0.),
fWSn(0.),
fSigmaHard(0.),
fSigmaNN(0.),
fA(0),
fBmin(0.),
fBmax(0.),
fEllDef(0),
fName()
{
gl.Copy(*this);
fXY[0] = fXY[1] = 0;
fI0I1[0] = fI0I1[1] = 0;
}
AliFastGlauber* AliFastGlauber::Instance()
{
if (fgGlauber) {
return fgGlauber;
} else {
fgGlauber = new AliFastGlauber();
return fgGlauber;
}
}
AliFastGlauber::~AliFastGlauber()
{
for(Int_t k=0; k<40; k++) delete fgWAlmondFixedB[k];
}
void AliFastGlauber::SetAuAuRhic()
{
SetWoodSaxonParametersAu();
SetHardCrossSection();
SetNNCrossSection(42);
SetNucleus(197);
SetFileName("$(ALICE_ROOT)/FASTSIM/data/glauberAuAu.root");
}
void AliFastGlauber::SetPbPbLHC()
{
SetWoodSaxonParametersPb();
SetHardCrossSection();
SetNNCrossSection();
SetNucleus();
SetFileName();
}
void AliFastGlauber::Init(Int_t mode)
{
Reset();
fgWSb = new TF1("WSb", WSb, 0, fgBMax, 4);
fgWSb->SetParameter(0, fWSr0);
fgWSb->SetParameter(1, fWSd);
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);
fgWSbz->SetParameter(2, fWSw);
fgWSbz->SetParameter(3, fWSn);
fgWSz = new TF1("WSz", WSz, 0, fgBMax, 5);
fgWSz->SetParameter(0, fWSr0);
fgWSz->SetParameter(1, fWSd);
fgWSz->SetParameter(2, fWSw);
fgWSz->SetParameter(3, fWSn);
fgWSta = new TF1("WSta", WSta, 0., fgBMax, 0);
fgWStarfi = new TF2("WStarfi", WStarfi, 0., fgBMax, 0., TMath::Pi(), 1);
fgWStarfi->SetParameter(0, 0.);
fgWStarfi->SetNpx(200);
fgWStarfi->SetNpy(20);
fgWKParticipants = new TF2("WKParticipants", WKParticipants, 0., fgBMax, 0., TMath::Pi(), 3);
fgWKParticipants->SetParameter(0, 0.);
fgWKParticipants->SetParameter(1, fSigmaNN);
fgWKParticipants->SetParameter(2, fA);
fgWKParticipants->SetNpx(200);
fgWKParticipants->SetNpy(20);
if (!mode) {
fgWStaa = new TF1("WStaa", WStaa, 0., fgBMax, 1);
fgWStaa->SetNpx(100);
fgWStaa->SetParameter(0,fA);
fgWStaa->SetNpx(100);
fgWParticipants = new TF1("WParticipants", WParticipants, 0., fgBMax, 2);
fgWParticipants->SetParameter(0, fSigmaNN);
fgWParticipants->SetParameter(1, fA);
fgWParticipants->SetNpx(100);
} else {
Info("Init","Reading overlap function from file %s",fName.Data());
TFile* f = new TFile(fName.Data());
if(!f->IsOpen()){
Fatal("Init", "Could not open file %s",fName.Data());
}
fgWStaa = (TF1*) f->Get("WStaa");
fgWParticipants = (TF1*) f->Get("WParticipants");
delete f;
}
fgWEnergyDensity = new TF1("WEnergyDensity", WEnergyDensity, 0., 2. * fWSr0, 1);
fgWEnergyDensity->SetParameter(0, fWSr0 + 1.);
fgWSgeo = new TF1("WSgeo", WSgeo, 0., fgBMax, 1);
fgWSgeo->SetParameter(0,fSigmaNN);
fgWSgeo->SetNpx(100);
fgWSbinary = new TF1("WSbinary", WSbinary, 0., fgBMax, 1);
fgWSbinary->SetParameter(0, fSigmaHard);
fgWSbinary->SetNpx(100);
fgWSN = new TF1("WSN", WSN, 0.01, fgBMax, 1);
fgWSN->SetNpx(100);
fgWAlmond = new TF2("WAlmond", WAlmond, -fgBMax, fgBMax, -fgBMax, fgBMax, 1);
fgWAlmond->SetParameter(0, 0.);
fgWAlmond->SetNpx(200);
fgWAlmond->SetNpy(200);
if(mode==2) {
Info("Init","Reading interaction almonds from file: %s",fName.Data());
Char_t almondName[100];
TFile* ff = new TFile(fName.Data());
for(Int_t k=0; k<40; k++) {
snprintf(almondName,100, "WAlmondFixedB%d",k);
fgWAlmondCurrent = (TF2*)ff->Get(almondName);
fgWAlmondFixedB[k] = fgWAlmondCurrent;
}
delete ff;
}
fgWIntRadius = new TF1("WIntRadius", WIntRadius, 0., fgBMax, 1);
fgWIntRadius->SetParameter(0, 0.);
fgWPathLength0 = new TF1("WPathLength0", WPathLength0, -TMath::Pi(), TMath::Pi(), 2);
fgWPathLength0->SetParameter(0, 0.);
fgWPathLength0->SetParameter(1, 0.);
fgWPathLength = new TF1("WPathLength", WPathLength, -TMath::Pi(), TMath::Pi(), 3);
fgWPathLength->SetParameter(0, 0.);
fgWPathLength->SetParameter(1, 1000.);
fgWPathLength->SetParameter(2, 0);
}
void AliFastGlauber::Reset() const
{
if(fgWSb) delete fgWSb;
if(fgRWSb) delete fgRWSb;
if(fgWSbz) delete fgWSbz;
if(fgWSz) delete fgWSz;
if(fgWSta) delete fgWSta;
if(fgWStarfi) delete fgWStarfi;
if(fgWAlmond) delete fgWAlmond;
if(fgWStaa) delete fgWStaa;
if(fgWSgeo) delete fgWSgeo;
if(fgWSbinary) delete fgWSbinary;
if(fgWSN) delete fgWSN;
if(fgWPathLength0) delete fgWPathLength0;
if(fgWPathLength) delete fgWPathLength;
if(fgWEnergyDensity) delete fgWEnergyDensity;
if(fgWIntRadius) delete fgWIntRadius;
if(fgWKParticipants) delete fgWKParticipants;
if(fgWParticipants) delete fgWParticipants;
}
void AliFastGlauber::DrawWSb() const
{
TCanvas *c1 = new TCanvas("c1","Wood Saxon",400,10,600,700);
c1->cd();
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]");
h2f->GetYaxis()->SetNoExponent(kTRUE);
h2f->GetYaxis()->SetTitle("#rho [fm^{-3}]");
h2f->Draw();
fgWSb->Draw("same");
TLegend *l1a = new TLegend(0.45,0.6,.90,0.8);
l1a->SetFillStyle(0);
l1a->SetBorderSize(0);
Char_t label[100];
snprintf(label,100, "r_{0} = %.2f fm",fWSr0);
l1a->AddEntry(fgWSb,label,"");
snprintf(label,100, "d = %.2f fm",fWSd);
l1a->AddEntry(fgWSb,label,"");
snprintf(label,100, "n = %.2e fm^{-3}",fWSn);
l1a->AddEntry(fgWSb,label,"");
snprintf(label,100, "#omega = %.2f",fWSw);
l1a->AddEntry(fgWSb,label,"");
l1a->Draw();
c1->Update();
}
void AliFastGlauber::DrawOverlap() const
{
TCanvas *c2 = new TCanvas("c2","Overlap",400,10,600,700);
c2->cd();
Double_t max=fgWStaa->GetMaximum(0,fgBMax)*1.01;
TH2F *h2f=new TH2F("h2ftaa","Overlap function: T_{AB} [mbarn^{-1}]",2,0,fgBMax,2,0, max);
h2f->SetStats(0);
h2f->GetXaxis()->SetTitle("b [fm]");
h2f->GetYaxis()->SetTitle("T_{AB} [mbarn^{-1}]");
h2f->Draw();
fgWStaa->Draw("same");
}
void AliFastGlauber::DrawParticipants() const
{
TCanvas *c3 = new TCanvas("c3","Participants",400,10,600,700);
c3->cd();
Double_t max=fgWParticipants->GetMaximum(0,fgBMax)*1.01;
TH2F *h2f=new TH2F("h2fpart","Number of Participants",2,0,fgBMax,2,0,max);
h2f->SetStats(0);
h2f->GetXaxis()->SetTitle("b [fm]");
h2f->GetYaxis()->SetTitle("N_{part}");
h2f->Draw();
fgWParticipants->Draw("same");
TLegend *l1a = new TLegend(0.50,0.75,.90,0.9);
l1a->SetFillStyle(0);
l1a->SetBorderSize(0);
Char_t label[100];
snprintf(label,100, "#sigma^{inel.}_{NN} = %.1f mbarn",fSigmaNN);
l1a->AddEntry(fgWParticipants,label,"");
l1a->Draw();
c3->Update();
}
void AliFastGlauber::DrawThickness() const
{
TCanvas *c4 = new TCanvas("c4","Thickness",400,10,600,700);
c4->cd();
Double_t max=fgWSta->GetMaximum(0,fgBMax)*1.01;
TH2F *h2f=new TH2F("h2fta","Thickness function: T_{A} [fm^{-2}]",2,0,fgBMax,2,0,max);
h2f->SetStats(0);
h2f->GetXaxis()->SetTitle("b [fm]");
h2f->GetYaxis()->SetTitle("T_{A} [fm^{-2}]");
h2f->Draw();
fgWSta->Draw("same");
}
void AliFastGlauber::DrawGeo() const
{
TCanvas *c5 = new TCanvas("c5","Geometrical Cross-Section",400,10,600,700);
c5->cd();
Double_t max=fgWSgeo->GetMaximum(0,fgBMax)*1.01;
TH2F *h2f=new TH2F("h2fgeo","Differential Geometrical Cross-Section: d#sigma^{geo}_{AB}/db [fm]",2,0,fgBMax,2,0,max);
h2f->SetStats(0);
h2f->GetXaxis()->SetTitle("b [fm]");
h2f->GetYaxis()->SetTitle("d#sigma^{geo}_{AB}/db [fm]");
h2f->Draw();
fgWSgeo->Draw("same");
TLegend *l1a = new TLegend(0.10,0.8,.40,0.9);
l1a->SetFillStyle(0);
l1a->SetBorderSize(0);
Char_t label[100];
snprintf(label,100, "#sigma_{NN}^{inel.} = %.1f mbarn",fSigmaNN);
l1a->AddEntry(fgWSgeo,label,"");
l1a->Draw();
c5->Update();
}
void AliFastGlauber::DrawBinary() const
{
TCanvas *c6 = new TCanvas("c6","Binary Cross-Section",400,10,600,700);
c6->cd();
Double_t max=fgWSbinary->GetMaximum(0,fgBMax)*1.01;
TH2F *h2f=new TH2F("h2fbinary","Differential Binary Cross-Section: #sigma^{hard}_{NN} dT_{AB}/db [fm]",2,0,fgBMax,2,0,max);
h2f->SetStats(0);
h2f->GetXaxis()->SetTitle("b [fm]");
h2f->GetYaxis()->SetTitle("d#sigma^{hard}_{AB}/db [fm]");
h2f->Draw();
fgWSbinary->Draw("same");
TLegend *l1a = new TLegend(0.50,0.8,.90,0.9);
l1a->SetFillStyle(0);
l1a->SetBorderSize(0);
Char_t label[100];
snprintf(label,100, "#sigma_{NN}^{hard} = %.1f mbarn",fSigmaHard);
l1a->AddEntry(fgWSb,label,"");
l1a->Draw();
c6->Update();
}
void AliFastGlauber::DrawN() const
{
TCanvas *c7 = new TCanvas("c7","Binaries per event",400,10,600,700);
c7->cd();
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]");
h2f->GetYaxis()->SetTitle("N_{coll}");
h2f->Draw();
fgWSN->Draw("same");
TLegend *l1a = new TLegend(0.50,0.75,.90,0.9);
l1a->SetFillStyle(0);
l1a->SetBorderSize(0);
Char_t label[100];
snprintf(label,100, "#sigma^{hard}_{NN} = %.1f mbarn",fSigmaHard);
l1a->AddEntry(fgWSN,label,"");
snprintf(label,100, "#sigma^{inel.}_{NN} = %.1f mbarn",fSigmaNN);
l1a->AddEntry(fgWSN,label,"");
l1a->Draw();
c7->Update();
}
void AliFastGlauber::DrawKernel(Double_t b) const
{
TCanvas *c8 = new TCanvas("c8","Kernel",400,10,600,700);
c8->cd();
fgWStarfi->SetParameter(0, b);
TH2F *h2f=new TH2F("h2fkernel","Kernel of Overlap function: d^{2}T_{AB}/dr/d#phi [fm^{-3}]",2,0,fgBMax,2,0,TMath::Pi());
h2f->SetStats(0);
h2f->GetXaxis()->SetTitle("r [fm]");
h2f->GetYaxis()->SetTitle("#phi [rad]");
h2f->Draw();
fgWStarfi->Draw("same");
TLegend *l1a = new TLegend(0.65,0.8,.90,0.9);
l1a->SetFillStyle(0);
l1a->SetBorderSize(0);
Char_t label[100];
snprintf(label, 100, "b = %.1f fm",b);
l1a->AddEntry(fgWStarfi,label,"");
l1a->Draw();
c8->Update();
}
void AliFastGlauber::DrawAlmond(Double_t b) const
{
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,-fgBMax, fgBMax, 2, -fgBMax, fgBMax);
h2f->SetStats(0);
h2f->GetXaxis()->SetTitle("x [fm]");
h2f->GetYaxis()->SetTitle("y [fm]");
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];
snprintf(label, 100, "b = %.1f fm",b);
l1a->AddEntry(fgWAlmond,label,"");
l1a->Draw();
c9->Update();
}
void AliFastGlauber::DrawEnergyDensity() const
{
TCanvas *c10 = new TCanvas("c10","Energy Density",400, 10, 600, 700);
c10->cd();
fgWEnergyDensity->SetMinimum(0.);
Double_t max=fgWEnergyDensity->GetMaximum(0,fgWEnergyDensity->GetParameter(0))*1.01;
TH2F *h2f=new TH2F("h2fenergydens","Energy density",2,0,fgBMax,2,0,max);
h2f->SetStats(0);
h2f->GetXaxis()->SetTitle("b [fm]");
h2f->GetYaxis()->SetTitle("fm^{-4}");
h2f->Draw();
fgWEnergyDensity->Draw("same");
c10->Update();
}
void AliFastGlauber::DrawPathLength0(Double_t b, Int_t iopt) const
{
TCanvas *c11 = new TCanvas("c11","Path Length",400,10,600,700);
c11->cd();
fgWPathLength0->SetParameter(0, b);
fgWPathLength0->SetParameter(1, Double_t(iopt));
fgWPathLength0->SetMinimum(0.);
fgWPathLength0->SetMaximum(10.);
TH2F *h2f=new TH2F("h2fpathlength0","Path length",2,-TMath::Pi(), TMath::Pi(),2,0,10.);
h2f->SetStats(0);
h2f->GetXaxis()->SetTitle("#phi [rad]");
h2f->GetYaxis()->SetTitle("l [fm]");
h2f->Draw();
fgWPathLength0->Draw("same");
}
void AliFastGlauber::DrawPathLength(Double_t b , Int_t ni, Int_t iopt) const
{
TCanvas *c12 = new TCanvas("c12","Path Length",400,10,600,700);
c12->cd();
fgWAlmond->SetParameter(0, b);
fgWPathLength->SetParameter(0, b);
fgWPathLength->SetParameter(1, Double_t (ni));
fgWPathLength->SetParameter(2, Double_t (iopt));
fgWPathLength->SetMinimum(0.);
fgWPathLength->SetMaximum(10.);
TH2F *h2f=new TH2F("h2fpathlength","Path length",2,-TMath::Pi(), TMath::Pi(),2,0,10.);
h2f->SetStats(0);
h2f->GetXaxis()->SetTitle("#phi [rad]");
h2f->GetYaxis()->SetTitle("l [fm]");
h2f->Draw();
fgWPathLength->Draw("same");
}
void AliFastGlauber::DrawIntRadius(Double_t b) const
{
TCanvas *c13 = new TCanvas("c13","Interaction Radius",400,10,600,700);
c13->cd();
fgWIntRadius->SetParameter(0, b);
fgWIntRadius->SetMinimum(0);
Double_t max=fgWIntRadius->GetMaximum(0,fgBMax)*1.01;
TH2F *h2f=new TH2F("h2fintradius","Interaction Density",2,0.,fgBMax,2,0,max);
h2f->SetStats(0);
h2f->GetXaxis()->SetTitle("r [fm]");
h2f->GetYaxis()->SetTitle("[fm^{-3}]");
h2f->Draw();
fgWIntRadius->Draw("same");
}
Double_t AliFastGlauber::WSb(const Double_t* x, const Double_t* par)
{
const Double_t kxx = x[0];
const Double_t kr0 = par[0];
const Double_t kd = par[1];
const Double_t kw = par[2];
const Double_t kn = par[3];
Double_t y = kn * (1.+kw*(kxx/kr0)*(kxx/kr0))/(1.+TMath::Exp((kxx-kr0)/kd));
return y;
}
Double_t AliFastGlauber::RWSb(const Double_t* x, const Double_t* par)
{
const Double_t kxx = x[0];
const Double_t kr0 = par[0];
const Double_t kd = par[1];
const Double_t kw = par[2];
const Double_t kn = par[3];
Double_t y = kxx * kxx * kn * (1.+kw*(kxx/kr0)*(kxx/kr0))/(1.+TMath::Exp((kxx-kr0)/kd));
return y;
}
Double_t AliFastGlauber::WSbz(const Double_t* x, const Double_t* par)
{
const Double_t kbb = x[0];
const Double_t kzz = x[1];
const Double_t kr0 = par[0];
const Double_t kd = par[1];
const Double_t kw = par[2];
const Double_t kn = par[3];
const Double_t kxx = TMath::Sqrt(kbb*kbb+kzz*kzz);
Double_t y = kn * (1.+kw*(kxx/kr0)*(kxx/kr0))/(1.+TMath::Exp((kxx-kr0)/kd));
return y;
}
Double_t AliFastGlauber::WSz(const Double_t* x, const Double_t* par)
{
const Double_t kzz = x[0];
const Double_t kr0 = par[0];
const Double_t kd = par[1];
const Double_t kw = par[2];
const Double_t kn = par[3];
const Double_t kbb = par[4];
const Double_t kxx = TMath::Sqrt(kbb*kbb+kzz*kzz);
Double_t y = kn * (1.+kw*(kxx/kr0)*(kxx/kr0))/(1.+TMath::Exp((kxx-kr0)/kd));
return y;
}
Double_t AliFastGlauber::WSta(const Double_t* x, const Double_t* )
{
const Double_t kb = x[0];
fgWSz->SetParameter(4, kb);
Double_t y = 2. * fgWSz->Integral(0., fgBMax);
return y;
}
Double_t AliFastGlauber::WStarfi(const Double_t* x, const Double_t* par)
{
const Double_t kr1 = x[0];
const Double_t kphi = x[1];
const Double_t kb = par[0];
const Double_t kr2 = TMath::Sqrt(kr1*kr1 + kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
Double_t y = kr1 * fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
return y;
}
Double_t AliFastGlauber::WStaa(const Double_t* x, const Double_t* par)
{
const Double_t kb = x[0];
const Double_t ka = par[0];
fgWStarfi->SetParameter(0, kb);
Double_t y = 0;
for (Int_t i = 0; i < fgkMCInts; i++)
{
const Double_t kphi = TMath::Pi() * gRandom->Rndm();
const Double_t kb1 = fgBMax * gRandom->Rndm();
y += fgWStarfi->Eval(kb1, kphi);
}
y *= 2. * TMath::Pi() * fgBMax / fgkMCInts;
y *= ka * ka * 0.1;
return y;
}
Double_t AliFastGlauber::WKParticipants(const Double_t* x, const Double_t* par)
{
const Double_t kr1 = x[0];
const Double_t kphi = x[1];
const Double_t kb = par[0];
const Double_t ksig = par[1];
const Double_t ka = par[2];
const Double_t kr2 = TMath::Sqrt(kr1*kr1 +kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
const Double_t kxsi = fgWSta->Eval(kr2) * ksig * 0.1;
Double_t a = ka;
Double_t sum = ka * kxsi;
Double_t y = sum;
for (Int_t i = 1; i <= ka; i++)
{
a--;
sum *= (-kxsi) * a / Float_t(i+1);
y += sum;
}
y *= kr1 * fgWSta->Eval(kr1);
return y;
}
Double_t AliFastGlauber::WParticipants(const Double_t* x, const Double_t* par)
{
const Double_t kb = x[0];
const Double_t ksig = par[0];
const Double_t ka = par[1];
fgWKParticipants->SetParameter(0, kb);
fgWKParticipants->SetParameter(1, ksig);
fgWKParticipants->SetParameter(2, ka);
Double_t y = 0;
for (Int_t i = 0; i < fgkMCInts; i++)
{
const Double_t kphi = TMath::Pi() * gRandom->Rndm();
const Double_t kb1 = fgBMax * gRandom->Rndm();
y += fgWKParticipants->Eval(kb1, kphi);
}
y *= 2. * ka * 2. * TMath::Pi() * fgBMax / fgkMCInts;
return y;
}
Double_t AliFastGlauber::WSgeo(const Double_t* x, const Double_t* par)
{
const Double_t kb = x[0];
const Double_t ksigNN = par[0];
const Double_t ktaa = fgWStaa->Eval(kb);
Double_t y = 2. * TMath::Pi() * kb * (1. - TMath::Exp(- ksigNN * ktaa));
return y;
}
Double_t AliFastGlauber::WSbinary(const Double_t* x, const Double_t* par)
{
const Double_t kb = x[0];
const Double_t ksig = par[0];
const Double_t ktaa = fgWStaa->Eval(kb);
Double_t y = 2. * TMath::Pi() * kb * ksig * ktaa;
return y;
}
Double_t AliFastGlauber::WSN(const Double_t* x, const Double_t* )
{
const Double_t kb = x[0];
Double_t y = fgWSbinary->Eval(kb)/fgWSgeo->Eval(kb);
return y;
}
Double_t AliFastGlauber::WEnergyDensity(const Double_t* x, const Double_t* par)
{
const Double_t kb = x[0];
const Double_t krA = par[0];
const Double_t krA2=krA*krA;
const Double_t kb2=kb*kb;
const Double_t ksaa = (TMath::Pi() - 2. * TMath::ASin(kb/ 2./ krA)) * krA2
- kb * TMath::Sqrt(krA2 - kb2/ 4.);
const Double_t ktaa = fgWStaa->Eval(kb);
Double_t y=ktaa/ksaa*10;
return y;
}
Double_t AliFastGlauber::WAlmond(const Double_t* x, const Double_t* par)
{
const Double_t kb = par[0];
const Double_t kxx = x[0] + kb/2.;
const Double_t kyy = x[1];
const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
const Double_t kphi = TMath::ATan2(kyy,kxx);
const Double_t kr2 = TMath::Sqrt(kr1*kr1 + kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
Double_t y = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
return y;
}
Double_t AliFastGlauber::WIntRadius(const Double_t* x, const Double_t* par)
{
const Double_t kr = x[0];
const Double_t kb = par[0];
fgWAlmond->SetParameter(0, kb);
const Double_t kdphi = 2. * TMath::Pi() / 100.;
Double_t phi = 0.;
Double_t y = 0.;
for (Int_t i = 0; i < 100; i++) {
const Double_t kxx = kr * TMath::Cos(phi);
const Double_t kyy = kr * TMath::Sin(phi);
y += fgWAlmond->Eval(kxx,kyy);
phi += kdphi;
}
y *= 2. * TMath::Pi() * kr / 100.;
return y;
}
Double_t AliFastGlauber::WPathLength0(const Double_t* x, const Double_t* par)
{
const Double_t kphi0 = x[0];
const Double_t kb = par[0];
const Int_t kiopt = Int_t(par[1]);
const Int_t kNp = 100;
const Double_t kDr = fgBMax/kNp;
Double_t r = 0.;
Double_t rw = 0.;
Double_t w = 0.;
for (Int_t i = 0; i < kNp; i++) {
const Double_t kxx = r * TMath::Cos(kphi0) + kb / 2.;
const Double_t kyy = r * TMath::Sin(kphi0);
const Double_t kphi = TMath::ATan2(kyy, kxx);
const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
const Double_t kr2 = TMath::Sqrt(kr1*kr1 + kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
const Double_t ky = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
rw += ky * r;
w += ky;
r += kDr;
}
Double_t y=0.;
if (!kiopt) {
if(w) y= 2. * rw / w;
} else {
const Double_t knorm=fgWSta->Eval(1e-4);
if(knorm) y = TMath::Sqrt(2. * rw * kDr / knorm / knorm);
}
return y;
}
Double_t AliFastGlauber::WPathLength(const Double_t* x, const Double_t* par)
{
const Double_t kphi0 = x[0];
const Double_t kb = par[0];
fgWAlmond->SetParameter(0, kb);
const Int_t kNpi = Int_t (par[1]);
const Int_t kiopt = Int_t(par[2]);
const Int_t kNp = 100;
const Double_t kDr = fgBMax/Double_t(kNp);
Double_t l = 0.;
for (Int_t in = 0; in < kNpi; in ++) {
Double_t rw = 0.;
Double_t w = 0.;
Double_t x0, y0;
fgWAlmond->GetRandom2(x0, y0);
const Double_t kr0 = TMath::Sqrt(x0*x0 + y0*y0);
const Int_t knps = Int_t ((fgBMax - kr0)/kDr) - 1;
Double_t r = 0.;
for (Int_t i = 0; (i < knps ); i++) {
const Double_t kxx = x0 + r * TMath::Cos(kphi0) + kb / 2.;
const Double_t kyy = y0 + r * TMath::Sin(kphi0);
const Double_t kphi = TMath::ATan2(kyy, kxx);
const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
const Double_t kr2 = TMath::Sqrt(kr1*kr1 + kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
const Double_t ky = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
rw += ky * r;
w += ky;
r += kDr;
}
if (!kiopt) {
if(w) l += (2. * rw / w);
} else {
const Double_t knorm=fgWSta->Eval(1e-4);
if(knorm) l+= 2. * rw * kDr / knorm / knorm;
}
}
Double_t ret=0;
if (!kiopt)
ret= l / kNpi;
else
ret=TMath::Sqrt( l / kNpi);
return ret;
}
Double_t AliFastGlauber::CrossSection(Double_t b1, Double_t b2) const
{
return fgWSgeo->Integral(b1, b2)*10.;
}
Double_t AliFastGlauber::HardCrossSection(Double_t b1, Double_t b2) const
{
return fgWSbinary->Integral(b1, b2)*10.;
}
Double_t AliFastGlauber::FractionOfHardCrossSection(Double_t b1, Double_t b2) const
{
return fgWSbinary->Integral(b1, b2)/fgWSbinary->Integral(0., 100.);
}
Double_t AliFastGlauber::NHard(const Double_t b1, const Double_t b2) const
{
const Double_t kshard=HardCrossSection(b1,b2);
const Double_t ksgeo=CrossSection(b1,b2);
if(ksgeo>0)
return kshard/ksgeo;
else return -1;
}
Double_t AliFastGlauber::Binaries(Double_t b) const
{
if(b < 1.e-4) b = 1e-4;
return fgWSN->Eval(b)/fgWSN->Eval(1e-4);
}
Double_t AliFastGlauber::MeanOverlap(Double_t b1, Double_t b2)
{
Double_t sum = 0.;
Double_t sumc = 0.;
Double_t b = b1;
while (b < b2-0.005) {
Double_t nc = GetNumberOfCollisions(b);
sum += 10. * fgWStaa->Eval(b) * fgWSgeo->Eval(b) * 0.01 / (1. - TMath::Exp(-nc));
sumc += 10. * fgWSgeo->Eval(b) * 0.01;
b += 0.01;
}
return (sum / CrossSection(b1, b2));
}
Double_t AliFastGlauber::MeanNumberOfCollisionsPerEvent(Double_t b1, Double_t b2)
{
Double_t sum = 0.;
Double_t sumc = 0.;
Double_t b = b1;
while (b < b2-0.005) {
Double_t nc = GetNumberOfCollisions(b);
sum += nc / (1. - TMath::Exp(-nc)) * 10. * fgWSgeo->Eval(b) * 0.01;
sumc += 10. * fgWSgeo->Eval(b) * 0.01;
b += 0.01;
}
return (sum / CrossSection(b1, b2));
}
Double_t AliFastGlauber::GetNumberOfBinaries(Double_t b) const
{
if(b<1.e-4) b=1e-4;
return fgWSN->Eval(b);
}
Double_t AliFastGlauber::Participants(Double_t b) const
{
if(b<1.e-4) b=1e-4;
return (fgWParticipants->Eval(b)/fgWParticipants->Eval(1e-4));
}
Double_t AliFastGlauber::GetNumberOfParticipants(Double_t b) const
{
if(b<1.e-4) b=1e-4;
return (fgWParticipants->Eval(b));
}
Double_t AliFastGlauber::GetNumberOfCollisions(Double_t b) const
{
if(b<1.e-4) b=1e-4;
return (fgWStaa->Eval(b)*fSigmaNN);
}
Double_t AliFastGlauber::GetNumberOfCollisionsPerEvent(Double_t b) const
{
Double_t n = GetNumberOfCollisions(b);
if (n > 0.) {
return (n / (1. - TMath::Exp(- n)));
} else {
return (0.);
}
}
void AliFastGlauber::SimulateTrigger(Int_t n)
{
TH1F* mbtH = new TH1F("mbtH", "MB Trigger b-Distribution", 100, 0., 20.);
TH1F* hdtH = new TH1F("hdtH", "Hard Trigger b-Distribution", 100, 0., 20.);
TH1F* mbmH = new TH1F("mbmH", "MB Trigger Multiplicity Distribution", 100, 0., 8000.);
TH1F* hdmH = new TH1F("hdmH", "Hard Trigger Multiplicity Distribution", 100, 0., 8000.);
mbtH->SetXTitle("b [fm]");
hdtH->SetXTitle("b [fm]");
mbmH->SetXTitle("Multiplicity");
hdmH->SetXTitle("Multiplicity");
TCanvas *c0 = new TCanvas("c0","Trigger Simulation",400,10,600,700);
c0->Divide(2,1);
TCanvas *c1 = new TCanvas("c1","Trigger Simulation",400,10,600,700);
c1->Divide(1,2);
Init(1);
for (Int_t iev = 0; iev < n; iev++)
{
Float_t b, p, mult;
GetRandom(b, p, mult);
mbtH->Fill(b,1.);
hdtH->Fill(b, p);
mbmH->Fill(mult, 1.);
hdmH->Fill(mult, p);
c0->cd(1);
mbtH->Draw();
c0->cd(2);
hdtH->Draw();
c0->Update();
c1->cd(1);
mbmH->Draw();
c1->cd(2);
hdmH->Draw();
c1->Update();
}
}
void AliFastGlauber::GetRandom(Float_t& b, Float_t& p, Float_t& mult)
{
b = fgWSgeo->GetRandom();
const Float_t kmu = fgWSN->Eval(b);
p = 1.-TMath::Exp(-kmu);
mult = 6000./fgWSN->Eval(1.) * kmu;
}
void AliFastGlauber::GetRandom(Int_t& bin, Bool_t& hard)
{
const Float_t kb = fgWSgeo->GetRandom();
const Float_t kmu = fgWSN->Eval(kb) * fSigmaHard;
const Float_t kp = 1.-TMath::Exp(-kmu);
if (kb < 5.) {
bin = 1;
} else if (kb < 8.6) {
bin = 2;
} else if (kb < 11.2) {
bin = 3;
} else if (kb < 13.2) {
bin = 4;
} else if (kb < 15.0) {
bin = 5;
} else {
bin = 6;
}
hard = kFALSE;
const Float_t kr = gRandom->Rndm();
if (kr < kp) hard = kTRUE;
}
Double_t AliFastGlauber::GetRandomImpactParameter(Double_t bmin, Double_t bmax)
{
Float_t b = -1.;
while(b < bmin || b > bmax)
b = fgWSgeo->GetRandom();
return b;
}
void AliFastGlauber::StoreFunctions() const
{
TFile* ff = new TFile(fName.Data(),"recreate");
fgWStaa->Write("WStaa");
fgWParticipants->Write("WParticipants");
ff->Close();
return;
}
void AliFastGlauber::StoreAlmonds() const
{
Char_t almondName[100];
TFile* ff = new TFile(fName.Data(),"update");
for(Int_t k=0; k<40; 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);
fgWAlmond->Write(almondName);
}
ff->Close();
return;
}
void AliFastGlauber::SetCentralityClass(Double_t xsecFrLow,Double_t xsecFrUp)
{
if(xsecFrLow>1. || xsecFrUp>1. || xsecFrLow>xsecFrUp) {
Error("SetCentralityClass", "Please set 0 <= xsecFrLow <= xsecFrUp <= 1\n");
return;
}
Double_t bLow=0.,bUp=0.;
Double_t xsecFr=0.;
const Double_t knorm=fgWSgeo->Integral(0.,100.);
while(xsecFr<xsecFrLow) {
xsecFr = fgWSgeo->Integral(0.,bLow)/knorm;
bLow += 0.1;
}
bUp = bLow;
while(xsecFr<xsecFrUp) {
xsecFr = fgWSgeo->Integral(0.,bUp)/knorm;
bUp += 0.1;
}
Info("SetCentralityClass", "Centrality class: %4.2f-%4.2f; %4.1f < b < %4.1f fm",
xsecFrLow,xsecFrUp,bLow,bUp);
fgWSbinary->SetRange(bLow,bUp);
fBmin=bLow;
fBmax=bUp;
return;
}
void AliFastGlauber::GetRandomBHard(Double_t& b)
{
b = fgWSbinary->GetRandom();
Int_t bin = 2*(Int_t)b;
if( (b-(Int_t)b) > 0.5) bin++;
fgWAlmondCurrent = fgWAlmondFixedB[bin];
return;
}
void AliFastGlauber::GetRandomXY(Double_t& x,Double_t& y)
{
fgWAlmondCurrent->GetRandom2(x,y);
return;
}
void AliFastGlauber::GetRandomPhi(Double_t& phi)
{
phi = 2.*TMath::Pi()*gRandom->Rndm();
return;
}
Double_t AliFastGlauber::CalculateLength(Double_t b,Double_t x0,Double_t y0,Double_t phi0)
{
const Int_t kNp = 100;
const Double_t kDl = fgBMax/Double_t(kNp);
if(fEllDef==1) {
const Double_t kr0 = TMath::Sqrt(x0*x0 + y0*y0);
const Int_t knps = Int_t ((fgBMax - kr0)/kDl) - 1;
Double_t l = 0.;
Double_t integral1 = 0.;
Double_t integral2 = 0.;
for (Int_t i = 0; i < knps; i++) {
const Double_t kxx = x0 + l * TMath::Cos(phi0) + b / 2.;
const Double_t kyy = y0 + l * TMath::Sin(phi0);
const Double_t kphi = TMath::ATan2(kyy, kxx);
const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
const Double_t kr2 = TMath::Sqrt(kr1*kr1 + b*b - 2.*kr1*b*TMath::Cos(kphi));
const Double_t kprodTATB = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
integral1 += kprodTATB * l * kDl;
integral2 += kprodTATB * kDl;
l += kDl;
}
Double_t ell=0.;
if(integral2)
ell = (2. * integral1 / integral2);
return ell;
} else if(fEllDef==2) {
const Double_t kr0 = TMath::Sqrt(x0*x0 + y0*y0);
const Int_t knps = Int_t ((fgBMax - kr0)/kDl) - 1;
const Double_t kprodTATBHalfMax = 0.5*fgWAlmondCurrent->Eval(0.,0.);
Double_t l = 0.;
Double_t integral = 0.;
for (Int_t i = 0; i < knps; i++) {
const Double_t kxx = x0 + l * TMath::Cos(phi0) + b / 2.;
const Double_t kyy = y0 + l * TMath::Sin(phi0);
const Double_t kphi = TMath::ATan2(kyy, kxx);
const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
const Double_t kr2 = TMath::Sqrt(kr1*kr1 + b*b - 2.*kr1*b*TMath::Cos(kphi));
const Double_t kprodTATB = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
if(kprodTATB>kprodTATBHalfMax) integral += kDl;
l += kDl;
}
Double_t ell = integral;
return ell;
} else {
Error("CalculateLength","Wrong length definition setting: %d !\n",fEllDef);
return -1.;
}
}
void AliFastGlauber::GetLengthAndPhi(Double_t& ell,Double_t& phi,Double_t b)
{
Double_t x0,y0,phi0;
if(b<0.) GetRandomBHard(b);
GetRandomXY(x0,y0);
GetRandomPhi(phi0);
phi = phi0;
ell = CalculateLength(b,x0,y0,phi0);
return;
}
void AliFastGlauber::GetLength(Double_t& ell,Double_t b)
{
Double_t phi;
GetLengthAndPhi(ell,phi,b);
return;
}
void AliFastGlauber::GetLengthsBackToBackAndPhi(Double_t& ell1,Double_t& ell2,Double_t &phi,Double_t b)
{
Double_t x0,y0,phi0;
if(b<0.) GetRandomBHard(b);
GetRandomXY(x0,y0);
GetRandomPhi(phi0);
const Double_t kphi0plusPi = phi0+TMath::Pi();
phi = phi0;
ell1 = CalculateLength(b,x0,y0,phi0);
ell2 = CalculateLength(b,x0,y0,kphi0plusPi);
return;
}
void AliFastGlauber::GetLengthsBackToBack(Double_t& ell1,Double_t& ell2,
Double_t b)
{
Double_t phi;
GetLengthsBackToBackAndPhi(ell1,ell2,phi,b);
return;
}
void AliFastGlauber::GetLengthsForPythia(Int_t n,Double_t* const phi,Double_t* ell, Double_t b)
{
Double_t x0, y0;
if(b < 0.) GetRandomBHard(b);
GetRandomXY(x0,y0);
for(Int_t i = 0; i< n; i++) ell[i] = CalculateLength(b,x0,y0,phi[i]);
return;
}
void AliFastGlauber::PlotBDistr(Int_t n)
{
Double_t b;
TH1F *hB = new TH1F("hB","dN/db",100,0,fgBMax);
hB->SetXTitle("b [fm]");
hB->SetYTitle("dN/db [a.u.]");
hB->SetFillColor(3);
for(Int_t i=0; i<n; i++) {
GetRandomBHard(b);
hB->Fill(b);
}
TCanvas *cB = new TCanvas("cB","Impact parameter distribution",0,0,500,500);
cB->cd();
hB->Draw();
return;
}
void AliFastGlauber::PlotLengthDistr(Int_t n,Bool_t save,const char *fname)
{
Double_t ell;
TH1F *hEll = new TH1F("hEll","Length distribution",64,-0.5,15);
hEll->SetXTitle("Transverse path length, L [fm]");
hEll->SetYTitle("Probability");
hEll->SetFillColor(2);
for(Int_t i=0; i<n; i++) {
GetLength(ell);
hEll->Fill(ell);
}
hEll->Scale(1/(Double_t)n);
TCanvas *cL = new TCanvas("cL","Length distribution",0,0,500,500);
cL->cd();
hEll->Draw();
if(save) {
TFile *f = new TFile(fname,"recreate");
hEll->Write();
f->Close();
}
return;
}
void AliFastGlauber::PlotLengthB2BDistr(Int_t n,Bool_t save,const char *fname)
{
Double_t ell1,ell2;
TH2F *hElls = new TH2F("hElls","Lengths back-to-back",100,0,15,100,0,15);
hElls->SetXTitle("Transverse path length, L [fm]");
hElls->SetYTitle("Transverse path length, L [fm]");
for(Int_t i=0; i<n; i++) {
GetLengthsBackToBack(ell1,ell2);
hElls->Fill(ell1,ell2);
}
hElls->Scale(1/(Double_t)n);
TCanvas *cLs = new TCanvas("cLs","Length back-to-back distribution",0,0,500,500);
gStyle->SetPalette(1,0);
cLs->cd();
hElls->Draw("col,Z");
if(save) {
TFile *f = new TFile(fname,"recreate");
hElls->Write();
f->Close();
}
return;
}
void AliFastGlauber::PlotAlmonds() const
{
TCanvas *c = new TCanvas("c","Almonds",0,0,500,500);
gStyle->SetPalette(1,0);
c->Divide(2,2);
c->cd(1);
fgWAlmondFixedB[0]->Draw("cont1");
c->cd(2);
fgWAlmondFixedB[10]->Draw("cont1");
c->cd(3);
fgWAlmondFixedB[20]->Draw("cont1");
c->cd(4);
fgWAlmondFixedB[30]->Draw("cont1");
return;
}
void AliFastGlauber::CalculateI0I1(Double_t& integral0,Double_t& integral1,
Double_t b,Double_t x0,Double_t y0,
Double_t phi0,Double_t ellCut) const
{
const Int_t kNp = 100;
const Double_t kDl = fgBMax/Double_t(kNp);
const Double_t kr0 = TMath::Sqrt(x0 * x0 + y0 * y0);
const Int_t knps = Int_t ((fgBMax - kr0)/kDl) - 1;
Double_t l = 0.;
integral0 = 0.;
integral1 = 0.;
Int_t i = 0;
while((i < knps) && (l < ellCut)) {
const Double_t kxx = x0 + l * TMath::Cos(phi0) + b / 2.;
const Double_t kyy = y0 + l * TMath::Sin(phi0);
const Double_t kphi = TMath::ATan2(kyy, kxx);
const Double_t kr1 = TMath::Sqrt(kxx*kxx + kyy*kyy);
const Double_t kr2 = TMath::Sqrt(kr1*kr1 + b*b - 2.*kr1*b*TMath::Cos(kphi));
const Double_t kprodTATB = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
integral0 += kprodTATB * kDl;
integral1 += kprodTATB * l * kDl;
l += kDl;
i++;
}
return;
}
void AliFastGlauber::GetI0I1AndPhi(Double_t& integral0,Double_t& integral1,
Double_t& phi,
Double_t ellCut,Double_t b)
{
Double_t x0,y0,phi0;
if(b<0.) GetRandomBHard(b);
GetRandomXY(x0,y0);
GetRandomPhi(phi0);
phi = phi0;
CalculateI0I1(integral0,integral1,b,x0,y0,phi0,ellCut);
return;
}
void AliFastGlauber::GetI0I1(Double_t& integral0,Double_t& integral1,
Double_t ellCut,Double_t b)
{
Double_t phi;
GetI0I1AndPhi(integral0,integral1,phi,ellCut,b);
return;
}
void AliFastGlauber::GetI0I1BackToBackAndPhi(Double_t& integral01,Double_t& integral11,
Double_t& integral02,Double_t& integral12,
Double_t& phi,
Double_t ellCut,Double_t b)
{
Double_t x0,y0,phi0;
if(b<0.) GetRandomBHard(b);
GetRandomXY(x0,y0);
GetRandomPhi(phi0);
phi = phi0;
const Double_t kphi0plusPi = phi0+TMath::Pi();
CalculateI0I1(integral01,integral11,b,x0,y0,phi0,ellCut);
CalculateI0I1(integral02,integral12,b,x0,y0,kphi0plusPi,ellCut);
return;
}
void AliFastGlauber::GetI0I1BackToBackAndPhiAndXY(Double_t& integral01,Double_t& integral11,
Double_t& integral02,Double_t& integral12,
Double_t& phi,Double_t &x,Double_t &y,
Double_t ellCut,Double_t b)
{
Double_t x0,y0,phi0;
if(b<0.) GetRandomBHard(b);
GetRandomXY(x0,y0);
GetRandomPhi(phi0);
phi = phi0; x=x0; y=y0;
const Double_t kphi0plusPi = phi0+TMath::Pi();
CalculateI0I1(integral01,integral11,b,x0,y0,phi0,ellCut);
CalculateI0I1(integral02,integral12,b,x0,y0,kphi0plusPi,ellCut);
return;
}
void AliFastGlauber::GetI0I1BackToBack(Double_t& integral01,Double_t& integral11,
Double_t& integral02,Double_t& integral12,
Double_t ellCut,Double_t b)
{
Double_t phi;
GetI0I1BackToBackAndPhi(integral01,integral11,integral02,integral12,
phi,ellCut,b);
return;
}
void AliFastGlauber::GetI0I1ForPythia(Int_t n,Double_t* phi,
Double_t* integral0,Double_t* integral1,
Double_t ellCut,Double_t b)
{
Double_t x0,y0;
if(b<0.) GetRandomBHard(b);
GetRandomXY(x0,y0);
for(Int_t i=0; i<n; i++)
CalculateI0I1(integral0[i],integral1[i],b,x0,y0,phi[i],ellCut);
return;
}
void AliFastGlauber::GetI0I1ForPythiaAndXY(Int_t n,Double_t* phi,
Double_t* integral0,Double_t* integral1,
Double_t &x,Double_t& y,
Double_t ellCut,Double_t b)
{
Double_t x0,y0;
if(b<0.) GetRandomBHard(b);
GetRandomXY(x0,y0);
for(Int_t i=0; i<n; i++)
CalculateI0I1(integral0[i],integral1[i],b,x0,y0,phi[i],ellCut);
x=x0;
y=y0;
return;
}
void AliFastGlauber::PlotI0I1Distr(Int_t n,Double_t ellCut,
Bool_t save,const char *fname)
{
Double_t i0,i1;
TH2F *hI0I1s = new TH2F("hI0I1s","I_{0} versus I_{1}",1000,0,0.001,1000,0,0.01);
hI0I1s->SetXTitle("I_{0} [fm^{-3}]");
hI0I1s->SetYTitle("I_{1} [fm^{-2}]");
TH1F *hI0 = new TH1F("hI0","I_{0} = #hat{q}L / k",
1000,0,0.001);
hI0->SetXTitle("I_{0} [fm^{-3}]");
hI0->SetYTitle("Probability");
hI0->SetFillColor(3);
TH1F *hI1 = new TH1F("hI1","I_{1} = #omega_{c} / k",
1000,0,0.01);
hI1->SetXTitle("I_{1} [fm^{-2}]");
hI1->SetYTitle("Probability");
hI1->SetFillColor(4);
TH1F *h2 = new TH1F("h2","2 I_{1}^{2}/I_{0} = R / k",
100,0,0.02);
h2->SetXTitle("2 I_{1}^{2}/I_{0} [fm^{-1}]");
h2->SetYTitle("Probability");
h2->SetFillColor(2);
TH1F *h3 = new TH1F("h3","2 I_{1}/I_{0} = L",
100,0,15);
h3->SetXTitle("2 I_{1}/I_{0} [fm]");
h3->SetYTitle("Probability");
h3->SetFillColor(5);
TH1F *h4 = new TH1F("h4","I_{0}^{2}/(2 I_{1}) = #hat{q} / k",
100,0,0.00015);
h4->SetXTitle("I_{0}^{2}/(2 I_{1}) [fm^{-4}]");
h4->SetYTitle("Probability");
h4->SetFillColor(7);
for(Int_t i=0; i<n; i++) {
GetI0I1(i0,i1,ellCut);
hI0I1s->Fill(i0,i1);
hI0->Fill(i0);
hI1->Fill(i1);
h2->Fill(2.*i1*i1/i0);
h3->Fill(2.*i1/i0);
h4->Fill(i0*i0/2./i1);
}
hI0->Scale(1/(Double_t)n);
hI1->Scale(1/(Double_t)n);
h2->Scale(1/(Double_t)n);
h3->Scale(1/(Double_t)n);
h4->Scale(1/(Double_t)n);
hI0I1s->Scale(1/(Double_t)n);
TCanvas *cI0I1 = new TCanvas("cI0I1","I0 and I1",0,0,900,700);
cI0I1->Divide(3,2);
cI0I1->cd(1);
hI0->Draw();
cI0I1->cd(2);
hI1->Draw();
cI0I1->cd(3);
h2->Draw();
cI0I1->cd(4);
h3->Draw();
cI0I1->cd(5);
h4->Draw();
cI0I1->cd(6);
gStyle->SetPalette(1,0);
hI0I1s->Draw("col,Z");
if(save) {
TFile *f = new TFile(fname,"recreate");
hI0I1s->Write();
hI0->Write();
hI1->Write();
h2->Write();
h3->Write();
h4->Write();
f->Close();
}
return;
}
void AliFastGlauber::PlotI0I1B2BDistr(Int_t n,Double_t ellCut,
Bool_t save,const char *fname)
{
Double_t i01,i11,i02,i12;
TH2F *hI0s = new TH2F("hI0s","I_{0}'s back-to-back",100,0,100,100,0,100);
hI0s->SetXTitle("I_{0} [fm^{-3}]");
hI0s->SetYTitle("I_{0} [fm^{-3}]");
TH2F *hI1s = new TH2F("hI1s","I_{1}'s back-to-back",100,0,100,100,0,100);
hI1s->SetXTitle("I_{1} [fm^{-2}]");
hI1s->SetYTitle("I_{1} [fm^{-2}]");
for(Int_t i=0; i<n; i++) {
GetI0I1BackToBack(i01,i11,i02,i12,ellCut);
hI0s->Fill(i01,i02);
hI1s->Fill(i11,i12);
}
hI0s->Scale(1/(Double_t)n);
hI1s->Scale(1/(Double_t)n);
TCanvas *cI0I1s = new TCanvas("cI0I1s","I0 and I1 back-to-back distributions",0,0,800,400);
gStyle->SetPalette(1,0);
cI0I1s->Divide(2,1);
cI0I1s->cd(1);
hI0s->Draw("col,Z");
cI0I1s->cd(2);
hI1s->Draw("col,Z");
if(save) {
TFile *f = new TFile(fname,"recreate");
hI0s->Write();
hI1s->Write();
f->Close();
}
return;
}
AliFastGlauber& AliFastGlauber::operator=(const AliFastGlauber& rhs)
{
rhs.Copy(*this);
return *this;
}
void AliFastGlauber::Copy(TObject&) const
{
Fatal("Copy","Not implemented!\n");
}