FASTSIM/AliFastGlauber.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17 //
  18 // Utility class to make simple Glauber type calculations
  19 //           for SYMMETRIC collision geometries (AA):
  20 // Impact parameter, production points, reaction plane dependence
  21 //
  22 // The SimulateTrigger method can be used for simple MB and hard-process
  23 // (binary scaling) trigger studies.
  24 //
  25 // Some basic quantities can be visualized directly.
  26 //
  27 // The default set-up for PbPb or AUAu collisions can be read from a file
  28 // calling Init(1) or Init(2) if you want to read Almonds too.
  29 //
  30 // ***** If you change settings dont forget to call init afterwards, *****
  31 // ***** in order to update the formulas with the new parameters.    *****
  32 //
  33 // Author: andreas.morsch@cern.ch
  34 //=================== Added by A. Dainese 11/02/04 ===========================
  35 // Calculate path length for a parton with production point (x0,y0)
  36 // and propagation direction (ux=cos(phi0),uy=sin(phi0))
  37 // in a collision with impact parameter b and functions that make use
  38 // of it.
  39 //=================== Added by A. Dainese 05/03/04 ===========================
  40 // Calculation of line integrals I0 and I1
  41 //  integral0 =    \int_0^ellCut dl*(T_A*T_B)(x0+l*ux,y0+l*uy)
  42 //  integral1 =    \int_0^ellCut dl*l*(T_A*T_B)(x0+l*ux,y0+l*uy)
  43 // mostly for use in the Quenching class
  44 //=================== Added by C. Loizdes 27/03/04 ===========================
  45 // Handling of AuAu collisions
  46 // More get/set functions
  47 // Comments, units and clearing of code
  48 //
  49
  50 // from AliRoot
  51 #include "AliFastGlauber.h"
  52 // from root
  53 #include <Riostream.h>
  54 #include <TCanvas.h>
  55 #include <TF1.h>
  56 #include <TF2.h>
  57 #include <TFile.h>
  58 #include <TH1F.h>
  59 #include <TH2F.h>
  60 #include <TLegend.h>
  61 #include <TMath.h>
  62 #include <TROOT.h>
  63 #include <TRandom.h>
  64 #include <TStyle.h>
  65
  66 ClassImp(AliFastGlauber)
  67
  68 Float_t AliFastGlauber::fgBMax           = 0.;
  69 TF1*    AliFastGlauber::fgWSb            = NULL;
  70 TF2*    AliFastGlauber::fgWSbz           = NULL;
  71 TF1*    AliFastGlauber::fgWSz            = NULL;
  72 TF1*    AliFastGlauber::fgWSta           = NULL;
  73 TF2*    AliFastGlauber::fgWStarfi        = NULL;
  74 TF2*    AliFastGlauber::fgWAlmond        = NULL;
  75 TF1*    AliFastGlauber::fgWStaa          = NULL;
  76 TF1*    AliFastGlauber::fgWSgeo          = NULL;
  77 TF1*    AliFastGlauber::fgWSbinary       = NULL;
  78 TF1*    AliFastGlauber::fgWSN            = NULL;
  79 TF1*    AliFastGlauber::fgWPathLength0   = NULL;
  80 TF1*    AliFastGlauber::fgWPathLength    = NULL;
  81 TF1*    AliFastGlauber::fgWEnergyDensity = NULL;
  82 TF1*    AliFastGlauber::fgWIntRadius     = NULL;
  83 TF2*    AliFastGlauber::fgWKParticipants = NULL;
  84 TF1*    AliFastGlauber::fgWParticipants  = NULL;
  85 TF2*    AliFastGlauber::fgWAlmondCurrent = NULL;
  86 TF2*    AliFastGlauber::fgWAlmondFixedB[40];
  87 const Int_t AliFastGlauber::fgkMCInts = 100000;
  88 AliFastGlauber* AliFastGlauber::fgGlauber = NULL;
  89
  90
  91 AliFastGlauber::AliFastGlauber():
  92     fWSr0(0.),
  93     fWSd(0.),
  94     fWSw(0.),
  95     fWSn(0.),
  96     fSigmaHard(0.),
  97     fSigmaNN(0.),
  98     fA(0),
  99     fBmin(0.),
 100     fBmax(0.),
 101     fEllDef(0),
 102     fName()
 103 {
 104   //  Default Constructor
 105   //  Defaults for Pb
 106   SetMaxImpact();
 107   SetLengthDefinition();
 108   SetPbPbLHC();
 109   fXY[0] = fXY[1] = 0;
 110   fI0I1[0] = fI0I1[1] = 0;
 111 }
 112
 113 AliFastGlauber::AliFastGlauber(const AliFastGlauber & gl)
 114     :TObject(gl),
 115      fWSr0(0.),
 116      fWSd(0.),
 117      fWSw(0.),
 118      fWSn(0.),
 119      fSigmaHard(0.),
 120      fSigmaNN(0.),
 121      fA(0),
 122      fBmin(0.),
 123      fBmax(0.),
 124      fEllDef(0),
 125      fName()
 126 {
 127 // Copy constructor
 128     gl.Copy(*this);
 129     fXY[0] = fXY[1] = 0;
 130     fI0I1[0] = fI0I1[1] = 0;
 131 }
 132
 133 AliFastGlauber* AliFastGlauber::Instance()
 134 {
 135 // Set random number generator
 136     if (fgGlauber) {
 137         return fgGlauber;
 138     } else {
 139         fgGlauber = new AliFastGlauber();
 140         return fgGlauber;
 141     }
 142 }
 143
 144 AliFastGlauber::~AliFastGlauber()
 145 {
 146 // Destructor
 147   for(Int_t k=0; k<40; k++) delete fgWAlmondFixedB[k];
 148 }
 149
 150 void AliFastGlauber::SetAuAuRhic()
 151 {
 152   //Set all parameters for RHIC
 153   SetWoodSaxonParametersAu();
 154   SetHardCrossSection();
 155   SetNNCrossSection(42);
 156   SetNucleus(197);
 157   SetFileName("$(ALICE_ROOT)/FASTSIM/data/glauberAuAu.root");
 158 }
 159
 160 void AliFastGlauber::SetPbPbLHC()
 161 {
 162   //Set all parameters for LHC
 163   SetWoodSaxonParametersPb();
 164   SetHardCrossSection();
 165   SetNNCrossSection();
 166   SetNucleus();
 167   SetFileName();
 168 }
 169
 170 void AliFastGlauber::Init(Int_t mode)
 171 {
 172   // Initialisation
 173   //    mode = 0; all functions are calculated
 174   //    mode = 1; overlap function is read from file (for Pb-Pb only)
 175   //    mode = 2; interaction almond functions are read from file
 176   //              USE THIS FOR PATH LENGTH CALC.!
 177   //
 178
 179   //
 180   Reset();
 181   //
 182
 183   //
 184   //  Wood-Saxon
 185   //
 186   fgWSb = new TF1("WSb", WSb, 0, fgBMax, 4);
 187   fgWSb->SetParameter(0, fWSr0);
 188   fgWSb->SetParameter(1, fWSd);
 189   fgWSb->SetParameter(2, fWSw);
 190   fgWSb->SetParameter(3, fWSn);
 191
 192   fgWSbz = new TF2("WSbz", WSbz, 0, fgBMax, 0, fgBMax, 4);
 193   fgWSbz->SetParameter(0, fWSr0);
 194   fgWSbz->SetParameter(1, fWSd);
 195   fgWSbz->SetParameter(2, fWSw);
 196   fgWSbz->SetParameter(3, fWSn);
 197
 198   fgWSz = new TF1("WSz", WSz, 0, fgBMax, 5);
 199   fgWSz->SetParameter(0, fWSr0);
 200   fgWSz->SetParameter(1, fWSd);
 201   fgWSz->SetParameter(2, fWSw);
 202   fgWSz->SetParameter(3, fWSn);
 203
 204   //
 205   //  Thickness
 206   //
 207   fgWSta = new TF1("WSta", WSta, 0., fgBMax, 0);
 208
 209   //
 210   //  Overlap Kernel
 211   //
 212   fgWStarfi = new TF2("WStarfi", WStarfi, 0., fgBMax, 0., TMath::Pi(), 1);
 213   fgWStarfi->SetParameter(0, 0.);
 214   fgWStarfi->SetNpx(200);
 215   fgWStarfi->SetNpy(20);
 216
 217   //
 218   //  Participants Kernel
 219   //
 220   fgWKParticipants = new TF2("WKParticipants", WKParticipants, 0., fgBMax, 0., TMath::Pi(), 3);
 221   fgWKParticipants->SetParameter(0, 0.);
 222   fgWKParticipants->SetParameter(1, fSigmaNN);
 223   fgWKParticipants->SetParameter(2, fA);
 224   fgWKParticipants->SetNpx(200);
 225   fgWKParticipants->SetNpy(20);
 226
 227   //
 228   //  Overlap and Participants
 229   //
 230   if (!mode) {
 231     fgWStaa = new TF1("WStaa", WStaa, 0., fgBMax, 1);
 232     fgWStaa->SetNpx(100);
 233     fgWStaa->SetParameter(0,fA);
 234     fgWStaa->SetNpx(100);
 235     fgWParticipants = new TF1("WParticipants", WParticipants, 0., fgBMax, 2);
 236     fgWParticipants->SetParameter(0, fSigmaNN);
 237     fgWParticipants->SetParameter(1, fA);
 238     fgWParticipants->SetNpx(100);
 239   } else {
 240     Info("Init","Reading overlap function from file %s",fName.Data());
 241     TFile* f = new TFile(fName.Data());
 242     if(!f->IsOpen()){
 243       Fatal("Init", "Could not open file %s",fName.Data());
 244     }
 245     fgWStaa = (TF1*) f->Get("WStaa");
 246     fgWParticipants = (TF1*) f->Get("WParticipants");
 247     delete f;
 248   }
 249
 250   //
 251   //  Energy Density
 252   //
 253   fgWEnergyDensity = new TF1("WEnergyDensity", WEnergyDensity, 0., 2. * fWSr0, 1);
 254   fgWEnergyDensity->SetParameter(0, fWSr0 + 1.);
 255
 256   //
 257   //  Geometrical Cross-Section
 258   //
 259   fgWSgeo = new TF1("WSgeo", WSgeo, 0., fgBMax, 1);
 260   fgWSgeo->SetParameter(0,fSigmaNN); //mbarn
 261   fgWSgeo->SetNpx(100);
 262
 263   //
 264   //  Hard cross section (binary collisions)
 265   //
 266   fgWSbinary = new TF1("WSbinary", WSbinary, 0., fgBMax, 1);
 267   fgWSbinary->SetParameter(0, fSigmaHard); //mbarn
 268   fgWSbinary->SetNpx(100);
 269
 270   //
 271   // Hard collisions per event
 272   //
 273   fgWSN = new TF1("WSN", WSN, 0., fgBMax, 1);
 274   fgWSN->SetNpx(100);
 275
 276   //
 277   //  Almond shaped interaction region
 278   //
 279   fgWAlmond = new TF2("WAlmond", WAlmond, -fgBMax, fgBMax, -fgBMax, fgBMax, 1);
 280   fgWAlmond->SetParameter(0, 0.);
 281   fgWAlmond->SetNpx(200);
 282   fgWAlmond->SetNpy(200);
 283
 284   if(mode==2) {
 285     Info("Init","Reading interaction almonds from file: %s",fName.Data());
 286     Char_t almondName[100];
 287     TFile* ff = new TFile(fName.Data());
 288     for(Int_t k=0; k<40; k++) {
 289       sprintf(almondName,"WAlmondFixedB%d",k);
 290       fgWAlmondCurrent = (TF2*)ff->Get(almondName);
 291       fgWAlmondFixedB[k] = fgWAlmondCurrent;
 292     }
 293     delete ff;
 294   }
 295
 296   fgWIntRadius = new TF1("WIntRadius", WIntRadius, 0., fgBMax, 1);
 297   fgWIntRadius->SetParameter(0, 0.);
 298
 299   //
 300   //  Path Length as a function of Phi
 301   //
 302   fgWPathLength0 = new TF1("WPathLength0", WPathLength0, -TMath::Pi(), TMath::Pi(), 2);
 303   fgWPathLength0->SetParameter(0, 0.);
 304   fgWPathLength0->SetParameter(1, 0.); //Pathlength definition
 305
 306   fgWPathLength = new TF1("WPathLength", WPathLength, -TMath::Pi(), TMath::Pi(), 3);
 307   fgWPathLength->SetParameter(0, 0.);    //Impact Parameter
 308   fgWPathLength->SetParameter(1, 1000.); //Number of interactions used for average
 309   fgWPathLength->SetParameter(2, 0);     //Pathlength definition
 310 }
 311
 312 void AliFastGlauber::Reset() const
 313 {
 314   //
 315   // Reset dynamic allocated formulas
 316   // in case init is called twice
 317
 318   if(fgWSb)            delete fgWSb;
 319   if(fgWSbz)           delete fgWSbz;
 320   if(fgWSz)            delete fgWSz;
 321   if(fgWSta)           delete fgWSta;
 322   if(fgWStarfi)        delete fgWStarfi;
 323   if(fgWAlmond)        delete fgWAlmond;
 324   if(fgWStaa)          delete fgWStaa;
 325   if(fgWSgeo)          delete fgWSgeo;
 326   if(fgWSbinary)       delete fgWSbinary;
 327   if(fgWSN)            delete fgWSN;
 328   if(fgWPathLength0)   delete fgWPathLength0;
 329   if(fgWPathLength)    delete fgWPathLength;
 330   if(fgWEnergyDensity) delete fgWEnergyDensity;
 331   if(fgWIntRadius)     delete fgWIntRadius;
 332   if(fgWKParticipants) delete fgWKParticipants;
 333   if(fgWParticipants)  delete fgWParticipants;
 334 }
 335
 336 void AliFastGlauber::DrawWSb() const
 337 {
 338   //
 339   //  Draw Wood-Saxon Nuclear Density Function
 340   //
 341   TCanvas *c1 = new TCanvas("c1","Wood Saxon",400,10,600,700);
 342   c1->cd();
 343   Double_t max=fgWSb->GetMaximum(0,fgBMax)*1.01;
 344   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);
 345   h2f->SetStats(0);
 346   h2f->GetXaxis()->SetTitle("r [fm]");
 347   h2f->GetYaxis()->SetNoExponent(kTRUE);
 348   h2f->GetYaxis()->SetTitle("#rho [fm^{-3}]");
 349   h2f->Draw();
 350   fgWSb->Draw("same");
 351   TLegend *l1a = new TLegend(0.45,0.6,.90,0.8);
 352   l1a->SetFillStyle(0);
 353   l1a->SetBorderSize(0);
 354   Char_t label[100];
 355   sprintf(label,"r_{0} = %.2f fm",fWSr0);
 356   l1a->AddEntry(fgWSb,label,"");
 357   sprintf(label,"d = %.2f fm",fWSd);
 358   l1a->AddEntry(fgWSb,label,"");
 359   sprintf(label,"n = %.2e fm^{-3}",fWSn);
 360   l1a->AddEntry(fgWSb,label,"");
 361   sprintf(label,"#omega = %.2f",fWSw);
 362   l1a->AddEntry(fgWSb,label,"");
 363   l1a->Draw();
 364   c1->Update();
 365 }
 366
 367 void AliFastGlauber::DrawOverlap() const
 368 {
 369   //
 370   //  Draw Overlap Function
 371   //
 372   TCanvas *c2 = new TCanvas("c2","Overlap",400,10,600,700);
 373   c2->cd();
 374   Double_t max=fgWStaa->GetMaximum(0,fgBMax)*1.01;
 375   TH2F *h2f=new TH2F("h2ftaa","Overlap function: T_{AB} [mbarn^{-1}]",2,0,fgBMax,2,0, max);
 376   h2f->SetStats(0);
 377   h2f->GetXaxis()->SetTitle("b [fm]");
 378   h2f->GetYaxis()->SetTitle("T_{AB} [mbarn^{-1}]");
 379   h2f->Draw();
 380   fgWStaa->Draw("same");
 381 }
 382
 383 void AliFastGlauber::DrawParticipants() const
 384 {
 385   //
 386   //  Draw Number of Participants Npart
 387   //
 388   TCanvas *c3 = new TCanvas("c3","Participants",400,10,600,700);
 389   c3->cd();
 390   Double_t max=fgWParticipants->GetMaximum(0,fgBMax)*1.01;
 391   TH2F *h2f=new TH2F("h2fpart","Number of Participants",2,0,fgBMax,2,0,max);
 392   h2f->SetStats(0);
 393   h2f->GetXaxis()->SetTitle("b [fm]");
 394   h2f->GetYaxis()->SetTitle("N_{part}");
 395   h2f->Draw();
 396   fgWParticipants->Draw("same");
 397   TLegend *l1a = new TLegend(0.50,0.75,.90,0.9);
 398   l1a->SetFillStyle(0);
 399   l1a->SetBorderSize(0);
 400   Char_t label[100];
 401   sprintf(label,"#sigma^{inel.}_{NN} = %.1f mbarn",fSigmaNN);
 402   l1a->AddEntry(fgWParticipants,label,"");
 403   l1a->Draw();
 404   c3->Update();
 405 }
 406
 407 void AliFastGlauber::DrawThickness() const
 408 {
 409   //
 410   //  Draw Thickness Function
 411   //
 412   TCanvas *c4 = new TCanvas("c4","Thickness",400,10,600,700);
 413   c4->cd();
 414   Double_t max=fgWSta->GetMaximum(0,fgBMax)*1.01;
 415   TH2F *h2f=new TH2F("h2fta","Thickness function: T_{A} [fm^{-2}]",2,0,fgBMax,2,0,max);
 416   h2f->SetStats(0);
 417   h2f->GetXaxis()->SetTitle("b [fm]");
 418   h2f->GetYaxis()->SetTitle("T_{A} [fm^{-2}]");
 419   h2f->Draw();
 420   fgWSta->Draw("same");
 421 }
 422
 423 void AliFastGlauber::DrawGeo() const
 424 {
 425   //
 426   //  Draw Geometrical Cross-Section
 427   //
 428   TCanvas *c5 = new TCanvas("c5","Geometrical Cross-Section",400,10,600,700);
 429   c5->cd();
 430   Double_t max=fgWSgeo->GetMaximum(0,fgBMax)*1.01;
 431   TH2F *h2f=new TH2F("h2fgeo","Differential Geometrical Cross-Section: d#sigma^{geo}_{AB}/db [fm]",2,0,fgBMax,2,0,max);
 432   h2f->SetStats(0);
 433   h2f->GetXaxis()->SetTitle("b [fm]");
 434   h2f->GetYaxis()->SetTitle("d#sigma^{geo}_{AB}/db [fm]");
 435   h2f->Draw();
 436   fgWSgeo->Draw("same");
 437   TLegend *l1a = new TLegend(0.10,0.8,.40,0.9);
 438   l1a->SetFillStyle(0);
 439   l1a->SetBorderSize(0);
 440   Char_t label[100];
 441   sprintf(label,"#sigma_{NN}^{inel.} = %.1f mbarn",fSigmaNN);
 442   l1a->AddEntry(fgWSgeo,label,"");
 443   l1a->Draw();
 444   c5->Update();
 445 }
 446
 447 void AliFastGlauber::DrawBinary() const
 448 {
 449   //
 450   //  Draw Binary Cross-Section
 451   //
 452   TCanvas *c6 = new TCanvas("c6","Binary Cross-Section",400,10,600,700);
 453   c6->cd();
 454   Double_t max=fgWSbinary->GetMaximum(0,fgBMax)*1.01;
 455   TH2F *h2f=new TH2F("h2fbinary","Differential Binary Cross-Section: #sigma^{hard}_{NN} dT_{AB}/db [fm]",2,0,fgBMax,2,0,max);
 456   h2f->SetStats(0);
 457   h2f->GetXaxis()->SetTitle("b [fm]");
 458   h2f->GetYaxis()->SetTitle("d#sigma^{hard}_{AB}/db [fm]");
 459   h2f->Draw();
 460   fgWSbinary->Draw("same");
 461   TLegend *l1a = new TLegend(0.50,0.8,.90,0.9);
 462   l1a->SetFillStyle(0);
 463   l1a->SetBorderSize(0);
 464   Char_t label[100];
 465   sprintf(label,"#sigma_{NN}^{hard} = %.1f mbarn",fSigmaHard);
 466   l1a->AddEntry(fgWSb,label,"");
 467   l1a->Draw();
 468   c6->Update();
 469 }
 470
 471 void AliFastGlauber::DrawN() const
 472 {
 473   //
 474   //  Draw Binaries per event (Ncoll)
 475   //
 476   TCanvas *c7 = new TCanvas("c7","Binaries per event",400,10,600,700);
 477   c7->cd();
 478   Double_t max=fgWSN->GetMaximum(0,fgBMax)*1.01;
 479   TH2F *h2f=new TH2F("h2fhardcols","Number of hard collisions: T_{AB} #sigma^{hard}_{NN}/#sigma_{AB}^{geo}",2,0,fgBMax,2,0,max);
 480   h2f->SetStats(0);
 481   h2f->GetXaxis()->SetTitle("b [fm]");
 482   h2f->GetYaxis()->SetTitle("N_{coll}");
 483   h2f->Draw();
 484   fgWSN->Draw("same");
 485   TLegend *l1a = new TLegend(0.50,0.75,.90,0.9);
 486   l1a->SetFillStyle(0);
 487   l1a->SetBorderSize(0);
 488   Char_t label[100];
 489   sprintf(label,"#sigma^{hard}_{NN} = %.1f mbarn",fSigmaHard);
 490   l1a->AddEntry(fgWSN,label,"");
 491   sprintf(label,"#sigma^{inel.}_{NN} = %.1f mbarn",fSigmaNN);
 492   l1a->AddEntry(fgWSN,label,"");
 493   l1a->Draw();
 494   c7->Update();
 495 }
 496
 497 void AliFastGlauber::DrawKernel(Double_t b) const
 498 {
 499   //
 500   //  Draw Kernel
 501   //
 502   TCanvas *c8 = new TCanvas("c8","Kernel",400,10,600,700);
 503   c8->cd();
 504   fgWStarfi->SetParameter(0, b);
 505   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());
 506   h2f->SetStats(0);
 507   h2f->GetXaxis()->SetTitle("r [fm]");
 508   h2f->GetYaxis()->SetTitle("#phi [rad]");
 509   h2f->Draw();
 510   fgWStarfi->Draw("same");
 511   TLegend *l1a = new TLegend(0.65,0.8,.90,0.9);
 512   l1a->SetFillStyle(0);
 513   l1a->SetBorderSize(0);
 514   Char_t label[100];
 515   sprintf(label,"b = %.1f fm",b);
 516   l1a->AddEntry(fgWStarfi,label,"");
 517   l1a->Draw();
 518   c8->Update();
 519 }
 520
 521 void AliFastGlauber::DrawAlmond(Double_t b) const
 522 {
 523   //
 524   //  Draw Interaction Almond
 525   //
 526   TCanvas *c9 = new TCanvas("c9","Almond",400,10,600,700);
 527   c9->cd();
 528   fgWAlmond->SetParameter(0, b);
 529   TH2F *h2f=new TH2F("h2falmond","Interaction Almond [fm^{-4}]",2,0,fgBMax,2,0,fgBMax);
 530   h2f->SetStats(0);
 531   h2f->GetXaxis()->SetTitle("x [fm]");
 532   h2f->GetYaxis()->SetTitle("y [fm]");
 533   h2f->Draw();
 534   fgWAlmond->Draw("same");
 535   TLegend *l1a = new TLegend(0.65,0.8,.90,0.9);
 536   l1a->SetFillStyle(0);
 537   l1a->SetBorderSize(0);
 538   Char_t label[100];
 539   sprintf(label,"b = %.1f fm",b);
 540   l1a->AddEntry(fgWAlmond,label,"");
 541   l1a->Draw();
 542   c9->Update();
 543 }
 544
 545 void AliFastGlauber::DrawEnergyDensity() const
 546 {
 547   //
 548   //  Draw energy density
 549   //
 550   TCanvas *c10 = new TCanvas("c10","Energy Density",400, 10, 600, 700);
 551   c10->cd();
 552   fgWEnergyDensity->SetMinimum(0.);
 553   Double_t max=fgWEnergyDensity->GetMaximum(0,fgWEnergyDensity->GetParameter(0))*1.01;
 554   TH2F *h2f=new TH2F("h2fenergydens","Energy density",2,0,fgBMax,2,0,max);
 555   h2f->SetStats(0);
 556   h2f->GetXaxis()->SetTitle("b [fm]");
 557   h2f->GetYaxis()->SetTitle("fm^{-4}");
 558   h2f->Draw();
 559   fgWEnergyDensity->Draw("same");
 560   c10->Update();
 561 }
 562
 563 void AliFastGlauber::DrawPathLength0(Double_t b, Int_t iopt) const
 564 {
 565   //
 566   //  Draw Path Length
 567   //
 568   TCanvas *c11 = new TCanvas("c11","Path Length",400,10,600,700);
 569   c11->cd();
 570   fgWPathLength0->SetParameter(0, b);
 571   fgWPathLength0->SetParameter(1, Double_t(iopt));
 572   fgWPathLength0->SetMinimum(0.);
 573   fgWPathLength0->SetMaximum(10.);
 574   TH2F *h2f=new TH2F("h2fpathlength0","Path length",2,-TMath::Pi(), TMath::Pi(),2,0,10.);
 575   h2f->SetStats(0);
 576   h2f->GetXaxis()->SetTitle("#phi [rad]");
 577   h2f->GetYaxis()->SetTitle("l [fm]");
 578   h2f->Draw();
 579   fgWPathLength0->Draw("same");
 580 }
 581
 582 void AliFastGlauber::DrawPathLength(Double_t b , Int_t ni, Int_t iopt) const
 583 {
 584   //
 585   //  Draw Path Length
 586   //
 587   TCanvas *c12 = new TCanvas("c12","Path Length",400,10,600,700);
 588   c12->cd();
 589   fgWAlmond->SetParameter(0, b);
 590   fgWPathLength->SetParameter(0, b);
 591   fgWPathLength->SetParameter(1, Double_t (ni));
 592   fgWPathLength->SetParameter(2, Double_t (iopt));
 593   fgWPathLength->SetMinimum(0.);
 594   fgWPathLength->SetMaximum(10.);
 595   TH2F *h2f=new TH2F("h2fpathlength","Path length",2,-TMath::Pi(), TMath::Pi(),2,0,10.);
 596   h2f->SetStats(0);
 597   h2f->GetXaxis()->SetTitle("#phi [rad]");
 598   h2f->GetYaxis()->SetTitle("l [fm]");
 599   h2f->Draw();
 600   fgWPathLength->Draw("same");
 601 }
 602
 603 void AliFastGlauber::DrawIntRadius(Double_t b) const
 604 {
 605   //
 606   //  Draw Interaction Radius
 607   //
 608   TCanvas *c13 = new TCanvas("c13","Interaction Radius",400,10,600,700);
 609   c13->cd();
 610   fgWIntRadius->SetParameter(0, b);
 611   fgWIntRadius->SetMinimum(0);
 612   Double_t max=fgWIntRadius->GetMaximum(0,fgBMax)*1.01;
 613   TH2F *h2f=new TH2F("h2fintradius","Interaction Density",2,0.,fgBMax,2,0,max);
 614   h2f->SetStats(0);
 615   h2f->GetXaxis()->SetTitle("r [fm]");
 616   h2f->GetYaxis()->SetTitle("[fm^{-3}]");
 617   h2f->Draw();
 618   fgWIntRadius->Draw("same");
 619 }
 620
 621 Double_t AliFastGlauber::WSb(Double_t* x, Double_t* par)
 622 {
 623   //
 624   //  Woods-Saxon Parameterisation
 625   //  as a function of radius (xx)
 626   //
 627   const Double_t kxx  = x[0];   //fm
 628   const Double_t kr0  = par[0]; //fm
 629   const Double_t kd   = par[1]; //fm
 630   const Double_t kw   = par[2]; //no units
 631   const Double_t kn   = par[3]; //fm^-3 (used to normalize integral to one)
 632   Double_t y   = kn * (1.+kw*(kxx/kr0)*(kxx/kr0))/(1.+TMath::Exp((kxx-kr0)/kd));
 633   return y; //fm^-3
 634 }
 635
 636 Double_t AliFastGlauber::WSbz(Double_t* x, Double_t* par)
 637 {
 638   //
 639   //  Wood Saxon Parameterisation
 640   //  as a function of z and  b
 641   //
 642   const Double_t kbb  = x[0];   //fm
 643   const Double_t kzz  = x[1];   //fm
 644   const Double_t kr0  = par[0]; //fm
 645   const Double_t kd   = par[1]; //fm
 646   const Double_t kw   = par[2]; //no units
 647   const Double_t kn   = par[3]; //fm^-3 (used to normalize integral to one)
 648   const Double_t kxx  = TMath::Sqrt(kbb*kbb+kzz*kzz);
 649   Double_t y  = kn * (1.+kw*(kxx/kr0)*(kxx/kr0))/(1.+TMath::Exp((kxx-kr0)/kd));
 650   return y; //fm^-3
 651 }
 652
 653 Double_t AliFastGlauber::WSz(Double_t* x, Double_t* par)
 654 {
 655   //
 656   //  Wood Saxon Parameterisation
 657   //  as a function of z for fixed b
 658   //
 659   const Double_t kzz  = x[0];   //fm
 660   const Double_t kr0  = par[0]; //fm
 661   const Double_t kd   = par[1]; //fm
 662   const Double_t kw   = par[2]; //no units
 663   const Double_t kn   = par[3]; //fm^-3 (used to normalize integral to one)
 664   const Double_t kbb  = par[4]; //fm
 665   const Double_t kxx  = TMath::Sqrt(kbb*kbb+kzz*kzz);
 666   Double_t y  = kn * (1.+kw*(kxx/kr0)*(kxx/kr0))/(1.+TMath::Exp((kxx-kr0)/kd));
 667   return y; //fm^-3
 668 }
 669
 670 Double_t AliFastGlauber::WSta(Double_t* x, Double_t* /*par*/)
 671 {
 672   //
 673   //  Thickness function T_A
 674   //  as a function of b
 675   //
 676   const Double_t kb  = x[0];
 677   fgWSz->SetParameter(4, kb);
 678   Double_t y  = 2. * fgWSz->Integral(0., fgBMax);
 679   return y; //fm^-2
 680 }
 681
 682 Double_t AliFastGlauber::WStarfi(Double_t* x, Double_t* par)
 683 {
 684   //
 685   //  Kernel for overlap function: T_A(s)*T_A(s-b)
 686   //  as a function of r and phi
 687   const Double_t kr1  = x[0];
 688   const Double_t kphi = x[1];
 689   const Double_t kb   = par[0];
 690   const Double_t kr2  = TMath::Sqrt(kr1*kr1 + kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
 691   Double_t y = kr1 * fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
 692   return y; //fm^-3
 693 }
 694
 695 Double_t AliFastGlauber::WStaa(Double_t* x, Double_t* par)
 696 {
 697   //
 698   //  Overlap function
 699   //  T_{AB}=Int d2s T_A(s)*T_B(s-b)
 700   //  as a function of b
 701   // (normalized to fA*fB)
 702   //
 703   const Double_t kb  = x[0];
 704   const Double_t ka = par[0];
 705   fgWStarfi->SetParameter(0, kb);
 706
 707   // root integration seems to fail
 708   /*
 709      Double_t al[2];
 710      Double_t bl[2];
 711      al[0] = 1e-6;
 712      al[1] = fgBMax;
 713      bl[0] = 0.;
 714      bl[1] = TMath::Pi();
 715      Double_t err;
 716
 717      Double_t y =  2. *  208. * 208. * fgWStarfi->IntegralMultiple(2, al, bl, 0.001, err);
 718      printf("WStaa: %.5e %.5e %.5e\n", b, y, err);
 719   */
 720
 721   //
 722   //  MC Integration
 723   //
 724   Double_t y = 0;
 725
 726
 727   for (Int_t i = 0; i < fgkMCInts; i++)
 728     {
 729
 730       const Double_t kphi = TMath::Pi() * gRandom->Rndm();
 731       const Double_t kb1  = fgBMax      * gRandom->Rndm();
 732       y += fgWStarfi->Eval(kb1, kphi);
 733     }
 734   y *= 2. * TMath::Pi() * fgBMax / fgkMCInts; //fm^-2
 735   y *= ka * ka * 0.1; //mbarn^-1
 736   return y;
 737 }
 738
 739 Double_t AliFastGlauber::WKParticipants(Double_t* x, Double_t* par)
 740 {
 741   //
 742   //  Kernel for number of participants
 743   //  as a function of r and phi
 744   //
 745   const Double_t kr1   = x[0];
 746   const Double_t kphi  = x[1];
 747   const Double_t kb    = par[0]; //fm
 748   const Double_t ksig  = par[1]; //mbarn
 749   const Double_t ka    = par[2]; //mass number
 750   const Double_t kr2   = TMath::Sqrt(kr1*kr1 +kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
 751   const Double_t kxsi  = fgWSta->Eval(kr2) * ksig * 0.1; //no units
 752   /*
 753     Double_t y=(1-TMath::Power((1-xsi),aa))
 754    */
 755   Double_t a   = ka;
 756   Double_t sum = ka * kxsi;
 757   Double_t y   = sum;
 758   for (Int_t i = 1; i <= ka; i++)
 759     {
 760       a--;
 761       sum *= (-kxsi) * a / Float_t(i+1);
 762       y  += sum;
 763     }
 764   y *= kr1 * fgWSta->Eval(kr1);
 765   return y; //fm^-1
 766 }
 767
 768 Double_t AliFastGlauber::WParticipants(Double_t* x, Double_t* par)
 769 {
 770   //
 771   //  Number of Participants as
 772   //  a function of b
 773   //
 774   const Double_t kb = x[0];
 775   const Double_t ksig  = par[0]; //mbarn
 776   const Double_t ka   = par[1];  //mass number
 777   fgWKParticipants->SetParameter(0, kb);
 778   fgWKParticipants->SetParameter(1, ksig);
 779   fgWKParticipants->SetParameter(2, ka);
 780
 781   //
 782   //  MC Integration
 783   //
 784   Double_t y = 0;
 785   for (Int_t i = 0; i < fgkMCInts; i++)
 786     {
 787       const Double_t kphi = TMath::Pi() * gRandom->Rndm();
 788       const Double_t kb1  = fgBMax      * gRandom->Rndm();
 789       y += fgWKParticipants->Eval(kb1, kphi);
 790     }
 791   y *= 2. *  ka * 2. * TMath::Pi() * fgBMax / fgkMCInts;
 792   return y; //no units
 793 }
 794
 795 Double_t AliFastGlauber::WSgeo(Double_t* x, Double_t* par)
 796 {
 797   //
 798   //  Geometrical Cross-Section
 799   //  as a function of b
 800   //
 801   const Double_t kb     = x[0];              //fm
 802   const Double_t ksigNN = par[0];            //mbarn
 803   const Double_t ktaa   = fgWStaa->Eval(kb); //mbarn^-1
 804   Double_t y     = 2. * TMath::Pi() * kb * (1. - TMath::Exp(- ksigNN * ktaa));
 805   return y; //fm
 806 }
 807
 808 Double_t AliFastGlauber::WSbinary(Double_t* x, Double_t* par)
 809 {
 810   //
 811   //  Number of binary hard collisions
 812   //  as a function of b
 813   //
 814   const Double_t kb     = x[0];              //fm
 815   const Double_t ksig   = par[0];            //mbarn
 816   const Double_t ktaa   = fgWStaa->Eval(kb); //mbarn^-1
 817   Double_t y = 2. * TMath::Pi() * kb * ksig * ktaa;
 818   return y; //fm
 819 }
 820
 821 Double_t AliFastGlauber::WSN(Double_t* x, Double_t* /*par*/)
 822 {
 823   //
 824   //  Number of hard processes per event
 825   //  as a function of b
 826   const Double_t kb = x[0];
 827   Double_t y = fgWSbinary->Eval(kb)/fgWSgeo->Eval(kb);
 828   return y; //no units
 829 }
 830
 831 Double_t AliFastGlauber::WEnergyDensity(Double_t* x, Double_t* par)
 832 {
 833   //
 834   //  Initial energy density
 835   //  as a function of the impact parameter
 836   //
 837   const Double_t kb     = x[0];
 838   const Double_t krA    = par[0];
 839   //
 840   //  Attention: area of transverse reaction zone in hard-sphere approximation !
 841   const Double_t krA2=krA*krA;
 842   const Double_t kb2=kb*kb;
 843   const Double_t ksaa = (TMath::Pi() - 2. * TMath::ASin(kb/ 2./ krA)) * krA2
 844                       - kb * TMath::Sqrt(krA2 - kb2/ 4.); //fm^2
 845   const Double_t ktaa = fgWStaa->Eval(kb); //mbarn^-1
 846   Double_t y=ktaa/ksaa*10;
 847   return y; //fm^-4
 848 }
 849
 850 Double_t AliFastGlauber::WAlmond(Double_t* x, Double_t* par)
 851 {
 852   //
 853   //  Almond shaped interaction region
 854   //  as a function of cartesian x,y.
 855   //
 856   const Double_t kb    = par[0];
 857   const Double_t kxx   = x[0] + kb/2.;
 858   const Double_t kyy   = x[1];
 859   const Double_t kr1   = TMath::Sqrt(kxx*kxx + kyy*kyy);
 860   const Double_t kphi  = TMath::ATan2(kyy,kxx);
 861   const Double_t kr2   = TMath::Sqrt(kr1*kr1 + kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
 862   //
 863   //  Interaction probability calculated as product of thicknesses
 864   //
 865   Double_t y    = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
 866   return y; //fm^-4
 867 }
 868
 869 Double_t AliFastGlauber::WIntRadius(Double_t* x, Double_t* par)
 870 {
 871   //
 872   //  Average interaction density over radius
 873   //  at which interaction takes place
 874   //  as a function of radius
 875   //
 876   const Double_t kr    = x[0];
 877   const Double_t kb    = par[0];
 878   fgWAlmond->SetParameter(0, kb);
 879   //  Average over phi in small steps
 880   const Double_t kdphi = 2. * TMath::Pi() / 100.;
 881   Double_t phi  = 0.;
 882   Double_t y    = 0.;
 883   for (Int_t i = 0; i < 100; i++) {
 884     const Double_t kxx = kr * TMath::Cos(phi);
 885     const Double_t kyy = kr * TMath::Sin(phi);
 886     y   += fgWAlmond->Eval(kxx,kyy);
 887     phi += kdphi;
 888   } // phi loop
 889   // Result multiplied by Jacobian (2 pi r)
 890   y *= 2. * TMath::Pi() * kr / 100.;
 891   return y; //fm^-3
 892 }
 893
 894 Double_t AliFastGlauber::WPathLength0(Double_t* x, Double_t* par)
 895 {
 896   //
 897   //  Path Length as a function of phi
 898   //  for interaction point fixed at (0,0)
 899   //  as a function of phi-direction
 900   //
 901   //  Phi direction in Almond
 902   const Double_t kphi0   = x[0];
 903   const Double_t kb      = par[0];
 904   //  Path Length definition
 905   const Int_t    kiopt   = Int_t(par[1]);
 906
 907   //  Step along radial direction phi
 908   const Int_t    kNp  = 100; // Steps in r
 909   const Double_t kDr  = fgBMax/kNp;
 910   Double_t r  = 0.;
 911   Double_t rw = 0.;
 912   Double_t w  = 0.;
 913   for (Int_t i = 0; i < kNp; i++) {
 914     //
 915     //  Transform into target frame
 916     //
 917     const Double_t kxx   = r * TMath::Cos(kphi0) + kb / 2.;
 918     const Double_t kyy   = r * TMath::Sin(kphi0);
 919     const Double_t kphi  = TMath::ATan2(kyy, kxx);
 920     const Double_t kr1   = TMath::Sqrt(kxx*kxx + kyy*kyy);
 921     // Radius in projectile frame
 922     const Double_t kr2   = TMath::Sqrt(kr1*kr1 + kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
 923     const Double_t ky    = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
 924
 925     rw += ky * r;
 926     w  += ky;
 927     r  += kDr;
 928   } // radial steps
 929
 930   Double_t y=0.;
 931   if (!kiopt) { // My length definition (is exact for hard disk)
 932       if(w) y= 2. * rw / w;
 933   } else {
 934       const Double_t knorm=fgWSta->Eval(1e-4);
 935       if(knorm) y =  TMath::Sqrt(2. * rw * kDr / knorm / knorm);
 936   }
 937   return y; //fm
 938 }
 939
 940 Double_t AliFastGlauber::WPathLength(Double_t* x, Double_t* par)
 941 {
 942   //
 943   //  Path Length as a function of phi
 944   //  Interaction point from random distribution
 945   //  as a function of the phi-direction
 946   const Double_t kphi0   = x[0];
 947   const Double_t kb      = par[0];
 948   fgWAlmond->SetParameter(0, kb);
 949   const Int_t    kNpi    = Int_t (par[1]); //Number of interactions
 950   const Int_t    kiopt   = Int_t(par[2]);  //Path Length definition
 951
 952   //
 953   //  r-steps
 954   //
 955   const Int_t    kNp   = 100;
 956   const Double_t kDr  = fgBMax/Double_t(kNp);
 957   Double_t l = 0.;  //  Path length
 958   for (Int_t in = 0; in < kNpi; in ++) {
 959     Double_t rw = 0.;
 960     Double_t w  = 0.;
 961     // Interaction point
 962     Double_t x0, y0;
 963     fgWAlmond->GetRandom2(x0, y0);
 964     // Initial radius
 965     const Double_t kr0  = TMath::Sqrt(x0*x0 + y0*y0);
 966     const Int_t    knps = Int_t ((fgBMax - kr0)/kDr) - 1;
 967
 968     // Radial steps
 969     Double_t r  = 0.;
 970     for (Int_t i = 0; (i < knps ); i++) {
 971       // Transform into target frame
 972       const Double_t kxx   = x0 + r * TMath::Cos(kphi0) + kb / 2.;
 973       const Double_t kyy   = y0 + r * TMath::Sin(kphi0);
 974       const Double_t kphi  = TMath::ATan2(kyy, kxx);
 975       const Double_t kr1   = TMath::Sqrt(kxx*kxx + kyy*kyy);
 976       // Radius in projectile frame
 977       const Double_t kr2   = TMath::Sqrt(kr1*kr1 + kb*kb - 2.*kr1*kb*TMath::Cos(kphi));
 978       const Double_t ky    = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
 979
 980       rw += ky * r;
 981       w  += ky;
 982       r  += kDr;
 983     } // steps
 984     // Average over interactions
 985     if (!kiopt) {
 986       if(w) l += (2. * rw / w);
 987     } else {
 988       const Double_t knorm=fgWSta->Eval(1e-4);
 989       if(knorm) l+= 2. * rw * kDr / knorm / knorm;
 990     }
 991   } // interactions
 992   Double_t ret=0;
 993   if (!kiopt)
 994     ret= l / kNpi;
 995   else
 996     ret=TMath::Sqrt( l / kNpi);
 997   return ret; //fm
 998 }
 999
1000 Double_t AliFastGlauber::CrossSection(Double_t b1, Double_t b2) const
1001 {
1002   //
1003   // Return the geometrical cross-section integrated from b1 to b2
1004   //
1005   return fgWSgeo->Integral(b1, b2)*10.; //mbarn
1006 }
1007
1008 Double_t AliFastGlauber::HardCrossSection(Double_t b1, Double_t b2) const
1009 {
1010   //
1011   // Return the hard cross-section integrated from b1 to b2
1012   //
1013   return fgWSbinary->Integral(b1, b2)*10.; //mbarn
1014 }
1015
1016 Double_t AliFastGlauber::FractionOfHardCrossSection(Double_t b1, Double_t b2) const
1017 {
1018   //
1019   // Return fraction of hard cross-section integrated from b1 to b2
1020   //
1021   return fgWSbinary->Integral(b1, b2)/fgWSbinary->Integral(0., 100.);
1022 }
1023
1024 Double_t AliFastGlauber::NHard(Double_t b1, Double_t b2) const
1025 {
1026   //
1027   //  Number of binary hard collisions
1028   //  as a function of b (nucl/ex/0302016 eq. 19)
1029   //
1030   const Double_t kshard=HardCrossSection(b1,b2);
1031   const Double_t ksgeo=CrossSection(b1,b2);
1032   if(ksgeo>0)
1033     return kshard/ksgeo;
1034   else return -1;
1035 }
1036
1037 Double_t AliFastGlauber::Binaries(Double_t b) const
1038 {
1039   //
1040   // Return number of binary hard collisions normalized to 1 at b=0
1041   //
1042   if(b==0) b=1e-4;
1043   return fgWSN->Eval(b)/fgWSN->Eval(1e-4);
1044 }
1045
1046 Double_t AliFastGlauber::MeanOverlap(Double_t b1, Double_t b2)
1047 {
1048 //
1049 // Calculate the mean overlap for impact parameter range b1 .. b2
1050 //
1051     Double_t sum  = 0.;
1052     Double_t sumc = 0.;
1053     Double_t b    = b1;
1054
1055     while (b < b2-0.005) {
1056         Double_t  nc = GetNumberOfCollisions(b);
1057         sum  += 10. * fgWStaa->Eval(b) * fgWSgeo->Eval(b) * 0.01 / (1. - TMath::Exp(-nc));
1058         sumc += 10. * fgWSgeo->Eval(b) * 0.01;
1059         b += 0.01;
1060     }
1061     return (sum / CrossSection(b1, b2));
1062 }
1063
1064
1065 Double_t AliFastGlauber::MeanNumberOfCollisionsPerEvent(Double_t b1, Double_t b2)
1066 {
1067 //
1068 // Calculate the mean number of collisions per event for impact parameter range b1 .. b2
1069 //
1070     Double_t sum  = 0.;
1071     Double_t sumc = 0.;
1072     Double_t b    = b1;
1073
1074     while (b < b2-0.005) {
1075         Double_t  nc = GetNumberOfCollisions(b);
1076         sum  += nc / (1. - TMath::Exp(-nc)) * 10. * fgWSgeo->Eval(b) * 0.01;
1077         sumc += 10. * fgWSgeo->Eval(b) * 0.01;
1078         b += 0.01;
1079     }
1080     return (sum / CrossSection(b1, b2));
1081 }
1082
1083
1084 Double_t AliFastGlauber::GetNumberOfBinaries(Double_t b) const
1085 {
1086   //
1087   // Return number of binary hard collisions at b
1088   //
1089   if(b==0) b=1e-4;
1090   return fgWSN->Eval(b);
1091 }
1092
1093 Double_t AliFastGlauber::Participants(Double_t  b) const
1094 {
1095   //
1096   // Return the number of participants normalized to 1 at b=0
1097   //
1098   if(b==0) b=1e-4;
1099   return (fgWParticipants->Eval(b)/fgWParticipants->Eval(1e-4));
1100 }
1101
1102 Double_t AliFastGlauber::GetNumberOfParticipants(Double_t  b) const
1103 {
1104   //
1105   // Return the number of participants for impact parameter b
1106   //
1107   if(b==0) b=1e-4;
1108   return (fgWParticipants->Eval(b));
1109 }
1110
1111 Double_t AliFastGlauber::GetNumberOfCollisions(Double_t  b) const
1112 {
1113   //
1114   // Return the number of collisions for impact parameter b
1115   //
1116   if(b==0) b=1e-4;
1117   return (fgWStaa->Eval(b)*fSigmaNN);
1118 }
1119
1120 Double_t AliFastGlauber::GetNumberOfCollisionsPerEvent(Double_t  b) const
1121 {
1122   //
1123   // Return the number of collisions per event (at least one collision)
1124   // for impact parameter b
1125   //
1126     Double_t n = GetNumberOfCollisions(b);
1127     if (n > 0.) {
1128         return (n / (1. - TMath::Exp(- n)));
1129     } else {
1130         return (0.);
1131     }
1132 }
1133
1134 void AliFastGlauber::SimulateTrigger(Int_t n)
1135 {
1136   //
1137   //  Simulates Trigger
1138   //
1139   TH1F* mbtH = new TH1F("mbtH", "MB Trigger b-Distribution",   100, 0., 20.);
1140   TH1F* hdtH = new TH1F("hdtH", "Hard Trigger b-Distribution", 100, 0., 20.);
1141   TH1F* mbmH = new TH1F("mbmH", "MB Trigger Multiplicity Distribution",   100, 0., 8000.);
1142   TH1F* hdmH = new TH1F("hdmH", "Hard Trigger Multiplicity Distribution", 100, 0., 8000.);
1143
1144   mbtH->SetXTitle("b [fm]");
1145   hdtH->SetXTitle("b [fm]");
1146   mbmH->SetXTitle("Multiplicity");
1147   hdmH->SetXTitle("Multiplicity");
1148
1149   TCanvas *c0 = new TCanvas("c0","Trigger Simulation",400,10,600,700);
1150   c0->Divide(2,1);
1151   TCanvas *c1 = new TCanvas("c1","Trigger Simulation",400,10,600,700);
1152   c1->Divide(1,2);
1153
1154   //
1155   //
1156   Init(1);
1157   for (Int_t iev = 0; iev < n; iev++)
1158     {
1159       Float_t b, p, mult;
1160       GetRandom(b, p, mult);
1161       mbtH->Fill(b,1.);
1162       hdtH->Fill(b, p);
1163       mbmH->Fill(mult, 1.);
1164       hdmH->Fill(mult, p);
1165
1166       c0->cd(1);
1167       mbtH->Draw();
1168       c0->cd(2);
1169       hdtH->Draw();
1170       c0->Update();
1171
1172       c1->cd(1);
1173       mbmH->Draw();
1174       c1->cd(2);
1175       hdmH->Draw();
1176       c1->Update();
1177     }
1178 }
1179
1180 void AliFastGlauber::GetRandom(Float_t& b, Float_t& p, Float_t& mult)
1181 {
1182   //
1183   // Gives back a random impact parameter, hard trigger probability and multiplicity
1184   //
1185   b = fgWSgeo->GetRandom();
1186   const Float_t kmu = fgWSN->Eval(b);
1187   p = 1.-TMath::Exp(-kmu);
1188   mult = 6000./fgWSN->Eval(1.) * kmu;
1189 }
1190
1191 void AliFastGlauber::GetRandom(Int_t& bin, Bool_t& hard)
1192 {
1193   //
1194   // Gives back a random impact parameter bin, and hard trigger decission
1195   //
1196   const Float_t kb  = fgWSgeo->GetRandom();
1197   const Float_t kmu = fgWSN->Eval(kb) * fSigmaHard;
1198   const Float_t kp  = 1.-TMath::Exp(-kmu);
1199   if (kb < 5.) {
1200     bin = 1;
1201   } else if (kb <  8.6) {
1202     bin = 2;
1203   } else if (kb < 11.2) {
1204     bin = 3;
1205   } else if (kb < 13.2) {
1206     bin = 4;
1207   } else if (kb < 15.0) {
1208     bin = 5;
1209   } else {
1210     bin = 6;
1211   }
1212   hard = kFALSE;
1213   const Float_t kr = gRandom->Rndm();
1214   if (kr < kp) hard = kTRUE;
1215 }
1216
1217 Double_t  AliFastGlauber::GetRandomImpactParameter(Double_t bmin, Double_t bmax)
1218 {
1219   //
1220   // Gives back a random impact parameter in the range bmin .. bmax
1221   //
1222   Float_t b = -1.;
1223   while(b < bmin || b > bmax)
1224     b = fgWSgeo->GetRandom();
1225   return b;
1226 }
1227
1228 void AliFastGlauber::StoreFunctions() const
1229 {
1230   //
1231   // Store in file functions
1232   //
1233   TFile* ff = new TFile(fName.Data(),"recreate");
1234   fgWStaa->Write("WStaa");
1235   fgWParticipants->Write("WParticipants");
1236   ff->Close();
1237   return;
1238 }
1239
1240 //=================== Added by A. Dainese 11/02/04 ===========================
1241
1242 void AliFastGlauber::StoreAlmonds() const
1243 {
1244   //
1245   // Store in file
1246   // 40 almonds for b = (0.25+k*0.5) fm (k=0->39)
1247   //
1248   Char_t almondName[100];
1249   TFile* ff = new TFile(fName.Data(),"update");
1250   for(Int_t k=0; k<40; k++) {
1251     sprintf(almondName,"WAlmondFixedB%d",k);
1252     Double_t b = 0.25+k*0.5;
1253     Info("StoreAlmonds"," b = %f\n",b);
1254     fgWAlmond->SetParameter(0,b);
1255     fgWAlmond->Write(almondName);
1256   }
1257   ff->Close();
1258   return;
1259 }
1260
1261 void AliFastGlauber::SetCentralityClass(Double_t xsecFrLow,Double_t xsecFrUp)
1262 {
1263   //
1264   // Set limits of centrality class as fractions
1265   // of the geomtrical cross section
1266   //
1267   if(xsecFrLow>1. || xsecFrUp>1. || xsecFrLow>xsecFrUp) {
1268     Error("SetCentralityClass", "Please set 0 <= xsecFrLow <= xsecFrUp <= 1\n");
1269     return;
1270   }
1271
1272   Double_t bLow=0.,bUp=0.;
1273   Double_t xsecFr=0.;
1274   const Double_t knorm=fgWSgeo->Integral(0.,100.);
1275   while(xsecFr<xsecFrLow) {
1276     xsecFr = fgWSgeo->Integral(0.,bLow)/knorm;
1277     bLow += 0.1;
1278   }
1279   bUp = bLow;
1280   while(xsecFr<xsecFrUp) {
1281     xsecFr = fgWSgeo->Integral(0.,bUp)/knorm;
1282     bUp += 0.1;
1283   }
1284
1285   Info("SetCentralityClass", "Centrality class: %4.2f-%4.2f; %4.1f < b < %4.1f fm",
1286          xsecFrLow,xsecFrUp,bLow,bUp);
1287   fgWSbinary->SetRange(bLow,bUp);
1288   fBmin=bLow;
1289   fBmax=bUp;
1290   return;
1291 }
1292
1293 void AliFastGlauber::GetRandomBHard(Double_t& b)
1294 {
1295   //
1296   // Get random impact parameter according to distribution of
1297   // hard (binary) cross-section, in the range defined by the centrality class
1298   //
1299   b = fgWSbinary->GetRandom();
1300   Int_t bin = 2*(Int_t)b;
1301   if( (b-(Int_t)b) > 0.5) bin++;
1302   fgWAlmondCurrent = fgWAlmondFixedB[bin];
1303   return;
1304 }
1305
1306 void AliFastGlauber::GetRandomXY(Double_t& x,Double_t& y)
1307 {
1308   //
1309   // Get random position of parton production point according to
1310   // product of thickness functions
1311   //
1312   fgWAlmondCurrent->GetRandom2(x,y);
1313   return;
1314 }
1315
1316 void AliFastGlauber::GetRandomPhi(Double_t& phi)
1317 {
1318   //
1319   // Get random parton azimuthal propagation direction
1320   //
1321   phi = 2.*TMath::Pi()*gRandom->Rndm();
1322   return;
1323 }
1324
1325 Double_t AliFastGlauber::CalculateLength(Double_t b,Double_t x0,Double_t y0,Double_t phi0)
1326 {
1327   //
1328   // Calculate path length for a parton with production point (x0,y0)
1329   // and propagation direction (ux=cos(phi0),uy=sin(phi0))
1330   // in a collision with impact parameter b
1331   //
1332
1333   // number of steps in l
1334   const Int_t    kNp  = 100;
1335   const Double_t kDl  = fgBMax/Double_t(kNp);
1336
1337   if(fEllDef==1) {
1338     //
1339     // Definition 1:
1340     //
1341     // ell = 2 * \int_0^\infty dl*l*(T_A*T_B)(x0+l*ux,y0+l*uy) /
1342     //           \int_0^\infty dl*(T_A*T_B)(x0+l*ux,y0+l*uy)
1343     //
1344
1345     // Initial radius
1346     const Double_t kr0 = TMath::Sqrt(x0*x0 + y0*y0);
1347     const Int_t knps = Int_t ((fgBMax - kr0)/kDl) - 1;
1348     Double_t l  = 0.;
1349     Double_t integral1 = 0.;
1350     Double_t integral2 = 0.;
1351     // Radial steps
1352     for (Int_t i = 0; i < knps; i++) {
1353
1354       // Transform into target frame
1355       const Double_t kxx   = x0 + l * TMath::Cos(phi0) + b / 2.;
1356       const Double_t kyy   = y0 + l * TMath::Sin(phi0);
1357       const Double_t kphi  = TMath::ATan2(kyy, kxx);
1358       const Double_t kr1   = TMath::Sqrt(kxx*kxx + kyy*kyy);
1359       // Radius in projectile frame
1360       const Double_t kr2   = TMath::Sqrt(kr1*kr1 + b*b - 2.*kr1*b*TMath::Cos(kphi));
1361       const Double_t kprodTATB = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
1362
1363       integral1 += kprodTATB * l * kDl;
1364       integral2 += kprodTATB * kDl;
1365       l  += kDl;
1366     } // steps
1367
1368     Double_t ell=0.;
1369     if(integral2)
1370       ell = (2. * integral1 / integral2);
1371     return ell;
1372   } else if(fEllDef==2) {
1373     //
1374     // Definition 2:
1375     //
1376     // ell = \int_0^\infty dl*
1377     //       \Theta((T_A*T_B)(x0+l*ux,y0+l*uy)-0.5*(T_A*T_B)(0,0))
1378     //
1379
1380     // Initial radius
1381     const Double_t kr0  = TMath::Sqrt(x0*x0 + y0*y0);
1382     const Int_t knps = Int_t ((fgBMax - kr0)/kDl) - 1;
1383     const Double_t kprodTATBHalfMax = 0.5*fgWAlmondCurrent->Eval(0.,0.);
1384     // Radial steps
1385     Double_t l  = 0.;
1386     Double_t integral = 0.;
1387     for (Int_t i = 0; i < knps; i++) {
1388       // Transform into target frame
1389       const Double_t kxx   = x0 + l * TMath::Cos(phi0) + b / 2.;
1390       const Double_t kyy   = y0 + l * TMath::Sin(phi0);
1391       const Double_t kphi  = TMath::ATan2(kyy, kxx);
1392       const Double_t kr1   = TMath::Sqrt(kxx*kxx + kyy*kyy);
1393       // Radius in projectile frame
1394       const Double_t kr2   = TMath::Sqrt(kr1*kr1 + b*b - 2.*kr1*b*TMath::Cos(kphi));
1395       const Double_t kprodTATB = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
1396       if(kprodTATB>kprodTATBHalfMax) integral += kDl;
1397       l  += kDl;
1398     } // steps
1399     Double_t ell = integral;
1400     return ell;
1401   } else {
1402     Error("CalculateLength","Wrong length definition setting: %d !\n",fEllDef);
1403     return -1.;
1404   }
1405 }
1406
1407 void AliFastGlauber::GetLengthAndPhi(Double_t& ell,Double_t& phi,Double_t b)
1408 {
1409   //
1410   // Return length from random b, x0, y0, phi0
1411   // Return also phi0
1412   //
1413   Double_t x0,y0,phi0;
1414   if(b<0.) GetRandomBHard(b);
1415   GetRandomXY(x0,y0);
1416   GetRandomPhi(phi0);
1417   phi = phi0;
1418   ell = CalculateLength(b,x0,y0,phi0);
1419   return;
1420 }
1421
1422 void AliFastGlauber::GetLength(Double_t& ell,Double_t b)
1423 {
1424   //
1425   // Return length from random b, x0, y0, phi0
1426   //
1427   Double_t phi;
1428   GetLengthAndPhi(ell,phi,b);
1429   return;
1430 }
1431
1432 void AliFastGlauber::GetLengthsBackToBackAndPhi(Double_t& ell1,Double_t& ell2,Double_t &phi,Double_t b)
1433 {
1434   //
1435   // Return 2 lengths back to back from random b, x0, y0, phi0
1436   // Return also phi0
1437  //
1438   Double_t x0,y0,phi0;
1439   if(b<0.) GetRandomBHard(b);
1440   GetRandomXY(x0,y0);
1441   GetRandomPhi(phi0);
1442   const Double_t kphi0plusPi = phi0+TMath::Pi();
1443   phi = phi0;
1444   ell1 = CalculateLength(b,x0,y0,phi0);
1445   ell2 = CalculateLength(b,x0,y0,kphi0plusPi);
1446   return;
1447 }
1448
1449 void AliFastGlauber::GetLengthsBackToBack(Double_t& ell1,Double_t& ell2,
1450                                           Double_t b)
1451 {
1452   //
1453   // Return 2 lengths back to back from random b, x0, y0, phi0
1454   //
1455   Double_t phi;
1456   GetLengthsBackToBackAndPhi(ell1,ell2,phi,b);
1457   return;
1458 }
1459
1460 void AliFastGlauber::GetLengthsForPythia(Int_t n,Double_t* phi,Double_t* ell, Double_t b)
1461 {
1462   //
1463   // Returns lenghts for n partons with azimuthal angles phi[n]
1464   // from random b, x0, y0
1465   //
1466   Double_t x0, y0;
1467   if(b < 0.) GetRandomBHard(b);
1468   GetRandomXY(x0,y0);
1469   for(Int_t i = 0; i< n; i++) ell[i] = CalculateLength(b,x0,y0,phi[i]);
1470   return;
1471 }
1472
1473 void AliFastGlauber::PlotBDistr(Int_t n)
1474 {
1475   //
1476   // Plot distribution of n impact parameters
1477   //
1478   Double_t b;
1479   TH1F *hB = new TH1F("hB","dN/db",100,0,fgBMax);
1480   hB->SetXTitle("b [fm]");
1481   hB->SetYTitle("dN/db [a.u.]");
1482   hB->SetFillColor(3);
1483   for(Int_t i=0; i<n; i++) {
1484     GetRandomBHard(b);
1485     hB->Fill(b);
1486   }
1487   TCanvas *cB = new TCanvas("cB","Impact parameter distribution",0,0,500,500);
1488   cB->cd();
1489   hB->Draw();
1490   return;
1491 }
1492
1493 void AliFastGlauber::PlotLengthDistr(Int_t n,Bool_t save,const char *fname)
1494 {
1495   //
1496   // Plot length distribution
1497   //
1498   Double_t ell;
1499   TH1F *hEll = new TH1F("hEll","Length distribution",64,-0.5,15);
1500   hEll->SetXTitle("Transverse path length, L [fm]");
1501   hEll->SetYTitle("Probability");
1502   hEll->SetFillColor(2);
1503   for(Int_t i=0; i<n; i++) {
1504     GetLength(ell);
1505     hEll->Fill(ell);
1506   }
1507   hEll->Scale(1/(Double_t)n);
1508   TCanvas *cL = new TCanvas("cL","Length distribution",0,0,500,500);
1509   cL->cd();
1510   hEll->Draw();
1511
1512   if(save) {
1513     TFile *f = new TFile(fname,"recreate");
1514     hEll->Write();
1515     f->Close();
1516   }
1517   return;
1518 }
1519
1520 void AliFastGlauber::PlotLengthB2BDistr(Int_t n,Bool_t save,const char *fname)
1521 {
1522   //
1523   // Plot lengths back-to-back distributions
1524   //
1525   Double_t ell1,ell2;
1526   TH2F *hElls = new TH2F("hElls","Lengths back-to-back",100,0,15,100,0,15);
1527   hElls->SetXTitle("Transverse path length, L [fm]");
1528   hElls->SetYTitle("Transverse path length, L [fm]");
1529   for(Int_t i=0; i<n; i++) {
1530     GetLengthsBackToBack(ell1,ell2);
1531     hElls->Fill(ell1,ell2);
1532   }
1533   hElls->Scale(1/(Double_t)n);
1534   TCanvas *cLs = new TCanvas("cLs","Length back-to-back distribution",0,0,500,500);
1535   gStyle->SetPalette(1,0);
1536   cLs->cd();
1537   hElls->Draw("col,Z");
1538   if(save) {
1539     TFile *f = new TFile(fname,"recreate");
1540     hElls->Write();
1541     f->Close();
1542   }
1543   return;
1544 }
1545
1546 void AliFastGlauber::PlotAlmonds() const
1547 {
1548   //
1549   // Plot almonds for some impact parameters
1550   //
1551   TCanvas *c = new TCanvas("c","Almonds",0,0,500,500);
1552   gStyle->SetPalette(1,0);
1553   c->Divide(2,2);
1554   c->cd(1);
1555   fgWAlmondFixedB[0]->Draw("cont1");
1556   c->cd(2);
1557   fgWAlmondFixedB[10]->Draw("cont1");
1558   c->cd(3);
1559   fgWAlmondFixedB[20]->Draw("cont1");
1560   c->cd(4);
1561   fgWAlmondFixedB[30]->Draw("cont1");
1562   return;
1563 }
1564
1565 //=================== Added by A. Dainese 05/03/04 ===========================
1566
1567 void AliFastGlauber::CalculateI0I1(Double_t& integral0,Double_t& integral1,
1568                                    Double_t b,Double_t x0,Double_t y0,
1569                                    Double_t phi0,Double_t ellCut) const
1570 {
1571   //
1572   // Calculate integrals:
1573   //  integral0 = \int_0^ellCut dl*(T_A*T_B)(x0+l*ux,y0+l*uy)
1574   //  integral1 = \int_0^ellCut dl*l*(T_A*T_B)(x0+l*ux,y0+l*uy)
1575   //
1576   // for a parton with production point (x0,y0)
1577   // and propagation direction (ux=cos(phi0),uy=sin(phi0))
1578   // in a collision with impact parameter b
1579   //
1580
1581   // number of steps in l
1582   const Int_t    kNp  = 100;
1583   const Double_t kDl  = fgBMax/Double_t(kNp);
1584
1585   // Initial radius
1586   const Double_t kr0  = TMath::Sqrt(x0 * x0 + y0 * y0);
1587   const Int_t knps = Int_t ((fgBMax - kr0)/kDl) - 1;
1588
1589   // Radial steps
1590   Double_t l  = 0.;
1591   integral0 = 0.;
1592   integral1 = 0.;
1593   Int_t i = 0;
1594   while((i < knps) && (l < ellCut)) {
1595     // Transform into target frame
1596     const Double_t kxx   = x0 + l * TMath::Cos(phi0) + b / 2.;
1597     const Double_t kyy   = y0 + l * TMath::Sin(phi0);
1598     const Double_t kphi  = TMath::ATan2(kyy, kxx);
1599     const Double_t kr1   = TMath::Sqrt(kxx*kxx + kyy*kyy);
1600     // Radius in projectile frame
1601     const Double_t kr2   = TMath::Sqrt(kr1*kr1 + b*b - 2.*kr1*b*TMath::Cos(kphi));
1602     const Double_t kprodTATB = fgWSta->Eval(kr1) * fgWSta->Eval(kr2);
1603     integral0 += kprodTATB * kDl;
1604     integral1 += kprodTATB * l * kDl;
1605     l  += kDl;
1606     i++;
1607   } // steps
1608   return;
1609 }
1610
1611 void AliFastGlauber::GetI0I1AndPhi(Double_t& integral0,Double_t& integral1,
1612                                    Double_t& phi,
1613                                    Double_t ellCut,Double_t b)
1614 {
1615   //
1616   // Return I0 and I1 from random b, x0, y0, phi0
1617   // Return also phi
1618   //
1619   Double_t x0,y0,phi0;
1620   if(b<0.) GetRandomBHard(b);
1621   GetRandomXY(x0,y0);
1622   GetRandomPhi(phi0);
1623   phi = phi0;
1624   CalculateI0I1(integral0,integral1,b,x0,y0,phi0,ellCut);
1625   return;
1626 }
1627
1628 void AliFastGlauber::GetI0I1(Double_t& integral0,Double_t& integral1,
1629                              Double_t ellCut,Double_t b)
1630 {
1631   //
1632   // Return I0 and I1 from random b, x0, y0, phi0
1633   //
1634   Double_t phi;
1635   GetI0I1AndPhi(integral0,integral1,phi,ellCut,b);
1636   return;
1637 }
1638
1639 void AliFastGlauber::GetI0I1BackToBackAndPhi(Double_t& integral01,Double_t& integral11,
1640                                              Double_t& integral02,Double_t& integral12,
1641                                              Double_t& phi,
1642                                              Double_t ellCut,Double_t b)
1643 {
1644   //
1645   // Return 2 pairs of I0 and I1 back to back from random b, x0, y0, phi0
1646   // Return also phi0
1647   //
1648   Double_t x0,y0,phi0;
1649   if(b<0.) GetRandomBHard(b);
1650   GetRandomXY(x0,y0);
1651   GetRandomPhi(phi0);
1652   phi = phi0;
1653   const Double_t kphi0plusPi = phi0+TMath::Pi();
1654   CalculateI0I1(integral01,integral11,b,x0,y0,phi0,ellCut);
1655   CalculateI0I1(integral02,integral12,b,x0,y0,kphi0plusPi,ellCut);
1656   return;
1657 }
1658
1659 void AliFastGlauber::GetI0I1BackToBackAndPhiAndXY(Double_t& integral01,Double_t& integral11,
1660                                                   Double_t& integral02,Double_t& integral12,
1661                                                   Double_t& phi,Double_t &x,Double_t &y,
1662                                                   Double_t ellCut,Double_t b)
1663 {
1664   //
1665   // Return 2 pairs of I0 and I1 back to back from random b, x0, y0, phi0
1666   // Return also phi0
1667   //
1668   Double_t x0,y0,phi0;
1669   if(b<0.) GetRandomBHard(b);
1670   GetRandomXY(x0,y0);
1671   GetRandomPhi(phi0);
1672   phi = phi0; x=x0; y=y0;
1673   const Double_t kphi0plusPi = phi0+TMath::Pi();
1674   CalculateI0I1(integral01,integral11,b,x0,y0,phi0,ellCut);
1675   CalculateI0I1(integral02,integral12,b,x0,y0,kphi0plusPi,ellCut);
1676   return;
1677 }
1678
1679 void AliFastGlauber::GetI0I1BackToBack(Double_t& integral01,Double_t& integral11,
1680                                        Double_t& integral02,Double_t& integral12,
1681                                        Double_t ellCut,Double_t b)
1682 {
1683   //
1684   // Return 2 pairs of I0 and I1 back to back from random b, x0, y0, phi0
1685   //
1686   Double_t phi;
1687   GetI0I1BackToBackAndPhi(integral01,integral11,integral02,integral12,
1688                           phi,ellCut,b);
1689   return;
1690 }
1691
1692 void AliFastGlauber::GetI0I1ForPythia(Int_t n,Double_t* phi,
1693                                       Double_t* integral0,Double_t* integral1,
1694                                       Double_t ellCut,Double_t b)
1695 {
1696   //
1697   // Returns I0 and I1 pairs for n partons with azimuthal angles phi[n]
1698   // from random b, x0, y0
1699   //
1700   Double_t x0,y0;
1701   if(b<0.) GetRandomBHard(b);
1702   GetRandomXY(x0,y0);
1703   for(Int_t i=0; i<n; i++)
1704     CalculateI0I1(integral0[i],integral1[i],b,x0,y0,phi[i],ellCut);
1705   return;
1706 }
1707
1708 void AliFastGlauber::GetI0I1ForPythiaAndXY(Int_t n,Double_t* phi,
1709                                       Double_t* integral0,Double_t* integral1,
1710                                       Double_t &x,Double_t& y,
1711                                       Double_t ellCut,Double_t b)
1712 {
1713   //
1714   // Returns I0 and I1 pairs for n partons with azimuthal angles phi[n]
1715   // from random b, x0, y0 and return x0,y0
1716   //
1717   Double_t x0,y0;
1718   if(b<0.) GetRandomBHard(b);
1719   GetRandomXY(x0,y0);
1720   for(Int_t i=0; i<n; i++)
1721     CalculateI0I1(integral0[i],integral1[i],b,x0,y0,phi[i],ellCut);
1722   x=x0;
1723   y=y0;
1724   return;
1725 }
1726
1727 void AliFastGlauber::PlotI0I1Distr(Int_t n,Double_t ellCut,
1728                                    Bool_t save,const char *fname)
1729 {
1730   //
1731   // Plot I0-I1 distribution
1732   //
1733   Double_t i0,i1;
1734   TH2F *hI0I1s = new TH2F("hI0I1s","I_{0} versus I_{1}",1000,0,0.001,1000,0,0.01);
1735   hI0I1s->SetXTitle("I_{0} [fm^{-3}]");
1736   hI0I1s->SetYTitle("I_{1} [fm^{-2}]");
1737
1738   TH1F *hI0 = new TH1F("hI0","I_{0} = #hat{q}L / k",
1739                        1000,0,0.001);
1740   hI0->SetXTitle("I_{0} [fm^{-3}]");
1741   hI0->SetYTitle("Probability");
1742   hI0->SetFillColor(3);
1743   TH1F *hI1 = new TH1F("hI1","I_{1} = #omega_{c} / k",
1744                        1000,0,0.01);
1745   hI1->SetXTitle("I_{1} [fm^{-2}]");
1746   hI1->SetYTitle("Probability");
1747   hI1->SetFillColor(4);
1748   TH1F *h2 = new TH1F("h2","2 I_{1}^{2}/I_{0} = R / k",
1749                       100,0,0.02);
1750   h2->SetXTitle("2 I_{1}^{2}/I_{0} [fm^{-1}]");
1751   h2->SetYTitle("Probability");
1752   h2->SetFillColor(2);
1753   TH1F *h3 = new TH1F("h3","2 I_{1}/I_{0} = L",
1754                       100,0,15);
1755   h3->SetXTitle("2 I_{1}/I_{0} [fm]");
1756   h3->SetYTitle("Probability");
1757   h3->SetFillColor(5);
1758   TH1F *h4 = new TH1F("h4","I_{0}^{2}/(2 I_{1}) = #hat{q} / k",
1759                       100,0,0.00015);
1760   h4->SetXTitle("I_{0}^{2}/(2 I_{1}) [fm^{-4}]");
1761   h4->SetYTitle("Probability");
1762   h4->SetFillColor(7);
1763
1764   for(Int_t i=0; i<n; i++) {
1765     GetI0I1(i0,i1,ellCut);
1766     hI0I1s->Fill(i0,i1);
1767     hI0->Fill(i0);
1768     hI1->Fill(i1);
1769     h2->Fill(2.*i1*i1/i0);
1770     h3->Fill(2.*i1/i0);
1771     h4->Fill(i0*i0/2./i1);
1772   }
1773   hI0->Scale(1/(Double_t)n);
1774   hI1->Scale(1/(Double_t)n);
1775   h2->Scale(1/(Double_t)n);
1776   h3->Scale(1/(Double_t)n);
1777   h4->Scale(1/(Double_t)n);
1778   hI0I1s->Scale(1/(Double_t)n);
1779
1780   TCanvas *cI0I1 = new TCanvas("cI0I1","I0 and I1",0,0,900,700);
1781   cI0I1->Divide(3,2);
1782   cI0I1->cd(1);
1783   hI0->Draw();
1784   cI0I1->cd(2);
1785   hI1->Draw();
1786   cI0I1->cd(3);
1787   h2->Draw();
1788   cI0I1->cd(4);
1789   h3->Draw();
1790   cI0I1->cd(5);
1791   h4->Draw();
1792   cI0I1->cd(6);
1793   gStyle->SetPalette(1,0);
1794   hI0I1s->Draw("col,Z");
1795
1796   if(save) {
1797     TFile *f = new TFile(fname,"recreate");
1798     hI0I1s->Write();
1799     hI0->Write();
1800     hI1->Write();
1801     h2->Write();
1802     h3->Write();
1803     h4->Write();
1804     f->Close();
1805   }
1806   return;
1807 }
1808
1809 void AliFastGlauber::PlotI0I1B2BDistr(Int_t n,Double_t ellCut,
1810                                       Bool_t save,const char *fname)
1811 {
1812   //
1813   // Plot I0-I1 back-to-back distributions
1814   //
1815   Double_t i01,i11,i02,i12;
1816   TH2F *hI0s = new TH2F("hI0s","I_{0}'s back-to-back",100,0,100,100,0,100);
1817   hI0s->SetXTitle("I_{0} [fm^{-3}]");
1818   hI0s->SetYTitle("I_{0} [fm^{-3}]");
1819   TH2F *hI1s = new TH2F("hI1s","I_{1}'s back-to-back",100,0,100,100,0,100);
1820   hI1s->SetXTitle("I_{1} [fm^{-2}]");
1821   hI1s->SetYTitle("I_{1} [fm^{-2}]");
1822
1823   for(Int_t i=0; i<n; i++) {
1824     GetI0I1BackToBack(i01,i11,i02,i12,ellCut);
1825     hI0s->Fill(i01,i02);
1826     hI1s->Fill(i11,i12);
1827   }
1828   hI0s->Scale(1/(Double_t)n);
1829   hI1s->Scale(1/(Double_t)n);
1830
1831   TCanvas *cI0I1s = new TCanvas("cI0I1s","I0 and I1 back-to-back distributions",0,0,800,400);
1832   gStyle->SetPalette(1,0);
1833   cI0I1s->Divide(2,1);
1834   cI0I1s->cd(1);
1835   hI0s->Draw("col,Z");
1836   cI0I1s->cd(2);
1837   hI1s->Draw("col,Z");
1838
1839   if(save) {
1840     TFile *f = new TFile(fname,"recreate");
1841     hI0s->Write();
1842     hI1s->Write();
1843     f->Close();
1844   }
1845   return;
1846 }
1847
1848 AliFastGlauber& AliFastGlauber::operator=(const  AliFastGlauber& rhs)
1849 {
1850 // Assignment operator
1851     rhs.Copy(*this);
1852     return *this;
1853 }
1854
1855 void AliFastGlauber::Copy(TObject&) const
1856 {
1857     //
1858     // Copy
1859     //
1860     Fatal("Copy","Not implemented!\n");
1861 }
1862