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