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