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