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