[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()
{
  //
  // Default constructor
  //
  fB = 0;
}

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

//_____________________________________________________________________________
void AliZDCFragment::GenerateIMF(Int_t* fZZ, Int_t &fNalpha)
{

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