Added support for I0I1 calculation need for k-Quenching scenario (A. Dainese).
authorloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 27 Mar 2004 17:13:00 +0000 (17:13 +0000)
committerloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 27 Mar 2004 17:13:00 +0000 (17:13 +0000)
Handling of AuAu collisions at Rhic
Reworked plots, comments and cleanup of code

FASTSIM/AliFastGlauber.cxx
FASTSIM/AliFastGlauber.h

index 6692039..b39b7cd 100644 (file)
 
 /* $Id$ */
 //
-// Utility class to make simple Glauber type calculations for collision geometries:
+// Utility class to make simple Glauber type calculations 
+//           for SYMMTRIC collision geometries (AA):
 // Impact parameter, production points, reaction plane dependence
+//
 // The SimulateTrigger method can be used for simple MB and hard-process
 // (binary scaling) trigger studies.
+// 
 // Some basic quantities can be visualized directly.
-// The default set-up for PbPb collisions can be read from a file calling Init(1).
 //
+// The default set-up for PbPb or AUAu collisions can be read from a file 
+// calling Init(1) or Init(2) if you want to read Almonds too.
+//
+// ***** If you change settings dont forget to call init afterwards, *****
+// ***** in order to update the formulas with the new parameters.    *****
 // 
 // Author: andreas.morsch@cern.ch
+//=================== Added by A. Dainese 11/02/04 ===========================
+// Calculate path length for a parton with production point (x0,y0)
+// and propagation direction (ux=cos(phi0),uy=sin(phi0)) 
+// in a collision with impact parameter b and functions that make use
+// of it.
+//=================== Added by A. Dainese 05/03/04 ===========================
+// Calculation of line integrals I0 and I1
+//  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)
+// mostly for use in the Quenching class
+//=================== Added by C. Loizdes 27/03/04 ===========================
+// Handling of AuAu collisions
+// More get/set functions
+// Comments, units and clearing of code
+//
 
 // from AliRoot
 #include "AliFastGlauber.h"
 #include <TH2F.h>
 #include <TF1.h>
 #include <TF2.h>
+#include <TLegend.h>
 #include <TCanvas.h>
 #include <TRandom.h>
 #include <TFile.h>
+#include <TROOT.h>
+#include <Riostream.h>
 
 ClassImp(AliFastGlauber)
 
+Float_t AliFastGlauber::fgBMax           = 0.;
 TF1*    AliFastGlauber::fgWSb            = NULL;     
 TF2*    AliFastGlauber::fgWSbz           = NULL;    
 TF1*    AliFastGlauber::fgWSz            = NULL;     
 TF1*    AliFastGlauber::fgWSta           = NULL;    
 TF2*    AliFastGlauber::fgWStarfi        = NULL; 
 TF2*    AliFastGlauber::fgWAlmond        = NULL; 
-//TF2*    AliFastGlauber::fWAlmondCurrent = NULL; 
 TF1*    AliFastGlauber::fgWStaa          = NULL;   
 TF1*    AliFastGlauber::fgWSgeo          = NULL;   
 TF1*    AliFastGlauber::fgWSbinary       = NULL;   
@@ -56,831 +81,1111 @@ TF1*    AliFastGlauber::fgWEnergyDensity = NULL;
 TF1*    AliFastGlauber::fgWIntRadius     = NULL;   
 TF2*    AliFastGlauber::fgWKParticipants = NULL; 
 TF1*    AliFastGlauber::fgWParticipants  = NULL; 
-Float_t AliFastGlauber::fgBMax           = 0.;
+TF2*    AliFastGlauber::fgWAlmondCurrent = NULL;    
+TF2     AliFastGlauber::fgWAlmondFixedB[40]; 
+const Int_t AliFastGlauber::fgkMCInts = 100000;
+Int_t AliFastGlauber::fgCounter = 0;       
 
-AliFastGlauber::AliFastGlauber()
+AliFastGlauber::AliFastGlauber() : fName()
 {
-//  Default Constructor
-//
-//  Defaults for Pb
-    SetWoodSaxonParameters(6.624, 0.549, 0.00, 7.69e-4);
-    SetHardCrossSection();
-    SetMaxImpact();
-    SetLengthDefinition();
+  //  Default Constructor
+  //
+  fgCounter++;
+  if(fgCounter>1)
+    Error("AliFastGlauber","More than more instance (%d) is not supported, check your code!",fgCounter);
+
+  //  Defaults for Pb
+  SetMaxImpact();
+  SetLengthDefinition();
+  SetPbPbLHC();
+}
+
+AliFastGlauber::~AliFastGlauber()
+{
+  fgCounter--;
+  //  if(fgCounter==0) Reset();
+}
+
+void AliFastGlauber::SetAuAuRhic()
+{
+  //Set all parameters for RHIC
+  SetWoodSaxonParametersAu();
+  SetHardCrossSection();
+  SetNNCrossSection(40);
+  SetNucleus(197);
+  SetFileName("$(ALICE_ROOT)/FASTSIM/data/glauberAuAu.root");
+}
+
+void AliFastGlauber::SetPbPbLHC()
+{
+  //Set all parameters for LHC
+  SetWoodSaxonParametersPb();
+  SetHardCrossSection();
+  SetNNCrossSection();
+  SetNucleus();
+  SetFileName();
 }
 
 void AliFastGlauber::Init(Int_t mode)
 {
-// Initialisation
-//    mode = 0; all functions are calculated 
-//    mode = 1; overlap function is read from file (for Pb-Pb only)
-//    mode = 2; interaction almond functions are read from file 
-//              (for Pb-Pb only); USE THIS FOR PATH LENGTH CALC.!
-//
+  // Initialisation
+  //    mode = 0; all functions are calculated 
+  //    mode = 1; overlap function is read from file (for Pb-Pb only)
+  //    mode = 2; interaction almond functions are read from file 
+  //              USE THIS FOR PATH LENGTH CALC.!
+  //
 
-//
-//  Wood-Saxon
-//
-    fgWSb = new TF1("WSb", WSb, 0, fgBMax, 4);
-    fgWSb->SetParameter(0, fWSr0);
-    fgWSb->SetParameter(1, fWSd);
-    fgWSb->SetParameter(2, fWSw);
-    fgWSb->SetParameter(3, fWSn);
-
-    fgWSbz = new TF2("WSbz", WSbz, 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);
+  // 
+  Reset(); 
+  //
 
-//
-//  Thickness
-//
-    fgWSta = new TF1("WSta", WSta, 0., fgBMax, 0);
+  //
+  //  Wood-Saxon
+  //
+  fgWSb = new TF1("WSb", WSb, 0, fgBMax, 4);
+  fgWSb->SetParameter(0, fWSr0);
+  fgWSb->SetParameter(1, fWSd);
+  fgWSb->SetParameter(2, fWSw);
+  fgWSb->SetParameter(3, fWSn);
+
+  fgWSbz = new TF2("WSbz", WSbz, 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);
+
+  //
+  //  Thickness
+  //
+  fgWSta = new TF1("WSta", WSta, 0., fgBMax, 0);
     
-//
-//  Overlap Kernel
-//
-    fgWStarfi = new TF2("WStarfi", WStarfi, 0., fgBMax, 0., TMath::Pi(), 1);
-    fgWStarfi->SetParameter(0, 0.);     
-    fgWStarfi->SetNpx(200);     
-    fgWStarfi->SetNpy(20);    
-//
-//  Participants Kernel
-//
-    fgWKParticipants = new TF2("WKParticipants", WKParticipants, 0., fgBMax, 0., TMath::Pi(), 1);
-    fgWKParticipants->SetParameter(0, 0.);     
-    fgWKParticipants->SetNpx(200);     
-    fgWKParticipants->SetNpy(20);      
-//
-//  Almond shaped interaction region
-//
-    fgWAlmond = new TF2("WAlmond", WAlmond, -fgBMax, fgBMax, -fgBMax, fgBMax, 1);
-    fgWAlmond->SetParameter(0, 0.);     
-    fgWAlmond->SetNpx(200);     
-    fgWAlmond->SetNpy(200);  
-  
-    if(mode==2) {
-      printf(" Reading interaction almonds from file\n");
-      Char_t almondName[100];
-      TFile* ff = new TFile("$(ALICE_ROOT)/FASTSIM/data/glauberPbPb.root");
-      for(Int_t k=0; k<40; k++) {
-       sprintf(almondName,"WAlmondFixedB%d",k);
-       fWAlmondCurrent = (TF2*)ff->Get(almondName);
-       new(&fWAlmondFixedB[k]) TF2(*fWAlmondCurrent);
-      }
-    }
-//
-//  Path Length as a function of Phi
-//    
-    fgWPathLength0 = new TF1("WPathLength0", WPathLength0, -TMath::Pi(), TMath::Pi(), 2);
-    fgWPathLength0->SetParameter(0, 0.);
-//  Pathlength definition     
-    fgWPathLength0->SetParameter(1, 0.);     
+  //
+  //  Overlap Kernel
+  //
+  fgWStarfi = new TF2("WStarfi", WStarfi, 0., fgBMax, 0., TMath::Pi(), 1);
+  fgWStarfi->SetParameter(0, 0.);     
+  fgWStarfi->SetNpx(200);     
+  fgWStarfi->SetNpy(20);    
 
-    fgWPathLength = new TF1("WPathLength", WPathLength, -TMath::Pi(), TMath::Pi(), 3);
-//  Impact Parameter
-    fgWPathLength->SetParameter(0, 0.);    
-//  Number of interactions used for average
-    fgWPathLength->SetParameter(1, 1000.);    
-//  Pathlength definition
-    fgWPathLength->SetParameter(2, 0);    
+  //
+  //  Participants Kernel
+  //
+  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);      
 
