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