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