BmnRoot 26.10.0-18-g4e3414f9
Software framework of the BM@N experiment
Loading...
Searching...
No Matches
SrcInnerTrackingRun7.cxx
Go to the documentation of this file.
1#include "SrcInnerTrackingRun7.h"
2
3#include "BmnCellDuet.h"
4#include "BmnFieldMap.h"
5#include "BmnGemStripHit.h"
6#include "BmnKalmanFilter.h"
7#include "BmnMatch.h"
8#include "BmnMath.h"
9#include "FairRunAna.h"
10#include "FairTrackParam.h"
11#include "TStyle.h"
12
13#include <TCanvas.h>
14#include <TH2F.h>
15#include <TLine.h>
16#include <TMath.h>
17#include <map>
18#include <vector>
19#if defined(_OPENMP)
20#include "omp.h"
21
22#include <cstdlib>
23#endif
24
25//-----------------------------------------
26static Double_t workTime = 0.0;
27//-----------------------------------------
28
29using namespace std;
30using namespace TMath;
31
32SrcInnerTrackingRun7::SrcInnerTrackingRun7(Int_t run, Bool_t target)
33 : fSteerFile("")
34 , fSteering(nullptr)
35{
36 fEventNo = 0;
37 fIsTarget = target;
38 fGemHitsArray = nullptr;
39 fHitsArray = nullptr;
40 fGlobTracksArray = nullptr;
41 fGemTracksArray = nullptr;
42 fField = nullptr;
43 fGemHitsBranchName = "BmnGemStripHit";
44 fGlobTracksBranchName = "BmnGlobalTrack";
45 fGemTracksBranchName = "BmnGemTrack";
46 fGemDetector = nullptr;
47 fHitXCutMin = nullptr;
48 fHitXCutMax = nullptr;
49 fHitYCutMin = nullptr;
50 fHitYCutMax = nullptr;
51 fNHitsCut = 0;
52}
53
54SrcInnerTrackingRun7::~SrcInnerTrackingRun7()
55{
56 delete fKalman;
57 delete fGemDetector;
58 delete[] fHitXCutMin;
59 delete[] fHitXCutMax;
60 delete[] fHitYCutMin;
61 delete[] fHitYCutMax;
62}
63
64InitStatus SrcInnerTrackingRun7::Init()
65{
66 if (fVerbose > 1)
67 cout << "======================== GEM tracking init started ====================" << endl;
68
69 fKalman = new BmnKalmanFilter();
70
71 // Get ROOT Manager
72 FairRootManager* ioman = FairRootManager::Instance();
73 if (!ioman)
74 Fatal("Init", "FairRootManager is not instantiated");
75
76 // MC tracks for drawing during debugging
77 fMCTracksArray = (TClonesArray*)ioman->GetObject("MCTrack"); // in
78
79 // GEM
80 fGemPointsArray = (TClonesArray*)ioman->GetObject("StsPoint"); // in
81 fGemHitsArray = (TClonesArray*)ioman->GetObject(fGemHitsBranchName); // in
82 fGemTracksArray = new TClonesArray("BmnTrack", 100); // out
83 ioman->Register(fGemTracksBranchName, "GEM", fGemTracksArray, kTRUE);
84
85 if (!fGemHitsArray) {
86 cout << "SrcInnerTrackingRun7::Init(): branch " << fGemHitsBranchName << " not found! Task will be deactivated"
87 << endl;
88 SetActive(kFALSE);
89 return kERROR;
90 }
91
92 fGlobTracksArray = new TClonesArray("BmnGlobalTrack", 100); // out
93 ioman->Register(fGlobTracksBranchName, "GLOBAL", fGlobTracksArray, kTRUE);
94 fHitsArray = new TClonesArray("BmnHit", 100); // out
95 ioman->Register("BmnInnerHits", "HITS", fHitsArray, kTRUE);
96
97 fField = FairRunAna::Instance()->GetField();
98 if (!fField)
99 Fatal("Init", "No Magnetic Field found");
100
101 TString gPathConfig = gSystem->Getenv("VMCWORKDIR");
102 TString gPathGemConfig = gPathConfig + "/parameters/gem/XMLConfigs/GemRunSRCSpring2018.xml";
103 fGemDetector = new BmnGemStripStationSet(gPathGemConfig);
104
105 fNStations = fGemDetector->GetNStations() - 4; // first four stations are in arms
106 fHitXCutMin = new Double_t[fNStations];
107 fHitXCutMax = new Double_t[fNStations];
108 fHitYCutMin = new Double_t[fNStations];
109 fHitYCutMax = new Double_t[fNStations];
110
111 BmnFieldMap* field = (BmnFieldMap*)FairRunAna::Instance()->GetField();
112 Bool_t isField = !(field->IsFieldOff());
113 fSteering = new BmnSteering(isField ? "SRC_run7_withField.dat" : "SRC_run7_noField.dat");
114
115 if (fVerbose > 1)
116 cout << "======================== GEM tracking init finished ===================" << endl;
117
118 return kSUCCESS;
119}
120
121void SrcInnerTrackingRun7::Exec(Option_t* opt)
122{
123 if (fVerbose > 1)
124 cout << "\n======================== GEM tracking exec started ====================" << endl;
125 if (fVerbose > 1)
126 cout << "\n Event number: " << fEventNo << endl;
127
128 if (!IsActive())
129 return;
130
131 clock_t tStart = clock();
132
133 fGemTracksArray->Delete();
134 fGlobTracksArray->Delete();
135 fHitsArray->Delete();
136
137 fEventNo++;
138 Int_t nHitsCut = 200; // fSteering->GetNHitsCutTotal();
139
140 for (Int_t iHit = 0; iHit < fGemHitsArray->GetEntriesFast(); ++iHit) {
141 BmnGemStripHit hit = *((BmnGemStripHit*)fGemHitsArray->At(iHit));
142 if (hit.GetStation() < 4)
143 continue;
144 if (Abs(hit.GetStripTotalSignalInUpperLayer() - hit.GetStripTotalSignalInLowerLayer())
145 / (hit.GetStripTotalSignalInUpperLayer() + hit.GetStripTotalSignalInLowerLayer())
146 > 0.5)
147 continue;
148 hit.SetStation(hit.GetStation() - 4); // shift for correct station numbering
149 hit.SetIndex(iHit);
150 hit.SetDetId(kGEM);
151 hit.SetDxyz(0.5, 0.5, 0.5);
152 new ((*fHitsArray)[fHitsArray->GetEntriesFast()]) BmnHit(hit);
153 }
154
155 if (fHitsArray->GetEntriesFast() > nHitsCut || fHitsArray->GetEntriesFast() == 0)
156 return;
157
158 for (Int_t iStat = 0; iStat < fNStations; iStat++) {
159 fHitXCutMin[iStat] = fSteering->GetHitXCutMin()[iStat];
160 fHitXCutMax[iStat] = fSteering->GetHitXCutMax()[iStat];
161 fHitYCutMin[iStat] = fSteering->GetHitYCutMin()[iStat];
162 fHitYCutMax[iStat] = fSteering->GetHitYCutMax()[iStat];
163 }
164
165 fNHitsCut = fSteering->GetNHitsCut();
166 fChiSquareCut = fSteering->GetChiSquareCut();
167
168 fDistCut = 0.45;
169 fNHitsCut = 4;
170 FindTracks_6of6();
171 // FindTracks_5of6();
172 FindTracks_4of6();
173 // FindTracks_3of6();
174 // cout << "FindTracks_4of4_OnLastStations: " << fGlobTracksArray->GetEntriesFast() << " tracks" << endl;
175
176 // DrawHits();
177
178 clock_t tFinish = clock();
179 if (fVerbose > 0)
180 cout << "SrcInnerTrackingRun7: " << fGlobTracksArray->GetEntriesFast() << " tracks" << endl;
181
182 workTime += ((Double_t)(tFinish - tStart)) / CLOCKS_PER_SEC;
183
184 if (fVerbose > 1)
185 cout << "\n======================== GEM tracking exec finished ===================" << endl;
186 return;
187}
188
189BmnStatus SrcInnerTrackingRun7::FindTracks_6of6()
190{
191 vector<BmnHit*> sortedHits[fNStations];
192 SortHits(sortedHits);
193
194 vector<BmnTrack> candidates;
195 vector<BmnTrack> sortedCandidates;
196
197 for (BmnHit* hit5 : sortedHits[5]) {
198 for (BmnHit* hit4 : sortedHits[4]) {
199 if (Abs(hit5->GetX() - hit4->GetX()) > 5 || Abs(hit5->GetY() - hit4->GetY()) > 2)
200 continue;
201 for (BmnHit* hit3 : sortedHits[3]) {
202 if (Abs(hit4->GetX() - hit3->GetX()) > 5 || Abs(hit4->GetY() - hit3->GetY()) > 2)
203 continue;
204 for (BmnHit* hit2 : sortedHits[2]) {
205 if (Abs(hit3->GetX() - hit2->GetX()) > 5 || Abs(hit3->GetY() - hit2->GetY()) > 2)
206 continue;
207 for (BmnHit* hit1 : sortedHits[1]) {
208 if (Abs(hit2->GetX() - hit1->GetX()) > 5 || Abs(hit2->GetY() - hit1->GetY()) > 2)
209 continue;
210 for (BmnHit* hit0 : sortedHits[0]) {
211 if (Abs(hit1->GetX() - hit0->GetX()) > 5 || Abs(hit1->GetY() - hit0->GetY()) > 2)
212 continue;
213 BmnTrack cand;
214 cand.AddHit(hit5->GetRefIndex(), hit5);
215 cand.AddHit(hit4->GetRefIndex(), hit4);
216 cand.AddHit(hit3->GetRefIndex(), hit3);
217 cand.AddHit(hit2->GetRefIndex(), hit2);
218 cand.AddHit(hit1->GetRefIndex(), hit1);
219 cand.AddHit(hit0->GetRefIndex(), hit0);
220 cand.SortHits();
221 CalculateTrackParams(&cand);
222 cand.SetNDF(cand.GetNHits() * 2 - 5);
223 candidates.push_back(cand);
224 }
225 }
226 }
227 }
228 }
229 }
230
231 BmnFieldMap* field = (BmnFieldMap*)FairRunAna::Instance()->GetField();
232 Bool_t isField = !(field->IsFieldOff());
233 if (isField)
234 TrackUpdateByKalman(candidates);
235 else
236 TrackUpdateByLine(candidates);
237 SortTracks(candidates, sortedCandidates);
238 TrackSelection(sortedCandidates);
239
240 return kBMNSUCCESS;
241}
242
243BmnStatus SrcInnerTrackingRun7::FindTracks_3of6()
244{
245 vector<BmnHit*> sortedHits[fNStations];
246 SortHits(sortedHits);
247
248 vector<BmnTrack> candidates;
249 vector<BmnTrack> sortedCandidates;
250
251 FindCandidateByThreeStations(0, 1, 2, candidates, sortedHits);
252 FindCandidateByThreeStations(1, 2, 3, candidates, sortedHits);
253 FindCandidateByThreeStations(2, 3, 4, candidates, sortedHits);
254 FindCandidateByThreeStations(3, 4, 5, candidates, sortedHits);
255
256 BmnFieldMap* field = (BmnFieldMap*)FairRunAna::Instance()->GetField();
257 Bool_t isField = !(field->IsFieldOff());
258 if (isField)
259 TrackUpdateByKalman(candidates);
260 else
261 TrackUpdateByLine(candidates);
262 SortTracks(candidates, sortedCandidates);
263 TrackSelection(sortedCandidates);
264
265 return kBMNSUCCESS;
266}
267
268BmnStatus SrcInnerTrackingRun7::FindTracks_5of6()
269{
270 vector<BmnHit*> sortedHits[fNStations];
271 SortHits(sortedHits);
272
273 vector<BmnTrack> candidates;
274 vector<BmnTrack> sortedCandidates;
275
276 FindCandidateByFiveStations(0, 1, 2, 3, 4, candidates, sortedHits);
277 FindCandidateByFiveStations(0, 1, 2, 3, 5, candidates, sortedHits);
278 FindCandidateByFiveStations(0, 1, 2, 4, 5, candidates, sortedHits);
279 FindCandidateByFiveStations(0, 1, 3, 4, 5, candidates, sortedHits);
280 FindCandidateByFiveStations(0, 2, 3, 4, 5, candidates, sortedHits);
281 FindCandidateByFiveStations(1, 2, 3, 4, 5, candidates, sortedHits);
282
283 BmnFieldMap* field = (BmnFieldMap*)FairRunAna::Instance()->GetField();
284 Bool_t isField = !(field->IsFieldOff());
285 if (isField)
286 TrackUpdateByKalman(candidates);
287 else
288 TrackUpdateByLine(candidates);
289 SortTracks(candidates, sortedCandidates);
290 TrackSelection(sortedCandidates);
291
292 return kBMNSUCCESS;
293}
294
295BmnStatus SrcInnerTrackingRun7::FindTracks_4of6()
296{
297 vector<BmnHit*> sortedHits[fNStations];
298 SortHits(sortedHits);
299
300 vector<BmnTrack> candidates;
301 vector<BmnTrack> sortedCandidates;
302
303 FindCandidateByFourStations(0, 1, 2, 3, candidates, sortedHits);
304 FindCandidateByFourStations(1, 2, 3, 4, candidates, sortedHits);
305 FindCandidateByFourStations(2, 3, 4, 5, candidates, sortedHits);
306
307 FindCandidateByFourStations(1, 3, 4, 5, candidates, sortedHits);
308 FindCandidateByFourStations(0, 2, 3, 4, candidates, sortedHits);
309 FindCandidateByFourStations(1, 2, 3, 5, candidates, sortedHits);
310 FindCandidateByFourStations(0, 1, 2, 4, candidates, sortedHits);
311
312 FindCandidateByFourStations(0, 3, 4, 5, candidates, sortedHits);
313 FindCandidateByFourStations(0, 1, 2, 5, candidates, sortedHits);
314 FindCandidateByFourStations(1, 2, 4, 5, candidates, sortedHits);
315 FindCandidateByFourStations(0, 2, 4, 5, candidates, sortedHits);
316 FindCandidateByFourStations(0, 2, 3, 5, candidates, sortedHits);
317 FindCandidateByFourStations(0, 1, 4, 5, candidates, sortedHits);
318 FindCandidateByFourStations(0, 1, 3, 5, candidates, sortedHits);
319 FindCandidateByFourStations(0, 1, 3, 4, candidates, sortedHits);
320
321 BmnFieldMap* field = (BmnFieldMap*)FairRunAna::Instance()->GetField();
322 Bool_t isField = !(field->IsFieldOff());
323 if (isField)
324 TrackUpdateByKalman(candidates);
325 else
326 TrackUpdateByLine(candidates);
327 SortTracks(candidates, sortedCandidates);
328 TrackSelection(sortedCandidates);
329
330 return kBMNSUCCESS;
331}
332
333BmnStatus SrcInnerTrackingRun7::FindCandidateByFiveStations(Short_t st0,
334 Short_t st1,
335 Short_t st2,
336 Short_t st3,
337 Short_t st4,
338 vector<BmnTrack>& candidates,
339 vector<BmnHit*>* sortedHits)
340{
341 for (BmnHit* hit4 : sortedHits[st4]) {
342 for (BmnHit* hit3 : sortedHits[st3]) {
343 if (Abs(hit4->GetX() - hit3->GetX()) > 5 || Abs(hit4->GetY() - hit3->GetY()) > 2)
344 continue;
345 for (BmnHit* hit2 : sortedHits[st2]) {
346 if (Abs(hit3->GetX() - hit2->GetX()) > 5 || Abs(hit3->GetY() - hit2->GetY()) > 2)
347 continue;
348 for (BmnHit* hit1 : sortedHits[st1]) {
349 if (Abs(hit2->GetX() - hit1->GetX()) > 5 || Abs(hit2->GetY() - hit1->GetY()) > 2)
350 continue;
351 for (BmnHit* hit0 : sortedHits[st0]) {
352 if (Abs(hit1->GetX() - hit0->GetX()) > 5 || Abs(hit1->GetY() - hit0->GetY()) > 2)
353 continue;
354 BmnTrack cand;
355 cand.AddHit(hit4->GetRefIndex(), hit4);
356 cand.AddHit(hit3->GetRefIndex(), hit3);
357 cand.AddHit(hit2->GetRefIndex(), hit2);
358 cand.AddHit(hit1->GetRefIndex(), hit1);
359 cand.AddHit(hit0->GetRefIndex(), hit0);
360 cand.SortHits();
361 CalculateTrackParams(&cand);
362 candidates.push_back(cand);
363 }
364 }
365 }
366 }
367 }
368
369 return kBMNSUCCESS;
370}
371
372BmnStatus SrcInnerTrackingRun7::FindCandidateByFourStations(Short_t st0,
373 Short_t st1,
374 Short_t st2,
375 Short_t st3,
376 vector<BmnTrack>& candidates,
377 vector<BmnHit*>* sortedHits)
378{
379 for (BmnHit* hit3 : sortedHits[st3]) {
380 for (BmnHit* hit2 : sortedHits[st2]) {
381 if (Abs(hit3->GetX() - hit2->GetX()) > 5 || Abs(hit3->GetY() - hit2->GetY()) > 2)
382 continue;
383 for (BmnHit* hit1 : sortedHits[st1]) {
384 if (Abs(hit2->GetX() - hit1->GetX()) > 5 || Abs(hit2->GetY() - hit1->GetY()) > 2)
385 continue;
386 for (BmnHit* hit0 : sortedHits[st0]) {
387 if (Abs(hit1->GetX() - hit0->GetX()) > 5 || Abs(hit1->GetY() - hit0->GetY()) > 2)
388 continue;
389 BmnTrack cand;
390 cand.AddHit(hit3->GetRefIndex(), hit3);
391 cand.AddHit(hit2->GetRefIndex(), hit2);
392 cand.AddHit(hit1->GetRefIndex(), hit1);
393 cand.AddHit(hit0->GetRefIndex(), hit0);
394 cand.SortHits();
395 CalculateTrackParams(&cand);
396 candidates.push_back(cand);
397 }
398 }
399 }
400 }
401 vector<Short_t> restSt(2);
402 for (Int_t i = 0; i < fNStations; ++i) {
403 if (i == st0 || i == st1 || i == st2 || i == st3)
404 continue;
405 restSt.push_back(i);
406 }
407 for (BmnTrack& cand : candidates) {
408 for (Int_t i = 0; i < 2; ++i) {
409 MatchHit(&cand, sortedHits[restSt[i]], restSt[i] > st3);
410 }
411
412 // NDF = (N counts in ZX plane + N counts in ZY plane) - 2 parameters of Line in ZY plane - 3 parameters of
413 // Circle in ZX plane Check it!!!
414 cand.SetNDF(cand.GetNHits() * 2 - 5);
415 }
416
417 return kBMNSUCCESS;
418}
419
420BmnStatus SrcInnerTrackingRun7::FindCandidateByThreeStations(Short_t st0,
421 Short_t st1,
422 Short_t st2,
423 vector<BmnTrack>& candidates,
424 vector<BmnHit*>* sortedHits)
425{
426 for (BmnHit* hit2 : sortedHits[st2]) {
427 for (BmnHit* hit1 : sortedHits[st1]) {
428 if (Abs(hit2->GetX() - hit1->GetX()) > 10 || Abs(hit2->GetY() - hit1->GetY()) > 5)
429 continue;
430 for (BmnHit* hit0 : sortedHits[st0]) {
431 if (Abs(hit1->GetX() - hit0->GetX()) > 10 || Abs(hit1->GetY() - hit0->GetY()) > 5)
432 continue;
433 BmnTrack cand;
434 cand.AddHit(hit2->GetRefIndex(), hit2);
435 cand.AddHit(hit1->GetRefIndex(), hit1);
436 cand.AddHit(hit0->GetRefIndex(), hit0);
437 cand.SortHits();
438 CalculateTrackParams(&cand);
439 candidates.push_back(cand);
440 }
441 }
442 }
443 vector<Short_t> restSt(3);
444 for (Int_t i = 0; i < fNStations; ++i) {
445 if (i == st0 || i == st1 || i == st2)
446 continue;
447 restSt.push_back(i);
448 }
449 for (BmnTrack& cand : candidates) {
450 for (Int_t i = 0; i < 3; ++i) {
451 MatchHit(&cand, sortedHits[restSt[i]], restSt[i] > st2);
452 }
453
454 // NDF = (N counts in ZX plane + N counts in ZY plane) - 2 parameters of Line in ZY plane - 3 parameters of
455 // Circle in ZX plane Check it!!!
456 cand.SetNDF(cand.GetNHits() * 2 - 5);
457 }
458
459 return kBMNSUCCESS;
460}
461
462BmnStatus SrcInnerTrackingRun7::MatchHit(BmnTrack* cand, vector<BmnHit*> sortedHits, Bool_t downstream)
463{
464 // Double_t cutX = 0.45;
465 // Double_t cutY = 0.45;
466 if (sortedHits.size() == 0)
467 return kBMNERROR;
468 FairTrackParam par = downstream ? *(cand->GetParamLast()) : *(cand->GetParamFirst());
469 Double_t minDist = DBL_MAX;
470 // Double_t minDX = DBL_MAX;
471 // Double_t minDY = DBL_MAX;
472 BmnHit* minHit = nullptr;
473 fKalman->TGeoTrackPropagate(&par, sortedHits.at(0)->GetZ(), 2212, nullptr, nullptr);
474 for (BmnHit* hit : sortedHits) {
475 // Double_t dX = Abs(par.GetX() - hit->GetX());
476 // Double_t dY = Abs(par.GetY() - hit->GetY());
477 Double_t d = Sqrt(Sq(par.GetX() - hit->GetX()) + Sq(par.GetY() - hit->GetY()));
478 // if (dX < cutX && dY < cutY && dX < minDX && dY < minDY) {
479 if (d < fDistCut && d < minDist) {
480 // minDX = dX;
481 // minDY = dY;
482 minDist = d;
483 minHit = hit;
484 }
485 }
486 if (minHit) {
487 cand->AddHit(minHit->GetRefIndex(), minHit);
488 cand->SortHits();
489 Double_t chi2 = 0.0;
490 fKalman->Update(&par, minHit, chi2);
491 cand->SetChi2(cand->GetChi2() + chi2);
492 }
493 if (downstream)
494 cand->SetParamLast(par);
495 else
496 cand->SetParamFirst(par);
497
498 return kBMNSUCCESS;
499}
500
501BmnStatus SrcInnerTrackingRun7::SortTracks(vector<BmnTrack>& inTracks, vector<BmnTrack>& sortedTracks)
502{
503 const Int_t n = fNStations - fNHitsCut + 1; // 6 for geometry 2018 (4, 5, 6, 7, 8, 9)
504 multimap<Float_t, Int_t> sortedMap[n]; // array of map<Chi2,trIdx>. Each element of array corresponds fixed number
505 // of hits on track (4, 5, 6)
506 for (size_t iTr = 0; iTr < inTracks.size(); ++iTr) {
507 if (inTracks.at(iTr).GetNHits() < fNHitsCut)
508 continue;
509 if (inTracks.at(iTr).GetChi2() / inTracks.at(iTr).GetNDF() > fChiSquareCut)
510 continue;
511 sortedMap[inTracks.at(iTr).GetNHits() - fNHitsCut].insert(
512 pair<Float_t, Int_t>(inTracks.at(iTr).GetChi2() / inTracks.at(iTr).GetNDF(), iTr));
513 }
514
515 for (Int_t i = n - 1; i >= 0; i--) {
516 for (auto it : sortedMap[i])
517 sortedTracks.push_back(inTracks.at(it.second));
518 sortedMap[i].clear();
519 }
520 return kBMNSUCCESS;
521}
522
523void SrcInnerTrackingRun7::Finish()
524{
525 ofstream outFile;
526 outFile.open("QA/timing.txt");
527 outFile << "Track Finder Time: " << workTime << endl;
528 cout << "Work time of the GEM tracking: " << workTime << endl;
529 return;
530}
531
532BmnStatus SrcInnerTrackingRun7::SortHits(vector<BmnHit*>* sortedHits)
533{
534 for (Int_t iHit = 0; iHit < fHitsArray->GetEntriesFast(); ++iHit) {
535 BmnHit* hit = (BmnHit*)fHitsArray->At(iHit);
536 if (!hit)
537 continue;
538 if (hit->IsUsed())
539 continue;
540 hit->SetRefIndex(iHit);
541 Int_t station = hit->GetStation();
542 if (station < 0)
543 continue;
544 // if (hit->GetX() > fHitXCutMax[station] || hit->GetX() < fHitXCutMin[station]) continue;
545 // if (hit->GetY() > fHitYCutMax[station] || hit->GetY() < fHitYCutMin[station]) continue;
546 sortedHits[station].push_back(hit);
547 }
548
549 return kBMNSUCCESS;
550}
551
552BmnStatus SrcInnerTrackingRun7::TrackUpdateByLine(vector<BmnTrack>& cands)
553{
554 for (BmnTrack& cand : cands) {
555 Double_t Tx = LineFit((BmnTrack*)&cand, fHitsArray, "ZX").X();
556 Double_t chiX = LineFit((BmnTrack*)&cand, fHitsArray, "ZX").Z();
557 Double_t Ty = LineFit((BmnTrack*)&cand, fHitsArray, "ZY").X();
558 Double_t chiY = LineFit((BmnTrack*)&cand, fHitsArray, "ZY").Z();
559
560 cand.SetChi2((chiX - chiY) > 0. ? chiX : chiY);
561
562 cand.GetParamFirst()->SetTx(Tx);
563 cand.GetParamFirst()->SetTy(Ty);
564
565 cand.GetParamLast()->SetTx(Tx);
566 cand.GetParamLast()->SetTy(Ty);
567 }
568
569 return kBMNSUCCESS;
570}
571
572BmnStatus SrcInnerTrackingRun7::TrackUpdateByKalman(vector<BmnTrack>& cands)
573{
574 for (BmnTrack& cand : cands) {
575 if (cand.GetNHits() < fNHitsCut)
576 continue;
577 FairTrackParam par = *(cand.GetParamFirst());
578 // par.SetQp(1.0/8.0);
579 Double_t chiTot = 0.0;
580 for (Int_t iHit = 0; iHit < cand.GetNHits(); ++iHit) {
581 BmnHit* hit = (BmnHit*)fHitsArray->At(cand.GetHitIndex(iHit));
582 Double_t chi = 0.0;
583 fKalman->TGeoTrackPropagate(&par, hit->GetZ(), 2212, nullptr, nullptr);
584 fKalman->Update(&par, hit, chi);
585 }
586 cand.SetParamLast(par);
587 for (Int_t iHit = cand.GetNHits() - 1; iHit >= 0; iHit--) {
588 BmnHit* hit = (BmnHit*)fHitsArray->At(cand.GetHitIndex(iHit));
589 Double_t chi = 0.0;
590 fKalman->TGeoTrackPropagate(&par, hit->GetZ(), 2212, nullptr, nullptr);
591 fKalman->Update(&par, hit, chi);
592 chiTot += chi;
593 }
594 cand.SetParamFirst(par);
595 cand.SetChi2(chiTot);
596 }
597
598 return kBMNSUCCESS;
599}
600
601BmnStatus SrcInnerTrackingRun7::TrackSelection(vector<BmnTrack>& sortedTracks)
602{
603 CheckSharedHits(sortedTracks);
604 for (size_t iTr = 0; iTr < sortedTracks.size(); ++iTr) {
605 BmnTrack tr = sortedTracks[iTr];
606 if (tr.GetFlag() == 666 || !IsParCorrect(tr.GetParamFirst()) || !IsParCorrect(tr.GetParamLast()))
607 continue;
608 BmnTrack gemTr;
609 BmnGlobalTrack globTr;
610
611 for (Int_t iHit = 0; iHit < tr.GetNHits(); ++iHit) {
612 BmnHit* hit = (BmnHit*)fHitsArray->At(tr.GetHitIndex(iHit));
613 globTr.AddHit(hit->GetRefIndex(), hit);
614 DetectorId detId = hit->GetDetId();
615 if (detId == kGEM)
616 gemTr.AddHit(hit->GetIndex(), hit);
617 else
618 continue;
619 }
620 globTr.SortHits();
621 gemTr.SortHits();
622
623 if (gemTr.GetNHits() > 0) {
624 gemTr.SetParamFirst(*(tr.GetParamFirst()));
625 gemTr.SetParamLast(*(tr.GetParamLast()));
626 new ((*fGemTracksArray)[fGemTracksArray->GetEntriesFast()]) BmnTrack(gemTr);
627 globTr.SetGemTrackIndex(fGemTracksArray->GetEntriesFast() - 1);
628 }
629 globTr.SetLength(CalculateLength(&tr));
630 globTr.SetParamFirst(*(tr.GetParamFirst()));
631 globTr.SetParamLast(*(tr.GetParamLast()));
632 globTr.SetNHits(tr.GetNHits());
633 globTr.SetNDF(tr.GetNDF());
634 globTr.SetChi2(tr.GetChi2());
635 new ((*fGlobTracksArray)[fGlobTracksArray->GetEntriesFast()]) BmnGlobalTrack(globTr);
636 SetHitsUsing(&tr, kTRUE);
637 if (!fIsTarget && iTr == 0)
638 break;
639 }
640
641 return kBMNSUCCESS;
642}
643
644void SrcInnerTrackingRun7::SetHitsUsing(BmnTrack* tr, Bool_t use)
645{
646 for (Int_t i = 0; i < tr->GetNHits(); ++i) {
647 BmnHit* hit = (BmnHit*)fHitsArray->At(tr->GetHitIndex(i));
648 if (!hit)
649 continue;
650 hit->SetUsing(use);
651 // for (Int_t iLink = 0; iLink < hit->GetDigitNumberMatch().GetNofLinks(); ++iLink) {
652 // BmnLink link = hit->GetDigitNumberMatch().GetLink(iLink);
653 // Int_t idxI = link.GetIndex();
654 // for (Int_t j = 0; j < fHitsArray->GetEntriesFast(); ++j) {
655 // BmnHit* hitJ = (BmnHit*)fHitsArray->At(j);
656 // if (!hitJ) continue;
657 // if (hitJ->IsUsed()) continue;
658 // for (Int_t iLink = 0; iLink < hit->GetDigitNumberMatch().GetNofLinks(); ++iLink) {
659 // BmnLink linkJ = hit->GetDigitNumberMatch().GetLink(iLink);
660 // Int_t idxJ = linkJ.GetIndex();
661 // if (idxI == idxJ) {
662 // hitJ->SetUsing(use);
663 // break;
664 // }
665 // }
666 // }
667 // }
668 }
669 return;
670}
671
672BmnStatus SrcInnerTrackingRun7::CalcCovMatrix(BmnTrack* tr)
673{
674 const UInt_t nHits = tr->GetNHits();
675 Double_t chi2circ = 0.0;
676 TVector3 CircParZX = CircleFit(tr, fHitsArray, chi2circ);
677
678 Double_t Xc = CircParZX.Y(); // x-coordinate of fit-circle center
679 Double_t Zc = CircParZX.X(); // z-coordinate of fit-circle center
680 fField = FairRunAna::Instance()->GetField();
681 const Double_t B = (LineFit(tr, fHitsArray, "ZY")).X(); // angle coefficient for helicoid
682
683 Double_t Q = (Xc > 0) ? +1 : -1;
684
685 Double_t sumX(0.0), sumY(0.0), sumTx(0.0), sumTy(0.0), sumQP(0.0);
686 Double_t sumXX(0.0), sumXY(0.0), sumXTx(0.0), sumXTy(0.0), sumXQP(0.0);
687 Double_t sumYY(0.0), sumYTx(0.0), sumYTy(0.0), sumYQP(0.0);
688 Double_t sumTxTx(0.0), sumTxTy(0.0), sumTxQP(0.0);
689 Double_t sumTyTy(0.0), sumTyQP(0.0);
690 Double_t sumQPQP(0.0);
691
692 for (UInt_t i = 0; i < nHits; ++i) {
693 BmnHit* hit = (BmnHit*)fHitsArray->At(tr->GetHitIndex(i));
694 if (!hit)
695 continue;
696 Double_t Xi = hit->GetX();
697 Double_t Yi = hit->GetY();
698 Double_t Zi = hit->GetZ();
699
700 Double_t Txi = -1.0 * (Zi - Zc) / (Xi - Xc);
701 Double_t Tyi = B;
702 Double_t Ri = Sqrt(Sq(Xi - Xc) + Sq(Zi - Zc));
703 Double_t QPi = Q / (0.0003 * Abs(fField->GetBy(Xi, Yi, Zi)) * Ri);
704
705 sumX += Xi;
706 sumY += Yi;
707 sumTx += Txi;
708 sumTy += Tyi;
709 sumQP += QPi;
710 sumXX += Xi * Xi;
711 sumXY += Xi * Yi;
712 sumXTx += Xi * Txi;
713 sumXTy += Xi * Tyi;
714 sumXQP += Xi * QPi;
715 sumYY += Yi * Yi;
716 sumYTx += Yi * Txi;
717 sumYTy += Yi * Tyi;
718 sumYQP += Yi * QPi;
719 sumTxTx += Txi * Txi;
720 sumTxTy += Txi * Tyi;
721 sumTxQP += Txi * QPi;
722 sumTyTy += Tyi * Tyi;
723 sumTyQP += Tyi * QPi;
724 sumQPQP += QPi * QPi;
725 }
726
727 Double_t meanX = sumX / nHits;
728 Double_t meanY = sumY / nHits;
729 Double_t meanTx = sumTx / nHits;
730 Double_t meanTy = sumTy / nHits;
731 Double_t meanQP = sumQP / nHits;
732
733 Double_t Cov_X_X = sumXX / nHits - Sq(meanX);
734 Double_t Cov_X_Y = sumXY / nHits - meanX * meanY;
735 Double_t Cov_X_Tx = sumXTx / nHits - meanX * meanTx;
736 Double_t Cov_X_Ty = sumXTy / nHits - meanX * meanTy;
737 Double_t Cov_X_Qp = sumXQP / nHits - meanX * meanQP;
738 Double_t Cov_Y_Y = sumYY / nHits - Sq(meanY);
739 Double_t Cov_Y_Tx = sumYTx / nHits - meanY * meanTx;
740 Double_t Cov_Y_Ty = sumYTy / nHits - meanY * meanTy;
741 Double_t Cov_Y_Qp = sumYQP / nHits - meanY * meanQP;
742 Double_t Cov_Tx_Tx = sumTxTx / nHits - Sq(meanTx);
743 Double_t Cov_Tx_Ty = sumTxTy / nHits - meanTx * meanTy;
744 Double_t Cov_Tx_Qp = sumTxQP / nHits - meanTx * meanQP;
745 Double_t Cov_Ty_Ty = sumTyTy / nHits - Sq(meanTy);
746 Double_t Cov_Ty_Qp = sumTyQP / nHits - meanTy * meanQP;
747 Double_t Cov_Qp_Qp = sumQPQP / nHits - Sq(meanQP);
748
749 FairTrackParam par;
750 par.SetCovariance(0, 0, Cov_X_X);
751 par.SetCovariance(0, 1, Cov_X_Y);
752 par.SetCovariance(0, 2, Cov_X_Tx);
753 par.SetCovariance(0, 3, Cov_X_Ty);
754 par.SetCovariance(0, 4, Cov_X_Qp);
755 par.SetCovariance(1, 1, Cov_Y_Y);
756 par.SetCovariance(1, 2, Cov_Y_Tx);
757 par.SetCovariance(1, 3, Cov_Y_Ty);
758 par.SetCovariance(1, 4, Cov_Y_Qp);
759 par.SetCovariance(2, 2, Cov_Tx_Tx);
760 par.SetCovariance(2, 3, Cov_Tx_Ty);
761 par.SetCovariance(2, 4, Cov_Tx_Qp);
762 par.SetCovariance(3, 3, Cov_Ty_Ty);
763 par.SetCovariance(3, 4, Cov_Ty_Qp);
764 par.SetCovariance(4, 4, Cov_Qp_Qp);
765
766 tr->SetParamFirst(par);
767 tr->SetParamLast(par);
768
769 return kBMNSUCCESS;
770}
771
772BmnStatus SrcInnerTrackingRun7::CalculateTrackParams(BmnTrack* tr)
773{
774
775 BmnFieldMap* field = (BmnFieldMap*)FairRunAna::Instance()->GetField();
776 Bool_t isField = !(field->IsFieldOff());
777
778 // Estimation of track parameters for events with magnetic field
779 const UInt_t nHits = tr->GetNHits();
780 if (nHits < 3)
781 return kBMNERROR;
782 TVector3 lineParZY = LineFit(tr, fHitsArray, "ZY");
783 if (lineParZY.Z() > 1)
784 return kBMNERROR; // cout << "chi2_lineFit = " << lineParZY.Z() << endl;
785 // tr->SetNDF(nHits - (isField ? 3 : 2));
786 const Double_t B = lineParZY.X(); // angle coefficient for helicoid
787
788 Double_t Tx_first = CalcTx((BmnHit*)fHitsArray->At(tr->GetHitIndex(0)), (BmnHit*)fHitsArray->At(tr->GetHitIndex(1)),
789 (BmnHit*)fHitsArray->At(tr->GetHitIndex(2)));
790 Double_t Tx_last =
791 CalcTx((BmnHit*)fHitsArray->At(tr->GetHitIndex(nHits - 1)), (BmnHit*)fHitsArray->At(tr->GetHitIndex(nHits - 2)),
792 (BmnHit*)fHitsArray->At(tr->GetHitIndex(nHits - 3)));
793
794 if (isField)
795 CalcCovMatrix(tr);
796 TVector3 firstPos;
797 TVector3 lastPos;
798 ((BmnHit*)fHitsArray->At(tr->GetHitIndex(0)))->Position(firstPos);
799 ((BmnHit*)fHitsArray->At(tr->GetHitIndex(nHits - 1)))->Position(lastPos);
800 tr->GetParamFirst()->SetPosition(firstPos);
801 tr->GetParamFirst()->SetTx(Tx_first);
802 tr->GetParamFirst()->SetTy(B);
803 tr->GetParamLast()->SetPosition(lastPos);
804 tr->GetParamLast()->SetTx(Tx_last);
805 tr->GetParamLast()->SetTy(B);
806 Bool_t doSimple = (nHits == 3) ? kTRUE : kFALSE;
807 Double_t QP = isField ? CalcQp(tr, doSimple) : 0.0;
808 tr->GetParamFirst()->SetQp(QP);
809 tr->GetParamLast()->SetQp(QP);
810 return kBMNSUCCESS;
811}
812
813TVector2 SrcInnerTrackingRun7::CalcMeanSigma(vector<Double_t> QpSegm)
814{
815 Double_t QpSum = 0.;
816 for (size_t iSegm = 0; iSegm < QpSegm.size(); iSegm++)
817 QpSum += QpSegm[iSegm];
818
819 Double_t QpMean = QpSum / QpSegm.size();
820
821 Double_t sqSigmaSum = 0.;
822 for (size_t iSegm = 0; iSegm < QpSegm.size(); iSegm++)
823 sqSigmaSum += Sq(QpSegm[iSegm] - QpMean);
824
825 return TVector2(QpMean, Sqrt(sqSigmaSum / QpSegm.size()));
826}
827
828Double_t SrcInnerTrackingRun7::CalcQp(BmnTrack* track, Bool_t doSimple)
829{
830 // Bool_t fRobustRefit = kTRUE;
831 // Bool_t fSimpleRefit = kFALSE;
832 vector<BmnHit*> hits;
833
834 for (Int_t iHit = 0; iHit < track->GetNHits(); iHit++)
835 hits.push_back((BmnHit*)fHitsArray->At(track->GetHitIndex(iHit)));
836
837 Int_t kNSegm = track->GetNHits() - 2;
838
839 Double_t QpRefit = 0.;
840 vector<Double_t> QpSegmBefore;
841
842 // Get q/p info from all track segments
843 for (Int_t iHit = 0; iHit < kNSegm; iHit++) {
844 BmnHit* first = hits[iHit];
845 BmnHit* second = hits[iHit + 1];
846 BmnHit* third = hits[iHit + 2];
847
848 TVector3 CircParZX = CircleBy3Hit(first, second, third);
849 Double_t R = CircParZX.Z();
850 Double_t Xc = CircParZX.Y();
851 Double_t Zc = CircParZX.X();
852
853 Double_t Q = (Xc > 0) ? +1. : -1.;
854 Double_t S = 0.0003
855 * (Abs(fField->GetBy(third->GetX(), third->GetY(), third->GetZ()))
856 + Abs(fField->GetBy(second->GetX(), second->GetY(), second->GetZ()))
857 + Abs(fField->GetBy(first->GetX(), first->GetY(), first->GetZ())))
858 / 3.;
859
860 Double_t Pt = S * R; // actually Pt/Q, but it doesn't matter
861 Double_t fX = first->GetX();
862 Double_t fZ = first->GetZ();
863
864 Double_t h = -1.0;
865
866 Double_t Tx_first = h * (fZ - Zc) / (fX - Xc);
867 TVector3 lineParZY = LineFit(track, fHitsArray, "ZY");
868 const Double_t B = lineParZY.X(); // angle coefficient for helicoid
869 Double_t Ty_first = B; // / (fX - Xc);
870
871 Double_t Pz = Pt / Sqrt(1 + Sq(Tx_first));
872 Double_t Px = Tx_first * Pz;
873 Double_t Py = Ty_first * Pz;
874 Double_t P = Sqrt(Sq(Px) + Sq(Py) + Sq(Pz));
875 Double_t QP = Q / P;
876
877 QpSegmBefore.push_back(QP);
878 }
879
880 // Non-robust (simple) refit when segments with bad q/p are not taken into account
881 if (doSimple) {
882 vector<Double_t> QpSegmAfter;
883 while (kTRUE) {
884 TVector2 meanSig = CalcMeanSigma(QpSegmBefore);
885 Double_t mean = meanSig.X();
886 Double_t sigma = meanSig.Y();
887 if (std::isnan(sigma) && fVerbose == 3) {
888 cout << "Bad refit convergence for track segment!!" << endl;
889 return kBMNERROR;
890 }
891
892 for (size_t iSegm = 0; iSegm < QpSegmBefore.size(); iSegm++)
893 if (Abs(QpSegmBefore[iSegm] - mean) - sigma <= 0.001) // Топорное сравнение FIXME
894 QpSegmAfter.push_back(QpSegmBefore[iSegm]);
895
896 if (QpSegmAfter.size() == QpSegmBefore.size()) {
897 QpRefit = mean;
898 break;
899 } else {
900 QpSegmBefore.clear();
901 QpSegmBefore.resize(0);
902
903 for (size_t iSegm = 0; iSegm < QpSegmAfter.size(); iSegm++)
904 QpSegmBefore.push_back(QpSegmAfter[iSegm]);
905
906 QpSegmAfter.clear();
907 QpSegmAfter.resize(0);
908 }
909 }
910 }
911
912 // Robust refit with use of Tukey weights calculation algorithm
913 // if (fRobustRefit) {
914 else
915 {
916 for (size_t iEle = 0; iEle < QpSegmBefore.size(); iEle++)
917 QpRefit += QpSegmBefore[iEle];
918
919 QpRefit /= QpSegmBefore.size();
920
921 vector<Double_t> d = dist(QpSegmBefore, QpRefit);
922
923 Double_t sigma = 0.;
924 for (size_t i = 0; i < QpSegmBefore.size(); i++)
925 sigma += (QpSegmBefore[i] - QpRefit) * (QpSegmBefore[i] - QpRefit);
926 sigma = Sqrt(sigma / QpSegmBefore.size());
927
928 vector<Double_t> w = W(d, sigma);
929 sigma = Sigma(d, w);
930
931 const Int_t kNIter = 20; // FIXME
932 for (Int_t iIter = 1; iIter < kNIter; iIter++) {
933 QpRefit = Mu(QpSegmBefore, w);
934 d = dist(QpSegmBefore, QpRefit);
935 w = W(d, sigma);
936 sigma = Sigma(d, w);
937 }
938 }
939
940 return QpRefit;
941}
942
943Double_t SrcInnerTrackingRun7::CalculateLength(BmnTrack* tr)
944{
945 if (!tr)
946 return 0.0;
947
948 vector<Double_t> X, Y, Z;
949 for (Int_t iGem = 0; iGem < tr->GetNHits(); iGem++) {
950 BmnHit* hit = (BmnHit*)fHitsArray->At(tr->GetHitIndex(iGem));
951 if (!hit)
952 continue;
953 X.push_back(hit->GetX());
954 Y.push_back(hit->GetY());
955 Z.push_back(hit->GetZ());
956 }
957 // Calculate distances between hits
958 Double_t length = 0.;
959 for (size_t i = 0; i < X.size() - 1; i++) {
960 Double_t dX = X[i] - X[i + 1];
961 Double_t dY = Y[i] - Y[i + 1];
962 Double_t dZ = Z[i] - Z[i + 1];
963 length += Sqrt(dX * dX + dY * dY + dZ * dZ);
964 }
965 tr->SetLength(length);
966 return length;
967}
968
969BmnStatus SrcInnerTrackingRun7::CheckSharedHits(vector<BmnTrack>& sortedTracks)
970{
971 set<Int_t> hitsId;
972
973 const Int_t kNSharedHits = 0; // fSteering->GetNSharedHits();
974
975 for (size_t iTr = 0; iTr < sortedTracks.size(); ++iTr) {
976 BmnTrack* tr = &(sortedTracks.at(iTr));
977 if (tr->GetFlag() == 666)
978 continue;
979
980 Int_t nofSharedHits = 0;
981 Int_t nofHits = tr->GetNHits();
982 for (Int_t iHit = 0; iHit < nofHits; iHit++)
983 if (hitsId.find(tr->GetHitIndex(iHit)) != hitsId.end()) {
984 nofSharedHits++;
985 if (nofSharedHits > kNSharedHits) {
986 tr->SetFlag(666);
987 break;
988 }
989 }
990 if (tr->GetFlag() == 666)
991 continue;
992
993 for (Int_t iHit = 0; iHit < nofHits; iHit++)
994 hitsId.insert(tr->GetHitIndex(iHit));
995 }
996 hitsId.clear();
997
998 return kBMNSUCCESS;
999}
1000
1001BmnStatus SrcInnerTrackingRun7::DrawHits()
1002{
1003 TH2F* h_HitsZX = new TH2F("h_HitsZX", "h_HitsZX", 400, 0.0, 200.0, 400, -100.0, 100.0);
1004 TH2F* h_HitsZY = new TH2F("h_HitsZY", "h_HitsZY", 400, 0.0, 200.0, 400, -10.0, 50.0);
1005 for (Int_t i = 0; i < fHitsArray->GetEntriesFast(); ++i) {
1006 BmnHit* hit = (BmnHit*)fHitsArray->At(i);
1007 h_HitsZX->Fill(hit->GetZ(), hit->GetX());
1008 h_HitsZY->Fill(hit->GetZ(), hit->GetY());
1009 }
1010
1011 TCanvas* c = new TCanvas("c", "c", 1000, 1000);
1012 c->Divide(1, 2);
1013 c->cd(1);
1014 h_HitsZX->SetMarkerStyle(8);
1015 h_HitsZX->SetMarkerSize(0.7);
1016 h_HitsZX->SetMarkerColor(kRed);
1017 h_HitsZX->Draw("P");
1018
1019 c->cd(2);
1020 h_HitsZY->SetMarkerStyle(8);
1021 h_HitsZY->SetMarkerSize(0.7);
1022 h_HitsZY->SetMarkerColor(kRed);
1023 h_HitsZY->Draw("P");
1024
1025 // if (fGlobTracksArray->GetEntriesFast() <= 0) return kBMNERROR;
1026
1027 for (Int_t iTr = 0; iTr < fGlobTracksArray->GetEntriesFast(); ++iTr) {
1028 BmnGlobalTrack* glTrack = (BmnGlobalTrack*)fGlobTracksArray->At(iTr);
1029
1030 BmnTrack* gemTr = nullptr;
1031 // BmnTrack* silTr = nullptr;
1032
1033 Double_t xPrev;
1034 Double_t yPrev;
1035 Double_t zPrev;
1036
1037 if (glTrack->GetGemTrackIndex() != -1) {
1038 gemTr = (BmnTrack*)fGemTracksArray->UncheckedAt(glTrack->GetGemTrackIndex());
1039
1040 Int_t startIdx = 0;
1041
1042 if (glTrack->GetSilTrackIndex() == -1) {
1043 xPrev = ((BmnHit*)fGemHitsArray->At(gemTr->GetHitIndex(0)))->GetX();
1044 yPrev = ((BmnHit*)fGemHitsArray->At(gemTr->GetHitIndex(0)))->GetY();
1045 zPrev = ((BmnHit*)fGemHitsArray->At(gemTr->GetHitIndex(0)))->GetZ();
1046 startIdx = 1;
1047 }
1048
1049 for (Int_t iHit = startIdx; iHit < gemTr->GetNHits(); iHit++) {
1050 BmnHit* hit = (BmnHit*)fGemHitsArray->UncheckedAt(gemTr->GetHitIndex(iHit));
1051 Double_t x = hit->GetX();
1052 Double_t y = hit->GetY();
1053 Double_t z = hit->GetZ();
1054 TLine* lineZX = new TLine(z, x, zPrev, xPrev);
1055 lineZX->SetLineColor(kRed);
1056 lineZX->SetLineWidth(2);
1057 TLine* lineZY = new TLine(z, y, zPrev, yPrev);
1058 lineZY->SetLineColor(kRed);
1059 lineZY->SetLineWidth(2);
1060 c->cd(1);
1061 lineZX->Draw("same");
1062 c->cd(2);
1063 lineZY->Draw("same");
1064 zPrev = z;
1065 yPrev = y;
1066 xPrev = x;
1067 }
1068 }
1069 }
1070
1071 // Int_t nMCtracks = 0;
1072 // for (Int_t iTr = 0; iTr < fMCTracksArray->GetEntriesFast(); ++iTr) {
1073 // CbmMCTrack* tr = (CbmMCTrack*)fMCTracksArray->At(iTr);
1074 // Int_t nPoints = 0; //tr->GetNPoints(kGEM);
1075 // vector<TLine*> vLineZX;
1076 // vector<TLine*> vLineZY;
1077 // Double_t z0 = tr->GetStartZ();
1078 // Double_t x0 = tr->GetStartX();
1079 // Double_t y0 = tr->GetStartY();
1080 // Double_t x1, y1, z1;
1081 // set<Int_t> vSt;
1082 // for (Int_t i = 0; i < fSilPointsArray->GetEntriesFast(); ++i) {
1083 // FairMCPoint* pnt = (FairMCPoint*)fSilPointsArray->At(i);
1084 // if (pnt->GetTrackID() == iTr) {
1085 // nPoints++;
1086 // z1 = pnt->GetZ();
1087 // x1 = pnt->GetX();
1088 // y1 = pnt->GetY();
1089 // vSt.insert(fSilDetector->GetPointStationOwnership(z1));
1090 // vLineZX.push_back(new TLine(z0, x0, z1, x1));
1091 // vLineZY.push_back(new TLine(z0, y0, z1, y1));
1092 // z0 = z1;
1093 // x0 = x1;
1094 // y0 = y1;
1095 // }
1096 // }
1097
1098 // for (Int_t i = 0; i < fGemPointsArray->GetEntriesFast(); ++i) {
1099 // FairMCPoint* pnt = (FairMCPoint*)fGemPointsArray->At(i);
1100 // if (pnt->GetTrackID() == iTr) {
1101 // nPoints++;
1102 // z1 = pnt->GetZ();
1103 // x1 = pnt->GetX();
1104 // y1 = pnt->GetY();
1105 // vSt.insert(fGemDetector->GetPointStationOwnership(z1) + fNSiliconStations);
1106 // vLineZX.push_back(new TLine(z0, x0, z1, x1));
1107 // vLineZY.push_back(new TLine(z0, y0, z1, y1));
1108 // z0 = z1;
1109 // x0 = x1;
1110 // y0 = y1;
1111 // }
1112 // }
1113
1114 // if (nPoints > 3 && vSt.size() > 3) {
1115 // nMCtracks++;
1116 // c->cd(1);
1117 // for (TLine* it : vLineZX) {
1118 // it->SetLineColor(kBlue);
1119 // it->SetLineStyle(9);
1120 // it->Draw("same");
1121 // }
1122 // c->cd(2);
1123 // for (TLine* it : vLineZY) {
1124 // it->SetLineColor(kBlue);
1125 // it->SetLineStyle(9);
1126 // it->Draw("same");
1127 // }
1128 // }
1129 // }
1130
1131 // printf("NMCTracks = %d GlobTracks = %d\n", nMCtracks, fGlobTracksArray->GetEntriesFast());
1132
1133 // arc->Draw("same");
1134 c->SaveAs("hits.png");
1135 getchar();
1136 delete h_HitsZX;
1137 delete h_HitsZY;
1138 delete c;
1139
1140 return kBMNSUCCESS;
1141}
LineFit
TVector3 LineFit(BmnTrack *track, const TClonesArray *arr, TString type)
Definition BmnMath.cxx:113
CircleBy3Hit
TVector3 CircleBy3Hit(BmnTrack *track, const TClonesArray *arr)
Definition BmnMath.cxx:470
Mu
Double_t Mu(vector< Double_t > qp, vector< Double_t > w)
Definition BmnMath.cxx:906
CircleFit
TVector3 CircleFit(BmnTrack *track, const TClonesArray *arr, Double_t &chi2)
Definition BmnMath.cxx:312
CalcTx
Double_t CalcTx(const BmnHit *h0, const BmnHit *h1, const BmnHit *h2)
Definition BmnMath.cxx:422
dist
vector< Double_t > dist(vector< Double_t > qp, Double_t mu)
Definition BmnMath.cxx:869
IsParCorrect
Bool_t IsParCorrect(const FairTrackParam *par, const Bool_t isField)
Definition BmnMath.cxx:59
W
vector< Double_t > W(vector< Double_t > dist, Double_t sig)
Definition BmnMath.cxx:878
workTime
static Double_t workTime
Definition SrcInnerTrackingRun7.cxx:26
i
int i
Definition P4_F32vec4.h:22
Bool_t
bool Bool_t
Int_t
int Int_t
Short_t
short Short_t
kFALSE
constexpr Bool_t kFALSE
Double_t
double Double_t
kTRUE
constexpr Bool_t kTRUE
Option_t
const char Option_t
kRed
kRed
Fatal
void Fatal(const char *location, const char *msgfmt,...)
gSystem
R__EXTERN TSystem * gSystem
kGEM
@ kGEM
Definition BmnDetectorList.h:15
BmnStatus
BmnStatus
Define status code.
Definition BmnEnums.h:93
kBMNERROR
@ kBMNERROR
Error.
Definition BmnEnums.h:95
kBMNSUCCESS
@ kBMNSUCCESS
Success.
Definition BmnEnums.h:94
BmnFieldMap::IsFieldOff
Bool_t IsFieldOff() const
Whether magnetic field is off.
Definition BmnFieldMap.h:93
BmnGemStripHit::GetStripTotalSignalInUpperLayer
Double_t GetStripTotalSignalInUpperLayer()
Definition BmnGemStripHit.cxx:30
BmnGemStripHit::GetStripTotalSignalInLowerLayer
Double_t GetStripTotalSignalInLowerLayer()
Strip layer characteristics of the hit (total signal, cluster size, strip position as the cluster mea...
Definition BmnGemStripHit.cxx:24
BmnGlobalTrack
BmnGlobalTrack.h.
Definition BmnGlobalTrack.h:14
BmnGlobalTrack::SetGemTrackIndex
void SetGemTrackIndex(Int_t iGem)
Modifiers.
Definition BmnGlobalTrack.h:83
BmnGlobalTrack::GetGemTrackIndex
Int_t GetGemTrackIndex() const
Accessors.
Definition BmnGlobalTrack.h:34
BmnHit::GetIndex
Int_t GetIndex() const
Definition BmnHit.h:32
BmnHit::GetDetId
DetectorId GetDetId() const
Definition BmnHit.h:34
BmnHit::SetUsing
void SetUsing(Bool_t use)
Definition BmnHit.h:47
BmnHit::SetDetId
void SetDetId(DetectorId det)
Definition BmnHit.h:51
BmnHit::SetIndex
void SetIndex(Int_t id)
Definition BmnHit.h:49
BmnHit::IsUsed
Bool_t IsUsed() const
Definition BmnHit.h:30
BmnHit::SetStation
void SetStation(Short_t st)
Definition BmnHit.h:52
BmnHit::GetStation
Short_t GetStation() const
Definition BmnHit.h:35
BmnTrack::GetNDF
Int_t GetNDF() const
Definition BmnTrack.h:57
BmnTrack::SetChi2
void SetChi2(Double_t chi2)
Definition BmnTrack.h:76
BmnTrack::SetParamLast
void SetParamLast(FairTrackParam &par)
Definition BmnTrack.h:72
BmnTrack::SetFlag
void SetFlag(Int_t flag)
Definition BmnTrack.h:74
BmnTrack::GetParamFirst
FairTrackParam * GetParamFirst()
Definition BmnTrack.h:63
BmnTrack::GetChi2
Float_t GetChi2() const
Definition BmnTrack.h:55
BmnTrack::GetNHits
Int_t GetNHits() const
Accessors.
Definition BmnTrack.h:49
BmnTrack::GetHitIndex
Int_t GetHitIndex(Int_t iHit) const
Definition BmnTrack.h:51
BmnTrack::AddHit
void AddHit(Int_t hitIndex, FairHit *Hit)
Associate a FairHit to the track.
Definition BmnTrack.cxx:28
BmnTrack::SetNHits
void SetNHits(Int_t n)
Definition BmnTrack.h:80
BmnTrack::SetParamFirst
void SetParamFirst(FairTrackParam &par)
Modifiers.
Definition BmnTrack.h:70
BmnTrack::GetParamLast
FairTrackParam * GetParamLast()
Definition BmnTrack.h:65
BmnTrack::SetLength
void SetLength(Double_t length)
Definition BmnTrack.h:84
BmnTrack::GetFlag
Int_t GetFlag() const
Definition BmnTrack.h:53
BmnTrack::SortHits
void SortHits()
Public method SortHits Sorts the hits in downstream direction and writes the hit indices into the mem...
Definition BmnTrack.cxx:45
BmnTrack::SetNDF
void SetNDF(Int_t ndf)
Definition BmnTrack.h:78
SrcInnerTrackingRun7::MatchHit
BmnStatus MatchHit(BmnTrack *cand, vector< BmnHit * > sortedHits, Bool_t downstream)
Definition SrcInnerTrackingRun7.cxx:462
SrcInnerTrackingRun7::fGlobTracksArray
TClonesArray * fGlobTracksArray
Definition SrcInnerTrackingRun7.h:96
SrcInnerTrackingRun7::SortTracks
BmnStatus SortTracks(vector< BmnTrack > &inTracks, vector< BmnTrack > &sortedTracks)
Definition SrcInnerTrackingRun7.cxx:501
SrcInnerTrackingRun7::fGemDetector
BmnGemStripStationSet * fGemDetector
Definition SrcInnerTrackingRun7.h:91
SrcInnerTrackingRun7::Init
virtual InitStatus Init()
Definition SrcInnerTrackingRun7.cxx:64
SrcInnerTrackingRun7::CalcCovMatrix
BmnStatus CalcCovMatrix(BmnTrack *tr)
Definition SrcInnerTrackingRun7.cxx:672
SrcInnerTrackingRun7::TrackUpdateByLine
BmnStatus TrackUpdateByLine(vector< BmnTrack > &cands)
Definition SrcInnerTrackingRun7.cxx:552
SrcInnerTrackingRun7::CalcMeanSigma
TVector2 CalcMeanSigma(vector< Double_t >)
Definition SrcInnerTrackingRun7.cxx:813
SrcInnerTrackingRun7::CalculateTrackParams
BmnStatus CalculateTrackParams(BmnTrack *tr)
Definition SrcInnerTrackingRun7.cxx:772
SrcInnerTrackingRun7::FindCandidateByThreeStations
BmnStatus FindCandidateByThreeStations(Short_t st0, Short_t st1, Short_t st2, vector< BmnTrack > &candidates, vector< BmnHit * > *sortedHits)
Definition SrcInnerTrackingRun7.cxx:420
SrcInnerTrackingRun7::SortHits
BmnStatus SortHits(vector< BmnHit * > *sortedHits)
Definition SrcInnerTrackingRun7.cxx:532
SrcInnerTrackingRun7::TrackSelection
BmnStatus TrackSelection(vector< BmnTrack > &cands)
Definition SrcInnerTrackingRun7.cxx:601
SrcInnerTrackingRun7::TrackUpdateByKalman
BmnStatus TrackUpdateByKalman(vector< BmnTrack > &cands)
Definition SrcInnerTrackingRun7.cxx:572
SrcInnerTrackingRun7::FindCandidateByFourStations
BmnStatus FindCandidateByFourStations(Short_t st0, Short_t st1, Short_t st2, Short_t st3, vector< BmnTrack > &candidates, vector< BmnHit * > *sortedHits)
Definition SrcInnerTrackingRun7.cxx:372
SrcInnerTrackingRun7::FindCandidateByFiveStations
BmnStatus FindCandidateByFiveStations(Short_t st0, Short_t st1, Short_t st2, Short_t st3, Short_t st4, vector< BmnTrack > &candidates, vector< BmnHit * > *sortedHits)
Definition SrcInnerTrackingRun7.cxx:333
SrcInnerTrackingRun7::SetHitsUsing
void SetHitsUsing(BmnTrack *tr, Bool_t use)
Definition SrcInnerTrackingRun7.cxx:644
SrcInnerTrackingRun7::CalculateLength
Double_t CalculateLength(BmnTrack *tr)
Definition SrcInnerTrackingRun7.cxx:943
SrcInnerTrackingRun7::Exec
virtual void Exec(Option_t *opt)
Definition SrcInnerTrackingRun7.cxx:121
SrcInnerTrackingRun7::CheckSharedHits
BmnStatus CheckSharedHits(vector< BmnTrack > &sortedTracks)
Definition SrcInnerTrackingRun7.cxx:969
SrcInnerTrackingRun7::CalcQp
Double_t CalcQp(BmnTrack *track, Bool_t doSimple)
Definition SrcInnerTrackingRun7.cxx:828
TAttLine::SetLineWidth
virtual void SetLineWidth(Width_t lwidth)
TAttLine::SetLineColor
virtual void SetLineColor(Color_t lcolor)
TAttMarker::SetMarkerColor
virtual void SetMarkerColor(Color_t mcolor=1)
TAttMarker::SetMarkerStyle
virtual void SetMarkerStyle(Style_t mstyle=1)
TAttMarker::SetMarkerSize
virtual void SetMarkerSize(Size_t msize=1)
TCanvas::cd
TVirtualPad * cd(Int_t subpadnumber=0) override
TH1::Draw
void Draw(Option_t *option="") override
TH2::Fill
virtual Int_t Fill(const char *namex, const char *namey, Double_t w)
TObject::Draw
virtual void Draw(Option_t *option="")
TPad::Divide
void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0) override
TPad::SaveAs
void SaveAs(const char *filename="", Option_t *option="") const override
TVector2::Y
Double_t Y() const
TVector2::X
Double_t X() const
TVector3::Z
Double_t Z() const
TVector3::Y
Double_t Y() const
TVector3::X
Double_t X() const
n
*brief user defined string literal for JSON values *sa std::size_t n
Definition json.hpp:20362
TMath::Sigma
constexpr Double_t Sigma()
TMath::Sq
Double_t Sq(Double_t x)
TMath::Sqrt
Double_t Sqrt(Double_t x)
TMath::R
constexpr Double_t R()
TMath::Abs
Double_t Abs(Double_t d)
std
Definition json.hpp:4446