coverity bug fix
[u/mrichter/AliRoot.git] / TUHKMgen / UHKM / InitialStateHydjet.cxx
1 //expanding localy equilibated fireball with volume hadron radiation
2 //thermal part: Blast wave model, Bjorken-like parametrization
3 //hyght-pt: PYTHIA + jet quenching model PYQUEN
4 //                                                                           
5 //         HYDJET++ 
6 //         version 1.0:  
7 //         InitialStateHydjet is the modified InitialStateBjorken 
8 //         The high-pt part related with PYTHIA-PYQUEN is included       
9 //         InitialStateBjorken (FASTMC) was used.
10 //
11 //
12 //         
13 //         InitialStateBjorken           
14 //         version 2.0: 
15 //         Ludmila Malinina  malinina@lav01.sinp.msu.ru,   SINP MSU/Moscow and JINR/Dubna
16 //         Ionut Arsene  i.c.arsene@fys.uio.no,            Oslo University                                                
17 //                     June 2007
18 //        
19 //         version 1.0:                                                               
20 //         Nikolai Amelin, Ludmila Malinina, Timur Pocheptsov (C) JINR/Dubna
21 //         amelin@sunhe.jinr.ru, malinina@sunhe.jinr.ru, pocheptsov@sunhe.jinr.ru 
22 //                           November. 2, 2005 
23
24                      
25
26
27 #include <iostream> 
28 #include <fstream>
29
30 #include <TLorentzVector.h>
31 #include <TVector3.h>
32 #include <TMath.h>
33
34 #include "InitialStateHydjet.h"
35 #include "RandArrayFunction.h"
36 #include "GrandCanonical.h"
37 #include "StrangePotential.h"
38 #include "Particle.h"
39 #include "ParticlePDG.h"
40 #include "HYJET_COMMONS.h"
41
42 extern "C" void  hyevnt_();
43 extern "C" void  myini_();
44 extern HYIPARCommon HYIPAR;
45 extern HYFPARCommon HYFPAR;
46 extern HYJPARCommon HYJPAR;
47 extern HYPARTCommon HYPART;
48 extern SERVICECommon SERVICE;
49
50 using std::cout;
51 using std::endl;
52
53 class ParticleAllocator;
54 class TRandom3;
55
56 // declaration of the static member fLastIndex
57 Int_t Particle::fgLastIndex;
58
59 //_________________________________________________________________________________
60 void InitialStateHydjet::Initialize(List_t &source, ParticleAllocator & allocator) {
61   // Generate initial particles from the soft and hard components
62
63   // Initialize the static "last index variable"
64   Particle::InitIndexing(); 
65
66   //----- high-pt part------------------------------
67   TLorentzVector partJMom, partJPos, zeroVec;
68
69   // run a HYDJET event
70   hyevnt_(); 
71
72   fBgen  =  HYFPAR.bgen * HYIPAR.RA;
73   fNpart =  HYFPAR.npart;
74   fNcoll =  HYFPAR.nbcol;
75
76   
77       
78   if(fParams.fNhsel != 0) {   
79     //get number of particles in jets
80     Int_t numbJetPart = HYPART.njp;
81
82     for(Int_t i = 0; i <numbJetPart; i++) {
83       Int_t pdg = Int_t(HYPART.ppart[i][1]);
84       Double_t px = HYPART.ppart[i][2];
85       Double_t py = HYPART.ppart[i][3];
86       Double_t pz = HYPART.ppart[i][4];
87       Double_t e =  HYPART.ppart[i][5];
88       Double_t vx = HYPART.ppart[i][6];
89       Double_t vy = HYPART.ppart[i][7];
90       Double_t vz = HYPART.ppart[i][8];
91       Double_t vt = HYPART.ppart[i][9];    
92       ParticlePDG *partDef = fDatabase->GetPDGParticle(pdg);
93       Int_t type =1;                //from jet
94       if(partDef) {
95         partJMom.SetXYZT(px, py, pz, e);
96         partJPos.SetXYZT(vx, vy, vz, vt);
97         Particle *particle=new Particle(partDef, partJPos, partJMom, 0, 0, type, -1, zeroVec, zeroVec);
98         particle->SetIndex();
99         allocator.AddParticle(*particle, source);
100         delete particle;
101       }
102     }
103   }       //nhsel !=0 not only hydro!             
104
105          
106   //----------HYDRO part------------------------------------------------
107   if(fParams.fNhsel < 3) {
108     const Double_t  weightMax = 2*TMath::CosH(fParams.fUmax);
109     const Int_t nBins = 100;
110     Double_t probList[nBins];
111     RandArrayFunction arrayFunctDistE(nBins);
112     RandArrayFunction arrayFunctDistR(nBins);
113     
114     TLorentzVector partPos, partMom, n1, p0;
115     TVector3 vec3;
116     //set maximal hadron energy
117     const Double_t eMax = 5.;  
118     //-------------------------------------
119     // get impact parameter    
120     
121     //effective volume for central     
122     double dYl= 2 * fParams.fYlmax; //uniform distr. [-Ylmax; Ylmax]  
123     if(fParams.fEtaType >0) dYl = TMath::Sqrt(2 * TMath::Pi()) * fParams.fYlmax ;  //Gaussian distr. 
124     Double_t volEffcent = 2 * TMath::Pi() * fParams.fTau * dYl * 
125     (fParams.fR * fParams.fR)/TMath::Power((fParams.fUmax),2)*
126     ((fParams.fUmax)*TMath::SinH((fParams.fUmax))-TMath::CosH((fParams.fUmax))+ 1);
127  
128     //effective volume for non-central Simpson2 
129     Double_t volEffnoncent = fParams.fTau * dYl * SimpsonIntegrator2(0., 2.*TMath::Pi());
130     fVolEff = volEffcent * HYFPAR.npart/HYFPAR.npart0;
131
132     Double_t coeffRB = TMath::Sqrt(volEffcent * HYFPAR.npart/HYFPAR.npart0/volEffnoncent);
133     Double_t coeffR1 = HYFPAR.npart/HYFPAR.npart0;
134     coeffR1 = TMath::Power(coeffR1, 0.333333);
135
136     double veff=fVolEff;
137     
138     //------------------------------------
139     //cycle on particles types
140     for(Int_t i = 0; i < fParams.fNPartTypes; ++i) {
141       Double_t mparam = fParams.fPartMult[2 * i] * veff;
142       Int_t multiplicity = gRandom->Poisson(mparam);
143       const Int_t encoding = fParams.fPartEnc[i];
144
145       if(multiplicity > 0) {
146         ParticlePDG *partDef = fDatabase->GetPDGParticle(encoding);
147         if(!partDef) {
148           Error("InitialStateHydjet::Initialize", "No particle with encoding %d", encoding);
149           continue;
150         }
151         //no charm now !
152         if(TMath::Abs(partDef->GetCharmQNumber())>0 || TMath::Abs(partDef->GetCharmAQNumber())>0){
153           continue;
154         }
155
156         //compute chemical potential for single f.o. mu==mu_ch
157         //compute chemical potential for thermal f.o.                
158         Double_t mu = fParams.fPartMu[2 * i];
159         
160         //choose Bose-Einstein or Fermi-Dirac statistics
161         const Double_t d    = !(Int_t(2*partDef->GetSpin()) & 1) ? -1 : 1;
162         const Double_t mass = partDef->GetMass();                        
163         
164         //prepare histogram to sample hadron energy: 
165         Double_t h = (eMax - mass) / nBins;
166         Double_t x = mass + 0.5 * h;
167         Int_t ii; 
168         for(ii = 0; ii < nBins; ++ii) {
169           if(x>=mu && fParams.fThFO>0)probList[ii] = x * TMath::Sqrt(x * x - mass * mass) / 
170                                         (TMath::Exp((x - mu) / (fParams.fThFO)) + d);
171           if(x>=mu && fParams.fThFO<=0)probList[ii] = x * TMath::Sqrt(x * x - mass * mass) / 
172                                          (TMath::Exp((x - mu) / (fParams.fT)) + d);                                                              
173           if(x<mu)probList[ii] = 0.; 
174           x += h;
175         }
176         arrayFunctDistE.PrepareTable(probList);
177
178         //prepare histogram to sample hadron transverse radius: 
179         h = (fParams.fR) / nBins;
180         x =  0.5 * h;
181         Double_t param = (fParams.fUmax) / (fParams.fR);
182         for(ii = 0; ii < nBins; ++ii) {
183           probList[ii] = x * TMath::CosH(param*x);
184           x += h;
185         }
186         arrayFunctDistR.PrepareTable(probList);
187
188         //loop over hadrons, assign hadron coordinates and momenta
189         Double_t weight = 0., yy = 0., px0 = 0., py0 = 0., pz0 = 0.;
190         Double_t e = 0., x0 = 0., y0 = 0., z0 = 0., t0 = 0., etaF = 0.; 
191         Double_t r, rB, phiF;
192       
193         for(Int_t j = 0; j < multiplicity; ++j) {               
194           do {
195             fParams.fEtaType <=0 ? etaF = fParams.fYlmax * (2. * gRandom->Rndm() - 1.) 
196               : etaF = (fParams.fYlmax) * (gRandom->Gaus());                                             
197             n1.SetXYZT(0.,0.,TMath::SinH(etaF),TMath::CosH(etaF));  
198             if(TMath::Abs(etaF)>5.)continue;
199             
200             rB = fParams.fR * coeffRB * coeffR1;
201                 
202             Double_t rho = TMath::Sqrt(gRandom->Rndm());
203             Double_t phi = TMath::TwoPi() * gRandom->Rndm();
204             Double_t rx =  TMath::Sqrt(1-fParams.fEpsilon)*rB; 
205             Double_t ry =  TMath::Sqrt(1+fParams.fEpsilon)*rB;
206             
207             x0 = rx * rho * TMath::Cos(phi);
208             y0 = ry * rho * TMath::Sin(phi);
209             r = TMath::Sqrt(x0*x0+y0*y0);
210             phiF = TMath::Abs(TMath::ATan(y0/x0));
211             
212             if(x0<0&&y0>0)phiF = TMath::Pi()-phiF;
213             if(x0<0&&y0<0)phiF = TMath::Pi()+phiF;
214             if(x0>0&&y0<0)phiF = 2.*TMath::Pi()-phiF;
215             
216             //proper time with emission duration                                                               
217             Double_t tau = coeffR1 * fParams.fTau +  sqrt(2.) * fParams.fSigmaTau * coeffR1 * (gRandom->Gaus());          
218             z0 = tau  * TMath::SinH(etaF);                                                                             
219             t0 = tau  * TMath::CosH(etaF);
220           
221             Double_t rhou = fParams.fUmax * r / rB;
222
223         
224             Double_t uxf = TMath::SinH(rhou)*TMath::Sqrt(1+fParams.fDelta)*TMath::Cos(phiF); 
225             Double_t uyf = TMath::SinH(rhou)*TMath::Sqrt(1-fParams.fDelta)*TMath::Sin(phiF);
226             Double_t utf = TMath::CosH(etaF) * TMath::CosH(rhou) * 
227               TMath::Sqrt(1+fParams.fDelta*TMath::Cos(2*phiF)*TMath::TanH(rhou)*TMath::TanH(rhou));
228             Double_t uzf = TMath::SinH(etaF) * TMath::CosH(rhou) * 
229               TMath::Sqrt(1+fParams.fDelta*TMath::Cos(2*phiF)*TMath::TanH(rhou)*TMath::TanH(rhou));
230         
231             vec3.SetXYZ(uxf / utf, uyf / utf, uzf / utf);
232             n1.Boost(-vec3); 
233     
234             yy = weightMax * gRandom->Rndm();        
235                
236             Double_t php0 = TMath::TwoPi() * gRandom->Rndm();
237             Double_t ctp0 = 2. * gRandom->Rndm() - 1.;
238             Double_t stp0 = TMath::Sqrt((1.-ctp0)*(1.+ctp0)); 
239             e = mass + (eMax - mass) * arrayFunctDistE(); 
240             Double_t pp0 = TMath::Sqrt((e-mass)*(e+mass));
241             px0 = pp0 * stp0 * TMath::Sin(php0); 
242             py0 = pp0 * stp0 * TMath::Cos(php0);
243             pz0 = pp0 * ctp0;
244             p0.SetXYZT(px0, py0, pz0, e);
245             
246             //weight for rdr          
247             weight = (n1 * p0) /e;  // weight for rdr gammar: weight = (n1 * p0) / n1[3] / e; 
248           } while(yy >= weight); 
249         
250           partMom.SetXYZT(px0, py0, pz0, e);
251           partPos.SetXYZT(x0, y0, z0, t0);
252           partMom.Boost(vec3);
253           Int_t type =0; //hydro
254
255           Particle *particle=new Particle(partDef, partPos, partMom, 0., 0, type, -1, zeroVec, zeroVec);
256           particle->SetIndex();
257           allocator.AddParticle(*particle, source);
258           delete particle;
259         } //nhsel==4 , no hydro part
260       }
261     } 
262   }
263   
264 }
265
266 //_________________________________________________________________________________
267 Bool_t InitialStateHydjet::ReadParams() {     
268   // Read parameters from an input file in ascii 
269  
270   Float_t par[200] = {0.};
271   Int_t i = 0; 
272   std::string s(40,' '); 
273   std::ifstream input("RunInputHydjet");
274   if (!input) {
275     Error("Ukm::ReadParams", "Cannot open RunInputHydjet");
276     return kFALSE;
277   }
278   
279   while (std::getline(input, s)) {
280     input>>par[i];
281     if (i < 140) 
282       std::cout<<s<<"     =  "<<par[i]<<std::endl;
283     ++i;
284     std::getline(input,s);
285   }
286   
287   std::cout<<"\nFor output use the files RunOutput.root  \n\n"<< std::endl; 
288   
289   fParams.fNevnt  = Int_t(par[0]); //number of events
290   fParams.fSqrtS  = par[1];        //cms energy per nucleon
291   fParams.fAw     = par[2];        // atomic number of colliding nuclei
292   fParams.fIfb    = Int_t(par[3]);      // flag of type of centrality generation (=0 is fixed by fBfix, not 0 
293                                          //impact parameter is generated in each event between fBfmin 
294                                           //and fBmax according with Glauber model (f-la 30)
295   fParams.fBmin = par[4];         //minimum impact parameter in units of nuclear radius RA 
296   fParams.fBmax = par[5];         //maximum impact parameter in units of nuclear radius RA
297   fParams.fBfix = par[6];         //fix impact parameter in units of nuclear radius RA
298   
299   fParams.fSeed = Int_t(par[7]);         //parameter to set the random nuber seed (=0 the current time is used
300                                    //to set the random generator seed, !=0 the value fSeed is 
301                                    //used to set the random generator seed and then the state of random
302                                    //number generator in PYTHIA MRPY(1)=fSeed
303        
304   fParams.fT         = par[8];     //chemical freeze-out temperature in GeV    
305   fParams.fMuB       = par[9];     //baryon potential 
306   fParams.fMuS       = par[10];    //strangeness potential 
307   fParams.fMuI3      = par[11];    //isospin potential   
308   fParams.fThFO      = par[12];    //thermal freeze-out temperature T^th in GeV
309   fParams.fMu_th_pip = par[13];    // effective chemical potential of positivly charged pions at thermal in GeV 
310
311        
312   fParams.fTau       = par[14];     //proper time value
313   fParams.fSigmaTau  = par[15];     //its standart deviation (emission duration)
314   fParams.fR         = par[16];     //maximal transverse radius 
315   fParams.fYlmax     = par[17];     //maximal longitudinal rapidity 
316   fParams.fUmax      = par[18];     //maximal transverse velocity multiplaed on \gamma_r 
317   fParams.fDelta     = par[19];     //momentum asymmetry parameter
318   fParams.fEpsilon   = par[20];     //coordinate asymmetry parameter
319   
320   fParams.fDecay      = Int_t(par[21]);    // flag to switch on/off hadron decays<0: decays off,>=0: decays on, (default: 0)
321   fParams.fWeakDecay  = Int_t(par[22]);    //flag to switch on/off weak hadron decays <0: decays off, >0: decays on, (default: 0)
322   
323   fParams.fEtaType   = Int_t(par[23]);     // flag to choose rapidity distribution, if fEtaType<=0, 
324                                      //then uniform rapidity distribution in [-fYlmax,fYlmax] if fEtaType>0,
325                                      //then Gaussian with dispertion = fYlmax 
326   
327   fParams.fTMuType   = Int_t(par[24]);     // flag to use calculated chemical freeze-out temperature,
328                                      //baryon potential and strangeness potential as a function of fSqrtS 
329
330   fParams.fCorrS     = par[25];     // flag and value to include strangeness supression factor    
331   fParams.fNhsel = Int_t(par[26]);         //flag to switch on/off jet and hydro-state production (0: jet
332                                      // production off and hydro on, 1: jet production on and jet quenching
333                                      // off and hydro on, 2: jet production on and jet quenching on and
334                                      // hydro on, 3: jet production on and jet quenching off and hydro
335                                      // off, 4: jet production on and jet quenching on and hydro off
336  
337   fParams.fIshad= Int_t(par[27]);         //flag to switch on/off impact parameter dependent nuclear
338                                     // shadowing for gluons and light sea quarks (u,d,s) (0: shadowing off,
339                                     // 1: shadowing on for fAw=207, 197, 110, 40, default: 1
340   
341   fParams.fPtmin = par[28];       //minimal transverse momentum transfer p_T of hard
342                                    // parton-parton scatterings in GeV (the PYTHIA parameter ckin(3)=fPtmin)
343    
344   //  PYQUEN energy loss model parameters:
345  
346   fParams.fT0 = par[29];          // initial temperature (in GeV) of QGP for
347                                    //central Pb+Pb collisions at mid-rapidity (initial temperature for other
348                                   //centralities and atomic numbers will be calculated automatically) (allowed range is 0.2<fT0<2) 
349   
350   fParams.fTau0= par[30];        //proper QGP formation time in fm/c (0.01<fTau0<10)
351   fParams.fNf= Int_t(par[31]);          //number of active quark flavours N_f in QGP fNf=0, 1,2 or 3 
352   fParams.fIenglu= Int_t(par[32]);      // flag to fix type of in-medium partonic energy loss 
353                                            //(0: radiative and collisional loss, 1: radiative loss only, 2:
354                                   //collisional loss only) (default: 0);
355   fParams.fIanglu= Int_t(par[33]);      //flag to fix type of angular distribution of in-medium emitted
356                                   // gluons (0: small-angular, 1: wide-angular, 2:collinear) (default: 0).
357
358
359   //PYTHIA parameters:
360   Int_t jj;
361   for (Int_t j = 0; j <25; ++j) {
362     jj= 35+j;
363     SERVICE.parPYTH[j]=par[jj];
364   } 
365
366   // Set Random Number seed 
367         
368   gRandom->SetSeed(fParams.fSeed); //Set 0 to use the current time
369   //to send seed in PYTHIA
370   SERVICE.iseed_fromC=gRandom->GetSeed(); 
371   std::cout<<"Seed for random number generation= "<<gRandom->GetSeed()<<std::endl;  
372   
373   fParams.fNPartTypes = 0;         //counter of hadron species
374   return kTRUE; 
375 }
376
377 //_________________________________________________________________________________
378 Bool_t InitialStateHydjet::MultIni() {
379   // Calculate average multiplicities, chemical potentials (if necessary),
380   // initialize pyquen 
381
382   //check and redefine input parameters
383   if(fParams.fTMuType>0 &&  fParams.fSqrtS > 2.24) {
384     if(fParams.fSqrtS < 2.24){
385       Error("InitialStateHydjet::MultIni", "SqrtS<2.24 not allowed with fParams.fTMuType>0");
386       return 0;
387     }
388     
389     //sqrt(s) = 2.24 ==> T_kin = 0.8 GeV
390     //see J. Cleymans, H. Oeschler, K. Redlich,S. Wheaton, Phys Rev. C73 034905 (2006)
391     fParams.fMuB = 1.308/(1. + fParams.fSqrtS*0.273);
392     fParams.fT = 0.166 - 0.139*fParams.fMuB*fParams.fMuB - 0.053*fParams.fMuB*fParams.fMuB*
393       fParams.fMuB*fParams.fMuB;
394     fParams.fMuI3 = 0.;
395     fParams.fMuS = 0.;
396     //create strange potential object and set strangeness density 0
397     StrangePotential psp(0., fDatabase);
398     psp.SetBaryonPotential(fParams.fMuB);
399     psp.SetTemperature(fParams.fT);
400     //compute strangeness potential
401     if(fParams.fMuB > 0.01)
402       fParams.fMuS = psp.CalculateStrangePotential();
403     
404     //if user choose fYlmax larger then allowed by kinematics at the specified beam energy sqrt(s)     
405     if(fParams.fYlmax > TMath::Log(fParams.fSqrtS/0.94)){
406       Error("InitialStateHydjet::MultIni", "fParams.fYlmax > TMath::Log(fParams.fSqrtS/0.94)!!! ");
407       return 0;
408     }
409     
410     
411     if(fParams.fCorrS <= 0.) {
412       //see F. Becattini, J. Mannien, M. Gazdzicki, Phys Rev. C73 044905 (2006)
413       fParams.fCorrS = 1. - 0.386* TMath::Exp(-1.23*fParams.fT/fParams.fMuB);
414       
415     }
416     std::cout<<"The phenomenological f-la J. Cleymans et al. PRC73 034905 (2006) for Tch mu_B was used." << std::endl;
417     std::cout<<"The simulation will be done with the calculated parameters:" << std::endl;
418     std::cout<<"Baryon chemical potential = "<<fParams.fMuB<< " [GeV]" << std::endl;
419     std::cout<<"Strangeness chemical potential = "<<fParams.fMuS<< " [GeV]" << std::endl;
420     std::cout<<"Isospin chemical potential = "<<fParams.fMuI3<< " [GeV]" << std::endl;
421     std::cout<<"Strangeness suppression parameter = "<<fParams.fCorrS << std::endl;
422     std::cout<<"Eta_max = "<<fParams.fYlmax<<  std::endl;
423     std::cout << std::endl;
424     
425   }
426   
427   
428   //initialisation of high-pt part 
429   
430   HYJPAR.nhsel = fParams.fNhsel;
431   HYJPAR.ptmin = fParams.fPtmin;
432   HYJPAR.ishad = fParams.fIshad;
433   HYIPAR.bminh = fParams.fBmin;
434   HYIPAR.bmaxh = fParams.fBmax;
435   HYIPAR.AW = fParams.fAw;
436   
437   HYPYIN.ifb = fParams.fIfb;
438   HYPYIN.bfix = fParams.fBfix;
439   HYPYIN.ene = fParams.fSqrtS;
440   
441   PYQPAR.T0 = fParams.fT0;
442   PYQPAR.tau0 = fParams.fTau0;
443   PYQPAR.nf = fParams.fNf;
444   PYQPAR.ienglu = fParams.fIenglu;
445   PYQPAR.ianglu = fParams.fIanglu;
446   
447   
448   myini_();  //
449   
450   
451   // calculation of  multiplicities of different particle species
452   // according to the grand canonical approach
453   GrandCanonical gc(15, fParams.fT, fParams.fMuB, fParams.fMuS, fParams.fMuI3);
454   GrandCanonical gcCh(15, fParams.fT, fParams.fMuB, fParams.fMuS, fParams.fMuI3);
455   GrandCanonical gcPiTh(15, fParams.fThFO, 0., 0., fParams.fMu_th_pip);
456   GrandCanonical gcTh0(15, fParams.fThFO, 0., 0., 0.);
457   
458   //effective volume for central     
459   double dYl= 2 * fParams.fYlmax; //uniform distr. [-Ylmax; Ylmax]  
460   if (fParams.fEtaType >0) dYl = TMath::Sqrt(2 * TMath::Pi()) * fParams.fYlmax ;  //Gaussian distr.                                                                            
461   fVolEff = 2 * TMath::Pi() * fParams.fTau * dYl * (fParams.fR * fParams.fR)/TMath::Power((fParams.fUmax),2) * 
462     ((fParams.fUmax)*TMath::SinH((fParams.fUmax))-TMath::CosH((fParams.fUmax))+ 1);
463   
464   Double_t particleDensityPiCh=0;
465   Double_t particleDensityPiTh=0;
466   
467   if(fParams.fThFO != fParams.fT && fParams.fThFO > 0){
468     particleDensityPiCh = gcCh.ParticleNumberDensity(fDatabase->GetPDGParticle(211));
469     particleDensityPiTh = gcPiTh.ParticleNumberDensity(fDatabase->GetPDGParticle(211));
470   }
471
472   if(fDatabase->GetNParticles()>=kNPartTypes) {
473     cout << "InitialStateHydjet::MultIni(): ERROR Particle definitions in the PDG database exceeds the hardcoded limit of " << kNPartTypes << endl;
474     cout << "      There is either an error with reading the particles file or you might need to increase the maximum allowed definitions" << endl;
475   }
476   for(Int_t particleIndex = 0; particleIndex<fDatabase->GetNParticles() && particleIndex<kNPartTypes; particleIndex++) {
477     ParticlePDG *currParticle = fDatabase->GetPDGParticleByIndex(particleIndex);
478     Int_t encoding = currParticle->GetPDG();
479     //strangeness supression
480     Double_t gammaS = 1;
481     Int_t s = Int_t(currParticle->GetStrangeness());
482     if(encoding == 333) 
483       s = 2;
484     if(fParams.fCorrS < 1. && s != 0)
485       gammaS = TMath::Power(fParams.fCorrS,-TMath::Abs(s));
486     //average densities      
487     Double_t particleDensity = gc.ParticleNumberDensity(currParticle)/gammaS;
488     
489     //compute chemical potential for single f.o. mu==mu_ch
490     Double_t mu = fParams.fMuB  * Int_t(currParticle->GetBaryonNumber()) + 
491       fParams.fMuS  * Int_t(currParticle->GetStrangeness()) +
492       fParams.fMuI3 * Int_t(currParticle->GetElectricCharge());
493     
494     //thermal f.o.
495     if(fParams.fThFO != fParams.fT && fParams.fThFO > 0){
496       Double_t particleDensityCh = gcCh.ParticleNumberDensity(currParticle);
497       Double_t particleDensityTh0 = gcTh0.ParticleNumberDensity(currParticle);
498       Double_t numbDensBolt = particleDensityPiTh*particleDensityCh/particleDensityPiCh;               
499       mu = fParams.fThFO*TMath::Log(numbDensBolt/particleDensityTh0);
500       if(abs(encoding)==211 || encoding==111)mu= fParams.fMu_th_pip; 
501       particleDensity = numbDensBolt;         
502     }
503     
504     // set particle number density to zero for some species
505     // photons
506     if(abs(encoding)==22)
507       particleDensity=0;
508     // K0L and K0S
509     if(abs(encoding)==130 || abs(encoding)==310) {
510       particleDensity=0;
511     }    
512     
513     if(particleDensity > 0.) {
514       if(fParams.fNPartTypes>=kNPartTypes) {
515         Error("in Bool_t MultIni:", "fNPartTypes exceeds the maximum allowed particle species of %d. Check it out!", kNPartTypes);
516         return kFALSE;
517       }
518       fParams.fPartEnc[fParams.fNPartTypes] = encoding;
519       fParams.fPartMult[2 * fParams.fNPartTypes] = particleDensity;
520       fParams.fPartMu[2 * fParams.fNPartTypes] = mu;
521       ++fParams.fNPartTypes;
522     }
523   }
524   return kTRUE;
525 }
526
527 //_________________________________________________________________________________
528 Double_t InitialStateHydjet::SimpsonIntegrator2(Double_t a, Double_t b) {
529   // Simpson integration
530   Int_t nsubIntervals=10000;
531   Double_t h = (b - a)/nsubIntervals; //0-pi, phi
532   Double_t s=0;
533   Double_t x = 0; //phi
534   for(Int_t j = 1; j < nsubIntervals; j++) {
535     x += h; // phi
536     Double_t e = fParams.fEpsilon;
537     Double_t rSB = fParams.fR; //test: podstavit' *coefff_RB
538     Double_t rB = rSB *(TMath::Sqrt(1-e*e)/TMath::Sqrt(1+e*TMath::Cos(2*x))); //f-la7 rB    
539     Double_t sr = SimpsonIntegrator(0,rB,x);
540     s += sr;
541   }
542   return s*h;
543   
544 }
545
546 //_________________________________________________________________________________
547 Double_t InitialStateHydjet::SimpsonIntegrator(Double_t a, Double_t b, Double_t phi) {
548   // Simpson integration
549   Int_t nsubIntervals=100;
550   Double_t h = (b - a)/nsubIntervals;
551   Double_t s = F2(phi,a + 0.5*h);
552   Double_t t = 0.5*(F2(phi,a) + F2(phi,b));
553   Double_t x = a;
554   Double_t y = a + 0.5*h;
555   for(Int_t i = 1; i < nsubIntervals; i++) {
556     x += h;
557     y += h;
558     s += F2(phi,y);
559     t += F2(phi,x);
560   }     
561   t += 2.0*s;
562   return t*h/3.0;
563 }
564
565
566 //f2=f(phi,r)
567 //_________________________________________________________________________________
568 Double_t InitialStateHydjet::F2(Double_t x, Double_t y) {
569   // formula
570   Double_t rSB = fParams.fR; //test: podstavit' *coefff_RB
571   Double_t rhou =  fParams.fUmax * y / rSB;
572   Double_t ff = y*TMath::CosH(rhou)*
573     TMath::Sqrt(1+fParams.fDelta*TMath::Cos(2*x)*TMath::TanH(rhou)*TMath::TanH(rhou));
574 //n_mu u^mu f-la 20
575   return ff;
576 }
577
578 //_________________________________________________________________________________
579 Double_t InitialStateHydjet::MidpointIntegrator2(Double_t a, Double_t b) {
580   // Perform integration through the mid-point method
581   Int_t nsubIntervals=2000; 
582   Int_t nsubIntervals2=1; 
583   Double_t h = (b - a)/nsubIntervals; //0-pi , phi
584   Double_t h2 = (fParams.fR)/nsubIntervals; //0-R maximal rB ?
585   Double_t x = a + 0.5*h;
586   Double_t y = 0;
587   Double_t t = F2(x,y);                    
588   Double_t e = fParams.fEpsilon;
589   for(Int_t j = 1; j < nsubIntervals; j++) {
590     x += h; // integr  phi
591     Double_t rSB = fParams.fR; //test: podstavit' *coefff_RB
592     Double_t  rB = rSB *(TMath::Sqrt(1-e*e)/TMath::Sqrt(1+e*TMath::Cos(2*x))); //f-la7 rB
593     nsubIntervals2 = Int_t(rB / h2)+1;
594     // integr R 
595     y=0;
596     for(Int_t i = 1; i < nsubIntervals2; i++) 
597       t += F2(x,(y += h2));
598   }
599   return t*h*h2;
600 }
601