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