-    fgWIntRadius = new TF1("WIntRadius", WIntRadius, 0., fgBMax, 1);
-    fgWIntRadius->SetParameter(0, 0.);    
+  //
+  //  Overlap and Participants
+  //
+  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;
+  }
 
+  //
+  //  Energy Density
+  //
+  fgWEnergyDensity = new TF1("WEnergyDensity", WEnergyDensity, 0., 2. * fWSr0, 1);
+  fgWEnergyDensity->SetParameter(0, fWSr0 + 1.);
+    
+  //
+  //  Geometrical Cross-Section
+  //
+  fgWSgeo = new TF1("WSgeo", WSgeo, 0., fgBMax, 1);
+  fgWSgeo->SetParameter(0,fSigmaNN); //mbarn
+  fgWSgeo->SetNpx(100);
 
-//
-//  Overlap
-//
-    if (! mode) {
-       fgWStaa = new TF1("WStaa", WStaa, 0., fgBMax, 0);
-       fgWStaa->SetNpx(100);
-       fgWParticipants = new TF1("WParticipants", WParticipants, 0., fgBMax, 0);
-       fgWParticipants->SetNpx(100);
+  //
+  //  Hard cross section (binary collisions)
+  //
+  fgWSbinary = new TF1("WSbinary", WSbinary, 0., fgBMax, 1);
+  fgWSbinary->SetParameter(0, fSigmaHard); //mbarn
+  fgWSbinary->SetNpx(100);
 
-    } else {
-       TFile* f = new TFile("$(ALICE_ROOT)/FASTSIM/data/glauberPbPb.root");
-       fgWStaa = (TF1*) f->Get("WStaa");
-       fgWParticipants = (TF1*) f->Get("WParticipants");
+  //
+  // Hard collisions per event
+  //
+  fgWSN = new TF1("WSN", WSN, 0., fgBMax, 1);
+  fgWSN->SetNpx(100);
+
+  //
+  //  Almond shaped interaction region
+  //
+  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++) {
+      sprintf(almondName,"WAlmondFixedB%d",k);
+      fgWAlmondCurrent = (TF2*)ff->Get(almondName);
+      new(&fgWAlmondFixedB[k]) TF2(*fgWAlmondCurrent);
     }
-    
-//
-//  Participants
-//
-    
-//
-    fgWEnergyDensity = new TF1("WEnergyDensity", WEnergyDensity, 0., 2. * fWSr0, 1);
-    fgWEnergyDensity->SetParameter(0, fWSr0 + 1.);
-    
-//
-//  Geometrical Cross-Section
-//
-    fgWSgeo = new TF1("WSgeo", WSgeo, 0., fgBMax, 0);
-    fgWSgeo->SetNpx(100);
-//
-//  Hard cross section (~ binary collisions)
-//
-    fgWSbinary = new TF1("WSbinary", WSbinary, 0., fgBMax, 1);
-    fgWSbinary->SetParameter(0, fSigmaHard); // mb
-    fgWSbinary->SetNpx(100);
-//
-// Hard collisions per event
-//
-    fgWSN = new TF1("WSN", WSN, 0., fgBMax, 1);
-    fgWSN->SetNpx(100);
+    delete ff;
+  }
+
+  fgWIntRadius = new TF1("WIntRadius", WIntRadius, 0., fgBMax, 1);
+  fgWIntRadius->SetParameter(0, 0.);    
+
+  //
+  //  Path Length as a function of Phi
+  //    
+  fgWPathLength0 = new TF1("WPathLength0", WPathLength0, -TMath::Pi(), TMath::Pi(), 2);
+  fgWPathLength0->SetParameter(0, 0.);
+  fgWPathLength0->SetParameter(1, 0.); //Pathlength definition     
+
+  fgWPathLength = new TF1("WPathLength", WPathLength, -TMath::Pi(), TMath::Pi(), 3);
+  fgWPathLength->SetParameter(0, 0.);    //Impact Parameter
+  fgWPathLength->SetParameter(1, 1000.); //Number of interactions used for average
+  fgWPathLength->SetParameter(2, 0);     //Pathlength definition
+}
+
+void AliFastGlauber::Reset()
+{
+  //
+  // Reset dynamic allocated formulas
+  // in case init is called twice
+
+  if(fgWSb)            delete fgWSb;     
+  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()
 {
-//
-//  Draw Wood-Saxon Nuclear Density Function
-//
-    TCanvas *c1 = new TCanvas("c1","Wood Saxon",400,10,600,700);
-    c1->cd();
-    fgWSb->Draw();
+  //
+  //  Draw Wood-Saxon Nuclear Density Function
+  //
+  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];
+  sprintf(label,"r_{0} = %.2f fm",fWSr0);
+  l1a->AddEntry(fgWSb,label,"");
+  sprintf(label,"d = %.2f fm",fWSd);
+  l1a->AddEntry(fgWSb,label,"");
+  sprintf(label,"n = %.2e fm^{-3}",fWSn);
+  l1a->AddEntry(fgWSb,label,"");
+  sprintf(label,"#omega = %.2f",fWSw);
+  l1a->AddEntry(fgWSb,label,"");
+  l1a->Draw();
+  c1->Update();
 }
 
 void AliFastGlauber::DrawOverlap()
 {
-//
-//  Draw Overlap Function
-//
-    TCanvas *c2 = new TCanvas("c2","Overlap",400,10,600,700);
-    c2->cd();
-    fgWStaa->Draw();
+  //
+  //  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);
+  h2f->SetStats(0);
+  h2f->GetXaxis()->SetTitle("b [fm]");
+  h2f->GetYaxis()->SetTitle("T_{AB} [mbarn^{-1}]");
+  h2f->Draw(); 
+  fgWStaa->Draw("same");
 }
 
 void AliFastGlauber::DrawParticipants()
 {
-//
-//  Draw Number of Participants
-//
-    TCanvas *c2 = new TCanvas("c2","Participants",400,10,600,700);
-    c2->cd();
-    fgWParticipants->Draw();
+  //
+  //  Draw Number of Participants Npart
+  //
+  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];
+  sprintf(label,"#sigma^{inel.}_{NN} = %.1f mbarn",fSigmaNN);
+  l1a->AddEntry(fgWParticipants,label,"");
+  l1a->Draw();
+  c3->Update();
 }
 
 void AliFastGlauber::DrawThickness()
 {
-//
-//  Draw Thickness Function
-//
-    TCanvas *c3 = new TCanvas("c3","Thickness",400,10,600,700);
-    c3->cd();
-    fgWSta->Draw();
+  //
+  //  Draw Thickness Function
+  //
+  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()
 {
-//
-//  Draw Geometrical Cross-Section
-//
-    TCanvas *c3 = new TCanvas("c3","Geometrical Cross-Section",400,10,600,700);
-    c3->cd();
-    fgWSgeo->Draw();
+  //
+  //  Draw Geometrical Cross-Section
+  //
+  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];
+  sprintf(label,"#sigma_{NN}^{inel.} = %.1f mbarn",fSigmaNN);
+  l1a->AddEntry(fgWSgeo,label,"");
+  l1a->Draw();
+  c5->Update();
 }
 
 void AliFastGlauber::DrawBinary()
 {
-//
-//  Draw Binary Cross-Section
-//
-    TCanvas *c4 = new TCanvas("c4","Binary Cross-Section",400,10,600,700);
-    c4->cd();
-    fgWSbinary->Draw();
+  //
+  //  Draw Binary Cross-Section
+  //
+  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];
+  sprintf(label,"#sigma_{NN}^{hard} = %.1f mbarn",fSigmaHard);
+  l1a->AddEntry(fgWSb,label,"");
+  l1a->Draw();
+  c6->Update();
 }
 
 void AliFastGlauber::DrawN()
 {
-//
-//  Draw Binaries per event
-//
-    TCanvas *c5 = new TCanvas("c5","Binaries per event",400,10,600,700);
-    c5->cd();
-    fgWSN->Draw();
+  //
+  //  Draw Binaries per event (Ncoll)
+  //
+  TCanvas *c7 = new TCanvas("c7","Binaries per event",400,10,600,700);
+  c7->cd();
+  Double_t max=fgWSN->GetMaximum(0,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];
+  sprintf(label,"#sigma^{hard}_{NN} = %.1f mbarn",fSigmaHard);
+  l1a->AddEntry(fgWSN,label,"");
+  sprintf(label,"#sigma^{inel.}_{NN} = %.1f mbarn",fSigmaNN);
+  l1a->AddEntry(fgWSN,label,"");
+  l1a->Draw();
+  c7->Update();
 }
 
 void AliFastGlauber::DrawKernel(Double_t b)
 {
-//
-//  Draw Kernel
-//
-    TCanvas *c6 = new TCanvas("c6","Kernel",400,10,600,700);
-    c6->cd();
-    fgWStarfi->SetParameter(0, b);
-    fgWStarfi->Draw();
+  //
+  //  Draw Kernel
+  //
+  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];
+  sprintf(label,"b = %.f fm",b);
+  l1a->AddEntry(fgWStarfi,label,"");
+  l1a->Draw();
+  c8->Update();
 }
 
 void AliFastGlauber::DrawAlmond(Double_t b)
 {
-//
-//  Draw Interaction Almond
-//
-    TCanvas *c7 = new TCanvas("c7","Almond",400,10,600,700);
-    c7->cd();
-    fgWAlmond->SetParameter(0, b);
-    fgWAlmond->Draw();
+  //
+  //  Draw Interaction Almond
+  //
+  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);
+  h2f->SetStats(0);
+  h2f->GetXaxis()->SetTitle("x [fm]");
+  h2f->GetYaxis()->SetTitle("y [fm]");
+  h2f->Draw(); 
+  fgWAlmond->Draw("same");
+  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);
+  l1a->AddEntry(fgWAlmond,label,"");
+  l1a->Draw();
+  c9->Update();
 }
 
