[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  ParamNimf[5]={0.011236,1.8364,56.572,-116.24,58.289}; 
   // Coefficients of polynomial for fluctuations on average number of IMF
   const Float_t  ParamFluctNimf[4]={-0.13176,2.9392,-5.2147,2.3092}; 
   // Coefficients of polynomial for average maximum Z of fragments
   const Float_t  ParamZmax[4]={0.16899,14.203,-2.8284,65.036}; 
   // Coefficients of polynomial for fluctuations on maximum Z of fragments
   const Float_t  ParamFluctZmax[5]={0.013782,-0.17282,1.5065,1.0654,-2.4317}; 
   // Coefficients of polynomial for exponent tau of fragments Z distribution
   const Float_t  ParamTau[3]={6.7233,-15.85,13.047};  
   //Coefficients of polynomial for average number of alphas
   const Float_t  ParamNalpha[4]={-0.68554,39.605,-68.311,30.165}; 
   // Coefficients of polynomial for fluctuations on average number of alphas
   const Float_t  ParamFluctNalpha[5]={0.283,6.2141,-17.113,17.394,-6.6084}; 
   // Coefficients of function for Pb nucleus skin
   const Float_t  ParamSkinPb[2]={0.93,11.05};
   
   // Thickness of nuclear surface
   const Float_t  NuclearThick = 0.52;
   // Maximum impact parameter for U [r0*A**(1/3)]
   const Float_t  bMaxU = 14.87;
   // Maximum impact parameter for Pb [r0*A**(1/3)]
   const Float_t  bMaxPb = 14.22;
   // Z of the projectile
   const Float_t  ZProj = 82.;
   
   // From b(Pb) to b(U)
   Float_t  bU = fB*bMaxU/bMaxPb;
    
   // 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 = ZProj/92.*ZbU;
   
   // Zbound is proportional to b**2 up to b < bMaxPb-2*NuclearThick
   // and then it is an increasing exponential, imposing that at 
   // b=bMaxPb-2NuclearThick the two functions have the same derivative
   Float_t bCore = bMaxPb-2*NuclearThick;
   if(fB>bCore){
     fZbAverage=ZProj*(1.-TMath::Exp(-ParamSkinPb[0]*(fB-ParamSkinPb[1])));
   }
   if(fZbAverage>ZProj) fZbAverage = ZProj;
   Float_t ZbNorm = fZbAverage/ZProj;
   Float_t bNorm = fB/bMaxPb;
   
   // From Zbound to <Nimf>,<Zmax>,tau
   // Polinomial fits to Aladin distribution
   // --- A.Schuttauf et al, Nuc.Phys. A607 (1996) 457.
   Float_t AverageNimf = ParamNimf[0]+ParamNimf[1]*ZbNorm+ParamNimf[2]*
           TMath::Power(ZbNorm,2)+ParamNimf[3]*TMath::Power(ZbNorm,3)+
	   ParamNimf[4]*TMath::Power(ZbNorm,4);
   
   // Add fluctuation: from Singh et al. 
   Float_t FluctNimf = ParamFluctNimf[0]+ParamFluctNimf[1]*ZbNorm+
           ParamFluctNimf[2]*TMath::Power(ZbNorm,2)+ParamFluctNimf[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 = ParamZmax[0]+ParamZmax[1]*ZbNorm+ParamZmax[2]*
           TMath::Power(ZbNorm,2)+ParamZmax[3]*TMath::Power(ZbNorm,3);
   fTau = ParamTau[0]+ParamTau[1]*ZbNorm+ParamTau[2]*TMath::Power(ZbNorm,2);
   
   // Add fluctuation to mean value of Zmax (see Hubele)
   Float_t FluctZmax = ParamFluctZmax[0]+ParamFluctZmax[1]*ZbNorm+
           ParamFluctZmax[2]*TMath::Power(ZbNorm,2)+ParamFluctZmax[3]*
	   TMath::Power(ZbNorm,3)+ParamFluctZmax[4]*TMath::Power(ZbNorm,4);
   FluctZmax = FluctZmax*ZProj/6.;
   Float_t xg = gRandom->Gaus(0.0,1.0);
   FluctZmax = FluctZmax*xg;
   fZmax = AverageZmax+FluctZmax;
   if(fZmax>ZProj) fZmax = ZProj;
   
//   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 = ParamNalpha[0]+ParamNalpha[1]*ZbNorm+
           ParamNalpha[2]*TMath::Power(ZbNorm,2)+ParamNalpha[3]*
	   TMath::Power(ZbNorm,3);
   Float_t FluctAlpha = ParamFluctNalpha[0]+ParamFluctNalpha[1]*
           ZbNorm+ParamFluctNalpha[2]*TMath::Power(ZbNorm,2)+
	   ParamFluctNalpha[3]*TMath::Power(ZbNorm,3)+
	   ParamFluctNalpha[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)
{
   const Float_t AIon[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 ZIon[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 = ZIon[j];
	if((fZZ[i]-iZ) == 0){
	  iA = Int_t(AIon[j]/0.93149432+0.5);
	  fNN[i] = iA - iZ;
          break;
	}
	else if((fZZ[i]-iZ) < 0){
	  fZZ[i] = ZIon[j-1];
	  iA = Int_t (AIon[j-1]/0.93149432+0.5);
	  fNN[i] = iA - ZIon[j-1];
          break;
	}
      }
      fZtot += fZZ[i];
      fNtot += fNN[i];
   }		     
   

}
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
5a881c97	70	const Float_t ParamNimf[5]={0.011236,1.8364,56.572,-116.24,58.289};
28e0901a	71	// Coefficients of polynomial for fluctuations on average number of IMF
5a881c97	72	const Float_t ParamFluctNimf[4]={-0.13176,2.9392,-5.2147,2.3092};
28e0901a	73	// Coefficients of polynomial for average maximum Z of fragments
5a881c97	74	const Float_t ParamZmax[4]={0.16899,14.203,-2.8284,65.036};
28e0901a	75	// Coefficients of polynomial for fluctuations on maximum Z of fragments
5a881c97	76	const Float_t ParamFluctZmax[5]={0.013782,-0.17282,1.5065,1.0654,-2.4317};
28e0901a	77	// Coefficients of polynomial for exponent tau of fragments Z distribution
5a881c97	78	const Float_t ParamTau[3]={6.7233,-15.85,13.047};
28e0901a	79	//Coefficients of polynomial for average number of alphas
5a881c97	80	const Float_t ParamNalpha[4]={-0.68554,39.605,-68.311,30.165};
28e0901a	81	// Coefficients of polynomial for fluctuations on average number of alphas
5a881c97	82	const Float_t ParamFluctNalpha[5]={0.283,6.2141,-17.113,17.394,-6.6084};
28e0901a	83	// Coefficients of function for Pb nucleus skin
5a881c97	84	const Float_t ParamSkinPb[2]={0.93,11.05};
28e0901a	85
28e0901a	86	// Thickness of nuclear surface
5a881c97	87	const Float_t NuclearThick = 0.52;
28e0901a	88	// Maximum impact parameter for U [r0A*(1/3)]
5a881c97	89	const Float_t bMaxU = 14.87;
28e0901a	90	// Maximum impact parameter for Pb [r0A*(1/3)]
5a881c97	91	const Float_t bMaxPb = 14.22;
28e0901a	92	// Z of the projectile
5a881c97	93	const Float_t ZProj = 82.;
28e0901a	94
	95	// From b(Pb) to b(U)
	96	Float_t bU = fB*bMaxU/bMaxPb;
	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).
	104	Float_t ZbU = bUbUTMath::Pi()/7.48;
	105
	106	// Rescale Zbound for Pb
	107	fZbAverage = ZProj/92.*ZbU;
	108
	109	// Zbound is proportional to b*2 up to b < bMaxPb-2NuclearThick
	110	// and then it is an increasing exponential, imposing that at
	111	// b=bMaxPb-2NuclearThick the two functions have the same derivative
	112	Float_t bCore = bMaxPb-2*NuclearThick;
	113	if(fB>bCore){
	114	fZbAverage=ZProj(1.-TMath::Exp(-ParamSkinPb[0](fB-ParamSkinPb[1])));
	115	}
	116	if(fZbAverage>ZProj) fZbAverage = ZProj;
	117	Float_t ZbNorm = fZbAverage/ZProj;
	118	Float_t bNorm = fB/bMaxPb;
	119
	120	// From Zbound to <Nimf>,<Zmax>,tau
	121	// Polinomial fits to Aladin distribution
	122	// --- A.Schuttauf et al, Nuc.Phys. A607 (1996) 457.
	123	Float_t AverageNimf = ParamNimf[0]+ParamNimf[1]ZbNorm+ParamNimf[2]
	124	TMath::Power(ZbNorm,2)+ParamNimf[3]*TMath::Power(ZbNorm,3)+
	125	ParamNimf[4]*TMath::Power(ZbNorm,4);
	126
	127	// Add fluctuation: from Singh et al.
	128	Float_t FluctNimf = ParamFluctNimf[0]+ParamFluctNimf[1]*ZbNorm+
	129	ParamFluctNimf[2]*TMath::Power(ZbNorm,2)+ParamFluctNimf[3]
	130	*TMath::Power(ZbNorm,3);
	131	Float_t xx = gRandom->Gaus(0.0,1.0);
	132	FluctNimf = FluctNimf*xx;
	133	fNimf = Int_t(AverageNimf+FluctNimf);
	134	Float_t y = gRandom->Rndm();
	135	if(y < ((AverageNimf+FluctNimf)-fNimf)) fNimf += 1;
	136	if(fNimf ==0 && ZbNorm>0.75) fNimf = 1;
28e0901a	137
	138	Float_t AverageZmax = ParamZmax[0]+ParamZmax[1]ZbNorm+ParamZmax[2]
	139	TMath::Power(ZbNorm,2)+ParamZmax[3]*TMath::Power(ZbNorm,3);
	140	fTau = ParamTau[0]+ParamTau[1]ZbNorm+ParamTau[2]TMath::Power(ZbNorm,2);
	141
	142	// Add fluctuation to mean value of Zmax (see Hubele)
	143	Float_t FluctZmax = ParamFluctZmax[0]+ParamFluctZmax[1]*ZbNorm+
	144	ParamFluctZmax[2]TMath::Power(ZbNorm,2)+ParamFluctZmax[3]
	145	TMath::Power(ZbNorm,3)+ParamFluctZmax[4]*TMath::Power(ZbNorm,4);
	146	FluctZmax = FluctZmax*ZProj/6.;
	147	Float_t xg = gRandom->Gaus(0.0,1.0);
	148	FluctZmax = FluctZmax*xg;
	149	fZmax = AverageZmax+FluctZmax;
	150	if(fZmax>ZProj) fZmax = ZProj;
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
	159	Float_t AverageAlpha = ParamNalpha[0]+ParamNalpha[1]*ZbNorm+
	160	ParamNalpha[2]TMath::Power(ZbNorm,2)+ParamNalpha[3]
	161	TMath::Power(ZbNorm,3);
	162	Float_t FluctAlpha = ParamFluctNalpha[0]+ParamFluctNalpha[1]*
	163	ZbNorm+ParamFluctNalpha[2]*TMath::Power(ZbNorm,2)+
	164	ParamFluctNalpha[3]*TMath::Power(ZbNorm,3)+
	165	ParamFluctNalpha[4]*TMath::Power(ZbNorm,4);
	166	Float_t xxx = gRandom->Gaus(0.0,1.0);
	167	FluctAlpha = FluctAlpha*xxx;
	168	fNalpha = Int_t(AverageAlpha+FluctAlpha);
	169	Float_t yy = gRandom->Rndm();
	170	if(yy < ((AverageAlpha+FluctAlpha)-fNalpha)) fNalpha += 1;
	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
5a881c97	176	Int_t Choice = 0;
	177	Float_t ZbFrag = 0, SumZ = 0.;
	178
28e0901a	179	if(bNorm<=0.9) {
	180	// remove alpha from zbound to find zbound for fragments (Z>=3)
	181	ZbFrag = fZbAverage-fNalpha*2;
	182	Choice = 1;
	183	}
	184	else {
	185	ZbFrag = fZbAverage;
	186	Choice = 0;
	187	}
	188	// printf("\n Choice = %d, fZbAverage = %f, ZbFrag = %f \n", Choice, fZbAverage, ZbFrag);
	189
	190
	191	// Check if ZbFrag < fZmax
	192	if(ZbFrag<=fZmax) {
	193	if(fNimf>0 && ZbFrag>=2){
	194	fNimf = 1;
	195	fZZ[0] = Int_t(ZbFrag);
	196	SumZ = ZbFrag;
	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
	248	for(Int_t ii=0; ii<fNimf; ii++){
	249	SumZ += fZZ[ii];
	250	}
	251
	252	Int_t k = 0;
	253	if(SumZ>ZbFrag){
583731c7	254	for(i=0; i< fNimf; i++){
28e0901a	255	k += 1;
	256	SumZ -= fZZ[i];
	257	if(SumZ<=ZbFrag){
	258	fNimf -= (i+1);
	259	break;
	260	}
	261	}
	262	}
	263	else {
	264	if(Choice == 1) return;
	265	Int_t iDiff = Int_t((ZbFrag-SumZ)/2);
	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
	281	SumZ=0;
583731c7	282	for(i=0; i<fNimf; i++){
28e0901a	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	{
5a881c97	291	const Float_t AIon[68]={1.87612,2.80943,3.7284,5.60305,6.53536,
28e0901a	292	6.53622,8.39479,9.32699,10.2551,11.17793,
	293	13.04378,14.89917,17.6969,18.62284,21.41483,
	294	22.34193,25.13314,26.06034,28.85188,29.7818,
	295	32.57328,33.50356,36.29447,37.22492,41.87617,
	296	44.66324,47.45401,48.38228,51.17447,52.10307,
	297	54.89593,53.96644,58.61856,59.54963,68.85715,
	298	74.44178,78.16309,81.88358,83.74571,91.19832,
	299	98.64997,106.10997,111.68821,122.86796,
	300	128.45793,
	301	130.32111,141.51236,
	302	141.55,146.477,148.033,152.699,153.631,
	303	155.802,157.357,162.022,162.984,166.2624,
	304	168.554,171.349,173.4536,177.198,179.0518,
	305	180.675,183.473,188.1345,190.77,193.729,
	306	221.74295};
5a881c97	307	const Int_t ZIon[68]={1,1,2,3,3,
28e0901a	308	4,4,5,5,6,
	309	7,8,9,10,11,
	310	12,13,14,15,16,
	311	17,18,19,20,21,
	312	22,23,24,25,26,
	313	27,28,29,30,32,
	314	34,36,38,40,42,
	315	46,48,50,54,56,
	316	58,62,
	317	63,64,65,66,67,
	318	68,69,70,71,72,
	319	73,74,75,76,77,
	320	78,79,80,81,82,
	321	92};
	322
	323	Int_t iZ, iA;
5a881c97	324	// printf("\n fNimf=%d\n",fNimf);
5a881c97	325
28e0901a	326	for(Int_t i=0; i<fNimf; i++) {
	327	for(Int_t j=0; j<68; j++) {
	328	iZ = ZIon[j];
	329	if((fZZ[i]-iZ) == 0){
	330	iA = Int_t(AIon[j]/0.93149432+0.5);
	331	fNN[i] = iA - iZ;
28e0901a	332	break;
	333	}
	334	else if((fZZ[i]-iZ) < 0){
	335	fZZ[i] = ZIon[j-1];
	336	iA = Int_t (AIon[j-1]/0.93149432+0.5);
	337	fNN[i] = iA - ZIon[j-1];
28e0901a	338	break;
	339	}
	340	}
	341	fZtot += fZZ[i];
	342	fNtot += fNN[i];
	343	}
	344
	345
	346	}