FASTSIM/AliFastGlauber.cxx

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