BmnRoot 26.10.0-18-g4e3414f9
Software framework of the BM@N experiment
Loading...
Searching...
No Matches
BmnInnerTrackingRun7.cxx
Go to the documentation of this file.
1#include "BmnInnerTrackingRun7.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 "BmnSiliconHit.h"
10#include "FairRunAna.h"
11#include "FairTrackParam.h"
12#include "TStyle.h"
13
14#include <TCanvas.h>
15#include <TH2F.h>
16#include <TLine.h>
17#include <TMath.h>
18#include <map>
19#include <vector>
20#if defined(_OPENMP)
21#include "omp.h"
22
23#include <cstdlib>
24#endif
25
26//-----------------------------------------
27static Double_t workTime = 0.0;
28static Double_t hitMatchTime = 0.0;
29static Double_t hitMatchTimeLoop = 0.0;
30static Double_t construct_4of4_time = 0.0;
31static Double_t construct_3of4_time = 0.0;
32static Double_t construct_2of4_time = 0.0;
33static Double_t CalculateTrackParTime = 0.0;
34static Double_t SortHitsTime = 0.0;
35static Double_t trackUpdateTime = 0.0;
36static Double_t sortTrackTime = 0.0;
37static Double_t trackSelectionTime = 0.0;
38//-----------------------------------------
39
40using namespace std;
41using namespace TMath;
42
43BmnInnerTrackingRun7::BmnInnerTrackingRun7(Int_t run, Bool_t target)
44 : fSteerFile("")
45 , fSteering(nullptr)
46{
47 fEventNo = 0;
48 fIsTarget = target;
49 fNSiliconStations = 0;
50 fGemHitsArray = nullptr;
51 fSilHitsArray = nullptr;
52 fHitsArray = nullptr;
53 fRoughVertex = TVector3(0.5, -4.6, -2.3);
54 fGlobTracksArray = nullptr;
55 fGemTracksArray = nullptr;
56 fSilTracksArray = nullptr;
57 fField = nullptr;
58 fGemHitsBranchName = "BmnGemStripHit";
59 fSilHitsBranchName = "BmnSiliconHit";
60 fGlobTracksBranchName = "BmnGlobalTrack";
61 fGemTracksBranchName = "BmnGemTrack";
62 fSilTracksBranchName = "BmnSiliconTrack";
63 fGemDetector = nullptr;
64 fSilDetector = nullptr;
65 fHitXCutMin = nullptr;
66 fHitXCutMax = nullptr;
67 fHitYCutMin = nullptr;
68 fHitYCutMax = nullptr;
69 fNHitsCut = 0;
70 fDoHitAsymFiltration = kTRUE;
71}
72
73BmnInnerTrackingRun7::~BmnInnerTrackingRun7()
74{
75 delete fKalman;
76 delete fGemDetector;
77 delete fSilDetector;
78 delete[] fHitXCutMin;
79 delete[] fHitXCutMax;
80 delete[] fHitYCutMin;
81 delete[] fHitYCutMax;
82}
83
84InitStatus BmnInnerTrackingRun7::Init()
85{
86 if (fVerbose > 1)
87 cout << "======================== GEM tracking init started ====================" << endl;
88
89 fKalman = new BmnKalmanFilter();
90
91 // Get ROOT Manager
92 FairRootManager* ioman = FairRootManager::Instance();
93 if (!ioman)
94 Fatal("Init", "FairRootManager is not instantiated");
95
96 // MC tracks for drawing during debugging
97 fMCTracksArray = (TClonesArray*)ioman->GetObject("MCTrack"); // in
98
99 // GEM
100 fGemPointsArray = (TClonesArray*)ioman->GetObject("StsPoint"); // in
101 fGemHitsArray = (TClonesArray*)ioman->GetObject(fGemHitsBranchName); // in
102 fGemTracksArray = new TClonesArray("BmnTrack", 100); // out
103 ioman->Register(fGemTracksBranchName, "GEM", fGemTracksArray, kTRUE);
104
105 // SILICON
106 fSilPointsArray = (TClonesArray*)ioman->GetObject("SiliconPoint"); // in
107 fSilHitsArray = (TClonesArray*)ioman->GetObject(fSilHitsBranchName); // in
108 fSilTracksArray = new TClonesArray("BmnTrack", 100); // out
109 ioman->Register(fSilTracksBranchName, "SILICON", fSilTracksArray, kTRUE);
110
111 if (!fGemHitsArray) {
112 cout << "BmnInnerTrackingRun7::Init(): branch " << fGemHitsBranchName << " not found! Task will be deactivated"
113 << endl;
114 SetActive(kFALSE);
115 return InitStatus::kERROR;
116 }
117
118 fGlobTracksArray = new TClonesArray("BmnGlobalTrack", 100); // out
119 ioman->Register(fGlobTracksBranchName, "GLOBAL", fGlobTracksArray, kTRUE);
120 fHitsArray = new TClonesArray("BmnHit", 100); // out
121 ioman->Register("BmnInnerHits", "HITS", fHitsArray, kTRUE);
122
123 fField = FairRunAna::Instance()->GetField();
124 if (!fField)
125 Fatal("Init", "No Magnetic Field found");
126
127 TString gPathConfig = gSystem->Getenv("VMCWORKDIR");
128
129 TString gPathGemConfig = gPathConfig + "/parameters/gem/XMLConfigs/";
130 TString confGem = "GemRunSpring2018.xml";
131 fGemDetector = new BmnGemStripStationSet(gPathGemConfig + confGem);
132
133 TString gPathSiConfig = gPathConfig + "/parameters/silicon/XMLConfigs/";
134 TString confSi = "SiliconRunSpring2018.xml";
135 fSilDetector = new BmnSiliconStationSet(gPathSiConfig + confSi);
136
137 Int_t nGemStations = fGemDetector->GetNStations();
138 Int_t nSilStations = fSilDetector->GetNStations();
139
140 fNSiliconStations = nSilStations;
141
142 fNStations = nGemStations + nSilStations;
143
144 fHitXCutMin = new Double_t[fNStations];
145 fHitXCutMax = new Double_t[fNStations];
146 fHitYCutMin = new Double_t[fNStations];
147 fHitYCutMax = new Double_t[fNStations];
148
149 BmnFieldMap* field = (BmnFieldMap*)FairRunAna::Instance()->GetField();
150 Bool_t isField = !(field->IsFieldOff());
151 fSteering = new BmnSteering(isField ? "BMN_run7_withField.dat" : "BMN_run7_noField.dat");
152
153 if (fVerbose > 1)
154 cout << "======================== GEM tracking init finished ===================" << endl;
155
156 return kSUCCESS;
157}
158
159void BmnInnerTrackingRun7::Exec(Option_t* opt)
160{
161 if (fVerbose > 1)
162 cout << "\n======================== GEM tracking exec started ====================" << endl;
163 if (fVerbose > 1)
164 cout << "\n Event number: " << fEventNo << endl;
165
166 if (!IsActive())
167 return;
168
169 clock_t tStart = clock();
170
171 fSilTracksArray->Delete();
172 fGemTracksArray->Delete();
173 fGlobTracksArray->Delete();
174 fHitsArray->Delete();
175
176 fEventNo++;
177
178 Int_t nHitsCut = fSteering->GetNHitsCutTotal();
179
180 Int_t nSilStations = fSilDetector->GetNStations();
181 // Int_t nGemStations = fGemDetector->GetNStations();
182
183 Double_t a[2][9] = {{1.29, 1.50, 1.33, 1.14, 1.17, 1.14, 1.00, 1.14, 1.14},
184 {1.33, 1.17, 1.00, 1.14, 1.33, 1.33, 1.14, 1.33, 1.33}};
185 Double_t b[2][9] = {{-645, -1500, -1330, -570, -585, -570, -500, -570, -570},
186 {500, 1000, 1000, 500, 500, 500, 500, 500, 500}};
187
188 for (Int_t iHit = 0; iHit < fSilHitsArray->GetEntriesFast(); ++iHit) {
189 BmnSiliconHit* hit = (BmnSiliconHit*)fSilHitsArray->At(iHit);
190 Double_t Sl = hit->GetStripTotalSignalInLowerLayer();
191 Double_t Su = hit->GetStripTotalSignalInUpperLayer();
192 Int_t iSt = hit->GetStation();
193 if (fDoHitAsymFiltration) {
194 if (Sl < a[0][iSt] * Su + b[0][iSt] || Sl > a[1][iSt] * Su + b[1][iSt]) {
195 hit->SetType(0);
196 continue;
197 } else
198 hit->SetType(1);
199 }
200 BmnHit innerHit = *hit;
201 innerHit.SetIndex(iHit);
202 innerHit.SetDetId(kSILICON);
203 innerHit.SetDxyz(0.5, 0.5, 0.5);
204 // innerHit.SetDxyz(hit->GetDx(), hit->GetDy(), hit->GetDz());
205 new ((*fHitsArray)[fHitsArray->GetEntriesFast()]) BmnHit(innerHit);
206 }
207 for (Int_t iHit = 0; iHit < fGemHitsArray->GetEntriesFast(); ++iHit) {
208 BmnGemStripHit* hit = (BmnGemStripHit*)fGemHitsArray->At(iHit);
209 Double_t Sl = hit->GetStripTotalSignalInLowerLayer();
210 Double_t Su = hit->GetStripTotalSignalInUpperLayer();
211 Int_t iSt = hit->GetStation() + nSilStations;
212 if (fDoHitAsymFiltration) {
213 if (Sl < a[0][iSt] * Su + b[0][iSt] || Sl > a[1][iSt] * Su + b[1][iSt]) {
214 hit->SetType(0);
215 continue;
216 } else
217 hit->SetType(1);
218 }
219 BmnHit innerHit = *hit;
220 innerHit.SetStation(iSt); // shift for correct station numbering
221 innerHit.SetIndex(iHit);
222 innerHit.SetDetId(kGEM);
223 innerHit.SetDxyz(0.5, 0.5, 0.5);
224 // innerHit.SetDxyz(hit->GetDx(), hit->GetDy(), hit->GetDz());
225 new ((*fHitsArray)[fHitsArray->GetEntriesFast()]) BmnHit(innerHit);
226 }
227
228 if (fHitsArray->GetEntriesFast() > nHitsCut || fHitsArray->GetEntriesFast() == 0)
229 return;
230
231 for (Int_t iStat = 0; iStat < fNStations; iStat++) {
232 fHitXCutMin[iStat] = fSteering->GetHitXCutMin()[iStat];
233 fHitXCutMax[iStat] = fSteering->GetHitXCutMax()[iStat];
234 fHitYCutMin[iStat] = fSteering->GetHitYCutMin()[iStat];
235 fHitYCutMax[iStat] = fSteering->GetHitYCutMax()[iStat];
236 }
237
238 fNHitsCut = fSteering->GetNHitsCut();
239 fChiSquareCut = fSteering->GetChiSquareCut();
240
241 fDistCut = 0.45;
242 fNHitsCut = 4;
243 FindTracks_4of4_OnLastGEMStations();
244 // cout << "FindTracks_4of4_OnLastStations: " << fGlobTracksArray->GetEntriesFast() << " tracks" << endl;
245 FindTracks_3of4_OnLastGEMStations();
246 // cout << "FindTracks_3of4_OnLastStations: " << fGlobTracksArray->GetEntriesFast() << " tracks" << endl;
247 fNHitsCut = 5;
248 FindTracks_2of2_OnFirstGEMStationsDownstream();
249 // cout << "FindTracks_2of2_OnFirstGEMStationsDownstream: " << fGlobTracksArray->GetEntriesFast() << " tracks" <<
250 // endl;
251 FindTracks_2of2_OnFirstGEMStationsUpstream();
252 // cout << "FindTracks_2of2_OnFirstGEMStationsUpstream: " << fGlobTracksArray->GetEntriesFast() << " tracks" <<
253 // endl;
254
255 // DrawHits();
256
257 clock_t tFinish = clock();
258 if (fVerbose > 0)
259 cout << "BmnInnerTrackingRun7: " << fGlobTracksArray->GetEntriesFast() << " tracks" << endl;
260
261 workTime += ((Double_t)(tFinish - tStart)) / CLOCKS_PER_SEC;
262
263 if (fVerbose > 1)
264 cout << "\n======================== GEM tracking exec finished ===================" << endl;
265}
266
267BmnStatus BmnInnerTrackingRun7::FindTracks_4of4_OnLastGEMStations()
268{
269 // Fit of dX vs X for different stations {p0, p1, sigma} (from qgsm simulation)
270 const Int_t nxRanges = 8;
271 const Int_t nyRanges = 5;
272 Double_t xRangeMin[nxRanges] = {-80, -60, -40, -20, 0, 20, 40, 60};
273 Double_t xRangeMax[nxRanges] = {-60, -40, -20, 0, 20, 40, 60, 80};
274 map<Short_t, vector<Double_t>> xRMS = {//<station, vector of RMS >
275 {6, {7.0, 5.3, 4.4, 3.5, 3.4, 3.9, 5.1, 5.3}},
276 {7, {6.8, 6.3, 4.8, 3.7, 3.7, 4.5, 5.3, 6.0}},
277 {8, {6.0, 5.9, 5.2, 4.3, 3.4, 4.2, 4.8, 5.1}}};
278 Double_t yRangeMin[nyRanges] = {-10, 0, 10, 20, 30};
279 Double_t yRangeMax[nyRanges] = {0, 10, 20, 30, 40};
280 map<Short_t, vector<Double_t>> yRMS = {//<station, vector of RMS >
281 {6, {0.8, 0.8, 0.9, 1.0, 0.9}},
282 {7, {0.9, 0.7, 1.0, 1.0, 0.9}},
283 {8, {0.5, 0.6, 0.8, 0.7, 0.7}}};
284
285 Double_t Y6[2] = {0.87, 0.16};
286 Double_t Y7[2] = {0.84, 0.15};
287 Double_t Y8[2] = {0.61, 0.13};
288 Double_t X6[2] = {1.04, 0.19};
289 Double_t X7[2] = {1.42, 0.18};
290 Double_t X8[2] = {1.41, 0.14};
291
292 vector<BmnHit*> sortedHits[fNStations];
293 SortHits(sortedHits);
294
295 vector<BmnTrack> candidates;
296 vector<BmnTrack> sortedCandidates;
297
298 clock_t t0 = 0; // clock();
299 // Int_t nThreads = 1;
300 // #if defined(_OPENMP)
301 // //char * number_of_threads = getenv("OMP_NUM_THREADS");
302 // //nThreads = atoi(number_of_threads);
303 // nThreads = 4;
304 // omp_set_num_threads(nThreads);
305 // #endif
306 // #pragma omp parallel for
307 for (BmnHit* hit8 : sortedHits[8]) {
308 for (BmnHit* hit7 : sortedHits[7]) {
309 Double_t x8 = hit8->GetX();
310 Double_t y8 = hit8->GetY();
311 Double_t dxEx = x8 - hit7->GetX();
312 Double_t dxTh = x8 * X8[1] + X8[0];
313 Double_t dyEx = y8 - hit7->GetY();
314 Double_t dyTh = y8 * Y8[1] + Y8[0];
315
316 Int_t iX = 0;
317 for (iX = 0; iX < nxRanges; ++iX)
318 if (x8 > xRangeMin[iX] && x8 < xRangeMax[iX])
319 break;
320 Int_t iY = 0;
321 for (iY = 0; iY < nyRanges; ++iY)
322 if (y8 > yRangeMin[iY] && y8 < yRangeMax[iY])
323 break;
324 if ((iX >= nxRanges) || (iY >= nyRanges)) {
325 if (fVerbose > 0)
326 cout << "WARNING: x8 or y8 is out of range, the hit is skipped (x = " << x8 << ", y = " << y8 << ")"
327 << endl;
328 continue;
329 }
330 Double_t rmsX = 2.0 * xRMS[8][iX];
331 Double_t rmsY = 3.0 * yRMS[8][iY]; // 2.0 * yRMS[8][iY];
332
333 if (Abs(dxEx - dxTh) > rmsX || Abs(dyEx - dyTh) > rmsY)
334 continue;
335 for (BmnHit* hit6 : sortedHits[6]) {
336 Double_t x7 = hit7->GetX();
337 Double_t y7 = hit7->GetY();
338 dxEx = x7 - hit6->GetX();
339 dxTh = x7 * X7[1] + X7[0];
340 dyEx = y7 - hit6->GetY();
341 dyTh = y7 * Y7[1] + Y7[0];
342
343 iX = 0;
344 for (iX = 0; iX < nxRanges; ++iX)
345 if (x7 > xRangeMin[iX] && x7 < xRangeMax[iX])
346 break;
347 iY = 0;
348 for (iY = 0; iY < nyRanges; ++iY)
349 if (y7 > yRangeMin[iY] && y7 < yRangeMax[iY])
350 break;
351 if ((iX >= nxRanges) || (iY >= nyRanges)) {
352 if (fVerbose > 0)
353 cout << "WARNING: x7 or y7 is out of range, the hit is skipped (x = " << x7 << ", y = " << y7
354 << ")" << endl;
355 continue;
356 }
357 rmsX = 2.0 * xRMS[7][iX];
358 rmsY = 3.0 * yRMS[7][iY]; // 2.0 * yRMS[7][iY];
359
360 if (Abs(dxEx - dxTh) > rmsX || Abs(dyEx - dyTh) > rmsY)
361 continue;
362 for (BmnHit* hit5 : sortedHits[5]) {
363 Double_t x6 = hit6->GetX();
364 Double_t y6 = hit6->GetY();
365 dxEx = x6 - hit5->GetX();
366 dxTh = x6 * X6[1] + X6[0];
367 dyEx = y6 - hit5->GetY();
368 dyTh = y6 * Y6[1] + Y6[0];
369
370 iX = 0;
371 for (iX = 0; iX < nxRanges; ++iX)
372 if (x6 > xRangeMin[iX] && x6 < xRangeMax[iX])
373 break;
374 iY = 0;
375 for (iY = 0; iY < nyRanges; ++iY)
376 if (y6 > yRangeMin[iY] && y6 < yRangeMax[iY])
377 break;
378 if ((iX >= nxRanges) || (iY >= nyRanges)) {
379 if (fVerbose > 0)
380 cout << "WARNING: x6 or y6 is out of range, the hit is skipped (x = " << x6
381 << ", y = " << y6 << ")" << endl;
382 continue;
383 }
384 rmsX = 2.0 * xRMS[6][iX];
385 rmsY = 3.0 * yRMS[6][iY]; // 2.0 * yRMS[6][iY];
386
387 if (Abs(dxEx - dxTh) > rmsX || Abs(dyEx - dyTh) > rmsY)
388 continue;
389 BmnTrack cand;
390 cand.AddHit(hit8->GetRefIndex(), hit8);
391 cand.AddHit(hit7->GetRefIndex(), hit7);
392 cand.AddHit(hit6->GetRefIndex(), hit6);
393 cand.AddHit(hit5->GetRefIndex(), hit5);
394 cand.SortHits();
395 candidates.push_back(cand);
396 }
397 }
398 }
399 }
400 clock_t t1 = 0; // clock();
401 construct_4of4_time += ((Double_t)(t1 - t0)) / CLOCKS_PER_SEC;
402
403 // cout << "FindTracks_4of4_OnLastGEMStations::candidates = " << candidates.size() << endl;
404
405 for (BmnTrack& cand : candidates) {
406 // for (Int_t iCand = 0; iCand < candidates.size(); ++iCand) {
407 // BmnTrack cand = candidates[iCand];
408 if (CalculateTrackParams(&cand) == kBMNERROR)
409 continue;
410 // First two GEM
411 MatchHit(&cand, sortedHits[4], kFALSE);
412 MatchHit(&cand, sortedHits[3], kFALSE);
413 // Silicon
414 MatchHit(&cand, sortedHits[2], kFALSE);
415 MatchHit(&cand, sortedHits[1], kFALSE);
416 MatchHit(&cand, sortedHits[0], kFALSE);
417
418 // NDF = (N counts in ZX plane + N counts in ZY plane) - 2 parameters of Line in ZY plane - 3 parameters of
419 // Circle in ZX plane Check it!!!
420 cand.SetNDF(cand.GetNHits() * 2 - 5);
421 }
422
423 // for (BmnTrack& cand : candidates) {
424 // cout << cand.GetNDF() << endl;
425 // cand.Print();
426 // }
427
428 BmnFieldMap* field = (BmnFieldMap*)FairRunAna::Instance()->GetField();
429 Bool_t isField = !(field->IsFieldOff());
430 if (isField)
431 TrackUpdateByKalman(candidates);
432 else
433 TrackUpdateByLine(candidates);
434 SortTracks(candidates, sortedCandidates);
435 TrackSelection(sortedCandidates);
436
437 return kBMNSUCCESS;
438}
439
440BmnStatus BmnInnerTrackingRun7::FindCandidateByThreeStations(Short_t st0,
441 Short_t st1,
442 Short_t st2,
443 vector<BmnTrack>& candidates,
444 vector<BmnHit*>* sortedHits)
445{
446 clock_t t0 = 0; // clock();
447 // Fit of dX vs X (dY vs Y) for different stations {p0, p1, sigma} (from qgsm simulation)
448 const Int_t nxRanges = 8;
449 const Int_t nyRanges = 5;
450 Double_t xRangeMin[nxRanges] = {-80, -60, -40, -20, 0, 20, 40, 60};
451 Double_t xRangeMax[nxRanges] = {-60, -40, -20, 0, 20, 40, 60, 80};
452
453 map<Short_t, vector<Double_t>> xRMS = {//<station, vector of RMS >
454 {6, {7.0, 5.3, 4.4, 3.5, 3.4, 3.9, 5.1, 5.3}},
455 {7, {6.8, 6.3, 4.8, 3.7, 3.7, 4.5, 5.3, 6.0}},
456 {8, {6.0, 5.9, 5.2, 4.3, 3.4, 4.2, 4.8, 5.1}}};
457 Double_t yRangeMin[nyRanges] = {-10, 0, 10, 20, 30};
458 Double_t yRangeMax[nyRanges] = {0, 10, 20, 30, 40};
459 map<Short_t, vector<Double_t>> yRMS = {//<station, vector of RMS >
460 {6, {0.8, 0.8, 0.9, 1.0, 0.9}},
461 {7, {0.9, 0.7, 1.0, 1.0, 0.9}},
462 {8, {0.5, 0.6, 0.8, 0.7, 0.7}}};
463 map<Short_t, vector<Double_t>> xMap = {//<station, <p0, p1> >
464 {6, {1.04, 0.19}},
465 {7, {1.42, 0.18}},
466 {8, {1.41, 0.14}}};
467 map<Short_t, vector<Double_t>> yMap = {//<station, <p0, p1> >
468 {6, {0.87, 0.16}},
469 {7, {0.84, 0.15}},
470 {8, {0.61, 0.13}}};
471
472 for (BmnHit* hit2 : sortedHits[st2]) {
473 for (BmnHit* hit1 : sortedHits[st1]) {
474 Double_t x2 = hit2->GetX();
475 Double_t y2 = hit2->GetY();
476 Double_t x1 = hit1->GetX();
477 Double_t y1 = hit1->GetY();
478 Double_t dxEx = Abs(x2 - x1);
479 Double_t dxTh = Abs(x2 * xMap[st2][1] + xMap[st2][0]);
480 Double_t dyEx = Abs(y2 - y1);
481 Double_t dyTh = Abs(y2 * yMap[st2][1] + yMap[st2][0]);
482
483 Int_t iX = 0;
484 for (iX = 0; iX < nxRanges; ++iX)
485 if (x2 > xRangeMin[iX] && x2 < xRangeMax[iX])
486 break;
487 Int_t iY = 0;
488 for (iY = 0; iY < nyRanges; ++iY)
489 if (y2 > yRangeMin[iY] && y2 < yRangeMax[iY])
490 break;
491 if ((iX >= nxRanges) || (iY >= nyRanges)) {
492 if (fVerbose > 0)
493 cout << "WARNING: x2 or y2 is out of range, the hit is skipped (x = " << x2 << ", y = " << y2 << ")"
494 << endl;
495 continue;
496 }
497 Double_t rmsX = 2.0 * xRMS[st2][iX];
498 Double_t rmsY = 3.0 * yRMS[st2][iY]; // 2.0 * yRMS[st2][iY];
499
500 if (Abs(dxEx - dxTh) > rmsX || Abs(dyEx - dyTh) > rmsY)
501 continue;
502 for (BmnHit* hit0 : sortedHits[st0]) {
503 dxEx = Abs(x1 - hit0->GetX());
504 dxTh = Abs(x1 * xMap[st1][1] + xMap[st1][0]);
505 dyEx = Abs(y1 - hit0->GetY());
506 dyTh = Abs(y1 * yMap[st1][1] + yMap[st1][0]);
507
508 iX = 0;
509 for (iX = 0; iX < nxRanges; ++iX)
510 if (x1 > xRangeMin[iX] && x1 < xRangeMax[iX])
511 break;
512 iY = 0;
513 for (iY = 0; iY < nyRanges; ++iY)
514 if (y1 > yRangeMin[iY] && y1 < yRangeMax[iY])
515 break;
516 if ((iX >= nxRanges) || (iY >= nyRanges)) {
517 if (fVerbose > 0)
518 cout << "WARNING: x1 or y1 is out of range, the hit is skipped (x = " << x1 << ", y = " << y1
519 << ")" << endl;
520 continue;
521 }
522 rmsX = 2.0 * xRMS[st1][iX];
523 rmsY = 3.0 * yRMS[st1][iY]; // 2.0 * yRMS[st1][iY];
524
525 if (Abs(dxEx - dxTh) > rmsX || Abs(dyEx - dyTh) > rmsY)
526 continue;
527 BmnTrack cand;
528 cand.AddHit(hit2->GetRefIndex(), hit2);
529 cand.AddHit(hit1->GetRefIndex(), hit1);
530 cand.AddHit(hit0->GetRefIndex(), hit0);
531 cand.SortHits();
532 CalculateTrackParams(&cand);
533 candidates.push_back(cand);
534 }
535 }
536 }
537 clock_t t1 = 0; // clock();
538 construct_3of4_time += ((Double_t)(t1 - t0)) / CLOCKS_PER_SEC;
539
540 return kBMNSUCCESS;
541}
542
543BmnStatus BmnInnerTrackingRun7::FindTracks_3of4_OnLastGEMStations()
544{
545 vector<BmnHit*> sortedHits[fNStations];
546 SortHits(sortedHits);
547
548 vector<BmnTrack> candidates;
549 vector<BmnTrack> sortedCandidates;
550
551 FindCandidateByThreeStations(5, 6, 7, candidates, sortedHits);
552 FindCandidateByThreeStations(5, 7, 8, candidates, sortedHits);
553 FindCandidateByThreeStations(5, 6, 8, candidates, sortedHits);
554 FindCandidateByThreeStations(6, 7, 8, candidates, sortedHits);
555
556 // cout << "FindTracks_3of4_OnLastGEMStations::candidates = " << candidates.size() << endl;
557
558 for (BmnTrack& cand : candidates) {
559 if (CalculateTrackParams(&cand) == kBMNERROR)
560 continue;
561 // First two GEM
562 MatchHit(&cand, sortedHits[4], kFALSE);
563 MatchHit(&cand, sortedHits[3], kFALSE);
564 // Silicon
565 MatchHit(&cand, sortedHits[2], kFALSE);
566 MatchHit(&cand, sortedHits[1], kFALSE);
567 MatchHit(&cand, sortedHits[0], kFALSE);
568
569 // NDF = (N counts in ZX plane + N counts in ZY plane) - 2 parameters of Line in ZY plane - 3 parameters of
570 // Circle in ZX plane Check it!!!
571 cand.SetNDF(cand.GetNHits() * 2 - 5);
572 }
573
574 BmnFieldMap* field = (BmnFieldMap*)FairRunAna::Instance()->GetField();
575 Bool_t isField = !(field->IsFieldOff());
576 if (isField)
577 TrackUpdateByKalman(candidates);
578 else
579 TrackUpdateByLine(candidates);
580 SortTracks(candidates, sortedCandidates);
581 TrackSelection(sortedCandidates);
582
583 return kBMNSUCCESS;
584}
585
586BmnStatus BmnInnerTrackingRun7::MatchHit(BmnTrack* cand, vector<BmnHit*> sortedHits, Bool_t downstream)
587{
588 clock_t t0 = 0; // clock();
589 // Double_t cutX = 0.45;
590 // Double_t cutY = 0.45;
591 if (sortedHits.size() == 0)
592 return kBMNERROR;
593 FairTrackParam par = downstream ? *(cand->GetParamLast()) : *(cand->GetParamFirst());
594 Double_t minDist = DBL_MAX;
595 // Double_t minDX = DBL_MAX;
596 // Double_t minDY = DBL_MAX;
597 BmnHit* minHit = nullptr;
598 clock_t t0_0 = 0; // clock();
599 fKalman->TGeoTrackPropagate(&par, sortedHits.at(0)->GetZ(), 2212, nullptr, nullptr);
600 clock_t t0_1 = 0; // clock();
601 for (BmnHit* hit : sortedHits) {
602 // Double_t dX = Abs(par.GetX() - hit->GetX());
603 // Double_t dY = Abs(par.GetY() - hit->GetY());
604 Double_t d = Sqrt(Sq(par.GetX() - hit->GetX()) + Sq(par.GetY() - hit->GetY()));
605 // if (dX < cutX && dY < cutY && dX < minDX && dY < minDY) {
606 if (d < fDistCut && d < minDist) {
607 // minDX = dX;
608 // minDY = dY;
609 minDist = d;
610 minHit = hit;
611 }
612 }
613 if (minHit) {
614 cand->AddHit(minHit->GetRefIndex(), minHit);
615 cand->SortHits();
616 Double_t chi2 = 0.0;
617 fKalman->Update(&par, minHit, chi2);
618 cand->SetChi2(cand->GetChi2() + chi2);
619 }
620 if (downstream)
621 cand->SetParamLast(par);
622 else
623 cand->SetParamFirst(par);
624 clock_t t1 = 0; // clock();
625 hitMatchTime += ((Double_t)(t1 - t0)) / CLOCKS_PER_SEC;
626 hitMatchTimeLoop += ((Double_t)(t0_1 - t0_0)) / CLOCKS_PER_SEC;
627
628 return kBMNSUCCESS;
629}
630
631BmnStatus BmnInnerTrackingRun7::FindCandidateByTwoStations(Short_t st0,
632 Short_t st1,
633 vector<BmnTrack>& candidates,
634 vector<BmnHit*>* sortedHits)
635{
636 for (BmnHit* hit1 : sortedHits[st1]) {
637 Double_t y1 = hit1->GetY();
638 Double_t z1 = hit1->GetZ();
639 for (BmnHit* hit0 : sortedHits[st0]) {
640 Double_t y0 = hit0->GetY();
641 Double_t z0 = hit0->GetZ();
642 if (Abs(hit1->GetX() - hit0->GetX()) > 40.0 || Abs(y1 - y0) > 20.0)
643 continue;
644 BmnTrack cand;
645 cand.AddHit(hit1->GetRefIndex(), hit1);
646 cand.AddHit(hit0->GetRefIndex(), hit0);
647 Double_t a = (y1 - y0) / (z1 - z0);
648 Double_t b = (z1 * y0 - y1 * z0) / (z1 - z0);
649 Double_t minDist = DBL_MAX;
650 BmnHit* minHit = nullptr;
651 for (BmnHit* hit : sortedHits[4]) {
652 Double_t d = Abs(hit->GetY() - (a * hit->GetZ() + b));
653 if (d < minDist && d < 0.5) {
654 minDist = d;
655 minHit = hit;
656 }
657 }
658 if (minHit) {
659 cand.AddHit(minHit->GetRefIndex(), minHit);
660 }
661 minDist = DBL_MAX;
662 minHit = nullptr;
663 for (BmnHit* hit : sortedHits[3]) {
664 Double_t d = Abs(hit->GetY() - (a * hit->GetZ() + b));
665 if (d < minDist && d < 0.5) {
666 minDist = d;
667 minHit = hit;
668 }
669 }
670 if (minHit) {
671 cand.AddHit(minHit->GetRefIndex(), minHit);
672 }
673 cand.SortHits();
674 candidates.push_back(cand);
675 }
676 }
677
678 return kBMNSUCCESS;
679}
680
681BmnStatus BmnInnerTrackingRun7::FindTracks_2of2_OnFirstGEMStationsDownstream()
682{
683 vector<BmnHit*> sortedHits[fNStations];
684 SortHits(sortedHits);
685 clock_t t0 = 0; // clock();
686
687 vector<BmnTrack> candidates;
688 vector<BmnTrack> sortedCandidates;
689 for (BmnHit* hit1 : sortedHits[4]) {
690 Double_t y1 = hit1->GetY();
691 Double_t z1 = hit1->GetZ();
692 for (BmnHit* hit0 : sortedHits[3]) {
693 Double_t y0 = hit0->GetY();
694 Double_t z0 = hit0->GetZ();
695 if (Abs(hit1->GetX() - hit0->GetX()) > 40.0 || Abs(y1 - y0) > 10.0)
696 continue;
697 BmnTrack cand;
698 cand.AddHit(hit1->GetRefIndex(), hit1);
699 cand.AddHit(hit0->GetRefIndex(), hit0);
700 Double_t a = (y1 - y0) / (z1 - z0);
701 Double_t b = (z1 * y0 - y1 * z0) / (z1 - z0);
702
703 Double_t minDist = DBL_MAX;
704 BmnHit* minHit = nullptr;
705 for (BmnHit* hit : sortedHits[5]) {
706 Double_t d = Abs(hit->GetY() - (a * hit->GetZ() + b));
707 if (d < minDist && d < fDistCut) {
708 minDist = d;
709 minHit = hit;
710 }
711 }
712 if (minHit) {
713 cand.AddHit(minHit->GetRefIndex(), minHit);
714 cand.SortHits();
715
716 if (CalculateTrackParams(&cand) == kBMNERROR)
717 continue;
718 // NDF = (N counts in ZX plane + N counts in ZY plane) - 2 parameters of Line in ZY plane - 3 parameters
719 // of Circle in ZX plane Check it!!!
720 cand.SetNDF(cand.GetNHits() * 2 - 5);
721
722 // Last GEMs
723 MatchHit(&cand, sortedHits[6], kTRUE);
724 MatchHit(&cand, sortedHits[7], kTRUE);
725 MatchHit(&cand, sortedHits[8], kTRUE);
726
727 // Silicon
728 MatchHit(&cand, sortedHits[2], kFALSE);
729 MatchHit(&cand, sortedHits[1], kFALSE);
730 MatchHit(&cand, sortedHits[0], kFALSE);
731
732 candidates.push_back(cand);
733 } else {
734 minDist = DBL_MAX;
735 minHit = nullptr;
736 for (BmnHit* hit : sortedHits[6]) {
737 Double_t d = Abs(hit->GetY() - (a * hit->GetZ() + b));
738 if (d < minDist && d < fDistCut) {
739 minDist = d;
740 minHit = hit;
741 }
742 }
743 if (minHit) {
744 cand.AddHit(minHit->GetRefIndex(), minHit);
745 cand.SortHits();
746
747 if (CalculateTrackParams(&cand) == kBMNERROR)
748 continue;
749 // NDF = (N counts in ZX plane + N counts in ZY plane) - 2 parameters of Line in ZY plane - 3
750 // parameters of Circle in ZX plane Check it!!!
751 cand.SetNDF(cand.GetNHits() * 2 - 5);
752
753 // Last GEMs
754 MatchHit(&cand, sortedHits[7], kTRUE);
755 MatchHit(&cand, sortedHits[8], kTRUE);
756
757 // Silicon
758 MatchHit(&cand, sortedHits[2], kFALSE);
759 MatchHit(&cand, sortedHits[1], kFALSE);
760 MatchHit(&cand, sortedHits[0], kFALSE);
761
762 candidates.push_back(cand);
763 } else {
764 minDist = DBL_MAX;
765 minHit = nullptr;
766 for (BmnHit* hit : sortedHits[7]) {
767 Double_t d = Abs(hit->GetY() - (a * hit->GetZ() + b));
768 if (d < minDist && d < fDistCut) {
769 minDist = d;
770 minHit = hit;
771 }
772 }
773 if (minHit) {
774 cand.AddHit(minHit->GetRefIndex(), minHit);
775 cand.SortHits();
776
777 if (CalculateTrackParams(&cand) == kBMNERROR)
778 continue;
779 // NDF = (N counts in ZX plane + N counts in ZY plane) - 2 parameters of Line in ZY plane - 3
780 // parameters of Circle in ZX plane Check it!!!
781 cand.SetNDF(cand.GetNHits() * 2 - 5);
782
783 // Last GEMs
784 MatchHit(&cand, sortedHits[8], kTRUE);
785
786 // Silicon
787 MatchHit(&cand, sortedHits[2], kFALSE);
788 MatchHit(&cand, sortedHits[1], kFALSE);
789 MatchHit(&cand, sortedHits[0], kFALSE);
790
791 candidates.push_back(cand);
792 } else {
793 minDist = DBL_MAX;
794 minHit = nullptr;
795 for (BmnHit* hit : sortedHits[8]) {
796 Double_t d = Abs(hit->GetY() - (a * hit->GetZ() + b));
797 if (d < minDist && d < fDistCut) {
798 minDist = d;
799 minHit = hit;
800 }
801 }
802 if (minHit) {
803 cand.AddHit(minHit->GetRefIndex(), minHit);
804 cand.SortHits();
805
806 if (CalculateTrackParams(&cand) == kBMNERROR)
807 continue;
808 // NDF = (N counts in ZX plane + N counts in ZY plane) - 2 parameters of Line in ZY plane -
809 // 3 parameters of Circle in ZX plane Check it!!!
810 cand.SetNDF(cand.GetNHits() * 2 - 5);
811
812 // Silicon
813 MatchHit(&cand, sortedHits[2], kFALSE);
814 MatchHit(&cand, sortedHits[1], kFALSE);
815 MatchHit(&cand, sortedHits[0], kFALSE);
816
817 candidates.push_back(cand);
818 } else {
819 minDist = DBL_MAX;
820 minHit = nullptr;
821 for (BmnHit* hit : sortedHits[2]) {
822 Double_t d = Abs(hit->GetY() - (a * hit->GetZ() + b));
823 if (d < minDist && d < fDistCut) {
824 minDist = d;
825 minHit = hit;
826 }
827 }
828 if (minHit) {
829 cand.AddHit(minHit->GetRefIndex(), minHit);
830 cand.SortHits();
831
832 if (CalculateTrackParams(&cand) == kBMNERROR)
833 continue;
834 // NDF = (N counts in ZX plane + N counts in ZY plane) - 2 parameters of Line in ZY
835 // plane - 3 parameters of Circle in ZX plane Check it!!!
836 cand.SetNDF(cand.GetNHits() * 2 - 5);
837
838 // Silicon
839 MatchHit(&cand, sortedHits[1], kFALSE);
840 MatchHit(&cand, sortedHits[0], kFALSE);
841
842 candidates.push_back(cand);
843 } else {
844 minDist = DBL_MAX;
845 minHit = nullptr;
846 for (BmnHit* hit : sortedHits[1]) {
847 Double_t d = Abs(hit->GetY() - (a * hit->GetZ() + b));
848 if (d < minDist && d < fDistCut) {
849 minDist = d;
850 minHit = hit;
851 }
852 }
853 if (minHit) {
854 cand.AddHit(minHit->GetRefIndex(), minHit);
855 cand.SortHits();
856
857 if (CalculateTrackParams(&cand) == kBMNERROR)
858 continue;
859 // NDF = (N counts in ZX plane + N counts in ZY plane) - 2 parameters of Line in ZY
860 // plane - 3 parameters of Circle in ZX plane Check it!!!
861 cand.SetNDF(cand.GetNHits() * 2 - 5);
862
863 // Silicon
864 MatchHit(&cand, sortedHits[0], kFALSE);
865
866 candidates.push_back(cand);
867 }
868 }
869 }
870 }
871 }
872 }
873 }
874 }
875 clock_t t1 = 0; // clock();
876 construct_2of4_time += ((Double_t)(t1 - t0)) / CLOCKS_PER_SEC;
877 // cout << "FindTracks_2of2_OnFirstGEMstStations::candidates = " << candidates.size() << endl;
878
879 BmnFieldMap* field = (BmnFieldMap*)FairRunAna::Instance()->GetField();
880 Bool_t isField = !(field->IsFieldOff());
881 if (isField)
882 TrackUpdateByKalman(candidates);
883 else
884 TrackUpdateByLine(candidates);
885 SortTracks(candidates, sortedCandidates);
886 TrackSelection(sortedCandidates);
887
888 return kBMNSUCCESS;
889}
890
891BmnStatus BmnInnerTrackingRun7::FindTracks_2of2_OnFirstGEMStationsUpstream()
892{
893 vector<BmnHit*> sortedHits[fNStations];
894 SortHits(sortedHits);
895 clock_t t0 = 0; // clock();
896 vector<BmnTrack> candidates;
897 vector<BmnTrack> sortedCandidates;
898 for (BmnHit* hit1 : sortedHits[4]) {
899 Double_t y1 = hit1->GetY();
900 Double_t z1 = hit1->GetZ();
901 for (BmnHit* hit0 : sortedHits[3]) {
902 Double_t y0 = hit0->GetY();
903 Double_t z0 = hit0->GetZ();
904 if (Abs(hit1->GetX() - hit0->GetX()) > 40.0 || Abs(y1 - y0) > 10.0)
905 continue;
906 BmnTrack cand;
907 cand.AddHit(hit1->GetRefIndex(), hit1);
908 cand.AddHit(hit0->GetRefIndex(), hit0);
909 Double_t a = (y1 - y0) / (z1 - z0);
910 Double_t b = (z1 * y0 - y1 * z0) / (z1 - z0);
911 Double_t minDist = DBL_MAX;
912 BmnHit* minHit = nullptr;
913 for (BmnHit* hit : sortedHits[2]) {
914 Double_t d = Abs(hit->GetY() - (a * hit->GetZ() + b));
915 if (d < minDist && d < fDistCut) {
916 minDist = d;
917 minHit = hit;
918 }
919 }
920 if (minHit)
921 cand.AddHit(minHit->GetRefIndex(), minHit);
922 else {
923 minDist = DBL_MAX;
924 minHit = nullptr;
925 for (BmnHit* hit : sortedHits[1]) {
926 Double_t d = Abs(hit->GetY() - (a * hit->GetZ() + b));
927 if (d < minDist && d < fDistCut) {
928 minDist = d;
929 minHit = hit;
930 }
931 }
932 if (minHit)
933 cand.AddHit(minHit->GetRefIndex(), minHit);
934 }
935 cand.SortHits();
936 candidates.push_back(cand);
937 }
938 }
939 clock_t t1 = 0; // clock();
940 construct_2of4_time += ((Double_t)(t1 - t0)) / CLOCKS_PER_SEC;
941 // cout << "FindTracks_2of2_OnFirstGEMstStations::candidates = " << candidates.size() << endl;
942
943 for (BmnTrack& cand : candidates) {
944 if (CalculateTrackParams(&cand) == kBMNERROR)
945 continue;
946
947 // //Last GEMs
948 MatchHit(&cand, sortedHits[5], kTRUE);
949 MatchHit(&cand, sortedHits[6], kTRUE);
950 MatchHit(&cand, sortedHits[7], kTRUE);
951 MatchHit(&cand, sortedHits[8], kTRUE);
952
953 // Silicon
954 MatchHit(&cand, sortedHits[2], kFALSE);
955 MatchHit(&cand, sortedHits[1], kFALSE);
956 MatchHit(&cand, sortedHits[0], kFALSE);
957
958 // NDF = (N counts in ZX plane + N counts in ZY plane) - 2 parameters of Line in ZY plane - 3 parameters of
959 // Circle in ZX plane Check it!!!
960 cand.SetNDF(cand.GetNHits() * 2 - 5);
961 }
962
963 BmnFieldMap* field = (BmnFieldMap*)FairRunAna::Instance()->GetField();
964 Bool_t isField = !(field->IsFieldOff());
965 if (isField)
966 TrackUpdateByKalman(candidates);
967 else
968 TrackUpdateByLine(candidates);
969 SortTracks(candidates, sortedCandidates);
970 TrackSelection(sortedCandidates);
971
972 return kBMNSUCCESS;
973}
974
975BmnStatus BmnInnerTrackingRun7::SortTracks(vector<BmnTrack>& inTracks, vector<BmnTrack>& sortedTracks)
976{
977 clock_t t0 = 0; // clock();
978 const Int_t n = fNStations - fNHitsCut + 1; // 6 for geometry 2018 (4, 5, 6, 7, 8, 9)
979 multimap<Float_t, Int_t> sortedMap[n]; // array of map<Chi2,trIdx>. Each element of array corresponds fixed number
980 // of hits on track (4, 5, 6)
981 for (size_t iTr = 0; iTr < inTracks.size(); ++iTr) {
982 if (inTracks.at(iTr).GetNHits() < fNHitsCut || inTracks.at(iTr).GetNHits() > fNStations)
983 continue;
984 if (inTracks.at(iTr).GetChi2() / inTracks.at(iTr).GetNDF() > fChiSquareCut)
985 continue;
986 sortedMap[inTracks.at(iTr).GetNHits() - fNHitsCut].insert(
987 pair<Float_t, Int_t>(inTracks.at(iTr).GetChi2() / inTracks.at(iTr).GetNDF(), iTr));
988 }
989
990 for (Int_t i = n - 1; i >= 0; i--) {
991 for (auto it : sortedMap[i])
992 sortedTracks.push_back(inTracks.at(it.second));
993 sortedMap[i].clear();
994 }
995
996 // multimap <Float_t, Int_t> sortedTracksMap; //map<Chi2,trIdx>
997 // for (Int_t iTr = 0; iTr < inTracks.size(); ++iTr)
998 // sortedTracksMap.insert(pair<Float_t, Int_t>(inTracks.at(iTr).GetChi2() / inTracks.at(iTr).GetNDF(), iTr));
999 //
1000 // for (auto it : sortedTracksMap)
1001 // sortedTracks.push_back(inTracks.at(it.second));
1002 clock_t t1 = 0; // clock();
1003 sortTrackTime += ((Double_t)(t1 - t0)) / CLOCKS_PER_SEC;
1004 return kBMNSUCCESS;
1005}
1006
1007void BmnInnerTrackingRun7::Finish()
1008{
1009 ofstream outFile;
1010 outFile.open("QA/timing.txt");
1011 outFile << "Track Finder Time: " << workTime << endl;
1012 if (fVerbose > 0) {
1013 cout << "hitMatchTime: " << hitMatchTime << endl;
1014 cout << "hitMatchTimeLoop: " << hitMatchTimeLoop << endl;
1015 cout << "construct_4of4_time: " << construct_4of4_time << endl;
1016 cout << "construct_3of4_time: " << construct_3of4_time << endl;
1017 cout << "construct_2of4_time: " << construct_2of4_time << endl;
1018 cout << "CalculateTrackParTime: " << CalculateTrackParTime << endl;
1019 cout << "SortHitsTime: " << SortHitsTime << endl;
1020 cout << "trackUpdateTime: " << trackUpdateTime << endl;
1021 cout << "sortTrackTime: " << sortTrackTime << endl;
1022 cout << "trackSelectionTime: " << trackSelectionTime << endl;
1023 }
1024 cout << "Work time of the GEM tracking: " << workTime << endl;
1025 return;
1026}
1027
1028BmnStatus BmnInnerTrackingRun7::SortHits(vector<BmnHit*>* sortedHits)
1029{
1030 clock_t t0 = 0; // clock();
1031 for (Int_t iHit = 0; iHit < fHitsArray->GetEntriesFast(); ++iHit) {
1032 BmnHit* hit = (BmnHit*)fHitsArray->At(iHit);
1033 if (!hit)
1034 continue;
1035 if (hit->IsUsed())
1036 continue;
1037 hit->SetRefIndex(iHit);
1038 Int_t station = hit->GetStation();
1039 // if (hit->GetX() > fHitXCutMax[station] || hit->GetX() < fHitXCutMin[station]) continue;
1040 if (hit->GetY() > fHitYCutMax[station] || hit->GetY() < fHitYCutMin[station])
1041 continue;
1042 sortedHits[station].push_back(hit);
1043 }
1044 clock_t t1 = 0; // clock();
1045 SortHitsTime += ((Double_t)(t1 - t0)) / CLOCKS_PER_SEC;
1046
1047 return kBMNSUCCESS;
1048}
1049
1050BmnStatus BmnInnerTrackingRun7::TrackUpdateByLine(vector<BmnTrack>& cands)
1051{
1052 for (BmnTrack& cand : cands) {
1053 Double_t Tx = LineFit((BmnTrack*)&cand, fHitsArray, "ZX").X();
1054 Double_t chiX = LineFit((BmnTrack*)&cand, fHitsArray, "ZX").Z();
1055 Double_t Ty = LineFit((BmnTrack*)&cand, fHitsArray, "ZY").X();
1056 Double_t chiY = LineFit((BmnTrack*)&cand, fHitsArray, "ZY").Z();
1057
1058 cand.SetChi2((chiX - chiY) > 0. ? chiX : chiY);
1059
1060 cand.GetParamFirst()->SetTx(Tx);
1061 cand.GetParamFirst()->SetTy(Ty);
1062
1063 cand.GetParamLast()->SetTx(Tx);
1064 cand.GetParamLast()->SetTy(Ty);
1065 }
1066
1067 return kBMNSUCCESS;
1068}
1069
1070BmnStatus BmnInnerTrackingRun7::TrackUpdateByKalman(vector<BmnTrack>& cands)
1071{
1072 clock_t t0 = 0; // clock();
1073 for (BmnTrack& cand : cands) {
1074 if (cand.GetNHits() < fNHitsCut)
1075 continue;
1076 FairTrackParam par = *(cand.GetParamFirst());
1077 // par.SetQp(1.0/8.0);
1078 Double_t chiTot = 0.0;
1079 for (Int_t iHit = 0; iHit < cand.GetNHits(); ++iHit) {
1080 BmnHit* hit = (BmnHit*)fHitsArray->At(cand.GetHitIndex(iHit));
1081 Double_t chi = 0.0;
1082 fKalman->TGeoTrackPropagate(&par, hit->GetZ(), 2212, nullptr, nullptr);
1083 fKalman->Update(&par, hit, chi);
1084 }
1085 cand.SetParamLast(par);
1086 for (Int_t iHit = cand.GetNHits() - 1; iHit >= 0; iHit--) {
1087 BmnHit* hit = (BmnHit*)fHitsArray->At(cand.GetHitIndex(iHit));
1088 Double_t chi = 0.0;
1089 fKalman->TGeoTrackPropagate(&par, hit->GetZ(), 2212, nullptr, nullptr);
1090 fKalman->Update(&par, hit, chi);
1091 chiTot += chi;
1092 }
1093 cand.SetParamFirst(par);
1094 cand.SetChi2(chiTot);
1095 }
1096 clock_t t1 = 0; // clock();
1097 trackUpdateTime += ((Double_t)(t1 - t0)) / CLOCKS_PER_SEC;
1098
1099 return kBMNSUCCESS;
1100}
1101
1102BmnStatus BmnInnerTrackingRun7::TrackSelection(vector<BmnTrack>& sortedTracks)
1103{
1104 clock_t t0 = 0; // clock();
1105
1106 // Switch on/off station skip in tracks
1107
1108 // for (Int_t iTr = 0; iTr < sortedTracks.size(); ++iTr) {
1109 // BmnTrack* tr = &(sortedTracks[iTr]);
1110 // if (tr->GetFlag() == 666 || !IsParCorrect(tr->GetParamFirst()) || !IsParCorrect(tr->GetParamLast()))
1111 // continue; Bool_t badTrack = kFALSE; BmnHit* hit0 = (BmnHit*)fHitsArray->At(tr->GetHitIndex(0)); Int_t stPrev
1112 // = hit0->GetStation(); for (Int_t iHit = 1; iHit < tr->GetNHits(); ++iHit) {
1113 // BmnHit* hit = (BmnHit*)fHitsArray->At(tr->GetHitIndex(iHit));
1114 // Int_t st = hit->GetStation();
1115 // Int_t dSt = st - stPrev - 1;
1116 // stPrev = st;
1117 // if (dSt == 1) { //number of skipped stations in track
1118 // badTrack = kTRUE;
1119 // break;
1120 // }
1121 // }
1122 // if (badTrack) tr->SetFlag(666);
1123 // }
1124 CheckSharedHits(sortedTracks);
1125 for (size_t iTr = 0; iTr < sortedTracks.size(); ++iTr) {
1126 BmnTrack tr = sortedTracks[iTr];
1127 if (tr.GetFlag() == 666 || !IsParCorrect(tr.GetParamFirst()) || !IsParCorrect(tr.GetParamLast()))
1128 continue;
1129 BmnTrack gemTr;
1130 BmnTrack silTr;
1131 BmnGlobalTrack globTr;
1132
1133 for (Int_t iHit = 0; iHit < tr.GetNHits(); ++iHit) {
1134 BmnHit* hit = (BmnHit*)fHitsArray->At(tr.GetHitIndex(iHit));
1135 globTr.AddHit(hit->GetRefIndex(), hit);
1136 DetectorId detId = hit->GetDetId();
1137 if (detId == kSILICON)
1138 silTr.AddHit(hit->GetIndex(), hit);
1139 else if (detId == kGEM)
1140 gemTr.AddHit(hit->GetIndex(), hit);
1141 else
1142 continue;
1143 }
1144 globTr.SortHits();
1145
1146 silTr.SortHits();
1147 gemTr.SortHits();
1148
1149 if (silTr.GetNHits() > 0) {
1150 new ((*fSilTracksArray)[fSilTracksArray->GetEntriesFast()]) BmnTrack(silTr);
1151 globTr.SetSilTrackIndex(fSilTracksArray->GetEntriesFast() - 1);
1152 }
1153 if (gemTr.GetNHits() > 0) {
1154 gemTr.SetParamFirst(*(tr.GetParamFirst()));
1155 gemTr.SetParamLast(*(tr.GetParamLast()));
1156 new ((*fGemTracksArray)[fGemTracksArray->GetEntriesFast()]) BmnTrack(gemTr);
1157 globTr.SetGemTrackIndex(fGemTracksArray->GetEntriesFast() - 1);
1158 }
1159 globTr.SetLength(CalculateLength(&tr));
1160 globTr.SetParamFirst(*(tr.GetParamFirst()));
1161 globTr.SetParamLast(*(tr.GetParamLast()));
1162 globTr.SetNHits(tr.GetNHits());
1163 globTr.SetNDF(tr.GetNDF());
1164 globTr.SetChi2(tr.GetChi2());
1165 new ((*fGlobTracksArray)[fGlobTracksArray->GetEntriesFast()]) BmnGlobalTrack(globTr);
1166 SetHitsUsing(&tr, kTRUE);
1167 if (!fIsTarget && iTr == 0)
1168 break;
1169 }
1170 clock_t t1 = 0; // clock();
1171 trackSelectionTime += ((Double_t)(t1 - t0)) / CLOCKS_PER_SEC;
1172
1173 return kBMNSUCCESS;
1174}
1175
1176void BmnInnerTrackingRun7::SetHitsUsing(BmnTrack* tr, Bool_t use)
1177{
1178 for (Int_t i = 0; i < tr->GetNHits(); ++i) {
1179 BmnHit* hit = (BmnHit*)fHitsArray->At(tr->GetHitIndex(i));
1180 if (!hit)
1181 continue;
1182 hit->SetUsing(use);
1183 // for (Int_t iLink = 0; iLink < hit->GetDigitNumberMatch().GetNofLinks(); ++iLink) {
1184 // BmnLink link = hit->GetDigitNumberMatch().GetLink(iLink);
1185 // Int_t idxI = link.GetIndex();
1186 // for (Int_t j = 0; j < fHitsArray->GetEntriesFast(); ++j) {
1187 // BmnHit* hitJ = (BmnHit*)fHitsArray->At(j);
1188 // if (!hitJ) continue;
1189 // if (hitJ->IsUsed()) continue;
1190 // for (Int_t iLink = 0; iLink < hit->GetDigitNumberMatch().GetNofLinks(); ++iLink) {
1191 // BmnLink linkJ = hit->GetDigitNumberMatch().GetLink(iLink);
1192 // Int_t idxJ = linkJ.GetIndex();
1193 // if (idxI == idxJ) {
1194 // hitJ->SetUsing(use);
1195 // break;
1196 // }
1197 // }
1198 // }
1199 // }
1200 }
1201}
1202
1203BmnStatus BmnInnerTrackingRun7::CalcCovMatrix(BmnTrack* tr)
1204{
1205 const UInt_t nHits = tr->GetNHits();
1206 Double_t chi2circ = 0.0;
1207 TVector3 CircParZX = CircleFit(tr, fHitsArray, chi2circ);
1208
1209 Double_t Xc = CircParZX.Y(); // x-coordinate of fit-circle center
1210 Double_t Zc = CircParZX.X(); // z-coordinate of fit-circle center
1211 fField = FairRunAna::Instance()->GetField();
1212 const Double_t B = (LineFit(tr, fHitsArray, "ZY")).X(); // angle coefficient for helicoid
1213
1214 Double_t Q = (Xc > 0) ? +1 : -1;
1215
1216 Double_t sumX(0.0), sumY(0.0), sumTx(0.0), sumTy(0.0), sumQP(0.0);
1217 Double_t sumXX(0.0), sumXY(0.0), sumXTx(0.0), sumXTy(0.0), sumXQP(0.0);
1218 Double_t sumYY(0.0), sumYTx(0.0), sumYTy(0.0), sumYQP(0.0);
1219 Double_t sumTxTx(0.0), sumTxTy(0.0), sumTxQP(0.0);
1220 Double_t sumTyTy(0.0), sumTyQP(0.0);
1221 Double_t sumQPQP(0.0);
1222
1223 for (UInt_t i = 0; i < nHits; ++i) {
1224 BmnHit* hit = (BmnHit*)fHitsArray->At(tr->GetHitIndex(i));
1225 if (!hit)
1226 continue;
1227 Double_t Xi = hit->GetX();
1228 Double_t Yi = hit->GetY();
1229 Double_t Zi = hit->GetZ();
1230
1231 Double_t Txi = -1.0 * (Zi - Zc) / (Xi - Xc);
1232 Double_t Tyi = B;
1233 Double_t Ri = Sqrt(Sq(Xi - Xc) + Sq(Zi - Zc));
1234 Double_t QPi = Q / (0.0003 * Abs(fField->GetBy(Xi, Yi, Zi)) * Ri);
1235
1236 sumX += Xi;
1237 sumY += Yi;
1238 sumTx += Txi;
1239 sumTy += Tyi;
1240 sumQP += QPi;
1241 sumXX += Xi * Xi;
1242 sumXY += Xi * Yi;
1243 sumXTx += Xi * Txi;
1244 sumXTy += Xi * Tyi;
1245 sumXQP += Xi * QPi;
1246 sumYY += Yi * Yi;
1247 sumYTx += Yi * Txi;
1248 sumYTy += Yi * Tyi;
1249 sumYQP += Yi * QPi;
1250 sumTxTx += Txi * Txi;
1251 sumTxTy += Txi * Tyi;
1252 sumTxQP += Txi * QPi;
1253 sumTyTy += Tyi * Tyi;
1254 sumTyQP += Tyi * QPi;
1255 sumQPQP += QPi * QPi;
1256 }
1257
1258 Double_t meanX = sumX / nHits;
1259 Double_t meanY = sumY / nHits;
1260 Double_t meanTx = sumTx / nHits;
1261 Double_t meanTy = sumTy / nHits;
1262 Double_t meanQP = sumQP / nHits;
1263
1264 Double_t Cov_X_X = sumXX / nHits - Sq(meanX);
1265 Double_t Cov_X_Y = sumXY / nHits - meanX * meanY;
1266 Double_t Cov_X_Tx = sumXTx / nHits - meanX * meanTx;
1267 Double_t Cov_X_Ty = sumXTy / nHits - meanX * meanTy;
1268 Double_t Cov_X_Qp = sumXQP / nHits - meanX * meanQP;
1269 Double_t Cov_Y_Y = sumYY / nHits - Sq(meanY);
1270 Double_t Cov_Y_Tx = sumYTx / nHits - meanY * meanTx;
1271 Double_t Cov_Y_Ty = sumYTy / nHits - meanY * meanTy;
1272 Double_t Cov_Y_Qp = sumYQP / nHits - meanY * meanQP;
1273 Double_t Cov_Tx_Tx = sumTxTx / nHits - Sq(meanTx);
1274 Double_t Cov_Tx_Ty = sumTxTy / nHits - meanTx * meanTy;
1275 Double_t Cov_Tx_Qp = sumTxQP / nHits - meanTx * meanQP;
1276 Double_t Cov_Ty_Ty = sumTyTy / nHits - Sq(meanTy);
1277 Double_t Cov_Ty_Qp = sumTyQP / nHits - meanTy * meanQP;
1278 Double_t Cov_Qp_Qp = sumQPQP / nHits - Sq(meanQP);
1279
1280 FairTrackParam par;
1281 par.SetCovariance(0, 0, Cov_X_X);
1282 par.SetCovariance(0, 1, Cov_X_Y);
1283 par.SetCovariance(0, 2, Cov_X_Tx);
1284 par.SetCovariance(0, 3, Cov_X_Ty);
1285 par.SetCovariance(0, 4, Cov_X_Qp);
1286 par.SetCovariance(1, 1, Cov_Y_Y);
1287 par.SetCovariance(1, 2, Cov_Y_Tx);
1288 par.SetCovariance(1, 3, Cov_Y_Ty);
1289 par.SetCovariance(1, 4, Cov_Y_Qp);
1290 par.SetCovariance(2, 2, Cov_Tx_Tx);
1291 par.SetCovariance(2, 3, Cov_Tx_Ty);
1292 par.SetCovariance(2, 4, Cov_Tx_Qp);
1293 par.SetCovariance(3, 3, Cov_Ty_Ty);
1294 par.SetCovariance(3, 4, Cov_Ty_Qp);
1295 par.SetCovariance(4, 4, Cov_Qp_Qp);
1296
1297 tr->SetParamFirst(par);
1298 tr->SetParamLast(par);
1299
1300 return kBMNSUCCESS;
1301}
1302
1303BmnStatus BmnInnerTrackingRun7::CalculateTrackParams(BmnTrack* tr)
1304{
1305
1306 BmnFieldMap* field = (BmnFieldMap*)FairRunAna::Instance()->GetField();
1307 Bool_t isField = !(field->IsFieldOff());
1308 clock_t t0 = 0; // clock();
1309 // Estimation of track parameters for events with magnetic field
1310 const UInt_t nHits = tr->GetNHits();
1311 if (nHits < 3)
1312 return kBMNERROR;
1313 TVector3 lineParZY = LineFit(tr, fHitsArray, "ZY");
1314 if (lineParZY.Z() > 1)
1315 return kBMNERROR; // cout << "chi2_lineFit = " << lineParZY.Z() << endl;
1316 // tr->SetNDF(nHits - (isField ? 3 : 2));
1317 const Double_t B = lineParZY.X(); // angle coefficient for helicoid
1318
1319 Double_t Tx_first = CalcTx((BmnHit*)fHitsArray->At(tr->GetHitIndex(0)), (BmnHit*)fHitsArray->At(tr->GetHitIndex(1)),
1320 (BmnHit*)fHitsArray->At(tr->GetHitIndex(2)));
1321 Double_t Tx_last =
1322 CalcTx((BmnHit*)fHitsArray->At(tr->GetHitIndex(nHits - 1)), (BmnHit*)fHitsArray->At(tr->GetHitIndex(nHits - 2)),
1323 (BmnHit*)fHitsArray->At(tr->GetHitIndex(nHits - 3)));
1324
1325 if (isField)
1326 CalcCovMatrix(tr);
1327 TVector3 firstPos;
1328 TVector3 lastPos;
1329 ((BmnHit*)fHitsArray->At(tr->GetHitIndex(0)))->Position(firstPos);
1330 ((BmnHit*)fHitsArray->At(tr->GetHitIndex(nHits - 1)))->Position(lastPos);
1331 tr->GetParamFirst()->SetPosition(firstPos);
1332 tr->GetParamFirst()->SetTx(Tx_first);
1333 tr->GetParamFirst()->SetTy(B);
1334 tr->GetParamLast()->SetPosition(lastPos);
1335 tr->GetParamLast()->SetTx(Tx_last);
1336 tr->GetParamLast()->SetTy(B);
1337 Bool_t doSimple = (nHits == 3) ? kTRUE : kFALSE;
1338 Double_t QP = isField ? CalcQp(tr, doSimple) : 0.0;
1339 tr->GetParamFirst()->SetQp(QP);
1340 tr->GetParamLast()->SetQp(QP);
1341 clock_t t1 = 0; // clock();
1342 CalculateTrackParTime += ((Double_t)(t1 - t0)) / CLOCKS_PER_SEC;
1343 return kBMNSUCCESS;
1344}
1345
1346TVector2 BmnInnerTrackingRun7::CalcMeanSigma(vector<Double_t> QpSegm)
1347{
1348 Double_t QpSum = 0.;
1349 for (size_t iSegm = 0; iSegm < QpSegm.size(); iSegm++)
1350 QpSum += QpSegm[iSegm];
1351
1352 Double_t QpMean = QpSum / QpSegm.size();
1353
1354 Double_t sqSigmaSum = 0.;
1355 for (size_t iSegm = 0; iSegm < QpSegm.size(); iSegm++)
1356 sqSigmaSum += Sq(QpSegm[iSegm] - QpMean);
1357
1358 return TVector2(QpMean, Sqrt(sqSigmaSum / QpSegm.size()));
1359}
1360
1361Double_t BmnInnerTrackingRun7::CalcQp(BmnTrack* track, Bool_t doSimple)
1362{
1363 // Bool_t fRobustRefit = kTRUE;
1364 // Bool_t fSimpleRefit = kFALSE;
1365 vector<BmnHit*> hits;
1366
1367 for (Int_t iHit = 0; iHit < track->GetNHits(); iHit++)
1368 hits.push_back((BmnHit*)fHitsArray->At(track->GetHitIndex(iHit)));
1369
1370 Int_t kNSegm = track->GetNHits() - 2;
1371
1372 Double_t QpRefit = 0.;
1373 vector<Double_t> QpSegmBefore;
1374
1375 // Get q/p info from all track segments
1376 for (Int_t iHit = 0; iHit < kNSegm; iHit++) {
1377 BmnHit* first = hits[iHit];
1378 BmnHit* second = hits[iHit + 1];
1379 BmnHit* third = hits[iHit + 2];
1380
1381 TVector3 CircParZX = CircleBy3Hit(first, second, third);
1382 Double_t R = CircParZX.Z();
1383 Double_t Xc = CircParZX.Y();
1384 Double_t Zc = CircParZX.X();
1385
1386 Double_t Q = (Xc > 0) ? +1. : -1.;
1387 Double_t S = 0.0003
1388 * (Abs(fField->GetBy(third->GetX(), third->GetY(), third->GetZ()))
1389 + Abs(fField->GetBy(second->GetX(), second->GetY(), second->GetZ()))
1390 + Abs(fField->GetBy(first->GetX(), first->GetY(), first->GetZ())))
1391 / 3.;
1392
1393 Double_t Pt = S * R; // actually Pt/Q, but it doesn't matter
1394 Double_t fX = first->GetX();
1395 Double_t fZ = first->GetZ();
1396
1397 Double_t h = -1.0;
1398
1399 Double_t Tx_first = h * (fZ - Zc) / (fX - Xc);
1400 TVector3 lineParZY = LineFit(track, fHitsArray, "ZY");
1401 const Double_t B = lineParZY.X(); // angle coefficient for helicoid
1402 Double_t Ty_first = B; // / (fX - Xc);
1403
1404 Double_t Pz = Pt / Sqrt(1 + Sq(Tx_first));
1405 Double_t Px = Tx_first * Pz;
1406 Double_t Py = Ty_first * Pz;
1407 Double_t P = Sqrt(Sq(Px) + Sq(Py) + Sq(Pz));
1408 Double_t QP = Q / P;
1409
1410 QpSegmBefore.push_back(QP);
1411 }
1412
1413 // Non-robust (simple) refit when segments with bad q/p are not taken into account
1414 if (doSimple) {
1415 vector<Double_t> QpSegmAfter;
1416 while (kTRUE) {
1417 TVector2 meanSig = CalcMeanSigma(QpSegmBefore);
1418 Double_t mean = meanSig.X();
1419 Double_t sigma = meanSig.Y();
1420 if (std::isnan(sigma) && fVerbose == 3) {
1421 cout << "Bad refit convergence for track segment!!" << endl;
1422 return kBMNERROR;
1423 }
1424
1425 for (size_t iSegm = 0; iSegm < QpSegmBefore.size(); iSegm++)
1426 if (Abs(QpSegmBefore[iSegm] - mean) - sigma <= 0.001) // Топорное сравнение FIXME
1427 QpSegmAfter.push_back(QpSegmBefore[iSegm]);
1428
1429 if (QpSegmAfter.size() == QpSegmBefore.size()) {
1430 QpRefit = mean;
1431 break;
1432 } else {
1433 QpSegmBefore.clear();
1434 QpSegmBefore.resize(0);
1435
1436 for (size_t iSegm = 0; iSegm < QpSegmAfter.size(); iSegm++)
1437 QpSegmBefore.push_back(QpSegmAfter[iSegm]);
1438
1439 QpSegmAfter.clear();
1440 QpSegmAfter.resize(0);
1441 }
1442 }
1443 }
1444
1445 // Robust refit with use of Tukey weights calculation algorithm
1446 // if (fRobustRefit) {
1447 else
1448 {
1449 for (size_t iEle = 0; iEle < QpSegmBefore.size(); iEle++)
1450 QpRefit += QpSegmBefore[iEle];
1451
1452 QpRefit /= QpSegmBefore.size();
1453
1454 vector<Double_t> d = dist(QpSegmBefore, QpRefit);
1455
1456 Double_t sigma = 0.;
1457 for (size_t i = 0; i < QpSegmBefore.size(); i++)
1458 sigma += (QpSegmBefore[i] - QpRefit) * (QpSegmBefore[i] - QpRefit);
1459 sigma = Sqrt(sigma / QpSegmBefore.size());
1460
1461 vector<Double_t> w = W(d, sigma);
1462 sigma = Sigma(d, w);
1463
1464 const Int_t kNIter = 20; // FIXME
1465 for (Int_t iIter = 1; iIter < kNIter; iIter++) {
1466 QpRefit = Mu(QpSegmBefore, w);
1467 d = dist(QpSegmBefore, QpRefit);
1468 w = W(d, sigma);
1469 sigma = Sigma(d, w);
1470 }
1471 }
1472
1473 return QpRefit;
1474}
1475
1476Double_t BmnInnerTrackingRun7::CalculateLength(BmnTrack* tr)
1477{
1478 if (!tr)
1479 return 0.0;
1480
1481 vector<Double_t> X, Y, Z;
1482 for (Int_t iGem = 0; iGem < tr->GetNHits(); iGem++) {
1483 BmnHit* hit = (BmnHit*)fHitsArray->At(tr->GetHitIndex(iGem));
1484 if (!hit)
1485 continue;
1486 X.push_back(hit->GetX());
1487 Y.push_back(hit->GetY());
1488 Z.push_back(hit->GetZ());
1489 }
1490 // Calculate distances between hits
1491 Double_t length = 0.;
1492 for (size_t i = 0; i < X.size() - 1; i++) {
1493 Double_t dX = X[i] - X[i + 1];
1494 Double_t dY = Y[i] - Y[i + 1];
1495 Double_t dZ = Z[i] - Z[i + 1];
1496 length += Sqrt(dX * dX + dY * dY + dZ * dZ);
1497 }
1498 tr->SetLength(length);
1499 return length;
1500}
1501
1502BmnStatus BmnInnerTrackingRun7::CheckSharedHits(vector<BmnTrack>& sortedTracks)
1503{
1504 set<Int_t> hitsId;
1505
1506 const Int_t kNSharedHits = 1; // fSteering->GetNSharedHits();
1507
1508 for (size_t iTr = 0; iTr < sortedTracks.size(); ++iTr) {
1509 BmnTrack* tr = &(sortedTracks.at(iTr));
1510 if (tr->GetFlag() == 666)
1511 continue;
1512
1513 Int_t nofSharedHits = 0;
1514 Int_t nofHits = tr->GetNHits();
1515 for (Int_t iHit = 0; iHit < nofHits; iHit++)
1516 if (hitsId.find(tr->GetHitIndex(iHit)) != hitsId.end()) {
1517 nofSharedHits++;
1518 if (nofSharedHits > kNSharedHits) {
1519 tr->SetFlag(666);
1520 break;
1521 }
1522 }
1523 if (tr->GetFlag() == 666)
1524 continue;
1525
1526 for (Int_t iHit = 0; iHit < nofHits; iHit++)
1527 hitsId.insert(tr->GetHitIndex(iHit));
1528 }
1529 hitsId.clear();
1530
1531 return kBMNSUCCESS;
1532}
1533
1534BmnStatus BmnInnerTrackingRun7::DrawHits()
1535{
1536 TH2F* h_HitsZX = new TH2F("h_HitsZX", "h_HitsZX", 400, 0.0, 200.0, 400, -100.0, 100.0);
1537 TH2F* h_HitsZY = new TH2F("h_HitsZY", "h_HitsZY", 400, 0.0, 200.0, 400, -10.0, 50.0);
1538 for (Int_t i = 0; i < fHitsArray->GetEntriesFast(); ++i) {
1539 BmnHit* hit = (BmnHit*)fHitsArray->At(i);
1540 h_HitsZX->Fill(hit->GetZ(), hit->GetX());
1541 h_HitsZY->Fill(hit->GetZ(), hit->GetY());
1542 }
1543
1544 TCanvas* c = new TCanvas("c", "c", 1000, 1000);
1545 c->Divide(1, 2);
1546 c->cd(1);
1547 h_HitsZX->SetMarkerStyle(8);
1548 h_HitsZX->SetMarkerSize(0.7);
1549 h_HitsZX->SetMarkerColor(kRed);
1550 h_HitsZX->Draw("P");
1551
1552 c->cd(2);
1553 h_HitsZY->SetMarkerStyle(8);
1554 h_HitsZY->SetMarkerSize(0.7);
1555 h_HitsZY->SetMarkerColor(kRed);
1556 h_HitsZY->Draw("P");
1557
1558 // if (fGlobTracksArray->GetEntriesFast() <= 0) return kBMNERROR;
1559
1560 for (Int_t iTr = 0; iTr < fGlobTracksArray->GetEntriesFast(); ++iTr) {
1561 BmnGlobalTrack* glTrack = (BmnGlobalTrack*)fGlobTracksArray->At(iTr);
1562
1563 BmnTrack* gemTr = nullptr;
1564 BmnTrack* silTr = nullptr;
1565
1566 Double_t xPrev;
1567 Double_t yPrev;
1568 Double_t zPrev;
1569
1570 if (glTrack->GetSilTrackIndex() != -1) {
1571 silTr = (BmnTrack*)fSilTracksArray->UncheckedAt(glTrack->GetSilTrackIndex());
1572
1573 xPrev = ((BmnHit*)fSilHitsArray->At(silTr->GetHitIndex(0)))->GetX();
1574 yPrev = ((BmnHit*)fSilHitsArray->At(silTr->GetHitIndex(0)))->GetY();
1575 zPrev = ((BmnHit*)fSilHitsArray->At(silTr->GetHitIndex(0)))->GetZ();
1576
1577 for (Int_t iHit = 1; iHit < silTr->GetNHits(); iHit++) {
1578 BmnHit* hit = (BmnHit*)fSilHitsArray->UncheckedAt(silTr->GetHitIndex(iHit));
1579 Double_t x = hit->GetX();
1580 Double_t y = hit->GetY();
1581 Double_t z = hit->GetZ();
1582 TLine* lineZX = new TLine(z, x, zPrev, xPrev);
1583 lineZX->SetLineColor(kRed);
1584 lineZX->SetLineWidth(2);
1585 TLine* lineZY = new TLine(z, y, zPrev, yPrev);
1586 lineZY->SetLineColor(kRed);
1587 lineZY->SetLineWidth(2);
1588 c->cd(1);
1589 lineZX->Draw("same");
1590 c->cd(2);
1591 lineZY->Draw("same");
1592 zPrev = z;
1593 yPrev = y;
1594 xPrev = x;
1595 }
1596 }
1597
1598 if (glTrack->GetGemTrackIndex() != -1) {
1599 gemTr = (BmnTrack*)fGemTracksArray->UncheckedAt(glTrack->GetGemTrackIndex());
1600
1601 Int_t startIdx = 0;
1602
1603 if (glTrack->GetSilTrackIndex() == -1) {
1604 xPrev = ((BmnHit*)fGemHitsArray->At(gemTr->GetHitIndex(0)))->GetX();
1605 yPrev = ((BmnHit*)fGemHitsArray->At(gemTr->GetHitIndex(0)))->GetY();
1606 zPrev = ((BmnHit*)fGemHitsArray->At(gemTr->GetHitIndex(0)))->GetZ();
1607 startIdx = 1;
1608 }
1609
1610 for (Int_t iHit = startIdx; iHit < gemTr->GetNHits(); iHit++) {
1611 BmnHit* hit = (BmnHit*)fGemHitsArray->UncheckedAt(gemTr->GetHitIndex(iHit));
1612 Double_t x = hit->GetX();
1613 Double_t y = hit->GetY();
1614 Double_t z = hit->GetZ();
1615 TLine* lineZX = new TLine(z, x, zPrev, xPrev);
1616 lineZX->SetLineColor(kRed);
1617 lineZX->SetLineWidth(2);
1618 TLine* lineZY = new TLine(z, y, zPrev, yPrev);
1619 lineZY->SetLineColor(kRed);
1620 lineZY->SetLineWidth(2);
1621 c->cd(1);
1622 lineZX->Draw("same");
1623 c->cd(2);
1624 lineZY->Draw("same");
1625 zPrev = z;
1626 yPrev = y;
1627 xPrev = x;
1628 }
1629 }
1630 }
1631
1632 // Int_t nMCtracks = 0;
1633 // for (Int_t iTr = 0; iTr < fMCTracksArray->GetEntriesFast(); ++iTr) {
1634 // CbmMCTrack* tr = (CbmMCTrack*)fMCTracksArray->At(iTr);
1635 // Int_t nPoints = 0; //tr->GetNPoints(kGEM);
1636 // vector<TLine*> vLineZX;
1637 // vector<TLine*> vLineZY;
1638 // Double_t z0 = tr->GetStartZ();
1639 // Double_t x0 = tr->GetStartX();
1640 // Double_t y0 = tr->GetStartY();
1641 // Double_t x1, y1, z1;
1642 // set<Int_t> vSt;
1643 // for (Int_t i = 0; i < fSilPointsArray->GetEntriesFast(); ++i) {
1644 // FairMCPoint* pnt = (FairMCPoint*)fSilPointsArray->At(i);
1645 // if (pnt->GetTrackID() == iTr) {
1646 // nPoints++;
1647 // z1 = pnt->GetZ();
1648 // x1 = pnt->GetX();
1649 // y1 = pnt->GetY();
1650 // vSt.insert(fSilDetector->GetPointStationOwnership(z1));
1651 // vLineZX.push_back(new TLine(z0, x0, z1, x1));
1652 // vLineZY.push_back(new TLine(z0, y0, z1, y1));
1653 // z0 = z1;
1654 // x0 = x1;
1655 // y0 = y1;
1656 // }
1657 // }
1658
1659 // for (Int_t i = 0; i < fGemPointsArray->GetEntriesFast(); ++i) {
1660 // FairMCPoint* pnt = (FairMCPoint*)fGemPointsArray->At(i);
1661 // if (pnt->GetTrackID() == iTr) {
1662 // nPoints++;
1663 // z1 = pnt->GetZ();
1664 // x1 = pnt->GetX();
1665 // y1 = pnt->GetY();
1666 // vSt.insert(fGemDetector->GetPointStationOwnership(z1) + fNSiliconStations);
1667 // vLineZX.push_back(new TLine(z0, x0, z1, x1));
1668 // vLineZY.push_back(new TLine(z0, y0, z1, y1));
1669 // z0 = z1;
1670 // x0 = x1;
1671 // y0 = y1;
1672 // }
1673 // }
1674
1675 // if (nPoints > 3 && vSt.size() > 3) {
1676 // nMCtracks++;
1677 // c->cd(1);
1678 // for (TLine* it : vLineZX) {
1679 // it->SetLineColor(kBlue);
1680 // it->SetLineStyle(9);
1681 // it->Draw("same");
1682 // }
1683 // c->cd(2);
1684 // for (TLine* it : vLineZY) {
1685 // it->SetLineColor(kBlue);
1686 // it->SetLineStyle(9);
1687 // it->Draw("same");
1688 // }
1689 // }
1690 // }
1691
1692 // printf("NMCTracks = %d GlobTracks = %d\n", nMCtracks, fGlobTracksArray->GetEntriesFast());
1693
1694 // arc->Draw("same");
1695 c->SaveAs("hits.png");
1696 getchar();
1697 delete h_HitsZX;
1698 delete h_HitsZY;
1699 delete c;
1700
1701 return kBMNSUCCESS;
1702}
CalculateTrackParTime
static Double_t CalculateTrackParTime
Definition BmnInnerTrackingRun7.cxx:33
construct_3of4_time
static Double_t construct_3of4_time
Definition BmnInnerTrackingRun7.cxx:31
sortTrackTime
static Double_t sortTrackTime
Definition BmnInnerTrackingRun7.cxx:36
construct_4of4_time
static Double_t construct_4of4_time
Definition BmnInnerTrackingRun7.cxx:30
trackUpdateTime
static Double_t trackUpdateTime
Definition BmnInnerTrackingRun7.cxx:35
construct_2of4_time
static Double_t construct_2of4_time
Definition BmnInnerTrackingRun7.cxx:32
hitMatchTimeLoop
static Double_t hitMatchTimeLoop
Definition BmnInnerTrackingRun7.cxx:29
SortHitsTime
static Double_t SortHitsTime
Definition BmnInnerTrackingRun7.cxx:34
trackSelectionTime
static Double_t trackSelectionTime
Definition BmnInnerTrackingRun7.cxx:37
hitMatchTime
static Double_t hitMatchTime
Definition BmnInnerTrackingRun7.cxx:28
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 BmnInnerTrackingRun7.cxx:27
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
iX
@ iX
Definition BmnAlignDefines.h:85
iY
@ iY
Definition BmnAlignDefines.h:86
kSILICON
@ kSILICON
Definition BmnDetectorList.h:25
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::SetType
void SetType(Int_t type)
Definition BmnHit.h:54
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
BmnInnerTrackingRun7::FindTracks_3of4_OnLastGEMStations
BmnStatus FindTracks_3of4_OnLastGEMStations()
Definition BmnInnerTrackingRun7.cxx:543
BmnInnerTrackingRun7::MatchHit
BmnStatus MatchHit(BmnTrack *cand, vector< BmnHit * > sortedHits, Bool_t downstream)
Definition BmnInnerTrackingRun7.cxx:586
BmnInnerTrackingRun7::FindCandidateByTwoStations
BmnStatus FindCandidateByTwoStations(Short_t st0, Short_t st1, vector< BmnTrack > &candidates, vector< BmnHit * > *sortedHits)
Definition BmnInnerTrackingRun7.cxx:631
BmnInnerTrackingRun7::fGlobTracksArray
TClonesArray * fGlobTracksArray
Definition BmnInnerTrackingRun7.h:90
BmnInnerTrackingRun7::FindCandidateByThreeStations
BmnStatus FindCandidateByThreeStations(Short_t st0, Short_t st1, Short_t st2, vector< BmnTrack > &candidates, vector< BmnHit * > *sortedHits)
Definition BmnInnerTrackingRun7.cxx:440
BmnInnerTrackingRun7::CalculateLength
Double_t CalculateLength(BmnTrack *tr)
Definition BmnInnerTrackingRun7.cxx:1476
BmnInnerTrackingRun7::SortHits
BmnStatus SortHits(vector< BmnHit * > *sortedHits)
Definition BmnInnerTrackingRun7.cxx:1028
BmnInnerTrackingRun7::fSilDetector
BmnSiliconStationSet * fSilDetector
Definition BmnInnerTrackingRun7.h:83
BmnInnerTrackingRun7::CalculateTrackParams
BmnStatus CalculateTrackParams(BmnTrack *tr)
Definition BmnInnerTrackingRun7.cxx:1303
BmnInnerTrackingRun7::Init
virtual InitStatus Init()
Definition BmnInnerTrackingRun7.cxx:84
BmnInnerTrackingRun7::SetHitsUsing
void SetHitsUsing(BmnTrack *tr, Bool_t use)
Definition BmnInnerTrackingRun7.cxx:1176
BmnInnerTrackingRun7::FindTracks_2of2_OnFirstGEMStationsDownstream
BmnStatus FindTracks_2of2_OnFirstGEMStationsDownstream()
Definition BmnInnerTrackingRun7.cxx:681
BmnInnerTrackingRun7::fGemDetector
BmnGemStripStationSet * fGemDetector
Definition BmnInnerTrackingRun7.h:82
BmnInnerTrackingRun7::CheckSharedHits
BmnStatus CheckSharedHits(vector< BmnTrack > &sortedTracks)
Definition BmnInnerTrackingRun7.cxx:1502
BmnInnerTrackingRun7::SortTracks
BmnStatus SortTracks(vector< BmnTrack > &inTracks, vector< BmnTrack > &sortedTracks)
Definition BmnInnerTrackingRun7.cxx:975
BmnInnerTrackingRun7::FindTracks_4of4_OnLastGEMStations
BmnStatus FindTracks_4of4_OnLastGEMStations()
Definition BmnInnerTrackingRun7.cxx:267
BmnInnerTrackingRun7::CalcQp
Double_t CalcQp(BmnTrack *track, Bool_t doSimple)
Definition BmnInnerTrackingRun7.cxx:1361
BmnInnerTrackingRun7::FindTracks_2of2_OnFirstGEMStationsUpstream
BmnStatus FindTracks_2of2_OnFirstGEMStationsUpstream()
Definition BmnInnerTrackingRun7.cxx:891
BmnInnerTrackingRun7::CalcMeanSigma
TVector2 CalcMeanSigma(vector< Double_t >)
Definition BmnInnerTrackingRun7.cxx:1346
BmnInnerTrackingRun7::TrackUpdateByKalman
BmnStatus TrackUpdateByKalman(vector< BmnTrack > &cands)
Definition BmnInnerTrackingRun7.cxx:1070
BmnInnerTrackingRun7::TrackUpdateByLine
BmnStatus TrackUpdateByLine(vector< BmnTrack > &cands)
Definition BmnInnerTrackingRun7.cxx:1050
BmnInnerTrackingRun7::TrackSelection
BmnStatus TrackSelection(vector< BmnTrack > &cands)
Definition BmnInnerTrackingRun7.cxx:1102
BmnInnerTrackingRun7::CalcCovMatrix
BmnStatus CalcCovMatrix(BmnTrack *tr)
Definition BmnInnerTrackingRun7.cxx:1203
BmnInnerTrackingRun7::Exec
virtual void Exec(Option_t *opt)
Definition BmnInnerTrackingRun7.cxx:159
BmnSiliconHit::GetStripTotalSignalInLowerLayer
Double_t GetStripTotalSignalInLowerLayer()
Strip layer characteristics of the hit (total signal, cluster size, strip position as the cluster mea...
Definition BmnSiliconHit.cxx:22
BmnSiliconHit::GetStripTotalSignalInUpperLayer
Double_t GetStripTotalSignalInUpperLayer()
Definition BmnSiliconHit.cxx:28
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
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