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