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