ThermalModelBase.h

Go to the documentation of this file.

/*
 * Thermal-FIST package
 * 
 * Copyright (c) 2014-2019 Volodymyr Vovchenko
 *
 * GNU General Public License (GPLv3 or later)
 */
#ifndef THERMALMODELBASE_H
#define THERMALMODELBASE_H
 
#include <string>
 
#include "HRGBase/ThermalParticleSystem.h"
#include "HRGBase/xMath.h"
#include "HRGBase/Broyden.h"
 
 
namespace thermalfist {

  class ThermalModelBase
  {
  public:
    enum ThermalModelEnsemble { 
      GCE = 0, 
      CE = 1,  
      SCE = 2, 
      CCE = 3 
    };

    enum ThermalModelInteraction { 
      Ideal = 0,        
      DiagonalEV = 1,   
      CrosstermsEV = 2, 
      QvdW = 3,         
      RealGas = 4,      
      MeanField = 5     
    };

    ThermalModelBase(ThermalParticleSystem *TPS, const ThermalModelParameters& params = ThermalModelParameters());
 
    virtual ~ThermalModelBase(void) { }

    int ComponentsNumber() const { return static_cast<int>(m_densities.size()); }

    virtual void FillVirial(const std::vector<double> & ri = std::vector<double>(0));

    bool UseWidth() const { return m_UseWidth; }

    void SetUseWidth(bool useWidth);

    void SetUseWidth(ThermalParticle::ResonanceWidthIntegration type);


    void SetNormBratio(bool normBratio);
    bool NormBratio() const { return m_NormBratio; }

    void SetOMP(bool openMP) { m_useOpenMP = openMP; }
     

    virtual void SetParameters(const ThermalModelParameters& params);
    const ThermalModelParameters& Parameters() const { return m_Parameters; }

    void UpdateParameters() { SetParameters(m_Parameters); }

    virtual void SetTemperature(double T);

    virtual void SetBaryonChemicalPotential(double muB);

    virtual void SetElectricChemicalPotential(double muQ);

    virtual void SetStrangenessChemicalPotential(double muS);

    virtual void SetCharmChemicalPotential(double muC);

    virtual void SetGammaq(double gammaq);

    virtual void SetGammaS(double gammaS);

    virtual void SetGammaC(double gammaC);

    virtual void SetBaryonCharge(int B);

    virtual void SetElectricCharge(int Q);

    virtual void SetStrangeness(int S);

    virtual void SetCharm(int C);

    virtual void SetRadius(double /*rad*/) { }

    virtual void SetRadius(int /*i*/, double /*rad*/) { }

    virtual void SetVirial(int /*i*/, int /*j*/, double /*b*/) { }

    virtual void SetRepulsion(int i, int j, double b) { SetVirial(i,j,b); }

    virtual void SetAttraction(int /*i*/, int /*j*/, double /*a*/) { }


    virtual void DisableMesonMesonVirial();
    void DisableMesonMesonRepulsion() { return DisableMesonMesonVirial(); }

    virtual void DisableMesonMesonAttraction();


    virtual void DisableMesonBaryonVirial();
    void DisableMesonBaryonRepulsion() { return DisableMesonBaryonVirial(); }

    virtual void DisableMesonBaryonAttraction();


    virtual void DisableBaryonBaryonVirial();
    void DisableBaryonBaryonRepulsion() { return DisableBaryonBaryonVirial(); }

    virtual void DisableBaryonBaryonAttraction();


    virtual void DisableBaryonAntiBaryonVirial();
    void DisableBaryonAntiBaryonRepulsion() { return DisableBaryonAntiBaryonVirial(); }

    virtual void DisableBaryonAntiBaryonAttraction();

    virtual void ReadInteractionParameters(const std::string & /*filename*/) { }

    virtual void WriteInteractionParameters(const std::string & /*filename*/) { }

    virtual void ChangeTPS(ThermalParticleSystem *TPS);


    virtual double VirialCoefficient(int /*i*/, int /*j*/) const { return 0.; }
    double RepulsionCoefficient(int i, int j) const { return VirialCoefficient(i,j); }

