CbmRoot
Loading...
Searching...
No Matches
CbmKFTrackFitQa.cxx
Go to the documentation of this file.
1/* Copyright (C) 2011-2019 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: Maksym Zyzak [committer] */
4
5/*
6 *====================================================================
7 *
8 * KF Fit performance
9 *
10 *====================================================================
11 */
12#include "CbmKFTrackFitQa.h"
13
14#include "CbmKF.h"
15#include "CbmKFMath.h"
16#include "CbmKFTrack.h"
17#include "CbmMCTrack.h"
18#include "CbmMatch.h"
19#include "CbmStsTrack.h"
20#include "CbmTrackMatch.h"
21#include "TBranch.h"
22#include "TDatabasePDG.h"
23#include "TDirectory.h"
24#include "TFile.h"
25#include "TTree.h"
26
27#include <FairMCPoint.h>
28#include <FairRootManager.h>
29
30#include <cmath>
31#include <iostream>
32using namespace std;
33using std::vector;
34
35ClassImp(CbmKFTrackFitQa)
36
37 CbmKFTrackFitQa::CbmKFTrackFitQa()
38 : listStsPts(nullptr)
39 , listMvdPts(nullptr)
40 , listMCTracks(nullptr)
41 , listStsTracksMatch(nullptr)
42 , listStsTracks(nullptr)
43 , listStsHits(nullptr)
44 , listMvdHits(nullptr)
45 , listMvdHitMatches(nullptr)
46 , listStsClusters(nullptr)
47 , listStsDigi(nullptr)
48 , listStsDigiMatch(nullptr)
49 ,
50
51
52 // Names of files
53 outfileName("outCbmTrackError.root")
54 ,
55
56 vStsHitMatch()
57 ,
58
59 res_STShit_x(nullptr)
60 , res_STShit_y(nullptr)
61 , pull_STShit_x(nullptr)
62 , pull_STShit_y(nullptr)
63 ,
64
65 res_MVDhit_x(nullptr)
66 , res_MVDhit_y(nullptr)
67 , pull_MVDhit_x(nullptr)
68 , pull_MVDhit_y(nullptr)
69 ,
70
71 res_AtPV_x(nullptr)
72 , res_AtPV_y(nullptr)
73 , res_AtPV_tx(nullptr)
74 , res_AtPV_ty(nullptr)
75 , res_AtPV_qp(nullptr)
76 ,
77
78 pull_AtPV_x(nullptr)
79 , pull_AtPV_y(nullptr)
80 , pull_AtPV_tx(nullptr)
81 , pull_AtPV_ty(nullptr)
82 , pull_AtPV_qp(nullptr)
83 ,
84
85
86 res_AtFP_x(nullptr)
87 , res_AtFP_y(nullptr)
88 , res_AtFP_tx(nullptr)
89 , res_AtFP_ty(nullptr)
90 , res_AtFP_qp(nullptr)
91 ,
92
93 pull_AtFP_x(nullptr)
94 , pull_AtFP_y(nullptr)
95 , pull_AtFP_tx(nullptr)
96 , pull_AtFP_ty(nullptr)
97 , pull_AtFP_qp(nullptr)
98 ,
99
100 q_QA(nullptr)
101 , dp_p(nullptr)
102 ,
103
104 ggg(nullptr)
105 ,
106
107 Nback(0)
108{
109 TDirectory* currentDir = gDirectory;
110 gDirectory->cd(outfileName.Data());
111
112 ggg = new TH1F("ggg", "ggg", 6, -0.5, 5.5);
113 q_QA = new TProfile("q_QA", "q quality", 100, 0.0, 1.0, 0.0, 100.0);
114 q_QA->GetXaxis()->SetTitle("p, GeV/c");
115 q_QA->GetYaxis()->SetTitle("Q determenition efficiency, %");
116
117 dp_p = new TProfile("dp_p", "dp_p", 100, 0.0, 1.0, 0.0, 5.0);
118 dp_p->GetXaxis()->SetTitle("p, GeV/c");
119 dp_p->GetYaxis()->SetTitle("#Deltap/p, %");
120
121 res_STShit_x = new TH1F("residual_STShit_x", "residual_STShit_x", 200, -10., 10.);
122 res_STShit_x->GetXaxis()->SetTitle("dX, um");
123 res_STShit_y = new TH1F("residual_STShit_y", "residual_STShit_y", 200, -10., 10.);
124 res_STShit_y->GetXaxis()->SetTitle("dY, um");
125 pull_STShit_x = new TH1F("pull_STShit_x", "pull_STShit_x", 100, -15., 15.);
126 pull_STShit_y = new TH1F("pull_STShit_y", "pull_STShit_y", 100, -15., 15.);
127
128 res_MVDhit_x = new TH1F("residual_MVDhit_x", "residual_MVDhit_x", 200, -30., 30.);
129 res_MVDhit_x->GetXaxis()->SetTitle("dX, um");
130 res_MVDhit_y = new TH1F("residual_MVDhit_y", "residual_MVDhit_y", 200, -30., 30.);
131 res_MVDhit_y->GetXaxis()->SetTitle("dY, um");
132 pull_MVDhit_x = new TH1F("pull_MVDhit_x", "pull_MVDhit_x", 100, -5., 5.);
133 pull_MVDhit_y = new TH1F("pull_MVDhit_y", "pull_MVDhit_y", 100, -5., 5.);
134
135
136 res_AtPV_x = new TH1F("residual_AtPV_x", "residual_AtPV_x", 2000, -1., 1.);
137 res_AtPV_x->GetXaxis()->SetTitle("dX, cm");
138 res_AtPV_y = new TH1F("residual_AtPV_y", "residual_AtPV_y", 2000, -1., 1.);
139 res_AtPV_y->GetXaxis()->SetTitle("dY, cm");
140 res_AtPV_tx = new TH1F("residual_AtPV_tx", "residual_AtPV_tx", 200, -0.004, 0.004);
141 res_AtPV_tx->GetXaxis()->SetTitle("dtx");
142 res_AtPV_ty = new TH1F("residual_AtPV_ty", "residual_AtPV_ty", 200, -0.004, 0.004);
143 res_AtPV_ty->GetXaxis()->SetTitle("dty");
144 res_AtPV_qp = new TH1F("residual_AtPV_qp", "residual_AtPV_qp", 200, -0.05, 0.05);
145 res_AtPV_qp->GetXaxis()->SetTitle("d(Q/P), c/GeV");
146
147 pull_AtPV_x = new TH1F("pull_AtPV_x", "pull_AtPV_x", 100, -5., 5.);
148 pull_AtPV_y = new TH1F("pull_AtPV_y", "pull_AtPV_y", 100, -5., 5.);
149 pull_AtPV_tx = new TH1F("pull_AtPV_tx", "pull_AtPV_tx", 100, -5., 5.);
150 pull_AtPV_ty = new TH1F("pull_AtPV_ty", "pull_AtPV_ty", 100, -5., 5.);
151 pull_AtPV_qp = new TH1F("pull_AtPV_qp", "pull_AtPV_qp", 100, -5., 5.);
152
153 res_AtFP_x = new TH1F("residual_AtFP_x", "residual_AtFP_x", 200, -50., 50.);
154 res_AtFP_x->GetXaxis()->SetTitle("dX, cm");
155 res_AtFP_y = new TH1F("residual_AtFP_y", "residual_AtFP_y", 200, -400., 400.);
156 res_AtFP_y->GetXaxis()->SetTitle("dY, cm");
157 res_AtFP_tx = new TH1F("residual_AtFP_tx", "residual_AtFP_tx", 200, -0.004, 0.004);
158 res_AtFP_tx->GetXaxis()->SetTitle("dtx, GeV/c");
159 res_AtFP_ty = new TH1F("residual_AtFP_ty", "residual_AtFP_ty", 200, -0.004, 0.004);
160 res_AtFP_ty->GetXaxis()->SetTitle("dty, GeV/c");
161 res_AtFP_qp = new TH1F("residual_AtFP_qp", "residual_AtFP_qp", 200, -0.05, 0.05);
162 res_AtFP_qp->GetXaxis()->SetTitle("d(Q/P), c/GeV");
163
164 pull_AtFP_x = new TH1F("pull_AtFP_x", "pull_AtFP_x", 100, -5., 5.);
165 pull_AtFP_y = new TH1F("pull_AtFP_y", "pull_AtFP_y", 100, -5., 5.);
166 pull_AtFP_tx = new TH1F("pull_AtFP_tx", "pull_AtFP_tx", 100, -5., 5.);
167 pull_AtFP_ty = new TH1F("pull_AtFP_ty", "pull_AtFP_ty", 100, -5., 5.);
168 pull_AtFP_qp = new TH1F("pull_AtFP_qp", "pull_AtFP_qp", 100, -5., 5.);
169
170 gDirectory = currentDir;
171}
172
173CbmKFTrackFitQa::~CbmKFTrackFitQa()
174{
175 delete res_STShit_x;
176 delete res_STShit_y;
177 delete pull_STShit_x;
178 delete pull_STShit_y;
179
180 delete res_MVDhit_x;
181 delete res_MVDhit_y;
182 delete pull_MVDhit_x;
183 delete pull_MVDhit_y;
184
185 delete res_AtPV_x;
186 delete res_AtPV_y;
187 delete res_AtPV_tx;
188 delete res_AtPV_ty;
189 delete res_AtPV_qp;
190
191 delete pull_AtPV_x;
192 delete pull_AtPV_y;
193 delete pull_AtPV_tx;
194 delete pull_AtPV_ty;
195 delete pull_AtPV_qp;
196
197 delete res_AtFP_x;
198 delete res_AtFP_y;
199 delete res_AtFP_tx;
200 delete res_AtFP_ty;
201 delete res_AtFP_qp;
202
203 delete pull_AtFP_x;
204 delete pull_AtFP_y;
205 delete pull_AtFP_tx;
206 delete pull_AtFP_ty;
207 delete pull_AtFP_qp;
208}
209
210InitStatus CbmKFTrackFitQa::ReInit() { return Init(); }
211
212InitStatus CbmKFTrackFitQa::Init()
213{
214 FairRootManager* fManger = FairRootManager::Instance();
215
216 listMCTracks = dynamic_cast<TClonesArray*>(fManger->GetObject("MCTrack"));
217 listStsPts = dynamic_cast<TClonesArray*>(fManger->GetObject("StsPoint"));
218 listStsTracks = dynamic_cast<TClonesArray*>(fManger->GetObject("StsTrack"));
219 listStsTracksMatch = dynamic_cast<TClonesArray*>(fManger->GetObject("StsTrackMatch"));
220 listStsHits = dynamic_cast<TClonesArray*>(fManger->GetObject("StsHit"));
221 listStsClusters = dynamic_cast<TClonesArray*>(fManger->GetObject("StsCluster"));
222 listStsDigi = dynamic_cast<TClonesArray*>(fManger->GetObject("StsDigi"));
223 listStsDigiMatch = dynamic_cast<TClonesArray*>(fManger->GetObject("StsDigiMatch"));
224
225 if (CbmKF::Instance()->vMvdMaterial.size() == 0) {
226 listMvdPts = nullptr;
227 listMvdHits = nullptr;
228 listMvdHitMatches = nullptr;
229 }
230 else {
231 listMvdPts = dynamic_cast<TClonesArray*>(fManger->GetObject("MvdPoint"));
232 listMvdHits = dynamic_cast<TClonesArray*>(fManger->GetObject("MvdHit"));
233 listMvdHitMatches = dynamic_cast<TClonesArray*>(fManger->GetObject("MvdHitMatch"));
234 }
235
236 return kSUCCESS;
237}
238
239void CbmKFTrackFitQa::Exec(Option_t* /*option*/)
240{
241 FillHitHistos();
242
243 CbmKFTrErrMCPoints* MCTrackSortedArray = new CbmKFTrErrMCPoints[listMCTracks->GetEntriesFast() + 1];
244 if (CbmKF::Instance()->vMvdMaterial.size() != 0) {
245 for (Int_t iMvd = 0; iMvd < listMvdPts->GetEntriesFast(); iMvd++) {
246 CbmMvdPoint* MvdPoint = (CbmMvdPoint*) listMvdPts->At(iMvd);
247 //MCTrackSortedArray[MvdPoint->GetTrackID()].SetMvdPoint(MvdPoint);
248 MCTrackSortedArray[MvdPoint->GetTrackID()].MvdArray.push_back(MvdPoint);
249 }
250 }
251 for (Int_t iSts = 0; iSts < listStsPts->GetEntriesFast(); iSts++) {
252 CbmStsPoint* StsPoint = (CbmStsPoint*) listStsPts->At(iSts);
253 //MCTrackSortedArray[StsPoint->GetTrackID()].SetStsPoint(StsPoint);
254 MCTrackSortedArray[StsPoint->GetTrackID()].StsArray.push_back(StsPoint);
255 }
256 for (Int_t itrack = 0; itrack < listStsTracks->GetEntriesFast(); itrack++) {
257 CbmStsTrack* StsTrack = (CbmStsTrack*) listStsTracks->At(itrack);
258 CbmTrackMatch* StsTrackMatch = (CbmTrackMatch*) listStsTracksMatch->At(itrack);
259 if (StsTrackMatch->GetNofMCTracks() == 0) {
260 continue;
261 }
262 CbmMCTrack* MCTrack = (CbmMCTrack*) listMCTracks->At(StsTrackMatch->GetMCTrackId());
263 CbmKFTrack KFTrack(*StsTrack);
264 FillHistoAtFirstPoint(&MCTrackSortedArray[StsTrackMatch->GetMCTrackId()], MCTrack, &KFTrack);
265 // if(MCTrack->GetP() >= 1 && MCTrack->GetMotherId()==-1 && MCTrackSortedArray[StsTrackMatch->GetMCTrackId()].MvdArray.size()>1)
266 // {
267 FillHistoAtParticleVertex(MCTrack, &KFTrack);
268 // }
269 }
270 delete[] MCTrackSortedArray;
271}
272void CbmKFTrackFitQa::Finish() { KFWrite(); }
273
274void CbmKFTrackFitQa::FillHistoAtParticleVertex(CbmMCTrack* track_mc, CbmKFTrack* track_kf)
275{
276 CbmKFTrackInterface* tr_help1 = track_kf;
277 CbmKFTrack* tr1 = new CbmKFTrack(*tr_help1);
278 CbmKFTrackInterface* tr_help = tr1;
279
280 tr_help->Extrapolate(track_mc->GetStartZ()); //extrapolating of the track to the prodaction point
281
282 Double_t* fT = tr_help->GetTrack();
283 Double_t* fC = tr_help->GetCovMatrix();
284
285 // getting the q (charge) of the mc track using the pdg information
286
287 Double_t qtrack = 0;
288 if (track_mc->GetPdgCode() < 9999999) {
289 qtrack = TDatabasePDG::Instance()->GetParticle(track_mc->GetPdgCode())->Charge() / 3.0;
290 }
291 else {
292 qtrack = 0;
293 }
294
295 Double_t PointPx = track_mc->GetPx();
296 Double_t PointPy = track_mc->GetPy();
297 Double_t PointPz = track_mc->GetPz();
298 Double_t P_mc = sqrt(PointPx * PointPx + PointPy * PointPy + PointPz * PointPz);
299
300 //differences of the KF track and MC track parameters calculation
301 Double_t ddx =
302 fT[0] - track_mc->GetStartX(); // The difference between x coordinates of the reconstructed and MC tracks
303 Double_t ddy =
304 fT[1] - track_mc->GetStartY(); // The difference between y coordinates of the reconstructed and MC tracks
305 Double_t ddtx = fT[2] - PointPx / PointPz; // The difference between tx of the reconstructed and MC tracks
306 Double_t ddty = fT[3] - PointPy / PointPz; // The difference between ty of the reconstructed and MC tracks
307 Double_t ddqp = (fabs(1. / fT[4]) - P_mc) / P_mc; // p resolution
308 Double_t ddqp_p = fT[4] - (qtrack / sqrt(PointPx * PointPx + PointPy * PointPy + PointPz * PointPz)); //p residual
309
310 //differences of the KF track and MC track parameters
311 res_AtPV_x->Fill(ddx);
312 res_AtPV_y->Fill(ddy);
313 res_AtPV_tx->Fill(ddtx);
314 res_AtPV_ty->Fill(ddty);
315 if (isfinite(fT[4]) && (fabs(fT[4]) > 1.e-20)) {
316 res_AtPV_qp->Fill(ddqp);
317 }
318 //pulls of the parameters
319 if (isfinite(fC[0]) && fC[0] > 0) {
320 pull_AtPV_x->Fill(ddx / sqrt(fC[0]));
321 }
322 if (isfinite(fC[2]) && fC[2] > 0) {
323 pull_AtPV_y->Fill(ddy / sqrt(fC[2]));
324 }
325 if (isfinite(fC[5]) && fC[5] > 0) {
326 pull_AtPV_tx->Fill(ddtx / sqrt(fC[5]));
327 }
328 if (isfinite(fC[9]) && fC[9] > 0) {
329 pull_AtPV_ty->Fill(ddty / sqrt(fC[9]));
330 }
331 if (isfinite(fC[14]) && fC[14] > 0) {
332 pull_AtPV_qp->Fill(ddqp_p / sqrt(fC[14]));
333 }
334
335 if (isfinite(fT[4]) && (fabs(fT[4]) > 1.e-20)) {
336 if (qtrack == (fabs(fT[4]) / fT[4])) {
337 q_QA->Fill(P_mc, 100.0);
338 }
339 else {
340 q_QA->Fill(P_mc, 0.0);
341 }
342
343 if (isfinite(fC[14]) && fC[14] > 0) {
344 dp_p->Fill(P_mc, fabs(1. / fT[4]) * sqrt(fC[14]) * 100, 1);
345 }
346 }
347}
348
349void CbmKFTrackFitQa::FillHistoAtFirstPoint(CbmKFTrErrMCPoints* mc_points, CbmMCTrack* track_mc, CbmKFTrack* track_kf)
350{
351 Double_t* fT = track_kf->GetTrack();
352 Double_t* fC = track_kf->GetCovMatrix();
353
354 //Double_t StartX = track_kf -> GetHit(0)[];
355 Bool_t nomvdpoint = 1;
356
357 FairMCPoint* MCFirstPoint;
358
359 // if (fabs(track_mc->GetP()) < 1.) return;
360 // if (fabs(track_mc->GetPz()) < 0.0001) return;
361 //if (fT[4]<0.) return;
362
363 if (!mc_points) {
364 return;
365 }
366 if ((mc_points->MvdArray.size() == 0) && (mc_points->StsArray.size() == 0)) {
367 return;
368 }
369
370 if (mc_points->MvdArray.size() > 0) {
371 MCFirstPoint = *mc_points->MvdArray.begin();
372 }
373 else {
374 MCFirstPoint = *mc_points->StsArray.begin();
375 }
376
377 for (vector<CbmMvdPoint*>::iterator imvd = mc_points->MvdArray.begin(); imvd != mc_points->MvdArray.end(); ++imvd) {
378 MCFirstPoint = *imvd;
379 //cout << "Z " << MCFirstPoint->GetZ() << " fT[5] " << fT[5] << " mctrackZ " << track_mc->GetStartZ() << " X " << fT[0] << " Xmc "<< track_mc->GetStartX()<<
380 //" MCID " << track_mc->GetMotherId()<< " MCIDPart " << MCFirstPoint->GetTrackID() << endl;
381 if (fabs(MCFirstPoint->GetZ() - fT[5]) < 1.) {
382 nomvdpoint = 0;
383 break;
384 }
385 }
386 if (nomvdpoint) {
387 for (vector<CbmStsPoint*>::iterator ists = mc_points->StsArray.begin(); ists != mc_points->StsArray.end(); ++ists) {
388 MCFirstPoint = *ists;
389 //cout << "Z " << MCFirstPoint->GetZ() << " fT[5] " << fT[5] << " mctrackZ " << track_mc->GetStartZ() << " X " << fT[0] << " Xmc "<< track_mc->GetStartX()<<
390 //" MCID " << track_mc->GetMotherId() << " MCIDPart " << MCFirstPoint->GetTrackID() << endl;
391 if (fabs(MCFirstPoint->GetZ() - fT[5]) < 1.) {
392 break;
393 }
394 }
395 }
396 track_kf->Extrapolate(MCFirstPoint->GetZ());
397
398
399 // getting of the q (charge) of the mc track using the pdg information
400 {
401 Double_t qtrack = 0;
402 if (track_mc->GetPdgCode() < 9999999) {
403 qtrack = TDatabasePDG::Instance()->GetParticle(track_mc->GetPdgCode())->Charge() / 3.0;
404 }
405 else {
406 qtrack = 0;
407 }
408
409 Double_t PointPx = MCFirstPoint->GetPx();
410 Double_t PointPy = MCFirstPoint->GetPy();
411 Double_t PointPz = MCFirstPoint->GetPz();
412 Double_t P_mc = sqrt(PointPx * PointPx + PointPy * PointPy + PointPz * PointPz);
413
414 //differences of the KF track and MC track parameters calculation
415
416 Double_t ddx =
417 fT[0] - MCFirstPoint->GetX(); // The difference between x coordinates of the reconstructed and MC tracks
418 Double_t ddy =
419 fT[1] - MCFirstPoint->GetY(); // The difference between y coordinates of the reconstructed and MC tracks
420 Double_t ddtx = fT[2] - PointPx / PointPz; // The difference between tx of the reconstructed and MC tracks
421 Double_t ddty = fT[3] - PointPy / PointPz; // The difference between ty of the reconstructed and MC tracks
422 Double_t ddqp = (fabs(1. / fT[4]) - P_mc) / P_mc; // p resolution
423 Double_t ddqp_p = fT[4] - (qtrack / sqrt(PointPx * PointPx + PointPy * PointPy + PointPz * PointPz)); //p residual
424
425 //differences of the KF track and MC track parameters
426 res_AtFP_x->Fill(ddx * 10000.);
427 res_AtFP_y->Fill(ddy * 10000.);
428 res_AtFP_tx->Fill(ddtx);
429 res_AtFP_ty->Fill(ddty);
430 if (isfinite(fT[4]) && (fabs(fT[4]) > 1.e-20)) {
431 res_AtFP_qp->Fill(ddqp);
432 }
433 //pulls of the parameters
434 if (isfinite(fC[0]) && fC[0] > 0) {
435 pull_AtFP_x->Fill(ddx / sqrt(fC[0]));
436 }
437 if (isfinite(fC[2]) && fC[2] > 0) {
438 pull_AtFP_y->Fill(ddy / sqrt(fC[2]));
439 }
440 if (isfinite(fC[5]) && fC[5] > 0) {
441 pull_AtFP_tx->Fill(ddtx / sqrt(fC[5]));
442 }
443 if (isfinite(fC[9]) && fC[9] > 0) {
444 pull_AtFP_ty->Fill(ddty / sqrt(fC[9]));
445 }
446 if (isfinite(fC[14]) && fC[14] > 0) {
447 pull_AtFP_qp->Fill(ddqp_p / sqrt(fC[14]));
448 }
449
450 if (isfinite(fT[4]) && (fabs(fT[4]) > 1.e-20)) {
451 if (qtrack == (fabs(fT[4]) / fT[4])) {
452 q_QA->Fill(P_mc, 100.0);
453 }
454 else {
455 q_QA->Fill(P_mc, 0.0);
456 }
457
458 if (isfinite(fC[14]) && fC[14] > 0) {
459 dp_p->Fill(P_mc, fabs(1. / fT[4]) * sqrt(fC[14]) * 100, 1);
460 }
461 //cout << ddx <<" "<< ddx/sqrt(fC[0]) << endl;
462 }
463 }
464}
465
466void CbmKFTrackFitQa::KFWrite()
467{
468 //This function writes obtained histograms to the root file
469 TDirectory* curr = gDirectory;
470 TFile* currentFile = gFile;
471 TFile* fout = new TFile(outfileName.Data(), "Recreate");
472
473 //differences of the KF track and MC track parameters
474 res_AtPV_x->Write();
475 res_AtPV_y->Write();
476 res_AtPV_tx->Write();
477 res_AtPV_ty->Write();
478 res_AtPV_qp->Write();
479 //pulls of the parameters
480 pull_AtPV_x->Write();
481 pull_AtPV_y->Write();
482 pull_AtPV_tx->Write();
483 pull_AtPV_ty->Write();
484 pull_AtPV_qp->Write();
485 //differences of the KF track and MC track parameters
486
487 res_AtFP_x->Write();
488 res_AtFP_y->Write();
489 res_AtFP_tx->Write();
490 res_AtFP_ty->Write();
491 res_AtFP_qp->Write();
492 //pulls of the parameters
493 pull_AtFP_x->Write();
494 pull_AtFP_y->Write();
495 pull_AtFP_tx->Write();
496 pull_AtFP_ty->Write();
497 pull_AtFP_qp->Write();
498
499 res_STShit_x->Write();
500 res_STShit_y->Write();
501 pull_STShit_x->Write();
502 pull_STShit_y->Write();
503
504 res_MVDhit_x->Write();
505 res_MVDhit_y->Write();
506 pull_MVDhit_x->Write();
507 pull_MVDhit_y->Write();
508
509 ggg->Write();
510 q_QA->Write();
511 dp_p->Write();
512
513 fout->Close();
514 fout->Delete();
515 gFile = currentFile;
516 gDirectory = curr;
517}
518
519void CbmKFTrackFitQa::Save()
520{
521 //This function writes obtained histograms to the root file
522
523 //differences of the KF track and MC track parameters
524 res_AtPV_x->GetXaxis()->SetTitle("dX, cm");
525 res_AtPV_x->SaveAs("res_AtPV_x.eps");
526 res_AtPV_y->GetXaxis()->SetTitle("dY, cm");
527 res_AtPV_y->SaveAs("res_AtPV_y.gif");
528 res_AtPV_tx->GetXaxis()->SetTitle("dtx");
529 res_AtPV_tx->SaveAs("res_AtPV_tx.gif");
530 res_AtPV_ty->GetXaxis()->SetTitle("dty");
531 res_AtPV_ty->SaveAs("res_AtPV_ty.gif");
532 res_AtPV_qp->GetXaxis()->SetTitle("d(Q/P), c/GeV");
533 res_AtPV_qp->SaveAs("res_AtPV_qp.gif");
534
535 //pulls of the parameters
536 pull_AtPV_x->SaveAs("pull_AtPV_x.gif");
537 pull_AtPV_y->SaveAs("pull_AtPV_y.gif");
538 pull_AtPV_tx->SaveAs("pull_AtPV_tx.gif");
539 pull_AtPV_ty->SaveAs("pull_AtPV_ty.gif");
540 pull_AtPV_qp->SaveAs("pull_AtPV_qp.gif");
541
542 //differences of the KF track and MC track parameters
543 res_AtFP_x->GetXaxis()->SetTitle("dX, cm");
544 res_AtFP_x->SaveAs("res_AtFP_x.gif");
545 res_AtFP_y->GetXaxis()->SetTitle("dY, cm");
546 res_AtFP_y->SaveAs("res_AtFP_y.gif");
547 res_AtFP_tx->GetXaxis()->SetTitle("dtx, GeV/c");
548 res_AtFP_tx->SaveAs("res_AtFP_tx.gif");
549 res_AtFP_ty->GetXaxis()->SetTitle("dty, GeV/c");
550 res_AtFP_ty->SaveAs("res_AtFP_ty.gif");
551 res_AtFP_qp->GetXaxis()->SetTitle("d(Q/P), c/GeV");
552 res_AtFP_qp->SaveAs("res_AtFP_qp.gif");
553
554 //pulls of the parameters
555 pull_AtFP_x->SaveAs("pull_AtFP_x.gif");
556 pull_AtFP_y->SaveAs("pull_AtFP_y.gif");
557 pull_AtFP_tx->SaveAs("pull_AtFP_tx.gif");
558 pull_AtFP_ty->SaveAs("pull_AtFP_ty.gif");
559 pull_AtFP_qp->SaveAs("pull_AtFP_qp.gif");
560}
561
562void CbmKFTrackFitQa::FillHitHistos()
563{
564 vStsHitMatch.resize(listStsHits->GetEntriesFast());
565 StsHitMatch();
566 for (Int_t iSts = 0; iSts < listStsHits->GetEntriesFast(); iSts++) {
567 CbmStsHit* StsHit = (CbmStsHit*) listStsHits->At(iSts);
568 if (vStsHitMatch[iSts] > -1) {
569 CbmStsPoint* StsP = (CbmStsPoint*) listStsPts->At(vStsHitMatch[iSts]);
570 Double_t dx = (StsHit->GetX() - 0.5 * (StsP->GetXIn() + StsP->GetXOut()));
571 Double_t dy = (StsHit->GetY() - 0.5 * (StsP->GetYIn() + StsP->GetYOut()));
572 res_STShit_x->Fill(dx);
573 res_STShit_y->Fill(dy);
574 pull_STShit_x->Fill(dx / (StsHit->GetDx()));
575 pull_STShit_y->Fill(dy / (StsHit->GetDy()));
576 }
577 }
578 // TODO correct the matching for MVD
579 // if(CbmKF::Instance()->vMvdMaterial.size() != 0)
580 // {
581 // for (Int_t iMvd=0; iMvd<listMvdHits->GetEntriesFast(); iMvd++)
582 // {
583 // CbmMvdHit* MvdHit = (CbmMvdHit*)listMvdHits->At(iMvd);
584 // CbmMvdHitMatch* MvdHitMatch = (CbmMvdHitMatch*)listMvdHitMatches->At(iMvd);
585 //
586 // if(MvdHitMatch->GetPointId() > -1)
587 // {
588 // CbmMvdPoint* MvdP = (CbmMvdPoint*) listMvdPts->At(MvdHitMatch->GetPointId());
589 // if(!MvdP) continue;
590 // Double_t dx = (MvdHit->GetX() - 0.5*(MvdP->GetX() + MvdP->GetXOut()));
591 // Double_t dy = (MvdHit->GetY() - 0.5*(MvdP->GetY() + MvdP->GetYOut()));
592 // res_MVDhit_x -> Fill( dx*10000. );
593 // res_MVDhit_y -> Fill( dy*10000. );
594 // pull_MVDhit_x -> Fill( dx/(MvdHit->GetDx()) );
595 // pull_MVDhit_y -> Fill( dy/(MvdHit->GetDy()) );
596 // }
597 // }
598 // }
599}
600
601void CbmKFTrackFitQa::StsHitMatch()
602{
603 const bool useLinks =
604 0; // 0 - use HitMatch, one_to_one; 1 - use FairLinks, many_to_many. Set 0 to switch to old definition of efficiency.
605 // TODO: fix trunk problem with links. Set useLinks = 1
606
607 for (int iH = 0; iH < listStsHits->GetEntriesFast(); iH++) {
608 CbmStsHit* stsHit = dynamic_cast<CbmStsHit*>(listStsHits->At(iH));
609 vStsHitMatch[iH] = -1;
610 int gggi = 0;
611 if (useLinks) {
612 if (listStsClusters) {
613 // if ( NLinks != 2 ) cout << "HitMatch: Error. Hit wasn't matched with 2 clusters." << endl;
614 // 1-st cluster
615 vector<int> stsPointIds; // save here all mc-points matched with first cluster
616 vector<int> stsPointIds_hit;
617 int iL = 0;
618 CbmLink link = stsHit->GetMatch()->GetLink(iL);
619 CbmStsCluster* stsCluster = dynamic_cast<CbmStsCluster*>(listStsClusters->At(link.GetIndex()));
620 int NLinks2 = stsCluster->GetMatch()->GetNofLinks();
621 for (int iL2 = 0; iL2 < NLinks2; iL2++) {
622 const CbmLink& link2 = stsCluster->GetMatch()->GetLink(iL2);
623 //CbmStsDigi *stsDigi = dynamic_cast<CbmStsDigi*>( listStsDigi->At( link2.GetIndex() ) );
624 // const int NLinks3 = stsDigi->GetNLinks();
625 CbmMatch* stsDigiMatch = dynamic_cast<CbmMatch*>(listStsDigiMatch->At(link2.GetIndex()));
626 const int NLinks3 = stsDigiMatch->GetNofLinks();
627 for (int iL3 = 0; iL3 < NLinks3; iL3++) {
628 CbmLink link3 = stsDigiMatch->GetLink(iL3);
629 // FairLink link3 = stsDigi->GetLink(iL3);
630 int stsPointId = link3.GetIndex();
631 stsPointIds.push_back(stsPointId);
632 } // for mcPoint
633 } // for digi
634 // 2-nd cluster
635 iL = 1;
636 link = stsHit->GetMatch()->GetLink(iL);
637 stsCluster = dynamic_cast<CbmStsCluster*>(listStsClusters->At(link.GetIndex()));
638 NLinks2 = stsCluster->GetMatch()->GetNofLinks();
639 for (int iL2 = 0; iL2 < NLinks2; iL2++) {
640 const CbmLink& link2 = stsCluster->GetMatch()->GetLink(iL2);
641 //CbmStsDigi *stsDigi = dynamic_cast<CbmStsDigi*>( listStsDigi->At( link2.GetIndex() ) );
642 CbmMatch* stsDigiMatch = dynamic_cast<CbmMatch*>(listStsDigiMatch->At(link2.GetIndex()));
643 // const int NLinks3 = stsDigi->GetNLinks();
644 const int NLinks3 = stsDigiMatch->GetNofLinks();
645 for (int iL3 = 0; iL3 < NLinks3; iL3++) {
646 // FairLink link3 = stsDigi->GetLink(iL3);
647 CbmLink link3 = stsDigiMatch->GetLink(iL3);
648 int stsPointId = link3.GetIndex();
649
650 if (!find(&(stsPointIds[0]), &(stsPointIds[stsPointIds.size()]), stsPointId)) {
651 continue; // check if first cluster matched with same mc-point
652 }
653 // CbmStsPoint *pt = dynamic_cast<CbmStsPoint*>( listStsPts->At(stsPointId) );
654 // if(!pt) continue;
655 // CbmMCTrack *MCTrack = dynamic_cast<CbmMCTrack*>( listMCTracks->At( pt->GetTrackID() ) );
656 // if ( !MCTrack ) continue;
657 // if ( abs(MCTrack->GetPdgCode()) >= 10000 ) continue;
658
659 stsPointIds_hit.push_back(stsPointId);
660 vStsHitMatch[iH] = stsPointId;
661 gggi++;
662 } // for mcPoint
663 } // for digi
664 } // if clusters
665 } // if useLinks
666 ggg->Fill(gggi);
667 } // for hits
668}
669
670
671void CbmKFTrackFitQa::FindBackTracks()
672{
673 /*
674 Nback=0;
675 int NallTracks=0;
676
677 for (Int_t ientr=0; ientr< treco->GetEntriesFast(); ientr++)
678// for (Int_t ientr=0; ientr< 1; ientr++)
679 {
680 treco->GetEntry(ientr);
681 tmc->GetEntry(ientr);
682 cout << "Entry " << ientr << endl;
683
684 CbmKFTrErrMCPoints MCTrackSortedArray[MCTrackArray->GetEntriesFast()+1];
685
686 for (Int_t iMvd=0; iMvd<MvdPointArray->GetEntriesFast(); iMvd++)
687 {
688 CbmMvdPoint* MvdPoint = (CbmMvdPoint*)MvdPointArray->At(iMvd);
689 //MCTrackSortedArray[MvdPoint->GetTrackID()].SetMvdPoint(MvdPoint);
690 MCTrackSortedArray[MvdPoint->GetTrackID()].MvdArray.push_back(MvdPoint);
691 }
692
693 for (Int_t iSts=0; iSts<StsPointArray->GetEntriesFast(); iSts++)
694 {
695 CbmStsPoint* StsPoint = (CbmStsPoint*)StsPointArray->At(iSts);
696 //MCTrackSortedArray[StsPoint->GetTrackID()].SetStsPoint(StsPoint);
697 MCTrackSortedArray[StsPoint->GetTrackID()].StsArray.push_back(StsPoint);
698 }
699
700 NallTracks+=StsTrackArray->GetEntriesFast();
701
702 for (Int_t itrack=0; itrack<StsTrackArray->GetEntriesFast(); itrack++)
703// for (Int_t itrack=0; itrack<1; itrack++)
704
705 {
706 CbmStsTrack* StsTrack = (CbmStsTrack*)StsTrackArray->At(itrack);
707 CbmTrackMatch* StsTrackMatch = (CbmTrackMatch*)StsTrackMatchArray->At(itrack);
708 Int_t NofMC = StsTrackMatch -> GetNofMCTracks();
709//cout << itrack << "!!!!!!!" << StsTrackArray->GetEntriesFast() << endl;
710// cout << StsTrackMatch -> GetNofMCTracks() << endl;
711 if(StsTrackMatch -> GetNofMCTracks() != 1) continue;
712 CbmMCTrack* MCTrack = (CbmMCTrack*)MCTrackArray->At(StsTrackMatch->GetMCTrackId());
713
714 CbmKFTrack *KFTrack;
715 KFTrack = new CbmKFTrack(*StsTrack);
716
717// FindBackTracks(&MCTrackSortedArray[StsTrackMatch->GetMCTrackId()], MCTrack, KFTrack, ientr);
718 FillHistoAtFirstPoint(&MCTrackSortedArray[StsTrackMatch->GetMCTrackId()], MCTrack, KFTrack);
719// FillHistoAtParticleVertex(MCTrack, KFTrack);
720 }
721 }*/
722}
723
724void CbmKFTrackFitQa::FindBackTracks(CbmKFTrErrMCPoints* /*mc_points*/, CbmMCTrack* /*track_mc*/,
725 CbmKFTrack* /*track_kf*/, int /*iEvent*/)
726{
727 /*
728 FILE *fBack = fopen("BackTracks.txt","a+");
729
730 Double_t *fT = track_kf -> GetTrack();
731 Double_t *fC = track_kf -> GetCovMatrix();
732
733 Bool_t back = 0;
734
735 FairMCPoint *MCPoint;
736 FairMCPoint *MCFirstPoint;
737
738 Bool_t nomvdpoint = 1;
739 Bool_t nostspoint = 1;
740
741 if(mc_points->MvdArray.size() > 0.)
742 {
743 track_kf -> Extrapolate(mc_points->MvdArray[0]->GetZ());
744 for( vector<CbmMvdPoint*>::iterator imvd = mc_points->MvdArray.begin(); imvd < mc_points->MvdArray.begin()+1; ++imvd)
745 {
746 MCFirstPoint = *imvd;
747 }
748 }
749 else
750 {
751 track_kf -> Extrapolate(mc_points->StsArray[0]->GetZ());
752 for( vector<CbmStsPoint*>::iterator ists = mc_points->StsArray.begin(); ists < mc_points->StsArray.begin()+1; ++ists)
753 {
754 MCFirstPoint = *ists;
755 }
756 }
757
758 for( vector<CbmMvdPoint*>::iterator imvd = mc_points->MvdArray.begin(); imvd != mc_points->MvdArray.end(); ++imvd)
759 {
760 MCPoint = *imvd;
761 if(MCPoint->GetPz()<0)
762 {
763 back=1;
764 }
765 }
766 for( vector<CbmStsPoint*>::iterator ists = mc_points->StsArray.begin(); ists != mc_points->StsArray.end(); ++ists)
767 {
768 MCPoint = *ists;
769 if(MCPoint->GetPz()<0)
770 {
771 back=1;
772 }
773 }
774 // getting of the q (charge) of the mc track using the pdg information
775 if (back)
776 {
777 cout <<"mc "<< sqrt(MCFirstPoint->GetPz()*MCFirstPoint->GetPz()+MCFirstPoint->GetPy()*MCFirstPoint->GetPy()+MCFirstPoint->GetPx()*MCFirstPoint->GetPx()) << endl;
778 cout <<"Px "<< MCFirstPoint->GetPx()<<" ";
779 cout <<"Py "<< MCFirstPoint->GetPy()<<" ";
780 cout <<"Pz "<< MCFirstPoint->GetPz()<<endl;
781 cout << "re " << fabs(1./fT[4])<<endl;
782 cout <<"PDG = "<<track_mc->GetPdgCode()<<endl;
783
784 Nback++;
785
786 Double_t qtrack=0;
787 if ( track_mc->GetPdgCode() < 9999999 )
788 qtrack = TDatabasePDG::Instance()->GetParticle(track_mc->GetPdgCode())->Charge()/3.0;
789 else qtrack = 0;
790
791 Double_t PointPx = MCFirstPoint->GetPx();
792 Double_t PointPy = MCFirstPoint->GetPy();
793 Double_t PointPz = MCFirstPoint->GetPz();
794
795 //differences of the KF track and MC track parameters calculation
796 Double_t ddx = fT[0] - MCFirstPoint->GetX(); // The difference between the x coordinates of the reconstructed and MC tracks
797 Double_t ddy = fT[1] - MCFirstPoint->GetY(); // The difference between the y coordinates of the reconstructed and MC tracks
798 Double_t ddtx = fT[2] - PointPx/PointPz; // The difference between the tx coordinates of the reconstructed and MC tracks
799 Double_t ddty = fT[3] - PointPy/PointPz; // The difference between the ty coordinates of the reconstructed and MC tracks
800 Double_t ddqp = fT[4] - (qtrack/sqrt(PointPx*PointPx+
801 PointPy*PointPy+
802 PointPz*PointPz)); // The difference between the qp coordinates of the reconstructed and MC tracks
803
804 fprintf(fBack,"Event %i\n",iEvent);
805 fprintf(fBack,"PDG = %i\n",track_mc->GetPdgCode());
806 fprintf(fBack,"MC P = %f, Px = %f, Py = %f, Pz = %f",
807 sqrt(MCPoint->GetPz()*MCPoint->GetPz()+MCPoint->GetPy()*MCPoint->GetPy()+MCPoint->GetPx()*MCPoint->GetPx()),
808 MCPoint->GetPx(),MCPoint->GetPy(),MCPoint->GetPz());
809 fprintf(fBack," mc track reco track\n");
810 fprintf(fBack,"X %f %f\n",fT[0],MCFirstPoint->GetX());
811 fprintf(fBack,"Y %f %f\n",fT[1],MCFirstPoint->GetY());
812 fprintf(fBack,"Z %f %f\n",fT[5],MCFirstPoint->GetZ());
813 fprintf(fBack,"Tx %f %f\n",fT[2],PointPx/PointPz);
814 fprintf(fBack,"Ty %f %f\n",fT[3],PointPy/PointPz);
815 fprintf(fBack,"qp %f %f\n\n",fT[4],qtrack/sqrt(PointPx*PointPx+PointPy*PointPy+PointPz*PointPz));
816
817 //differences of the KF track and MC track parameters
818 res_AtFP_x -> Fill(ddx*10000.);
819 res_AtFP_y -> Fill(ddy*10000.);
820 res_AtFP_tx -> Fill(ddtx);
821 res_AtFP_ty -> Fill(ddty);
822 res_AtFP_qp -> Fill(ddqp);
823
824 //pulls of the parameters
825 pull_AtFP_x -> Fill(ddx/sqrt(fC[0]));
826 pull_AtFP_y -> Fill(ddy/sqrt(fC[2]));
827 pull_AtFP_tx -> Fill(ddtx/sqrt(fC[5]));
828 pull_AtFP_ty -> Fill(ddty/sqrt(fC[9]));
829 pull_AtFP_qp -> Fill(ddqp/sqrt(fC[14]));
830 }
831*/
832}
ClassImp
ClassImp(CbmKFTrackFitQa) CbmKFTrackFitQa
Definition CbmKFTrackFitQa.cxx:35
CbmStsTrack.h
Data class for STS tracks.
sqrt
friend fvec sqrt(const fvec &a)
Definition KfSimdPseudo.h:118
CbmCluster::GetMatch
CbmMatch * GetMatch() const
Definition CbmCluster.h:91
CbmHit::GetMatch
CbmMatch * GetMatch() const
Definition CbmHit.h:75
CbmKFTrErrMCPoints::MvdArray
std::vector< CbmMvdPoint * > MvdArray
Definition CbmKFTrErrMCPoints.h:68
CbmKFTrErrMCPoints::StsArray
std::vector< CbmStsPoint * > StsArray
Definition CbmKFTrErrMCPoints.h:67
CbmKFTrackFitQa::listMvdPts
TClonesArray * listMvdPts
Definition CbmKFTrackFitQa.h:60
CbmKFTrackFitQa::listStsTracks
TClonesArray * listStsTracks
Definition CbmKFTrackFitQa.h:63
CbmKFTrackFitQa::vStsHitMatch
std::vector< int > vStsHitMatch
Definition CbmKFTrackFitQa.h:75
CbmKFTrackFitQa::FillHistoAtParticleVertex
void FillHistoAtParticleVertex(CbmMCTrack *track_mc, CbmKFTrack *track_kf)
Definition CbmKFTrackFitQa.cxx:274
CbmKFTrackFitQa::Init
virtual InitStatus Init()
Definition CbmKFTrackFitQa.cxx:212
CbmKFTrackFitQa::listStsDigi
TClonesArray * listStsDigi
Definition CbmKFTrackFitQa.h:68
CbmKFTrackFitQa::listStsHits
TClonesArray * listStsHits
Definition CbmKFTrackFitQa.h:64
CbmKFTrackFitQa::FillHistoAtFirstPoint
void FillHistoAtFirstPoint(CbmKFTrErrMCPoints *mc_points, CbmMCTrack *track_mc, CbmKFTrack *track_kf)
Definition CbmKFTrackFitQa.cxx:349
CbmKFTrackFitQa::listStsPts
TClonesArray * listStsPts
Definition CbmKFTrackFitQa.h:59
CbmKFTrackFitQa::listMvdHitMatches
TClonesArray * listMvdHitMatches
Definition CbmKFTrackFitQa.h:66
CbmKFTrackFitQa::listStsTracksMatch
TClonesArray * listStsTracksMatch
Definition CbmKFTrackFitQa.h:62
CbmKFTrackFitQa::listStsClusters
TClonesArray * listStsClusters
Definition CbmKFTrackFitQa.h:67
CbmKFTrackFitQa::listMCTracks
TClonesArray * listMCTracks
Definition CbmKFTrackFitQa.h:61
CbmKFTrackFitQa::listMvdHits
TClonesArray * listMvdHits
Definition CbmKFTrackFitQa.h:65
CbmKFTrackFitQa::Exec
void Exec(Option_t *option)
Definition CbmKFTrackFitQa.cxx:239
CbmKFTrackFitQa::listStsDigiMatch
TClonesArray * listStsDigiMatch
Definition CbmKFTrackFitQa.h:69
CbmKFTrackFitQa::ReInit
virtual InitStatus ReInit()
Definition CbmKFTrackFitQa.cxx:210
CbmKFTrackInterface::GetTrack
virtual Double_t * GetTrack()
Is it electron.
Definition CbmKFTrackInterface.cxx:35
CbmKFTrackInterface::GetCovMatrix
virtual Double_t * GetCovMatrix()
array[6] of track parameters(x,y,tx,ty,qp,z)
Definition CbmKFTrackInterface.cxx:36
CbmKFTrackInterface::Extrapolate
Int_t Extrapolate(Double_t z, Double_t *QP0=nullptr)
Access to i-th hit.
Definition CbmKFTrackInterface.cxx:41
CbmKFTrack::GetCovMatrix
Double_t * GetCovMatrix() override
array[6] of track parameters(x,y,tx,ty,qp,z)
Definition CbmKFTrack.h:78
CbmKFTrack::GetTrack
Double_t * GetTrack() override
Is it electron.
Definition CbmKFTrack.h:77
CbmKF::Instance
static CbmKF * Instance()
Definition CbmKF.h:40
CbmMCTrack::GetPz
double GetPz() const
Definition CbmMCTrack.h:72
CbmMCTrack::GetPx
double GetPx() const
Definition CbmMCTrack.h:70
CbmMCTrack::GetStartZ
double GetStartZ() const
Definition CbmMCTrack.h:75
CbmMCTrack::GetStartX
double GetStartX() const
Definition CbmMCTrack.h:73
CbmMCTrack::GetPdgCode
int32_t GetPdgCode() const
Definition CbmMCTrack.h:68
CbmMCTrack::GetStartY
double GetStartY() const
Definition CbmMCTrack.h:74
CbmMCTrack::GetPy
double GetPy() const
Definition CbmMCTrack.h:71
CbmMatch::GetLink
const CbmLink & GetLink(int32_t i) const
Definition CbmMatch.h:39
CbmMatch::GetNofLinks
int32_t GetNofLinks() const
Definition CbmMatch.h:42
CbmStsCluster
Data class for STS clusters.
Definition CbmStsCluster.h:35
CbmStsHit
data class for a reconstructed 3-d hit in the STS
Definition CbmStsHit.h:35
CbmStsPoint::GetXOut
double GetXOut() const
Definition CbmStsPoint.h:76
CbmStsPoint::GetYIn
double GetYIn() const
Definition CbmStsPoint.h:74
CbmStsPoint::GetXIn
double GetXIn() const
Definition CbmStsPoint.h:73
CbmStsPoint::GetYOut
double GetYOut() const
Definition CbmStsPoint.h:77
CbmTrackMatch::GetMCTrackId
int32_t GetMCTrackId() const
Definition CbmTrackMatch.h:45
CbmTrackMatch::GetNofMCTracks
int32_t GetNofMCTracks() const
Definition CbmTrackMatch.h:57
std
Hash for CbmL1LinkKey.
Definition algo/base/Options.cxx:21