ThermalModelFit.h

Go to the documentation of this file.

#ifndef THERMALMODELFIT_H
#define THERMALMODELFIT_H
 
#include "ThermalModelBase.h"
#include <cmath>
#include <string>
 
struct FitParameter {
  bool toFit;
  double value;
  double error;
  double xmin, xmax;
  std::string name;
  FitParameter(std::string name_ = "",
               bool toFit_       = true,
               double val        = 1.,
               double err        = 0.,
               double vmin       = -2.,
               double vmax       = 2.)
    : toFit(toFit_)
    , value(val)
    , error(err)
    , xmin(vmin)
    , xmax(vmax)
    , name(name_) {}
};
 
struct ThermalModelFitParametersExtended {
  FitParameter T, muB, muS, muQ, gammaS, V;
  double chi2ndf;
  FitParameter nB, rhoB, rhoQ, en, entropy, pressure, eta;
  ThermalModelFitParametersExtended()
    : T()
    , muB()
    , muS()
    , muQ()
    , gammaS()
    , V()
    , chi2ndf()
    , nB()
    , rhoB()
    , rhoQ()
    , en()
    , entropy()
    , pressure()
    , eta() {}
  ThermalModelFitParametersExtended(ThermalModelBase* model)
    : T()
    , muB()
    , muS()
    , muQ()
    , gammaS()
    , V()
    , chi2ndf()
    , nB()
    , rhoB()
    , rhoQ()
    , en()
    , entropy()
    , pressure()
    , eta() {
    T.value        = model->Parameters.T;
    muB.value      = model->Parameters.muB;
    muS.value      = model->Parameters.muS;
    muQ.value      = model->Parameters.muQ;
    gammaS.value   = model->Parameters.gammaS;
    V.value        = model->Parameters.V;
    nB.value       = model->CalculateHadronDensity();
    rhoB.value     = model->CalculateBaryonDensity();
    rhoQ.value     = model->CalculateChargeDensity();
    en.value       = model->CalculateEnergyDensity();
    entropy.value  = model->CalculateEntropyDensity();
    pressure.value = model->CalculatePressure();
    eta.value =
      model->CalculateShearViscosity() / model->CalculateEntropyDensity();
  }
};
 
struct ThermalModelFitParameters {
  FitParameter T, muB, muS, muQ, gammaS, V, R;
  double chi2ndf;
  ThermalModelFitParameters(double T_      = 0.12,
                            double muB_    = 0.5,
                            double muS_    = 0.1,
                            double muQ_    = -0.01,
                            double gammaS_ = 1.,
                            double V_      = 4000.,
                            double R_      = 1.)
    : T(FitParameter("T", true, T_, 0.05, 0.05, 0.18))
    , muB(FitParameter("muB", true, muB_, 0.05, 0., 1.))
    , muS(FitParameter("muS", false, muS_))
    , muQ(FitParameter("muQ", false, muQ_))
    , gammaS(FitParameter("gammaS", false, gammaS_, 0.2, 0.5, 1.5))
    , V(FitParameter("V", true, V_, 2000., 1., 20000.))
    , R(FitParameter("R", false, R_, 0., 0., 100.))
    , chi2ndf(1.) {
    // T = FitParameter("T", true, T_, 0.05, 0.05, 0.18);
    // muB = FitParameter("muB", true, muB_, 0.05, 0., 1.);
    // muS = FitParameter("muS", false, muS_);
    // muQ = FitParameter("muQ", false, muQ_);
    // gammaS = FitParameter("gammaS", false, gammaS_, 0.2, 0.5, 1.5);
    // V = FitParameter("V", true, V_, 2000., 1., 20000.);
    // R = FitParameter("R", false, R_, 0., 0., 100.);
  }
  FitParameter GetParamater(const std::string& name) const {
    if (T.name == name) return T;
    if (muB.name == name) return muB;
    if (muS.name == name) return muS;
    if (muQ.name == name) return muQ;
    if (gammaS.name == name) return gammaS;
    if (V.name == name) return V;
    if (R.name == name) return R;
    return FitParameter();
  }
  void SetParameter(const std::string& name, const FitParameter& param) {
    if (T.name == name) T = param;
    if (muB.name == name) muB = param;
    if (muS.name == name) muS = param;
    if (muQ.name == name) muQ = param;
    if (gammaS.name == name) gammaS = param;
    if (V.name == name) V = param;
    if (R.name == name) R = param;
  }
  void SetParameter(const std::string& name,
                    double val,
                    double err,
                    double xmin,
                    double xmax) {
    if (T.name == name) T = FitParameter(name, true, val, err, xmin, xmax);
    if (muB.name == name) muB = FitParameter(name, true, val, err, xmin, xmax);
    if (muS.name == name) muS = FitParameter(name, true, val, err, xmin, xmax);
    if (muQ.name == name) muQ = FitParameter(name, true, val, err, xmin, xmax);
    if (gammaS.name == name)
      gammaS = FitParameter(name, true, val, err, xmin, xmax);
    if (V.name == name) V = FitParameter(name, true, val, err, xmin, xmax);
    if (R.name == name) R = FitParameter(name, true, val, err, xmin, xmax);
  }
  void SetParameterFitFlag(const std::string& name, bool toFit) {
    if (T.name == name) T.toFit = toFit;
    if (muB.name == name) muB.toFit = toFit;
    if (muS.name == name) muS.toFit = toFit;
    if (muQ.name == name) muQ.toFit = toFit;
    if (gammaS.name == name) gammaS.toFit = toFit;
    if (V.name == name) V.toFit = toFit;
    if (R.name == name) R.toFit = toFit;
  }
  ThermalModelParameters GetThermalModelParameters() {
    return ThermalModelParameters(
      T.value, muB.value, muS.value, muQ.value, gammaS.value, V.value, R.value);
  }
};
 