-void AliFastGlauber::DrawPathLength0(Double_t b, Int_t iopt)
+void AliFastGlauber::DrawEnergyDensity()
 {
-//
-//  Draw Path Length
-//
-    TCanvas *c8 = new TCanvas("c8","Path Length",400,10,600,700);
-    c8->cd();
-    fgWPathLength0->SetParameter(0, b);
-    fgWPathLength0->SetParameter(1, Double_t(iopt));
-    fgWPathLength0->SetMinimum(0.); 
-    fgWPathLength0->SetMaximum(10.); 
-    fgWPathLength0->Draw();
+  //
+  //  Draw energy density
+  //
+  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::DrawPathLength(Double_t b , Int_t ni, Int_t iopt)
+void AliFastGlauber::DrawPathLength0(Double_t b, Int_t iopt)
 {
-//
-//  Draw Path Length
-//
-    TCanvas *c9 = new TCanvas("c9","Path Length",400,10,600,700);
-    c9->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.); 
-    fgWPathLength->Draw();
+  //
+  //  Draw Path Length
+  //
+  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::DrawIntRadius(Double_t b)
+void AliFastGlauber::DrawPathLength(Double_t b , Int_t ni, Int_t iopt)
 {
-//
-//  Draw Interaction Radius
-//
-    TCanvas *c10 = new TCanvas("c10","Interaction Radius",400,10,600,700);
-    c10->cd();
-    fgWIntRadius->SetParameter(0, b);
-    fgWIntRadius->SetMinimum(0.);
-    fgWIntRadius->Draw();
+  //
+  //  Draw Path Length
+  //
+  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::DrawEnergyDensity()
+void AliFastGlauber::DrawIntRadius(Double_t b)
 {
-//
-//  Draw energy density
-//
-    TCanvas *c11 = new TCanvas("c11","Energy Density",400, 10, 600, 700);
-    c11->cd();
-    fgWEnergyDensity->SetMinimum(0.);
-    fgWEnergyDensity->Draw();
+  //
+  //  Draw Interaction Radius
+  //
+  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(Double_t* x, Double_t* par)
 {
-//
-//  Woods-Saxon Parameterisation
-//  as a function of radius
-//
-    Double_t xx  = x[0];
-    Double_t r0  = par[0];
-    Double_t d   = par[1];
-    Double_t w   = par[2];
-    Double_t n   = par[3];
-    
-    Double_t y  = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
-
-    return y;
+  //
+  //  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));
+  return y; //fm^-3
 }
 
 Double_t AliFastGlauber::WSbz(Double_t* x, Double_t* par)
 {
-//
-//  Wood Saxon Parameterisation
-//  as a function of z and  b
-//
-    Double_t bb  = x[0];
-    Double_t zz  = x[1];
-    Double_t r0  = par[0];
-    Double_t d   = par[1];
-    Double_t w   = par[2];
-    Double_t n   = par[3];
-    Double_t xx  = TMath::Sqrt(bb*bb+zz*zz);
-    Double_t y  = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
-
-    return y;
+  //
+  //  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));
+  return y; //fm^-3
 }
 
 Double_t AliFastGlauber::WSz(Double_t* x, Double_t* par)
 {
-//
-//  Wood Saxon Parameterisation
-//  as a function of z for fixed b
-//
-    Double_t bb  = par[4];
-    Double_t zz  = x[0];
-    Double_t r0  = par[0];
-    Double_t d   = par[1];
-    Double_t w   = par[2];
-    Double_t n   = par[3];
-    Double_t xx  = TMath::Sqrt(bb*bb+zz*zz);
-    Double_t y  = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
-
-    return y;
+  //
+  //  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));
+  return y; //fm^-3
 }
 
 Double_t AliFastGlauber::WSta(Double_t* x, Double_t* /*par*/)
 {
-//
-//  Thickness function 
-//
-    Double_t b  = x[0];
-    fgWSz->SetParameter(4, b);
-    Double_t y  = 2. * fgWSz->Integral(0., fgBMax);
-    return y;
+  //
+  //  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);
+  return y; //fm^-2
 }
 
-
-
 Double_t AliFastGlauber::WStarfi(Double_t* x, Double_t* par)
 {
-//
-//  Kernel for overlap function
-//
-    Double_t b    = par[0];
-    Double_t r1   = x[0];
-    Double_t phi  = x[1];
-    Double_t r2   = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi)); 
-    Double_t y    = r1 * fgWSta->Eval(r1) * fgWSta->Eval(r2);
-    return y;
+  //
+  //  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);
+  return y; //fm^-3
+}
+
+Double_t AliFastGlauber::WStaa(Double_t* x, Double_t* par)
+{
+  //
+  //  Overlap function 
+  //  T_{AB}=Int d2s T_A(s)*T_B(s-b)
+  //  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);
+
+  // root integration seems to fail
+  /* 
+     Double_t al[2];
+     Double_t bl[2];
+     al[0] = 1e-6;
+     al[1] = fgBMax;
+     bl[0] = 0.;
+     bl[1] = TMath::Pi();
+     Double_t err;
+     
+     Double_t y =  2. *  208. * 208. * fgWStarfi->IntegralMultiple(2, al, bl, 0.001, err);
+     printf("WStaa: %.5e %.5e %.5e\n", b, y, err);
+  */
+
+  //
+  //  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);
+    }
+  y *= 2. * TMath::Pi() * fgBMax / fgkMCInts; //fm^-2
+  y *= A * A * 0.1; //mbarn^-1
+  return y;
 }
 
 Double_t AliFastGlauber::WKParticipants(Double_t* x, Double_t* par)
 {
-//
-//  Kernel for number of participants
-//
-    Double_t b    = par[0];
-    Double_t r1   = x[0];
-    Double_t phi  = x[1];
-    Double_t r2   = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi)); 
-    Double_t xsi  = fgWSta->Eval(r2) * 5.6;
-    Double_t a = 208;
-    Double_t sum = a * xsi;
-    Double_t y   = sum;
-    for (Int_t i = 1; i <= 208; i++)
+  //
+  //  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;
+  /*
+    Double_t y=(1-TMath::Power((1-xsi),A))
+   */
+  Double_t a = A;
+  Double_t sum = a * xsi;
+  Double_t y   = sum;
+  for (Int_t i = 1; i <= A; i++)
     {
-       a--;
-       sum *= (-xsi) * a / Float_t(i+1);
-       y  += sum;
+      a--;
+      sum *= (-xsi) * a / Float_t(i+1);
+      y  += sum;
     }
