Averaging for pathlength corrected.
[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 // from AliRoot
19 #include "AliFastGlauber.h"
20 // from root
21 #include <TH1F.h>
22 #include <TF1.h>
23 #include <TF2.h>
24 #include <TCanvas.h>
25 #include <TRandom.h>
26 #include <TFile.h>
27
28 ClassImp(AliFastGlauber)
29
30 TF1*    AliFastGlauber::fWSb            = NULL;     
31 TF2*    AliFastGlauber::fWSbz           = NULL;    
32 TF1*    AliFastGlauber::fWSz            = NULL;     
33 TF1*    AliFastGlauber::fWSta           = NULL;    
34 TF2*    AliFastGlauber::fWStarfi        = NULL; 
35 TF2*    AliFastGlauber::fWAlmond        = NULL; 
36 TF1*    AliFastGlauber::fWStaa          = NULL;   
37 TF1*    AliFastGlauber::fWSgeo          = NULL;   
38 TF1*    AliFastGlauber::fWSbinary       = NULL;   
39 TF1*    AliFastGlauber::fWSN            = NULL;   
40 TF1*    AliFastGlauber::fWPathLength0   = NULL;   
41 TF1*    AliFastGlauber::fWPathLength    = NULL;
42 TF1*    AliFastGlauber::fWEnergyDensity = NULL;   
43 TF1*    AliFastGlauber::fWIntRadius     = NULL;   
44 Float_t AliFastGlauber::fbMax           = 0.;
45
46 AliFastGlauber::AliFastGlauber()
47 {
48 //  Default Constructor
49 //
50 //  Defaults for Pb
51     SetWoodSaxonParameters(6.624, 0.549, 0.00, 7.69e-4);
52     SetHardCrossSection();
53     SetMaxImpact();
54 }
55
56 void AliFastGlauber::Init(Int_t mode)
57 {
58 // Initialisation
59 //
60 //  Wood-Saxon
61 //
62     fWSb = new TF1("WSb", WSb, 0, fbMax, 4);
63     fWSb->SetParameter(0, fWSr0);
64     fWSb->SetParameter(1, fWSd);
65     fWSb->SetParameter(2, fWSw);
66     fWSb->SetParameter(3, fWSn);
67
68     fWSbz = new TF2("WSbz", WSbz, 0, fbMax, 4);
69     fWSbz->SetParameter(0, fWSr0);
70     fWSbz->SetParameter(1, fWSd);
71     fWSbz->SetParameter(2, fWSw);
72     fWSbz->SetParameter(3, fWSn);
73
74     fWSz = new TF1("WSz", WSz, 0, fbMax, 5);
75     fWSz->SetParameter(0, fWSr0);
76     fWSz->SetParameter(1, fWSd);
77     fWSz->SetParameter(2, fWSw);
78     fWSz->SetParameter(3, fWSn);
79
80 //
81 //  Thickness
82 //
83     fWSta = new TF1("WSta", WSta, 0., fbMax, 0);
84     
85 //
86 //  Overlap Kernel
87 //
88     fWStarfi = new TF2("WStarfi", WStarfi, 0., fbMax, 0., TMath::Pi(), 1);
89     fWStarfi->SetParameter(0, 0.);     
90     fWStarfi->SetNpx(200);     
91     fWStarfi->SetNpy(20);     
92 //
93 //  Almond shaped interaction region
94 //
95     fWAlmond = new TF2("WAlmond", WAlmond, -fbMax, fbMax, -fbMax, fbMax, 1);
96     fWAlmond->SetParameter(0, 0.);     
97     fWAlmond->SetNpx(200);     
98     fWAlmond->SetNpy(200);    
99 //
100 //  Path Length as a function of Phi
101 //    
102     fWPathLength0 = new TF1("WPathLength0", WPathLength0, -TMath::Pi(), TMath::Pi(), 2);
103     fWPathLength0->SetParameter(0, 0.);
104 //  Pathlength definition     
105     fWPathLength0->SetParameter(1, 0.);     
106
107     fWPathLength = new TF1("WPathLength", WPathLength, -TMath::Pi(), TMath::Pi(), 3);
108 //  Impact Parameter
109     fWPathLength->SetParameter(0, 0.);    
110 //  Number of interactions used for average
111     fWPathLength->SetParameter(1, 1000.);    
112 //  Pathlength definition
113     fWPathLength->SetParameter(2, 0);    
114
115     fWIntRadius = new TF1("WIntRadius", WIntRadius, 0., fbMax, 1);
116     fWIntRadius->SetParameter(0, 0.);    
117
118
119 //
120 //  Overlap
121 //
122     if (! mode) {
123         fWStaa = new TF1("WStaa", WStaa, 0., fbMax, 0);
124         fWStaa->SetNpx(100);
125     } else {
126         TFile* f = new TFile("$(ALICE_ROOT)/FASTSIM/data/glauberPbPb.root");
127         fWStaa = (TF1*) f->Get("WStaa");
128     }
129     
130 //
131     fWEnergyDensity = new TF1("WEnergyDensity", WEnergyDensity, 0., 2. * fWSr0, 1);
132     fWEnergyDensity->SetParameter(0, fWSr0 + 1.);
133     
134 //
135 //  Geometrical Cross-Section
136 //
137     fWSgeo = new TF1("WSgeo", WSgeo, 0., fbMax, 0);
138     fWSgeo->SetNpx(100);
139 //
140 //  Hard cross section (~ binary collisions)
141 //
142     fWSbinary = new TF1("WSbinary", WSbinary, 0., fbMax, 1);
143     fWSbinary->SetParameter(0, fSigmaHard); // mb
144     fWSbinary->SetNpx(100);
145 //
146 // Hard collisions per event
147 //
148     fWSN = new TF1("WSN", WSN, 0., fbMax, 1);
149     fWSN->SetNpx(100);
150 }
151
152 void AliFastGlauber::DrawWSb()
153 {
154 //
155 //  Draw Wood-Saxon Nuclear Density Function
156 //
157     TCanvas *c1 = new TCanvas("c1","Wood Saxon",400,10,600,700);
158     c1->cd();
159     fWSb->Draw();
160 }
161
162 void AliFastGlauber::DrawOverlap()
163 {
164 //
165 //  Draw Overlap Function
166 //
167     TCanvas *c2 = new TCanvas("c2","Overlap",400,10,600,700);
168     c2->cd();
169     fWStaa->Draw();
170 }
171
172 void AliFastGlauber::DrawThickness()
173 {
174 //
175 //  Draw Thickness Function
176 //
177     TCanvas *c3 = new TCanvas("c3","Thickness",400,10,600,700);
178     c3->cd();
179     fWSta->Draw();
180 }
181
182 void AliFastGlauber::DrawGeo()
183 {
184 //
185 //  Draw Geometrical Cross-Section
186 //
187     TCanvas *c3 = new TCanvas("c3","Geometrical Cross-Section",400,10,600,700);
188     c3->cd();
189     fWSgeo->Draw();
190 }
191
192 void AliFastGlauber::DrawBinary()
193 {
194 //
195 //  Draw Wood-Saxon Nuclear Density Function
196 //
197     TCanvas *c4 = new TCanvas("c4","Binary Cross-Section",400,10,600,700);
198     c4->cd();
199     fWSbinary->Draw();
200 }
201
202 void AliFastGlauber::DrawN()
203 {
204 //
205 //  Draw Binaries per event
206 //
207     TCanvas *c5 = new TCanvas("c5","Binaries per event",400,10,600,700);
208     c5->cd();
209     fWSN->Draw();
210 }
211
212 void AliFastGlauber::DrawKernel(Double_t b)
213 {
214 //
215 //  Draw Kernel
216 //
217     TCanvas *c6 = new TCanvas("c6","Kernel",400,10,600,700);
218     c6->cd();
219     fWStarfi->SetParameter(0, b);
220     fWStarfi->Draw();
221 }
222
223 void AliFastGlauber::DrawAlmond(Double_t b)
224 {
225 //
226 //  Draw Interaction Almond
227 //
228     TCanvas *c7 = new TCanvas("c7","Almond",400,10,600,700);
229     c7->cd();
230     fWAlmond->SetParameter(0, b);
231     fWAlmond->Draw();
232 }
233
234 void AliFastGlauber::DrawPathLength0(Double_t b, Int_t iopt)
235 {
236 //
237 //  Draw Path Length
238 //
239     TCanvas *c8 = new TCanvas("c8","Path Length",400,10,600,700);
240     c8->cd();
241     fWPathLength0->SetParameter(0, b);
242     fWPathLength0->SetParameter(1, Double_t(iopt));
243     fWPathLength0->SetMinimum(0.); 
244     fWPathLength0->SetMaximum(10.); 
245     fWPathLength0->Draw();
246 }
247
248 void AliFastGlauber::DrawPathLength(Double_t b , Int_t ni, Int_t iopt)
249 {
250 //
251 //  Draw Path Length
252 //
253     TCanvas *c9 = new TCanvas("c9","Path Length",400,10,600,700);
254     c9->cd();
255     fWAlmond->SetParameter(0, b);
256
257     fWPathLength->SetParameter(0, b);
258     fWPathLength->SetParameter(1, Double_t (ni));
259     fWPathLength->SetParameter(2, Double_t (iopt));
260     fWPathLength->SetMinimum(0.); 
261     fWPathLength->SetMaximum(10.); 
262     fWPathLength->Draw();
263 }
264
265 void AliFastGlauber::DrawIntRadius(Double_t b)
266 {
267 //
268 //  Draw Interaction Radius
269 //
270     TCanvas *c10 = new TCanvas("c10","Interaction Radius",400,10,600,700);
271     c10->cd();
272     fWIntRadius->SetParameter(0, b);
273     fWIntRadius->SetMinimum(0.);
274     fWIntRadius->Draw();
275 }
276
277 void AliFastGlauber::DrawEnergyDensity()
278 {
279 //
280 //  Draw energy density
281 //
282     TCanvas *c11 = new TCanvas("c11","Energy Density",400, 10, 600, 700);
283     c11->cd();
284     fWEnergyDensity->SetMinimum(0.);
285     fWEnergyDensity->Draw();
286 }
287
288 Double_t AliFastGlauber::WSb(Double_t* x, Double_t* par)
289 {
290 //
291 //  Woods-Saxon Parameterisation
292 //  as a function of radius
293 //
294     Double_t xx  = x[0];
295     Double_t r0  = par[0];
296     Double_t d   = par[1];
297     Double_t w   = par[2];
298     Double_t n   = par[3];
299     
300     Double_t y  = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
301
302     return y;
303 }
304
305 Double_t AliFastGlauber::WSbz(Double_t* x, Double_t* par)
306 {
307 //
308 //  Wood Saxon Parameterisation
309 //  as a function of z and  b
310 //
311     Double_t bb  = x[0];
312     Double_t zz  = x[1];
313     Double_t r0  = par[0];
314     Double_t d   = par[1];
315     Double_t w   = par[2];
316     Double_t n   = par[3];
317     Double_t xx  = TMath::Sqrt(bb*bb+zz*zz);
318     Double_t y  = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
319
320     return y;
321 }
322
323 Double_t AliFastGlauber::WSz(Double_t* x, Double_t* par)
324 {
325 //
326 //  Wood Saxon Parameterisation
327 //  as a function of z for fixed b
328 //
329     Double_t bb  = par[4];
330     Double_t zz  = x[0];
331     Double_t r0  = par[0];
332     Double_t d   = par[1];
333     Double_t w   = par[2];
334     Double_t n   = par[3];
335     Double_t xx  = TMath::Sqrt(bb*bb+zz*zz);
336     Double_t y  = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
337
338     return y;
339 }
340
341 Double_t AliFastGlauber::WSta(Double_t* x, Double_t* /*par*/)
342 {
343 //
344 //  Thickness function 
345 //
346     Double_t b  = x[0];
347     fWSz->SetParameter(4, b);
348     Double_t y  = 2. * fWSz->Integral(0., fbMax);
349     return y;
350 }
351
352
353
354 Double_t AliFastGlauber::WStarfi(Double_t* x, Double_t* par)
355 {
356 //
357 //  Kernel for overlap function
358 //
359     Double_t b    = par[0];
360     Double_t r1   = x[0];
361     Double_t phi  = x[1];
362     Double_t r2   = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi)); 
363     Double_t y    = r1 * fWSta->Eval(r1) * fWSta->Eval(r2);
364     return y;
365 }
366
367
368 Double_t AliFastGlauber::WAlmond(Double_t* x, Double_t* par)
369 {
370 //
371 //  Almond shaped interaction region
372 //
373     Double_t b    = par[0];
374     Double_t xx   = x[0] + b/2.;
375     Double_t yy   = x[1];
376     Double_t r1   = TMath::Sqrt(xx * xx + yy * yy);
377     Double_t phi  = TMath::ATan2(yy,xx);
378     
379     Double_t r2   = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi)); 
380 //
381 //  Interaction probability calculated as product of thicknesses
382 //
383     Double_t y    = fWSta->Eval(r1) * fWSta->Eval(r2);
384     return y;
385 }
386
387 Double_t AliFastGlauber::WIntRadius(Double_t* x, Double_t* par)
388 {
389 //
390 //  Average radius at which interaction takes place
391 //
392 //  Radius in the Almond
393     Double_t r    = x[0];
394 //  Impact parameter
395     Double_t b    = par[0];
396     fWAlmond->SetParameter(0, b);
397 //  Steps in phi
398     Double_t dphi = 2. * TMath::Pi() / 100.;
399 //  Average over phi    
400     Double_t phi  = 0.;
401     Double_t y    = 0.;
402
403     for (Int_t i = 0; i < 100; i++) {
404         Double_t xx = r * TMath::Cos(phi);
405         Double_t yy = r * TMath::Sin(phi);
406         y   += fWAlmond->Eval(xx,yy);
407         phi += dphi;
408     } // phi loop
409 // Result multiplied by Jacobian (2 pi r)     
410     return (2. * TMath::Pi() * y * r / 100.);
411 }
412
413 Double_t AliFastGlauber::WPathLength0(Double_t* x, Double_t* par)
414 {
415 //
416 //  Path Length as a function of phi for interaction point fixed at (0,0)
417 //
418 //
419 //  Steps in r 
420     const Int_t    np  = 100;
421     const Double_t dr  = fbMax/Double_t(np);
422 //  Impact parameter    
423     Double_t b      = par[0];
424 //  Path Length definition
425     Int_t    iopt   = Int_t(par[1]);
426     
427 //  Phi direction in Almond
428     Double_t phi0   = x[0];
429     Double_t r  = 0.;
430     Double_t rw = 0.;
431     Double_t w  = 0.;
432 //  Step along radial direction phi   
433     for (Int_t i = 0; i < np; i++) {
434 //
435 //  Transform into target frame
436 //
437         Double_t xx   = r * TMath::Cos(phi0) + b / 2.;
438         Double_t yy   = r * TMath::Sin(phi0);
439         Double_t phi  = TMath::ATan2(yy, xx);
440         
441         Double_t r1   = TMath::Sqrt(xx * xx + yy * yy);
442 // Radius in projectile frame
443         Double_t r2   = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi)); 
444         Double_t y    = fWSta->Eval(r1) * fWSta->Eval(r2);
445
446         rw += y * r;
447         w  += y;
448         r  += dr;
449     } // radial steps
450 //
451 //  My length definition (is exact for hard disk)
452     if (!iopt) {
453         return (2. * rw / w);
454     } else {
455         return TMath::Sqrt(2. * rw * dr / fWSta->Eval(0.01) / fWSta->Eval(0.01));
456     }
457 }
458
459 Double_t AliFastGlauber::WPathLength(Double_t* x, Double_t* par)
460 {
461 //
462 //  Path Length as a function of phi 
463 //  Interaction point from random distribution
464 //
465 //
466 //  r-steps
467 // 
468     const Int_t    np   = 100;
469     const Double_t dr  = fbMax/Double_t(np);
470 //  Number of interactions
471     const Int_t    npi  = Int_t (par[1]);
472
473 //
474 //  Impact parameter    
475     Double_t b      = par[0];
476 //  Path Length definition 
477     Int_t    iopt   = Int_t(par[2]);
478 //  Phi direction
479     Double_t phi0   = x[0];
480
481     printf("phi0 %f \n", phi0);
482     
483 //  Path length 
484     Double_t l = 0.;
485     
486     for (Int_t in = 0; in < npi; in ++) {
487         Double_t rw = 0.;
488         Double_t w  = 0.;
489         
490         // Interaction point
491         Double_t x0, y0;
492         fWAlmond->GetRandom2(x0, y0);
493 // Initial radius
494         Double_t r0  = TMath::Sqrt(x0 * x0 + y0 * y0);
495         Int_t    nps = Int_t ((fbMax - r0)/dr) - 1;
496         
497         Double_t r  = 0.;
498 // Radial steps
499         for (Int_t i = 0; (i < nps ); i++) {
500             
501 // Transform into target frame
502             Double_t xx   = x0 + r * TMath::Cos(phi0) + b / 2.;
503             Double_t yy   = y0 + r * TMath::Sin(phi0);
504             Double_t phi  = TMath::ATan2(yy, xx);
505             Double_t r1   = TMath::Sqrt(xx * xx + yy * yy);
506 // Radius in projectile frame
507             Double_t r2   = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi)); 
508             Double_t y    = fWSta->Eval(r1) * fWSta->Eval(r2);
509             
510             rw += y * r;
511             w  += y;
512             r  += dr;
513         } // steps
514 // Average over interactions
515         if (!iopt) {
516             l += (2. * rw / w);
517         } else {
518             l+= 2. * rw * dr / fWSta->Eval(0.01) / fWSta->Eval(0.01);
519         }
520     } // interactions
521     if (!iopt) 
522         return (l / Double_t(npi));
523     else 
524         return (TMath::Sqrt(l / Double_t(npi)));
525 }
526
527 Double_t AliFastGlauber::WStaa(Double_t* x, Double_t* /*par*/)
528 {
529 //
530 //  Overlap function
531 //
532     Double_t b    = x[0];
533     fWStarfi->SetParameter(0, b);
534 /*
535     Double_t al[2];
536     Double_t bl[2];
537     al[0] = 0.;
538     al[1] = 0.;
539     bl[0] = 6.6;
540     bl[1] = TMath::Pi();
541     Double_t err;
542     
543     Double_t y =  2. * fWStarfi->IntegralMultiple(2, al, bl, 0.001, err);
544     printf("WStaa: %f %f %f\n", b, y, err);
545 */
546 //
547 //  MC Integration
548 //
549     Double_t y = 0;
550     for (Int_t i = 0; i < 100000; i++)
551     {
552         Double_t phi = TMath::Pi() * gRandom->Rndm();
553         Double_t b1  = fbMax       * gRandom->Rndm();   
554         y += fWStarfi->Eval(b1, phi);
555     }
556     y *= 2. * 0.1 *  208. * 208. * TMath::Pi() * fbMax / 100000.;
557     return y;
558 }
559
560 Double_t AliFastGlauber::WSgeo(Double_t* x, Double_t* /*par*/)
561 {
562 //
563 //  Geometrical Cross-Section
564 //
565     Double_t b    = x[0];
566     Double_t taa  = fWStaa->Eval(b);
567     const Double_t sigma = 55.6; // mbarn
568     
569     Double_t y    = 2. * TMath::Pi() * b * (1. - TMath::Exp(- sigma * taa)); // fm
570     return y;
571 }
572
573
574 Double_t AliFastGlauber::WSbinary(Double_t* x, Double_t* par)
575 {
576 //
577 //  Number of binary collisions
578 //
579     Double_t b     = x[0];
580     Double_t sigma = par[0];
581     Double_t taa   = fWStaa->Eval(b);
582     
583     Double_t y    = 2. * TMath::Pi() * b * sigma * taa; // fm
584     return y;
585 }
586
587 Double_t AliFastGlauber::WSN(Double_t* x, Double_t* /*par*/)
588 {
589 //
590 //  Number of hard processes per event
591 //
592     Double_t b     = x[0];
593     Double_t y     = fWSbinary->Eval(b)/fWSgeo->Eval(b);
594     return y;
595 }
596
597 Double_t AliFastGlauber::WEnergyDensity(Double_t* x, Double_t* par)
598 {
599 //
600 //  Initial energy density as a function of the impact parameter
601 //
602     Double_t b     = x[0];
603     Double_t rA    = par[0];
604 //
605 //  Attention: area of transverse reaction zone in hard-sphere approximation !     
606     Double_t saa   = (TMath::Pi() - 2. * TMath::ASin(b/ 2./ rA)) * rA * rA 
607         - b * TMath::Sqrt(rA * rA - b * b/ 4.);
608     Double_t taa   = fWStaa->Eval(b);
609     
610     return (taa/saa);
611 }
612
613 void AliFastGlauber::SimulateTrigger(Int_t n)
614 {
615     //
616     //  Simulates Trigger
617     //
618     TH1F* mbtH = new TH1F("mbtH", "MB Trigger b-Distribution",   100, 0., 20.);
619     TH1F* hdtH = new TH1F("hdtH", "Hard Trigger b-Distribution", 100, 0., 20.);   
620     TH1F* mbmH = new TH1F("mbmH", "MB Trigger Multiplicity Distribution",   100, 0., 8000.);
621     TH1F* hdmH = new TH1F("hdmH", "Hard Trigger Multiplicity Distribution", 100, 0., 8000.);   
622
623     mbtH->SetXTitle("b [fm]");
624     hdtH->SetXTitle("b [fm]");    
625     mbmH->SetXTitle("Multiplicity");
626     hdmH->SetXTitle("Multiplicity");    
627
628     TCanvas *c0 = new TCanvas("c0","Trigger Simulation",400,10,600,700);    
629     c0->Divide(2,1);
630     TCanvas *c1 = new TCanvas("c1","Trigger Simulation",400,10,600,700);    
631     c1->Divide(1,2);
632
633     //
634     //
635     Init(1);
636     for (Int_t iev = 0; iev < n; iev++)
637     {
638         Float_t b, p, mult;
639         GetRandom(b, p, mult);
640         mbtH->Fill(b,1.);
641         hdtH->Fill(b, p);
642         mbmH->Fill(mult, 1.);
643         hdmH->Fill(mult, p);
644
645         c0->cd(1);
646         mbtH->Draw();
647         c0->cd(2);
648         hdtH->Draw();   
649         c0->Update();
650
651         c1->cd(1);
652         mbmH->Draw();
653         c1->cd(2);
654         hdmH->Draw();   
655         c1->Update();
656     }
657 }
658
659 void AliFastGlauber::GetRandom(Float_t& b, Float_t& p, Float_t& mult)
660 {
661     //
662     // Gives back a random impact parameter, hard trigger probability and multiplicity
663     //
664         b = fWSgeo->GetRandom();
665         Float_t mu = fWSN->Eval(b);
666         p = 1.-TMath::Exp(-mu);
667         mult = 6000./fWSN->Eval(1.) * mu;
668 }
669
670 void AliFastGlauber::GetRandom(Int_t& bin, Bool_t& hard)
671 {
672     //
673     // Gives back a random impact parameter bin, and hard trigger decission
674     //
675         Float_t b  = fWSgeo->GetRandom();
676         Float_t mu = fWSN->Eval(b) * fSigmaHard;
677         Float_t p  = 1.-TMath::Exp(-mu);
678         if (b < 5.) {
679             bin = 1;
680         } else if (b <  8.6) {
681             bin = 2;
682         } else if (b < 11.2) {
683             bin = 3;
684         } else if (b < 13.2) {
685             bin = 4;
686         } else if (b < 15.0) {
687             bin = 5;
688         } else {
689             bin = 6;
690         }
691         
692         hard = kFALSE;
693         
694         Float_t r = gRandom->Rndm();
695         
696         if (r < p) hard = kTRUE;
697 }
698
699
700 Float_t  AliFastGlauber::GetRandomImpactParameter(Float_t bmin, Float_t bmax)
701 {
702     //
703     // Gives back a random impact parameter in the range bmin .. bmax
704     //
705
706     Float_t b = -1.;
707     while(b < bmin || b > bmax)
708         b = fWSgeo->GetRandom();
709     return b;
710 }
711
712 Double_t AliFastGlauber::CrossSection(Double_t b1, Double_t b2)
713 {
714     //
715     // Return cross-section integrated from b1 to b2 
716     //
717     
718     return fWSgeo->Integral(b1, b2)/100.;
719 }
720
721 Double_t AliFastGlauber::FractionOfHardCrossSection(Double_t b1, Double_t b2)
722 {
723     //
724     // Return raction of hard cross-section integrated from b1 to b2 
725     //
726     
727     return fWSbinary->Integral(b1, b2)/fWSbinary->Integral(0., 100.);
728 }
729
730
731 Double_t AliFastGlauber::Binaries(Double_t b)
732 {
733     //
734     // Return number of binary collisions normalized to 1 at b=0
735     //
736     
737     return fWSN->Eval(b)/fWSN->Eval(0.001);
738 }