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