27 #include "TClonesArray.h"
28 #include "TDirectory.h"
33 #include "TGraphErrors.h"
42 #include "TStopwatch.h"
54 #include "CbmTofHit.h"
56 #include "TDatabasePDG.h"
63 #include "KFPTrackVector.h"
64 #include "KFParticleTopoReconstructor.h"
84 KFParticleTopoReconstructor* tr,
97 , fRecoLevel(recoLevel)
99 , fEventStats(EventStats)
117 , histo1DIntervals(0)
131 , fMass(TDatabasePDG::Instance()->GetParticle(fPDGID)->Mass())
135 , paramGlobalInterval(0)
136 , param2GlobalInterval(0)
138 , paramLocalInterval(0)
141 , totalGlobalInterval(0)
142 , totalLocalInterval(0)
170 ycm = 0.5 *
log((1. + betacm) / (1. - betacm));
174 TDirectory* currentDir = gDirectory;
176 gDirectory->cd(
"Models");
178 histodir = gDirectory;
181 sprintf(ccc,
"BlastWave %s", name.Data());
182 gDirectory->mkdir(ccc);
184 gDirectory->mkdir(fModeName);
185 gDirectory->cd(fModeName);
186 TString tname =
"PerEvent";
187 if (fEventStats != 1)
189 TString(
"Each ") + TString::Itoa(fEventStats, 10) + TString(
" events");
190 gDirectory->mkdir(tname);
191 gDirectory->cd(tname);
192 int CurrentIndex = 0;
194 IndexT = CurrentIndex;
195 histo1D[CurrentIndex] =
new TH1F(
"T",
"Event-by-event T", 100, 0., 0.4);
196 histo1D[CurrentIndex]->SetXTitle(
"T (GeV)");
197 histo1D[CurrentIndex]->SetYTitle(
"Entries");
200 IndexMt = CurrentIndex;
201 histo1D[CurrentIndex] =
new TH1F(
"f(m_{T})",
"mt distribution", 200, 0., 2.5);
202 histo1D[CurrentIndex]->SetXTitle(
"m_{T} - m_{0} (GeV)");
203 histo1D[CurrentIndex]->SetYTitle(
"Entries");
206 IndexModelMt = CurrentIndex;
207 histo1D[CurrentIndex] =
208 new TH1F(
"Model f(m_{T})",
"Model mt distribution", 200, 0., 2.5);
209 histo1D[CurrentIndex]->SetXTitle(
"m_{T} - m_{0} (GeV)");
210 histo1D[CurrentIndex]->SetYTitle(
"Entries");
213 IndexMt2 = CurrentIndex;
214 histo1D[CurrentIndex] =
215 new TH1F(
"f2(m_{T})",
"mt2 distribution", 40, 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}]");
221 IndexModelMt2 = CurrentIndex;
222 histo1D[CurrentIndex] =
223 new TH1F(
"Model f2(m_{T})",
"Model mt2 distribution", 200, 0., 2.5);
224 histo1D[CurrentIndex]->SetXTitle(
"m_{T} - m_{0} (GeV)");
225 histo1D[CurrentIndex]->SetYTitle(
226 "#frac{1}{N} #frac{dN}{m_{T} dm_{T}} [(GeV)^{-2}]");
229 IndexModelMt4Pi = CurrentIndex;
230 histo1D[CurrentIndex] =
231 new TH1F(
"Model f2(m_{T}) 4Pi",
"Model mt2 distribution", 200, 0., 2.5);
232 histo1D[CurrentIndex]->SetXTitle(
"m_{T} - m_{0} (GeV)");
233 histo1D[CurrentIndex]->SetYTitle(
234 "#frac{1}{N} #frac{dN}{m_{T} dm_{T}} [(GeV)^{-2}]");
238 new TH1F(
"dN/dy",
"Rapidity distribution", 40, -2. * ycm, 2. * ycm);
239 histodndy->SetXTitle(
"y");
240 histodndy->SetYTitle(
"dN/dy");
242 histodndymodel =
new TH1F(
243 "dN/dy model",
"Rapidity distribution from model", 40, -2. * ycm, 2. * ycm);
244 histodndymodel->SetXTitle(
"y");
245 histodndymodel->SetYTitle(
"dN/dy");
247 histoeta =
new TH1F(
"etamax",
"Event-by-event etamax", 50, 0., ycm);
248 histoeta->SetXTitle(
"#eta_{max}");
249 histoeta->SetYTitle(
"Entries");
251 histodndymodel2 =
new TH1F(
"dN/dy model 2",
252 "Rapidity distribution from model",
256 histodndymodel2->SetXTitle(
"y");
257 histodndymodel2->SetYTitle(
"dN/dy");
260 new TH1F(
"Multiplicity",
"Event-by-event multiplicity", 1000, 0.5, 1000.5);
261 histomult->SetXTitle(
"Event");
262 histomult->SetYTitle(
"N");
264 histomultmodel =
new TH1F(
"Multiplicity from model",
265 "Event-by-event multiplicity",
269 histomultmodel->SetXTitle(
"Event");
270 histomultmodel->SetYTitle(
"N");
272 grdndyReco =
new TGraphErrors();
273 grdndyReco->SetTitle(
"Rapidity distribution");
274 grdndyReco->GetXaxis()->SetTitle(
"y");
275 grdndyReco->GetYaxis()->SetTitle(
"dN/dy");
276 grdndyReco->SetName(TString(
"dNdy-") + fModeName);
277 gDirectory->Add(grdndyReco);
280 if (fYminv.size() > 0) {
281 histo1DIntervals =
new TH1F**[fYminv.size()];
282 for (
unsigned int ind = 0; ind < fYminv.size(); ++ind) {
283 histo1DIntervals[ind] =
new TH1F*[nHisto1D];
284 char cc[200], cc2[200];
287 sprintf(cc,
"Event-by-event T, %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
288 sprintf(cc2,
"T, %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
289 histo1DIntervals[ind][CurrentIndex] =
new TH1F(cc2, cc, 100, 0., 1.);
290 histo1DIntervals[ind][CurrentIndex]->SetXTitle(
"T (GeV)");
291 histo1DIntervals[ind][CurrentIndex]->SetYTitle(
"Entries");
294 sprintf(cc,
"mt distribution, %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
295 sprintf(cc2,
"f(m_{T}, %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
296 histo1DIntervals[ind][CurrentIndex] =
new TH1F(cc2, cc, 200, 0., 2.5);
297 histo1DIntervals[ind][CurrentIndex]->SetXTitle(
"m_{T} - m_{0} (GeV)");
298 histo1DIntervals[ind][CurrentIndex]->SetYTitle(
"Entries");
302 cc,
"Model mt distribution, %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
303 sprintf(cc2,
"Model f(m_{T}), %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
304 histo1DIntervals[ind][CurrentIndex] =
new TH1F(cc2, cc, 200, 0., 2.5);
305 histo1DIntervals[ind][CurrentIndex]->SetXTitle(
"m_{T} - m_{0} (GeV)");
306 histo1DIntervals[ind][CurrentIndex]->SetYTitle(
"Entries");
309 sprintf(cc,
"mt2 distribution, %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
310 sprintf(cc2,
"f2(m_{T}, %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
311 histo1DIntervals[ind][CurrentIndex] =
new TH1F(cc2, cc, 40, 0., 2.5);
312 histo1DIntervals[ind][CurrentIndex]->SetXTitle(
"m_{T} - m_{0} (GeV)");
313 histo1DIntervals[ind][CurrentIndex]->SetYTitle(
314 "#frac{1}{N} #frac{dN}{m_{T} dm_{T}} [(GeV)^{-2}]");
318 cc,
"Model mt2 distribution, %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
319 sprintf(cc2,
"Model f2(m_{T}), %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
320 histo1DIntervals[ind][CurrentIndex] =
new TH1F(cc2, cc, 200, 0., 2.5);
321 histo1DIntervals[ind][CurrentIndex]->SetXTitle(
"m_{T} - m_{0} (GeV)");
322 histo1DIntervals[ind][CurrentIndex]->SetYTitle(
323 "#frac{1}{N} #frac{dN}{m_{T} dm_{T}} [(GeV)^{-2}]");
327 "Model mt2 distribution 4Pi, %.2lf<y<%.2lf",
331 cc2,
"Model f2(m_{T}) 4Pi, %.2lf<y<%.2lf", fYminv[ind], fYmaxv[ind]);
332 histo1DIntervals[ind][CurrentIndex] =
new TH1F(cc2, cc, 200, 0., 2.5);
333 histo1DIntervals[ind][CurrentIndex]->SetXTitle(
"m_{T} - m_{0} (GeV)");
334 histo1DIntervals[ind][CurrentIndex]->SetYTitle(
335 "#frac{1}{N} #frac{dN}{m_{T} dm_{T}} [(GeV)^{-2}]");
341 gDirectory->cd(
"..");
342 gDirectory->cd(
"..");
343 gDirectory->cd(
"..");
345 gDirectory = currentDir;
351 std::cout <<
"ekin = " << ekin <<
"\n";
352 std::cout <<
"ycm = " << ycm <<
"\n";
356 paramGlobal = paramLocal = 0.;
357 y2Global = y2Local = 0.;
358 y4Global = y4Local = 0.;
359 totalGlobal = totalLocal = 0;
360 paramGlobalInterval.resize(0);
361 param2GlobalInterval.resize(0);
362 paramLocalInterval.resize(0);
363 totalGlobalInterval.resize(0);
364 totalLocalInterval.resize(0);
368 if (fRecoLevel == -1 || fRecoLevel > 10)
369 model =
new BlastWave(fMass, fPDGID,
false, -2. * ycm, 2. * ycm, ycm, 0.);
371 model =
new BlastWave(fMass, fPDGID,
true, -2. * ycm, 2. * ycm, ycm, 0.);
372 modelmc =
new BlastWave(fMass, fPDGID,
true, -2. * ycm, 2. * ycm, ycm, 0.);
374 if (fRecoLevel == -1 || fRecoLevel > 10)
376 fMass, fPDGID,
false, -2. * ycm, 2. * ycm, ycm, Tlong);
379 fMass, fPDGID,
true, -2. * ycm, 2. * ycm, ycm, Tlong);
381 fMass, fPDGID,
false, -2. * ycm, 2. * ycm, ycm, Tlong);
388 for (
unsigned int ind = 0; ind <
fYminv.size(); ++ind)
407 TDirectory* currentDir = gDirectory;
409 gDirectory->cd(
"Models");
414 sprintf(ccc,
"BlastWave %s",
name.Data());
419 TString tname =
"PerEvent";
422 TString(
"Each ") + TString::Itoa(
fEventStats, 10) + TString(
" events");
424 gDirectory->cd(tname);
425 int CurrentIndex = 0;
429 for (
unsigned int ind = 0; ind <
fYminv.size(); ++ind) {
431 sprintf(cc3,
"%.2lf<y<%.2lf",
fYminv[ind],
fYmaxv[ind]);
432 gDirectory->mkdir(cc3);
436 char cc[200], cc2[200];
439 sprintf(cc,
"Event-by-event T, %.2lf<y<%.2lf",
fYminv[ind],
fYmaxv[ind]);
440 sprintf(cc2,
"T, %.2lf<y<%.2lf",
fYminv[ind],
fYmaxv[ind]);
446 sprintf(cc,
"mt distribution, %.2lf<y<%.2lf",
fYminv[ind],
fYmaxv[ind]);
447 sprintf(cc2,
"f(m_{T}, %.2lf<y<%.2lf",
fYminv[ind],
fYmaxv[ind]);
454 cc,
"Model mt distribution, %.2lf<y<%.2lf",
fYminv[ind],
fYmaxv[ind]);
455 sprintf(cc2,
"Model f(m_{T}), %.2lf<y<%.2lf",
fYminv[ind],
fYmaxv[ind]);
461 sprintf(cc,
"mt2 distribution, %.2lf<y<%.2lf",
fYminv[ind],
fYmaxv[ind]);
462 sprintf(cc2,
"f2(m_{T}, %.2lf<y<%.2lf",
fYminv[ind],
fYmaxv[ind]);
466 "#frac{1}{N} #frac{dN}{m_{T} dm_{T}} [(GeV)^{-2}]");
470 cc,
"Model mt2 distribution, %.2lf<y<%.2lf",
fYminv[ind],
fYmaxv[ind]);
471 sprintf(cc2,
"Model f2(m_{T}), %.2lf<y<%.2lf",
fYminv[ind],
fYmaxv[ind]);
475 "#frac{1}{N} #frac{dN}{m_{T} dm_{T}} [(GeV)^{-2}]");
479 "Model mt2 distribution 4Pi, %.2lf<y<%.2lf",
483 cc2,
"Model f2(m_{T}) 4Pi, %.2lf<y<%.2lf",
fYminv[ind],
fYmaxv[ind]);
487 "#frac{1}{N} #frac{dN}{m_{T} dm_{T}} [(GeV)^{-2}]");
490 gDirectory->cd(
"..");
492 grTy =
new TGraphErrors();
493 grTy->SetTitle(
"T(y)");
494 grTy->GetXaxis()->SetTitle(
"y");
495 grTy->GetYaxis()->SetTitle(
"T(y) (MeV)");
497 gDirectory->Add(
grTy);
501 gDirectory->cd(
"..");
502 gDirectory->cd(
"..");
503 gDirectory->cd(
"..");
505 gDirectory = currentDir;
509 flistMCTracks =
dynamic_cast<TClonesArray*
>(fManger->GetObject(
"MCTrack"));
550 * TMath::Sqrt((s2 - s1 * s1 / Numb) / (Numb - 1.) / Numb);
567 for (
unsigned int i = 0;
i <
fYminv.size(); ++
i) {
618 std::vector<double> ys(0), dndys(0), dndyerrs(0);
620 for (
unsigned int ind = 0; ind <
fYminv.size(); ++ind) {
685 std::cout <<
fYminv[ind] <<
"<y<" <<
fYmaxv[ind] <<
"\tT = " << T
688 <<
"\t<m_T>_2 = " << avmt <<
"\n";
694 * TMath::Sqrt((s2 - s1 * s1 / Numb) / (Numb - 1.) / Numb);
696 if (T != T || errT != errT || errT > T / 3.)
continue;
699 grindex, 0. * 0.5 * (
fYmaxv[ind] -
fYminv[ind]), errT * 1.e3);
704 double errA =
modelsY[ind]->GetAerror(
709 std::cout <<
"A = " << A <<
" error = " << errA <<
"\n";
724 grTy->GetXaxis()->SetLimits(-3., 3.);
725 grTy->GetYaxis()->SetLimits(0., 300.);
728 grTy->GetXaxis()->SetTitle(
"y");
729 grTy->GetYaxis()->SetTitle(
"T (MeV)");
752 std::cout <<
"eta = " << params.
eta.
value <<
" "
753 <<
"Error = " << params.
eta.
error <<
"\n";
760 * TMath::Sqrt((s2 - s1 * s1 / Numb) / (Numb - 1.) / Numb);
762 std::cout <<
"etamax = " << etamax <<
" "
763 <<
"Error = " << erretamax <<
"\n";
773 std::cout <<
"<" <<
name <<
"> = "
788 if (RecoLevel == -1) {
789 vector<CbmMCTrack> vRTracksMC;
791 std::cout <<
"MC tracks: " << nTracksMC <<
"\n";
792 vRTracksMC.resize(nTracksMC);
793 for (
int iTr = 0; iTr < nTracksMC; iTr++)
797 for (
int iTr = 0; iTr < nTracksMC; iTr++) {
798 if (vRTracksMC[iTr].GetPdgCode() ==
fPDGID
799 && vRTracksMC[iTr].GetMotherId() == -1) {
801 double pt = vRTracksMC[iTr].GetPt();
802 double ty = vRTracksMC[iTr].GetRapidity() -
ycm;
803 double mt = TMath::Sqrt(
804 vRTracksMC[iTr].GetMass() * vRTracksMC[iTr].GetMass() + pt * pt);
818 for (
unsigned int ind = 0; ind <
fYminv.size(); ++ind) {
834 for (
int itype = 2; itype <= 3; ++itype) {
836 const kfvector_int& pdgs = tr.PDG();
837 for (
unsigned int ind = 0; ind < pdgs.size(); ++ind) {
838 int iPDG = pdgs[ind];
842 double pt = tr.Pt(ind);
843 double p = tr.P(ind);
844 double m = TDatabasePDG::Instance()->GetParticle(iPDG)->Mass();
845 double p0 = TMath::Sqrt(
m *
m + p * p);
846 double pz = TMath::Sqrt(p * p - pt * pt);
847 double ty = 0.5 *
log((p0 + pz) / (p0 - pz)) -
ycm;
848 pz = TMath::Sqrt(
m *
m + pt * pt) * TMath::SinH(ty);
849 double mt = TMath::Sqrt(
m *
m + pt * pt);
863 for (
unsigned int ind2 = 0; ind2 <
fYminv.size(); ++ind2) {