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 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 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 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 CalculateTwoParticleCorrelations();

    virtual void CalculateTwoParticleFluctuationsDecays();

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

    virtual double TwoParticleSusceptibilityPrimordialByPdg(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 void CalculateSusceptibilityMatrix();

    virtual void CalculateProxySusceptibilityMatrix();

    virtual void CalculateParticleChargeCorrelationMatrix();

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

    //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(); }

    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 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 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; }

    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; }

    void ResetCalculatedFlags();

    bool IsCalculated() const { return m_Calculated; }

    bool IsFluctuationsCalculated() const { return m_FluctuationsCalculated; }

    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 Susc(ConservedCharge::Name i, ConservedCharge::Name j) const { return m_Susc[i][j]; }

    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; }

  protected:
    ThermalModelParameters m_Parameters;
    ThermalParticleSystem* m_TPS;

    bool   m_LastCalculationSuccessFlag;
    double m_MaxDiff;

    bool m_Calculated;
    bool m_FeeddownCalculated;
    bool m_FluctuationsCalculated;
    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;

    // 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;

    // Susceptibility matrix of net-p, net-Q, and net-K
    std::vector< std::vector<double> > m_ProxySusc;

    // 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