    virtual double AttractionCoefficient(int /*i*/, int /*j*/) const { return 0.; }

    bool QuantumStatistics() const { return m_QuantumStats; }

    virtual void SetStatistics(bool stats);

    virtual void SetCalculationType(IdealGasFunctions::QStatsCalculationType type) { m_TPS->SetCalculationType(type); }
    
    virtual void SetClusterExpansionOrder(int order) { m_TPS->SetClusterExpansionOrder(order); }
    
    void SetResonanceWidthShape(ThermalParticle::ResonanceWidthShape shape) { m_TPS->SetResonanceWidthShape(shape); }
    
    void SetResonanceWidthIntegrationType(ThermalParticle::ResonanceWidthIntegration type);// { m_TPS->SetResonanceWidthIntegrationType(type); }

    virtual void FillChemicalPotentials();

    virtual void SetChemicalPotentials(const std::vector<double> & chem = std::vector<double>(0));

    const std::vector<double>& ChemicalPotentials() const { return m_Chem; }

    double ChemicalPotential(int i) const;

    virtual void SetChemicalPotential(int i, double chem);

    virtual double FullIdealChemicalPotential(int i) const;
 

    bool ConstrainMuB() const { return m_ConstrainMuB; }

    void ConstrainMuB(bool constrain) { m_ConstrainMuB = constrain; }

    bool ConstrainMuQ() const { return m_ConstrainMuQ; }

    void ConstrainMuQ(bool constrain) { m_ConstrainMuQ = constrain; }

    bool ConstrainMuS() const { return m_ConstrainMuS; }

    void ConstrainMuS(bool constrain) { m_ConstrainMuS = constrain; }

    bool ConstrainMuC() const { return m_ConstrainMuC; }

    void ConstrainMuC(bool constrain) { m_ConstrainMuC = constrain; }

    void UsePartialChemicalEquilibrium(bool usePCE) { m_PCE = usePCE; }

    bool UsePartialChemicalEquilibrium() { return m_PCE; }


    void SetSoverB(double SB) { m_SBgoal = SB; }
    double SoverB() const { return m_SBgoal; }


    void SetQoverB(double QB) { m_QBgoal = QB; }
    double QoverB() const { return m_QBgoal; }

    void SetVolume(double Volume) { m_Volume = Volume; m_Parameters.V = Volume; }
    double Volume() const { return m_Parameters.V; }
    
    void SetVolumeRadius(double R) { m_Volume = 4. / 3.*xMath::Pi() * R * R * R; m_Parameters.V = m_Volume; }

    double CanonicalVolume() const { return m_Parameters.SVc; }
 
    void SetCanonicalVolume(double Volume) { m_Parameters.SVc = Volume; }

    void SetCanonicalVolumeRadius(double Rc) { m_Parameters.SVc = 4. / 3. * xMath::Pi() * Rc * Rc * Rc; }
 
 
    //double StrangenessCanonicalVolume() const { return CanonicalVolume(); }
    //void SetStrangenessCanonicalVolume(double Volume) { SetCanonicalVolume(Volume); }
    //void SetStrangenessCanonicalVolumeRadius(double Radius) { SetCanonicalVolumeRadius(Radius); }
 
    // Same as FixParameters but with a more clear name on what is actually does
    void ConstrainChemicalPotentials(bool resetInitialValues = true);
 
    
    virtual void FixParameters();

    virtual void FixParametersNoReset();

    virtual bool SolveChemicalPotentials(double totB = 0., double totQ = 0., double totS = 0., double totC = 0.,
      double muBinit = 0., double muQinit = 0., double muSinit = 0., double muCinit = 0.,
      bool ConstrMuB = true, bool ConstrMuQ = true, bool ConstrMuS = true, bool ConstrMuC = true);

    virtual bool FixChemicalPotentialsThroughDensities(double rhoB = 0., double rhoQ = 0., double rhoS = 0., double rhoC = 0.,
      double muBinit = 0., double muQinit = 0., double muSinit = 0., double muCinit = 0.,
      bool ConstrMuB = true, bool ConstrMuQ = true, bool ConstrMuS = true, bool ConstrMuC = true);

