28e0901a |
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 | // --- Standard libraries |
17 | #include <stdlib.h> |
28e0901a |
18 | |
19 | // --- ROOT system |
20 | #include <TRandom.h> |
21 | #include <TF1.h> |
22 | |
23 | // --- AliRoot classes |
24 | #include "AliZDCFragment.h" |
28e0901a |
25 | |
5a881c97 |
26 | ClassImp(AliZDCFragment) |
27 | |
28 | int comp(const void *i,const void *j) {return *(int *)i - *(int *)j;} |
29 | |
30 | |
28e0901a |
31 | //_____________________________________________________________________________ |
32 | AliZDCFragment::AliZDCFragment() |
33 | { |
34 | // |
35 | // Default constructor |
36 | // |
37 | fB = 0; |
38 | } |
39 | |
40 | //_____________________________________________________________________________ |
41 | AliZDCFragment::AliZDCFragment(Float_t b) |
42 | : TNamed(" "," ") |
43 | { |
44 | // |
45 | // Standard constructor |
46 | // |
47 | fB = b; |
48 | fZbAverage = 0; |
49 | fNimf = 0; |
50 | fZmax = 0; |
51 | fTau = 0; |
28e0901a |
52 | for(Int_t i=0; i<=99; i++){ |
53 | fZZ[i] = 0; |
54 | fNN[i] = 0; |
55 | } |
5a881c97 |
56 | fNalpha = 0; |
57 | fZtot = 0; |
58 | fNtot = 0; |
28e0901a |
59 | |
60 | } |
61 | |
62 | //_____________________________________________________________________________ |
63 | void AliZDCFragment::GenerateIMF(Int_t* fZZ, Int_t &fNalpha) |
64 | { |
65 | // Coefficients of polynomial for average number of IMF |
5a881c97 |
66 | const Float_t ParamNimf[5]={0.011236,1.8364,56.572,-116.24,58.289}; |
28e0901a |
67 | // Coefficients of polynomial for fluctuations on average number of IMF |
5a881c97 |
68 | const Float_t ParamFluctNimf[4]={-0.13176,2.9392,-5.2147,2.3092}; |
28e0901a |
69 | // Coefficients of polynomial for average maximum Z of fragments |
5a881c97 |
70 | const Float_t ParamZmax[4]={0.16899,14.203,-2.8284,65.036}; |
28e0901a |
71 | // Coefficients of polynomial for fluctuations on maximum Z of fragments |
5a881c97 |
72 | const Float_t ParamFluctZmax[5]={0.013782,-0.17282,1.5065,1.0654,-2.4317}; |
28e0901a |
73 | // Coefficients of polynomial for exponent tau of fragments Z distribution |
5a881c97 |
74 | const Float_t ParamTau[3]={6.7233,-15.85,13.047}; |
28e0901a |
75 | //Coefficients of polynomial for average number of alphas |
5a881c97 |
76 | const Float_t ParamNalpha[4]={-0.68554,39.605,-68.311,30.165}; |
28e0901a |
77 | // Coefficients of polynomial for fluctuations on average number of alphas |
5a881c97 |
78 | const Float_t ParamFluctNalpha[5]={0.283,6.2141,-17.113,17.394,-6.6084}; |
28e0901a |
79 | // Coefficients of function for Pb nucleus skin |
5a881c97 |
80 | const Float_t ParamSkinPb[2]={0.93,11.05}; |
28e0901a |
81 | |
82 | // Thickness of nuclear surface |
5a881c97 |
83 | const Float_t NuclearThick = 0.52; |
28e0901a |
84 | // Maximum impact parameter for U [r0*A**(1/3)] |
5a881c97 |
85 | const Float_t bMaxU = 14.87; |
28e0901a |
86 | // Maximum impact parameter for Pb [r0*A**(1/3)] |
5a881c97 |
87 | const Float_t bMaxPb = 14.22; |
28e0901a |
88 | // Z of the projectile |
5a881c97 |
89 | const Float_t ZProj = 82.; |
28e0901a |
90 | |
91 | // From b(Pb) to b(U) |
92 | Float_t bU = fB*bMaxU/bMaxPb; |
93 | |
94 | // From b(U) to Zbound(U) |
95 | // --- A.Schuttauf et al, Nuc.Phys. A607 (1996) 457 --------------- |
96 | // From geometrical consideration and from dsigma/dZbound for U+U, |
97 | // which is approx. constant, the constant value is found |
98 | // integrating the nucleus cross surface from 0 to bmax=R1+R2 where |
99 | // R = 1.2*A**(1/3). This value has been measured in Aladin (U+U). |
100 | Float_t ZbU = bU*bU*TMath::Pi()/7.48; |
101 | |
102 | // Rescale Zbound for Pb |
103 | fZbAverage = ZProj/92.*ZbU; |
104 | |
105 | // Zbound is proportional to b**2 up to b < bMaxPb-2*NuclearThick |
106 | // and then it is an increasing exponential, imposing that at |
107 | // b=bMaxPb-2NuclearThick the two functions have the same derivative |
108 | Float_t bCore = bMaxPb-2*NuclearThick; |
109 | if(fB>bCore){ |
110 | fZbAverage=ZProj*(1.-TMath::Exp(-ParamSkinPb[0]*(fB-ParamSkinPb[1]))); |
111 | } |
112 | if(fZbAverage>ZProj) fZbAverage = ZProj; |
113 | Float_t ZbNorm = fZbAverage/ZProj; |
114 | Float_t bNorm = fB/bMaxPb; |
115 | |
116 | // From Zbound to <Nimf>,<Zmax>,tau |
117 | // Polinomial fits to Aladin distribution |
118 | // --- A.Schuttauf et al, Nuc.Phys. A607 (1996) 457. |
119 | Float_t AverageNimf = ParamNimf[0]+ParamNimf[1]*ZbNorm+ParamNimf[2]* |
120 | TMath::Power(ZbNorm,2)+ParamNimf[3]*TMath::Power(ZbNorm,3)+ |
121 | ParamNimf[4]*TMath::Power(ZbNorm,4); |
122 | |
123 | // Add fluctuation: from Singh et al. |
124 | Float_t FluctNimf = ParamFluctNimf[0]+ParamFluctNimf[1]*ZbNorm+ |
125 | ParamFluctNimf[2]*TMath::Power(ZbNorm,2)+ParamFluctNimf[3] |
126 | *TMath::Power(ZbNorm,3); |
127 | Float_t xx = gRandom->Gaus(0.0,1.0); |
128 | FluctNimf = FluctNimf*xx; |
129 | fNimf = Int_t(AverageNimf+FluctNimf); |
130 | Float_t y = gRandom->Rndm(); |
131 | if(y < ((AverageNimf+FluctNimf)-fNimf)) fNimf += 1; |
132 | if(fNimf ==0 && ZbNorm>0.75) fNimf = 1; |
28e0901a |
133 | |
134 | Float_t AverageZmax = ParamZmax[0]+ParamZmax[1]*ZbNorm+ParamZmax[2]* |
135 | TMath::Power(ZbNorm,2)+ParamZmax[3]*TMath::Power(ZbNorm,3); |
136 | fTau = ParamTau[0]+ParamTau[1]*ZbNorm+ParamTau[2]*TMath::Power(ZbNorm,2); |
137 | |
138 | // Add fluctuation to mean value of Zmax (see Hubele) |
139 | Float_t FluctZmax = ParamFluctZmax[0]+ParamFluctZmax[1]*ZbNorm+ |
140 | ParamFluctZmax[2]*TMath::Power(ZbNorm,2)+ParamFluctZmax[3]* |
141 | TMath::Power(ZbNorm,3)+ParamFluctZmax[4]*TMath::Power(ZbNorm,4); |
142 | FluctZmax = FluctZmax*ZProj/6.; |
143 | Float_t xg = gRandom->Gaus(0.0,1.0); |
144 | FluctZmax = FluctZmax*xg; |
145 | fZmax = AverageZmax+FluctZmax; |
146 | if(fZmax>ZProj) fZmax = ZProj; |
5a881c97 |
147 | |
148 | // printf("\n\n ------------------------------------------------------------"); |
149 | // printf("\n Generation of nuclear fragments\n"); |
150 | // printf("\n fNimf = %d\n", fNimf); |
151 | // printf("\n fZmax = %f\n", fZmax); |
28e0901a |
152 | |
153 | // Find the number of alpha particles |
154 | // from Singh et al. : Pb+emulsion |
155 | Float_t AverageAlpha = ParamNalpha[0]+ParamNalpha[1]*ZbNorm+ |
156 | ParamNalpha[2]*TMath::Power(ZbNorm,2)+ParamNalpha[3]* |
157 | TMath::Power(ZbNorm,3); |
158 | Float_t FluctAlpha = ParamFluctNalpha[0]+ParamFluctNalpha[1]* |
159 | ZbNorm+ParamFluctNalpha[2]*TMath::Power(ZbNorm,2)+ |
160 | ParamFluctNalpha[3]*TMath::Power(ZbNorm,3)+ |
161 | ParamFluctNalpha[4]*TMath::Power(ZbNorm,4); |
162 | Float_t xxx = gRandom->Gaus(0.0,1.0); |
163 | FluctAlpha = FluctAlpha*xxx; |
164 | fNalpha = Int_t(AverageAlpha+FluctAlpha); |
165 | Float_t yy = gRandom->Rndm(); |
166 | if(yy < ((AverageAlpha+FluctAlpha)-fNalpha)) fNalpha += 1; |
167 | |
28e0901a |
168 | // 2 possibilities: |
169 | // 1) for bNorm < 0.9 ==> first remove alphas, then fragments |
170 | // 2) for bNorm > 0.9 ==> first remove fragments, then alphas |
171 | |
5a881c97 |
172 | Int_t Choice = 0; |
173 | Float_t ZbFrag = 0, SumZ = 0.; |
174 | |
28e0901a |
175 | if(bNorm<=0.9) { |
176 | // remove alpha from zbound to find zbound for fragments (Z>=3) |
177 | ZbFrag = fZbAverage-fNalpha*2; |
178 | Choice = 1; |
179 | } |
180 | else { |
181 | ZbFrag = fZbAverage; |
182 | Choice = 0; |
183 | } |
184 | // printf("\n Choice = %d, fZbAverage = %f, ZbFrag = %f \n", Choice, fZbAverage, ZbFrag); |
185 | |
186 | |
187 | // Check if ZbFrag < fZmax |
188 | if(ZbFrag<=fZmax) { |
189 | if(fNimf>0 && ZbFrag>=2){ |
190 | fNimf = 1; |
191 | fZZ[0] = Int_t(ZbFrag); |
192 | SumZ = ZbFrag; |
193 | } |
194 | else { |
195 | fNimf = 0; |
196 | } |
197 | return; |
198 | } |
199 | |
200 | // Prepare the exponential charge distribution dN/dZ |
201 | if(fZmax <= 0.01) { |
202 | fNimf = 0; |
203 | return; |
204 | } |
205 | if(fNimf == 0) { |
206 | fNimf = 0; |
207 | return; |
208 | } |
209 | |
210 | TF1 *funTau = new TF1("funTau","1./(x**[0])",0.01,fZmax); |
211 | funTau->SetParameter(0,fTau); |
212 | |
213 | // Extract randomly the charge of the fragments from the distribution |
214 | |
4b016189 |
215 | Float_t * zz = new Float_t[fNimf]; |
28e0901a |
216 | for(Int_t j=0; j<fNimf; j++){ |
217 | zz[j] =0; |
218 | } |
219 | for(Int_t i=0; i<fNimf; i++){ |
220 | zz[i] = Float_t(funTau->GetRandom()); |
221 | // printf("\n zz[%d] = %f \n",i,zz[i]); |
222 | } |
223 | delete funTau; |
224 | |
5a881c97 |
225 | // Sorting vector in ascending order with C function QSORT |
4b016189 |
226 | qsort((void*)zz,fNimf,sizeof(Float_t),comp); |
28e0901a |
227 | |
5a881c97 |
228 | |
229 | // for(Int_t i=0; i<fNimf; i++){ |
28e0901a |
230 | // printf("\n After sorting -> zz[%d] = %f \n",i,zz[i]); |
5a881c97 |
231 | // } |
28e0901a |
232 | |
233 | // Rescale the maximum charge to fZmax |
234 | for(Int_t j=0; j<fNimf; j++){ |
235 | fZZ[j] = Int_t (zz[j]*fZmax/zz[fNimf-1]); |
236 | if(fZZ[j]<3) fZZ[j] = 3; |
5a881c97 |
237 | // printf("\n fZZ[%d] = %d \n",j,fZZ[j]); |
28e0901a |
238 | } |
4b016189 |
239 | |
240 | delete[] zz; |
28e0901a |
241 | |
5a881c97 |
242 | // Check that the sum of the bound charges is not > than Zbound-Zalfa |
28e0901a |
243 | |
244 | for(Int_t ii=0; ii<fNimf; ii++){ |
245 | SumZ += fZZ[ii]; |
246 | } |
247 | |
248 | Int_t k = 0; |
249 | if(SumZ>ZbFrag){ |
250 | for(Int_t i=0; i< fNimf; i++){ |
251 | k += 1; |
252 | SumZ -= fZZ[i]; |
253 | if(SumZ<=ZbFrag){ |
254 | fNimf -= (i+1); |
255 | break; |
256 | } |
257 | } |
258 | } |
259 | else { |
260 | if(Choice == 1) return; |
261 | Int_t iDiff = Int_t((ZbFrag-SumZ)/2); |
262 | if(iDiff<fNalpha){ |
263 | fNalpha=iDiff; |
264 | return; |
265 | } |
266 | else{ |
267 | return; |
268 | } |
269 | } |
270 | |
271 | fNimf += k; |
272 | for(Int_t i=0; i<fNimf; i++){ |
273 | fZZ[i] = fZZ[i+k]; |
274 | } |
275 | fNimf -= k; |
276 | |
277 | SumZ=0; |
278 | for(Int_t i=0; i<fNimf; i++){ |
279 | SumZ += fZZ[i]; |
280 | } |
28e0901a |
281 | |
282 | } |
283 | |
28e0901a |
284 | //_____________________________________________________________________________ |
285 | void AliZDCFragment::AttachNeutrons(Int_t *fZZ, Int_t *fNN, Int_t &fZtot,Int_t &fNtot) |
286 | { |
5a881c97 |
287 | const Float_t AIon[68]={1.87612,2.80943,3.7284,5.60305,6.53536, |
28e0901a |
288 | 6.53622,8.39479,9.32699,10.2551,11.17793, |
289 | 13.04378,14.89917,17.6969,18.62284,21.41483, |
290 | 22.34193,25.13314,26.06034,28.85188,29.7818, |
291 | 32.57328,33.50356,36.29447,37.22492,41.87617, |
292 | 44.66324,47.45401,48.38228,51.17447,52.10307, |
293 | 54.89593,53.96644,58.61856,59.54963,68.85715, |
294 | 74.44178,78.16309,81.88358,83.74571,91.19832, |
295 | 98.64997,106.10997,111.68821,122.86796, |
296 | 128.45793, |
297 | 130.32111,141.51236, |
298 | 141.55,146.477,148.033,152.699,153.631, |
299 | 155.802,157.357,162.022,162.984,166.2624, |
300 | 168.554,171.349,173.4536,177.198,179.0518, |
301 | 180.675,183.473,188.1345,190.77,193.729, |
302 | 221.74295}; |
5a881c97 |
303 | const Int_t ZIon[68]={1,1,2,3,3, |
28e0901a |
304 | 4,4,5,5,6, |
305 | 7,8,9,10,11, |
306 | 12,13,14,15,16, |
307 | 17,18,19,20,21, |
308 | 22,23,24,25,26, |
309 | 27,28,29,30,32, |
310 | 34,36,38,40,42, |
311 | 46,48,50,54,56, |
312 | 58,62, |
313 | 63,64,65,66,67, |
314 | 68,69,70,71,72, |
315 | 73,74,75,76,77, |
316 | 78,79,80,81,82, |
317 | 92}; |
318 | |
319 | Int_t iZ, iA; |
5a881c97 |
320 | // printf("\n fNimf=%d\n",fNimf); |
321 | |
28e0901a |
322 | for(Int_t i=0; i<fNimf; i++) { |
323 | for(Int_t j=0; j<68; j++) { |
324 | iZ = ZIon[j]; |
325 | if((fZZ[i]-iZ) == 0){ |
326 | iA = Int_t(AIon[j]/0.93149432+0.5); |
327 | fNN[i] = iA - iZ; |
28e0901a |
328 | break; |
329 | } |
330 | else if((fZZ[i]-iZ) < 0){ |
331 | fZZ[i] = ZIon[j-1]; |
332 | iA = Int_t (AIon[j-1]/0.93149432+0.5); |
333 | fNN[i] = iA - ZIon[j-1]; |
28e0901a |
334 | break; |
335 | } |
336 | } |
337 | fZtot += fZZ[i]; |
338 | fNtot += fNN[i]; |
339 | } |
340 | |
341 | |
342 | } |