]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - ZDC/AliZDCFragment.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Fix by Ruben
[u/mrichter/AliRoot.git] / ZDC / AliZDCFragment.cxx
... / ...
CommitLineData
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
33ClassImp(AliZDCFragment)
34
35int comp(const void *i,const void *j) {return *(int *)i - *(int *)j;}
36
37
38//_____________________________________________________________________________
39AliZDCFragment::AliZDCFragment():
40 fB(0),
41 fZbAverage(0),
42 fNimf(0),
43 fZmax(0),
44 fTau(0),
45 fNalpha(0),
46 fZtot(0),
47 fNtot(0)
48{
49 //
50 // Default constructor
51 //
52 for(Int_t i=0; i<=99; i++){
53 fZZ[i] = 0;
54 fNN[i] = 0;
55 }
56}
57
58//_____________________________________________________________________________
59AliZDCFragment::AliZDCFragment(Float_t b):
60 TNamed(" "," "),
61 fB(b),
62 fZbAverage(0),
63 fNimf(0),
64 fZmax(0),
65 fTau(0),
66 fNalpha(0),
67 fZtot(0),
68 fNtot(0)
69{
70 //
71 // Standard constructor
72 //
73 for(Int_t i=0; i<=99; i++){
74 fZZ[i] = 0;
75 fNN[i] = 0;
76 }
77
78}
79
80//_____________________________________________________________________________
81void 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 if(fB>kbMaxPb) fB = 2*kbMaxPb-fB;
117
118 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 Float_t zbU = bU*bU*TMath::Pi()/7.48;
127
128 // Rescale Zbound for Pb
129 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 if(fB>kbMaxPb){
136 fZbAverage = TMath::Exp(-kParamSkinPb[0]*(fB-kParamSkinPb[1]));
137 //printf(" b = %1.2f fm Z_bound %1.2f\n", fB, fZbAverage);
138 }
139 if(fZbAverage>kZProj) fZbAverage = kZProj;
140 Float_t zbNorm = fZbAverage/kZProj;
141 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 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);
149
150 // Add fluctuation: from Singh et al.
151 Float_t fluctNimf = kParamFluctNimf[0]+kParamFluctNimf[1]*zbNorm+
152 kParamFluctNimf[2]*TMath::Power(zbNorm,2)+kParamFluctNimf[3]
153 *TMath::Power(zbNorm,3);
154 Float_t xx = gRandom->Gaus(0.0,1.0);
155 fluctNimf = fluctNimf*xx;
156 fNimf = Int_t(averageNimf+fluctNimf);
157 Float_t y = gRandom->Rndm();
158 if(y < ((averageNimf+fluctNimf)-fNimf)) fNimf += 1;
159 if(fNimf ==0 && zbNorm>0.75) fNimf = 1;
160
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);
164
165 // Add fluctuation to mean value of Zmax (see Hubele)
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.;
170 Float_t xg = gRandom->Gaus(0.0,1.0);
171 fluctZmax = fluctZmax*xg;
172 fZmax = (averageZmax+fluctZmax);
173 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 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);
189 Float_t xxx = gRandom->Gaus(0.0,1.0);
190 fluctAlpha = fluctAlpha*xxx;
191 fNalpha = Int_t(averageAlpha+fluctAlpha);
192 Float_t yy = gRandom->Rndm();
193 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 if(bNorm<=0.9) {
203 // remove alpha from zbound to find zbound for fragments (Z>=3)
204 zbFrag = fZbAverage-fNalpha*2;
205 choice = 1;
206 }
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 if(zbFrag<=fZmax) {
216 if(fNimf>0 && zbFrag>=2){
217 fNimf = 1;
218 fZZ[0] = Int_t(zbFrag);
219 sumZ = zbFrag;
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
242 Float_t * zz = new Float_t[fNimf];
243 for(j=0; j<fNimf; j++){
244 zz[j] =0;
245 }
246 for(i=0; i<fNimf; i++){
247 zz[i] = Float_t(funTau->GetRandom());
248// printf("\n zz[%d] = %f \n",i,zz[i]);
249 }
250 delete funTau;
251
252 // Sorting vector in ascending order with C function QSORT
253 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 for(j=0; j<fNimf; j++){
262 fZZ[j] = Int_t (zz[j]*fZmax/zz[fNimf-1]);
263 if(fZZ[j]<3) fZZ[j] = 3;
264// printf("\n fZZ[%d] = %d \n",j,fZZ[j]);
265 }
266
267 delete[] zz;
268
269 // Check that the sum of the bound charges is not > than Zbound-Zalfa
270
271 for(Int_t ii=0; ii<fNimf; ii++){
272 sumZ += fZZ[ii];
273 }
274
275 Int_t k = 0;
276 if(sumZ>zbFrag){
277 for(i=0; i< fNimf; i++){
278 k += 1;
279 sumZ -= fZZ[i];
280 if(sumZ<=zbFrag){
281 fNimf -= (i+1);
282 break;
283 }
284 }
285 }
286 else {
287 if(choice == 1) return;
288 Int_t iDiff = Int_t((zbFrag-sumZ)/2);
289 if(iDiff<fNalpha){
290 fNalpha=iDiff;
291 return;
292 }
293 else{
294 return;
295 }
296 }
297
298 fNimf += k;
299 for(i=0; i<fNimf; i++){
300 fZZ[i] = fZZ[i+k];
301 }
302 fNimf -= k;
303
304 sumZ=0;
305 for(i=0; i<fNimf; i++){
306 sumZ += fZZ[i];
307 }
308
309}
310
311//_____________________________________________________________________________
312void 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 for(Int_t i=0; i<fNimf; i++) {
353 for(Int_t j=0; j<68; j++) {
354 iZ = kZIon[j];
355 if((fZZ[i]-iZ) == 0){
356 iA = Int_t(kAIon[j]/0.93149432+0.5);
357 fNN[i] = iA - iZ;
358 break;
359 }
360 else if((fZZ[i]-iZ) < 0){
361 fZZ[i] = kZIon[j-1];
362 iA = Int_t (kAIon[j-1]/0.93149432+0.5);
363 fNN[i] = iA - kZIon[j-1];
364 break;
365 }
366 }
367 fZtot += fZZ[i];
368 fNtot += fNN[i];
369 }
370
371
372}
373
374//_____________________________________________________________________________
375Float_t AliZDCFragment::DeuteronNumber()
376{
377 // Calculates the fraction of deuterum nucleus produced
378 //
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;
383}