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