    virtual void CalculatePrimordialDensities() = 0;

    virtual void CalculateDensities();

    virtual void ValidateCalculation();

    std::string ValidityCheckLog() const { return m_ValidityLog; }

    virtual void CalculateDensitiesGCE() { CalculateDensities(); m_GCECalculated = true; }

    virtual void CalculateFeeddown();

    virtual void CalculateFluctuations();

    virtual void CalculateTemperatureDerivatives();

    virtual void CalculateTwoParticleCorrelations();

    virtual void CalculateTwoParticleFluctuationsDecays();

    virtual double TwoParticleSusceptibilityPrimordial(int i, int j) const;

    virtual double TwoParticleSusceptibilityPrimordialByPdg(long long id1, long long id2);

    virtual double TwoParticleSusceptibilityTemperatureDerivativePrimordial(int i, int j) const;

    virtual double TwoParticleSusceptibilityTemperatureDerivativePrimordialByPdg(long long id1, long long id2);

    virtual double NetParticleSusceptibilityPrimordialByPdg(long long id1, long long id2);

    virtual double TwoParticleSusceptibilityFinal(int i, int j) const;

    virtual double TwoParticleSusceptibilityFinalByPdg(long long id1, long long id2);

    virtual double NetParticleSusceptibilityFinalByPdg(long long id1, long long id2);

    virtual double PrimordialParticleChargeSusceptibility(int i, ConservedCharge::Name chg) const;

    virtual double PrimordialParticleChargeSusceptibilityByPdg(long long id1, ConservedCharge::Name chg);

    virtual double PrimordialNetParticleChargeSusceptibilityByPdg(long long id1, ConservedCharge::Name chg);
 

    virtual double FinalParticleChargeSusceptibility(int i, ConservedCharge::Name chg) const;

    virtual double FinalParticleChargeSusceptibilityByPdg(long long id1, ConservedCharge::Name chg);

    virtual double FinalNetParticleChargeSusceptibilityByPdg(long long id1, ConservedCharge::Name chg);
 

    virtual double SusceptibilityDimensionfull(ConservedCharge::Name i, ConservedCharge::Name j) const;

    virtual void CalculateSusceptibilityMatrix();

    virtual void CalculateProxySusceptibilityMatrix();

    virtual void CalculateParticleChargeCorrelationMatrix();

    virtual std::vector<double> CalculateChargeFluctuations(const std::vector<double> &chgs, int order = 4);

    virtual std::vector<double> CalculateGeneralizedSusceptibilities(const std::vector<std::vector<double>> &chgs);
 
 
    //virtual double GetParticlePrimordialDensity(unsigned int);
    //virtual double GetParticleTotalDensity(unsigned int);
 
    // Equation of state etc.

    double Pressure() { return CalculatePressure(); }

    double EnergyDensity() { return CalculateEnergyDensity(); }

    double EntropyDensity() { return CalculateEntropyDensity(); }

    double HadronDensity() { return CalculateHadronDensity(); }

    double BaryonDensity() { return CalculateBaryonDensity(); }

    double ElectricChargeDensity() { return CalculateChargeDensity(); }

    double StrangenessDensity() { return CalculateStrangenessDensity(); }

    double CharmDensity() { return CalculateCharmDensity(); }

    double AbsoluteBaryonDensity() { return CalculateAbsoluteBaryonDensity(); }

    double AbsoluteElectricChargeDensity() { return CalculateAbsoluteChargeDensity(); }

    double AbsoluteStrangenessDensity() { return CalculateAbsoluteStrangenessDensity(); }

    double AbsoluteCharmDensity() { return CalculateAbsoluteCharmDensity(); }

    double SpecificHeatChem() { return CalculateSpecificHeatChem(); }

    double cs2(bool rhoBconst = true, bool rhoQconst = true, bool rhoSconst = true, bool rhoCconst = true) { 
      return CalculateAdiabaticSpeedOfSoundSquared(rhoBconst, rhoQconst, rhoSconst, rhoCconst); 
    }

