[u/mrichter/AliRoot.git] / ZDC / AliZDCFragment.cxx

/**************************************************************************
 * 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.                  *
 **************************************************************************/


// ******************************************************************
//
//  	Class for nuclear fragments formation
//
// ******************************************************************

// --- Standard libraries
#include <stdlib.h>

// --- ROOT system
#include <TRandom.h>
#include <TF1.h>

// --- AliRoot classes
#include "AliZDCFragment.h"
 
ClassImp(AliZDCFragment)
   
int comp(const void *i,const void *j) {return *(int *)i - *(int *)j;}


//_____________________________________________________________________________
AliZDCFragment::AliZDCFragment():
  fB(0),
  fZbAverage(0),
  fNimf(0),
  fZmax(0),
  fTau(0),
  fNalpha(0),
  fZtot(0),
  fNtot(0)
{
  //
  // Default constructor
  //
}

//_____________________________________________________________________________
AliZDCFragment::AliZDCFragment(Float_t b): 
  TNamed(" "," "),
  fB(b),
  fZbAverage(0),
  fNimf(0),
  fZmax(0),
  fTau(0),
  fNalpha(0),
  fZtot(0),
  fNtot(0)
{
  //
  // Standard constructor
  //
  for(Int_t i=0; i<=99; i++){
     fZZ[i] = 0;
     fNN[i] = 0;
  }
  
}

