98ab726f33a920449ef9f3e14c9fe1ede754df55
[u/mrichter/AliRoot.git] / TUHKMgen / UHKM / GrandCanonical.cxx
1 /*                                                                            
2                                                                             
3         Nikolai Amelin, Ludmila Malinina, Timur Pocheptsov (C) JINR/Dubna
4       amelin@sunhe.jinr.ru, malinina@sunhe.jinr.ru, pocheptsov@sunhe.jinr.ru 
5                            November. 2, 2005                                
6
7 */
8 #include <iostream>
9
10 #include <TError.h>
11 #include <TMath.h>
12 #ifndef GRANDCANONICAL_INCLUDED
13 #include "GrandCanonical.h"
14 #endif
15 #ifndef HANKELFUNCTION_INCLUDED
16 #include "HankelFunction.h"
17 #endif
18 #ifndef UKUTILITY_INCLUDED
19 #include "UKUtility.h"
20 #endif
21
22 GrandCanonical::GrandCanonical():
23   fTemperature(-1111),
24   fBaryonPotential(-1111),
25   fStrangePotential(-1111),
26   fElectroPotential(-1111),
27   fNMax(-1111),
28   fInitialized(kFALSE)
29 {
30 }
31
32 GrandCanonical::GrandCanonical(Int_t nmax, Double_t temperature, Double_t baryonPotential, Double_t strangePotential, Double_t electroPotential):
33   fTemperature(temperature),
34   fBaryonPotential(baryonPotential),
35   fStrangePotential(strangePotential),
36   fElectroPotential(electroPotential),
37   fNMax(nmax),
38   fInitialized(kTRUE)
39 {
40 }
41
42 GrandCanonical::~GrandCanonical() {}
43
44
45 void GrandCanonical::Temperature(Double_t value) {
46   fTemperature = value;
47   if(fNMax!=-1111 && fBaryonPotential!=-1111 && fStrangePotential!=-1111 && fElectroPotential!=-1111)
48     fInitialized = kTRUE;
49 }
50
51 void GrandCanonical::BaryonPotential(Double_t value) {
52   fBaryonPotential = value;
53   if(fNMax!=-1111 && fTemperature!=-1111 && fStrangePotential!=-1111 && fElectroPotential!=-1111)
54     fInitialized = kTRUE;
55 }
56
57 void GrandCanonical::StrangePotential(Double_t value) {
58   fStrangePotential = value;
59   if(fNMax!=-1111 && fTemperature!=-1111 && fBaryonPotential!=-1111 && fElectroPotential!=-1111)
60     fInitialized = kTRUE;
61 }
62
63 void GrandCanonical::ElectroPotential(Double_t value) {
64   fElectroPotential = value;
65   if(fNMax!=-1111 && fTemperature!=-1111 && fBaryonPotential!=-1111 && fStrangePotential!=-1111)
66     fInitialized = kTRUE;
67 }
68
69 void GrandCanonical::NMax(Int_t value) {
70   fNMax = value;
71   if(fTemperature!=-1111 && fBaryonPotential!=-1111 && fStrangePotential!=-1111 && fElectroPotential!=-1111)
72     fInitialized = kTRUE;
73 }
74
75 Double_t GrandCanonical::ParticleEnergyDensity(ParticlePDG* particle) {
76   // Check if all the thermodinamic parameters are set
77   if(!fInitialized)
78     Fatal("GrandCanonical::ParticleEnergyDensity", "GrandCanonical object not fully initialized!!");
79   
80   // Compute the particle energy density
81   Double_t degFactor = 2.*particle->GetSpin() + 1.;                                    // degeneracy factor
82   Double_t mass = particle->GetMass();                                                // PDG table mass
83   Double_t d = Int_t(2.*particle->GetSpin()) & 1 ? 1. : -1;                                   // Bose-Einstein/Fermi-Dirac factor
84   Double_t preFactor = (degFactor*mass*mass*fTemperature*fTemperature/hbarc/hbarc/hbarc)/(2.*TMath::Pi()*TMath::Pi()); 
85
86   Double_t postFactor = 0.;
87   //compute chemical potential
88   Double_t potential = fBaryonPotential * particle->GetBaryonNumber() + 
89                        fStrangePotential * particle->GetStrangeness() +
90                        fElectroPotential * particle->GetElectricCharge();
91
92   for(Int_t n = 1; n <= fNMax; ++n) {
93     postFactor += TMath::Power(-d, n+1)/(n*n) *
94                   TMath::Exp(n*potential/fTemperature) *
95                   (3.*HankelKn(2, n*mass/fTemperature) + n*mass/fTemperature*HankelK1(n*mass/fTemperature));
96   }
97   return preFactor * postFactor;
98 }
99
100 Double_t GrandCanonical::ParticleNumberDensity(ParticlePDG* particle) {
101   // Check if all the thermodinamic parameters are set
102   if(!fInitialized)
103     Fatal("GrandCanonical::ParticleNumberDensity", "GrandCanonical object not fully initialized!!");
104
105   Double_t degFactor = 2.*particle->GetSpin() + 1.;
106   Double_t mass = particle->GetMass();     
107   Double_t d = Int_t(2*particle->GetSpin()) & 1 ? 1. : -1.;
108   Double_t preFactor = (degFactor*mass*mass*fTemperature/hbarc/hbarc/hbarc)/(2.*TMath::Pi()*TMath::Pi());
109
110   Double_t postFactor = 0.;
111   Double_t potential = fBaryonPotential * particle->GetBaryonNumber() + 
112                        fStrangePotential * particle->GetStrangeness() +
113                        fElectroPotential * particle->GetElectricCharge();
114                            
115   for(Int_t n = 1; n <= fNMax; ++n) {
116     postFactor += TMath::Power(-d, n+1)/n * 
117                   TMath::Exp(n*potential/fTemperature) *
118                   HankelKn(2, n*mass/fTemperature);        
119   }
120   return preFactor * postFactor;
121 }
122
123
124 Double_t GrandCanonical::EnergyDensity(DatabasePDG* database) {
125   // Check if all the thermodinamic parameters are set
126   if(!fInitialized)
127     Fatal("GrandCanonical::EnergyDensity", "GrandCanonical object not fully initialized!!");
128
129   Double_t meanEnergyDensity = 0.;
130
131   for(Int_t currParticle = 0; currParticle<database->GetNParticles(); currParticle++) {
132     ParticlePDG *particle = database->GetPDGParticleByIndex(currParticle);
133     meanEnergyDensity += ParticleEnergyDensity(particle);
134   }
135
136   return meanEnergyDensity;
137 }
138
139 Double_t GrandCanonical::BaryonDensity(DatabasePDG* database) {
140   // Check if all the thermodinamic parameters are set
141   if(!fInitialized)
142     Fatal("GrandCanonical::BaryonDensity", "GrandCanonical object not fully initialized!!");
143
144   Double_t meanBaryonDensity = 0.;
145
146   for(Int_t currParticle = 0; currParticle<database->GetNParticles(); currParticle++) {
147     ParticlePDG *particle = database->GetPDGParticleByIndex(currParticle);
148     meanBaryonDensity += ParticleNumberDensity(particle)*particle->GetBaryonNumber();
149   }
150   return meanBaryonDensity;
151 }
152
153 Double_t GrandCanonical::StrangeDensity(DatabasePDG* database) {
154   // Check if all the thermodinamic parameters are set
155   if(!fInitialized)
156     Fatal("GrandCanonical::StrangeDensity", "GrandCanonical object not fully initialized!!");
157
158   Double_t meanStrangeDensity = 0.;
159
160   for(Int_t currParticle = 0; currParticle<database->GetNParticles(); currParticle++) {
161     ParticlePDG *particle = database->GetPDGParticleByIndex(currParticle);
162     meanStrangeDensity += ParticleNumberDensity(particle)*particle->GetStrangeness();
163   }
164
165   return meanStrangeDensity;
166 }
167
168 Double_t GrandCanonical::ElectroDensity(DatabasePDG* database) {
169   // Check if all the thermodinamic parameters are set
170   if(!fInitialized)
171     Fatal("GrandCanonical::ElectroDensity", "GrandCanonical object not fully initialized!!");
172
173   Double_t meanElectroDensity = 0.;
174   
175   //hadrons
176   for(Int_t currParticle = 0; currParticle<database->GetNParticles(); currParticle++) {
177     ParticlePDG *particle = database->GetPDGParticleByIndex(currParticle);
178     meanElectroDensity += ParticleNumberDensity(particle)*particle->GetElectricCharge();
179   }
180
181   return meanElectroDensity;
182 }