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 cbm::algo::ca::McHitInfo& 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 const auto& rActSetup = fpAlgo->GetParameters().GetActiveSetup();
727
728 static int NEvents = 0;
729 if (NEvents > 0) {
730 h_reg_MCmom->Scale(NEvents);
731 h_acc_MCmom->Scale(NEvents);
732 h_reco_MCmom->Scale(NEvents);
733 h_ghost_Rmom->Scale(NEvents);
734 h_reg_prim_MCmom->Scale(NEvents);
735 h_acc_prim_MCmom->Scale(NEvents);
736 h_reco_prim_MCmom->Scale(NEvents);
737 h_reg_sec_MCmom->Scale(NEvents);
738 h_acc_sec_MCmom->Scale(NEvents);
739 h_reco_sec_MCmom->Scale(NEvents);
740 }
741 NEvents++;
742
743 // //hit density calculation: h_hit_density[10]
744 //
745 // for (vector<CbmL1HitDebugInfo>::iterator hIt = fvHitDebugInfo.begin(); hIt != fvHitDebugInfo.end(); ++hIt){
746 // float x = hIt->x;
747 // float y = hIt->y;
748 // float r = sqrt(x*x+y*y);
749 // h_hit_density[hIt->iStation]->Fill(r, 1.0/(2.0*3.1415*r));
750 // }
751
752 //
753 for (vector<CbmL1Track>::iterator rtraIt = fvRecoTracks.begin(); rtraIt != fvRecoTracks.end(); ++rtraIt) {
754 CbmL1Track* prtra = &(*rtraIt);
755 if ((prtra->Hits).size() < 1) continue;
756 { // fill histos
757 if (fabs(prtra->GetQp()) > 1.e-10)
758 h_reco_mom->Fill(
759 fabs(1.0 / prtra->GetQp())); // TODO: Is it a right precision? In FairTrackParam it is 1.e-4 (S.Zharko)
760 // NOTE: p = (TMath::Abs(fQp) > 1.e-4) ? 1. / TMath::Abs(fQp) : 1.e4; // FairTrackParam::Momentum(TVector3)
761 // 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
762
763 h_reco_phi->Fill(prtra->GetPhi()); // TODO: What is precision?
764 h_reco_nhits->Fill((prtra->Hits).size());
765 cbm::algo::ca::McHitInfo& mh = fvHitDebugInfo[prtra->Hits[0]];
766 h_reco_station->Fill(mh.iStation);
767
768 int iFstHit = prtra->GetFirstHitIndex();
769 auto& fstHit = fvHitDebugInfo[iFstHit];
770 h2_fst_hit_yz->Fill(rActSetup.GetActiveLayer(fstHit.iStation).GetZref()[0], fstHit.GetY());
771
772 int iLstHit = prtra->GetLastHitIndex();
773 auto& lstHit = fvHitDebugInfo[iLstHit];
774 h2_lst_hit_yz->Fill(rActSetup.GetActiveLayer(lstHit.iStation).GetZref()[0], lstHit.GetY());
775
776 for (int iH : prtra->Hits) {
777 const auto& hit = fvHitDebugInfo[iH];
778 int y = hit.GetY();
779 int z = rActSetup.GetActiveLayer(hit.iStation).GetZref()[0];
780 h2_all_hit_yz->Fill(z, y);
781 }
782 }
783
784
785 h_reco_purity->Fill(100 * prtra->GetMaxPurity());
786
787 if (prtra->GetNdf() > 0) {
788 if (prtra->IsGhost()) {
789 h_ghost_chi2->Fill(prtra->GetChiSq() / prtra->GetNdf());
790 h_ghost_prob->Fill(TMath::Prob(prtra->GetChiSq(), prtra->GetNdf()));
791 }
792 else if (prtra->GetNofMCTracks() > 0) {
793 const auto& mcTrk = fMCData.GetTrack(prtra->GetMatchedMCTrackIndex());
794 if (mcTrk.IsReconstructable()) {
795 h_reco_chi2->Fill(prtra->GetChiSq() / prtra->GetNdf());
796 h_reco_prob->Fill(TMath::Prob(prtra->GetChiSq(), prtra->GetNdf()));
797 }
798 else {
799 // h_rest_chi2->Fill(prtra->GetChiSq()/prtra->NDF);
800 h_rest_prob->Fill(TMath::Prob(prtra->GetChiSq(), prtra->GetNdf()));
801 }
802 }
803 }
804
805
806 // fill ghost histos
807 if (prtra->IsGhost()) {
808 //fMonitor.IncrementCounter(EMonitorKey::kGhostTrack);
809 h_ghost_purity->Fill(100 * prtra->GetMaxPurity());
810 if (fabs(prtra->GetQp()) > 1.e-10) {
811 h_ghost_mom->Fill(prtra->GetP());
812 h_ghost_phi->Fill(prtra->GetPhi()); // phi = atan(py / px) = atan(ty / tx)
813 h_ghost_Rmom->Fill(prtra->GetP());
814 }
815 h_ghost_nhits->Fill((prtra->Hits).size());
816 cbm::algo::ca::McHitInfo& h1 = fvHitDebugInfo[prtra->Hits[0]];
817 cbm::algo::ca::McHitInfo& h2 = fvHitDebugInfo[prtra->Hits[1]];
818 h_ghost_fstation->Fill(h1.iStation);
819 h_ghost_r->Fill(sqrt(fabs(h1.x * h1.x + h1.y * h1.y)));
820 double z1 = rActSetup.GetActiveLayer(h1.iStation).GetZref()[0];
821 double z2 = rActSetup.GetActiveLayer(h2.iStation).GetZref()[0];
822 if (fabs(z2 - z1) > 1.e-4) {
823 h_ghost_tx->Fill((h2.x - h1.x) / (z2 - z1));
824 h_ghost_ty->Fill((h2.y - h1.y) / (z2 - z1));
825 }
826
827 cbm::algo::ca::McHitInfo& hf = fvHitDebugInfo[prtra->Hits[0]];
828 cbm::algo::ca::McHitInfo& hl = fvHitDebugInfo[prtra->Hits[(prtra->Hits).size() - 1]];
829 if (fabs(prtra->GetQp()) > 1.e-10) {
830 h2_ghost_nhits_vs_mom->Fill(prtra->GetP(), prtra->Hits.size());
831 h2_ghost_fstation_vs_mom->Fill(prtra->GetP(), hf.iStation + 1);
832 }
833 if (hl.iStation >= hf.iStation) h2_ghost_lstation_vs_fstation->Fill(hf.iStation + 1, hl.iStation + 1);
834 }
835
836 } // for reco tracks
837
838 int mc_total = 0; // total amount of reconstructable mcTracks
839 for (const auto& mcTrk : fMCData.GetTrackContainer()) {
840 // if( !( mcTrk.GetPdgCode() == -11 && mcTrk.GetPdgCode() == -1 ) ) continue; // electrons only
841
842 // No Sts hits?
843 int nmchits = mcTrk.GetNofHits();
844 if (nmchits == 0) continue;
845
846 double momentum = mcTrk.GetP();
847 double pt = mcTrk.GetPt();
848 double theta = mcTrk.GetTheta() * 180 / 3.1415;
849 double phi = mcTrk.GetPhi();
850
851 h_mcp->Fill(momentum);
852 h_nmchits->Fill(nmchits);
853
854 int nSta = mcTrk.GetTotNofStationsWithHit();
855
856 cbm::algo::ca::McHitInfo& fh = fvHitDebugInfo[*(mcTrk.GetHitIndexes().begin())];
857 cbm::algo::ca::McHitInfo& lh = fvHitDebugInfo[*(mcTrk.GetHitIndexes().rbegin())];
858
859 h_reg_MCmom->Fill(momentum);
860 if (mcTrk.IsPrimary()) {
861 h_reg_mom_prim->Fill(momentum);
862 h_reg_phi_prim->Fill(phi);
863 h_reg_prim_MCmom->Fill(momentum);
864 h_reg_nhits_prim->Fill(nSta);
865 h2_reg_nhits_vs_mom_prim->Fill(momentum, nSta);
866 h2_reg_fstation_vs_mom_prim->Fill(momentum, fh.iStation + 1);
867 if (lh.iStation >= fh.iStation) h2_reg_lstation_vs_fstation_prim->Fill(fh.iStation + 1, lh.iStation + 1);
868 }
869 else {
870 h_reg_mom_sec->Fill(momentum);
871 h_reg_phi_sec->Fill(phi);
872 h_reg_sec_MCmom->Fill(momentum);
873 h_reg_nhits_sec->Fill(nSta);
874 if (momentum > 0.01) {
875 h2_reg_nhits_vs_mom_sec->Fill(momentum, nSta);
876 h2_reg_fstation_vs_mom_sec->Fill(momentum, fh.iStation + 1);
877 if (lh.iStation >= fh.iStation) h2_reg_lstation_vs_fstation_sec->Fill(fh.iStation + 1, lh.iStation + 1);
878 }
879 }
880
881 if (mcTrk.IsAdditional()) {
882 h_acc_mom_short123s->Fill(momentum);
883 }
884
885 if (!mcTrk.IsReconstructable()) continue;
886 mc_total++;
887
888 h_acc_MCmom->Fill(momentum);
889 if (mcTrk.IsPrimary()) {
890 h_acc_mom_prim->Fill(momentum);
891 h_acc_phi_prim->Fill(phi);
892 h_acc_prim_MCmom->Fill(momentum);
893 h_acc_nhits_prim->Fill(nSta);
894 h2_acc_nhits_vs_mom_prim->Fill(momentum, nSta);
895 h2_acc_fstation_vs_mom_prim->Fill(momentum, fh.iStation + 1);
896 if (lh.iStation >= fh.iStation) h2_acc_lstation_vs_fstation_prim->Fill(fh.iStation + 1, lh.iStation + 1);
897 }
898 else {
899 h_acc_mom_sec->Fill(momentum);
900 h_acc_phi_sec->Fill(phi);
901 h_acc_sec_MCmom->Fill(momentum);
902 h_acc_nhits_sec->Fill(nSta);
903 if (momentum > 0.01) {
904 h2_acc_nhits_vs_mom_sec->Fill(momentum, nSta);
905 h2_acc_fstation_vs_mom_sec->Fill(momentum, fh.iStation + 1);
906 if (lh.iStation >= fh.iStation) h2_acc_lstation_vs_fstation_sec->Fill(fh.iStation + 1, lh.iStation + 1);
907 }
908 }
909
910
911 // vertex distribution of primary and secondary tracks
912 h2_vertex->Fill(mcTrk.GetStartZ(), mcTrk.GetStartY());
913 //h3_vertex->Fill(mcTrk.GetStartZ(), mcTrk.GetStartX(), mcTrk.GetStartY());
914 if (mcTrk.GetMotherId() < 0) { // primary
915 h2_prim_vertex->Fill(mcTrk.GetStartZ(), mcTrk.GetStartY());
916 //h3_prim_vertex->Fill(mcTrk.GetStartZ(), mcTrk.GetStartX(), mcTrk.GetStartY());
917 }
918 else {
919 h2_sec_vertex->Fill(mcTrk.GetStartZ(), mcTrk.GetStartY());
920 //h3_sec_vertex->Fill(mcTrk.GetStartZ(), mcTrk.GetStartX(), mcTrk.GetStartY());
921 }
922
923
924 int iph = mcTrk.GetHitIndexes()[0];
925 cbm::algo::ca::McHitInfo& ph = fvHitDebugInfo[iph];
926
927 h_sec_r->Fill(sqrt(fabs(ph.x * ph.x + ph.y * ph.y)));
928 if (fabs(mcTrk.GetPz()) > 1.e-8) {
929 h_tx->Fill(mcTrk.GetTx());
930 h_ty->Fill(mcTrk.GetTy());
931 }
932
933 // reconstructed tracks
934 bool reco = (mcTrk.IsReconstructed());
935 int nMCHits = mcTrk.GetTotNofStationsWithHit();
936
937 if (reco) {
938 p_eff_all_vs_mom->Fill(momentum, 100.0);
939 p_eff_all_vs_nhits->Fill(nMCHits, 100.0);
940 p_eff_all_vs_pt->Fill(pt, 100.0);
941 h_reco_MCmom->Fill(momentum);
942 if (mcTrk.GetMotherId() < 0) { // primary
943 p_eff_prim_vs_mom->Fill(momentum, 100.0);
944 p_eff_prim_vs_nhits->Fill(nMCHits, 100.0);
945 p_eff_prim_vs_pt->Fill(pt, 100.0);
946 p_eff_prim_vs_theta->Fill(theta, 100.0);
947 }
948 else {
949 p_eff_sec_vs_mom->Fill(momentum, 100.0);
950 p_eff_sec_vs_nhits->Fill(nMCHits, 100.0);
951 }
952 if (mcTrk.GetMotherId() < 0) { // primary
953 h_reco_mom_prim->Fill(momentum);
954 h_reco_phi_prim->Fill(phi);
955 h_reco_prim_MCmom->Fill(momentum);
956 h_reco_nhits_prim->Fill(nSta);
957 h2_reco_nhits_vs_mom_prim->Fill(momentum, nSta);
958 h2_reco_fstation_vs_mom_prim->Fill(momentum, fh.iStation + 1);
959 if (lh.iStation >= fh.iStation) h2_reco_lstation_vs_fstation_prim->Fill(fh.iStation + 1, lh.iStation + 1);
960 }
961 else {
962 h_reco_mom_sec->Fill(momentum);
963 h_reco_phi_sec->Fill(phi);
964 h_reco_sec_MCmom->Fill(momentum);
965 h_reco_nhits_sec->Fill(nSta);
966 if (momentum > 0.01) {
967 h2_reco_nhits_vs_mom_sec->Fill(momentum, nSta);
968 h2_reco_fstation_vs_mom_sec->Fill(momentum, fh.iStation + 1);
969 if (lh.iStation >= fh.iStation) h2_reco_lstation_vs_fstation_sec->Fill(fh.iStation + 1, lh.iStation + 1);
970 }
971 }
972 }
973 else {
974 h_notfound_mom->Fill(momentum);
975 h_notfound_phi->Fill(phi);
976 p_eff_all_vs_mom->Fill(momentum, 0.0);
977 p_eff_all_vs_nhits->Fill(nMCHits, 0.0);
978 p_eff_all_vs_pt->Fill(pt, 0.0);
979 if (mcTrk.GetMotherId() < 0) { // primary
980 p_eff_prim_vs_mom->Fill(momentum, 0.0);
981 p_eff_prim_vs_nhits->Fill(nMCHits, 0.0);
982 p_eff_prim_vs_pt->Fill(pt, 0.0);
983 p_eff_prim_vs_theta->Fill(theta, 0.0);
984 }
985 else {
986 p_eff_sec_vs_mom->Fill(momentum, 0.0);
987 p_eff_sec_vs_nhits->Fill(nMCHits, 0.0);
988 }
989 if (mcTrk.GetMotherId() < 0) { // primary
990 h_rest_mom_prim->Fill(momentum);
991 h_rest_phi_prim->Fill(phi);
992 h_rest_nhits_prim->Fill(nSta);
993 h2_rest_nhits_vs_mom_prim->Fill(momentum, nSta);
994 h2_rest_fstation_vs_mom_prim->Fill(momentum, fh.iStation + 1);
995 if (lh.iStation >= fh.iStation) h2_rest_lstation_vs_fstation_prim->Fill(fh.iStation + 1, lh.iStation + 1);
996 }
997 else {
998 h_rest_mom_sec->Fill(momentum);
999 h_rest_phi_sec->Fill(phi);
1000 h_rest_nhits_sec->Fill(nSta);
1001 if (momentum > 0.01) {
1002 h2_rest_nhits_vs_mom_sec->Fill(momentum, nSta);
1003 h2_rest_fstation_vs_mom_sec->Fill(momentum, fh.iStation + 1);
1004 if (lh.iStation >= fh.iStation) h2_rest_lstation_vs_fstation_sec->Fill(fh.iStation + 1, lh.iStation + 1);
1005 }
1006 }
1007 }
1008
1009 // killed tracks. At least one hit of they belong to some recoTrack
1010 // bool killed = 0;
1011 if (!reco) {
1012 h_notfound_nhits->Fill(nmchits);
1013 h_notfound_station->Fill(ph.iStation);
1014 h_notfound_r->Fill(sqrt(fabs(ph.x * ph.x + ph.y * ph.y)));
1015
1016 if (mcTrk.GetNofPoints() > 0) {
1017 const auto& pMC = fMCData.GetPoint(mcTrk.GetPointIndexes()[0]);
1018 h_notfound_tx->Fill(pMC.GetTx());
1019 h_notfound_ty->Fill(pMC.GetTy());
1020 }
1021
1022 // CbmL1HitDebugInfo &ph21 = fvHitDebugInfo[mtra.Hits[0]];
1023 // CbmL1HitDebugInfo &ph22 = fvHitDebugInfo[mtra.Hits[1]];
1024
1025 // double z21 = rActSetup.GetActiveLayer(ph21.iStation).GetZref()[0];
1026 // double z22 = rActSetup.GetActiveLayer(ph22.iStation).GetZref()[0];
1027 // if( fabs(z22-z21)>1.e-4 ){
1028 // h_notfound_tx->Fill((ph22.x-ph21.x)/(z22-z21));
1029 // h_notfound_ty->Fill((ph22.y-ph21.y)/(z22-z21));
1030 // }
1031
1032 // if( mtra.IsDisturbed() ) killed = 1;
1033 }
1034
1035
1036 if (mcTrk.IsSignal()) { // D0
1037 h_reco_d0_mom->Fill(momentum);
1038 if (reco)
1039 p_eff_d0_vs_mom->Fill(momentum, 100.0);
1040 else
1041 p_eff_d0_vs_mom->Fill(momentum, 0.0);
1042 }
1043
1044 } // for mcTracks
1045
1046 int NFakes = 0;
1047 for (unsigned int ih = 0; ih < fpAlgo->GetInputData().GetNhits(); ih++) {
1048 int iMC = fvHitDebugInfo[ih].GetBestMcPointId(); // TODO2: adapt to linking
1049 if (iMC < 0) NFakes++;
1050 }
1051
1052 h_reco_time->Fill(fTrackingTime);
1053 h_reco_timeNtr->Fill(mc_total, fTrackingTime);
1054 h_reco_timeNhit->Fill(fpAlgo->GetInputData().GetNhits(), fTrackingTime);
1055
1056 h_reco_fakeNtr->Fill(mc_total, NFakes);
1057 h_reco_fakeNhit->Fill(fpAlgo->GetInputData().GetNhits() - NFakes, NFakes);
1058
1059
1060 // NOTE: Must be called once
1061 h_reg_MCmom->Scale(1.f / NEvents);
1062 h_acc_MCmom->Scale(1.f / NEvents);
1063 h_reco_MCmom->Scale(1.f / NEvents);
1064 h_ghost_Rmom->Scale(1.f / NEvents);
1065 h_reg_prim_MCmom->Scale(1.f / NEvents);
1066 h_acc_prim_MCmom->Scale(1.f / NEvents);
1067 h_reco_prim_MCmom->Scale(1.f / NEvents);
1068 h_reg_sec_MCmom->Scale(1.f / NEvents);
1069 h_acc_sec_MCmom->Scale(1.f / NEvents);
1070 h_reco_sec_MCmom->Scale(1.f / NEvents);
1071
1072} // void CbmL1::HistoPerformance()
1073
1074
1075void CbmL1::TrackFitPerformance()
1076{
1077 const auto& rActSetup = fpAlgo->GetParameters().GetActiveSetup();
1078
1079 const int Nh_fit = 17;
1080 static TH1F *h_fit[Nh_fit], *h_fitL[Nh_fit], *h_fitSV[Nh_fit], *h_fitPV[Nh_fit], *h_fit_chi2;
1081
1082 static TH2F *PRes2D, *PRes2DPrim, *PRes2DSec;
1083
1084 static TH2F* pion_res_pt_fstDca;
1085 static TH2F* kaon_res_pt_fstDca;
1086 static TH2F* pton_res_pt_fstDca;
1087
1088 static TH2F* pion_res_pt_vtxDca;
1089 static TH2F* kaon_res_pt_vtxDca;
1090 static TH2F* pton_res_pt_vtxDca;
1091
1092 static TH2F* pion_res_p_fstDca;
1093 static TH2F* kaon_res_p_fstDca;
1094 static TH2F* pton_res_p_fstDca;
1095
1096 static TH2F* pion_res_p_vtxDca;
1097 static TH2F* kaon_res_p_vtxDca;
1098 static TH2F* pton_res_p_vtxDca;
1099
1100 static bool first_call = 1;
1101
1102 kf::TrackKalmanFilter<fvec> fit;
1103 fit.SetParticleMass(fpAlgo->GetDefaultParticleMass());
1104 fit.SetMask(fmask::One());
1105 //fit.SetMaxExtrapolationStep(10.);
1106
1107 if (first_call) {
1108 first_call = 0;
1109
1110 TDirectory* currentDir = gDirectory;
1111 gDirectory = fHistoDir;
1112 gDirectory->cd("Fit");
1113 {
1114 PRes2D = new TH2F("PRes2D", "Resolution P [%] vs P", 100, 0., 20., 100, -15., 15.);
1115 PRes2DPrim = new TH2F("PRes2DPrim", "Resolution P [%] vs P", 100, 0., 20., 100, -15., 15.);
1116 PRes2DSec = new TH2F("PRes2DSec", "Resolution P [%] vs P", 100, 0., 20., 100, -15., 15.);
1117
1118 pion_res_pt_fstDca = new TH2F("pion_res_pt_fstDca", "", 100, 0, 10, 1000, -500, 500);
1119 kaon_res_pt_fstDca = new TH2F("kaon_res_pt_fstDca", "", 100, 0, 10, 1000, -500, 500);
1120 pton_res_pt_fstDca = new TH2F("pton_res_pt_fstDca", "", 100, 0, 10, 1000, -500, 500);
1121
1122 pion_res_pt_vtxDca = new TH2F("pion_res_pt_vtxDca", "", 100, 0, 10, 1000, -5000, 5000);
1123 kaon_res_pt_vtxDca = new TH2F("kaon_res_pt_vtxDca", "", 100, 0, 10, 1000, -5000, 5000);
1124 pton_res_pt_vtxDca = new TH2F("pton_res_pt_vtxDca", "", 100, 0, 10, 1000, -5000, 5000);
1125
1126 pion_res_p_fstDca = new TH2F("pion_res_p_fstDca", "", 100, 0, 10, 1000, -500, 500);
1127 kaon_res_p_fstDca = new TH2F("kaon_res_p_fstDca", "", 100, 0, 10, 1000, -500, 500);
1128 pton_res_p_fstDca = new TH2F("pton_res_p_fstDca", "", 100, 0, 10, 1000, -500, 500);
1129
1130 pion_res_p_vtxDca = new TH2F("pion_res_p_vtxDca", "", 100, 0, 10, 1000, -5000, 5000);
1131 kaon_res_p_vtxDca = new TH2F("kaon_res_p_vtxDca", "", 100, 0, 10, 1000, -5000, 5000);
1132 pton_res_p_vtxDca = new TH2F("pton_res_p_vtxDca", "", 100, 0, 10, 1000, -5000, 5000);
1133
1134 struct {
1135 const char* name;
1136 const char* title;
1137 Int_t n;
1138 Double_t l, r;
1139 } Table[Nh_fit] = {{"x", "Residual X [#mum]", 140, -40., 40.},
1140 {"y", "Residual Y [#mum]", 100, -450., 450.},
1141 {"tx", "Residual Tx [mrad]", 100, -4., 4.},
1142 {"ty", "Residual Ty [mrad]", 100, -3.5, 3.5},
1143 {"P", "Resolution P/Q [%]", 100, -15., 15.},
1144 {"px", "Pull X [residual/estimated_error]", 100, -6., 6.},
1145 {"py", "Pull Y [residual/estimated_error]", 100, -6., 6.},
1146 {"ptx", "Pull Tx [residual/estimated_error]", 100, -6., 6.},
1147 {"pty", "Pull Ty [residual/estimated_error]", 100, -6., 6.},
1148 {"pQP", "Pull Q/P [residual/estimated_error]", 100, -6., 6.},
1149 {"QPreco", "Reco Q/P ", 100, -10., 10.},
1150 {"QPmc", "MC Q/P ", 100, -10., 10.},
1151 {"time", "Residual Time [ns]", 100, -10., 10.},
1152 {"ptime", "Pull Time [residual/estimated_error]", 100, -10., 10.},
1153 {"vi", "Residual Vi [1/c]", 100, -4., 4.},
1154 {"pvi", "Pull Vi [residual/estimated_error]", 100, -10., 10.},
1155 {"distrVi", "Vi distribution [1/c]", 100, 0., 4.}};
1156
1157 if (ca::TrackingMode::kGlobal == fpAlgo->GetTrackingMode()) {
1158 Table[4].l = -100.;
1159 Table[4].r = 100.;
1160 }
1161
1162 struct Tab {
1163 const char* name;
1164 const char* title;
1165 Int_t n;
1166 Double_t l, r;
1167 };
1168 Tab TableVertex[Nh_fit] = {//{"x", "Residual X [cm]", 200, -0.01, 0.01},
1169 {"x", "Residual X [cm]", 2000, -1., 1.},
1170 //{"y", "Residual Y [cm]", 200, -0.01, 0.01},
1171 {"y", "Residual Y [cm]", 2000, -1., 1.},
1172 //{"tx", "Residual Tx [mrad]", 100, -3., 3.},
1173 {"tx", "Residual Tx [mrad]", 100, -2., 2.},
1174 //{"ty", "Residual Ty [mrad]", 100, -3., 3.},
1175 {"ty", "Residual Ty [mrad]", 100, -2., 2.},
1176 {"P", "Resolution P/Q [%]", 100, -15., 15.},
1177 {"px", "Pull X [residual/estimated_error]", 100, -10., 10.},
1178 {"py", "Pull Y [residual/estimated_error]", 100, -10., 10.},
1179 {"ptx", "Pull Tx [residual/estimated_error]", 100, -10., 10.},
1180 {"pty", "Pull Ty [residual/estimated_error]", 100, -10., 10.},
1181 {"pQP", "Pull Q/P [residual/estimated_error]", 100, -10., 10.},
1182 {"QPreco", "Reco Q/P ", 100, -10., 10.},
1183 {"QPmc", "MC Q/P ", 100, -10., 10.},
1184 {"time", "Residual Time [ns]", 100, -10., 10.},
1185 {"ptime", "Pull Time [residual/estimated_error]", 100, -10., 10.},
1186 {"vi", "Residual Vi [1/c]", 100, -10., 10.},
1187 {"pvi", "Pull Vi [residual/estimated_error]", 100, -10., 10.},
1188 {"distrVi", "Vi distribution [1/c]", 100, -10., 10.}};
1189
1190 if (ca::TrackingMode::kGlobal == fpAlgo->GetTrackingMode()) {
1191 TableVertex[4].l = -1.;
1192 TableVertex[4].r = 1.;
1193 }
1194
1195 for (int i = 0; i < Nh_fit; i++) {
1196 char n[225], t[255];
1197 sprintf(n, "fst_%s", Table[i].name);
1198 sprintf(t, "First point %s", Table[i].title);
1199 h_fit[i] = new TH1F(n, t, Table[i].n, Table[i].l, Table[i].r);
1200 sprintf(n, "lst_%s", Table[i].name);
1201 sprintf(t, "Last point %s", Table[i].title);
1202 h_fitL[i] = new TH1F(n, t, Table[i].n, Table[i].l, Table[i].r);
1203 sprintf(n, "svrt_%s", TableVertex[i].name);
1204 sprintf(t, "Secondary vertex point %s", TableVertex[i].title);
1205 h_fitSV[i] = new TH1F(n, t, TableVertex[i].n, TableVertex[i].l, TableVertex[i].r);
1206 sprintf(n, "pvrt_%s", TableVertex[i].name);
1207 sprintf(t, "Primary vertex point %s", TableVertex[i].title);
1208 h_fitPV[i] = new TH1F(n, t, TableVertex[i].n, TableVertex[i].l, TableVertex[i].r);
1209
1210 for (int j = 0; j < 12; j++) {
1211 sprintf(n, "smth_%s_sta_%i", Table[i].name, j);
1212 sprintf(t, "Station %i %s", j, Table[i].title);
1213 // h_smoothed[j][i] = new TH1F(n,t, Table[i].n, Table[i].l, Table[i].r);
1214 }
1215 }
1216 h_fit_chi2 = new TH1F("h_fit_chi2", "Chi2/NDF", 50, -0.5, 10.0);
1217 }
1218 gDirectory = currentDir;
1219 } // if first call
1220
1221
1222 for (vector<CbmL1Track>::iterator it = fvRecoTracks.begin(); it != fvRecoTracks.end(); ++it) {
1223
1224 if (it->IsGhost()) continue;
1225
1226 bool isTimeFitted = 0;
1227
1228 {
1229 int nTimeMeasurenments = 0;
1230 for (unsigned int ih = 0; ih < it->Hits.size(); ih++) {
1231 int ista = fvHitDebugInfo[it->Hits[ih]].iStation;
1232 if (rActSetup.GetActiveLayer(ista).IsTimeMeasured()) {
1233 nTimeMeasurenments++;
1234 }
1235 }
1236 isTimeFitted = (nTimeMeasurenments > 1);
1237 }
1238
1239 do { // first hit
1240
1241 if (it->GetNofMCTracks() < 1) {
1242 break;
1243 }
1244 const auto& mcTrk = fMCData.GetTrack(it->GetMCTrackIndexes()[0]);
1245 int imcPoint = fvHitDebugInfo[it->Hits.front()].GetBestMcPointId();
1246
1247 // extrapolate to the first mc point of the mc track !
1248
1249 imcPoint = mcTrk.GetPointIndexes()[0];
1250
1251 if (imcPoint < 0) break;
1252
1253 const auto& mcP = fMCData.GetPoint(imcPoint);
1254
1255 TrackParamV tr(*it);
1256 FillFitHistos(tr, mcP, isTimeFitted, h_fit);
1257
1258 double dx = tr.GetX()[0] - mcP.GetXIn();
1259 double dy = tr.GetY()[0] - mcP.GetYIn();
1260 double dca = sqrt(dx * dx + dy * dy);
1261 // make dca distance negative for the half of the tracks to ease the gaussian fit and the rms calculation
1262 if (mcTrk.GetId() % 2) dca = -dca;
1263 double pt = mcP.GetPt();
1264 double dP = (fabs(1. / tr.GetQp()[0]) - mcP.GetP()) / mcP.GetP();
1265
1266 PRes2D->Fill(mcP.GetP(), 100. * dP);
1267
1268 if (mcTrk.IsPrimary()) {
1269
1270 h_fit[4]->Fill(100. * dP);
1271
1272 PRes2DPrim->Fill(mcP.GetP(), 100. * dP);
1273
1274 if (abs(mcTrk.GetPdgCode()) == 211) {
1275 pion_res_p_fstDca->Fill(mcP.GetP(), dca * 1.e4);
1276 pion_res_pt_fstDca->Fill(pt, dca * 1.e4);
1277 }
1278 if (abs(mcTrk.GetPdgCode()) == 321) {
1279 kaon_res_p_fstDca->Fill(mcP.GetP(), dca * 1.e4);
1280 kaon_res_pt_fstDca->Fill(pt, dca * 1.e4);
1281 }
1282 if (abs(mcTrk.GetPdgCode()) == 2212) {
1283 pton_res_p_fstDca->Fill(mcP.GetP(), dca * 1.e4);
1284 pton_res_pt_fstDca->Fill(pt, dca * 1.e4);
1285 }
1286 }
1287 else {
1288 PRes2DSec->Fill(mcP.GetP(), 100. * dP);
1289 }
1290 } while (0);
1291
1292
1293 { // last hit
1294
1295 int iMC = fvHitDebugInfo[it->Hits.back()].GetBestMcPointId(); // TODO2: adapt to linking
1296 if (iMC < 0) continue;
1297 const auto& mcP = fMCData.GetPoint(iMC);
1298 TrackParamV tr(it->TLast);
1299 FillFitHistos(tr, mcP, isTimeFitted, h_fitL);
1300 }
1301
1302 do { // vertex
1303
1304 if (it->GetNofMCTracks() < 1) {
1305 break;
1306 }
1307
1308 const auto& mc = fMCData.GetTrack(it->GetMatchedMCTrackIndex());
1309 fit.SetTrack(*it);
1310
1311 const TrackParamV& tr = fit.Tr();
1312
1313 // if (mc.mother_ID != -1) { // secondary
1314 if (!mc.IsPrimary()) { // secondary
1315
1316 { // extrapolate to vertex
1317 kf::FieldRegion<fvec> fld(kf::GlobalField::fgOriginalFieldType, kf::GlobalField::fgOriginalField);
1318 fit.Extrapolate(mc.GetStartZ(), fld);
1319 // add material
1320 const int fSta = fvHitDebugInfo[it->Hits[0]].iStation;
1321 const int dir = int((mc.GetStartZ() - rActSetup.GetActiveLayer(fSta).GetZref()[0])
1322 / fabs(mc.GetStartZ() - rActSetup.GetActiveLayer(fSta).GetZref()[0]));
1323 // if (abs(mc.z - fpAlgo->GetParameters().GetStation(fSta).fZ[0]) > 10.) continue; // can't extrapolate on large distance
1324 for (int iSta = fSta /*+dir*/; (iSta >= 0) && (iSta < fNStations)
1325 && (dir * (mc.GetStartZ() - rActSetup.GetActiveLayer(iSta).GetZref()[0]) > 0);
1326 iSta += dir) {
1327 // LOG(info) << iSta << " " << dir;
1328 auto radThick = rActSetup.GetMaterial(iSta).GetThicknessX0(fit.Tr().GetX(), fit.Tr().GetY());
1329 fit.MultipleScattering(radThick);
1330 fit.EnergyLossCorrection(radThick, kf::FitDirection::kUpstream);
1331 }
1332 }
1333 if (mc.GetStartZ() != tr.GetZ()[0]) continue;
1334
1335 // static int good = 0;
1336 // static int bad = 0;
1337 // if (mc.z != tr.GetZ()[0]){
1338 // bad++;
1339 // continue;
1340 // }
1341 // else good++;
1342 // LOG(info) << "bad\\good" << bad << " " << good;
1343
1344 // calculate pulls
1345 //h_fitSV[0]->Fill( (mc.x-tr.GetX()[0]) *1.e4);
1346 //h_fitSV[1]->Fill( (mc.y-tr.GetY()[0]) *1.e4);
1347 h_fitSV[0]->Fill((tr.GetX()[0] - mc.GetStartX()));
1348 h_fitSV[1]->Fill((tr.GetY()[0] - mc.GetStartY()));
1349 h_fitSV[2]->Fill((tr.GetTx()[0] - mc.GetTx()) * 1.e3);
1350 h_fitSV[3]->Fill((tr.GetTy()[0] - mc.GetTy()) * 1.e3);
1351 h_fitSV[4]->Fill(100. * (fabs(1. / tr.GetQp()[0]) / mc.GetP() - 1.));
1352 if (std::isfinite((fscal) tr.C00()[0]) && tr.C00()[0] > 0) {
1353 h_fitSV[5]->Fill((tr.GetX()[0] - mc.GetStartX()) / sqrt(tr.C00()[0]));
1354 }
1355 if (std::isfinite((fscal) tr.C11()[0]) && tr.C11()[0] > 0)
1356 h_fitSV[6]->Fill((tr.GetY()[0] - mc.GetStartY()) / sqrt(tr.C11()[0]));
1357 if (std::isfinite((fscal) tr.C22()[0]) && tr.C22()[0] > 0)
1358 h_fitSV[7]->Fill((tr.GetTx()[0] - mc.GetTx()) / sqrt(tr.C22()[0]));
1359 if (std::isfinite((fscal) tr.C33()[0]) && tr.C33()[0] > 0)
1360 h_fitSV[8]->Fill((tr.GetTy()[0] - mc.GetTy()) / sqrt(tr.C33()[0]));
1361 if (std::isfinite((fscal) tr.C44()[0]) && tr.C44()[0] > 0)
1362 h_fitSV[9]->Fill((tr.GetQp()[0] - mc.GetCharge() / mc.GetP()) / sqrt(tr.C44()[0]));
1363 h_fitSV[10]->Fill(tr.GetQp()[0]);
1364 h_fitSV[11]->Fill(mc.GetCharge() / mc.GetP());
1365 if (isTimeFitted) {
1366 h_fitSV[12]->Fill(tr.GetTime()[0] - mc.GetStartT());
1367 if (std::isfinite((fscal) tr.C55()[0]) && tr.C55()[0] > 0.) {
1368 h_fitSV[13]->Fill((tr.GetTime()[0] - mc.GetStartT()) / sqrt(tr.C55()[0]));
1369 }
1370 double dvi =
1371 tr.GetVi()[0]
1372 - sqrt(1. + mc.GetMass() * mc.GetMass() / mc.GetP() / mc.GetP()) * constants::phys::SpeedOfLightInvD;
1373 h_fitSV[14]->Fill(dvi * constants::phys::SpeedOfLightD);
1374 if (std::isfinite((fscal) tr.C66()[0]) && tr.C66()[0] > 0.) {
1375 h_fitSV[15]->Fill(dvi / sqrt(tr.C66()[0]));
1376 }
1377 h_fitSV[16]->Fill(tr.GetVi()[0] * constants::phys::SpeedOfLightD);
1378 }
1379 }
1380 else { // primary
1381
1382 { // extrapolate to vertex
1383 kf::FieldRegion<fvec> fld(kf::GlobalField::fgOriginalFieldType, kf::GlobalField::fgOriginalField);
1384
1385 // add material
1386 const int fSta = fvHitDebugInfo[it->Hits[0]].iStation;
1387
1388 const int dir = (mc.GetStartZ() - rActSetup.GetActiveLayer(fSta).GetZref()[0])
1389 / abs(mc.GetStartZ() - rActSetup.GetActiveLayer(fSta).GetZref()[0]);
1390 // if (abs(mc.z - fpAlgo->GetParameters().GetStation(fSta].fZ[0]) > 10.) continue; // can't extrapolate on large distance
1391
1392 for (int iSta = fSta + dir; (iSta >= 0) && (iSta < fNStations)
1393 && (dir * (mc.GetStartZ() - rActSetup.GetActiveLayer(fSta).GetZref()[0]) > 0);
1394 iSta += dir) {
1395
1396 fit.Extrapolate(rActSetup.GetActiveLayer(iSta).GetZref(), fld);
1397 auto radThick = rActSetup.GetMaterial(iSta).GetThicknessX0(fit.Tr().GetX(), 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::algo::ca::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.SetTrack(track);
1530 kf::FieldRegion<fvec> fld(kf::EFieldType::Normal, kf::GlobalField::fgOriginalField);
1531 fit.Extrapolate(mcP.GetZOut(), fld);
1532 track = fit.Tr();
1533
1534 const TrackParamV& tr = track;
1535
1536 h[0]->Fill((tr.GetX()[0] - mcP.GetXOut()) * 1.e4);
1537 h[1]->Fill((tr.GetY()[0] - mcP.GetYOut()) * 1.e4);
1538 h[2]->Fill((tr.GetTx()[0] - mcP.GetTxOut()) * 1.e3);
1539 h[3]->Fill((tr.GetTy()[0] - mcP.GetTyOut()) * 1.e3);
1540 h[4]->Fill(100. * (fabs(1. / tr.GetQp()[0]) / mcP.GetP() - 1.));
1541
1542 if (std::isfinite((fscal) tr.C00()[0]) && tr.C00()[0] > 0)
1543 h[5]->Fill((tr.GetX()[0] - mcP.GetXOut()) / sqrt(tr.C00()[0]));
1544 if (std::isfinite((fscal) tr.C11()[0]) && tr.C11()[0] > 0)
1545 h[6]->Fill((tr.GetY()[0] - mcP.GetYOut()) / sqrt(tr.C11()[0]));
1546 if (std::isfinite((fscal) tr.C22()[0]) && tr.C22()[0] > 0)
1547 h[7]->Fill((tr.GetTx()[0] - mcP.GetTxOut()) / sqrt(tr.C22()[0]));
1548 if (std::isfinite((fscal) tr.C33()[0]) && tr.C33()[0] > 0)
1549 h[8]->Fill((tr.GetTy()[0] - mcP.GetTyOut()) / sqrt(tr.C33()[0]));
1550 if (std::isfinite((fscal) tr.C44()[0]) && tr.C44()[0] > 0)
1551 h[9]->Fill((tr.GetQp()[0] - mcP.GetQp()) / sqrt(tr.C44()[0]));
1552 h[10]->Fill(tr.GetQp()[0]);
1553 h[11]->Fill(mcP.GetQp());
1554 if (isTimeFitted) {
1555 h[12]->Fill(tr.GetTime()[0] - mcP.GetTime());
1556 if (std::isfinite((fscal) tr.C55()[0]) && tr.C55()[0] > 0.) {
1557 h[13]->Fill((tr.GetTime()[0] - mcP.GetTime()) / sqrt(tr.C55()[0]));
1558 }
1559
1560 double dvi = tr.GetVi()[0] - mcP.GetInvSpeed();
1561 h[14]->Fill(dvi * constants::phys::SpeedOfLightD);
1562 if (std::isfinite((fscal) tr.C66()[0]) && tr.C66()[0] > 0.) {
1563 h[15]->Fill(dvi / sqrt(tr.C66()[0]));
1564 }
1565 h[16]->Fill(tr.GetVi()[0] * constants::phys::SpeedOfLightD);
1566 }
1567}
1568
1569
1570void CbmL1::FieldApproxCheck()
1571{
1572 TDirectory* curr = gDirectory;
1573 TFile* currentFile = gFile;
1574 TFile* fout = new TFile("CaFieldApproxQa.root", "RECREATE");
1575 fout->cd();
1576
1577 const auto& rActSetup = fpAlgo->GetParameters().GetActiveSetup();
1578
1579 assert(FairRunAna::Instance());
1580 FairField* MF = FairRunAna::Instance()->GetField();
1581 assert(MF);
1582
1583 for (int ist = 0; ist < rActSetup.GetNofLayers(); ist++) {
1584
1585 const auto& st = rActSetup.GetActiveLayer(ist);
1586
1587 double z = st.GetZref()[0];
1588 double Xmax = st.GetXmax()[0];
1589 double Ymax = st.GetYmax()[0];
1590
1591 int NbinsX = 101;
1592 int NbinsY = 101;
1593
1594 TH2F* hdB[4] = {
1595 new TH2F(Form("station_%i_dB", ist), Form("station %i, dB, z = %0.f cm", ist, z), //
1596 NbinsX, -Xmax, Xmax, NbinsY, -Ymax, Ymax),
1597 new TH2F(Form("station_%i_dBx", ist), Form("station %i, dBx, z = %0.f cm", ist, z), //
1598 NbinsX, -Xmax, Xmax, NbinsY, -Ymax, Ymax),
1599 new TH2F(Form("station_%i_dBy", ist), Form("station %i, dBy, z = %0.f cm", ist, z), //
1600 NbinsX, -Xmax, Xmax, NbinsY, -Ymax, Ymax),
1601 new TH2F(Form("station_%i_dBz", ist), Form("station %i, dBz, z = %0.f cm", ist, z), //
1602 NbinsX, -Xmax, Xmax, NbinsY, -Ymax, Ymax) //
1603 };
1604
1605 TH2F* hB[4] = {
1606 new TH2F(Form("station_%i_B", ist), Form("station %i, B, z = %0.f cm", ist, z), //
1607 NbinsX, -Xmax, Xmax, NbinsY, -Ymax, Ymax),
1608 new TH2F(Form("station_%i_Bx", ist), Form("station %i, Bx, z = %0.f cm", ist, z), //
1609 NbinsX, -Xmax, Xmax, NbinsY, -Ymax, Ymax),
1610 new TH2F(Form("station_%i_By", ist), Form("station %i, By, z = %0.f cm", ist, z), //
1611 NbinsX, -Xmax, Xmax, NbinsY, -Ymax, Ymax),
1612 new TH2F(Form("station_%i_Bz", ist), Form("station %i, Bz, z = %0.f cm", ist, z), //
1613 NbinsX, -Xmax, Xmax, NbinsY, -Ymax, Ymax) //
1614 };
1615
1616 for (int i = 0; i < 4; i++) {
1617 hdB[i]->GetXaxis()->SetTitle("X, cm");
1618 hdB[i]->GetYaxis()->SetTitle("Y, cm");
1619 hdB[i]->GetXaxis()->SetTitleOffset(1);
1620 hdB[i]->GetYaxis()->SetTitleOffset(1);
1621 hdB[i]->GetZaxis()->SetTitleOffset(1.3);
1622 hB[i]->GetXaxis()->SetTitle("X, cm");
1623 hB[i]->GetYaxis()->SetTitle("Y, cm");
1624 hB[i]->GetXaxis()->SetTitleOffset(1);
1625 hB[i]->GetYaxis()->SetTitleOffset(1);
1626 hB[i]->GetZaxis()->SetTitleOffset(1.3);
1627 }
1628
1629 hdB[0]->GetZaxis()->SetTitle("B_map - B_L1, kGauss");
1630 hdB[1]->GetZaxis()->SetTitle("Bx_map - Bx_L1, kGauss");
1631 hdB[2]->GetZaxis()->SetTitle("By_map - By_L1, kGauss");
1632 hdB[3]->GetZaxis()->SetTitle("Bz_map - Bz_L1, kGauss");
1633
1634 hdB[0]->GetZaxis()->SetTitle("B_map, kGauss");
1635 hdB[1]->GetZaxis()->SetTitle("Bx_map, kGauss");
1636 hdB[2]->GetZaxis()->SetTitle("By_map, kGauss");
1637 hdB[3]->GetZaxis()->SetTitle("Bz_map, kGauss");
1638
1639
1640 for (int i = 1; i <= hdB[0]->GetXaxis()->GetNbins(); i++) {
1641 double x = hdB[0]->GetXaxis()->GetBinCenter(i);
1642 for (int j = 1; j <= hdB[0]->GetYaxis()->GetNbins(); j++) {
1643 double y = hdB[0]->GetYaxis()->GetBinCenter(j);
1644 double r[] = {x, y, z};
1645 double B[] = {0., 0., 0.};
1646 MF->GetFieldValue(r, B);
1647
1648 double bbb = sqrt(B[0] * B[0] + B[1] * B[1] + B[2] * B[2]);
1649
1650 auto B_L1 = rActSetup.GetField().GetFieldValue(ist, x, y);
1651
1652 hdB[0]->SetBinContent(i, j, (bbb - B_L1.GetAbs()[0]));
1653 hdB[1]->SetBinContent(i, j, (B[0] - B_L1.GetBx()[0]));
1654 hdB[2]->SetBinContent(i, j, (B[1] - B_L1.GetBy()[0]));
1655 hdB[3]->SetBinContent(i, j, (B[2] - B_L1.GetBz()[0]));
1656
1657 hB[0]->SetBinContent(i, j, bbb);
1658 hB[1]->SetBinContent(i, j, B[0]);
1659 hB[2]->SetBinContent(i, j, B[1]);
1660 hB[3]->SetBinContent(i, j, B[2]);
1661 }
1662 }
1663
1664 for (int i = 0; i < 4; i++) {
1665 hdB[i]->Write();
1666 hB[i]->Write();
1667 }
1668
1669 } // for ista
1670
1671 fout->Close();
1672 fout->Delete();
1673 gFile = currentFile;
1674 gDirectory = curr;
1675} // void CbmL1::FieldApproxCheck()
1676
1677
1678#include "TMath.h" // ????????????????????? FIXME
1679void CbmL1::FieldIntegralCheck()
1680{
1681 TDirectory* curr = gDirectory;
1682 TFile* currentFile = gFile;
1683 TFile* fout = new TFile("FieldApprox.root", "RECREATE");
1684 fout->cd();
1685
1686 assert(FairRunAna::Instance());
1687 FairField* MF = FairRunAna::Instance()->GetField();
1688 assert(MF);
1689
1690 int nPointsZ = 1000;
1691 int nPointsPhi = 100;
1692 int nPointsTheta = 100;
1693 double startZ = 0, endZ = 100.;
1694 double startPhi = 0, endPhi = 2 * TMath::Pi();
1695 double startTheta = -30. / 180. * TMath::Pi(), endTheta = 30. / 180. * TMath::Pi();
1696
1697 double DZ = endZ - startZ;
1698 double DP = endPhi - startPhi;
1699 double DT = endTheta - startTheta;
1700
1701 float ddp = endPhi / nPointsPhi;
1702 float ddt = 2 * endTheta / nPointsTheta;
1703
1704 TH2F* hSb =
1705 new TH2F("Field Integral", "Field Integral", static_cast<int>(nPointsPhi), -(startPhi + ddp / 2.),
1706 (endPhi + ddp / 2.), static_cast<int>(nPointsTheta), (startTheta - ddt / 2.), (endTheta + ddt / 2.));
1707
1708 for (int iP = 0; iP < nPointsPhi; iP++) {
1709 double phi = startPhi + iP * DP / nPointsPhi;
1710 for (int iT = 0; iT < nPointsTheta; iT++) {
1711 double theta = startTheta + iT * DT / nPointsTheta;
1712
1713 double Sb = 0;
1714 for (int iZ = 1; iZ < nPointsZ; iZ++) {
1715 double z = startZ + DZ * iZ / nPointsZ;
1716 double x = z * TMath::Tan(theta) * TMath::Cos(phi);
1717 double y = z * TMath::Tan(theta) * TMath::Sin(phi);
1718 double r[3] = {x, y, z};
1719 double b[3];
1720 MF->GetFieldValue(r, b);
1721 double B = sqrt(b[0] * b[0] + b[1] * b[1] + b[2] * b[2]);
1722 Sb += B * DZ / nPointsZ / 100. / 10.;
1723 }
1724 hSb->SetBinContent(iP + 1, iT + 1, Sb);
1725 }
1726 }
1727
1728 hSb->GetXaxis()->SetTitle("#phi [rad]");
1729 hSb->GetYaxis()->SetTitle("#theta [rad]");
1730 hSb->GetXaxis()->SetTitleOffset(1);
1731 hSb->GetYaxis()->SetTitleOffset(1);
1732 hSb->GetZaxis()->SetTitle("Field Integral [T#dotm]");
1733 hSb->GetZaxis()->SetTitleOffset(1.3);
1734
1735 hSb->Write();
1736
1737
1738 fout->Close();
1739 fout->Delete();
1740 gFile = currentFile;
1741 gDirectory = curr;
1742} // void CbmL1::FieldApproxCheck()
1743
1744// ---------------------------------------------------------------------------------------------------------------------
1745//
1746void CbmL1::DumpMCTripletsToTree()
1747{
1748 if (!fpMcTripletsTree) {
1749 TDirectory* currentDir = gDirectory;
1750 TFile* currentFile = gFile;
1751
1752
1753 // Get prefix and directory
1754 boost::filesystem::path p = (FairRunAna::Instance()->GetUserOutputFileName()).Data();
1755 std::string dir = p.parent_path().string();
1756 if (dir.empty()) dir = ".";
1757 std::string filename = dir + "/" + fsMcTripletsOutputFilename + "." + p.filename().string();
1758 LOG(info) << "\033[1;31mSAVING A TREE: " << filename << "\033[0m";
1759
1760 fpMcTripletsOutFile = new TFile(filename.c_str(), "RECREATE");
1761 fpMcTripletsOutFile->cd();
1762 fpMcTripletsTree = new TTree("t", "MC Triplets");
1763 // motherId - id of mother particle (< 0, if primary)
1764 // pdg - PDG code of particle
1765 // processId - id of the process (ROOT::TMCProcessID)
1766 // pq - absolute value of track momentum [GeV/c], divided by charge [e]
1767 // zv - z-component of track vertex [cm]
1768 // s - global index of station
1769 // x0, y0, z0 - position of the left MC point in a triplet [cm]
1770 // x1, y1, z1 - position of the middle MC point in a triplet [cm]
1771 // x2, y2, z2 - position of the right MC point in a triplet [cm]
1772
1773 gFile = currentFile;
1774 gDirectory = currentDir;
1775 }
1776
1777 // Variables for tree branches
1778 int brMotherId = -1;
1779 int brPdg = -1;
1780 unsigned int brProcId = (unsigned int) -1;
1781 float brP = 0.f;
1782 int brQ = 0;
1783 float brVertexZ = 0.f;
1784 int brStation = -1;
1785
1786 float brX0 = 0.f;
1787 float brY0 = 0.f;
1788 float brZ0 = 0.f;
1789 float brX1 = 0.f;
1790 float brY1 = 0.f;
1791 float brZ1 = 0.f;
1792 float brX2 = 0.f;
1793 float brY2 = 0.f;
1794 float brZ2 = 0.f;
1795
1796 // Register branches
1797 fpMcTripletsTree->Branch("brMotherId", &brMotherId, "motherId/I");
1798 fpMcTripletsTree->Branch("brPdg", &brPdg, "pdg/I");
1799 fpMcTripletsTree->Branch("brProcId", &brProcId, "processId/i");
1800 fpMcTripletsTree->Branch("brP", &brP, "p/F");
1801 fpMcTripletsTree->Branch("brQ", &brQ, "q/I");
1802 fpMcTripletsTree->Branch("brVertexZ", &brVertexZ, "zv/F");
1803 fpMcTripletsTree->Branch("brStation", &brStation, "s/I");
1804 fpMcTripletsTree->Branch("brX0", &brX0, "x0/F");
1805 fpMcTripletsTree->Branch("brY0", &brY0, "y0/F");
1806 fpMcTripletsTree->Branch("brZ0", &brZ0, "z0/F");
1807 fpMcTripletsTree->Branch("brX1", &brX1, "x1/F");
1808 fpMcTripletsTree->Branch("brY1", &brY1, "y1/F");
1809 fpMcTripletsTree->Branch("brZ1", &brZ1, "z1/F");
1810 fpMcTripletsTree->Branch("brX2", &brX2, "x2/F");
1811 fpMcTripletsTree->Branch("brY2", &brY2, "y2/F");
1812 fpMcTripletsTree->Branch("brZ2", &brZ2, "z2/F");
1813
1814 struct ReducedMcPoint {
1815 int s;
1816 float x;
1817 float y;
1818 float z;
1819 };
1820
1821 for (const auto& tr : fMCData.GetTrackContainer()) {
1822 // Use only reconstructable tracks
1823 if (!tr.IsReconstructable()) {
1824 continue;
1825 }
1826
1827 std::vector<ReducedMcPoint> vPoints;
1828 vPoints.reserve(tr.GetNofPoints());
1829 for (unsigned int iP : tr.GetPointIndexes()) {
1830 const auto& point = fMCData.GetPoint(iP);
1831 vPoints.emplace_back(
1832 ReducedMcPoint{point.GetActiveStationId(), float(point.GetX()), float(point.GetY()), float(point.GetZ())});
1833 }
1834
1835 std::sort(vPoints.begin(), vPoints.end(),
1836 [](const ReducedMcPoint& lhs, const ReducedMcPoint& rhs) { return lhs.s < rhs.s; });
1837
1838 for (unsigned int i = 0; i + 2 < vPoints.size(); ++i) {
1839 // Condition to collect only triplets without gaps in stations
1840 // TODO: SZh 20.10.2022 Add cases for jump iterations
1841 if (vPoints[i + 1].s == vPoints[i].s + 1 && vPoints[i + 2].s == vPoints[i].s + 2) {
1842 // Fill MC-triplets tree
1843 brMotherId = tr.GetMotherId();
1844 brPdg = tr.GetPdgCode();
1845 brProcId = tr.GetProcessId();
1846 brP = tr.GetP();
1847 brQ = tr.GetCharge();
1848 brVertexZ = tr.GetStartZ();
1849 brStation = vPoints[i].s;
1850 brX0 = vPoints[i].x;
1851 brY0 = vPoints[i].y;
1852 brZ0 = vPoints[i].z;
1853 brX1 = vPoints[i + 1].x;
1854 brY1 = vPoints[i + 1].y;
1855 brZ1 = vPoints[i + 1].z;
1856 brX2 = vPoints[i + 2].x;
1857 brY2 = vPoints[i + 2].y;
1858 brZ2 = vPoints[i + 2].z;
1859
1860 fpMcTripletsTree->Fill();
1861 }
1862 }
1863 }
1864}
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.
Int_t
int Int_t
Definition RootTypesDef.h:16
B
Some class B.
Definition _GTestStructLooper.cxx:36
h
Generates beam ions for transport simulation.
CbmL1Track::GetMaxPurity
double GetMaxPurity() const
Gets max purity.
Definition CbmL1Track.h:68
CbmL1Track::Hits
std::vector< int > Hits
Indexes of hits of this track.
Definition CbmL1Track.h:101
CbmL1Track::IsGhost
bool IsGhost() const
Definition CbmL1Track.h:52
CbmL1Track::GetMatchedMCTrackIndex
int GetMatchedMCTrackIndex() const
Gets index of matched MC track.
Definition CbmL1Track.h:73
CbmL1Track::GetFirstHitIndex
int GetFirstHitIndex() const
Gets first hit index.
Definition CbmL1Track.h:55
CbmL1Track::GetNofMCTracks
int GetNofMCTracks() const
Gets number of associated MC tracks.
Definition CbmL1Track.h:79
CbmL1Track::GetLastHitIndex
int GetLastHitIndex() const
Gets last hit index.
Definition CbmL1Track.h:58
CbmL1::fNStations
int fNStations
number of total active detector stations
Definition CbmL1.h:435
CbmL1::fHistoDir
TDirectory * fHistoDir
Definition CbmL1.h:495
CbmL1::fTrackingTime
double fTrackingTime
time of track finding procedure
Definition CbmL1.h:438
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:505
CbmL1::fsMcTripletsOutputFilename
std::string fsMcTripletsOutputFilename
Name of file to save MC-triplets tree.
Definition CbmL1.h:506
CbmL1::fMCData
cbm::algo::ca::McData fMCData
MC Data object.
Definition CbmL1.h:410
CbmL1::fpAlgo
ca::Framework * fpAlgo
Pointer to the L1 track finder algorithm.
Definition CbmL1.h:418
CbmL1::FillFitHistos
void FillFitHistos(cbm::algo::kf::TrackParamV &tr, const cbm::algo::ca::McPoint &mc, bool isTimeFitted, TH1F *h[])
Definition CbmL1Performance.cxx:1523
CbmL1::fTableDir
TDirectory * fTableDir
Definition CbmL1.h:496
CbmL1::FieldIntegralCheck
void FieldIntegralCheck()
Definition CbmL1Performance.cxx:1679
CbmL1::fvRecoTracks
ca::Vector< CbmL1Track > fvRecoTracks
Reconstructed tracks container.
Definition CbmL1.h:422
CbmL1::TrackFitPerformance
void TrackFitPerformance()
Definition CbmL1Performance.cxx:1075
CbmL1::DumpMCTripletsToTree
void DumpMCTripletsToTree()
Definition CbmL1Performance.cxx:1746
CbmL1::HistoPerformance
void HistoPerformance()
Fills performance histograms.
Definition CbmL1Performance.cxx:432
CbmL1::fvHitDebugInfo
ca::Vector< cbm::algo::ca::McHitInfo > fvHitDebugInfo
Container of hits with extended information.
Definition CbmL1.h:482
CbmL1::FieldApproxCheck
void FieldApproxCheck()
Definition CbmL1Performance.cxx:1570
CbmL1::fpMcTripletsOutFile
TFile * fpMcTripletsOutFile
File to save MC-triplets tree.
Definition CbmL1.h:504
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
Data
Class to handle QA-objects in the online reconstruction.
Definition QaData.h:27
Table
Table element in the report.
Definition CbmQaReportTable.h:21
cbm::algo::ca::Hit
ca::Hit class describes a generic hit for the CA tracker
Definition CaHit.h:32
cbm::algo::ca::McHitInfo::dy
double dy
y coordinate error [cm]
Definition CaMcHitInfo.h:125
cbm::algo::ca::McHitInfo::y
double y
y coordinate of position [cm]
Definition CaMcHitInfo.h:121
cbm::algo::ca::McHitInfo::iStation
int iStation
index of station in active stations array
Definition CaMcHitInfo.h:118
cbm::algo::ca::McHitInfo::dx
double dx
x coordinate error [cm]
Definition CaMcHitInfo.h:124
cbm::algo::ca::McHitInfo::x
double x
x coordinate of position [cm]
Definition CaMcHitInfo.h:120
cbm::algo::ca::McPoint
Class describes a unified MC-point, used in CA tracking QA analysis.
Definition CaMcPoint.h:33
cbm::algo::ca::McPoint::GetInvSpeed
int GetInvSpeed() const
Gets inverse speed at reference z of station [ns/cm].
Definition CaMcPoint.h:90
cbm::algo::ca::McPoint::GetTxOut
double GetTxOut() const
Gets slope along x-axis at exit of station.
Definition CaMcPoint.h:219
cbm::algo::ca::McPoint::GetTime
double GetTime() const
Gets time [ns].
Definition CaMcPoint.h:207
cbm::algo::ca::McPoint::GetP
double GetP() const
Gets track momentum absolute value at reference z of station [GeV/c].
Definition CaMcPoint.h:120
cbm::algo::ca::McPoint::GetYOut
double GetYOut() const
Gets y coordinate at exit of station [cm].
Definition CaMcPoint.h:246
cbm::algo::ca::McPoint::GetXOut
double GetXOut() const
Gets x coordinate at exit of station [cm].
Definition CaMcPoint.h:237
cbm::algo::ca::McPoint::GetQp
double GetQp() const
Gets track charge over momentum at reference z of station [ec/GeV].
Definition CaMcPoint.h:189
cbm::algo::ca::McPoint::GetTyOut
double GetTyOut() const
Gets slope along x-axis at exit of station.
Definition CaMcPoint.h:228
cbm::algo::ca::McPoint::GetZOut
double GetZOut() const
Gets z coordinate at exit of station [cm].
Definition CaMcPoint.h:255
cbm::algo::ca::utils::Debugger::Instance
static Debugger & Instance()
Instance.
Definition CaDebugger.cxx:49
cbm::algo::ca::utils::Debugger::FillNtuple
virtual void FillNtuple(const char *name, float v[])=0
Add an entry to ntuple.
cbm::algo::ca::utils::Debugger::AddNtuple
virtual void AddNtuple(const char *name, const char *varlist)=0
Set new ntuple.
cbm::algo::kf::FieldRegion
Magnetic field region, corresponding to a hit triplet.
Definition KfFieldRegion.h:307
cbm::algo::kf::GlobalField::fgOriginalFieldType
static EFieldType fgOriginalFieldType
Global field type.
Definition KfFieldRegion.h:43
cbm::algo::kf::GlobalField::fgOriginalField
static FieldFn_t fgOriginalField
Global variable to store the field function (x,y,z)->(Bx,By,Bz)
Definition KfFieldRegion.h:40
cbm::algo::kf::TrackKalmanFilter::EnergyLossCorrection
void EnergyLossCorrection(DataT radThick, FitDirection direction)
Definition KfTrackKalmanFilter.cxx:1235
cbm::algo::kf::TrackKalmanFilter::Extrapolate
void Extrapolate(DataT z, const kf::FieldRegion< DataT > &F)
Definition KfTrackKalmanFilter.h:358
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:1157
cbm::algo::kf::TrackKalmanFilter::SetParticleMass
void SetParticleMass(DataT mass)
set particle mass for the fit
Definition KfTrackKalmanFilter.h:137
cbm::algo::kf::TrackKalmanFilter::Tr
kf::TrackParam< DataT > & Tr()
Definition KfTrackKalmanFilter.h:125
cbm::algo::kf::TrackKalmanFilter::SetTrack
void SetTrack(const kf::TrackParam< T > &t)
Definition KfTrackKalmanFilter.h:102
cbm::algo::kf::TrackParamBase::GetP
T GetP() const
Gets momentum [GeV/ec]. For the straight tracks returns 1.e4 [GeV/ec].
Definition KfTrackParam.h:290
cbm::algo::kf::TrackParamBase::GetPhi
T GetPhi() const
Gets azimuthal angle [rad].
Definition KfTrackParam.h:278
cbm::algo::kf::TrackParamBase::GetQp
T GetQp() const
Gets charge over momentum [ec/GeV].
Definition KfTrackParam.h:236
cbm::algo::kf::TrackParamBase::GetNdf
T GetNdf() const
Gets NDF of track fit model.
Definition KfTrackParam.h:266
cbm::algo::kf::TrackParamBase::GetY
T GetY() const
Gets y position [cm].
Definition KfTrackParam.h:218
cbm::algo::kf::TrackParamBase::GetChiSq
T GetChiSq() const
Gets Chi-square of track fit model.
Definition KfTrackParam.h:263
cbm::algo::kf::TrackParamBase::GetX
T GetX() const
Gets x position [cm].
Definition KfTrackParam.h:212
cbm::algo::kf::fmask::One
static fmask One()
Definition KfSimdPseudo.h:60
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:86
cbm::algo::ca::constants::phys::SpeedOfLightInvD
constexpr double SpeedOfLightInvD
Inverse speed of light [ns/cm].
Definition CaDefs.h:87
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
cbm::algo::kf::EFieldType::Normal
@ Normal
Field near the tracker subsystem.
Definition KfDefs.h:127
cbm::algo::kf::TrackParamV
TrackParam< fvec > TrackParamV
Definition KfTrackParam.h:675
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::counters
cacore::Vector< T > counters
Definition CbmL1Counters.h:122