RapidityFitBlastWave.cxx

Go to the documentation of this file.

#include "RapidityFitBlastWave.h"
 
#include "Minuit2/FCNBase.h"
#include "Minuit2/FunctionMinimum.h"
#include "Minuit2/MnHesse.h"
#include "Minuit2/MnMigrad.h"
#include "Minuit2/MnMinos.h"
#include "Minuit2/MnPrint.h"
#include "Minuit2/MnUserParameterState.h"
#include "Minuit2/SimplexMinimizer.h"
#include "TMath.h"
 
#include <fstream>
 
 
using namespace ROOT::Minuit2;
 
 
namespace RapidityFitBlastWaveNamespace {
  double dndy(double y, double A, double eta, double T, double mass) {
    double ret  = 0.;
    int iters   = 100;
    double deta = 2. * eta / iters;
    for (int i = 0; i < iters; ++i) {
      double et = -eta + (0.5 + i) * deta;
      ret += (mass * mass / T / T + mass / T * 2. / TMath::CosH(y - et)
              + 2. / TMath::CosH(y - et) / TMath::CosH(y - et))
             * TMath::Exp(-mass / T * TMath::CosH(y - et));
    }
    return ret * deta * A * T * T * T;
  }
 
  class dndyFCN : public FCNBase {
 
  public:
    dndyFCN(const std::vector<double>& y_,
            const std::vector<double>& dndy_,
            const std::vector<double>& err_,
            double T,
            double mass)
      : fPositions(y_)
      , fMeasurements(dndy_)
      , fErrors(err_)
      , fT(T)
      , fMass(mass)
    //    fPositions(y_), fMeasurements(dndy_), fErrors(err_), fT(T), fMass(mass), iter(0)
    {}
 
    ~dndyFCN() {}
 
    double operator()(const std::vector<double>& par) const {
      double chi2 = 0.;
      for (unsigned int i = 0; i < fPositions.size(); ++i) {
        double val = dndy(fPositions[i], par[0], par[1], fT, fMass);
        chi2 += (fMeasurements[i] - val) * (fMeasurements[i] - val) / fErrors[i]
                / fErrors[i];
      }
      return chi2;
    }
 
    double Up() const { return 1.; }
 
  private:
    std::vector<double> fPositions;
    std::vector<double> fMeasurements;
    std::vector<double> fErrors;
    double fT, fMass;
    //    int iter;
  };
}  // namespace RapidityFitBlastWaveNamespace
 
using namespace RapidityFitBlastWaveNamespace;
 
RapidityFitBlastWave::~RapidityFitBlastWave(void) {}
 
RapidityFitBlastWaveParameters RapidityFitBlastWave::PerformFit() {
  dndyFCN mfunc(y, dndy, err, fT, fMass);
  std::vector<double> params(2, 0.);
  params[0] = Parameters.A.value;
  params[1] = Parameters.eta.value;
  //params[2] = Parameters.yav.value;
  MnUserParameters upar;
  upar.Add("A",
           Parameters.A.value,
           Parameters.A.error,
           Parameters.A.xmin,
           Parameters.A.xmax);
  upar.Add("eta",
           Parameters.eta.value,
           Parameters.eta.error,
           Parameters.eta.xmin,
           Parameters.eta.xmax);
  //upar.Add("yav", Parameters.yav.value, Parameters.yav.error, Parameters.yav.xmin, Parameters.yav.xmax);
  int nparams = 2;
 
  MnMigrad migrad(mfunc, upar);
  //std::cout<<"start migrad "<<std::endl;
  FunctionMinimum min = migrad();
  MnHesse hess;
  hess(mfunc, min);
  RapidityFitBlastWaveParameters ret = Parameters;
  ret.A.value                        = (min.UserParameters()).Params()[0];
  ret.A.error                        = (min.UserParameters()).Errors()[0];
  ret.eta.value                      = (min.UserParameters()).Params()[1];
  ret.eta.error                      = (min.UserParameters()).Errors()[1];
  //ret.yav.value = (min.UserParameters()).Params()[2];
  //ret.yav.error = (min.UserParameters()).Errors()[2];
 
  params[0] = ret.A.value;
  params[1] = ret.eta.value;
  //params[2] = ret.yav.value;
  std::cout << "A = " << params[0] << " "
            << " eta = " << params[1] << " "
            << " chi2/ndf = " << mfunc(params) / (y.size() - nparams)
            << std::endl;
 
  ret.chi2ndf = mfunc(params) / (y.size() - nparams);
 
  return ret;
}