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