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