[u/mrichter/AliRoot.git] / ZDC / AliZDCFragment.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/


// ******************************************************************
//
//  	Class for nuclear fragments formation
//
// ******************************************************************

// --- Standard libraries
#include <stdlib.h>

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

// --- AliRoot classes
#include "AliZDCFragment.h"
 
ClassImp(AliZDCFragment)
   
int comp(const void *i,const void *j) {return *(int *)i - *(int *)j;}


//_____________________________________________________________________________
AliZDCFragment::AliZDCFragment():
  fB(0),
  fZbAverage(0),
  fNimf(0),
  fZmax(0),
  fTau(0),
  fNalpha(0),
  fZtot(0),
  fNtot(0)
{
  //
  // Default constructor
  //
}

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

//_____________________________________________________________________________
void AliZDCFragment::GenerateIMF()
{

   // Loop variables
  Int_t i,j;

   // Coefficients of polynomial for average number of IMF
   const Float_t  kParamNimf[5]={0.011236,1.8364,56.572,-116.24,58.289}; 
   // Coefficients of polynomial for fluctuations on average number of IMF
   const Float_t  kParamFluctNimf[4]={-0.13176,2.9392,-5.2147,2.3092}; 
   // Coefficients of polynomial for average maximum Z of fragments
   //const Float_t  kParamZmax[4]={0.16899,14.203,-2.8284,65.036}; 
   const Float_t  kParamZmax[4]={0.16899,14.203,-2.8284,70.5}; 
   // 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+4*kNuclearThick;
   // Z of the projectile
   const Float_t  kZProj = 82.;
   
   // From b(Pb) to b(U)
   if(fB>kbMaxPb) fB = 2*kbMaxPb-fB;
   
   Float_t  bU = fB*kbMaxU/kbMaxPb;
    
   // From b(U) to Zbound(U) 
   // --- A.Schuttauf et al, Nuc.Phys. A607 (1996) 457 ---------------
   // From geometrical consideration and from dsigma/dZbound for U+U,
   // which is approx. constant, the constant value is found  
   // integrating the nucleus cross surface from 0 to bmax=R1+R2 where 
   // R = 1.2*A**(1/3). This value has been measured in Aladin (U+U).
   Float_t  zbU = bU*bU*TMath::Pi()/7.48;
   
   //  Rescale Zbound for Pb
   fZbAverage = kZProj/92.*zbU;
   
   // Zbound is proportional to b**2 up to b < kbMaxPb-2*kNuclearThick
   // and then it is an increasing exponential, imposing that at 
   // b=kbMaxPb-2kNuclearThick the two functions have the same derivative
   /*Float_t bCore = kbMaxPb-2*kNuclearThick;
   if(fB>bCore){
     fZbAverage=kZProj*(1.-TMath::Exp(-kParamSkinPb[0]*(fB-kParamSkinPb[1])));
   }*/
   if(fZbAverage>kZProj) fZbAverage = kZProj;
   Float_t zbNorm = fZbAverage/kZProj;
   Float_t bNorm = fB/kbMaxPb;
   
   // From Zbound to <Nimf>,<Zmax>,tau
   // Polinomial fits to Aladin distribution
   // --- A.Schuttauf et al, Nuc.Phys. A607 (1996) 457.
   Float_t averageNimf = kParamNimf[0]+kParamNimf[1]*zbNorm+kParamNimf[2]*
           TMath::Power(zbNorm,2)+kParamNimf[3]*TMath::Power(zbNorm,3)+
	   kParamNimf[4]*TMath::Power(zbNorm,4);
   
   // Add fluctuation: from Singh et al. 
   Float_t fluctNimf = kParamFluctNimf[0]+kParamFluctNimf[1]*zbNorm+
           kParamFluctNimf[2]*TMath::Power(zbNorm,2)+kParamFluctNimf[3]
	   *TMath::Power(zbNorm,3);
   Float_t xx = gRandom->Gaus(0.0,1.0);
   fluctNimf = fluctNimf*xx;
   fNimf = Int_t(averageNimf+fluctNimf);
   Float_t y = gRandom->Rndm();
   if(y < ((averageNimf+fluctNimf)-fNimf)) fNimf += 1;
   if(fNimf ==0 && zbNorm>0.75) fNimf = 1;
   
   Float_t averageZmax = kParamZmax[0]+kParamZmax[1]*zbNorm+kParamZmax[2]*
           TMath::Power(zbNorm,2)+kParamZmax[3]*TMath::Power(zbNorm,3);
   fTau = kParamTau[0]+kParamTau[1]*zbNorm+kParamTau[2]*TMath::Power(zbNorm,2);
   
   // Add fluctuation to mean value of Zmax (see Hubele)
   Float_t fluctZmax = kParamFluctZmax[0]+kParamFluctZmax[1]*zbNorm+
           kParamFluctZmax[2]*TMath::Power(zbNorm,2)+kParamFluctZmax[3]*
	   TMath::Power(zbNorm,3)+kParamFluctZmax[4]*TMath::Power(zbNorm,4);
   fluctZmax = fluctZmax*kZProj/6.;
   Float_t xg = gRandom->Gaus(0.0,1.0);
   fluctZmax = fluctZmax*xg;
   fZmax = (averageZmax+fluctZmax);
   if(fZmax>kZProj) fZmax = kZProj;
   
//   printf("\n\n ------------------------------------------------------------");   
//   printf("\n Generation of nuclear fragments\n");   
//   printf("\n fNimf = %d\n", fNimf);   
//   printf("\n fZmax = %f\n", fZmax); 

   // Find the number of alpha particles 
   // from Singh et al. : Pb+emulsion
   Float_t averageAlpha = kParamNalpha[0]+kParamNalpha[1]*zbNorm+
           kParamNalpha[2]*TMath::Power(zbNorm,2)+kParamNalpha[3]*
	   TMath::Power(zbNorm,3);
   Float_t fluctAlpha = kParamFluctNalpha[0]+kParamFluctNalpha[1]*
           zbNorm+kParamFluctNalpha[2]*TMath::Power(zbNorm,2)+
	   kParamFluctNalpha[3]*TMath::Power(zbNorm,3)+
	   kParamFluctNalpha[4]*TMath::Power(zbNorm,4);
   Float_t xxx = gRandom->Gaus(0.0,1.0);
   fluctAlpha = fluctAlpha*xxx;
   fNalpha = Int_t(averageAlpha+fluctAlpha);
   Float_t yy = gRandom->Rndm();
   if(yy < ((averageAlpha+fluctAlpha)-fNalpha)) fNalpha += 1;

   // 2 possibilities:
   // 1) for bNorm < 0.9 ==> first remove alphas, then fragments
   // 2) for bNorm > 0.9 ==> first remove fragments, then alphas

   Int_t choice = 0;
   Float_t zbFrag = 0, sumZ = 0.;

   if(bNorm<=0.9) {
   // remove alpha from zbound to find zbound for fragments  (Z>=3)
     zbFrag = fZbAverage-fNalpha*2;
     choice = 1;
   }
   else {
     zbFrag = fZbAverage;
     choice = 0;
   }
//   printf("\n choice = %d, fZbAverage = %f, zbFrag = %f \n", choice, fZbAverage, zbFrag);
   
   
   // Check if zbFrag < fZmax
   if(zbFrag<=fZmax) {
     if(fNimf>0 && zbFrag>=2){
       fNimf = 1;
       fZZ[0] = Int_t(zbFrag);
       sumZ = zbFrag;
     }
     else {
       fNimf = 0;
     }
     return;
   }
   
   // Prepare the exponential charge distribution dN/dZ
   if(fZmax <= 0.01) {
     fNimf = 0;
     return;
   }
   if(fNimf == 0) {
     fNimf = 0;
     return;
   }
   
   TF1 *funTau = new TF1("funTau","1./(x**[0])",0.01,fZmax);
   funTau->SetParameter(0,fTau);

   // Extract randomly the charge of the fragments from the distribution
 
   Float_t * zz = new Float_t[fNimf];
   for(j=0; j<fNimf; j++){
      zz[j] =0;
   }
   for(i=0; i<fNimf; i++){
      zz[i] = Float_t(funTau->GetRandom());
//      printf("\n	zz[%d] = %f \n",i,zz[i]);
   }
   delete funTau;
   
   // Sorting vector in ascending order with C function QSORT 
   qsort((void*)zz,fNimf,sizeof(Float_t),comp);

   
//   for(Int_t i=0; i<fNimf; i++){
//      printf("\n After sorting -> zz[%d] = %f \n",i,zz[i]);
//   }
   
   // Rescale the maximum charge to fZmax
   for(j=0; j<fNimf; j++){
     fZZ[j] = Int_t (zz[j]*fZmax/zz[fNimf-1]);
     if(fZZ[j]<3) fZZ[j] = 3;
//     printf("\n 	fZZ[%d] = %d \n",j,fZZ[j]);
   }

   delete[] zz;
   
   // Check that the sum of the bound charges is not > than Zbound-Zalfa
   
   for(Int_t ii=0; ii<fNimf; ii++){
     sumZ += fZZ[ii];
   }
   
   Int_t k = 0;
   if(sumZ>zbFrag){
     for(i=0; i< fNimf; i++){
       k += 1;
       sumZ -= fZZ[i];
       if(sumZ<=zbFrag){
         fNimf -= (i+1);
         break;
       }
     }
   }
   else {
     if(choice == 1) return;
     Int_t iDiff = Int_t((zbFrag-sumZ)/2);
     if(iDiff<fNalpha){
       fNalpha=iDiff;
       return;
     }
     else{
       return;
     }
   }

   fNimf += k;
   for(i=0; i<fNimf; i++){
     fZZ[i] = fZZ[i+k];
   }
   fNimf -= k;
   
   sumZ=0;
   for(i=0; i<fNimf; i++){
     sumZ += fZZ[i];
   }
   
}

