GrandCanonical.cc

Go to the documentation of this file.

/*                                                                            
                                                                            
Nikolai Amelin, Ludmila Malinina, Timur Pocheptsov (C) JINR/Dubna
amelin@sunhe.jinr.ru, malinina@sunhe.jinr.ru, pocheptsov@sunhe.jinr.ru 
November. 2, 2005                                

*/
#include <iostream>
#include <TMath.h>

#include "GeneratorInterface/Hydjet2Interface/interface/GrandCanonical.h"
#include "GeneratorInterface/Hydjet2Interface/interface/HankelFunction.h"
#include "GeneratorInterface/Hydjet2Interface/interface/UKUtility.h"

GrandCanonical::GrandCanonical() {
  fInitialized = kFALSE;
  fNMax = -1111;
  fTemperature = -1111;
  fBaryonPotential = -1111;
  fStrangePotential = -1111;
  fElectroPotential = -1111;
  fCharmPotential = -1111;
}

GrandCanonical::GrandCanonical(int nmax,
                               double temperature,
                               double baryonPotential,
                               double strangePotential,
                               double electroPotential,
                               double charmPotential) {
  fNMax = nmax;
  fTemperature = temperature;
  fBaryonPotential = baryonPotential;
  fStrangePotential = strangePotential;
  fElectroPotential = electroPotential;
  fCharmPotential = charmPotential;
  fInitialized = kTRUE;
}

GrandCanonical::~GrandCanonical() {}

void GrandCanonical::Temperature(double value) {
  fTemperature = value;
  if (fNMax != -1111 && fBaryonPotential != -1111 && fStrangePotential != -1111 && fElectroPotential != -1111 &&
      fCharmPotential != -1111)
    fInitialized = kTRUE;
}

void GrandCanonical::BaryonPotential(double value) {
  fBaryonPotential = value;
  if (fNMax != -1111 && fTemperature != -1111 && fStrangePotential != -1111 && fElectroPotential != -1111 &&
      fCharmPotential != -1111)
    fInitialized = kTRUE;
}

void GrandCanonical::StrangePotential(double value) {
  fStrangePotential = value;
  if (fNMax != -1111 && fTemperature != -1111 && fBaryonPotential != -1111 && fElectroPotential != -1111 &&
      fCharmPotential != -1111)
    fInitialized = kTRUE;
}

void GrandCanonical::ElectroPotential(double value) {
  fElectroPotential = value;
  if (fNMax != -1111 && fTemperature != -1111 && fBaryonPotential != -1111 && fStrangePotential != -1111 &&
      fCharmPotential != -1111)
    fInitialized = kTRUE;
}

void GrandCanonical::CharmPotential(double value) {
  fCharmPotential = value;
  if (fNMax != -1111 && fTemperature != -1111 && fBaryonPotential != -1111 && fStrangePotential != -1111 &&
      fElectroPotential != -1111)
    fInitialized = kTRUE;
}

void GrandCanonical::NMax(int value) {
  fNMax = value;
  if (fTemperature != -1111 && fBaryonPotential != -1111 && fStrangePotential != -1111 && fElectroPotential != -1111 &&
      fCharmPotential != -1111)
    fInitialized = kTRUE;
}

double GrandCanonical::ParticleEnergyDensity(ParticlePDG* particle) {
  // Check if all the thermodinamic parameters are set
  if (!fInitialized)
    Fatal("GrandCanonical::ParticleEnergyDensity", "GrandCanonical object not fully initialized!!");

  // Compute the particle energy density
  double degFactor = 2. * particle->GetSpin() + 1.;        // degeneracy factor
  double mass = particle->GetMass();                       // PDG table mass
  double d = int(2. * particle->GetSpin()) & 1 ? 1. : -1;  // Bose-Einstein/Fermi-Dirac factor
  double preFactor = (degFactor * mass * mass * fTemperature * fTemperature / hbarc / hbarc / hbarc) /
                     (2. * TMath::Pi() * TMath::Pi());

  double postFactor = 0.;
  //compute chemical potential
  double potential = fBaryonPotential * particle->GetBaryonNumber() + fStrangePotential * particle->GetStrangeness() +
                     fElectroPotential * particle->GetElectricCharge() + fCharmPotential * particle->GetCharmness();

  for (int n = 1; n <= fNMax; ++n) {
    postFactor +=
        TMath::Power(-d, n + 1) / (n * n) * TMath::Exp(n * potential / fTemperature) *
        (3. * HankelKn(2, n * mass / fTemperature) + n * mass / fTemperature * HankelK1(n * mass / fTemperature));
  }
  return preFactor * postFactor;
}

