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