Changes to the DEtaDPhi corr fun (making 0 in the center of the bin in the histogram...
[u/mrichter/AliRoot.git] / FASTSIM / AliFastGlauber.cxx
index 201207b..a792609 100644 (file)
@@ -16,7 +16,7 @@
 /* $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)
 
 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;    
@@ -82,28 +82,68 @@ TF1*    AliFastGlauber::fgWIntRadius     = NULL;
 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()
@@ -148,7 +188,13 @@ void AliFastGlauber::Init(Int_t mode)
   fgWSb->SetParameter(2, fWSw);
   fgWSb->SetParameter(3, fWSn);
 
-  fgWSbz = new TF2("WSbz", WSbz, 0, fgBMax, 4);
+  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);
@@ -229,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);
 
   //
@@ -245,9 +291,9 @@ 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);
-      new(&fgWAlmondFixedB[k]) TF2(*fgWAlmondCurrent);
+      fgWAlmondFixedB[k] = fgWAlmondCurrent;
     }
     delete ff;
   }
@@ -268,13 +314,14 @@ void AliFastGlauber::Init(Int_t mode)
   fgWPathLength->SetParameter(2, 0);     //Pathlength definition
 }
 
-void AliFastGlauber::Reset()
+void AliFastGlauber::Reset() const
 {
   //
   // Reset dynamic allocated formulas
   // 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;
@@ -299,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]");
@@ -311,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();
@@ -331,7 +378,7 @@ 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);
+  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}]");
@@ -357,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();
@@ -397,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();
@@ -421,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();
@@ -434,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]");
@@ -445,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();
@@ -471,7 +518,7 @@ void AliFastGlauber::DrawKernel(Double_t b) const
   l1a->SetFillStyle(0);
   l1a->SetBorderSize(0);
   Char_t label[100];
-  sprintf(label,"b = %.f fm",b);
+  snprintf(label, 100, "b = %.1f fm",b);
   l1a->AddEntry(fgWStarfi,label,"");
   l1a->Draw();
   c8->Update();
@@ -485,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 = %.f fm",b);
+  snprintf(label, 100, "b = %.1f fm",b);
   l1a->AddEntry(fgWAlmond,label,"");
   l1a->Draw();
   c9->Update();
@@ -577,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
@@ -592,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
@@ -609,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
@@ -626,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
@@ -638,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)
@@ -651,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 
@@ -681,8 +745,11 @@ Double_t AliFastGlauber::WStaa(Double_t* x, Double_t* par)
   //  MC Integration
   //
   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);
@@ -692,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
@@ -721,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 
@@ -748,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
@@ -761,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
@@ -774,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
@@ -784,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 
@@ -803,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
@@ -822,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 
@@ -847,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 
@@ -884,16 +951,16 @@ Double_t AliFastGlauber::WPathLength0(Double_t* x, Double_t* par)
   } // radial steps
 
   Double_t y=0.;
-  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);
+  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 
@@ -972,26 +1039,77 @@ Double_t AliFastGlauber::HardCrossSection(Double_t b1, Double_t b2) const
 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::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 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::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 hard collisions at b
   //
-  if(b==0) b=1e-4;
+  if(b<1.e-4) b=1e-4;
   return fgWSN->Eval(b);
 }
 
@@ -1000,7 +1118,7 @@ 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));
 }
 
@@ -1009,7 +1127,7 @@ 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));
 }
 
@@ -1018,10 +1136,24 @@ Double_t AliFastGlauber::GetNumberOfCollisions(Double_t  b) const
   //
   // 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);
 }
 
+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)
 {
   //
@@ -1139,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);
@@ -1190,7 +1322,7 @@ 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]; 
+  fgWAlmondCurrent = fgWAlmondFixedB[bin]; 
   return;
 }
 
@@ -1295,35 +1427,60 @@ Double_t AliFastGlauber::CalculateLength(Double_t b,Double_t x0,Double_t y0,Doub
   }
 }
 
-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 2 lengths back to back from random b, x0, y0, phi0 
+  // 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 kphi0plusPi = phi0+TMath::Pi();
+  phi = phi0;
   ell1 = CalculateLength(b,x0,y0,phi0);
   ell2 = CalculateLength(b,x0,y0,kphi0plusPi);
   return;
 }
 
-void AliFastGlauber::GetLengthsForPythia(Int_t n,Double_t* phi,Double_t* ell, Double_t b)
+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;
+}
+
+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] 
@@ -1356,7 +1513,7 @@ void AliFastGlauber::PlotBDistr(Int_t n)
   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
@@ -1383,7 +1540,7 @@ void AliFastGlauber::PlotLengthDistr(Int_t n,Bool_t save,Char_t *fname)
   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
@@ -1418,13 +1575,13 @@ void AliFastGlauber::PlotAlmonds() const
   gStyle->SetPalette(1,0);
   c->Divide(2,2);
   c->cd(1);
-  fgWAlmondFixedB[0].Draw("cont1");
+  fgWAlmondFixedB[0]->Draw("cont1");
   c->cd(2);
-  fgWAlmondFixedB[10].Draw("cont1");
+  fgWAlmondFixedB[10]->Draw("cont1");
   c->cd(3);
-  fgWAlmondFixedB[20].Draw("cont1");
+  fgWAlmondFixedB[20]->Draw("cont1");
   c->cd(4);
-  fgWAlmondFixedB[30].Draw("cont1");
+  fgWAlmondFixedB[30]->Draw("cont1");
   return;
 }
 
@@ -1436,8 +1593,8 @@ void AliFastGlauber::CalculateI0I1(Double_t& integral0,Double_t& integral1,
 {
   // 
   // 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)) 
@@ -1474,37 +1631,87 @@ void AliFastGlauber::CalculateI0I1(Double_t& integral0,Double_t& integral1,
   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);
+  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)
+{
+  //
+  // 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;
+}
+
 void AliFastGlauber::GetI0I1ForPythia(Int_t n,Double_t* phi,
                                      Double_t* integral0,Double_t* integral1,
                                      Double_t ellCut,Double_t b)
@@ -1521,20 +1728,43 @@ void AliFastGlauber::GetI0I1ForPythia(Int_t n,Double_t* phi,
   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);
@@ -1556,6 +1786,7 @@ void AliFastGlauber::PlotI0I1Distr(Int_t n,Double_t ellCut,
 
   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);
@@ -1564,6 +1795,10 @@ void AliFastGlauber::PlotI0I1Distr(Int_t n,Double_t ellCut,
   }
   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);
@@ -1577,9 +1812,13 @@ void AliFastGlauber::PlotI0I1Distr(Int_t n,Double_t ellCut,
   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();
@@ -1591,10 +1830,10 @@ void AliFastGlauber::PlotI0I1Distr(Int_t n,Double_t ellCut,
 }
 
 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); 
@@ -1629,3 +1868,18 @@ void AliFastGlauber::PlotI0I1B2BDistr(Int_t n,Double_t ellCut,
   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");
+}
+