Line data Source code
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 : // ******************************************************************
18 : //
19 : // Class for nuclear fragments formation
20 : //
21 : // ******************************************************************
22 :
23 : // --- Standard libraries
24 : #include <stdlib.h>
25 :
26 : // --- ROOT system
27 : #include <TRandom.h>
28 : #include <TF1.h>
29 :
30 : // --- AliRoot classes
31 : #include "AliZDCFragment.h"
32 :
33 12 : ClassImp(AliZDCFragment)
34 :
35 0 : int comp(const void *i,const void *j) {return *(int *)i - *(int *)j;}
36 :
37 :
38 : //_____________________________________________________________________________
39 0 : AliZDCFragment::AliZDCFragment():
40 0 : fB(0),
41 0 : fZbAverage(0),
42 0 : fNimf(0),
43 0 : fZmax(0),
44 0 : fTau(0),
45 0 : fNalpha(0),
46 0 : fZtot(0),
47 0 : fNtot(0)
48 0 : {
49 : //
50 : // Default constructor
51 : //
52 0 : for(Int_t i=0; i<=99; i++){
53 0 : fZZ[i] = 0;
54 0 : fNN[i] = 0;
55 : }
56 0 : }
57 :
58 : //_____________________________________________________________________________
59 : AliZDCFragment::AliZDCFragment(Float_t b):
60 0 : TNamed(" "," "),
61 0 : fB(b),
62 0 : fZbAverage(0),
63 0 : fNimf(0),
64 0 : fZmax(0),
65 0 : fTau(0),
66 0 : fNalpha(0),
67 0 : fZtot(0),
68 0 : fNtot(0)
69 0 : {
70 : //
71 : // Standard constructor
72 : //
73 0 : for(Int_t i=0; i<=99; i++){
74 0 : fZZ[i] = 0;
75 0 : fNN[i] = 0;
76 : }
77 :
78 0 : }
79 :
80 : //_____________________________________________________________________________
81 : void AliZDCFragment::GenerateIMF()
82 : {
83 :
84 : // Loop variables
85 : Int_t i,j;
86 :
87 : // Coefficients of polynomial for average number of IMF
88 : const Float_t kParamNimf[5]={0.011236,1.8364,56.572,-116.24,58.289};
89 : // Coefficients of polynomial for fluctuations on average number of IMF
90 : const Float_t kParamFluctNimf[4]={-0.13176,2.9392,-5.2147,2.3092};
91 : // Coefficients of polynomial for average maximum Z of fragments
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};
94 : // Coefficients of polynomial for fluctuations on maximum Z of fragments
95 : const Float_t kParamFluctZmax[5]={0.013782,-0.17282,1.5065,1.0654,-2.4317};
96 : // Coefficients of polynomial for exponent tau of fragments Z distribution
97 : const Float_t kParamTau[3]={6.7233,-15.85,13.047};
98 : //Coefficients of polynomial for average number of alphas
99 : const Float_t kParamNalpha[4]={-0.68554,39.605,-68.311,30.165};
100 : // Coefficients of polynomial for fluctuations on average number of alphas
101 : const Float_t kParamFluctNalpha[5]={0.283,6.2141,-17.113,17.394,-6.6084};
102 : // Coefficients of function for Pb nucleus skin
103 : const Float_t kParamSkinPb[2]={0.762408, 20.};
104 :
105 : // Thickness of nuclear surface
106 : //const Float_t kNuclearThick = 0.52;
107 : // Maximum impact parameter for U [r0*A**(1/3)]
108 : const Float_t kbMaxU = 14.87;
109 : // Maximum impact parameter for Pb [r0*A**(1/3)]
110 : //const Float_t kbMaxPb = 14.22+4*kNuclearThick;
111 : const Float_t kbMaxPb = 14.22;
112 : // Z of the projectile
113 : const Float_t kZProj = 82.;
114 :
115 : // From b(Pb) to b(U)
116 0 : if(fB>kbMaxPb) fB = 2*kbMaxPb-fB;
117 :
118 0 : Float_t bU = fB*kbMaxU/kbMaxPb;
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).
126 0 : Float_t zbU = bU*bU*TMath::Pi()/7.48;
127 :
128 : // Rescale Zbound for Pb
129 0 : fZbAverage = kZProj/92.*zbU;
130 :
131 : // Zbound is proportional to b**2 up to b < kbMaxPb-2*kNuclearThick
132 : // and then it is an increasing exponential, imposing that at
133 : // b=kbMaxPb-2kNuclearThick the two functions have the same derivative
134 : //Float_t bCore = kbMaxPb-2*kNuclearThick;
135 0 : if(fB>kbMaxPb){
136 0 : fZbAverage = TMath::Exp(-kParamSkinPb[0]*(fB-kParamSkinPb[1]));
137 : //printf(" b = %1.2f fm Z_bound %1.2f\n", fB, fZbAverage);
138 0 : }
139 0 : if(fZbAverage>kZProj) fZbAverage = kZProj;
140 0 : Float_t zbNorm = fZbAverage/kZProj;
141 0 : Float_t bNorm = fB/kbMaxPb;
142 :
143 : // From Zbound to <Nimf>,<Zmax>,tau
144 : // Polinomial fits to Aladin distribution
145 : // --- A.Schuttauf et al, Nuc.Phys. A607 (1996) 457.
146 0 : Float_t averageNimf = kParamNimf[0]+kParamNimf[1]*zbNorm+kParamNimf[2]*
147 0 : TMath::Power(zbNorm,2)+kParamNimf[3]*TMath::Power(zbNorm,3)+
148 0 : kParamNimf[4]*TMath::Power(zbNorm,4);
149 :
150 : // Add fluctuation: from Singh et al.
151 0 : Float_t fluctNimf = kParamFluctNimf[0]+kParamFluctNimf[1]*zbNorm+
152 0 : kParamFluctNimf[2]*TMath::Power(zbNorm,2)+kParamFluctNimf[3]
153 0 : *TMath::Power(zbNorm,3);
154 0 : Float_t xx = gRandom->Gaus(0.0,1.0);
155 0 : fluctNimf = fluctNimf*xx;
156 0 : fNimf = Int_t(averageNimf+fluctNimf);
157 0 : Float_t y = gRandom->Rndm();
158 0 : if(y < ((averageNimf+fluctNimf)-fNimf)) fNimf += 1;
159 0 : if(fNimf ==0 && zbNorm>0.75) fNimf = 1;
160 :
161 0 : Float_t averageZmax = kParamZmax[0]+kParamZmax[1]*zbNorm+kParamZmax[2]*
162 0 : TMath::Power(zbNorm,2)+kParamZmax[3]*TMath::Power(zbNorm,3);
163 0 : fTau = kParamTau[0]+kParamTau[1]*zbNorm+kParamTau[2]*TMath::Power(zbNorm,2);
164 :
165 : // Add fluctuation to mean value of Zmax (see Hubele)
166 0 : Float_t fluctZmax = kParamFluctZmax[0]+kParamFluctZmax[1]*zbNorm+
167 0 : kParamFluctZmax[2]*TMath::Power(zbNorm,2)+kParamFluctZmax[3]*
168 0 : TMath::Power(zbNorm,3)+kParamFluctZmax[4]*TMath::Power(zbNorm,4);
169 0 : fluctZmax = fluctZmax*kZProj/6.;
170 0 : Float_t xg = gRandom->Gaus(0.0,1.0);
171 0 : fluctZmax = fluctZmax*xg;
172 0 : fZmax = (averageZmax+fluctZmax);
173 0 : if(fZmax>kZProj) fZmax = kZProj;
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);
179 :
180 : // Find the number of alpha particles
181 : // from Singh et al. : Pb+emulsion
182 0 : Float_t averageAlpha = kParamNalpha[0]+kParamNalpha[1]*zbNorm+
183 0 : kParamNalpha[2]*TMath::Power(zbNorm,2)+kParamNalpha[3]*
184 0 : TMath::Power(zbNorm,3);
185 0 : Float_t fluctAlpha = kParamFluctNalpha[0]+kParamFluctNalpha[1]*
186 0 : zbNorm+kParamFluctNalpha[2]*TMath::Power(zbNorm,2)+
187 0 : kParamFluctNalpha[3]*TMath::Power(zbNorm,3)+
188 0 : kParamFluctNalpha[4]*TMath::Power(zbNorm,4);
189 0 : Float_t xxx = gRandom->Gaus(0.0,1.0);
190 0 : fluctAlpha = fluctAlpha*xxx;
191 0 : fNalpha = Int_t(averageAlpha+fluctAlpha);
192 0 : Float_t yy = gRandom->Rndm();
193 0 : if(yy < ((averageAlpha+fluctAlpha)-fNalpha)) fNalpha += 1;
194 :
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 :
199 : Int_t choice = 0;
200 : Float_t zbFrag = 0, sumZ = 0.;
201 :
202 0 : if(bNorm<=0.9) {
203 : // remove alpha from zbound to find zbound for fragments (Z>=3)
204 0 : zbFrag = fZbAverage-fNalpha*2;
205 : choice = 1;
206 0 : }
207 : else {
208 : zbFrag = fZbAverage;
209 : choice = 0;
210 : }
211 : // printf("\n choice = %d, fZbAverage = %f, zbFrag = %f \n", choice, fZbAverage, zbFrag);
212 :
213 :
214 : // Check if zbFrag < fZmax
215 0 : if(zbFrag<=fZmax) {
216 0 : if(fNimf>0 && zbFrag>=2){
217 0 : fNimf = 1;
218 0 : fZZ[0] = Int_t(zbFrag);
219 : sumZ = zbFrag;
220 0 : }
221 : else {
222 0 : fNimf = 0;
223 : }
224 0 : return;
225 : }
226 :
227 : // Prepare the exponential charge distribution dN/dZ
228 0 : if(fZmax <= 0.01) {
229 0 : fNimf = 0;
230 0 : return;
231 : }
232 0 : if(fNimf == 0) {
233 0 : fNimf = 0;
234 0 : return;
235 : }
236 :
237 0 : TF1 *funTau = new TF1("funTau","1./(x**[0])",0.01,fZmax);
238 0 : funTau->SetParameter(0,fTau);
239 :
240 : // Extract randomly the charge of the fragments from the distribution
241 :
242 0 : Float_t * zz = new Float_t[fNimf];
243 0 : for(j=0; j<fNimf; j++){
244 0 : zz[j] =0;
245 : }
246 0 : for(i=0; i<fNimf; i++){
247 0 : zz[i] = Float_t(funTau->GetRandom());
248 : // printf("\n zz[%d] = %f \n",i,zz[i]);
249 : }
250 0 : delete funTau;
251 :
252 : // Sorting vector in ascending order with C function QSORT
253 0 : qsort((void*)zz,fNimf,sizeof(Float_t),comp);
254 :
255 :
256 : // for(Int_t i=0; i<fNimf; i++){
257 : // printf("\n After sorting -> zz[%d] = %f \n",i,zz[i]);
258 : // }
259 :
260 : // Rescale the maximum charge to fZmax
261 0 : for(j=0; j<fNimf; j++){
262 0 : fZZ[j] = Int_t (zz[j]*fZmax/zz[fNimf-1]);
263 0 : if(fZZ[j]<3) fZZ[j] = 3;
264 : // printf("\n fZZ[%d] = %d \n",j,fZZ[j]);
265 : }
266 :
267 0 : delete[] zz;
268 :
269 : // Check that the sum of the bound charges is not > than Zbound-Zalfa
270 :
271 0 : for(Int_t ii=0; ii<fNimf; ii++){
272 0 : sumZ += fZZ[ii];
273 : }
274 :
275 : Int_t k = 0;
276 0 : if(sumZ>zbFrag){
277 0 : for(i=0; i< fNimf; i++){
278 0 : k += 1;
279 0 : sumZ -= fZZ[i];
280 0 : if(sumZ<=zbFrag){
281 0 : fNimf -= (i+1);
282 0 : break;
283 : }
284 : }
285 : }
286 : else {
287 0 : if(choice == 1) return;
288 0 : Int_t iDiff = Int_t((zbFrag-sumZ)/2);
289 0 : if(iDiff<fNalpha){
290 0 : fNalpha=iDiff;
291 0 : return;
292 : }
293 : else{
294 0 : return;
295 : }
296 : }
297 :
298 0 : fNimf += k;
299 0 : for(i=0; i<fNimf; i++){
300 0 : fZZ[i] = fZZ[i+k];
301 : }
302 0 : fNimf -= k;
303 :
304 : sumZ=0;
305 0 : for(i=0; i<fNimf; i++){
306 0 : sumZ += fZZ[i];
307 : }
308 :
309 0 : }
310 :
311 : //_____________________________________________________________________________
312 : void AliZDCFragment::AttachNeutrons()
313 : {
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,
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};
333 : const Int_t kZIon[68]={1,1,2,3,3,
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;
350 : // printf("\n fNimf=%d\n",fNimf);
351 :
352 0 : for(Int_t i=0; i<fNimf; i++) {
353 0 : for(Int_t j=0; j<68; j++) {
354 0 : iZ = kZIon[j];
355 0 : if((fZZ[i]-iZ) == 0){
356 0 : iA = Int_t(kAIon[j]/0.93149432+0.5);
357 0 : fNN[i] = iA - iZ;
358 0 : break;
359 : }
360 0 : else if((fZZ[i]-iZ) < 0){
361 0 : fZZ[i] = kZIon[j-1];
362 0 : iA = Int_t (kAIon[j-1]/0.93149432+0.5);
363 0 : fNN[i] = iA - kZIon[j-1];
364 0 : break;
365 : }
366 : }
367 0 : fZtot += fZZ[i];
368 0 : fNtot += fNN[i];
369 : }
370 :
371 :
372 0 : }
373 :
374 : //_____________________________________________________________________________
375 : Float_t AliZDCFragment::DeuteronNumber()
376 : {
377 : // Calculates the fraction of deuterum nucleus produced
378 : //
379 : Float_t deuteronProdPar[2] = {-0.068,0.0385};
380 0 : Float_t deutNum = deuteronProdPar[0] + deuteronProdPar[1]*fB;
381 0 : if(deutNum<0.) deutNum = 0.;
382 0 : return deutNum;
383 : }
|
