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