double GrandCanonical::ParticleNumberDensity(ParticlePDG* particle) {
  // Check if all the thermodinamic parameters are set
  if (!fInitialized)
    Fatal("GrandCanonical::ParticleNumberDensity", "GrandCanonical object not fully initialized!!");

  double degFactor = 2. * particle->GetSpin() + 1.;
  double mass = particle->GetMass();
  double d = int(2 * particle->GetSpin()) & 1 ? 1. : -1.;
  double preFactor =
      (degFactor * mass * mass * fTemperature / hbarc / hbarc / hbarc) / (2. * TMath::Pi() * TMath::Pi());

  double postFactor = 0.;
  double potential = fBaryonPotential * particle->GetBaryonNumber() + fStrangePotential * particle->GetStrangeness() +
                     fElectroPotential * particle->GetElectricCharge() + fCharmPotential * particle->GetCharmness();

  for (int n = 1; n <= fNMax; ++n) {
    postFactor +=
        TMath::Power(-d, n + 1) / n * TMath::Exp(n * potential / fTemperature) * HankelKn(2, n * mass / fTemperature);
  }
  return preFactor * postFactor;
}

double GrandCanonical::EnergyDensity(DatabasePDG* database) {
  // Check if all the thermodinamic parameters are set
  if (!fInitialized)
    Fatal("GrandCanonical::EnergyDensity", "GrandCanonical object not fully initialized!!");

  double meanEnergyDensity = 0.;

  for (int currParticle = 0; currParticle < database->GetNParticles(); currParticle++) {
    ParticlePDG* particle = database->GetPDGParticleByIndex(currParticle);
    meanEnergyDensity += ParticleEnergyDensity(particle);
  }

  return meanEnergyDensity;
}

double GrandCanonical::BaryonDensity(DatabasePDG* database) {
  // Check if all the thermodinamic parameters are set
  if (!fInitialized)
    Fatal("GrandCanonical::BaryonDensity", "GrandCanonical object not fully initialized!!");

  double meanBaryonDensity = 0.;

  for (int currParticle = 0; currParticle < database->GetNParticles(); currParticle++) {
    ParticlePDG* particle = database->GetPDGParticleByIndex(currParticle);
    meanBaryonDensity += ParticleNumberDensity(particle) * particle->GetBaryonNumber();
  }
  return meanBaryonDensity;
}

double GrandCanonical::StrangeDensity(DatabasePDG* database) {
  // Check if all the thermodinamic parameters are set
  if (!fInitialized)
    Fatal("GrandCanonical::StrangeDensity", "GrandCanonical object not fully initialized!!");

  double meanStrangeDensity = 0.;

  for (int currParticle = 0; currParticle < database->GetNParticles(); currParticle++) {
    ParticlePDG* particle = database->GetPDGParticleByIndex(currParticle);
    meanStrangeDensity += ParticleNumberDensity(particle) * particle->GetStrangeness();
  }

  return meanStrangeDensity;
}

double GrandCanonical::ElectroDensity(DatabasePDG* database) {
  // Check if all the thermodinamic parameters are set
  if (!fInitialized)
    Fatal("GrandCanonical::ElectroDensity", "GrandCanonical object not fully initialized!!");

  double meanElectroDensity = 0.;

  //hadrons
  for (int currParticle = 0; currParticle < database->GetNParticles(); currParticle++) {
    ParticlePDG* particle = database->GetPDGParticleByIndex(currParticle);
    meanElectroDensity += ParticleNumberDensity(particle) * particle->GetElectricCharge();
  }

  return meanElectroDensity;
}

double GrandCanonical::CharmDensity(DatabasePDG* database) {
  // Check if all the thermodinamic parameters are set
  if (!fInitialized)
    Fatal("GrandCanonical::CharmDensity", "GrandCanonical object not fully initialized!!");

  double meanCharmDensity = 0.;

  for (int currParticle = 0; currParticle < database->GetNParticles(); currParticle++) {
    ParticlePDG* particle = database->GetPDGParticleByIndex(currParticle);
    meanCharmDensity += ParticleNumberDensity(particle) * particle->GetCharmness();
  }

  return meanCharmDensity;
}