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