Thermal-FIST 1.5
Package for hadron resonance gas model applications
Loading...
Searching...
No Matches
ThermalModelFit.cpp
Go to the documentation of this file.
1/*
2 * Thermal-FIST package
3 *
4 * Copyright (c) 2014-2019 Volodymyr Vovchenko
5 *
6 * GNU General Public License (GPLv3 or later)
7 */
8#include "HRGFit/ThermalModelFit.h"
9
10#include <fstream>
11#include <cstdio>
12#include <cstdlib>
13#include <ctime>
14#include <iostream>
15#include <iomanip>
16#include <sstream>
17
18#ifdef USE_MINUIT
19#include "Minuit2/FCNBase.h"
20#include "Minuit2/FunctionMinimum.h"
21#include "Minuit2/MnMigrad.h"
22#include "Minuit2/MnMinos.h"
23#include "Minuit2/MnHesse.h"
24#include "Minuit2/MnUserParameterState.h"
25//#include "Minuit2/MnPrint.h"
26//#include "Minuit2/SimplexMinimizer.h"
27#endif
28
29#include "HRGBase/ThermalModelBase.h"
30#include "HRGBase/Utility.h"
31#include "HRGPCE/ThermalModelPCEAnnihilation.h"
32
33namespace thermalfist {
34
35 #ifdef USE_MINUIT
36
37 using namespace ROOT::Minuit2;
38
39 namespace {
40
41 class FitFCN : public FCNBase {
42
43 public:
44
45 FitFCN(ThermalModelFit *thmfit_, bool verbose_ = true) : m_THMFit(thmfit_), m_verbose(verbose_) {
46 }
47
48 ~FitFCN() {}
49
50 double operator()(const std::vector<double>& par) const {
51 m_THMFit->Increment();
52 double chi2 = 0.;
53 if (par[2]<0.) return 1e12;
54 if (par[3]<0.) return 1e12;
55
56 m_THMFit->SetParameterValue("T", par[0]);
57 m_THMFit->SetParameterValue("muB", par[1]);
58 m_THMFit->SetParameterValue("gammaS", par[2]);
59 m_THMFit->SetParameterValue("R", par[3]);
60 m_THMFit->SetParameterValue("Rc", par[4]);
61 m_THMFit->SetParameterValue("gammaq", par[5]);
62 m_THMFit->SetParameterValue("muQ", par[6]);
63 m_THMFit->SetParameterValue("muS", par[7]);
64 m_THMFit->SetParameterValue("muC", par[8]);
65 m_THMFit->SetParameterValue("gammac", par[9]);
66 m_THMFit->SetParameterValue("Tkin", par[10]);
67
68 m_THMFit->model()->SetTemperature(par[0]);
69
70 if (!m_THMFit->model()->ConstrainMuB())
71 m_THMFit->model()->SetBaryonChemicalPotential(par[1]);
72
73 m_THMFit->model()->SetGammaS(par[2]);
74
75 m_THMFit->model()->SetVolumeRadius(par[3]);
76
77 if (m_THMFit->FixVcOverV())
78 m_THMFit->model()->SetCanonicalVolume(m_THMFit->model()->Volume() * m_THMFit->VcOverV());
79 else
80 m_THMFit->model()->SetCanonicalVolumeRadius(par[4]);
81
82
83 m_THMFit->model()->SetGammaq(par[5]);
84
85 m_THMFit->model()->SetGammaC(par[9]);
86
87 if (m_THMFit->model()->ConstrainMuQ())
88 m_THMFit->model()->SetElectricChemicalPotential(-par[1] / 50.);
89 else
90 m_THMFit->model()->SetElectricChemicalPotential(par[6]);
91
92 if (m_THMFit->model()->ConstrainMuS())
93 m_THMFit->model()->SetStrangenessChemicalPotential(par[1] / 5.);
94 else
95 m_THMFit->model()->SetStrangenessChemicalPotential(par[7]);
96
97 if (m_THMFit->model()->ConstrainMuC())
98 m_THMFit->model()->SetCharmChemicalPotential(par[1] / 5.);
99 else
100 m_THMFit->model()->SetCharmChemicalPotential(par[8]);
101
102 m_THMFit->model()->SetParameters(m_THMFit->model()->Parameters());
103
104 //m_THMFit->model()->SetQoverB(m_THMFit->QoverB());
105
106 IdealGasFunctions::calculationHadBECIssue = false;
107
108 m_THMFit->model()->ConstrainChemicalPotentials();
109
110 if (m_THMFit->UseTkin()) {
111 m_THMFit->modelpce()->SetChemicalFreezeout(m_THMFit->model()->Parameters(), m_THMFit->model()->ChemicalPotentials());
112 m_THMFit->modelpce()->CalculatePCE(par[10]);
113 }
114 else {
115 m_THMFit->model()->CalculateDensities();
116 }
117
118 // If current chemical potentials lead to
119 // Bose-Einstein function divergence (\mu > m),
120 // then effectively discard parameter of the current iteration by setting chi^2 to 10^12
121 if (IdealGasFunctions::calculationHadBECIssue) {
122 std::cerr << m_THMFit->Iters();
123 std::cerr << "**WARNING** Issue with Bose-Einstein condensation, discarding this iteration..." << std::endl;
124 return m_THMFit->Chi2() = chi2 = 1.e12;
125 }
126
127 // Ratios first
128 for (size_t i = 0; i < m_THMFit->FittedQuantities().size(); ++i) {
129 if (m_THMFit->FittedQuantities()[i].type == FittedQuantity::Ratio) {
130 const ExperimentRatio &ratio = m_THMFit->FittedQuantities()[i].ratio;
131 double dens1 = m_THMFit->model()->GetDensity(ratio.fPDGID1, ratio.fFeedDown1);
132 double dens2 = m_THMFit->model()->GetDensity(ratio.fPDGID2, ratio.fFeedDown2);
133 double ModelRatio = dens1 / dens2;
134 m_THMFit->ModelData(i) = ModelRatio;
135 if (m_THMFit->FittedQuantities()[i].toFit)
136 chi2 += (ModelRatio - ratio.fValue) * (ModelRatio - ratio.fValue) / ratio.fError / ratio.fError;
137 }
138 }
139
140 // Yields second
141 for (size_t i = 0; i < m_THMFit->FittedQuantities().size(); ++i) {
142 if (m_THMFit->FittedQuantities()[i].type == FittedQuantity::Multiplicity) {
143 const ExperimentMultiplicity &multiplicity = m_THMFit->FittedQuantities()[i].mult;
144 double dens = m_THMFit->model()->GetDensity(multiplicity.fPDGID, multiplicity.fFeedDown);
145 double ModelMult = dens * m_THMFit->model()->Parameters().V;
146 m_THMFit->ModelData(i) = ModelMult;
147 if (m_THMFit->FittedQuantities()[i].toFit)
148 chi2 += (ModelMult - multiplicity.fValue) * (ModelMult - multiplicity.fValue) / multiplicity.fError / multiplicity.fError;
149 }
150 }
151
152 if (m_verbose) {
153 printf("%15d ", m_THMFit->Iters());
154 printf("%15lf ", chi2);
155 if (m_THMFit->Parameters().T.toFit)
156 printf("%15lf ", par[0]);
157 if (m_THMFit->Parameters().muB.toFit)
158 printf("%15lf ", m_THMFit->model()->Parameters().muB);
159 if (m_THMFit->Parameters().muQ.toFit)
160 printf("%15lf ", m_THMFit->model()->Parameters().muQ);
161 if (m_THMFit->Parameters().muS.toFit)
162 printf("%15lf ", m_THMFit->model()->Parameters().muS);
163 if (m_THMFit->Parameters().muC.toFit)
164 printf("%15lf ", m_THMFit->model()->Parameters().muC);
165 if (m_THMFit->Parameters().R.toFit)
166 printf("%15lf ", par[3]);
167 if (m_THMFit->Parameters().Rc.toFit)
168 printf("%15lf ", par[4]);
169 if (m_THMFit->Parameters().gammaq.toFit)
170 printf("%15lf ", m_THMFit->model()->Parameters().gammaq);
171 if (m_THMFit->Parameters().gammaS.toFit)
172 printf("%15lf ", m_THMFit->model()->Parameters().gammaS);
173 if (m_THMFit->Parameters().gammaC.toFit)
174 printf("%15lf ", m_THMFit->model()->Parameters().gammaC);
175 if (m_THMFit->Parameters().Tkin.toFit)
176 printf("%15lf ", par[10]);
177 printf("\n");
178
179 if (m_THMFit->model()->Ensemble() == ThermalModelBase::CE)
180 printf("B = %10.5lf\tQ = %10.5lf\tS = %10.5lf\tC = %10.5lf\n",
181 m_THMFit->model()->CalculateBaryonDensity() * m_THMFit->model()->Parameters().V,
182 m_THMFit->model()->CalculateChargeDensity() * m_THMFit->model()->Parameters().V,
183 m_THMFit->model()->CalculateStrangenessDensity() * m_THMFit->model()->Parameters().V,
184 m_THMFit->model()->CalculateCharmDensity() * m_THMFit->model()->Parameters().V);
185 }
186
187 m_THMFit->Chi2() = chi2;
188 if (m_THMFit->model()->Ensemble() == ThermalModelBase::CE) {
189 m_THMFit->BT() = m_THMFit->model()->CalculateBaryonDensity() * m_THMFit->model()->Parameters().V;
190 m_THMFit->QT() = m_THMFit->model()->CalculateChargeDensity() * m_THMFit->model()->Parameters().V;
191 m_THMFit->ST() = m_THMFit->model()->CalculateStrangenessDensity() * m_THMFit->model()->Parameters().V;
192 m_THMFit->CT() = m_THMFit->model()->CalculateCharmDensity() * m_THMFit->model()->Parameters().V;
193 }
194
195 if (chi2!=chi2) {
196 chi2 = 1.e12;
197 printf("**WARNING** chi2 evaluated to NaN\n");
198 }
199
200 return chi2;
201 }
202
203 double Up() const {return 1.;}
204
205 private:
206 ThermalModelFit *m_THMFit;
207 int m_iter;
208 bool m_verbose;
209 };
210 }
211
212 #endif
213
214 ThermalModelFit::ThermalModelFit(ThermalModelBase *model_):
215 m_model(model_), m_modelpce(NULL), m_Parameters(model_->Parameters()), m_FixVcToV(true), m_VcOverV(1.),
216 m_YieldsAtTkin(false), m_SahaForNuclei(true), m_PCEFreezeLongLived(false), m_PCEWidthCut(0.015),
217 m_PCEAnnihilation(false), m_PCEPionAnnihilationNumber(5.)
218 {
219 }
220
221
225
226 ThermalModelFitParameters ThermalModelFit::PerformFit(bool verbose, bool AsymmErrors) {
227 #ifdef USE_MINUIT
228 m_ModelData.resize(m_Quantities.size(), 0.);
229
230 m_Parameters.B = m_model->Parameters().B;
231 m_Parameters.Q = m_model->Parameters().Q;
232 m_Parameters.S = m_model->Parameters().S;
233 m_Parameters.C = m_model->Parameters().C;
234
235 if (UseTkin()) {
236 if (m_PCEAnnihilation) {
237 m_modelpce = new ThermalModelPCEAnnihilation(m_model, m_PCEFreezeLongLived, m_PCEWidthCut);
238 static_cast<ThermalModelPCEAnnihilation*>(m_modelpce)->SetPionAnnihilationNumber(m_PCEPionAnnihilationNumber);
239 m_modelpce->UseSahaForNuclei(m_SahaForNuclei);
240 m_modelpce->SetStabilityFlags(static_cast<ThermalModelPCEAnnihilation*>(m_modelpce)->RecalculateStabilityFlags());
241 }
242 else {
243 m_modelpce = new ThermalModelPCE(m_model, m_PCEFreezeLongLived, m_PCEWidthCut);
244 if (!m_SahaForNuclei)
245 m_modelpce->SetStabilityFlags(m_modelpce->ComputePCEStabilityFlags(m_model->TPS(), m_SahaForNuclei, m_PCEFreezeLongLived, m_PCEWidthCut));
246 }
247 }
248
249 m_Iters = 0;
250 FitFCN mfunc(this, verbose);
251 std::vector<double> params(11, 0.);
252 params[0] = m_Parameters.T.value;
253 params[1] = m_Parameters.muB.value;
254 params[2] = m_Parameters.gammaS.value;
255 params[3] = m_Parameters.R.value;
256 params[4] = m_Parameters.Rc.value;
257 params[5] = m_Parameters.gammaq.value;
258 params[6] = m_Parameters.muQ.value;
259 params[7] = m_Parameters.muS.value;
260 params[8] = m_Parameters.muC.value;
261 params[9] = m_Parameters.gammaC.value;
262 params[10] = m_Parameters.Tkin.value;
263
264 MnUserParameters upar;
265 upar.Add("T", m_Parameters.T.value, m_Parameters.T.error, m_Parameters.T.xmin, m_Parameters.T.xmax);
266 upar.Add("muB", m_Parameters.muB.value, m_Parameters.muB.error, m_Parameters.muB.xmin, m_Parameters.muB.xmax);
267 upar.Add("gammaS", m_Parameters.gammaS.value, m_Parameters.gammaS.error, m_Parameters.gammaS.xmin, m_Parameters.gammaS.xmax);
268 upar.Add("R", m_Parameters.R.value, m_Parameters.R.error, m_Parameters.R.xmin, m_Parameters.R.xmax);
269 upar.Add("Rc", m_Parameters.Rc.value, m_Parameters.Rc.error, m_Parameters.Rc.xmin, m_Parameters.Rc.xmax);
270 upar.Add("gammaq", m_Parameters.gammaq.value, m_Parameters.gammaq.error, m_Parameters.gammaq.xmin, m_Parameters.gammaq.xmax);
271 upar.Add("muQ", m_Parameters.muQ.value, m_Parameters.muQ.error, m_Parameters.muQ.xmin, m_Parameters.muQ.xmax);
272 upar.Add("muS", m_Parameters.muS.value, m_Parameters.muS.error, m_Parameters.muS.xmin, m_Parameters.muS.xmax);
273 upar.Add("muC", m_Parameters.muC.value, m_Parameters.muC.error, m_Parameters.muC.xmin, m_Parameters.muC.xmax);
274 upar.Add("gammaC", m_Parameters.gammaC.value, m_Parameters.gammaC.error, m_Parameters.gammaC.xmin, m_Parameters.gammaC.xmax);
275 upar.Add("Tkin", m_Parameters.Tkin.value, m_Parameters.Tkin.error, m_Parameters.Tkin.xmin, m_Parameters.Tkin.xmax);
276
277 int nparams = 11;
278
279 // If only ratios fitted then volume drops out
280 //if (m_Multiplicities.size() == 0 && m_model->Ensemble() == ThermalModelBase::GCE)
281 if (m_Multiplicities.size() == 0)
282 m_Parameters.SetParameterFitFlag("R", false);
283
284 // If GCE, or Vc fixed to V, then correlation volume drops out
285 if (m_model->Ensemble() == ThermalModelBase::GCE || this->FixVcOverV())
286 m_Parameters.SetParameterFitFlag("Rc", false);
287
288 // If full CE, or S/B fixed, then muB drops out
289 if (m_model->IsConservedChargeCanonical(ConservedCharge::BaryonCharge)
290 || m_model->ConstrainMuB()
291 || !m_model->TPS()->hasBaryons())
292 m_Parameters.SetParameterFitFlag("muB", false);
293
294 // If full CE, or Q/B fixed, or no charged particles, then muQ drops out
295 if (m_model->IsConservedChargeCanonical(ConservedCharge::CharmCharge)
296 || m_model->ConstrainMuQ()
297 || !m_model->TPS()->hasCharged())
298 m_Parameters.SetParameterFitFlag("muQ", false);
299
300 // If full CE, SCE, or S fixed to zero, or no strange particles then muS drops out
301 if (m_model->IsConservedChargeCanonical(ConservedCharge::StrangenessCharge)
302 || m_model->ConstrainMuS()
303 || !m_model->TPS()->hasStrange())
304 m_Parameters.SetParameterFitFlag("muS", false);
305
306 // If not GCE, or C fixed to zero, or no charm particles then muC drops out
307 if (m_model->IsConservedChargeCanonical(ConservedCharge::CharmCharge)
308 || m_model->ConstrainMuC()
309 || !m_model->TPS()->hasCharmed())
310 m_Parameters.SetParameterFitFlag("muC", false);
311
312 // If no strangeness, then gammaS drops out (but beware of neutral particles with strange quarks!)
313 if (!m_model->TPS()->hasStrange())
314 m_Parameters.SetParameterFitFlag("gammaS", false);
315
316 // If no charm, then gammaC drops out (but beware of neutral particles with charm quarks!)
317 if (!m_model->TPS()->hasCharmed())
318 m_Parameters.SetParameterFitFlag("gammaC", false);
319
320 // If not using PCE, no point using Tkin
321 if (!UseTkin())
322 m_Parameters.SetParameterFitFlag("Tkin", false);
323
324 if (!m_Parameters.T.toFit) { upar.Fix("T"); nparams--; }
325 if (!m_Parameters.R.toFit) { upar.Fix("R"); nparams--; }
326 if (!m_Parameters.Rc.toFit) { upar.Fix("Rc"); nparams--; }
327 if (!m_Parameters.muB.toFit) { upar.Fix("muB"); nparams--; }
328 if (!m_Parameters.muQ.toFit) { upar.Fix("muQ"); nparams--; }
329 if (!m_Parameters.muS.toFit) { upar.Fix("muS"); nparams--; }
330 if (!m_Parameters.muC.toFit) { upar.Fix("muC"); nparams--; }
331 if (!m_Parameters.gammaq.toFit) { upar.Fix("gammaq"); nparams--; }
332 if (!m_Parameters.gammaS.toFit) { upar.Fix("gammaS"); nparams--; }
333 if (!m_Parameters.gammaC.toFit) { upar.Fix("gammaC"); nparams--; }
334 if (!m_Parameters.Tkin.toFit) { upar.Fix("Tkin"); nparams--; }
335
336
337 m_Ndf = GetNdf();
338
339 bool repeat = 0;
340
341 ThermalModelFitParameters ret;
342
343 if (nparams>0) {
344
345 if (verbose) {
346 printf("Starting a thermal fit...\n\n");
347 printf("%15s ", "Iteration");
348 printf("%15s ", "chi2");
349 if (m_Parameters.T.toFit)
350 printf("%15s ", "T [GeV]");
351 if (m_Parameters.muB.toFit)
352 printf("%15s ", "muB [GeV]");
353 if (m_Parameters.muQ.toFit)
354 printf("%15s ", "muQ [GeV]");
355 if (m_Parameters.muS.toFit)
356 printf("%15s ", "muS [GeV]");
357 if (m_Parameters.muC.toFit)
358 printf("%15s ", "muC [GeV]");
359 if (m_Parameters.R.toFit)
360 printf("%15s ", "R [fm]");
361 if (m_Parameters.Rc.toFit)
362 printf("%15s ", "Rc [fm]");
363 if (m_Parameters.gammaq.toFit)
364 printf("%15s ", "gammaq");
365 if (m_Parameters.gammaS.toFit)
366 printf("%15s ", "gammaS");
367 if (m_Parameters.gammaC.toFit)
368 printf("%15s ", "gammaC");
369 if (m_Parameters.Tkin.toFit)
370 printf("%15s ", "Tkin [GeV]");
371 printf("\n");
372 }
373
374 MnMigrad migrad(mfunc, upar);
375
376 FunctionMinimum min = migrad();
377
378 if (verbose)
379 printf("\nMinimum found! Now calculating the error matrix...\n\n");
380
381 MnHesse hess;
382 hess(mfunc, min);
383
384 ret = m_Parameters;
385
386 if (AsymmErrors) {
387 MnMinos mino(mfunc, min);
388 std::pair<double, double> errs;
389 if (m_Parameters.T.toFit) { errs = mino(0); ret.T.errm = abs(errs.first); ret.T.errp = abs(errs.second); }
390 if (m_Parameters.muB.toFit) { errs = mino(1); ret.muB.errm = abs(errs.first); ret.muB.errp = abs(errs.second); }
391 if (m_Parameters.gammaS.toFit) { errs = mino(2); ret.gammaS.errm = abs(errs.first); ret.gammaS.errp = abs(errs.second); }
392 if (m_Parameters.R.toFit) { errs = mino(3); ret.R.errm = abs(errs.first); ret.R.errp = abs(errs.second); }
393 if (m_Parameters.Rc.toFit) { errs = mino(4); ret.Rc.errm = abs(errs.first); ret.Rc.errp = abs(errs.second); }
394 if (m_Parameters.gammaq.toFit) { errs = mino(5); ret.gammaq.errm = abs(errs.first); ret.gammaq.errp = abs(errs.second); }
395 }
396
397 ret.T.value = (min.UserParameters()).Params()[0];
398 ret.T.error = (min.UserParameters()).Errors()[0];
399 ret.muB.value = (min.UserParameters()).Params()[1];
400 ret.muB.error = (min.UserParameters()).Errors()[1];
401 ret.gammaS.value = (min.UserParameters()).Params()[2];
402 ret.gammaS.error = (min.UserParameters()).Errors()[2];
403 ret.R.value = (min.UserParameters()).Params()[3];
404 ret.R.error = (min.UserParameters()).Errors()[3];
405 ret.Rc.value = (min.UserParameters()).Params()[4];
406 ret.Rc.error = (min.UserParameters()).Errors()[4];
407 ret.gammaq.value = (min.UserParameters()).Params()[5];
408 ret.gammaq.error = (min.UserParameters()).Errors()[5];
409 ret.muQ.value = (min.UserParameters()).Params()[6];
410 ret.muQ.error = (min.UserParameters()).Errors()[6];
411 ret.muS.value = (min.UserParameters()).Params()[7];
412 ret.muS.error = (min.UserParameters()).Errors()[7];
413 ret.muC.value = (min.UserParameters()).Params()[8];
414 ret.muC.error = (min.UserParameters()).Errors()[8];
415 ret.gammaC.value = (min.UserParameters()).Params()[9];
416 ret.gammaC.error = (min.UserParameters()).Errors()[9];
417 ret.Tkin.value = (min.UserParameters()).Params()[10];
418 ret.Tkin.error = (min.UserParameters()).Errors()[10];
419
420 if (!m_Parameters.Rc.toFit && FixVcOverV()) {
421 ret.Rc.value = ret.R.value * pow(VcOverV(), 1./3.);
422 ret.Rc.error = 0.;
423 }
424
425 //if (repeat) {
426 // m_model->SetUseWidth(true);
427
428 // upar.SetValue("T", (min.UserParameters()).Params()[0]);
429 // upar.SetValue("muB", (min.UserParameters()).Params()[1]);
430 // upar.SetValue("gammaS", (min.UserParameters()).Params()[2]);
431 // upar.SetValue("R", (min.UserParameters()).Params()[3]);
432 // upar.SetValue("Rc", (min.UserParameters()).Params()[4]);
433 // upar.SetValue("gammaq", (min.UserParameters()).Params()[5]);
434 // upar.SetValue("muQ", (min.UserParameters()).Params()[6]);
435 // upar.SetValue("muS", (min.UserParameters()).Params()[7]);
436 // upar.SetValue("muC", (min.UserParameters()).Params()[8]);
437 // upar.SetValue("gammaC", (min.UserParameters()).Params()[9]);
438 // upar.SetValue("Tkin", (min.UserParameters()).Params()[10]);
439
440 // MnMigrad migradd(mfunc, upar);
441 // min = migradd();
442
443 // ret = m_Parameters;
444 // ret.T.value = (min.UserParameters()).Params()[0];
445 // ret.T.error = (min.UserParameters()).Errors()[0];
446 // ret.muB.value = (min.UserParameters()).Params()[1];
447 // ret.muB.error = (min.UserParameters()).Errors()[1];
448 // ret.gammaS.value = (min.UserParameters()).Params()[2];
449 // ret.gammaS.error = (min.UserParameters()).Errors()[2];
450 // ret.R.value = (min.UserParameters()).Params()[3];
451 // ret.R.error = (min.UserParameters()).Errors()[3];
452 // ret.Rc.value = (min.UserParameters()).Params()[4];
453 // ret.Rc.error = (min.UserParameters()).Errors()[4];
454 // ret.gammaq.value = (min.UserParameters()).Params()[5];
455 // ret.gammaq.error = (min.UserParameters()).Errors()[5];
456 // ret.muQ.value = (min.UserParameters()).Params()[6];
457 // ret.muQ.error = (min.UserParameters()).Errors()[6];
458 // ret.muS.value = (min.UserParameters()).Params()[7];
459 // ret.muS.error = (min.UserParameters()).Errors()[7];
460 // ret.muC.value = (min.UserParameters()).Params()[8];
461 // ret.muC.error = (min.UserParameters()).Errors()[8];
462 // ret.gammaC.value = (min.UserParameters()).Params()[9];
463 // ret.gammaC.error = (min.UserParameters()).Errors()[9];
464 // ret.Tkin.value = (min.UserParameters()).Params()[10];
465 // ret.Tkin.error = (min.UserParameters()).Errors()[10];
466 //}
467 }
468 else {
469 ret = m_Parameters;
470
471 ret.T.value = upar.Params()[0];
472 ret.T.error = upar.Errors()[0];
473 ret.muB.value = upar.Params()[1];
474 ret.muB.error = upar.Errors()[1];
475 ret.gammaS.value = upar.Params()[2];
476 ret.gammaS.error = upar.Errors()[2];
477 ret.R.value = upar.Params()[3];
478 ret.R.error = upar.Errors()[3];
479 ret.Rc.value = upar.Params()[4];
480 ret.Rc.error = upar.Errors()[4];
481 ret.gammaq.value = upar.Params()[5];
482 ret.gammaq.error = upar.Errors()[5];
483 ret.muQ.value = upar.Params()[6];
484 ret.muQ.error = upar.Errors()[6];
485 ret.muS.value = upar.Params()[7];
486 ret.muS.error = upar.Errors()[7];
487 ret.muC.value = upar.Params()[8];
488 ret.muC.error = upar.Errors()[8];
489 ret.gammaC.value = upar.Params()[9];
490 ret.gammaC.error = upar.Errors()[9];
491 ret.Tkin.value = upar.Params()[10];
492 ret.Tkin.error = upar.Errors()[10];
493 }
494
495
496 ThermalModelParameters parames = ret.GetThermalModelParameters();
497
498 parames.B = m_model->Parameters().B;
499 parames.S = m_model->Parameters().S;
500 parames.Q = m_model->Parameters().Q;
501 parames.C = m_model->Parameters().C;
502 //m_model->SetParameters(parames);
503 //m_model->FixParameters();
504 if (!ret.muQ.toFit && m_model->ConstrainMuQ()) {
505 ret.muQ.value = m_model->Parameters().muQ;
506 ret.muQ.error = 0.;
507 }
508 if (!ret.muS.toFit && m_model->ConstrainMuS()) {
509 ret.muS.value = m_model->Parameters().muS;
510 ret.muS.error = 0.;
511 }
512 if (!ret.muC.toFit && m_model->ConstrainMuC()) {
513 ret.muC.value = m_model->Parameters().muC;
514 ret.muC.error = 0.;
515 }
516 m_Parameters = ret;
517
518 params[0] = ret.T.value;
519 params[1] = ret.muB.value;
520 params[2] = ret.gammaS.value;
521 params[3] = ret.R.value;
522 params[4] = ret.Rc.value;
523 params[5] = ret.gammaq.value;
524 params[6] = ret.muQ.value;
525 params[7] = ret.muS.value;
526 params[8] = ret.muC.value;
527 params[9] = ret.gammaC.value;
528 params[10] = ret.Tkin.value;
529
530 ret.chi2 = mfunc(params);
531 int ndf = 0;
532 for (size_t i = 0; i < m_Quantities.size(); ++i)
533 if (m_Quantities[i].toFit) ndf++;
534 ndf = ndf - nparams;
535
536 ret.chi2ndf = ret.chi2 / ndf;
537 ret.ndf = ndf;
538
539 if (!AsymmErrors) {
540 for (size_t i = 0; i < ret.ParameterList.size(); ++i) {
541 ret.GetParameter(i).errm = ret.GetParameter(i).errp = ret.GetParameter(i).error;
542 }
543 }
544
545
546 m_Parameters = ret;
547
548
549 m_ExtendedParameters = ThermalModelFitParametersExtended(m_model);
550
551 if (m_modelpce != NULL) {
552 delete m_modelpce;
553 m_modelpce = NULL;
554 }
555
556 if (verbose)
557 printf("Thermal fit finished\n\n");
558 return ret;
559 #else
560 throw std::runtime_error("Cannot fit without MINUIT2 library!");
561 #endif
562 }
563
564 void ThermalModelFit::PrintRatios() {
565 m_model->SetParameters(m_Parameters.GetThermalModelParameters());
566 m_model->FixParameters();
567 m_model->CalculateDensities();
568
569 for (size_t i = 0; i < m_Quantities.size(); ++i) {
570 if (m_Quantities[i].type == FittedQuantity::Ratio) {
571 const ExperimentRatio &ratio = m_Quantities[i].ratio;
572 int ind1 = m_model->TPS()->PdgToId(ratio.fPDGID1);
573 int ind2 = m_model->TPS()->PdgToId(ratio.fPDGID2);
574 double dens1 = m_model->GetDensity(ratio.fPDGID1, ratio.fFeedDown1);
575 double dens2 = m_model->GetDensity(ratio.fPDGID2, ratio.fFeedDown2);
576 std::cout << m_model->TPS()->Particles()[ind1].Name() << "/" << m_model->TPS()->Particles()[ind2].Name() << " = " <<
577 dens1 / dens2 << " " << ratio.fValue << " " << ratio.fError << "\n";
578 }
579 }
580 }
581
582 void ThermalModelFit::PrintParameters() {
583 printf("%20s\n", "Fit parameters");
584 for(int i=0;i<6;++i) {
585 FitParameter param = m_Parameters.GetParameter(i);
586 if (param.name == "Tkin" && !m_YieldsAtTkin)
587 continue;
588 if (param.toFit) {
589 printf("%10s = %10lf %2s %lf\n", param.name.c_str(), param.value, "+-", param.error);
590 }
591 }
592 if (m_model->Ensemble() == ThermalModelBase::CE) {
593 printf("%10s = %10lf\n", "B", m_model->CalculateBaryonDensity() * m_model->Volume());
594 printf("%10s = %10lf\n", "S", m_model->CalculateStrangenessDensity() * m_model->Volume());
595 printf("%10s = %10lf\n", "Q", m_model->CalculateChargeDensity() * m_model->Volume());
596 }
597 printf("\n");
598 }
599
600
601 void ThermalModelFit::PrintMultiplicities() {
602 // No longer needed, done at the end of PerformFit
603 //m_model->SetParameters(m_Parameters.GetThermalModelParameters());
604 //m_model->SetQoverB(m_QBgoal);
605 //m_model->FixParameters();
606 //m_model->CalculateDensities();
607
608 for (size_t i = 0; i < m_Quantities.size(); ++i) {
609 if (m_Quantities[i].type == FittedQuantity::Multiplicity) {
610 const ExperimentMultiplicity &mult = m_Quantities[i].mult;
611
612 double dens1 = m_model->GetDensity(mult.fPDGID, mult.fFeedDown);
613
614 std::string tname = m_model->TPS()->GetNameFromPDG(mult.fPDGID);
615
616 if (mult.fPDGID == 1)
617 printf("%10s\t%11s\t%6lf %2s %lf\n", "Npart", "Experiment:",
618 mult.fValue, "+-", mult.fError);
619 else if (mult.fPDGID == 33340)
620 printf("%10s\t%11s\t%6lf %2s %lf\n", (m_model->TPS()->ParticleByPDG(3334).Name() + " + " + m_model->TPS()->ParticleByPDG(-3334).Name()).c_str(), "Experiment:",
621 mult.fValue, "+-", mult.fError);
622 else
623 printf("%10s\t%11s\t%6lf %2s %lf\n", tname.c_str(), "Experiment:", mult.fValue, "+-", mult.fError);
624 printf("%10s\t%11s\t%6lf %2s %lf\n", "", "Model:", dens1 * m_model->Parameters().V, "+-", 0.);
625 printf("%10s\t%11s\t%6lf\n", "", "Deviation:", (dens1 * m_model->Parameters().V - mult.fValue) / mult.fError);
626 printf("\n");
627
628 }
629 }
630
631 }
632
633 void ThermalModelFit::PrintYields() {
634 //m_model->SetParameters(m_Parameters.GetThermalModelParameters());
635 //m_model->SetQoverB(m_QBgoal);
636 //m_model->FixParameters();
637 //m_model->CalculateDensities();
638
639 for (size_t i = 0; i < m_Quantities.size(); ++i) {
640 if (m_Quantities[i].type == FittedQuantity::Multiplicity) {
641 const ExperimentMultiplicity &mult = m_Quantities[i].mult;
642
643 double dens1 = m_model->GetDensity(mult.fPDGID, mult.fFeedDown);
644
645 std::string tname = m_model->TPS()->GetNameFromPDG(mult.fPDGID);
646
647 if (mult.fPDGID == 1)
648 printf("%10s\t%11s\t%6lf %2s %lf\n", "Npart", "Experiment:",
649 mult.fValue, "+-", mult.fError);
650 else if (mult.fPDGID == 33340)
651 printf("%10s\t%11s\t%6lf %2s %lf\n", (m_model->TPS()->ParticleByPDG(3334).Name() + " + " + m_model->TPS()->ParticleByPDG(-3334).Name()).c_str(), "Experiment:",
652 mult.fValue, "+-", mult.fError);
653 else
654 printf("%10s\t%11s\t%6lf %2s %lf\n", tname.c_str(), "Experiment:", mult.fValue, "+-", mult.fError);
655 printf("%10s\t%11s\t%6lf %2s %lf\n", "", "Model:", dens1 * m_model->Parameters().V, "+-", 0.);
656 printf("%10s\t%11s\t%6lf\n", "", "Deviation:", (dens1 * m_model->Parameters().V - mult.fValue) / mult.fError);
657 printf("\n");
658
659 }
660 }
661
662 for (size_t i = 0; i < m_Quantities.size(); ++i) {
663 if (m_Quantities[i].type == FittedQuantity::Ratio) {
664 const ExperimentRatio &ratio = m_Quantities[i].ratio;
665 double dens1 = m_model->GetDensity(ratio.fPDGID1, ratio.fFeedDown1);
666 double dens2 = m_model->GetDensity(ratio.fPDGID2, ratio.fFeedDown2);
667 std::string name1 = m_model->TPS()->GetNameFromPDG(ratio.fPDGID1);
668 std::string name2 = m_model->TPS()->GetNameFromPDG(ratio.fPDGID2);
669 if (ratio.fPDGID1 == 1)
670 name1 = "Npart";
671 if (ratio.fPDGID2 == 1)
672 name2 = "Npart";
673 if (ratio.fPDGID1 == 33340)
674 name1 = m_model->TPS()->ParticleByPDG(3334).Name() + " + " + m_model->TPS()->ParticleByPDG(-3334).Name();
675 if (ratio.fPDGID2 == 33340)
676 name2 = m_model->TPS()->ParticleByPDG(3334).Name() + " + " + m_model->TPS()->ParticleByPDG(-3334).Name();
677 printf("%10s\t%11s\t%6lf %2s %lf\n", (std::string(name1 + "/" + name2)).c_str(), "Experiment:",
678 ratio.fValue, "+-", ratio.fError);
679 printf("%10s\t%11s\t%6lf %2s %lf\n", "", "Model:", dens1 / dens2, "+-", 0.);
680 printf("%10s\t%11s\t%6lf\n", "", "Deviation:", (dens1 / dens2 - ratio.fValue) / ratio.fError);
681 printf("\n");
682 }
683 }
684 }
685
686 void ThermalModelFit::PrintYieldsTable(std::string filename) {
687 FILE *f = fopen(filename.c_str(), "w");
688
689 fprintf(f, "%15s\t%25s\t%15s\t%15s\t%15s\t%15s\t%15s\t%15s\t%15s\t%15s\t%15s\n",
690 "N", "Name", "Data", "Error", "Fit", "Deviation", "Deviation2",
691 "Residual", "ResidualError", "Data/Model", "Data/ModelError");
692
693 //m_model->SetParameters(m_Parameters.GetThermalModelParameters());
694 //m_model->SetQoverB(m_QBgoal);
695 //m_model->FixParameters();
696 //m_model->CalculateDensities();
697 for(size_t i=0;i<m_Multiplicities.size();++i) {
698 double dens1 = m_model->GetDensity(m_Multiplicities[i].fPDGID, m_Multiplicities[i].fFeedDown);
699
700 std::string tname = m_model->TPS()->GetNameFromPDG(m_Multiplicities[i].fPDGID);
701 if (m_Multiplicities[i].fPDGID==1) tname = "Npart";
702 else if (m_Multiplicities[i].fPDGID==33340) tname = m_model->TPS()->ParticleByPDG(3334).Name() + " + " + m_model->TPS()->ParticleByPDG(-3334).Name();
703 fprintf(f, "%15d\t%25s\t%15lf\t%15lf\t%15lf\t%15lf\t%15lf\t%15lf\t%15lf\t%15lf\t%15lf\n", static_cast<int>(i)+1,
704 tname.c_str(), m_Multiplicities[i].fValue, m_Multiplicities[i].fError, dens1 * m_model->Parameters().V,
705 (m_Multiplicities[i].fValue-dens1 * m_model->Parameters().V)/m_Multiplicities[i].fError,
706 -(m_Multiplicities[i].fValue-dens1 * m_model->Parameters().V)/m_Multiplicities[i].fError,
707 (dens1 * m_model->Parameters().V-m_Multiplicities[i].fValue)/(dens1 * m_model->Parameters().V),
708 m_Multiplicities[i].fError/(dens1 * m_model->Parameters().V),
709 m_Multiplicities[i].fValue/(dens1 * m_model->Parameters().V),
710 m_Multiplicities[i].fError/(dens1 * m_model->Parameters().V));
711 }
712 for(size_t i=0;i<m_Ratios.size();++i) {
713 double dens1 = m_model->GetDensity(m_Ratios[i].fPDGID1, m_Ratios[i].fFeedDown1);
714 double dens2 = m_model->GetDensity(m_Ratios[i].fPDGID2, m_Ratios[i].fFeedDown2);
715 std::string name1 = m_model->TPS()->GetNameFromPDG(m_Ratios[i].fPDGID1);
716 std::string name2 = m_model->TPS()->GetNameFromPDG(m_Ratios[i].fPDGID2);
717 if (m_Ratios[i].fPDGID1==1) name1 = "Npart";
718 if (m_Ratios[i].fPDGID2==1) name2 = "Npart";
719 if (m_Ratios[i].fPDGID1==33340) name1 = m_model->TPS()->ParticleByPDG(3334).Name() + " + " + m_model->TPS()->ParticleByPDG(-3334).Name();
720 if (m_Ratios[i].fPDGID2==33340) name2 = m_model->TPS()->ParticleByPDG(3334).Name() + " + " + m_model->TPS()->ParticleByPDG(-3334).Name();
721
722 fprintf(f, "%15d\t%25s\t%15lf\t%15lf\t%15lf\t%15lf\t%15lf\t%15lf\t%15lf\t%15lf\t%15lf\n", static_cast<int>(i)+1,
723 (std::string(name1 + "/" + name2)).c_str(), m_Ratios[i].fValue, m_Ratios[i].fError, dens1 / dens2,
724 (m_Ratios[i].fValue - dens1 / dens2)/m_Ratios[i].fError,
725 -(m_Ratios[i].fValue - dens1 / dens2)/m_Ratios[i].fError,
726 (dens1 / dens2-m_Ratios[i].fValue)/(dens1 / dens2),
727 m_Ratios[i].fError/(dens1 / dens2),
728 m_Ratios[i].fValue/(dens1 / dens2),
729 m_Ratios[i].fError/(dens1 / dens2));
730 }
731
732 fclose(f);
733 }
734
735 void ThermalModelFit::PrintYieldsTable2(std::string filename) {
736 FILE *f = fopen(filename.c_str(), "w");
737
738 fprintf(f, "%15s\t%25s\t%15s\t%15s\t%15s\t%15s\n", "N", "Name", "Data", "Error", "Fit", "Deviation");
739 //m_model->SetParameters(m_Parameters.GetThermalModelParameters());
740 //m_model->SetQoverB(m_QBgoal);
741 //m_model->FixParameters();
742 //m_model->CalculateDensities();
743 for(size_t i=0;i<m_Multiplicities.size();++i) {
744
745 double dens1 = m_model->GetDensity(m_Multiplicities[i].fPDGID, m_Multiplicities[i].fFeedDown);
746 std::string tname = m_model->TPS()->GetNameFromPDG(m_Multiplicities[i].fPDGID);
747 if (m_Multiplicities[i].fPDGID==1) tname = "Npart";
748 else if (m_Multiplicities[i].fPDGID==33340) tname = m_model->TPS()->ParticleByPDG(3334).Name() + " + " + m_model->TPS()->ParticleByPDG(-3334).Name();
749 fprintf(f, "%15lf\t%25s\t%15lf\t%15lf\t%15lf\t%15lf\n", i+1 - 0.3, tname.c_str(), m_Multiplicities[i].fValue, m_Multiplicities[i].fError, dens1 * m_model->Parameters().V, (m_Multiplicities[i].fValue-dens1 * m_model->Parameters().V)/m_Multiplicities[i].fError);
750 fprintf(f, "%15lf\t%25s\t%15lf\t%15lf\t%15lf\t%15lf\n", i+1 + 0.3, tname.c_str(), m_Multiplicities[i].fValue, m_Multiplicities[i].fError, dens1 * m_model->Parameters().V, (m_Multiplicities[i].fValue-dens1 * m_model->Parameters().V)/m_Multiplicities[i].fError);
751 }
752 for(size_t i=0;i<m_Ratios.size();++i) {
753 double dens1 = m_model->GetDensity(m_Ratios[i].fPDGID1, m_Ratios[i].fFeedDown1);
754 double dens2 = m_model->GetDensity(m_Ratios[i].fPDGID2, m_Ratios[i].fFeedDown2);
755
756 std::string name1 = m_model->TPS()->GetNameFromPDG(m_Ratios[i].fPDGID1);
757 std::string name2 = m_model->TPS()->GetNameFromPDG(m_Ratios[i].fPDGID2);
758 if (m_Ratios[i].fPDGID1==1) name1 = "Npart";
759 if (m_Ratios[i].fPDGID2==1) name2 = "Npart";
760 if (m_Ratios[i].fPDGID1==33340) name1 = m_model->TPS()->ParticleByPDG(3334).Name() + " + " + m_model->TPS()->ParticleByPDG(-3334).Name();
761 if (m_Ratios[i].fPDGID2==33340) name2 = m_model->TPS()->ParticleByPDG(3334).Name() + " + " + m_model->TPS()->ParticleByPDG(-3334).Name();
762
763 fprintf(f, "%15lf\t%25s\t%15lf\t%15lf\t%15lf\t%15lf\n", i+1 - 0.3, (std::string(name1 + "/" + name2)).c_str(), m_Ratios[i].fValue, m_Ratios[i].fError, dens1 / dens2, (m_Ratios[i].fValue - dens1 / dens2)/m_Ratios[i].fError);
764 fprintf(f, "%15lf\t%25s\t%15lf\t%15lf\t%15lf\t%15lf\n", i+1 + 0.3, (std::string(name1 + "/" + name2)).c_str(), m_Ratios[i].fValue, m_Ratios[i].fError, dens1 / dens2, (m_Ratios[i].fValue - dens1 / dens2)/m_Ratios[i].fError);
765 }
766
767 fclose(f);
768 }
769
770 void ThermalModelFit::PrintYieldsLatex(std::string filename, std::string name) {
771 FILE *f = fopen(filename.c_str(), "w");
772 fprintf(f, "\\begin{tabular}{|c|c|c|c|}\n");
773 fprintf(f, "\\hline\n");
774 fprintf(f, "\\multicolumn{4}{|c|}{%s} \\\\\n", name.c_str());
775 fprintf(f, "\\hline\n");
776 fprintf(f, "Yield & Measurement & Fit & Deviation \\\\\n");
777 fprintf(f, "\\hline\n");
778 //m_model->SetParameters(m_Parameters.GetThermalModelParameters());
779 //m_model->SetQoverB(m_QBgoal);
780 //m_model->FixParameters();
781 //m_model->CalculateDensities();
782 for(size_t i=0;i<m_Multiplicities.size();++i) {
783 double dens1 = m_model->GetDensity(m_Multiplicities[i].fPDGID, m_Multiplicities[i].fFeedDown);
784
785 std::string tname = m_model->TPS()->GetNameFromPDG(m_Multiplicities[i].fPDGID);
786
787 if (m_Multiplicities[i].fPDGID==1) tname = "Npart";
788 else if (m_Multiplicities[i].fPDGID==33340) tname = m_model->TPS()->ParticleByPDG(3334).Name() + " + " + m_model->TPS()->ParticleByPDG(-3334).Name();
789 fprintf(f, "$%s$ & $%.4lf \\pm %.4lf$ & $%.4lf$ & $%.4lf$ \\\\\n", tname.c_str(), m_Multiplicities[i].fValue, m_Multiplicities[i].fError, dens1 * m_model->Parameters().V, (dens1 * m_model->Parameters().V-m_Multiplicities[i].fValue)/m_Multiplicities[i].fError);
790 }
791 for(size_t i=0;i<m_Ratios.size();++i) {
792 int ind1 = m_model->TPS()->PdgToId(m_Ratios[i].fPDGID1);
793 int ind2 = m_model->TPS()->PdgToId(m_Ratios[i].fPDGID2);
794 double dens1 = m_model->GetDensity(m_Ratios[i].fPDGID1, m_Ratios[i].fFeedDown1);
795 double dens2 = m_model->GetDensity(m_Ratios[i].fPDGID2, m_Ratios[i].fFeedDown2);
796
797 std::string name1 = m_model->TPS()->Particles()[ind1].Name();
798 std::string name2 = m_model->TPS()->Particles()[ind2].Name();
799 if (m_Ratios[i].fPDGID1==1) name1 = "Npart";
800 if (m_Ratios[i].fPDGID2==1) name2 = "Npart";
801 if (m_Ratios[i].fPDGID1==33340) name1 = m_model->TPS()->Particles()[m_model->TPS()->PdgToId(3334)].Name() + " + " + m_model->TPS()->Particles()[m_model->TPS()->PdgToId(-3334)].Name();
802 if (m_Ratios[i].fPDGID2==33340) name2 = m_model->TPS()->Particles()[m_model->TPS()->PdgToId(3334)].Name() + " + " + m_model->TPS()->Particles()[m_model->TPS()->PdgToId(-3334)].Name();
803 fprintf(f, "$%s$ & $%.4lf \\pm %.4lf$ & $%.4lf$ & $%.4lf$ \\\\\n", (std::string(name1 + "/" + name2)).c_str(), m_Ratios[i].fValue, m_Ratios[i].fError, dens1 / dens2 , (dens1 / dens2 - m_Ratios[i].fValue)/m_Ratios[i].fError);
804 }
805 fprintf(f, "\\hline\n");
806 fprintf(f, "\\end{tabular}\n");
807
808 fprintf(f, "\\begin{tabular}{|c|c|}\n");
809 fprintf(f, "\\hline\n");
810 fprintf(f, "\\multicolumn{2}{|c|}{%s} \\\\\n", name.c_str());
811 fprintf(f, "\\hline\n");
812 fprintf(f, "Parameter & Fit result \\\\\n");
813 fprintf(f, "\\hline\n");
814 for(int i=0;i<6;++i) {
815 FitParameter param = m_Parameters.GetParameter(i);
816 if (param.toFit) {
817 double tval = param.value, terr = param.error;
818 if (param.name=="T" || param.name.substr(0,2)=="mu") { tval *= 1000.; terr *= 1000.; }
819 fprintf(f, "$%s$ & $%.2lf \\pm %.2lf$ \\\\\n", param.name.c_str(), tval, terr);
820 }
821 }
822 fprintf(f, "$%s$ & $%.2lf/%d$ \\\\\n", "\\chi^2/N_{\\rm df}", m_Parameters.chi2ndf * GetNdf(), GetNdf());
823 fprintf(f, "\\hline\n");
824 fprintf(f, "\\end{tabular}\n");
825
826 fclose(f);
827 }
828
829 void ThermalModelFit::PrintYieldsLatexAll(std::string filename, std::string name, bool asymm) {
830 FILE *f = fopen(filename.c_str(), "w");
831 fprintf(f, "\\begin{tabular}{|c|c|c|c|}\n");
832 fprintf(f, "\\hline\n");
833 fprintf(f, "\\multicolumn{4}{|c|}{%s} \\\\\n", name.c_str());
834 fprintf(f, "\\hline\n");
835 fprintf(f, "Yield & Measurement & Fit \\\\\n");
836 fprintf(f, "\\hline\n");
837 //m_model->SetParameters(m_Parameters.GetThermalModelParameters());
838 //m_model->SetQoverB(m_QBgoal);
839 //m_model->FixParameters();
840 //m_model->CalculateDensities();
841
842 std::vector<std::string> prt(0);
843 std::vector<long long> pdgs(0);
844 std::vector<int> fl(m_Multiplicities.size(), 0);
845 prt.push_back("\\pi^+"); pdgs.push_back(211);
846 prt.push_back("\\pi^-"); pdgs.push_back(-211);
847 prt.push_back("K^+"); pdgs.push_back(321);
848 prt.push_back("K^-"); pdgs.push_back(-321);
849 prt.push_back("p"); pdgs.push_back(2212);
850 prt.push_back("\\bar{p}"); pdgs.push_back(-2212);
851 prt.push_back("\\Lambda"); pdgs.push_back(3122);
852 prt.push_back("\\bar{\\Lambda}"); pdgs.push_back(-3122);
853 prt.push_back("\\Sigma^+"); pdgs.push_back(3222);
854 prt.push_back("\\bar{\\Sigma}^+"); pdgs.push_back(-3222);
855 prt.push_back("\\Sigma^-"); pdgs.push_back(3112);
856 prt.push_back("\\bar{\\Sigma}^-"); pdgs.push_back(-3112);
857 prt.push_back("\\Xi^0"); pdgs.push_back(3322);
858 prt.push_back("\\bar{\\Xi}^0"); pdgs.push_back(-3322);
859 prt.push_back("\\Xi^-"); pdgs.push_back(3312);
860 prt.push_back("\\bar{\\Xi}^-"); pdgs.push_back(-3312);
861 prt.push_back("\\Omega"); pdgs.push_back(3334);
862 prt.push_back("\\bar{\\Omega}"); pdgs.push_back(-3334);
863 prt.push_back("\\pi^0"); pdgs.push_back(111);
864 prt.push_back("K^0_S"); pdgs.push_back(310);
865 prt.push_back("\\eta"); pdgs.push_back(221);
866 prt.push_back("\\omega"); pdgs.push_back(223);
867 prt.push_back("K^{*+}"); pdgs.push_back(323);
868 prt.push_back("K^{*-}"); pdgs.push_back(-323);
869 prt.push_back("K^{*0}"); pdgs.push_back(313);
870 prt.push_back("\\bar{K^{*0}}"); pdgs.push_back(-313);
871 //prt.push_back("K^*(892)"); pdgs.push_back(323);
872 prt.push_back("\\rho^+"); pdgs.push_back(213);
873 prt.push_back("\\rho^-"); pdgs.push_back(-213);
874 prt.push_back("\\rho^0"); pdgs.push_back(113);
875 prt.push_back("\\eta'"); pdgs.push_back(331);
876 prt.push_back("\\phi"); pdgs.push_back(333);
877 prt.push_back("\\Lambda(1520)"); pdgs.push_back(3124);
878
879
880
881 prt.push_back("d"); pdgs.push_back(1000010020);
882 prt.push_back("^3He"); pdgs.push_back(1000020030);
883 prt.push_back("^3_{\\Lambda}H"); pdgs.push_back(1010010030);
884 prt.push_back("D^+"); pdgs.push_back(411);
885 prt.push_back("D^-"); pdgs.push_back(-411);
886 prt.push_back("D^0"); pdgs.push_back(421);
887 prt.push_back("\\bar{D}^0"); pdgs.push_back(-421);
888
889 for(size_t i=0;i<prt.size();++i) {
890 bool isexp = false;
891 int tj = -1;
892 for(size_t j=0;j<m_Multiplicities.size();++j)
893 if (pdgs[i]==m_Multiplicities[j].fPDGID) { isexp = true; fl[j] = 1; tj = j; break; }
894
895 if (m_model->TPS()->PdgToId(pdgs[i]) == -1) continue;
896
897 double dens1 = 0.;
898 dens1 = m_model->GetDensity(pdgs[i], Feeddown::StabilityFlag);
899
900 if (isexp) dens1 = m_model->GetDensity(pdgs[i], m_Multiplicities[tj].fFeedDown);
901
902 std::string tname = prt[i];
903
904 if (!isexp) {
905 if (dens1 * m_model->Parameters().V>1.e-3) fprintf(f, "$%s$ & & $%.3g$ \\\\\n", tname.c_str(), dens1 * m_model->Parameters().V);
906 else {
907 char tnum[400];
908 sprintf(tnum, "%.2E", dens1 * m_model->Parameters().V);
909 double val = 0.;
910 int step = 0;
911 char valst[6];
912 char stepst[6];
913 for(int ii=0;ii<4;++ii) valst[ii] = tnum[ii];
914 valst[4] = 0;
915 for(int ii=0;ii<4;++ii) stepst[ii] = tnum[5+ii];
916 stepst[4] = 0;
917 val = atof(valst);
918 step = atoi(stepst);
919 fprintf(f, "$%s$ & & $%.2lf \\cdot 10^{%d}$ \\\\\n", tname.c_str(), val, step);
920 }
921 }
922 else {
923 if (m_Multiplicities[tj].fValue>0.999) fprintf(f, "$%s$ & $%.4g \\pm %.4g$ & $%.4g$ \\\\\n", tname.c_str(), m_Multiplicities[tj].fValue, m_Multiplicities[tj].fError, dens1 * m_model->Parameters().V);//, (dens1 * m_model->Parameters().V-m_Multiplicities[tj].fValue)/m_Multiplicities[tj].fError);
924 else if (dens1 * m_model->Parameters().V>1.e-3) fprintf(f, "$%s$ & $%.3g \\pm %.3g$ & $%.3g$ \\\\\n", tname.c_str(), m_Multiplicities[tj].fValue, m_Multiplicities[tj].fError, dens1 * m_model->Parameters().V);//, (dens1 * m_model->Parameters().V-m_Multiplicities[tj].fValue)/m_Multiplicities[tj].fError);
925 else fprintf(f, "$%s$ & $%.6g \\pm %.6g$ & $%.3g$ \\\\\n", tname.c_str(), m_Multiplicities[tj].fValue, m_Multiplicities[tj].fError, dens1 * m_model->Parameters().V);//, (dens1 * m_model->Parameters().V-m_Multiplicities[tj].fValue)/m_Multiplicities[tj].fError);
926 }
927 }
928
929 for(size_t i=0;i<m_Multiplicities.size();++i)
930 if (!fl[i]) {
931 double dens1 = m_model->GetDensity(m_Multiplicities[i].fPDGID, m_Multiplicities[i].fFeedDown);
932
933 std::string tname = m_model->TPS()->GetNameFromPDG(m_Multiplicities[i].fPDGID);
934 if (m_Multiplicities[i].fPDGID==1) tname = "Npart";
935 else if (m_Multiplicities[i].fPDGID==33340) tname = m_model->TPS()->Particles()[m_model->TPS()->PdgToId(3334)].Name() + " + " + m_model->TPS()->Particles()[m_model->TPS()->PdgToId(-3334)].Name();
936 fprintf(f, "$%s$ & $%.4lf \\pm %.4lf$ & $%.4lf$ & $%.4lf$ \\\\\n", tname.c_str(), m_Multiplicities[i].fValue, m_Multiplicities[i].fError, dens1 * m_model->Parameters().V, (dens1 * m_model->Parameters().V-m_Multiplicities[i].fValue)/m_Multiplicities[i].fError);
937 }
938 for(size_t i=0;i<m_Ratios.size();++i) {
939 double dens1 = m_model->GetDensity(m_Ratios[i].fPDGID1, m_Ratios[i].fFeedDown1);
940 double dens2 = m_model->GetDensity(m_Ratios[i].fPDGID2, m_Ratios[i].fFeedDown2);
941
942 std::string name1 = m_model->TPS()->GetNameFromPDG(m_Ratios[i].fPDGID1);
943 std::string name2 = m_model->TPS()->GetNameFromPDG(m_Ratios[i].fPDGID2);
944 if (m_Ratios[i].fPDGID1==1) name1 = "Npart";
945 if (m_Ratios[i].fPDGID2==1) name2 = "Npart";
946 if (m_Ratios[i].fPDGID1==33340) name1 = m_model->TPS()->Particles()[m_model->TPS()->PdgToId(3334)].Name() + " + " + m_model->TPS()->Particles()[m_model->TPS()->PdgToId(-3334)].Name();
947 if (m_Ratios[i].fPDGID2==33340) name2 = m_model->TPS()->Particles()[m_model->TPS()->PdgToId(3334)].Name() + " + " + m_model->TPS()->Particles()[m_model->TPS()->PdgToId(-3334)].Name();
948 fprintf(f, "$%s$ & $%.4lf \\pm %.4lf$ & $%.4lf$ & $%.4lf$ \\\\\n", (std::string(name1 + "/" + name2)).c_str(), m_Ratios[i].fValue, m_Ratios[i].fError, dens1 / dens2 , (dens1 / dens2 - m_Ratios[i].fValue)/m_Ratios[i].fError);
949 }
950 fprintf(f, "\\hline\n");
951 fprintf(f, "\\end{tabular}\n");
952
953 fprintf(f, "\\begin{tabular}{|c|c|}\n");
954 fprintf(f, "\\hline\n");
955 fprintf(f, "\\multicolumn{2}{|c|}{%s} \\\\\n", name.c_str());
956 fprintf(f, "\\hline\n");
957 fprintf(f, "Parameter & Fit result \\\\\n");
958 fprintf(f, "\\hline\n");
959 for(int i=0;i<9;++i) {
960 FitParameter param = m_Parameters.GetParameter(i);
961 if (param.toFit) {
962 double tval = param.value, terr = param.error, terrp = param.errp, terrm = param.errm;
963 if (param.name=="T" || param.name.substr(0,2)=="mu") { tval *= 1000.; terr *= 1000.; terrp *= 1000.; terrm *= 1000.; }
964 if (!asymm) fprintf(f, "$%s$ & $%.3lf \\pm %.3lf$ \\\\\n", param.name.c_str(), tval, terr);
965 else fprintf(f, "$%s$ & $%.3lf^{+%.3lf}_{-%.3lf}$ \\\\\n", param.name.c_str(), tval, terrp, terrm);
966 }
967 }
968 fprintf(f, "$%s$ & $%.2lf/%d$ \\\\\n", "\\chi^2/N_{\\rm df}", m_Parameters.chi2ndf * GetNdf(), GetNdf());
969 fprintf(f, "\\hline\n");
970 fprintf(f, "\\end{tabular}\n");
971
972 fclose(f);
973 }
974
975 std::string ThermalModelFit::GetCurrentTime() {
976 time_t rawtime;
977 struct tm * timeinfo;
978 char buffer[1000];
979
980 time(&rawtime);
981 timeinfo = localtime(&rawtime);
982
983 strftime(buffer, sizeof(buffer), "%d-%m-%Y at %I:%M:%S", timeinfo);
984 std::string str(buffer);
985
986 return str;
987 }
988
989 void ThermalModelFit::PrintFitLog(std::string filename, std::string comment, bool asymm)
990 {
991 std::ostream *fout;
992 if (filename != "")
993 fout = new std::ofstream(filename.c_str());
994 else
995 fout = &std::cout;
996
997 //FILE *f;
998 //if (filename != "") {
999 // f = fopen(filename.c_str(), "w");
1000 //}
1001 //else {
1002 // f = stdout;
1003 //}
1004
1005 if (comment != "")
1006 *fout << comment << std::endl << std::endl;
1007 //fprintf(f, "%s\n\n", comment.c_str());
1008
1009 //fprintf(f, "Performed on %s\n\n", GetCurrentTime().c_str());
1010 *fout << "Performed on " << GetCurrentTime() << std::endl << std::endl;
1011
1012 //fprintf(f, "chi2/dof = %lf/%d = %lf\n\n", m_Parameters.chi2, m_Parameters.ndf, m_Parameters.chi2ndf);
1013 *fout << "chi2/dof = " << m_Parameters.chi2 << "/" << m_Parameters.ndf << " = " << m_Parameters.chi2ndf << std::endl << std::endl;
1014
1015 std::pair<double, double> accuracy = ModelDescriptionAccuracy();
1016 //fprintf(f, "Data description accuracy = (%.2lf +- %.2lf) %%\n\n", accuracy.first * 100., accuracy.second * 100.);
1017 *fout << "Data description accuracy = ("
1018 << std::fixed
1019 << std::setprecision(2)
1020 << accuracy.first * 100. << " +- " << accuracy.second * 100. << ") %" << std::endl << std::endl;
1021
1022 fout->unsetf(std::ios_base::floatfield);
1023 *fout << std::setprecision(6);
1024
1025 //fprintf(f, "Extracted parameters:\n");
1026 *fout << "Extracted parameters:" << std::endl;
1027 for (int i = 0; i < 10 ; ++i) {
1028 if (i == 6 && m_model->Ensemble() != ThermalModelBase::SCE && m_model->Ensemble() != ThermalModelBase::CE)
1029 continue;
1030 FitParameter param = m_Parameters.GetParameter(i);
1031 std::string sunit = "";
1032 std::string tname = "";
1033 double tmn = 1.;
1034 if (param.name == "T" || param.name.substr(0, 2) == "mu") {
1035 sunit = "[MeV]";
1036 tmn = 1.e3;
1037 }
1038 else if (param.name == "R" || param.name == "Rc") {
1039 sunit = "[fm]";
1040 }
1041 tname = param.name + sunit;
1042
1043 if (param.name == "Rc" && m_model->Ensemble() == ThermalModelBase::GCE)
1044 continue;
1045
1046 if (param.name == "muS" && (m_model->Ensemble() == ThermalModelBase::SCE || m_model->Ensemble() == ThermalModelBase::CE))
1047 continue;
1048
1049 if (param.name == "muB" && m_model->Ensemble() == ThermalModelBase::CE)
1050 continue;
1051
1052 if (param.name == "muQ" && m_model->Ensemble() == ThermalModelBase::CE)
1053 continue;
1054
1055 if ((param.name == "muS" || param.name == "gammaS") && !m_model->TPS()->hasStrange())
1056 continue;
1057 if (param.name == "muQ" && !m_model->TPS()->hasCharged())
1058 continue;
1059 if ((param.name == "muC" || param.name == "gammaC") && !m_model->TPS()->hasCharmed())
1060 continue;
1061
1062 if (param.toFit) {
1063 double tval = param.value, terr = param.error, terrp = param.errp, terrm = param.errm;
1064
1065 if (param.name == "Tkin" && !m_YieldsAtTkin)
1066 continue;
1067
1068 if (param.name != "R" && param.name != "Rc") {
1069 if (!asymm)
1070 *fout << std::setw(15) << tname
1071 << " = "
1072 << std::setw(15) << tval * tmn
1073 << " +- "
1074 << std::left << std::setw(15) << terr * tmn << std::right
1075 << std::endl;
1076 //fprintf(f, "%15s = %15lf +- %-15lf\n", tname.c_str(), tval * tmn, terr * tmn);
1077 else
1078 *fout << std::setw(15) << tname
1079 << " = "
1080 << std::setw(15) << tval * tmn
1081 << " + "
1082 << std::left << std::setw(15) << terrp * tmn << std::right
1083 << " - "
1084 << std::left << std::setw(15) << terrm * tmn << std::right
1085 << std::endl;
1086 //fprintf(f, "%15s = %15lf + %-15lf - %-15lf\n", tname.c_str(), tval * tmn, terrp * tmn, terrm * tmn);
1087 }
1088 else {
1089 if (!asymm) {
1090 *fout << std::setw(15) << tname
1091 << " = "
1092 << std::setw(15) << tval * tmn
1093 << " +- "
1094 << std::left << std::setw(15) << terr * tmn << std::right
1095 << "\t";
1096 //fprintf(f, "%15s = %15lf +- %-15lf\t", tname.c_str(), tval * tmn, terr * tmn);
1097 std::string tname2 = "V[fm3]";
1098 if (param.name == "Rc")
1099 tname2 = "Vc[fm3]";
1100 *fout << std::setw(15) << tname2
1101 << " = "
1102 << std::setw(15) << 4. / 3. * xMath::Pi() * pow(tval * tmn, 3)
1103 << " +- "
1104 << std::left << std::setw(15) << 4. * xMath::Pi() * pow(tval * tmn, 2) * terr * tmn << std::right
1105 << std::endl;
1106 //fprintf(f, "%15s = %15lf +- %-15lf\n", tname2.c_str(), 4. / 3. * xMath::Pi() * pow(tval * tmn, 3), 4. * xMath::Pi() * pow(tval * tmn, 2) * terr * tmn);
1107 }
1108 else {
1109 *fout << std::setw(15) << tname
1110 << " = "
1111 << std::setw(15) << tval * tmn
1112 << " + "
1113 << std::left << std::setw(15) << terrp * tmn << std::right
1114 << " - "
1115 << std::left << std::setw(15) << terrm * tmn << std::right
1116 << "\t";
1117 //fprintf(f, "%15s = %15lf + %-15lf - %-15lf\t", tname.c_str(), tval * tmn, terrp * tmn, terrm * tmn);
1118 std::string tname2 = "V[fm3]";
1119 if (param.name == "Rc")
1120 tname2 = "Vc[fm3]";
1121 *fout << std::setw(15) << tname2
1122 << " = "
1123 << std::setw(15) << 4. / 3. * xMath::Pi() * pow(tval * tmn, 3)
1124 << " + "
1125 << std::left << std::setw(15) << 4. * xMath::Pi() * pow(tval * tmn, 2) * terrp * tmn << std::right
1126 << " - "
1127 << std::left << std::setw(15) << 4. * xMath::Pi() * pow(tval * tmn, 2) * terrm * tmn << std::right
1128 << std::endl;
1129 //fprintf(f, "%15s = %15lf + %-15lf - %-15lf\n", tname2.c_str(), 4. / 3. * xMath::Pi() * pow(tval * tmn, 3), 4. * xMath::Pi() * pow(tval * tmn, 2) * terrp * tmn, 4. * xMath::Pi() * pow(tval * tmn, 2) * terrm * tmn);
1130 }
1131 }
1132 }
1133 else {
1134 double tval = param.value;
1135 if (param.name != "R" && param.name != "Rc")
1136 *fout << std::setw(15) << tname
1137 << " = "
1138 << std::setw(15) << tval * tmn
1139 << " (FIXED)" << std::endl;
1140 //fprintf(f, "%15s = %15lf (FIXED)\n", tname.c_str(), tval * tmn);
1141 else {
1142 //fprintf(f, "%15s = %15lf (FIXED)\t", tname.c_str(), tval * tmn);
1143 *fout << std::setw(15) << tname
1144 << " = "
1145 << std::setw(15) << tval * tmn
1146 << " (FIXED)" << "\t";
1147 std::string tname2 = "V[fm3]";
1148 if (param.name == "Rc")
1149 tname2 = "Vc[fm3]";
1150 *fout << std::setw(15) << tname2
1151 << " = "
1152 << std::setw(15) << 4. / 3. * xMath::Pi() * pow(tval * tmn, 3)
1153 << " (FIXED)" << std::endl;
1154 //fprintf(f, "%15s = %15lf (FIXED)\n", tname2.c_str(), 4. / 3. * xMath::Pi() * pow(tval * tmn, 3));
1155 }
1156 }
1157 }
1158
1159 //fprintf(f, "\n\n");
1160 *fout << std::endl << std::endl;
1161
1162 *fout << "Yields:" << std::endl;
1163 //fprintf(f, "Yields:\n");
1164 for (size_t i = 0; i<m_Multiplicities.size(); ++i) {
1165 double dens1 = m_model->GetDensity(m_Multiplicities[i].fPDGID, m_Multiplicities[i].fFeedDown);
1166
1167 std::string tname = m_model->TPS()->GetNameFromPDG(m_Multiplicities[i].fPDGID);
1168 if (m_Multiplicities[i].fPDGID == 1) tname = "Npart";
1169 else if (m_Multiplicities[i].fPDGID == 33340) tname = m_model->TPS()->Particles()[m_model->TPS()->PdgToId(3334)].Name() + " + " + m_model->TPS()->Particles()[m_model->TPS()->PdgToId(-3334)].Name();
1170
1171 *fout << std::setw(25) << tname
1172 << " Experiment: ";
1173
1174 if (m_Multiplicities[i].fValue > 1.e-5 && m_Multiplicities[i].fError > 1.e-5)
1175 fout->unsetf(std::ios_base::floatfield);
1176 else
1177 *fout << std::scientific;
1178
1179 *fout << std::setw(15) << m_Multiplicities[i].fValue;
1180 *fout << " +- ";
1181 *fout << std::left << std::setw(15) << m_Multiplicities[i].fError << std::right << " ";
1182 *fout << "Model: ";
1183 *fout << std::left << std::setw(15) << dens1 * m_model->Parameters().V << std::right << " ";
1184
1185 fout->unsetf(std::ios_base::floatfield);
1186
1187 *fout << "Std.dev.: ";
1188 *fout << std::left << std::setw(15) << (dens1* m_model->Parameters().V - m_Multiplicities[i].fValue) / m_Multiplicities[i].fError << std::right << " ";
1189 *fout << std::endl;
1190
1191 //if (m_Multiplicities[i].fValue > 1.e-5 && m_Multiplicities[i].fError > 1.e-5)
1192 // fprintf(f, "%25s Experiment: %15lf +- %-15lf Model: %-15lf Std.dev.: %-15lf \n", tname.c_str(), m_Multiplicities[i].fValue, m_Multiplicities[i].fError, dens1 * m_model->Parameters().V, (dens1 * m_model->Parameters().V - m_Multiplicities[i].fValue) / m_Multiplicities[i].fError);
1193 //else
1194 // fprintf(f, "%25s Experiment: %15E +- %-15E Model: %-15E Std.dev.: %-15lf \n", tname.c_str(), m_Multiplicities[i].fValue, m_Multiplicities[i].fError, dens1 * m_model->Parameters().V, (dens1 * m_model->Parameters().V - m_Multiplicities[i].fValue) / m_Multiplicities[i].fError);
1195 }
1196
1197 for (size_t i = 0; i<m_Ratios.size(); ++i) {
1198 double dens1 = m_model->GetDensity(m_Ratios[i].fPDGID1, m_Ratios[i].fFeedDown1);
1199 double dens2 = m_model->GetDensity(m_Ratios[i].fPDGID2, m_Ratios[i].fFeedDown2);
1200
1201 std::string name1 = m_model->TPS()->GetNameFromPDG(m_Ratios[i].fPDGID1);
1202 std::string name2 = m_model->TPS()->GetNameFromPDG(m_Ratios[i].fPDGID2);
1203 if (m_Ratios[i].fPDGID1 == 1) name1 = "Npart";
1204 if (m_Ratios[i].fPDGID2 == 1) name2 = "Npart";
1205 if (m_Ratios[i].fPDGID1 == 33340) name1 = m_model->TPS()->Particles()[m_model->TPS()->PdgToId(3334)].Name() + " + " + m_model->TPS()->Particles()[m_model->TPS()->PdgToId(-3334)].Name();
1206 if (m_Ratios[i].fPDGID2 == 33340) name2 = m_model->TPS()->Particles()[m_model->TPS()->PdgToId(3334)].Name() + " + " + m_model->TPS()->Particles()[m_model->TPS()->PdgToId(-3334)].Name();
1207
1208 *fout << std::setw(25) << std::string(name1 + "/" + name2)
1209 << " Experiment: ";
1210
1211 if (m_Ratios[i].fValue > 1.e-5 && m_Ratios[i].fError > 1.e-5)
1212 fout->unsetf(std::ios_base::floatfield);
1213 else
1214 *fout << std::scientific;
1215
1216 *fout << std::setw(15) << m_Ratios[i].fValue;
1217 *fout << " +- ";
1218 *fout << std::left << std::setw(15) << m_Ratios[i].fError << std::right << " ";
1219 *fout << "Model: ";
1220 *fout << std::left << std::setw(15) << dens1 / dens2 << std::right << " ";
1221
1222 fout->unsetf(std::ios_base::floatfield);
1223
1224 *fout << "Std.dev.: ";
1225 *fout << std::left << std::setw(15) << (dens1 / dens2 - m_Ratios[i].fValue) / m_Ratios[i].fError << std::right << " ";
1226 *fout << std::endl;
1227
1228 //if (m_Ratios[i].fValue > 1.e-5 && m_Ratios[i].fError > 1.e-5)
1229 // fprintf(f, "%25s Experiment: %15lf +- %-15lf Model: %-15lf Std.dev.: %-15lf \n", (std::string(name1 + "/" + name2)).c_str(), m_Ratios[i].fValue, m_Ratios[i].fError, dens1 / dens2, (dens1 / dens2 - m_Ratios[i].fValue) / m_Ratios[i].fError);
1230 //else
1231 // fprintf(f, "%25s Experiment: %15E +- %-15E Model: %-15E Std.dev.: %-15lf \n", (std::string(name1 + "/" + name2)).c_str(), m_Ratios[i].fValue, m_Ratios[i].fError, dens1 / dens2, (dens1 / dens2 - m_Ratios[i].fValue) / m_Ratios[i].fError);
1232 }
1233
1234 //if (f != stdout)
1235 // fclose(f);
1236 if (filename != "" && fout != NULL) {
1237 delete fout;
1238 }
1239 }
1240
1241 using namespace std;
1242
1243 std::pair<double, double> ThermalModelFit::ModelDescriptionAccuracy() const
1244 {
1245 double chi2total = 0.;
1246 std::vector<double> weights;
1247 std::vector<double> vals, errors;
1248 for (int i = 0; i < ModelDataSize(); ++i) {
1249 const FittedQuantity &quantity = FittedQuantities()[i];
1250 if (quantity.toFit) {
1251 vals.push_back(fabs(m_ModelData[i] / quantity.Value() - 1));
1252 errors.push_back(quantity.ValueError() / quantity.Value() * m_ModelData[i] / quantity.Value());
1253
1254 double chi2contrib = (m_ModelData[i] - quantity.Value()) * (m_ModelData[i] - quantity.Value()) / quantity.ValueError() / quantity.ValueError();
1255
1256 chi2total += chi2contrib;
1257 weights.push_back(chi2contrib);
1258 }
1259 }
1260
1261 for (size_t i = 0; i < weights.size(); ++i) {
1262 weights[i] /= chi2total;
1263 }
1264
1265 double mean = 0., error = 0.;
1266 for (size_t i = 0; i < vals.size(); ++i) {
1267 mean += vals[i] * weights[i];
1268 error += errors[i] * weights[i];
1269 }
1270
1271 return make_pair(mean, error);
1272 }
1273
1274 std::vector<FittedQuantity> ThermalModelFit::loadExpDataFromFile(const std::string & filename) {
1275 std::vector<FittedQuantity> ret(0);
1276 fstream fin;
1277 fin.open(filename.c_str());
1278 if (fin.is_open()) {
1279 char tmpc[2000];
1280 fin.getline(tmpc, 2000);
1281 string tmp = string(tmpc);
1282 vector<string> elems = CuteHRGHelper::split(tmp, '#');
1283
1284 int flnew = 0;
1285 if (elems.size() < 1 || CuteHRGHelper::split(elems[0], ';').size() < 4)
1286 flnew = 1;
1287 else
1288 flnew = 0;
1289
1290 fin.clear();
1291 fin.seekg(0, ios::beg);
1292
1293 if (flnew == 1)
1294 ret = loadExpDataFromFile_NewFormat(fin);
1295 else
1296 ret = loadExpDataFromFile_OldFormat(fin);
1297
1298 fin.close();
1299 }
1300 return ret;
1301 }
1302
1303 std::vector<FittedQuantity> ThermalModelFit::loadExpDataFromFile_OldFormat(std::fstream & fin)
1304 {
1305 std::vector<FittedQuantity> ret(0);
1306 if (fin.is_open()) {
1307 string tmp;
1308 char tmpc[2000];
1309 fin.getline(tmpc, 2000);
1310 tmp = string(tmpc);
1311 while (1) {
1312 vector<string> fields = CuteHRGHelper::split(tmp, ';');
1313 if (fields.size()<4) break;
1314 int type = atoi(fields[0].c_str());
1315 long long pdgid1 = stringToLongLong(fields[1]), pdgid2 = 0;
1316 double value = 0., error = 0., error2 = 0.;
1317 int feeddown1 = 1, feeddown2 = 1;
1318 if (type) {
1319 pdgid2 = stringToLongLong(fields[2]);
1320 value = atof(fields[3].c_str());
1321 if (fields.size() >= 5) error = atof(fields[4].c_str());
1322 if (fields.size() >= 6) {
1323 error2 = atof(fields[5].c_str());
1324 error = sqrt(error*error + error2*error2);
1325 }
1326 if (fields.size() >= 7) feeddown1 = atoi(fields[6].c_str());
1327 if (fields.size() >= 8) feeddown2 = atoi(fields[7].c_str());
1328 }
1329 else {
1330 value = atof(fields[2].c_str());
1331 error = atof(fields[3].c_str());
1332 if (fields.size() >= 5) {
1333 error2 = atof(fields[4].c_str());
1334 error = sqrt(error*error + error2*error2);
1335 }
1336 if (fields.size() >= 6) feeddown1 = atoi(fields[5].c_str());
1337 }
1338
1339 if (type) ret.push_back(FittedQuantity(ExperimentRatio(pdgid1, pdgid2, value, error, static_cast<Feeddown::Type>(feeddown1), static_cast<Feeddown::Type>(feeddown2))));
1340 else ret.push_back(FittedQuantity(ExperimentMultiplicity(pdgid1, value, error, static_cast<Feeddown::Type>(feeddown1))));
1341
1342 fin.getline(tmpc, 500);
1343 tmp = string(tmpc);
1344 }
1345 fin.close();
1346 }
1347 return ret;
1348 }
1349
1350 std::vector<FittedQuantity> ThermalModelFit::loadExpDataFromFile_NewFormat(std::fstream & fin)
1351 {
1352 std::vector<FittedQuantity> ret(0);
1353 if (fin.is_open()) {
1354 char cc[2000];
1355 while (!fin.eof()) {
1356 fin.getline(cc, 2000);
1357 string tmp = string(cc);
1358 vector<string> elems = CuteHRGHelper::split(tmp, '#');
1359 if (elems.size() < 1)
1360 continue;
1361
1362 istringstream iss(elems[0]);
1363
1364 int fitflag;
1365 long long pdgid1, pdgid2;
1366 int feeddown1, feeddown2;
1367 double value, error;
1368
1369 if (iss >> fitflag >> pdgid1 >> pdgid2 >> feeddown1 >> feeddown2 >> value >> error) {
1370 if (pdgid2 != 0) ret.push_back(FittedQuantity(ExperimentRatio(pdgid1, pdgid2, value, error, static_cast<Feeddown::Type>(feeddown1), static_cast<Feeddown::Type>(feeddown2))));
1371 else ret.push_back(FittedQuantity(ExperimentMultiplicity(pdgid1, value, error, static_cast<Feeddown::Type>(feeddown1))));
1372
1373 if (fitflag != 0)
1374 ret[ret.size() - 1].toFit = true;
1375 else
1376 ret[ret.size() - 1].toFit = false;
1377 }
1378 }
1379 }
1380 return ret;
1381 }
1382
1383 void ThermalModelFit::saveExpDataToFile(const std::vector<FittedQuantity>& outQuantities, const std::string & filename)
1384 {
1385 std::ofstream fout(filename.c_str());
1386 if (fout.is_open()) {
1387 fout << "#"
1388 << std::setw(14) << "is_fitted"
1389 << std::setw(15) << "pdg1"
1390 << std::setw(15) << "pdg2"
1391 << std::setw(15) << "feeddown1"
1392 << std::setw(15) << "feeddown2"
1393 << std::setw(15) << "value"
1394 << std::setw(15) << "error"
1395 << std::endl;
1396 for (size_t i = 0; i < outQuantities.size(); ++i) {
1397 if (outQuantities[i].type == FittedQuantity::Multiplicity) {
1398 fout << std::setw(15) << static_cast<int>(outQuantities[i].toFit)
1399 << std::setw(15) << outQuantities[i].mult.fPDGID
1400 << std::setw(15) << 0
1401 << std::setw(15) << outQuantities[i].mult.fFeedDown
1402 << std::setw(15) << 0
1403 << std::setw(15) << outQuantities[i].mult.fValue
1404 << std::setw(15) << outQuantities[i].mult.fError
1405 << std::endl;
1406 }
1407 else {
1408 fout << std::setw(15) << static_cast<int>(outQuantities[i].toFit)
1409 << std::setw(15) << outQuantities[i].ratio.fPDGID1
1410 << std::setw(15) << outQuantities[i].ratio.fPDGID2
1411 << std::setw(15) << outQuantities[i].ratio.fFeedDown1
1412 << std::setw(15) << outQuantities[i].ratio.fFeedDown2
1413 << std::setw(15) << outQuantities[i].ratio.fValue
1414 << std::setw(15) << outQuantities[i].ratio.fError
1415 << std::endl;
1416 }
1417 }
1418 fout.close();
1419 }
1420 else {
1421 printf("**WARNING** ThermalModelFit::saveExpDataToFile: Cannot open file for writing data!");
1422 }
1423 }
1424
1425 int ThermalModelFit::GetNdf() const {
1426 int nparams = 11;
1427 if (!m_Parameters.T.toFit) nparams--;
1428 if (!m_Parameters.muB.toFit) nparams--;
1429 if (!m_Parameters.muS.toFit) nparams--;
1430 if (!m_Parameters.muQ.toFit) nparams--;
1431 if (!m_Parameters.muC.toFit) nparams--;
1432 if (!m_Parameters.gammaq.toFit) nparams--;
1433 if (!m_Parameters.gammaS.toFit) nparams--;
1434 if (!m_Parameters.gammaC.toFit) nparams--;
1435 if (!m_Parameters.R.toFit || (m_Multiplicities.size()==0 && m_model->Ensemble() == ThermalModelBase::GCE)) nparams--;
1436 if (!m_Parameters.Rc.toFit || (m_model->Ensemble() != ThermalModelBase::CE && m_model->Ensemble() != ThermalModelBase::SCE && m_model->Ensemble() != ThermalModelBase::CCE)) nparams--;
1437 if (!m_Parameters.Tkin.toFit || !UseTkin()) nparams--;
1438 int ndof = 0;
1439 for (size_t i = 0; i < m_Quantities.size(); ++i)
1440 if (m_Quantities[i].toFit) ndof++;
1441 return (ndof - nparams);
1442 }
1443
1444
1445} // namespace thermalfist
params
map< string, double > params
Definition NeutronStars-CSHRG.cpp:34
ThermalModelFit.h
ThermalModelFit class.
Utility.h
Contains some helper functions.
thermalfist::ThermalModelBase
Abstract base class for an HRG model implementation.
Definition ThermalModelBase.h:29
thermalfist::ThermalModelBase::SCE
@ SCE
Strangeness-canonical ensemble.
Definition ThermalModelBase.h:38
thermalfist::ThermalModelBase::GCE
@ GCE
Grand canonical ensemble.
Definition ThermalModelBase.h:36
thermalfist::ThermalModelBase::CCE
@ CCE
Charm-canonical ensemble.
Definition ThermalModelBase.h:39
thermalfist::ThermalModelFit::FittedQuantities
const std::vector< FittedQuantity > & FittedQuantities() const
Return a vector of the fitted quantities (data)
Definition ThermalModelFit.h:266
thermalfist::ThermalModelFit::PrintYieldsLatexAll
void PrintYieldsLatexAll(std::string filename="Yield.dat", std::string name="A+A", bool asymm=false)
Definition ThermalModelFit.cpp:829
thermalfist::ThermalModelFit::PrintYieldsLatex
void PrintYieldsLatex(std::string filename="Yield.dat", std::string name="A+A")
Definition ThermalModelFit.cpp:770
thermalfist::ThermalModelFit::~ThermalModelFit
~ThermalModelFit(void)
Destroy the Thermal Model Fit object.
Definition ThermalModelFit.cpp:222
thermalfist::ThermalModelFit::PrintFitLog
void PrintFitLog(std::string filename="", std::string comment="Thermal fit", bool asymm=false)
Prints the result of the fitting procedure to a file.
Definition ThermalModelFit.cpp:989
thermalfist::ThermalModelFit::ThermalModelFit
ThermalModelFit(ThermalModelBase *model)
Construct a new ThermalModelFit object.
Definition ThermalModelFit.cpp:214
thermalfist::ThermalModelFit::PrintParameters
void PrintParameters()
Print function.
Definition ThermalModelFit.cpp:582
thermalfist::ThermalModelFit::loadExpDataFromFile
static std::vector< FittedQuantity > loadExpDataFromFile(const std::string &filename)
Load the experimental data from a file.
Definition ThermalModelFit.cpp:1274
thermalfist::ThermalModelFit::UseTkin
void UseTkin(bool YieldsAtTkin)
Sets whether the yields should be evaluated at Tkin using partial chemical equilibrium.
Definition ThermalModelFit.h:309
thermalfist::ThermalModelFit::GetNdf
int GetNdf() const
Number of degrees of freedom in the fit.
Definition ThermalModelFit.cpp:1425
thermalfist::ThermalModelFit::saveExpDataToFile
static void saveExpDataToFile(const std::vector< FittedQuantity > &outQuantities, const std::string &filename)
Definition ThermalModelFit.cpp:1383
thermalfist::ThermalModelFit::ModelDescriptionAccuracy
std::pair< double, double > ModelDescriptionAccuracy() const
Returns a relative error of the data description (and its uncertainty estimate)
Definition ThermalModelFit.cpp:1243
thermalfist::ThermalModelFit::Parameters
const ThermalModelFitParameters & Parameters() const
Definition ThermalModelFit.h:208
thermalfist::ThermalModelFit::PerformFit
ThermalModelFitParameters PerformFit(bool verbose=true, bool AsymmErrors=false)
The thermal fit procedure.
Definition ThermalModelFit.cpp:226
thermalfist::ThermalModelFit::PrintYieldsTable2
void PrintYieldsTable2(std::string filename="Yield.dat")
Definition ThermalModelFit.cpp:735
thermalfist::ThermalModelFit::PrintYieldsTable
void PrintYieldsTable(std::string filename="Yield.dat")
Definition ThermalModelFit.cpp:686
thermalfist::ThermalModelPCEAnnihilation
Class implementing HRG in partial chemical equilibrium with baryon annihilation.
Definition ThermalModelPCEAnnihilation.h:26
thermalfist::ThermalModelPCE
Class implementing HRG in partial chemical equilibrium.
Definition ThermalModelPCE.h:37
thermalfist::ThermalModelPCE::UseSahaForNuclei
void UseSahaForNuclei(bool flag)
Whether the nuclear abundances are evaluated through the Saha equation.
Definition ThermalModelPCE.h:77
thermalfist::CuteHRGHelper::split
std::vector< std::string > split(const std::string &s, char delim)
Definition ThermalParticleSystem.cpp:1652
thermalfist::xMath::Pi
constexpr double Pi()
Pi constant.
Definition xMath.h:23
thermalfist
The main namespace where all classes and functions of the Thermal-FIST library reside.
Definition CosmicEoS.h:9
thermalfist::stringToLongLong
long long stringToLongLong(const std::string &str)
thermalfist::ExperimentMultiplicity
Structure containing the experimental yield (multiplicity) to be fitted.
Definition ThermalModelFitQuantities.h:28
thermalfist::ExperimentMultiplicity::fFeedDown
Feeddown::Type fFeedDown
The feeddown contributions to be included.
Definition ThermalModelFitQuantities.h:39
thermalfist::ExperimentMultiplicity::fPDGID
long long fPDGID
PDG code of the particle yield.
Definition ThermalModelFitQuantities.h:30
thermalfist::ExperimentRatio
Structure containing the experimental ratio of yields to be fitted.
Definition ThermalModelFitQuantities.h:63
thermalfist::ExperimentRatio::fFeedDown2
Feeddown::Type fFeedDown2
The feeddown contributions to be included for the yield in the denominator.
Definition ThermalModelFitQuantities.h:80
thermalfist::ExperimentRatio::fPDGID1
long long fPDGID1
PDG code of the particle yield in the numerator.
Definition ThermalModelFitQuantities.h:65
thermalfist::ExperimentRatio::fValue
double fValue
Experimental value of the yield ratio.
Definition ThermalModelFitQuantities.h:71
thermalfist::ExperimentRatio::fError
double fError
Experimental error of the yield ratio.
Definition ThermalModelFitQuantities.h:74
thermalfist::ExperimentRatio::fPDGID2
long long fPDGID2
PDG code of the particle yield in the denominator.
Definition ThermalModelFitQuantities.h:68
thermalfist::ExperimentRatio::fFeedDown1
Feeddown::Type fFeedDown1
The feeddown contributions to be included for the yield in the numerator.
Definition ThermalModelFitQuantities.h:77
thermalfist::Feeddown::StabilityFlag
@ StabilityFlag
Feeddown from all particles marked as unstable.
Definition ParticleDecay.h:41
thermalfist::FitParameter
Structure for an arbitrary fit parameter.
Definition ThermalModelFitParameters.h:20
thermalfist::FitParameter::toFit
bool toFit
Whether the parameter is fitted or fixed.
Definition ThermalModelFitParameters.h:28
thermalfist::FitParameter::error
double error
Parameter error (symmetric)
Definition ThermalModelFitParameters.h:25
thermalfist::FitParameter::name
std::string name
Name of the parameter.
Definition ThermalModelFitParameters.h:29
thermalfist::FittedQuantity
Structure describing the measurement to be fitted or compared to model.
Definition ThermalModelFitQuantities.h:122
thermalfist::FittedQuantity::toFit
bool toFit
Whether this quantity contributes to the of a fit.
Definition ThermalModelFitQuantities.h:133
thermalfist::FittedQuantity::Value
double Value() const
Value of the measurement.
Definition ThermalModelFitQuantities.h:165
thermalfist::FittedQuantity::ValueError
double ValueError() const
Error of the measurement.
Definition ThermalModelFitQuantities.h:173
thermalfist::ThermalModelFitParametersExtended
Extended structure for calculating uncertainties in non-fit quantities resulting from uncertanties in...
Definition ThermalModelFitParameters.h:43
thermalfist::ThermalModelFitParameters
Structure holding information about parameters of a thermal fit.
Definition ThermalModelFitParameters.h:69
thermalfist::ThermalModelFitParameters::GetThermalModelParameters
ThermalModelParameters GetThermalModelParameters()
Return ThermalModelParameters corresponding to the this ThermalModelFitParameters.
Definition ThermalModelFitParameters.cpp:182
thermalfist::ThermalModelFitParameters::Tkin
FitParameter Tkin
Since version 1.3: Kinetic freeze-out temperature.
Definition ThermalModelFitParameters.h:78
thermalfist::ThermalModelFitParameters::T
FitParameter T
All the fit parameters.
Definition ThermalModelFitParameters.h:75
thermalfist::ThermalModelFitParameters::ParameterList
std::vector< FitParameter * > ParameterList
Vector of pointer to all the parameters.
Definition ThermalModelFitParameters.h:81
thermalfist::ThermalModelFitParameters::GetParameter
FitParameter GetParameter(const std::string &name) const
Get the FitParameter by its name.
Definition ThermalModelFitParameters.cpp:120
thermalfist::ThermalModelParameters
Structure containing all thermal parameters of the model.
Definition ThermalModelParameters.h:18
Generated by doxygen 1.13.2