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