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