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 }
  53
  54 //_____________________________________________________________________________
  55 AliZDCFragment::AliZDCFragment(Float_t b):
  56   TNamed(" "," "),
  57   fB(b),
  58   fZbAverage(0),
  59   fNimf(0),
  60   fZmax(0),
  61   fTau(0),
  62   fNalpha(0),
  63   fZtot(0),
  64   fNtot(0)
  65 {
  66   //
  67   // Standard constructor
  68   //
  69   for(Int_t i=0; i<=99; i++){
  70      fZZ[i] = 0;
  71      fNN[i] = 0;
  72   }
  73
  74 }
  75
  76 //_____________________________________________________________________________
  77 void AliZDCFragment::GenerateIMF()
  78 {
  79
  80    // Loop variables
  81   Int_t i,j;
  82
  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.};
 100
 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.;
 110
 111    // From b(Pb) to b(U)
 112    if(fB>kbMaxPb) fB = 2*kbMaxPb-fB;
 113
 114    Float_t  bU = fB*kbMaxU/kbMaxPb;
 115
 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;
 123
 124    //  Rescale Zbound for Pb
 125    fZbAverage = kZProj/92.*zbU;
 126
 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;
 131    if(fB>kbMaxPb){
 132      fZbAverage = TMath::Exp(-kParamSkinPb[0]*(fB-kParamSkinPb[1]));
 133      //printf(" b = %1.2f fm   Z_bound %1.2f\n", fB, fZbAverage);
 134    }
 135    if(fZbAverage>kZProj) fZbAverage = kZProj;
 136    Float_t zbNorm = fZbAverage/kZProj;
 137    Float_t bNorm = fB/kbMaxPb;
 138
 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);
 145
 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;
 156
 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);
 160
 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;
 170
 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);
 175
 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;
 190
 191    // 2 possibilities:
 192    // 1) for bNorm < 0.9 ==> first remove alphas, then fragments
 193    // 2) for bNorm > 0.9 ==> first remove fragments, then alphas
 194
 195    Int_t choice = 0;
 196    Float_t zbFrag = 0, sumZ = 0.;
 197
 198    if(bNorm<=0.9) {
 199    // remove alpha from zbound to find zbound for fragments  (Z>=3)
 200      zbFrag = fZbAverage-fNalpha*2;
 201      choice = 1;
 202    }
 203    else {
 204      zbFrag = fZbAverage;
 205      choice = 0;
 206    }
 207 //   printf("\n choice = %d, fZbAverage = %f, zbFrag = %f \n", choice, fZbAverage, zbFrag);
 208
 209
 210    // Check if zbFrag < fZmax
 211    if(zbFrag<=fZmax) {
 212      if(fNimf>0 && zbFrag>=2){
 213        fNimf = 1;
 214        fZZ[0] = Int_t(zbFrag);
 215        sumZ = zbFrag;
 216      }
 217      else {
 218        fNimf = 0;
 219      }
 220      return;
 221    }
 222
 223    // Prepare the exponential charge distribution dN/dZ
 224    if(fZmax <= 0.01) {
 225      fNimf = 0;
 226      return;
 227    }
 228    if(fNimf == 0) {
 229      fNimf = 0;
 230      return;
 231    }
 232
 233    TF1 *funTau = new TF1("funTau","1./(x**[0])",0.01,fZmax);
 234    funTau->SetParameter(0,fTau);
 235
 236    // Extract randomly the charge of the fragments from the distribution
 237
 238    Float_t * zz = new Float_t[fNimf];
 239    for(j=0; j<fNimf; j++){
 240       zz[j] =0;
 241    }
 242    for(i=0; i<fNimf; i++){
 243       zz[i] = Float_t(funTau->GetRandom());
 244 //      printf("\n      zz[%d] = %f \n",i,zz[i]);
 245    }
 246    delete funTau;
 247
 248    // Sorting vector in ascending order with C function QSORT
 249    qsort((void*)zz,fNimf,sizeof(Float_t),comp);
 250
 251
 252 //   for(Int_t i=0; i<fNimf; i++){
 253 //      printf("\n After sorting -> zz[%d] = %f \n",i,zz[i]);
 254 //   }
 255
 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]);
 261    }
 262
 263    delete[] zz;
 264
 265    // Check that the sum of the bound charges is not > than Zbound-Zalfa
 266
 267    for(Int_t ii=0; ii<fNimf; ii++){
 268      sumZ += fZZ[ii];
 269    }
 270
 271    Int_t k = 0;
 272    if(sumZ>zbFrag){
 273      for(i=0; i< fNimf; i++){
 274        k += 1;
 275        sumZ -= fZZ[i];
 276        if(sumZ<=zbFrag){
 277          fNimf -= (i+1);
 278          break;
 279        }
 280      }
 281    }
 282    else {
 283      if(choice == 1) return;
 284      Int_t iDiff = Int_t((zbFrag-sumZ)/2);
 285      if(iDiff<fNalpha){
 286        fNalpha=iDiff;
 287        return;
 288      }
 289      else{
 290        return;
 291      }
 292    }
 293
 294    fNimf += k;
 295    for(i=0; i<fNimf; i++){
 296      fZZ[i] = fZZ[i+k];
 297    }
 298    fNimf -= k;
 299
 300    sumZ=0;
 301    for(i=0; i<fNimf; i++){
 302      sumZ += fZZ[i];
 303    }
 304
 305 }
 306
 307 //_____________________________________________________________________________
 308 void AliZDCFragment::AttachNeutrons()
 309 {
 310 //
 311 // Prepare nuclear fragment by attaching a suitable number of neutrons
 312 //
 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,
 322                      128.45793,
 323                      130.32111,141.51236,
 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,
 328                      221.74295};
 329    const Int_t kZIon[68]={1,1,2,3,3,
 330                      4,4,5,5,6,
 331                      7,8,9,10,11,
 332                      12,13,14,15,16,
 333                      17,18,19,20,21,
 334                      22,23,24,25,26,
 335                      27,28,29,30,32,
 336                      34,36,38,40,42,
 337                      46,48,50,54,56,
 338                      58,62,
 339                      63,64,65,66,67,
 340                      68,69,70,71,72,
 341                      73,74,75,76,77,
 342                      78,79,80,81,82,
 343                      92};
 344
 345    Int_t iZ, iA;
 346 //   printf("\n fNimf=%d\n",fNimf);
 347
 348    for(Int_t i=0; i<fNimf; i++) {
 349       for(Int_t j=0; j<68; j++) {
 350         iZ = kZIon[j];
 351         if((fZZ[i]-iZ) == 0){
 352           iA = Int_t(kAIon[j]/0.93149432+0.5);
 353           fNN[i] = iA - iZ;
 354           break;
 355         }
 356         else if((fZZ[i]-iZ) < 0){
 357           fZZ[i] = kZIon[j-1];
 358           iA = Int_t (kAIon[j-1]/0.93149432+0.5);
 359           fNN[i] = iA - kZIon[j-1];
 360           break;
 361         }
 362       }
 363       fZtot += fZZ[i];
 364       fNtot += fNN[i];
 365    }
 366
 367
 368 }
 369
 370 //_____________________________________________________________________________
 371 Float_t AliZDCFragment::DeuteronNumber()
 372 {
 373     // Calculates the fraction of deuterum nucleus produced
 374     //
 375     Float_t deuteronProdPar[2] = {-0.068,0.0385};
 376     Float_t deutNum = deuteronProdPar[0] + deuteronProdPar[1]*fB;
 377     if(deutNum<0.) deutNum = 0.;
 378     return deutNum;
 379 }