fixed bug in AliSTARTTrigger
[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 // --- Standard libraries
17 #include <stdlib.h>
18
19 // --- ROOT system
20 #include <TRandom.h>
21 #include <TF1.h>
22
23 // --- AliRoot classes
24 #include "AliZDCFragment.h"
25  
26 ClassImp(AliZDCFragment)
27    
28 int comp(const void *i,const void *j) {return *(int *)i - *(int *)j;}
29
30
31 //_____________________________________________________________________________
32 AliZDCFragment::AliZDCFragment()
33 {
34   //
35   // Default constructor
36   //
37   fB = 0;
38 }
39
40 //_____________________________________________________________________________
41 AliZDCFragment::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;
52   for(Int_t i=0; i<=99; i++){
53      fZZ[i] = 0;
54      fNN[i] = 0;
55   }
56   fNalpha = 0;
57   fZtot = 0;
58   fNtot = 0;
59   
60 }
61
62 //_____________________________________________________________________________
63 void AliZDCFragment::GenerateIMF(Int_t* fZZ, Int_t &fNalpha)
64 {
65
66    // Loop variables
67   Int_t i,j;
68
69    // Coefficients of polynomial for average number of IMF
70    const Float_t  kParamNimf[5]={0.011236,1.8364,56.572,-116.24,58.289}; 
71    // Coefficients of polynomial for fluctuations on average number of IMF
72    const Float_t  kParamFluctNimf[4]={-0.13176,2.9392,-5.2147,2.3092}; 
73    // Coefficients of polynomial for average maximum Z of fragments
74    const Float_t  kParamZmax[4]={0.16899,14.203,-2.8284,65.036}; 
75    // Coefficients of polynomial for fluctuations on maximum Z of fragments
76    const Float_t  kParamFluctZmax[5]={0.013782,-0.17282,1.5065,1.0654,-2.4317}; 
77    // Coefficients of polynomial for exponent tau of fragments Z distribution
78    const Float_t  kParamTau[3]={6.7233,-15.85,13.047};  
79    //Coefficients of polynomial for average number of alphas
80    const Float_t  kParamNalpha[4]={-0.68554,39.605,-68.311,30.165}; 
81    // Coefficients of polynomial for fluctuations on average number of alphas
82    const Float_t  kParamFluctNalpha[5]={0.283,6.2141,-17.113,17.394,-6.6084}; 
83    // Coefficients of function for Pb nucleus skin
84    const Float_t  kParamSkinPb[2]={0.93,11.05};
85    
86    // Thickness of nuclear surface
87    const Float_t  kNuclearThick = 0.52;
88    // Maximum impact parameter for U [r0*A**(1/3)]
89    const Float_t  kbMaxU = 14.87;
90    // Maximum impact parameter for Pb [r0*A**(1/3)]
91    const Float_t  kbMaxPb = 14.22;
92    // Z of the projectile
93    const Float_t  kZProj = 82.;
94    
95    // From b(Pb) to b(U)
96    Float_t  bU = fB*kbMaxU/kbMaxPb;
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).
104    Float_t  zbU = bU*bU*TMath::Pi()/7.48;
105    
106    //  Rescale Zbound for Pb
107    fZbAverage = kZProj/92.*zbU;
108    
109    // Zbound is proportional to b**2 up to b < kbMaxPb-2*kNuclearThick
110    // and then it is an increasing exponential, imposing that at 
111    // b=kbMaxPb-2kNuclearThick the two functions have the same derivative
112    Float_t bCore = kbMaxPb-2*kNuclearThick;
113    if(fB>bCore){
114      fZbAverage=kZProj*(1.-TMath::Exp(-kParamSkinPb[0]*(fB-kParamSkinPb[1])));
115    }
116    if(fZbAverage>kZProj) fZbAverage = kZProj;
117    Float_t zbNorm = fZbAverage/kZProj;
118    Float_t bNorm = fB/kbMaxPb;
119    
120    // From Zbound to <Nimf>,<Zmax>,tau
121    // Polinomial fits to Aladin distribution
122    // --- A.Schuttauf et al, Nuc.Phys. A607 (1996) 457.
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);
126    
127    // Add fluctuation: from Singh et al. 
128    Float_t fluctNimf = kParamFluctNimf[0]+kParamFluctNimf[1]*zbNorm+
129            kParamFluctNimf[2]*TMath::Power(zbNorm,2)+kParamFluctNimf[3]
130            *TMath::Power(zbNorm,3);
131    Float_t xx = gRandom->Gaus(0.0,1.0);
132    fluctNimf = fluctNimf*xx;
133    fNimf = Int_t(averageNimf+fluctNimf);
134    Float_t y = gRandom->Rndm();
135    if(y < ((averageNimf+fluctNimf)-fNimf)) fNimf += 1;
136    if(fNimf ==0 && zbNorm>0.75) fNimf = 1;
137    
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);
141    
142    // Add fluctuation to mean value of Zmax (see Hubele)
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.;
147    Float_t xg = gRandom->Gaus(0.0,1.0);
148    fluctZmax = fluctZmax*xg;
149    fZmax = averageZmax+fluctZmax;
150    if(fZmax>kZProj) fZmax = kZProj;
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); 
156
157    // Find the number of alpha particles 
158    // from Singh et al. : Pb+emulsion
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);
166    Float_t xxx = gRandom->Gaus(0.0,1.0);
167    fluctAlpha = fluctAlpha*xxx;
168    fNalpha = Int_t(averageAlpha+fluctAlpha);
169    Float_t yy = gRandom->Rndm();
170    if(yy < ((averageAlpha+fluctAlpha)-fNalpha)) fNalpha += 1;
171
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
176    Int_t choice = 0;
177    Float_t zbFrag = 0, sumZ = 0.;
178
179    if(bNorm<=0.9) {
180    // remove alpha from zbound to find zbound for fragments  (Z>=3)
181      zbFrag = fZbAverage-fNalpha*2;
182      choice = 1;
183    }
184    else {
185      zbFrag = fZbAverage;
186      choice = 0;
187    }
188 //   printf("\n choice = %d, fZbAverage = %f, zbFrag = %f \n", choice, fZbAverage, zbFrag);
189    
190    
191    // Check if zbFrag < fZmax
192    if(zbFrag<=fZmax) {
193      if(fNimf>0 && zbFrag>=2){
194        fNimf = 1;
195        fZZ[0] = Int_t(zbFrag);
196        sumZ = zbFrag;
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  
219    Float_t * zz = new Float_t[fNimf];
220    for(j=0; j<fNimf; j++){
221       zz[j] =0;
222    }
223    for(i=0; i<fNimf; i++){
224       zz[i] = Float_t(funTau->GetRandom());
225 //      printf("\n      zz[%d] = %f \n",i,zz[i]);
226    }
227    delete funTau;
228    
229    // Sorting vector in ascending order with C function QSORT 
230    qsort((void*)zz,fNimf,sizeof(Float_t),comp);
231
232    
233 //   for(Int_t i=0; i<fNimf; i++){
234 //      printf("\n After sorting -> zz[%d] = %f \n",i,zz[i]);
235 //   }
236    
237    // Rescale the maximum charge to fZmax
238    for(j=0; j<fNimf; j++){
239      fZZ[j] = Int_t (zz[j]*fZmax/zz[fNimf-1]);
240      if(fZZ[j]<3) fZZ[j] = 3;
241 //     printf("\n       fZZ[%d] = %d \n",j,fZZ[j]);
242    }
243
244    delete[] zz;
245    
246    // Check that the sum of the bound charges is not > than Zbound-Zalfa
247    
248    for(Int_t ii=0; ii<fNimf; ii++){
249      sumZ += fZZ[ii];
250    }
251    
252    Int_t k = 0;
253    if(sumZ>zbFrag){
254      for(i=0; i< fNimf; i++){
255        k += 1;
256        sumZ -= fZZ[i];
257        if(sumZ<=zbFrag){
258          fNimf -= (i+1);
259          break;
260        }
261      }
262    }
263    else {
264      if(choice == 1) return;
265      Int_t iDiff = Int_t((zbFrag-sumZ)/2);
266      if(iDiff<fNalpha){
267        fNalpha=iDiff;
268        return;
269      }
270      else{
271        return;
272      }
273    }
274
275    fNimf += k;
276    for(i=0; i<fNimf; i++){
277      fZZ[i] = fZZ[i+k];
278    }
279    fNimf -= k;
280    
281    sumZ=0;
282    for(i=0; i<fNimf; i++){
283      sumZ += fZZ[i];
284    }
285    
286 }
287
288 //_____________________________________________________________________________
289 void AliZDCFragment::AttachNeutrons(Int_t *fZZ, Int_t *fNN, Int_t &fZtot,Int_t &fNtot)
290 {
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,
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};
310    const Int_t kZIon[68]={1,1,2,3,3,
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;  
327 //   printf("\n fNimf=%d\n",fNimf);  
328
329    for(Int_t i=0; i<fNimf; i++) {
330       for(Int_t j=0; j<68; j++) {
331         iZ = kZIon[j];
332         if((fZZ[i]-iZ) == 0){
333           iA = Int_t(kAIon[j]/0.93149432+0.5);
334           fNN[i] = iA - iZ;
335           break;
336         }
337         else if((fZZ[i]-iZ) < 0){
338           fZZ[i] = kZIon[j-1];
339           iA = Int_t (kAIon[j-1]/0.93149432+0.5);
340           fNN[i] = iA - kZIon[j-1];
341           break;
342         }
343       }
344       fZtot += fZZ[i];
345       fNtot += fNN[i];
346    }                 
347    
348
349 }
350
351 //_____________________________________________________________________________
352 Float_t AliZDCFragment::DeuteronFraction()
353 {
354     Float_t DeuteronProdPar[2] = {1.068,0.0385};
355     Float_t DeutFrac = DeuteronProdPar[0]-DeuteronProdPar[1]*fB;
356     if(DeutFrac>1.) DeutFrac=1.;
357     return DeutFrac;
358 }