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