Calculation of number of participants added.
[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 <TStyle.h>
32 #include <TH1F.h>
33 #include <TH2F.h>
34 #include <TF1.h>
35 #include <TF2.h>
36 #include <TCanvas.h>
37 #include <TRandom.h>
38 #include <TFile.h>
39
40 ClassImp(AliFastGlauber)
41
42 TF1*    AliFastGlauber::fgWSb            = NULL;     
43 TF2*    AliFastGlauber::fgWSbz           = NULL;    
44 TF1*    AliFastGlauber::fgWSz            = NULL;     
45 TF1*    AliFastGlauber::fgWSta           = NULL;    
46 TF2*    AliFastGlauber::fgWStarfi        = NULL; 
47 TF2*    AliFastGlauber::fgWAlmond        = NULL; 
48 //TF2*    AliFastGlauber::fWAlmondCurrent = NULL; 
49 TF1*    AliFastGlauber::fgWStaa          = NULL;   
50 TF1*    AliFastGlauber::fgWSgeo          = NULL;   
51 TF1*    AliFastGlauber::fgWSbinary       = NULL;   
52 TF1*    AliFastGlauber::fgWSN            = NULL;   
53 TF1*    AliFastGlauber::fgWPathLength0   = NULL;   
54 TF1*    AliFastGlauber::fgWPathLength    = NULL;
55 TF1*    AliFastGlauber::fgWEnergyDensity = NULL;   
56 TF1*    AliFastGlauber::fgWIntRadius     = NULL;   
57 TF2*    AliFastGlauber::fgWKParticipants = NULL; 
58 TF1*    AliFastGlauber::fgWParticipants  = NULL; 
59 Float_t AliFastGlauber::fgBMax           = 0.;
60
61 AliFastGlauber::AliFastGlauber()
62 {
63 //  Default Constructor
64 //
65 //  Defaults for Pb
66     SetWoodSaxonParameters(6.624, 0.549, 0.00, 7.69e-4);
67     SetHardCrossSection();
68     SetMaxImpact();
69     SetLengthDefinition();
70 }
71
72 void AliFastGlauber::Init(Int_t mode)
73 {
74 // Initialisation
75 //    mode = 0; all functions are calculated 
76 //    mode = 1; overlap function is read from file (for Pb-Pb only)
77 //    mode = 2; interaction almond functions are read from file 
78 //              (for Pb-Pb only); USE THIS FOR PATH LENGTH CALC.!
79 //
80
81 //
82 //  Wood-Saxon
83 //
84     fgWSb = new TF1("WSb", WSb, 0, fgBMax, 4);
85     fgWSb->SetParameter(0, fWSr0);
86     fgWSb->SetParameter(1, fWSd);
87     fgWSb->SetParameter(2, fWSw);
88     fgWSb->SetParameter(3, fWSn);
89
90     fgWSbz = new TF2("WSbz", WSbz, 0, fgBMax, 4);
91     fgWSbz->SetParameter(0, fWSr0);
92     fgWSbz->SetParameter(1, fWSd);
93     fgWSbz->SetParameter(2, fWSw);
94     fgWSbz->SetParameter(3, fWSn);
95
96     fgWSz = new TF1("WSz", WSz, 0, fgBMax, 5);
97     fgWSz->SetParameter(0, fWSr0);
98     fgWSz->SetParameter(1, fWSd);
99     fgWSz->SetParameter(2, fWSw);
100     fgWSz->SetParameter(3, fWSn);
101
102 //
103 //  Thickness
104 //
105     fgWSta = new TF1("WSta", WSta, 0., fgBMax, 0);
106     
107 //
108 //  Overlap Kernel
109 //
110     fgWStarfi = new TF2("WStarfi", WStarfi, 0., fgBMax, 0., TMath::Pi(), 1);
111     fgWStarfi->SetParameter(0, 0.);     
112     fgWStarfi->SetNpx(200);     
113     fgWStarfi->SetNpy(20);    
114 //
115 //  Participants Kernel
116 //
117     fgWKParticipants = new TF2("WKParticipants", WKParticipants, 0., fgBMax, 0., TMath::Pi(), 1);
118     fgWKParticipants->SetParameter(0, 0.);     
119     fgWKParticipants->SetNpx(200);     
120     fgWKParticipants->SetNpy(20);      
121 //
122 //  Almond shaped interaction region
123 //
124     fgWAlmond = new TF2("WAlmond", WAlmond, -fgBMax, fgBMax, -fgBMax, fgBMax, 1);
125     fgWAlmond->SetParameter(0, 0.);     
126     fgWAlmond->SetNpx(200);     
127     fgWAlmond->SetNpy(200);  
128   
129     if(mode==2) {
130       printf(" Reading interaction almonds from file\n");
131       Char_t almondName[100];
132       TFile* ff = new TFile("$(ALICE_ROOT)/FASTSIM/data/glauberPbPb.root");
133       for(Int_t k=0; k<40; k++) {
134         sprintf(almondName,"WAlmondFixedB%d",k);
135         fWAlmondCurrent = (TF2*)ff->Get(almondName);
136         new(&fWAlmondFixedB[k]) TF2(*fWAlmondCurrent);
137       }
138     }
139 //
140 //  Path Length as a function of Phi
141 //    
142     fgWPathLength0 = new TF1("WPathLength0", WPathLength0, -TMath::Pi(), TMath::Pi(), 2);
143     fgWPathLength0->SetParameter(0, 0.);
144 //  Pathlength definition     
145     fgWPathLength0->SetParameter(1, 0.);     
146
147     fgWPathLength = new TF1("WPathLength", WPathLength, -TMath::Pi(), TMath::Pi(), 3);
148 //  Impact Parameter
149     fgWPathLength->SetParameter(0, 0.);    
150 //  Number of interactions used for average
151     fgWPathLength->SetParameter(1, 1000.);    
152 //  Pathlength definition
153     fgWPathLength->SetParameter(2, 0);    
154
155     fgWIntRadius = new TF1("WIntRadius", WIntRadius, 0., fgBMax, 1);
156     fgWIntRadius->SetParameter(0, 0.);    
157
158
159 //
160 //  Overlap
161 //
162     if (! mode) {
163         fgWStaa = new TF1("WStaa", WStaa, 0., fgBMax, 0);
164         fgWStaa->SetNpx(100);
165         fgWParticipants = new TF1("WParticipants", WParticipants, 0., fgBMax, 0);
166         fgWParticipants->SetNpx(100);
167
168     } else {
169         TFile* f = new TFile("$(ALICE_ROOT)/FASTSIM/data/glauberPbPb.root");
170         fgWStaa = (TF1*) f->Get("WStaa");
171         fgWParticipants = (TF1*) f->Get("WParticipants");
172     }
173     
174 //
175 //  Participants
176 //
177     
178 //
179     fgWEnergyDensity = new TF1("WEnergyDensity", WEnergyDensity, 0., 2. * fWSr0, 1);
180     fgWEnergyDensity->SetParameter(0, fWSr0 + 1.);
181     
182 //
183 //  Geometrical Cross-Section
184 //
185     fgWSgeo = new TF1("WSgeo", WSgeo, 0., fgBMax, 0);
186     fgWSgeo->SetNpx(100);
187 //
188 //  Hard cross section (~ binary collisions)
189 //
190     fgWSbinary = new TF1("WSbinary", WSbinary, 0., fgBMax, 1);
191     fgWSbinary->SetParameter(0, fSigmaHard); // mb
192     fgWSbinary->SetNpx(100);
193 //
194 // Hard collisions per event
195 //
196     fgWSN = new TF1("WSN", WSN, 0., fgBMax, 1);
197     fgWSN->SetNpx(100);
198 }
199
200 void AliFastGlauber::DrawWSb()
201 {
202 //
203 //  Draw Wood-Saxon Nuclear Density Function
204 //
205     TCanvas *c1 = new TCanvas("c1","Wood Saxon",400,10,600,700);
206     c1->cd();
207     fgWSb->Draw();
208 }
209
210 void AliFastGlauber::DrawOverlap()
211 {
212 //
213 //  Draw Overlap Function
214 //
215     TCanvas *c2 = new TCanvas("c2","Overlap",400,10,600,700);
216     c2->cd();
217     fgWStaa->Draw();
218 }
219
220 void AliFastGlauber::DrawParticipants()
221 {
222 //
223 //  Draw Number of Participants
224 //
225     TCanvas *c2 = new TCanvas("c2","Participants",400,10,600,700);
226     c2->cd();
227     fgWParticipants->Draw();
228 }
229
230 void AliFastGlauber::DrawThickness()
231 {
232 //
233 //  Draw Thickness Function
234 //
235     TCanvas *c3 = new TCanvas("c3","Thickness",400,10,600,700);
236     c3->cd();
237     fgWSta->Draw();
238 }
239
240 void AliFastGlauber::DrawGeo()
241 {
242 //
243 //  Draw Geometrical Cross-Section
244 //
245     TCanvas *c3 = new TCanvas("c3","Geometrical Cross-Section",400,10,600,700);
246     c3->cd();
247     fgWSgeo->Draw();
248 }
249
250 void AliFastGlauber::DrawBinary()
251 {
252 //
253 //  Draw Binary Cross-Section
254 //
255     TCanvas *c4 = new TCanvas("c4","Binary Cross-Section",400,10,600,700);
256     c4->cd();
257     fgWSbinary->Draw();
258 }
259
260 void AliFastGlauber::DrawN()
261 {
262 //
263 //  Draw Binaries per event
264 //
265     TCanvas *c5 = new TCanvas("c5","Binaries per event",400,10,600,700);
266     c5->cd();
267     fgWSN->Draw();
268 }
269
270 void AliFastGlauber::DrawKernel(Double_t b)
271 {
272 //
273 //  Draw Kernel
274 //
275     TCanvas *c6 = new TCanvas("c6","Kernel",400,10,600,700);
276     c6->cd();
277     fgWStarfi->SetParameter(0, b);
278     fgWStarfi->Draw();
279 }
280
281 void AliFastGlauber::DrawAlmond(Double_t b)
282 {
283 //
284 //  Draw Interaction Almond
285 //
286     TCanvas *c7 = new TCanvas("c7","Almond",400,10,600,700);
287     c7->cd();
288     fgWAlmond->SetParameter(0, b);
289     fgWAlmond->Draw();
290 }
291
292 void AliFastGlauber::DrawPathLength0(Double_t b, Int_t iopt)
293 {
294 //
295 //  Draw Path Length
296 //
297     TCanvas *c8 = new TCanvas("c8","Path Length",400,10,600,700);
298     c8->cd();
299     fgWPathLength0->SetParameter(0, b);
300     fgWPathLength0->SetParameter(1, Double_t(iopt));
301     fgWPathLength0->SetMinimum(0.); 
302     fgWPathLength0->SetMaximum(10.); 
303     fgWPathLength0->Draw();
304 }
305
306 void AliFastGlauber::DrawPathLength(Double_t b , Int_t ni, Int_t iopt)
307 {
308 //
309 //  Draw Path Length
310 //
311     TCanvas *c9 = new TCanvas("c9","Path Length",400,10,600,700);
312     c9->cd();
313     fgWAlmond->SetParameter(0, b);
314
315     fgWPathLength->SetParameter(0, b);
316     fgWPathLength->SetParameter(1, Double_t (ni));
317     fgWPathLength->SetParameter(2, Double_t (iopt));
318     fgWPathLength->SetMinimum(0.); 
319     fgWPathLength->SetMaximum(10.); 
320     fgWPathLength->Draw();
321 }
322
323 void AliFastGlauber::DrawIntRadius(Double_t b)
324 {
325 //
326 //  Draw Interaction Radius
327 //
328     TCanvas *c10 = new TCanvas("c10","Interaction Radius",400,10,600,700);
329     c10->cd();
330     fgWIntRadius->SetParameter(0, b);
331     fgWIntRadius->SetMinimum(0.);
332     fgWIntRadius->Draw();
333 }
334
335 void AliFastGlauber::DrawEnergyDensity()
336 {
337 //
338 //  Draw energy density
339 //
340     TCanvas *c11 = new TCanvas("c11","Energy Density",400, 10, 600, 700);
341     c11->cd();
342     fgWEnergyDensity->SetMinimum(0.);
343     fgWEnergyDensity->Draw();
344 }
345
346 Double_t AliFastGlauber::WSb(Double_t* x, Double_t* par)
347 {
348 //
349 //  Woods-Saxon Parameterisation
350 //  as a function of radius
351 //
352     Double_t xx  = x[0];
353     Double_t r0  = par[0];
354     Double_t d   = par[1];
355     Double_t w   = par[2];
356     Double_t n   = par[3];
357     
358     Double_t y  = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
359
360     return y;
361 }
362
363 Double_t AliFastGlauber::WSbz(Double_t* x, Double_t* par)
364 {
365 //
366 //  Wood Saxon Parameterisation
367 //  as a function of z and  b
368 //
369     Double_t bb  = x[0];
370     Double_t zz  = x[1];
371     Double_t r0  = par[0];
372     Double_t d   = par[1];
373     Double_t w   = par[2];
374     Double_t n   = par[3];
375     Double_t xx  = TMath::Sqrt(bb*bb+zz*zz);
376     Double_t y  = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
377
378     return y;
379 }
380
381 Double_t AliFastGlauber::WSz(Double_t* x, Double_t* par)
382 {
383 //
384 //  Wood Saxon Parameterisation
385 //  as a function of z for fixed b
386 //
387     Double_t bb  = par[4];
388     Double_t zz  = x[0];
389     Double_t r0  = par[0];
390     Double_t d   = par[1];
391     Double_t w   = par[2];
392     Double_t n   = par[3];
393     Double_t xx  = TMath::Sqrt(bb*bb+zz*zz);
394     Double_t y  = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
395
396     return y;
397 }
398
399 Double_t AliFastGlauber::WSta(Double_t* x, Double_t* /*par*/)
400 {
401 //
402 //  Thickness function 
403 //
404     Double_t b  = x[0];
405     fgWSz->SetParameter(4, b);
406     Double_t y  = 2. * fgWSz->Integral(0., fgBMax);
407     return y;
408 }
409
410
411
412 Double_t AliFastGlauber::WStarfi(Double_t* x, Double_t* par)
413 {
414 //
415 //  Kernel for overlap function
416 //
417     Double_t b    = par[0];
418     Double_t r1   = x[0];
419     Double_t phi  = x[1];
420     Double_t r2   = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi)); 
421     Double_t y    = r1 * fgWSta->Eval(r1) * fgWSta->Eval(r2);
422     return y;
423 }
424
425 Double_t AliFastGlauber::WKParticipants(Double_t* x, Double_t* par)
426 {
427 //
428 //  Kernel for number of participants
429 //
430     Double_t b    = par[0];
431     Double_t r1   = x[0];
432     Double_t phi  = x[1];
433     Double_t r2   = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi)); 
434     Double_t xsi  = fgWSta->Eval(r2) * 5.6;
435     Double_t a = 208;
436     Double_t sum = a * xsi;
437     Double_t y   = sum;
438     for (Int_t i = 1; i <= 208; i++)
439     {
440         a--;
441         sum *= (-xsi) * a / Float_t(i+1);
442         y  += sum;
443     }
444     
445     y    = r1 * fgWSta->Eval(r1) * y;
446     return y;
447 }
448
449
450 Double_t AliFastGlauber::WAlmond(Double_t* x, Double_t* par)
451 {
452 //
453 //  Almond shaped interaction region
454 //
455     Double_t b    = par[0];
456     Double_t xx   = x[0] + b/2.;
457     Double_t yy   = x[1];
458     Double_t r1   = TMath::Sqrt(xx * xx + yy * yy);
459     Double_t phi  = TMath::ATan2(yy,xx);
460     
461     Double_t r2   = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi)); 
462 //
463 //  Interaction probability calculated as product of thicknesses
464 //
465     Double_t y    = fgWSta->Eval(r1) * fgWSta->Eval(r2);
466     return y;
467 }
468
469 Double_t AliFastGlauber::WIntRadius(Double_t* x, Double_t* par)
470 {
471 //
472 //  Average radius at which interaction takes place
473 //
474 //  Radius in the Almond
475     Double_t r    = x[0];
476 //  Impact parameter
477     Double_t b    = par[0];
478     fgWAlmond->SetParameter(0, b);
479 //  Steps in phi
480     Double_t dphi = 2. * TMath::Pi() / 100.;
481 //  Average over phi    
482     Double_t phi  = 0.;
483     Double_t y    = 0.;
484
485     for (Int_t i = 0; i < 100; i++) {
486         Double_t xx = r * TMath::Cos(phi);
487         Double_t yy = r * TMath::Sin(phi);
488         y   += fgWAlmond->Eval(xx,yy);
489         phi += dphi;
490     } // phi loop
491 // Result multiplied by Jacobian (2 pi r)     
492     return (2. * TMath::Pi() * y * r / 100.);
493 }
494
495 Double_t AliFastGlauber::WPathLength0(Double_t* x, Double_t* par)
496 {
497 //
498 //  Path Length as a function of phi for interaction point fixed at (0,0)
499 //
500 //
501 //  Steps in r 
502     const Int_t    kNp  = 100;
503     const Double_t kDr  = fgBMax/Double_t(kNp);
504 //  Impact parameter    
505     Double_t b      = par[0];
506 //  Path Length definition
507     Int_t    iopt   = Int_t(par[1]);
508     
509 //  Phi direction in Almond
510     Double_t phi0   = x[0];
511     Double_t r  = 0.;
512     Double_t rw = 0.;
513     Double_t w  = 0.;
514 //  Step along radial direction phi   
515     for (Int_t i = 0; i < kNp; i++) {
516 //
517 //  Transform into target frame
518 //
519         Double_t xx   = r * TMath::Cos(phi0) + b / 2.;
520         Double_t yy   = r * TMath::Sin(phi0);
521         Double_t phi  = TMath::ATan2(yy, xx);
522         
523         Double_t r1   = TMath::Sqrt(xx * xx + yy * yy);
524 // Radius in projectile frame
525         Double_t r2   = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi)); 
526         Double_t y    = fgWSta->Eval(r1) * fgWSta->Eval(r2);
527
528         rw += y * r;
529         w  += y;
530         r  += kDr;
531     } // radial steps
532 //
533 //  My length definition (is exact for hard disk)
534     if (!iopt) {
535         return (2. * rw / w);
536     } else {
537         return TMath::Sqrt(2. * rw * kDr / fgWSta->Eval(0.01) / fgWSta->Eval(0.01));
538     }
539 }
540
541 Double_t AliFastGlauber::WPathLength(Double_t* x, Double_t* par)
542 {
543 //
544 //  Path Length as a function of phi 
545 //  Interaction point from random distribution
546 //
547 //
548 //  r-steps
549 // 
550     const Int_t    kNp   = 100;
551     const Double_t kDr  = fgBMax/Double_t(kNp);
552 //  Number of interactions
553     const Int_t    kNpi  = Int_t (par[1]);
554
555 //
556 //  Impact parameter    
557     Double_t b      = par[0];
558 //  Path Length definition 
559     Int_t    iopt   = Int_t(par[2]);
560 //  Phi direction
561     Double_t phi0   = x[0];
562
563     printf("phi0 %f \n", phi0);
564     
565 //  Path length 
566     Double_t l = 0.;
567     
568     for (Int_t in = 0; in < kNpi; in ++) {
569         Double_t rw = 0.;
570         Double_t w  = 0.;
571         
572         // Interaction point
573         Double_t x0, y0;
574         fgWAlmond->GetRandom2(x0, y0);
575 // Initial radius
576         Double_t r0  = TMath::Sqrt(x0 * x0 + y0 * y0);
577         Int_t    nps = Int_t ((fgBMax - r0)/kDr) - 1;
578         
579         Double_t r  = 0.;
580 // Radial steps
581         for (Int_t i = 0; (i < nps ); i++) {
582             
583 // Transform into target frame
584             Double_t xx   = x0 + r * TMath::Cos(phi0) + b / 2.;
585             Double_t yy   = y0 + r * TMath::Sin(phi0);
586             Double_t phi  = TMath::ATan2(yy, xx);
587             Double_t r1   = TMath::Sqrt(xx * xx + yy * yy);
588 // Radius in projectile frame
589             Double_t r2   = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi)); 
590             Double_t y    = fgWSta->Eval(r1) * fgWSta->Eval(r2);
591             
592             rw += y * r;
593             w  += y;
594             r  += kDr;
595         } // steps
596 // Average over interactions
597         if (!iopt) {
598             l += (2. * rw / w);
599         } else {
600             l+= 2. * rw * kDr / fgWSta->Eval(0.01) / fgWSta->Eval(0.01);
601         }
602     } // interactions
603     if (!iopt) 
604         return (l / Double_t(kNpi));
605     else 
606         return (TMath::Sqrt(l / Double_t(kNpi)));
607 }
608
609 Double_t AliFastGlauber::WStaa(Double_t* x, Double_t* /*par*/)
610 {
611 //
612 //  Overlap function
613 //
614     Double_t b    = x[0];
615     fgWStarfi->SetParameter(0, b);
616 /*
617     Double_t al[2];
618     Double_t bl[2];
619     al[0] = 0.;
620     al[1] = 0.;
621     bl[0] = 6.6;
622     bl[1] = TMath::Pi();
623     Double_t err;
624     
625     Double_t y =  2. * fgWStarfi->IntegralMultiple(2, al, bl, 0.001, err);
626     printf("WStaa: %f %f %f\n", b, y, err);
627 */
628 //
629 //  MC Integration
630 //
631     Double_t y = 0;
632     for (Int_t i = 0; i < 100000; i++)
633     {
634         Double_t phi = TMath::Pi() * gRandom->Rndm();
635         Double_t b1  = fgBMax       * gRandom->Rndm();  
636         y += fgWStarfi->Eval(b1, phi);
637     }
638     y *= 2. * 0.1 *  208. * 208. * TMath::Pi() * fgBMax / 100000.;
639     return y;
640 }
641
642 Double_t AliFastGlauber::WParticipants(Double_t* x, Double_t* /*par*/)
643 {
644 //
645 //  Overlap function
646 //
647     Double_t b    = x[0];
648     fgWKParticipants->SetParameter(0, b);
649 //
650 //  MC Integration
651 //
652     Double_t y = 0;
653     for (Int_t i = 0; i < 100000; i++)
654     {
655         Double_t phi = TMath::Pi() * gRandom->Rndm();
656         Double_t b1  = fgBMax      * gRandom->Rndm();   
657         y += fgWKParticipants->Eval(b1, phi);
658     }
659     y *= 4. *  208. * TMath::Pi() * fgBMax / 100000.;
660     return y;
661 }
662
663
664 Double_t AliFastGlauber::WSgeo(Double_t* x, Double_t* /*par*/)
665 {
666 //
667 //  Geometrical Cross-Section
668 //
669     Double_t b    = x[0];
670     Double_t taa  = fgWStaa->Eval(b);
671     const Double_t kSigma = 55.6; // mbarn
672     
673     Double_t y    = 2. * TMath::Pi() * b * (1. - TMath::Exp(- kSigma * taa)); // fm
674     return y;
675 }
676
677
678 Double_t AliFastGlauber::WSbinary(Double_t* x, Double_t* par)
679 {
680 //
681 //  Number of binary collisions
682 //
683     Double_t b     = x[0];
684     Double_t sigma = par[0];
685     Double_t taa   = fgWStaa->Eval(b);
686     
687     Double_t y    = 2. * TMath::Pi() * b * sigma * taa; // fm
688     return y;
689 }
690
691 Double_t AliFastGlauber::WSN(Double_t* x, Double_t* /*par*/)
692 {
693 //
694 //  Number of hard processes per event
695 //
696     Double_t b     = x[0];
697     Double_t y     = fgWSbinary->Eval(b)/fgWSgeo->Eval(b);
698     return y;
699 }
700
701 Double_t AliFastGlauber::WEnergyDensity(Double_t* x, Double_t* par)
702 {
703 //
704 //  Initial energy density as a function of the impact parameter
705 //
706     Double_t b     = x[0];
707     Double_t rA    = par[0];
708 //
709 //  Attention: area of transverse reaction zone in hard-sphere approximation !     
710     Double_t saa   = (TMath::Pi() - 2. * TMath::ASin(b/ 2./ rA)) * rA * rA 
711         - b * TMath::Sqrt(rA * rA - b * b/ 4.);
712     Double_t taa   = fgWStaa->Eval(b);
713     
714     return (taa/saa);
715 }
716
717 void AliFastGlauber::SimulateTrigger(Int_t n)
718 {
719     //
720     //  Simulates Trigger
721     //
722     TH1F* mbtH = new TH1F("mbtH", "MB Trigger b-Distribution",   100, 0., 20.);
723     TH1F* hdtH = new TH1F("hdtH", "Hard Trigger b-Distribution", 100, 0., 20.);   
724     TH1F* mbmH = new TH1F("mbmH", "MB Trigger Multiplicity Distribution",   100, 0., 8000.);
725     TH1F* hdmH = new TH1F("hdmH", "Hard Trigger Multiplicity Distribution", 100, 0., 8000.);   
726
727     mbtH->SetXTitle("b [fm]");
728     hdtH->SetXTitle("b [fm]");    
729     mbmH->SetXTitle("Multiplicity");
730     hdmH->SetXTitle("Multiplicity");    
731
732     TCanvas *c0 = new TCanvas("c0","Trigger Simulation",400,10,600,700);    
733     c0->Divide(2,1);
734     TCanvas *c1 = new TCanvas("c1","Trigger Simulation",400,10,600,700);    
735     c1->Divide(1,2);
736
737     //
738     //
739     Init(1);
740     for (Int_t iev = 0; iev < n; iev++)
741     {
742         Float_t b, p, mult;
743         GetRandom(b, p, mult);
744         mbtH->Fill(b,1.);
745         hdtH->Fill(b, p);
746         mbmH->Fill(mult, 1.);
747         hdmH->Fill(mult, p);
748
749         c0->cd(1);
750         mbtH->Draw();
751         c0->cd(2);
752         hdtH->Draw();   
753         c0->Update();
754
755         c1->cd(1);
756         mbmH->Draw();
757         c1->cd(2);
758         hdmH->Draw();   
759         c1->Update();
760     }
761 }
762
763 void AliFastGlauber::GetRandom(Float_t& b, Float_t& p, Float_t& mult)
764 {
765     //
766     // Gives back a random impact parameter, hard trigger probability and multiplicity
767     //
768         b = fgWSgeo->GetRandom();
769         Float_t mu = fgWSN->Eval(b);
770         p = 1.-TMath::Exp(-mu);
771         mult = 6000./fgWSN->Eval(1.) * mu;
772 }
773
774 void AliFastGlauber::GetRandom(Int_t& bin, Bool_t& hard)
775 {
776     //
777     // Gives back a random impact parameter bin, and hard trigger decission
778     //
779         Float_t b  = fgWSgeo->GetRandom();
780         Float_t mu = fgWSN->Eval(b) * fSigmaHard;
781         Float_t p  = 1.-TMath::Exp(-mu);
782         if (b < 5.) {
783             bin = 1;
784         } else if (b <  8.6) {
785             bin = 2;
786         } else if (b < 11.2) {
787             bin = 3;
788         } else if (b < 13.2) {
789             bin = 4;
790         } else if (b < 15.0) {
791             bin = 5;
792         } else {
793             bin = 6;
794         }
795         
796         hard = kFALSE;
797         
798         Float_t r = gRandom->Rndm();
799         
800         if (r < p) hard = kTRUE;
801 }
802
803
804
805 Float_t  AliFastGlauber::GetRandomImpactParameter(Float_t bmin, Float_t bmax)
806 {
807     //
808     // Gives back a random impact parameter in the range bmin .. bmax
809     //
810
811     Float_t b = -1.;
812     while(b < bmin || b > bmax)
813         b = fgWSgeo->GetRandom();
814     return b;
815 }
816
817 Double_t AliFastGlauber::CrossSection(Double_t b1, Double_t b2)
818 {
819     //
820     // Return cross-section integrated from b1 to b2 
821     //
822     
823     return fgWSgeo->Integral(b1, b2)/100.;
824 }
825
826 Float_t AliFastGlauber::GetNumberOfParticipants(Float_t  b)
827 {
828     //
829     // Gives back the number of participants for impact parameter b
830     //
831     
832     return (fgWParticipants->Eval(b));
833 }
834
835 Double_t AliFastGlauber::FractionOfHardCrossSection(Double_t b1, Double_t b2)
836 {
837     //
838     // Return raction of hard cross-section integrated from b1 to b2 
839     //
840     
841     return fgWSbinary->Integral(b1, b2)/fgWSbinary->Integral(0., 100.);
842 }
843
844
845 Double_t AliFastGlauber::Binaries(Double_t b)
846 {
847     //
848     // Return number of binary collisions normalized to 1 at b=0
849     //
850     
851     return fgWSN->Eval(b)/fgWSN->Eval(0.001);
852 }
853
854 //=================== Added by A. Dainese 11/02/04 ===========================
855 void AliFastGlauber::SetCentralityClass(Double_t xsecFrLow,Double_t xsecFrUp)
856 {
857   //
858   // Set limits of centrality class as fractions of the geomtrical 
859   // cross section  
860   //
861   if(xsecFrLow>1. || xsecFrUp>1. || xsecFrLow>xsecFrUp) {
862     printf(" Please set 0 <= xsecFrLow <= xsecFrUp <= 1\n");
863     return;
864   }
865
866   Double_t bLow=0.,bUp=0.;
867   Double_t xsecFr=0.;
868
869   while(xsecFr<xsecFrLow) {
870     xsecFr = fgWSgeo->Integral(0.,bLow)/fgWSgeo->Integral(0.,100.);
871     bLow += 0.1;
872   }
873   bUp = bLow;
874   while(xsecFr<xsecFrUp) {
875     xsecFr = fgWSgeo->Integral(0.,bUp)/fgWSgeo->Integral(0.,100.);
876     bUp += 0.1;
877   }
878
879   printf(" Centrality class: %4.2f-%4.2f; %4.1f < b < %4.1f fm\n",
880          xsecFrLow,xsecFrUp,bLow,bUp);
881
882   fgWSbinary->SetRange(bLow,bUp);
883
884   return;
885 }
886
887 void AliFastGlauber::GetRandomBHard(Double_t& b)
888 {
889   //
890   // Get random impact parameter according to distribution of 
891   // hard (binary) cross-section, in the range defined by the centrality class
892   //
893   b = fgWSbinary->GetRandom();
894   Int_t bin = 2*(Int_t)b;
895   if( (b-(Int_t)b) > 0.5) bin++;
896   fWAlmondCurrent = &fWAlmondFixedB[bin]; 
897   return;
898 }
899
900 void AliFastGlauber::GetRandomXY(Double_t& x,Double_t& y)
901 {
902   //
903   // Get random position of parton production point according to 
904   // product of thickness functions
905   //
906   fWAlmondCurrent->GetRandom2(x,y);
907   return;
908 }
909
910 void AliFastGlauber::GetRandomPhi(Double_t& phi) 
911 {
912   //
913   // Get random parton azimuthal propagation direction
914   //
915   phi = 2.*TMath::Pi()*gRandom->Rndm();
916   return;
917 }
918
919 Double_t AliFastGlauber::CalculateLength(Double_t b,
920                                          Double_t x0,Double_t y0,Double_t phi0)
921 {
922   // 
923   // Calculate path length for a parton with production point (x0,y0)
924   // and propagation direction (ux=cos(phi0),uy=sin(phi0)) 
925   // in a collision with impact parameter b
926   //
927
928   // number of steps in l
929   const Int_t    kNp  = 100;
930   const Double_t kDl  = fgBMax/Double_t(kNp);
931
932   Double_t l,integral,integral1,integral2,r0,xx,yy,phi,r1,r2,ell;
933   Double_t rhocoll,rhocollHalfMax;
934   Int_t i,nps;
935
936   if(fEllDef==1) {
937     //
938     // Definition 1:
939     // 
940     // ell = 2 * \int_0^\infty dl*l*\rho_coll(x0+l*ux,y0+l*uy) /
941     //           \int_0^\infty dl*\rho_coll(x0+l*ux,y0+l*uy)
942     //
943     
944
945     l  = 0.;
946     integral1 = 0.;
947     integral2 = 0.;
948     
949     // Initial radius
950     r0  = TMath::Sqrt(x0 * x0 + y0 * y0);
951     nps = Int_t ((fgBMax - r0)/kDl) - 1;
952     
953     // Radial steps
954     for (i = 0; i < nps; i++) {
955       
956       // Transform into target frame
957       xx   = x0 + l * TMath::Cos(phi0) + b / 2.;
958       yy   = y0 + l * TMath::Sin(phi0);
959       phi  = TMath::ATan2(yy, xx);
960       r1   = TMath::Sqrt(xx * xx + yy * yy);
961       // Radius in projectile frame
962       r2   = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi)); 
963       rhocoll = fgWSta->Eval(r1) * fgWSta->Eval(r2);
964       
965       integral1 += rhocoll * l * kDl;
966       integral2 += rhocoll * kDl;
967       l  += kDl;
968     } // steps
969     
970     ell = (2. * integral1 / integral2);
971     return ell;
972     
973   } else if(fEllDef==2) {
974     //
975     // Definition 2:
976     // 
977     // ell = \int_0^\infty dl*
978     //       \Theta(\rho_coll(x0+l*ux,y0+l*uy)-0.5*\rho_coll(0,0))
979     //
980
981     l  = 0.;
982     integral = 0.;
983     
984     // Initial radius
985     r0  = TMath::Sqrt(x0 * x0 + y0 * y0);
986     nps = Int_t ((fgBMax - r0)/kDl) - 1;
987
988     rhocollHalfMax = 0.5*fWAlmondCurrent->Eval(0.,0.);
989
990     // Radial steps
991     for (i = 0; i < nps; i++) {
992       
993       // Transform into target frame
994       xx   = x0 + l * TMath::Cos(phi0) + b / 2.;
995       yy   = y0 + l * TMath::Sin(phi0);
996       phi  = TMath::ATan2(yy, xx);
997       r1   = TMath::Sqrt(xx * xx + yy * yy);
998       // Radius in projectile frame
999       r2   = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi)); 
1000       rhocoll = fgWSta->Eval(r1) * fgWSta->Eval(r2);
1001       
1002       if(rhocoll>rhocollHalfMax) integral += kDl;
1003
1004       l  += kDl;
1005     } // steps
1006     
1007     ell = integral;
1008     return ell;
1009
1010   } else {
1011     printf(" Wrong length definition setting!\n");
1012     return -1.;
1013   }
1014 }
1015
1016 void AliFastGlauber::GetLength(Double_t& ell,Double_t b)
1017 {
1018   //
1019   // Return length from random b, x0, y0, phi0 
1020   //
1021   Double_t x0,y0,phi0;
1022   if(b<0.) GetRandomBHard(b);
1023   GetRandomXY(x0,y0);
1024   GetRandomPhi(phi0);
1025
1026   ell = CalculateLength(b,x0,y0,phi0);
1027
1028   return;
1029 }
1030
1031 void AliFastGlauber::GetLengthsBackToBack(Double_t& ell1,Double_t& ell2,
1032                                           Double_t b)
1033 {
1034   //
1035   // Return 2 lengths back to back from random b, x0, y0, phi0 
1036   //
1037   Double_t x0,y0,phi0;
1038   if(b<0.) GetRandomBHard(b);
1039   GetRandomXY(x0,y0);
1040   GetRandomPhi(phi0);
1041   Double_t phi0plusPi = phi0+TMath::Pi();
1042
1043   ell1 = CalculateLength(b,x0,y0,phi0);
1044   ell2 = CalculateLength(b,x0,y0,phi0plusPi);
1045
1046   return;
1047 }
1048
1049 void AliFastGlauber::GetLengthsForPythia(Int_t n,Double_t* phi,Double_t* ell, Double_t b)
1050 {
1051   //
1052   // Returns lenghts for n partons with azimuthal angles phi[n] 
1053   // from random b, x0, y0
1054   //
1055   Double_t x0, y0;
1056   if(b < 0.) GetRandomBHard(b);
1057   GetRandomXY(x0,y0);
1058
1059   for(Int_t i = 0; i< n; i++) ell[i] = CalculateLength(b,x0,y0,phi[i]);
1060
1061   return;
1062 }
1063
1064 void AliFastGlauber::PlotBDistr(Int_t n)
1065 {
1066   // 
1067   // Plot distribution of n impact parameters
1068   //
1069   Double_t b;
1070   TH1F *hB = new TH1F("hB","dN/db",100,0,fgBMax); 
1071   hB->SetXTitle("b [fm]");
1072   hB->SetYTitle("dN/db [a.u.]");
1073   hB->SetFillColor(3);
1074
1075   for(Int_t i=0; i<n; i++) {
1076     GetRandomBHard(b);
1077     hB->Fill(b);
1078   }
1079
1080   TCanvas *cB = new TCanvas("cB","Impact parameter distribution",0,0,500,500);
1081   cB->cd();
1082   hB->Draw();
1083
1084   return;
1085 }
1086
1087 void AliFastGlauber::PlotLengthDistr(Int_t n,Bool_t save,Char_t *fname)
1088 {
1089   //
1090   // Plot length distribution
1091   //
1092   Double_t ell;
1093   TH1F *hEll = new TH1F("hEll","Length distribution",16,-0.5,15.5); 
1094   hEll->SetXTitle("Transverse path length, L [fm]");
1095   hEll->SetYTitle("Probability");
1096   hEll->SetFillColor(2);
1097
1098   for(Int_t i=0; i<n; i++) {
1099     GetLength(ell);
1100     hEll->Fill(ell);
1101   }
1102   hEll->Scale(1/(Double_t)n);
1103
1104   TCanvas *cL = new TCanvas("cL","Length distribution",0,0,500,500);
1105   cL->cd();
1106   hEll->Draw();
1107
1108   if(save) {
1109     TFile *f = new TFile(fname,"recreate");
1110     hEll->Write();
1111     f->Close();
1112   }
1113   return;
1114 }
1115
1116 void AliFastGlauber::PlotLengthB2BDistr(Int_t n,Bool_t save,Char_t *fname)
1117 {
1118   //
1119   // Plot lengths back-to-back distributions
1120   //
1121   Double_t ell1,ell2;
1122   TH2F *hElls = new TH2F("hElls","Lengths back-to-back",100,0,15,100,0,15); 
1123   hElls->SetXTitle("Transverse path length, L [fm]");
1124   hElls->SetYTitle("Transverse path length, L [fm]");
1125
1126   for(Int_t i=0; i<n; i++) {
1127     GetLengthsBackToBack(ell1,ell2);
1128     hElls->Fill(ell1,ell2);
1129   }
1130   hElls->Scale(1/(Double_t)n);
1131
1132
1133   TCanvas *cLs = new TCanvas("cLs","Length back-to-back distribution",0,0,500,500);
1134   gStyle->SetPalette(1,0);
1135   cLs->cd();
1136   hElls->Draw("col,Z");
1137
1138   if(save) {
1139     TFile *f = new TFile(fname,"recreate");
1140     hElls->Write();
1141     f->Close();
1142   }
1143   return;
1144 }
1145
1146 void AliFastGlauber::StoreAlmonds()
1147 {
1148   //
1149   // Store in glauberPbPb.root 40 almonds for b = (0.25+k*0.5) fm (k=0->39)
1150   //
1151
1152   Char_t almondName[100];
1153   TFile* ff = new TFile("$(ALICE_ROOT)/FASTSIM/data/glauberPbPb.root","update");
1154   for(Int_t k=0; k<40; k++) {
1155     sprintf(almondName,"WAlmondFixedB%d",k);
1156     Double_t b = 0.25+k*0.5;
1157     printf(" b = %f\n",b);
1158     fgWAlmond->SetParameter(0,b);
1159
1160     fgWAlmond->Write(almondName);
1161   }
1162   ff->Close();
1163
1164   return;
1165 }
1166
1167 void AliFastGlauber::PlotAlmonds()
1168 {
1169   //
1170   // Plot almonds for some impact parameters
1171   //
1172
1173   TCanvas *c = new TCanvas("c","Almonds",0,0,500,500);
1174   gStyle->SetPalette(1,0);
1175   c->Divide(2,2);
1176   c->cd(1);
1177   fWAlmondFixedB[0].Draw("cont1");
1178   c->cd(2);
1179   fWAlmondFixedB[10].Draw("cont1");
1180   c->cd(3);
1181   fWAlmondFixedB[20].Draw("cont1");
1182   c->cd(4);
1183   fWAlmondFixedB[30].Draw("cont1");
1184
1185   return;
1186 }