//_____________________________________________________________________________
void AliZDCFragment::AttachNeutrons()
{
//
// Prepare nuclear fragment by attaching a suitable number of neutrons
//
   const Float_t kAIon[68]={1.87612,2.80943,3.7284,5.60305,6.53536,
     		     6.53622,8.39479,9.32699,10.2551,11.17793,
     		     13.04378,14.89917,17.6969,18.62284,21.41483,
     		     22.34193,25.13314,26.06034,28.85188,29.7818,
     		     32.57328,33.50356,36.29447,37.22492,41.87617,
     		     44.66324,47.45401,48.38228,51.17447,52.10307,
     		     54.89593,53.96644,58.61856,59.54963,68.85715,
     		     74.44178,78.16309,81.88358,83.74571,91.19832,
     		     98.64997,106.10997,111.68821,122.86796,
     		     128.45793,
     		     130.32111,141.51236,
     		     141.55,146.477,148.033,152.699,153.631,
     		     155.802,157.357,162.022,162.984,166.2624,
     		     168.554,171.349,173.4536,177.198,179.0518,
     		     180.675,183.473,188.1345,190.77,193.729,
     		     221.74295};
   const Int_t kZIon[68]={1,1,2,3,3,
     		     4,4,5,5,6,
     		     7,8,9,10,11,
     		     12,13,14,15,16,
     		     17,18,19,20,21,
     		     22,23,24,25,26,
     		     27,28,29,30,32,
     		     34,36,38,40,42,
     		     46,48,50,54,56,
     		     58,62,
     		     63,64,65,66,67,
     		     68,69,70,71,72,
     		     73,74,75,76,77,
     		     78,79,80,81,82,
      		     92};
    
   Int_t iZ, iA;  
//   printf("\n fNimf=%d\n",fNimf);  

   for(Int_t i=0; i<fNimf; i++) {
      for(Int_t j=0; j<68; j++) {
        iZ = kZIon[j];
	if((fZZ[i]-iZ) == 0){
	  iA = Int_t(kAIon[j]/0.93149432+0.5);
	  fNN[i] = iA - iZ;
          break;
	}
	else if((fZZ[i]-iZ) < 0){
	  fZZ[i] = kZIon[j-1];
	  iA = Int_t (kAIon[j-1]/0.93149432+0.5);
	  fNN[i] = iA - kZIon[j-1];
          break;
	}
      }
      fZtot += fZZ[i];
      fNtot += fNN[i];
   }		     
   

}

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