FASTSIM/AliFastGlauber.cxx

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