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