-    
-    y    = r1 * fgWSta->Eval(r1) * y;
-    return y;
+  y    = r1 * fgWSta->Eval(r1) * y;
+  return y; //fm^-1
 }
 
+Double_t AliFastGlauber::WParticipants(Double_t* x, 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);
+
+  //
+  //  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);
+    }
+  y *= 2. *  A * 2. * TMath::Pi() * fgBMax / fgkMCInts;
+  return y; //no units
+}
+
+Double_t AliFastGlauber::WSgeo(Double_t* x, 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)); 
+  return y; //fm
+}
+
+Double_t AliFastGlauber::WSbinary(Double_t* x, Double_t* par)
+{
+  //
+  //  Number of binary 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; 
+  return y; //fm
+}
+
+Double_t AliFastGlauber::WSN(Double_t* x, 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);
+  return y; //no units
+}
+
+Double_t AliFastGlauber::WEnergyDensity(Double_t* x, 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];
+  //
+  //  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;
+  return y; //fm^-4
+}
 
 Double_t AliFastGlauber::WAlmond(Double_t* x, Double_t* par)
 {
-//
-//  Almond shaped interaction region
-//
-    Double_t b    = par[0];
-    Double_t xx   = x[0] + b/2.;
-    Double_t yy   = x[1];
-    Double_t r1   = TMath::Sqrt(xx * xx + yy * yy);
-    Double_t phi  = TMath::ATan2(yy,xx);
-    
-    Double_t r2   = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi)); 
-//
-//  Interaction probability calculated as product of thicknesses
-//
-    Double_t y    = fgWSta->Eval(r1) * fgWSta->Eval(r2);
-    return y;
+  //
+  //  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)); 
+  //
+  //  Interaction probability calculated as product of thicknesses
+  //
+  const Double_t y    = fgWSta->Eval(r1) * fgWSta->Eval(r2);
+  return y; //fm^-4
 }
 
 Double_t AliFastGlauber::WIntRadius(Double_t* x, Double_t* par)
 {
-//
-//  Average radius at which interaction takes place
-//
-//  Radius in the Almond
-    Double_t r    = x[0];
-//  Impact parameter
-    Double_t b    = par[0];
-    fgWAlmond->SetParameter(0, b);
-//  Steps in phi
-    Double_t dphi = 2. * TMath::Pi() / 100.;
-//  Average over phi    
-    Double_t phi  = 0.;
-    Double_t y    = 0.;
-
-    for (Int_t i = 0; i < 100; i++) {
-       Double_t xx = r * TMath::Cos(phi);
-       Double_t yy = r * TMath::Sin(phi);
-       y   += fgWAlmond->Eval(xx,yy);
-       phi += dphi;
-    } // phi loop
-// Result multiplied by Jacobian (2 pi r)     
-    return (2. * TMath::Pi() * y * r / 100.);
+  //
+  //  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);
+  //  Average over phi in small steps   
+  const Double_t dphi = 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;
+  } // phi loop
+  // Result multiplied by Jacobian (2 pi r)     
+  y *= 2. * TMath::Pi() * r / 100.;
+  return y; //fm^-3
 }
 
 Double_t AliFastGlauber::WPathLength0(Double_t* x, Double_t* par)
 {
-//
-//  Path Length as a function of phi for interaction point fixed at (0,0)
-//
-//
-//  Steps in r 
-    const Int_t    kNp  = 100;
-    const Double_t kDr  = fgBMax/Double_t(kNp);
-//  Impact parameter    
-    Double_t b      = par[0];
-//  Path Length definition
-    Int_t    iopt   = Int_t(par[1]);
-    
-//  Phi direction in Almond
-    Double_t phi0   = x[0];
-    Double_t r  = 0.;
+  //
+  //  Path Length as a function of phi 
+  //  for interaction point fixed at (0,0)
+  //  as a function of phi-direction
+  //
+  //  Phi direction in Almond
+  const Double_t phi0   = x[0];
+  const Double_t b      = par[0];
+  //  Path Length definition
+  const Int_t    iopt   = Int_t(par[1]);
+
+  //  Step along radial direction phi   
+  const Int_t    kNp  = 100; // Steps in r 
+  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++) {
+    //
+    //  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);
+    // 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);
+
+    rw += y * r;
+    w  += y;
+    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);
+  }
+  return y; //fm
+}
+
+Double_t AliFastGlauber::WPathLength(Double_t* x, 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 
+
+  //
+  //  r-steps
+  // 
+  const Int_t    kNp   = 100;
+  const Double_t kDr  = fgBMax/Double_t(kNp);
+  Double_t l = 0.;  //  Path length 
+  for (Int_t in = 0; in < kNpi; in ++) {
     Double_t rw = 0.;
     Double_t w  = 0.;
-//  Step along radial direction phi   
-    for (Int_t i = 0; i < kNp; i++) {
-//
-//  Transform into target frame
-//
-       Double_t xx   = r * TMath::Cos(phi0) + b / 2.;
-       Double_t yy   = r * TMath::Sin(phi0);
-       Double_t phi  = TMath::ATan2(yy, xx);
+    // Interaction point
+    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;
        
-       Double_t r1   = TMath::Sqrt(xx * xx + yy * yy);
-// Radius in projectile frame
-       Double_t r2   = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi)); 
-       Double_t y    = fgWSta->Eval(r1) * fgWSta->Eval(r2);
-
-       rw += y * r;
-       w  += y;
-       r  += kDr;
-    } // radial steps
-//
-//  My length definition (is exact for hard disk)
+    // Radial steps
+    Double_t r  = 0.;
+    for (Int_t i = 0; (i < nps ); 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);
+      // 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);
+           
+      rw += y * r;
+      w  += y;
+      r  += kDr;
+    } // steps
+    // Average over interactions
     if (!iopt) {
-       return (2. * rw / w);
+      if(w) l += (2. * rw / w);
     } else {
-       return TMath::Sqrt(2. * rw * kDr / fgWSta->Eval(0.01) / fgWSta->Eval(0.01));
+      const Double_t norm=fgWSta->Eval(0.01);
+      if(norm) l+= 2. * rw * kDr / norm / norm;
     }
+  } // interactions
+  Double_t ret=0;
+  if (!iopt) 
+    ret= l / kNpi;
+  else 
+    ret=TMath::Sqrt( l / kNpi);
+  return ret; //fm
 }
 
-Double_t AliFastGlauber::WPathLength(Double_t* x, Double_t* par)
+Double_t AliFastGlauber::CrossSection(Double_t b1, Double_t b2)
 {
-//
-//  Path Length as a function of phi 
-//  Interaction point from random distribution
-//
-//
-//  r-steps
-// 
-    const Int_t    kNp   = 100;
-    const Double_t kDr  = fgBMax/Double_t(kNp);
-//  Number of interactions
-    const Int_t    kNpi  = Int_t (par[1]);
+  //
+  // Return the geometrical cross-section integrated from b1 to b2 
+  //
+  return fgWSgeo->Integral(b1, b2)*10.; //mbarn
+}
 
