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