CbmRoot
ThermalModelFit.cxx
Go to the documentation of this file.
1 #include "ThermalModelFit.h"
2 #include "ThermalModelBase.h"
3 
4 #include "Minuit2/FCNBase.h"
5 #include "Minuit2/FunctionMinimum.h"
6 #include "Minuit2/MnHesse.h"
7 #include "Minuit2/MnMigrad.h"
8 #include "Minuit2/MnMinos.h"
9 #include "Minuit2/MnPrint.h"
10 #include "Minuit2/MnUserParameterState.h"
11 #include "Minuit2/SimplexMinimizer.h"
12 
13 //#include "DebugText.h"
14 #include <fstream>
15 
16 //#include <QDebug>
17 
18 
19 using namespace ROOT::Minuit2;
20 
21 int migraditer = 0;
22 
23 namespace ThermalModelFitNamespace {
24  class RatiosFCN : public FCNBase {
25 
26  public:
27  RatiosFCN(ThermalModelFit* thmfit_) : THMFit(thmfit_), iter(0) {}
28 
29  ~RatiosFCN() {}
30 
31  double operator()(const std::vector<double>& par) const {
32  //std::cout << ++iter << "\n";
33  ++migraditer;
34  double chi2 = 0.;
35  THMFit->model->Parameters.T = par[0];
36  THMFit->model->Parameters.muB = par[1];
37  THMFit->model->Parameters.gammaS = par[2];
38  THMFit->model->Parameters.V = par[3];
39  THMFit->model->Parameters.muQ = -par[1] / 50.;
40  THMFit->model->Parameters.muS = par[1] / 5.;
41  THMFit->model->SetQBgoal(THMFit->QBgoal);
42  THMFit->model->FixParameters();
43  THMFit->Parameters.muQ.value = THMFit->model->Parameters.muQ;
44  THMFit->Parameters.muS.value = THMFit->model->Parameters.muS;
45  THMFit->model->CalculateDensities();
46  for (unsigned int i = 0; i < THMFit->Ratios.size(); ++i) {
47  int id1 = THMFit->model->TPS->PDGtoID[THMFit->Ratios[i].PDGID1];
48  int id2 = THMFit->model->TPS->PDGtoID[THMFit->Ratios[i].PDGID2];
49  double ModelRatio = THMFit->model->densitiestotal[id1]
50  / THMFit->model->densitiestotal[id2];
51  chi2 += (ModelRatio - THMFit->Ratios[i].value)
52  * (ModelRatio - THMFit->Ratios[i].value)
53  / THMFit->Ratios[i].error / THMFit->Ratios[i].error;
54  }
55  for (unsigned int i = 0; i < THMFit->Multiplicities.size(); ++i) {
56  int id1 = THMFit->model->TPS->PDGtoID[THMFit->Multiplicities[i].fPDGID];
57  double ModelMult =
58  THMFit->model->densitiestotal[id1] * THMFit->model->Parameters.V;
59  chi2 += (ModelMult - THMFit->Multiplicities[i].fValue)
60  * (ModelMult - THMFit->Multiplicities[i].fValue)
61  / THMFit->Multiplicities[i].fError
62  / THMFit->Multiplicities[i].fError;
63  }
64  return chi2;
65  }
66 
67  double Up() const { return 1.; }
68 
69  RatiosFCN(const RatiosFCN&);
70  RatiosFCN& operator=(const RatiosFCN&);
71 
72  private:
73  ThermalModelFit* THMFit;
74  int iter;
75  };
76 } // namespace ThermalModelFitNamespace
77 
78 using namespace ThermalModelFitNamespace;
79 
80 ThermalModelFit::~ThermalModelFit(void) {}
81 
82 ThermalModelFitParameters ThermalModelFit::PerformFit() {
83  RatiosFCN mfunc(this);
84  std::vector<double> params(4, 0.);
85  params[0] = Parameters.T.value;
86  params[1] = Parameters.muB.value;
87  params[2] = Parameters.gammaS.value;
88  params[3] = Parameters.V.value;
89  MnUserParameters upar;
90  upar.Add("T",
91  Parameters.T.value,
92  Parameters.T.error,
93  Parameters.T.xmin,
94  Parameters.T.xmax);
95  upar.Add("muB",
96  Parameters.muB.value,
97  Parameters.muB.error,
98  Parameters.muB.xmin,
99  Parameters.muB.xmax);
100  upar.Add("gammaS",
101  Parameters.gammaS.value,
102  Parameters.gammaS.error,
103  Parameters.gammaS.xmin,
104  Parameters.gammaS.xmax);
105  upar.Add("V",
106  Parameters.V.value,
107  Parameters.V.error,
108  Parameters.V.xmin,
109  Parameters.V.xmax);
110  int nparams = 4;
111  if (!Parameters.gammaS.toFit) {
112  upar.Fix("gammaS");
113  nparams--;
114  }
115  if (!Parameters.V.toFit || Multiplicities.size() == 0) {
116  upar.Fix("V");
117  nparams--;
118  }
119 
120  bool repeat = 0;
121  if (model->fUseWidth) {
122  repeat = 1;
123  model->SetUseWidth(false);
124  }
125 
126  MnMigrad migrad(mfunc, upar);
127  FunctionMinimum min = migrad();
128  MnHesse hess;
129  hess(mfunc, min);
130  ThermalModelFitParameters ret = Parameters;
131  ret.T.value = (min.UserParameters()).Params()[0];
132  ret.T.error = (min.UserParameters()).Errors()[0];
133  ret.muB.value = (min.UserParameters()).Params()[1];
134  ret.muB.error = (min.UserParameters()).Errors()[1];
135  ret.gammaS.value = (min.UserParameters()).Params()[2];
136  ret.gammaS.error = (min.UserParameters()).Errors()[2];
137  ret.V.value = (min.UserParameters()).Params()[3];
138  ret.V.error = (min.UserParameters()).Errors()[3];
139 
140  if (repeat) {
141  model->SetUseWidth(true);
142 
143  upar.SetValue("T", (min.UserParameters()).Params()[0]);
144  upar.SetValue("muB", (min.UserParameters()).Params()[1]);
145  upar.SetValue("gammaS", (min.UserParameters()).Params()[2]);
146  upar.SetValue("V", (min.UserParameters()).Params()[3]);
147 
148  MnMigrad migradd(mfunc, upar);
149  min = migradd();
150 
151  ret = Parameters;
152  ret.T.value = (min.UserParameters()).Params()[0];
153  ret.T.error = (min.UserParameters()).Errors()[0];
154  ret.muB.value = (min.UserParameters()).Params()[1];
155  ret.muB.error = (min.UserParameters()).Errors()[1];
156  ret.gammaS.value = (min.UserParameters()).Params()[2];
157  ret.gammaS.error = (min.UserParameters()).Errors()[2];
158  ret.V.value = (min.UserParameters()).Params()[3];
159  ret.V.error = (min.UserParameters()).Errors()[3];
160  }
161 
162  model->SetParameters(ret.GetThermalModelParameters());
163  model->SetQBgoal(QBgoal);
164  model->FixParameters();
165  ret.muQ.value = model->Parameters.muQ;
166  ret.muS.value = model->Parameters.muS;
167  Parameters = ret;
168  params[0] = ret.T.value;
169  params[1] = ret.muB.value;
170  params[2] = ret.gammaS.value;
171  params[3] = ret.V.value;
172  std::cout << "T= " << params[0] * 1000. << " MeV "
173  << "muB= " << params[1] * 1000. << " MeV "
174  << " gammaS = " << params[2] << " V = " << params[3] << " fm^3"
175  << " chi2/ndf = "
176  << mfunc(params) / (Multiplicities.size() + Ratios.size() - nparams)
177  << std::endl;
178 
179  ret.chi2ndf =
180  mfunc(params) / (Multiplicities.size() + Ratios.size() - nparams);
181 
182  ExtendedParameters = ThermalModelFitParametersExtended(model);
183 
184  ExtendedParameters.T.error = ret.T.error;
185  ExtendedParameters.muB.error = ret.muB.error;
186  ExtendedParameters.gammaS.error = ret.gammaS.error;
187  ExtendedParameters.V.error = ret.V.error;
188 
189  double Tinit = model->Parameters.T;
190  double muBinit = model->Parameters.muB;
191 
192  double dT = 0.0005, dmu = 0.0005;
193  model->Parameters.T = Tinit + dT;
194  model->Parameters.muB = muBinit;
195  model->fCalculated = false;
196  model->FixParameters();
197  model->CalculateDensities();
198  ThermalModelFitParametersExtended ParametersdT =
199  ThermalModelFitParametersExtended(model);
200  model->Parameters.T = Tinit;
201  model->Parameters.muB = muBinit + dmu;
202  model->fCalculated = false;
203  model->FixParameters();
204  model->CalculateDensities();
205  ThermalModelFitParametersExtended ParametersdMu =
206  ThermalModelFitParametersExtended(model);
207  model->Parameters.T = Tinit;
208  model->Parameters.muB = muBinit;
209  model->FixParameters();
210  model->CalculateDensities();
211 
212 
213  double err1 = (min.UserCovariance()).operator()(0, 0);
214  double err2 = (min.UserCovariance()).operator()(1, 1);
215  double cov = (min.UserCovariance()).operator()(0, 1);
216  double deriv1 = (ParametersdT.muQ.value - ExtendedParameters.muQ.value) / dT;
217  double deriv2 =
218  (ParametersdMu.muQ.value - ExtendedParameters.muQ.value) / dmu;
219  ExtendedParameters.muQ.error =
220  sqrt(deriv1 * deriv1 * err1 + deriv2 * deriv2 * err2
221  + 2. * deriv1 * deriv2 * cov);
222  deriv1 = (ParametersdT.muS.value - ExtendedParameters.muS.value) / dT;
223  deriv2 = (ParametersdMu.muS.value - ExtendedParameters.muS.value) / dmu;
224  ExtendedParameters.muS.error =
225  sqrt(deriv1 * deriv1 * err1 + deriv2 * deriv2 * err2
226  + 2. * deriv1 * deriv2 * cov);
227  deriv1 = (ParametersdT.nB.value - ExtendedParameters.nB.value) / dT;
228  deriv2 = (ParametersdMu.nB.value - ExtendedParameters.nB.value) / dmu;
229  ExtendedParameters.nB.error =
230  sqrt(deriv1 * deriv1 * err1 + deriv2 * deriv2 * err2
231  + 2. * deriv1 * deriv2 * cov);
232  deriv1 = (ParametersdT.rhoB.value - ExtendedParameters.rhoB.value) / dT;
233  deriv2 = (ParametersdMu.rhoB.value - ExtendedParameters.rhoB.value) / dmu;
234  ExtendedParameters.rhoB.error =
235  sqrt(deriv1 * deriv1 * err1 + deriv2 * deriv2 * err2
236  + 2. * deriv1 * deriv2 * cov);
237  deriv1 = (ParametersdT.rhoQ.value - ExtendedParameters.rhoQ.value) / dT;
238  deriv2 = (ParametersdMu.rhoQ.value - ExtendedParameters.rhoQ.value) / dmu;
239  ExtendedParameters.rhoQ.error =
240  sqrt(deriv1 * deriv1 * err1 + deriv2 * deriv2 * err2
241  + 2. * deriv1 * deriv2 * cov);
242  deriv1 = (ParametersdT.en.value - ExtendedParameters.en.value) / dT;
243  deriv2 = (ParametersdMu.en.value - ExtendedParameters.en.value) / dmu;
244  ExtendedParameters.en.error =
245  sqrt(deriv1 * deriv1 * err1 + deriv2 * deriv2 * err2
246  + 2. * deriv1 * deriv2 * cov);
247  deriv1 = (ParametersdT.entropy.value - ExtendedParameters.entropy.value) / dT;
248  deriv2 =
249  (ParametersdMu.entropy.value - ExtendedParameters.entropy.value) / dmu;
250  ExtendedParameters.entropy.error =
251  sqrt(deriv1 * deriv1 * err1 + deriv2 * deriv2 * err2
252  + 2. * deriv1 * deriv2 * cov);
253  deriv1 =
254  (ParametersdT.pressure.value - ExtendedParameters.pressure.value) / dT;
255  deriv2 =
256  (ParametersdMu.pressure.value - ExtendedParameters.pressure.value) / dmu;
257  ExtendedParameters.pressure.error =
258  sqrt(deriv1 * deriv1 * err1 + deriv2 * deriv2 * err2
259  + 2. * deriv1 * deriv2 * cov);
260  deriv1 = (ParametersdT.eta.value - ExtendedParameters.eta.value) / dT;
261  deriv2 = (ParametersdMu.eta.value - ExtendedParameters.eta.value) / dmu;
262  ExtendedParameters.eta.error =
263  sqrt(deriv1 * deriv1 * err1 + deriv2 * deriv2 * err2
264  + 2. * deriv1 * deriv2 * cov);
265 
266  return ret;
267 }
268 
269 void ThermalModelFit::PrintRatios() {
270  model->SetParameters(Parameters.GetThermalModelParameters());
271  model->SetQBgoal(QBgoal);
272  model->FixParameters();
273  model->CalculateDensities();
274  for (unsigned int i = 0; i < Ratios.size(); ++i) {
275  int ind1 = model->TPS->PDGtoID[Ratios[i].PDGID1];
276  int ind2 = model->TPS->PDGtoID[Ratios[i].PDGID2];
277  std::cout << model->TPS->fParticles[ind1].fName << "/"
278  << model->TPS->fParticles[ind2].fName << " = "
279  << model->densitiestotal[ind1] / model->densitiestotal[ind2]
280  << " " << Ratios[i].value << " " << Ratios[i].error << "\n";
281  }
282 }
283 
284 void ThermalModelFit::PrintMultiplicities() {
285  model->SetParameters(Parameters.GetThermalModelParameters());
286  model->SetQBgoal(QBgoal);
287  model->FixParameters();
288  model->CalculateDensities();
289  for (unsigned int i = 0; i < Multiplicities.size(); ++i) {
290  int ind1 = model->TPS->PDGtoID[Multiplicities[i].fPDGID];
291  std::cout << "<" << model->TPS->fParticles[ind1].fName
292  << "> = " << model->densitiestotal[ind1] * model->Parameters.V
293  << " " << Multiplicities[i].fValue << " "
294  << Multiplicities[i].fError << "\n";
295  }
296 }
297 
298 using namespace std;
299 
300 std::vector<FittedQuantity> loadExpDataFromFile(const std::string& filename) {
301  std::vector<FittedQuantity> ret(0);
302  ifstream fin(filename.c_str());
303  //fin.open(filename);
304  if (fin.is_open()) {
305  string tmp;
306  char tmpc[500];
307  fin.getline(tmpc, 500);
308  tmp = string(tmpc);
309  while (1) {
310  vector<string> fields = split(tmp, ';');
311  if (fields.size() < 4) break;
312  int type = atoi(fields[0].c_str());
313  int pdgid1 = atoi(fields[1].c_str()), pdgid2 = 0;
314  double value, error;
315  if (type) {
316  pdgid2 = atoi(fields[2].c_str());
317  value = atof(fields[3].c_str());
318  error = 1.;
319  if (fields.size() >= 4) error = atof(fields[4].c_str());
320  } else {
321  value = atof(fields[2].c_str());
322  error = atof(fields[3].c_str());
323  }
324 
325  if (type)
326  ret.push_back(
327  FittedQuantity(ExperimentRatio(pdgid1, pdgid2, value, error)));
328  else
329  ret.push_back(
330  FittedQuantity(ExperimentMultiplicity(pdgid1, value, error)));
331 
332  fin.getline(tmpc, 500);
333  tmp = string(tmpc);
334  }
335  fin.close();
336  }
337  return ret;
338 }
339 
340 double ThermalModelFit::chi2Ndf(double T, double muB) {
341  double Tinit = T, muBinit = muB;
342  RatiosFCN mfunc(this);
343  int nparams = 4;
344  if (!Parameters.gammaS.toFit) nparams--;
345  if (!Parameters.V.toFit || Multiplicities.size() == 0) nparams--;
346  std::vector<double> params(4, 0.);
347  params[0] = T;
348  params[1] = muB;
349  params[2] = model->Parameters.gammaS;
350  params[3] = model->Parameters.V;
351  model->Parameters.muS = muB / 5.;
352  model->Parameters.muQ = -muB / 50.;
353  double ret = mfunc(params);
354  model->Parameters.T = Tinit;
355  model->Parameters.muB = muBinit;
356  return ret / (Multiplicities.size() + Ratios.size() - nparams);
357 }
FitParameter::value
double value
Definition: ThermalModelFit.h:10
ThermalModelFit::PerformFit
ThermalModelFitParameters PerformFit()
Definition: ThermalModelFit.cxx:82
split
std::vector< std::string > & split(const std::string &s, char delim, std::vector< std::string > &elems)
Definition: ThermalParticleSystem.cxx:144
ThermalModelFit::PrintMultiplicities
void PrintMultiplicities()
Definition: ThermalModelFit.cxx:284
migraditer
int migraditer
Definition: ThermalModelFit.cxx:21
ThermalModelFitParameters::GetThermalModelParameters
ThermalModelParameters GetThermalModelParameters()
Definition: ThermalModelFit.h:147
ThermalModelFit::~ThermalModelFit
~ThermalModelFit(void)
Definition: ThermalModelFit.cxx:80
ThermalModelFitNamespace::RatiosFCN::RatiosFCN
RatiosFCN(const RatiosFCN &)
FitParameter::error
double error
Definition: ThermalModelFit.h:11
ThermalModelFitNamespace::RatiosFCN::operator()
double operator()(const std::vector< double > &par) const
Definition: ThermalModelFit.cxx:31
ThermalModelFitParameters::muQ
FitParameter muQ
Definition: ThermalModelFit.h:80
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
ThermalModelFitNamespace::RatiosFCN::THMFit
ThermalModelFit * THMFit
Definition: ThermalModelFit.cxx:73
ThermalModelFit::chi2Ndf
double chi2Ndf(double T, double muB)
Definition: ThermalModelFit.cxx:340
ThermalModelBase.h
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
ThermalModelFitNamespace::RatiosFCN::~RatiosFCN
~RatiosFCN()
Definition: ThermalModelFit.cxx:29
ThermalModelFitParameters::muB
FitParameter muB
Definition: ThermalModelFit.h:80
FittedQuantity
Definition: ThermalModelFit.h:194
ThermalModelFitParameters::chi2ndf
double chi2ndf
Definition: ThermalModelFit.h:81
ThermalModelFitParameters::muS
FitParameter muS
Definition: ThermalModelFit.h:80
ExperimentMultiplicity
Definition: ThermalModelFit.h:153
ThermalModelFitParametersExtended::nB
FitParameter nB
Definition: ThermalModelFit.h:31
ThermalModelFitParameters::V
FitParameter V
Definition: ThermalModelFit.h:80
ThermalModelFitParametersExtended::rhoB
FitParameter rhoB
Definition: ThermalModelFit.h:31
loadExpDataFromFile
std::vector< FittedQuantity > loadExpDataFromFile(const std::string &filename)
Definition: ThermalModelFit.cxx:300
ThermalModelFitParametersExtended::eta
FitParameter eta
Definition: ThermalModelFit.h:31
ThermalModelFit
Definition: ThermalModelFit.h:216
min
friend F32vec4 min(const F32vec4 &a, const F32vec4 &b)
Definition: L1/vectors/P4_F32vec4.h:33
ThermalModelFitParametersExtended::rhoQ
FitParameter rhoQ
Definition: ThermalModelFit.h:31
ThermalModelFitNamespace::RatiosFCN::RatiosFCN
RatiosFCN(ThermalModelFit *thmfit_)
Definition: ThermalModelFit.cxx:27
ThermalModelFitParameters
Definition: ThermalModelFit.h:79
ThermalModelFitNamespace::RatiosFCN::Up
double Up() const
Definition: ThermalModelFit.cxx:67
ThermalModelFitParametersExtended::en
FitParameter en
Definition: ThermalModelFit.h:31
ThermalModelFitParametersExtended::muS
FitParameter muS
Definition: ThermalModelFit.h:29
ThermalModelFitNamespace::RatiosFCN::operator=
RatiosFCN & operator=(const RatiosFCN &)
ExperimentRatio
Definition: ThermalModelFit.h:165
ThermalModelFitNamespace::RatiosFCN::iter
int iter
Definition: ThermalModelFit.cxx:74
ThermalModelFit::PrintRatios
void PrintRatios()
Definition: ThermalModelFit.cxx:269
ThermalModelFitParametersExtended
Definition: ThermalModelFit.h:28
ThermalModelFitParametersExtended::entropy
FitParameter entropy
Definition: ThermalModelFit.h:31
ThermalModelFitParametersExtended::muQ
FitParameter muQ
Definition: ThermalModelFit.h:29
ThermalModelFitParameters::gammaS
FitParameter gammaS
Definition: ThermalModelFit.h:80
ThermalModelFitNamespace
Definition: ThermalModelFit.cxx:23
ThermalModelFitParametersExtended::pressure
FitParameter pressure
Definition: ThermalModelFit.h:31
ThermalModelFitParameters::T
FitParameter T
Definition: ThermalModelFit.h:80
ThermalModelFit.h
ThermalModelFitNamespace::RatiosFCN
Definition: ThermalModelFit.cxx:24