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