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 #include "GrandCanonical.h"
  13 #include "HankelFunction.h"
  14 #include "UKUtility.h"
  15 #include "ParticlePDG.h"
  16 #include "DatabasePDG.h"
  17
  18
  19 //_______________________________________________________________________________
  20 GrandCanonical::GrandCanonical():
  21   fTemperature(-1111),
  22   fBaryonPotential(-1111),
  23   fStrangePotential(-1111),
  24   fElectroPotential(-1111),
  25   fNMax(-1111),
  26   fInitialized(kFALSE)
  27 {
  28   //
  29   // default constructor
  30   //
  31 }
  32
  33 //_______________________________________________________________________________
  34 GrandCanonical::GrandCanonical(Int_t nmax, Double_t temperature, Double_t baryonPotential, Double_t strangePotential, Double_t electroPotential):
  35   fTemperature(temperature),
  36   fBaryonPotential(baryonPotential),
  37   fStrangePotential(strangePotential),
  38   fElectroPotential(electroPotential),
  39   fNMax(nmax),
  40   fInitialized(kTRUE)
  41 {
  42   //
  43   // constructor
  44   //
  45 }
  46
  47 //_______________________________________________________________________________
  48 GrandCanonical::~GrandCanonical() {
  49 //
  50 // destructor
  51 //
  52 }
  53
  54 //_______________________________________________________________________________
  55 void GrandCanonical::Temperature(Double_t value) {
  56   //
  57   // set temperature
  58   //
  59   fTemperature = value;
  60   if(fNMax!=-1111 && TMath::Abs(fBaryonPotential+1111)>1.e-10 &&
  61     TMath::Abs(fStrangePotential+1111)>1.e-10 && TMath::Abs(fElectroPotential+1111)>1.e-10)
  62     fInitialized = kTRUE;
  63 }
  64
  65 //_______________________________________________________________________________
  66 void GrandCanonical::BaryonPotential(Double_t value) {
  67   //
  68   // set baryo chemical potential
  69   //
  70   fBaryonPotential = value;
  71   if(fNMax!=-1111 && TMath::Abs(fTemperature+1111)>1.e-10 &&
  72     TMath::Abs(fStrangePotential+1111)>1.e-10 && TMath::Abs(fElectroPotential+1111)>1.e-10)
  73     fInitialized = kTRUE;
  74 }
  75
  76 //_______________________________________________________________________________
  77 void GrandCanonical::StrangePotential(Double_t value) {
  78   //
  79   // set strange potential
  80   //
  81   fStrangePotential = value;
  82   if(fNMax!=-1111 && TMath::Abs(fTemperature+1111)>1.e-10 &&
  83      TMath::Abs(fBaryonPotential+1111)>1.e-10 && TMath::Abs(fElectroPotential+1111)>1.e-10)
  84     fInitialized = kTRUE;
  85 }
  86
  87 //_______________________________________________________________________________
  88 void GrandCanonical::ElectroPotential(Double_t value) {
  89   //
  90   // set electro chemical potential
  91   //
  92   fElectroPotential = value;
  93   if(fNMax!=-1111 && TMath::Abs(fTemperature+1111)>1.e-10 &&
  94      TMath::Abs(fBaryonPotential+1111)>1.e-10 && TMath::Abs(fStrangePotential+1111)>1.e-10)
  95     fInitialized = kTRUE;
  96 }
  97
  98 //_______________________________________________________________________________
  99 void GrandCanonical::NMax(Int_t value) {
 100   //
 101   // set the number of iterations
 102   //
 103   fNMax = value;
 104   if(fTemperature!=-1111 && fBaryonPotential!=-1111 && fStrangePotential!=-1111 && fElectroPotential!=-1111)
 105     if(fNMax!=-1111 && TMath::Abs(fBaryonPotential+1111)>1.e-10 &&
 106     TMath::Abs(fStrangePotential+1111)>1.e-10 && TMath::Abs(fElectroPotential+1111)>1.e-10)
 107     fInitialized = kTRUE;
 108 }
 109
 110 //_______________________________________________________________________________
 111 Double_t GrandCanonical::ParticleEnergyDensity(ParticlePDG *const particle) {
 112   //
 113   // compute the energy density for a given particle
 114   //
 115   // Check if all the thermodinamic parameters are set
 116   if(!fInitialized)
 117     Fatal("GrandCanonical::ParticleEnergyDensity", "GrandCanonical object not fully initialized!!");
 118
 119   // Compute the particle energy density
 120   Double_t degFactor = 2.*particle->GetSpin() + 1.;                                    // degeneracy factor
 121   Double_t mass = particle->GetMass();                                                // PDG table mass
 122   Double_t d = Int_t(2.*particle->GetSpin()) & 1 ? 1. : -1;                                   // Bose-Einstein/Fermi-Dirac factor
 123   Double_t preFactor = (degFactor*mass*mass*fTemperature*fTemperature/kHbarc/kHbarc/kHbarc)/(2.*TMath::Pi()*TMath::Pi());
 124
 125   Double_t postFactor = 0.;
 126   //compute chemical potential
 127   Double_t potential = fBaryonPotential * particle->GetBaryonNumber() +
 128                        fStrangePotential * particle->GetStrangeness() +
 129                        fElectroPotential * particle->GetElectricCharge();
 130
 131   for(Int_t n = 1; n <= fNMax; ++n) {
 132     postFactor += TMath::Power(-d, n+1)/(n*n) *
 133                   TMath::Exp(n*potential/fTemperature) *
 134                   (3.*HankelKn(2, n*mass/fTemperature) + n*mass/fTemperature*HankelK1(n*mass/fTemperature));
 135   }
 136   return preFactor * postFactor;
 137 }
 138
 139 //_______________________________________________________________________________
 140 Double_t GrandCanonical::ParticleNumberDensity(ParticlePDG *const particle) {
 141   //
 142   // compute the particle number density
 143   //
 144   // Check if all the thermodinamic parameters are set
 145   if(!fInitialized)
 146     Fatal("GrandCanonical::ParticleNumberDensity", "GrandCanonical object not fully initialized!!");
 147
 148   Double_t degFactor = 2.*particle->GetSpin() + 1.;
 149   Double_t mass = particle->GetMass();
 150   Double_t d = Int_t(2*particle->GetSpin()) & 1 ? 1. : -1.;
 151   Double_t preFactor = (degFactor*mass*mass*fTemperature/kHbarc/kHbarc/kHbarc)/(2.*TMath::Pi()*TMath::Pi());
 152
 153   Double_t postFactor = 0.;
 154   Double_t potential = fBaryonPotential * particle->GetBaryonNumber() +
 155                        fStrangePotential * particle->GetStrangeness() +
 156                        fElectroPotential * particle->GetElectricCharge();
 157
 158   for(Int_t n = 1; n <= fNMax; ++n) {
 159     postFactor += TMath::Power(-d, n+1)/n *
 160                   TMath::Exp(n*potential/fTemperature) *
 161                   HankelKn(2, n*mass/fTemperature);
 162   }
 163   return preFactor * postFactor;
 164 }
 165
 166 //_______________________________________________________________________________
 167 Double_t GrandCanonical::EnergyDensity(DatabasePDG *const database) {
 168   //
 169   // compute the total energy density
 170   //
 171   // Check if all the thermodinamic parameters are set
 172   if(!fInitialized)
 173     Fatal("GrandCanonical::EnergyDensity", "GrandCanonical object not fully initialized!!");
 174
 175   Double_t meanEnergyDensity = 0.;
 176
 177   for(Int_t currParticle = 0; currParticle<database->GetNParticles(); currParticle++) {
 178     ParticlePDG *particle = database->GetPDGParticleByIndex(currParticle);
 179     meanEnergyDensity += ParticleEnergyDensity(particle);
 180   }
 181
 182   return meanEnergyDensity;
 183 }
 184
 185 //_______________________________________________________________________________
 186 Double_t GrandCanonical::BaryonDensity(DatabasePDG *const database) {
 187   //
 188   // compute the baryon density
 189   //
 190   // Check if all the thermodinamic parameters are set
 191   if(!fInitialized)
 192     Fatal("GrandCanonical::BaryonDensity", "GrandCanonical object not fully initialized!!");
 193
 194   Double_t meanBaryonDensity = 0.;
 195
 196   for(Int_t currParticle = 0; currParticle<database->GetNParticles(); currParticle++) {
 197     ParticlePDG *particle = database->GetPDGParticleByIndex(currParticle);
 198     meanBaryonDensity += ParticleNumberDensity(particle)*particle->GetBaryonNumber();
 199   }
 200   return meanBaryonDensity;
 201 }
 202
 203 //_______________________________________________________________________________
 204 Double_t GrandCanonical::StrangeDensity(DatabasePDG *const database) {
 205   //
 206   // compute the strangeness density
 207   //
 208   // Check if all the thermodinamic parameters are set
 209   if(!fInitialized)
 210     Fatal("GrandCanonical::StrangeDensity", "GrandCanonical object not fully initialized!!");
 211
 212   Double_t meanStrangeDensity = 0.;
 213
 214   for(Int_t currParticle = 0; currParticle<database->GetNParticles(); currParticle++) {
 215     ParticlePDG *particle = database->GetPDGParticleByIndex(currParticle);
 216     meanStrangeDensity += ParticleNumberDensity(particle)*particle->GetStrangeness();
 217   }
 218
 219   return meanStrangeDensity;
 220 }
 221
 222 //_______________________________________________________________________________
 223 Double_t GrandCanonical::ElectroDensity(DatabasePDG *const database) {
 224   //
 225   // compute the electro number density
 226   //
 227   // Check if all the thermodinamic parameters are set
 228   if(!fInitialized)
 229     Fatal("GrandCanonical::ElectroDensity", "GrandCanonical object not fully initialized!!");
 230
 231   Double_t meanElectroDensity = 0.;
 232
 233   //hadrons
 234   for(Int_t currParticle = 0; currParticle<database->GetNParticles(); currParticle++) {
 235     ParticlePDG *particle = database->GetPDGParticleByIndex(currParticle);
 236     meanElectroDensity += ParticleNumberDensity(particle)*particle->GetElectricCharge();
 237   }
 238
 239   return meanElectroDensity;
 240 }