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

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