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