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