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