    double cT2(bool rhoBconst = true, bool rhoQconst = true, bool rhoSconst = true, bool rhoCconst = true) { 
      return CalculateIsothermalSpeedOfSoundSquared(rhoBconst, rhoQconst, rhoSconst, rhoCconst); 
    }

    double HeatCapacity(bool rhoBconst = true, bool rhoQconst = true, bool rhoSconst = true, bool rhoCconst = true) { 
      return CalculateHeatCapacity(rhoBconst, rhoQconst, rhoSconst, rhoCconst); 
    }


    virtual double CalculatePressure() = 0;
    virtual double CalculateEnergyDensity() = 0;
    virtual double CalculateEntropyDensity() = 0;
    virtual double CalculateHadronDensity();
    virtual double CalculateBaryonDensity();
    virtual double CalculateChargeDensity();
    virtual double CalculateStrangenessDensity();
    virtual double CalculateCharmDensity();
    virtual double CalculateAbsoluteBaryonDensity();
    virtual double CalculateAbsoluteChargeDensity();
    virtual double CalculateAbsoluteStrangenessDensity();
    virtual double CalculateAbsoluteCharmDensity();
    virtual double CalculateAbsoluteStrangenessDensityModulo();
    virtual double CalculateAbsoluteCharmDensityModulo();
    virtual double CalculateEnergyDensityDerivativeT() = 0;
    virtual double CalculateSpecificHeatChem();
    virtual double CalculateAdiabaticSpeedOfSoundSquared(bool rhoBconst = true, bool rhoQconst = true, bool rhoSconst = true, bool rhoCconst = true);
    virtual double CalculateIsothermalSpeedOfSoundSquared(bool rhoBconst = true, bool rhoQconst = true, bool rhoSconst = true, bool rhoCconst = true);
    virtual double CalculateHeatCapacity(bool rhoBconst = true, bool rhoQconst = true, bool rhoSconst = true, bool rhoCconst = true);

    virtual double CalculateArbitraryChargeDensity();
    
    virtual double CalculateBaryonMatterEntropyDensity() { return 0.; }

    virtual double CalculateMesonMatterEntropyDensity() { return 0.; }

    virtual double ParticleScaledVariance(int) { return 1.; }

    virtual double ParticleSkewness(int) { return 1.; }

    virtual double ParticleKurtosis(int) { return 1.; }

    virtual double CalculateEigenvolumeFraction() { return 0.; }

    virtual double ParticleScalarDensity(int i) = 0;
 
    virtual double GetMaxDiff() const { return m_MaxDiff; }

    virtual bool   IsLastSolutionOK() const { return m_LastCalculationSuccessFlag; }

    double GetdndT(int i) const;

    double GetDensity(long long PDGID, int feeddown) { return GetDensity(PDGID, static_cast<Feeddown::Type>(feeddown)); }
    
    double GetDensity(long long PDGID, Feeddown::Type feeddown);

    double GetYield(long long PDGID, Feeddown::Type feeddown) { return GetDensity(PDGID, feeddown) * Volume(); }
 
    std::vector<double> GetIdealGasDensities() const;

    ThermalParticleSystem* TPS() { return m_TPS; }

    const std::vector<double>& Densities()      const { return m_densities; }
    std::vector<double>& Densities() { return m_densities; }

    const std::vector<double>& TotalDensities() const { return m_densitiestotal; }

    const std::vector< std::vector<double> >& AllDensities()  const { return m_densitiesbyfeeddown; }

    void SetTAG(const std::string & tag) { m_TAG = tag; }

    const std::string& TAG() const { return m_TAG; }

    int  PdgToId(long long pdgid) const { return m_TPS->PdgToId(pdgid); }

    void ResetCalculatedFlags();

    bool IsCalculated() const { return m_Calculated; }

    bool IsFluctuationsCalculated() const { return m_FluctuationsCalculated; }

    bool IsSusceptibilitiesCalculated() const { return m_SusceptibilitiesCalculated; }

    bool IsTemperatureDerivativesCalculated() const { return m_TemperatureDerivativesCalculated; }

    bool IsGCECalculated() const { return m_GCECalculated; }

