CbmRoot
Loading...
Searching...
No Matches
CbmL1Performance.cxx
Go to the documentation of this file.
1/* Copyright (C) 2006-2020 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: Ivan Kisel, Sergey Gorbunov [committer], Igor Kulakov, Valentina Akishina, Grigory Kozlov */
4
5/*
6 *====================================================================
7 *
8 * CBM Level 1 Reconstruction
9 *
10 * Authors: I.Kisel, S.Gorbunov
11 *
12 * e-mail : ikisel@kip.uni-heidelberg.de
13 *
14 *====================================================================
15 *
16 * L1 Fit performance
17 *
18 *====================================================================
19 */
20
21#include "CaFramework.h"
22#include "CaToolsDebugger.h"
23#include "CbmL1.h"
24#include "CbmL1Constants.h"
25#include "CbmL1Counters.h"
26#include "CbmMatch.h"
27#include "CbmMuchPixelHit.h"
28#include "CbmMuchPoint.h"
29#include "CbmQaTable.h"
30#include "CbmStsDigi.h"
31#include "CbmStsHit.h"
32#include "CbmStsPoint.h"
33#include "CbmStsSetup.h"
34#include "CbmStsStation.h"
35#include "CbmTofHit.h"
36#include "CbmTofPoint.h"
37#include "CbmTrdHit.h"
38#include "CbmTrdPoint.h"
39#include "FairField.h"
40#include "FairRunAna.h"
41#include "FairTrackParam.h" // for vertex pulls
42#include "KfTrackKalmanFilter.h"
43#include "KfTrackParam.h"
44#include "TFile.h"
45#include "TH1.h"
46#include "TH2.h"
47#include "TMath.h"
48#include "TNtuple.h"
49#include "TProfile.h"
50
51#include <Logger.h>
52
53#include <boost/filesystem.hpp>
54
55#include <cmath>
56#include <list>
57#include <map>
58#include <sstream>
59#include <vector>
60
61using namespace cbm::algo::ca;
62using namespace cbm::algo;
63
64using cbm::algo::kf::TrackParamV;
65using cbm::ca::tools::Debugger;
66using std::map;
67using std::vector;
68
69struct TL1PerfEfficiencies : public TL1Efficiencies {
70
71 TL1PerfEfficiencies()
72 : TL1Efficiencies()
73 , ratio_killed()
74 , ratio_clone()
75 , ratio_length()
76 , ratio_fakes()
77 , killed()
78 , clone()
79 , reco_length()
80 , reco_fakes()
81 , mc_length()
82 , mc_length_hits()
83 {
84 // add total efficiency
85 AddCounter("long_fast_prim", "LongRPrim efficiency");
86 AddCounter("fast_prim", "FastPrim efficiency");
87 AddCounter("fast_sec", "FastSec efficiency");
88 AddCounter("fast", "Fastset efficiency");
89 AddCounter("total", "Allset efficiency");
90 AddCounter("slow_prim", "SlowPrim efficiency");
91 AddCounter("slow_sec", "SlowSec efficiency");
92 AddCounter("slow", "Slow efficiency");
93 AddCounter("d0", "D0 efficiency");
94 AddCounter("short", "Short123s efficiency");
95 AddCounter("shortPion", "Short123s pion eff");
96 AddCounter("shortProton", "Short123s proton eff");
97 AddCounter("shortKaon", "Short123s kaon eff");
98 AddCounter("shortE", "Short123s e eff");
99 AddCounter("shortRest", "Short123s rest eff");
100
101 AddCounter("fast_sec_e", "FastSec E efficiency");
102 AddCounter("fast_e", "Fastset E efficiency");
103 AddCounter("total_e", "Allset E efficiency");
104 AddCounter("slow_sec_e", "SlowSec E efficiency");
105 AddCounter("slow_e", "Slow E efficiency");
106 }
107
108 virtual ~TL1PerfEfficiencies(){};
109
110 virtual void AddCounter(const TString& shortname, const TString& name)
111 {
112 TL1Efficiencies::AddCounter(shortname, name);
113 ratio_killed.AddCounter();
114 ratio_clone.AddCounter();
115 ratio_length.AddCounter();
116 ratio_fakes.AddCounter();
117 killed.AddCounter();
118 clone.AddCounter();
119 reco_length.AddCounter();
120 reco_fakes.AddCounter();
121 mc_length.AddCounter();
122 mc_length_hits.AddCounter();
123 };
124
136
137 void CalcEff()
138 {
139 TL1Efficiencies::CalcEff();
140 ratio_killed = killed / mc;
141 ratio_clone = clone / mc;
142 TL1TracksCatCounters<int> allReco = reco + clone;
143 ratio_length = reco_length / allReco;
144 ratio_fakes = reco_fakes / allReco;
145 };
146
147 void Inc(bool isReco, bool isKilled, double _ratio_length, double _ratio_fakes, int _nclones, int _mc_length,
148 int _mc_length_hits, const TString& name)
149 {
150 TL1Efficiencies::Inc(isReco, name);
151
152 const int index = indices[name];
153
154 if (isKilled) killed.counters[index]++;
155 reco_length.counters[index] += _ratio_length;
156 reco_fakes.counters[index] += _ratio_fakes;
157 clone.counters[index] += _nclones;
158 mc_length.counters[index] += _mc_length;
159 mc_length_hits.counters[index] += _mc_length_hits;
160 };
161
162 void PrintEff(bool ifPrintTableToLog = false, TDirectory* outDir = nullptr,
163 const std::string& nameOfTable = "efficiency_table")
164 {
165 int NCounters = mc.GetNcounters();
166 std::vector<std::string> rowNames(NCounters + 2);
167 for (int iC = 0; iC < NCounters; ++iC) {
168 rowNames[iC] = std::string(names[iC].Data());
169 }
170 rowNames[NCounters] = "Ghost prob.";
171 rowNames[NCounters + 1] = "N ghosts";
172 std::vector<std::string> colNames = {"Eff.", "Killed", "Length", "Fakes", "Clones",
173 "All Reco", "All MC", "MCl(hits)", "MCl(MCps)"};
174
175 CbmQaTable* aTable = new CbmQaTable(nameOfTable.c_str(), "Track Efficiency", NCounters + 2, 9);
176 aTable->SetColWidth(20);
177 aTable->SetNamesOfRows(rowNames);
178 aTable->SetNamesOfCols(colNames);
179 for (int iC = 0; iC < NCounters; iC++) {
180 aTable->SetCell(iC, 0, ratio_reco.counters[iC]);
181 aTable->SetCell(iC, 1, ratio_killed.counters[iC]);
182 aTable->SetCell(iC, 2, ratio_length.counters[iC]);
183 aTable->SetCell(iC, 3, ratio_fakes.counters[iC]);
184 aTable->SetCell(iC, 4, ratio_clone.counters[iC]);
185 if (nEvents > 0) {
186 aTable->SetCell(iC, 5, reco.counters[iC] / double(nEvents));
187 aTable->SetCell(iC, 6, mc.counters[iC] / double(nEvents));
188 }
189 else {
190 aTable->SetCell(iC, 5, -1.);
191 aTable->SetCell(iC, 6, -1.);
192 }
193 if (mc.counters[iC] > 0) {
194 aTable->SetCell(iC, 7, mc_length_hits.counters[iC] / double(mc.counters[iC]));
195 aTable->SetCell(iC, 8, mc_length.counters[iC] / double(mc.counters[iC]));
196 }
197 else {
198 aTable->SetCell(iC, 7, -1.);
199 aTable->SetCell(iC, 8, -1.);
200 }
201 }
202 aTable->SetCell(NCounters, 0, ratio_ghosts);
203 aTable->SetCell(NCounters + 1, 0, ghosts);
204
205 if (ifPrintTableToLog) {
206 LOG(info) << *aTable; // print a table to log
207 }
208 if (outDir != nullptr) {
209 aTable->SetDirectory(outDir);
210 }
211 else {
212 aTable->SetDirectory(nullptr);
213 delete aTable;
214 }
215 };
216
217 TL1TracksCatCounters<double> ratio_killed;
218 TL1TracksCatCounters<double> ratio_clone;
219 TL1TracksCatCounters<double> ratio_length;
220 TL1TracksCatCounters<double> ratio_fakes;
221
222 TL1TracksCatCounters<int> killed;
223 TL1TracksCatCounters<int> clone;
224 TL1TracksCatCounters<double> reco_length;
225 TL1TracksCatCounters<double> reco_fakes;
226 TL1TracksCatCounters<int> mc_length;
227 TL1TracksCatCounters<int> mc_length_hits;
228};
229
230
231void CbmL1::EfficienciesPerformance(bool doFinish)
232{
233 static TL1PerfEfficiencies L1_NTRA; // average efficiencies
234
235 static int L1_NEVENTS = 0;
236 static double L1_CATIME = 0.0;
237
238 if (doFinish) {
239 L1_NTRA.CalcEff();
240 L1_NTRA.PrintEff(false, fTableDir);
241 return;
242 }
243
244 TL1PerfEfficiencies ntra; // efficiencies for current event
245
246 cbm::ca::tools::Debugger::Instance().AddNtuple("ghost", "it:ih:p:x:y:z:t:dx:dy");
247 static int statNghost = 0;
248
249 for (vector<CbmL1Track>::iterator rtraIt = fvRecoTracks.begin(); rtraIt != fvRecoTracks.end(); ++rtraIt) {
250 ntra.ghosts += rtraIt->IsGhost();
251 if (0 && rtraIt->IsGhost()) { // Debug.
252 TrackParamV tr(*rtraIt);
253 std::stringstream ss;
254 ss << " L1: ghost track: nhits " << rtraIt->GetNofHits() << " p " << 1. / rtraIt->GetQp() << " purity "
255 << rtraIt->GetMaxPurity() << " | hits ";
256 for (map<int, int>::iterator posIt = rtraIt->hitMap.begin(); posIt != rtraIt->hitMap.end(); posIt++) {
257 ss << " (" << posIt->second << " " << posIt->first << ") ";
258 }
259 LOG(info) << ss.str();
260 for (map<int, int>::iterator posIt = rtraIt->hitMap.begin(); posIt != rtraIt->hitMap.end(); posIt++) {
261 const auto& mcTrk = fMCData.GetTrack(posIt->first);
262 LOG(info) << "mc " << posIt->first << " pdg " << mcTrk.GetPdgCode() << " mother: " << mcTrk.GetMotherId()
263 << " chain " << mcTrk.GetChainId() << " n mc stations: " << mcTrk.GetNofConsStationsWithPoint();
264 }
265 for (unsigned int i = 0; i < rtraIt->Hits.size(); i++) {
266 const ca::Hit& h = fpAlgo->GetInputData().GetHit(rtraIt->Hits[i]);
267 const CbmL1HitDebugInfo& s = fvHitDebugInfo[rtraIt->Hits[i]];
268 LOG(info) << " x y z t " << s.x << " " << s.y << " " << h.Z() << " dx " << s.dx << " dy " << s.dy;
269 cbm::ca::tools::Debugger::Instance().FillNtuple("ghost", statNghost, i, fabs(1. / tr.GetQp()[0]), s.x, s.y,
270 h.Z(), h.T(), s.dx, s.dy);
271 }
272 LOG(info) << tr.ToString(0);
273 statNghost++;
274 }
275 }
276
277 int sta_nhits[fpAlgo->GetParameters().GetNstationsActive()];
278 int sta_nfakes[fpAlgo->GetParameters().GetNstationsActive()];
279
280 for (int i = 0; i < fpAlgo->GetParameters().GetNstationsActive(); i++) {
281 sta_nhits[i] = 0;
282 sta_nfakes[i] = 0;
283 }
284
285 for (const auto& mcTrk : fMCData.GetTrackContainer()) {
286 // if( !( mcTrk.GetPdgCode() == -11 && mcTrk.GetMotherId() == -1 ) ) continue; // electrons only
287 if (!mcTrk.IsReconstructable() && !mcTrk.IsAdditional()) continue;
288
289 // -- find used constans --
290 // is track reconstructed
291 const bool reco = (mcTrk.IsReconstructed());
292 // is track killed. At least one hit of it belong to some recoTrack
293 const bool killed = !reco && mcTrk.IsDisturbed();
294 // ration length for current mcTrack
295 double ratio_length = 0;
296 double ratio_fakes = 0;
297 double mc_length_hits = mcTrk.GetTotNofStationsWithHit();
298
299
300 int mc_length = mcTrk.GetTotNofStationsWithPoint();
301 if (reco) {
302 for (unsigned int irt : mcTrk.GetRecoTrackIndexes()) {
303 const auto& rTrk = fvRecoTracks[irt];
304 ratio_length += static_cast<double>(rTrk.GetNofHits()) * rTrk.GetMaxPurity() / mc_length_hits;
305 ratio_fakes += 1 - rTrk.GetMaxPurity();
306 }
307 }
308
309
310 // number of clones
311 int nclones = 0;
312 if (reco) nclones = mcTrk.GetNofClones();
313 // if (nclones){ // Debug. Look at clones
314 // for (int irt = 0; irt < rTracks.size(); irt++){
315 // const int ista = fvHitDebugInfo[rTracks[irt]->Hits[0]].iStation;
316 // LOG(info) << rTracks[irt]->GetNOfHits() << "(" << ista << ") ";
317 // }
318 // LOG(info) << mtra.NStations();
319 // }
320
321 if (mcTrk.IsAdditional()) { // short
322 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "short");
323 switch (mcTrk.GetPdgCode()) {
324 case 11:
325 case -11:
326 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "shortE");
327 break;
328 case 211:
329 case -211:
330 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "shortPion");
331 break;
332 case 321:
333 case -321:
334 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "shortKaon");
335 break;
336 case 2212:
337 case -2212:
338 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "shortProton");
339 break;
340 default: ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "shortRest");
341 }
342 }
343 else { // separate all efficiecies from short eff
344
345
346 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "total");
347
348 if (mcTrk.IsSignal()) { // D0
349 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "d0");
350 }
351
352 if (mcTrk.GetP() > CbmL1Constants::MinFastMom) { // fast tracks
353 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "fast");
354
355 if (mcTrk.IsPrimary()) { // fast primary
356 if (mcTrk.GetTotNofStationsWithHit() == fNStations) { // long fast primary
357 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "long_fast_prim");
358 }
359 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "fast_prim");
360 }
361 else { // fast secondary
362 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "fast_sec");
363 }
364 }
365 else { // slow set of tracks
366 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "slow");
367
368 if (mcTrk.IsPrimary()) { // slow primary
369 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "slow_prim");
370 }
371 else {
372 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "slow_sec");
373 }
374 } // if slow
375
376 if (mcTrk.GetPdgCode() == 11 || mcTrk.GetPdgCode() == -11) {
377 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "total_e");
378
379 if (mcTrk.GetP() > CbmL1Constants::MinFastMom) { // fast tracks
380 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "fast_e");
381
382 if (mcTrk.IsPrimary()) { // fast primary
383 }
384 else { // fast secondary
385 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "fast_sec_e");
386 }
387 }
388 else { // slow set of tracks
389 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "slow_e");
390
391 if (mcTrk.IsPrimary()) { // slow primary
392 }
393 else {
394 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "slow_sec_e");
395 }
396 } // if slow
397 }
398 }
399
400 } // for mcTracks
401
402 L1_CATIME += fTrackingTime;
403 L1_NEVENTS++;
404 ntra.IncNEvents();
405 L1_NTRA += ntra;
406
407 ntra.CalcEff();
408 L1_NTRA.CalcEff();
409
410 if (fVerbose) {
411 if (fVerbose > 1) {
412 ntra.PrintEff(true);
413 std::stringstream ss;
414 ss << "Number of true and fake hits in stations: \n";
415 for (int i = 0; i < fpAlgo->GetParameters().GetNstationsActive(); i++) {
416 ss << sta_nhits[i] - sta_nfakes[i] << "+" << sta_nfakes[i] << " ";
417 }
418 LOG(info) << ss.str();
419 } // fVerbose > 1
420 LOG(info) << "\n"
421 << "L1 ACCUMULATED STAT : " << L1_NEVENTS << " EVENTS \n";
422 L1_NTRA.PrintEff(true);
423 LOG(info) << "Reconstructible MC tracks/event: "
424 << (double(L1_NTRA.mc.counters[L1_NTRA.indices["total"]]) / double(L1_NEVENTS));
425 LOG(info) << "Reconstructed MC tracks/event: "
426 << (double(L1_NTRA.reco.counters[L1_NTRA.indices["total"]]) / double(L1_NEVENTS));
427 LOG(info) << "\nCA Track Finder: " << L1_CATIME / L1_NEVENTS << " s/time slice \n";
428 }
429} // void CbmL1::Performance()
430
431
432void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitFast on match data in CbmL1**Track classes
433{
434
435 static TProfile *p_eff_all_vs_mom, *p_eff_prim_vs_mom, *p_eff_sec_vs_mom, *p_eff_d0_vs_mom, *p_eff_prim_vs_theta,
436 *p_eff_all_vs_pt, *p_eff_prim_vs_pt, *p_eff_all_vs_nhits, *p_eff_prim_vs_nhits, *p_eff_sec_vs_nhits;
437
438 static TH1F *h_reg_MCmom, *h_acc_MCmom, *h_reco_MCmom, *h_ghost_Rmom;
439 static TH1F *h_reg_prim_MCmom, *h_acc_prim_MCmom, *h_reco_prim_MCmom;
440 static TH1F *h_reg_sec_MCmom, *h_acc_sec_MCmom, *h_reco_sec_MCmom;
441
442 static TH1F* h_acc_mom_short123s;
443
444 static TH1F *h_reg_mom_prim, *h_reg_mom_sec, *h_reg_nhits_prim, *h_reg_nhits_sec;
445 static TH1F *h_acc_mom_prim, *h_acc_mom_sec, *h_acc_nhits_prim, *h_acc_nhits_sec;
446 static TH1F *h_reco_mom_prim, *h_reco_mom_sec, *h_reco_nhits_prim, *h_reco_nhits_sec;
447 static TH1F *h_rest_mom_prim, *h_rest_mom_sec, *h_rest_nhits_prim, *h_rest_nhits_sec;
448
449 //static TH1F *h_hit_density[10];
450
451 static TH1F *h_ghost_purity, *h_ghost_mom, *h_ghost_nhits, *h_ghost_fstation, *h_ghost_chi2, *h_ghost_prob,
452 *h_ghost_tx, *h_ghost_ty;
453 static TH1F *h_reco_mom, *h_reco_d0_mom, *h_reco_nhits, *h_reco_station, *h_reco_chi2, *h_reco_prob, *h_rest_prob,
454 *h_reco_purity, *h_reco_time;
455 static TProfile *h_reco_timeNtr, *h_reco_timeNhit;
456 static TProfile *h_reco_fakeNtr, *h_reco_fakeNhit;
457 static TH1F *h_tx, *h_ty, *h_sec_r, *h_ghost_r;
458
459 static TH1F *h_notfound_mom, *h_notfound_nhits, *h_notfound_station, *h_notfound_r, *h_notfound_tx, *h_notfound_ty;
460
461 static TH1F *h_mcp, *h_nmchits;
462 // static TH1F *h_chi2, *h_prob, *MC_vx, *MC_vy, *MC_vz, *VtxChiPrim, *VtxChiSec;
463
464 // static TH2F *P_vs_P ;
465
466 static TH2F *h2_vertex, *h2_prim_vertex, *h2_sec_vertex;
467 //static TH3F *h3_vertex, *h3_prim_vertex, *h3_sec_vertex;
468
469 static TH2F *h2_reg_nhits_vs_mom_prim, *h2_reg_nhits_vs_mom_sec, *h2_reg_fstation_vs_mom_prim,
470 *h2_reg_fstation_vs_mom_sec, *h2_reg_lstation_vs_fstation_prim, *h2_reg_lstation_vs_fstation_sec;
471 static TH2F *h2_acc_nhits_vs_mom_prim, *h2_acc_nhits_vs_mom_sec, *h2_acc_fstation_vs_mom_prim,
472 *h2_acc_fstation_vs_mom_sec, *h2_acc_lstation_vs_fstation_prim, *h2_acc_lstation_vs_fstation_sec;
473 static TH2F *h2_reco_nhits_vs_mom_prim, *h2_reco_nhits_vs_mom_sec, *h2_reco_fstation_vs_mom_prim,
474 *h2_reco_fstation_vs_mom_sec, *h2_reco_lstation_vs_fstation_prim, *h2_reco_lstation_vs_fstation_sec;
475 static TH2F *h2_ghost_nhits_vs_mom, *h2_ghost_fstation_vs_mom, *h2_ghost_lstation_vs_fstation;
476 static TH2F *h2_rest_nhits_vs_mom_prim, *h2_rest_nhits_vs_mom_sec, *h2_rest_fstation_vs_mom_prim,
477 *h2_rest_fstation_vs_mom_sec, *h2_rest_lstation_vs_fstation_prim, *h2_rest_lstation_vs_fstation_sec;
478
479 static TH1F* h_reg_phi_prim;
480 static TH1F* h_reg_phi_sec;
481 static TH1F* h_acc_phi_prim;
482 static TH1F* h_acc_phi_sec;
483 static TH1F* h_reco_phi_prim;
484 static TH1F* h_reco_phi_sec;
485 static TH1F* h_rest_phi_prim;
486 static TH1F* h_rest_phi_sec;
487 static TH1F* h_ghost_phi;
488 static TH1F* h_reco_phi;
489 static TH1F* h_notfound_phi;
490
491 static TH2F* h2_fst_hit_yz; // occupancy of track first hit in y-z plane
492 static TH2F* h2_lst_hit_yz; // occupancy of track last hit in y-z plane
493 static TH2F* h2_all_hit_yz; // occupancy of track all hits in y-z plane
494
495 static bool first_call = 1;
496
497 if (first_call) {
498 first_call = 0;
499
500 TDirectory* curdir = gDirectory;
501 gDirectory = fHistoDir;
502
503 h2_fst_hit_yz =
504 new TH2F("h2_fst_hit_yz", "Track first hit occupancy in y-z plane;z [cm];y [cm]", 80, 0, 0, 80, 0, 0);
505 h2_lst_hit_yz =
506 new TH2F("h2_lst_hit_yz", "Track last hit occupancy in y-z plane;z[ cm];y [cm]", 80, 0, 0, 80, 0, 0);
507 h2_all_hit_yz = new TH2F("h2_all_hit_yz", "Track hit occupancy in y-z plane;z[ cm];y [cm]", 80, 0, 0, 80, 0, 0);
508
509 p_eff_all_vs_mom = new TProfile("p_eff_all_vs_mom", "AllSet Efficiency vs Momentum", 100, 0.0, 5.0, 0.0, 100.0);
510 p_eff_prim_vs_mom =
511 new TProfile("p_eff_prim_vs_mom", "Primary Set Efficiency vs Momentum", 100, 0.0, 5.0, 0.0, 100.0);
512 p_eff_sec_vs_mom =
513 new TProfile("p_eff_sec_vs_mom", "Secondary Set Efficiency vs Momentum", 100, 0.0, 5.0, 0.0, 100.0);
514 p_eff_d0_vs_mom =
515 new TProfile("p_eff_d0_vs_mom", "D0 Secondary Tracks Efficiency vs Momentum", 150, 0.0, 15.0, 0.0, 100.0);
516 p_eff_prim_vs_theta =
517 new TProfile("p_eff_prim_vs_theta", "All Primary Set Efficiency vs Theta", 100, 0.0, 30.0, 0.0, 100.0);
518 p_eff_all_vs_pt = new TProfile("p_eff_all_vs_pt", "AllSet Efficiency vs Pt", 100, 0.0, 5.0, 0.0, 100.0);
519 p_eff_prim_vs_pt = new TProfile("p_eff_prim_vs_pt", "Primary Set Efficiency vs Pt", 100, 0.0, 5.0, 0.0, 100.0);
520
521 p_eff_all_vs_nhits = new TProfile("p_eff_all_vs_nhits", "AllSet Efficiency vs NMCHits", 8, 3.0, 11.0, 0.0, 100.0);
522 p_eff_prim_vs_nhits =
523 new TProfile("p_eff_prim_vs_nhits", "PrimSet Efficiency vs NMCHits", 8, 3.0, 11.0, 0.0, 100.0);
524 p_eff_sec_vs_nhits = new TProfile("p_eff_sec_vs_nhits", "SecSet Efficiency vs NMCHits", 8, 3.0, 11.0, 0.0, 100.0);
525
526 h_reg_MCmom = new TH1F("h_reg_MCmom", "Momentum of registered tracks", 100, 0.0, 5.0);
527 h_acc_MCmom = new TH1F("h_acc_MCmom", "Reconstructable tracks", 100, 0.0, 5.0);
528 h_reco_MCmom = new TH1F("h_reco_MCmom", "Reconstructed tracks", 100, 0.0, 5.0);
529
530 h_ghost_Rmom = new TH1F("h_ghost_Rmom", "Ghost tracks", 100, 0.0, 5.0);
531 h_reg_prim_MCmom = new TH1F("h_reg_prim_MCmom", "Momentum of registered tracks", 100, 0.0, 5.0);
532 h_acc_prim_MCmom = new TH1F("h_acc_prim_MCmom", "Reconstructable tracks", 100, 0.0, 5.0);
533 h_reco_prim_MCmom = new TH1F("h_reco_prim_MCmom", "Reconstructed tracks", 100, 0.0, 5.0);
534 h_reg_sec_MCmom = new TH1F("h_reg_sec_MCmom", "Momentum of registered tracks", 100, 0.0, 5.0);
535 h_acc_sec_MCmom = new TH1F("h_acc_sec_MCmom", "Reconstructable tracks", 100, 0.0, 5.0);
536 h_reco_sec_MCmom = new TH1F("h_reco_sec_MCmom", "Reconstructed tracks", 100, 0.0, 5.0);
537
538 h_acc_mom_short123s =
539 new TH1F("h_acc_mom_short123s", "Momentum of accepted tracks with 3 hits on first stations", 500, 0.0, 5.0);
540
541 h_reg_mom_prim = new TH1F("h_reg_mom_prim", "Momentum of registered primary tracks", 500, 0.0, 5.0);
542 h_reg_mom_sec = new TH1F("h_reg_mom_sec", "Momentum of registered secondary tracks", 500, 0.0, 5.0);
543 h_acc_mom_prim = new TH1F("h_acc_mom_prim", "Momentum of accepted primary tracks", 500, 0.0, 5.0);
544 h_acc_mom_sec = new TH1F("h_acc_mom_sec", "Momentum of accepted secondary tracks", 500, 0.0, 5.0);
545 h_reco_mom_prim = new TH1F("h_reco_mom_prim", "Momentum of reconstructed primary tracks", 500, 0.0, 5.0);
546 h_reco_mom_sec = new TH1F("h_reco_mom_sec", "Momentum of reconstructed secondary tracks", 500, 0.0, 5.0);
547 h_rest_mom_prim = new TH1F("h_rest_mom_prim", "Momentum of not found primary tracks", 500, 0.0, 5.0);
548 h_rest_mom_sec = new TH1F("h_rest_mom_sec", "Momentum of not found secondary tracks", 500, 0.0, 5.0);
549
550 h_reg_phi_prim = new TH1F("h_reg_phi_prim", "Azimuthal angle of registered primary tracks", 1000, -3.15, 3.15);
551 h_reg_phi_sec = new TH1F("h_reg_phi_sec", "Azimuthal angle of registered secondary tracks", 1000, -3.15, 3.15);
552 h_acc_phi_prim = new TH1F("h_acc_phi_prim", "Azimuthal angle of accepted primary tracks", 1000, -3.15, 3.15);
553 h_acc_phi_sec = new TH1F("h_acc_phi_sec", "Azimuthal angle of accepted secondary tracks", 1000, -3.15, 3.15);
554 h_reco_phi_prim = new TH1F("h_reco_phi_prim", "Azimuthal angle of reconstructed primary tracks", 1000, -3.15, 3.15);
555 h_reco_phi_sec = new TH1F("h_reco_phi_sec", "Azimuthal angle of reconstructed secondary tracks", 1000, -3.15, 3.15);
556 h_rest_phi_prim = new TH1F("h_rest_phi_prim", "Azimuthal angle of not found primary tracks", 1000, -3.15, 3.15);
557 h_rest_phi_sec = new TH1F("h_rest_phi_sec", "Azimuthal angle of not found secondary tracks", 1000, -3.15, 3.15);
558
559 h_reg_nhits_prim = new TH1F("h_reg_nhits_prim", "Number of hits in registered primary track", 51, -0.1, 10.1);
560 h_reg_nhits_sec = new TH1F("h_reg_nhits_sec", "Number of hits in registered secondary track", 51, -0.1, 10.1);
561 h_acc_nhits_prim = new TH1F("h_acc_nhits_prim", "Number of hits in accepted primary track", 51, -0.1, 10.1);
562 h_acc_nhits_sec = new TH1F("h_acc_nhits_sec", "Number of hits in accepted secondary track", 51, -0.1, 10.1);
563 h_reco_nhits_prim = new TH1F("h_reco_nhits_prim", "Number of hits in reconstructed primary track", 51, -0.1, 10.1);
564 h_reco_nhits_sec = new TH1F("h_reco_nhits_sec", "Number of hits in reconstructed secondary track", 51, -0.1, 10.1);
565 h_rest_nhits_prim = new TH1F("h_rest_nhits_prim", "Number of hits in not found primary track", 51, -0.1, 10.1);
566 h_rest_nhits_sec = new TH1F("h_rest_nhits_sec", "Number of hits in not found secondary track", 51, -0.1, 10.1);
567
568 h_ghost_mom = new TH1F("h_ghost_mom", "Momentum of ghost tracks", 500, 0.0, 5.0);
569 h_ghost_phi = new TH1F("h_ghost_phi", "Azimuthal angle of ghost tracks", 1000, -3.15, 3.15);
570 h_ghost_nhits = new TH1F("h_ghost_nhits", "Number of hits in ghost track", 51, -0.1, 10.1);
571 h_ghost_fstation = new TH1F("h_ghost_fstation", "First station of ghost track", 50, -0.5, 10.0);
572 h_ghost_chi2 = new TH1F("h_ghost_chi2", "Chi2/NDF of ghost track", 50, -0.5, 10.0);
573 h_ghost_prob = new TH1F("h_ghost_prob", "Prob of ghost track", 505, 0., 1.01);
574 h_ghost_r = new TH1F("h_ghost_r", "R of ghost track at the first hit", 50, 0.0, 150.0);
575 h_ghost_tx = new TH1F("h_ghost_tx", "tx of ghost track at the first hit", 50, -5.0, 5.0);
576 h_ghost_ty = new TH1F("h_ghost_ty", "ty of ghost track at the first hit", 50, -1.0, 1.0);
577 h_ghost_purity = new TH1F("h_ghost_purity", "Ghost: percentage of correct hits", 100, -0.5, 100.5);
578
579 h_reco_mom = new TH1F("h_reco_mom", "Momentum of reco track", 50, 0.0, 5.0);
580 h_reco_phi = new TH1F("h_reco_phi", "Azimuthal angle of reco track", 1000, -3.15, 3.15);
581 h_reco_nhits = new TH1F("h_reco_nhits", "Number of hits in reco track", 50, 0.0, 10.0);
582 h_reco_station = new TH1F("h_reco_station", "First station of reco track", 50, -0.5, 10.0);
583 h_reco_chi2 = new TH1F("h_reco_chi2", "Chi2/NDF of reco track", 50, -0.5, 10.0);
584 h_reco_prob = new TH1F("h_reco_prob", "Prob of reco track", 505, 0., 1.01);
585 h_rest_prob = new TH1F("h_rest_prob", "Prob of reco rest track", 505, 0., 1.01);
586 h_reco_purity = new TH1F("h_reco_purity", "Percentage of correct hits", 100, -0.5, 100.5);
587 h_reco_time = new TH1F("h_reco_time", "CA Track Finder Time (s/ev)", 20, 0.0, 20.0);
588 h_reco_timeNtr = new TProfile("h_reco_timeNtr", "CA Track Finder Time (s/ev) vs N Tracks", 200, 0.0, 1000.0);
589 h_reco_timeNhit = new TProfile("h_reco_timeNhit", "CA Track Finder Time (s/ev) vs N Hits", 200, 0.0, 30000.0);
590 h_reco_fakeNtr = new TProfile("h_reco_fakeNtr", "N Fake hits vs N Tracks", 200, 0.0, 1000.0);
591 h_reco_fakeNhit = new TProfile("h_reco_fakeNhit", "N Fake hits vs N Hits", 200, 0.0, 30000.0);
592
593 h_reco_d0_mom = new TH1F("h_reco_d0_mom", "Momentum of reco D0 track", 150, 0.0, 15.0);
594
595 // h_hit_density[0] = new TH1F("h_hit_density0", "Hit density station 1", 50, 0.0, 5.00);
596 // h_hit_density[1] = new TH1F("h_hit_density1", "Hit density station 2", 100, 0.0, 10.00);
597 // h_hit_density[2] = new TH1F("h_hit_density2", "Hit density station 3", 140, 0.0, 13.99);
598 // h_hit_density[3] = new TH1F("h_hit_density3", "Hit density station 4", 180, 0.0, 18.65);
599 // h_hit_density[4] = new TH1F("h_hit_density4", "Hit density station 5", 230, 0.0, 23.32);
600 // h_hit_density[5] = new TH1F("h_hit_density5", "Hit density station 6", 270, 0.0, 27.98);
601 // h_hit_density[6] = new TH1F("h_hit_density6", "Hit density station 7", 340, 0.0, 34.97);
602 // h_hit_density[7] = new TH1F("h_hit_density7", "Hit density station 8", 460, 0.0, 46.63);
603 // h_hit_density[8] = new TH1F("h_hit_density8", "Hit density station 9", 500, 0.0, 50.00);
604 // h_hit_density[9] = new TH1F("h_hit_density8", "Hit density station 9", 500, 0.0, 50.00);
605 // h_hit_density[10] = new TH1F("h_hit_density8", "Hit density station 9", 500, 0.0, 50.00);
606
607 h_tx = new TH1F("h_tx", "tx of MC track at the vertex", 50, -0.5, 0.5);
608 h_ty = new TH1F("h_ty", "ty of MC track at the vertex", 50, -0.5, 0.5);
609
610 h_sec_r = new TH1F("h_sec_r", "R of sec MC track at the first hit", 50, 0.0, 15.0);
611
612 h_notfound_mom = new TH1F("h_notfound_mom", "Momentum of not found track", 50, 0.0, 5.0);
613 h_notfound_phi = new TH1F("h_notfound_phi", "Momentum of not found track", 50, 0.0, 5.0);
614 h_notfound_nhits = new TH1F("h_notfound_nhits", "Num hits of not found track", 50, 0.0, 10.0);
615 h_notfound_station = new TH1F("h_notfound_station", "First station of not found track", 50, 0.0, 10.0);
616 h_notfound_r = new TH1F("h_notfound_r", "R at first hit of not found track", 50, 0.0, 15.0);
617 h_notfound_tx = new TH1F("h_notfound_tx", "tx of not found track at the first hit", 50, -5.0, 5.0);
618 h_notfound_ty = new TH1F("h_notfound_ty", "ty of not found track at the first hit", 50, -1.0, 1.0);
619
620 /*
621 h_chi2 = new TH1F("chi2", "Chi^2", 100, 0.0, 10.);
622 h_prob = new TH1F("prob", "Prob", 100, 0.0, 1.01);
623 VtxChiPrim = new TH1F("vtxChiPrim", "Chi to primary vtx for primary tracks", 100, 0.0, 10.);
624 VtxChiSec = new TH1F("vtxChiSec", "Chi to primary vtx for secondary tracks", 100, 0.0, 10.);
625*/
626 h_mcp = new TH1F("h_mcp", "MC P ", 500, 0.0, 5.0);
627 /*
628 MC_vx = new TH1F("MC_vx", "MC Vertex X", 100, -.05, .05);
629 MC_vy = new TH1F("MC_vy", "MC Vertex Y", 100, -.05, .05);
630 MC_vz = new TH1F("MC_vz", "MC Vertex Z", 100, -.05, .05);
631*/
632 h_nmchits = new TH1F("h_nmchits", "N STS hits on MC track", 50, 0.0, 10.0);
633
634 // P_vs_P = new TH2F("P_vs_P", "Resolution P/Q vs P", 20, 0., 20.,100, -.05, .05);
635
636 h2_vertex = new TH2F("h2_vertex", "2D vertex distribution", 105, -5., 100., 100, -50., 50.);
637 h2_prim_vertex = new TH2F("h2_primvertex", "2D primary vertex distribution", 105, -5., 100., 100, -50., 50.);
638 h2_sec_vertex = new TH2F("h2_sec_vertex", "2D secondary vertex distribution", 105, -5., 100., 100, -50., 50.);
639
640 //h3_vertex = new TH3F("h3_vertex", "3D vertex distribution", 20, -5., 100., 100, -50., 50., 100, -50., 50.);
641 //h3_prim_vertex = new TH3F("h3_primvertex", "3D vertex distribution", 20, -5., 100., 100, -50., 50., 100, -50., 50.);
642 //h3_sec_vertex = new TH3F("h3_sec_vertex", "3D vertex distribution", 20, -5., 100., 100, -50., 50., 100, -50., 50.);
643
644 h2_reg_nhits_vs_mom_prim =
645 new TH2F("h2_reg_nhits_vs_mom_prim", "NHits vs. Momentum (Reg. Primary Tracks)", 51, -0.05, 5.05, 11, -0.5, 10.5);
646 h2_reg_nhits_vs_mom_sec = new TH2F("h2_reg_nhits_vs_mom_sec", "NHits vs. Momentum (Reg. Secondary Tracks)", 51,
647 -0.05, 5.05, 11, -0.5, 10.5);
648 h2_acc_nhits_vs_mom_prim =
649 new TH2F("h2_acc_nhits_vs_mom_prim", "NHits vs. Momentum (Acc. Primary Tracks)", 51, -0.05, 5.05, 11, -0.5, 10.5);
650 h2_acc_nhits_vs_mom_sec = new TH2F("h2_acc_nhits_vs_mom_sec", "NHits vs. Momentum (Acc. Secondary Tracks)", 51,
651 -0.05, 5.05, 11, -0.5, 10.5);
652 h2_reco_nhits_vs_mom_prim = new TH2F("h2_reco_nhits_vs_mom_prim", "NHits vs. Momentum (Reco Primary Tracks)", 51,
653 -0.05, 5.05, 11, -0.5, 10.5);
654 h2_reco_nhits_vs_mom_sec = new TH2F("h2_reco_nhits_vs_mom_sec", "NHits vs. Momentum (Reco Secondary Tracks)", 51,
655 -0.05, 5.05, 11, -0.5, 10.5);
656 h2_ghost_nhits_vs_mom =
657 new TH2F("h2_ghost_nhits_vs_mom", "NHits vs. Momentum (Ghost Tracks)", 51, -0.05, 5.05, 11, -0.5, 10.5);
658 h2_rest_nhits_vs_mom_prim = new TH2F("h2_rest_nhits_vs_mom_prim", "NHits vs. Momentum (!Found Primary Tracks)", 51,
659 -0.05, 5.05, 11, -0.5, 10.5);
660 h2_rest_nhits_vs_mom_sec = new TH2F("h2_rest_nhits_vs_mom_sec", "NHits vs. Momentum (!Found Secondary Tracks)", 51,
661 -0.05, 5.05, 11, -0.5, 10.5);
662
663 h2_reg_fstation_vs_mom_prim =
664 new TH2F("h2_reg_fstation_vs_mom_prim", "First Station vs. Momentum (Reg. Primary Tracks)", 51, -0.05, 5.05, 11,
665 -0.5, 10.5);
666 h2_reg_fstation_vs_mom_sec =
667 new TH2F("h2_reg_fstation_vs_mom_sec", "First Station vs. Momentum (Reg. Secondary Tracks)", 51, -0.05, 5.05, 11,
668 -0.5, 10.5);
669 h2_acc_fstation_vs_mom_prim =
670 new TH2F("h2_acc_fstation_vs_mom_prim", "First Station vs. Momentum (Acc. Primary Tracks)", 51, -0.05, 5.05, 11,
671 -0.5, 10.5);
672 h2_acc_fstation_vs_mom_sec =
673 new TH2F("h2_acc_fstation_vs_mom_sec", "First Station vs. Momentum (Acc. Secondary Tracks)", 51, -0.05, 5.05, 11,
674 -0.5, 10.5);
675 h2_reco_fstation_vs_mom_prim =
676 new TH2F("h2_reco_fstation_vs_mom_prim", "First Station vs. Momentum (Reco Primary Tracks)", 51, -0.05, 5.05, 11,
677 -0.5, 10.5);
678 h2_reco_fstation_vs_mom_sec =
679 new TH2F("h2_reco_fstation_vs_mom_sec", "First Station vs. Momentum (Reco Secondary Tracks)", 51, -0.05, 5.05, 11,
680 -0.5, 10.5);
681 h2_ghost_fstation_vs_mom = new TH2F("h2_ghost_fstation_vs_mom", "First Station vs. Momentum (Ghost Tracks)", 51,
682 -0.05, 5.05, 11, -0.5, 10.5);
683 h2_rest_fstation_vs_mom_prim =
684 new TH2F("h2_rest_fstation_vs_mom_prim", "First Station vs. Momentum (!Found Primary Tracks)", 51, -0.05, 5.05,
685 11, -0.5, 10.5);
686 h2_rest_fstation_vs_mom_sec =
687 new TH2F("h2_rest_fstation_vs_mom_sec", "First Station vs. Momentum (!Found Secondary Tracks)", 51, -0.05, 5.05,
688 11, -0.5, 10.5);
689
690 h2_reg_lstation_vs_fstation_prim =
691 new TH2F("h2_reg_lstation_vs_fstation_prim", "Last vs. First Station (Reg. Primary Tracks)", 11, -0.5, 10.5, 11,
692 -0.5, 10.5);
693 h2_reg_lstation_vs_fstation_sec =
694 new TH2F("h2_reg_lstation_vs_fstation_sec", "Last vs. First Station (Reg. Secondary Tracks)", 11, -0.5, 10.5, 11,
695 -0.5, 10.5);
696 h2_acc_lstation_vs_fstation_prim =
697 new TH2F("h2_acc_lstation_vs_fstation_prim", "Last vs. First Station (Acc. Primary Tracks)", 11, -0.5, 10.5, 11,
698 -0.5, 10.5);
699 h2_acc_lstation_vs_fstation_sec =
700 new TH2F("h2_acc_lstation_vs_fstation_sec", "Last vs. First Station (Acc. Secondary Tracks)", 11, -0.5, 10.5, 11,
701 -0.5, 10.5);
702 h2_reco_lstation_vs_fstation_prim =
703 new TH2F("h2_reco_lstation_vs_fstation_prim", "Last vs. First Station (Reco Primary Tracks)", 11, -0.5, 10.5, 11,
704 -0.5, 10.5);
705 h2_reco_lstation_vs_fstation_sec =
706 new TH2F("h2_reco_lstation_vs_fstation_sec", "Last vs. First Station (Reco Secondary Tracks)", 11, -0.5, 10.5, 11,
707 -0.5, 10.5);
708 h2_ghost_lstation_vs_fstation = new TH2F("h2_ghost_lstation_vs_fstation", "Last vs. First Station (Ghost Tracks)",
709 11, -0.5, 10.5, 11, -0.5, 10.5);
710 h2_rest_lstation_vs_fstation_prim =
711 new TH2F("h2_rest_lstation_vs_fstation_prim", "Last vs. First Station (!Found Primary Tracks)", 11, -0.5, 10.5,
712 11, -0.5, 10.5);
713 h2_rest_lstation_vs_fstation_sec =
714 new TH2F("h2_rest_lstation_vs_fstation_sec", "Last vs. First Station (!Found Secondary Tracks)", 11, -0.5, 10.5,
715 11, -0.5, 10.5);
716
717 //maindir->cd();
718
719 // ----- Create list of all histogram pointers
720
721 gDirectory = curdir;
722
723 } // first_call
724
725
726 static int NEvents = 0;
727 if (NEvents > 0) {
728 h_reg_MCmom->Scale(NEvents);
729 h_acc_MCmom->Scale(NEvents);
730 h_reco_MCmom->Scale(NEvents);
731 h_ghost_Rmom->Scale(NEvents);
732 h_reg_prim_MCmom->Scale(NEvents);
733 h_acc_prim_MCmom->Scale(NEvents);
734 h_reco_prim_MCmom->Scale(NEvents);
735 h_reg_sec_MCmom->Scale(NEvents);
736 h_acc_sec_MCmom->Scale(NEvents);
737 h_reco_sec_MCmom->Scale(NEvents);
738 }
739 NEvents++;
740
741 // //hit density calculation: h_hit_density[10]
742 //
743 // for (vector<CbmL1HitDebugInfo>::iterator hIt = fvHitDebugInfo.begin(); hIt != fvHitDebugInfo.end(); ++hIt){
744 // float x = hIt->x;
745 // float y = hIt->y;
746 // float r = sqrt(x*x+y*y);
747 // h_hit_density[hIt->iStation]->Fill(r, 1.0/(2.0*3.1415*r));
748 // }
749
750 //
751 for (vector<CbmL1Track>::iterator rtraIt = fvRecoTracks.begin(); rtraIt != fvRecoTracks.end(); ++rtraIt) {
752 CbmL1Track* prtra = &(*rtraIt);
753 if ((prtra->Hits).size() < 1) continue;
754 { // fill histos
755 if (fabs(prtra->GetQp()) > 1.e-10)
756 h_reco_mom->Fill(
757 fabs(1.0 / prtra->GetQp())); // TODO: Is it a right precision? In FairTrackParam it is 1.e-4 (S.Zharko)
758 // NOTE: p = (TMath::Abs(fQp) > 1.e-4) ? 1. / TMath::Abs(fQp) : 1.e4; // FairTrackParam::Momentum(TVector3)
759 // h_reco_mom->Fill(TMath::Abs(prtra->T[4] > 1.e-4) ? 1. / TMath::Abs(prtra->T[4]) : 1.e+4); // this should be correct
760
761 h_reco_phi->Fill(prtra->GetPhi()); // TODO: What is precision?
762 h_reco_nhits->Fill((prtra->Hits).size());
763 CbmL1HitDebugInfo& mh = fvHitDebugInfo[prtra->Hits[0]];
764 h_reco_station->Fill(mh.iStation);
765
766 int iFstHit = prtra->GetFirstHitIndex();
767 auto& fstHit = fvHitDebugInfo[iFstHit];
768 h2_fst_hit_yz->Fill(fpAlgo->GetParameters().GetStation(fstHit.iStation).fZ[0], fstHit.GetY());
769
770 int iLstHit = prtra->GetLastHitIndex();
771 auto& lstHit = fvHitDebugInfo[iLstHit];
772 h2_lst_hit_yz->Fill(fpAlgo->GetParameters().GetStation(lstHit.iStation).fZ[0], lstHit.GetY());
773
774 for (int iH : prtra->Hits) {
775 const auto& hit = fvHitDebugInfo[iH];
776 int y = hit.GetY();
777 int z = fpAlgo->GetParameters().GetStation(hit.iStation).fZ[0];
778 h2_all_hit_yz->Fill(z, y);
779 }
780 }
781
782
783 h_reco_purity->Fill(100 * prtra->GetMaxPurity());
784
785 if (prtra->GetNdf() > 0) {
786 if (prtra->IsGhost()) {
787 h_ghost_chi2->Fill(prtra->GetChiSq() / prtra->GetNdf());
788 h_ghost_prob->Fill(TMath::Prob(prtra->GetChiSq(), prtra->GetNdf()));
789 }
790 else if (prtra->GetNofMCTracks() > 0) {
791 const auto& mcTrk = fMCData.GetTrack(prtra->GetMatchedMCTrackIndex());
792 if (mcTrk.IsReconstructable()) {
793 h_reco_chi2->Fill(prtra->GetChiSq() / prtra->GetNdf());
794 h_reco_prob->Fill(TMath::Prob(prtra->GetChiSq(), prtra->GetNdf()));
795 }
796 else {
797 // h_rest_chi2->Fill(prtra->GetChiSq()/prtra->NDF);
798 h_rest_prob->Fill(TMath::Prob(prtra->GetChiSq(), prtra->GetNdf()));
799 }
800 }
801 }
802
803
804 // fill ghost histos
805 if (prtra->IsGhost()) {
806 //fMonitor.IncrementCounter(EMonitorKey::kGhostTrack);
807 h_ghost_purity->Fill(100 * prtra->GetMaxPurity());
808 if (fabs(prtra->GetQp()) > 1.e-10) {
809 h_ghost_mom->Fill(prtra->GetP());
810 h_ghost_phi->Fill(prtra->GetPhi()); // phi = atan(py / px) = atan(ty / tx)
811 h_ghost_Rmom->Fill(prtra->GetP());
812 }
813 h_ghost_nhits->Fill((prtra->Hits).size());
814 CbmL1HitDebugInfo& h1 = fvHitDebugInfo[prtra->Hits[0]];
815 CbmL1HitDebugInfo& h2 = fvHitDebugInfo[prtra->Hits[1]];
816 h_ghost_fstation->Fill(h1.iStation);
817 h_ghost_r->Fill(sqrt(fabs(h1.x * h1.x + h1.y * h1.y)));
818 double z1 = fpAlgo->GetParameters().GetStation(h1.iStation).fZ[0];
819 double z2 = fpAlgo->GetParameters().GetStation(h2.iStation).fZ[0];
820 if (fabs(z2 - z1) > 1.e-4) {
821 h_ghost_tx->Fill((h2.x - h1.x) / (z2 - z1));
822 h_ghost_ty->Fill((h2.y - h1.y) / (z2 - z1));
823 }
824
825 CbmL1HitDebugInfo& hf = fvHitDebugInfo[prtra->Hits[0]];
826 CbmL1HitDebugInfo& hl = fvHitDebugInfo[prtra->Hits[(prtra->Hits).size() - 1]];
827 if (fabs(prtra->GetQp()) > 1.e-10) {
828 h2_ghost_nhits_vs_mom->Fill(prtra->GetP(), prtra->Hits.size());
829 h2_ghost_fstation_vs_mom->Fill(prtra->GetP(), hf.iStation + 1);
830 }
831 if (hl.iStation >= hf.iStation) h2_ghost_lstation_vs_fstation->Fill(hf.iStation + 1, hl.iStation + 1);
832 }
833
834 } // for reco tracks
835
836 int mc_total = 0; // total amount of reconstructable mcTracks
837 for (const auto& mcTrk : fMCData.GetTrackContainer()) {
838 // if( !( mcTrk.GetPdgCode() == -11 && mcTrk.GetPdgCode() == -1 ) ) continue; // electrons only
839
840 // No Sts hits?
841 int nmchits = mcTrk.GetNofHits();
842 if (nmchits == 0) continue;
843
844 double momentum = mcTrk.GetP();
845 double pt = mcTrk.GetPt();
846 double theta = mcTrk.GetTheta() * 180 / 3.1415;
847 double phi = mcTrk.GetPhi();
848
849 h_mcp->Fill(momentum);
850 h_nmchits->Fill(nmchits);
851
852 int nSta = mcTrk.GetTotNofStationsWithHit();
853
854 CbmL1HitDebugInfo& fh = fvHitDebugInfo[*(mcTrk.GetHitIndexes().begin())];
855 CbmL1HitDebugInfo& lh = fvHitDebugInfo[*(mcTrk.GetHitIndexes().rbegin())];
856
857 h_reg_MCmom->Fill(momentum);
858 if (mcTrk.IsPrimary()) {
859 h_reg_mom_prim->Fill(momentum);
860 h_reg_phi_prim->Fill(phi);
861 h_reg_prim_MCmom->Fill(momentum);
862 h_reg_nhits_prim->Fill(nSta);
863 h2_reg_nhits_vs_mom_prim->Fill(momentum, nSta);
864 h2_reg_fstation_vs_mom_prim->Fill(momentum, fh.iStation + 1);
865 if (lh.iStation >= fh.iStation) h2_reg_lstation_vs_fstation_prim->Fill(fh.iStation + 1, lh.iStation + 1);
866 }
867 else {
868 h_reg_mom_sec->Fill(momentum);
869 h_reg_phi_sec->Fill(phi);
870 h_reg_sec_MCmom->Fill(momentum);
871 h_reg_nhits_sec->Fill(nSta);
872 if (momentum > 0.01) {
873 h2_reg_nhits_vs_mom_sec->Fill(momentum, nSta);
874 h2_reg_fstation_vs_mom_sec->Fill(momentum, fh.iStation + 1);
875 if (lh.iStation >= fh.iStation) h2_reg_lstation_vs_fstation_sec->Fill(fh.iStation + 1, lh.iStation + 1);
876 }
877 }
878
879 if (mcTrk.IsAdditional()) {
880 h_acc_mom_short123s->Fill(momentum);
881 }
882
883 if (!mcTrk.IsReconstructable()) continue;
884 mc_total++;
885
886 h_acc_MCmom->Fill(momentum);
887 if (mcTrk.IsPrimary()) {
888 h_acc_mom_prim->Fill(momentum);
889 h_acc_phi_prim->Fill(phi);
890 h_acc_prim_MCmom->Fill(momentum);
891 h_acc_nhits_prim->Fill(nSta);
892 h2_acc_nhits_vs_mom_prim->Fill(momentum, nSta);
893 h2_acc_fstation_vs_mom_prim->Fill(momentum, fh.iStation + 1);
894 if (lh.iStation >= fh.iStation) h2_acc_lstation_vs_fstation_prim->Fill(fh.iStation + 1, lh.iStation + 1);
895 }
896 else {
897 h_acc_mom_sec->Fill(momentum);
898 h_acc_phi_sec->Fill(phi);
899 h_acc_sec_MCmom->Fill(momentum);
900 h_acc_nhits_sec->Fill(nSta);
901 if (momentum > 0.01) {
902 h2_acc_nhits_vs_mom_sec->Fill(momentum, nSta);
903 h2_acc_fstation_vs_mom_sec->Fill(momentum, fh.iStation + 1);
904 if (lh.iStation >= fh.iStation) h2_acc_lstation_vs_fstation_sec->Fill(fh.iStation + 1, lh.iStation + 1);
905 }
906 }
907
908
909 // vertex distribution of primary and secondary tracks
910 h2_vertex->Fill(mcTrk.GetStartZ(), mcTrk.GetStartY());
911 //h3_vertex->Fill(mcTrk.GetStartZ(), mcTrk.GetStartX(), mcTrk.GetStartY());
912 if (mcTrk.GetMotherId() < 0) { // primary
913 h2_prim_vertex->Fill(mcTrk.GetStartZ(), mcTrk.GetStartY());
914 //h3_prim_vertex->Fill(mcTrk.GetStartZ(), mcTrk.GetStartX(), mcTrk.GetStartY());
915 }
916 else {
917 h2_sec_vertex->Fill(mcTrk.GetStartZ(), mcTrk.GetStartY());
918 //h3_sec_vertex->Fill(mcTrk.GetStartZ(), mcTrk.GetStartX(), mcTrk.GetStartY());
919 }
920
921
922 int iph = mcTrk.GetHitIndexes()[0];
923 CbmL1HitDebugInfo& ph = fvHitDebugInfo[iph];
924
925 h_sec_r->Fill(sqrt(fabs(ph.x * ph.x + ph.y * ph.y)));
926 if (fabs(mcTrk.GetPz()) > 1.e-8) {
927 h_tx->Fill(mcTrk.GetTx());
928 h_ty->Fill(mcTrk.GetTy());
929 }
930
931 // reconstructed tracks
932 bool reco = (mcTrk.IsReconstructed());
933 int nMCHits = mcTrk.GetTotNofStationsWithHit();
934
935 if (reco) {
936 p_eff_all_vs_mom->Fill(momentum, 100.0);
937 p_eff_all_vs_nhits->Fill(nMCHits, 100.0);
938 p_eff_all_vs_pt->Fill(pt, 100.0);
939 h_reco_MCmom->Fill(momentum);
940 if (mcTrk.GetMotherId() < 0) { // primary
941 p_eff_prim_vs_mom->Fill(momentum, 100.0);
942 p_eff_prim_vs_nhits->Fill(nMCHits, 100.0);
943 p_eff_prim_vs_pt->Fill(pt, 100.0);
944 p_eff_prim_vs_theta->Fill(theta, 100.0);
945 }
946 else {
947 p_eff_sec_vs_mom->Fill(momentum, 100.0);
948 p_eff_sec_vs_nhits->Fill(nMCHits, 100.0);
949 }
950 if (mcTrk.GetMotherId() < 0) { // primary
951 h_reco_mom_prim->Fill(momentum);
952 h_reco_phi_prim->Fill(phi);
953 h_reco_prim_MCmom->Fill(momentum);
954 h_reco_nhits_prim->Fill(nSta);
955 h2_reco_nhits_vs_mom_prim->Fill(momentum, nSta);
956 h2_reco_fstation_vs_mom_prim->Fill(momentum, fh.iStation + 1);
957 if (lh.iStation >= fh.iStation) h2_reco_lstation_vs_fstation_prim->Fill(fh.iStation + 1, lh.iStation + 1);
958 }
959 else {
960 h_reco_mom_sec->Fill(momentum);
961 h_reco_phi_sec->Fill(phi);
962 h_reco_sec_MCmom->Fill(momentum);
963 h_reco_nhits_sec->Fill(nSta);
964 if (momentum > 0.01) {
965 h2_reco_nhits_vs_mom_sec->Fill(momentum, nSta);
966 h2_reco_fstation_vs_mom_sec->Fill(momentum, fh.iStation + 1);
967 if (lh.iStation >= fh.iStation) h2_reco_lstation_vs_fstation_sec->Fill(fh.iStation + 1, lh.iStation + 1);
968 }
969 }
970 }
971 else {
972 h_notfound_mom->Fill(momentum);
973 h_notfound_phi->Fill(phi);
974 p_eff_all_vs_mom->Fill(momentum, 0.0);
975 p_eff_all_vs_nhits->Fill(nMCHits, 0.0);
976 p_eff_all_vs_pt->Fill(pt, 0.0);
977 if (mcTrk.GetMotherId() < 0) { // primary
978 p_eff_prim_vs_mom->Fill(momentum, 0.0);
979 p_eff_prim_vs_nhits->Fill(nMCHits, 0.0);
980 p_eff_prim_vs_pt->Fill(pt, 0.0);
981 p_eff_prim_vs_theta->Fill(theta, 0.0);
982 }
983 else {
984 p_eff_sec_vs_mom->Fill(momentum, 0.0);
985 p_eff_sec_vs_nhits->Fill(nMCHits, 0.0);
986 }
987 if (mcTrk.GetMotherId() < 0) { // primary
988 h_rest_mom_prim->Fill(momentum);
989 h_rest_phi_prim->Fill(phi);
990 h_rest_nhits_prim->Fill(nSta);
991 h2_rest_nhits_vs_mom_prim->Fill(momentum, nSta);
992 h2_rest_fstation_vs_mom_prim->Fill(momentum, fh.iStation + 1);
993 if (lh.iStation >= fh.iStation) h2_rest_lstation_vs_fstation_prim->Fill(fh.iStation + 1, lh.iStation + 1);
994 }
995 else {
996 h_rest_mom_sec->Fill(momentum);
997 h_rest_phi_sec->Fill(phi);
998 h_rest_nhits_sec->Fill(nSta);
999 if (momentum > 0.01) {
1000 h2_rest_nhits_vs_mom_sec->Fill(momentum, nSta);
1001 h2_rest_fstation_vs_mom_sec->Fill(momentum, fh.iStation + 1);
1002 if (lh.iStation >= fh.iStation) h2_rest_lstation_vs_fstation_sec->Fill(fh.iStation + 1, lh.iStation + 1);
1003 }
1004 }
1005 }
1006
1007 // killed tracks. At least one hit of they belong to some recoTrack
1008 // bool killed = 0;
1009 if (!reco) {
1010 h_notfound_nhits->Fill(nmchits);
1011 h_notfound_station->Fill(ph.iStation);
1012 h_notfound_r->Fill(sqrt(fabs(ph.x * ph.x + ph.y * ph.y)));
1013
1014 if (mcTrk.GetNofPoints() > 0) {
1015 const auto& pMC = fMCData.GetPoint(mcTrk.GetPointIndexes()[0]);
1016 h_notfound_tx->Fill(pMC.GetTx());
1017 h_notfound_ty->Fill(pMC.GetTy());
1018 }
1019
1020 // CbmL1HitDebugInfo &ph21 = fvHitDebugInfo[mtra.Hits[0]];
1021 // CbmL1HitDebugInfo &ph22 = fvHitDebugInfo[mtra.Hits[1]];
1022
1023 // double z21 = fpAlgo->GetParameters().GetStation(ph21.iStation).fZ[0];
1024 // double z22 = fpAlgo->GetParameters().GetStation(ph22.iStation).fZ[0];
1025 // if( fabs(z22-z21)>1.e-4 ){
1026 // h_notfound_tx->Fill((ph22.x-ph21.x)/(z22-z21));
1027 // h_notfound_ty->Fill((ph22.y-ph21.y)/(z22-z21));
1028 // }
1029
1030 // if( mtra.IsDisturbed() ) killed = 1;
1031 }
1032
1033
1034 if (mcTrk.IsSignal()) { // D0
1035 h_reco_d0_mom->Fill(momentum);
1036 if (reco)
1037 p_eff_d0_vs_mom->Fill(momentum, 100.0);
1038 else
1039 p_eff_d0_vs_mom->Fill(momentum, 0.0);
1040 }
1041
1042 } // for mcTracks
1043
1044 int NFakes = 0;
1045 for (unsigned int ih = 0; ih < fpAlgo->GetInputData().GetNhits(); ih++) {
1046 int iMC = fvHitDebugInfo[ih].GetBestMcPointId(); // TODO2: adapt to linking
1047 if (iMC < 0) NFakes++;
1048 }
1049
1050 h_reco_time->Fill(fTrackingTime);
1051 h_reco_timeNtr->Fill(mc_total, fTrackingTime);
1052 h_reco_timeNhit->Fill(fpAlgo->GetInputData().GetNhits(), fTrackingTime);
1053
1054 h_reco_fakeNtr->Fill(mc_total, NFakes);
1055 h_reco_fakeNhit->Fill(fpAlgo->GetInputData().GetNhits() - NFakes, NFakes);
1056
1057
1058 // NOTE: Must be called once
1059 h_reg_MCmom->Scale(1.f / NEvents);
1060 h_acc_MCmom->Scale(1.f / NEvents);
1061 h_reco_MCmom->Scale(1.f / NEvents);
1062 h_ghost_Rmom->Scale(1.f / NEvents);
1063 h_reg_prim_MCmom->Scale(1.f / NEvents);
1064 h_acc_prim_MCmom->Scale(1.f / NEvents);
1065 h_reco_prim_MCmom->Scale(1.f / NEvents);
1066 h_reg_sec_MCmom->Scale(1.f / NEvents);
1067 h_acc_sec_MCmom->Scale(1.f / NEvents);
1068 h_reco_sec_MCmom->Scale(1.f / NEvents);
1069
1070} // void CbmL1::HistoPerformance()
1071
1072
1073void CbmL1::TrackFitPerformance()
1074{
1075 const int Nh_fit = 17;
1076 static TH1F *h_fit[Nh_fit], *h_fitL[Nh_fit], *h_fitSV[Nh_fit], *h_fitPV[Nh_fit], *h_fit_chi2;
1077
1078 static TH2F *PRes2D, *PRes2DPrim, *PRes2DSec;
1079
1080 static TH2F* pion_res_pt_fstDca;
1081 static TH2F* kaon_res_pt_fstDca;
1082 static TH2F* pton_res_pt_fstDca;
1083
1084 static TH2F* pion_res_pt_vtxDca;
1085 static TH2F* kaon_res_pt_vtxDca;
1086 static TH2F* pton_res_pt_vtxDca;
1087
1088 static TH2F* pion_res_p_fstDca;
1089 static TH2F* kaon_res_p_fstDca;
1090 static TH2F* pton_res_p_fstDca;
1091
1092 static TH2F* pion_res_p_vtxDca;
1093 static TH2F* kaon_res_p_vtxDca;
1094 static TH2F* pton_res_p_vtxDca;
1095
1096 static bool first_call = 1;
1097
1098 kf::TrackKalmanFilter<fvec> fit;
1099 fit.SetParticleMass(fpAlgo->GetDefaultParticleMass());
1100 fit.SetMask(fmask::One());
1101 //fit.SetMaxExtrapolationStep(10.);
1102 fit.SetDoFitVelocity(true);
1103
1104 if (first_call) {
1105 first_call = 0;
1106
1107 TDirectory* currentDir = gDirectory;
1108 gDirectory = fHistoDir;
1109 gDirectory->cd("Fit");
1110 {
1111 PRes2D = new TH2F("PRes2D", "Resolution P [%] vs P", 100, 0., 20., 100, -15., 15.);
1112 PRes2DPrim = new TH2F("PRes2DPrim", "Resolution P [%] vs P", 100, 0., 20., 100, -15., 15.);
1113 PRes2DSec = new TH2F("PRes2DSec", "Resolution P [%] vs P", 100, 0., 20., 100, -15., 15.);
1114
1115 pion_res_pt_fstDca = new TH2F("pion_res_pt_fstDca", "", 100, 0, 10, 1000, -500, 500);
1116 kaon_res_pt_fstDca = new TH2F("kaon_res_pt_fstDca", "", 100, 0, 10, 1000, -500, 500);
1117 pton_res_pt_fstDca = new TH2F("pton_res_pt_fstDca", "", 100, 0, 10, 1000, -500, 500);
1118
1119 pion_res_pt_vtxDca = new TH2F("pion_res_pt_vtxDca", "", 100, 0, 10, 1000, -5000, 5000);
1120 kaon_res_pt_vtxDca = new TH2F("kaon_res_pt_vtxDca", "", 100, 0, 10, 1000, -5000, 5000);
1121 pton_res_pt_vtxDca = new TH2F("pton_res_pt_vtxDca", "", 100, 0, 10, 1000, -5000, 5000);
1122
1123 pion_res_p_fstDca = new TH2F("pion_res_p_fstDca", "", 100, 0, 10, 1000, -500, 500);
1124 kaon_res_p_fstDca = new TH2F("kaon_res_p_fstDca", "", 100, 0, 10, 1000, -500, 500);
1125 pton_res_p_fstDca = new TH2F("pton_res_p_fstDca", "", 100, 0, 10, 1000, -500, 500);
1126
1127 pion_res_p_vtxDca = new TH2F("pion_res_p_vtxDca", "", 100, 0, 10, 1000, -5000, 5000);
1128 kaon_res_p_vtxDca = new TH2F("kaon_res_p_vtxDca", "", 100, 0, 10, 1000, -5000, 5000);
1129 pton_res_p_vtxDca = new TH2F("pton_res_p_vtxDca", "", 100, 0, 10, 1000, -5000, 5000);
1130
1131 struct {
1132 const char* name;
1133 const char* title;
1134 Int_t n;
1135 Double_t l, r;
1136 } Table[Nh_fit] = {{"x", "Residual X [#mum]", 140, -40., 40.},
1137 {"y", "Residual Y [#mum]", 100, -450., 450.},
1138 {"tx", "Residual Tx [mrad]", 100, -4., 4.},
1139 {"ty", "Residual Ty [mrad]", 100, -3.5, 3.5},
1140 {"P", "Resolution P/Q [%]", 100, -15., 15.},
1141 {"px", "Pull X [residual/estimated_error]", 100, -6., 6.},
1142 {"py", "Pull Y [residual/estimated_error]", 100, -6., 6.},
1143 {"ptx", "Pull Tx [residual/estimated_error]", 100, -6., 6.},
1144 {"pty", "Pull Ty [residual/estimated_error]", 100, -6., 6.},
1145 {"pQP", "Pull Q/P [residual/estimated_error]", 100, -6., 6.},
1146 {"QPreco", "Reco Q/P ", 100, -10., 10.},
1147 {"QPmc", "MC Q/P ", 100, -10., 10.},
1148 {"time", "Residual Time [ns]", 100, -10., 10.},
1149 {"ptime", "Pull Time [residual/estimated_error]", 100, -10., 10.},
1150 {"vi", "Residual Vi [1/c]", 100, -4., 4.},
1151 {"pvi", "Pull Vi [residual/estimated_error]", 100, -10., 10.},
1152 {"distrVi", "Vi distribution [1/c]", 100, 0., 4.}};
1153
1154 if (ca::TrackingMode::kGlobal == fpAlgo->GetTrackingMode()) {
1155 Table[4].l = -100.;
1156 Table[4].r = 100.;
1157 }
1158
1159 struct Tab {
1160 const char* name;
1161 const char* title;
1162 Int_t n;
1163 Double_t l, r;
1164 };
1165 Tab TableVertex[Nh_fit] = {//{"x", "Residual X [cm]", 200, -0.01, 0.01},
1166 {"x", "Residual X [cm]", 2000, -1., 1.},
1167 //{"y", "Residual Y [cm]", 200, -0.01, 0.01},
1168 {"y", "Residual Y [cm]", 2000, -1., 1.},
1169 //{"tx", "Residual Tx [mrad]", 100, -3., 3.},
1170 {"tx", "Residual Tx [mrad]", 100, -2., 2.},
1171 //{"ty", "Residual Ty [mrad]", 100, -3., 3.},
1172 {"ty", "Residual Ty [mrad]", 100, -2., 2.},
1173 {"P", "Resolution P/Q [%]", 100, -15., 15.},
1174 {"px", "Pull X [residual/estimated_error]", 100, -10., 10.},
1175 {"py", "Pull Y [residual/estimated_error]", 100, -10., 10.},
1176 {"ptx", "Pull Tx [residual/estimated_error]", 100, -10., 10.},
1177 {"pty", "Pull Ty [residual/estimated_error]", 100, -10., 10.},
1178 {"pQP", "Pull Q/P [residual/estimated_error]", 100, -10., 10.},
1179 {"QPreco", "Reco Q/P ", 100, -10., 10.},
1180 {"QPmc", "MC Q/P ", 100, -10., 10.},
1181 {"time", "Residual Time [ns]", 100, -10., 10.},
1182 {"ptime", "Pull Time [residual/estimated_error]", 100, -10., 10.},
1183 {"vi", "Residual Vi [1/c]", 100, -10., 10.},
1184 {"pvi", "Pull Vi [residual/estimated_error]", 100, -10., 10.},
1185 {"distrVi", "Vi distribution [1/c]", 100, -10., 10.}};
1186
1187 if (ca::TrackingMode::kGlobal == fpAlgo->GetTrackingMode()) {
1188 TableVertex[4].l = -1.;
1189 TableVertex[4].r = 1.;
1190 }
1191
1192 for (int i = 0; i < Nh_fit; i++) {
1193 char n[225], t[255];
1194 sprintf(n, "fst_%s", Table[i].name);
1195 sprintf(t, "First point %s", Table[i].title);
1196 h_fit[i] = new TH1F(n, t, Table[i].n, Table[i].l, Table[i].r);
1197 sprintf(n, "lst_%s", Table[i].name);
1198 sprintf(t, "Last point %s", Table[i].title);
1199 h_fitL[i] = new TH1F(n, t, Table[i].n, Table[i].l, Table[i].r);
1200 sprintf(n, "svrt_%s", TableVertex[i].name);
1201 sprintf(t, "Secondary vertex point %s", TableVertex[i].title);
1202 h_fitSV[i] = new TH1F(n, t, TableVertex[i].n, TableVertex[i].l, TableVertex[i].r);
1203 sprintf(n, "pvrt_%s", TableVertex[i].name);
1204 sprintf(t, "Primary vertex point %s", TableVertex[i].title);
1205 h_fitPV[i] = new TH1F(n, t, TableVertex[i].n, TableVertex[i].l, TableVertex[i].r);
1206
1207 for (int j = 0; j < 12; j++) {
1208 sprintf(n, "smth_%s_sta_%i", Table[i].name, j);
1209 sprintf(t, "Station %i %s", j, Table[i].title);
1210 // h_smoothed[j][i] = new TH1F(n,t, Table[i].n, Table[i].l, Table[i].r);
1211 }
1212 }
1213 h_fit_chi2 = new TH1F("h_fit_chi2", "Chi2/NDF", 50, -0.5, 10.0);
1214 }
1215 gDirectory = currentDir;
1216 } // if first call
1217
1218
1219 for (vector<CbmL1Track>::iterator it = fvRecoTracks.begin(); it != fvRecoTracks.end(); ++it) {
1220
1221 if (it->IsGhost()) continue;
1222
1223 bool isTimeFitted = 0;
1224
1225 {
1226 int nTimeMeasurenments = 0;
1227 for (unsigned int ih = 0; ih < it->Hits.size(); ih++) {
1228 int ista = fvHitDebugInfo[it->Hits[ih]].iStation;
1229 if (fpAlgo->GetParameters().GetStation(ista).timeInfo) {
1230 nTimeMeasurenments++;
1231 }
1232 }
1233 isTimeFitted = (nTimeMeasurenments > 1);
1234 }
1235
1236 do { // first hit
1237
1238 if (it->GetNofMCTracks() < 1) {
1239 break;
1240 }
1241 const auto& mcTrk = fMCData.GetTrack(it->GetMCTrackIndexes()[0]);
1242 int imcPoint = fvHitDebugInfo[it->Hits.front()].GetBestMcPointId();
1243
1244 // extrapolate to the first mc point of the mc track !
1245
1246 imcPoint = mcTrk.GetPointIndexes()[0];
1247
1248 if (imcPoint < 0) break;
1249
1250 const auto& mcP = fMCData.GetPoint(imcPoint);
1251
1252 TrackParamV tr(*it);
1253 FillFitHistos(tr, mcP, isTimeFitted, h_fit);
1254
1255 double dx = tr.GetX()[0] - mcP.GetXIn();
1256 double dy = tr.GetY()[0] - mcP.GetYIn();
1257 double dca = sqrt(dx * dx + dy * dy);
1258 // make dca distance negative for the half of the tracks to ease the gaussian fit and the rms calculation
1259 if (mcTrk.GetId() % 2) dca = -dca;
1260 double pt = mcP.GetPt();
1261 double dP = (fabs(1. / tr.GetQp()[0]) - mcP.GetP()) / mcP.GetP();
1262
1263 PRes2D->Fill(mcP.GetP(), 100. * dP);
1264
1265 if (mcTrk.IsPrimary()) {
1266
1267 h_fit[4]->Fill(100. * dP);
1268
1269 PRes2DPrim->Fill(mcP.GetP(), 100. * dP);
1270
1271 if (abs(mcTrk.GetPdgCode()) == 211) {
1272 pion_res_p_fstDca->Fill(mcP.GetP(), dca * 1.e4);
1273 pion_res_pt_fstDca->Fill(pt, dca * 1.e4);
1274 }
1275 if (abs(mcTrk.GetPdgCode()) == 321) {
1276 kaon_res_p_fstDca->Fill(mcP.GetP(), dca * 1.e4);
1277 kaon_res_pt_fstDca->Fill(pt, dca * 1.e4);
1278 }
1279 if (abs(mcTrk.GetPdgCode()) == 2212) {
1280 pton_res_p_fstDca->Fill(mcP.GetP(), dca * 1.e4);
1281 pton_res_pt_fstDca->Fill(pt, dca * 1.e4);
1282 }
1283 }
1284 else {
1285 PRes2DSec->Fill(mcP.GetP(), 100. * dP);
1286 }
1287 } while (0);
1288
1289
1290 { // last hit
1291
1292 int iMC = fvHitDebugInfo[it->Hits.back()].GetBestMcPointId(); // TODO2: adapt to linking
1293 if (iMC < 0) continue;
1294 const auto& mcP = fMCData.GetPoint(iMC);
1295 TrackParamV tr(it->TLast);
1296 FillFitHistos(tr, mcP, isTimeFitted, h_fitL);
1297 }
1298
1299 do { // vertex
1300
1301 if (it->GetNofMCTracks() < 1) {
1302 break;
1303 }
1304
1305 const auto& mc = fMCData.GetTrack(it->GetMatchedMCTrackIndex());
1306 fit.SetTrack(*it);
1307
1308 const TrackParamV& tr = fit.Tr();
1309
1310 // if (mc.mother_ID != -1) { // secondary
1311 if (!mc.IsPrimary()) { // secondary
1312
1313 { // extrapolate to vertex
1314 kf::FieldRegion<fvec> fld(kf::GlobalField::fgOriginalFieldType, kf::GlobalField::fgOriginalField);
1315 fit.Extrapolate(mc.GetStartZ(), fld);
1316 // add material
1317 const int fSta = fvHitDebugInfo[it->Hits[0]].iStation;
1318 const int dir = int((mc.GetStartZ() - fpAlgo->GetParameters().GetStation(fSta).fZ[0])
1319 / fabs(mc.GetStartZ() - fpAlgo->GetParameters().GetStation(fSta).fZ[0]));
1320 // if (abs(mc.z - fpAlgo->GetParameters().GetStation(fSta).fZ[0]) > 10.) continue; // can't extrapolate on large distance
1321 for (int iSta = fSta /*+dir*/;
1322 (iSta >= 0) && (iSta < fNStations)
1323 && (dir * (mc.GetStartZ() - fpAlgo->GetParameters().GetStation(iSta).fZ[0]) > 0);
1324 iSta += dir) {
1325 // LOG(info) << iSta << " " << dir;
1326 auto radThick = fpAlgo->GetParameters().GetActiveSetup().GetMaterial(iSta).GetThicknessX0(fit.Tr().GetX(),
1327 fit.Tr().GetY());
1328 fit.MultipleScattering(radThick);
1329 fit.EnergyLossCorrection(radThick, kf::FitDirection::kUpstream);
1330 }
1331 }
1332 if (mc.GetStartZ() != tr.GetZ()[0]) continue;
1333
1334 // static int good = 0;
1335 // static int bad = 0;
1336 // if (mc.z != tr.GetZ()[0]){
1337 // bad++;
1338 // continue;
1339 // }
1340 // else good++;
1341 // LOG(info) << "bad\\good" << bad << " " << good;
1342
1343 // calculate pulls
1344 //h_fitSV[0]->Fill( (mc.x-tr.GetX()[0]) *1.e4);
1345 //h_fitSV[1]->Fill( (mc.y-tr.GetY()[0]) *1.e4);
1346 h_fitSV[0]->Fill((tr.GetX()[0] - mc.GetStartX()));
1347 h_fitSV[1]->Fill((tr.GetY()[0] - mc.GetStartY()));
1348 h_fitSV[2]->Fill((tr.GetTx()[0] - mc.GetTx()) * 1.e3);
1349 h_fitSV[3]->Fill((tr.GetTy()[0] - mc.GetTy()) * 1.e3);
1350 h_fitSV[4]->Fill(100. * (fabs(1. / tr.GetQp()[0]) / mc.GetP() - 1.));
1351 if (std::isfinite((fscal) tr.C00()[0]) && tr.C00()[0] > 0) {
1352 h_fitSV[5]->Fill((tr.GetX()[0] - mc.GetStartX()) / sqrt(tr.C00()[0]));
1353 }
1354 if (std::isfinite((fscal) tr.C11()[0]) && tr.C11()[0] > 0)
1355 h_fitSV[6]->Fill((tr.GetY()[0] - mc.GetStartY()) / sqrt(tr.C11()[0]));
1356 if (std::isfinite((fscal) tr.C22()[0]) && tr.C22()[0] > 0)
1357 h_fitSV[7]->Fill((tr.GetTx()[0] - mc.GetTx()) / sqrt(tr.C22()[0]));
1358 if (std::isfinite((fscal) tr.C33()[0]) && tr.C33()[0] > 0)
1359 h_fitSV[8]->Fill((tr.GetTy()[0] - mc.GetTy()) / sqrt(tr.C33()[0]));
1360 if (std::isfinite((fscal) tr.C44()[0]) && tr.C44()[0] > 0)
1361 h_fitSV[9]->Fill((tr.GetQp()[0] - mc.GetCharge() / mc.GetP()) / sqrt(tr.C44()[0]));
1362 h_fitSV[10]->Fill(tr.GetQp()[0]);
1363 h_fitSV[11]->Fill(mc.GetCharge() / mc.GetP());
1364 if (isTimeFitted) {
1365 h_fitSV[12]->Fill(tr.GetTime()[0] - mc.GetStartT());
1366 if (std::isfinite((fscal) tr.C55()[0]) && tr.C55()[0] > 0.) {
1367 h_fitSV[13]->Fill((tr.GetTime()[0] - mc.GetStartT()) / sqrt(tr.C55()[0]));
1368 }
1369 double dvi =
1370 tr.GetVi()[0]
1371 - sqrt(1. + mc.GetMass() * mc.GetMass() / mc.GetP() / mc.GetP()) * constants::phys::SpeedOfLightInvD;
1372 h_fitSV[14]->Fill(dvi * constants::phys::SpeedOfLightD);
1373 if (std::isfinite((fscal) tr.C66()[0]) && tr.C66()[0] > 0.) {
1374 h_fitSV[15]->Fill(dvi / sqrt(tr.C66()[0]));
1375 }
1376 h_fitSV[16]->Fill(tr.GetVi()[0] * constants::phys::SpeedOfLightD);
1377 }
1378 }
1379 else { // primary
1380
1381 { // extrapolate to vertex
1382 kf::FieldRegion<fvec> fld(kf::GlobalField::fgOriginalFieldType, kf::GlobalField::fgOriginalField);
1383
1384 // add material
1385 const int fSta = fvHitDebugInfo[it->Hits[0]].iStation;
1386
1387 const int dir = (mc.GetStartZ() - fpAlgo->GetParameters().GetStation(fSta).fZ[0])
1388 / abs(mc.GetStartZ() - fpAlgo->GetParameters().GetStation(fSta).fZ[0]);
1389 // if (abs(mc.z - fpAlgo->GetParameters().GetStation(fSta].fZ[0]) > 10.) continue; // can't extrapolate on large distance
1390
1391 for (int iSta = fSta + dir; (iSta >= 0) && (iSta < fNStations)
1392 && (dir * (mc.GetStartZ() - fpAlgo->GetParameters().GetStation(iSta).fZ[0]) > 0);
1393 iSta += dir) {
1394
1395 fit.Extrapolate(fpAlgo->GetParameters().GetStation(iSta).fZ, fld);
1396 auto radThick = fpAlgo->GetParameters().GetActiveSetup().GetMaterial(iSta).GetThicknessX0(fit.Tr().GetX(),
1397 fit.Tr().GetY());
1398 fit.MultipleScattering(radThick);
1399 fit.EnergyLossCorrection(radThick, kf::FitDirection::kUpstream);
1400 }
1401 fit.Extrapolate(mc.GetStartZ(), fld);
1402 }
1403 if (mc.GetStartZ() != tr.GetZ()[0]) continue;
1404
1405 double dx = tr.GetX()[0] - mc.GetStartX();
1406 double dy = tr.GetY()[0] - mc.GetStartY();
1407 double dt = sqrt(dx * dx + dy * dy);
1408 // make dt distance negative for the half of the tracks to ease the gaussian fit and the rms calculation
1409 if (mc.GetId() % 2) dt = -dt;
1410
1411 double pt = mc.GetPt();
1412
1413 if (abs(mc.GetPdgCode()) == 211) {
1414 pion_res_p_vtxDca->Fill(mc.GetP(), dt * 1.e4);
1415 pion_res_pt_vtxDca->Fill(pt, dt * 1.e4);
1416 }
1417 if (abs(mc.GetPdgCode()) == 321) {
1418 kaon_res_p_vtxDca->Fill(mc.GetPdgCode(), dt * 1.e4);
1419 kaon_res_pt_vtxDca->Fill(mc.GetPt(), dt * 1.e4);
1420 }
1421 if (abs(mc.GetPdgCode()) == 2212) {
1422 pton_res_p_vtxDca->Fill(mc.GetP(), dt * 1.e4);
1423 pton_res_pt_vtxDca->Fill(pt, dt * 1.e4);
1424 }
1425
1426 // calculate pulls
1427 h_fitPV[0]->Fill((mc.GetStartX() - tr.GetX()[0]));
1428 h_fitPV[1]->Fill((mc.GetStartY() - tr.GetY()[0]));
1429 h_fitPV[2]->Fill((mc.GetTx() - tr.GetTx()[0]) * 1.e3);
1430 h_fitPV[3]->Fill((mc.GetTy() - tr.GetTy()[0]) * 1.e3);
1431 h_fitPV[4]->Fill(100. * (fabs(1 / tr.GetQp()[0]) / mc.GetP() - 1.));
1432 if (std::isfinite((fscal) tr.C00()[0]) && tr.C00()[0] > 0)
1433 h_fitPV[5]->Fill((mc.GetStartX() - tr.GetX()[0]) / sqrt(tr.C00()[0]));
1434 if (std::isfinite((fscal) tr.C11()[0]) && tr.C11()[0] > 0)
1435 h_fitPV[6]->Fill((mc.GetStartY() - tr.GetY()[0]) / sqrt(tr.C11()[0]));
1436 if (std::isfinite((fscal) tr.C22()[0]) && tr.C22()[0] > 0)
1437 h_fitPV[7]->Fill((mc.GetTx() - tr.GetTx()[0]) / sqrt(tr.C22()[0]));
1438 if (std::isfinite((fscal) tr.C33()[0]) && tr.C33()[0] > 0)
1439 h_fitPV[8]->Fill((mc.GetTy() - tr.GetTy()[0]) / sqrt(tr.C33()[0]));
1440 if (std::isfinite((fscal) tr.C44()[0]) && tr.C44()[0] > 0)
1441 h_fitPV[9]->Fill((mc.GetCharge() / mc.GetP() - tr.GetQp()[0]) / sqrt(tr.C44()[0]));
1442 h_fitPV[10]->Fill(tr.GetQp()[0]);
1443 h_fitPV[11]->Fill(mc.GetCharge() / mc.GetP());
1444 if (isTimeFitted) {
1445 h_fitPV[12]->Fill(tr.GetTime()[0] - mc.GetStartT());
1446 if (std::isfinite((fscal) tr.C55()[0]) && tr.C55()[0] > 0.) {
1447 h_fitPV[13]->Fill((tr.GetTime()[0] - mc.GetStartT()) / sqrt(tr.C55()[0]));
1448 }
1449 double dvi =
1450 tr.GetVi()[0]
1451 - sqrt(1. + mc.GetMass() * mc.GetMass() / mc.GetP() / mc.GetP()) * constants::phys::SpeedOfLightInvD;
1452 h_fitPV[14]->Fill(dvi * constants::phys::SpeedOfLightD);
1453 if (std::isfinite((fscal) tr.C66()[0]) && tr.C66()[0] > 0.) {
1454 h_fitPV[15]->Fill(dvi / sqrt(tr.C66()[0]));
1455 }
1456 h_fitPV[16]->Fill(tr.GetVi()[0] * constants::phys::SpeedOfLightD);
1457 }
1458 }
1459 } while (0);
1460
1461 h_fit_chi2->Fill(it->GetChiSq() / it->GetNdf());
1462
1463 // last TRD point
1464 /*
1465 do {
1466 break; // only produce these plots in debug mode
1467 if (it->GetNMCTracks() < 1) { break; }
1468
1469 if (!fpTrdPoints) break;
1470 const CbmL1MCTrack& mcTrack = *(it->GetMCTracks()[0]);
1471 int nTrdPoints = fpTrdPoints->Size(mcTrack.iFile, mcTrack.iEvent);
1472 int lastP = -1;
1473 double lastZ = -1.e8;
1474 for (int iPoint = 0; iPoint < nTrdPoints; iPoint++) {
1475 const CbmTrdPoint* pt = dynamic_cast<CbmTrdPoint*>(fpTrdPoints->Get(mcTrack.iFile, mcTrack.iEvent, iPoint));
1476 if (!pt) { continue; }
1477 if (pt->GetTrackID() != mcTrack.ID) { continue; }
1478 if (lastP < 0 || lastZ < pt->GetZ()) {
1479 lastP = iPoint;
1480 lastZ = pt->GetZ();
1481 }
1482 }
1483 CbmL1MCPoint mcP;
1484 if (CbmL1::ReadMCPoint(&mcP, lastP, mcTrack.iFile, mcTrack.iEvent, 3)) { break; }
1485 TrackParamV tr;
1486 tr.copyFromArrays(it->TLast, it->CLast);
1487 FillFitHistos(tr, mcP, isTimeFitted, h_fitTrd);
1488 } while (0); // Trd point
1489
1490
1491 // last TOF point
1492 do {
1493 break; // only produce these plots in debug mode
1494
1495 if (it->GetNMCTracks() < 1) { break; }
1496
1497 if (!fpTofPoints) break;
1498 const CbmL1MCTrack& mcTrack = *(it->GetMCTracks()[0]);
1499 int nTofPoints = fpTofPoints->Size(mcTrack.iFile, mcTrack.iEvent);
1500 int lastP = -1;
1501 double lastZ = -1.e8;
1502 for (int iPoint = 0; iPoint < nTofPoints; iPoint++) {
1503 const CbmTofPoint* pt = dynamic_cast<CbmTofPoint*>(fpTofPoints->Get(mcTrack.iFile, mcTrack.iEvent, iPoint));
1504 if (!pt) { continue; }
1505 if (pt->GetTrackID() != mcTrack.ID) { continue; }
1506 if (lastP < 0 || lastZ < pt->GetZ()) {
1507 lastP = iPoint;
1508 lastZ = pt->GetZ();
1509 }
1510 }
1511 CbmL1MCPoint mcP;
1512 if (CbmL1::ReadMCPoint(&mcP, lastP, mcTrack.iFile, mcTrack.iEvent, 4)) { break; }
1513 TrackParamV tr;
1514 tr.copyFromArrays(it->TLast, it->CLast);
1515 FillFitHistos(tr, mcP, isTimeFitted, h_fitTof);
1516 } while (0); // tof point
1517*/
1518 } // reco track
1519
1520} // void CbmL1::TrackFitPerformance()
1521
1522
1523void CbmL1::FillFitHistos(TrackParamV& track, const cbm::ca::tools::MCPoint& mcP, bool isTimeFitted, TH1F* h[])
1524{
1525 kf::TrackKalmanFilter<fvec> fit;
1526 fit.SetParticleMass(fpAlgo->GetDefaultParticleMass());
1527 fit.SetMask(fmask::One());
1528 //fit.SetMaxExtrapolationStep(10.);
1529 fit.SetDoFitVelocity(true);
1530 fit.SetTrack(track);
1531 kf::FieldRegion<fvec> fld(kf::EFieldType::Normal, kf::GlobalField::fgOriginalField);
1532 fit.Extrapolate(mcP.GetZOut(), fld);
1533 track = fit.Tr();
1534
1535 const TrackParamV& tr = track;
1536
1537 h[0]->Fill((tr.GetX()[0] - mcP.GetXOut()) * 1.e4);
1538 h[1]->Fill((tr.GetY()[0] - mcP.GetYOut()) * 1.e4);
1539 h[2]->Fill((tr.GetTx()[0] - mcP.GetTxOut()) * 1.e3);
1540 h[3]->Fill((tr.GetTy()[0] - mcP.GetTyOut()) * 1.e3);
1541 h[4]->Fill(100. * (fabs(1. / tr.GetQp()[0]) / mcP.GetP() - 1.));
1542
1543 if (std::isfinite((fscal) tr.C00()[0]) && tr.C00()[0] > 0)
1544 h[5]->Fill((tr.GetX()[0] - mcP.GetXOut()) / sqrt(tr.C00()[0]));
1545 if (std::isfinite((fscal) tr.C11()[0]) && tr.C11()[0] > 0)
1546 h[6]->Fill((tr.GetY()[0] - mcP.GetYOut()) / sqrt(tr.C11()[0]));
1547 if (std::isfinite((fscal) tr.C22()[0]) && tr.C22()[0] > 0)
1548 h[7]->Fill((tr.GetTx()[0] - mcP.GetTxOut()) / sqrt(tr.C22()[0]));
1549 if (std::isfinite((fscal) tr.C33()[0]) && tr.C33()[0] > 0)
1550 h[8]->Fill((tr.GetTy()[0] - mcP.GetTyOut()) / sqrt(tr.C33()[0]));
1551 if (std::isfinite((fscal) tr.C44()[0]) && tr.C44()[0] > 0)
1552 h[9]->Fill((tr.GetQp()[0] - mcP.GetQp()) / sqrt(tr.C44()[0]));
1553 h[10]->Fill(tr.GetQp()[0]);
1554 h[11]->Fill(mcP.GetQp());
1555 if (isTimeFitted) {
1556 h[12]->Fill(tr.GetTime()[0] - mcP.GetTime());
1557 if (std::isfinite((fscal) tr.C55()[0]) && tr.C55()[0] > 0.) {
1558 h[13]->Fill((tr.GetTime()[0] - mcP.GetTime()) / sqrt(tr.C55()[0]));
1559 }
1560
1561 double dvi = tr.GetVi()[0] - mcP.GetInvSpeed();
1562 h[14]->Fill(dvi * constants::phys::SpeedOfLightD);
1563 if (std::isfinite((fscal) tr.C66()[0]) && tr.C66()[0] > 0.) {
1564 h[15]->Fill(dvi / sqrt(tr.C66()[0]));
1565 }
1566 h[16]->Fill(tr.GetVi()[0] * constants::phys::SpeedOfLightD);
1567 }
1568}
1569
1570
1571void CbmL1::FieldApproxCheck()
1572{
1573 TDirectory* curr = gDirectory;
1574 TFile* currentFile = gFile;
1575 TFile* fout = new TFile("CaFieldApproxQa.root", "RECREATE");
1576 fout->cd();
1577
1578 assert(FairRunAna::Instance());
1579 FairField* MF = FairRunAna::Instance()->GetField();
1580 assert(MF);
1581
1582 for (int ist = 0; ist < fpAlgo->GetParameters().GetNstationsActive(); ist++) {
1583
1584 const ca::Station<fvec>& st = fpAlgo->GetParameters().GetStation(ist);
1585
1586 double z = st.fZ[0];
1587 double Xmax = st.Xmax[0];
1588 double Ymax = st.Ymax[0];
1589
1590 int NbinsX = 101;
1591 int NbinsY = 101;
1592
1593 TH2F* hdB[4] = {
1594 new TH2F(Form("station_%i_dB", ist), Form("station %i, dB, z = %0.f cm", ist, z), //
1595 NbinsX, -Xmax, Xmax, NbinsY, -Ymax, Ymax),
1596 new TH2F(Form("station_%i_dBx", ist), Form("station %i, dBx, z = %0.f cm", ist, z), //
1597 NbinsX, -Xmax, Xmax, NbinsY, -Ymax, Ymax),
1598 new TH2F(Form("station_%i_dBy", ist), Form("station %i, dBy, z = %0.f cm", ist, z), //
1599 NbinsX, -Xmax, Xmax, NbinsY, -Ymax, Ymax),
1600 new TH2F(Form("station_%i_dBz", ist), Form("station %i, dBz, z = %0.f cm", ist, z), //
1601 NbinsX, -Xmax, Xmax, NbinsY, -Ymax, Ymax) //
1602 };
1603
1604 TH2F* hB[4] = {
1605 new TH2F(Form("station_%i_B", ist), Form("station %i, B, z = %0.f cm", ist, z), //
1606 NbinsX, -Xmax, Xmax, NbinsY, -Ymax, Ymax),
1607 new TH2F(Form("station_%i_Bx", ist), Form("station %i, Bx, z = %0.f cm", ist, z), //
1608 NbinsX, -Xmax, Xmax, NbinsY, -Ymax, Ymax),
1609 new TH2F(Form("station_%i_By", ist), Form("station %i, By, z = %0.f cm", ist, z), //
1610 NbinsX, -Xmax, Xmax, NbinsY, -Ymax, Ymax),
1611 new TH2F(Form("station_%i_Bz", ist), Form("station %i, Bz, z = %0.f cm", ist, z), //
1612 NbinsX, -Xmax, Xmax, NbinsY, -Ymax, Ymax) //
1613 };
1614
1615 for (int i = 0; i < 4; i++) {
1616 hdB[i]->GetXaxis()->SetTitle("X, cm");
1617 hdB[i]->GetYaxis()->SetTitle("Y, cm");
1618 hdB[i]->GetXaxis()->SetTitleOffset(1);
1619 hdB[i]->GetYaxis()->SetTitleOffset(1);
1620 hdB[i]->GetZaxis()->SetTitleOffset(1.3);
1621 hB[i]->GetXaxis()->SetTitle("X, cm");
1622 hB[i]->GetYaxis()->SetTitle("Y, cm");
1623 hB[i]->GetXaxis()->SetTitleOffset(1);
1624 hB[i]->GetYaxis()->SetTitleOffset(1);
1625 hB[i]->GetZaxis()->SetTitleOffset(1.3);
1626 }
1627
1628 hdB[0]->GetZaxis()->SetTitle("B_map - B_L1, kGauss");
1629 hdB[1]->GetZaxis()->SetTitle("Bx_map - Bx_L1, kGauss");
1630 hdB[2]->GetZaxis()->SetTitle("By_map - By_L1, kGauss");
1631 hdB[3]->GetZaxis()->SetTitle("Bz_map - Bz_L1, kGauss");
1632
1633 hdB[0]->GetZaxis()->SetTitle("B_map, kGauss");
1634 hdB[1]->GetZaxis()->SetTitle("Bx_map, kGauss");
1635 hdB[2]->GetZaxis()->SetTitle("By_map, kGauss");
1636 hdB[3]->GetZaxis()->SetTitle("Bz_map, kGauss");
1637
1638
1639 for (int i = 1; i <= hdB[0]->GetXaxis()->GetNbins(); i++) {
1640 double x = hdB[0]->GetXaxis()->GetBinCenter(i);
1641 for (int j = 1; j <= hdB[0]->GetYaxis()->GetNbins(); j++) {
1642 double y = hdB[0]->GetYaxis()->GetBinCenter(j);
1643 double r[] = {x, y, z};
1644 double B[] = {0., 0., 0.};
1645 MF->GetFieldValue(r, B);
1646
1647 double bbb = sqrt(B[0] * B[0] + B[1] * B[1] + B[2] * B[2]);
1648
1649 kf::FieldValue<fvec> B_L1 = st.fieldSlice.GetFieldValue(x, y);
1650
1651 hdB[0]->SetBinContent(i, j, (bbb - B_L1.GetAbs()[0]));
1652 hdB[1]->SetBinContent(i, j, (B[0] - B_L1.GetBx()[0]));
1653 hdB[2]->SetBinContent(i, j, (B[1] - B_L1.GetBy()[0]));
1654 hdB[3]->SetBinContent(i, j, (B[2] - B_L1.GetBz()[0]));
1655
1656 hB[0]->SetBinContent(i, j, bbb);
1657 hB[1]->SetBinContent(i, j, B[0]);
1658 hB[2]->SetBinContent(i, j, B[1]);
1659 hB[3]->SetBinContent(i, j, B[2]);
1660 }
1661 }
1662
1663 for (int i = 0; i < 4; i++) {
1664 hdB[i]->Write();
1665 hB[i]->Write();
1666 }
1667
1668 } // for ista
1669
1670 fout->Close();
1671 fout->Delete();
1672 gFile = currentFile;
1673 gDirectory = curr;
1674} // void CbmL1::FieldApproxCheck()
1675
1676
1677#include "TMath.h"
1678void CbmL1::FieldIntegralCheck()
1679{
1680 TDirectory* curr = gDirectory;
1681 TFile* currentFile = gFile;
1682 TFile* fout = new TFile("FieldApprox.root", "RECREATE");
1683 fout->cd();
1684
1685 assert(FairRunAna::Instance());
1686 FairField* MF = FairRunAna::Instance()->GetField();
1687 assert(MF);
1688
1689 int nPointsZ = 1000;
1690 int nPointsPhi = 100;
1691 int nPointsTheta = 100;
1692 double startZ = 0, endZ = 100.;
1693 double startPhi = 0, endPhi = 2 * TMath::Pi();
1694 double startTheta = -30. / 180. * TMath::Pi(), endTheta = 30. / 180. * TMath::Pi();
1695
1696 double DZ = endZ - startZ;
1697 double DP = endPhi - startPhi;
1698 double DT = endTheta - startTheta;
1699
1700 float ddp = endPhi / nPointsPhi;
1701 float ddt = 2 * endTheta / nPointsTheta;
1702
1703 TH2F* hSb =
1704 new TH2F("Field Integral", "Field Integral", static_cast<int>(nPointsPhi), -(startPhi + ddp / 2.),
1705 (endPhi + ddp / 2.), static_cast<int>(nPointsTheta), (startTheta - ddt / 2.), (endTheta + ddt / 2.));
1706
1707 for (int iP = 0; iP < nPointsPhi; iP++) {
1708 double phi = startPhi + iP * DP / nPointsPhi;
1709 for (int iT = 0; iT < nPointsTheta; iT++) {
1710 double theta = startTheta + iT * DT / nPointsTheta;
1711
1712 double Sb = 0;
1713 for (int iZ = 1; iZ < nPointsZ; iZ++) {
1714 double z = startZ + DZ * iZ / nPointsZ;
1715 double x = z * TMath::Tan(theta) * TMath::Cos(phi);
1716 double y = z * TMath::Tan(theta) * TMath::Sin(phi);
1717 double r[3] = {x, y, z};
1718 double b[3];
1719 MF->GetFieldValue(r, b);
1720 double B = sqrt(b[0] * b[0] + b[1] * b[1] + b[2] * b[2]);
1721 Sb += B * DZ / nPointsZ / 100. / 10.;
1722 }
1723 hSb->SetBinContent(iP + 1, iT + 1, Sb);
1724 }
1725 }
1726
1727 hSb->GetXaxis()->SetTitle("#phi [rad]");
1728 hSb->GetYaxis()->SetTitle("#theta [rad]");
1729 hSb->GetXaxis()->SetTitleOffset(1);
1730 hSb->GetYaxis()->SetTitleOffset(1);
1731 hSb->GetZaxis()->SetTitle("Field Integral [T#dotm]");
1732 hSb->GetZaxis()->SetTitleOffset(1.3);
1733
1734 hSb->Write();
1735
1736
1737 fout->Close();
1738 fout->Delete();
1739 gFile = currentFile;
1740 gDirectory = curr;
1741} // void CbmL1::FieldApproxCheck()
1742
1743// ---------------------------------------------------------------------------------------------------------------------
1744//
1745void CbmL1::DumpMCTripletsToTree()
1746{
1747 if (!fpMcTripletsTree) {
1748 TDirectory* currentDir = gDirectory;
1749 TFile* currentFile = gFile;
1750
1751
1752 // Get prefix and directory
1753 boost::filesystem::path p = (FairRunAna::Instance()->GetUserOutputFileName()).Data();
1754 std::string dir = p.parent_path().string();
1755 if (dir.empty()) dir = ".";
1756 std::string filename = dir + "/" + fsMcTripletsOutputFilename + "." + p.filename().string();
1757 LOG(info) << "\033[1;31mSAVING A TREE: " << filename << "\033[0m";
1758
1759 fpMcTripletsOutFile = new TFile(filename.c_str(), "RECREATE");
1760 fpMcTripletsOutFile->cd();
1761 fpMcTripletsTree = new TTree("t", "MC Triplets");
1762 // motherId - id of mother particle (< 0, if primary)
1763 // pdg - PDG code of particle
1764 // processId - id of the process (ROOT::TMCProcessID)
1765 // pq - absolute value of track momentum [GeV/c], divided by charge [e]
1766 // zv - z-component of track vertex [cm]
1767 // s - global index of station
1768 // x0, y0, z0 - position of the left MC point in a triplet [cm]
1769 // x1, y1, z1 - position of the middle MC point in a triplet [cm]
1770 // x2, y2, z2 - position of the right MC point in a triplet [cm]
1771
1772 gFile = currentFile;
1773 gDirectory = currentDir;
1774 }
1775
1776 // Variables for tree branches
1777 int brMotherId = -1;
1778 int brPdg = -1;
1779 unsigned int brProcId = (unsigned int) -1;
1780 float brP = 0.f;
1781 int brQ = 0;
1782 float brVertexZ = 0.f;
1783 int brStation = -1;
1784
1785 float brX0 = 0.f;
1786 float brY0 = 0.f;
1787 float brZ0 = 0.f;
1788 float brX1 = 0.f;
1789 float brY1 = 0.f;
1790 float brZ1 = 0.f;
1791 float brX2 = 0.f;
1792 float brY2 = 0.f;
1793 float brZ2 = 0.f;
1794
1795 // Register branches
1796 fpMcTripletsTree->Branch("brMotherId", &brMotherId, "motherId/I");
1797 fpMcTripletsTree->Branch("brPdg", &brPdg, "pdg/I");
1798 fpMcTripletsTree->Branch("brProcId", &brProcId, "processId/i");
1799 fpMcTripletsTree->Branch("brP", &brP, "p/F");
1800 fpMcTripletsTree->Branch("brQ", &brQ, "q/I");
1801 fpMcTripletsTree->Branch("brVertexZ", &brVertexZ, "zv/F");
1802 fpMcTripletsTree->Branch("brStation", &brStation, "s/I");
1803 fpMcTripletsTree->Branch("brX0", &brX0, "x0/F");
1804 fpMcTripletsTree->Branch("brY0", &brY0, "y0/F");
1805 fpMcTripletsTree->Branch("brZ0", &brZ0, "z0/F");
1806 fpMcTripletsTree->Branch("brX1", &brX1, "x1/F");
1807 fpMcTripletsTree->Branch("brY1", &brY1, "y1/F");
1808 fpMcTripletsTree->Branch("brZ1", &brZ1, "z1/F");
1809 fpMcTripletsTree->Branch("brX2", &brX2, "x2/F");
1810 fpMcTripletsTree->Branch("brY2", &brY2, "y2/F");
1811 fpMcTripletsTree->Branch("brZ2", &brZ2, "z2/F");
1812
1813 struct ReducedMcPoint {
1814 int s;
1815 float x;
1816 float y;
1817 float z;
1818 };
1819
1820 for (const auto& tr : fMCData.GetTrackContainer()) {
1821 // Use only reconstructable tracks
1822 if (!tr.IsReconstructable()) {
1823 continue;
1824 }
1825
1826 std::vector<ReducedMcPoint> vPoints;
1827 vPoints.reserve(tr.GetNofPoints());
1828 for (unsigned int iP : tr.GetPointIndexes()) {
1829 const auto& point = fMCData.GetPoint(iP);
1830 vPoints.emplace_back(
1831 ReducedMcPoint{point.GetStationId(), float(point.GetX()), float(point.GetY()), float(point.GetZ())});
1832 }
1833
1834 std::sort(vPoints.begin(), vPoints.end(),
1835 [](const ReducedMcPoint& lhs, const ReducedMcPoint& rhs) { return lhs.s < rhs.s; });
1836
1837 for (unsigned int i = 0; i + 2 < vPoints.size(); ++i) {
1838 // Condition to collect only triplets without gaps in stations
1839 // TODO: SZh 20.10.2022 Add cases for jump iterations
1840 if (vPoints[i + 1].s == vPoints[i].s + 1 && vPoints[i + 2].s == vPoints[i].s + 2) {
1841 // Fill MC-triplets tree
1842 brMotherId = tr.GetMotherId();
1843 brPdg = tr.GetPdgCode();
1844 brProcId = tr.GetProcessId();
1845 brP = tr.GetP();
1846 brQ = tr.GetCharge();
1847 brVertexZ = tr.GetStartZ();
1848 brStation = vPoints[i].s;
1849 brX0 = vPoints[i].x;
1850 brY0 = vPoints[i].y;
1851 brZ0 = vPoints[i].z;
1852 brX1 = vPoints[i + 1].x;
1853 brY1 = vPoints[i + 1].y;
1854 brZ1 = vPoints[i + 1].z;
1855 brX2 = vPoints[i + 2].x;
1856 brY2 = vPoints[i + 2].y;
1857 brZ2 = vPoints[i + 2].z;
1858
1859 fpMcTripletsTree->Fill();
1860 }
1861 }
1862 }
1863}
CaToolsDebugger.h
Tracking Debugger class (implementation)
CbmMuchPixelHit.h
Class for pixel hits in MUCH detector.
CbmQaTable.h
Definition of CbmQaTable class.
CbmStsHit.h
Data class for a reconstructed hit in the STS.
CbmTrdHit.h
Class for hits in TRD detector.
sqrt
friend fvec sqrt(const fvec &a)
Definition KfSimdPseudo.h:118
size
static constexpr size_t size()
Definition KfSimdPseudo.h:2
KfTrackKalmanFilter.h
Track fit utilities for the CA tracking based on the Kalman filter.
h
Generates beam ions for transport simulation.
CbmL1HitDebugInfo::dy
double dy
y coordinate error [cm]
Definition CbmL1Hit.h:151
CbmL1HitDebugInfo::y
double y
y coordinate of position [cm]
Definition CbmL1Hit.h:147
CbmL1HitDebugInfo::iStation
int iStation
index of station in active stations array
Definition CbmL1Hit.h:144
CbmL1HitDebugInfo::dx
double dx
x coordinate error [cm]
Definition CbmL1Hit.h:150
CbmL1HitDebugInfo::x
double x
x coordinate of position [cm]
Definition CbmL1Hit.h:146
CbmL1Track::GetMaxPurity
double GetMaxPurity() const
Gets max purity.
Definition CbmL1Track.h:71
CbmL1Track::Hits
std::vector< int > Hits
Indexes of hits of this track.
Definition CbmL1Track.h:104
CbmL1Track::IsGhost
bool IsGhost() const
Definition CbmL1Track.h:55
CbmL1Track::GetMatchedMCTrackIndex
int GetMatchedMCTrackIndex() const
Gets index of matched MC track.
Definition CbmL1Track.h:76
CbmL1Track::GetFirstHitIndex
int GetFirstHitIndex() const
Gets first hit index.
Definition CbmL1Track.h:58
CbmL1Track::GetNofMCTracks
int GetNofMCTracks() const
Gets number of associated MC tracks.
Definition CbmL1Track.h:82
CbmL1Track::GetLastHitIndex
int GetLastHitIndex() const
Gets last hit index.
Definition CbmL1Track.h:61
CbmL1::fNStations
int fNStations
number of total active detector stations
Definition CbmL1.h:441
CbmL1::fHistoDir
TDirectory * fHistoDir
Definition CbmL1.h:500
CbmL1::fTrackingTime
double fTrackingTime
time of track finding procedure
Definition CbmL1.h:444
CbmL1::EfficienciesPerformance
void EfficienciesPerformance(bool doFinish=kFALSE)
Calculates tracking efficiencies (counters)
Definition CbmL1Performance.cxx:231
CbmL1::fpMcTripletsTree
TTree * fpMcTripletsTree
Tree to save MC-triplets.
Definition CbmL1.h:510
CbmL1::fsMcTripletsOutputFilename
std::string fsMcTripletsOutputFilename
Name of file to save MC-triplets tree.
Definition CbmL1.h:511
CbmL1::fMCData
cbm::ca::tools::MCData fMCData
MC Data object.
Definition CbmL1.h:417
CbmL1::fpAlgo
ca::Framework * fpAlgo
Pointer to the L1 track finder algorithm.
Definition CbmL1.h:425
CbmL1::fTableDir
TDirectory * fTableDir
Definition CbmL1.h:501
CbmL1::FieldIntegralCheck
void FieldIntegralCheck()
Definition CbmL1Performance.cxx:1678
CbmL1::fvRecoTracks
ca::Vector< CbmL1Track > fvRecoTracks
Reconstructed tracks container.
Definition CbmL1.h:429
CbmL1::FillFitHistos
void FillFitHistos(cbm::algo::kf::TrackParamV &tr, const cbm::ca::tools::MCPoint &mc, bool isTimeFitted, TH1F *h[])
Definition CbmL1Performance.cxx:1523
CbmL1::TrackFitPerformance
void TrackFitPerformance()
Definition CbmL1Performance.cxx:1073
CbmL1::fvHitDebugInfo
ca::Vector< CbmL1HitDebugInfo > fvHitDebugInfo
Container of hits with extended information.
Definition CbmL1.h:487
CbmL1::DumpMCTripletsToTree
void DumpMCTripletsToTree()
Definition CbmL1Performance.cxx:1745
CbmL1::HistoPerformance
void HistoPerformance()
Fills performance histograms.
Definition CbmL1Performance.cxx:432
CbmL1::FieldApproxCheck
void FieldApproxCheck()
Definition CbmL1Performance.cxx:1571
CbmL1::fpMcTripletsOutFile
TFile * fpMcTripletsOutFile
File to save MC-triplets tree.
Definition CbmL1.h:509
CbmQaTable
TODO: SZh, 30.01.2023: Override THistPainter::PaintText() to add zeroes in tables.
Definition CbmQaTable.h:24
CbmQaTable::SetNamesOfCols
void SetNamesOfCols(const std::vector< std::string > &names)
Sets the names of table columns.
Definition CbmQaTable.cxx:72
CbmQaTable::SetNamesOfRows
void SetNamesOfRows(const std::vector< std::string > &names)
Sets the names of table rows.
Definition CbmQaTable.cxx:83
CbmQaTable::SetCell
void SetCell(Int_t iRow, Int_t iCol, Double_t content, Double_t error=0.)
Definition CbmQaTable.cxx:64
CbmQaTable::SetColWidth
void SetColWidth(Int_t width)
Definition CbmQaTable.h:74
cbm::algo::ca::Framework::GetInputData
const InputData & GetInputData() const
Gets pointer to input data object for external access.
Definition CaFramework.h:96
cbm::algo::ca::Framework::GetTrackingMode
TrackingMode GetTrackingMode()
Definition CaFramework.h:136
cbm::algo::ca::Framework::GetDefaultParticleMass
fscal GetDefaultParticleMass() const
Definition CaFramework.h:110
cbm::algo::ca::Framework::GetParameters
const Parameters< fvec > & GetParameters() const
Gets a pointer to the Framework parameters object.
Definition CaFramework.h:87
cbm::algo::ca::Hit
ca::Hit class describes a generic hit for the CA tracker
Definition CaHit.h:32
cbm::algo::ca::InputData::GetNhits
HitIndex_t GetNhits() const
Gets number of hits in the hits vector.
Definition CaInputData.h:61
cbm::algo::ca::InputData::GetHit
const Hit & GetHit(HitIndex_t index) const
Gets reference to hit by its index.
Definition CaInputData.h:55
cbm::algo::ca::Station::Xmax
DataT Xmax
min radius of the station [cm]
Definition CaStation.h:30
cbm::algo::ca::Station::fZ
DataT fZ
z position of station [cm]
Definition CaStation.h:29
cbm::algo::ca::Station::Ymax
DataT Ymax
max radius of the station [cm]
Definition CaStation.h:31
cbm::algo::ca::Station::fieldSlice
kf::FieldSlice< DataT > fieldSlice
Magnetic field near the station.
Definition CaStation.h:33
cbm::algo::kf::FieldRegion
Magnetic field region, corresponding to a hit triplet.
Definition KfFieldRegion.h:257
cbm::algo::kf::FieldSlice::GetFieldValue
constexpr FieldValue< T > GetFieldValue(const T &x, const T &y) const
Gets field value at a point on the transverse plane.
Definition KfFieldSlice.h:73
cbm::algo::kf::FieldValue
Magnetic flux density vector.
Definition KfFieldValue.h:31
cbm::algo::kf::FieldValue::GetBz
T GetBz() const
Gets magnetic flux density z-component [kG].
Definition KfFieldValue.h:108
cbm::algo::kf::FieldValue::GetBx
T GetBx() const
Gets magnetic flux x-component [kG].
Definition KfFieldValue.h:102
cbm::algo::kf::FieldValue::GetBy
T GetBy() const
Gets magnetic flux density y-component [kG].
Definition KfFieldValue.h:105
cbm::algo::kf::FieldValue::GetAbs
T GetAbs() const
Gets absolute magnetic flux [kG].
Definition KfFieldValue.h:95
cbm::algo::kf::GlobalField::fgOriginalFieldType
static EFieldType fgOriginalFieldType
Global field type.
Definition KfFieldRegion.h:39
cbm::algo::kf::GlobalField::fgOriginalField
static FieldFn_t fgOriginalField
Global variable to store the fielf funciton (x,y,z)->(Bx,By,Bz)
Definition KfFieldRegion.h:36
cbm::algo::kf::TrackKalmanFilter::Extrapolate
void Extrapolate(DataT z, const kf::FieldRegion< DataT > &F)
Definition KfTrackKalmanFilter.cxx:601
cbm::algo::kf::TrackKalmanFilter::Tr
kf::TrackParam< DataT > & Tr()
Definition KfTrackKalmanFilter.h:70
cbm::algo::kf::TrackKalmanFilter::MultipleScattering
void MultipleScattering(DataT radThick, DataT tx0, DataT ty0, DataT qp0)
apply multiple scattering correction to the track with the given Qp0
Definition KfTrackKalmanFilter.cxx:980
cbm::algo::kf::TrackKalmanFilter::EnergyLossCorrection
void EnergyLossCorrection(DataT radThick, FitDirection direction)
Definition KfTrackKalmanFilter.cxx:1057
cbm::algo::kf::TrackKalmanFilter::SetParticleMass
void SetParticleMass(DataT mass)
set particle mass for the fit
Definition KfTrackKalmanFilter.h:84
cbm::algo::kf::TrackKalmanFilter::SetTrack
void SetTrack(const kf::TrackParam< T > &t)
Definition KfTrackKalmanFilter.h:62
cbm::algo::kf::TrackParamBase::C00
T C00() const
Individual getters for covariance matrix elements.
Definition KfTrackParam.h:143
cbm::algo::kf::TrackParamBase::GetVi
T GetVi() const
Gets inverse velocity [ns/cm] in downstream direction.
Definition KfTrackParam.h:243
cbm::algo::kf::TrackParamBase::GetTy
T GetTy() const
Gets slope along y-axis.
Definition KfTrackParam.h:225
cbm::algo::kf::TrackParamBase::GetZ
T GetZ() const
Gets z position [cm].
Definition KfTrackParam.h:204
cbm::algo::kf::TrackParamBase::GetP
T GetP() const
Gets momentum [GeV/ec]. For the straight tracks returns 1.e4 [GeV/ec].
Definition KfTrackParam.h:285
cbm::algo::kf::TrackParamBase::GetPhi
T GetPhi() const
Gets azimuthal angle [rad].
Definition KfTrackParam.h:273
cbm::algo::kf::TrackParamBase::GetTime
T GetTime() const
Gets time [ns].
Definition KfTrackParam.h:237
cbm::algo::kf::TrackParamBase::GetTx
T GetTx() const
Gets slope along x-axis.
Definition KfTrackParam.h:219
cbm::algo::kf::TrackParamBase::GetQp
T GetQp() const
Gets charge over momentum [ec/GeV].
Definition KfTrackParam.h:231
cbm::algo::kf::TrackParamBase::GetNdf
T GetNdf() const
Gets NDF of track fit model.
Definition KfTrackParam.h:261
cbm::algo::kf::TrackParamBase::GetY
T GetY() const
Gets y position [cm].
Definition KfTrackParam.h:213
cbm::algo::kf::TrackParamBase::GetChiSq
T GetChiSq() const
Gets Chi-square of track fit model.
Definition KfTrackParam.h:258
cbm::algo::kf::TrackParamBase::ToString
std::string ToString(int i=-1) const
Prints parameters to a string.
Definition KfTrackParam.cxx:15
cbm::algo::kf::TrackParamBase::GetX
T GetX() const
Gets x position [cm].
Definition KfTrackParam.h:207
cbm::algo::kf::fmask::One
static fmask One()
Definition KfSimdPseudo.h:60
cbm::ca::tools::Debugger::Instance
static Debugger & Instance()
Instance.
Definition CaToolsDebugger.cxx:29
cbm::ca::tools::Debugger::AddNtuple
void AddNtuple(const char *name, const char *varlist)
Set new ntuple.
Definition CaToolsDebugger.cxx:67
cbm::ca::tools::Debugger::FillNtuple
void FillNtuple(const char *name, float v[])
Add an entry to ntuple.
Definition CaToolsDebugger.cxx:94
cbm::ca::tools::MCData::GetTrack
const auto & GetTrack(int idx) const
Gets a reference to MC track by its internal index.
Definition CaToolsMCData.h:149
cbm::ca::tools::MCData::GetTrackContainer
const auto & GetTrackContainer() const
Gets a reference to the vector of tracks.
Definition CaToolsMCData.h:155
cbm::ca::tools::MCData::GetPoint
const auto & GetPoint(int idx) const
Gets a reference to MC point by its index.
Definition CaToolsMCData.h:138
cbm::ca::tools::MCPoint
Class describes a unified MC-point, used in CA tracking QA analysis.
Definition CaToolsMCPoint.h:33
cbm::ca::tools::MCPoint::GetTyOut
double GetTyOut() const
Gets slope along x-axis at exit of station.
Definition CaToolsMCPoint.h:222
cbm::ca::tools::MCPoint::GetTime
double GetTime() const
Gets time [ns].
Definition CaToolsMCPoint.h:201
cbm::ca::tools::MCPoint::GetYOut
double GetYOut() const
Gets y coordinate at exit of station [cm].
Definition CaToolsMCPoint.h:240
cbm::ca::tools::MCPoint::GetInvSpeed
int GetInvSpeed() const
Gets inverse speed at reference z of station [ns/cm].
Definition CaToolsMCPoint.h:84
cbm::ca::tools::MCPoint::GetZOut
double GetZOut() const
Gets z coordinate at exit of station [cm].
Definition CaToolsMCPoint.h:249
cbm::ca::tools::MCPoint::GetQp
double GetQp() const
Gets track charge over momentum at reference z of station [ec/GeV].
Definition CaToolsMCPoint.h:180
cbm::ca::tools::MCPoint::GetXOut
double GetXOut() const
Gets x coordinate at exit of station [cm].
Definition CaToolsMCPoint.h:231
cbm::ca::tools::MCPoint::GetTxOut
double GetTxOut() const
Gets slope along x-axis at exit of station.
Definition CaToolsMCPoint.h:213
cbm::ca::tools::MCPoint::GetP
double GetP() const
Gets track momentum absolute value at reference z of station [GeV/c].
Definition CaToolsMCPoint.h:111
CbmL1Constants::MinFastMom
const double MinFastMom
Definition CbmL1Constants.h:12
cbm::algo::ca::constants::phys::SpeedOfLightD
constexpr double SpeedOfLightD
Speed of light [cm/ns].
Definition CaDefs.h:79
cbm::algo::ca::constants::phys::SpeedOfLightInvD
constexpr double SpeedOfLightInvD
Inverse speed of light [ns/cm].
Definition CaDefs.h:80
cbm::algo::ca
TODO: SZh 8.11.2022: add selection of parameterisation.
Definition CaBranch.h:14
cbm::algo::ca::fscal
kf::fscal fscal
Definition CaSimd.h:14
TL1Efficiencies::AddCounter
virtual void AddCounter(const TString &shortname, const TString &name)
Definition CbmL1Counters.h:167
TL1Efficiencies::indices
std::map< TString, int > indices
Definition CbmL1Counters.h:154
TL1Efficiencies::Inc
void Inc(bool isReco, const TString &name)
Definition CbmL1Counters.h:202
TL1Efficiencies::reco
TL1TracksCatCounters< int > reco
Definition CbmL1Counters.h:161
TL1Efficiencies::names
std::vector< TString > names
Definition CbmL1Counters.h:153
TL1Efficiencies::mc
TL1TracksCatCounters< int > mc
Definition CbmL1Counters.h:160
TL1Efficiencies::operator+=
TL1Efficiencies & operator+=(TL1Efficiencies &a)
Definition CbmL1Counters.h:191
TL1Efficiencies::ratio_reco
TL1TracksCatCounters< double > ratio_reco
Definition CbmL1Counters.h:156
TL1PerfEfficiencies::ratio_clone
TL1TracksCatCounters< double > ratio_clone
Definition CbmL1Performance.cxx:218
TL1PerfEfficiencies::mc_length_hits
TL1TracksCatCounters< int > mc_length_hits
Definition CbmL1Performance.cxx:227
TL1PerfEfficiencies::ratio_length
TL1TracksCatCounters< double > ratio_length
Definition CbmL1Performance.cxx:219
TL1PerfEfficiencies::ratio_fakes
TL1TracksCatCounters< double > ratio_fakes
Definition CbmL1Performance.cxx:220
TL1PerfEfficiencies::PrintEff
void PrintEff(bool ifPrintTableToLog=false, TDirectory *outDir=nullptr, const std::string &nameOfTable="efficiency_table")
Definition CbmL1Performance.cxx:162
TL1PerfEfficiencies::Inc
void Inc(bool isReco, bool isKilled, double _ratio_length, double _ratio_fakes, int _nclones, int _mc_length, int _mc_length_hits, const TString &name)
Definition CbmL1Performance.cxx:147
TL1PerfEfficiencies::reco_length
TL1TracksCatCounters< double > reco_length
Definition CbmL1Performance.cxx:224
TL1PerfEfficiencies::mc_length
TL1TracksCatCounters< int > mc_length
Definition CbmL1Performance.cxx:226
TL1PerfEfficiencies::reco_fakes
TL1TracksCatCounters< double > reco_fakes
Definition CbmL1Performance.cxx:225
TL1PerfEfficiencies::killed
TL1TracksCatCounters< int > killed
Definition CbmL1Performance.cxx:222
TL1PerfEfficiencies::ratio_killed
TL1TracksCatCounters< double > ratio_killed
Definition CbmL1Performance.cxx:217
TL1PerfEfficiencies::AddCounter
virtual void AddCounter(const TString &shortname, const TString &name)
Definition CbmL1Performance.cxx:110
TL1PerfEfficiencies::clone
TL1TracksCatCounters< int > clone
Definition CbmL1Performance.cxx:223
TL1PerfEfficiencies::operator+=
TL1PerfEfficiencies & operator+=(TL1PerfEfficiencies &a)
Definition CbmL1Performance.cxx:125
TL1TracksCatCounters
counters used for efficiency calculation
Definition CbmL1Counters.h:26
TL1TracksCatCounters::GetNcounters
int GetNcounters() const
Definition CbmL1Counters.h:31
TL1TracksCatCounters::counters
cacore::Vector< T > counters
Definition CbmL1Counters.h:122