struct ExperimentMultiplicity {
  int fPDGID;
  double fValue, fError;
  ExperimentMultiplicity(int PDGID    = -211,
                         double value = 300.,
                         double error = 20.)
    : fPDGID(PDGID), fValue(value), fError(error) {}
  void addSystematicError(double percent) {
    fError = sqrt(fError * fError + percent * percent * fValue * fValue);
  }
};
 
struct ExperimentRatio {
  int PDGID1, PDGID2;
  double value, error;
  ExperimentRatio(int PDGID1_   = 211,
                  int PDGID2_   = -211,
                  double value_ = 1.,
                  double error_ = 0.1)
    : PDGID1(PDGID1_), PDGID2(PDGID2_), value(value_), error(error_) {}
  ExperimentRatio(int PDGID1_,
                  int PDGID2_,
                  double value1,
                  double error1,
                  double value2,
                  double error2)
    : PDGID1(PDGID1_)
    , PDGID2(PDGID2_)
    , value(value1 / value2)
    , error(sqrt(error1 * error1 / value2 / value2
                 + value1 * value1 / value2 / value2 / value2 / value2 * error2
                     * error2)) {
    // value = value1 / value2;
    // error = sqrt(error1*error1 / value2 / value2 + value1 * value1 / value2  / value2  / value2  / value2  * error2 * error2);
    //error += value * 0.05;
  }
  void addSystematicError(double percent) {
    error = sqrt(error * error + percent * percent * value * value);
  }
};
 
struct FittedQuantity {
  int type;  // 0 - Multiplicity, 1 - Ratio
  ExperimentMultiplicity mult;
  ExperimentRatio ratio;
  FittedQuantity()
    : type(0), mult(ExperimentMultiplicity(-211, 10., 1.)), ratio() {
    // type = 0;
    // mult = ExperimentMultiplicity(-211, 10., 1.);
  }
  FittedQuantity(const ExperimentMultiplicity& op)
    : type(0), mult(op), ratio() {
    // type = 0;
    // mult = op;
  }
  FittedQuantity(const ExperimentRatio& op) : type(1), mult(), ratio(op) {
    // type = 1;
    // ratio = op;
  }
};
 
class ThermalModelBase;
 
class ThermalModelFit {
public:
  ThermalModelBase* model;
  ThermalModelFitParameters Parameters;
  ThermalModelFitParametersExtended ExtendedParameters;
  double QBgoal;
  std::vector<ExperimentMultiplicity> Multiplicities;
  std::vector<ExperimentRatio> Ratios;
  ThermalModelFit(ThermalModelBase* model_,
                  double T      = 0.12,
                  double muB    = 0.5,
                  double muS    = 0.,
                  double muQ    = 0.,
                  double gammaS = 1.,
                  double V      = 4000.,
                  double R      = 1.)
    : model(model_)
    , Parameters(T, muB, muS, muQ, gammaS, V, R)
    , ExtendedParameters()
    , QBgoal(0.4)
    , Multiplicities()
    , Ratios() {
    Multiplicities.resize(0);
    Ratios.resize(0);
    model->SetParameters(Parameters.GetThermalModelParameters());
  }
  ~ThermalModelFit(void);
  void SetFitFlag(const std::string& name, bool flag) {
    Parameters.SetParameterFitFlag(name, flag);
    QBgoal = 0.4;
  }
  void SetQBConstraint(double QB) { QBgoal = QB; }
  void SetRatios(const std::vector<ExperimentRatio>& Ratios_) {
    Ratios = Ratios_;
  }
  void AddRatios(const std::vector<ExperimentRatio>& Ratios_) {
    Ratios.insert(Ratios.end(), Ratios_.begin(), Ratios_.end());
  }
  void AddRatio(const ExperimentRatio& Ratio_) { Ratios.push_back(Ratio_); }
  void ClearRatios() { Ratios.resize(0); }
  void PrintRatios();
 
  void SetMultiplicities(
    const std::vector<ExperimentMultiplicity>& Multiplicities_) {
    Multiplicities = Multiplicities_;
  }
  void AddMultiplicities(
    const std::vector<ExperimentMultiplicity>& Multiplicities_) {
    Multiplicities.insert(
      Multiplicities.end(), Multiplicities_.begin(), Multiplicities_.end());
  }
  void AddMultiplicity(const ExperimentMultiplicity& Multiplicity_) {
    Multiplicities.push_back(Multiplicity_);
  }
  void ClearMultiplicities() { Multiplicities.resize(0); }
  void PrintMultiplicities();
  double chi2Ndf(double T, double muB);
 
  ThermalModelFitParameters PerformFit();
 
  ThermalModelFit(const ThermalModelFit&);
  ThermalModelFit& operator=(const ThermalModelFit&);
};
 
std::vector<FittedQuantity> loadExpDataFromFile(const std::string& filename);
 
#endif