    double ScaledVariancePrimordial(int id) const { return (id >= 0 && id < static_cast<int>(m_wprim.size())) ? m_wprim[id] : 1.; }
    
    double ScaledVarianceTotal(int id) const { return (id >= 0 && id < static_cast<int>(m_wtot.size())) ? m_wtot[id] : 1.; }
    
    double SkewnessPrimordial(int id) const { return (id >= 0 && id < static_cast<int>(m_skewprim.size())) ? m_skewprim[id] : 1.; }
    
    double SkewnessTotal(int id) const { return (id >= 0 && id < static_cast<int>(m_skewtot.size())) ? m_skewtot[id] : 1.; }
    
    double KurtosisPrimordial(int id) const { return (id >= 0 && id < static_cast<int>(m_kurtprim.size())) ? m_kurtprim[id] : 1.; }
    
    double KurtosisTotal(int id) const { return (id >= 0 && id < static_cast<int>(m_kurttot.size())) ? m_kurttot[id] : 1.; }

    double ConservedChargeDensity(ConservedCharge::Name chg);

    double ConservedChargeDensitydT(ConservedCharge::Name chg);

    double Susc(ConservedCharge::Name i, ConservedCharge::Name j) const;// { return m_Susc[i][j]; }

    double dSuscdT(ConservedCharge::Name i, ConservedCharge::Name j) const;

    double ProxySusc(ConservedCharge::Name i, ConservedCharge::Name j) const;// { return m_ProxySusc[i][j]; }

    double ChargedMultiplicity(int type = 0);

    double ChargedScaledVariance(int type = 0);

    double ChargedMultiplicityFinal(int type = 0);

    double ChargedScaledVarianceFinal(int type = 0);

    ThermalModelEnsemble Ensemble() { return m_Ensemble; }
 

    virtual bool IsConservedChargeCanonical(ConservedCharge::Name charge) const { return 0; }

    ThermalModelInteraction InteractionModel() { return m_InteractionModel; }
 

    void SetDensityModelForParticleSpecies(int i, GeneralizedDensity* density_model = NULL);

    void SetDensityModelForParticleSpeciesByPdg(long long PDGID, GeneralizedDensity* density_model = NULL);

    void ClearDensityModels();
 
    const IdealGasFunctions::IdealGasFunctionsExtraConfig& GetIdealGasFunctionsExtraConfig() const { return m_IGFExtraConfig; }
 

     void SetMagneticField(double B = 0.0, int lmax = -1);
    void RecomputeThresholdsDueToMagneticField();

    std::vector<double> PartialPressures();

    void ClearMagneticField();
 
  protected:
    ThermalModelParameters m_Parameters;
    ThermalParticleSystem* m_TPS;

    IdealGasFunctions::IdealGasFunctionsExtraConfig m_IGFExtraConfig;
 
    bool   m_LastCalculationSuccessFlag;
    double m_MaxDiff;
 
    bool m_Calculated;
    bool m_FeeddownCalculated;
    bool m_FluctuationsCalculated;
    bool m_SusceptibilitiesCalculated;
    bool m_GCECalculated;
    bool m_UseWidth;
    bool m_NormBratio;
    bool m_QuantumStats;
    double m_QBgoal;
    double m_SBgoal;
    double m_Volume;
 
    bool m_ConstrainMuB;
    bool m_ConstrainMuQ;
    bool m_ConstrainMuS;
    bool m_ConstrainMuC;
 
    bool m_PCE;
 
    bool m_useOpenMP;
 
    std::vector<double> m_densities;
    std::vector<double> m_densitiestotal;
    //std::vector<double> m_densitiestotalweak;
    std::vector< std::vector<double> > m_densitiesbyfeeddown;
    std::vector<double> m_Chem;
 
    // Scaled variance
    std::vector<double> m_wprim;
    std::vector<double> m_wtot;
 
    // Skewness
    std::vector<double> m_skewprim;
    std::vector<double> m_skewtot;
 
    // Kurtosis
    std::vector<double> m_kurtprim;
    std::vector<double> m_kurttot;
 
