ThermalModelEVDiagonal.cpp

Go to the documentation of this file.

/*
 * Thermal-FIST package
 * 
 * Copyright (c) 2014-2019 Volodymyr Vovchenko
 *
 * GNU General Public License (GPLv3 or later)
 */
#include "HRGEV/ThermalModelEVDiagonal.h"
 
#ifdef USE_OPENMP
#include <omp.h>
#endif
 
#include <vector>
#include <cmath>
#include <cassert>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>
 
#include "HRGBase/xMath.h"
#include "HRGEV/ExcludedVolumeHelper.h"
 
#include <Eigen/Dense>
 
using namespace Eigen;
 
using namespace std;
 
namespace thermalfist {
 
  ThermalModelEVDiagonal::ThermalModelEVDiagonal(ThermalParticleSystem *TPS, const ThermalModelParameters& params) :
    ThermalModelBase(TPS, params)
  {
    m_densitiesid.resize(m_TPS->Particles().size());
    m_v.resize(m_TPS->Particles().size());
    m_Volume = params.V;
    m_TAG = "ThermalModelEVDiagonal";
 
    m_Ensemble = GCE;
    m_InteractionModel = DiagonalEV;
  }
 
 
  ThermalModelEVDiagonal::~ThermalModelEVDiagonal(void)
  {
  }
 
 
  void ThermalModelEVDiagonal::SetRadius(double rad) {
    if (m_v.size() != m_TPS->Particles().size())
      m_v.resize(m_TPS->Particles().size());
    for (int i = 0; i < m_TPS->ComponentsNumber(); ++i) {
      m_v[i] = CuteHRGHelper::vr(rad);
    }
  }
 
  void ThermalModelEVDiagonal::FillVirial(const std::vector<double>& ri)
  {
    if (ri.size() != m_TPS->Particles().size()) {
      std::cerr << "**WARNING** " << m_TAG << "::FillVirial(const std::vector<double> & ri): size of ri does not match number of hadrons in the list" << std::endl;
      return;
    }
    m_v.resize(m_TPS->Particles().size());
    for (int i = 0; i < m_TPS->ComponentsNumber(); ++i)
      m_v[i] = CuteHRGHelper::vr(ri[i]);
  }
 
  void ThermalModelEVDiagonal::FillVirialEV(const std::vector<double>& vi)
  {
    if (vi.size() != m_TPS->Particles().size()) {
      std::cerr << "**WARNING** " << m_TAG << "::FillVirialEV(const std::vector<double> & vi): size of vi does not match number of hadrons in the list" << std::endl;
      return;
    }
    m_v = vi;
  }
 
  void ThermalModelEVDiagonal::ReadInteractionParameters(const std::string & filename)
  {
    m_v = std::vector<double>(m_TPS->Particles().size(), 0.);
 
    ifstream fin(filename.c_str());
    char cc[2000];
    while (!fin.eof()) {
      fin.getline(cc, 2000);
      string tmp = string(cc);
      vector<string> elems = CuteHRGHelper::split(tmp, '#');
      if (elems.size() < 1)
        continue;
      istringstream iss(elems[0]);
      long long pdgid;
      double b;
      if (iss >> pdgid >> b) {
        int ind = m_TPS->PdgToId(pdgid);
        if (ind != -1)
          m_v[ind] = b;
      }
    }
    fin.close();
  }
 
  void ThermalModelEVDiagonal::WriteInteractionParameters(const std::string & filename)
  {
    ofstream fout(filename.c_str());
    fout << "# List of eigenvolume parameters to be used in the Diagonal excluded-volume HRG model"
      << std::endl;
    fout << "# Only particles with a non-zero eigenvolume parameter are listed here"
      << std::endl;
    fout << "#" << std::setw(14) << "pdg"
      << std::setw(15) << "v_i[fm^3]"
      << std::endl;
    for (int i = 0; i < m_TPS->ComponentsNumber(); ++i) {
      if (m_v[i] != 0.) {
        fout << std::setw(15) << m_TPS->Particle(i).PdgId();
        fout << std::setw(15) << m_v[i];
        fout << std::endl;
      }
    }
    fout.close();
  }
 
  double ThermalModelEVDiagonal::ExcludedVolume(int i) const
  {
    if (i < 0 || i >= static_cast<int>(m_v.size()))
      return 0.;
    return m_v[i];
  }
 
 
  void ThermalModelEVDiagonal::ChangeTPS(ThermalParticleSystem *TPS) {
    ThermalModelBase::ChangeTPS(TPS);
    m_densitiesid.resize(m_TPS->Particles().size());
  }
 
 
  double ThermalModelEVDiagonal::DensityId(int i, double Pressure) {
    double dMu = -m_v[i] * Pressure;
 
    return m_TPS->Particles()[i].Density(m_Parameters, IdealGasFunctions::ParticleDensity, m_UseWidth, m_Chem[i] + dMu);
  }
 
  double ThermalModelEVDiagonal::PressureId(int i, double Pressure) {
    double dMu = -m_v[i] * Pressure;
 
    return m_TPS->Particles()[i].Density(m_Parameters, IdealGasFunctions::Pressure, m_UseWidth, m_Chem[i] + dMu);
  }
 
  double ThermalModelEVDiagonal::ScaledVarianceId(int i, double Pressure) {
    double dMu = -m_v[i] * Pressure;
 
    return m_TPS->Particles()[i].ScaledVariance(m_Parameters, m_UseWidth, m_Chem[i] + dMu);
  }
 
  double ThermalModelEVDiagonal::Pressure(double P) {
    double ret = 0.;
 
    for (int i = 0; i < m_TPS->ComponentsNumber(); ++i) {
      double dMu = -m_v[i] * P;
      ret += m_TPS->Particles()[i].Density(m_Parameters, IdealGasFunctions::Pressure, m_UseWidth, m_Chem[i] + dMu);
    }
 
    return ret;
  }
 
  void ThermalModelEVDiagonal::SolvePressure() {
    BroydenEquationsDEV eqs(this);
    BroydenJacobianDEV jac(this);
    Broyden broydn(&eqs, &jac);
    BroydenSolutionCriteriumDEV crit(this, 1.0E-8);
 
    m_Pressure = Pressure(0.);
    double mnc = pow(m_Parameters.T, 4.) * pow(xMath::GeVtoifm(), 3.);
    if (m_Parameters.T == 0.0)
      mnc = pow(xMath::mnucleon(), 4.) * pow(xMath::GeVtoifm(), 3.);
    eqs.SetMnc(mnc);
    jac.SetMnc(mnc);
    double x0 = log(m_Pressure / mnc);
    std::vector<double> x(1, x0);
 
    x = broydn.Solve(x, &crit);
 
    
 
    double PressureNew = mnc * exp(x[0]);
 
    // UPDATE: Currently not used
    //double current_precision = Broyden::TOL;
    //while (abs(PressureNew) < current_precision && current_precision > 1.e-50 && abs(PressureNew /m_Pressure) < 1.e-5) {
    //  current_precision *= 1.e-10;
    //  x = broydn.Solve(x, &Broyden::BroydenSolutionCriterium(current_precision));
    //  PressureNew = mnc * exp(x[0]);
    //}
 
    m_Pressure = PressureNew;
 
    if (broydn.Iterations() == broydn.MaxIterations())
      m_LastCalculationSuccessFlag = false;
    else m_LastCalculationSuccessFlag = true;
 
    m_MaxDiff = broydn.MaxDifference();
  }
 
  void ThermalModelEVDiagonal::CalculatePrimordialDensities() {
    m_FluctuationsCalculated = false;
 
    SolvePressure();
 
    m_wnSum = 0.;
    m_Densityid = 0.;
    m_TotalDensity = 0.;
    m_Suppression = 0.;
    double densityid = 0., suppression = 0.;
 
    m_densitiesidnoshift = m_densitiesid;
 
#pragma omp parallel for reduction(+:densityid) reduction(+:suppression) if(useOpenMP)
    for (int i = 0; i < m_TPS->ComponentsNumber(); ++i) {
      double dMu = -m_v[i] * m_Pressure;
      m_densitiesid[i] = m_TPS->Particles()[i].Density(m_Parameters, IdealGasFunctions::ParticleDensity, m_UseWidth, m_Chem[i] + dMu);
      densityid += m_densitiesid[i];
      suppression += m_v[i] * m_densitiesid[i];
 
      m_densitiesidnoshift[i] = m_TPS->Particles()[i].Density(m_Parameters, IdealGasFunctions::ParticleDensity, m_UseWidth, m_Chem[i]);
    }
 
    m_Densityid = densityid;
    m_Suppression = suppression;
 
    m_Suppression = 1. / (1. + m_Suppression);
 
    for (int i = 0; i < m_TPS->ComponentsNumber(); ++i) {
      m_densities[i] = m_densitiesid[i] * m_Suppression;
      m_TotalDensity += m_densities[i];
      m_wnSum += m_densities[i] * m_TPS->Particles()[i].ScaledVariance(m_Parameters, m_UseWidth, m_Chem[i] - m_v[i] * m_Pressure);
    }
 
    CalculateFeeddown();
 
    m_Calculated = true;
    ValidateCalculation();
  }
 
  void ThermalModelEVDiagonal::CalculateTwoParticleCorrelations() {
    int NN = m_densities.size();
    vector<double> tN(NN);
    for (int i = 0; i < NN; ++i) tN[i] = DensityId(i, m_Pressure);
 
    vector<double> chi2id(NN);
    for (int i = 0; i < NN; ++i) {
      chi2id[i] = m_TPS->Particles()[i].chiDimensionfull(2, m_Parameters, m_UseWidth, m_Chem[i] - m_v[i] * m_Pressure) * xMath::GeVtoifm3();
      if (tN[i] > 0.0)
        chi2id[i] *= m_densities[i] / tN[i];
    }
 
    m_PrimCorrel.resize(NN);
    for (int i = 0; i < NN; ++i) m_PrimCorrel[i].resize(NN);
    m_dmusdmu = m_PrimCorrel;
    m_TotalCorrel = m_PrimCorrel;
 
    for (int i = 0; i < NN; ++i) {
      for (int j = 0; j < NN; ++j) {
        m_PrimCorrel[i][j] = 0.;
        if (i == j) m_PrimCorrel[i][j] += chi2id[i];
        m_PrimCorrel[i][j] += -m_v[i] * m_densities[j] * chi2id[i];
        m_PrimCorrel[i][j] += -m_v[j] * m_densities[i] * chi2id[j];
        double tmp = 0.;
        for (size_t k = 0; k < m_densities.size(); ++k)
          tmp += m_v[k] * m_v[k] * chi2id[k];
        m_PrimCorrel[i][j] += m_densities[i] * m_densities[j] * tmp;
 
        m_dmusdmu[i][j] = (i == j ? 1. : 0.) - m_v[i] * m_densities[j];
      }
    }
 
    for (int i = 0; i < NN; ++i) {
      m_wprim[i] = m_PrimCorrel[i][i];
      if (m_densities[i] > 0.) m_wprim[i] *= m_Parameters.T / m_densities[i];
      else m_wprim[i] = 1.;
      if (m_wprim[i] != m_wprim[i]) m_wprim[i] = 1.;
    }
 
    m_TwoParticleCorrelationsCalculated = true;
  }
 
  double ThermalModelEVDiagonal::CalculateEnergyDensityDerivativeT() {
    if (!IsTemperatureDerivativesCalculated())
      CalculateTemperatureDerivatives();
 
    // Compute dE/dT
    double ret = 0.;
 
    double eid = 0., dTbn = 0., bn = 0., dTeidterm = 0., dmueidterm = 0.;
    for (int i = 0; i < m_TPS->ComponentsNumber(); ++i) {
      const ThermalParticle &part = m_TPS->Particles()[i];
      eid += part.Density(m_Parameters, IdealGasFunctions::EnergyDensity, m_UseWidth, m_Chem[i] - m_v[i] * m_Pressure);
      dTbn += m_v[i] * m_dndT[i];
      bn += m_v[i] * m_densities[i];
      dTeidterm += part.Density(m_Parameters, IdealGasFunctions::dedT, m_UseWidth, m_Chem[i] - m_v[i] * m_Pressure);
      dmueidterm += -m_v[i] * EntropyDensity() * part.Density(m_Parameters, IdealGasFunctions::dedmu, m_UseWidth, m_Chem[i] - m_v[i] * m_Pressure);
    }
 
    ret = -dTbn * eid + (1. - bn) * (dTeidterm + dmueidterm);
 
    return ret;
  }
 
  double ThermalModelEVDiagonal::CalculateEntropyDensityDerivativeT() {
    if (!IsTemperatureDerivativesCalculated())
      CalculateTemperatureDerivatives();
 
    // Compute ds/dT
    double ret = 0.;
 
    double sid = 0., dTbn = 0., bn = 0., dTsidterm = 0., dmusidterm = 0.;
    for (int i = 0; i < m_TPS->ComponentsNumber(); ++i) {
      const ThermalParticle &part = m_TPS->Particles()[i];
      sid += part.Density(m_Parameters, IdealGasFunctions::EntropyDensity, m_UseWidth, m_Chem[i] - m_v[i] * m_Pressure);
      dTbn += m_v[i] * m_dndT[i];
      bn += m_v[i] * m_densities[i];
      dTsidterm += part.Density(m_Parameters, IdealGasFunctions::dsdT, m_UseWidth, m_Chem[i] - m_v[i] * m_Pressure);
      // ds/dmu = dn/dT (Maxwell relation), dmu*/dT = -b_i * s
      dmusidterm += -m_v[i] * EntropyDensity() * part.Density(m_Parameters, IdealGasFunctions::dndT, m_UseWidth, m_Chem[i] - m_v[i] * m_Pressure);
    }
 
    ret = -dTbn * sid + (1. - bn) * (dTsidterm + dmusidterm);
 
    return ret;
  }
 
  void ThermalModelEVDiagonal::CalculateTemperatureDerivatives() {
    int N = m_TPS->ComponentsNumber();
    m_dndT = vector<double>(N, 0.);
    m_dmusdT = vector<double>(N, 0.);
    m_PrimChi2sdT = vector<vector<double>>(N, vector<double>(N, 0.));
 
    double T = m_Parameters.T;
 
    vector<double> nids(N, 0.), dniddTs(N, 0.), chi2ids(N, 0.), dchi2idsdT(N, 0.), chi3ids(N, 0.);
    for (int i = 0; i < N; ++i) {
      nids[i] = m_TPS->Particles()[i].Density(m_Parameters, IdealGasFunctions::ParticleDensity, m_UseWidth, m_Chem[i] - m_v[i] * m_Pressure);
      dniddTs[i] = m_TPS->Particles()[i].Density(m_Parameters, IdealGasFunctions::dndT, m_UseWidth, m_Chem[i] - m_v[i] * m_Pressure);
      chi2ids[i] = m_TPS->Particles()[i].chiDimensionfull(2, m_Parameters, m_UseWidth, m_Chem[i] - m_v[i] * m_Pressure) * xMath::GeVtoifm3();
      // dchi2idsdT = d(dn^id/dmu*)/dT = T^2 * dchi2/dT + 2 * chi2ids / T
      // The two terms individually diverge as 1/T but their sum is O(T).
      // At T=0 the result is exactly zero (Sommerfeld: dn/dmu = const + O(T^2)).
      if (T > 0.)
        dchi2idsdT[i] = (T * T * m_TPS->Particles()[i].Density(m_Parameters, IdealGasFunctions::dchi2dT, m_UseWidth, m_Chem[i] - m_v[i] * m_Pressure) * xMath::GeVtoifm3() + 2. * chi2ids[i] / T);
      else
        dchi2idsdT[i] = 0.;
      chi3ids[i] = m_TPS->Particles()[i].chiDimensionfull(3, m_Parameters, m_UseWidth, m_Chem[i] - m_v[i] * m_Pressure) * xMath::GeVtoifm3();
    }
 
    double s = EntropyDensity();
    double sum_chi2idb2 = 0., sum_bns= 0., sum_bn = 0., sum_dTbns = 0., sum_dchi2dTidb2 = 0., sum_chi3idb3 = 0.;
    for (int i = 0; i < N; ++i) {
      m_dmusdT[i] = -m_v[i] * s;
      sum_chi2idb2 += m_v[i] * m_v[i] * chi2ids[i];
      sum_bns += m_v[i] * nids[i];
      sum_bn += m_v[i] * m_densities[i];
      sum_dTbns += m_v[i] * dniddTs[i];
      sum_dchi2dTidb2 += m_v[i] * m_v[i] * dchi2idsdT[i];
      sum_chi3idb3 += m_v[i] * m_v[i] * m_v[i] * chi3ids[i];
    }
 
    double sum_dbndT = (1. - sum_bn) * (1. - sum_bn) * (sum_dTbns - sum_chi2idb2 * s);
 
    for (int i = 0; i < N; ++i) {
      m_dndT[i] = -sum_dbndT * nids[i] + (1. - sum_bn) * (dniddTs[i] + chi2ids[i] * m_dmusdT[i]);
    }
 
    if (IsSusceptibilitiesCalculated()) {
      for (int i = 0; i < N; ++i) {
        for (int j = 0; j < N; ++j) {
          m_PrimChi2sdT[i][j] = 0.;
 
          double tval = 0.;
          tval += -chi2ids[j] * m_densities[i] * m_v[j];
          tval += sum_chi2idb2 * m_densities[j] * m_densities[i];
          if (i == j)
            tval += chi2ids[i];
          tval += -chi2ids[i] * m_densities[j] * m_v[i];
          m_PrimChi2sdT[i][j] += -sum_dbndT * tval;
 
          tval = 0.;
          tval += -dchi2idsdT[j] * m_densities[i] * m_v[j];
          tval += -chi3ids[j] * m_dmusdT[j] * m_densities[i] * m_v[j];
          tval += -chi2ids[j] * m_dndT[i] * m_v[j];
          tval += sum_dchi2dTidb2 * m_densities[j] * m_densities[i];
          tval += -sum_chi3idb3 * s * m_densities[j] * m_densities[i];
          tval += sum_chi2idb2 * (m_dndT[j] * m_densities[i] + m_densities[j] * m_dndT[i]);
          if (i == j) {
            tval += dchi2idsdT[i];
            tval += chi3ids[i] * m_dmusdT[i];
          }
          tval += -dchi2idsdT[i] * m_densities[j] * m_v[i];
          tval += -chi3ids[i] * m_dmusdT[i] * m_densities[j] * m_v[i];
          tval += -chi2ids[i] * m_dndT[j] * m_v[i];
 
          m_PrimChi2sdT[i][j] += (1. - sum_bn) * tval;
 
          // dchi2/dT = dchi2til/dT / T^2  - 2 * chi2til / T^3
          // Both terms diverge as T -> 0; guard against division by zero.
          if (T > 0.)
            m_PrimChi2sdT[i][j] = m_PrimChi2sdT[i][j] / (xMath::GeVtoifm3() * T * T) - 2. * m_PrimCorrel[i][j] / T / T / T / xMath::GeVtoifm3();
          else
            m_PrimChi2sdT[i][j] = 0.;
        }
      }
 
      m_TemperatureDerivativesCalculated = true;
      CalculateSusceptibilityMatrix();
    }
 
    m_TemperatureDerivativesCalculated = true;
  }
 
  // TODO include correlations
  void ThermalModelEVDiagonal::CalculateFluctuations() {
    CalculateTwoParticleCorrelations();
    CalculateSusceptibilityMatrix();
    CalculateTwoParticleFluctuationsDecays();
    CalculateProxySusceptibilityMatrix();
    CalculateParticleChargeCorrelationMatrix();
 
    m_FluctuationsCalculated = true;
 
    for (size_t i = 0; i < m_wprim.size(); ++i) {
      m_skewprim[i] = ParticleSkewness(i);
      m_kurtprim[i] = ParticleKurtosis(i);
    }
 
    for (size_t i = 0; i < m_wprim.size(); ++i) {
      m_skewtot[i] = 1.;
      m_kurttot[i] = 1.;
    }
 
    if (IsTemperatureDerivativesCalculated()) {
      m_TemperatureDerivativesCalculated = false;
      CalculateTemperatureDerivatives();
    }
  }
 
 
  std::vector<double> ThermalModelEVDiagonal::CalculateChargeFluctuations(const std::vector<double>& chgs, int order, bool dimensionfull)
  {
    vector<double> ret(order + 1, 0.);
 
    // chi1
    for (size_t i = 0; i < m_densities.size(); ++i)
      ret[0] += chgs[i] * m_densities[i];
 
    if (!dimensionfull)
      ret[0] /= pow(m_Parameters.T * xMath::GeVtoifm(), 3);
    else
      ret[0] /= pow(xMath::GeVtoifm(), 3);
 
    if (order < 2) return ret;
 
 
 
    // Preparing matrix for system of linear equations
    int NN = m_densities.size();
 
    vector<double> MuStar(NN, 0.);
    for (int i = 0; i < NN; ++i) {
      MuStar[i] = m_Chem[i] + MuShift(i);
    }
 
 
    MatrixXd densMatrix(2 * NN, 2 * NN);
    VectorXd solVector(2 * NN), xVector(2 * NN);
 
    vector<double> DensitiesId(m_densities.size()), chi2id(m_densities.size());
    for (int i = 0; i < NN; ++i) {
      DensitiesId[i] = m_TPS->Particles()[i].Density(m_Parameters, IdealGasFunctions::ParticleDensity, m_UseWidth, MuStar[i]);
      chi2id[i] = m_TPS->Particles()[i].chiDimensionfull(2, m_Parameters, m_UseWidth, MuStar[i]) * pow(xMath::GeVtoifm(), 3);
    }
 
    for (int i = 0; i < NN; ++i)
      for (int j = 0; j < NN; ++j) {
        densMatrix(i, j) = m_v[j] * DensitiesId[i];
        if (i == j) densMatrix(i, j) += 1.;
      }
 
    for (int i = 0; i < NN; ++i)
      for (int j = 0; j < NN; ++j)
        densMatrix(i, NN + j) = 0.;
 
    for (int i = 0; i < NN; ++i) {
      densMatrix(i, NN + i) = 0.;
      for (int k = 0; k < NN; ++k) {
        densMatrix(i, NN + i) += m_v[k] * m_densities[k];
      }
      densMatrix(i, NN + i) = (densMatrix(i, NN + i) - 1.) * chi2id[i];
    }
 
    for (int i = 0; i < NN; ++i)
      for (int j = 0; j < NN; ++j) {
        densMatrix(NN + i, j) = 0.;
        densMatrix(NN + i, NN + j) = m_v[i] * DensitiesId[j];
        if (i == j) densMatrix(NN + i, NN + j) += 1.;
      }
 
 
    PartialPivLU<MatrixXd> decomp(densMatrix);
 
    // chi2
    vector<double> dni(NN, 0.), dmus(NN, 0.);
 
    for (int i = 0; i < NN; ++i) {
      xVector[i] = 0.;
      xVector[NN + i] = chgs[i];
    }
 
    solVector = decomp.solve(xVector);
 
    for (int i = 0; i < NN; ++i) {
      dni[i] = solVector[i];
      dmus[i] = solVector[NN + i];
    }
 
    for (int i = 0; i < NN; ++i)
      ret[1] += chgs[i] * dni[i];
 
    if (!dimensionfull)
      ret[1] /= pow(m_Parameters.T, 2) * pow(xMath::GeVtoifm(), 3);
    else
      ret[1] /= pow(xMath::GeVtoifm(), 3);
 
    if (order < 3) return ret;
    // chi3
    vector<double> d2ni(NN, 0.), d2mus(NN, 0.);
 
    vector<double> chi3id(m_densities.size());
    for (int i = 0; i < NN; ++i)
      chi3id[i] = m_TPS->Particles()[i].chiDimensionfull(3, m_Parameters, m_UseWidth, MuStar[i]) * pow(xMath::GeVtoifm(), 3);
 
    for (int i = 0; i < NN; ++i) {
      xVector[i] = 0.;
 
      double tmp = 0.;
      for (int j = 0; j < NN; ++j) tmp += m_v[j] * dni[j];
      tmp = -2. * tmp * chi2id[i] * dmus[i];
      xVector[i] += tmp;
 
      tmp = 0.;
      for (int j = 0; j < NN; ++j) tmp += m_v[j] * m_densities[j];
      tmp = -(tmp - 1.) * chi3id[i] * dmus[i] * dmus[i];
      xVector[i] += tmp;
    }
    for (int i = 0; i < NN; ++i) {
      xVector[NN + i] = 0.;
 
      double tmp = 0.;
      for (int j = 0; j < NN; ++j) tmp += -m_v[i] * dmus[j] * chi2id[j] * dmus[j];
 
      xVector[NN + i] = tmp;
    }
 
    solVector = decomp.solve(xVector);
 
    for (int i = 0; i < NN; ++i) {
      d2ni[i] = solVector[i];
      d2mus[i] = solVector[NN + i];
    }
 
    for (int i = 0; i < NN; ++i)
      ret[2] += chgs[i] * d2ni[i];
 
    if (!dimensionfull)
      ret[2] /= m_Parameters.T * pow(xMath::GeVtoifm(), 3);
    else
      ret[2] /= pow(xMath::GeVtoifm(), 3);
 
 
    if (order < 4) return ret;
 
    // chi4
    vector<double> d3ni(NN, 0.), d3mus(NN, 0.);
 
    vector<double> chi4id(m_densities.size());
    for (int i = 0; i < NN; ++i)
      chi4id[i] = m_TPS->Particles()[i].chiDimensionfull(4, m_Parameters, m_UseWidth, MuStar[i]) * pow(xMath::GeVtoifm(), 3);
 
    vector<double> dnis(NN, 0.);
    for (int i = 0; i < NN; ++i) {
      dnis[i] = chi2id[i] * dmus[i];
    }
 
    vector<double> d2nis(NN, 0.);
    for (int i = 0; i < NN; ++i) {
      d2nis[i] = chi3id[i] * dmus[i] * dmus[i] +
        chi2id[i] * d2mus[i];
    }
 
    for (int i = 0; i < NN; ++i) {
      xVector[i] = 0.;
 
      double tmp = 0.;
      for (int j = 0; j < NN; ++j) tmp += m_v[j] * dni[j];
      tmp = -3. * tmp * d2nis[i];
      xVector[i] += tmp;
 
      tmp = 0.;
      for (int j = 0; j < NN; ++j) tmp += m_v[j] * d2ni[j];
      tmp = -3. * tmp * dnis[i];
      xVector[i] += tmp;
 
      double tmps = 0.;
      for (int j = 0; j < NN; ++j) tmps += m_v[j] * m_densities[j];
 
      tmp = -(tmps - 1.) * chi3id[i] * d2mus[i] * 3. * dmus[i];
      xVector[i] += tmp;
 
      tmp = -(tmps - 1.) * chi4id[i] * dmus[i] * dmus[i] * dmus[i];
      xVector[i] += tmp;
    }
    for (int i = 0; i < NN; ++i) {
      xVector[NN + i] = 0.;
 
      double tmp = 0.;
      for (int j = 0; j < NN; ++j) tmp += -2. * m_v[i] * d2mus[j] * dnis[j];
      xVector[NN + i] += tmp;
 
      tmp = 0.;
      for (int j = 0; j < NN; ++j) tmp += -m_v[i] * dmus[j] * d2nis[j];
      xVector[NN + i] += tmp;
    }
 
    solVector = decomp.solve(xVector);
 
    for (int i = 0; i < NN; ++i) {
      d3ni[i] = solVector[i];
      d3mus[i] = solVector[NN + i];
    }
 
    for (int i = 0; i < NN; ++i)
      ret[3] += chgs[i] * d3ni[i];
 
    ret[3] /= pow(xMath::GeVtoifm(), 3);
 
    return ret;
  }
 
 
  double ThermalModelEVDiagonal::CalculateEnergyDensity() {
    if (!m_Calculated) CalculateDensities();
    double ret = 0.;
    double dMu = 0.;
    for (int i = 0; i < m_TPS->ComponentsNumber(); ++i) {
      dMu = -m_v[i] * m_Pressure;
      ret += m_TPS->Particles()[i].Density(m_Parameters, IdealGasFunctions::EnergyDensity, m_UseWidth, m_Chem[i] + dMu);
    }
    return ret * m_Suppression;
  }
 
  double ThermalModelEVDiagonal::CalculateEntropyDensity() {
    if (!m_Calculated) CalculateDensities();
    double ret = 0.;
    double dMu = 0.;
    for (int i = 0; i < m_TPS->ComponentsNumber(); ++i) {
      dMu = -m_v[i] * m_Pressure;
      ret += m_TPS->Particles()[i].Density(m_Parameters, IdealGasFunctions::EntropyDensity, m_UseWidth, m_Chem[i] + dMu);
    }
    return ret * m_Suppression;
  }
 
  double ThermalModelEVDiagonal::CalculatePressure() {
    if (!m_Calculated) CalculateDensities();
    return m_Pressure;
  }
 
  double ThermalModelEVDiagonal::ParticleScalarDensity(int part) {
    if (!m_Calculated) CalculateDensities();
 
    double dMu = -m_v[part] * m_Pressure;
    double ret = m_TPS->Particles()[part].Density(m_Parameters, IdealGasFunctions::ScalarDensity, m_UseWidth, m_Chem[part] + dMu);
    return ret * m_Suppression;
  }
 
  double ThermalModelEVDiagonal::CommonSuppressionFactor()
  {
    if (!m_Calculated)
      CalculateDensities();
    return m_Suppression;
  }
 
  double ThermalModelEVDiagonal::MuShift(int id) const
  {
    if (id >= 0. && id < static_cast<int>(m_v.size()))
      return -m_v[id] * m_Pressure;
    else
      return 0.0;
  }
 
  double ThermalModelEVDiagonal::CalculateEigenvolumeFraction() {
    double tEV = 0.;
    for (int i = 0; i < m_TPS->ComponentsNumber(); ++i) {
      tEV += m_v[i] * m_densities[i] / 4.;
    }
    return tEV;
  }
 
  void ThermalModelEVDiagonal::SetRadius(int i, double rad)
  {
    if (i >= 0 && i < static_cast<int>(m_v.size()))
      m_v[i] = CuteHRGHelper::vr(rad);
  }
 
  double ThermalModelEVDiagonal::VirialCoefficient(int i, int /*j*/) const
  {
    return m_v[i];
  }
 
  void ThermalModelEVDiagonal::SetVirial(int i, int j, double b)
  {
    if (i == j)
      m_v[i] = b;
  }
 
  std::vector<double> ThermalModelEVDiagonal::BroydenEquationsDEV::Equations(const std::vector<double>& x)
  {
    std::vector<double> ret(1);
    double pressure = m_mnc * exp(x[0]);
 
    ret[0] = pressure - m_THM->Pressure(pressure);
 
    return ret;
  }
 
  std::vector<double> ThermalModelEVDiagonal::BroydenJacobianDEV::Jacobian(const std::vector<double>& x)
  {
    double pressure = m_mnc * exp(x[0]);
 
    double ret = 0.;
    for (size_t i = 0; i < m_THM->Densities().size(); ++i)
      ret += m_THM->VirialCoefficient(i, i) * m_THM->DensityId(i, pressure);
    ret += 1.;
    ret *= pressure;
 
    return std::vector<double>(1, ret);
  }
 
  std::vector<double> ThermalModelEVDiagonal::BroydenEquationsDEVOrig::Equations(const std::vector<double>& x)
  {
    std::vector<double> ret(1);
    ret[0] = x[0] - m_THM->Pressure(x[0]);
    return ret;
  }
 
  std::vector<double> ThermalModelEVDiagonal::BroydenJacobianDEVOrig::Jacobian(const std::vector<double>& x)
  {
    const double &pressure = x[0];
 
    double ret = 0.;
    for (size_t i = 0; i < m_THM->Densities().size(); ++i)
      ret += m_THM->VirialCoefficient(i, i) * m_THM->DensityId(i, pressure);
    ret += 1.;
 
    return std::vector<double>(1, ret);
  }
 
  bool ThermalModelEVDiagonal::BroydenSolutionCriteriumDEV::IsSolved(const std::vector<double>& x, const std::vector<double>& f, const std::vector<double>& /*xdelta*/) const
  {
    double maxdiff = 0.;
    for (size_t i = 0; i < x.size(); ++i) {
      maxdiff = std::max(maxdiff, fabs(f[i]) / m_THM->m_Pressure);
    }
    return (maxdiff < m_MaximumError);
  }
} // namespace thermalfist