//_____________________________________________________________________________
void AliZDCFragment::GenerateIMF()
{

   // Loop variables
  Int_t i,j;

   // Coefficients of polynomial for average number of IMF
   const Float_t  kParamNimf[5]={0.011236,1.8364,56.572,-116.24,58.289}; 
   // Coefficients of polynomial for fluctuations on average number of IMF
   const Float_t  kParamFluctNimf[4]={-0.13176,2.9392,-5.2147,2.3092}; 
   // Coefficients of polynomial for average maximum Z of fragments
   //const Float_t  kParamZmax[4]={0.16899,14.203,-2.8284,65.036}; 
   const Float_t  kParamZmax[4]={0.16899,14.203,-2.8284,70.5}; 
   // Coefficients of polynomial for fluctuations on maximum Z of fragments
   const Float_t  kParamFluctZmax[5]={0.013782,-0.17282,1.5065,1.0654,-2.4317}; 
   // Coefficients of polynomial for exponent tau of fragments Z distribution
   const Float_t  kParamTau[3]={6.7233,-15.85,13.047};  
   //Coefficients of polynomial for average number of alphas
   const Float_t  kParamNalpha[4]={-0.68554,39.605,-68.311,30.165}; 
   // Coefficients of polynomial for fluctuations on average number of alphas
   const Float_t  kParamFluctNalpha[5]={0.283,6.2141,-17.113,17.394,-6.6084}; 
   // Coefficients of function for Pb nucleus skin
   const Float_t  kParamSkinPb[2]={0.762408, 20.};
   
   // Thickness of nuclear surface
   const Float_t  kNuclearThick = 0.52;
   // Maximum impact parameter for U [r0*A**(1/3)]
   const Float_t  kbMaxU = 14.87;
   // Maximum impact parameter for Pb [r0*A**(1/3)]
   //const Float_t  kbMaxPb = 14.22+4*kNuclearThick;
   const Float_t  kbMaxPb = 14.22;
   // Z of the projectile
   const Float_t  kZProj = 82.;
   
   // From b(Pb) to b(U)
   if(fB>kbMaxPb) fB = 2*kbMaxPb-fB;
   
   Float_t  bU = fB*kbMaxU/kbMaxPb;
    
   // 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 = kZProj/92.*zbU;
   
   // Zbound is proportional to b**2 up to b < kbMaxPb-2*kNuclearThick
   // and then it is an increasing exponential, imposing that at 
   // b=kbMaxPb-2kNuclearThick the two functions have the same derivative
   //Float_t bCore = kbMaxPb-2*kNuclearThick;
   if(fB>kbMaxPb){
     fZbAverage = TMath::Exp(-kParamSkinPb[0]*(fB-kParamSkinPb[1]));
     //printf(" b = %1.2f fm   Z_bound %1.2f\n", fB, fZbAverage);
   }
   if(fZbAverage>kZProj) fZbAverage = kZProj;
   Float_t zbNorm = fZbAverage/kZProj;
   Float_t bNorm = fB/kbMaxPb;
   
   // From Zbound to <Nimf>,<Zmax>,tau
   // Polinomial fits to Aladin distribution
   // --- A.Schuttauf et al, Nuc.Phys. A607 (1996) 457.
   Float_t averageNimf = kParamNimf[0]+kParamNimf[1]*zbNorm+kParamNimf[2]*
           TMath::Power(zbNorm,2)+kParamNimf[3]*TMath::Power(zbNorm,3)+
	   kParamNimf[4]*TMath::Power(zbNorm,4);
   
   // Add fluctuation: from Singh et al. 
   Float_t fluctNimf = kParamFluctNimf[0]+kParamFluctNimf[1]*zbNorm+
           kParamFluctNimf[2]*TMath::Power(zbNorm,2)+kParamFluctNimf[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;
   
   Float_t averageZmax = kParamZmax[0]+kParamZmax[1]*zbNorm+kParamZmax[2]*
           TMath::Power(zbNorm,2)+kParamZmax[3]*TMath::Power(zbNorm,3);
   fTau = kParamTau[0]+kParamTau[1]*zbNorm+kParamTau[2]*TMath::Power(zbNorm,2);
   
   // Add fluctuation to mean value of Zmax (see Hubele)
   Float_t fluctZmax = kParamFluctZmax[0]+kParamFluctZmax[1]*zbNorm+
           kParamFluctZmax[2]*TMath::Power(zbNorm,2)+kParamFluctZmax[3]*
	   TMath::Power(zbNorm,3)+kParamFluctZmax[4]*TMath::Power(zbNorm,4);
   fluctZmax = fluctZmax*kZProj/6.;
   Float_t xg = gRandom->Gaus(0.0,1.0);
   fluctZmax = fluctZmax*xg;
   fZmax = (averageZmax+fluctZmax);
   if(fZmax>kZProj) fZmax = kZProj;
   
//   printf("\n\n ------------------------------------------------------------");   
//   printf("\n Generation of nuclear fragments\n");   
//   printf("\n fNimf = %d\n", fNimf);   
//   printf("\n fZmax = %f\n", fZmax); 

   // Find the number of alpha particles 
   // from Singh et al. : Pb+emulsion
   Float_t averageAlpha = kParamNalpha[0]+kParamNalpha[1]*zbNorm+
           kParamNalpha[2]*TMath::Power(zbNorm,2)+kParamNalpha[3]*
	   TMath::Power(zbNorm,3);
   Float_t fluctAlpha = kParamFluctNalpha[0]+kParamFluctNalpha[1]*
           zbNorm+kParamFluctNalpha[2]*TMath::Power(zbNorm,2)+
	   kParamFluctNalpha[3]*TMath::Power(zbNorm,3)+
	   kParamFluctNalpha[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;

   // 2 possibilities:
   // 1) for bNorm < 0.9 ==> first remove alphas, then fragments
   // 2) for bNorm > 0.9 ==> first remove fragments, then alphas

   Int_t choice = 0;
   Float_t zbFrag = 0, sumZ = 0.;

   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 = new Float_t[fNimf];
   for(j=0; j<fNimf; j++){
      zz[j] =0;
   }
   for(i=0; i<fNimf; i++){
      zz[i] = Float_t(funTau->GetRandom());
//      printf("\n	zz[%d] = %f \n",i,zz[i]);
   }
   delete funTau;
   
   // Sorting vector in ascending order with C function QSORT 
   qsort((void*)zz,fNimf,sizeof(Float_t),comp);

   
//   for(Int_t i=0; i<fNimf; i++){
//      printf("\n After sorting -> zz[%d] = %f \n",i,zz[i]);
//   }
   
   // Rescale the maximum charge to fZmax
   for(j=0; j<fNimf; j++){
     fZZ[j] = Int_t (zz[j]*fZmax/zz[fNimf-1]);
     if(fZZ[j]<3) fZZ[j] = 3;
//     printf("\n 	fZZ[%d] = %d \n",j,fZZ[j]);
   }

   delete[] zz;
   
   // Check that the sum of the bound charges is not > than Zbound-Zalfa
   
   for(Int_t ii=0; ii<fNimf; ii++){
     sumZ += fZZ[ii];
   }
   
   Int_t k = 0;
   if(sumZ>zbFrag){
     for(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<fNalpha){
       fNalpha=iDiff;
       return;
     }
     else{
       return;
     }
   }

   fNimf += k;
   for(i=0; i<fNimf; i++){
     fZZ[i] = fZZ[i+k];
   }
   fNimf -= k;
   
   sumZ=0;
   for(i=0; i<fNimf; i++){
     sumZ += fZZ[i];
   }
   
}

//_____________________________________________________________________________
void AliZDCFragment::AttachNeutrons()
{
//
// Prepare nuclear fragment by attaching a suitable number of neutrons
//
   const Float_t kAIon[68]={1.87612,2.80943,3.7284,5.60305,6.53536,
     		     6.53622,8.39479,9.32699,10.2551,11.17793,
     		     13.04378,14.89917,17.6969,18.62284,21.41483,
     		     22.34193,25.13314,26.06034,28.85188,29.7818,
     		     32.57328,33.50356,36.29447,37.22492,41.87617,
     		     44.66324,47.45401,48.38228,51.17447,52.10307,
     		     54.89593,53.96644,58.61856,59.54963,68.85715,
     		     74.44178,78.16309,81.88358,83.74571,91.19832,
     		     98.64997,106.10997,111.68821,122.86796,
     		     128.45793,
     		     130.32111,141.51236,
     		     141.55,146.477,148.033,152.699,153.631,
     		     155.802,157.357,162.022,162.984,166.2624,
     		     168.554,171.349,173.4536,177.198,179.0518,
     		     180.675,183.473,188.1345,190.77,193.729,
     		     221.74295};
   const Int_t kZIon[68]={1,1,2,3,3,
     		     4,4,5,5,6,
     		     7,8,9,10,11,
     		     12,13,14,15,16,
     		     17,18,19,20,21,
     		     22,23,24,25,26,
     		     27,28,29,30,32,
     		     34,36,38,40,42,
     		     46,48,50,54,56,
     		     58,62,
     		     63,64,65,66,67,
     		     68,69,70,71,72,
     		     73,74,75,76,77,
     		     78,79,80,81,82,
      		     92};
    
   Int_t iZ, iA;  
//   printf("\n fNimf=%d\n",fNimf);  

   for(Int_t i=0; i<fNimf; i++) {
      for(Int_t j=0; j<68; j++) {
        iZ = kZIon[j];
	if((fZZ[i]-iZ) == 0){
	  iA = Int_t(kAIon[j]/0.93149432+0.5);
	  fNN[i] = iA - iZ;
          break;
	}
	else if((fZZ[i]-iZ) < 0){
	  fZZ[i] = kZIon[j-1];
	  iA = Int_t (kAIon[j-1]/0.93149432+0.5);
	  fNN[i] = iA - kZIon[j-1];
          break;
	}
      }
      fZtot += fZZ[i];
      fNtot += fNN[i];
   }		     
   

}

//_____________________________________________________________________________
Float_t AliZDCFragment::DeuteronNumber()
{
    // Calculates the fraction of deuterum nucleus produced
    //
    Float_t deuteronProdPar[2] = {-0.068,0.0385};
    Float_t deutNum = deuteronProdPar[0] + deuteronProdPar[1]*fB;
    if(deutNum<0.) deutNum = 0.;
    return deutNum;
}
Commit	Line	Data
28e0901a	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
8a2624cc	16
	17	// ******************************************************************
	18	//
	19	// Class for nuclear fragments formation
	20	//
	21	// ******************************************************************
	22
28e0901a	23	// --- Standard libraries
28e0901a	24	#include <stdlib.h>
28e0901a	25
	26	// --- ROOT system
	27	#include <TRandom.h>
	28	#include <TF1.h>
	29
	30	// --- AliRoot classes
	31	#include "AliZDCFragment.h"
28e0901a	32
5a881c97	33	ClassImp(AliZDCFragment)
	34
	35	int comp(const void i,const void j) {return (int )i - (int )j;}
	36
	37
28e0901a	38	//_____________________________________________________________________________
a718c993	39	AliZDCFragment::AliZDCFragment():
ab28a71c	40	fB(0),
	41	fZbAverage(0),
	42	fNimf(0),
	43	fZmax(0),
	44	fTau(0),
	45	fNalpha(0),
	46	fZtot(0),
	47	fNtot(0)
28e0901a	48	{
	49	//
	50	// Default constructor
	51	//
28e0901a	52	}
	53
	54	//_____________________________________________________________________________
ab28a71c	55	AliZDCFragment::AliZDCFragment(Float_t b):
ab28a71c	56	TNamed(" "," "),
a718c993	57	fB(b),
	58	fZbAverage(0),
	59	fNimf(0),
	60	fZmax(0),
	61	fTau(0),
	62	fNalpha(0),
	63	fZtot(0),
	64	fNtot(0)
28e0901a	65	{
	66	//
	67	// Standard constructor
	68	//
28e0901a	69	for(Int_t i=0; i<=99; i++){
	70	fZZ[i] = 0;
	71	fNN[i] = 0;
	72	}
	73
	74	}
	75
	76	//_____________________________________________________________________________
92339a90	77	void AliZDCFragment::GenerateIMF()
28e0901a	78	{
583731c7	79
	80	// Loop variables
	81	Int_t i,j;
	82
28e0901a	83	// Coefficients of polynomial for average number of IMF
86d81a3c	84	const Float_t kParamNimf[5]={0.011236,1.8364,56.572,-116.24,58.289};
28e0901a	85	// Coefficients of polynomial for fluctuations on average number of IMF
86d81a3c	86	const Float_t kParamFluctNimf[4]={-0.13176,2.9392,-5.2147,2.3092};
28e0901a	87	// Coefficients of polynomial for average maximum Z of fragments
1f359fbe	88	//const Float_t kParamZmax[4]={0.16899,14.203,-2.8284,65.036};
1f359fbe	89	const Float_t kParamZmax[4]={0.16899,14.203,-2.8284,70.5};
28e0901a	90	// Coefficients of polynomial for fluctuations on maximum Z of fragments
86d81a3c	91	const Float_t kParamFluctZmax[5]={0.013782,-0.17282,1.5065,1.0654,-2.4317};
28e0901a	92	// Coefficients of polynomial for exponent tau of fragments Z distribution
86d81a3c	93	const Float_t kParamTau[3]={6.7233,-15.85,13.047};
28e0901a	94	//Coefficients of polynomial for average number of alphas
86d81a3c	95	const Float_t kParamNalpha[4]={-0.68554,39.605,-68.311,30.165};
28e0901a	96	// Coefficients of polynomial for fluctuations on average number of alphas
86d81a3c	97	const Float_t kParamFluctNalpha[5]={0.283,6.2141,-17.113,17.394,-6.6084};
28e0901a	98	// Coefficients of function for Pb nucleus skin
c15e4ed3	99	const Float_t kParamSkinPb[2]={0.762408, 20.};
28e0901a	100
28e0901a	101	// Thickness of nuclear surface
86d81a3c	102	const Float_t kNuclearThick = 0.52;
28e0901a	103	// Maximum impact parameter for U [r0A*(1/3)]
86d81a3c	104	const Float_t kbMaxU = 14.87;
28e0901a	105	// Maximum impact parameter for Pb [r0A*(1/3)]
c15e4ed3	106	//const Float_t kbMaxPb = 14.22+4*kNuclearThick;
c15e4ed3	107	const Float_t kbMaxPb = 14.22;
28e0901a	108	// Z of the projectile
86d81a3c	109	const Float_t kZProj = 82.;
28e0901a	110
28e0901a	111	// From b(Pb) to b(U)
1f359fbe	112	if(fB>kbMaxPb) fB = 2*kbMaxPb-fB;
1f359fbe	113
86d81a3c	114	Float_t bU = fB*kbMaxU/kbMaxPb;
28e0901a	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.2A*(1/3). This value has been measured in Aladin (U+U).
86d81a3c	122	Float_t zbU = bUbUTMath::Pi()/7.48;
28e0901a	123
28e0901a	124	// Rescale Zbound for Pb
86d81a3c	125	fZbAverage = kZProj/92.*zbU;
28e0901a	126
86d81a3c	127	// Zbound is proportional to b*2 up to b < kbMaxPb-2kNuclearThick
28e0901a	128	// and then it is an increasing exponential, imposing that at
86d81a3c	129	// b=kbMaxPb-2kNuclearThick the two functions have the same derivative
c15e4ed3	130	//Float_t bCore = kbMaxPb-2*kNuclearThick;
	131	if(fB>kbMaxPb){
	132	fZbAverage = TMath::Exp(-kParamSkinPb[0]*(fB-kParamSkinPb[1]));
30158a90	133	//printf(" b = %1.2f fm Z_bound %1.2f\n", fB, fZbAverage);
c15e4ed3	134	}
86d81a3c	135	if(fZbAverage>kZProj) fZbAverage = kZProj;
	136	Float_t zbNorm = fZbAverage/kZProj;
	137	Float_t bNorm = fB/kbMaxPb;
28e0901a	138
	139	// From Zbound to <Nimf>,<Zmax>,tau
	140	// Polinomial fits to Aladin distribution
	141	// --- A.Schuttauf et al, Nuc.Phys. A607 (1996) 457.
86d81a3c	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);
28e0901a	145
28e0901a	146	// Add fluctuation: from Singh et al.
86d81a3c	147	Float_t fluctNimf = kParamFluctNimf[0]+kParamFluctNimf[1]*zbNorm+
	148	kParamFluctNimf[2]*TMath::Power(zbNorm,2)+kParamFluctNimf[3]
	149	*TMath::Power(zbNorm,3);
28e0901a	150	Float_t xx = gRandom->Gaus(0.0,1.0);
86d81a3c	151	fluctNimf = fluctNimf*xx;
86d81a3c	152	fNimf = Int_t(averageNimf+fluctNimf);
28e0901a	153	Float_t y = gRandom->Rndm();
86d81a3c	154	if(y < ((averageNimf+fluctNimf)-fNimf)) fNimf += 1;
86d81a3c	155	if(fNimf ==0 && zbNorm>0.75) fNimf = 1;
28e0901a	156
86d81a3c	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);
28e0901a	160
28e0901a	161	// Add fluctuation to mean value of Zmax (see Hubele)
86d81a3c	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.;
28e0901a	166	Float_t xg = gRandom->Gaus(0.0,1.0);
86d81a3c	167	fluctZmax = fluctZmax*xg;
1f359fbe	168	fZmax = (averageZmax+fluctZmax);
86d81a3c	169	if(fZmax>kZProj) fZmax = kZProj;
5a881c97	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);
28e0901a	175
	176	// Find the number of alpha particles
	177	// from Singh et al. : Pb+emulsion
