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