doxy: consider comments after last line
[u/mrichter/AliRoot.git] / EVGEN / glauber.C
1 static TF1* f_wdsx_z;
2 static TF1* f_ta;
3 static TF2* f_ta_rfi;
4
5 void glauber()
6 {
7 //
8 //  Calculates some geometrical properties of PbPb collisions 
9 //  in the Glauber Model
10 //
11 //  Wood-Saxon nuclear density function
12 //
13     TCanvas *c1 = new TCanvas("c1","Wood Saxon",400,10,600,700);
14     TF1* f_wdsx = new TF1("f_wdsx", wdsx, 0, 15., 4);
15     f_wdsx->SetParameter(0,6.624);
16     f_wdsx->SetParameter(1,0.549);
17     f_wdsx->SetParameter(2,0.000);
18     f_wdsx->SetParameter(3,7.69e-4);
19     f_wdsx->Draw();
20 //
21 //  Wood Saxon-nuclear density (b-z)
22 //
23     TCanvas *c2 = new TCanvas("c2","Wood Saxon",400,10,600,700);
24     TF2* f_wdsx_bz = new TF2("f_wdsx_bz", wdsx_bz, 0, 15., 0., 15., 4);
25     f_wdsx_bz->SetParameter(0,6.624);
26     f_wdsx_bz->SetParameter(1,0.549);
27     f_wdsx_bz->SetParameter(2,0.000);
28     f_wdsx_bz->SetParameter(3,7.69e-4);
29     f_wdsx_bz->Draw();
30 //
31 //  Wood Saxon-nuclear density (z, for fixed b)
32 //
33     TCanvas *c3 = new TCanvas("c3","Wood Saxon",400,10,600,700);
34     f_wdsx_z = new TF1("f_wdsx_z", wdsx_z, 0, 15., 5);
35     f_wdsx_z->SetParameter(0,6.624);
36     f_wdsx_z->SetParameter(1,0.549);
37     f_wdsx_z->SetParameter(2,0.000);
38     f_wdsx_z->SetParameter(3,7.69e-4);
39     f_wdsx_z->SetParameter(4,0.);
40     f_wdsx_z->Draw();
41 //
42 //  Thickness function 
43 //
44     TCanvas *c4 = new TCanvas("c4","T_A",400,10,600,700);
45     f_ta = new TF1("f_ta", ta, 0, 15., 0);
46     f_ta->Draw();
47 //
48 //  Kernel of overlap function
49 //    
50
51     TCanvas *c5 = new TCanvas("c5","T_A",400,10,600,700);
52     f_ta_rfi = new TF2("f_ta_rfi", ta_rfi, 0, 15., 0., TMath::Pi(), 1);
53     f_ta_rfi->SetParameter(0,0.);
54     f_ta_rfi->Draw();
55 //
56 //  Overlap Function
57 //
58     TCanvas *c6 = new TCanvas("c6","T_AA",400,10,600,700);
59     TF1* f_taa = new TF1("f_taa", taa,0.,15., 0);
60     f_taa->Draw();
61
62 }
63
64
65
66 Double_t taa(Double_t* x, Double_t* dum)
67 {
68     printf("taa %f\n", x[0]);
69     
70     Double_t b    = x[0];
71     Double_t y    = 10./(208*208)*hijing->Profile((Float_t)b);
72     return y;
73 }
74
75
76 Double_t wdsx(Double_t* x, Double_t* par)
77 {
78 //
79 //  Wood Saxon Parameterisation
80 //  as a function of radius
81 //
82     Double_t xx  = x[0];
83     Double_t r0  = par[0];
84     Double_t d   = par[1];
85     Double_t w   = par[2];
86     Double_t n   = par[3];
87     
88     Double_t y  = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
89     return y;
90 }
91
92 Double_t wdsx_bz(Double_t* x, Double_t* par)
93 {
94 //
95 //  Wood Saxon Parameterisation
96 //  as a function of z and  b
97 //
98     Double_t bb  = x[0];
99     Double_t zz  = x[1];
100     Double_t r0  = par[0];
101     Double_t d   = par[1];
102     Double_t w   = par[2];
103     Double_t n   = par[3];
104     Double_t xx  = TMath::Sqrt(bb*bb+zz*zz);
105     Double_t y  = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
106     return y;
107 }
108
109 Double_t wdsx_z(Double_t* x, Double_t* par)
110 {
111 //
112 //  Wood Saxon Parameterisation
113 //  as a function of z for fixed b
114 //
115     Double_t bb  = par[4];
116     Double_t zz  = x[0];
117     Double_t r0  = par[0];
118     Double_t d   = par[1];
119     Double_t w   = par[2];
120     Double_t n   = par[3];
121     Double_t xx  = TMath::Sqrt(bb*bb+zz*zz);
122     Double_t y  = n * (1.+w*(xx/r0)*(xx/r0))/(1.+TMath::Exp((xx-r0)/d));
123     return y;
124 }
125
126 Double_t ta(Double_t* x, Double_t* par)
127 {
128 //
129 //  Thickness function 
130 //
131     Double_t b  = x[0];
132     f_wdsx_z->SetParameter(4,b);
133     Double_t y  = 2.*f_wdsx_z->Integral(0.,15.);
134     return y;
135 }
136
137 Double_t ta_rfi(Double_t* x, Double_t* par)
138 {
139 //
140 //  Kernel for overlap function
141 //
142     Double_t b    = par[0];
143     Double_t r1   = x[0];
144     Double_t phi  = x[1];
145     Double_t r2   = TMath::Sqrt(r1*r1+b*b-2.*r1*b*TMath::Cos(phi));
146     Double_t y    = r1*f_ta->Eval(r1)*f_ta->Eval(r2);
147 //    Double_t y    = r1*f_ta->Eval(r1);
148     return y;
149 }
150
151
152 Double_t taa(Double_t* x, Double_t* par)
153 {
154 //
155 //  Overlap function
156 //
157
158     
159     Double_t b    = x[0];
160     f_ta_rfi->SetParameter(0,b);
161     f_ta_rfi->Update();
162     
163     Double_t y    = 2.*
164         f_ta_rfi->Integral(0.,15., 0., TMath::Pi(), 0.001);
165     printf("taa %f %f\n", x[0], y);
166     return y;
167 }
168
169
170
171
172
173
174
175
176
177