86d81a3c	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);
28e0901a	185	Float_t xxx = gRandom->Gaus(0.0,1.0);
86d81a3c	186	fluctAlpha = fluctAlpha*xxx;
86d81a3c	187	fNalpha = Int_t(averageAlpha+fluctAlpha);
28e0901a	188	Float_t yy = gRandom->Rndm();
86d81a3c	189	if(yy < ((averageAlpha+fluctAlpha)-fNalpha)) fNalpha += 1;
28e0901a	190
28e0901a	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
86d81a3c	195	Int_t choice = 0;
86d81a3c	196	Float_t zbFrag = 0, sumZ = 0.;
5a881c97	197
28e0901a	198	if(bNorm<=0.9) {
28e0901a	199	// remove alpha from zbound to find zbound for fragments (Z>=3)
86d81a3c	200	zbFrag = fZbAverage-fNalpha*2;
86d81a3c	201	choice = 1;
28e0901a	202	}
28e0901a	203	else {
86d81a3c	204	zbFrag = fZbAverage;
86d81a3c	205	choice = 0;
28e0901a	206	}
86d81a3c	207	// printf("\n choice = %d, fZbAverage = %f, zbFrag = %f \n", choice, fZbAverage, zbFrag);
28e0901a	208
28e0901a	209
86d81a3c	210	// Check if zbFrag < fZmax
	211	if(zbFrag<=fZmax) {
	212	if(fNimf>0 && zbFrag>=2){
28e0901a	213	fNimf = 1;
86d81a3c	214	fZZ[0] = Int_t(zbFrag);
86d81a3c	215	sumZ = zbFrag;
28e0901a	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
4b016189	238	Float_t * zz = new Float_t[fNimf];
583731c7	239	for(j=0; j<fNimf; j++){
28e0901a	240	zz[j] =0;
28e0901a	241	}
583731c7	242	for(i=0; i<fNimf; i++){
28e0901a	243	zz[i] = Float_t(funTau->GetRandom());
	244	// printf("\n zz[%d] = %f \n",i,zz[i]);
	245	}
	246	delete funTau;
	247
5a881c97	248	// Sorting vector in ascending order with C function QSORT
4b016189	249	qsort((void*)zz,fNimf,sizeof(Float_t),comp);
28e0901a	250
5a881c97	251
5a881c97	252	// for(Int_t i=0; i<fNimf; i++){
28e0901a	253	// printf("\n After sorting -> zz[%d] = %f \n",i,zz[i]);
5a881c97	254	// }
28e0901a	255
28e0901a	256	// Rescale the maximum charge to fZmax
583731c7	257	for(j=0; j<fNimf; j++){
28e0901a	258	fZZ[j] = Int_t (zz[j]*fZmax/zz[fNimf-1]);
28e0901a	259	if(fZZ[j]<3) fZZ[j] = 3;
5a881c97	260	// printf("\n fZZ[%d] = %d \n",j,fZZ[j]);
28e0901a	261	}
4b016189	262
4b016189	263	delete[] zz;
28e0901a	264
5a881c97	265	// Check that the sum of the bound charges is not > than Zbound-Zalfa
28e0901a	266
28e0901a	267	for(Int_t ii=0; ii<fNimf; ii++){
86d81a3c	268	sumZ += fZZ[ii];
28e0901a	269	}
	270
	271	Int_t k = 0;
86d81a3c	272	if(sumZ>zbFrag){
583731c7	273	for(i=0; i< fNimf; i++){
28e0901a	274	k += 1;
86d81a3c	275	sumZ -= fZZ[i];
86d81a3c	276	if(sumZ<=zbFrag){
28e0901a	277	fNimf -= (i+1);
	278	break;
	279	}
	280	}
	281	}
	282	else {
86d81a3c	283	if(choice == 1) return;
86d81a3c	284	Int_t iDiff = Int_t((zbFrag-sumZ)/2);
28e0901a	285	if(iDiff<fNalpha){
	286	fNalpha=iDiff;
	287	return;
	288	}
	289	else{
	290	return;
	291	}
	292	}
	293
	294	fNimf += k;
583731c7	295	for(i=0; i<fNimf; i++){
28e0901a	296	fZZ[i] = fZZ[i+k];
	297	}
	298	fNimf -= k;
	299
86d81a3c	300	sumZ=0;
583731c7	301	for(i=0; i<fNimf; i++){
86d81a3c	302	sumZ += fZZ[i];
28e0901a	303	}
28e0901a	304
	305	}
	306
28e0901a	307	//_____________________________________________________________________________
92339a90	308	void AliZDCFragment::AttachNeutrons()
28e0901a	309	{
86d81a3c	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,
28e0901a	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};
86d81a3c	329	const Int_t kZIon[68]={1,1,2,3,3,
28e0901a	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;
5a881c97	346	// printf("\n fNimf=%d\n",fNimf);
5a881c97	347
28e0901a	348	for(Int_t i=0; i<fNimf; i++) {
28e0901a	349	for(Int_t j=0; j<68; j++) {
86d81a3c	350	iZ = kZIon[j];
28e0901a	351	if((fZZ[i]-iZ) == 0){
86d81a3c	352	iA = Int_t(kAIon[j]/0.93149432+0.5);
28e0901a	353	fNN[i] = iA - iZ;
28e0901a	354	break;
	355	}
	356	else if((fZZ[i]-iZ) < 0){
86d81a3c	357	fZZ[i] = kZIon[j-1];
	358	iA = Int_t (kAIon[j-1]/0.93149432+0.5);
	359	fNN[i] = iA - kZIon[j-1];
28e0901a	360	break;
	361	}
	362	}
	363	fZtot += fZZ[i];
	364	fNtot += fNN[i];
	365	}
	366
	367
	368	}
39202b85	369
39202b85	370	//_____________________________________________________________________________
3848c666	371	Float_t AliZDCFragment::DeuteronNumber()
39202b85	372	{
8a2624cc	373	// Calculates the fraction of deuterum nucleus produced
8a2624cc	374	//
3848c666	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;
39202b85	379	}