-//
-//  Impact parameter    
-    Double_t b      = par[0];
-//  Path Length definition 
-    Int_t    iopt   = Int_t(par[2]);
-//  Phi direction
-    Double_t phi0   = x[0];
-
-    printf("phi0 %f \n", phi0);
-    
-//  Path length 
-    Double_t l = 0.;
-    
-    for (Int_t in = 0; in < kNpi; in ++) {
-       Double_t rw = 0.;
-       Double_t w  = 0.;
-       
-       // Interaction point
-       Double_t x0, y0;
-       fgWAlmond->GetRandom2(x0, y0);
-// Initial radius
-       Double_t r0  = TMath::Sqrt(x0 * x0 + y0 * y0);
-       Int_t    nps = Int_t ((fgBMax - r0)/kDr) - 1;
-       
-       Double_t r  = 0.;
-// Radial steps
-       for (Int_t i = 0; (i < nps ); i++) {
-           
-// Transform into target frame
-           Double_t xx   = x0 + r * TMath::Cos(phi0) + b / 2.;
-           Double_t yy   = y0 + r * TMath::Sin(phi0);
-           Double_t phi  = TMath::ATan2(yy, xx);
-           Double_t r1   = TMath::Sqrt(xx * xx + yy * yy);
-// Radius in projectile frame
-           Double_t r2   = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi)); 
-           Double_t y    = fgWSta->Eval(r1) * fgWSta->Eval(r2);
-           
-           rw += y * r;
-           w  += y;
-           r  += kDr;
-       } // steps
-// Average over interactions
-       if (!iopt) {
-           l += (2. * rw / w);
-       } else {
-           l+= 2. * rw * kDr / fgWSta->Eval(0.01) / fgWSta->Eval(0.01);
-       }
-    } // interactions
-    if (!iopt) 
-       return (l / Double_t(kNpi));
-    else 
-       return (TMath::Sqrt(l / Double_t(kNpi)));
-}
-
-Double_t AliFastGlauber::WStaa(Double_t* x, Double_t* /*par*/)
+Double_t AliFastGlauber::HardCrossSection(Double_t b1, Double_t b2)
 {
-//
-//  Overlap function
-//
-    Double_t b    = x[0];
-    fgWStarfi->SetParameter(0, b);
-/*
-    Double_t al[2];
-    Double_t bl[2];
-    al[0] = 0.;
-    al[1] = 0.;
-    bl[0] = 6.6;
-    bl[1] = TMath::Pi();
-    Double_t err;
-    
-    Double_t y =  2. * fgWStarfi->IntegralMultiple(2, al, bl, 0.001, err);
-    printf("WStaa: %f %f %f\n", b, y, err);
-*/
-//
-//  MC Integration
-//
-    Double_t y = 0;
-    for (Int_t i = 0; i < 100000; i++)
-    {
-       Double_t phi = TMath::Pi() * gRandom->Rndm();
-       Double_t b1  = fgBMax       * gRandom->Rndm();  
-       y += fgWStarfi->Eval(b1, phi);
-    }
-    y *= 2. * 0.1 *  208. * 208. * TMath::Pi() * fgBMax / 100000.;
-    return y;
+  //
+  // Return the hard cross-section integrated from b1 to b2 
+  //
+  return fgWSbinary->Integral(b1, b2)*10.; //mbarn
 }
 
-Double_t AliFastGlauber::WParticipants(Double_t* x, Double_t* /*par*/)
+Double_t AliFastGlauber::FractionOfHardCrossSection(Double_t b1, Double_t b2)
 {
-//
-//  Overlap function
-//
-    Double_t b    = x[0];
-    fgWKParticipants->SetParameter(0, b);
-//
-//  MC Integration
-//
-    Double_t y = 0;
-    for (Int_t i = 0; i < 100000; i++)
-    {
-       Double_t phi = TMath::Pi() * gRandom->Rndm();
-       Double_t b1  = fgBMax      * gRandom->Rndm();   
-       y += fgWKParticipants->Eval(b1, phi);
-    }
-    y *= 4. *  208. * TMath::Pi() * fgBMax / 100000.;
-    return y;
+  //
+  // Return raction of hard cross-section integrated from b1 to b2 
+  //
+  return fgWSbinary->Integral(b1, b2)/fgWSbinary->Integral(0., 100.);
 }
 
-
-Double_t AliFastGlauber::WSgeo(Double_t* x, Double_t* /*par*/)
+Double_t AliFastGlauber::Binaries(Double_t b)
 {
-//
-//  Geometrical Cross-Section
-//
-    Double_t b    = x[0];
-    Double_t taa  = fgWStaa->Eval(b);
-    const Double_t kSigma = 55.6; // mbarn
-    
-    Double_t y    = 2. * TMath::Pi() * b * (1. - TMath::Exp(- kSigma * taa)); // fm
-    return y;
+  //
+  // Return number of binary collisions normalized to 1 at b=0
+  //
+  if(b==0) b=1e-4;
+  return fgWSN->Eval(b)/fgWSN->Eval(1e-4);
 }
 
-
-Double_t AliFastGlauber::WSbinary(Double_t* x, Double_t* par)
+Double_t AliFastGlauber::GetNumberofBinaries(Double_t b)
 {
-//
-//  Number of binary collisions
-//
-    Double_t b     = x[0];
-    Double_t sigma = par[0];
-    Double_t taa   = fgWStaa->Eval(b);
-    
-    Double_t y    = 2. * TMath::Pi() * b * sigma * taa; // fm
-    return y;
+  //
+  // Return number of binary collisions at b
+  //
+  if(b==0) b=1e-4;
+  return fgWSN->Eval(b);
 }
 
-Double_t AliFastGlauber::WSN(Double_t* x, Double_t* /*par*/)
+Double_t AliFastGlauber::Participants(Double_t  b)
 {
-//
-//  Number of hard processes per event
-//
-    Double_t b     = x[0];
-    Double_t y     = fgWSbinary->Eval(b)/fgWSgeo->Eval(b);
-    return y;
+  //
+  // Return the number of participants normalized to 1 at b=0
+  //
+  if(b==0) b=1e-4;
+  return (fgWParticipants->Eval(b)/fgWParticipants->Eval(1e-4));
 }
 
-Double_t AliFastGlauber::WEnergyDensity(Double_t* x, Double_t* par)
+Double_t AliFastGlauber::GetNumberOfParticipants(Double_t  b)
 {
-//
-//  Initial energy density as a function of the impact parameter
-//
-    Double_t b     = x[0];
-    Double_t rA    = par[0];
-//
-//  Attention: area of transverse reaction zone in hard-sphere approximation !     
-    Double_t saa   = (TMath::Pi() - 2. * TMath::ASin(b/ 2./ rA)) * rA * rA 
-       - b * TMath::Sqrt(rA * rA - b * b/ 4.);
-    Double_t taa   = fgWStaa->Eval(b);
-    
-    return (taa/saa);
+  //
+  // Return the number of participants for impact parameter b
+  //
+  if(b==0) b=1e-4;
+  return (fgWParticipants->Eval(b));
 }
 
 void AliFastGlauber::SimulateTrigger(Int_t n)
 {
-    //
-    //  Simulates Trigger
-    //
-    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.);   
+  //
+  //  Simulates Trigger
+  //
+  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");    
+  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);
+  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++)
+  //
+  //
+  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();
+      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)
 {
-    //
-    // Gives back a random impact parameter, hard trigger probability and multiplicity
-    //
-       b = fgWSgeo->GetRandom();
-       Float_t mu = fgWSN->Eval(b);
-       p = 1.-TMath::Exp(-mu);
-       mult = 6000./fgWSN->Eval(1.) * mu;
+  //
+  // 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;
 }
 
 void AliFastGlauber::GetRandom(Int_t& bin, Bool_t& hard)
 {
-    //
-    // Gives back a random impact parameter bin, and hard trigger decission
-    //
-       Float_t b  = fgWSgeo->GetRandom();
-       Float_t mu = fgWSN->Eval(b) * fSigmaHard;
-       Float_t p  = 1.-TMath::Exp(-mu);
-       if (b < 5.) {
-           bin = 1;
-       } else if (b <  8.6) {
-           bin = 2;
-       } else if (b < 11.2) {
-           bin = 3;
-       } else if (b < 13.2) {
-           bin = 4;
-       } else if (b < 15.0) {
-           bin = 5;
-       } else {
-           bin = 6;
-       }
-       
-       hard = kFALSE;
-       
-       Float_t r = gRandom->Rndm();
-       
-       if (r < p) hard = kTRUE;
+  //
+  // 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.) {
+    bin = 1;
+  } else if (b <  8.6) {
+    bin = 2;
+  } else if (b < 11.2) {
+    bin = 3;
+  } else if (b < 13.2) {
+    bin = 4;
+  } else if (b < 15.0) {
+    bin = 5;
+  } else {
+    bin = 6;
+  }
+  hard = kFALSE;
+  const Float_t r = gRandom->Rndm();
+  if (r < p) hard = kTRUE;
 }
 