    // 2nd order correlations of primordial and total numbers
    std::vector< std::vector<double> > m_PrimCorrel;
    std::vector< std::vector<double> > m_TotalCorrel;
 
    // Shifted chemical potential derivatives
    std::vector< std::vector<double> > m_dmusdmu;
 
 
    // Particle number-conserved charge correlators
    std::vector< std::vector<double> > m_PrimChargesCorrel;
    std::vector< std::vector<double> > m_FinalChargesCorrel;
 
    // Conserved charges susceptibility matrix
    std::vector< std::vector<double> > m_Susc;
 
    // Temperature derivative of the conserved charges susceptibility matrix (GeV^-1)
    std::vector< std::vector<double> > m_dSuscdT;
 
    // Susceptibility matrix of net-p, net-Q, and net-K
    std::vector< std::vector<double> > m_ProxySusc;
 
    // Temperature derivatives
    bool m_TemperatureDerivativesCalculated;
    // Densities
    std::vector<double> m_dndT;
    // (Shifted) chemical potentials
    std::vector<double> m_dmusdT;
    // Primordial mixed susceptibilities
    std::vector< std::vector<double> > m_PrimChi2sdT; // GeV^-1
 
    // Cumulants of arbitrary charge calculation
    //std::vector< std::vector<double> > m_chi;
 
    // Contains log of possible errors when checking the calculation
    std::string m_ValidityLog;
 
    double m_wnSum;
 
    std::string m_TAG;
 
    ThermalModelEnsemble m_Ensemble;
    ThermalModelInteraction m_InteractionModel;

    virtual double MuShift(int /*id*/) const { return 0.; }
    
 
  private:
    void ResetChemicalPotentials();
 
    double GetDensity(long long PDGID, const std::vector<double> *dens);
 
    class BroydenEquationsChem : public BroydenEquations
    {
    public:
      BroydenEquationsChem(ThermalModelBase *model) : BroydenEquations(), m_THM(model) { m_N = 2; }
      std::vector<double> Equations(const std::vector<double> &x);
    private:
      ThermalModelBase *m_THM;
    };
 
    class BroydenJacobianChem : public BroydenJacobian
    {
    public:
      BroydenJacobianChem(ThermalModelBase *model) : BroydenJacobian(), m_THM(model) { }
      std::vector<double> Jacobian(const std::vector<double> &x);
    private:
      ThermalModelBase *m_THM;
    };
 
    class BroydenChem : public Broyden
    {
    public:
      BroydenChem(ThermalModelBase *THM, BroydenEquations *eqs = NULL, BroydenJacobian *jaco = NULL) : Broyden(eqs, jaco) { m_THM = THM; }
      ~BroydenChem(void) { }
      std::vector<double> Solve(const std::vector<double> &x0, BroydenSolutionCriterium *solcrit = NULL, int max_iterations = MAX_ITERS);
    private:
      ThermalModelBase *m_THM;
    };
 
 
    class BroydenEquationsChemTotals : public BroydenEquations
    {
    public:
      BroydenEquationsChemTotals(const std::vector<int> & vConstr, const std::vector<int> & vType, const std::vector<double> & vTotals, ThermalModelBase *model);// : BroydenEquations(), m_THM(model) { m_N = 3; }
      std::vector<double> Equations(const std::vector<double> &x);
    private:
      std::vector<int> m_Constr;
      std::vector<int> m_Type;
      std::vector<double> m_Totals;
      ThermalModelBase *m_THM;
    };
 
    class BroydenJacobianChemTotals : public BroydenJacobian
    {
    public:
      BroydenJacobianChemTotals(const std::vector<int> & vConstr, const std::vector<int> & vType, const std::vector<double> & vTotals, ThermalModelBase *model) : BroydenJacobian(), m_Constr(vConstr), m_Type(vType), m_Totals(vTotals), m_THM(model) { }
      std::vector<double> Jacobian(const std::vector<double> &x);
    private:
      std::vector<int> m_Constr;
      std::vector<int> m_Type;
      std::vector<double> m_Totals;
      ThermalModelBase *m_THM;
    };
  };
 
} // namespace thermalfist
 
#endif // THERMALMODELBASE_H