ThermalModelPCE.h

Go to the documentation of this file.

#ifndef THERMALMODELPCE_H
#define THERMALMODELPCE_H
#include <map>

#include "HRGBase/ThermalModelBase.h"
#include "HRGBase/Broyden.h"

namespace thermalfist {

  class ThermalModelPCE
  {
  public:
    enum PCEMode {
      AtFixedTemperature = 0, 
      AtFixedVolume = 1
    };


    ThermalModelPCE(ThermalModelBase *THMbase, bool FreezeLonglived = false, double LonglivedResoWidthCut = 0.015);

    virtual ~ThermalModelPCE(void) { }


    void UseSahaForNuclei(bool flag) { m_UseSahaForNuclei = true; m_StabilityFlagsSet = false; }
    bool UseSahaForNuclei() const { return m_UseSahaForNuclei; }


    void FreezeLonglivedResonances(bool flag) { m_FreezeLonglivedResonances = true; m_StabilityFlagsSet = false; }
    bool FreezeLonglivedResonances() const { return m_FreezeLonglivedResonances; }


    void SetLonglivedResonanceWidthCut(double width_cut) { m_ResoWidthCut = width_cut; FreezeLonglivedResonances(true); }
    double LonglivedResonanceWidthCut() const { return m_ResoWidthCut; }


    virtual void SetStabilityFlags(const std::vector<int>& StabilityFlags);
    const std::vector<int>& StabilityFlags() const { return m_StabilityFlags; }


    void SetChemicalFreezeout(const ThermalModelParameters& params, const std::vector<double>& ChemInit = std::vector<double>(0));


    void SetEntropyDensityChem(double sinit) { m_EntropyDensityInit = sinit; }
    double EntropyDensityChem() const { return m_EntropyDensityInit; }

    ThermalModelBase* ThermalModel() const { return m_model; }

    virtual void CalculatePCE(double param, PCEMode mode = AtFixedTemperature);

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

    double Volume() const { return m_ParametersCurrent.V; }

    static void PrepareNucleiForPCE(ThermalParticleSystem *TPS);

    static std::vector<int> ComputePCEStabilityFlags(
      const ThermalParticleSystem* TPS,
      bool SahaEquationForNuclei = true,
      bool FreezeLongLived = false,
      double WidthCut = 0.015);

    void ApplyFixForBoseCondensation();
    

  protected:
    ThermalModelBase *m_model;

    bool m_UseSahaForNuclei;

    bool m_FreezeLonglivedResonances;
    double m_ResoWidthCut;

    bool m_ChemicalFreezeoutSet;

    bool m_IsCalculated;

    bool m_StabilityFlagsSet;
    std::vector<int> m_StabilityFlags;
    int m_StableComponentsNumber;
    std::vector< std::vector<double> > m_EffectiveCharges;
    std::vector<int> m_StableMapTo;

    ThermalModelParameters m_ParametersInit;
    std::vector<double> m_ChemInit;
    std::vector<double> m_DensitiesInit;
    std::vector<double> m_StableDensitiesInit;
    double m_EntropyDensityInit;
    double m_ParticleDensityInit;

    ThermalModelParameters m_ParametersCurrent;
    std::vector<double> m_ChemCurrent;

    class BroydenEquationsPCE : public BroydenEquations
    {
    public:
      BroydenEquationsPCE(ThermalModelPCE *model, int mode = 0) : BroydenEquations(), m_THM(model), m_Mode(mode) { m_N = m_THM->m_StableComponentsNumber + 1; }
      std::vector<double> Equations(const std::vector<double> &x);
    private:
      ThermalModelPCE *m_THM;
      int m_Mode;
    };

  };

} // namespace thermalfist

#endif