-
-
-Float_t  AliFastGlauber::GetRandomImpactParameter(Float_t bmin, Float_t bmax)
+Double_t  AliFastGlauber::GetRandomImpactParameter(Double_t bmin, Double_t bmax)
 {
-    //
-    // Gives back a random impact parameter in the range bmin .. bmax
-    //
-
-    Float_t b = -1.;
-    while(b < bmin || b > bmax)
-       b = fgWSgeo->GetRandom();
-    return b;
+  //
+  // Gives back a random impact parameter in the range bmin .. bmax
+  //
+  Float_t b = -1.;
+  while(b < bmin || b > bmax)
+    b = fgWSgeo->GetRandom();
+  return b;
 }
 
-Double_t AliFastGlauber::CrossSection(Double_t b1, Double_t b2)
+void AliFastGlauber::StoreFunctions()
 {
-    //
-    // Return cross-section integrated from b1 to b2 
-    //
-    
-    return fgWSgeo->Integral(b1, b2)/100.;
+  //
+  // Store in file functions
+  //
+  TFile* ff = new TFile(fName.Data(),"recreate");
+  fgWStaa->Write("WStaa");
+  fgWParticipants->Write("WParticipants");
+  ff->Close();
+  return;
 }
 
-Float_t AliFastGlauber::GetNumberOfParticipants(Float_t  b)
-{
-    //
-    // Gives back the number of participants for impact parameter b
-    //
-    
-    return (fgWParticipants->Eval(b));
-}
+//=================== Added by A. Dainese 11/02/04 ===========================
 
-Double_t AliFastGlauber::FractionOfHardCrossSection(Double_t b1, Double_t b2)
+void AliFastGlauber::PlotAlmonds()
 {
-    //
-    // Return raction of hard cross-section integrated from b1 to b2 
-    //
-    
-    return fgWSbinary->Integral(b1, b2)/fgWSbinary->Integral(0., 100.);
+  //
+  // 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;
 }
 
-
-Double_t AliFastGlauber::Binaries(Double_t b)
+void AliFastGlauber::StoreAlmonds()
 {
-    //
-    // Return number of binary collisions normalized to 1 at b=0
-    //
-    
-    return fgWSN->Eval(b)/fgWSN->Eval(0.001);
+  //
+  // Store in file 
+  // 40 almonds for b = (0.25+k*0.5) fm (k=0->39)
+  //
+  Char_t almondName[100];
+  TFile* ff = new TFile(fName.Data(),"update");
+  for(Int_t k=0; k<40; k++) {
+    sprintf(almondName,"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;
 }
 
-//=================== Added by A. Dainese 11/02/04 ===========================
 void AliFastGlauber::SetCentralityClass(Double_t xsecFrLow,Double_t xsecFrUp)
 {
   //
-  // Set limits of centrality class as fractions of the geomtrical 
-  // cross section  
+  // Set limits of centrality class as fractions 
+  // of the geomtrical cross section  
   //
   if(xsecFrLow>1. || xsecFrUp>1. || xsecFrLow>xsecFrUp) {
-    printf(" Please set 0 <= xsecFrLow <= xsecFrUp <= 1\n");
+    Error("SetCentralityClass", "Please set 0 <= xsecFrLow <= xsecFrUp <= 1\n");
     return;
   }
 
   Double_t bLow=0.,bUp=0.;
   Double_t xsecFr=0.;
-
+  const Double_t norm=fgWSgeo->Integral(0.,100.);
   while(xsecFr<xsecFrLow) {
-    xsecFr = fgWSgeo->Integral(0.,bLow)/fgWSgeo->Integral(0.,100.);
+    xsecFr = fgWSgeo->Integral(0.,bLow)/norm;
     bLow += 0.1;
   }
   bUp = bLow;
   while(xsecFr<xsecFrUp) {
-    xsecFr = fgWSgeo->Integral(0.,bUp)/fgWSgeo->Integral(0.,100.);
+    xsecFr = fgWSgeo->Integral(0.,bUp)/norm;
     bUp += 0.1;
   }
 
-  printf(" 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\n",
         xsecFrLow,xsecFrUp,bLow,bUp);
-
   fgWSbinary->SetRange(bLow,bUp);
-
   return;
 }
 
@@ -893,7 +1198,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++;
-  fWAlmondCurrent = &fWAlmondFixedB[bin]; 
+  fgWAlmondCurrent = &fgWAlmondFixedB[bin]; 
   return;
 }
 
@@ -903,7 +1208,7 @@ void AliFastGlauber::GetRandomXY(Double_t& x,Double_t& y)
   // Get random position of parton production point according to 
   // product of thickness functions
   //
-  fWAlmondCurrent->GetRandom2(x,y);
+  fgWAlmondCurrent->GetRandom2(x,y);
   return;
 }
 
@@ -916,8 +1221,7 @@ void AliFastGlauber::GetRandomPhi(Double_t& phi)
   return;
 }
 
-Double_t AliFastGlauber::CalculateLength(Double_t b,
-                                        Double_t x0,Double_t y0,Double_t phi0)
+Double_t AliFastGlauber::CalculateLength(Double_t b,Double_t x0,Double_t y0,Double_t phi0)
 {
   // 
   // Calculate path length for a parton with production point (x0,y0)
@@ -929,86 +1233,72 @@ Double_t AliFastGlauber::CalculateLength(Double_t b,
   const Int_t    kNp  = 100;
   const Double_t kDl  = fgBMax/Double_t(kNp);
 
-  Double_t l,integral,integral1,integral2,r0,xx,yy,phi,r1,r2,ell;
-  Double_t rhocoll,rhocollHalfMax;
-  Int_t i,nps;
-
   if(fEllDef==1) {
     //
     // Definition 1:
     // 
-    // ell = 2 * \int_0^\infty dl*l*\rho_coll(x0+l*ux,y0+l*uy) /
-    //           \int_0^\infty dl*\rho_coll(x0+l*ux,y0+l*uy)
+    // ell = 2 * \int_0^\infty dl*l*(T_A*T_B)(x0+l*ux,y0+l*uy) /
+    //           \int_0^\infty dl*(T_A*T_B)(x0+l*ux,y0+l*uy)
     //
-    
 
-    l  = 0.;
-    integral1 = 0.;
-    integral2 = 0.;
-    
     // Initial radius
-    r0  = TMath::Sqrt(x0 * x0 + y0 * y0);
-    nps = Int_t ((fgBMax - r0)/kDl) - 1;
-    
+    const Double_t r0  = TMath::Sqrt(x0 * x0 + y0 * y0);
+    const Int_t nps = Int_t ((fgBMax - r0)/kDl) - 1;
+    Double_t l  = 0.;
+    Double_t integral1 = 0.;
+    Double_t integral2 = 0.;
     // Radial steps
-    for (i = 0; i < nps; i++) {
+    for (Int_t i = 0; i < nps; i++) {
       
       // Transform into target frame
-      xx   = x0 + l * TMath::Cos(phi0) + b / 2.;
-      yy   = y0 + l * TMath::Sin(phi0);
-      phi  = TMath::ATan2(yy, xx);
-      r1   = TMath::Sqrt(xx * xx + yy * yy);
+      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);
       // Radius in projectile frame
-      r2   = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi)); 
-      rhocoll = fgWSta->Eval(r1) * fgWSta->Eval(r2);
+      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);
       
-      integral1 += rhocoll * l * kDl;
-      integral2 += rhocoll * kDl;
+      integral1 += prodTATB * l * kDl;
+      integral2 += prodTATB * kDl;
       l  += kDl;
     } // steps
     
-    ell = (2. * integral1 / integral2);
+    Double_t ell=0.;
+    if(integral2)
+      ell = (2. * integral1 / integral2);
     return ell;
-    
   } else if(fEllDef==2) {
     //
     // Definition 2:
     // 
     // ell = \int_0^\infty dl*
-    //       \Theta(\rho_coll(x0+l*ux,y0+l*uy)-0.5*\rho_coll(0,0))
+    //       \Theta((T_A*T_B)(x0+l*ux,y0+l*uy)-0.5*(T_A*T_B)(0,0))
     //
 
-    l  = 0.;
-    integral = 0.;
-    
     // Initial radius
-    r0  = TMath::Sqrt(x0 * x0 + y0 * y0);
-    nps = Int_t ((fgBMax - r0)/kDl) - 1;
-
-    rhocollHalfMax = 0.5*fWAlmondCurrent->Eval(0.,0.);
-
+    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.);
     // Radial steps
