From: coppedis Date: Fri, 23 Feb 2001 16:46:42 +0000 (+0000) Subject: Class for spectator nucleons fragmentation X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=commitdiff_plain;h=28e0901a81c086a85a31ae310a93b4ca23d1767b;ds=inline Class for spectator nucleons fragmentation --- diff --git a/ZDC/AliZDCFragment.cxx b/ZDC/AliZDCFragment.cxx new file mode 100644 index 00000000000..c56ad5511dd --- /dev/null +++ b/ZDC/AliZDCFragment.cxx @@ -0,0 +1,344 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +// --- Standard libraries +#include +#include + +// --- ROOT system +#include +#include + +// --- AliRoot classes +#include "AliZDCFragment.h" + +ClassImp(AliZDCFragment) + +//_____________________________________________________________________________ +AliZDCFragment::AliZDCFragment() +{ + // + // Default constructor + // + fB = 0; +} + +//_____________________________________________________________________________ +AliZDCFragment::AliZDCFragment(Float_t b) + : TNamed(" "," ") +{ + // + // Standard constructor + // + fB = b; + fZbAverage = 0; + fNimf = 0; + fZmax = 0; + fTau = 0; + fNalpha = 0; + fZtot = 0; + fNtot = 0; + for(Int_t i=0; i<=99; i++){ + fZZ[i] = 0; + fNN[i] = 0; + } + +} + +//_____________________________________________________________________________ +void AliZDCFragment::GenerateIMF(Int_t* fZZ, Int_t &fNalpha) +{ + // Coefficients of polynomial for average number of IMF + Float_t ParamNimf[5]={0.011236,1.8364,56.572,-116.24,58.289}; + // Coefficients of polynomial for fluctuations on average number of IMF + Float_t ParamFluctNimf[4]={-0.13176,2.9392,-5.2147,2.3092}; + // Coefficients of polynomial for average maximum Z of fragments + Float_t ParamZmax[4]={0.16899,14.203,-2.8284,65.036}; + // Coefficients of polynomial for fluctuations on maximum Z of fragments + Float_t ParamFluctZmax[5]={0.013782,-0.17282,1.5065,1.0654,-2.4317}; + // Coefficients of polynomial for exponent tau of fragments Z distribution + Float_t ParamTau[3]={6.7233,-15.85,13.047}; + //Coefficients of polynomial for average number of alphas + Float_t ParamNalpha[4]={-0.68554,39.605,-68.311,30.165}; + // Coefficients of polynomial for fluctuations on average number of alphas + Float_t ParamFluctNalpha[5]={0.283,6.2141,-17.113,17.394,-6.6084}; + // Coefficients of function for Pb nucleus skin + Float_t ParamSkinPb[2]={0.93,11.05}; + + // Thickness of nuclear surface + Float_t NuclearThick = 0.52; + // Maximum impact parameter for U [r0*A**(1/3)] + Float_t bMaxU = 14.87; + // Maximum impact parameter for Pb [r0*A**(1/3)] + Float_t bMaxPb = 14.22; + // Z of the projectile + Float_t ZProj = 82.; + + // From b(Pb) to b(U) + Float_t bU = fB*bMaxU/bMaxPb; + + // From b(U) to Zbound(U) + // --- A.Schuttauf et al, Nuc.Phys. A607 (1996) 457 --------------- + // From geometrical consideration and from dsigma/dZbound for U+U, + // which is approx. constant, the constant value is found + // integrating the nucleus cross surface from 0 to bmax=R1+R2 where + // R = 1.2*A**(1/3). This value has been measured in Aladin (U+U). + Float_t ZbU = bU*bU*TMath::Pi()/7.48; + + // Rescale Zbound for Pb + fZbAverage = ZProj/92.*ZbU; + + // Zbound is proportional to b**2 up to b < bMaxPb-2*NuclearThick + // and then it is an increasing exponential, imposing that at + // b=bMaxPb-2NuclearThick the two functions have the same derivative + Float_t bCore = bMaxPb-2*NuclearThick; + if(fB>bCore){ + fZbAverage=ZProj*(1.-TMath::Exp(-ParamSkinPb[0]*(fB-ParamSkinPb[1]))); + } + if(fZbAverage>ZProj) fZbAverage = ZProj; + Float_t ZbNorm = fZbAverage/ZProj; + Float_t bNorm = fB/bMaxPb; + + // From Zbound to ,,tau + // Polinomial fits to Aladin distribution + // --- A.Schuttauf et al, Nuc.Phys. A607 (1996) 457. + Float_t AverageNimf = ParamNimf[0]+ParamNimf[1]*ZbNorm+ParamNimf[2]* + TMath::Power(ZbNorm,2)+ParamNimf[3]*TMath::Power(ZbNorm,3)+ + ParamNimf[4]*TMath::Power(ZbNorm,4); + + // Add fluctuation: from Singh et al. + Float_t FluctNimf = ParamFluctNimf[0]+ParamFluctNimf[1]*ZbNorm+ + ParamFluctNimf[2]*TMath::Power(ZbNorm,2)+ParamFluctNimf[3] + *TMath::Power(ZbNorm,3); + Float_t xx = gRandom->Gaus(0.0,1.0); + FluctNimf = FluctNimf*xx; + fNimf = Int_t(AverageNimf+FluctNimf); + Float_t y = gRandom->Rndm(); + if(y < ((AverageNimf+FluctNimf)-fNimf)) fNimf += 1; + if(fNimf ==0 && ZbNorm>0.75) fNimf = 1; +// printf("\n fNimf = %d\n", fNimf); + + Float_t AverageZmax = ParamZmax[0]+ParamZmax[1]*ZbNorm+ParamZmax[2]* + TMath::Power(ZbNorm,2)+ParamZmax[3]*TMath::Power(ZbNorm,3); + fTau = ParamTau[0]+ParamTau[1]*ZbNorm+ParamTau[2]*TMath::Power(ZbNorm,2); + + // Add fluctuation to mean value of Zmax (see Hubele) + Float_t FluctZmax = ParamFluctZmax[0]+ParamFluctZmax[1]*ZbNorm+ + ParamFluctZmax[2]*TMath::Power(ZbNorm,2)+ParamFluctZmax[3]* + TMath::Power(ZbNorm,3)+ParamFluctZmax[4]*TMath::Power(ZbNorm,4); + FluctZmax = FluctZmax*ZProj/6.; + Float_t xg = gRandom->Gaus(0.0,1.0); + FluctZmax = FluctZmax*xg; + fZmax = AverageZmax+FluctZmax; + if(fZmax>ZProj) fZmax = ZProj; +// printf("\n fZmax = %f\n", fZmax); + + // Find the number of alpha particles + // from Singh et al. : Pb+emulsion + Float_t AverageAlpha = ParamNalpha[0]+ParamNalpha[1]*ZbNorm+ + ParamNalpha[2]*TMath::Power(ZbNorm,2)+ParamNalpha[3]* + TMath::Power(ZbNorm,3); + Float_t FluctAlpha = ParamFluctNalpha[0]+ParamFluctNalpha[1]* + ZbNorm+ParamFluctNalpha[2]*TMath::Power(ZbNorm,2)+ + ParamFluctNalpha[3]*TMath::Power(ZbNorm,3)+ + ParamFluctNalpha[4]*TMath::Power(ZbNorm,4); + Float_t xxx = gRandom->Gaus(0.0,1.0); + FluctAlpha = FluctAlpha*xxx; + fNalpha = Int_t(AverageAlpha+FluctAlpha); + Float_t yy = gRandom->Rndm(); + if(yy < ((AverageAlpha+FluctAlpha)-fNalpha)) fNalpha += 1; + + Int_t Choice = 0; + Float_t ZbFrag = 0, SumZ = 0.; + // 2 possibilities: + // 1) for bNorm < 0.9 ==> first remove alphas, then fragments + // 2) for bNorm > 0.9 ==> first remove fragments, then alphas + + if(bNorm<=0.9) { + // remove alpha from zbound to find zbound for fragments (Z>=3) + ZbFrag = fZbAverage-fNalpha*2; + Choice = 1; + } + else { + ZbFrag = fZbAverage; + Choice = 0; + } +// printf("\n Choice = %d, fZbAverage = %f, ZbFrag = %f \n", Choice, fZbAverage, ZbFrag); + + + // Check if ZbFrag < fZmax + if(ZbFrag<=fZmax) { + if(fNimf>0 && ZbFrag>=2){ + fNimf = 1; + fZZ[0] = Int_t(ZbFrag); + SumZ = ZbFrag; + } + else { + fNimf = 0; + } + return; + } + + // Prepare the exponential charge distribution dN/dZ + if(fZmax <= 0.01) { + fNimf = 0; + return; + } + if(fNimf == 0) { + fNimf = 0; + return; + } + + TF1 *funTau = new TF1("funTau","1./(x**[0])",0.01,fZmax); + funTau->SetParameter(0,fTau); + + // Extract randomly the charge of the fragments from the distribution + + Float_t zz[fNimf]; + for(Int_t j=0; jGetRandom()); +// printf("\n zz[%d] = %f \n",i,zz[i]); + } + delete funTau; + + // Sorting vector in ascending order + +// int comp(const void*, const void*); + qsort((void*)zz,fNimf,sizeof(float),comp); + + for(Int_t i=0; i zz[%d] = %f \n",i,zz[i]); + } + + // Rescale the maximum charge to fZmax + for(Int_t j=0; jZbFrag){ + for(Int_t i=0; i< fNimf; i++){ + k += 1; + SumZ -= fZZ[i]; + if(SumZ<=ZbFrag){ + fNimf -= (i+1); + break; + } + } + } + else { + if(Choice == 1) return; + Int_t iDiff = Int_t((ZbFrag-SumZ)/2); + if(iDiff fNimf = %d, SumZ = %f, fZmax = %f\n", +// fNimf, SumZ, fZmax); +// printf("\n fNalpha = %d\n", fNalpha); +// for(Int_t j=0; j + +int comp(const void*, const void*); + +class AliZDCFragment : public TNamed { + +public: + AliZDCFragment(); + AliZDCFragment(Float_t b); + virtual ~AliZDCFragment() {} + void GenerateIMF(Int_t* fZZ, Int_t &fNalpha); + void AttachNeutrons(Int_t* fZZ, Int_t* fNN, Int_t &Ztot, Int_t &Ntot); + + // Setting parameters + virtual void SetImpactParameter(Float_t b) {fB=b;}; + + // Getting parameters + Int_t GetFragmentNum() {return fNimf;}; + + +protected: + + Float_t fB; // Impact parameter + Float_t fZbAverage ; // Mean value of Z bound + Int_t fNimf; // Number of IMF + Float_t fZmax; // Mean value of maximum Z of fragment + Float_t fTau; // Exponent of charge distribution: dN/dZ = Z*exp(-fTau) + Int_t fZZ[100]; // Array of atomic numbers of fragments + Int_t fNalpha; // Number of alpha particles + + Int_t fNN[100]; // Array of number of neutrons of fragments + Int_t fZtot; // Total number of bound protons + Int_t fNtot; // Total number of bound neutrons + + + ClassDef(AliZDCFragment,1) // Generator for AliZDC fragment class +}; + +#endif