[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
  //
  for(Int_t i=0; i<=99; i++){
     fZZ[i] = 0;
     fNN[i] = 0;
  }
}

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