Added Eventplane Dependence in dPhi Correlations code

[u/mrichter/AliRoot.git] / FASTSIM / AliFastGlauber.cxx

diff --git a/FASTSIM/AliFastGlauber.cxx b/FASTSIM/AliFastGlauber.cxx
index 672df6ba83ef3ae65c467a8e2d6909c0a22755da..a792609dcccd2b3600bb482d12c3dd1fdcc67434 100644 (file)
--- a/FASTSIM/AliFastGlauber.cxx
+++ b/FASTSIM/AliFastGlauber.cxx
@@ -65,6 +65,7 @@ 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;    
@@ -187,6 +188,12 @@ void AliFastGlauber::Init(Int_t mode)
   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);
@@ -268,7 +275,7 @@ void AliFastGlauber::Init(Int_t mode)
   //
   // Hard collisions per event
   //
-  fgWSN = new TF1("WSN", WSN, 0., fgBMax, 1);
+  fgWSN = new TF1("WSN", WSN, 0.01, fgBMax, 1);
   fgWSN->SetNpx(100);
 
   //
@@ -284,7 +291,7 @@ void AliFastGlauber::Init(Int_t mode)
     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;
     }
@@ -314,6 +321,7 @@ void AliFastGlauber::Reset() const
   // 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;
@@ -338,7 +346,7 @@ 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;
+  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]");
@@ -350,13 +358,13 @@ void AliFastGlauber::DrawWSb() const
   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();
@@ -396,7 +404,7 @@ void AliFastGlauber::DrawParticipants() const
   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();
@@ -436,7 +444,7 @@ void AliFastGlauber::DrawGeo() const
   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();
@@ -460,7 +468,7 @@ void AliFastGlauber::DrawBinary() const
   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();
@@ -473,7 +481,7 @@ void AliFastGlauber::DrawN() const
   //
   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]");
@@ -484,9 +492,9 @@ void AliFastGlauber::DrawN() const
   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();
@@ -510,7 +518,7 @@ void AliFastGlauber::DrawKernel(Double_t b) const
   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();
@@ -524,17 +532,18 @@ 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,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();
@@ -616,7 +625,7 @@ void AliFastGlauber::DrawIntRadius(Double_t b) const
   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
@@ -631,7 +640,23 @@ Double_t AliFastGlauber::WSb(Double_t* x, Double_t* par)
   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
@@ -648,7 +673,7 @@ Double_t AliFastGlauber::WSbz(Double_t* x, Double_t* par)
   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
@@ -665,7 +690,7 @@ Double_t AliFastGlauber::WSz(Double_t* x, Double_t* par)
   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
@@ -677,7 +702,7 @@ Double_t AliFastGlauber::WSta(Double_t* x, Double_t* /*par*/)
   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)
@@ -690,7 +715,7 @@ Double_t AliFastGlauber::WStarfi(Double_t* x, Double_t* par)
   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 
@@ -734,7 +759,7 @@ Double_t AliFastGlauber::WStaa(Double_t* x, Double_t* par)
   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
@@ -763,7 +788,7 @@ Double_t AliFastGlauber::WKParticipants(Double_t* x, Double_t* par)
   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 
@@ -790,7 +815,7 @@ Double_t AliFastGlauber::WParticipants(Double_t* x, Double_t* par)
   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
@@ -803,7 +828,7 @@ Double_t AliFastGlauber::WSgeo(Double_t* x, Double_t* par)
   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
@@ -816,7 +841,7 @@ Double_t AliFastGlauber::WSbinary(Double_t* x, Double_t* par)
   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
@@ -826,7 +851,7 @@ Double_t AliFastGlauber::WSN(Double_t* x, Double_t* /*par*/)
   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 
@@ -845,7 +870,7 @@ Double_t AliFastGlauber::WEnergyDensity(Double_t* x, Double_t* par)
   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
@@ -864,7 +889,7 @@ Double_t AliFastGlauber::WAlmond(Double_t* x, Double_t* par)
   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 
@@ -889,7 +914,7 @@ Double_t AliFastGlauber::WIntRadius(Double_t* x, Double_t* par)
   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 
@@ -935,7 +960,7 @@ Double_t AliFastGlauber::WPathLength0(Double_t* x, Double_t* par)
   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 
@@ -1019,7 +1044,7 @@ Double_t AliFastGlauber::FractionOfHardCrossSection(Double_t b1, Double_t b2) co
   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 
@@ -1246,7 +1271,7 @@ void AliFastGlauber::StoreAlmonds() const
   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);
@@ -1455,7 +1480,7 @@ void AliFastGlauber::GetLengthsBackToBack(Double_t& ell1,Double_t& ell2,
   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]