CbmRoot
CbmThermalModelNoFlow.cxx
Go to the documentation of this file.
1 /*
2  *====================================================================
3  *
4  * CBM Thermal Model, No Flow. This module is obsolete, to be completely replaced by CbmBoltzmannDistribution
5  *
6  * Authors: V.Vovchenko
7  *
8  * e-mail : new
9  *
10  *====================================================================
11  *
12  * Boltzmann distribution calculations
13  *
14  *====================================================================
15  */
16 
17 #include "CbmThermalModelNoFlow.h"
18 #include "CbmL1Def.h"
19 
20 
21 #include "CbmKFVertex.h"
22 #include "CbmStsTrack.h"
23 
24 #include "TClonesArray.h"
25 #include "TDirectory.h"
26 #include "TFile.h"
27 #include "TMath.h"
28 
29 #include "TH1F.h"
30 #include "TH2F.h"
31 #include "TSpline.h"
32 
33 #include "CbmVertex.h"
34 
35 #include "L1Field.h"
36 
37 #include "TStopwatch.h"
38 #include <algorithm>
39 #include <cmath>
40 #include <cstdio>
41 #include <iostream>
42 #include <map>
43 
44 #include "CbmMCTrack.h"
45 #include "CbmTrackMatch.h"
46 //for particle ID from global track
47 #include "CbmGlobalTrack.h"
48 #include "CbmRichRing.h"
49 #include "CbmTofHit.h"
50 #include "CbmTrdTrack.h"
51 #include "TDatabasePDG.h"
52 //for RICH identification
53 #include "TSystem.h"
54 // #include "CbmRichElectronIdAnn.h"
55 
56 #include "CbmL1PFFitter.h"
57 
58 
59 using std::ios;
60 using std::vector;
61 
62 namespace ThermalModelNoFlowNamespace {
63 
64 #include "NumericalIntegration.h"
65 
66  const Double_t kProtonMass = 0.938271998;
67  const int p_sz1 = 2;
68  //const int p_sz = 4;
69  const TString p_names[p_sz1] = {"pi-", "pi+"}; //, "K-", "K+"};//, "p"};
70  const TString p_names_gr[p_sz1] = {
71  "#pi^{-}",
72  "#pi^{+}"}; //, "K^{-}", "K^{+}"};//, "p"};
73  const int pdgIds[p_sz1] = {-211, 211}; //, -321, 321};//, 2212};
74  const Double_t GeVtoifm = 5.06423;
75 
76  const int recoLevels = 4;
77  const TString LevelNames[recoLevels] = {"Level 0",
78  "Level 1",
79  "Level 2",
80  "Level 3"};
81 
82  Double_t AcceptanceFunction::getAcceptance(const Double_t& y,
83  const Double_t& pt) const {
84  double ret = sfunc.Eval(y, pt);
85  if (ret < 0.) ret = 0.;
86  if (ret > 1.) ret = 1.;
87  return ret;
88  }
89 
90  class ThermalDistributionFunction {
91 
92  public:
93  ThermalDistributionFunction(int part,
94  double T_,
95  double R_,
96  double ekin_,
97  AcceptanceFunction* af_ = NULL,
98  ReconstructionEfficiencyFunction* rf_ = NULL)
99  : fT(T_)
100  , fV(4. / 3. * TMath::Pi() * R_ * R_ * R_)
101  , mass(0.)
102  , ekin(ekin_)
103  , ycm(0.)
104  , af(af_)
105  , rf(rf_) {
106  mass = TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass();
107  double v =
108  sqrt(2. * kProtonMass * ekin + ekin * ekin) / (2. * kProtonMass + ekin);
109  ycm = 0.5 * log((1 + v) / (1 - v)); //TMath::ATanH(v);
110  }
111 
112  ~ThermalDistributionFunction() {}
113 
114  double dndycm(double y) const {
115  double ret = 0.;
116  double dpt = 0.05;
117  double tmpmt = 0.;
118  double tmpr = 0.;
119  for (double pt = 0.5 * dpt; pt < 6.; pt += dpt) {
120  tmpmt = sqrt(pt * pt + mass * mass);
121  tmpr = pt * tmpmt * TMath::CosH(y)
122  * TMath::Exp(-tmpmt * TMath::CosH(y) / fT);
123  if (af != NULL) tmpr *= af->getAcceptance(y + ycm, pt);
124  if (rf != NULL) {
125  double tp = tmpmt * TMath::CosH(y + ycm);
126  tp = sqrt(tp * tp - mass * mass);
127  tmpr *= rf->f(tp);
128  }
129  ret += tmpr;
130  }
131  ret *= dpt * GeVtoifm * GeVtoifm * GeVtoifm * fV / 4. / TMath::Pi()
132  / TMath::Pi();
133  return ret;
134  }
135 
136  double dndylab(double y) const {
137  double dpt = 0.025;
138  double tmpmt = 0.;
139  double tmpr = 0.;
140  double ret = 0.;
141  double ss2 = 0.5 * sqrt(2 * kProtonMass * (ekin + 2 * kProtonMass));
142  double v = sqrt(1. - kProtonMass * kProtonMass / ss2 / ss2);
143  y = y - TMath::ATanH(v);
144 
145  for (double pt = 0.5 * dpt; pt < 6.; pt += dpt) {
146  tmpmt = sqrt(pt * pt + mass * mass);
147  tmpr = pt * tmpmt * TMath::CosH(y)
148  * TMath::Exp(-tmpmt * TMath::CosH(y) / fT);
149  if (af != NULL) tmpr *= af->getAcceptance(y + ycm, pt);
150  if (rf != NULL) {
151  double tp = tmpmt * TMath::CosH(y + ycm);
152  tp = sqrt(tp * tp - mass * mass);
153  tmpr *= rf->f(tp);
154  }
155  ret += tmpr;
156  }
157  ret *= dpt * GeVtoifm * GeVtoifm * GeVtoifm * fV / 4. / TMath::Pi()
158  / TMath::Pi();
159  return ret;
160  }
161 
162  double dndybinlab(double ymin, double ymax, int itery) const {
163  double dy = (ymax - ymin) / itery;
164  double ty = 0., ret = 0., ret2 = 0.;
165  for (int iy = 0; iy < itery; ++iy) {
166  ty = ymin + dy * iy + dy * 0.5;
167  ret2 += 1.;
168  ret += dndylab(ty);
169  }
170  return ret / ret2;
171  }
172 
173  double dndydptcm(double y, double pt) const {
174  double tmpmt = sqrt(pt * pt + mass * mass);
175  double ret = GeVtoifm * GeVtoifm * GeVtoifm * fV / 4. / TMath::Pi()
176  / TMath::Pi() * pt * tmpmt * TMath::CosH(y)
177  * TMath::Exp(-tmpmt * TMath::CosH(y) / fT);
178  if (af != NULL) ret *= af->getAcceptance(y + ycm, pt);
179  if (rf != NULL) {
180  double tp = tmpmt * TMath::CosH(y + ycm);
181  tp = sqrt(tp * tp - mass * mass);
182  ret *= rf->f(tp);
183  }
184  return ret;
185  }
186 
187  double dndydptlab(double y, double pt) const {
188  double tmpmt = sqrt(pt * pt + mass * mass);
189  double ss2 = 0.5 * sqrt(2 * kProtonMass * (ekin + 2 * kProtonMass));
190  double v = sqrt(1. - kProtonMass * kProtonMass / ss2 / ss2);
191  double ret = 1.;
192  y = y - TMath::ATanH(v);
193  ret *= GeVtoifm * GeVtoifm * GeVtoifm * fV / 4. / TMath::Pi()
194  / TMath::Pi() * pt * tmpmt * TMath::CosH(y)
195  * TMath::Exp(-tmpmt * TMath::CosH(y) / fT);
196  if (af != NULL) ret *= af->getAcceptance(y + ycm, pt);
197  if (rf != NULL) {
198  double tp = tmpmt * TMath::CosH(y + ycm);
199  tp = sqrt(tp * tp - mass * mass);
200  ret *= rf->f(tp);
201  }
202  return ret;
203  }
204 
205  double dndydptbin(double ymin,
206  double ymax,
207  int itery,
208  double ptmin,
209  double ptmax,
210  int iterpt) const {
211  double dpt = (ptmax - ptmin) / iterpt;
212  double dy = (ymax - ymin) / itery;
213  double ty = 0., tpt = 0., ret = 0., ret2 = 0.;
214  double tmpmt = 0.;
215  double tmpr = 0.;
216  double ss2 = 0.5 * sqrt(2 * kProtonMass * (ekin + 2 * kProtonMass));
217  double v = sqrt(1. - kProtonMass * kProtonMass / ss2 / ss2);
218  ymin = ymin - TMath::ATanH(v);
219  for (int ipt = 0; ipt < iterpt; ++ipt) {
220  tpt = ptmin + dpt * ipt + dpt * 0.5;
221  tmpmt = sqrt(tpt * tpt + mass * mass);
222  for (int iy = 0; iy < itery; ++iy) {
223  ty = ymin + dy * iy + dy * 0.5;
224  tmpr = GeVtoifm * GeVtoifm * GeVtoifm * fV / 4. / TMath::Pi()
225  / TMath::Pi() * tpt * tmpmt * TMath::CosH(ty)
226  * TMath::Exp(-tmpmt * TMath::CosH(ty) / fT);
227  ret2 += 1.;
228  if (af != NULL) tmpr *= af->getAcceptance(ty + ycm, tpt);
229  if (rf != NULL) {
230  double tp = tmpmt * TMath::CosH(ty + ycm);
231  tp = sqrt(tp * tp - mass * mass);
232  tmpr *= rf->f(tp);
233  }
234  ret += tmpr;
235  }
236  }
237  return ret / ret2;
238  }
239 
240  static double chi2func(double /*y*/, double pt, double m) {
241  return sqrt(pt * pt + m * m); //y*y;//pt*pt + m*m;
242  }
243 
244  // Function to determine temperature
245  static double tempCrit(double /*y*/, double pt, double m) {
246  return sqrt(pt * pt + m * m); //*cosh(y);
247  }
248 
249  double tempCritAv() {
250  if (af == NULL) {
251  double ret1 = 0., ret2 = 0.;
252  double ymin = -3., ymax = 6.;
253  int itery = 500;
254  double dy = (ymax - ymin) / itery;
255  double ptmin = 0., ptmax = 3.;
256  int iterpt = 400;
257  double dpt = (ptmax - ptmin) / iterpt;
258  double ty = 0., tpt = 0., tmp = 0.;
259  double tp = 0.;
260  for (int iy = 0; iy < itery; ++iy) {
261  ty = ymin + (iy + 0.5) * dy;
262  for (int ipt = 0; ipt < iterpt; ++ipt) {
263  tpt = ptmin + (ipt + 0.5) * dpt;
264  tmp = tpt * sqrt(tpt * tpt + mass * mass) * TMath::CosH(ty - ycm)
265  * TMath::Exp(-sqrt(tpt * tpt + mass * mass)
266  * TMath::CosH(ty - ycm) / fT);
267  ret1 += tmp;
268  ret2 += tmp * tempCrit(ty, tpt, mass);
269  if (rf != NULL) {
270  tp = sqrt(mass * mass * TMath::SinH(ty) * TMath::SinH(ty)
271  + tpt * tpt * TMath::CosH(ty) * TMath::CosH(ty));
272  ret1 *= rf->f(tp);
273  ret2 *= rf->f(tp);
274  }
275  }
276  }
277  return ret2 / ret1;
278  } else {
279  if (af->ys.size() == 0) return -1.;
280  double ret1 = 0., ret2 = 0.;
281  double tmp = 0., tp = 0.;
282  for (unsigned int i = 0; i < af->ys.size(); ++i) {
283  tp = sqrt(af->pts[i] * af->pts[i] + mass * mass) * cosh(af->ys[i]);
284  tp = sqrt(tp * tp - mass * mass);
285  tmp = af->pts[i] * sqrt(af->pts[i] * af->pts[i] + mass * mass)
286  * cosh(af->ys[i] - ycm)
287  * exp(-sqrt(af->pts[i] * af->pts[i] + mass * mass)
288  * cosh(af->ys[i] - ycm) / fT)
289  * af->probs[i];
290  if (rf != NULL) tmp *= rf->f(tp);
291  ret1 += tmp;
292  ret2 += tmp * tempCrit(af->ys[i], af->pts[i], mass);
293  }
294  return ret2 / ret1;
295  }
296  }
297 
298  double chi2FuncAv() {
299  if (af == NULL) {
300  double ret1 = 0., ret2 = 0.;
301  double ymin = -3., ymax = 6.;
302  int itery = 500;
303  double dy = (ymax - ymin) / itery;
304  double ptmin = 0., ptmax = 3.;
305  int iterpt = 400;
306  double dpt = (ptmax - ptmin) / iterpt;
307  double ty = 0., tpt = 0., tmp = 0.;
308  double tp = 0.;
309  for (int iy = 0; iy < itery; ++iy) {
310  ty = ymin + (iy + 0.5) * dy;
311  for (int ipt = 0; ipt < iterpt; ++ipt) {
312  tpt = ptmin + (ipt + 0.5) * dpt;
313  tmp = tpt * sqrt(tpt * tpt + mass * mass) * TMath::CosH(ty - ycm)
314  * TMath::Exp(-sqrt(tpt * tpt + mass * mass)
315  * TMath::CosH(ty - ycm) / fT);
316  ret1 += tmp;
317  ret2 += tmp * chi2func(ty, tpt, mass);
318  if (rf != NULL) {
319  tp = sqrt(mass * mass * TMath::SinH(ty) * TMath::SinH(ty)
320  + tpt * tpt * TMath::CosH(ty) * TMath::CosH(ty));
321  ret1 *= rf->f(tp);
322  ret2 *= rf->f(tp);
323  }
324  }
325  }
326  return ret2 / ret1;
327  } else {
328  if (af->ys.size() == 0) return -1.;
329  double ret1 = 0., ret2 = 0.;
330  double tmp = 0., tp = 0.;
331  for (unsigned int i = 0; i < af->ys.size(); ++i) {
332  tp = sqrt(af->pts[i] * af->pts[i] + mass * mass) * cosh(af->ys[i]);
333  tp = sqrt(tp * tp - mass * mass);
334  tmp = af->pts[i] * sqrt(af->pts[i] * af->pts[i] + mass * mass)
335  * cosh(af->ys[i] - ycm)
336  * exp(-sqrt(af->pts[i] * af->pts[i] + mass * mass)
337  * cosh(af->ys[i] - ycm) / fT)
338  * af->probs[i];
339  if (rf != NULL) tmp *= rf->f(tp);
340  ret1 += tmp;
341  ret2 += tmp * chi2func(af->ys[i], af->pts[i], mass);
342  }
343  return ret2 / ret1;
344  }
345  }
346 
347  double T() const { return fT; }
348  double V() const { return fV; }
349 
350  ThermalDistributionFunction(const ThermalDistributionFunction&);
351  ThermalDistributionFunction& operator=(const ThermalDistributionFunction&);
352 
353  private:
354  double fT, fV;
355  double mass;
356  double ekin;
357  double ycm;
358  AcceptanceFunction* af;
359  ReconstructionEfficiencyFunction* rf;
360  };
361 
362  class ThermalChi2Func {
363 
364  public:
365  ThermalChi2Func(TH1F* dndyexp, TH2F* dndydptexp, double Norm_)
366  : Norm(Norm_), dndyHist(0), dndydptHist(0) {
367  dndyHist = dndyexp;
368  dndydptHist = dndydptexp;
369  }
370 
371  ~ThermalChi2Func() {}
372 
373  double chi2dndy(int part,
374  double T_,
375  double R_,
376  double ekin_,
377  AcceptanceFunction* af_ = NULL,
378  ReconstructionEfficiencyFunction* rf_ = NULL,
379  double systerr = 0.) const {
380  ThermalDistributionFunction pl(part, T_, R_, ekin_, af_, rf_);
381  double chi2 = 0.;
382  double tmpval = 0., tmpval2 = 0., tmpvar = 0.;
383 
384  int iters = 0;
385 
386 
387  for (int n = 0; n < dndyHist->GetNbinsX(); n++) {
388  tmpval = dndyHist->GetBinContent(n);
389  if (tmpval > 100.) {
390  iters++;
391  tmpval /= Norm * dndyHist->GetXaxis()->GetBinWidth(n);
392  tmpvar = dndyHist->GetBinError(n) / Norm
393  / dndyHist->GetXaxis()->GetBinWidth(n);
394  tmpvar = sqrt(tmpvar * tmpvar + systerr * systerr * tmpval * tmpval);
395  tmpval2 = pl.dndybinlab(dndyHist->GetXaxis()->GetBinLowEdge(n),
396  dndyHist->GetXaxis()->GetBinUpEdge(n),
397  10);
398  chi2 += (tmpval - tmpval2) * (tmpval - tmpval2) / tmpvar
399  / tmpvar;
400  }
401  }
402 
403  chi2 /= iters;
404 
405  return chi2;
406  }
407 
408  double chi2ypt(int part,
409  double T_,
410  double R_,
411  double ekin_,
412  AcceptanceFunction* af_ = NULL,
413  ReconstructionEfficiencyFunction* rf_ = NULL,
414  double systerr = 0.) const {
415  ThermalDistributionFunction pl(part, T_, R_, ekin_, af_, rf_);
416  double chi2 = 0.;
417  int iters = 0;
418  double tmpval = 0., tmpval2 = 0., tmpvar = 0.;
419  for (int nx = 0; nx < dndydptHist->GetNbinsX(); nx++) {
420  for (int ny = 0; ny < dndydptHist->GetNbinsY(); ny++) {
421  tmpval = dndydptHist->GetBinContent(nx, ny);
422  if (tmpval > 5.) {
423  iters++;
424  tmpval /= Norm * dndydptHist->GetXaxis()->GetBinWidth(nx)
425  * dndydptHist->GetYaxis()->GetBinWidth(ny);
426  tmpvar = dndydptHist->GetBinError(nx, ny) / Norm
427  / dndydptHist->GetXaxis()->GetBinWidth(nx)
428  / dndydptHist->GetYaxis()->GetBinWidth(ny);
429  tmpvar =
430  sqrt(tmpvar * tmpvar + systerr * systerr * tmpval * tmpval);
431  tmpval2 =
432  pl.dndydptbin(dndydptHist->GetXaxis()->GetBinLowEdge(nx),
433  dndydptHist->GetXaxis()->GetBinUpEdge(nx),
434  10,
435  dndydptHist->GetYaxis()->GetBinLowEdge(ny),
436  dndydptHist->GetYaxis()->GetBinUpEdge(ny),
437  10); //pl.dndyexp(dndyHist->GetBinCenter(n));
438  chi2 += (tmpval - tmpval2) * (tmpval - tmpval2) / tmpvar
439  / tmpvar;
440  }
441  }
442  }
443  return chi2 / iters;
444  }
445 
446  ThermalChi2Func(const ThermalChi2Func&);
447  ThermalChi2Func& operator=(const ThermalChi2Func&);
448 
449  private:
450  double Norm;
451  TH1F* dndyHist;
452  TH2F* dndydptHist;
453  };
454 } // namespace ThermalModelNoFlowNamespace
455 
456 
457 using namespace ThermalModelNoFlowNamespace;
458 using namespace std;
459 
460 ClassImp(CbmThermalModelNoFlow)
461 
462  CbmThermalModelNoFlow::CbmThermalModelNoFlow(Float_t ekin_,
463  Int_t recoLevel,
464  Int_t usePID,
465  Int_t trackNumber,
466  Int_t /*iVerbose*/)
467  : ekin(ekin_)
468  , ycm(1.)
469  , fUpdate(true)
470  , fusePID(usePID)
471  , fRecoLevel(recoLevel)
472  , fTrackNumber(trackNumber)
473  , flistStsPts(0)
474  , flistTofPts(0)
475  , flistStsTracksMatch(0)
476  , flistStsTracks(0)
477  , fPrimVtx(0)
478  , outfileName("")
479  , fChiToPrimVtx()
480  , NPrimGlobalMC(0)
481  , NPrimGlobalReco(0)
482  , flistMCTracks(0)
483  , flsitGlobalTracks(0)
484  , flistTofHits(0)
485  , flistTofPoints(0)
486  , histodir(0)
487  , events(0)
488  , AcceptanceSTS()
489  , AcceptanceSTSTOF() {
490  NPrimGlobalMC = 0;
491  NPrimGlobalReco = 0;
492 
493  double pbeam = sqrt((kProtonMass + ekin) * (kProtonMass + ekin)
494  - kProtonMass * kProtonMass);
495  double betacm = pbeam / (2. * kProtonMass + ekin);
496  ycm = 0.5 * log((1. + betacm) / (1. - betacm));
497 
498  for (int part = 0; part < p_sz; ++part)
499  ComputeThermalDependence(part);
500 
501  /*globalmtavMC = new Double_t[p_sz];
502  globalmt2avMC = new Double_t[p_sz];
503  globalfavMC = new Double_t[p_sz];
504  globalf2avMC = new Double_t[p_sz];
505  globalmtmomerrMC = new Double_t[p_sz];
506  globalnTracksMC = new Int_t[p_sz];*/
507  for (int part = 0; part < p_sz; ++part) {
508  globalmtavMC[part] = 0.;
509  globalmt2avMC[part] = 0.;
510  globalfavMC[part] = 0.;
511  globalf2avMC[part] = 0.;
512  globalmtmomerrMC[part] = 0.;
513  globalnTracksMC[part] = 0;
514  }
515 
516  //for(Int_t i=0;i<recoLevels;++i) {
517  /*globalmtavReco = new Double_t[p_sz];
518  globalmt2avReco = new Double_t[p_sz];
519  globalfavReco = new Double_t[p_sz];
520  globalf2avReco = new Double_t[p_sz];
521  globalnTracksReco = new Int_t[p_sz];
522  globalmtmomerrReco = new Double_t[p_sz];*/
523  for (int part = 0; part < p_sz; ++part) {
524  globalmtavReco[part] = 0.;
525  globalmt2avReco[part] = 0.;
526  globalfavReco[part] = 0.;
527  globalf2avReco[part] = 0.;
528  //globalmt4avReco[part] = 0.;
529  globalmtmomerrReco[part] = 0.;
530  globalnTracksReco[part] = 0;
531  }
532  //}
533 
534 
535  TDirectory* currentDir = gDirectory;
536 
537  gDirectory->cd("Models");
538 
539  histodir = gDirectory;
540 
541  gDirectory->mkdir("BoltzmannDistribution");
542  gDirectory->cd("BoltzmannDistribution");
543  gDirectory->mkdir("All tracks");
544  gDirectory->cd("All tracks");
545  gDirectory->mkdir("PerEvent");
546  gDirectory->cd("PerEvent");
547  /*hTempFullMC = new TH1F*[1];
548  hTempErrStatFullMC = new TH1F*[1];
549  hTempErrMomFullMC = new TH1F*[1];
550  hRadFullMC = new TH1F*[1];
551  hRadErrStatFullMC = new TH1F*[1];
552  hRadErrMomFullMC = new TH1F*[1];*/
553  for (int part = 0; part < p_sz; ++part) {
554  gDirectory->mkdir(p_names[part]);
555  gDirectory->cd(p_names[part]);
556  hTempFullMC[part] = new TH1F("Temperature " + p_names_gr[part],
557  "Event-by-event temperature",
558  300,
559  0.,
560  0.3);
561  hTempFullMC[part]->SetXTitle("T (GeV)");
562  hTempFullMC[part]->SetYTitle("Entries");
563 
564  hTempErrStatFullMC[part] =
565  new TH1F("Temperature statistical error" + p_names_gr[part],
566  "Event-by-event temperature statistical error",
567  100,
568  0.,
569  0.1);
570  hTempErrStatFullMC[part]->SetXTitle("T (GeV)");
571  hTempErrStatFullMC[part]->SetYTitle("Entries");
572 
573  hTempErrMomFullMC[part] =
574  new TH1F("Temperature momentum error" + p_names_gr[part],
575  "Event-by-event temperature momentum error",
576  100,
577  0.,
578  0.01);
579  hTempErrMomFullMC[part]->SetXTitle("T (GeV)");
580  hTempErrMomFullMC[part]->SetYTitle("Entries");
581 
582  hRadFullMC[part] = new TH1F("Fireball radius " + p_names_gr[part],
583  "Event-by-event radius",
584  200,
585  0.,
586  100.);
587  hRadFullMC[part]->SetXTitle("R (fm)");
588  hRadFullMC[part]->SetYTitle("Entries");
589 
590  hRadErrStatFullMC[part] =
591  new TH1F("Fireball radius statistical error " + p_names_gr[part],
592  "Event-by-event radius statistical error",
593  100,
594  0.,
595  10.);
596  hRadErrStatFullMC[part]->SetXTitle("R (fm)");
597  hRadErrStatFullMC[part]->SetYTitle("Entries");
598 
599  hRadErrMomFullMC[part] =
600  new TH1F("Fireball radius momentum error " + p_names_gr[part],
601  "Event-by-event radius momentum error",
602  100,
603  0.,
604  10.);
605  hRadErrMomFullMC[part]->SetXTitle("R (fm)");
606  hRadErrMomFullMC[part]->SetYTitle("Entries");
607 
608  gDirectory->cd("..");
609  }
610  gDirectory->cd("..");
611  gDirectory->cd("..");
612  gDirectory->cd("..");
613 
614  gDirectory->cd("BoltzmannDistribution");
615  gDirectory->cd("All tracks");
616  gDirectory->mkdir("Global");
617  gDirectory->cd("Global");
618  /*hfyMC = new TH1F*[1];
619  hfyMCmodel = new TH1F*[1];
620  hfdndydptMC = new TH2F*[1];
621  hfdndydptMCmodel = new TH2F*[1];*/
622  for (int part = 0; part < p_sz; ++part) {
623  gDirectory->mkdir(p_names[part]);
624  gDirectory->cd(p_names[part]);
625 
626  hfyMC[part] = new TH1F(
627  "dndy " + p_names_gr[part], "Rapidity distribution", 100, -5., 5.);
628  hfyMCmodel[part] = new TH1F("dndy " + p_names_gr[part] + " model",
629  "Rapidity distribution",
630  100,
631  -5.,
632  5.);
633  hfdndydptMC[part] = new TH2F("dndydpt " + p_names_gr[part],
634  "YPt distribution",
635  50,
636  -5.,
637  5.,
638  25,
639  0.,
640  5.);
641  hfdndydptMCmodel[part] = new TH2F("dndydpt " + p_names_gr[part] + " model",
642  "YPt distribution",
643  50,
644  -5.,
645  5.,
646  25,
647  0.,
648  5.);
649 
650  gDirectory->cd("..");
651  }
652  gDirectory->cd("..");
653  gDirectory->cd("..");
654  gDirectory->cd("..");
655 
656  gDirectory->cd("BoltzmannDistribution");
657  gDirectory->mkdir("Reconstructible tracks");
658  gDirectory->cd("Reconstructible tracks");
659  gDirectory->mkdir("PerEvent");
660  gDirectory->cd("PerEvent");
661 
662  for (int part = 0; part < p_sz; ++part) {
663  gDirectory->mkdir(p_names[part]);
664  gDirectory->cd(p_names[part]);
665  hTempFullReco[part] = new TH1F("Temperature " + p_names_gr[part],
666  "Event-by-event temperature",
667  300,
668  0.,
669  0.3);
670  hTempFullReco[part]->SetXTitle("T (GeV)");
671  hTempFullReco[part]->SetYTitle("Entries");
672 
673  hTempErrStatFullReco[part] =
674  new TH1F("Temperature statistical error " + p_names_gr[part],
675  "Event-by-event temperature statistical error ",
676  100,
677  0.,
678  0.1);
679  hTempErrStatFullReco[part]->SetXTitle("T (GeV)");
680  hTempErrStatFullReco[part]->SetYTitle("Entries");
681 
682  hTempErrMomFullReco[part] =
683  new TH1F("Temperature momentum error " + p_names_gr[part],
684  "Event-by-event temperature momentum error",
685  100,
686  0.,
687  0.01);
688  hTempErrMomFullReco[part]->SetXTitle("T (GeV)");
689  hTempErrMomFullReco[part]->SetYTitle("Entries");
690 
691 
692  hRadFullReco[part] = new TH1F("Fireball radius " + p_names_gr[part],
693  "Event-by-event radius",
694  200,
695  0.,
696  100.);
697  hRadFullReco[part]->SetXTitle("R (fm)");
698  hRadFullReco[part]->SetYTitle("Entries");
699 
700  hRadErrStatFullReco[part] =
701  new TH1F("Fireball radius statistical error " + p_names_gr[part],
702  "Event-by-event radius statistical error ",
703  10,
704  0.,
705  10.);
706  hRadErrStatFullReco[part]->SetXTitle("R (fm)");
707  hRadErrStatFullReco[part]->SetYTitle("Entries");
708 
709  hRadErrMomFullReco[part] =
710  new TH1F("Fireball radius momentum error " + p_names_gr[part],
711  "Event-by-event radius momentum error",
712  100,
713  0.,
714  10.);
715  hRadErrMomFullReco[part]->SetXTitle("R (fm)");
716  hRadErrMomFullReco[part]->SetYTitle("Entries");
717 
718  gDirectory->cd("..");
719  }
720 
721  gDirectory->cd("..");
722  gDirectory->cd("..");
723  gDirectory->cd("..");
724 
725  gDirectory->cd("BoltzmannDistribution");
726  gDirectory->cd("Reconstructible tracks");
727  gDirectory->mkdir("Global");
728  gDirectory->cd("Global");
729  for (int part = 0; part < p_sz; ++part) {
730  gDirectory->mkdir(p_names[part]);
731  gDirectory->cd(p_names[part]);
732 
733  hfyReco[part] = new TH1F(
734  "dndy " + p_names_gr[part], "Rapidity distribution", 100, -5., 5.);
735  hfyRecomodel[part] = new TH1F("dndy " + p_names_gr[part] + " model",
736  "Rapidity distribution",
737  100,
738  -5.,
739  5.);
740  hfdndydptReco[part] = new TH2F("dndydpt " + p_names_gr[part],
741  "YPt distribution",
742  50,
743  -5.,
744  5.,
745  25,
746  0.,
747  5.);
748  hfdndydptRecomodel[part] =
749  new TH2F("dndydpt " + p_names_gr[part] + " model",
750  "YPt distribution",
751  50,
752  -5.,
753  5.,
754  25,
755  0.,
756  5.);
757 
758  gDirectory->cd("..");
759  }
760  gDirectory->cd("..");
761  gDirectory->cd("..");
762  gDirectory->cd("..");
763 
764  gDirectory->cd("BoltzmannDistribution");
765  gDirectory->mkdir("Reconstructible tracks corrected");
766  gDirectory->cd("Reconstructible tracks corrected");
767  gDirectory->mkdir("PerEvent");
768  gDirectory->cd("PerEvent");
769  for (int part = 0; part < p_sz; ++part) {
770  gDirectory->mkdir(p_names[part]);
771  gDirectory->cd(p_names[part]);
772  hTempFullRecoCor[part] = new TH1F("Temperature " + p_names_gr[part],
773  "Event-by-event temperature",
774  300,
775  0.,
776  0.3);
777  hTempFullRecoCor[part]->SetXTitle("T (GeV)");
778  hTempFullRecoCor[part]->SetYTitle("Entries");
779 
780  hTempErrStatFullRecoCor[part] =
781  new TH1F("Temperature statistical error " + p_names_gr[part],
782  "Event-by-event temperature statistical error ",
783  100,
784  0.,
785  0.1);
786  hTempErrStatFullRecoCor[part]->SetXTitle("T (GeV)");
787  hTempErrStatFullRecoCor[part]->SetYTitle("Entries");
788 
789  hTempErrMomFullRecoCor[part] =
790  new TH1F("Temperature momentum error " + p_names_gr[part],
791  "Event-by-event temperature momentum error",
792  100,
793  0.,
794  0.01);
795  hTempErrMomFullRecoCor[part]->SetXTitle("T (GeV)");
796  hTempErrMomFullRecoCor[part]->SetYTitle("Entries");
797 
798  hRadFullRecoCor[part] = new TH1F("Fireball radius " + p_names_gr[part],
799  "Event-by-event radius",
800  200,
801  0.,
802  100.);
803  hRadFullRecoCor[part]->SetXTitle("R (fm)");
804  hRadFullRecoCor[part]->SetYTitle("Entries");
805 
806  hRadErrStatFullRecoCor[part] =
807  new TH1F("Fireball radius statistical error " + p_names_gr[part],
808  "Event-by-event radius statistical error ",
809  100,
810  0.,
811  10.);
812  hRadErrStatFullRecoCor[part]->SetXTitle("R (fm)");
813  hRadErrStatFullRecoCor[part]->SetYTitle("Entries");
814 
815  hRadErrMomFullRecoCor[part] =
816  new TH1F("Fireball radius momentum error " + p_names_gr[part],
817  "Event-by-event radius momentum error",
818  100,
819  0.,
820  10.);
821  hRadErrMomFullRecoCor[part]->SetXTitle("R (fm)");
822  hRadErrMomFullRecoCor[part]->SetYTitle("Entries");
823 
824  gDirectory->cd("..");
825  }
826  gDirectory->cd("..");
827  gDirectory->cd("..");
828  gDirectory->cd("..");
829 
830  gDirectory->cd("BoltzmannDistribution");
831  gDirectory->cd("Reconstructible tracks corrected");
832  gDirectory->mkdir("Global");
833  gDirectory->cd("Global");
834  for (int part = 0; part < p_sz; ++part) {
835  gDirectory->mkdir(p_names[part]);
836  gDirectory->cd(p_names[part]);
837 
838  hfyRecoCor[part] = new TH1F(
839  "dndy " + p_names_gr[part], "Rapidity distribution", 100, -5., 5.);
840  hfyRecoCormodel[part] = new TH1F("dndy " + p_names_gr[part] + " model",
841  "Rapidity distribution",
842  100,
843  -5.,
844  5.);
845  hfdndydptRecoCor[part] = new TH2F("dndydpt " + p_names_gr[part],
846  "YPt distribution",
847  50,
848  -5.,
849  5.,
850  25,
851  0.,
852  5.);
853  hfdndydptRecoCormodel[part] =
854  new TH2F("dndydpt " + p_names_gr[part] + " model",
855  "YPt distribution",
856  50,
857  -5.,
858  5.,
859  25,
860  0.,
861  5.);
862 
863  gDirectory->cd("..");
864  }
865  gDirectory->cd("..");
866  gDirectory->cd("..");
867  gDirectory->cd("..");
868 
869  gDirectory = currentDir;
870  //cout << "ycm = " << ycm;
871 
872  events = 0;
873 }
874 
875 CbmThermalModelNoFlow::~CbmThermalModelNoFlow() {}
876 
877 void CbmThermalModelNoFlow::ReadAcceptanceFunction(AcceptanceFunction& func,
878  TString filename) {
879  double ymin = 0., ymax = 6.;
880  double ptmin = 0., ptmax = 2.5;
881  func.ys.resize(0);
882  func.pts.resize(0);
883  func.probs.resize(0);
884  ifstream fin(filename.Data());
885  fin >> func.dy >> func.dpt;
886  double ty, tpt, prob;
887  func.ys.resize(0);
888  func.pts.resize(0);
889  func.probs.resize(0);
890  while (fin >> ty >> tpt >> prob) {
891  if (tpt < ptmin || tpt > ptmax || ty < ymin || ty > ymax) continue;
892  func.ys.push_back(ty);
893  func.pts.push_back(tpt);
894  func.probs.push_back(prob);
895  }
896  func.setSpline();
897  fin.close();
898 }
899 
900 void CbmThermalModelNoFlow::ReInit(FairRootManager* fManger) {
901  //FairRootManager *fManger = FairRootManager::Instance();
902 
903  // flistStsTracks = (TClonesArray *) fManger->GetObject("StsTrack");
904  // flistTofPts = (TClonesArray *) fManger->GetObject("TofPoint");
905  // fPrimVtx = (CbmVertex*) fManger->GetObject("PrimaryVertex");
906  flistStsTracks = dynamic_cast<TClonesArray*>(fManger->GetObject("StsTrack"));
907  flistTofPts = dynamic_cast<TClonesArray*>(fManger->GetObject("TofPoint"));
908 
909  // Get pointer to PrimaryVertex object from IOManager if it exists
910  // The old name for the object is "PrimaryVertex" the new one
911  // "PrimaryVertex." Check first for the new name
912  fPrimVtx = dynamic_cast<CbmVertex*>(fManger->GetObject("PrimaryVertex."));
913  if (nullptr == fPrimVtx) {
914  fPrimVtx = dynamic_cast<CbmVertex*>(fManger->GetObject("PrimaryVertex"));
915  }
916  if (nullptr == fPrimVtx) {
917  Error("CbmThermalModelNoFlow::ReInit", "vertex not found!");
918  }
919  //fPrimVtx = dynamic_cast<CbmVertex*>( fManger->GetObject("PrimaryVertex") );
920  //fPrimVtx = new CbmVertex();
921  //for the particle id
922  flistStsTracksMatch =
923  dynamic_cast<TClonesArray*>(fManger->GetObject("StsTrackMatch"));
924  flistMCTracks = dynamic_cast<TClonesArray*>(fManger->GetObject("MCTrack"));
925  // flistStsPts = (TClonesArray *) fManger->GetObject("StsPoint");
926  // flistTofPts = (TClonesArray *) fManger->GetObject("TofPoint");
927  flistStsPts = dynamic_cast<TClonesArray*>(fManger->GetObject("StsPoint"));
928  flistTofPts = dynamic_cast<TClonesArray*>(fManger->GetObject("TofPoint"));
929 
930 
931  if (fusePID == 2) {
932  flsitGlobalTracks =
933  dynamic_cast<TClonesArray*>(fManger->GetObject("GlobalTrack"));
934  flistTofHits = dynamic_cast<TClonesArray*>(fManger->GetObject("TofHit"));
935  }
936 }
937 
938 void CbmThermalModelNoFlow::Init() {
939  //ReInit();
940 }
941 
942 void CbmThermalModelNoFlow::Exec() {
943  if (!flistStsTracks) return;
944 
945  if (!flistStsTracksMatch) return;
946  if (!flistStsPts) return;
947  if (!flistMCTracks) return;
948 
949  vector<CbmStsTrack> vRTracks;
950  vector<CbmMCTrack> vRTracksMC;
951  int nTracks = flistStsTracks->GetEntries();
952  int nTracksMC = flistMCTracks->GetEntries();
953  vRTracks.resize(nTracks);
954  for (int iTr = 0; iTr < nTracks; iTr++)
955  // vRTracks[iTr] = *( (CbmStsTrack*) flistStsTracks->At(iTr));
956  vRTracks[iTr] = *(static_cast<CbmStsTrack*>(flistStsTracks->At(iTr)));
957  vRTracksMC.resize(nTracksMC);
958  for (int iTr = 0; iTr < nTracksMC; iTr++)
959  // vRTracksMC[iTr] = *( (CbmMCTrack*) flistMCTracks->At(iTr));
960  vRTracksMC[iTr] = *(static_cast<CbmMCTrack*>(flistMCTracks->At(iTr)));
961 
962  CbmKFVertex kfVertex;
963  if (fPrimVtx) kfVertex = CbmKFVertex(*fPrimVtx);
964 
965  vector<int> vTrackPDG(vRTracks.size(), -1);
966  if (fusePID == 1) {
967  for (int iTr = 0; iTr < nTracks; iTr++) {
968  // CbmTrackMatch* stsTrackMatch = (CbmTrackMatch*)flistStsTracksMatch->At(iTr);
969  CbmTrackMatch* stsTrackMatch =
970  static_cast<CbmTrackMatch*>(flistStsTracksMatch->At(iTr));
971  if (stsTrackMatch->GetNofMCTracks() == 0) continue;
972  const int mcTrackId = stsTrackMatch->GetMCTrackId();
973  // CbmMCTrack* mcTrack = (CbmMCTrack*)flistMCTracks->At(mcTrackId);
974  CbmMCTrack* mcTrack =
975  static_cast<CbmMCTrack*>(flistMCTracks->At(mcTrackId));
976  vTrackPDG[iTr] = mcTrack->GetPdgCode();
977  }
978  }
979 
980  if (fusePID == 2) {
981  if (NULL == flsitGlobalTracks) {
982  Fatal("KF Particle Finder", "No GlobalTrack array!");
983  }
984  if (NULL == flistTofHits) {
985  Fatal("KF Particle Finder", "No TOF hits array!");
986  }
987 
988  const Double_t m2P = 0.885;
989  const Double_t m2K = 0.245;
990  const Double_t m2Pi = 0.019479835;
991 
992  Double_t sP[3][5] = {
993  {0.0337428, -0.013939, 0.00567602, -0.000202229, 4.07531e-06},
994  {0.00717827, -0.00257353, 0.00389851, -9.83097e-05, 1.33011e-06},
995  {0.001348, 0.00220126, 0.0023619, 7.35395e-05, -4.06706e-06}}; //SIS-300
996  // Double_t sP[3][5] = { {0.056908,-0.0470572,0.0216465,-0.0021016,8.50396e-05},
997  // {0.00943075,-0.00635429,0.00998695,-0.00111527,7.77811e-05},
998  // {0.00176298,0.00367263,0.00308013,0.000844013,-0.00010423} }; //SIS-100
999 
1000  // sP[0][0] = 0.0618927;
1001  // sP[0][1] = -0.0719277;
1002  // sP[0][2] = 0.0396255;
1003  // sP[0][3] = -0.00583356;
1004  // sP[0][4] = 0.000317072;
1005  //
1006  // sP[1][0] = 0.0291881;
1007  // sP[1][1] = -0.0904978;
1008  // sP[1][2] = 0.086161;
1009  // sP[1][3] = -0.0229688;
1010  // sP[1][4] = 0.00199382;
1011  //
1012  // sP[2][0] = 0.162171;
1013  // sP[2][1] = -0.194007;
1014  // sP[2][2] = 0.0893264;
1015  // sP[2][3] = -0.014626;
1016  // sP[2][4] = 0.00088203;
1017 
1018  const Int_t PdgHypo[4] = {2212, 321, 211, -11};
1019 
1020  for (Int_t igt = 0; igt < flsitGlobalTracks->GetEntriesFast(); igt++) {
1021  const CbmGlobalTrack* globalTrack =
1022  static_cast<const CbmGlobalTrack*>(flsitGlobalTracks->At(igt));
1023 
1024  Int_t stsTrackIndex = globalTrack->GetStsTrackIndex();
1025  if (stsTrackIndex < 0) continue;
1026 
1027  // Bool_t isElectronTRD = 0;
1028  // Bool_t isElectronRICH = 0;
1029  // Bool_t isElectron = 0;
1030 
1031  const FairTrackParam* stsPar = vRTracks[stsTrackIndex].GetParamFirst();
1032  TVector3 mom;
1033  stsPar->Momentum(mom);
1034 
1035  Double_t p = mom.Mag();
1036  Int_t q = stsPar->GetQp() > 0 ? 1 : -1;
1037 
1038  Double_t l = globalTrack->GetLength();
1039  if (!((l > 800.) && (l < 1400.))) continue; //SIS 300
1040  // if( !((l>580.) && (l<700.)) ) continue;//SIS 100
1041 
1042  Double_t time;
1043  Int_t tofHitIndex = globalTrack->GetTofHitIndex();
1044  if (tofHitIndex >= 0) {
1045  vTrackPDG[stsTrackIndex] = -211;
1046  continue;
1047  const CbmTofHit* tofHit =
1048  static_cast<const CbmTofHit*>(flistTofHits->At(tofHitIndex));
1049  if (!tofHit) continue;
1050  // time = tofHit->GetTime()-0.021;
1051  time = tofHit->GetTime();
1052  } else
1053  continue;
1054 
1055  if (!((time > 29.) && (time < 50.))) continue; //SIS 300
1056  // if( !((time>19.) && (time<42.)) ) continue; //SIS 100
1057 
1058  Double_t m2 =
1059  p * p * (1. / ((l / time / 29.9792458) * (l / time / 29.9792458)) - 1.);
1060 
1061  Double_t sigma[3];
1062  sigma[0] = sP[0][0] + sP[0][1] * p + sP[0][2] * p * p
1063  + sP[0][3] * p * p * p + sP[0][4] * p * p * p * p;
1064  sigma[1] = sP[1][0] + sP[1][1] * p + sP[1][2] * p * p
1065  + sP[1][3] * p * p * p + sP[1][4] * p * p * p * p;
1066  sigma[2] = sP[2][0] + sP[2][1] * p + sP[2][2] * p * p
1067  + sP[2][3] * p * p * p + sP[2][4] * p * p * p * p;
1068 
1069  Double_t dm2[3];
1070  dm2[0] = fabs(m2 - m2P) / sigma[0];
1071  dm2[1] = fabs(m2 - m2K) / sigma[1];
1072  dm2[2] = fabs(m2 - m2Pi) / sigma[2];
1073 
1074  int iPdg = 2;
1075  Double_t dm2min = dm2[2];
1076 
1077  // if(isElectronRICH && isElectronTRD)
1078  // {
1079  // if (p >= 1.) {
1080  // if (m2 < (0.01 + (p - 1.) * 0.09))
1081  // isElectron = 1;
1082  // }
1083  // else {
1084  // if (m2 < 0.0)
1085  // isElectron = 1;
1086  // }
1087  // }
1088  //
1089  // if(!isElectron)
1090  {
1091  if (p > 12.) continue;
1092  if (q > 0) {
1093  if (dm2[1] < dm2min) {
1094  iPdg = 1;
1095  dm2min = dm2[1];
1096  }
1097  if (dm2[0] < dm2min) {
1098  iPdg = 0;
1099  dm2min = dm2[0];
1100  }
1101 
1102  if (dm2min > 2) iPdg = -1;
1103  } else {
1104  if (dm2[1] < dm2min) {
1105  iPdg = 1;
1106  dm2min = dm2[1];
1107  }
1108  if ((dm2min > 3) && (dm2[0] < dm2min)) {
1109  iPdg = 0;
1110  dm2min = dm2[0];
1111  }
1112 
1113  if (dm2min > 2) iPdg = -1;
1114  }
1115 
1116  if (iPdg > -1) vTrackPDG[stsTrackIndex] = q * PdgHypo[iPdg];
1117  }
1118  // else
1119  // vTrackPDG[stsTrackIndex] = q*PdgHypo[3];
1120  }
1121  }
1122 
1123  static int NEv = 0;
1124  NEv++;
1125  events++;
1126 
1127  //cout << NEv << endl;
1128 
1129 
1130  CbmKFTrErrMCPoints* MCTrackSortedArray =
1131  new CbmKFTrErrMCPoints[flistMCTracks->GetEntriesFast() + 1];
1132  for (Int_t iSts = 0; iSts < flistStsPts->GetEntriesFast(); iSts++) {
1133  // CbmStsPoint* StsPoint = (CbmStsPoint*)flistStsPts->At(iSts);
1134  CbmStsPoint* StsPoint = static_cast<CbmStsPoint*>(flistStsPts->At(iSts));
1135  MCTrackSortedArray[StsPoint->GetTrackID()].StsArray.push_back(StsPoint);
1136  }
1137 
1138  for (Int_t iTof = 0; iTof < flistTofPts->GetEntriesFast(); iTof++) {
1139  // CbmTofPoint* TofPoint = (CbmTofPoint*)flistTofPts->At(iTof);
1140  CbmTofPoint* TofPoint = static_cast<CbmTofPoint*>(flistTofPts->At(iTof));
1141  MCTrackSortedArray[TofPoint->GetTrackID()].TofArray.push_back(TofPoint);
1142  }
1143 
1144  //int allrec;
1145  //allrec = 0;
1146  int* TrRecons = new int[flistMCTracks->GetEntriesFast() + 1];
1147  int iTr = 0;
1148 
1149  double mtavMC[p_sz];
1150  double mtavReco[p_sz];
1151  double mt2avMC[p_sz];
1152  double mt2avReco[p_sz];
1153  double mtmomerrMC[p_sz];
1154  double mtmomerrReco[p_sz];
1155  double tmpmt;
1156  int nTracksPaAllMC[p_sz];
1157  int nTracksReco[p_sz];
1158  for (int part = 0; part < p_sz; ++part) {
1159  mtavMC[part] = 0.;
1160  mt2avMC[part] = 0.;
1161  mtmomerrMC[part] = 0.;
1162  mtavReco[part] = 0.;
1163  mt2avReco[part] = 0.;
1164  mtmomerrReco[part] = 0.;
1165  nTracksPaAllMC[part] = 0;
1166  nTracksReco[part] = 0;
1167  }
1168 
1169  for (iTr = 0; iTr < nTracksMC; iTr++) {
1170  //TrRecons[iTr] = (int)(checkIfReconstructable(&MCTrackSortedArray[iTr]));
1171  TrRecons[iTr] =
1172  MCTrackSortedArray[iTr].IsReconstructable(&vRTracksMC[iTr], 4, 3, 0.);
1173  }
1174 
1175  //fParticles.clear();
1176  fChiToPrimVtx.clear();
1177 
1178  CbmL1PFFitter fitter;
1179 
1180  // fitter.Fit(vRTracks); //assumed, that the initial fit should be fixed and must return good results!!!
1181  vector<FairTrackParam> paramsinit;
1182  for (int iTrs = 0; iTrs < nTracks; iTrs++) {
1183  paramsinit.push_back(*(vRTracks[iTrs].GetParamFirst()));
1184  }
1185 
1186  vector<L1FieldRegion> vField;
1187  fitter.GetChiToVertex(vRTracks, vField, fChiToPrimVtx, kfVertex, 3000000);
1188 
1189  // Level 0
1190  {
1191  for (iTr = 0; iTr < nTracksMC; iTr++) {
1192  //iTr = ((CbmTrackMatch*)flistStsTracksMatch->At(iTrs))->GetMCTrackId();
1193  for (int part = 0; part < p_sz; ++part) {
1194  if (vRTracksMC[iTr].GetPdgCode() == pdgIds[part]
1195  && vRTracksMC[iTr].GetMotherId() == -1) {
1196  tmpmt = sqrt(vRTracksMC[iTr].GetMass() * vRTracksMC[iTr].GetMass()
1197  + vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1198  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy());
1199  nTracksPaAllMC[part]++;
1200  //mtavMC[part] += tmpmt;
1201  //mt2avMC[part] += tmpmt * tmpmt;
1202  mtavMC[part] +=
1203  tempCrit(vRTracksMC[iTr].GetRapidity(),
1204  sqrt(vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1205  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy()),
1206  vRTracksMC[iTr].GetMass());
1207  mt2avMC[part] +=
1208  tempCrit(vRTracksMC[iTr].GetRapidity(),
1209  sqrt(vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1210  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy()),
1211  vRTracksMC[iTr].GetMass())
1212  * tempCrit(
1213  vRTracksMC[iTr].GetRapidity(),
1214  sqrt(vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1215  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy()),
1216  vRTracksMC[iTr].GetMass());
1217  mtmomerrMC[part] += 0.;
1218 
1219  hfyMC[part]->Fill(vRTracksMC[iTr].GetRapidity());
1220  hfdndydptMC[part]->Fill(
1221  vRTracksMC[iTr].GetRapidity(),
1222  sqrt(vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1223  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy()));
1224 
1225  globalnTracksMC[part]++;
1226  //globalmtavMC[part] += tmpmt;
1227  //globalmt2avMC[part] += tmpmt * tmpmt;
1228  globalmtavMC[part] +=
1229  tempCrit(vRTracksMC[iTr].GetRapidity(),
1230  sqrt(vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1231  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy()),
1232  vRTracksMC[iTr].GetMass());
1233  globalmt2avMC[part] +=
1234  tempCrit(vRTracksMC[iTr].GetRapidity(),
1235  sqrt(vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1236  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy()),
1237  vRTracksMC[iTr].GetMass())
1238  * tempCrit(
1239  vRTracksMC[iTr].GetRapidity(),
1240  sqrt(vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1241  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy()),
1242  vRTracksMC[iTr].GetMass());
1243  globalfavMC[part] +=
1244  chi2func(vRTracksMC[iTr].GetRapidity(),
1245  sqrt(vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1246  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy()),
1247  vRTracksMC[iTr].GetMass());
1248  globalf2avMC[part] +=
1249  chi2func(vRTracksMC[iTr].GetRapidity(),
1250  sqrt(vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1251  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy()),
1252  vRTracksMC[iTr].GetMass())
1253  * chi2func(
1254  vRTracksMC[iTr].GetRapidity(),
1255  sqrt(vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1256  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy()),
1257  vRTracksMC[iTr].GetMass());
1258  //globalmt4avMC[part] += tmpmt * tmpmt * tmpmt * tmpmt;
1259  globalmtmomerrMC[part] += 0.;
1260  if (fRecoLevel == 0 && TrRecons[iTr]) {
1261  //if (1) {
1262  nTracksReco[part]++;
1263  //mtavReco[part] += tmpmt;
1264  //mt2avReco[part] += tmpmt * tmpmt;
1265  mtavReco[part] += tempCrit(
1266  vRTracksMC[iTr].GetRapidity(),
1267  sqrt(vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1268  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy()),
1269  vRTracksMC[iTr].GetMass());
1270  mt2avReco[part] +=
1271  tempCrit(
1272  vRTracksMC[iTr].GetRapidity(),
1273  sqrt(vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1274  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy()),
1275  vRTracksMC[iTr].GetMass())
1276  * tempCrit(
1277  vRTracksMC[iTr].GetRapidity(),
1278  sqrt(vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1279  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy()),
1280  vRTracksMC[iTr].GetMass());
1281  mtmomerrReco[part] += 0.;
1282  globalnTracksReco[part]++;
1283  //globalmtavReco[part] += tmpmt;
1284  //globalmt2avReco[part] += tmpmt * tmpmt;
1285  //globalmt4avReco[part] += tmpmt * tmpmt * tmpmt * tmpmt;
1286  globalmtavReco[part] += tempCrit(
1287  vRTracksMC[iTr].GetRapidity(),
1288  sqrt(vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1289  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy()),
1290  vRTracksMC[iTr].GetMass());
1291  globalmt2avReco[part] +=
1292  tempCrit(
1293  vRTracksMC[iTr].GetRapidity(),
1294  sqrt(vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1295  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy()),
1296  vRTracksMC[iTr].GetMass())
1297  * tempCrit(
1298  vRTracksMC[iTr].GetRapidity(),
1299  sqrt(vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1300  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy()),
1301  vRTracksMC[iTr].GetMass());
1302  globalfavReco[part] += chi2func(
1303  vRTracksMC[iTr].GetRapidity(),
1304  sqrt(vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1305  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy()),
1306  vRTracksMC[iTr].GetMass());
1307  globalf2avReco[part] +=
1308  chi2func(
1309  vRTracksMC[iTr].GetRapidity(),
1310  sqrt(vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1311  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy()),
1312  vRTracksMC[iTr].GetMass())
1313  * chi2func(
1314  vRTracksMC[iTr].GetRapidity(),
1315  sqrt(vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1316  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy()),
1317  vRTracksMC[iTr].GetMass());
1318 
1319  globalmtmomerrReco[part] += 0.;
1320 
1321  hfyReco[part]->Fill(vRTracksMC[iTr].GetRapidity());
1322  hfdndydptReco[part]->Fill(
1323  vRTracksMC[iTr].GetRapidity(),
1324  sqrt(vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1325  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy()));
1326 
1327  hfyRecoCor[part]->Fill(vRTracksMC[iTr].GetRapidity());
1328  hfdndydptRecoCor[part]->Fill(
1329  vRTracksMC[iTr].GetRapidity(),
1330  sqrt(vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1331  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy()));
1332  }
1333  }
1334  }
1335  }
1336  }
1337 
1338  //std::map<int, int> trs;
1339  //int clcnt = 0, totsize = 0;
1340 
1341  // Level 1
1342  if (fRecoLevel == 1) {
1343  for (int iTrs = 0; iTrs < nTracks; iTrs++) {
1344  if (vTrackPDG[iTrs] != -1) {
1345  // iTr = ((CbmTrackMatch*)flistStsTracksMatch->At(iTrs))->GetMCTrackId();
1346  iTr = (static_cast<CbmTrackMatch*>(flistStsTracksMatch->At(iTrs)))
1347  ->GetMCTrackId();
1348 
1349  for (int part = 0; part < p_sz; ++part) {
1350  if (vTrackPDG[iTrs] == pdgIds[part]
1351  && vRTracksMC[iTr].GetMotherId() == -1) {
1352  //if (!TrRecons[iTr]) break;
1353  //if (trs.count(iTr)) { trs[iTr]++; /*std::cout << trs[iTr] << " "; */ clcnt++; }
1354  //else trs[iTr] = 1;
1355  //totsize++;
1356  tmpmt = sqrt(vRTracksMC[iTr].GetMass() * vRTracksMC[iTr].GetMass()
1357  + vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1358  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy());
1359  nTracksReco[part]++;
1360  //mtavReco[part] += tmpmt;
1361  //mt2avReco[part] += tmpmt * tmpmt;
1362  mtavReco[part] += tempCrit(
1363  vRTracksMC[iTr].GetRapidity(),
1364  sqrt(vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1365  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy()),
1366  vRTracksMC[iTr].GetMass());
1367  mt2avReco[part] +=
1368  tempCrit(
1369  vRTracksMC[iTr].GetRapidity(),
1370  sqrt(vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1371  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy()),
1372  vRTracksMC[iTr].GetMass())
1373  * tempCrit(
1374  vRTracksMC[iTr].GetRapidity(),
1375  sqrt(vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1376  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy()),
1377  vRTracksMC[iTr].GetMass());
1378  mtmomerrReco[part] += 0.;
1379  globalnTracksReco[part]++;
1380  //globalmtavReco[part] += tmpmt;
1381  //globalmt2avReco[part] += tmpmt * tmpmt;
1382  globalmtavReco[part] += tempCrit(
1383  vRTracksMC[iTr].GetRapidity(),
1384  sqrt(vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1385  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy()),
1386  vRTracksMC[iTr].GetMass());
1387  globalmt2avReco[part] +=
1388  tempCrit(
1389  vRTracksMC[iTr].GetRapidity(),
1390  sqrt(vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1391  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy()),
1392  vRTracksMC[iTr].GetMass())
1393  * tempCrit(
1394  vRTracksMC[iTr].GetRapidity(),
1395  sqrt(vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1396  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy()),
1397  vRTracksMC[iTr].GetMass());
1398  //globalmt4avReco[part] += tmpmt * tmpmt * tmpmt * tmpmt;
1399  globalfavReco[part] += chi2func(
1400  vRTracksMC[iTr].GetRapidity(),
1401  sqrt(vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1402  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy()),
1403  vRTracksMC[iTr].GetMass());
1404  globalf2avReco[part] +=
1405  chi2func(
1406  vRTracksMC[iTr].GetRapidity(),
1407  sqrt(vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1408  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy()),
1409  vRTracksMC[iTr].GetMass())
1410  * chi2func(
1411  vRTracksMC[iTr].GetRapidity(),
1412  sqrt(vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1413  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy()),
1414  vRTracksMC[iTr].GetMass());
1415  globalmtmomerrReco[part] += 0.;
1416 
1417  hfyReco[part]->Fill(vRTracksMC[iTr].GetRapidity());
1418  hfdndydptReco[part]->Fill(
1419  vRTracksMC[iTr].GetRapidity(),
1420  sqrt(vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1421  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy()));
1422 
1423  hfyRecoCor[part]->Fill(vRTracksMC[iTr].GetRapidity());
1424  hfdndydptRecoCor[part]->Fill(
1425  vRTracksMC[iTr].GetRapidity(),
1426  sqrt(vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx()
1427  + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy()));
1428  }
1429  }
1430  }
1431  }
1432  //std::cout << (double)(clcnt) / totsize << endl;
1433  }
1434 
1435  double tmprapid;
1436 
1437  // Level 2
1438  if (fRecoLevel == 2) {
1439  for (int iTrs = 0; iTrs < nTracks; iTrs++) {
1440  if (vTrackPDG[iTrs] != -1) {
1441  // iTr = ((CbmTrackMatch*)flistStsTracksMatch->At(iTrs))->GetMCTrackId();
1442  iTr = (static_cast<CbmTrackMatch*>(flistStsTracksMatch->At(iTrs)))
1443  ->GetMCTrackId();
1444  TVector3 tmpv;
1445  vRTracks[iTrs].GetParamFirst()->Momentum(tmpv);
1446  for (int part = 0; part < p_sz; ++part) {
1447  tmpmt =
1448  sqrt(TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()
1449  * TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()
1450  + tmpv.Pt() * tmpv.Pt());
1451  tmprapid =
1452  0.5
1453  * log(
1454  (sqrt(
1455  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()
1456  * TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()
1457  + tmpv.Mag() * tmpv.Mag())
1458  + tmpv.Pz())
1459  / (sqrt(
1460  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()
1461  * TDatabasePDG::Instance()
1462  ->GetParticle(pdgIds[part])
1463  ->Mass()
1464  + tmpv.Mag() * tmpv.Mag())
1465  - tmpv.Pz()));
1466  if (!(tmpmt > 0. && tmpmt < 10.)) break;
1467  if (vTrackPDG[iTrs] == pdgIds[part]
1468  && vRTracksMC[iTr].GetMotherId() == -1) {
1469  //if (!TrRecons[iTr]) break;
1470  Double_t covMatrix[15];
1471  vRTracks[iTrs].GetParamFirst()->CovMatrix(covMatrix);
1472  double tmperr = getMtErrorSquare(
1473  vRTracks[iTrs].GetParamFirst()->GetQp(),
1474  vRTracks[iTrs].GetParamFirst()->GetTx(),
1475  vRTracks[iTrs].GetParamFirst()->GetTy(),
1476  covMatrix,
1477  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass());
1478  //tmperr = sqrt(tmperr);
1479  if (!(tmperr > 0. && tmperr < 0.1)) break;
1480  //if (tmperr<0. || sqrt(tmperr)>0.001) break;
1481  //tmpmt = sqrt(vRTracksMC[iTr].GetMass()*vRTracksMC[iTr].GetMass() + vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx() + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy());
1482  nTracksReco[part]++;
1483  //mtavReco[part] += tmpmt;
1484  //mt2avReco[part] += tmpmt * tmpmt;
1485  mtavReco[part] += tempCrit(
1486  tmprapid,
1487  sqrt(tmpmt * tmpmt
1488  - TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()
1489  * TDatabasePDG::Instance()
1490  ->GetParticle(pdgIds[part])
1491  ->Mass()),
1492  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass());
1493  mt2avReco[part] +=
1494  tempCrit(
1495  tmprapid,
1496  sqrt(
1497  tmpmt * tmpmt
1498  - TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()
1499  * TDatabasePDG::Instance()
1500  ->GetParticle(pdgIds[part])
1501  ->Mass()),
1502  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass())
1503  * tempCrit(
1504  tmprapid,
1505  sqrt(
1506  tmpmt * tmpmt
1507  - TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()
1508  * TDatabasePDG::Instance()
1509  ->GetParticle(pdgIds[part])
1510  ->Mass()),
1511  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass());
1512 
1513 
1514  mtmomerrReco[part] += getMtErrorSquare(
1515  vRTracks[iTrs].GetParamFirst()->GetQp(),
1516  vRTracks[iTrs].GetParamFirst()->GetTx(),
1517  vRTracks[iTrs].GetParamFirst()->GetTy(),
1518  covMatrix,
1519  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass());
1520  globalnTracksReco[part]++;
1521  //globalmtavReco[part] += tmpmt;
1522  //globalmt2avReco[part] += tmpmt * tmpmt;
1523  globalmtavReco[part] += tempCrit(
1524  tmprapid,
1525  sqrt(tmpmt * tmpmt
1526  - TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()
1527  * TDatabasePDG::Instance()
1528  ->GetParticle(pdgIds[part])
1529  ->Mass()),
1530  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass());
1531  globalmt2avReco[part] +=
1532  tempCrit(
1533  tmprapid,
1534  sqrt(
1535  tmpmt * tmpmt
1536  - TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()
1537  * TDatabasePDG::Instance()
1538  ->GetParticle(pdgIds[part])
1539  ->Mass()),
1540  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass())
1541  * tempCrit(
1542  tmprapid,
1543  sqrt(
1544  tmpmt * tmpmt
1545  - TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()
1546  * TDatabasePDG::Instance()
1547  ->GetParticle(pdgIds[part])
1548  ->Mass()),
1549  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass());
1550  //globalmt4avReco[part] += tmpmt * tmpmt * tmpmt * tmpmt;
1551  globalfavReco[part] += chi2func(
1552  tmprapid,
1553  sqrt(tmpmt * tmpmt
1554  - TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()
1555  * TDatabasePDG::Instance()
1556  ->GetParticle(pdgIds[part])
1557  ->Mass()),
1558  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass());
1559  globalf2avReco[part] +=
1560  chi2func(
1561  tmprapid,
1562  sqrt(
1563  tmpmt * tmpmt
1564  - TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()
1565  * TDatabasePDG::Instance()
1566  ->GetParticle(pdgIds[part])
1567  ->Mass()),
1568  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass())
1569  * chi2func(
1570  tmprapid,
1571  sqrt(
1572  tmpmt * tmpmt
1573  - TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()
1574  * TDatabasePDG::Instance()
1575  ->GetParticle(pdgIds[part])
1576  ->Mass()),
1577  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass());
1578  globalmtmomerrReco[part] += getMtErrorSquare(
1579  vRTracks[iTrs].GetParamFirst()->GetQp(),
1580  vRTracks[iTrs].GetParamFirst()->GetTx(),
1581  vRTracks[iTrs].GetParamFirst()->GetTy(),
1582  covMatrix,
1583  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass());
1584 
1585  hfyReco[part]->Fill(tmprapid);
1586  hfdndydptReco[part]->Fill(
1587  tmprapid,
1588  sqrt(tmpmt * tmpmt
1589  - TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()
1590  * TDatabasePDG::Instance()
1591  ->GetParticle(pdgIds[part])
1592  ->Mass()));
1593 
1594  hfyRecoCor[part]->Fill(tmprapid);
1595  hfdndydptRecoCor[part]->Fill(
1596  tmprapid,
1597  sqrt(tmpmt * tmpmt
1598  - TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()
1599  * TDatabasePDG::Instance()
1600  ->GetParticle(pdgIds[part])
1601  ->Mass()));
1602  }
1603  }
1604  }
1605  }
1606  }
1607 
1608  // Level 3
1609  if (fRecoLevel == 3) {
1610  for (int iTrs = 0; iTrs < nTracks; iTrs++) {
1611  if (vTrackPDG[iTrs] != -1) {
1612  // iTr = ((CbmTrackMatch*)flistStsTracksMatch->At(iTrs))->GetMCTrackId();
1613  iTr = (static_cast<CbmTrackMatch*>(flistStsTracksMatch->At(iTrs)))
1614  ->GetMCTrackId();
1615  TVector3 tmpv;
1616  vRTracks[iTrs].GetParamFirst()->Momentum(tmpv);
1617  for (int part = 0; part < p_sz; ++part) {
1618  tmpmt =
1619  sqrt(TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()
1620  * TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()
1621  + tmpv.Pt() * tmpv.Pt());
1622  tmprapid =
1623  0.5
1624  * log(
1625  (sqrt(
1626  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()
1627  * TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()
1628  + tmpv.Mag() * tmpv.Mag())
1629  + tmpv.Pz())
1630  / (sqrt(
1631  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()
1632  * TDatabasePDG::Instance()
1633  ->GetParticle(pdgIds[part])
1634  ->Mass()
1635  + tmpv.Mag() * tmpv.Mag())
1636  - tmpv.Pz()));
1637  if (!(tmpmt > 0. && tmpmt < 10.)) break;
1638  //std::cout << tmprapid << " ";
1639  //if (tmprapid<0.) break;
1640  //if (tmprapid>6.) break;
1641  //if (tmprapid<0. || tmprapid>6. || tmpv.Pt()<0. || tmpv.Pt()>2.5) continue;
1642  if (vTrackPDG[iTrs] == pdgIds[part]
1643  && fChiToPrimVtx[iTrs]
1644  < 3.) { //vRTracksMC[iTr].GetMotherId()==-1) {
1645  //tmpmt = sqrt(vRTracksMC[iTr].GetMass()*vRTracksMC[iTr].GetMass() + vRTracksMC[iTr].GetPx() * vRTracksMC[iTr].GetPx() + vRTracksMC[iTr].GetPy() * vRTracksMC[iTr].GetPy());
1646  //if (!TrRecons[iTr]) break;
1647  Double_t covMatrix[15];
1648  vRTracks[iTrs].GetParamFirst()->CovMatrix(covMatrix);
1649  double tmperr = getMtErrorSquare(
1650  vRTracks[iTrs].GetParamFirst()->GetQp(),
1651  vRTracks[iTrs].GetParamFirst()->GetTx(),
1652  vRTracks[iTrs].GetParamFirst()->GetTy(),
1653  covMatrix,
1654  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass());
1655  //tmperr = sqrt(tmperr);
1656  if (!(tmperr > 0. && tmperr < 0.1)) break;
1657 
1658  nTracksReco[part]++;
1659  //mtavReco[part] += tmpmt;
1660  //mt2avReco[part] += tmpmt * tmpmt;
1661  mtavReco[part] += tempCrit(
1662  tmprapid,
1663  sqrt(tmpmt * tmpmt
1664  - TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()
1665  * TDatabasePDG::Instance()
1666  ->GetParticle(pdgIds[part])
1667  ->Mass()),
1668  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass());
1669  mt2avReco[part] +=
1670  tempCrit(
1671  tmprapid,
1672  sqrt(
1673  tmpmt * tmpmt
1674  - TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()
1675  * TDatabasePDG::Instance()
1676  ->GetParticle(pdgIds[part])
1677  ->Mass()),
1678  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass())
1679  * tempCrit(
1680  tmprapid,
1681  sqrt(
1682  tmpmt * tmpmt
1683  - TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()
1684  * TDatabasePDG::Instance()
1685  ->GetParticle(pdgIds[part])
1686  ->Mass()),
1687  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass());
1688 
1689 
1690  mtmomerrReco[part] += getMtErrorSquare(
1691  vRTracks[iTrs].GetParamFirst()->GetQp(),
1692  vRTracks[iTrs].GetParamFirst()->GetTx(),
1693  vRTracks[iTrs].GetParamFirst()->GetTy(),
1694  covMatrix,
1695  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass());
1696 
1697  globalnTracksReco[part]++;
1698  //globalmtavReco[part] += tmpmt;
1699  //globalmt2avReco[part] += tmpmt * tmpmt;
1700  globalmtavReco[part] += tempCrit(
1701  tmprapid,
1702  sqrt(tmpmt * tmpmt
1703  - TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()
1704  * TDatabasePDG::Instance()
1705  ->GetParticle(pdgIds[part])
1706  ->Mass()),
1707  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass());
1708  globalmt2avReco[part] +=
1709  tempCrit(
1710  tmprapid,
1711  sqrt(
1712  tmpmt * tmpmt
1713  - TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()
1714  * TDatabasePDG::Instance()
1715  ->GetParticle(pdgIds[part])
1716  ->Mass()),
1717  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass())
1718  * tempCrit(
1719  tmprapid,
1720  sqrt(
1721  tmpmt * tmpmt
1722  - TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()
1723  * TDatabasePDG::Instance()
1724  ->GetParticle(pdgIds[part])
1725  ->Mass()),
1726  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass());
1727  globalfavReco[part] += chi2func(
1728  tmprapid,
1729  sqrt(tmpmt * tmpmt
1730  - TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()
1731  * TDatabasePDG::Instance()
1732  ->GetParticle(pdgIds[part])
1733  ->Mass()),
1734  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass());
1735  globalf2avReco[part] +=
1736  chi2func(
1737  tmprapid,
1738  sqrt(
1739  tmpmt * tmpmt
1740  - TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()
1741  * TDatabasePDG::Instance()
1742  ->GetParticle(pdgIds[part])
1743  ->Mass()),
1744  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass())
1745  * chi2func(
1746  tmprapid,
1747  sqrt(
1748  tmpmt * tmpmt
1749  - TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()
1750  * TDatabasePDG::Instance()
1751  ->GetParticle(pdgIds[part])
1752  ->Mass()),
1753  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass());
1754  globalmtmomerrReco[part] += getMtErrorSquare(
1755  vRTracks[iTrs].GetParamFirst()->GetQp(),
1756  vRTracks[iTrs].GetParamFirst()->GetTx(),
1757  vRTracks[iTrs].GetParamFirst()->GetTy(),
1758  covMatrix,
1759  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass());
1760 
1761  hfyReco[part]->Fill(tmprapid);
1762  hfdndydptReco[part]->Fill(
1763  tmprapid,
1764  sqrt(tmpmt * tmpmt
1765  - TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()
1766  * TDatabasePDG::Instance()
1767  ->GetParticle(pdgIds[part])
1768  ->Mass()));
1769 
1770  hfyRecoCor[part]->Fill(tmprapid);
1771  hfdndydptRecoCor[part]->Fill(
1772  tmprapid,
1773  sqrt(tmpmt * tmpmt
1774  - TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()
1775  * TDatabasePDG::Instance()
1776  ->GetParticle(pdgIds[part])
1777  ->Mass()));
1778  }
1779  }
1780  }
1781  }
1782  }
1783 
1784 
1785  for (int part = 0; part < p_sz; ++part) {
1786  mtavMC[part] /= nTracksPaAllMC[part];
1787  mt2avMC[part] /= nTracksPaAllMC[part];
1788  }
1789 
1790  for (int part = 0; part < p_sz; ++part) {
1791  //if (part==0)
1792  {
1793  double tmpT = getTemperatureAll(mtavMC[part], part, 20);
1794  hTempFullMC[part]->Fill(tmpT);
1795  hTempErrStatFullMC[part]->Fill(
1796  getTemperatureDerivAll(mtavMC[part], part, 20)
1797  * sqrt((mt2avMC[part] - mtavMC[part] * mtavMC[part])
1798  / (nTracksPaAllMC[part] - 1)));
1799  hTempErrMomFullMC[part]->Fill(
1800  getTemperatureDerivAll(mtavMC[part], part, 20) * sqrt(mtmomerrMC[part])
1801  / nTracksPaAllMC[part]);
1802 
1803  hRadFullMC[part]->Fill(
1804  getRadius(tmpT,
1805  nTracksPaAllMC[part],
1806  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()));
1807  double tmpTerr = getTemperatureDerivAll(mtavMC[part], part, 20)
1808  * sqrt((mt2avMC[part] - mtavMC[part] * mtavMC[part])
1809  / (nTracksPaAllMC[part] - 1));
1810  double tmp1 = 0., tmp2 = 0.;
1811  double tmpNerr = sqrt(nTracksPaAllMC[part]);
1812  tmp1 = getRadiusDerivT(
1813  tmpT,
1814  nTracksPaAllMC[part],
1815  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass());
1816  tmp2 = getRadiusDerivN(
1817  tmpT,
1818  nTracksPaAllMC[part],
1819  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass());
1820  double tmpRerr =
1821  sqrt(tmp1 * tmp1 * tmpTerr * tmpTerr + tmp2 * tmp2 * tmpNerr * tmpNerr);
1822  hRadErrStatFullMC[part]->Fill(tmpRerr);
1823 
1824  tmpTerr = getTemperatureDerivAll(mtavMC[part], part, 20)
1825  * sqrt(mtmomerrMC[part]) / nTracksPaAllMC[part];
1826  tmpRerr = sqrt(tmp1 * tmp1 * tmpTerr * tmpTerr
1827  + tmp2 * tmp2 * tmpNerr * tmpNerr * 0.);
1828  hRadErrMomFullMC[part]->Fill(tmpRerr);
1829  }
1830  }
1831 
1832  //for(int rl=0;rl<recoLevels;++rl) {
1833  for (int part = 0; part < p_sz; ++part) {
1834  mtavReco[part] /= nTracksReco[part];
1835  mt2avReco[part] /= nTracksReco[part];
1836  }
1837 
1838  for (int part = 0; part < p_sz; ++part) {
1839  //if (part==0)
1840  {
1841  mtTacc[part] = mtTsts[part];
1842  Npart[part] = Npartsts[part];
1843  if (fusePID == 2 && fRecoLevel >= 1) {
1844  mtTacc[part] = mtTststof[part];
1845  Npart[part] = Npartststof[part];
1846  }
1847  double tmpT = getTemperatureAll(mtavMC[part], part, 20);
1848  //hTempFull[part]->Fill(tmpT);
1849  //hRadFull[part]->Fill(getRadius(tmpT, nTracksPaAll[part], TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()));
1850  tmpT = getTemperatureAll(mtavReco[part], part, 20);
1851  hTempFullReco[part]->Fill(tmpT);
1852  hRadFullReco[part]->Fill(
1853  getRadius(tmpT,
1854  nTracksReco[part],
1855  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()));
1856 
1857  double tmpTerr =
1858  getTemperatureDerivAll(mtavReco[part], part, 20)
1859  * sqrt((mt2avReco[part] - mtavReco[part] * mtavReco[part])
1860  / (nTracksReco[part] - 1));
1861  double tmp1 = 0., tmp2 = 0.;
1862  double tmpNerr = sqrt(nTracksReco[part]);
1863  tmp1 = getRadiusDerivT(
1864  tmpT,
1865  nTracksReco[part],
1866  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass());
1867  tmp2 = getRadiusDerivN(
1868  tmpT,
1869  nTracksReco[part],
1870  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass());
1871  double tmpRerr =
1872  sqrt(tmp1 * tmp1 * tmpTerr * tmpTerr + tmp2 * tmp2 * tmpNerr * tmpNerr);
1873 
1874  hTempErrStatFullReco[part]->Fill(tmpTerr);
1875  hRadErrStatFullReco[part]->Fill(tmpRerr);
1876 
1877  tmpTerr = getTemperatureDerivAll(mtavReco[part], part, 20)
1878  * sqrt(mtmomerrReco[part]) / nTracksReco[part];
1879  tmpRerr =
1880  sqrt(tmp1 * tmp1 * tmpTerr * tmpTerr + tmp2 * tmp2 * tmpNerr * tmpNerr);
1881 
1882  hTempErrMomFullReco[part]->Fill(tmpTerr);
1883  hRadErrMomFullReco[part]->Fill(tmpRerr);
1884 
1885  tmpT = getTemperatureAllCor(mtavReco[part], part, 20, mtTacc[part]);
1886  hTempFullRecoCor[part]->Fill(tmpT);
1887  hRadFullRecoCor[part]->Fill(getRadiusCor(
1888  tmpT,
1889  nTracksReco[part],
1890  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass(),
1891  Npart[part]));
1892 
1893  tmpTerr =
1894  getTemperatureDerivAllCor(mtavReco[part], part, 20, mtTacc[part])
1895  * sqrt((mt2avReco[part] - mtavReco[part] * mtavReco[part])
1896  / (nTracksReco[part] - 1));
1897  tmp1 = 0., tmp2 = 0.;
1898  tmpNerr = sqrt(nTracksReco[part]);
1899  tmp1 = getRadiusDerivTCor(
1900  tmpT,
1901  nTracksReco[part],
1902  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass(),
1903  Npart[part]);
1904  tmp2 = getRadiusDerivNCor(
1905  tmpT,
1906  nTracksReco[part],
1907  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass(),
1908  Npart[part]);
1909  tmpRerr =
1910  sqrt(tmp1 * tmp1 * tmpTerr * tmpTerr + tmp2 * tmp2 * tmpNerr * tmpNerr);
1911 
1912  hTempErrStatFullRecoCor[part]->Fill(tmpTerr);
1913  hRadErrStatFullRecoCor[part]->Fill(tmpRerr);
1914 
1915  tmpTerr =
1916  getTemperatureDerivAllCor(mtavReco[part], part, 20, mtTacc[part])
1917  * sqrt(mtmomerrReco[part]) / nTracksReco[part];
1918  tmpRerr = sqrt(tmp1 * tmp1 * tmpTerr * tmpTerr
1919  + tmp2 * tmp2 * tmpNerr * tmpNerr * 0.);
1920 
1921  hTempErrMomFullRecoCor[part]->Fill(tmpTerr);
1922  hRadErrMomFullRecoCor[part]->Fill(tmpRerr);
1923  }
1924  }
1925  //}
1926 
1927  delete[] MCTrackSortedArray;
1928  delete[] TrRecons;
1929 }
1930 
1931 // Function to test method of moments
1932 double CbmThermalModelNoFlow::chi2func(double /*y*/, double pt, double m) {
1933  return sqrt(pt * pt + m * m); //y*y;//pt*pt + m*m;
1934 }
1935 
1936 // Function to determine temperature
1937 double CbmThermalModelNoFlow::tempCrit(double /*y*/, double pt, double m) {
1938  return sqrt(pt * pt + m * m); //*cosh(y);
1939 }
1940 
1941 void CbmThermalModelNoFlow::ComputeThermalDependence(Int_t part) {
1942  TString work = getenv("VMCWORKDIR");
1943  TString dir =
1944  work + "/KF/KFModelParameters/"; //ThermalModel.mT-t.sts_v13.txt";
1945  ReadAcceptanceFunction(AcceptanceSTS, dir + "pty_acc_sts_v13.txt");
1946  ReadAcceptanceFunction(AcceptanceSTSTOF, dir + "pty_acc_ststof_v13.txt");
1947 
1948  std::vector<double> mtall, Tvecall;
1949 
1950  std::cout << "Computing dependencies of ThermalModel...";
1951 
1952  fflush(stdout);
1953 
1954  std::vector<double> Tvec, mtacc, NPart;
1955  mtall.resize(0);
1956  double dT = 0.002;
1957  double Tmax = 0.5;
1958  for (int i = 0;; ++i) {
1959  double tmpT = dT * (0.5 + i);
1960  if (tmpT > Tmax) break;
1961  //cout << tmpT << "\n";
1962  Tvec.push_back(tmpT);
1963  mtall.push_back(ThermalMt(
1964  tmpT, TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()));
1965  mtacc.push_back(
1966  ThermalMtAcc(tmpT,
1967  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass(),
1968  AcceptanceSTS));
1969  NPart.push_back(
1970  NFracAcc(tmpT,
1971  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass(),
1972  AcceptanceSTS));
1973  //cout << tmpT << " " << mtall[mtall.size()-1] << "\n";
1974  }
1975 
1976  //mtTall = new TSpline3("mtTall", &Tvecall[0], &mtall[0], Tvecall.size());
1977  mtTall[part] = new TSpline3("mtTall", &Tvec[0], &mtall[0], Tvec.size());
1978  mtTsts[part] = new TSpline3("mtTacc", &Tvec[0], &mtacc[0], Tvec.size());
1979  Npartsts[part] = new TSpline3("Npart", &Tvec[0], &NPart[0], Tvec.size());
1980 
1981  Tvec.resize(0);
1982  //mtall.resize(0);
1983  mtacc.resize(0);
1984  NPart.resize(0);
1985 
1986  dT = 0.002;
1987  Tmax = 0.5;
1988  for (int i = 0;; ++i) {
1989  double tmpT = dT * (0.5 + i);
1990  if (tmpT > Tmax) break;
1991  //cout << tmpT << "\n";
1992  Tvec.push_back(tmpT);
1993  //mtall.push_back(ThermalMt(tmpT, TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()));
1994  mtacc.push_back(
1995  ThermalMtAcc(tmpT,
1996  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass(),
1997  AcceptanceSTSTOF));
1998  NPart.push_back(
1999  NFracAcc(tmpT,
2000  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass(),
2001  AcceptanceSTSTOF));
2002  }
2003 
2004  mtTststof[part] = new TSpline3("mtTacc", &Tvec[0], &mtacc[0], Tvec.size());
2005  Npartststof[part] = new TSpline3("Npart", &Tvec[0], &NPart[0], Tvec.size());
2006 
2007  if (fusePID == 1) {
2008  mtTacc[part] = mtTsts[part];
2009  Npart[part] = Npartsts[part];
2010  } else {
2011  mtTacc[part] = mtTststof[part];
2012  Npart[part] = Npartststof[part];
2013  }
2014 
2015  std::cout << " Done!\n";
2016 }
2017 
2018 Double_t CbmThermalModelNoFlow::ThermalMt(double T, double m) {
2019  double ret1 = 0., ret2 = 0.;
2020  vector<double> xlag, wlag, xleg, wleg;
2021  //GetCoefs2DLaguerre32Legendre32(0., 6., xlag, wlag, xleg, wleg);
2022  GetCoefs2DLegendre32Legendre32(0., 3., -3., 6., xlag, wlag, xleg, wleg);
2023  for (Int_t i = 0; i < 32; i++) {
2024  for (Int_t j = 0; j < 32; j++) {
2025  //if (xlag[i]>3.) break;
2026  //cout << xlag[i] << " " << xleg[j] << " " << accfunc.getAcceptance(xleg[j], xlag[i]) << "\t";
2027  double tmpf = 0.;
2028  tmpf = xlag[i] * sqrt(xlag[i] * xlag[i] + m * m)
2029  * TMath::CosH(xleg[j] - ycm)
2030  * exp(-sqrt(xlag[i] * xlag[i] + m * m) * TMath::CosH(xleg[j] - ycm)
2031  / T);
2032  ret1 += wlag[i] * wleg[j] * tmpf;
2033  ret2 += wlag[i] * wleg[j] * tmpf * tempCrit(xleg[j], xlag[i], m);
2034  }
2035  }
2036  return ret2 / ret1;
2037 }
2038 
2039 Double_t CbmThermalModelNoFlow::ThermalMt2(double T, double m) {
2040  double ret1 = 0., ret2 = 0.;
2041  double ymin = -3., ymax = 6.;
2042  int itery = 1500;
2043  double dy = (ymax - ymin) / itery;
2044  double ptmin = 0., ptmax = 3.;
2045  int iterpt = 1200;
2046  double dpt = (ptmax - ptmin) / iterpt;
2047  double ty = 0., tpt = 0., tmp = 0.;
2048  for (int iy = 0; iy < itery; ++iy) {
2049  ty = ymin + (iy + 0.5) * dy;
2050  for (int ipt = 0; ipt < iterpt; ++ipt) {
2051  tpt = ptmin + (ipt + 0.5) * dpt;
2052  tmp = tpt * sqrt(tpt * tpt + m * m) * cosh(ty - ycm)
2053  * exp(-sqrt(tpt * tpt + m * m) * cosh(ty - ycm) / T);
2054  ret1 += tmp;
2055  //ret2 += tmp * (tpt*tpt + m*m);
2056  ret2 += tmp * chi2func(ty, tpt, m);
2057  }
2058  }
2059  //std::cout << "T = " << T << "\tmT = " << ret2 / ret1 << endl;
2060  return ret2 / ret1;
2061 }
2062 
2063 Double_t CbmThermalModelNoFlow::RecEfficiency(Double_t p) {
2064  if (fTrackNumber == 0)
2065  return 0.99 - 0.98 * exp(-p * p / 2. / 0.135 / 0.135);
2066  else
2067  return 0.98 - 0.88 * exp(-p * p / 2. / 0.2 / 0.2);
2068 }
2069 
2070 Double_t
2071 CbmThermalModelNoFlow::ThermalMtAcc(double T,
2072  double m,
2073  const AcceptanceFunction& accfunc) {
2074  double ret1 = 0., ret2 = 0.;
2075  vector<double> xlag, wlag, xleg, wleg;
2076  //GetCoefs2DLaguerre32Legendre32(0., 6., xlag, wlag, xleg, wleg);
2077  GetCoefs2DLegendre32Legendre32(0., 2.5, 0., 6., xlag, wlag, xleg, wleg);
2078  //cout << ycm << "\n";
2079  for (Int_t i = 0; i < 32; i++) {
2080  for (Int_t j = 0; j < 32; j++) {
2081  //if (xlag[i]>3.) break;
2082  //cout << xlag[i] << " " << xleg[j] << " " << accfunc.getAcceptance(xleg[j], xlag[i]) << "\t";
2083  double tmpf = 0.;
2084  double tp = sqrt(xlag[i] * xlag[i] + m * m) * TMath::CosH(xleg[j]);
2085  tp = sqrt(tp * tp - m * m);
2086  tmpf =
2087  xlag[i] * sqrt(xlag[i] * xlag[i] + m * m) * TMath::CosH(xleg[j] - ycm)
2088  * exp(-sqrt(xlag[i] * xlag[i] + m * m) * TMath::CosH(xleg[j] - ycm) / T)
2089  * accfunc.getAcceptance(xleg[j], xlag[i]) * RecEfficiency(tp);
2090  ret1 += wlag[i] * wleg[j] * tmpf;
2091  ret2 += wlag[i] * wleg[j] * tmpf * tempCrit(xleg[j], xlag[i], m);
2092  }
2093  }
2094  return ret2 / ret1;
2095 }
2096 
2097 Double_t
2098 CbmThermalModelNoFlow::ThermalMt2Acc(double T,
2099  double m,
2100  const AcceptanceFunction& accfunc) {
2101  if (accfunc.ys.size() == 0) return -1.;
2102  double ret1 = 0., ret2 = 0.;
2103  double tmp = 0., tp = 0.;
2104  for (unsigned int i = 0; i < accfunc.ys.size(); ++i) {
2105  tp = sqrt(accfunc.pts[i] * accfunc.pts[i] + m * m) * cosh(accfunc.ys[i]);
2106  tp = sqrt(tp * tp - m * m);
2107  tmp = accfunc.pts[i] * sqrt(accfunc.pts[i] * accfunc.pts[i] + m * m)
2108  * cosh(accfunc.ys[i] - ycm)
2109  * exp(-sqrt(accfunc.pts[i] * accfunc.pts[i] + m * m)
2110  * cosh(accfunc.ys[i] - ycm) / T)
2111  * accfunc.probs[i] * RecEfficiency(tp);
2112  ret1 += tmp;
2113  //ret2 += tmp * (accfunc.pts[i]*accfunc.pts[i] + m*m);
2114  ret2 += tmp * chi2func(accfunc.ys[i], accfunc.pts[i], m);
2115  }
2116  //std::cout << "T = " << T << "\tmT = " << ret2 / ret1 << endl;
2117  return ret2 / ret1;
2118 }
2119 
2120 Double_t CbmThermalModelNoFlow::NFracAcc(double T,
2121  double m,
2122  const AcceptanceFunction& accfunc) {
2123  double ret1 = 0., ret2 = 0.;
2124  vector<double> xlag, wlag, xleg, wleg;
2125  //GetCoefs2DLaguerre32Legendre32(0., 6., xlag, wlag, xleg, wleg);
2126  GetCoefs2DLegendre32Legendre32(0., 2.5, 0., 6., xlag, wlag, xleg, wleg);
2127  for (Int_t i = 0; i < 32; i++) {
2128  for (Int_t j = 0; j < 32; j++) {
2129  //if (xlag[i]>3.) break;
2130  //cout << xlag[i] << " " << xleg[j] << " " << accfunc.getAcceptance(xleg[j], xlag[i]) << "\t";
2131  double tmpf = 0.;
2132  double tp = sqrt(xlag[i] * xlag[i] + m * m) * TMath::CosH(xleg[j]);
2133  tp = sqrt(tp * tp - m * m);
2134  tmpf =
2135  xlag[i] * sqrt(xlag[i] * xlag[i] + m * m) * TMath::CosH(xleg[j] - ycm)
2136  * exp(-sqrt(xlag[i] * xlag[i] + m * m) * TMath::CosH(xleg[j] - ycm) / T)
2137  * accfunc.getAcceptance(xleg[j], xlag[i]) * RecEfficiency(tp);
2138  ret2 += wlag[i] * wleg[j] * tmpf;
2139  //ret2 += wlag[i]*wleg[j]*tmpf * tempCrit(xleg[j], xlag[i], m);
2140  }
2141  }
2142  ret2 *= 1. / 2. / 2. / TMath::Pi() / TMath::Pi();
2143  ret1 = T * m * m * TMath::BesselK(2, m / T) / 2. / TMath::Pi() / TMath::Pi();
2144  return ret2 / ret1;
2145 }
2146 
2147 bool CbmThermalModelNoFlow::checkIfReconstructable(
2148  CbmKFTrErrMCPoints* inTrack) {
2149  //bool ret = 0;
2150  int stsHits = inTrack->GetNStsPoints();
2151  int tofHits = inTrack->GetNTofPoints();
2152  if (fusePID == 2 && tofHits == 0) return 0;
2153  vector<int> hitz(0);
2154  for (int hit = 0; hit < stsHits; ++hit) {
2155  hitz.push_back(
2156  static_cast<int>((inTrack->GetStsPoint(hit)->GetZ() + 5.) / 10.));
2157  }
2158  sort(hitz.begin(), hitz.end());
2159  for (int hit1 = 0; hit1 < stsHits - 3; ++hit1) {
2160  int station1 = hitz[hit1];
2161  int hitsConsecutive = 1;
2162  for (int hit2 = hit1 + 1; hit2 < stsHits && hitsConsecutive < 4; ++hit2) {
2163  int station2 = hitz[hit2];
2164  if (station2 == station1) continue;
2165  if (station2 - station1 > 1) break;
2166  if (station2 - station1 == 1) {
2167  hitsConsecutive++;
2168  station1 = station2;
2169  }
2170  }
2171  if (hitsConsecutive == 4) return 1;
2172  }
2173  return 0;
2174 }
2175 
2176 double CbmThermalModelNoFlow::getRadius(double T, double N, double mo) {
2177  //std::cout << T << "\t" << N << endl;
2178  double V = (2 * TMath::Pi() * TMath::Pi() * N) / T / mo / mo
2179  / TMath::BesselK(2, mo / T);
2180  return TMath::Power(3. * V / 4. / TMath::Pi(), 1. / 3.) / GeVtoifm;
2181 }
2182 
2183 double CbmThermalModelNoFlow::getRadiusDerivT(double T, double N, double mo) {
2184  double dT = 0.005;
2185  return fabs(getRadius(T + dT, N, mo) - getRadius(T, N, mo)) / dT;
2186 }
2187 
2188 double CbmThermalModelNoFlow::getRadiusDerivN(double T, double N, double mo) {
2189  double dN = 0.1;
2190  return fabs(getRadius(T, N + dN, mo) - getRadius(T, N, mo)) / dN;
2191 }
2192 
2193 double CbmThermalModelNoFlow::getRadiusCor(double T,
2194  double N,
2195  double mo,
2196  TSpline3* Np) {
2197  //std::cout << T << "\t" << N << "\t" << N/Npart->Eval(T) << endl ;
2198  double V = (2 * TMath::Pi() * TMath::Pi() * N / Np->Eval(T)) / T / mo / mo
2199  / TMath::BesselK(2, mo / T);
2200  return TMath::Power(3. * V / 4. / TMath::Pi(), 1. / 3.) / GeVtoifm;
2201 }
2202 
2203 double CbmThermalModelNoFlow::getRadiusDerivTCor(double T,
2204  double N,
2205  double mo,
2206  TSpline3* Np) {
2207  double dT = 0.005;
2208  return fabs(getRadiusCor(T + dT, N, mo, Np) - getRadiusCor(T, N, mo, Np))
2209  / dT;
2210 }
2211 
2212 double CbmThermalModelNoFlow::getRadiusDerivNCor(double T,
2213  double N,
2214  double mo,
2215  TSpline3* Np) {
2216  double dN = 0.1;
2217  return fabs(getRadiusCor(T, N + dN, mo, Np) - getRadiusCor(T, N, mo, Np))
2218  / dN;
2219 }
2220 
2221 double CbmThermalModelNoFlow::getTemperature(double mt, double mo, int iters) {
2222  double left = 0., right = 1., center;
2223  double valleft =
2224  tempFunction(left, mt, mo); //, valright = tempFunction(right, mt, mo)
2225  double valcenter;
2226  for (int i = 0; i < iters; ++i) {
2227  center = (left + right) / 2.;
2228  valcenter = tempFunction(center, mt, mo);
2229  if (valleft * valcenter > 0)
2230  left = center;
2231  else
2232  right = center;
2233  }
2234  return (left + right) / 2.;
2235 }
2236 
2237 double
2238 CbmThermalModelNoFlow::getTemperatureAll(double mt, int part, int iters) {
2239  double left = 0., right = 0.45, center;
2240  double valleft =
2241  mtTall[part]->Eval(left) - mt; //, valright = mtTall[part]->Eval(right)-mt,
2242  double valcenter;
2243  for (int i = 0; i < iters; ++i) {
2244  center = (left + right) / 2.;
2245  valcenter = mtTall[part]->Eval(center) - mt;
2246  if (valleft * valcenter > 0)
2247  left = center;
2248  else
2249  right = center;
2250  }
2251  return (left + right) / 2.;
2252 }
2253 
2254 double
2255 CbmThermalModelNoFlow::getTemperatureDerivAll(double mt, int part, int iters) {
2256  double dmt = 0.005;
2257  return fabs((getTemperatureAll(mt + dmt, part, iters)
2258  - getTemperatureAll(mt, part, iters)))
2259  / dmt;
2260 }
2261 
2262 double CbmThermalModelNoFlow::getTemperatureAllCor(double mt,
2263  int /*part*/,
2264  int iters,
2265  TSpline3* mtT) {
2266  double left = 0., right = 0.45, center;
2267  double valleft = mtT->Eval(left) - mt; //, valright = mtT->Eval(right)-mt,
2268  double valcenter;
2269  for (int i = 0; i < iters; ++i) {
2270  center = (left + right) / 2.;
2271  valcenter = mtT->Eval(center) - mt;
2272  if (valleft * valcenter > 0)
2273  left = center;
2274  else
2275  right = center;
2276  }
2277  return (left + right) / 2.;
2278 }
2279 
2280 double CbmThermalModelNoFlow::getTemperatureDerivAllCor(double mt,
2281  int part,
2282  int iters,
2283  TSpline3* mtT) {
2284  double dmt = 0.005;
2285  return fabs((getTemperatureAllCor(mt + dmt, part, iters, mtT)
2286  - getTemperatureAllCor(mt, part, iters, mtT)))
2287  / dmt;
2288 }
2289 
2290 double CbmThermalModelNoFlow::getTemperatureRapidity(double mt,
2291  double mo,
2292  int iters,
2293  double y) {
2294  double left = 0., right = 1., center;
2295  double valleft =
2296  tempFunction(left, mt, mo); //, valright = tempFunction(right, mt, mo),
2297  double valcenter;
2298  for (int i = 0; i < iters; ++i) {
2299  center = (left + right) / 2.;
2300  valcenter = tempFunction(center, mt, mo);
2301  if (valleft * valcenter > 0)
2302  left = center;
2303  else
2304  right = center;
2305  }
2306  return (left + right) * 0.5 * cosh(y);
2307 }
2308 
2309 double CbmThermalModelNoFlow::getInverseSlope(double mt, double mo, int
2310  /*iters*/) {
2311  return 0.5
2312  * (0.5 * mt - mo
2313  + sqrt((0.5 * mt - mo) * (0.5 * mt - mo) + 2 * mo * (mt - mo)));
2314 }
2315 
2316 double CbmThermalModelNoFlow::getDerivTemperature(double mt,
2317  double mo,
2318  double T,
2319  int iters) {
2320  double left = 0., right = 10., center;
2321  double valleft = tempDerivFunction(
2322  T, left, mt, mo); //, valright = tempDerivFunction(T, right, mt, mo),
2323  double valcenter;
2324  for (int i = 0; i < iters; ++i) {
2325  center = (left + right) / 2.;
2326  valcenter = tempDerivFunction(T, center, mt, mo);
2327  if (valleft * valcenter > 0)
2328  left = center;
2329  else
2330  right = center;
2331  }
2332  return (left + right) / 2.;
2333 }
2334 
2335 double CbmThermalModelNoFlow::getMtErrorSquare(double qp,
2336  double Tx,
2337  double Ty,
2338  double covMatrix[],
2339  double mo) {
2340  double tmpTxy = Tx * Tx + Ty * Ty;
2341  double mt = sqrt(mo * mo + tmpTxy / qp / qp / (tmpTxy + 1.));
2342  double ddqp = -2. / qp / qp / qp * tmpTxy / (tmpTxy + 1.);
2343  double ddTx = 2. * Tx / qp / qp / (tmpTxy + 1.) * (tmpTxy + 1.);
2344  double ddTy = 2. * Ty / qp / qp / (tmpTxy + 1.) * (tmpTxy + 1.);
2345  //return (1./4./mt/mt) * ( ddqp*ddqp*covMatrix[14] + ddTx*ddTx*covMatrix[5] + ddTy*ddTy*covMatrix[9]
2346  // + 2. * ( ddqp*ddTx*covMatrix[12] + ddqp*ddTy*covMatrix[13] + ddTx*ddTy*covMatrix[8] ) );
2347  return (1. / 4. / mt / mt)
2348  * (ddqp * ddqp * covMatrix[14] + ddTx * ddTx * covMatrix[9]
2349  + ddTy * ddTy * covMatrix[12]
2350  + 2.
2351  * (ddqp * ddTx * covMatrix[11] + ddqp * ddTy * covMatrix[13]
2352  + ddTx * ddTy * covMatrix[10]));
2353 }
2354 
2355 double CbmThermalModelNoFlow::tempFunction(double T, double mt, double mo) {
2356  return mt * (mo * mo + 2. * mo * T + 2 * T * T)
2357  - (mo * mo * mo + 3 * mo * mo * T + 6 * mo * T * T + 6 * T * T * T);
2358 }
2359 
2360 double CbmThermalModelNoFlow::tempDerivFunction(double T,
2361  double dT,
2362  double mt,
2363  double mo) {
2364  //return mt*(mo*mo+2.*mo*T+2*T*T) - (mo*mo*mo+3*mo*mo*T+6*mo*T*T+6*T*T*T);
2365  return mt * (2 * mo * dT + 4 * T * dT) + mo * mo + 2 * mo * T + 2 * T * T
2366  - (3. * mo * mo * dT + 12 * mo * T * dT + 18. * T * T * dT);
2367 }
2368 
2369 void CbmThermalModelNoFlow::Finish() {
2370  ofstream fout("ThermalModel.txt");
2371  for (int part = 0; part < p_sz; ++part) {
2372  globalmtavMC[part] /= globalnTracksMC[part];
2373  globalmt2avMC[part] /= globalnTracksMC[part];
2374  globalfavMC[part] /= globalnTracksMC[part];
2375  globalf2avMC[part] /= globalnTracksMC[part];
2376  globalmtavReco[part] /= globalnTracksReco[part];
2377  globalmt2avReco[part] /= globalnTracksReco[part];
2378  globalfavReco[part] /= globalnTracksReco[part];
2379  globalf2avReco[part] /= globalnTracksReco[part];
2380  }
2381  //if (!calcAcceptance) {
2382  //printf("Type\tAll\tReco\tReco cor\n");
2383  printf("%10s\t%10s\t%10s\t%10s\n", "Type", "All", "Reco", "Reco cor");
2384  fout << "MC Tracks:" << endl;
2385  fout << "Particle\tTemperature\tError\tRadius\tError\tchi2/ndf\n";
2386  for (int part = 0; part < p_sz; ++part) {
2387  //if (part==0)
2388  {
2389  double tmpTMC = getTemperatureAll(globalmtavMC[part], part, 20);
2390  double tmpRMC =
2391  getRadius(tmpTMC,
2392  globalnTracksMC[part] / static_cast<double>(events),
2393  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass());
2394  double tmpTMCerr =
2395  getTemperatureDerivAll(globalmtavMC[part], part, 20)
2396  * sqrt((globalmt2avMC[part] - globalmtavMC[part] * globalmtavMC[part])
2397  / (globalnTracksMC[part] - 1.));
2398  double tmp1MC = 0., tmp2MC = 0.;
2399  double tmpNMCerr = sqrt(globalnTracksMC[part]) / events;
2400  tmp1MC = getRadiusDerivT(
2401  tmpTMC,
2402  static_cast<double>(globalnTracksMC[part]) / events,
2403  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass());
2404  tmp2MC = getRadiusDerivN(
2405  tmpTMC,
2406  static_cast<double>(globalnTracksMC[part]) / events,
2407  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass());
2408  double tmpRMCerr = sqrt(tmp1MC * tmp1MC * tmpTMCerr * tmpTMCerr
2409  + tmp2MC * tmp2MC * tmpNMCerr * tmpNMCerr);
2410  double tmpmt2th = ThermalMt2(
2411  tmpTMC, TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass());
2412  double tmpchi2ndf =
2413  (globalfavMC[part] - tmpmt2th) * (globalfavMC[part] - tmpmt2th)
2414  / ((globalf2avMC[part] - globalfavMC[part] * globalfavMC[part])
2415  / (globalnTracksMC[part] - 1.));
2416  fout << p_names[part].Data() << "\t" << tmpTMC << "\t" << tmpTMCerr
2417  << "\t" << tmpRMC << "\t" << tmpRMCerr << "\t" << tmpchi2ndf << "\n";
2418  //cout << globalmt2avMC[part] << " " << tmpmt2th << " " << (globalmt4avMC[part]-globalmt2avMC[part]*globalmt2avMC[part]) / (globalnTracksMC[part]-1.) << "\n";
2419  }
2420  }
2421  //for(int rl=0;rl<recoLevels;++rl) {
2422  //printf("%s\n", LevelNames[fRecoLevel].Data());
2423  fout << "Reco level " << fRecoLevel << ":" << endl;
2424  fout << "Particle\tTemperature\tError stat.\tError mom.\tError "
2425  "tot.\tRadius\tError stat.\tError mom.\tError tot.\tchi2/ndf\n";
2426  for (int part = 0; part < p_sz; ++part) {
2427  //if (part==0)
2428  {
2429 
2430  mtTacc[part] = mtTsts[part];
2431  Npart[part] = Npartsts[part];
2432  if (fusePID == 2 && fRecoLevel >= 1) {
2433  mtTacc[part] = mtTststof[part];
2434  Npart[part] = Npartststof[part];
2435  }
2436  printf(
2437  "%s Temperature\t%10lf\t%10lf\t%10lf\n",
2438  p_names[part].Data(),
2439  getTemperatureAll(globalmtavMC[part], part, 20),
2440  getTemperatureAll(globalmtavReco[part], part, 20),
2441  getTemperatureAllCor(globalmtavReco[part], part, 20, mtTacc[part]));
2442  double tmpT1 = getTemperatureAll(globalmtavMC[part], part, 20);
2443  double tmpT2 = getTemperatureAll(globalmtavReco[part], part, 20);
2444  double tmpT3 =
2445  getTemperatureAllCor(globalmtavReco[part], part, 20, mtTacc[part]);
2446  printf(
2447  "%s Radius\t%10lf\t%10lf\t%10lf\n",
2448  p_names[part].Data(),
2449  getRadius(tmpT1,
2450  globalnTracksMC[part] / static_cast<double>(events),
2451  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()),
2452  getRadius(tmpT2,
2453  globalnTracksReco[part] / static_cast<double>(events),
2454  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()),
2455  getRadiusCor(
2456  tmpT3,
2457  globalnTracksReco[part] / static_cast<double>(events),
2458  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass(),
2459  Npart[part]));
2460  printf("%lf\t%lf\n",
2461  Npart[part]->Eval(tmpT3),
2462  globalnTracksReco[part]
2463  / static_cast<double>(globalnTracksMC[part]));
2464  printf("%lf\t%lf\n", globalmtavMC[part], globalmtavReco[part]);
2465 
2466  double tmpmt2thMC = ThermalMt2(
2467  tmpT1, TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass());
2468  double tmpchi2ndfMC =
2469  (globalfavMC[part] - tmpmt2thMC) * (globalfavMC[part] - tmpmt2thMC)
2470  / ((globalf2avMC[part] - globalfavMC[part] * globalfavMC[part])
2471  / (globalnTracksMC[part] - 1.));
2472 
2473  double tmpTRe =
2474  getTemperatureAllCor(globalmtavReco[part], part, 20, mtTacc[part]);
2475  double tmpRRe = getRadiusCor(
2476  tmpTRe,
2477  globalnTracksReco[part] / static_cast<double>(events),
2478  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass(),
2479  Npart[part]);
2480  double tmpTReerr =
2481  getTemperatureDerivAllCor(globalmtavReco[part], part, 20, mtTacc[part])
2482  * sqrt(
2483  (globalmt2avReco[part] - globalmtavReco[part] * globalmtavReco[part])
2484  / (globalnTracksReco[part] - 1.));
2485  double tmp1Re = 0., tmp2Re = 0.;
2486  double tmpNReerr = sqrt(globalnTracksReco[part]) / events;
2487  tmp1Re = getRadiusDerivTCor(
2488  tmpTRe,
2489  static_cast<double>(globalnTracksReco[part]) / events,
2490  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass(),
2491  Npart[part]);
2492  tmp2Re = getRadiusDerivNCor(
2493  tmpTRe,
2494  static_cast<double>(globalnTracksReco[part]) / events,
2495  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass(),
2496  Npart[part]);
2497  double tmpRReerr = sqrt(tmp1Re * tmp1Re * tmpTReerr * tmpTReerr
2498  + tmp2Re * tmp2Re * tmpNReerr * tmpNReerr);
2499 
2500  if (0)
2501  cout << sqrt(globalmtmomerrReco[part]) / globalnTracksReco[part]
2502  << endl;
2503  double tmpTReerrmom =
2504  getTemperatureDerivAllCor(globalmtavReco[part], part, 20, mtTacc[part])
2505  * sqrt(globalmtmomerrReco[part]) / globalnTracksReco[part];
2506  double tmpRReerrmom = sqrt(tmp1Re * tmp1Re * tmpTReerrmom * tmpTReerrmom);
2507 
2508  double tmpmt2thRecoUn = ThermalMt2(
2509  tmpT2, TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass());
2510  double tmpchi2ndfRecoUn =
2511  (globalfavReco[part] - tmpmt2thRecoUn)
2512  * (globalfavReco[part] - tmpmt2thRecoUn)
2513  / ((globalf2avReco[part] - globalfavReco[part] * globalfavReco[part])
2514  / (globalnTracksReco[part] - 1.));
2515 
2516  double tmpmt2th = ThermalMt2Acc(
2517  tmpTRe,
2518  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass(),
2519  AcceptanceSTS);
2520  double tmpchi2ndf =
2521  (globalfavReco[part] - tmpmt2th) * (globalfavReco[part] - tmpmt2th)
2522  / ((globalf2avReco[part] - globalfavReco[part] * globalfavReco[part])
2523  / (globalnTracksReco[part] - 1.));
2524 
2525  fout << p_names[part].Data() << "\t" << tmpTRe << "\t" << tmpTReerr
2526  << "\t" << tmpTReerrmom << "\t"
2527  << sqrt(tmpTReerr * tmpTReerr + tmpTReerrmom * tmpTReerrmom) << "\t"
2528  << tmpRRe << "\t" << tmpRReerr << "\t" << tmpRReerrmom << "\t"
2529  << sqrt(tmpRReerr * tmpRReerr + tmpRReerrmom * tmpRReerrmom) << "\t"
2530  << tmpchi2ndf << "\n";
2531 
2532  printf("%lf\t%lf\t%lf\n", tmpchi2ndfMC, tmpchi2ndfRecoUn, tmpchi2ndf);
2533 
2534 
2535  ThermalDistributionFunction pl(
2536  part,
2537  tmpT1,
2538  getRadius(tmpT1,
2539  globalnTracksMC[part] / static_cast<double>(events),
2540  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()),
2541  ekin,
2542  NULL,
2543  NULL);
2544 
2545  for (int n = 0; n < hfyMC[part]->GetNbinsX(); n++) {
2546  hfyMCmodel[part]->SetBinContent(
2547  n,
2548  pl.dndylab(hfyMC[part]->GetBinCenter(n)) * static_cast<double>(events)
2549  * (hfyMC[part]->GetXaxis()->GetBinUpEdge(n)
2550  - hfyMC[part]->GetXaxis()->GetBinLowEdge(n)));
2551  }
2552 
2553  for (int nx = 0; nx < hfdndydptMC[part]->GetNbinsX(); nx++) {
2554  for (int ny = 0; ny < hfdndydptMC[part]->GetNbinsY(); ny++) {
2555  //hfdndydptMCmodel[part]->SetBinContent(nx, ny, pl.dndydpt(hfdndydptMC[part]->GetXaxis()->GetBinCenter(nx), hfdndydptMC[part]->GetYaxis()->GetBinCenter(ny)) * globalmtavMC[0]);
2556  hfdndydptMCmodel[part]->SetBinContent(
2557  nx,
2558  ny,
2559  pl.dndydptbin(hfdndydptMC[part]->GetXaxis()->GetBinLowEdge(nx),
2560  hfdndydptMC[part]->GetXaxis()->GetBinUpEdge(nx),
2561  10,
2562  hfdndydptMC[part]->GetYaxis()->GetBinLowEdge(ny),
2563  hfdndydptMC[part]->GetYaxis()->GetBinUpEdge(ny),
2564  10)
2565  * (double) events * hfdndydptMC[part]->GetXaxis()->GetBinWidth(nx)
2566  * hfdndydptMC[part]->GetYaxis()->GetBinWidth(ny));
2567  }
2568  }
2569 
2570  ThermalDistributionFunction plR(
2571  part,
2572  tmpT2,
2573  getRadius(tmpT2,
2574  globalnTracksReco[part] / static_cast<double>(events),
2575  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()),
2576  ekin,
2577  NULL,
2578  NULL);
2579 
2580  for (int n = 0; n < hfyReco[part]->GetNbinsX(); n++) {
2581  hfyRecomodel[part]->SetBinContent(
2582  n,
2583  plR.dndylab(hfyReco[part]->GetBinCenter(n))
2584  * static_cast<double>(events)
2585  * (hfyReco[part]->GetXaxis()->GetBinUpEdge(n)
2586  - hfyReco[part]->GetXaxis()->GetBinLowEdge(n)));
2587  }
2588 
2589  for (int nx = 0; nx < hfdndydptReco[part]->GetNbinsX(); nx++) {
2590  for (int ny = 0; ny < hfdndydptReco[part]->GetNbinsY(); ny++) {
2591  //hfdndydptMCmodel[part]->SetBinContent(nx, ny, pl.dndydpt(hfdndydptMC[part]->GetXaxis()->GetBinCenter(nx), hfdndydptMC[part]->GetYaxis()->GetBinCenter(ny)) * globalmtavMC[0]);
2592  hfdndydptRecomodel[part]->SetBinContent(
2593  nx,
2594  ny,
2595  plR.dndydptbin(hfdndydptReco[part]->GetXaxis()->GetBinLowEdge(nx),
2596  hfdndydptReco[part]->GetXaxis()->GetBinUpEdge(nx),
2597  10,
2598  hfdndydptReco[part]->GetYaxis()->GetBinLowEdge(ny),
2599  hfdndydptReco[part]->GetYaxis()->GetBinUpEdge(ny),
2600  10)
2601  * static_cast<double>(events)
2602  * hfdndydptReco[part]->GetXaxis()->GetBinWidth(nx)
2603  * hfdndydptReco[part]->GetYaxis()->GetBinWidth(ny));
2604  }
2605  }
2606 
2607  ReconstructionEfficiencyFunction reff;
2608  if (fTrackNumber == 0)
2609  reff = ReconstructionEfficiencyFunction(0.99, 0.98, 0.135);
2610  else
2611  reff = ReconstructionEfficiencyFunction(0.98, 0.88, 0.2);
2612 
2613  ThermalDistributionFunction plRC(
2614  part,
2615  tmpT3,
2616  getRadiusCor(
2617  tmpT3,
2618  globalnTracksReco[part] / static_cast<double>(events),
2619  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass(),
2620  Npart[part]),
2621  ekin,
2622  &AcceptanceSTS,
2623  &reff);
2624 
2625  for (int n = 0; n < hfyRecoCor[part]->GetNbinsX(); n++) {
2626  hfyRecoCormodel[part]->SetBinContent(
2627  n,
2628  plRC.dndylab(hfyRecoCor[part]->GetBinCenter(n))
2629  * static_cast<double>(events)
2630  * (hfyRecoCor[part]->GetXaxis()->GetBinUpEdge(n)
2631  - hfyRecoCor[part]->GetXaxis()->GetBinLowEdge(n)));
2632  //hfyRecoCormodel[part]->SetBinContent(n, plRC.dndybinlab(hfyRecoCor[part]->GetXaxis()->GetBinLowEdge(n), hfyRecoCor[part]->GetXaxis()->GetBinUpEdge(n), 10) * (double)events * (hfyRecoCor[part]->GetXaxis()->GetBinUpEdge(n) - hfyRecoCor[part]->GetXaxis()->GetBinLowEdge(n)));
2633  }
2634 
2635  for (int nx = 0; nx < hfdndydptRecoCor[part]->GetNbinsX(); nx++) {
2636  for (int ny = 0; ny < hfdndydptRecoCor[part]->GetNbinsY(); ny++) {
2637  //hfdndydptMCmodel[part]->SetBinContent(nx, ny, pl.dndydpt(hfdndydptMC[part]->GetXaxis()->GetBinCenter(nx), hfdndydptMC[part]->GetYaxis()->GetBinCenter(ny)) * globalmtavMC[0]);
2638  hfdndydptRecoCormodel[part]->SetBinContent(
2639  nx,
2640  ny,
2641  plRC.dndydptbin(
2642  hfdndydptRecoCor[part]->GetXaxis()->GetBinLowEdge(nx),
2643  hfdndydptRecoCor[part]->GetXaxis()->GetBinUpEdge(nx),
2644  10,
2645  hfdndydptRecoCor[part]->GetYaxis()->GetBinLowEdge(ny),
2646  hfdndydptRecoCor[part]->GetYaxis()->GetBinUpEdge(ny),
2647  10)
2648  * static_cast<double>(events)
2649  * hfdndydptRecoCor[part]->GetXaxis()->GetBinWidth(nx)
2650  * hfdndydptRecoCor[part]->GetYaxis()->GetBinWidth(ny));
2651  }
2652  }
2653 
2654  ThermalChi2Func fFCN2(
2655  hfyMC[part], hfdndydptMC[part], static_cast<double>(events));
2656  cout << "MC chi2/ndf = "
2657  << fFCN2.chi2dndy(
2658  part,
2659  tmpT1,
2660  getRadius(
2661  tmpT1,
2662  globalnTracksMC[part] / static_cast<double>(events),
2663  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()),
2664  ekin,
2665  NULL,
2666  NULL)
2667  << "\n";
2668  //fflush(stdout);
2669  cout << "MC YPtchi2/ndf = "
2670  << fFCN2.chi2ypt(
2671  part,
2672  tmpT1,
2673  getRadius(
2674  tmpT1,
2675  globalnTracksMC[part] / static_cast<double>(events),
2676  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()),
2677  ekin,
2678  NULL,
2679  NULL)
2680  << "\n";
2681 
2682  ThermalChi2Func fFCN2Reco(
2683  hfyReco[part], hfdndydptReco[part], static_cast<double>(events));
2684  cout << "Reco uncor chi2/ndf = "
2685  << fFCN2Reco.chi2dndy(
2686  part,
2687  tmpT2,
2688  getRadius(
2689  tmpT2,
2690  globalnTracksReco[part] / static_cast<double>(events),
2691  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()),
2692  ekin,
2693  NULL,
2694  NULL,
2695  0.05)
2696  << "\n";
2697  //fflush(stdout);
2698  cout << "Reco uncor YPtchi2/ndf = "
2699  << fFCN2Reco.chi2ypt(
2700  part,
2701  tmpT2,
2702  getRadius(
2703  tmpT2,
2704  globalnTracksReco[part] / static_cast<double>(events),
2705  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass()),
2706  ekin,
2707  NULL,
2708  NULL,
2709  0.05)
2710  << "\n";
2711 
2712  ThermalChi2Func fFCN2RecoCor(
2713  hfyRecoCor[part], hfdndydptRecoCor[part], static_cast<double>(events));
2714  cout << "Reco chi2/ndf = "
2715  << fFCN2RecoCor.chi2dndy(
2716  part,
2717  tmpT3,
2718  getRadiusCor(
2719  tmpT3,
2720  globalnTracksReco[part] / static_cast<double>(events),
2721  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass(),
2722  Npart[part]),
2723  ekin,
2724  &AcceptanceSTS,
2725  &reff,
2726  0.05)
2727  << "\n";
2728  //fflush(stdout);
2729  cout << "Reco YPtchi2/ndf = "
2730  << fFCN2RecoCor.chi2ypt(
2731  part,
2732  tmpT3,
2733  getRadiusCor(
2734  tmpT3,
2735  globalnTracksReco[part] / static_cast<double>(events),
2736  TDatabasePDG::Instance()->GetParticle(pdgIds[part])->Mass(),
2737  Npart[part]),
2738  ekin,
2739  &AcceptanceSTS,
2740  &reff,
2741  0.05)
2742  << "\n";
2743  }
2744  }
2745  //}
2746  //printf("Type\tRapidity\tAll\tReco\tReco cor\n");
2747  //}
2748 }
CbmThermalModelNoFlow::globalf2avReco
Double_t globalf2avReco[p_sz]
Definition: CbmThermalModelNoFlow.h:173
exp
friend F32vec4 exp(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:134
ThermalModelNoFlowNamespace::ThermalDistributionFunction::dndydptlab
double dndydptlab(double y, double pt) const
Definition: CbmThermalModelNoFlow.cxx:187
CbmThermalModelNoFlow::mtTststof
TSpline3 * mtTststof[p_sz]
Definition: CbmThermalModelNoFlow.h:149
CbmThermalModelNoFlow::hfyMC
TH1F * hfyMC[p_sz]
Definition: CbmThermalModelNoFlow.h:211
CbmThermalModelNoFlow::fTrackNumber
Int_t fTrackNumber
Definition: CbmThermalModelNoFlow.h:132
ThermalModelNoFlowNamespace::ThermalDistributionFunction::operator=
ThermalDistributionFunction & operator=(const ThermalDistributionFunction &)
CbmThermalModelNoFlow::ThermalMt
Double_t ThermalMt(double T, double m)
Definition: CbmThermalModelNoFlow.cxx:2018
CbmThermalModelNoFlow::flistStsTracksMatch
TClonesArray * flistStsTracksMatch
Definition: CbmThermalModelNoFlow.h:137
ThermalModelNoFlowNamespace::ThermalDistributionFunction
Definition: CbmThermalModelNoFlow.cxx:90
CbmVertex.h
CbmThermalModelNoFlow::ThermalMtAcc
Double_t ThermalMtAcc(double T, double m, const ThermalModelNoFlowNamespace::AcceptanceFunction &func)
Definition: CbmThermalModelNoFlow.cxx:2071
CbmThermalModelNoFlow::getInverseSlope
Double_t getInverseSlope(double mt, double mo, int iters)
Definition: CbmThermalModelNoFlow.cxx:2309
CbmThermalModelNoFlow::hfdndydptRecomodel
TH2F * hfdndydptRecomodel[p_sz]
Definition: CbmThermalModelNoFlow.h:219
ThermalModelNoFlowNamespace::ThermalChi2Func::ThermalChi2Func
ThermalChi2Func(const ThermalChi2Func &)
ThermalModelNoFlowNamespace::ThermalChi2Func::~ThermalChi2Func
~ThermalChi2Func()
Definition: CbmThermalModelNoFlow.cxx:371
CbmThermalModelNoFlow::hTempErrMomFullRecoCor
TH1F * hTempErrMomFullRecoCor[p_sz]
Definition: CbmThermalModelNoFlow.h:200
CbmThermalModelNoFlow::ReadAcceptanceFunction
void ReadAcceptanceFunction(ThermalModelNoFlowNamespace::AcceptanceFunction &func, TString filename)
Definition: CbmThermalModelNoFlow.cxx:877
ThermalModelNoFlowNamespace::ThermalDistributionFunction::ekin
double ekin
Definition: CbmThermalModelNoFlow.cxx:356
CbmThermalModelNoFlow::ReInit
virtual void ReInit(FairRootManager *fManger)
Definition: CbmThermalModelNoFlow.cxx:900
CbmThermalModelNoFlow::globalnTracksReco
Int_t globalnTracksReco[p_sz]
Definition: CbmThermalModelNoFlow.h:175
ThermalModelNoFlowNamespace::ThermalChi2Func::Norm
double Norm
Definition: CbmThermalModelNoFlow.cxx:450
ThermalModelNoFlowNamespace::ThermalDistributionFunction::ThermalDistributionFunction
ThermalDistributionFunction(int part, double T_, double R_, double ekin_, AcceptanceFunction *af_=NULL, ReconstructionEfficiencyFunction *rf_=NULL)
Definition: CbmThermalModelNoFlow.cxx:93
CbmThermalModelNoFlow::getRadiusDerivT
Double_t getRadiusDerivT(double T, double N, double mo)
Definition: CbmThermalModelNoFlow.cxx:2183
L1Field.h
CbmL1PFFitter
Definition: CbmL1PFFitter.h:31
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
ThermalModelNoFlowNamespace::ThermalDistributionFunction::dndylab
double dndylab(double y) const
Definition: CbmThermalModelNoFlow.cxx:136
CbmThermalModelNoFlow::ycm
Float_t ycm
Definition: CbmThermalModelNoFlow.h:128
CbmThermalModelNoFlow::Exec
virtual void Exec()
Definition: CbmThermalModelNoFlow.cxx:942
CbmThermalModelNoFlow::hfdndydptMC
TH2F * hfdndydptMC[p_sz]
Definition: CbmThermalModelNoFlow.h:213
CbmThermalModelNoFlow::hfdndydptReco
TH2F * hfdndydptReco[p_sz]
Definition: CbmThermalModelNoFlow.h:218
CbmThermalModelNoFlow::~CbmThermalModelNoFlow
~CbmThermalModelNoFlow()
Definition: CbmThermalModelNoFlow.cxx:875
CbmThermalModelNoFlow::tempDerivFunction
Double_t tempDerivFunction(double T, double dT, double mt, double mo)
Definition: CbmThermalModelNoFlow.cxx:2360
CbmThermalModelNoFlow::Npartststof
TSpline3 * Npartststof[p_sz]
Definition: CbmThermalModelNoFlow.h:149
CbmThermalModelNoFlow::hTempFullReco
TH1F * hTempFullReco[p_sz]
Definition: CbmThermalModelNoFlow.h:193
CbmThermalModelNoFlow::getDerivTemperature
Double_t getDerivTemperature(double mt, double mo, double T, int iters)
Definition: CbmThermalModelNoFlow.cxx:2316
ThermalModelNoFlowNamespace::ReconstructionEfficiencyFunction::f
Double_t f(double p)
Definition: CbmThermalModelNoFlow.h:51
CbmThermalModelNoFlow::getRadiusDerivNCor
Double_t getRadiusDerivNCor(double T, double N, double mo, TSpline3 *Np)
Definition: CbmThermalModelNoFlow.cxx:2212
CbmThermalModelNoFlow::hTempFullRecoCor
TH1F * hTempFullRecoCor[p_sz]
Definition: CbmThermalModelNoFlow.h:194
CbmThermalModelNoFlow::getTemperatureAllCor
Double_t getTemperatureAllCor(double mt, int part, int iters, TSpline3 *Np)
Definition: CbmThermalModelNoFlow.cxx:2262
CbmMCTrack::GetPdgCode
Int_t GetPdgCode() const
Definition: CbmMCTrack.h:70
CbmThermalModelNoFlow::hTempFullMC
TH1F * hTempFullMC[p_sz]
Definition: CbmThermalModelNoFlow.h:192
CbmThermalModelNoFlow::mtTacc
TSpline3 * mtTacc[p_sz]
Definition: CbmThermalModelNoFlow.h:148
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
ThermalModelNoFlowNamespace::p_names
const TString p_names[p_sz1]
Definition: CbmThermalModelNoFlow.cxx:69
CbmTrackMatch
Definition: CbmTrackMatch.h:18
CbmThermalModelNoFlow::globalnTracksMC
Int_t globalnTracksMC[p_sz]
Definition: CbmThermalModelNoFlow.h:169
CbmThermalModelNoFlow::flistTofHits
TClonesArray * flistTofHits
Definition: CbmThermalModelNoFlow.h:183
ThermalModelNoFlowNamespace::ThermalDistributionFunction::mass
double mass
Definition: CbmThermalModelNoFlow.cxx:355
CbmThermalModelNoFlow::getTemperatureDerivAllCor
Double_t getTemperatureDerivAllCor(double mt, int part, int iters, TSpline3 *Np)
Definition: CbmThermalModelNoFlow.cxx:2280
CbmGlobalTrack.h
ThermalModelNoFlowNamespace::ThermalDistributionFunction::T
double T() const
Definition: CbmThermalModelNoFlow.cxx:347
ThermalModelNoFlowNamespace::AcceptanceFunction::probs
std::vector< Double_t > probs
Definition: CbmThermalModelNoFlow.h:38
ThermalModelNoFlowNamespace
Definition: CbmThermalModelNoFlow.cxx:62
CbmGlobalTrack::GetLength
Double_t GetLength() const
Definition: CbmGlobalTrack.h:50
CbmThermalModelNoFlow::globalmt2avReco
Double_t globalmt2avReco[p_sz]
Definition: CbmThermalModelNoFlow.h:171
ThermalModelNoFlowNamespace::AcceptanceFunction::dy
Double_t dy
Definition: CbmThermalModelNoFlow.h:37
ThermalModelNoFlowNamespace::AcceptanceFunction
Definition: CbmThermalModelNoFlow.h:36
CbmThermalModelNoFlow::hRadErrStatFullReco
TH1F * hRadErrStatFullReco[p_sz]
Definition: CbmThermalModelNoFlow.h:205
CbmThermalModelNoFlow::RecEfficiency
Double_t RecEfficiency(Double_t p)
Definition: CbmThermalModelNoFlow.cxx:2063
CbmRichRing.h
CbmThermalModelNoFlow
Definition: CbmThermalModelNoFlow.h:57
CbmStsPoint
Definition: CbmStsPoint.h:27
ThermalModelNoFlowNamespace::ThermalChi2Func
Definition: CbmThermalModelNoFlow.cxx:362
CbmTrackMatch.h
CbmThermalModelNoFlow::hfyRecoCormodel
TH1F * hfyRecoCormodel[p_sz]
Definition: CbmThermalModelNoFlow.h:222
CbmThermalModelNoFlow::hfdndydptRecoCor
TH2F * hfdndydptRecoCor[p_sz]
Definition: CbmThermalModelNoFlow.h:223
CbmThermalModelNoFlow::checkIfReconstructable
bool checkIfReconstructable(CbmKFTrErrMCPoints *inTrack)
Definition: CbmThermalModelNoFlow.cxx:2147
CbmThermalModelNoFlow::AcceptanceSTS
ThermalModelNoFlowNamespace::AcceptanceFunction AcceptanceSTS
Definition: CbmThermalModelNoFlow.h:226
CbmThermalModelNoFlow::globalf2avMC
Double_t globalf2avMC[p_sz]
Definition: CbmThermalModelNoFlow.h:167
CbmKFTrErrMCPoints::TofArray
std::vector< CbmTofPoint * > TofArray
Definition: CbmKFTrErrMCPoints.h:79
CbmThermalModelNoFlow.h
ThermalModelNoFlowNamespace::ThermalDistributionFunction::fV
double fV
Definition: CbmThermalModelNoFlow.cxx:354
CbmThermalModelNoFlow::getTemperatureDerivAll
Double_t getTemperatureDerivAll(double mt, int part, int iters)
Definition: CbmThermalModelNoFlow.cxx:2255
CbmThermalModelNoFlow::hTempErrStatFullMC
TH1F * hTempErrStatFullMC[p_sz]
Definition: CbmThermalModelNoFlow.h:195
CbmThermalModelNoFlow::getTemperatureRapidity
Double_t getTemperatureRapidity(double mt, double mo, int iters, double y=0.)
Definition: CbmThermalModelNoFlow.cxx:2290
CbmThermalModelNoFlow::fChiToPrimVtx
std::vector< float > fChiToPrimVtx
Definition: CbmThermalModelNoFlow.h:145
ThermalModelNoFlowNamespace::ThermalDistributionFunction::fT
double fT
Definition: CbmThermalModelNoFlow.cxx:354
ThermalModelNoFlowNamespace::ThermalDistributionFunction::chi2FuncAv
double chi2FuncAv()
Definition: CbmThermalModelNoFlow.cxx:298
CbmL1Def.h
CbmThermalModelNoFlow::globalmtmomerrMC
Double_t globalmtmomerrMC[p_sz]
Definition: CbmThermalModelNoFlow.h:168
ThermalModelNoFlowNamespace::ThermalChi2Func::operator=
ThermalChi2Func & operator=(const ThermalChi2Func &)
CbmStsTrack.h
Data class for STS tracks.
CbmGlobalTrack::GetStsTrackIndex
Int_t GetStsTrackIndex() const
Definition: CbmGlobalTrack.h:38
CbmThermalModelNoFlow::events
Int_t events
Definition: CbmThermalModelNoFlow.h:190
CbmThermalModelNoFlow::hfyRecoCor
TH1F * hfyRecoCor[p_sz]
Definition: CbmThermalModelNoFlow.h:221
ThermalModelNoFlowNamespace::pdgIds
const int pdgIds[p_sz1]
Definition: CbmThermalModelNoFlow.cxx:73
CbmThermalModelNoFlow::hRadFullReco
TH1F * hRadFullReco[p_sz]
Definition: CbmThermalModelNoFlow.h:202
ThermalModelNoFlowNamespace::ThermalChi2Func::chi2ypt
double chi2ypt(int part, double T_, double R_, double ekin_, AcceptanceFunction *af_=NULL, ReconstructionEfficiencyFunction *rf_=NULL, double systerr=0.) const
Definition: CbmThermalModelNoFlow.cxx:408
CbmThermalModelNoFlow::Finish
virtual void Finish()
Definition: CbmThermalModelNoFlow.cxx:2369
CbmL1PFFitter::GetChiToVertex
void GetChiToVertex(std::vector< CbmStsTrack > &Tracks, std::vector< L1FieldRegion > &field, std::vector< float > &chiToVtx, CbmKFVertex &primVtx, float chiPrim=-1)
Definition: CbmL1PFFitter.cxx:403
ThermalModelNoFlowNamespace::ThermalDistributionFunction::V
double V() const
Definition: CbmThermalModelNoFlow.cxx:348
CbmKFTrErrMCPoints::IsReconstructable
Bool_t IsReconstructable(CbmMCTrack *mcTr, int MinNStations, int PerformanceMode, float MinRecoMom)
Definition: CbmKFTrErrMCPoints.cxx:168
CbmThermalModelNoFlow::fPrimVtx
CbmVertex * fPrimVtx
Definition: CbmThermalModelNoFlow.h:139
CbmTrackMatch::GetNofMCTracks
Int_t GetNofMCTracks() const
Definition: CbmTrackMatch.h:56
CbmThermalModelNoFlow::globalmtavReco
Double_t globalmtavReco[p_sz]
Definition: CbmThermalModelNoFlow.h:170
CbmThermalModelNoFlow::globalmtmomerrReco
Double_t globalmtmomerrReco[p_sz]
Definition: CbmThermalModelNoFlow.h:174
ThermalModelNoFlowNamespace::p_sz1
const int p_sz1
Definition: CbmThermalModelNoFlow.cxx:67
CbmHit::GetTime
Double_t GetTime() const
Definition: CbmHit.h:75
CbmThermalModelNoFlow::getTemperature
Double_t getTemperature(double mt, double mo, int iters)
Definition: CbmThermalModelNoFlow.cxx:2221
ThermalModelNoFlowNamespace::p_names_gr
const TString p_names_gr[p_sz1]
Definition: CbmThermalModelNoFlow.cxx:70
BilinearSplineFunction::Eval
double Eval(double x, double y) const
Definition: CbmBilinearSplineFunction.h:55
CbmKFTrErrMCPoints::GetStsPoint
CbmStsPoint * GetStsPoint(Int_t i)
Definition: CbmKFTrErrMCPoints.h:39
log
friend F32vec4 log(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:135
CbmVertex
Definition: CbmVertex.h:26
CbmThermalModelNoFlow::globalmtavMC
Double_t globalmtavMC[p_sz]
Definition: CbmThermalModelNoFlow.h:164
CbmThermalModelNoFlow::NFracAcc
Double_t NFracAcc(double T, double m, const ThermalModelNoFlowNamespace::AcceptanceFunction &func)
Definition: CbmThermalModelNoFlow.cxx:2120
CbmThermalModelNoFlow::hRadFullRecoCor
TH1F * hRadFullRecoCor[p_sz]
Definition: CbmThermalModelNoFlow.h:203
CbmKFTrErrMCPoints
Definition: CbmKFTrErrMCPoints.h:32
CbmThermalModelNoFlow::hfyRecomodel
TH1F * hfyRecomodel[p_sz]
Definition: CbmThermalModelNoFlow.h:217
CbmThermalModelNoFlow::flsitGlobalTracks
TClonesArray * flsitGlobalTracks
Definition: CbmThermalModelNoFlow.h:182
CbmL1PFFitter.h
CbmThermalModelNoFlow::tempCrit
Double_t tempCrit(double y, double pt, double m)
Definition: CbmThermalModelNoFlow.cxx:1937
CbmThermalModelNoFlow::tempFunction
Double_t tempFunction(double T, double mt, double mo)
Definition: CbmThermalModelNoFlow.cxx:2355
CbmThermalModelNoFlow::getRadiusDerivN
Double_t getRadiusDerivN(double T, double N, double mo)
Definition: CbmThermalModelNoFlow.cxx:2188
ThermalModelNoFlowNamespace::ThermalDistributionFunction::ycm
double ycm
Definition: CbmThermalModelNoFlow.cxx:357
CbmThermalModelNoFlow::ThermalMt2Acc
Double_t ThermalMt2Acc(double T, double m, const ThermalModelNoFlowNamespace::AcceptanceFunction &func)
Definition: CbmThermalModelNoFlow.cxx:2098
CbmThermalModelNoFlow::hTempErrStatFullRecoCor
TH1F * hTempErrStatFullRecoCor[p_sz]
Definition: CbmThermalModelNoFlow.h:197
CbmThermalModelNoFlow::Npart
TSpline3 * Npart[p_sz]
Definition: CbmThermalModelNoFlow.h:148
CbmKFTrErrMCPoints::GetNStsPoints
int GetNStsPoints() const
Definition: CbmKFTrErrMCPoints.h:50
xMath::Pi
double Pi()
Definition: xMath.h:5
ThermalModelNoFlowNamespace::recoLevels
const int recoLevels
Definition: CbmThermalModelNoFlow.cxx:76
CbmThermalModelNoFlow::getRadiusCor
Double_t getRadiusCor(double T, double N, double mo, TSpline3 *Np)
Definition: CbmThermalModelNoFlow.cxx:2193
CbmThermalModelNoFlow::hRadErrMomFullRecoCor
TH1F * hRadErrMomFullRecoCor[p_sz]
Definition: CbmThermalModelNoFlow.h:209
ThermalModelNoFlowNamespace::ThermalDistributionFunction::rf
ReconstructionEfficiencyFunction * rf
Definition: CbmThermalModelNoFlow.cxx:359
CbmThermalModelNoFlow::mtTall
TSpline3 * mtTall[p_sz]
Definition: CbmThermalModelNoFlow.h:148
ThermalModelNoFlowNamespace::ThermalDistributionFunction::tempCrit
static double tempCrit(double, double pt, double m)
Definition: CbmThermalModelNoFlow.cxx:245
CbmThermalModelNoFlow::hfdndydptRecoCormodel
TH2F * hfdndydptRecoCormodel[p_sz]
Definition: CbmThermalModelNoFlow.h:224
CbmThermalModelNoFlow::Init
virtual void Init()
Definition: CbmThermalModelNoFlow.cxx:938
CbmThermalModelNoFlow::flistTofPts
TClonesArray * flistTofPts
Definition: CbmThermalModelNoFlow.h:135
CbmThermalModelNoFlow::hRadErrStatFullMC
TH1F * hRadErrStatFullMC[p_sz]
Definition: CbmThermalModelNoFlow.h:204
ThermalModelNoFlowNamespace::ThermalDistributionFunction::dndycm
double dndycm(double y) const
Definition: CbmThermalModelNoFlow.cxx:114
CbmThermalModelNoFlow::getTemperatureAll
Double_t getTemperatureAll(double mt, int part, int iters)
Definition: CbmThermalModelNoFlow.cxx:2238
ThermalModelNoFlowNamespace::ThermalDistributionFunction::dndybinlab
double dndybinlab(double ymin, double ymax, int itery) const
Definition: CbmThermalModelNoFlow.cxx:162
CbmThermalModelNoFlow::flistStsTracks
TClonesArray * flistStsTracks
Definition: CbmThermalModelNoFlow.h:138
CbmThermalModelNoFlow::hTempErrMomFullReco
TH1F * hTempErrMomFullReco[p_sz]
Definition: CbmThermalModelNoFlow.h:199
CbmGlobalTrack
Definition: CbmGlobalTrack.h:26
CbmThermalModelNoFlow::ekin
Float_t ekin
Definition: CbmThermalModelNoFlow.h:127
GetCoefs2DLegendre32Legendre32
void GetCoefs2DLegendre32Legendre32(double ay, double by, double a2y, double b2y, std::vector< double > &xlag, std::vector< double > &wlag, std::vector< double > &xleg, std::vector< double > &wleg)
Definition: HRGModel/NumericalIntegration.h:318
ClassImp
ClassImp(CbmThermalModelNoFlow) CbmThermalModelNoFlow
Definition: CbmThermalModelNoFlow.cxx:460
ThermalModelNoFlowNamespace::ThermalChi2Func::chi2dndy
double chi2dndy(int part, double T_, double R_, double ekin_, AcceptanceFunction *af_=NULL, ReconstructionEfficiencyFunction *rf_=NULL, double systerr=0.) const
Definition: CbmThermalModelNoFlow.cxx:373
CbmThermalModelNoFlow::mtTsts
TSpline3 * mtTsts[p_sz]
Definition: CbmThermalModelNoFlow.h:148
ThermalModelNoFlowNamespace::GeVtoifm
const Double_t GeVtoifm
Definition: CbmThermalModelNoFlow.cxx:74
CbmThermalModelNoFlow::ComputeThermalDependence
void ComputeThermalDependence(Int_t part=0)
Definition: CbmThermalModelNoFlow.cxx:1941
CbmMCTrack.h
CbmThermalModelNoFlow::hfdndydptMCmodel
TH2F * hfdndydptMCmodel[p_sz]
Definition: CbmThermalModelNoFlow.h:214
ThermalModelNoFlowNamespace::AcceptanceFunction::setSpline
void setSpline()
Definition: CbmThermalModelNoFlow.h:40
ThermalModelNoFlowNamespace::ThermalDistributionFunction::dndydptcm
double dndydptcm(double y, double pt) const
Definition: CbmThermalModelNoFlow.cxx:173
v
__m128 v
Definition: L1/vectors/P4_F32vec4.h:1
CbmThermalModelNoFlow::ThermalMt2
Double_t ThermalMt2(double T, double m)
Definition: CbmThermalModelNoFlow.cxx:2039
CbmThermalModelNoFlow::getRadiusDerivTCor
Double_t getRadiusDerivTCor(double T, double N, double mo, TSpline3 *Np)
Definition: CbmThermalModelNoFlow.cxx:2203
ThermalModelNoFlowNamespace::AcceptanceFunction::sfunc
BilinearSplineFunction sfunc
Definition: CbmThermalModelNoFlow.h:39
ThermalModelNoFlowNamespace::ThermalDistributionFunction::dndydptbin
double dndydptbin(double ymin, double ymax, int itery, double ptmin, double ptmax, int iterpt) const
Definition: CbmThermalModelNoFlow.cxx:205
xMath::BesselK
double BesselK(int n, double x)
Definition: xMath.cxx:180
CbmMCTrack
Definition: CbmMCTrack.h:34
CbmThermalModelNoFlow::hTempErrMomFullMC
TH1F * hTempErrMomFullMC[p_sz]
Definition: CbmThermalModelNoFlow.h:198
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:60
ThermalModelNoFlowNamespace::ReconstructionEfficiencyFunction
Definition: CbmThermalModelNoFlow.h:45
CbmTrackMatch::GetMCTrackId
Int_t GetMCTrackId() const
Definition: CbmTrackMatch.h:44
CbmKFTrErrMCPoints::GetNTofPoints
int GetNTofPoints() const
Definition: CbmKFTrErrMCPoints.h:51
m
__m128 m
Definition: L1/vectors/P4_F32vec4.h:26
y
Double_t y
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmKFTrErrMCPoints::StsArray
std::vector< CbmStsPoint * > StsArray
Definition: CbmKFTrErrMCPoints.h:77
CbmThermalModelNoFlow::globalfavReco
Double_t globalfavReco[p_sz]
Definition: CbmThermalModelNoFlow.h:172
ThermalModelNoFlowNamespace::ThermalDistributionFunction::tempCritAv
double tempCritAv()
Definition: CbmThermalModelNoFlow.cxx:249
CbmThermalModelNoFlow::hRadErrMomFullReco
TH1F * hRadErrMomFullReco[p_sz]
Definition: CbmThermalModelNoFlow.h:208
CbmThermalModelNoFlow::Npartsts
TSpline3 * Npartsts[p_sz]
Definition: CbmThermalModelNoFlow.h:149
ThermalModelNoFlowNamespace::AcceptanceFunction::ys
std::vector< Double_t > ys
Definition: CbmThermalModelNoFlow.h:38
ThermalModelNoFlowNamespace::AcceptanceFunction::dpt
Double_t dpt
Definition: CbmThermalModelNoFlow.h:37
CbmThermalModelNoFlow::fusePID
Int_t fusePID
Definition: CbmThermalModelNoFlow.h:130
ThermalModelNoFlowNamespace::ThermalChi2Func::dndydptHist
TH2F * dndydptHist
Definition: CbmThermalModelNoFlow.cxx:452
CbmThermalModelNoFlow::getMtErrorSquare
Double_t getMtErrorSquare(double qp, double Tx, double Ty, double covMatrix[], double mo)
Definition: CbmThermalModelNoFlow.cxx:2335
CbmTofPoint
Geometric intersection of a MC track with a TOFb detector.
Definition: CbmTofPoint.h:40
CbmTofHit
Definition: core/data/tof/CbmTofHit.h:26
NumericalIntegration.h
CbmGlobalTrack::GetTofHitIndex
Int_t GetTofHitIndex() const
Definition: CbmGlobalTrack.h:42
CbmThermalModelNoFlow::hRadErrStatFullRecoCor
TH1F * hRadErrStatFullRecoCor[p_sz]
Definition: CbmThermalModelNoFlow.h:206
ThermalModelNoFlowNamespace::ThermalDistributionFunction::ThermalDistributionFunction
ThermalDistributionFunction(const ThermalDistributionFunction &)
CbmTrdTrack.h
CbmThermalModelNoFlow::globalmt2avMC
Double_t globalmt2avMC[p_sz]
Definition: CbmThermalModelNoFlow.h:165
ThermalModelNoFlowNamespace::ThermalChi2Func::ThermalChi2Func
ThermalChi2Func(TH1F *dndyexp, TH2F *dndydptexp, double Norm_)
Definition: CbmThermalModelNoFlow.cxx:365
CbmThermalModelNoFlow::globalfavMC
Double_t globalfavMC[p_sz]
Definition: CbmThermalModelNoFlow.h:166
kProtonMass
const Double_t kProtonMass
Definition: CbmPhsdGenerator.cxx:31
CbmStsTrack
Definition: CbmStsTrack.h:37
CbmThermalModelNoFlow::fRecoLevel
Int_t fRecoLevel
Definition: CbmThermalModelNoFlow.h:131
ThermalModelNoFlowNamespace::ThermalDistributionFunction::~ThermalDistributionFunction
~ThermalDistributionFunction()
Definition: CbmThermalModelNoFlow.cxx:112
ThermalModelNoFlowNamespace::ThermalDistributionFunction::af
AcceptanceFunction * af
Definition: CbmThermalModelNoFlow.cxx:358
CbmThermalModelNoFlow::p_sz
static const int p_sz
Definition: CbmThermalModelNoFlow.h:59
CbmThermalModelNoFlow::hTempErrStatFullReco
TH1F * hTempErrStatFullReco[p_sz]
Definition: CbmThermalModelNoFlow.h:196
CbmThermalModelNoFlow::chi2func
Double_t chi2func(double y, double pt, double m)
Definition: CbmThermalModelNoFlow.cxx:1932
ThermalModelNoFlowNamespace::ThermalChi2Func::dndyHist
TH1F * dndyHist
Definition: CbmThermalModelNoFlow.cxx:451
ThermalModelNoFlowNamespace::AcceptanceFunction::getAcceptance
Double_t getAcceptance(const Double_t &y, const Double_t &pt) const
Definition: CbmThermalModelNoFlow.cxx:82
CbmThermalModelNoFlow::hRadErrMomFullMC
TH1F * hRadErrMomFullMC[p_sz]
Definition: CbmThermalModelNoFlow.h:207
CbmKFVertex.h
CbmThermalModelNoFlow::AcceptanceSTSTOF
ThermalModelNoFlowNamespace::AcceptanceFunction AcceptanceSTSTOF
Definition: CbmThermalModelNoFlow.h:227
ThermalModelNoFlowNamespace::AcceptanceFunction::pts
std::vector< Double_t > pts
Definition: CbmThermalModelNoFlow.h:38
ThermalModelNoFlowNamespace::kProtonMass
const Double_t kProtonMass
Definition: CbmThermalModelNoFlow.cxx:66
CbmThermalModelNoFlow::flistStsPts
TClonesArray * flistStsPts
Definition: CbmThermalModelNoFlow.h:134
CbmThermalModelNoFlow::getRadius
Double_t getRadius(double T, double N, double mo)
Definition: CbmThermalModelNoFlow.cxx:2176
CbmThermalModelNoFlow::CbmThermalModelNoFlow
CbmThermalModelNoFlow(Float_t ekin_=24.08, Int_t recoLevel=3, Int_t usePID=1, Int_t trackNumber=1, Int_t iVerbose=1)
ThermalModelNoFlowNamespace::ThermalDistributionFunction::chi2func
static double chi2func(double, double pt, double m)
Definition: CbmThermalModelNoFlow.cxx:240
CbmThermalModelNoFlow::flistMCTracks
TClonesArray * flistMCTracks
Definition: CbmThermalModelNoFlow.h:181
CbmThermalModelNoFlow::hRadFullMC
TH1F * hRadFullMC[p_sz]
Definition: CbmThermalModelNoFlow.h:201
ThermalModelNoFlowNamespace::LevelNames
const TString LevelNames[recoLevels]
Definition: CbmThermalModelNoFlow.cxx:77
CbmThermalModelNoFlow::hfyMCmodel
TH1F * hfyMCmodel[p_sz]
Definition: CbmThermalModelNoFlow.h:212
CbmKFVertex
Definition: CbmKFVertex.h:6
CbmThermalModelNoFlow::hfyReco
TH1F * hfyReco[p_sz]
Definition: CbmThermalModelNoFlow.h:216