cffcdabedf0a22b6650a0a644faa9227b8e110a6
[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 TF1*    AliFastGlauber::fWStaa    = NULL;   
36 TF1*    AliFastGlauber::fWSgeo    = NULL;   
37 TF1*    AliFastGlauber::fWSbinary = NULL;   
38 TF1*    AliFastGlauber::fWSN      = NULL;   
39 Float_t AliFastGlauber::fbMax     = 0.;
40
41 AliFastGlauber::AliFastGlauber()
42 {
43 //  Default Constructor
44 //
45 //  Defaults for Pb
46     SetWoodSaxonParameters(6.624, 0.549, 0.00, 7.69e-4);
47     SetHardCrossSection();
48     SetMaxImpact();
49 }
50
51 void AliFastGlauber::Init(Int_t mode)
52 {
53 // Initialisation
54 //
55 //  Wood-Saxon
56 //
57     fWSb = new TF1("WSb", WSb, 0, fbMax, 4);
58     fWSb->SetParameter(0, fWSr0);
59     fWSb->SetParameter(1, fWSd);
60     fWSb->SetParameter(2, fWSw);
61     fWSb->SetParameter(3, fWSn);
62
63     fWSbz = new TF2("WSbz", WSbz, 0, fbMax, 4);
64     fWSbz->SetParameter(0, fWSr0);
65     fWSbz->SetParameter(1, fWSd);
66     fWSbz->SetParameter(2, fWSw);
67     fWSbz->SetParameter(3, fWSn);
68
69     fWSz = new TF1("WSz", WSz, 0, fbMax, 5);
70     fWSz->SetParameter(0, fWSr0);
71     fWSz->SetParameter(1, fWSd);
72     fWSz->SetParameter(2, fWSw);
73     fWSz->SetParameter(3, fWSn);
74
75 //
76 //  Thickness
77 //
78     fWSta = new TF1("WSta", WSta, 0., fbMax, 0);
79     
80 //
81 //  Overlap Kernel
82 //
83     fWStarfi = new TF2("WStarfi", WStarfi, 0., fbMax, 0., TMath::Pi(), 1);
84     fWStarfi->SetParameter(0, 0.);     
85     fWStarfi->SetNpx(200);     
86     fWStarfi->SetNpy(20);     
87 //
88 //  Overlap
89 //
90     if (! mode) {
91         fWStaa = new TF1("WStaa", WStaa, 0., fbMax, 0);
92         fWStaa->SetNpx(100);
93     } else {
94         TFile* f = new TFile("$(ALICE_ROOT)/FASTSIM/data/glauberPbPb.root");
95         fWStaa = (TF1*) f->Get("WStaa");
96     }
97     
98 //
99 //  Geometrical Cross-Section
100 //
101     fWSgeo = new TF1("WSgeo", WSgeo, 0., fbMax, 0);
102     fWSgeo->SetNpx(100);
103 //
104 //  Hard cross section (~ binary collisions)
105 //
106     fWSbinary = new TF1("WSbinary", WSbinary, 0., fbMax, 1);
107     fWSbinary->SetParameter(0, fSigmaHard); // mb
108     fWSbinary->SetNpx(100);
109 //
110 // Hard collisions per event
111 //
112     fWSN = new TF1("WSN", WSN, 0., fbMax, 1);
113     fWSN->SetNpx(100);
114 }
115
116 void AliFastGlauber::DrawWSb()
117 {
118 //
119 //  Draw Wood-Saxon Nuclear Density Function
120 //
121     TCanvas *c1 = new TCanvas("c1","Wood Saxon",400,10,600,700);
122     c1->cd();
123     fWSb->Draw();
124 }
125
126 void AliFastGlauber::DrawOverlap()
127 {
128 //
129 //  Draw Overlap Function
130 //
131     TCanvas *c2 = new TCanvas("c2","Overlap",400,10,600,700);
132     c2->cd();
133     fWStaa->Draw();
134 }
135
136 void AliFastGlauber::DrawThickness()
137 {
138 //
139 //  Draw Thickness Function
140 //
141     TCanvas *c3 = new TCanvas("c3","Thickness",400,10,600,700);
142     c3->cd();
143     fWSta->Draw();
144 }
145
146 void AliFastGlauber::DrawGeo()
147 {
148 //
149 //  Draw Geometrical Cross-Section
150 //
151     TCanvas *c3 = new TCanvas("c3","Geometrical Cross-Section",400,10,600,700);
152     c3->cd();
153     fWSgeo->Draw();
154 }
155
156 void AliFastGlauber::DrawBinary()
157 {
158 //
159 //  Draw Wood-Saxon Nuclear Density Function
160 //
161     TCanvas *c4 = new TCanvas("c4","Binary Cross-Section",400,10,600,700);
162     c4->cd();
163     fWSbinary->Draw();
164 }
165
166 void AliFastGlauber::DrawN()
167 {
168 //
169 //  Draw Binaries per event
170 //
171     TCanvas *c5 = new TCanvas("c5","Binaries per event",400,10,600,700);
172     c5->cd();
173     fWSN->Draw();
174 }
175
176 void AliFastGlauber::DrawKernel()
177 {
178 //
179 //  Draw Kernel
180 //
181     TCanvas *c6 = new TCanvas("c6","Kernel",400,10,600,700);
182     c6->cd();
183     fWStarfi->Draw();
184 }
185
186 Double_t AliFastGlauber::WSb(Double_t* x, Double_t* par)
187 {
188 //
189 //  Woods-Saxon Parameterisation
190 //  as a function of radius
191 //
192     Double_t xx  = x[0];
193     Double_t r0  = par[0];
194     Double_t d   = par[1];
195     Double_t w   = par[2];
196     Double_t n   = par[3];
197     
198     Double_t y  = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
199     return y;
200 }
201
202 Double_t AliFastGlauber::WSbz(Double_t* x, Double_t* par)
203 {
204 //
205 //  Wood Saxon Parameterisation
206 //  as a function of z and  b
207 //
208     Double_t bb  = x[0];
209     Double_t zz  = x[1];
210     Double_t r0  = par[0];
211     Double_t d   = par[1];
212     Double_t w   = par[2];
213     Double_t n   = par[3];
214     Double_t xx  = TMath::Sqrt(bb*bb+zz*zz);
215     Double_t y  = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
216     return y;
217 }
218
219 Double_t AliFastGlauber::WSz(Double_t* x, Double_t* par)
220 {
221 //
222 //  Wood Saxon Parameterisation
223 //  as a function of z for fixed b
224 //
225     Double_t bb  = par[4];
226     Double_t zz  = x[0];
227     Double_t r0  = par[0];
228     Double_t d   = par[1];
229     Double_t w   = par[2];
230     Double_t n   = par[3];
231     Double_t xx  = TMath::Sqrt(bb*bb+zz*zz);
232     Double_t y  = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
233 //    if (y < 1.e-6) y = 0;
234     return y;
235 }
236
237 Double_t AliFastGlauber::WSta(Double_t* x, Double_t* /*par*/)
238 {
239 //
240 //  Thickness function 
241 //
242     Double_t b  = x[0];
243     fWSz->SetParameter(4, b);
244     Double_t y  = 2. * fWSz->Integral(0., fbMax);
245     return y;
246 }
247
248
249
250 Double_t AliFastGlauber::WStarfi(Double_t* x, Double_t* par)
251 {
252 //
253 //  Kernel for overlap function
254 //
255     Double_t b    = par[0];
256     Double_t r1   = x[0];
257     Double_t phi  = x[1];
258     Double_t r2   = TMath::Sqrt(r1 * r1 + b * b - 2. * r1 * b * TMath::Cos(phi)); 
259     Double_t y    = r1 * fWSta->Eval(r1) * fWSta->Eval(r2);
260     return y;
261 }
262
263
264 Double_t AliFastGlauber::WStaa(Double_t* x, Double_t* /*par*/)
265 {
266 //
267 //  Overlap function
268 //
269     Double_t b    = x[0];
270     fWStarfi->SetParameter(0, b);
271 /*
272     Double_t al[2];
273     Double_t bl[2];
274     al[0] = 0.;
275     al[1] = 0.;
276     bl[0] = 6.6;
277     bl[1] = TMath::Pi();
278     Double_t err;
279     
280     Double_t y =  2. * fWStarfi->IntegralMultiple(2, al, bl, 0.001, err);
281     printf("WStaa: %f %f %f\n", b, y, err);
282 */
283 //
284 //  MC Integration
285 //
286     Double_t y = 0;
287     for (Int_t i = 0; i < 100000; i++)
288     {
289         Double_t phi = TMath::Pi() * gRandom->Rndm();
290         Double_t b1  = fbMax       * gRandom->Rndm();   
291         y += fWStarfi->Eval(b1, phi);
292     }
293     y *= 2. * 0.1 *  208. * 208. * TMath::Pi() * fbMax / 100000.;
294     return y;
295 }
296
297 Double_t AliFastGlauber::WSgeo(Double_t* x, Double_t* /*par*/)
298 {
299 //
300 //  Geometrical Cross-Section
301 //
302     Double_t b    = x[0];
303     Double_t taa  = fWStaa->Eval(b);
304     const Double_t sigma = 55.6; // mbarn
305     
306     Double_t y    = 2. * TMath::Pi() * b * (1. - TMath::Exp(- sigma * taa)); // fm
307     return y;
308 }
309
310
311 Double_t AliFastGlauber::WSbinary(Double_t* x, Double_t* par)
312 {
313 //
314 //  Number of binary collisions
315 //
316     Double_t b     = x[0];
317     Double_t sigma = par[0];
318     Double_t taa   = fWStaa->Eval(b);
319     
320     Double_t y    = 2. * TMath::Pi() * b * sigma * taa; // fm
321     return y;
322 }
323
324 Double_t AliFastGlauber::WSN(Double_t* x, Double_t* /*par*/)
325 {
326 //
327 //  Number of hard processes per event
328 //
329     Double_t b     = x[0];
330     Double_t y     = fWSbinary->Eval(b)/fWSgeo->Eval(b);
331     return y;
332 }
333
334 void AliFastGlauber::SimulateTrigger(Int_t n)
335 {
336     //
337     //  Simulates Trigger
338     //
339     TH1F* mbtH = new TH1F("mbtH", "MB Trigger b-Distribution",   100, 0., 20.);
340     TH1F* hdtH = new TH1F("hdtH", "Hard Trigger b-Distribution", 100, 0., 20.);   
341     TH1F* mbmH = new TH1F("mbmH", "MB Trigger Multiplicity Distribution",   100, 0., 8000.);
342     TH1F* hdmH = new TH1F("hdmH", "Hard Trigger Multiplicity Distribution", 100, 0., 8000.);   
343
344     mbtH->SetXTitle("b [fm]");
345     hdtH->SetXTitle("b [fm]");    
346     mbmH->SetXTitle("Multiplicity");
347     hdmH->SetXTitle("Multiplicity");    
348
349     TCanvas *c0 = new TCanvas("c0","Trigger Simulation",400,10,600,700);    
350     c0->Divide(2,1);
351     TCanvas *c1 = new TCanvas("c1","Trigger Simulation",400,10,600,700);    
352     c1->Divide(1,2);
353
354     //
355     //
356     Init(1);
357     for (Int_t iev = 0; iev < n; iev++)
358     {
359         Float_t b, p, mult;
360         GetRandom(b, p, mult);
361         mbtH->Fill(b,1.);
362         hdtH->Fill(b, p);
363         mbmH->Fill(mult, 1.);
364         hdmH->Fill(mult, p);
365
366         c0->cd(1);
367         mbtH->Draw();
368         c0->cd(2);
369         hdtH->Draw();   
370         c0->Update();
371
372         c1->cd(1);
373         mbmH->Draw();
374         c1->cd(2);
375         hdmH->Draw();   
376         c1->Update();
377     }
378 }
379
380 void AliFastGlauber::GetRandom(Float_t& b, Float_t& p, Float_t& mult)
381 {
382     //
383     // Gives back a random impact parameter, hard trigger probability and multiplicity
384     //
385         b = fWSgeo->GetRandom();
386         Float_t mu = fWSN->Eval(b);
387         p = 1.-TMath::Exp(-mu);
388         mult = 6000./fWSN->Eval(1.) * mu;
389 }
390
391 void AliFastGlauber::GetRandom(Int_t& bin, Bool_t& hard)
392 {
393     //
394     // Gives back a random impact parameter bin, and hard trigger decission
395     //
396         Float_t b  = fWSgeo->GetRandom();
397         Float_t mu = fWSN->Eval(b) * fSigmaHard;
398         Float_t p  = 1.-TMath::Exp(-mu);
399         if (b < 5.) {
400             bin = 1;
401         } else if (b <  8.6) {
402             bin = 2;
403         } else if (b < 11.2) {
404             bin = 3;
405         } else if (b < 13.2) {
406             bin = 4;
407         } else if (b < 15.0) {
408             bin = 5;
409         } else {
410             bin = 6;
411         }
412         
413         hard = kFALSE;
414         
415         Float_t r = gRandom->Rndm();
416         
417         if (r < p) hard = kTRUE;
418 }
419
420
421 Float_t  AliFastGlauber::GetRandomImpactParameter(Float_t bmin, Float_t bmax)
422 {
423     //
424     // Gives back a random impact parameter in the range bmin .. bmax
425     //
426
427     Float_t b = -1.;
428     while(b < bmin || b > bmax)
429         b = fWSgeo->GetRandom();
430     return b;
431 }
432
433 Double_t AliFastGlauber::CrossSection(Double_t b1, Double_t b2)
434 {
435     //
436     // Return cross-section integrated from b1 to b2 
437     //
438     
439     return fWSgeo->Integral(b1, b2)/100.;
440 }
441
442 Double_t AliFastGlauber::FractionOfHardCrossSection(Double_t b1, Double_t b2)
443 {
444     //
445     // Return raction of hard cross-section integrated from b1 to b2 
446     //
447     
448     return fWSbinary->Integral(b1, b2)/fWSbinary->Integral(0., 100.);
449 }
450
451
452 Double_t AliFastGlauber::Binaries(Double_t b)
453 {
454     //
455     // Return number of binary collisions normalized to 1 at b=0
456     //
457     
458     return fWSN->Eval(b)/fWSN->Eval(0.);
459 }