/* $Id$ */
//
// Utility class to make simple Glauber type calculations
-// for SYMMTRIC collision geometries (AA):
+// for SYMMETRIC collision geometries (AA):
// Impact parameter, production points, reaction plane dependence
//
// The SimulateTrigger method can be used for simple MB and hard-process
// from AliRoot
#include "AliFastGlauber.h"
// from root
-#include <TStyle.h>
-#include <TH1F.h>
-#include <TH2F.h>
+#include <TCanvas.h>
#include <TF1.h>
#include <TF2.h>
+#include <TFile.h>
+#include <TH1F.h>
+#include <TH2F.h>
#include <TLegend.h>
-#include <TCanvas.h>
+#include <TMath.h>
#include <TRandom.h>
-#include <TFile.h>
-#include <TROOT.h>
-#include <Riostream.h>
+#include <TStyle.h>
ClassImp(AliFastGlauber)
TF2* AliFastGlauber::fgWKParticipants = NULL;
TF1* AliFastGlauber::fgWParticipants = NULL;
TF2* AliFastGlauber::fgWAlmondCurrent = NULL;
-TF2 AliFastGlauber::fgWAlmondFixedB[40];
+TF2* AliFastGlauber::fgWAlmondFixedB[40];
const Int_t AliFastGlauber::fgkMCInts = 100000;
-Int_t AliFastGlauber::fgCounter = 0;
-
-AliFastGlauber::AliFastGlauber() : fName()
-{
- // Default Constructor
- //
- fgCounter++;
- if(fgCounter>1)
- Error("AliFastGlauber","More than more instance (%d) is not supported, check your code!",fgCounter);
-
+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()
+{
+ // Default Constructor
// Defaults for Pb
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()
+{
+// 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()
{
- fgCounter--;
- // if(fgCounter==0) Reset();
+// Destructor
+ for(Int_t k=0; k<40; k++) delete fgWAlmondFixedB[k];
}
void AliFastGlauber::SetAuAuRhic()
//Set all parameters for RHIC
SetWoodSaxonParametersAu();
SetHardCrossSection();
- SetNNCrossSection(40);
+ SetNNCrossSection(42);
SetNucleus(197);
SetFileName("$(ALICE_ROOT)/FASTSIM/data/glauberAuAu.root");
}
fgWSb->SetParameter(2, fWSw);
fgWSb->SetParameter(3, fWSn);
- fgWSbz = new TF2("WSbz", WSbz, 0, fgBMax, 4);
+ fgWSbz = new TF2("WSbz", WSbz, 0, fgBMax, 0, fgBMax, 4);
fgWSbz->SetParameter(0, fWSr0);
fgWSbz->SetParameter(1, fWSd);
fgWSbz->SetParameter(2, fWSw);
for(Int_t k=0; k<40; k++) {
sprintf(almondName,"WAlmondFixedB%d",k);
fgWAlmondCurrent = (TF2*)ff->Get(almondName);
- new(&fgWAlmondFixedB[k]) TF2(*fgWAlmondCurrent);
+ fgWAlmondFixedB[k] = fgWAlmondCurrent;
}
delete ff;
}
fgWPathLength->SetParameter(2, 0); //Pathlength definition
}
-void AliFastGlauber::Reset()
+void AliFastGlauber::Reset() const
{
//
// Reset dynamic allocated formulas
if(fgWParticipants) delete fgWParticipants;
}
-void AliFastGlauber::DrawWSb()
+void AliFastGlauber::DrawWSb() const
{
//
// Draw Wood-Saxon Nuclear Density Function
c1->Update();
}
-void AliFastGlauber::DrawOverlap()
+void AliFastGlauber::DrawOverlap() const
{
//
// Draw Overlap Function
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);
+ 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}]");
fgWStaa->Draw("same");
}
-void AliFastGlauber::DrawParticipants()
+void AliFastGlauber::DrawParticipants() const
{
//
// Draw Number of Participants Npart
c3->Update();
}
-void AliFastGlauber::DrawThickness()
+void AliFastGlauber::DrawThickness() const
{
//
// Draw Thickness Function
fgWSta->Draw("same");
}
-void AliFastGlauber::DrawGeo()
+void AliFastGlauber::DrawGeo() const
{
//
// Draw Geometrical Cross-Section
c5->Update();
}
-void AliFastGlauber::DrawBinary()
+void AliFastGlauber::DrawBinary() const
{
//
// Draw Binary Cross-Section
c6->Update();
}
-void AliFastGlauber::DrawN()
+void AliFastGlauber::DrawN() const
{
//
// Draw Binaries per event (Ncoll)
c7->Update();
}
-void AliFastGlauber::DrawKernel(Double_t b)
+void AliFastGlauber::DrawKernel(Double_t b) const
{
//
// Draw Kernel
l1a->SetFillStyle(0);
l1a->SetBorderSize(0);
Char_t label[100];
- sprintf(label,"b = %.f fm",b);
+ sprintf(label,"b = %.1f fm",b);
l1a->AddEntry(fgWStarfi,label,"");
l1a->Draw();
c8->Update();
}
-void AliFastGlauber::DrawAlmond(Double_t b)
+void AliFastGlauber::DrawAlmond(Double_t b) const
{
//
// Draw Interaction Almond
l1a->SetFillStyle(0);
l1a->SetBorderSize(0);
Char_t label[100];
- sprintf(label,"b = %.f fm",b);
+ sprintf(label,"b = %.1f fm",b);
l1a->AddEntry(fgWAlmond,label,"");
l1a->Draw();
c9->Update();
}
-void AliFastGlauber::DrawEnergyDensity()
+void AliFastGlauber::DrawEnergyDensity() const
{
//
// Draw energy density
c10->Update();
}
-void AliFastGlauber::DrawPathLength0(Double_t b, Int_t iopt)
+void AliFastGlauber::DrawPathLength0(Double_t b, Int_t iopt) const
{
//
// Draw Path Length
fgWPathLength0->Draw("same");
}
-void AliFastGlauber::DrawPathLength(Double_t b , Int_t ni, Int_t iopt)
+void AliFastGlauber::DrawPathLength(Double_t b , Int_t ni, Int_t iopt) const
{
//
// Draw Path Length
fgWPathLength->Draw("same");
}
-void AliFastGlauber::DrawIntRadius(Double_t b)
+void AliFastGlauber::DrawIntRadius(Double_t b) const
{
//
// Draw Interaction Radius
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
// as a function of radius (xx)
//
- const Double_t xx = x[0]; //fm
- const Double_t r0 = par[0]; //fm
- const Double_t d = par[1]; //fm
- const Double_t w = par[2]; //no units
- const Double_t n = par[3]; //fm^-3 (used to normalize integral to one)
- const Double_t y = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
+ 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 = kn * (1.+kw*(kxx/kr0)*(kxx/kr0))/(1.+TMath::Exp((kxx-kr0)/kd));
return y; //fm^-3
}
-Double_t AliFastGlauber::WSbz(Double_t* x, Double_t* par)
+Double_t AliFastGlauber::WSbz(const Double_t* x, const Double_t* par)
{
//
// Wood Saxon Parameterisation
// as a function of z and b
//
- const Double_t bb = x[0]; //fm
- const Double_t zz = x[1]; //fm
- const Double_t r0 = par[0]; //fm
- const Double_t d = par[1]; //fm
- const Double_t w = par[2]; //no units
- const Double_t n = par[3]; //fm^-3 (used to normalize integral to one)
- const Double_t xx = TMath::Sqrt(bb*bb+zz*zz);
- const Double_t y = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
+ const Double_t kbb = x[0]; //fm
+ const Double_t kzz = x[1]; //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)
+ 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; //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
// as a function of z for fixed b
//
- const Double_t zz = x[0]; //fm
- const Double_t r0 = par[0]; //fm
- const Double_t d = par[1]; //fm
- const Double_t w = par[2]; //no units
- const Double_t n = par[3]; //fm^-3 (used to normalize integral to one)
- const Double_t bb = par[4]; //fm
- const Double_t xx = TMath::Sqrt(bb*bb+zz*zz);
- const Double_t y = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
+ const Double_t kzz = 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)
+ const Double_t kbb = par[4]; //fm
+ 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; //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
// as a function of b
//
- const Double_t b = x[0];
- fgWSz->SetParameter(4, b);
- const Double_t y = 2. * fgWSz->Integral(0., fgBMax);
+ const Double_t kb = x[0];
+ fgWSz->SetParameter(4, kb);
+ Double_t y = 2. * fgWSz->Integral(0., fgBMax);
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)
// as a function of r and phi
- const Double_t r1 = x[0];
- const Double_t phi = x[1];
- const Double_t b = par[0];
- const Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
- const Double_t y = r1 * fgWSta->Eval(r1) * fgWSta->Eval(r2);
+ 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; //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
// as a function of b
// (normalized to fA*fB)
//
- const Double_t b = x[0];
- const Double_t A = par[0];
- fgWStarfi->SetParameter(0, b);
+ const Double_t kb = x[0];
+ const Double_t ka = par[0];
+ fgWStarfi->SetParameter(0, kb);
// root integration seems to fail
/*
// MC Integration
//
Double_t y = 0;
+
+
for (Int_t i = 0; i < fgkMCInts; i++)
{
- const Double_t phi = TMath::Pi() * gRandom->Rndm();
- const Double_t b1 = fgBMax * gRandom->Rndm();
- y += fgWStarfi->Eval(b1, phi);
+
+ 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; //fm^-2
- y *= A * A * 0.1; //mbarn^-1
+ y *= ka * ka * 0.1; //mbarn^-1
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
// as a function of r and phi
//
- const Double_t r1 = x[0];
- const Double_t phi = x[1];
- const Double_t b = par[0]; //fm
- const Double_t sig = par[1]; //mbarn
- const Double_t A = par[2]; //mass number
- const Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
- const Double_t xsi = fgWSta->Eval(r2) * sig * 0.1;
+ const Double_t kr1 = x[0];
+ const Double_t kphi = x[1];
+ const Double_t kb = par[0]; //fm
+ const Double_t ksig = par[1]; //mbarn
+ const Double_t ka = par[2]; //mass number
+ 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; //no units
/*
- Double_t y=(1-TMath::Power((1-xsi),A))
+ Double_t y=(1-TMath::Power((1-xsi),aa))
*/
- Double_t a = A;
- Double_t sum = a * xsi;
+ Double_t a = ka;
+ Double_t sum = ka * kxsi;
Double_t y = sum;
- for (Int_t i = 1; i <= A; i++)
+ for (Int_t i = 1; i <= ka; i++)
{
a--;
- sum *= (-xsi) * a / Float_t(i+1);
+ sum *= (-kxsi) * a / Float_t(i+1);
y += sum;
}
- y = r1 * fgWSta->Eval(r1) * y;
+ y *= kr1 * fgWSta->Eval(kr1);
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
// a function of b
//
- const Double_t b = x[0];
- const Double_t sig = par[0]; //mbarn
- const Double_t A = par[1]; //mass number
- fgWKParticipants->SetParameter(0, b);
- fgWKParticipants->SetParameter(1, sig);
- fgWKParticipants->SetParameter(2, A);
+ const Double_t kb = x[0];
+ const Double_t ksig = par[0]; //mbarn
+ const Double_t ka = par[1]; //mass number
+ fgWKParticipants->SetParameter(0, kb);
+ fgWKParticipants->SetParameter(1, ksig);
+ fgWKParticipants->SetParameter(2, ka);
//
// MC Integration
Double_t y = 0;
for (Int_t i = 0; i < fgkMCInts; i++)
{
- Double_t phi = TMath::Pi() * gRandom->Rndm();
- Double_t b1 = fgBMax * gRandom->Rndm();
- y += fgWKParticipants->Eval(b1, phi);
+ const Double_t kphi = TMath::Pi() * gRandom->Rndm();
+ const Double_t kb1 = fgBMax * gRandom->Rndm();
+ y += fgWKParticipants->Eval(kb1, kphi);
}
- y *= 2. * A * 2. * TMath::Pi() * fgBMax / fgkMCInts;
+ y *= 2. * ka * 2. * TMath::Pi() * fgBMax / fgkMCInts;
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
// as a function of b
//
- const Double_t b = x[0]; //fm
- const Double_t sigNN = par[0]; //mbarn
- const Double_t taa = fgWStaa->Eval(b); //mbarn^-1
- Double_t y = 2. * TMath::Pi() * b * (1. - TMath::Exp(- sigNN * taa));
+ const Double_t kb = x[0]; //fm
+ const Double_t ksigNN = par[0]; //mbarn
+ const Double_t ktaa = fgWStaa->Eval(kb); //mbarn^-1
+ Double_t y = 2. * TMath::Pi() * kb * (1. - TMath::Exp(- ksigNN * ktaa));
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 collisions
+ // Number of binary hard collisions
// as a function of b
//
- const Double_t b = x[0]; //fm
- const Double_t sig = par[0]; //mbarn
- const Double_t taa = fgWStaa->Eval(b); //mbarn^-1
- const Double_t y = 2. * TMath::Pi() * b * sig * taa;
+ const Double_t kb = x[0]; //fm
+ const Double_t ksig = par[0]; //mbarn
+ const Double_t ktaa = fgWStaa->Eval(kb); //mbarn^-1
+ Double_t y = 2. * TMath::Pi() * kb * ksig * ktaa;
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
// as a function of b
- const Double_t b = x[0];
- const Double_t y = fgWSbinary->Eval(b)/fgWSgeo->Eval(b);
+ const Double_t kb = x[0];
+ Double_t y = fgWSbinary->Eval(kb)/fgWSgeo->Eval(kb);
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
// as a function of the impact parameter
//
- const Double_t b = x[0];
- const Double_t rA = par[0];
+ const Double_t kb = x[0];
+ const Double_t krA = par[0];
//
// Attention: area of transverse reaction zone in hard-sphere approximation !
- const Double_t rA2=rA*rA;
- const Double_t b2=b*b;
- Double_t saa = (TMath::Pi() - 2. * TMath::ASin(b/ 2./ rA)) * rA2
- - b * TMath::Sqrt(rA2 - b2/ 4.); //fm^2
- Double_t taa = fgWStaa->Eval(b); //mbarn^-1
- Double_t y=taa/saa*10;
+ 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.); //fm^2
+ const Double_t ktaa = fgWStaa->Eval(kb); //mbarn^-1
+ Double_t y=ktaa/ksaa*10;
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
// as a function of cartesian x,y.
//
- const Double_t b = par[0];
- const Double_t xx = x[0] + b/2.;
- const Double_t yy = x[1];
- const Double_t r1 = TMath::Sqrt(xx * xx + yy * yy);
- const Double_t phi = TMath::ATan2(yy,xx);
- const Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
+ 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));
//
// Interaction probability calculated as product of thicknesses
//
- const Double_t y = fgWSta->Eval(r1) * fgWSta->Eval(r2);
+ Double_t y = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
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
// at which interaction takes place
// as a function of radius
//
- const Double_t r = x[0];
- const Double_t b = par[0];
- fgWAlmond->SetParameter(0, b);
+ const Double_t kr = x[0];
+ const Double_t kb = par[0];
+ fgWAlmond->SetParameter(0, kb);
// Average over phi in small steps
- const Double_t dphi = 2. * TMath::Pi() / 100.;
+ 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 xx = r * TMath::Cos(phi);
- const Double_t yy = r * TMath::Sin(phi);
- y += fgWAlmond->Eval(xx,yy);
- phi += dphi;
+ const Double_t kxx = kr * TMath::Cos(phi);
+ const Double_t kyy = kr * TMath::Sin(phi);
+ y += fgWAlmond->Eval(kxx,kyy);
+ phi += kdphi;
} // phi loop
// Result multiplied by Jacobian (2 pi r)
- y *= 2. * TMath::Pi() * r / 100.;
+ y *= 2. * TMath::Pi() * kr / 100.;
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
// as a function of phi-direction
//
// Phi direction in Almond
- const Double_t phi0 = x[0];
- const Double_t b = par[0];
+ const Double_t kphi0 = x[0];
+ const Double_t kb = par[0];
// Path Length definition
- const Int_t iopt = Int_t(par[1]);
+ const Int_t kiopt = Int_t(par[1]);
// Step along radial direction phi
const Int_t kNp = 100; // Steps in r
//
// Transform into target frame
//
- const Double_t xx = r * TMath::Cos(phi0) + b / 2.;
- const Double_t yy = r * TMath::Sin(phi0);
- const Double_t phi = TMath::ATan2(yy, xx);
- const Double_t r1 = TMath::Sqrt(xx * xx + yy * yy);
+ 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);
// Radius in projectile frame
- const Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
- const Double_t y = fgWSta->Eval(r1) * fgWSta->Eval(r2);
+ 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 += y * r;
- w += y;
+ rw += ky * r;
+ w += ky;
r += kDr;
} // radial steps
Double_t y=0.;
- if (!iopt) // My length definition (is exact for hard disk)
- if(w) y= 2. * rw / w;
- else {
- const Double_t norm=fgWSta->Eval(0.01);
- if(norm) y = TMath::Sqrt(2. * rw * kDr /norm /norm);
+ if (!kiopt) { // My length definition (is exact for hard disk)
+ 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; //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
// Interaction point from random distribution
// as a function of the phi-direction
- const Double_t phi0 = x[0];
- const Double_t b = par[0];
- fgWAlmond->SetParameter(0, b);
- const Int_t kNpi = Int_t (par[1]); //Number of interactions
- const Int_t iopt = Int_t(par[2]); //Path Length definition
+ const Double_t kphi0 = x[0];
+ const Double_t kb = par[0];
+ fgWAlmond->SetParameter(0, kb);
+ const Int_t kNpi = Int_t (par[1]); //Number of interactions
+ const Int_t kiopt = Int_t(par[2]); //Path Length definition
//
// r-steps
Double_t x0, y0;
fgWAlmond->GetRandom2(x0, y0);
// Initial radius
- const Double_t r0 = TMath::Sqrt(x0 * x0 + y0 * y0);
- const Int_t nps = Int_t ((fgBMax - r0)/kDr) - 1;
+ const Double_t kr0 = TMath::Sqrt(x0*x0 + y0*y0);
+ const Int_t knps = Int_t ((fgBMax - kr0)/kDr) - 1;
// Radial steps
Double_t r = 0.;
- for (Int_t i = 0; (i < nps ); i++) {
+ for (Int_t i = 0; (i < knps ); i++) {
// Transform into target frame
- const Double_t xx = x0 + r * TMath::Cos(phi0) + b / 2.;
- const Double_t yy = y0 + r * TMath::Sin(phi0);
- const Double_t phi = TMath::ATan2(yy, xx);
- const Double_t r1 = TMath::Sqrt(xx * xx + yy * yy);
+ 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);
// Radius in projectile frame
- const Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
- const Double_t y = fgWSta->Eval(r1) * fgWSta->Eval(r2);
+ 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 += y * r;
- w += y;
+ rw += ky * r;
+ w += ky;
r += kDr;
} // steps
// Average over interactions
- if (!iopt) {
+ if (!kiopt) {
if(w) l += (2. * rw / w);
} else {
- const Double_t norm=fgWSta->Eval(0.01);
- if(norm) l+= 2. * rw * kDr / norm / norm;
+ const Double_t knorm=fgWSta->Eval(1e-4);
+ if(knorm) l+= 2. * rw * kDr / knorm / knorm;
}
} // interactions
Double_t ret=0;
- if (!iopt)
+ if (!kiopt)
ret= l / kNpi;
else
ret=TMath::Sqrt( l / kNpi);
return ret; //fm
}
-Double_t AliFastGlauber::CrossSection(Double_t b1, Double_t b2)
+Double_t AliFastGlauber::CrossSection(Double_t b1, Double_t b2) const
{
//
// Return the geometrical cross-section integrated from b1 to b2
return fgWSgeo->Integral(b1, b2)*10.; //mbarn
}
-Double_t AliFastGlauber::HardCrossSection(Double_t b1, Double_t b2)
+Double_t AliFastGlauber::HardCrossSection(Double_t b1, Double_t b2) const
{
//
// Return the hard cross-section integrated from b1 to b2
return fgWSbinary->Integral(b1, b2)*10.; //mbarn
}
-Double_t AliFastGlauber::FractionOfHardCrossSection(Double_t b1, Double_t b2)
+Double_t AliFastGlauber::FractionOfHardCrossSection(Double_t b1, Double_t b2) const
{
//
- // Return raction of hard cross-section integrated from b1 to b2
+ // Return fraction of hard cross-section integrated from b1 to b2
//
return fgWSbinary->Integral(b1, b2)/fgWSbinary->Integral(0., 100.);
}
-Double_t AliFastGlauber::Binaries(Double_t b)
+Double_t AliFastGlauber::NHard(const Double_t b1, const Double_t b2) const
+{
+ //
+ // Number of binary hard collisions
+ // as a function of b (nucl/ex/0302016 eq. 19)
+ //
+ 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
{
//
- // Return number of binary collisions normalized to 1 at b=0
+ // 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);
}
-Double_t AliFastGlauber::GetNumberofBinaries(Double_t b)
+Double_t AliFastGlauber::MeanOverlap(Double_t b1, Double_t b2)
+{
+//
+// Calculate the mean overlap for impact parameter range b1 .. 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)
+{
+//
+// Calculate the mean number of collisions per event for impact parameter range b1 .. 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
{
//
- // Return number of binary collisions at b
+ // 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);
}
-Double_t AliFastGlauber::Participants(Double_t b)
+Double_t AliFastGlauber::Participants(Double_t b) const
{
//
// 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));
}
-Double_t AliFastGlauber::GetNumberOfParticipants(Double_t b)
+Double_t AliFastGlauber::GetNumberOfParticipants(Double_t b) const
{
//
// 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));
}
+Double_t AliFastGlauber::GetNumberOfCollisions(Double_t b) const
+{
+ //
+ // Return the number of collisions for impact parameter b
+ //
+ if(b<1.e-4) b=1e-4;
+ return (fgWStaa->Eval(b)*fSigmaNN);
+}
+
+Double_t AliFastGlauber::GetNumberOfCollisionsPerEvent(Double_t b) const
+{
+ //
+ // Return the number of collisions per event (at least one collision)
+ // for impact parameter b
+ //
+ Double_t n = GetNumberOfCollisions(b);
+ if (n > 0.) {
+ return (n / (1. - TMath::Exp(- n)));
+ } else {
+ return (0.);
+ }
+}
+
void AliFastGlauber::SimulateTrigger(Int_t n)
{
//
// Gives back a random impact parameter, hard trigger probability and multiplicity
//
b = fgWSgeo->GetRandom();
- const Float_t mu = fgWSN->Eval(b);
- p = 1.-TMath::Exp(-mu);
- mult = 6000./fgWSN->Eval(1.) * mu;
+ 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)
//
// Gives back a random impact parameter bin, and hard trigger decission
//
- const Float_t b = fgWSgeo->GetRandom();
- const Float_t mu = fgWSN->Eval(b) * fSigmaHard;
- const Float_t p = 1.-TMath::Exp(-mu);
- if (b < 5.) {
+ 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 (b < 8.6) {
+ } else if (kb < 8.6) {
bin = 2;
- } else if (b < 11.2) {
+ } else if (kb < 11.2) {
bin = 3;
- } else if (b < 13.2) {
+ } else if (kb < 13.2) {
bin = 4;
- } else if (b < 15.0) {
+ } else if (kb < 15.0) {
bin = 5;
} else {
bin = 6;
}
hard = kFALSE;
- const Float_t r = gRandom->Rndm();
- if (r < p) hard = kTRUE;
+ const Float_t kr = gRandom->Rndm();
+ if (kr < kp) hard = kTRUE;
}
Double_t AliFastGlauber::GetRandomImpactParameter(Double_t bmin, Double_t bmax)
return b;
}
-void AliFastGlauber::StoreFunctions()
+void AliFastGlauber::StoreFunctions() const
{
//
// Store in file functions
//=================== Added by A. Dainese 11/02/04 ===========================
-void AliFastGlauber::PlotAlmonds()
-{
- //
- // Plot almonds for some impact parameters
- //
- 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::StoreAlmonds()
+void AliFastGlauber::StoreAlmonds() const
{
//
// Store in file
Double_t bLow=0.,bUp=0.;
Double_t xsecFr=0.;
- const Double_t norm=fgWSgeo->Integral(0.,100.);
+ const Double_t knorm=fgWSgeo->Integral(0.,100.);
while(xsecFr<xsecFrLow) {
- xsecFr = fgWSgeo->Integral(0.,bLow)/norm;
+ xsecFr = fgWSgeo->Integral(0.,bLow)/knorm;
bLow += 0.1;
}
bUp = bLow;
while(xsecFr<xsecFrUp) {
- xsecFr = fgWSgeo->Integral(0.,bUp)/norm;
+ xsecFr = fgWSgeo->Integral(0.,bUp)/knorm;
bUp += 0.1;
}
- Info("SetCentralityClass", "Centrality class: %4.2f-%4.2f; %4.1f < b < %4.1f fm\n",
+ 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;
}
b = fgWSbinary->GetRandom();
Int_t bin = 2*(Int_t)b;
if( (b-(Int_t)b) > 0.5) bin++;
- fgWAlmondCurrent = &fgWAlmondFixedB[bin];
+ fgWAlmondCurrent = fgWAlmondFixedB[bin];
return;
}
//
// Initial radius
- const Double_t r0 = TMath::Sqrt(x0 * x0 + y0 * y0);
- const Int_t nps = Int_t ((fgBMax - r0)/kDl) - 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.;
// Radial steps
- for (Int_t i = 0; i < nps; i++) {
+ for (Int_t i = 0; i < knps; i++) {
// Transform into target frame
- const Double_t xx = x0 + l * TMath::Cos(phi0) + b / 2.;
- const Double_t yy = y0 + l * TMath::Sin(phi0);
- const Double_t phi = TMath::ATan2(yy, xx);
- const Double_t r1 = TMath::Sqrt(xx * xx + yy * yy);
+ 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);
// Radius in projectile frame
- const Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
- const Double_t prodTATB = fgWSta->Eval(r1) * fgWSta->Eval(r2);
+ 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 += prodTATB * l * kDl;
- integral2 += prodTATB * kDl;
+ integral1 += kprodTATB * l * kDl;
+ integral2 += kprodTATB * kDl;
l += kDl;
} // steps
//
// Initial radius
- const Double_t r0 = TMath::Sqrt(x0 * x0 + y0 * y0);
- const Int_t nps = Int_t ((fgBMax - r0)/kDl) - 1;
- const Double_t prodTATBHalfMax = 0.5*fgWAlmondCurrent->Eval(0.,0.);
+ 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.);
// Radial steps
Double_t l = 0.;
Double_t integral = 0.;
- for (Int_t i = 0; i < nps; i++) {
+ for (Int_t i = 0; i < knps; i++) {
// Transform into target frame
- const Double_t xx = x0 + l * TMath::Cos(phi0) + b / 2.;
- const Double_t yy = y0 + l * TMath::Sin(phi0);
- const Double_t phi = TMath::ATan2(yy, xx);
- const Double_t r1 = TMath::Sqrt(xx * xx + yy * yy);
+ 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);
// Radius in projectile frame
- const Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
- const Double_t prodTATB = fgWSta->Eval(r1) * fgWSta->Eval(r2);
- if(prodTATB>prodTATBHalfMax) integral += kDl;
+ 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;
} // steps
Double_t ell = integral;
}
}
-void AliFastGlauber::GetLength(Double_t& ell,Double_t b)
+void AliFastGlauber::GetLengthAndPhi(Double_t& ell,Double_t& phi,Double_t b)
{
//
// Return length from random b, x0, y0, phi0
+ // Return also phi0
//
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::GetLengthsBackToBack(Double_t& ell1,Double_t& ell2,Double_t b)
+void AliFastGlauber::GetLength(Double_t& ell,Double_t b)
+{
+ //
+ // Return length from random b, x0, y0, phi0
+ //
+ Double_t phi;
+ GetLengthAndPhi(ell,phi,b);
+ return;
+}
+
+void AliFastGlauber::GetLengthsBackToBackAndPhi(Double_t& ell1,Double_t& ell2,Double_t &phi,Double_t b)
{
//
// Return 2 lengths back to back from random b, x0, y0, phi0
- //
+ // Return also phi0
+ //
Double_t x0,y0,phi0;
if(b<0.) GetRandomBHard(b);
GetRandomXY(x0,y0);
GetRandomPhi(phi0);
- const Double_t phi0plusPi = phi0+TMath::Pi();
+ const Double_t kphi0plusPi = phi0+TMath::Pi();
+ phi = phi0;
ell1 = CalculateLength(b,x0,y0,phi0);
- ell2 = CalculateLength(b,x0,y0,phi0plusPi);
+ ell2 = CalculateLength(b,x0,y0,kphi0plusPi);
+ return;
+}
+
+void AliFastGlauber::GetLengthsBackToBack(Double_t& ell1,Double_t& ell2,
+ Double_t b)
+{
+ //
+ // Return 2 lengths back to back from random b, x0, y0, phi0
+ //
+ Double_t phi;
+ GetLengthsBackToBackAndPhi(ell1,ell2,phi,b);
return;
}
return;
}
-void AliFastGlauber::PlotLengthDistr(Int_t n,Bool_t save,Char_t *fname)
+void AliFastGlauber::PlotLengthDistr(Int_t n,Bool_t save,const char *fname)
{
//
// Plot length distribution
//
Double_t ell;
- TH1F *hEll = new TH1F("hEll","Length distribution",16,-0.5,15.5);
+ TH1F *hEll = new TH1F("hEll","Length distribution",64,-0.5,15);
hEll->SetXTitle("Transverse path length, L [fm]");
hEll->SetYTitle("Probability");
hEll->SetFillColor(2);
return;
}
-void AliFastGlauber::PlotLengthB2BDistr(Int_t n,Bool_t save,Char_t *fname)
+void AliFastGlauber::PlotLengthB2BDistr(Int_t n,Bool_t save,const char *fname)
{
//
// Plot lengths back-to-back distributions
return;
}
+void AliFastGlauber::PlotAlmonds() const
+{
+ //
+ // Plot almonds for some impact parameters
+ //
+ 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;
+}
+
//=================== Added by A. Dainese 05/03/04 ===========================
void AliFastGlauber::CalculateI0I1(Double_t& integral0,Double_t& integral1,
Double_t b,Double_t x0,Double_t y0,
- Double_t phi0,Double_t ellCut)
+ Double_t phi0,Double_t ellCut) const
{
//
// Calculate integrals:
- // integral0 = \int_0^ellCut dl*(T_A*T_B)(x0+l*ux,y0+l*uy)
- // integral1 = \int_0^ellCut dl*l*(T_A*T_B)(x0+l*ux,y0+l*uy)
+ // integral0 = \int_0^ellCut dl*(T_A*T_B)(x0+l*ux,y0+l*uy)
+ // integral1 = \int_0^ellCut dl*l*(T_A*T_B)(x0+l*ux,y0+l*uy)
//
// for a parton with production point (x0,y0)
// and propagation direction (ux=cos(phi0),uy=sin(phi0))
const Double_t kDl = fgBMax/Double_t(kNp);
// Initial radius
- const Double_t r0 = TMath::Sqrt(x0 * x0 + y0 * y0);
- const Int_t nps = Int_t ((fgBMax - r0)/kDl) - 1;
+ const Double_t kr0 = TMath::Sqrt(x0 * x0 + y0 * y0);
+ const Int_t knps = Int_t ((fgBMax - kr0)/kDl) - 1;
// Radial steps
Double_t l = 0.;
integral0 = 0.;
integral1 = 0.;
Int_t i = 0;
- while((i < nps) && (l < ellCut)) {
+ while((i < knps) && (l < ellCut)) {
// Transform into target frame
- const Double_t xx = x0 + l * TMath::Cos(phi0) + b / 2.;
- const Double_t yy = y0 + l * TMath::Sin(phi0);
- const Double_t phi = TMath::ATan2(yy, xx);
- const Double_t r1 = TMath::Sqrt(xx * xx + yy * yy);
+ 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);
// Radius in projectile frame
- const Double_t r2 = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi));
- const Double_t prodTATB = fgWSta->Eval(r1) * fgWSta->Eval(r2);
- integral0 += prodTATB * kDl;
- integral1 += prodTATB * l * kDl;
+ 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++;
} // steps
return;
}
+void AliFastGlauber::GetI0I1AndPhi(Double_t& integral0,Double_t& integral1,
+ Double_t& phi,
+ Double_t ellCut,Double_t b)
+{
+ //
+ // Return I0 and I1 from random b, x0, y0, phi0
+ // Return also phi
+ //
+ 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)
{
//
// Return I0 and I1 from random b, x0, y0, phi0
//
+ 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)
+{
+ //
+ // Return 2 pairs of I0 and I1 back to back from random b, x0, y0, phi0
+ // Return also phi0
+ //
Double_t x0,y0,phi0;
if(b<0.) GetRandomBHard(b);
GetRandomXY(x0,y0);
GetRandomPhi(phi0);
- CalculateI0I1(integral0,integral1,b,x0,y0,phi0,ellCut);
+ 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::GetI0I1BackToBack(Double_t& integral01,Double_t& integral11,
- Double_t& integral02,Double_t& integral12,
- Double_t ellCut,Double_t b)
+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)
{
//
// Return 2 pairs of I0 and I1 back to back from random b, x0, y0, phi0
+ // Return also phi0
//
Double_t x0,y0,phi0;
if(b<0.) GetRandomBHard(b);
GetRandomXY(x0,y0);
GetRandomPhi(phi0);
- const Double_t phi0plusPi = phi0+TMath::Pi();
+ 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,phi0plusPi,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)
+{
+ //
+ // Return 2 pairs of I0 and I1 back to back from random b, x0, y0, phi0
+ //
+ Double_t phi;
+ GetI0I1BackToBackAndPhi(integral01,integral11,integral02,integral12,
+ phi,ellCut,b);
return;
}
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)
+{
+ //
+ // Returns I0 and I1 pairs for n partons with azimuthal angles phi[n]
+ // from random b, x0, y0 and return x0,y0
+ //
+ 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,Char_t *fname)
+ Bool_t save,const char *fname)
{
//
- // Plot length distribution
+ // Plot I0-I1 distribution
//
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",
- 100,0,0.001);
+ 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",
- 100,0,0.01);
+ 1000,0,0.01);
hI1->SetXTitle("I_{1} [fm^{-2}]");
hI1->SetYTitle("Probability");
hI1->SetFillColor(4);
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);
}
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);
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();
}
void AliFastGlauber::PlotI0I1B2BDistr(Int_t n,Double_t ellCut,
- Bool_t save,Char_t *fname)
+ Bool_t save,const char *fname)
{
//
- // Plot lengths back-to-back distributions
+ // Plot I0-I1 back-to-back distributions
//
Double_t i01,i11,i02,i12;
TH2F *hI0s = new TH2F("hI0s","I_{0}'s back-to-back",100,0,100,100,0,100);
}
return;
}
+
+AliFastGlauber& AliFastGlauber::operator=(const AliFastGlauber& rhs)
+{
+// Assignment operator
+ rhs.Copy(*this);
+ return *this;
+}
+
+void AliFastGlauber::Copy(TObject&) const
+{
+ //
+ // Copy
+ //
+ Fatal("Copy","Not implemented!\n");
+}
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff