ThermalModelCanonical.h

Go to the documentation of this file.

/*
 * Thermal-FIST package
 *
 * Copyright (c) 2015-2019 Volodymyr Vovchenko
 *
 * GNU General Public License (GPLv3 or later)
 */
#ifndef THERMALMODELCANONICAL_H
#define THERMALMODELCANONICAL_H
 
#include <map>
 
 
#include "HRGBase/ThermalModelIdeal.h"
 
namespace thermalfist {

  struct QuantumNumbers
  {
    int B; 
    int Q; 
    int S; 
    int C; 

    QuantumNumbers(int iB = 0, int iQ = 0, int iS = 0, int iC = 0) :
      B(iB), Q(iQ), S(iS), C(iC) { }
    const bool operator < (const QuantumNumbers &r) const {
      if (B != r.B)
        return (B < r.B);
      else if (Q != r.Q)
        return (Q < r.Q);
      else if (S != r.S)
        return (S < r.S);
      else
        return (C < r.C);
    }
  };

  enum CanonicalMethod {
    GaussLegendre   = 0,  
    SaddlePoint     = 1,  
    SaddlePointNLO  = 2,  
    SaddlePointLO   = 3   
  };

  class ThermalModelCanonical :
    public ThermalModelBase
  {
  public:
    ThermalModelCanonical(ThermalParticleSystem *TPS, const ThermalModelParameters& params = ThermalModelParameters());

    virtual ~ThermalModelCanonical(void);

    virtual void CalculateQuantumNumbersRange(bool computeFluctuations = false);

    virtual void CalculatePartitionFunctions(double Vc = -1.);

    virtual bool IsParticleCanonical(const ThermalParticle &part);

    virtual void ConserveBaryonCharge(bool conserve = true) { m_BCE = static_cast<int>(conserve); m_PartialZCalculated = false;  }

    virtual void ConserveElectricCharge(bool conserve = true) { m_QCE = static_cast<int>(conserve); m_PartialZCalculated = false; }

    virtual void ConserveStrangeness(bool conserve = true) { m_SCE = static_cast<int>(conserve); m_PartialZCalculated = false; }

    virtual void ConserveCharm(bool conserve = true) { m_CCE = static_cast<int>(conserve); m_PartialZCalculated = false; }
 
    virtual bool IsConservedChargeCanonical(ConservedCharge::Name charge) const;

    virtual double GetGCEDensity(int i) const;

    int IntegrationIterationsMultiplier() const { return m_IntegrationIterationsMultiplier; }

    void SetIntegrationIterationsMultiplier(int multiplier) { (multiplier > 0 ? m_IntegrationIterationsMultiplier = multiplier : m_IntegrationIterationsMultiplier = 1); }

    CanonicalMethod Method() const { return m_Method; }

    void SetMethod(CanonicalMethod method) { m_Method = method; m_PartialZCalculated = false; }

    void SetModelGCE(ThermalModelBase* model);
 
    /*
     * \brief Reset all flags which correspond to a calculation status
     */
    virtual void ResetCalculatedFlags() {
      ThermalModelBase::ResetCalculatedFlags();
      ResetPartialZCalculated();
    }
    /*
     * \brief Reset the flags indicating whether the partial partition functions are calculated
    */
    void ResetPartialZCalculated() { m_PartialZCalculated = false; m_PartialZCalculatedFlucts = false; }
 
    /*
     * \brief Whether the partial partition functions are calculated
     */
    bool IsPartialZCalculated() const { return m_PartialZCalculated; }
 
    /*
     * \brief Whether the partial partition functions are calculated with fluctuations
     */
    bool IsPartialZCalculatedFlucts() const { return m_PartialZCalculatedFlucts; }
 
    // Override functions begin
 
    void ChangeTPS(ThermalParticleSystem *TPS);
 
    virtual void SetStatistics(bool stats);
 
    virtual void FixParameters();
 
    virtual void FixParametersNoReset();

    void SolveChemicalPotentialsGCE(bool resetInitialValues);
 
    virtual void CalculatePrimordialDensities();
 
    virtual void ValidateCalculation();
 
    virtual double ParticleScaledVariance(int part);
 
    virtual void CalculateTwoParticleCorrelations();

    virtual void CalculateFluctuations();
 
    virtual double CalculateEnergyDensity();
 
    virtual double CalculatePressure();
 
    virtual double CalculateEntropyDensity();
 
    virtual double CalculateEnergyDensityDerivativeT() { throw std::runtime_error("CalculateEnergyDensityDerivativeT not implemented"); return 0.; } // Not implemented
 
    virtual double CalculateEntropyDensityDerivativeT() { throw std::runtime_error("CalculateEntropyDensityDerivativeT not implemented"); return 0.; } // Not implemented

    virtual double ParticleScalarDensity(int /*part*/) { return 0.; }
 
    // Override functions end
 
  private:

    virtual void PrepareModelGCE();
    void CleanModelGCE();

    void EnforceSaddlePointLOForNonIdeal();


    void ComputeSaddlePointPartitionFunctions();  
    void SolveSaddlePointEquations();             
    void SetupSaddlePointChargeIndices();          

    double ComputeSigmaLogDetFromModel();

    void ComputeAnalyticCumulants(bool computeCovariance = false);

 
  protected:

    std::vector<QuantumNumbers> m_QNvec;

    std::map<QuantumNumbers, int> m_QNMap;

    std::vector<double> m_Corr;

    std::vector<double> m_PartialZ;

    int m_IntegrationIterationsMultiplier;
 
    CanonicalMethod m_Method; 
 
    int m_BMAX; 
    int m_QMAX; 
    int m_SMAX; 
    int m_CMAX; 
    int m_BMAX_list, m_QMAX_list, m_SMAX_list, m_CMAX_list;
 
    double m_MultExp; 
    double m_MultExpBanalyt; 

    ThermalModelBase *m_modelgce;
    ThermalModelBase *m_modelgce_default; 
    ThermalModelBase *m_modelgce_ext;     
 
    int m_BCE; 
    int m_QCE; 
    int m_SCE; 
    int m_CCE; 

    bool m_Banalyt;
 
    bool m_PartialZCalculated;
    bool m_PartialZCalculatedFlucts;


    int m_SaddlePointDim;    
    std::vector<int> m_SaddlePointChargeIndices;  
    std::vector<double> m_SaddlePointMu;          
    std::vector<double> m_SaddlePointMuStar;      
    double m_SaddlePointLogDetSigma;              
    std::vector<double> m_W2;                      

    double m_NLOEnergyDensity;   
    double m_NLOPressure;        
    double m_NLOEntropyDensity;  
  };
 
} // namespace thermalfist
 
#endif