CbmRoot
CbmBoltzmannDistribution.cxx
Go to the documentation of this file.
1 /*
2  *====================================================================
3  *
4  * Extraction of temperature from Boltzmann distribution
5  *
6  * Authors: V.Vovchenko
7  *
8  * e-mail : new
9  *
10  *====================================================================
11  *
12  * Extraction of temperature from Boltzmann distribution
13  *
14  *====================================================================
15  */
16 
17 #include "CbmBoltzmannDistribution.h"
18 #include "CbmL1Def.h"
19 
20 
21 //#include "KFParticleFinder.h"
22 //#include "KFParticleSIMD.h"
23 #include "CbmKFTrack.h"
24 #include "CbmKFVertex.h"
25 #include "CbmStsTrack.h"
26 
27 #include "TClonesArray.h"
28 #include "TDirectory.h"
29 #include "TFile.h"
30 #include "TMath.h"
31 
32 #include "TCanvas.h"
33 #include "TGraphErrors.h"
34 #include "TH1F.h"
35 #include "TH2F.h"
36 #include "TSpline.h"
37 
38 #include "L1Field.h"
39 
40 #include "CbmVertex.h"
41 
42 #include "TStopwatch.h"
43 #include <cmath>
44 #include <cstdio>
45 #include <iostream>
46 #include <map>
47 
48 //#include "CbmTrackMatch.h"
49 #include "CbmMCTrack.h"
50 #include "CbmTrackMatchNew.h"
51 //for particle ID from global track
52 #include "CbmGlobalTrack.h"
53 #include "CbmRichRing.h"
54 #include "CbmTofHit.h"
55 #include "CbmTrdTrack.h"
56 #include "TDatabasePDG.h"
57 //for RICH identification
58 #include "TSystem.h"
59 // #include "CbmRichElectronIdAnn.h"
60 
61 #include "CbmL1PFFitter.h"
62 
63 #include "KFPTrackVector.h"
64 #include "KFParticleTopoReconstructor.h"
65 
66 #include "BoltzmannDistribution.h"
67 
68 
69 using std::ios;
70 using std::vector;
71 
72 using namespace std;
73 
74 ClassImp(CbmBoltzmannDistribution)
75 
76  CbmBoltzmannDistribution::CbmBoltzmannDistribution(
77  Int_t recoLevel,
78  Int_t /*iVerbose*/,
79  TString Mode,
80  Int_t PDG,
81  TString pname,
82  Int_t EventStats,
83  KFParticleTopoReconstructor* tr,
84  Float_t ekin_)
85  : CbmModelBase(tr)
86  , name(pname)
87  ,
88  //ekin(ekin_),
89  //fusePID(usePID),
90  //fTrackNumber(trackNumber),
91  //flistStsTracks(0),
92  //flistStsTracksMatch(0),
93  //fPrimVtx(0),
94  ekin(ekin_)
95  , p0cm(5.)
96  , ycm(2.)
97  , fUpdate(true)
98  , fusePID(true)
99  , fRecoLevel(recoLevel)
100  , fTrackNumber(1)
101  , fEventStats(EventStats)
102  , events(0)
103  , fModeName(Mode)
104  , outfileName("")
105  ,
106  //flsitGlobalTracks(0),
107  //flistTofHits(0),
108  histodir(0)
109  , flistMCTracks(0)
110  , IndexT(0)
111  , IndexMt(0)
112  , IndexModelMt(0)
113  , IndexMt2(0)
114  , IndexModelMt2(0)
115  , IndexModelMt4Pi(0)
116  , histodndy(0)
117  , histodndymodel(0)
118  , histo1DIntervals(0)
119  , grTy(0)
120  , grdndyReco(0)
121  , pullT(0)
122  , Ts()
123  , kProtonMass(0.938271998)
124  , fPDGID(PDG)
125  , fMass(TDatabasePDG::Instance()->GetParticle(fPDGID)->Mass())
126  , fYminv()
127  , fYmaxv()
128  , paramGlobal(0)
129  , paramGlobalInterval(0)
130  , param2GlobalInterval(0)
131  , paramLocal(0)
132  , paramLocalInterval(0)
133  , totalLocal(0)
134  , totalGlobal(0)
135  , totalGlobalInterval(0)
136  , totalLocalInterval(0)
137  , totalEvents(0)
138  , model(0)
139  , modelmc(0)
140  , modelsY(0)
141 
142 // flistRichRings(0),
143 // flistTrdTracks(0),
144 
145 {
146  // fModeName = Mode;
147  // fEventStats = EventStats;
148  // fMass = TDatabasePDG::Instance()->GetParticle(fPDGID)->Mass();
149 
150  // events = 0;
151  Ts.resize(0);
152 
153  //PPDG = 2212;
154  // kProtonMass = 0.938271998;
155 
156  //PDGtoIndex.clear();
157 
158  TDirectory* currentDir = gDirectory;
159 
160  //std::cout << 0.128 << "\t" << mtTall->Eval(0.128) << endl;
161  //gDirectory->mkdir("KFModelParameters");
162  gDirectory->cd("Models");
163 
164  histodir = gDirectory;
165 
166  char ccc[200];
167  sprintf(ccc, "BoltzmannDistribution %s", name.Data());
168  gDirectory->mkdir(ccc);
169  gDirectory->cd(ccc);
170  gDirectory->mkdir(fModeName);
171  gDirectory->cd(fModeName);
172  TString tname = "PerEvent";
173  if (fEventStats != 1)
174  tname =
175  TString("Each ") + TString::Itoa(fEventStats, 10) + TString(" events");
176  gDirectory->mkdir(tname);
177  gDirectory->cd(tname);
178  //for(int part=0;part<p_sz;++part) {
179  //gDirectory->mkdir(p_names[part]);
180  //gDirectory->cd(p_names[part]);
181  int CurrentIndex = 0;
182 
183  IndexT = CurrentIndex;
184  histo1D[CurrentIndex] = new TH1F("T", "Event-by-event T", 100, 0., 0.4);
185  histo1D[CurrentIndex]->SetXTitle("T (GeV)");
186  histo1D[CurrentIndex]->SetYTitle("Entries");
187  CurrentIndex++;
188 
189  IndexMt = CurrentIndex;
190  histo1D[CurrentIndex] = new TH1F("f(m_{T})", "mt distribution", 200, 0., 2.5);
191  histo1D[CurrentIndex]->SetXTitle("m_{T} - m_{0} (GeV)");
192  histo1D[CurrentIndex]->SetYTitle("Entries");
193  CurrentIndex++;
194 
195  IndexModelMt = CurrentIndex;
196  histo1D[CurrentIndex] =
197  new TH1F("Model f(m_{T})", "Model mt distribution", 200, 0., 2.5);
198  histo1D[CurrentIndex]->SetXTitle("m_{T} - m_{0} (GeV)");
199  histo1D[CurrentIndex]->SetYTitle("Entries");
200  CurrentIndex++;
201 
202  IndexMt2 = CurrentIndex;
203  histo1D[CurrentIndex] =
204  new TH1F("f2(m_{T})", "mt2 distribution", 40, 0., 2.5);
205  histo1D[CurrentIndex]->SetXTitle("m_{T} - m_{0} (GeV)");
206  histo1D[CurrentIndex]->SetYTitle(
207  "#frac{1}{N} #frac{dN}{m_{T} dm_{T}} [(GeV)^{-2}]");
208  CurrentIndex++;
209 
210  IndexModelMt2 = CurrentIndex;
211  histo1D[CurrentIndex] =
212  new TH1F("Model f2(m_{T})", "Model mt2 distribution", 200, 0., 2.5);
213  histo1D[CurrentIndex]->SetXTitle("m_{T} - m_{0} (GeV)");
214  histo1D[CurrentIndex]->SetYTitle(
215  "#frac{1}{N} #frac{dN}{m_{T} dm_{T}} [(GeV)^{-2}]");
216  CurrentIndex++;
217 
218  IndexModelMt4Pi = CurrentIndex;
219  histo1D[CurrentIndex] =
220  new TH1F("Model f2(m_{T}) 4Pi", "Model mt2 distribution", 200, 0., 2.5);
221  histo1D[CurrentIndex]->SetXTitle("m_{T} - m_{0} (GeV)");
222  histo1D[CurrentIndex]->SetYTitle(
223  "#frac{1}{N} #frac{dN}{m_{T} dm_{T}} [(GeV)^{-2}]");
224  CurrentIndex++;
225 
226  histodndy = new TH1F("dN/dy", "Rapidity distribution", 40, -3., 3.);
227  histodndy->SetXTitle("y");
228  histodndy->SetYTitle("dN/dy");
229 
230  histodndymodel =
231  new TH1F("dN/dy model", "Rapidity distribution from model", 40, -3., 3.);
232  histodndymodel->SetXTitle("y");
233  histodndymodel->SetYTitle("dN/dy");
234 
235  grdndyReco = new TGraphErrors();
236  grdndyReco->SetTitle("Rapidity distribution");
237  grdndyReco->GetXaxis()->SetTitle("y");
238  grdndyReco->GetYaxis()->SetTitle("dN/dy");
239  grdndyReco->SetName(TString("dNdy-") + fModeName);
240  gDirectory->Add(grdndyReco);
241 
242 
243  if (fYminv.size() > 0) {
244  histo1DIntervals = new TH1F**[fYminv.size()];
245  for (unsigned int ind = 0; ind < fYminv.size(); ++ind) {
246  histo1DIntervals[ind] = new TH1F*[nHisto1D];
247  char cc[200], cc2[200];
248  CurrentIndex = 0;
249 
250  sprintf(cc, "Event-by-event T, %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
251  sprintf(cc2, "T, %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
252  histo1DIntervals[ind][CurrentIndex] = new TH1F(cc2, cc, 100, 0., 1.);
253  histo1DIntervals[ind][CurrentIndex]->SetXTitle("T (GeV)");
254  histo1DIntervals[ind][CurrentIndex]->SetYTitle("Entries");
255  CurrentIndex++;
256 
257  sprintf(cc, "mt distribution, %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
258  sprintf(cc2, "f(m_{T}, %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
259  histo1DIntervals[ind][CurrentIndex] = new TH1F(cc2, cc, 200, 0., 2.5);
260  histo1DIntervals[ind][CurrentIndex]->SetXTitle("m_{T} - m_{0} (GeV)");
261  histo1DIntervals[ind][CurrentIndex]->SetYTitle("Entries");
262  CurrentIndex++;
263 
264  sprintf(
265  cc, "Model mt distribution, %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
266  sprintf(cc2, "Model f(m_{T}), %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
267  histo1DIntervals[ind][CurrentIndex] = new TH1F(cc2, cc, 200, 0., 2.5);
268  histo1DIntervals[ind][CurrentIndex]->SetXTitle("m_{T} - m_{0} (GeV)");
269  histo1DIntervals[ind][CurrentIndex]->SetYTitle("Entries");
270  CurrentIndex++;
271 
272  sprintf(cc, "mt2 distribution, %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
273  sprintf(cc2, "f2(m_{T}, %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
274  histo1DIntervals[ind][CurrentIndex] = new TH1F(cc2, cc, 40, 0., 2.5);
275  histo1DIntervals[ind][CurrentIndex]->SetXTitle("m_{T} - m_{0} (GeV)");
276  histo1DIntervals[ind][CurrentIndex]->SetYTitle(
277  "#frac{1}{N} #frac{dN}{m_{T} dm_{T}} [(GeV)^{-2}]");
278  CurrentIndex++;
279 
280  sprintf(
281  cc, "Model mt2 distribution, %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
282  sprintf(cc2, "Model f2(m_{T}), %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
283  histo1DIntervals[ind][CurrentIndex] = new TH1F(cc2, cc, 200, 0., 2.5);
284  histo1DIntervals[ind][CurrentIndex]->SetXTitle("m_{T} - m_{0} (GeV)");
285  histo1DIntervals[ind][CurrentIndex]->SetYTitle(
286  "#frac{1}{N} #frac{dN}{m_{T} dm_{T}} [(GeV)^{-2}]");
287  CurrentIndex++;
288 
289  sprintf(cc,
290  "Model mt2 distribution 4Pi, %.2lf<y<%.2lf",
291  fYminv[ind],
292  fYmaxv[ind]);
293  sprintf(
294  cc2, "Model f2(m_{T}) 4Pi, %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
295  histo1DIntervals[ind][CurrentIndex] = new TH1F(cc2, cc, 200, 0., 2.5);
296  histo1DIntervals[ind][CurrentIndex]->SetXTitle("m_{T} - m_{0} (GeV)");
297  histo1DIntervals[ind][CurrentIndex]->SetYTitle(
298  "#frac{1}{N} #frac{dN}{m_{T} dm_{T}} [(GeV)^{-2}]");
299  CurrentIndex++;
300  }
301  }
302 
303 
304  gDirectory->cd("..");
305  gDirectory->cd("..");
306  gDirectory->cd("..");
307 
308  gDirectory = currentDir;
309 
310  double pbeam = sqrt((kProtonMass + ekin) * (kProtonMass + ekin)
311  - kProtonMass * kProtonMass);
312  double betacm = pbeam / (2. * kProtonMass + ekin);
313  ycm = 0.5 * log((1. + betacm) / (1. - betacm));
314 
315  double ssqrt = sqrt(2. * kProtonMass * (ekin + 2. * kProtonMass));
316  p0cm = sqrt(0.25 * ssqrt * ssqrt - kProtonMass * kProtonMass);
317 
318  std::cout << "ekin = " << ekin << "\n";
319  std::cout << "ycm = " << ycm << "\n";
320 
321  events = 0;
322 
323  paramGlobal = paramLocal = 0.;
324  totalGlobal = totalLocal = 0;
325  paramGlobalInterval.resize(0);
326  param2GlobalInterval.resize(0);
327  paramLocalInterval.resize(0);
328  totalGlobalInterval.resize(0);
329  totalLocalInterval.resize(0);
330  modelsY.resize(0);
331 
332  if (fRecoLevel == -1 || fRecoLevel > 10)
333  model = new BoltzmannDistribution(fMass);
334  else
335  model = new BoltzmannDistribution(fMass, fPDGID, true, -ycm, ycm, ycm, 1.);
336  modelmc = new BoltzmannDistribution(fMass);
337 }
338 
339 CbmBoltzmannDistribution::~CbmBoltzmannDistribution() {
340  if (model != NULL) delete model;
341  for (unsigned int ind = 0; ind < fYminv.size(); ++ind)
342  if (modelsY[ind] != NULL) delete modelsY[ind];
343 }
344 
345 void CbmBoltzmannDistribution::AddRapidityInterval(double ymin, double ymax) {
346  fYminv.push_back(ymin);
347  fYmaxv.push_back(ymax);
348  paramGlobalInterval.push_back(0.);
349  param2GlobalInterval.push_back(0.);
350  paramLocalInterval.push_back(0.);
351  totalGlobalInterval.push_back(0);
352  totalLocalInterval.push_back(0);
353  if (fRecoLevel == -1 || fRecoLevel > 10)
354  modelsY.push_back(new BoltzmannDistribution(fMass));
355  else
356  modelsY.push_back(
357  new BoltzmannDistribution(fMass, fPDGID, true, ymin, ymax, ycm, 1.));
358 }
359 
360 void CbmBoltzmannDistribution::AddHistos() {
361  TDirectory* currentDir = gDirectory;
362 
363  gDirectory->cd("Models");
364 
365  histodir = gDirectory;
366 
367  char ccc[200];
368  sprintf(ccc, "BoltzmannDistribution %s", name.Data());
369  //gDirectory->mkdir(ccc);
370  gDirectory->cd(ccc);
371  //gDirectory->mkdir(fModeName);
372  gDirectory->cd(fModeName);
373  TString tname = "PerEvent";
374  if (fEventStats != 1)
375  tname =
376  TString("Each ") + TString::Itoa(fEventStats, 10) + TString(" events");
377  gDirectory->cd(tname);
378  int CurrentIndex = 0;
379 
380  if (fYminv.size() > 0) {
381  histo1DIntervals = new TH1F**[fYminv.size()];
382  for (unsigned int ind = 0; ind < fYminv.size(); ++ind) {
383  char cc3[200];
384  sprintf(cc3, "%.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
385  gDirectory->mkdir(cc3);
386  gDirectory->cd(cc3);
387 
388  histo1DIntervals[ind] = new TH1F*[nHisto1D];
389  char cc[200], cc2[200];
390  CurrentIndex = 0;
391 
392  sprintf(cc, "Event-by-event T, %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
393  sprintf(cc2, "T, %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
394  histo1DIntervals[ind][CurrentIndex] = new TH1F(cc2, cc, 100, 0., 1.);
395  histo1DIntervals[ind][CurrentIndex]->SetXTitle("T (GeV)");
396  histo1DIntervals[ind][CurrentIndex]->SetYTitle("Entries");
397  CurrentIndex++;
398 
399  sprintf(cc, "mt distribution, %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
400  sprintf(cc2, "f(m_{T}, %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
401  histo1DIntervals[ind][CurrentIndex] = new TH1F(cc2, cc, 200, 0., 2.5);
402  histo1DIntervals[ind][CurrentIndex]->SetXTitle("m_{T} - m_{0} (GeV)");
403  histo1DIntervals[ind][CurrentIndex]->SetYTitle("Entries");
404  CurrentIndex++;
405 
406  sprintf(
407  cc, "Model mt distribution, %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
408  sprintf(cc2, "Model f(m_{T}), %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
409  histo1DIntervals[ind][CurrentIndex] = new TH1F(cc2, cc, 200, 0., 2.5);
410  histo1DIntervals[ind][CurrentIndex]->SetXTitle("m_{T} - m_{0} (GeV)");
411  histo1DIntervals[ind][CurrentIndex]->SetYTitle("Entries");
412  CurrentIndex++;
413 
414  sprintf(cc, "mt2 distribution, %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
415  sprintf(cc2, "f2(m_{T}, %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
416  histo1DIntervals[ind][CurrentIndex] = new TH1F(cc2, cc, 40, 0., 2.5);
417  histo1DIntervals[ind][CurrentIndex]->SetXTitle("m_{T} - m_{0} (GeV)");
418  histo1DIntervals[ind][CurrentIndex]->SetYTitle(
419  "#frac{1}{N} #frac{dN}{m_{T} dm_{T}} [(GeV)^{-2}]");
420  CurrentIndex++;
421 
422  sprintf(
423  cc, "Model mt2 distribution, %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
424  sprintf(cc2, "Model f2(m_{T}), %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
425  histo1DIntervals[ind][CurrentIndex] = new TH1F(cc2, cc, 200, 0., 2.5);
426  histo1DIntervals[ind][CurrentIndex]->SetXTitle("m_{T} - m_{0} (GeV)");
427  histo1DIntervals[ind][CurrentIndex]->SetYTitle(
428  "#frac{1}{N} #frac{dN}{m_{T} dm_{T}} [(GeV)^{-2}]");
429  CurrentIndex++;
430 
431  sprintf(cc,
432  "Model mt2 distribution 4Pi, %.2lf<y<%.2lf",
433  fYminv[ind],
434  fYmaxv[ind]);
435  sprintf(
436  cc2, "Model f2(m_{T}) 4Pi, %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
437  histo1DIntervals[ind][CurrentIndex] = new TH1F(cc2, cc, 200, 0., 2.5);
438  histo1DIntervals[ind][CurrentIndex]->SetXTitle("m_{T} - m_{0} (GeV)");
439  histo1DIntervals[ind][CurrentIndex]->SetYTitle(
440  "#frac{1}{N} #frac{dN}{m_{T} dm_{T}} [(GeV)^{-2}]");
441  CurrentIndex++;
442 
443  gDirectory->cd("..");
444  }
445  grTy = new TGraphErrors();
446  grTy->SetTitle("T(y)");
447  grTy->GetXaxis()->SetTitle("y");
448  grTy->GetYaxis()->SetTitle("T(y) (MeV)");
449  grTy->SetName(TString("T(y)-") + fModeName);
450  gDirectory->Add(grTy);
451  }
452 
453 
454  gDirectory->cd("..");
455  gDirectory->cd("..");
456  gDirectory->cd("..");
457 
458  gDirectory = currentDir;
459 }
460 
461 void CbmBoltzmannDistribution::ReInit(FairRootManager* fManger) {
462  flistMCTracks = dynamic_cast<TClonesArray*>(fManger->GetObject("MCTrack"));
463 }
464 
465 void CbmBoltzmannDistribution::Init() {
466  //ReInit();
467 }
468 
469 void CbmBoltzmannDistribution::Exec() {
470  if (fRecoLevel == -1 && !flistMCTracks) return;
471  if (fRecoLevel != -1 && !fTopoReconstructor) return;
472 
473  CalculateAveragesInEvent(fRecoLevel, 1);
474 
475  events++;
476  totalEvents++;
477 
478  //if (events%10==0) std::cout << "Event # " << events << "\n";
479 
480  if (events % fEventStats == 0) {
481 
482  double T = model->GetT(paramLocal / totalLocal);
483 
484  histo1D[IndexT]->Fill(T);
485 
486  events = 0;
487  paramLocal = 0.;
488  totalLocal = 0;
489 
490  for (unsigned int i = 0; i < fYminv.size(); ++i) {
491  if (totalLocalInterval[i] > 0)
492  histo1DIntervals[i][IndexT]->Fill(
493  modelsY[i]->GetT(paramLocalInterval[i] / totalLocalInterval[i]));
494  //std::cout << events << " " << paramLocalInterval[i] / totalLocalInterval[i] << "\n";
495  paramLocalInterval[i] = 0.;
496  totalLocalInterval[i] = 0;
497  }
498  }
499 }
500 
501 void CbmBoltzmannDistribution::Finish() {
502 
503  double T = model->GetT(paramGlobal / totalGlobal);
504  std::cout << "T = " << T << "\t<m_T> = " << paramGlobal / totalGlobal << "\n";
505 
506  histo1D[IndexMt]->Sumw2();
507  histo1D[IndexMt2]->Sumw2();
508  histo1D[IndexMt]->Scale(1. / totalGlobal
509  / histo1D[IndexMt]->GetXaxis()->GetBinWidth(1));
510  histo1D[IndexMt2]->Scale(1. / totalGlobal
511  / histo1D[IndexMt2]->GetXaxis()->GetBinWidth(1));
512  for (int n = 1; n < histo1D[IndexMt2]->GetNbinsX(); n++) {
513  histo1D[IndexMt2]->SetBinContent(
514  n,
515  histo1D[IndexMt2]->GetBinContent(n)
516  / (histo1D[IndexMt2]->GetXaxis()->GetBinCenter(n) + fMass));
517  histo1D[IndexMt2]->SetBinError(
518  n,
519  histo1D[IndexMt2]->GetBinError(n)
520  / (histo1D[IndexMt2]->GetXaxis()->GetBinCenter(n) + fMass));
521  }
522  histodndy->Scale(1. / totalEvents / histodndy->GetXaxis()->GetBinWidth(1));
523 
524  for (int n = 1; n < histo1D[IndexModelMt]->GetNbinsX(); n++) {
525  histo1D[IndexModelMt]->SetBinContent(
526  n,
527  model->fmt(histo1D[IndexModelMt]->GetXaxis()->GetBinCenter(n) + fMass,
528  T));
529  histo1D[IndexModelMt2]->SetBinContent(
530  n,
531  model->fmt(histo1D[IndexModelMt2]->GetXaxis()->GetBinCenter(n) + fMass, T)
532  / (histo1D[IndexModelMt2]->GetXaxis()->GetBinCenter(n) + fMass));
533  histo1D[IndexModelMt4Pi]->SetBinContent(
534  n,
535  modelmc->fmt(
536  histo1D[IndexModelMt4Pi]->GetXaxis()->GetBinCenter(n) + fMass, T)
537  / (histo1D[IndexModelMt4Pi]->GetXaxis()->GetBinCenter(n) + fMass));
538  }
539 
540  std::vector<double> ys(0), dndys(0), dndyerrs(0);
541  int grindex = 0;
542  for (unsigned int ind = 0; ind < fYminv.size(); ++ind) {
543  if (totalGlobalInterval[ind] > 0)
544  T =
545  modelsY[ind]->GetT(paramGlobalInterval[ind] / totalGlobalInterval[ind]);
546  else
547  continue;
548  histo1DIntervals[ind][IndexMt]->Sumw2();
549  histo1DIntervals[ind][IndexMt2]->Sumw2();
550  histo1DIntervals[ind][IndexMt]->Scale(
551  1. / totalGlobalInterval[ind]
552  / histo1DIntervals[ind][IndexMt]->GetXaxis()->GetBinWidth(1));
553  histo1DIntervals[ind][IndexMt2]->Scale(
554  1. / totalGlobalInterval[ind]
555  / histo1DIntervals[ind][IndexMt2]->GetXaxis()->GetBinWidth(1));
556  for (int n = 1; n < histo1DIntervals[ind][IndexMt2]->GetNbinsX(); n++) {
557  histo1DIntervals[ind][IndexMt2]->SetBinContent(
558  n,
559  histo1DIntervals[ind][IndexMt2]->GetBinContent(n)
560  / (histo1DIntervals[ind][IndexMt2]->GetXaxis()->GetBinCenter(n)
561  + fMass));
562  histo1DIntervals[ind][IndexMt2]->SetBinError(
563  n,
564  histo1DIntervals[ind][IndexMt2]->GetBinError(n)
565  / (histo1DIntervals[ind][IndexMt2]->GetXaxis()->GetBinCenter(n)
566  + fMass));
567  }
568  double avmt = 0.;
569  if (T == T && totalGlobalInterval[ind] > 0) {
570  for (int n = 1; n < histo1DIntervals[ind][IndexModelMt]->GetNbinsX();
571  n++) {
572  histo1DIntervals[ind][IndexModelMt]->SetBinContent(
573  n,
574  modelsY[ind]->fmt(
575  histo1DIntervals[ind][IndexModelMt]->GetXaxis()->GetBinCenter(n)
576  + fMass,
577  T));
578  histo1DIntervals[ind][IndexModelMt2]->SetBinContent(
579  n,
580  modelsY[ind]->fmt(
581  histo1DIntervals[ind][IndexModelMt2]->GetXaxis()->GetBinCenter(n)
582  + fMass,
583  T)
584  / (histo1DIntervals[ind][IndexModelMt2]->GetXaxis()->GetBinCenter(n)
585  + fMass));
586  histo1DIntervals[ind][IndexModelMt4Pi]->SetBinContent(
587  n,
588  modelmc->fmt(
589  histo1DIntervals[ind][IndexModelMt4Pi]->GetXaxis()->GetBinCenter(n)
590  + fMass,
591  T)
592  / (histo1DIntervals[ind][IndexModelMt4Pi]->GetXaxis()->GetBinCenter(
593  n)
594  + fMass));
595  avmt +=
596  (histo1DIntervals[ind][IndexModelMt]->GetXaxis()->GetBinCenter(n)
597  + fMass)
598  * modelsY[ind]->fmt(
599  histo1DIntervals[ind][IndexModelMt]->GetXaxis()->GetBinCenter(n)
600  + fMass,
601  T);
602  }
603  }
604  avmt *= histo1DIntervals[ind][IndexModelMt]->GetXaxis()->GetBinWidth(1);
605  std::cout << fYminv[ind] << "<y<" << fYmaxv[ind] << "\tT = " << T
606  << "\t<m_T> = "
607  << paramGlobalInterval[ind] / totalGlobalInterval[ind]
608  << "\t<m_T>_2 = " << avmt << "\n";
609  double s1 = paramGlobalInterval[ind];
610  double s2 = param2GlobalInterval[ind];
611  int Numb = totalGlobalInterval[ind];
612  double errT =
613  modelsY[ind]->dTdmt(paramGlobalInterval[ind] / totalGlobalInterval[ind])
614  * TMath::Sqrt((s2 - s1 * s1 / Numb) / (Numb - 1.) / Numb);
615  //std::cout << T << " " << errT << " " << (s2 - s1*s1 / Numb) << "\n";
616  if (T != T || errT != errT || errT > T / 3.) continue;
617  grTy->SetPoint(grindex, 0.5 * (fYminv[ind] + fYmaxv[ind]), T * 1.e3);
618  grTy->SetPointError(
619  grindex, 0. * 0.5 * (fYmaxv[ind] - fYminv[ind]), errT * 1.e3);
620  //std::cout << T << " " << errT << "\n";
621 
622  double A = modelsY[ind]->GetA(
623  totalGlobalInterval[ind] / static_cast<double>(totalEvents), T);
624  double errA = modelsY[ind]->GetAerror(
625  totalGlobalInterval[ind] / static_cast<double>(totalEvents),
626  T,
627  sqrt(totalGlobalInterval[ind]) / static_cast<double>(totalEvents),
628  errT);
629  std::cout << "A = " << A << " error = " << errA << "\n";
630  grdndyReco->SetPoint(grindex,
631  0.5 * (fYminv[ind] + fYmaxv[ind]),
632  A / (fYmaxv[ind] - fYminv[ind]));
633  grdndyReco->SetPointError(grindex,
634  0. * 0.5 * (fYmaxv[ind] - fYminv[ind]),
635  errA / (fYmaxv[ind] - fYminv[ind]));
636 
637  ys.push_back(0.5 * (fYminv[ind] + fYmaxv[ind]));
638  dndys.push_back(A / (fYmaxv[ind] - fYminv[ind]));
639  dndyerrs.push_back(errA / (fYmaxv[ind] - fYminv[ind]));
640 
641  grindex++;
642  }
643 
644  T = model->GetT(paramGlobal / totalGlobal);
645 
646  for (int n = 1; n < histodndymodel->GetNbinsX(); n++) {
647  histodndymodel->SetBinContent(
648  n,
649  model->dndy(
650  histodndymodel->GetXaxis()->GetBinCenter(n),
651  model->GetA(totalGlobal / static_cast<double>(totalEvents), T),
652  T));
653  }
654 }
655 
656 void CbmBoltzmannDistribution::CalculateAveragesInEvent(int RecoLevel,
657  bool UpdateGlobal) {
658  if (RecoLevel == -1) {
659  vector<CbmMCTrack> vRTracksMC;
660  int nTracksMC = flistMCTracks->GetEntries();
661  std::cout << "MC tracks: " << nTracksMC << "\n";
662  vRTracksMC.resize(nTracksMC);
663  for (int iTr = 0; iTr < nTracksMC; iTr++)
664  // vRTracksMC[iTr] = *( (CbmMCTrack*) flistMCTracks->At(iTr));
665  vRTracksMC[iTr] = *(static_cast<CbmMCTrack*>(flistMCTracks->At(iTr)));
666 
667  for (int iTr = 0; iTr < nTracksMC; iTr++) {
668  if (vRTracksMC[iTr].GetPdgCode() == fPDGID
669  && vRTracksMC[iTr].GetMotherId() == -1) {
670  totalLocal++;
671  double pt = vRTracksMC[iTr].GetPt();
672  double ty = vRTracksMC[iTr].GetRapidity() - ycm;
673  double mt = TMath::Sqrt(
674  vRTracksMC[iTr].GetMass() * vRTracksMC[iTr].GetMass() + pt * pt);
675  paramLocal += mt;
676  if (UpdateGlobal) {
677  totalGlobal++;
678  paramGlobal += mt;
679  }
680  histo1D[IndexMt]->Fill(mt - fMass);
681  histo1D[IndexMt2]->Fill(mt - fMass);
682  histodndy->Fill(ty);
683  for (unsigned int ind = 0; ind < fYminv.size(); ++ind) {
684  if (ty >= fYminv[ind] && ty <= fYmaxv[ind]) {
685  if (UpdateGlobal) totalGlobalInterval[ind]++;
686  totalLocalInterval[ind]++;
687  if (UpdateGlobal) {
688  paramGlobalInterval[ind] += mt;
689  param2GlobalInterval[ind] += mt * mt;
690  }
691  paramLocalInterval[ind] += mt;
692  histo1DIntervals[ind][IndexMt]->Fill(mt - fMass);
693  histo1DIntervals[ind][IndexMt2]->Fill(mt - fMass);
694  }
695  }
696  }
697  }
698  } else {
699  for (int itype = 2; itype <= 3; ++itype) {
700  const KFPTrackVector& tr = fTopoReconstructor->GetTracks()[itype];
701  const kfvector_int& pdgs = tr.PDG();
702  for (unsigned int ind = 0; ind < pdgs.size(); ++ind) {
703  int iPDG = pdgs[ind];
704  //for(unsigned int part=0;part<ParticlePDGsTrack.size();++part) {
705  if (iPDG == fPDGID) {
706  totalLocal++;
707  double pt = tr.Pt(ind);
708  double p = tr.P(ind);
709  double m = TDatabasePDG::Instance()->GetParticle(iPDG)->Mass();
710  double p0 = TMath::Sqrt(m * m + p * p);
711  double pz = TMath::Sqrt(p * p - pt * pt);
712  double ty = 0.5 * log((p0 + pz) / (p0 - pz)) - ycm;
713  pz = TMath::Sqrt(m * m + pt * pt) * TMath::SinH(ty);
714  double mt = TMath::Sqrt(m * m + pt * pt);
715 
716  paramLocal += mt;
717  if (UpdateGlobal) {
718  totalGlobal++;
719  paramGlobal += mt;
720  }
721  histo1D[IndexMt]->Fill(mt - fMass);
722  histo1D[IndexMt2]->Fill(mt - fMass);
723  histodndy->Fill(ty);
724  for (unsigned int ind2 = 0; ind2 < fYminv.size(); ++ind2) {
725  if (ty >= fYminv[ind2] && ty <= fYmaxv[ind2]) {
726  if (UpdateGlobal) totalGlobalInterval[ind2]++;
727  totalLocalInterval[ind2]++;
728  if (UpdateGlobal) {
729  paramGlobalInterval[ind2] += mt;
730  param2GlobalInterval[ind2] += mt * mt;
731  }
732  paramLocalInterval[ind2] += mt;
733  histo1DIntervals[ind2][IndexMt]->Fill(mt - fMass);
734  histo1DIntervals[ind2][IndexMt2]->Fill(mt - fMass);
735  }
736  }
737  }
738  //}
739  }
740  }
741  }
742 }
CbmBoltzmannDistribution::totalLocal
int totalLocal
Definition: CbmBoltzmannDistribution.h:112
BoltzmannDistribution::GetT
double GetT(double amt)
Definition: BoltzmannDistribution.h:37
CbmModelBase::fTopoReconstructor
KFParticleTopoReconstructor * fTopoReconstructor
Definition: CbmModelBase.h:42
CbmBoltzmannDistribution::paramLocal
double paramLocal
Definition: CbmBoltzmannDistribution.h:110
CbmVertex.h
CbmBoltzmannDistribution::~CbmBoltzmannDistribution
~CbmBoltzmannDistribution()
Definition: CbmBoltzmannDistribution.cxx:339
CbmBoltzmannDistribution::flistMCTracks
TClonesArray * flistMCTracks
Definition: CbmBoltzmannDistribution.h:85
BoltzmannDistribution::fmt
double fmt(double amt, double T)
Definition: BoltzmannDistribution.cxx:110
L1Field.h
CbmBoltzmannDistribution::modelsY
std::vector< BoltzmannDistribution * > modelsY
Definition: CbmBoltzmannDistribution.h:121
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
CbmBoltzmannDistribution::Finish
virtual void Finish()
Definition: CbmBoltzmannDistribution.cxx:501
CbmBoltzmannDistribution::totalGlobalInterval
std::vector< int > totalGlobalInterval
Definition: CbmBoltzmannDistribution.h:113
CbmBoltzmannDistribution::paramGlobal
double paramGlobal
Definition: CbmBoltzmannDistribution.h:108
CbmBoltzmannDistribution::ReInit
virtual void ReInit(FairRootManager *fManger)
Definition: CbmBoltzmannDistribution.cxx:461
CbmBoltzmannDistribution::paramLocalInterval
std::vector< double > paramLocalInterval
Definition: CbmBoltzmannDistribution.h:111
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmBoltzmannDistribution::fYmaxv
std::vector< double > fYmaxv
Definition: CbmBoltzmannDistribution.h:107
CbmGlobalTrack.h
BoltzmannDistribution.h
CbmBoltzmannDistribution::fYminv
std::vector< double > fYminv
Definition: CbmBoltzmannDistribution.h:107
CbmBoltzmannDistribution::totalLocalInterval
std::vector< int > totalLocalInterval
Definition: CbmBoltzmannDistribution.h:113
CbmRichRing.h
CbmBoltzmannDistribution::fEventStats
Int_t fEventStats
Definition: CbmBoltzmannDistribution.h:76
CbmL1Def.h
CbmBoltzmannDistribution::IndexMt2
int IndexMt2
Definition: CbmBoltzmannDistribution.h:90
BoltzmannDistribution::GetA
double GetA(double multiplicity, double T)
Definition: BoltzmannDistribution.h:39
CbmStsTrack.h
Data class for STS tracks.
CbmBoltzmannDistribution::param2GlobalInterval
std::vector< double > param2GlobalInterval
Definition: CbmBoltzmannDistribution.h:109
CbmBoltzmannDistribution::grdndyReco
TGraphErrors * grdndyReco
Definition: CbmBoltzmannDistribution.h:97
CbmBoltzmannDistribution::Init
virtual void Init()
Definition: CbmBoltzmannDistribution.cxx:465
BoltzmannDistribution::dndy
double dndy(double y, double A, double T)
Definition: BoltzmannDistribution.cxx:170
CbmBoltzmannDistribution::Exec
virtual void Exec()
Definition: CbmBoltzmannDistribution.cxx:469
ClassImp
ClassImp(CbmBoltzmannDistribution) CbmBoltzmannDistribution
Definition: CbmBoltzmannDistribution.cxx:74
CbmBoltzmannDistribution::CalculateAveragesInEvent
void CalculateAveragesInEvent(int RecoLevel, bool UpdateGlobal=0)
Definition: CbmBoltzmannDistribution.cxx:656
CbmBoltzmannDistribution::IndexT
int IndexT
Definition: CbmBoltzmannDistribution.h:90
CbmBoltzmannDistribution::fPDGID
int fPDGID
Definition: CbmBoltzmannDistribution.h:105
CbmBoltzmannDistribution::histodir
TDirectory * histodir
Definition: CbmBoltzmannDistribution.h:83
CbmBoltzmannDistribution::ycm
Float_t ycm
Definition: CbmBoltzmannDistribution.h:71
CbmBoltzmannDistribution::IndexModelMt2
int IndexModelMt2
Definition: CbmBoltzmannDistribution.h:90
log
friend F32vec4 log(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:135
CbmTrackMatchNew.h
CbmL1PFFitter.h
CbmBoltzmannDistribution::histo1DIntervals
TH1F *** histo1DIntervals
Definition: CbmBoltzmannDistribution.h:94
CbmBoltzmannDistribution::totalEvents
int totalEvents
Definition: CbmBoltzmannDistribution.h:114
CbmBoltzmannDistribution::name
TString name
Definition: CbmBoltzmannDistribution.h:41
CbmKFTrack.h
CbmBoltzmannDistribution::histodndy
TH1F * histodndy
Definition: CbmBoltzmannDistribution.h:93
CbmModelBase
Definition: CbmModelBase.h:25
CbmBoltzmannDistribution::modelmc
BoltzmannDistribution * modelmc
Definition: CbmBoltzmannDistribution.h:120
CbmBoltzmannDistribution::histodndymodel
TH1F * histodndymodel
Definition: CbmBoltzmannDistribution.h:93
CbmBoltzmannDistribution::AddRapidityInterval
void AddRapidityInterval(double ymin, double ymax)
Definition: CbmBoltzmannDistribution.cxx:345
CbmBoltzmannDistribution::fRecoLevel
Int_t fRecoLevel
Definition: CbmBoltzmannDistribution.h:74
CbmBoltzmannDistribution::paramGlobalInterval
std::vector< double > paramGlobalInterval
Definition: CbmBoltzmannDistribution.h:109
CbmBoltzmannDistribution::nHisto1D
static const int nHisto1D
Definition: CbmBoltzmannDistribution.h:89
CbmMCTrack.h
BoltzmannDistribution
Definition: BoltzmannDistribution.h:12
CbmBoltzmannDistribution::fModeName
TString fModeName
Definition: CbmBoltzmannDistribution.h:80
CbmMCTrack
Definition: CbmMCTrack.h:34
CbmBoltzmannDistribution::model
BoltzmannDistribution * model
Definition: CbmBoltzmannDistribution.h:120
CbmBoltzmannDistribution::IndexModelMt4Pi
int IndexModelMt4Pi
Definition: CbmBoltzmannDistribution.h:90
CbmBoltzmannDistribution::IndexMt
int IndexMt
Definition: CbmBoltzmannDistribution.h:90
CbmBoltzmannDistribution::totalGlobal
int totalGlobal
Definition: CbmBoltzmannDistribution.h:112
m
__m128 m
Definition: L1/vectors/P4_F32vec4.h:26
CbmBoltzmannDistribution.h
CbmBoltzmannDistribution
Definition: CbmBoltzmannDistribution.h:38
CbmBoltzmannDistribution::fMass
double fMass
Definition: CbmBoltzmannDistribution.h:106
CbmTrdTrack.h
kProtonMass
const Double_t kProtonMass
Definition: CbmPhsdGenerator.cxx:31
CbmBoltzmannDistribution::CbmBoltzmannDistribution
CbmBoltzmannDistribution(Int_t recoLevel=-1, Int_t iVerbose=1, TString Mode="MC", Int_t PDG=-211, TString pname="pi-", Int_t EventStats=1, KFParticleTopoReconstructor *tr=0, Float_t ekin_=25.)
CbmBoltzmannDistribution::histo1D
TH1F * histo1D[nHisto1D]
Definition: CbmBoltzmannDistribution.h:92
CbmKFVertex.h
CbmBoltzmannDistribution::grTy
TGraphErrors * grTy
Definition: CbmBoltzmannDistribution.h:96
CbmBoltzmannDistribution::AddHistos
void AddHistos()
Definition: CbmBoltzmannDistribution.cxx:360
CbmBoltzmannDistribution::IndexModelMt
int IndexModelMt
Definition: CbmBoltzmannDistribution.h:90
CbmBoltzmannDistribution::events
Int_t events
Definition: CbmBoltzmannDistribution.h:77