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