-    for (i = 0; i < nps; i++) {
-      
+    Double_t l  = 0.;
+    Double_t integral = 0.;
+    for (Int_t i = 0; i < nps; i++) {
       // Transform into target frame
-      xx   = x0 + l * TMath::Cos(phi0) + b / 2.;
-      yy   = y0 + l * TMath::Sin(phi0);
-      phi  = TMath::ATan2(yy, xx);
-      r1   = TMath::Sqrt(xx * xx + yy * yy);
+      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);
       // Radius in projectile frame
-      r2   = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi)); 
-      rhocoll = fgWSta->Eval(r1) * fgWSta->Eval(r2);
-      
-      if(rhocoll>rhocollHalfMax) integral += kDl;
-
+      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;
       l  += kDl;
     } // steps
-    
-    ell = integral;
+    Double_t ell = integral;
     return ell;
-
   } else {
-    printf(" Wrong length definition setting!\n");
+    Error("CalculateLength","Wrong length definition setting: %d !\n",fEllDef);
     return -1.;
   }
 }
@@ -1022,14 +1312,11 @@ void AliFastGlauber::GetLength(Double_t& ell,Double_t b)
   if(b<0.) GetRandomBHard(b);
   GetRandomXY(x0,y0);
   GetRandomPhi(phi0);
-
   ell = CalculateLength(b,x0,y0,phi0);
-
   return;
 }
 
-void AliFastGlauber::GetLengthsBackToBack(Double_t& ell1,Double_t& ell2,
-                                         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 
@@ -1038,11 +1325,9 @@ void AliFastGlauber::GetLengthsBackToBack(Double_t& ell1,Double_t& ell2,
   if(b<0.) GetRandomBHard(b);
   GetRandomXY(x0,y0);
   GetRandomPhi(phi0);
-  Double_t phi0plusPi = phi0+TMath::Pi();
-
+  const Double_t phi0plusPi = phi0+TMath::Pi();
   ell1 = CalculateLength(b,x0,y0,phi0);
   ell2 = CalculateLength(b,x0,y0,phi0plusPi);
-
   return;
 }
 
@@ -1055,14 +1340,12 @@ void AliFastGlauber::GetLengthsForPythia(Int_t n,Double_t* phi,Double_t* ell, Do
   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)
-{
+{  
   // 
   // Plot distribution of n impact parameters
   //
@@ -1071,16 +1354,13 @@ void AliFastGlauber::PlotBDistr(Int_t n)
   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;
 }
 
@@ -1094,13 +1374,11 @@ void AliFastGlauber::PlotLengthDistr(Int_t n,Bool_t save,Char_t *fname)
   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();
@@ -1122,19 +1400,15 @@ void AliFastGlauber::PlotLengthB2BDistr(Int_t n,Bool_t save,Char_t *fname)
   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();
@@ -1143,44 +1417,203 @@ void AliFastGlauber::PlotLengthB2BDistr(Int_t n,Bool_t save,Char_t *fname)
   return;
 }
 
-void AliFastGlauber::StoreAlmonds()
+//=================== 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)
 {
+  // 
+  // 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)
   //
-  // Store in glauberPbPb.root 40 almonds for b = (0.25+k*0.5) fm (k=0->39)
+  // for a parton with production point (x0,y0)
+  // and propagation direction (ux=cos(phi0),uy=sin(phi0)) 
+  // in a collision with impact parameter b
+  // 
+
+  // number of steps in l
+  const Int_t    kNp  = 100;
+  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;
+    
+  // Radial steps
+  Double_t l  = 0.;
+  integral0 = 0.;
+  integral1 = 0.;
+  Int_t i = 0;
+  while((i < nps) && (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);
+    // 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;
+    l  += kDl;
+    i++;
+  } // steps
+  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 x0,y0,phi0;
+  if(b<0.) GetRandomBHard(b);
+  GetRandomXY(x0,y0);
+  GetRandomPhi(phi0);
+  CalculateI0I1(integral0,integral1,b,x0,y0,phi0,ellCut);
+  return;
+}
 
-  Char_t almondName[100];
-  TFile* ff = new TFile("$(ALICE_ROOT)/FASTSIM/data/glauberPbPb.root","update");
-  for(Int_t k=0; k<40; k++) {
-    sprintf(almondName,"WAlmondFixedB%d",k);
-    Double_t b = 0.25+k*0.5;
-    printf(" b = %f\n",b);
-    fgWAlmond->SetParameter(0,b);
+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 x0,y0,phi0;
+  if(b<0.) GetRandomBHard(b);
+  GetRandomXY(x0,y0);
+  GetRandomPhi(phi0);
+  const Double_t phi0plusPi = phi0+TMath::Pi();
+  CalculateI0I1(integral01,integral11,b,x0,y0,phi0,ellCut);
+  CalculateI0I1(integral02,integral12,b,x0,y0,phi0plusPi,ellCut);
+  return;
+}
 
-    fgWAlmond->Write(almondName);
+void AliFastGlauber::GetI0I1ForPythia(Int_t n,Double_t* phi,
+                                     Double_t* integral0,Double_t* integral1,
+                                     Double_t ellCut,Double_t b)
+{
+  //
+  // Returns I0 and I1 pairs for n partons with azimuthal angles phi[n] 
+  // from random b, 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);
+  return;
+}
+
+void AliFastGlauber::PlotI0I1Distr(Int_t n,Double_t ellCut,
+                                  Bool_t save,Char_t *fname)
+{
+  //
+  // Plot length distribution
+  //
+  Double_t i0,i1;
+  TH1F *hI0 = new TH1F("hI0","I_{0} = #hat{q}L / k",
+                      100,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); 
+  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);
+    hI0->Fill(i0);
+    hI1->Fill(i1);
+    h2->Fill(2.*i1*i1/i0);
+    h3->Fill(2.*i1/i0);
+    h4->Fill(i0*i0/2./i1);
   }
-  ff->Close();
+  hI0->Scale(1/(Double_t)n);
+  hI1->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();
 
+  if(save) {
+    TFile *f = new TFile(fname,"recreate");
+    hI0->Write();
+    hI1->Write();
+    h2->Write();
+    h3->Write();
+    h4->Write();
+    f->Close();
+  }
   return;
 }
 
-void AliFastGlauber::PlotAlmonds()
+void AliFastGlauber::PlotI0I1B2BDistr(Int_t n,Double_t ellCut,
+                                     Bool_t save,Char_t *fname)
 {
   //
-  // Plot almonds for some impact parameters
+  // Plot lengths 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); 
+  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}]");
 
-  TCanvas *c = new TCanvas("c","Almonds",0,0,500,500);
+  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);
-  c->Divide(2,2);
-  c->cd(1);
-  fWAlmondFixedB[0].Draw("cont1");
-  c->cd(2);
-  fWAlmondFixedB[10].Draw("cont1");
-  c->cd(3);
-  fWAlmondFixedB[20].Draw("cont1");
-  c->cd(4);
-  fWAlmondFixedB[30].Draw("cont1");
+  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;
 }
index 27f6696..1da9490 100644 (file)
 // Author: andreas.morsch@cern.ch
 
 #include <TObject.h>
+#include <TString.h>
 #include <TF2.h>
 class TF1;
