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