[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>
#include <search.h> 

// --- ROOT system
#include <TRandom.h>
#include <TF1.h>

// --- AliRoot classes
#include "AliZDCFragment.h"
   
ClassImp(AliZDCFragment)
 
//_____________________________________________________________________________
AliZDCFragment::AliZDCFragment()
{
  //
  // Default constructor
  //
  fB = 0;
}

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

//_____________________________________________________________________________
void AliZDCFragment::GenerateIMF(Int_t* fZZ, Int_t &fNalpha)
{
   // Coefficients of polynomial for average number of IMF
   Float_t  ParamNimf[5]={0.011236,1.8364,56.572,-116.24,58.289}; 
   // Coefficients of polynomial for fluctuations on average number of IMF
   Float_t  ParamFluctNimf[4]={-0.13176,2.9392,-5.2147,2.3092}; 
   // Coefficients of polynomial for average maximum Z of fragments
   Float_t  ParamZmax[4]={0.16899,14.203,-2.8284,65.036}; 
   // Coefficients of polynomial for fluctuations on maximum Z of fragments
   Float_t  ParamFluctZmax[5]={0.013782,-0.17282,1.5065,1.0654,-2.4317}; 
   // Coefficients of polynomial for exponent tau of fragments Z distribution
   Float_t  ParamTau[3]={6.7233,-15.85,13.047};  
   //Coefficients of polynomial for average number of alphas
   Float_t  ParamNalpha[4]={-0.68554,39.605,-68.311,30.165}; 
   // Coefficients of polynomial for fluctuations on average number of alphas
   Float_t  ParamFluctNalpha[5]={0.283,6.2141,-17.113,17.394,-6.6084}; 
   // Coefficients of function for Pb nucleus skin
   Float_t  ParamSkinPb[2]={0.93,11.05};
   
   // Thickness of nuclear surface
   Float_t  NuclearThick = 0.52;
   // Maximum impact parameter for U [r0*A**(1/3)]
   Float_t  bMaxU = 14.87;
   // Maximum impact parameter for Pb [r0*A**(1/3)]
   Float_t  bMaxPb = 14.22;
   // Z of the projectile
   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;
//   printf("\n fNimf = %d\n", fNimf);   
   
   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 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;

   Int_t Choice = 0;
   Float_t ZbFrag = 0, SumZ = 0.;
   // 2 possibilities:
   // 1) for bNorm < 0.9 ==> first remove alphas, then fragments
   // 2) for bNorm > 0.9 ==> first remove fragments, then alphas

   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[fNimf];
   for(Int_t j=0; j<fNimf; j++){
      zz[j] =0;
   }
   for(Int_t 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
 
//   int comp(const void*, const void*);
   qsort((void*)zz,fNimf,sizeof(float),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(Int_t 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]);
   }
   
   // Check that the sum of the bound charges is not greater than Zbound-Zalfa
   
   for(Int_t ii=0; ii<fNimf; ii++){
     SumZ += fZZ[ii];
   }
   
   Int_t k = 0;
   if(SumZ>ZbFrag){
     for(Int_t 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(Int_t i=0; i<fNimf; i++){
     fZZ[i] = fZZ[i+k];
   }
   fNimf -= k;
   
   SumZ=0;
   for(Int_t i=0; i<fNimf; i++){
     SumZ += fZZ[i];
   }
//   printf("\n	The END ->  fNimf = %d, SumZ = %f, fZmax = %f\n", 
//           fNimf, SumZ, fZmax);
//   printf("\n fNalpha = %d\n", fNalpha);   
//   for(Int_t j=0; j<fNimf; j++){
//     printf("\n 	fZZ[%d] = %d \n",j,fZZ[j]);
//   }
   
}

int comp(const void *i,const void *j){return (int*)j-(int*)i;}

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