-class TF2;
 
 class AliFastGlauber : public TObject {
  public:
     AliFastGlauber();
-    virtual ~AliFastGlauber(){;}
+    virtual ~AliFastGlauber();
+    void Init(Int_t mode = 0);
+
     void SetWoodSaxonParameters(Double_t r0, Double_t d, Double_t w, Double_t n)
        {fWSr0 = r0; fWSd = d; fWSw = w; fWSn = n;}
+    void SetWoodSaxonParametersAu()
+       {fWSr0 = 6.38; fWSd = 0.535; fWSw = 0.; fWSn = 8.59e-4;}
+    void SetWoodSaxonParametersPb()
+       {fWSr0 = 6.624; fWSd = 0.549; fWSw = 0.; fWSn = 7.69e-4;}
     void SetMaxImpact(Float_t bmax = 20.) {fgBMax = bmax;};
-    void SetHardCrossSection(Float_t xs = 6.6) {fSigmaHard = xs;}
+    void SetHardCrossSection(Float_t xs = 1.0) {fSigmaHard = xs;}
+    void SetNNCrossSection  (Float_t xs = 55.6) {fSigmaNN = xs;}
+    void SetNucleus(Int_t n=208) {fA=n;}
+    void SetAuAuRhic();
+    void SetPbPbLHC();
+    void SetFileName(TString &fn){fName=fn;}
+    void SetFileName(Char_t *fn="$(ALICE_ROOT)/FASTSIM/data/glauberPbPb.root"){fName=fn;}
 
     static Double_t WSb            (Double_t *xx, Double_t *par);
     static Double_t WSbz           (Double_t *xx, Double_t *par);
     static Double_t WSz            (Double_t *xx, Double_t *par);
     static Double_t WSta           (Double_t *xx, Double_t *par);
     static Double_t WStarfi        (Double_t *xx, Double_t *par);
+    static Double_t WStaa          (Double_t *xx, Double_t *par);
     static Double_t WKParticipants (Double_t *xx, Double_t *par);
     static Double_t WParticipants  (Double_t *xx, Double_t *par);    
-    static Double_t WStaa          (Double_t *xx, Double_t *par);
     static Double_t WSgeo          (Double_t *xx, Double_t *par);
     static Double_t WSbinary       (Double_t *xx, Double_t *par);
     static Double_t WSN            (Double_t *xx, Double_t *par);
@@ -40,8 +52,27 @@ class AliFastGlauber : public TObject {
     static Double_t WPathLength    (Double_t *xx, Double_t *par);
     static Double_t WIntRadius     (Double_t *xx, Double_t *par);
     static Double_t WEnergyDensity (Double_t *xx, Double_t *par);
+
+    const TF2* Kernel()            const {return fgWStarfi;}
+    const TF1* Overlap()           const {return fgWStaa;}
+    const TF1* GetWSB()            const {return fgWSb;}
+    const TF2* GetWSbz()           const {return fgWSbz;}
+    const TF1* GetWSz()            const {return fgWSz;} 
+    const TF1* GetWSta()           const {return fgWSta;}
+    const TF2* GetWStarfi()        const {return fgWStarfi;}
+    const TF2* GetWKParticipants() const {return fgWKParticipants;}
+    const TF1* GetWParticipants()  const {return fgWParticipants;} 
+    const TF1* GetWStaa()          const {return fgWStaa;} 
+    const TF2* GetWAlmond()        const {return fgWAlmond;}
+    const TF1* GetWPathLength0()   const {return fgWPathLength0;} 
+    const TF1* GetWPathLength()    const {return fgWPathLength;}
+    const TF1* GetWIntRadius()     const {return fgWIntRadius;}
+    const TF1* GetWSgeo()          const {return fgWSgeo;}
+    const TF1* GetWSbinary()       const {return fgWSbinary;}
+    const TF1* GetWSN()            const {return fgWSN;}     
+    const TF1* GetWEnergyDensity() const {return fgWEnergyDensity;} 
+    const TF2* GetWAlmondFixedB(Int_t i) const {return &fgWAlmondFixedB[i];}
     
-    void Init(Int_t mode = 0);
     void DrawWSb();
     void DrawThickness();
     void DrawOverlap();
@@ -57,18 +88,23 @@ class AliFastGlauber : public TObject {
     void DrawEnergyDensity();
     
     Double_t CrossSection(Double_t b1, Double_t b2);
+    Double_t HardCrossSection(Double_t b1, Double_t b2);
     Double_t FractionOfHardCrossSection(Double_t b1, Double_t b2);
     Double_t Binaries(Double_t b);
-    TF2* Kernel()  {return fgWStarfi;}
-    TF1* Overlap() {return fgWStaa;}
+    Double_t GetNumberofBinaries(Double_t b);
+    Double_t Participants(Double_t b);
+    Double_t GetNumberOfParticipants(Double_t  b);
+
     void SimulateTrigger(Int_t n);
     void GetRandom(Float_t& b, Float_t& p, Float_t& mult);
     void GetRandom(Int_t& bin, Bool_t& hard);
-    Float_t GetRandomImpactParameter(Float_t bmin, Float_t bmax);
-    Float_t GetNumberOfParticipants(Float_t b);
-    void SetLengthDefinition(Int_t def=1) { fEllDef=def; }
-    void SetCentralityClass(Double_t xsecFrLow=0.0,Double_t xsecFrUp=0.1);    
+    Double_t GetRandomImpactParameter(Double_t bmin, Double_t bmax);
+    void StoreFunctions();
+
     void StoreAlmonds();
+    void PlotAlmonds();
+    void SetLengthDefinition(Int_t def=1) {fEllDef=def;}
+    void SetCentralityClass(Double_t xsecFrLow=0.0,Double_t xsecFrUp=0.1);    
     void GetRandomBHard(Double_t& b);
     void GetRandomXY(Double_t& x,Double_t& y);
     void GetRandomPhi(Double_t& phi);
@@ -83,9 +119,30 @@ class AliFastGlauber : public TObject {
                         Char_t *fname="length.root");
     void PlotLengthB2BDistr(Int_t n=1000,Bool_t save=kFALSE,
                            Char_t *fname="lengthB2B.root");
-    void PlotAlmonds();
+    void CalculateI0I1(Double_t& integral0,Double_t& integral1,
+                      Double_t b=0.,
+                      Double_t x0=0.,Double_t y0=0.,Double_t phi0=0.,
+                      Double_t ellCut=20.);
+    void GetI0I1(Double_t& integral0,Double_t& integral1,
+                Double_t ellCut=20.,Double_t b=-1.);
+    void GetI0I1BackToBack(Double_t& integral01,Double_t& integral11,
+                          Double_t& integral02,Double_t& integral12,
+                          Double_t ellCut=20.,Double_t b=-1.);
+    void GetI0I1ForPythia(Int_t n,Double_t* phi,
+                         Double_t* integral0,Double_t* integral1,
+                         Double_t ellCut=20.,Double_t b=-1.);
+    void PlotI0I1Distr(Int_t n=1000,Double_t ellCut=20.,Bool_t save=kFALSE,
+                      Char_t *fname="i0i1.root");
+    void PlotI0I1B2BDistr(Int_t n=1000,Double_t ellCut=20.,Bool_t save=kFALSE,
+                         Char_t *fname="i0i1B2B.root");
 
  protected:
+    void Reset();
+
+    static Float_t fgBMax;        // Maximum Impact Parameter
+    static Int_t fgCounter;       // Counter to protect double instantiation
+    static const Int_t fgkMCInts; // Number of MC integrations
+
     static TF1*    fgWSb;            // Wood-Saxon Function (b)
     static TF2*    fgWSbz;           // Wood-Saxon Function (b, z)
     static TF1*    fgWSz;            // Wood-Saxon Function (b = b0, z)
@@ -102,18 +159,20 @@ class AliFastGlauber : public TObject {
     static TF1*    fgWSbinary;       // dSigma/db binary
     static TF1*    fgWSN;            // dN/db binary
     static TF1*    fgWEnergyDensity; // Energy density as a function of impact parameter
-    TF2  fWAlmondFixedB[40]; // Interaction Almonds read from file
-    TF2*    fWAlmondCurrent; // Interaction Almond used for length
-    
-    Float_t fWSr0;           // Wood-Saxon Parameter r0
-    Float_t fWSd;            // Wood-Saxon Parameter d
-    Float_t fWSw;            // Wood-Saxon Parameter w
-    Float_t fWSn;            // Wood-Saxon Parameter n
-    Float_t fSigmaHard;      // Hard Cross Section
-    static Float_t fgBMax;   // Maximum Impact Parameter
+    static TF2     fgWAlmondFixedB[40]; // Interaction Almonds read from file
+    static TF2*    fgWAlmondCurrent;    // Interaction Almond used for length
     
-    Int_t fEllDef;           // definition of length (see CalculateLength())
+    Float_t fWSr0;      // Wood-Saxon Parameter r0
+    Float_t fWSd;       // Wood-Saxon Parameter d
+    Float_t fWSw;       // Wood-Saxon Parameter w
+    Float_t fWSn;       // Wood-Saxon Parameter n
+                        // (chosen such that integral is one)
+    Float_t fSigmaHard; // Hard Cross Section [mbarn]
+    Float_t fSigmaNN;   // NN Cross Section [mbarn]   
+    Int_t fA;           // Nucleon number of nucleus A
 
+    Int_t fEllDef;      // definition of length (see CalculateLength())
+    TString fName;     // filename of stored distributions
     ClassDef(AliFastGlauber,1) // Event geometry simulation in the Glauber Model
 };