Fixing coverities
[u/mrichter/AliRoot.git] / ZDC / AliZDCFragment.cxx
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
17 // ******************************************************************
18 //
19 //      Class for nuclear fragments formation
20 //
21 // ******************************************************************
22
23 // --- Standard libraries
24 #include <stdlib.h>
25
26 // --- ROOT system
27 #include <TRandom.h>
28 #include <TF1.h>
29
30 // --- AliRoot classes
31 #include "AliZDCFragment.h"
32  
33 ClassImp(AliZDCFragment)
34    
35 int comp(const void *i,const void *j) {return *(int *)i - *(int *)j;}
36
37
38 //_____________________________________________________________________________
39 AliZDCFragment::AliZDCFragment():
40   fB(0),
41   fZbAverage(0),
42   fNimf(0),
43   fZmax(0),
44   fTau(0),
45   fNalpha(0),
46   fZtot(0),
47   fNtot(0)
48 {
49   //
50   // Default constructor
51   //
52   for(Int_t i=0; i<=99; i++){
53      fZZ[i] = 0;
54      fNN[i] = 0;
55   }
56 }
57
58 //_____________________________________________________________________________
59 AliZDCFragment::AliZDCFragment(Float_t b): 
60   TNamed(" "," "),
61   fB(b),
62   fZbAverage(0),
63   fNimf(0),
64   fZmax(0),
65   fTau(0),
66   fNalpha(0),
67   fZtot(0),
68   fNtot(0)
69 {
70   //
71   // Standard constructor
72   //
73   for(Int_t i=0; i<=99; i++){
74      fZZ[i] = 0;
75      fNN[i] = 0;
76   }
77   
78 }
79
80 //_____________________________________________________________________________
81 void AliZDCFragment::GenerateIMF()
82 {
83
84    // Loop variables
85   Int_t i,j;
86
87    // Coefficients of polynomial for average number of IMF
88    const Float_t  kParamNimf[5]={0.011236,1.8364,56.572,-116.24,58.289}; 
89    // Coefficients of polynomial for fluctuations on average number of IMF
90    const Float_t  kParamFluctNimf[4]={-0.13176,2.9392,-5.2147,2.3092}; 
91    // Coefficients of polynomial for average maximum Z of fragments
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}; 
94    // Coefficients of polynomial for fluctuations on maximum Z of fragments
95    const Float_t  kParamFluctZmax[5]={0.013782,-0.17282,1.5065,1.0654,-2.4317}; 
96    // Coefficients of polynomial for exponent tau of fragments Z distribution
97    const Float_t  kParamTau[3]={6.7233,-15.85,13.047};  
98    //Coefficients of polynomial for average number of alphas
99    const Float_t  kParamNalpha[4]={-0.68554,39.605,-68.311,30.165}; 
100    // Coefficients of polynomial for fluctuations on average number of alphas
101    const Float_t  kParamFluctNalpha[5]={0.283,6.2141,-17.113,17.394,-6.6084}; 
102    // Coefficients of function for Pb nucleus skin
103    const Float_t  kParamSkinPb[2]={0.762408, 20.};
104    
105    // Thickness of nuclear surface
106    //const Float_t  kNuclearThick = 0.52;
107    // Maximum impact parameter for U [r0*A**(1/3)]
108    const Float_t  kbMaxU = 14.87;
109    // Maximum impact parameter for Pb [r0*A**(1/3)]
110    //const Float_t  kbMaxPb = 14.22+4*kNuclearThick;
111    const Float_t  kbMaxPb = 14.22;
112    // Z of the projectile
113    const Float_t  kZProj = 82.;
114    
115    // From b(Pb) to b(U)
116    if(fB>kbMaxPb) fB = 2*kbMaxPb-fB;
117    
118    Float_t  bU = fB*kbMaxU/kbMaxPb;
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).
126    Float_t  zbU = bU*bU*TMath::Pi()/7.48;
127    
128    //  Rescale Zbound for Pb
129    fZbAverage = kZProj/92.*zbU;
130    
131    // Zbound is proportional to b**2 up to b < kbMaxPb-2*kNuclearThick
132    // and then it is an increasing exponential, imposing that at 
133    // b=kbMaxPb-2kNuclearThick the two functions have the same derivative
134    //Float_t bCore = kbMaxPb-2*kNuclearThick;
135    if(fB>kbMaxPb){
136      fZbAverage = TMath::Exp(-kParamSkinPb[0]*(fB-kParamSkinPb[1]));
137      //printf(" b = %1.2f fm   Z_bound %1.2f\n", fB, fZbAverage);
138    }
139    if(fZbAverage>kZProj) fZbAverage = kZProj;
140    Float_t zbNorm = fZbAverage/kZProj;
141    Float_t bNorm = fB/kbMaxPb;
142    
143    // From Zbound to <Nimf>,<Zmax>,tau
144    // Polinomial fits to Aladin distribution
145    // --- A.Schuttauf et al, Nuc.Phys. A607 (1996) 457.
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);
149    
150    // Add fluctuation: from Singh et al. 
151    Float_t fluctNimf = kParamFluctNimf[0]+kParamFluctNimf[1]*zbNorm+
152            kParamFluctNimf[2]*TMath::Power(zbNorm,2)+kParamFluctNimf[3]
153            *TMath::Power(zbNorm,3);
154    Float_t xx = gRandom->Gaus(0.0,1.0);
155    fluctNimf = fluctNimf*xx;
156    fNimf = Int_t(averageNimf+fluctNimf);
157    Float_t y = gRandom->Rndm();
158    if(y < ((averageNimf+fluctNimf)-fNimf)) fNimf += 1;
159    if(fNimf ==0 && zbNorm>0.75) fNimf = 1;
160    
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);
164    
165    // Add fluctuation to mean value of Zmax (see Hubele)
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.;
170    Float_t xg = gRandom->Gaus(0.0,1.0);
171    fluctZmax = fluctZmax*xg;
172    fZmax = (averageZmax+fluctZmax);
173    if(fZmax>kZProj) fZmax = kZProj;
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); 
179
180    // Find the number of alpha particles 
181    // from Singh et al. : Pb+emulsion
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);
189    Float_t xxx = gRandom->Gaus(0.0,1.0);
190    fluctAlpha = fluctAlpha*xxx;
191    fNalpha = Int_t(averageAlpha+fluctAlpha);
192    Float_t yy = gRandom->Rndm();
193    if(yy < ((averageAlpha+fluctAlpha)-fNalpha)) fNalpha += 1;
194
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
199    Int_t choice = 0;
200    Float_t zbFrag = 0, sumZ = 0.;
201
202    if(bNorm<=0.9) {
203    // remove alpha from zbound to find zbound for fragments  (Z>=3)
204      zbFrag = fZbAverage-fNalpha*2;
205      choice = 1;
206    }
207    else {
208      zbFrag = fZbAverage;
209      choice = 0;
210    }
211 //   printf("\n choice = %d, fZbAverage = %f, zbFrag = %f \n", choice, fZbAverage, zbFrag);
212    
213    
214    // Check if zbFrag < fZmax
215    if(zbFrag<=fZmax) {
216      if(fNimf>0 && zbFrag>=2){
217        fNimf = 1;
218        fZZ[0] = Int_t(zbFrag);
219        sumZ = zbFrag;
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  
242    Float_t * zz = new Float_t[fNimf];
243    for(j=0; j<fNimf; j++){
244       zz[j] =0;
245    }
246    for(i=0; i<fNimf; i++){
247       zz[i] = Float_t(funTau->GetRandom());
248 //      printf("\n      zz[%d] = %f \n",i,zz[i]);
249    }
250    delete funTau;
251    
252    // Sorting vector in ascending order with C function QSORT 
253    qsort((void*)zz,fNimf,sizeof(Float_t),comp);
254
255    
256 //   for(Int_t i=0; i<fNimf; i++){
257 //      printf("\n After sorting -> zz[%d] = %f \n",i,zz[i]);
258 //   }
259    
260    // Rescale the maximum charge to fZmax
261    for(j=0; j<fNimf; j++){
262      fZZ[j] = Int_t (zz[j]*fZmax/zz[fNimf-1]);
263      if(fZZ[j]<3) fZZ[j] = 3;
264 //     printf("\n       fZZ[%d] = %d \n",j,fZZ[j]);
265    }
266
267    delete[] zz;
268    
269    // Check that the sum of the bound charges is not > than Zbound-Zalfa
270    
271    for(Int_t ii=0; ii<fNimf; ii++){
272      sumZ += fZZ[ii];
273    }
274    
275    Int_t k = 0;
276    if(sumZ>zbFrag){
277      for(i=0; i< fNimf; i++){
278        k += 1;
279        sumZ -= fZZ[i];
280        if(sumZ<=zbFrag){
281          fNimf -= (i+1);
282          break;
283        }
284      }
285    }
286    else {
287      if(choice == 1) return;
288      Int_t iDiff = Int_t((zbFrag-sumZ)/2);
289      if(iDiff<fNalpha){
290        fNalpha=iDiff;
291        return;
292      }
293      else{
294        return;
295      }
296    }
297
298    fNimf += k;
299    for(i=0; i<fNimf; i++){
300      fZZ[i] = fZZ[i+k];
301    }
302    fNimf -= k;
303    
304    sumZ=0;
305    for(i=0; i<fNimf; i++){
306      sumZ += fZZ[i];
307    }
308    
309 }
310
311 //_____________________________________________________________________________
312 void AliZDCFragment::AttachNeutrons()
313 {
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,
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};
333    const Int_t kZIon[68]={1,1,2,3,3,
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;  
350 //   printf("\n fNimf=%d\n",fNimf);  
351
352    for(Int_t i=0; i<fNimf; i++) {
353       for(Int_t j=0; j<68; j++) {
354         iZ = kZIon[j];
355         if((fZZ[i]-iZ) == 0){
356           iA = Int_t(kAIon[j]/0.93149432+0.5);
357           fNN[i] = iA - iZ;
358           break;
359         }
360         else if((fZZ[i]-iZ) < 0){
361           fZZ[i] = kZIon[j-1];
362           iA = Int_t (kAIon[j-1]/0.93149432+0.5);
363           fNN[i] = iA - kZIon[j-1];
364           break;
365         }
366       }
367       fZtot += fZZ[i];
368       fNtot += fNN[i];
369    }                 
370    
371
372 }
373
374 //_____________________________________________________________________________
375 Float_t AliZDCFragment::DeuteronNumber()
376 {
377     // Calculates the fraction of deuterum nucleus produced
378     //
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;
383 }