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