]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TUHKMgen/UHKM/InitialStateHydjet.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fix for Coverity defect 14006: RESOURCE_LEAK
[u/mrichter/AliRoot.git] / TUHKMgen / UHKM / InitialStateHydjet.cxx
CommitLineData
3fa37a65 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
03896fc4 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
3fa37a65 23
03896fc4 24
b1c2e580 25
b1c2e580 26
03896fc4 27#include <iostream>
28#include <fstream>
b1c2e580 29
b1c2e580 30#include <TLorentzVector.h>
31#include <TVector3.h>
b1c2e580 32#include <TMath.h>
33
b1c2e580 34#include "InitialStateHydjet.h"
b1c2e580 35#include "RandArrayFunction.h"
b1c2e580 36#include "GrandCanonical.h"
b1c2e580 37#include "StrangePotential.h"
b1c2e580 38#include "Particle.h"
b1c2e580 39#include "ParticlePDG.h"
b1c2e580 40#include "HYJET_COMMONS.h"
03896fc4 41
b1c2e580 42extern "C" void hyevnt_();
43extern "C" void myini_();
b1c2e580 44extern HYIPARCommon HYIPAR;
45extern HYFPARCommon HYFPAR;
46extern HYJPARCommon HYJPAR;
47extern HYPARTCommon HYPART;
48extern SERVICECommon SERVICE;
49
50using std::cout;
51using std::endl;
52
3fa37a65 53class ParticleAllocator;
54class TRandom3;
55
56// declaration of the static member fLastIndex
03896fc4 57Int_t Particle::fgLastIndex;
3fa37a65 58
03896fc4 59//_________________________________________________________________________________
b1c2e580 60void InitialStateHydjet::Initialize(List_t &source, ParticleAllocator & allocator) {
7b7936e9 61 // Generate initial particles from the soft and hard components
3fa37a65 62
63 // Initialize the static "last index variable"
64 Particle::InitIndexing();
65
b1c2e580 66 //----- high-pt part------------------------------
b1c2e580 67 TLorentzVector partJMom, partJPos, zeroVec;
b1c2e580 68
7b7936e9 69 // run a HYDJET event
70 hyevnt_();
77b9c3bb 71
72 fBgen = HYFPAR.bgen * HYIPAR.RA;
73 fNpart = HYFPAR.npart;
74 fNcoll = HYFPAR.nbcol;
75
76
77
b1c2e580 78 if(fParams.fNhsel != 0) {
3fa37a65 79 //get number of particles in jets
b1c2e580 80 Int_t numbJetPart = HYPART.njp;
3fa37a65 81
82 for(Int_t i = 0; i <numbJetPart; i++) {
83 Int_t pdg = Int_t(HYPART.ppart[i][1]);
b1c2e580 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];
b1c2e580 92 ParticlePDG *partDef = fDatabase->GetPDGParticle(pdg);
93 Int_t type =1; //from jet
3fa37a65 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 }
b1c2e580 102 }
7b7936e9 103 } //nhsel !=0 not only hydro!
b1c2e580 104
b1c2e580 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;
b1c2e580 116 //set maximal hadron energy
117 const Double_t eMax = 5.;
118 //-------------------------------------
119 // get impact parameter
b1c2e580 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.
3eeebf49 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);
b1c2e580 127
128 //effective volume for non-central Simpson2
7b7936e9 129 Double_t volEffnoncent = fParams.fTau * dYl * SimpsonIntegrator2(0., 2.*TMath::Pi());
130 fVolEff = volEffcent * HYFPAR.npart/HYFPAR.npart0;
b1c2e580 131
7b7936e9 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);
b1c2e580 135
7b7936e9 136 double veff=fVolEff;
77b9c3bb 137
b1c2e580 138 //------------------------------------
139 //cycle on particles types
140 for(Int_t i = 0; i < fParams.fNPartTypes; ++i) {
7b7936e9 141 Double_t mparam = fParams.fPartMult[2 * i] * veff;
142 Int_t multiplicity = gRandom->Poisson(mparam);
b1c2e580 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 !
03896fc4 152 if(TMath::Abs(partDef->GetCharmQNumber())>0 || TMath::Abs(partDef->GetCharmAQNumber())>0){
b1c2e580 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;
7b7936e9 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) /
b1c2e580 170 (TMath::Exp((x - mu) / (fParams.fThFO)) + d);
7b7936e9 171 if(x>=mu && fParams.fThFO<=0)probList[ii] = x * TMath::Sqrt(x * x - mass * mass) /
b1c2e580 172 (TMath::Exp((x - mu) / (fParams.fT)) + d);
7b7936e9 173 if(x<mu)probList[ii] = 0.;
b1c2e580 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);
7b7936e9 182 for(ii = 0; ii < nBins; ++ii) {
183 probList[ii] = x * TMath::CosH(param*x);
b1c2e580 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.;
7b7936e9 191 Double_t r, rB, phiF;
b1c2e580 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
7b7936e9 200 rB = fParams.fR * coeffRB * coeffR1;
b1c2e580 201
202 Double_t rho = TMath::Sqrt(gRandom->Rndm());
203 Double_t phi = TMath::TwoPi() * gRandom->Rndm();
7b7936e9 204 Double_t rx = TMath::Sqrt(1-fParams.fEpsilon)*rB;
205 Double_t ry = TMath::Sqrt(1+fParams.fEpsilon)*rB;
b1c2e580 206
7b7936e9 207 x0 = rx * rho * TMath::Cos(phi);
208 y0 = ry * rho * TMath::Sin(phi);
b1c2e580 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
7b7936e9 217 Double_t tau = coeffR1 * fParams.fTau + sqrt(2.) * fParams.fSigmaTau * coeffR1 * (gRandom->Gaus());
b1c2e580 218 z0 = tau * TMath::SinH(etaF);
219 t0 = tau * TMath::CosH(etaF);
220
7b7936e9 221 Double_t rhou = fParams.fUmax * r / rB;
b1c2e580 222
b1c2e580 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.;
bfd20868 238 Double_t stp0 = TMath::Sqrt((1.-ctp0)*(1.+ctp0));
b1c2e580 239 e = mass + (eMax - mass) * arrayFunctDistE();
bfd20868 240 Double_t pp0 = TMath::Sqrt((e-mass)*(e+mass));
b1c2e580 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
b1c2e580 250 partMom.SetXYZT(px0, py0, pz0, e);
251 partPos.SetXYZT(x0, y0, z0, t0);
252 partMom.Boost(vec3);
253 Int_t type =0; //hydro
7b7936e9 254
3fa37a65 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;
b1c2e580 259 } //nhsel==4 , no hydro part
260 }
261 }
262 }
263
b1c2e580 264}
265
03896fc4 266//_________________________________________________________________________________
b1c2e580 267Bool_t InitialStateHydjet::ReadParams() {
7b7936e9 268 // Read parameters from an input file in ascii
269
b1c2e580 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
03896fc4 377//_________________________________________________________________________________
b1c2e580 378Bool_t InitialStateHydjet::MultIni() {
7b7936e9 379 // Calculate average multiplicities, chemical potentials (if necessary),
380 // initialize pyquen
381
b1c2e580 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
712797b3 397 StrangePotential psp(0., fDatabase);
398 psp.SetBaryonPotential(fParams.fMuB);
399 psp.SetTemperature(fParams.fT);
b1c2e580 400 //compute strangeness potential
401 if(fParams.fMuB > 0.01)
712797b3 402 fParams.fMuS = psp.CalculateStrangePotential();
b1c2e580 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);
b1c2e580 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
b1c2e580 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
7b7936e9 448 myini_(); //
b1c2e580 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);
3fa37a65 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.);
b1c2e580 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);
b1c2e580 463
3fa37a65 464 Double_t particleDensityPiCh=0;
465 Double_t particleDensityPiTh=0;
b1c2e580 466
467 if(fParams.fThFO != fParams.fT && fParams.fThFO > 0){
3fa37a65 468 particleDensityPiCh = gcCh.ParticleNumberDensity(fDatabase->GetPDGParticle(211));
469 particleDensityPiTh = gcPiTh.ParticleNumberDensity(fDatabase->GetPDGParticle(211));
b1c2e580 470 }
471
472 for(Int_t particleIndex = 0; particleIndex < fDatabase->GetNParticles(); particleIndex++) {
473 ParticlePDG *currParticle = fDatabase->GetPDGParticleByIndex(particleIndex);
474 Int_t encoding = currParticle->GetPDG();
475 //strangeness supression
476 Double_t gammaS = 1;
7b7936e9 477 Int_t s = Int_t(currParticle->GetStrangeness());
478 if(encoding == 333)
479 s = 2;
480 if(fParams.fCorrS < 1. && s != 0)
481 gammaS = TMath::Power(fParams.fCorrS,-TMath::Abs(s));
b1c2e580 482 //average densities
7b7936e9 483 Double_t particleDensity = gc.ParticleNumberDensity(currParticle)/gammaS;
b1c2e580 484
485 //compute chemical potential for single f.o. mu==mu_ch
486 Double_t mu = fParams.fMuB * Int_t(currParticle->GetBaryonNumber()) +
487 fParams.fMuS * Int_t(currParticle->GetStrangeness()) +
488 fParams.fMuI3 * Int_t(currParticle->GetElectricCharge());
489
490 //thermal f.o.
491 if(fParams.fThFO != fParams.fT && fParams.fThFO > 0){
3fa37a65 492 Double_t particleDensityCh = gcCh.ParticleNumberDensity(currParticle);
493 Double_t particleDensityTh0 = gcTh0.ParticleNumberDensity(currParticle);
494 Double_t numbDensBolt = particleDensityPiTh*particleDensityCh/particleDensityPiCh;
495 mu = fParams.fThFO*TMath::Log(numbDensBolt/particleDensityTh0);
b1c2e580 496 if(abs(encoding)==211 || encoding==111)mu= fParams.fMu_th_pip;
3fa37a65 497 particleDensity = numbDensBolt;
b1c2e580 498 }
499
7b7936e9 500 // set particle number density to zero for some species
501 // photons
502 if(abs(encoding)==22)
503 particleDensity=0;
504 // K0L and K0S
505 if(abs(encoding)==130 || abs(encoding)==310) {
506 particleDensity=0;
507 }
b1c2e580 508
509 if(particleDensity > 0.) {
510 // outMult<<encoding<< " " <<particleDensity<< " "<<mu<<std::endl;
511 fParams.fPartEnc[fParams.fNPartTypes] = encoding;
512 fParams.fPartMult[2 * fParams.fNPartTypes] = particleDensity;
513 fParams.fPartMu[2 * fParams.fNPartTypes] = mu;
514 ++fParams.fNPartTypes;
515 if(fParams.fNPartTypes > 1000)
516 Error("in Bool_t MultIni:", "fNPartTypes is too large %d", fParams.fNPartTypes);
517 }
518 }
519 return kTRUE;
520}
521
03896fc4 522//_________________________________________________________________________________
b1c2e580 523Double_t InitialStateHydjet::SimpsonIntegrator2(Double_t a, Double_t b) {
7b7936e9 524 // Simpson integration
b1c2e580 525 Int_t nsubIntervals=10000;
526 Double_t h = (b - a)/nsubIntervals; //0-pi, phi
527 Double_t s=0;
b1c2e580 528 Double_t x = 0; //phi
529 for(Int_t j = 1; j < nsubIntervals; j++) {
530 x += h; // phi
531 Double_t e = fParams.fEpsilon;
3fa37a65 532 Double_t rSB = fParams.fR; //test: podstavit' *coefff_RB
533 Double_t rB = rSB *(TMath::Sqrt(1-e*e)/TMath::Sqrt(1+e*TMath::Cos(2*x))); //f-la7 rB
7b7936e9 534 Double_t sr = SimpsonIntegrator(0,rB,x);
b1c2e580 535 s += sr;
536 }
537 return s*h;
538
539}
540
03896fc4 541//_________________________________________________________________________________
b1c2e580 542Double_t InitialStateHydjet::SimpsonIntegrator(Double_t a, Double_t b, Double_t phi) {
7b7936e9 543 // Simpson integration
b1c2e580 544 Int_t nsubIntervals=100;
545 Double_t h = (b - a)/nsubIntervals;
03896fc4 546 Double_t s = F2(phi,a + 0.5*h);
547 Double_t t = 0.5*(F2(phi,a) + F2(phi,b));
b1c2e580 548 Double_t x = a;
549 Double_t y = a + 0.5*h;
550 for(Int_t i = 1; i < nsubIntervals; i++) {
551 x += h;
552 y += h;
03896fc4 553 s += F2(phi,y);
554 t += F2(phi,x);
b1c2e580 555 }
556 t += 2.0*s;
557 return t*h/3.0;
558}
559
560
561//f2=f(phi,r)
03896fc4 562//_________________________________________________________________________________
563Double_t InitialStateHydjet::F2(Double_t x, Double_t y) {
7b7936e9 564 // formula
3fa37a65 565 Double_t rSB = fParams.fR; //test: podstavit' *coefff_RB
566 Double_t rhou = fParams.fUmax * y / rSB;
b1c2e580 567 Double_t ff = y*TMath::CosH(rhou)*
568 TMath::Sqrt(1+fParams.fDelta*TMath::Cos(2*x)*TMath::TanH(rhou)*TMath::TanH(rhou));
569//n_mu u^mu f-la 20
570 return ff;
571}
572
03896fc4 573//_________________________________________________________________________________
b1c2e580 574Double_t InitialStateHydjet::MidpointIntegrator2(Double_t a, Double_t b) {
7b7936e9 575 // Perform integration through the mid-point method
b1c2e580 576 Int_t nsubIntervals=2000;
577 Int_t nsubIntervals2=1;
578 Double_t h = (b - a)/nsubIntervals; //0-pi , phi
7b7936e9 579 Double_t h2 = (fParams.fR)/nsubIntervals; //0-R maximal rB ?
b1c2e580 580 Double_t x = a + 0.5*h;
581 Double_t y = 0;
03896fc4 582 Double_t t = F2(x,y);
b1c2e580 583 Double_t e = fParams.fEpsilon;
b1c2e580 584 for(Int_t j = 1; j < nsubIntervals; j++) {
585 x += h; // integr phi
3fa37a65 586 Double_t rSB = fParams.fR; //test: podstavit' *coefff_RB
587 Double_t rB = rSB *(TMath::Sqrt(1-e*e)/TMath::Sqrt(1+e*TMath::Cos(2*x))); //f-la7 rB
7b7936e9 588 nsubIntervals2 = Int_t(rB / h2)+1;
b1c2e580 589 // integr R
590 y=0;
591 for(Int_t i = 1; i < nsubIntervals2; i++)
03896fc4 592 t += F2(x,(y += h2));
b1c2e580 593 }
594 return t*h*h2;
595}
596