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