1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 // ******************************************************************
19 // Class for nuclear fragments formation
21 // ******************************************************************
23 // --- Standard libraries
30 // --- AliRoot classes
31 #include "AliZDCFragment.h"
33 ClassImp(AliZDCFragment)
35 int comp(const void *i,const void *j) {return *(int *)i - *(int *)j;}
38 //_____________________________________________________________________________
39 AliZDCFragment::AliZDCFragment():
50 // Default constructor
54 //_____________________________________________________________________________
55 AliZDCFragment::AliZDCFragment(Float_t b):
67 // Standard constructor
69 for(Int_t i=0; i<=99; i++){
76 //_____________________________________________________________________________
77 void AliZDCFragment::GenerateIMF()
83 // Coefficients of polynomial for average number of IMF
84 const Float_t kParamNimf[5]={0.011236,1.8364,56.572,-116.24,58.289};
85 // Coefficients of polynomial for fluctuations on average number of IMF
86 const Float_t kParamFluctNimf[4]={-0.13176,2.9392,-5.2147,2.3092};
87 // Coefficients of polynomial for average maximum Z of fragments
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};
90 // Coefficients of polynomial for fluctuations on maximum Z of fragments
91 const Float_t kParamFluctZmax[5]={0.013782,-0.17282,1.5065,1.0654,-2.4317};
92 // Coefficients of polynomial for exponent tau of fragments Z distribution
93 const Float_t kParamTau[3]={6.7233,-15.85,13.047};
94 //Coefficients of polynomial for average number of alphas
95 const Float_t kParamNalpha[4]={-0.68554,39.605,-68.311,30.165};
96 // Coefficients of polynomial for fluctuations on average number of alphas
97 const Float_t kParamFluctNalpha[5]={0.283,6.2141,-17.113,17.394,-6.6084};
98 // Coefficients of function for Pb nucleus skin
99 const Float_t kParamSkinPb[2]={0.762408, 20.};
101 // Thickness of nuclear surface
102 //const Float_t kNuclearThick = 0.52;
103 // Maximum impact parameter for U [r0*A**(1/3)]
104 const Float_t kbMaxU = 14.87;
105 // Maximum impact parameter for Pb [r0*A**(1/3)]
106 //const Float_t kbMaxPb = 14.22+4*kNuclearThick;
107 const Float_t kbMaxPb = 14.22;
108 // Z of the projectile
109 const Float_t kZProj = 82.;
111 // From b(Pb) to b(U)
112 if(fB>kbMaxPb) fB = 2*kbMaxPb-fB;
114 Float_t bU = fB*kbMaxU/kbMaxPb;
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).
122 Float_t zbU = bU*bU*TMath::Pi()/7.48;
124 // Rescale Zbound for Pb
125 fZbAverage = kZProj/92.*zbU;
127 // Zbound is proportional to b**2 up to b < kbMaxPb-2*kNuclearThick
128 // and then it is an increasing exponential, imposing that at
129 // b=kbMaxPb-2kNuclearThick the two functions have the same derivative
130 //Float_t bCore = kbMaxPb-2*kNuclearThick;
132 fZbAverage = TMath::Exp(-kParamSkinPb[0]*(fB-kParamSkinPb[1]));
133 //printf(" b = %1.2f fm Z_bound %1.2f\n", fB, fZbAverage);
135 if(fZbAverage>kZProj) fZbAverage = kZProj;
136 Float_t zbNorm = fZbAverage/kZProj;
137 Float_t bNorm = fB/kbMaxPb;
139 // From Zbound to <Nimf>,<Zmax>,tau
140 // Polinomial fits to Aladin distribution
141 // --- A.Schuttauf et al, Nuc.Phys. A607 (1996) 457.
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);
146 // Add fluctuation: from Singh et al.
147 Float_t fluctNimf = kParamFluctNimf[0]+kParamFluctNimf[1]*zbNorm+
148 kParamFluctNimf[2]*TMath::Power(zbNorm,2)+kParamFluctNimf[3]
149 *TMath::Power(zbNorm,3);
150 Float_t xx = gRandom->Gaus(0.0,1.0);
151 fluctNimf = fluctNimf*xx;
152 fNimf = Int_t(averageNimf+fluctNimf);
153 Float_t y = gRandom->Rndm();
154 if(y < ((averageNimf+fluctNimf)-fNimf)) fNimf += 1;
155 if(fNimf ==0 && zbNorm>0.75) fNimf = 1;
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);
161 // Add fluctuation to mean value of Zmax (see Hubele)
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.;
166 Float_t xg = gRandom->Gaus(0.0,1.0);
167 fluctZmax = fluctZmax*xg;
168 fZmax = (averageZmax+fluctZmax);
169 if(fZmax>kZProj) fZmax = kZProj;
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);
176 // Find the number of alpha particles
177 // from Singh et al. : Pb+emulsion
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);
185 Float_t xxx = gRandom->Gaus(0.0,1.0);
186 fluctAlpha = fluctAlpha*xxx;
187 fNalpha = Int_t(averageAlpha+fluctAlpha);
188 Float_t yy = gRandom->Rndm();
189 if(yy < ((averageAlpha+fluctAlpha)-fNalpha)) fNalpha += 1;
192 // 1) for bNorm < 0.9 ==> first remove alphas, then fragments
193 // 2) for bNorm > 0.9 ==> first remove fragments, then alphas
196 Float_t zbFrag = 0, sumZ = 0.;
199 // remove alpha from zbound to find zbound for fragments (Z>=3)
200 zbFrag = fZbAverage-fNalpha*2;
207 // printf("\n choice = %d, fZbAverage = %f, zbFrag = %f \n", choice, fZbAverage, zbFrag);
210 // Check if zbFrag < fZmax
212 if(fNimf>0 && zbFrag>=2){
214 fZZ[0] = Int_t(zbFrag);
223 // Prepare the exponential charge distribution dN/dZ
233 TF1 *funTau = new TF1("funTau","1./(x**[0])",0.01,fZmax);
234 funTau->SetParameter(0,fTau);
236 // Extract randomly the charge of the fragments from the distribution
238 Float_t * zz = new Float_t[fNimf];
239 for(j=0; j<fNimf; j++){
242 for(i=0; i<fNimf; i++){
243 zz[i] = Float_t(funTau->GetRandom());
244 // printf("\n zz[%d] = %f \n",i,zz[i]);
248 // Sorting vector in ascending order with C function QSORT
249 qsort((void*)zz,fNimf,sizeof(Float_t),comp);
252 // for(Int_t i=0; i<fNimf; i++){
253 // printf("\n After sorting -> zz[%d] = %f \n",i,zz[i]);
256 // Rescale the maximum charge to fZmax
257 for(j=0; j<fNimf; j++){
258 fZZ[j] = Int_t (zz[j]*fZmax/zz[fNimf-1]);
259 if(fZZ[j]<3) fZZ[j] = 3;
260 // printf("\n fZZ[%d] = %d \n",j,fZZ[j]);
265 // Check that the sum of the bound charges is not > than Zbound-Zalfa
267 for(Int_t ii=0; ii<fNimf; ii++){
273 for(i=0; i< fNimf; i++){
283 if(choice == 1) return;
284 Int_t iDiff = Int_t((zbFrag-sumZ)/2);
295 for(i=0; i<fNimf; i++){
301 for(i=0; i<fNimf; i++){
307 //_____________________________________________________________________________
308 void AliZDCFragment::AttachNeutrons()
311 // Prepare nuclear fragment by attaching a suitable number of neutrons
313 const Float_t kAIon[68]={1.87612,2.80943,3.7284,5.60305,6.53536,
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,
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,
329 const Int_t kZIon[68]={1,1,2,3,3,
346 // printf("\n fNimf=%d\n",fNimf);
348 for(Int_t i=0; i<fNimf; i++) {
349 for(Int_t j=0; j<68; j++) {
351 if((fZZ[i]-iZ) == 0){
352 iA = Int_t(kAIon[j]/0.93149432+0.5);
356 else if((fZZ[i]-iZ) < 0){
358 iA = Int_t (kAIon[j-1]/0.93149432+0.5);
359 fNN[i] = iA - kZIon[j-1];
370 //_____________________________________________________________________________
371 Float_t AliZDCFragment::DeuteronNumber()
373 // Calculates the fraction of deuterum nucleus produced
375 Float_t deuteronProdPar[2] = {-0.068,0.0385};
376 Float_t deutNum = deuteronProdPar[0] + deuteronProdPar[1]*fB;
377 if(deutNum<0.) deutNum = 0.;