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