BmnRoot 26.10.0-210-gc52d34b0
Software framework of the BM@N experiment
Loading...
Searching...
No Matches
BmnStsVectorFinder.cxx
Go to the documentation of this file.
1// -------------------------------------------------------------------------
2// ----- BmnStsVectorFinder source file -----
3// ----- Created 3/12/21 by A.Zinchenko, D.Zinchenko, R.Zinchenko -----
4// -------------------------------------------------------------------------
5#include "BmnStsVectorFinder.h"
6
7#include "CbmKF.h"
8#include "CbmKFTrack.h"
9#include "CbmStsCluster.h"
10#include "CbmStsDigiMatch.h"
11#include "CbmStsDigiScheme.h"
12#include "CbmStsHit.h"
13#include "CbmStsPoint.h"
14#include "CbmStsSector.h"
15#include "CbmStsStation.h"
16#include "CbmStsTrack.h"
17#include "FairRootManager.h"
18#include "FairRunAna.h"
19#include "TMVA/Tools.h" //AZ-300922
20
21#include <TClonesArray.h>
22#include <TGraph.h>
23#include <TMath.h>
24#include <iostream>
25#include <map>
26#include <set>
27#include <vector>
28
29using std::cout;
30using std::endl;
31using std::map;
32using std::multimap;
33using std::pair;
34using std::set;
35using std::vector;
36
37#include <TStopwatch.h>
38
39// std::ofstream outf("log.txt");
40
41static Double_t workTime = 0.0;
42
43// ----- Default constructor -------------------------------------------
44BmnStsVectorFinder::BmnStsVectorFinder()
45 : FairTask("STS Vector Finder")
46 , fHitArray(NULL)
47 , fTrackArray(NULL)
48 , fVectorArray(NULL)
49 , fPass(0)
50 , fNsta(0)
51 , fNbranches(5)
52 , fExact(0)
53 , fExactSel(-9)
54 , discarded(0)
55 , fOut(0)
56 , // 0,2,3,6
57 fNhitsMin(NULL)
58 , fMatBudgetFileName("")
59 , fUseTMVA(0)
60 , // 0,2,3 //AZ-300922
61 fNSi(4)
62// fGemConfigFile("GemRun8.xml"),
63// fSilConfigFile("SiliconRun8.xml")
64{
65 for (int j = 0; j < 4; ++j)
66 fClusArray[j] = nullptr;
67}
68// -------------------------------------------------------------------------
69
70// ----- Destructor ----------------------------------------------------
71BmnStsVectorFinder::~BmnStsVectorFinder()
72{
73 if (fVectorArray) {
74 fVectorArray->Delete();
75 delete fVectorArray;
76 }
77}
78// -------------------------------------------------------------------------
79
80// ----- Public method Init --------------------------------------------
81InitStatus BmnStsVectorFinder::Init()
82{
83
84 // Get RootManager
85 FairRootManager* ioman = FairRootManager::Instance();
86 if (!ioman) {
87 cout << "-E- BmnStsVectorFinder::Init: "
88 << "RootManager not instantiated!" << endl;
89 return kFATAL;
90 }
91
92 // Get input arrays
93 fClusArray[0] = (TClonesArray*)ioman->GetObject("StsCluster");
94 if (!fClusArray[0]) {
95 // cout << "-W- BmnStsVectorFinder::Init: "
96 // << "No StsCluster array!" << endl;
97 // return kERROR;
98 fClusArray[0] = (TClonesArray*)ioman->GetObject("BmnSiliconUpperCluster");
99 fClusArray[1] = (TClonesArray*)ioman->GetObject("BmnSiliconLowerCluster");
100 fClusArray[2] = (TClonesArray*)ioman->GetObject("BmnGemUpperCluster");
101 fClusArray[3] = (TClonesArray*)ioman->GetObject("BmnGemLowerCluster");
102 for (int j = 0; j < 4; ++j) {
103 // AZ-241222 if ( ! fClusArray[j] ) {
104 if (!fClusArray[0] && !fClusArray[2]) { // AZ-241222
105 cout << "-W- BmnStsVectorFinder::Init: "
106 << "No BmnCluster array!" << j << endl;
107 return kERROR;
108 }
109 }
110 }
111
112 fHitArray = (TClonesArray*)ioman->GetObject("StsHit");
113 if (!fHitArray) {
114 // cout << "-W- BmnStsVectorFinder::Init: "
115 // << "No StsHit array!" << endl;
116 return kERROR;
117 }
118
119 fTrackArray = (TClonesArray*)ioman->GetObject("StsTrack");
120 if (!fTrackArray) {
121 // cout << "-W- BmnStsVectorFinder::Init: "
122 // << "No StsTrack array!" << endl;
123 // AZ-100722 return kERROR;
124 }
125
126 fDigiMatches = (TClonesArray*)ioman->GetObject("StsDigiMatch");
127 if (!fDigiMatches) {
128 // cout << "-W- BmnStsVectorFinder::Init: "
129 // << "No StsDigiMatch array!" << endl;
130 // return kERROR;
131 }
132
133 fStsPoints = (TClonesArray*)ioman->GetObject("StsPoint");
134 if (!fStsPoints) {
135 // cout << "-W- BmnStsVectorFinder::Init: "<< "No StsPoint array!" << endl;
136 // AZ-181222 return kERROR;
137 }
138
139 fSilPoints = (TClonesArray*)ioman->GetObject("SiliconPoint");
140 if (!fSilPoints) {
141 // cout << "-E- BmnStsVectorFinder::Init: No SiliconPoint array!" << endl;
142 }
143
144 // Create and register output array
145 fVectorArray = new TClonesArray("CbmStsTrack");
146 ioman->Register("StsVector", "STS", fVectorArray, kTRUE);
147
148 if (fOut / 3 > 0) { // AZ-260922
149 // Store triplets
150 fMap3OutPtr = &fMap3Out;
151 ioman->RegisterAny("Triplets", fMap3OutPtr, kTRUE);
152 }
153 if (fOut && fOut % 2 == 0) { // AZ-260922
154 // Store doublets
155 fMap2OutPtr = &fMap2Out;
156 ioman->RegisterAny("Doublets", fMap2OutPtr, kTRUE);
157 }
158
159 fXyzv[0] = 0.0;
160 // fXyzv[1] = 0.86;
161 // fXyzv[2] = -24.1;
162 fXyzv[1] = 0.0;
163 fXyzv[2] = 0.0;
164 // fXyzv[2] = -3.0; //AZ-271222
165
166 fitter.Init();
167
168 // Define logic
169 static Int_t nHitsMin[20] = {7, 5, 4, 4}; // min number of hits per track vs iteration - Run 7
170 // static Int_t nHitsMin[20] = {7, 5, 4, 3}; // min number of hits per track vs iteration - Run 7
171 // static Double_t dTanX[20] = {0.001, 0.003, 0.002, 0.008}; // window size in TanX vs iteration
172 static Double_t dTanX[20] = {0.002, 0.003, 0.002, 0.008}; // window size in TanX vs iteration
173 static Double_t dTanXTMVA[20] = {0.003, 0.0045, 0.003, 0.012}; // window size in TanX vs iteration - for TMVA
174 static Double_t dTanXB0[20] = {0.001, 0.002, 0.002, 0.002}; // window size in TanX vs iteration - for zero field
175 // static Double_t dTanY[20] = {0.01, 0.02, 0.01, 0.02}; // window size in TanY vs iteration
176 static Double_t dTanY[20] = {0.01, 0.02, 0.02, 0.03}; // window size in TanY vs iteration
177 static Double_t dTanYTMVA[20] = {0.015, 0.03, 0.03, 0.045}; // window size in TanY vs iteration - for TMVA
178 static Double_t dTanYB0[20] = {0.01, 0.02, 0.02, 0.02}; // window size in TanY vs iteration - for zero field
179 static Double_t dTanY3[20] = {0.01, 0.02, 0.02, 0.03}; // window size in TanY vs iteration
180 static Double_t dTanY3TMVA[20] = {0.015, 0.03, 0.03, 0.045}; // window size in TanY vs iteration - for TMVA
181 static Double_t dTanY3B0[20] = {0.01, 0.02, 0.02, 0.02}; // window size in TanY vs iteration - for zero field
182 static Double_t dTanXzero[20] = {0.01, 0.02, 0.02, 0.02}; // window size in TanX vs iteration - for zero field
183 // static Double_t dVarX[20] = {0.004, 0.01, 0.006, 0.02}; // window size in VarX vs iteration - coherent with
184 // ptCut
185 static Double_t ptCut[20] = {1.0, 0.3, 0.7, 0.1}; // min pxz
186 static Double_t ptCutB0[20] = {1., 1., 1., 1.}; // AZ-010123: min pxz - zero field
187 // static Double_t ptCutB0[20] = {0.5, 0.5, 0.5, 0.5 }; //AZ-010123: min pxz - zero field
188 static Double_t curvSta[20] = {
189 12.0, 12.0, 12.0, 7.5, 5.5,
190 4.5, 3.6, 3.0, 3.0, 3.0}; // AZ-310722 - max curvature in stations (for tracks starting in Si)
191 // static Double_t curvSta4[20] = {12.0, 12.0, 12.0, 10.0, 7.5, 5.5, 4.5, 3.6, 3.0, 3.0, 3.0}; //AZ-051022 -
192 // max curvature in stations (for tracks starting in Si) - 4 Si stations
193 static Double_t phiXZ[20] = {0.76, 0.69, 0.69, 0.95, 0.74, 0.60, 0.51, 0.44, 0.39, 0.35}; // Max phi_XZ in
194 // stations
195 // static Double_t phiXZ4[20] = {0.76, 0.69, 0.69, 0.66, 0.95, 0.74, 0.60, 0.51, 0.44, 0.39, 0.35}; //Max
196 // phi_XZ in stations - 4 Si
197 static Double_t phiXZ4[20] = {0.60, 0.63, 0.61, 0.62, 0.95, 0.74,
198 0.60, 0.51, 0.44, 0.39, 0.35}; // AZ-051122 - 4 Si 6_10_14_18
199 static Double_t tanXmax[20] = {0.95, 0.83, 0.83, 1.40, 0.91,
200 0.68, 0.56, 0.47, 0.41, 0.37}; // Max tan_X in stations
201 // static Double_t tanXmax4[20] = {0.95, 0.83, 0.83, 0.77, 1.40, 0.91, 0.68, 0.56, 0.47, 0.41, 0.37}; //Max
202 // tan_X in stations - 4 Si
203 static Double_t tanXmax4[20] = {0.68, 0.73, 0.69, 0.71, 1.40, 0.91,
204 0.68, 0.56, 0.47, 0.41, 0.37}; // AZ-051122 - 4 Si 6_10_14_18
205 static Double_t zMean[20] = {18.5, 27.2, 35.8, 60.9, 91.7,
206 123.4, 154.0, 185.8, 216.5, 248.0}; // Mean Z of stations
207 // static Double_t zMean4[20] = {18.5, 27.2, 35.8, 44.4, 60.9, 91.7, 123.4, 154.0, 185.8, 216.5, 248.0}; //Mean Z
208 // of stations - 4 Si AZ-301222 static Double_t zMean4[20] = {14.9, 23.5, 32.1, 40.8, 60.9, 91.7, 123.4, 154.0,
209 // 185.8, 216.5, 248.0}; //AZ-051122 - 4Si 6_10_14_18
210 static Double_t zMean4[20] = {17.6, 26.4, 35.1, 43.7, 63.6, 94.2,
211 126.5, 156.8, 189.0, 219.2, 251.4}; // AZ-301222 - 4Si 6_10_14_18 - beam data
212 //*/
213
214 FairField* magField = FairRunAna::Instance()->GetField();
215 fBy = TMath::Abs(magField->GetBy(0.0, 0.0, 135.0)) / 10; // max. field in Tesla
216
217 fNhitsMin = nHitsMin;
218 fdTanX = (fUseTMVA) ? dTanXTMVA : dTanX;
219 fdTanY = (fUseTMVA) ? dTanYTMVA : dTanY;
220 fdTanY3 = (fUseTMVA) ? dTanY3TMVA : dTanY3;
221 fPTcut = ptCut;
222 // AZ-311222
223 if (fBy < 0.2) {
224 fdTanX = dTanXB0;
225 fdTanY = dTanYB0;
226 fdTanY3 = dTanY3B0;
227 fPTcut = ptCutB0;
228 } // AZ-311222
229 fdTanXB0 = dTanXzero; // AZ-010123
230 fPhiXZ = (fNSi == 4) ? phiXZ4 : phiXZ;
231 fTanXmax = (fNSi == 4) ? tanXmax4 : tanXmax;
232 fZmean = (fNSi == 4) ? zMean4 : zMean;
233 // fVarx = dVarX;
234 fCurvSta = curvSta;
235
236 for (int j = 0; j < 20; ++j)
237 if (fPTcut[j] > 0.001)
238 fPTcut[j] = 1 / fPTcut[j]; // AZ-170722 - curvature
239
240 // fExact = 1; // for debug // removed DZ 13.10.2021
241 // fExactSel = 13; // for debug
242
243 if (fExact) {
244 // if (1) {
245 // Open windows
246 // Double_t scale = 3;
247 Double_t scale = 1;
248 for (Int_t j = 0; j < 20; ++j) {
249 dTanX[j] *= scale;
250 dTanY[j] *= scale;
251 // dX[j] *= scale;
252 }
253 }
254
255 // AZ-271222 for (int j = 0; j < 20; ++j) fCurvSta[j] *= (0.8 / fBy); //AZ-290822 - max curvature vs field
256 if (fBy > 0.1)
257 for (int j = 0; j < 20; ++j)
258 fCurvSta[j] *= (0.8 / fBy);
259 // AZ-311222 else fBy = 0.1; // zero mag. field
260
261 if (fMatBudgetFileName != "")
262 fitter.ReadMatBudget(fMatBudgetFileName);
263
264 if (fUseTMVA == 2)
265 InitTMVA2(); // AZ-221022 - for doublets
266 if (fUseTMVA)
267 InitTMVA3(); // AZ-300922 - for triplets
268
269 // Create tracking detector setups------------------------------------------------------
270 // TString pathConfig = gSystem->Getenv("VMCWORKDIR");
271 // fGemStationSet = new BmnGemStripStationSet(pathConfig + "/parameters/gem/XMLConfigs/" + fGemConfigFile);
272 // fSilStationSet = new BmnSiliconStationSet(pathConfig + "/parameters/silicon/XMLConfigs/" + fSilConfigFile);
273
274 // cout << "-I- BmnStsVectorFinder: Intialisation successfull" << endl;
275 return kSUCCESS;
276}
277// -------------------------------------------------------------------------
278
279// ----- Public method Exec --------------------------------------------
280void BmnStsVectorFinder::Exec(Option_t* opt)
281{
282
283 TStopwatch sw;
284 sw.Start();
285
286 //*
287 // Reset output array
288 if (!fHitArray)
289 Fatal("Exec", "No StsHit array");
290
291 fVectorArray->Delete();
292 fMap2Out.clear(); // AZ-210922
293 fMap3Out.clear(); // AZ-260922
294
295 for (Int_t j = 0; j < 19; ++j) {
296 fCandSet[j].clear();
297 // fTripleCodes[j].clear();
298 fCandCodes[j].clear();
299 }
300
301 // AZ-100722 Int_t nTracks = fTrackArray->GetEntriesFast();
302 // AZ-100722 if (nTracks == 0) return;
303 Int_t nTracks = 0;
304 if (fTrackArray) {
305 nTracks = fTrackArray->GetEntriesFast();
306 if (nTracks == 0)
307 return;
308 }
309
310 const Int_t nPass = 4; // 4; //5; //8;
311 Int_t nsta = CbmKF::Instance()->GetNStsStations();
312 for (Int_t j = 0; j < nPass; ++j)
313 fNskips[j] = nsta - fNhitsMin[j];
314
315 Int_t nHits0 = fHitArray->GetEntriesFast(), idmaxP = 0;
316
317 // Fill cluster-to-hit maps
318 fHit2id.clear();
319 fClusMaps[0].clear();
320 fClusMaps[1].clear();
321 fmapHits.clear();
322
323 for (Int_t ihit = 0; ihit < nHits0; ++ihit) {
324 CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(ihit);
325 if (hit->GetZ() > 400.0)
326 continue; // AZ-290323 - exclude small GEM downstream
327 Int_t iclusF = hit->GetDigi(0);
328 Int_t iclusB = hit->GetDigi(1);
329 fClusMaps[0].insert(pair<Int_t, Int_t>(iclusF, ihit));
330 fClusMaps[1].insert(pair<Int_t, Int_t>(iclusB, ihit));
331 Int_t ista = hit->GetStationNr() - 0;
332 fNsta = TMath::Max(fNsta, ista);
333 // Code below for debugging
334 if (fVerbose > 0) {
335 set<Int_t> idset = GetHitId(hit, idmaxP);
336 for (set<Int_t>::iterator sit = idset.begin(); sit != idset.end(); ++sit)
337 fHit2id.insert(pair<Int_t, Int_t>(ihit, *sit));
338 }
339 }
340
341 Int_t nHitsOut = 0;
342 discarded = 0;
343 BuildTrackCand();
344
345 for (Int_t ipass = 0; ipass < nPass; ++ipass) {
346 // for (Int_t ipass = nPass-1; ipass < nPass; ++ipass) {
347 fPass = ipass;
348 // Int_t minHits = 5; //4; // minimum number of hits on tracks to accept
349 // Int_t minHits = (ipass == 0) ? 15 : 5; // minimum number of hits on tracks to accept
350 Int_t minHits = (ipass == 0) ? 15 : -fNhitsMin[ipass - 1]; // minimum number of hits on tracks to accept
351 TClonesArray* trArray = (ipass == 0) ? fTrackArray : fVectorArray;
352 nHitsOut += ExcludeHits(minHits, trArray);
353
354 if (fVerbose > 0) {
355 cout << "-I- BmnStsVectorFinder: start - " << nHits0 << ", end - " << nHits0 - nHitsOut << endl;
356 }
357
358 Int_t ntr0 = fVectorArray->GetEntriesFast();
359 // BuildTrackCand();
360 /*
361 //ExtendTracks(fNsta-1);
362 FitTracks();
363 RemoveFakes();
364 RemoveDoubles();
365 cout << "\n ***** Pass " << ipass << ": Number of found tracks = "
366 << fVectorArray->GetEntriesFast() - ntr0 << " " << fVectorArray->GetEntriesFast() << "\n" << endl;
367 //if (fNhitsMin[ipass] <= 4) RemoveFakes();
368 */
369 BuildDoublets();
370 BuildTriplets();
371 fTracks.clear(); //
372 BuildTracks();
373 FitTracks();
374 RemoveDoubles();
375 if (fVerbose > 0) {
376 cout << "\n ***** Pass " << ipass << ": Number of found tracks = " << fVectorArray->GetEntriesFast() - ntr0
377 << " " << fVectorArray->GetEntriesFast() << "\n"
378 << endl;
379 }
380 }
381
382 if (fVerbose > 0)
383 cout << "discarded " << discarded << " track candidates" << endl;
384
385 // Post-processing - try to exclude fake tracks (with too many shared clusters)
386 ExcludeFakes();
387
388 sw.Stop();
389 workTime += sw.RealTime();
390}
391
392// -------------------------------------------------------------------------
393
394Int_t BmnStsVectorFinder::ExcludeHits(Int_t minHits, TClonesArray* trArray)
395{
396 // Exclude hits used for tracking
397
398 // AZ-100722 Int_t nTracks = trArray->GetEntriesFast();
399 Int_t nTracks = 0; // AZ-100722
400 if (trArray)
401 nTracks = trArray->GetEntriesFast(); // AZ-100722
402
403 if (nTracks == 0)
404 return 0;
405 Int_t nHitsOut = 0;
406 // AZ-130722 multimap<Int_t,Int_t>::iterator mit;
407 // AZ-130722 pair<multimap<Int_t,Int_t>::iterator,multimap<Int_t,Int_t>::iterator> ret;
408 unordered_multimap<Int_t, Int_t>::iterator mit;
409 pair<unordered_multimap<Int_t, Int_t>::iterator, unordered_multimap<Int_t, Int_t>::iterator> ret; // AZ-130722
410
411 for (Int_t itra = 0; itra < nTracks; ++itra) {
412 CbmStsTrack* track = (CbmStsTrack*)trArray->UncheckedAt(itra);
413 if (track->GetNStsHits() < minHits) {
414 trArray->Remove(track);
415 continue;
416 }
417
418 Int_t nHitsTr = track->GetNStsHits();
419
420 for (Int_t ihit = 0; ihit < nHitsTr; ++ihit) {
421 int indx = track->GetStsHitIndex(ihit); // AZ-130722
422 // AZ-130722 CbmStsHit *hit = (CbmStsHit*) fHitArray->UncheckedAt(track->GetStsHitIndex(ihit));
423 CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(indx); // AZ-130722
424 if (!hit)
425 continue;
426
427 if (fNhitsMin[fPass] < -5) {
428 // if (fNhitsMin[fPass] > 6) { //AZ-180722
429 // Exclude all hits created from given clusters
430 Int_t iclusF = hit->GetDigi(0);
431 ret = fClusMaps[0].equal_range(iclusF);
432 for (mit = ret.first; mit != ret.second; ++mit) {
433 CbmStsHit* hit1 = (CbmStsHit*)fHitArray->UncheckedAt(mit->second);
434 if (hit1->GetUniqueID() == 0)
435 ++nHitsOut; // excluded hit - used in track
436 hit1->SetUniqueID(1); // exclude hit - used in track
437 // AZ-130722 fmapHits[ihit].used = 1;
438 fmapHits[mit->second].used = 1; // AZ-130722
439 }
440 Int_t iclusB = hit->GetDigi(1);
441 ret = fClusMaps[1].equal_range(iclusB);
442 for (mit = ret.first; mit != ret.second; ++mit) {
443 CbmStsHit* hit1 = (CbmStsHit*)fHitArray->UncheckedAt(mit->second);
444 if (hit1->GetUniqueID() == 0)
445 ++nHitsOut; // excluded hit - used in track
446 hit1->SetUniqueID(1); // exclude hit - used in track
447 // AZ-130722 fmapHits[ihit].used = 1;
448 fmapHits[mit->second].used = 1; // AZ-130722
449 }
450 }
451 if (hit->GetUniqueID() == 0)
452 ++nHitsOut;
453 hit->SetUniqueID(1);
454 // AZ-130722 fmapHits[ihit].used = 1;
455 fmapHits[indx].used = 1; // AZ-130722
456
457 // AZ-130722 Remove used doublets
458 //*
459 int ista = hit->GetStationNr() - 1;
460 ret = fMap2[ista].equal_range(indx);
461 fMap2[ista].erase(ret.first, ret.second);
462 // AZ-140722 Remove used triplets
463 ret = fMap3[ista].equal_range(indx);
464 for (mit = ret.first; mit != ret.second; ++mit) {
465 candvec cand = fCandVec3[ista][mit->second];
466 string code = MakeCode(cand);
467 fMapCode3[ista].erase(code);
468 }
469 fMap3[ista].erase(ret.first, ret.second);
470 //*/
471 }
472 }
473
474 trArray->Compress();
475 return nHitsOut;
476}
477
478// -------------------------------------------------------------------------
479
480set<Int_t> BmnStsVectorFinder::GetHitId(CbmStsHit* hit, Int_t& idmaxP)
481{
482 // Get IDs contributing to given hit (for debugging purposes)
483
484 // AZ-241222 if (fClusArray[1]) return GetHitIdBmn (hit, idmaxP); // BmnRoot data structures
485 if (1)
486 return GetHitIdBmn(hit, idmaxP); // BmnRoot data structures - AZ-241222
487
488 set<Int_t> ids;
489 // return ids; // FIXME
490
491 Int_t iclusF = hit->GetDigi(0);
492 Int_t iclusB = hit->GetDigi(1);
493 CbmStsCluster* clusF = (CbmStsCluster*)fClusArray[0]->UncheckedAt(iclusF);
494 CbmStsCluster* clusB = (CbmStsCluster*)fClusArray[0]->UncheckedAt(iclusB);
495
496 map<Int_t, Double_t> indF, indB;
497 // set<Int_t> ids;
498 Int_t nDigis = clusF->GetNDigis();
499 Double_t pmax = 0.0;
500 TVector3 mom3;
501
502 for (Int_t j = 0; j < nDigis; ++j) {
503 CbmStsDigiMatch* digiMatch = (CbmStsDigiMatch*)fDigiMatches->UncheckedAt(clusF->GetDigi(j));
504 for (Int_t j1 = 0; j1 < 3; ++j1) {
505 Int_t ip = digiMatch->GetRefIndex(j1);
506 if (ip < 0)
507 break;
508 CbmStsPoint* p = (CbmStsPoint*)fStsPoints->UncheckedAt(ip);
509 p->Momentum(mom3);
510 indF[p->GetTrackID()] = mom3.Mag();
511 }
512 }
513
514 nDigis = clusB->GetNDigis();
515
516 for (Int_t j = 0; j < nDigis; ++j) {
517 CbmStsDigiMatch* digiMatch = (CbmStsDigiMatch*)fDigiMatches->UncheckedAt(clusB->GetDigi(j));
518 for (Int_t j1 = 0; j1 < 3; ++j1) {
519 Int_t ip = digiMatch->GetRefIndex(j1);
520 if (ip < 0)
521 break;
522 CbmStsPoint* p = (CbmStsPoint*)fStsPoints->UncheckedAt(ip);
523 // Not needed p->Momentum(mom3);
524 indB[p->GetTrackID()] = 1.0; // mom3.Mag();
525 }
526 }
527
528 idmaxP = -1;
529
530 for (map<Int_t, Double_t>::iterator mit = indF.begin(); mit != indF.end(); ++mit) {
531 Int_t id = (indB.count(mit->first) == 0) ? -mit->first : mit->first;
532 ids.insert(id); // negative ID for ghost crossings
533 if (id >= 0 && mit->second > pmax) {
534 pmax = mit->second;
535 idmaxP = id;
536 }
537 }
538 // This is to be consistent with CBM hit matching procedure
539 // (hit is matched with the closest point)
540 if (hit->GetRefIndex() >= 0) {
541 CbmStsPoint* p = (CbmStsPoint*)fStsPoints->UncheckedAt(hit->GetRefIndex());
542 idmaxP = p->GetTrackID();
543 } else
544 idmaxP = -1;
545
546 for (map<Int_t, Double_t>::iterator mit = indB.begin(); mit != indB.end(); ++mit) {
547 if (ids.count(mit->first) > 0)
548 continue;
549 ids.insert(-mit->first); // negative ID for ghost crossings
550 }
551
552 return ids;
553}
554
555// -------------------------------------------------------------------------
556
557set<Int_t> BmnStsVectorFinder::GetHitIdBmn(CbmStsHit* hit, Int_t& idmaxP)
558{
559 // Get IDs contributing to given hit (for debugging purposes) - for BmnRoot data structures
560
561 set<Int_t> ids;
562 // int ind = 0; // Si hit
563 // if (hit->GetSystemId() == kGEM) ind = 2; // GEM hit
564
565 if (hit->GetRefIndex() == -1 || fStsPoints == nullptr || fSilPoints == nullptr)
566 ids.insert(-1); // fake hit
567 else {
568 // AZ-150722 CbmStsPoint *p = (CbmStsPoint*) fStsPoints->UncheckedAt(hit->GetRefIndex());
569 FairMCPoint* p = nullptr;
570 if (hit->GetSystemId() == kGEM)
571 p = (FairMCPoint*)fStsPoints->UncheckedAt(hit->GetRefIndex());
572 else
573 p = (FairMCPoint*)fSilPoints->UncheckedAt(hit->GetRefIndex());
574 ids.insert(p->GetTrackID());
575 }
576 idmaxP = *ids.begin(); // AZ-290922
577 return ids;
578}
579
580// -------------------------------------------------------------------------
581
582void BmnStsVectorFinder::BuildTrackCand()
583{
584 // Build track candidates
585
586 for (Int_t ist = 0; ist < fNsta; ++ist) {
587 fmapPhx[ist].clear();
588 fmapTy[ist].clear();
589 fmapX[ist].clear();
590 fmapY[ist].clear();
591 fSeedVec[ist].clear();
592 fCandVec[ist].clear();
593 // fCandMap2[ist].clear();
594 // fCandMap3[ist].clear();
595 fCandVec2[ist].clear();
596 fMap2[ist].clear();
597 fCandVec3[ist].clear();
598 fMap3[ist].clear();
599 fMapCode3[ist].clear();
600 }
601 // fmapHits.clear();
602
603 Int_t nHits = fHitArray->GetEntriesFast(), idmaxP = 0;
604 if (fVerbose > 0)
605 cout << "nHits " << nHits << endl;
606 TVector3 pos;
607
608 for (Int_t ih = 0; ih < nHits; ++ih) {
609 CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(ih);
610 if (hit->GetUniqueID())
611 continue; // used hit
612 if (hit->GetZ() > 400.0)
613 continue; // AZ-290323 - exclude small GEM downstream
614 // if (hit->GetSystemId() == kGEM) continue; //!!!AZ-100123 - exclude GEMs
615 // if (hit->GetSystemId() != kGEM) continue; //!!!AZ-100123 - exclude Si
616 // else if (hit->GetSystemId() != kGEM && hit->GetStationNr() == 1) continue; //!!!AZ-110123 - exclude station 1
617 // else if (hit->GetSystemId() != kGEM && hit->GetStationNr() == 4) continue; //!!!AZ-110123 - exclude station 1
618 // if (hit->GetSystemId() == kGEM && hit->GetStationNr() > 7 && TMath::Abs(hit->GetX()) < 10.0 &&
619 // TMath::Abs(hit->GetY()) < 15.0) continue; //!!!AZ-130123 - exclude beam region
620 /*if (hit->GetSystemId() == kGEM && hit->GetSectorNr() > 4) continue; //300123 - exclude lower GEMs
621 else if (hit->GetSystemId() != kGEM) {
622 //AZ-300123 - exclude lower Si
623 if (hit->GetStationNr() == 1 && hit->GetSectorNr() > 3) continue;
624 if (hit->GetStationNr() == 2 && hit->GetSectorNr() > 5) continue;
625 if (hit->GetStationNr() == 3 && hit->GetSectorNr() > 7) continue;
626 if (hit->GetStationNr() == 4 && hit->GetSectorNr() > 9) continue;
627 }*/
628 // AZ-260423 - apply cluster charge correlation cut
629 // if (hit->GetSystemId() == kGEM && TMath::Abs(hit->GetSignalDiv()+0.08) > 0.5) continue;
630 // else if (hit->GetSystemId() != kGEM && TMath::Abs(hit->GetSignalDiv()-0.11) > 0.5) continue;
631 if (fExactSel >= 0) {
632 // Accept only hits with given ID - for debug
633 set<Int_t> ids = GetHitId(hit, idmaxP);
634 if (ids.count(fExactSel) == 0)
635 continue;
636 }
637 if (fExact) {
638 // Accept only true hits
639 set<Int_t> ids = GetHitId(hit, idmaxP);
640 if (idmaxP < 0)
641 continue;
642 }
643 Int_t ista = hit->GetStationNr() - 1;
644 //!!! For debug
645 // if (ista > 2) continue;
646 // if (ista >= 4) continue; //AZ-201222 - Si only
647 //!!!
648 hit->Position(pos);
649 Double_t dx = hit->GetX() - fXyzv[0];
650 Double_t dy = hit->GetY() - fXyzv[1];
651 Double_t dz = hit->GetZ() - fXyzv[2];
652 Double_t by = (fBy < 0.2) ? 0.8 : fBy; // AZ-271222
653 // AZ-250822 if (fPass >= 0) fmapHits[ih] = hitinfo(pos,dx/dz,dy/dz);
654 // 271222 if (fPass >= 0) fmapHits[ih] = hitinfo (pos, TMath::ATan2(dx,dz)/fPhiXZ[ista]/fZmean[ista]/fBy,
655 // dy/dz/fTanXmax[ista]); //AZ-280822 AZ-010123 if (fPass >= 0) fmapHits[ih] = hitinfo (pos,
656 // TMath::ATan2(dx,dz)/fPhiXZ[ista]/fZmean[ista]/by, dy/dz/fTanXmax[ista]); //AZ-271222
657 if (fPass >= 0)
658 fmapHits[ih] = hitinfo(pos, TMath::ATan2(dx, dz) / fPhiXZ[ista] / fZmean[ista] / by,
659 dx / dz / fTanXmax[ista], dy / dz / fTanXmax[ista]); // AZ-010123
660 fmapX[ista].insert(pair<Double_t, Int_t>(pos[0], ih));
661 fmapY[ista].insert(pair<Double_t, Int_t>(pos[1], ih));
662 fmapPhx[ista].insert(pair<Double_t, Int_t>(fmapHits[ih].phx, ih));
663 fmapTy[ista].insert(pair<Double_t, Int_t>(fmapHits[ih].ty, ih));
664 }
665
666 // const Int_t stastop[9] = { fNsta-1, fNsta-2, fNsta-4}; // stop station vs iteration
667 // const Int_t stastop[9] = { fNsta-3, fNsta-2, fNsta-3}; // stop station vs iteration
668 Int_t stastop = fNsta - 1 - fNskips[fPass];
669 stastop = 0;
670
671 if (fVerbose > 0)
672 cout << "fNskips fPass " << fNskips[fPass] << " " << fPass << " stastop " << stastop << " fNsta-1 " << fNsta - 1
673 << endl;
674
675 for (Int_t ista = fNsta - 1; ista >= stastop; --ista) {
676 if (fVerbose > 0)
677 cout << "xmapsize " << fmapX[ista].size() << endl;
678 for (multimap<Double_t, Int_t>::iterator mit = fmapX[ista].begin(); mit != fmapX[ista].end(); ++mit) {
679 candvec aaa;
680 // AZ-280323 aaa.nskips = fNsta - 1 - ista;
681 aaa.nskips = 0; // AZ-280323
682 aaa.momxz = 0.0;
683 aaa.stahit[ista] = mit->second;
684 aaa.code = "-" + to_string(aaa.stahit[ista]) + "-"; // hit index coded
685 // if (fExact) {
686
687 set<Int_t> ids = GetHitId(mit->second, idmaxP);
688 aaa.idmaxP = idmaxP;
689 //}
690 // cout << "aaa.z " <<((CbmStsHit*) fHitArray->UncheckedAt(aaa.stahit[ista]))->GetZ() << " ";
691 fSeedVec[ista].push_back(aaa);
692 } // for (multimap<Double_t,Int_t>::iterator mit
693 if (fVerbose > 0) {
694 Int_t ncand = fSeedVec[ista].size();
695 cout << " Vector stat: " << ista << " " << ncand; // << endl;
696 // AZ-130722 pair<multimap<Int_t, Int_t>::iterator, multimap<Int_t, Int_t>::iterator> ret;
697 pair<unordered_multimap<Int_t, Int_t>::iterator, unordered_multimap<Int_t, Int_t>::iterator>
698 ret; // AZ-130722
699
700 for (Int_t j = 0; j < ncand; ++j) {
701 Int_t ih = fSeedVec[ista][j].stahit[ista];
702 ret = fHit2id.equal_range(ih);
703 cout << " (" << ih << "*";
704 // AZ-130722 for (multimap<Int_t, Int_t>::iterator mit = ret.first; mit != ret.second; ++mit) {
705 for (unordered_multimap<Int_t, Int_t>::iterator mit = ret.first; mit != ret.second; ++mit)
706 { // AZ-130722
707 if (mit != ret.first)
708 cout << ":";
709 cout << mit->second;
710 }
711 cout << ")";
712 }
713 cout << "\n";
714 }
715
716 // Extend track candidates
717 // if (fCandVec[ista].size()) ExtendTracks(ista);
718 } // for (Int_t ista = fNsta-1;
719}
720
721// -------------------------------------------------------------------------
722
723void BmnStsVectorFinder::BuildDoublets()
724{
725 // Build doublets
726
727 int rebuild2 = 1; // AZ-140722
728
729 for (Int_t ist = 0; ist < fNsta; ++ist) {
730 fCandVec[ist].clear();
731 fCandVec2[ist].clear();
732 fMap2[ist].clear();
733 fCandVec3[ist].clear();
734 fMap3[ist].clear();
735 fMapCode3[ist].clear();
736 // AZ-140722 - Does not seem to speed up
737 /*
738 if ( (fPass == 0) || (TMath::Abs(fdTanX[fPass]-fdTanX[fPass-1]) > 0.001) || //AZ-140722
739 (TMath::Abs(fdTanY[fPass]-fdTanY[fPass-1]) > 0.001) ) {
740 fCandVec2[ist].clear();
741 fMap2[ist].clear();
742 rebuild2 = 1;
743 }
744 */
745 }
746 if (rebuild2 == 0)
747 return; // AZ-140722 - do not rebuild doublets
748
749 Int_t idmaxP = 0;
750
751 for (Int_t ista = 0; ista < fNsta - 1; ++ista) {
752 Int_t nTra = fSeedVec[ista].size();
753 // AZ-280323 Int_t istanext = ista + 1;
754 if (nTra == 0)
755 continue;
756
757 for (int istanext = ista + 1; istanext < ista + 3; ++istanext) {
758 // for (int istanext = ista + 1; istanext < ista + 2; ++istanext) { // no jumps over station
759 if (istanext > fNsta - 1)
760 continue; // AZ-280323
761 if (fSeedVec[istanext].size() == 0)
762 continue;
763 Double_t dty = fdTanY[fPass], dphx = fdTanX[fPass];
764 // dty *= 5; //AZ-040822
765 // dtx *= 2; //AZ-040822
766 int nskips = istanext - (ista + 1);
767
768 for (Int_t itra = 0; itra < nTra; ++itra) {
769 candvec& aaa = fSeedVec[ista][itra];
770 Int_t ih = aaa.stahit[ista];
771 if (fmapHits[ih].used)
772 continue; // used hit
773 // if (fMap2[ista].count
774 Double_t phx = fmapHits[ih].phx;
775 phx *= (fPhiXZ[ista] * fZmean[ista]); // AZ-280822
776 phx /= (fPhiXZ[istanext] * fZmean[istanext]); // AZ-280822
777 Double_t ty = fmapHits[ih].ty;
778 ty *= (fTanXmax[ista] / fTanXmax[istanext]); // AZ-280822
779 // AZ-250822 Double_t distxz = TMath::Sqrt (fmapHits[ih].tx*fmapHits[ih].tx + 1) *
780 // (fmapHits[ih].xyz.Z() - fXyzv[2]); //AZ-310722 - XZ-distance to PV
781 Double_t tx = fmapHits[ih].tx;
782 tx *= (fTanXmax[ista] / fTanXmax[istanext]); // AZ-010123
783 multimap<Double_t, candvec> mapDoublets; // AZ-251022 TMVA output -> doublets
784
785 // Get hits on the downstream station
786 multimap<Double_t, Int_t>::iterator mityb = fmapTy[istanext].lower_bound(ty - dty);
787 multimap<Double_t, Int_t>::iterator mitye = fmapTy[istanext].upper_bound(ty + dty);
788 multimap<Double_t, Int_t>::iterator mitxb = fmapPhx[istanext].lower_bound(phx - dphx);
789 multimap<Double_t, Int_t>::iterator mitxe = fmapPhx[istanext].upper_bound(phx + dphx);
790 // Get hits from the acceptance window
791 set<Int_t> setTx, setTy, intersect;
792 for (multimap<Double_t, Int_t>::iterator mit = mitxb; mit != mitxe; ++mit)
793 if (fmapHits[mit->second].used == 0)
794 setTx.insert(mit->second);
795 for (multimap<Double_t, Int_t>::iterator mit = mityb; mit != mitye; ++mit)
796 if (fmapHits[mit->second].used == 0)
797 setTy.insert(mit->second);
798 set_intersection(setTx.begin(), setTx.end(), setTy.begin(), setTy.end(),
799 std::inserter(intersect, intersect.begin()));
800 // cout << " Intersect: " << ista << " " << setTx.size() << " " << setTy.size() << " " <<
801 // intersect.size() << endl;
802 Double_t xp = fmapHits[ih].xyz[0]; // AZ-280722
803 Double_t zp = fmapHits[ih].xyz[2]; // AZ-280722
804 // int n0 = fMap2[ista].size(); //AZ-241022 - for debug
805
806 for (set<Int_t>::iterator sit = intersect.begin(); sit != intersect.end(); ++sit) {
807 CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(*sit);
808 if (hit->GetUniqueID())
809 continue; // used hit
810 if (fmapHits[*sit].used)
811 continue;
812 if (fBy < 0.2 && TMath::Abs(tx - fmapHits[*sit].tx) > fdTanXB0[fPass])
813 continue; // AZ-010123 - zero mag. field
814 if (fExact) {
815 // Exact track ID match
816 set<Int_t> ids = GetHitId(*sit, idmaxP);
817 // if (ids.count(aaa.idmaxP) == 0) continue;
818 if (idmaxP != aaa.idmaxP)
819 continue;
820 }
821 // map<Int_t,Int_t> aaa1 = aaa;
822 candvec aaa1 = aaa;
823 aaa1.stahit[istanext] = *sit; // second hit of the doublet
824 // aaa1.nextind = *sit; // second hit of the doublet
825 Double_t dx = xp - fmapHits[*sit].xyz[0]; // AZ-280722
826 Double_t dz = zp - fmapHits[*sit].xyz[2]; // AZ-280722
827 aaa1.lengxz = TMath::Sqrt(dx * dx + dz * dz); // AZ-280722 - doublet length in XZ-plane
828 // AZ-250822 Double_t varx = 0.0;
829 // AZ-250822 if (!CheckVarx(dx, dz, tx, distxz, varx)) continue; //AZ-310722 - apply cut on
830 // X-variable
831 aaa1.momxz =
832 Curv3(aaa1, aaa1, aaa1, 0); // AZ-310722 - pxz-estimate with the origin at the primary vertex
833 if (fBy > 0.2
834 && (TMath::Abs(aaa1.momxz) > fPTcut[fPass] || TMath::Abs(aaa1.momxz) > fCurvSta[ista]))
835 continue; // AZ-010123
836 else if (fBy < 0.2 && TMath::Abs(aaa1.momxz) / fZmean[1] * fZmean[istanext] > fPTcut[fPass])
837 continue; // AZ-010123 - zero mag. field
838 // AZ-300922 aaa1.ty = fmapHits[*sit].ty; //AZ-290822
839 aaa1.ty = ty - fmapHits[*sit].ty; // AZ-300922
840 aaa1.tx = tx - fmapHits[*sit].tx; // AZ-010123 - for xero mag. field
841 // AZ-250822 aaa1.varx = varx;
842 // cout << " Doublet: " << aaa1.idmaxP << " " << aaa1.stahit.begin()->first << endl; // debug output
843
844 if (fOut && fOut % 2 == 0) { // AZ-210922 for writing out
845 string code = to_string(ih) + "-" + to_string(*sit);
846 if (fMap2Out.find(code) == fMap2Out.end())
847 fMap2Out[code] = fPass;
848 }
849
850 // AZ-221022
851 Double_t tmvaOut = 9.9;
852 if (fUseTMVA == 2)
853 tmvaOut = TMVAOutput2(aaa1);
854 // if (fUseTMVA) outf << tmvaOut << " " << fPass << endl;
855 // if (tmvaOut < 0.02) continue;
856 if (tmvaOut < 0.01)
857 continue;
858 aaa1.nskips = nskips;
859
860 mapDoublets.insert(pair<Double_t, candvec>(-tmvaOut, aaa1)); // AZ-251022
861 } // for (set<Int_t>::iterator sit = intersect.begin();
862
863 // AZ-251022
864 int n2tot = 0, n2max = (fUseTMVA == 2) ? fNbranches * 2 : 999;
865
866 for (auto ait = mapDoublets.begin(); ait != mapDoublets.end(); ++ait) {
867 fCandVec2[ista].push_back(ait->second);
868 fMap2[ista].insert(pair<int, int>(ih, fCandVec2[ista].size() - 1));
869 ++n2tot;
870 if (n2tot >= n2max)
871 break;
872 }
873 // cout << " *** Doublets : " << ista << " " << ih << " " << fMap2[ista].size()-n0 << endl;
874 // if (extendOK) break; // track has been extended - do not jump over station
875 } // for (Int_t itra = 0;
876 } // for (int istanext = ista + 1;
877
878 Int_t ncand = fCandVec2[ista].size();
879 if (fVerbose >= 1)
880 cout << " Doublet stat: " << ista << " " << ncand << endl;
881 } // for (Int_t ista = fNsta-1;
882}
883
884// -------------------------------------------------------------------------
885
886Bool_t BmnStsVectorFinder::CheckVarx(Double_t dx, Double_t dz, Double_t tx, Double_t distxz, Double_t& varx)
887{
888 // AZ-310722 - Check Varx
889
890 varx = (dx / dz - tx) / distxz;
891 // AZ-250822 if (TMath::Abs(varx) > fVarx[fPass]) return kFALSE;
892 return kTRUE;
893}
894
895// -------------------------------------------------------------------------
896
897void BmnStsVectorFinder::BuildTriplets()
898{
899 // Build triplets
900
901 // Int_t idmaxP = 0;
902
903 // AZ-130722 pair<multimap<Int_t,int>::iterator,multimap<Int_t,int>::iterator> ret;
904 // AZ-130722 multimap<Int_t,int>::iterator mit, mit1;
905 pair<unordered_multimap<Int_t, int>::iterator, unordered_multimap<Int_t, int>::iterator> ret; // AZ-130722
906 unordered_multimap<Int_t, int>::iterator mit, mit1; // AZ-130722
907 Double_t pxzCut = (fUseTMVA) ? 0.3 : 0.2; // AZ-041022 - pxz difference cut
908 // if (fBy < 0.2) pxzCut = 1.0; //AZ-261222 - zero magnetic field
909
910 for (Int_t ista = 0; ista < fNsta - 2; ++ista) {
911 discarded = 0;
912 Int_t nTra = fCandVec2[ista].size();
913 // AZ-280323 Int_t istanext = ista + 1;
914 // AZ-280323 if (nTra == 0 || fCandVec2[istanext].size() == 0) continue;
915 if (nTra == 0)
916 continue;
917
918 for (mit = fMap2[ista].begin(); mit != fMap2[ista].end(); ++mit) {
919 candvec& aaa = fCandVec2[ista][mit->second]; // first doublet
920 Int_t istanext = ista + 1 + aaa.nskips; // AZ-280323
921 if (fCandVec2[istanext].size() == 0)
922 continue; // AZ-280323
923 // AZ-300922 Double_t tany = fmapHits[aaa.stahit.begin()->second].ty; //AZ-2908822
924 // AZ-300922 tany *= (fTanXmax[ista] / fTanXmax[istanext+1]); //AZ-290822
925 Double_t tany = aaa.ty; // AZ-300922
926 Double_t tanx = aaa.tx; // AZ-010123
927
928 ret = fMap2[istanext].equal_range(aaa.stahit.rbegin()->second);
929 multimap<Double_t, candvec> mapDoublets; // AZ-110722 second doublets, attached to one first doublet
930 int nbr = 0, newtr = 0, newtr3 = 0;
931 Double_t c2max = 999.0;
932
933 for (mit1 = ret.first; mit1 != ret.second; ++mit1) {
934 candvec aaa1;
935 // AZ-150722 Int_t newtr = (mit1 == ret.first) ? 1 : 0;
936
937 // Int_t nhits = 3;
938 aaa1 = aaa;
939 candvec& aaa2 = fCandVec2[istanext][mit1->second]; // second doublet
940 if (fExact) {
941 // Exact track ID match
942 if (aaa.idmaxP != aaa2.idmaxP)
943 continue;
944 }
945 // AZ-310722 - Check doublet matching criteria
946 // if (TMath::Abs(aaa.varx - aaa2.varx) > 0.002) continue; //AZ-050822
947 // AZ-070822 if (TMath::Abs(aaa.ty - aaa2.ty) > 0.01) continue; //AZ-040822
948 // if (TMath::Abs(aaa.ty - aaa2.ty) > fdTanY3[fPass]) continue; //AZ-250822
949 // cout << TMath::Abs(tany - aaa2.ty) << " " << fdTanY3[fPass] << " "
950 // << TMath::Abs(tanx - aaa2.tx) << " " << fdTanXB0[fPass] << endl;
951 if (TMath::Abs(tany - aaa2.ty) > fdTanY3[fPass])
952 continue; // AZ-290822
953 if (fBy < 0.2 && TMath::Abs(tanx - aaa2.tx) > fdTanXB0[fPass])
954 continue; // AZ-010123 - for zero mag. field
955
956 if (fOut / 3 > 0) { // AZ-260922 for writing out
957 string code = to_string(aaa.stahit.begin()->second) + "-" + to_string(aaa.stahit.rbegin()->second)
958 + "-" + to_string(aaa2.stahit.rbegin()->second);
959 if (fMap3Out.find(code) == fMap3Out.end())
960 fMap3Out[code] = fPass;
961 }
962
963 if (fBy > 0.2
964 && TMath::Abs(aaa.momxz - aaa2.momxz)
965 / (TMath::Abs(aaa.momxz) + TMath::Abs(aaa2.momxz))
966 //// * TMath::Sqrt(aaa.lengxz+aaa2.lengxz) / 10 > 0.2) continue; //AZ-250822
967 * TMath::Sqrt(aaa.lengxz + aaa2.lengxz) / 10
968 > pxzCut)
969 continue; // AZ-041022
970 // if (fBy < 0.2 && TMath::Abs(aaa.momxz-aaa2.momxz) / TMath::Sqrt(fZmean[ista]) > pxzCut) continue;
971 // //AZ-311222
972
973 aaa1.stahit[aaa2.stahit.rbegin()->first] = aaa2.stahit.rbegin()->second; // third hit of the triplet
974
975 Double_t tmvaOut = 9.9; // AZ-300922
976 if (fUseTMVA)
977 tmvaOut = TMVAOutput3(aaa1, aaa2); // AZ-300922
978 // if (fUseTMVA) outf << tmvaOut << " " << fPass << endl; //AZ-300922
979 if (tmvaOut < 0.02)
980 continue; // AZ-300922 - use TMVA MLPBNN method
981
982 map<Int_t, Int_t>& hitMap = aaa1.stahit;
983
984 CbmStsHit* hit = nullptr;
985 CbmStsTrack track;
986 TArrayI& hits = *track.GetStsHits();
987 hits.Set(3);
988 Int_t indx = 0 /*, iok = 1*/;
989
990 /// hitMap ---- aaa1.stahit
991 for (map<Int_t, Int_t>::iterator mit2 = hitMap.begin(); mit2 != hitMap.end(); ++mit2) {
992 hits[indx++] = mit2->second;
993 hit = (CbmStsHit*)fHitArray->UncheckedAt(mit2->second);
994 // cout << " hit: " << mit2->first << " " << mit2->second << endl;
995 // hit->SetDx(0.08/TMath::Sqrt(12.0));
996 hit->SetDx(0.08 / TMath::Sqrt(12.0) * 1.2); // AZ-150823 - old
997 /// hit->SetDx(0.05); //AZ-110123 - early estimate
998 // hit->SetDx(0.015);
999 hit->SetDy(0.1234); // AZ-150823 - old
1000 // AZ-030123 if (hit->GetStationNr() <= 3) hit->SetDx(0.02/TMath::Sqrt(12.0));
1001 if (hit->GetStationNr() <= fNSi)
1002 hit->SetDx(0.02 / TMath::Sqrt(12.0)); // AZ-150823 - old
1003 // if (hit->GetStationNr() <= fNSi) hit->SetDx(0.01); //AZ-040123 - early estimate
1004 /// if (hit->GetStationNr() <= fNSi) hit->SetDx(0.02); //AZ-040123 - early estimate
1005 // if (hit->GetStationNr() % 2 != 0) hit->SetCovXY(1.968e-3); // this is for 2-D hits
1006 // else hit->SetCovXY(-1.968e-3);
1007 // hit->SetDy(0.03); // for test
1008 if (indx == 3)
1009 break;
1010 }
1011
1012 // fitter.DoFit(&track);
1013 // float chi2 = track.GetChi2();
1014 // AZ-310722 Double_t ty = 0.0;
1015 Double_t ty = aaa.ty; // AZ-310722
1016 // AZ-280722 float chi2 = LinearFit (&track, newtr, ty) / hit->GetDy() / hit->GetDy();
1017 float chi2 = LinearFit(aaa, aaa2, &track, newtr, ty) / hit->GetDy() / hit->GetDy(); // AZ-280722
1018 ++newtr; // AZ-150722
1019 // float chi2 = LinearFit (&track, newtr, ty) / 0.09 / 0.09;
1020 // cout << "3hit track Chi2 " << track.GetChi2() << " " << nhits << endl;
1021 // if Chi2 > 10.0 exclude track cand
1022 // TODO DZ add index so I know how many candidates were discarded
1023 // if (aaa1.idmaxP == 78 || aaa1.idmaxP == 222 || aaa1.idmaxP == 127)
1024 // cout << " Triplet: " << aaa1.idmaxP << " " << aaa1.stahit.begin()->first << " " << chi2 << endl; //
1025 // debug output
1026
1027 if (chi2 > 10.0) {
1028 // if (chi2 > 20.0) { //AZ-010123
1029 discarded++;
1030 continue;
1031 }
1032 // Triplet curvature
1033 aaa1.momxz = Curv3(aaa, aaa2, aaa1, newtr3);
1034 newtr3 = 1;
1035 // AZ-150722 if (TMath::Abs(aaa1.momxz) < 0.2) {
1036 // Double_t ptcut = (fPass < 2) ? 0.3 : 0.2; //AZ-150722
1037 // if (TMath::Abs(aaa1.momxz) < ptcut) { //AZ-150722
1038 if (TMath::Abs(aaa1.momxz) > fPTcut[fPass]) { // AZ-170722 - curvature
1039 // if (TMath::Abs(aaa1.momxz) > fPTcut[fPass] || TMath::Abs(aaa1.momxz) > fCurvSta[ista]) {
1040 // //AZ-310722 - curvature
1041 discarded++;
1042 continue;
1043 }
1044
1045 aaa1.ty = ty;
1046 /*AZ-110722
1047 fCandVec3[ista].push_back(aaa1);
1048 fMap3[ista].insert(pair<Int_t,int>(mit->first,fCandVec3[ista].size()-1));
1049 string code = "-";
1050 for (map<Int_t,Int_t>::iterator mitr = hitMap.begin(); mitr != hitMap.end(); ++mitr)
1051 code += (to_string(mitr->second) + "-");
1052 fMapCode3[ista][code] = fCandVec3[ista].size() - 1;
1053 */
1054 if (fUseTMVA) {
1055 // AZ-300922 Keep fNbranches with the highest NN response
1056 if (-tmvaOut > c2max && nbr >= fNbranches)
1057 continue;
1058 } else if (chi2 > c2max && nbr >= fNbranches)
1059 continue;
1060 aaa1.chi2 = chi2; // AZ-140722
1061 Double_t qual = (fUseTMVA) ? -tmvaOut : chi2; // AZ-300922
1062 aaa1.nskips = aaa.nskips + aaa2.nskips; // AZ-280323
1063 // AZ-300922 mapDoublets.insert(pair<Double_t,candvec>(chi2,aaa1)); //AZ-110722
1064 mapDoublets.insert(pair<Double_t, candvec>(qual, aaa1)); // AZ-300922
1065 // mapDoublets.insert(pair<Double_t,candvec>(TMath::Abs(aaa1.momxz),aaa1)); //AZ-280722
1066 nbr = mapDoublets.size();
1067 // AZ-310722 if (nbr > fNbranches) mapDoublets.erase(mapDoublets.rbegin()->first);
1068 if (nbr > fNbranches) {
1069 multimap<Double_t, candvec>::iterator mit0 = mapDoublets.end();
1070 --mit0;
1071 mapDoublets.erase(mit0);
1072 }
1073 c2max = mapDoublets.rbegin()->first;
1074 } // for (mit1 = ret.first;
1075
1076 // Take only first fNbranches with lowest chi2
1077 int nok = 0;
1078 for (multimap<Double_t, candvec>::iterator it = mapDoublets.begin(); it != mapDoublets.end(); ++it) {
1079 if (nok >= fNbranches)
1080 break;
1081 ++nok;
1082 fCandVec3[ista].push_back(it->second);
1083 fMap3[ista].insert(pair<Int_t, int>(mit->first, fCandVec3[ista].size() - 1));
1084 map<Int_t, Int_t>& hitMap = it->second.stahit;
1085 string code = "-";
1086 for (map<Int_t, Int_t>::iterator mitr = hitMap.begin(); mitr != hitMap.end(); ++mitr)
1087 code += (to_string(mitr->second) + "-");
1088 fMapCode3[ista][code] = fCandVec3[ista].size() - 1;
1089 // cout << ista << " " << code << endl;
1090 }
1091 } // for (mit = fCandMap2[ista].begin();
1092
1093 Int_t ncand = fCandVec3[ista].size();
1094 if (fVerbose > 0)
1095 cout << " Triplet stat: " << ista << " " << ncand << " " << discarded << endl;
1096 } // for (Int_t ista = fNsta-1;
1097}
1098
1099// -------------------------------------------------------------------------
1100
1101void BmnStsVectorFinder::BuildTracks()
1102{
1103 // Build tracks starting from triplets
1104
1105 int istaEnd = fNsta - fNhitsMin[fPass] + 1;
1106 istaEnd = TMath::Min(istaEnd, fNsta - 3); // do not allow simple triplets on last 3 stations
1107 // std::multimap<Int_t,int> &candMap = fMap3[ista];
1108 // vector<candvec> &candVec = fCandVec3[ista];
1109 // std::multimap<Int_t,int> *candMap = &fMap2[ista];
1110 // vector<candvec> *candVec = &fCandVec2[ista];
1111
1112 for (Int_t ista = 0; ista < istaEnd; ++ista) {
1113 // std::multimap<Int_t,int> *candMap = &fMap2[ista];
1114 // vector<candvec> *candVec = &fCandVec2[ista];
1115 // AZ-130722 std::multimap<Int_t,int> *candMap = &fMap3[ista];
1116 std::unordered_multimap<Int_t, int>* candMap = &fMap3[ista]; // AZ-130722
1117 vector<candvec>* candVec = &fCandVec3[ista];
1118 Int_t nTra = candMap->size();
1119 // cout << " BuildTracks: " << ista << " " << nTra << " " << istaEnd << " " << fNsta << endl;
1120 if (nTra == 0)
1121 continue;
1122 // Int_t istaup = ista - 2;
1123
1124 // AZ-130722 for (multimap<Int_t,int>::iterator mit = candMap->begin(); mit != candMap->end(); ++mit) {
1125 for (unordered_multimap<Int_t, int>::iterator mit = candMap->begin(); mit != candMap->end(); ++mit) {
1126 candvec& cand = (*candVec)[mit->second];
1127 if (fPass == 0 && cand.nskips > 0)
1128 continue; //!!! AZ-280323 - no jumps over stations at 1'st pass !!!
1129 map<int, int>::iterator it = cand.stahit.begin();
1130 int ista0 = it->first;
1131 ++it;
1132 int ista1 = it->first;
1133 if (ista1 - ista0 > 1)
1134 continue; // AZ-250423 - isolated first hit is not allowed
1135 ExtendTrack(cand);
1136 }
1137 }
1138}
1139
1140// -------------------------------------------------------------------------
1141
1142void BmnStsVectorFinder::ExtendTrack(candvec cand)
1143{
1144 // Extend the track by triplets or doublets
1145
1146 const Double_t c2ndfMax = 30.0; // cut on chi2/NDF
1147 Double_t pxzCut = 0.3; // AZ-261222 - pxz difference cut
1148 // if (fBy < 0.2) pxzCut = 1.0; //AZ-261222 - zero magnetic field
1149
1150 Int_t nTra = 0, ista = cand.stahit.rbegin()->first;
1151 // cout << ista << " " << cand.stahit.size() << endl;
1152 /*
1153 //!!! Check - step back
1154 if (ista == fNsta-2) {
1155 cand.stahit.erase(cand.stahit.rbegin()->first);
1156 ista = cand.stahit.rbegin()->first;
1157 }
1158 */
1159 // std::multimap<Int_t,int> &candMap = (ista == fNsta-2) ? fMap2[ista] : fMap3[ista]; // doublets or triplets
1160 // vector<candvec> &candVec = (ista == fNsta-2) ? fCandVec2[ista] : fCandVec3[ista]; // doublets or triplets
1161 // std::multimap<Int_t,int> *pCandMap = (ista == fNsta-2) ? &fMap2[ista] : &fMap2[ista]; // doublets or triplets
1162 // vector<candvec> *pCandVec = (ista == fNsta-2) ? &fCandVec2[ista] : &fCandVec2[ista]; // doublets or triplets
1163 // AZ-130722 std::multimap<Int_t,int> *pCandMap = (ista == fNsta-2) ? &fMap2[ista] : &fMap2[ista]; // doublets or
1164 // triplets
1165 std::unordered_multimap<Int_t, int>* pCandMap =
1166 (ista == fNsta - 2) ? &fMap2[ista] : &fMap2[ista]; // AZ-130722 doublets or triplets
1167 vector<candvec>* pCandVec = (ista == fNsta - 2) ? &fCandVec2[ista] : &fCandVec2[ista]; // doublets or triplets
1168 // AZ-130722 multimap<Int_t,int> &candMap = *pCandMap;
1169 unordered_multimap<Int_t, int>& candMap = *pCandMap; // AZ-130722
1170 vector<candvec>& candVec = *pCandVec;
1171
1172 if (ista < fNsta - 1) {
1173 /*
1174 if (ista == fNsta-2) {
1175 // Doublets
1176 candMap = fMap2[ista];
1177 candVec = fCandVec2[ista];
1178 }
1179 */
1180 nTra = candMap.count(cand.stahit.rbegin()->second);
1181 }
1182
1183 // cout << " ExtendTrack " << ista << " " << nTra << endl;
1184 if (nTra == 0) {
1185 // No extension found
1186 if ((Int_t)cand.stahit.size() < fNhitsMin[fPass])
1187 return; // too short track
1188 // cout << " ExtendTrack - store " << cand.stahit.size() << " " << cand.idmaxP << endl;
1189 // FitTrack(cand);
1190 string code("-");
1191 for (map<int, int>::iterator it = cand.stahit.begin(); it != cand.stahit.end(); ++it)
1192 code += (to_string(it->second) + "-");
1193 // cout << cand.stahit.rbegin()->first << " " << ista << endl;
1194 if (fCandCodes[cand.stahit.size()].find(code) == fCandCodes[cand.stahit.size()].end()) {
1195 fCandCodes[cand.stahit.size()].insert(code);
1196 // cout << " 1 " << endl;
1197 // AZ-170122 if (cand.stahit.size() >= 4 && FitTrack(cand) > c2ndfMax) return;
1198 // AZ-110123 if (cand.stahit.size() >= 4 && cand.stahit.rbegin()->first != ista && FitTrack(cand) >
1199 // c2ndfMax) return;
1200 if ((Int_t)cand.stahit.size() >= fNhitsMin[fPass] && cand.stahit.rbegin()->first != ista
1201 && FitTrack(cand) > c2ndfMax)
1202 return; // AZ-110123
1203 // cout << " ExtendTrack - store " << cand.stahit.size() << " " << cand.idmaxP << endl;
1204 fCandVec[cand.stahit.begin()->first].push_back(cand);
1205 }
1206 return;
1207 }
1208
1209 // AZ-130722 multimap<Int_t,int>::iterator mit;
1210 // AZ-130722 pair<multimap<Int_t,int>::iterator,multimap<Int_t,int>::iterator> ret;
1211 unordered_multimap<Int_t, int>::iterator mit; // AZ-130722
1212 pair<unordered_multimap<Int_t, int>::iterator, unordered_multimap<Int_t, int>::iterator> ret; // AZ-130722
1213 ret = candMap.equal_range(cand.stahit.rbegin()->second);
1214 string code = "-";
1215 int istaMid = 0;
1216
1217 // Use Kalman track constraints around extrapolated position
1218 /*
1219 set<int> kalmWind;
1220 if (nTra > 5)
1221 if (cand.stahit.size() >= 4) kalmWind = KalmanWindow(cand, candVec[ret.first->second].stahit.begin()->second);
1222 */
1223
1224 for (mit = ret.first; mit != ret.second; ++mit) {
1225 // candvec &aaa = candVec[mit->second];
1226 candvec aaa = candVec[mit->second];
1227 if (fExact && aaa.idmaxP != cand.idmaxP)
1228 continue;
1229 if (fPass == 0 && aaa.nskips > 0)
1230 continue; //!!!AZ-289323 - no jumps over station at 1'st pass
1231
1232 // if (nTra > 5 && cand.stahit.size() >= 4 &&
1233 // kalmWind.find(aaa.stahit.begin()->second) == kalmWind.end()) continue; // outside Kalman window
1234
1235 // Check if there is a triplet started from the last but one point of the current one (cand)
1236 map<Int_t, Int_t>::iterator mitr = aaa.stahit.begin();
1237 ++mitr;
1238 if (code.size() == 1) {
1239 map<Int_t, Int_t>::reverse_iterator mit1 = cand.stahit.rbegin();
1240 ++mit1;
1241 istaMid = mit1->first;
1242 code += (to_string(mit1->second) + "-");
1243 --mit1;
1244 code += (to_string(mit1->second) + "-");
1245 code += (to_string(mitr->second) + "-");
1246 } else {
1247 // Replace last index
1248 int pos = code.rfind("-", code.size() - 2);
1249 code.erase(code.begin() + pos + 1, code.end());
1250 code += (to_string(mitr->second) + "-");
1251 }
1252 // cout << " Check code: " << code << " " << istaMid << endl;
1253 // Check extension possibility
1254 int ok = 1;
1255 if (fMapCode3[istaMid].find(code) == fMapCode3[istaMid].end())
1256 ok = 0; // no triplet found
1257 //!!! Check triplet slopes and curvatures
1258 //*
1259 candvec* bbb = nullptr;
1260 if (ok) {
1261 bbb = &fCandVec3[istaMid][fMapCode3[istaMid][code]];
1262 // Double_t dty = TMath::Abs (cand.ty - bbb.ty);
1263 // if (dty * 2 / (TMath::Abs(cand.ty) + TMath::Abs(bbb.ty)) > 0.1) continue; // different slopes
1264 // if (dty > 0.005) continue; // different slopes
1265 // Curvatures
1266 // cout << cand.momxz << " " << bbb.momxz << endl;
1267 if ((fPass == 0 && bbb->nskips > 0) || (cand.nskips + bbb->nskips > 2))
1268 ok = 0; //!!!AZ-280323 - jumps over stations
1269 if (ok && cand.stahit.size() > 2) {
1270 Double_t dp = TMath::Abs(cand.momxz - bbb->momxz);
1271 // if (dp * 2 / (TMath::Abs (cand.momxz) + TMath::Abs(bbb->momxz)) > 0.3) {
1272 // AZ-060822 if (dp * 2 / (TMath::Abs (cand.momxz) + TMath::Abs(bbb->momxz)) > 0.2) {
1273 // AZ-261222 if (dp / (TMath::Abs (cand.momxz) + TMath::Abs(bbb->momxz)) > 0.3) { //AZ-060822
1274 if (TMath::Abs(fBy) > 0.2 && dp / (TMath::Abs(cand.momxz) + TMath::Abs(bbb->momxz)) > pxzCut)
1275 { // AZ-311222
1276 // Different curvatures
1277 ok = 0;
1278 }
1279 }
1280 }
1281 if (ok) {
1282 if ((fPass == 0 && aaa.nskips > 0) || (cand.nskips + aaa.nskips > 2))
1283 ok = 0; //!!!AZ-280323 - jumps over stations
1284 if (aaa.stahit.size() == 3) {
1285 // Extension by a triplet
1286 if (ok) {
1287 Double_t dp = TMath::Abs(cand.momxz - aaa.momxz);
1288 // if (dp * 2 / (TMath::Abs (cand.momxz) + TMath::Abs(aaa.momxz)) > 0.3) {
1289 // AZ-070822 if (dp * 2 / (TMath::Abs (cand.momxz) + TMath::Abs(aaa.momxz)) > 0.2) {
1290 // AZ-261222 if (dp / (TMath::Abs (cand.momxz) + TMath::Abs(aaa.momxz)) > 0.3) { //AZ-070822
1291 if (TMath::Abs(fBy) > 0.2 && dp / (TMath::Abs(cand.momxz) + TMath::Abs(aaa.momxz)) > pxzCut)
1292 { // AZ-020123
1293 // different curvatures
1294 ok = 0;
1295 }
1296 }
1297 }
1298 }
1299 if (!ok) {
1300 // No extension found
1301 if ((Int_t)cand.stahit.size() < fNhitsMin[fPass])
1302 continue; // too short track
1303 // cout << " ExtendTrack - store " << cand.stahit.size() << " " << cand.idmaxP << endl;
1304 string codeVec("-");
1305 for (map<int, int>::iterator it = cand.stahit.begin(); it != cand.stahit.end(); ++it)
1306 codeVec += (to_string(it->second) + "-");
1307 if (fCandCodes[cand.stahit.size()].find(codeVec) == fCandCodes[cand.stahit.size()].end()) {
1308 fCandCodes[cand.stahit.size()].insert(codeVec);
1309 if (cand.stahit.size() >= 4) {
1310 // cout << " 2 " << endl;
1311 if (cand.stahit.rbegin()->first == ista)
1312 continue; // AZ-171121
1313 Double_t chi2 = FitTrack(cand);
1314 // cout << " --- Chi2: " << chi2 << endl;
1315 if (chi2 > c2ndfMax)
1316 continue;
1317 }
1318 // cout << " ExtendTrack - store " << cand.stahit.size() << " " << cand.idmaxP << endl;
1319 fCandVec[cand.stahit.begin()->first].push_back(cand);
1320 }
1321 continue; // AZ-280323
1322 }
1323 //*/
1324 //!!! Extension found
1325 candvec candext = cand;
1326 if (cand.stahit.size() < 3)
1327 candext.momxz = bbb->momxz; // update momentum
1328 // candext.ty = bbb.ty;
1329 // Add triplet (or doublet)
1330 for (mitr = aaa.stahit.begin(); mitr != aaa.stahit.end(); ++mitr)
1331 candext.stahit[mitr->first] = mitr->second;
1332 candext.nskips += aaa.nskips; //!!!AZ-280323 - update number of jumps
1333 // cout << " ids: " << aaa.idmaxP << " " << cand.idmaxP << " " << candext.stahit.size() << " " << candext.nskips
1334 // << " " << fPass << endl; if (candext.nskips > 2) exit(0); AZ-120722 if (aaa.idmaxP != cand.idmaxP)
1335 // cand.idmaxP = -1;
1336 if (aaa.idmaxP != cand.idmaxP)
1337 candext.idmaxP = -1; // AZ-120722
1338 // Fit track to fight with combinatorics
1339 if (candext.stahit.size() >= 4) {
1340 Double_t chi2ndf = FitTrack(candext);
1341 // cout << " -3- " << chi2ndf << endl;
1342 if (chi2ndf > c2ndfMax)
1343 continue;
1344 /* This slows down
1345 if (chi2ndf > c2ndfMax) {
1346 string codext = MakeCode(candext);
1347 fCandCodes[candext.stahit.size()].insert(codext);
1348 continue;
1349 }
1350 */
1351 }
1352 ExtendTrack(candext);
1353 } // for (mit = ret.first;
1354}
1355
1356// -------------------------------------------------------------------------
1357
1358Double_t BmnStsVectorFinder::FitTrack(candvec& cand)
1359{
1360 // Fit track candidate
1361
1362 const Int_t gkChi2Cut = 5.0; // 10.0;
1363 // Int_t gkChi2Cut = 10.0; //AZ-010123
1364
1365 CbmStsTrack track;
1366 string hitcode;
1367
1368 MakeStsTrack(cand, hitcode, track);
1369 int nhits = track.GetNStsHits();
1370 if (fCandSet[nhits].find(hitcode) != fCandSet[nhits].end())
1371 return 0.0; // track has been saved already
1372
1373 // AZ-160122 fitter.DoFit(&track);
1374 if (nhits <= 4)
1375 fitter.DoFit(&track);
1376 else
1377 FilterHit(cand, track);
1378 /*
1379 else {
1380 // Debug
1381 CbmStsTrack track1(track);
1382 FilterHit (cand, track);
1383 fitter.DoFit(&track1);
1384 cout << " debug " << track.GetChi2() << " " << track1.GetChi2() << endl;
1385 }
1386 */
1387 // track.GetParamFirst()->Print();
1388 // cout << " aaaaaa " << track.GetChi2() << " " << track.GetNStsHits() << endl;
1389 int ndf = TMath::Max(track.GetNDF(), 1);
1390 Double_t chi2ndf = track.GetChi2() / ndf;
1391 cand.param = *(track.GetParamLast());
1392 cand.chi2 = track.GetChi2();
1393 if (chi2ndf < gkChi2Cut && nhits >= fNhitsMin[fPass]) {
1394 // Refit with material budget
1395 // fitter.DoFit(&track);
1396 fitter.Fit(&track); // AZ
1397 if (track.GetChi2() / ndf > gkChi2Cut) {
1398 // cout << " yyyyyy " << " " << track.GetChi2() / ndf << endl;
1399 // fitter.Fit(&track); //AZ - for debugging
1400 return chi2ndf;
1401 }
1402 // Save good tracks
1403 Double_t qual = nhits;
1404 qual += (100.0 - TMath::Min(TMath::Abs(track.GetChi2()), 100.0)) / 101.0;
1405 fTracks.insert(pair<Double_t, CbmStsTrack>(-qual, track));
1406 fCandSet[nhits].insert(hitcode);
1407 // cout << " FitTrack: " << hitcode << " " << track.GetChi2() << endl;
1408 }
1409 // fCandSet[nhits].insert(hitcode); - don't use this !!!
1410 return chi2ndf;
1411}
1412
1413// -------------------------------------------------------------------------
1414//*
1415Double_t BmnStsVectorFinder::FilterHit(candvec& cand, CbmStsTrack& track)
1416{
1417 // Filter last hit on track candidate
1418
1419 // const Int_t gkChi2Cut = 10.0;
1420 // AZ const Int_t gkChi2Cut = 20.0; //AZ-010123
1421 /*
1422 CbmStsTrack track;
1423 string hitcode;
1424
1425 MakeStsTrack(cand, hitcode, track);
1426 */
1427 Int_t nhits = cand.stahit.size();
1428 // CbmStsHit *hit = (CbmStsHit*) fHitArray->UncheckedAt(track.GetStsHitIndex(nhits-1));
1429 /*cout << " !!! parameters: " << cand.stahit.size() << " " << cand.param.GetX() << " " << cand.param.GetY()
1430 << " " << cand.param.GetZ() << " "
1431 << hit->GetZ() << " " << cand.chi2 << " " << cand.code << endl;*/
1432 // AZ-050222 fitter.Extrapolate(&track,hit->GetZ(),track.GetParamLast());
1433 // cout << hit->GetZ() << " " << track.GetParamLast()->GetZ() << endl;
1434
1435 CbmKFTrack kftr;
1436 kftr.SetStsTrack(track, kFALSE); // last parameter
1437 fitter.SetKFHits(kftr, &track);
1438 Bool_t downstream = kTRUE;
1439 // Double_t qp0 = 0.0;
1440 Double_t qp0 = track.GetParamLast()->GetQp();
1441 // Double_t *cov = kftr.GetCovMatrix();
1442 // for (Int_t j = 0; j < 15; ++j) cov[j] = 0.0;
1443 // cov[0] = cov[2] = cov[5] = cov[9] = cov[14] = 10000.0;
1444 // kftr.GetRefChi2() = 0.0;
1445 // Double_t chi20 = 0.0, dchi2 = 0.0;
1446 // Double_t chi20 = kftr.GetRefChi2(), dchi2 = 0.0;
1447
1448 /*CbmKF *KF = */ CbmKF::Instance();
1449 CbmKFHit* h = NULL;
1450 // Int_t istold = 0, imax = 0;
1451 Bool_t err = kFALSE;
1452
1453 // for (Int_t i = 0; i < nhits; i++) {
1454 h = kftr.GetHit(nhits - 1);
1455 err = err || h->Filter(kftr, downstream, qp0);
1456 // dchi2 = kftr.GetRefChi2() - chi20;
1457 // cout << err << " " << dchi2 << " " << kftr.GetRefChi2() << " " << cand.chi2+dchi2 << endl;
1458 kftr.GetTrackParam(*track.GetParamLast());
1459 track.SetChi2(kftr.GetRefChi2() + err * 1000);
1460 track.SetNDF(2 * nhits - 5);
1461 /*
1462 if (fCandSet[track.GetNStsHits()].find(hitcode) != fCandSet[track.GetNStsHits()].end()) return 0.0; // track has
1463 been saved already
1464
1465 fitter.DoFit(&track);
1466 //track.GetParamFirst()->Print();
1467 //cout << " aaaaaa " << track.GetChi2() << " " << track.GetNStsHits() << endl;
1468 int ndf = TMath::Max (track.GetNDF(), 1);
1469 Double_t chi2ndf = track.GetChi2() / ndf;
1470 if (chi2ndf < gkChi2Cut && track.GetNStsHits() >= fNhitsMin[fPass]) {
1471 // Save good tracks
1472 Double_t qual = track.GetNStsHits();
1473 qual += (100.0 - TMath::Min(TMath::Abs(track.GetChi2()),100.0)) / 101.0;
1474 fTracks.insert(pair<Double_t,CbmStsTrack>(-qual,track));
1475 fCandSet[track.GetNStsHits()].insert(hitcode);
1476 //cout << " FitTrack: " << hitcode << " " << track.GetChi2() << endl;
1477 }
1478 return chi2ndf;
1479 */
1480 return 0.0;
1481}
1482//*/
1483// -------------------------------------------------------------------------
1484
1485void BmnStsVectorFinder::MakeStsTrack(candvec& cand, string& hitcode, CbmStsTrack& track)
1486{
1487 // Fill STS track information
1488
1489 map<Int_t, Int_t>& hitMap = cand.stahit;
1490 CbmStsHit* hit = nullptr;
1491 TArrayI& hits = *track.GetStsHits();
1492 hits.Set(hitMap.size());
1493 Int_t indx = 0 /*, iok = 1*/;
1494 hitcode = "-";
1495
1496 /// hitMap ---- cand.stahit
1497 for (map<Int_t, Int_t>::iterator mit2 = hitMap.begin(); mit2 != hitMap.end(); ++mit2) {
1498 hits[indx++] = mit2->second;
1499 hit = (CbmStsHit*)fHitArray->UncheckedAt(mit2->second);
1500 // cout << " hit: " << mit2->first << " " << mit2->second << endl;
1501 hit->SetDx(0.08 / TMath::Sqrt(12.0)); // AZ-150823 - old
1502 /// hit->SetDx(0.05); //AZ-110123 - early estimate
1503 // hit->SetDx(0.015);
1504 hit->SetDy(0.1234); // AZ-150823 - old
1505 // AZ-030123 if (hit->GetStationNr() <= 3) hit->SetDx(0.02/TMath::Sqrt(12.0));
1506 if (hit->GetStationNr() <= fNSi)
1507 hit->SetDx(0.02 / TMath::Sqrt(12.0)); // AZ-150823 - old
1508 // if (hit->GetStationNr() % 2 != 0) hit->SetCovXY(1.968e-3); // this is for 2-D hits
1509 // else hit->SetCovXY(-1.968e-3);
1510 // hit->SetDy(0.03); // for test
1511 CbmStsStation* stat = CbmStsDigiScheme::Instance()->GetStationByNr(hit->GetStationNr());
1512 Double_t dx = stat->GetSector(0)->GetDx();
1513 // if (dx < 0.02) dx *= 2; // scale up for Si
1514 // AZ-120722 if (dx > 0.01) dx *= 1.2; // for GEMs
1515 if (dx > 0.02)
1516 dx *= 1.2; // AZ-120722 for GEMs
1517 hit->SetDx(dx / TMath::Sqrt(12.0)); // AZ-150823 - old
1518 // if (dx < 0.02) hit->SetDx(0.01); //AZ-040123 - early estimate for Si
1519 /// if (dx < 0.02) hit->SetDx(0.02); //AZ-040123 - early estimate for Si
1520 hitcode += (to_string(mit2->second) + "-");
1521 }
1522 track.SetParamFirst(cand.param);
1523 track.SetParamLast(cand.param);
1524 track.SetChi2(cand.chi2);
1525 cand.code = hitcode;
1526}
1527
1528// -------------------------------------------------------------------------
1529
1530void BmnStsVectorFinder::ExtendTracks(Int_t ista)
1531{
1532 // Extend tracks by going upstream
1533
1534 // const Double_t zsta[19] = {1.132, 32.750, 64.150, 96.750, 129.300, 160.800, 192.100}; // average hit Z-position -
1535 // Run 6 !!! const Double_t zsta[19] = {16.21, 22.79, 29.22, 42.56, 67.63, 116.10, 137.80, 164.10, 186.20}; //
1536 // average hit Z-position - Run 7 !!!
1537 /*const Double_t zsta[19] = {18.55, 27.16, 35.78, 61.05, 91.50, 123.50, 153.90, 185.90, 216.40, 248.20}; // average
1538 hit Z-position - Run 8 !!! const Double_t pars[19][4] = { {1.30358, 20.8953, 0.857369, 0.56339}, {1.64202, 7.18813,
1539 0.596002, 0.694601}, {3.10765, 3.21564, 0.291252, 1.08062}, {2.25484, 3.81810, 0.430724, 0.840519},
1540 {3.38309, 4.33504, 0.349704, 1.07882}, {3.99252, 3.35587, 0.302402, 1.65050},
1541 {10.3321, 4.24742, 0.168995, 1.71113} };*/
1542}
1543
1544// -------------------------------------------------------------------------
1545
1546void BmnStsVectorFinder::FinishEvent()
1547{
1548 // Finish event
1549 /*
1550 FairRootManager* ioman = FairRootManager::Instance();
1551 */
1552}
1553
1554// -------------------------------------------------------------------------
1555
1556void BmnStsVectorFinder::Finish()
1557{
1558 printf("Work time of BmnStsVectorFinder: %4.2f sec.\n", workTime);
1559}
1560
1561//__________________________________________________________________________
1562
1563TVector3 BmnStsVectorFinder::Parabola(TVector3& pos0, TVector3& pos1, TVector3& pos2)
1564{
1565 // Get parabolic track approximation from 3 points (X vs Z)
1566 // y = a*x^2 + bx + c
1567
1568 Double_t x[3] = {pos0[2] - fXyzv[2], pos1[2] - fXyzv[2], pos2[2] - fXyzv[2]};
1569 Double_t y[3] = {pos0[0] - fXyzv[0], pos1[0] - fXyzv[0], pos2[0] - fXyzv[0]};
1570
1571 Double_t denom = (x[0] - x[1]) * (x[0] - x[2]) * (x[1] - x[2]);
1572 Double_t dy10 = y[1] - y[0];
1573 Double_t dy02 = y[0] - y[2];
1574 Double_t dy21 = y[2] - y[1];
1575 // Double_t a = x[2] * (y[1] - y[0]) + x[1] * (y[0] - y[2]) + x[0] * (y[2] - y[1]);
1576 Double_t a = x[2] * dy10 + x[1] * dy02 + x[0] * dy21;
1577 a /= denom;
1578 // Double_t b = x[0]*x[0] * (y[1] - y[2]) + x[2]*x[2] * (y[0] - y[1]) +
1579 // x[1]*x[1] * (y[2] - y[0]);
1580 Double_t b = -x[0] * x[0] * dy21 - x[2] * x[2] * dy10 - x[1] * x[1] * dy02;
1581 b /= denom;
1582 Double_t c = x[1] * x[1] * (x[2] * y[0] - x[0] * y[2]) + x[1] * (x[0] * x[0] * y[2] - x[2] * x[2] * y[0])
1583 + x[0] * x[2] * (x[2] - x[0]) * y[1];
1584 c /= denom;
1585 return TVector3(c, b, a);
1586 // 𝐴=𝑥3(𝑦2−𝑦1)+𝑥2(𝑦1−𝑦3)+𝑥1(𝑦3−𝑦2) / (𝑥1−𝑥2)(𝑥1−𝑥3)(𝑥2−𝑥3)
1587 // 𝐵=𝑥21(𝑦2−𝑦3)+𝑥23(𝑦1−𝑦2)+𝑥22(𝑦3−𝑦1) / (𝑥1−𝑥2)(𝑥1−𝑥3)(𝑥2−𝑥3)
1588 // 𝐶=𝑥22(𝑥3𝑦1−𝑥1𝑦3)+𝑥2(𝑥21𝑦3−𝑥23𝑦1)+𝑥1𝑥3(𝑥3−𝑥1)𝑦2 / (𝑥1−𝑥2)(𝑥1−𝑥3)(𝑥2−𝑥3)
1589}
1590
1591// -------------------------------------------------------------------------
1592
1593void BmnStsVectorFinder::FitTracks()
1594{
1595 // Fit tracks
1596
1597 const Double_t gkChi2Cut = 5.0; // 10.0;
1598 // const Double_t gkChi2Cut = 10.0; //AZ-010123
1599 // fTracks.clear();
1600
1601 // for (Int_t ista = 0; ista < fNsta-3; ++ista) {
1602 for (Int_t ista = 0; ista < fNsta; ++ista) {
1603 Int_t ntra = fCandVec[ista].size();
1604 // cout << " FitTracks: " << ista << " " << ntra << endl;
1605
1606 for (Int_t itra = 0; itra < ntra; ++itra) {
1607 candvec& cand = fCandVec[ista][itra];
1608 // if (cand.nextind >= 0) continue; // track has an extension already processed
1609
1610 map<Int_t, Int_t>& hitMap = cand.stahit;
1611 Int_t nhits = hitMap.size();
1612 // if (nhits < 3) continue; // too short track
1613 if (nhits < fNhitsMin[fPass])
1614 continue; // too short track
1615 CbmStsTrack track;
1616 TArrayI& hits = *track.GetStsHits();
1617 hits.Set(nhits);
1618 Int_t indx = 0, iok = 1;
1619
1620 // Check if the track candidate has been checked already
1621 string hitcode("-");
1622 for (map<Int_t, Int_t>::iterator mit = hitMap.begin(); mit != hitMap.end(); ++mit)
1623 hitcode += (to_string(mit->second) + "-");
1624 // cout << itra << " " << hitcode << endl;
1625
1626 // Int_t sta0 = hitMap.begin()->first, same = 0;
1627
1628 // same = fCandSet[sta0].count(hitcode);
1629 // same = fCandSet[nhits].count(hitcode);
1630 // if (same) cout << " same: " << sta0 << " " << same << " " << endl;
1631
1632 /////AZ if (same) continue; // used combination - skip it ???
1633 // fCandSet[sta0].insert(hitcode);
1634 fCandSet[nhits].insert(hitcode);
1635
1636 for (map<Int_t, Int_t>::iterator mit = hitMap.begin(); mit != hitMap.end(); ++mit) {
1637 hits[indx++] = mit->second;
1638 // Debug
1639 CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(mit->second);
1640 // cout << " hit: " << mit->first << " " << hit->GetX() << " " << hit->GetY() << " " << hit->GetZ() << "
1641 // "
1642 // << hit->GetDx() << " " << hit->GetDy() << endl;
1643 // cout << " hit: " << mit->first << " " << mit->second << endl;
1644 // hit->SetDx(0.08/TMath::Sqrt(12.0));
1645 /// hit->SetDx(0.05); //AZ-110123 - early estimate
1646 // hit->SetDx(0.015);
1647 /// hit->SetDy(0.1234);
1648 // AZ-030123 if (hit->GetStationNr() <= 3) hit->SetDx(0.02/TMath::Sqrt(12.0));
1649 // if (hit->GetStationNr() <= fNSi) hit->SetDx(0.02/TMath::Sqrt(12.0)); //AZ-030123
1650 // if (hit->GetStationNr() <= 3) hit->SetDx(0.0025);
1651 // if (hit->GetStationNr() % 2 != 0) hit->SetCovXY(1.968e-3); // this is for 2-D hits (CbmFitter)
1652 // else hit->SetCovXY(-1.968e-3);
1653 // hit->SetDy(0.03); // for test
1654 CbmStsStation* stat = CbmStsDigiScheme::Instance()->GetStationByNr(hit->GetStationNr());
1655 Double_t dx = stat->GetSector(0)->GetDx();
1656 // if (dx < 0.02) dx *= 2; // scale up for Si
1657 // AZ-120722 if (dx > 0.01) dx *= 1.2; // scale up for GEMs
1658 if (dx > 0.02)
1659 dx *= 1.2; // AZ-120722 scale up for GEMs
1660 hit->SetDx(dx / TMath::Sqrt(12.0)); // AZ-150823 - old
1661 // if (dx < 0.02) hit->SetDx(0.01); //AZ-040123 - early estimate for Si
1662 /// if (dx < 0.02) hit->SetDx(0.02); //AZ-040123 - early estimate for Si
1663 /*
1664 if (fMatBudgetFileName != "" && hit->GetDz() < 1.e-7) {
1665 // Store material budget for the hit (along the beam)
1666 Double_t zshift = (dx < 0.02) ? 0.5 : 1.5; // Si or GEM
1667 TProfile2D *prof = fMatHistos.lower_bound(hit->GetZ()-zshift)->second;
1668 int ix = prof->FindBin(hit->GetX());
1669 int iy = prof->FindBin(hit->GetY());
1670 Double_t radThick = prof->GetBinContent(ix,iy);
1671 hit->SetDz(radThick/100);
1672 }
1673 */
1674 }
1675
1676 // cout << " before " << track.GetParamFirst()->GetZ() << endl;
1677 // fitter.DoFit(&track);
1678 fitter.Fit(&track); // AZ
1679 // cout << " after " << track.GetParamFirst()->GetZ() << " " << track.GetNStsHits() << " " <<
1680 // track.GetChi2() << endl;
1681 /*
1682 CbmKFTrack kftr = CbmKFTrack(track);
1683 fitter.SetKFHits(kftr, &track);
1684 kftr.Fit(kTRUE); // downstream
1685 kftr.Fit(kFALSE); // upstream
1686 */
1687 if (fVerbose > 0) {
1688 cout << " aaaaaaa " << endl;
1689 // track.GetParamFirst()->Print();
1690 cout << track.GetChi2() << " " << nhits << " " << endl;
1691 }
1692 // AZ-111221 if (track.GetChi2() / track.GetNDF() > gkChi2Cut) {
1693 if (track.GetChi2() / track.GetNDF() > gkChi2Cut && track.GetNStsHits() > 3) {
1694 /*
1695 // Too high chi2/NDF - try to remove one outlier with the largest dChi2
1696 // Based on the code from CbmKFTrackInterface::Fit()
1697 CbmKFTrack kftr;
1698 //kftr.SetStsTrack(track, kFALSE); // last parameter
1699 //Bool_t downstream = kFALSE;
1700 kftr.SetStsTrack(track, kTRUE); // first parameter
1701 fitter.SetKFHits(kftr,&track);
1702 Bool_t downstream = kTRUE;
1703 Double_t qp0 = 0.0;
1704 Double_t *cov = kftr.GetCovMatrix();
1705 for (Int_t j = 0; j < 15; ++j) cov[j] = 0.0;
1706 cov[0] = cov[2] = cov[5] = cov[9] = cov[14] = 10000.0;
1707 kftr.GetRefChi2() = 0.0;
1708 Double_t chi20 = 0.0, dchi2 = 0.0, dchi2max = 0.0;
1709
1710 CbmKF *KF = CbmKF::Instance();
1711 CbmKFHit *h = NULL;
1712 Int_t istold = 0, imax = 0;
1713 Bool_t err = kFALSE;
1714 set<Int_t> indok;
1715
1716 //for (Int_t i = nhits - 1; i >= 0; i--) {
1717 // h = kftr.GetHit( i );
1718 // Int_t ist = h->MaterialIndex;
1719 // if (i < nhits - 1) {
1720 // for (Int_t j = istold - 1; j > ist; j--)
1721 // err = err || KF->vMaterial[j]->Pass( kftr, downstream, qp0 );
1722 // }
1723
1724 for (Int_t i = 0; i < nhits; i++) {
1725 h = kftr.GetHit( i );
1726 Int_t ist = h->MaterialIndex;
1727 if (i > 0) {
1728 for (Int_t j = istold + 1; j < ist; j++)
1729 err = err || KF->vMaterial[j]->Pass( kftr, downstream, qp0 );
1730 }
1731 err = err || h->Filter( kftr, downstream, qp0 );
1732 dchi2 = kftr.GetRefChi2() - chi20;
1733 cout << dchi2 << " " << kftr.GetRefChi2() << endl;
1734 //if (dchi2 <= gkChi2Cut) indok.insert(i);
1735 //if (dchi2 <= 2*gkChi2Cut) indok.insert(i);
1736 indok.insert(i);
1737 if (dchi2 > dchi2max) {
1738 dchi2max = dchi2;
1739 imax = i;
1740 }
1741 chi20 = kftr.GetRefChi2();
1742 istold = ist;
1743 }
1744 */
1745 // Remove last hit and refit
1746 set<Int_t> indok;
1747 for (Int_t i = 0; i < nhits; i++)
1748 indok.insert(i);
1749 int imax = nhits - 1;
1750
1751 if (indok.size() > 3) {
1752 // if (indok.size() >= fNhitsMin[fPass]) {
1753 // if (indok.size()-1 >= fNhitsMin[fPass]) {
1754 indok.erase(imax);
1755 // Good short track - refit
1756 indx = 0;
1757 hitcode = "-";
1758 for (set<Int_t>::iterator sit = indok.begin(); sit != indok.end(); ++sit) {
1759 hits[indx++] = hits[*sit];
1760 // hitcode += hits[*sit];
1761 hitcode += (to_string(hits[*sit]) + "-");
1762 }
1763 // if (fCandSet[sta0].count(hitcode)) continue;
1764 // fCandSet[sta0].insert(hitcode);
1765 // cout << " failed " << hitcode << endl;
1766 if (fCandSet[indok.size()].find(hitcode) != fCandSet[indok.size()].end())
1767 continue;
1768 fCandSet[indok.size()].insert(hitcode);
1769
1770 hits.Set(indok.size());
1771 // fitter.DoFit(&track);
1772 fitter.Fit(&track); // AZ
1773 int ndf = (track.GetFlag() == 0) ? track.GetNDF() : 1;
1774 // cout << " bbbbbb " << endl;
1775 // track.GetParamFirst()->Print();
1776 if (fVerbose > 0)
1777 cout << " bbbbbb " << track.GetChi2() << " " << hits.GetSize() << " " << track.GetChi2() / ndf
1778 << endl;
1779 if (track.GetChi2() / track.GetNDF() > gkChi2Cut)
1780 iok = 0;
1781 } else
1782 iok = 0;
1783
1784 } // if (track.GetChi2() / track.GetNDF() > gkChi2Cut)
1785
1786 // Good track - extra printout
1787 if (iok) {
1788 Double_t qual = track.GetNStsHits();
1789 // qual += (100.0 - TMath::Min(TMath::Abs(track.GetChi2()),100.0)) / 101.0;
1790 // if (track.GetNStsHits() <= 4) qual += TMath::Min(1/TMath::Abs(track.GetParamFirst()->GetQp()),10.0)
1791 // / 11.0;
1792 if (track.GetNStsHits() < 4) {
1793 // Vertex constraint
1794 /*
1795 TVector3 mom3, pos3;
1796 track.GetParamFirst()->Momentum(mom3);
1797 track.GetParamFirst()->Position(pos3);
1798 //Double_t ang = mom3.Angle(pos3);
1799 Double_t phi1 = TMath::ATan2 (track.GetParamFirst()->GetY(), track.GetParamFirst()->GetX());
1800 Double_t phi2 = TMath::ATan2 (track.GetParamFirst()->GetTy(), track.GetParamFirst()->GetTx());
1801 Double_t dphi = phi1 - Proxim(phi1,phi2);
1802 //qual += (2.1 - TMath::Min(TMath::Abs(dphi),2.1)) / 2.2;
1803 //Double_t www = TMath::Exp(-TMath::Abs(dphi)/2.1);
1804 Double_t www = TMath::Min (TMath::Exp(-TMath::Abs(dphi)/0.2), 0.999); // 0.2 - rms
1805 if (TMath::Abs(dphi) > 2.1) www = 0.0;
1806 */
1807 // Linear extrapolation to Zvert
1808 Double_t xextr; // = pos3.X() - pos3.Z() * track.GetParamFirst()->GetTx() + 0.04;
1809 Double_t yextr; // = pos3.Y() - pos3.Z() * track.GetParamFirst()->GetTy();
1810 // Extrapolation to Zvert
1811 FairTrackParam param = *track.GetParamFirst();
1812 fitter.Extrapolate(&param, fXyzv[2], &param);
1813 xextr = param.GetX() - fXyzv[0];
1814 yextr = param.GetY() - fXyzv[1];
1815 Double_t rad =
1816 TMath::Sqrt(xextr * xextr / 0.035 / 0.035 + yextr * yextr / 0.2 / 0.2); // 0.035, 0.2 - sigmas
1817 Double_t www = TMath::Min(TMath::Exp(-rad), 0.999);
1818 qual += TMath::Exp(-TMath::Abs(track.GetChi2())) * www;
1819 } else
1820 qual += (100.0 - TMath::Min(TMath::Abs(track.GetChi2()), 100.0)) / 101.0;
1821 fTracks.insert(pair<Double_t, CbmStsTrack>(-qual, track));
1822 // AZ-130722 pair<multimap<Int_t,Int_t>::iterator,multimap<Int_t,Int_t>::iterator> ret;
1823 pair<unordered_multimap<Int_t, Int_t>::iterator, unordered_multimap<Int_t, Int_t>::iterator>
1824 ret; // AZ-130722
1825 nhits = hits.GetSize();
1826 if (fVerbose > 0) {
1827 cout << " Good track: " << endl;
1828 // track.GetParamFirst()->Print();
1829 cout << track.GetChi2() << " " << nhits << " " << track.GetParamFirst()->GetZ() << endl;
1830
1831 for (Int_t j = nhits - 1; j >= 0; --j) {
1832 ret = fHit2id.equal_range(hits[j]);
1833 cout << " (";
1834 // AZ-130722 for (multimap<Int_t,Int_t>::iterator mit = ret.first; mit != ret.second; ++mit) {
1835 for (unordered_multimap<Int_t, Int_t>::iterator mit = ret.first; mit != ret.second; ++mit) {
1836 if (mit == ret.first)
1837 cout << mit->first << "*";
1838 else
1839 cout << ":";
1840 cout << mit->second;
1841 }
1842 cout << ")";
1843 }
1844 cout << "\n";
1845 }
1846 }
1847
1848 } // for (Int_t itra = 0;
1849 }
1850}
1851
1852// -------------------------------------------------------------------------
1853
1854void BmnStsVectorFinder::RemoveDoubles()
1855{
1856 // Remove double tracks (keep the ones with better quality)
1857
1858 multimap<Double_t, CbmStsTrack>::iterator mit = fTracks.begin(), mit1;
1859 Int_t nOK = fTracks.size(), nOut = fVectorArray->GetEntriesFast();
1860
1861 for (; mit != fTracks.end(); ++mit) {
1862 if (mit->second.GetFlag())
1863 continue; // discarded track
1864 mit1 = mit;
1865 ++mit1;
1866
1867 for (; mit1 != fTracks.end(); ++mit1) {
1868 if (mit1->second.GetFlag())
1869 continue; // discarded track
1870 if (!AreTracksDoubles(mit1->second, mit->second))
1871 continue;
1872 mit1->second.SetFlag(1); // double - discard it
1873 --nOK;
1874 /*
1875 // Debug printout
1876 TArrayI &hits = *mit1->second.GetStsHits();
1877 Int_t nhits = mit1->second.GetNStsHits();
1878 pair<multimap<Int_t,Int_t>::iterator,multimap<Int_t,Int_t>::iterator> ret;
1879 cout << " Double: " << mit->first << " " << mit1->first << endl;
1880 for (Int_t j = nhits-1; j >= 0; --j) {
1881 ret = fHit2id.equal_range(hits[j]);
1882 cout << " (";
1883 for (multimap<Int_t,Int_t>::iterator mit2 = ret.first; mit2 != ret.second; ++mit2) {
1884 if (mit2 == ret.first) cout << mit2->first << "*";
1885 else cout << ":";
1886 cout << mit2->second;
1887 }
1888 cout << ")";
1889 }
1890 cout << "\n";
1891 */
1892 }
1893
1894 new ((*fVectorArray)[nOut++]) CbmStsTrack(mit->second);
1895 // cout << " yyyyy " << mit->second.GetParamFirst()->GetZ() << " " << mit->second.GetNStsHits() << endl;
1896 }
1897 if (fVerbose > 0)
1898 cout << " RemoveDoubles: " << fTracks.size() << " " << nOK << endl;
1899}
1900
1901// -------------------------------------------------------------------------
1902
1903Bool_t BmnStsVectorFinder::AreTracksDoubles(CbmStsTrack& tr1, CbmStsTrack& tr2)
1904{
1905 /// Searching common hits in 2 tracks to determine doubles
1906 /// tr1 contains fewer hits than tr2
1907
1908 Int_t limCommonPoint = (tr1.GetNStsHits() + 1) / 2; // at least so many common hits should be found
1909 // limCommonPoint = TMath::Min (limCommonPoint, 2); // 2 common hits at max
1910 limCommonPoint = TMath::Min(limCommonPoint, 1); // 1 common hit at max
1911 TArrayI &hits1 = *tr1.GetStsHits(), &hits2 = *tr2.GetStsHits();
1912 Int_t nh1 = hits1.GetSize(), nh2 = hits2.GetSize(), nHitsCommon = 0, j = nh2 - 1;
1913 if (nh1 == 3)
1914 limCommonPoint = 0;
1915
1916 for (Int_t i = nh1 - 1; i >= 0; i--) {
1917 // for (Int_t i = 0; i <= nh1 - 1; ++i){
1918 CbmStsHit* hit1 = (CbmStsHit*)fHitArray->UncheckedAt(tr1.GetStsHitIndex(i));
1919
1920 for (; j >= 0; j--) {
1921 // for ( ; j <= nh2 - 1; ++j){
1922 CbmStsHit* hit2 = (CbmStsHit*)fHitArray->UncheckedAt(tr2.GetStsHitIndex(j));
1923
1924 // Is the hit common for two tracks?
1925 if (hit1 == hit2) {
1926 nHitsCommon++;
1927 // AZ-281121 if (nHitsCommon == limCommonPoint) return kTRUE; // already enough common hits
1928 if (nHitsCommon > limCommonPoint)
1929 return kTRUE; // already enough common hits
1930 break;
1931 }
1932
1933 if (hit2->GetStationNr() < hit1->GetStationNr())
1934 break; // already closer to target
1935 }
1936
1937 // AZ-281121 if ((nh1 - i - 1) + limCommonPoint - nHitsCommon > nh1) return kFALSE; // there'll be not enough
1938 // common hits already
1939 if (i + nHitsCommon <= limCommonPoint)
1940 return kFALSE; // there'll be not enough common hits already
1941 }
1942
1943 // AZ-281121 if (nHitsCommon < limCommonPoint) return kFALSE; // not too many common hits
1944 if (nHitsCommon <= limCommonPoint)
1945 return kFALSE; // not too many common hits
1946
1947 return kTRUE;
1948}
1949
1950// -------------------------------------------------------------------------
1951/*
1952void BmnStsVectorFinder::RemoveFakes()
1953{
1954 // Remove fake tracks (according to some empirical rules)
1955
1956 Int_t ncand = fVectorArray->GetEntriesFast();
1957
1958 for (Int_t j = 0; j < ncand; ++j) {
1959 CbmStsTrack *tr = (CbmStsTrack*) fVectorArray->UncheckedAt(j);
1960 Int_t nhits = tr->GetNStsHits();
1961 if (1/TMath::Abs(tr->GetParamFirst()->GetQp()) < 0.2) { fVectorArray->Remove(tr); continue; }
1962 if (nhits < 4) {
1963 CbmStsHit *hit = (CbmStsHit*) fHitArray->UncheckedAt(tr->GetStsHitIndex(0));
1964 if (hit->GetStationNr() <= 3) continue;
1965 fVectorArray->Remove(tr); // short track starts in GEM
1966 }
1967 }
1968 fVectorArray->Compress();
1969
1970 cout << " Remove fakes: " << ncand << " " << fVectorArray->GetEntriesFast() << endl;
1971}
1972*/
1973// -------------------------------------------------------------------------
1974void BmnStsVectorFinder::RemoveFakes()
1975{
1976 // Remove fake tracks (according to some empirical rules)
1977
1978 multimap<Double_t, CbmStsTrack>::iterator mit = fTracks.begin();
1979 Int_t nOK = fTracks.size(), nOut = nOK;
1980
1981 for (; mit != fTracks.end(); ++mit) {
1982 if (mit->second.GetFlag())
1983 continue; // discarded track
1984 CbmStsTrack* tr = &mit->second;
1985 Int_t nhits = tr->GetNStsHits();
1986 if (1 / TMath::Abs(tr->GetParamFirst()->GetQp()) < 0.2) {
1987 tr->SetFlag(1);
1988 --nOut;
1989 continue;
1990 }
1991 if (nhits < 4) {
1992 CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(tr->GetStsHitIndex(0));
1993 if (hit->GetStationNr() <= 3)
1994 continue;
1995 tr->SetFlag(1); // short track starts in GEM
1996 --nOut;
1997 }
1998 }
1999
2000 if (fVerbose > 0)
2001 cout << " Remove fakes: " << nOK << " " << nOut << endl;
2002}
2003
2004// -------------------------------------------------------------------------
2005
2006void BmnStsVectorFinder::ExcludeFakes()
2007{
2008 // Exclude fake tracks (with too many shared clusters)
2009
2010 Int_t ncand = fVectorArray->GetEntriesFast();
2011
2012 map<Int_t, set<Int_t>> mClusTr;
2013 map<Int_t, Double_t> mWeight;
2014
2015 for (Int_t j = 0; j < ncand; ++j) {
2016 CbmStsTrack* tr = (CbmStsTrack*)fVectorArray->UncheckedAt(j);
2017 Int_t nhits = tr->GetNStsHits();
2018
2019 for (Int_t jh = 0; jh < nhits; ++jh) {
2020 Int_t indx = tr->GetStsHitIndex(jh);
2021 CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(indx);
2022
2023 for (Int_t side = 0; side < 2; ++side) {
2024 Int_t iclus = hit->GetDigi(side);
2025 if (mClusTr.count(iclus) == 0) {
2026 set<Int_t> aaa;
2027 mClusTr[iclus] = aaa;
2028 }
2029 mClusTr[iclus].insert(j);
2030 }
2031 }
2032
2033 Double_t www = TMath::Exp(-tr->GetChi2() / tr->GetNDF());
2034 // Extrapolation to Zvert
2035 if (nhits == 3) {
2036 FairTrackParam param = *tr->GetParamFirst();
2037 fitter.Extrapolate(&param, fXyzv[2], &param);
2038 Double_t xextr = param.GetX() - fXyzv[0];
2039 Double_t yextr = param.GetY() - fXyzv[1];
2040 Double_t rad =
2041 TMath::Sqrt(xextr * xextr / 0.035 / 0.035 + yextr * yextr / 0.2 / 0.2); // 0.035, 0.2 - sigmas
2042 www *= TMath::Exp(-rad);
2043 }
2044 mWeight[j] = www;
2045 }
2046
2047 while (1) {
2048 multimap<Double_t, pair<Int_t, Int_t>> mapQual;
2049 set<Int_t> set2kill, set2kill1;
2050
2051 for (Int_t j = 0; j < ncand; ++j) {
2052 CbmStsTrack* tr = (CbmStsTrack*)fVectorArray->UncheckedAt(j);
2053 if (tr == NULL)
2054 continue; // removed track
2055 Int_t nhits = tr->GetNStsHits();
2056 Int_t nover = 0;
2057
2058 for (Int_t jh = 0; jh < nhits; ++jh) {
2059 Int_t indx = tr->GetStsHitIndex(jh);
2060 CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(indx);
2061 // if (mClusTr[hit->GetDigi(0)].size() > 1 || mClusTr[hit->GetDigi(1)].size() > 1) ++nover;
2062
2063 for (Int_t side = 0; side < 2; ++side) {
2064 Int_t iclus = hit->GetDigi(side);
2065 if (mClusTr[iclus].size() > 1)
2066 ++nover;
2067 }
2068 }
2069 // Double_t quality = nhits * 1000 + (1 - Double_t(nover) / nhits) * 100
2070 //+ (1 - TMath::Min(tr->GetChi2(),100.0) / 100);
2071 Double_t quality = nhits * 1000;
2072 Double_t www = TMath::Exp(-Double_t(nover) / nhits / 2) * mWeight[j];
2073 mapQual.insert(pair<Double_t, pair<Int_t, Int_t>>(quality + www, pair<Int_t, Int_t>(j, nover)));
2074 // if ((nover == nhits) || (nhits == 3 && nover == 2)) set2kill.insert(j);
2075 // if ((nhits <= 4 && nover == nhits) || (nhits == 3 && nover == 2)) set2kill.insert(j);
2076 // if (nhits <= 4 && Double_t(nover)/nhits/2 > 0.81) set2kill.insert(j);
2077 if (nhits <= 5 && Double_t(nover) / nhits / 2 > 0.6)
2078 set2kill.insert(j);
2079 else if (!fExact && nhits <= 3 && www < 1.e-8)
2080 set2kill1.insert(j); // 1e-8 - empirical value
2081 }
2082
2083 if (set2kill.size())
2084 set2kill1.clear(); // first kill in set2kill
2085 if (set2kill.size() || set2kill1.size()) {
2086 // Remove track and update cluster-to-track info
2087
2088 for (multimap<Double_t, pair<Int_t, Int_t>>::iterator mit = mapQual.begin(); mit != mapQual.end(); ++mit) {
2089 Int_t itr = mit->second.first;
2090 CbmStsTrack* tr = (CbmStsTrack*)fVectorArray->UncheckedAt(itr);
2091 Int_t nhits = tr->GetNStsHits() /*, nover = mit->second.second*/;
2092 if (set2kill.count(itr) == 0 && set2kill1.count(itr) == 0)
2093 continue;
2094
2095 for (Int_t jh = 0; jh < nhits; ++jh) {
2096 Int_t indx = tr->GetStsHitIndex(jh);
2097 CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(indx);
2098 if (mClusTr[hit->GetDigi(0)].count(itr) == 0 || mClusTr[hit->GetDigi(1)].count(itr) == 0)
2099 // AZ-300822 { cout << " !!! No track found: " << itr << endl; exit(0); }
2100 {
2101 if (fVerbose > 0)
2102 cout << " !!! No track found: " << itr << endl; /*exit(0);*/
2103 } // AZ-300822
2104 mClusTr[hit->GetDigi(0)].erase(itr);
2105 mClusTr[hit->GetDigi(1)].erase(itr);
2106 }
2107 fVectorArray->Remove(tr);
2108 break;
2109 }
2110 } else
2111 break;
2112 } // while (1)
2113
2114 fVectorArray->Compress();
2115
2116 if (fVerbose > 0)
2117 cout << " Exclude fakes: " << ncand << " " << fVectorArray->GetEntriesFast() << endl;
2118
2119 // AZ-310723 - compute number of explicitly owned track hits and number of missing stations
2120 /*
2121 ncand = fVectorArray->GetEntriesFast();
2122 map<int,int> usedHits;
2123
2124 for (int itr = 0; itr < ncand; ++itr) {
2125 CbmStsTrack *tr = (CbmStsTrack*) fVectorArray->UncheckedAt(itr);
2126 Int_t nhits = tr->GetNStsHits(), stat0 = 0, stat1 = 0, nmiss = 0;
2127 Double_t edgeXY[2] = {9, 9};
2128
2129 for (Int_t jh = 0; jh < nhits; ++jh) {
2130 Int_t indx = tr->GetStsHitIndex(jh);
2131 usedHits[indx]++;
2132 CbmStsHit *hit = (CbmStsHit*) fHitArray->UncheckedAt(indx);
2133 if (jh == 0) {
2134 stat0 = hit->GetStationNr() - 1;
2135 if (hit->GetStationNr() <= fNSi) {
2136 // Silicon
2137 BmnSiliconStation *stat = fSilStationSet->GetStation(stat0);
2138 //cout << hit->GetZ() << " jjj " << tr->GetParamFirst()->GetZ() << endl;
2139
2140 for (int imod = 0; imod < stat->GetNModules(); ++imod) {
2141 BmnSiliconModule *mod = stat->GetModule(imod);
2142 if (mod->IsPointInsideModule(hit->GetX(),hit->GetY())) {
2143 //cout << " +++ " << hit->GetX() << " " << hit->GetY() << " "
2144 // << mod->GetXMinModule() << " " << mod->GetXMaxModule() << " "
2145 // << mod->GetYMinModule() << " " << mod->GetYMaxModule() << endl;
2146 edgeXY[0] = hit->GetX() - mod->GetXMinModule();
2147 Double_t ddx = hit->GetX() - mod->GetXMaxModule();
2148 if (TMath::Abs(ddx) < TMath::Abs(edgeXY[0])) edgeXY[0] = ddx;
2149 edgeXY[1] = hit->GetY() - mod->GetYMinModule();
2150 Double_t ddy = hit->GetY() - mod->GetYMaxModule();
2151 if (TMath::Abs(ddy) < TMath::Abs(edgeXY[1])) edgeXY[1] = ddy;
2152 break;
2153 }
2154 }
2155 }
2156 } else if (jh == nhits - 1) stat1 = hit->GetStationNr() - 1;
2157 }
2158 nmiss = stat1 - stat0 - nhits + 1;
2159 //tr->SetEdges (edgeXY[0], edgeXY[1]);
2160
2161 if (stat0 > 0) {
2162 // Check upstream stations
2163 FairTrackParam param = *tr->GetParamFirst();
2164
2165 for (int ista = 0; ista < stat0; ++ista) {
2166 int igem = stat0 - fNSi - 1 - ista;
2167 if (igem >= 0) {
2168 BmnGemStripStation *stat = fGemStationSet->GetStation(igem);
2169 fitter.Extrapolate (&param, stat->GetZPosition(), &param);
2170 //if (param.GetX() < stat->GetXMinStation() || param.GetX() > stat->GetXMaxStation()) continue;
2171 //if (param.GetY() < stat->GetYMinStation() || param.GetY() > stat->GetYMaxStation()) continue;
2172
2173 for (int imod = 0; imod < stat->GetNModules(); ++imod) {
2174 BmnGemStripModule *mod = stat->GetModule(imod);
2175 if (mod->IsPointInsideModule(param.GetX(),param.GetY())) {
2176 ++nmiss;
2177 //cout << " +++ " << param.GetX() << " " << param.GetY() << endl;
2178 break;
2179 } //else cout << " --- " << param.GetX() << " " << param.GetY() << endl;
2180 }
2181 } else {
2182 BmnSiliconStation *stat = fSilStationSet->GetStation(stat0-1-ista);
2183 fitter.Extrapolate (&param, stat->GetZPosition(), &param);
2184 //if (param.GetX() < stat->GetXMinStation() || param.GetX() > stat->GetXMaxStation()) continue;
2185 //if (param.GetY() < stat->GetYMinStation() || param.GetY() > stat->GetYMaxStation()) continue;
2186
2187 for (int imod = 0; imod < stat->GetNModules(); ++imod) {
2188 BmnSiliconModule *mod = stat->GetModule(imod);
2189 if (mod->IsPointInsideModule(param.GetX(),param.GetY())) {
2190 ++nmiss;
2191 //cout << " +++ " << param.GetX() << " " << param.GetY() << endl;
2192 break;
2193 } //else cout << " --- " << param.GetX() << " " << param.GetY() << endl;
2194 }
2195 }
2196 }
2197 } else if (stat1 < fNsta - 1) {
2198 // Check downstream stations
2199 FairTrackParam param = *tr->GetParamLast();
2200
2201 for (int ista = stat1 + 1; ista < fNsta; ++ista) {
2202 int igem = ista - fNSi;
2203 if (igem >= 0) {
2204 BmnGemStripStation *stat = fGemStationSet->GetStation(igem);
2205 fitter.Extrapolate (&param, stat->GetZPosition(), &param);
2206 //if (param.GetX() < stat->GetXMinStation() || param.GetX() > stat->GetXMaxStation()) continue;
2207 //if (param.GetY() < stat->GetYMinStation() || param.GetY() > stat->GetYMaxStation()) continue;
2208
2209 for (int imod = 0; imod < stat->GetNModules(); ++imod) {
2210 BmnGemStripModule *mod = stat->GetModule(imod);
2211 if (mod->IsPointInsideModule(param.GetX(),param.GetY())) {
2212 ++nmiss;
2213 //cout << " +++ " << param.GetX() << " " << param.GetY() << endl;
2214 break;
2215 } //else cout << " --- " << param.GetX() << " " << param.GetY() << endl;
2216 }
2217 } else {
2218 BmnSiliconStation *stat = fSilStationSet->GetStation(ista);
2219 fitter.Extrapolate (&param, stat->GetZPosition(), &param);
2220 //if (param.GetX() < stat->GetXMinStation() || param.GetX() > stat->GetXMaxStation()) continue;
2221 //if (param.GetY() < stat->GetYMinStation() || param.GetY() > stat->GetYMaxStation()) continue;
2222
2223 for (int imod = 0; imod < stat->GetNModules(); ++imod) {
2224 BmnSiliconModule *mod = stat->GetModule(imod);
2225 if (mod->IsPointInsideModule(param.GetX(),param.GetY())) {
2226 ++nmiss;
2227 //cout << " +++ " << param.GetX() << " " << param.GetY() << endl;
2228 break;
2229 } //else cout << " --- " << param.GetX() << " " << param.GetY() << endl;
2230 }
2231 }
2232 }
2233 } // if (stat0 > 0)
2234 //tr->SetFlag(nmiss); //AZ - fFlag is used (in global tracking?)
2235 tr->SetStsEv(nmiss);
2236 } // for (int itr = 0; itr < ncand;
2237
2238 for (int itr = 0; itr < ncand; ++itr) {
2239 CbmStsTrack *tr = (CbmStsTrack*) fVectorArray->UncheckedAt(itr);
2240 Int_t nhits = tr->GetNStsHits(), nown = 0;
2241 Double_t weight = 0.0;
2242
2243 for (Int_t jh = 0; jh < nhits; ++jh) {
2244 Int_t indx = tr->GetStsHitIndex(jh);
2245 if (usedHits[indx] == 1) ++nown;
2246 CbmStsHit *hit = (CbmStsHit*) fHitArray->UncheckedAt(indx);
2247 if (hit->GetStationNr() <= fNSi) {
2248 if (TMath::Abs(hit->GetSignalDiv() - 0.1238) > 0.1245*2) weight += 0.5;
2249 else weight += 1.0;
2250 } else {
2251 if (TMath::Abs(hit->GetSignalDiv() + 0.0766) > 0.0749*3) weight += 0.5;
2252 else weight += 1.0;
2253 }
2254 }
2255 //tr->SetFlag (100 * nown + tr->GetFlag());
2256 tr->SetStsEv (100 * nown + tr->GetStsEv());
2257 tr->SetB (weight);
2258 }
2259 */
2260}
2261
2262// -------------------------------------------------------------------------
2263
2264Double_t BmnStsVectorFinder::DxVsMom(Int_t ista, candvec& vec)
2265{
2266 // Compute momentum correction for parabolic extrapolation
2267
2268 Double_t pars[19][4] = {{0, 0, 0, 0},
2269 {0, 0, 0, 0},
2270 {8.28698e-06, 13.3515, 0, 0},
2271 {0.031189, 4.56933, 0.031189, 4.56933},
2272 {0.395468, 3.77954, 0.190848, 3.41027},
2273 {0.128467, 4.26652, 0.128467, 4.26652},
2274 {0.160979, 4.95095, 0.160979, 4.95095}};
2275
2276 Double_t pxz = vec.momxz;
2277 Int_t endsta = vec.stahit.rbegin()->first;
2278 Double_t p0 = pars[ista][0];
2279 Double_t p1 = pars[ista][1];
2280 if (endsta > ista + 2) {
2281 p0 = pars[ista][2];
2282 p1 = pars[ista][3];
2283 }
2284
2285 Double_t dx = p0 / TMath::Power(TMath::Log10(10 * pxz), p1);
2286 return dx;
2287}
2288
2289// -------------------------------------------------------------------------
2290
2291Double_t BmnStsVectorFinder::Proxim(Double_t phi0, Double_t phi)
2292{
2293 /// Adjust angle phi to be "around" phi0 - to avoid discontinuity around +- Pi
2294
2295 Double_t dPhi = phi0 - phi;
2296 if (TMath::Abs(dPhi) > TMath::Pi())
2297 phi += TMath::Pi() * 2 * TMath::Sign(1., dPhi);
2298 return phi;
2299}
2300
2301// -------------------------------------------------------------------------
2302// AZ-280722 Double_t BmnStsVectorFinder::LinearFit(CbmStsTrack *tr, Int_t newtr, Double_t& ty)
2303Double_t BmnStsVectorFinder::LinearFit(candvec& cand,
2304 candvec& cand2,
2305 CbmStsTrack* tr,
2306 Int_t newtr,
2307 Double_t& ty) // AZ-280722
2308{
2309 /// Linear LS fit of 3 points (in R-Y plane): y = a + b * l, l is a track projection length on X-Z plane
2310
2311 // static Double_t x[3], y[3], z[3], xs[3], ys[3], xys[3], x2s[3];
2312 static Double_t /*x[3], */ y[3], /*z[3], */ l[3] = {0}, lens[3] = {0}, ys[3], lys[3], l2s[3];
2313
2314 // AZ-150722 Int_t i0 = (newtr == 1) ? 0 : 2;
2315 Int_t i0 = (newtr == 0) ? 0 : 2;
2316
2317 for (Int_t i = i0; i < 3; ++i) {
2318 CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(tr->GetStsHitIndex(i));
2319 // AZ-280722 x[i] = hit->GetX();
2320 // AZ-280722 z[i] = hit->GetZ();
2321 if (i) {
2322 /*AZ-280722
2323 Double_t dx = x[i] - x[i-1];
2324 Double_t dz = z[i] - z[i-1];
2325 l[i] = lens[i] = TMath::Sqrt (dx * dx + dz * dz) + l[i-1];
2326 */
2327 if (i == 1)
2328 l[i] = lens[i] = cand.lengxz; // AZ-280722
2329 else
2330 l[i] = lens[i] = cand2.lengxz + l[i - 1]; // AZ-280722
2331 }
2332 y[i] = ys[i] = hit->GetY();
2333 // cout << " ily " << i << " " << l[i] << " " << y[i] << endl;
2334 }
2335
2336 for (Int_t i = i0; i < 3; ++i) {
2337 lys[i] = l[i] * y[i];
2338 l2s[i] = l[i] * l[i];
2339 if (i) {
2340 lens[i] += lens[i - 1];
2341 ys[i] += ys[i - 1];
2342 lys[i] += lys[i - 1];
2343 l2s[i] += l2s[i - 1];
2344 }
2345 }
2346 Double_t b = 3 * lys[2] - lens[2] * ys[2];
2347 b /= (3 * l2s[2] - lens[2] * lens[2]);
2348 Double_t a = (ys[2] - b * lens[2]) / 3;
2349
2350 Double_t chi2 = 0.0;
2351
2352 for (int i = 0; i < 3; ++i) {
2353 Double_t dy = a + b * l[i] - y[i];
2354 chi2 += (dy * dy);
2355 }
2356 // Check
2357 // TGraph gr(3,l,y);
2358 // gr.Fit("pol1");
2359 // cout << " LS fit: " << a << " " << b << " " << chi2 << endl;
2360 ty = b;
2361 return chi2;
2362}
2363
2364// -------------------------------------------------------------------------
2365
2366Double_t BmnStsVectorFinder::Curv3(candvec& cand1, candvec& cand2, candvec& cand3, int newtr)
2367{
2368 // Compute triplet curvature in ZX plane
2369
2370 static TVector3 points3[3], midPoint;
2371 int indx = 1;
2372 static Double_t by = 0.0;
2373 static int first = 99;
2374
2375 // if (first != fPass) {
2376 // first = fPass;
2377 if (first) {
2378 first = 0;
2379 points3[0].SetXYZ(fXyzv[0], fXyzv[1], fXyzv[2]);
2380 }
2381
2382 map<int, int>::iterator mit = cand3.stahit.begin();
2383 if (0) { // fPass > 1) { //0 - assuming tracks coming from (0,0,0)
2384 // if (newtr == 0 && cand3.stahit.size() > 2) {
2385 points3[0] = fmapHits[mit->second].xyz; // AZ-280722
2386 points3[0].SetY(0.0); // AZ-280722
2387 }
2388 // AZ-310722 ++mit;
2389 if (cand3.stahit.size() > 2)
2390 ++mit; // AZ-310722
2391
2392 for (; mit != cand3.stahit.end(); ++mit) {
2393 // if (newtr > 0 && indx < 2) { ++indx; continue; }
2394 points3[indx] = fmapHits[mit->second].xyz;
2395 if (indx == 1)
2396 midPoint = points3[indx];
2397 points3[indx] -= points3[0];
2398 points3[indx++].SetY(0.0);
2399 }
2400
2401 TVector3 vec21 = points3[2] - points3[1];
2402 // AZ-280722 Double_t cosAlpha = points3[2] * vec21 / points3[2].Mag() / vec21.Mag();
2403 // Double_t cosAlpha = points3[1] * vec21 / points3[1].Mag() / vec21.Mag(); //AZ-280722
2404 Double_t cosAlpha = points3[1] * vec21 / points3[1].Mag() / cand2.lengxz; // AZ-280722
2405 // AZ-280722 Double_t rad = points3[1].Mag() / 2. / TMath::Sin(TMath::ACos(cosAlpha));
2406 Double_t rad = points3[2].Mag() / 2. / TMath::Sin(TMath::ACos(cosAlpha)); // AZ-280722
2407 // Double_t bz = FairRunAna::Instance()->GetField()->GetBz(0.,0.,0.);
2408 // Double_t factor = 0.003 * bz / 10.; // 0.3 * 0.01 * 5kG / 10
2409 // AZ-310722 Double_t sign = TMath::Sign (1.0, points3[1].Cross(points3[2]).Y());
2410 Double_t sign = TMath::Sign(3e-4, points3[1].Cross(points3[2]).Y()); // AZ-310722
2411 if (newtr == 0) {
2412 FairField* magField = FairRunAna::Instance()->GetField();
2413 by = magField->GetBy(midPoint[0], midPoint[1], midPoint[2]);
2414 if (TMath::Abs(by) < 0.1)
2415 by = -10.0; // AZ-261222 - zero magnetic field
2416 }
2417 Double_t pxz = by * rad * sign;
2418 // cout << " pt " << pxz << endl;
2419 // AZ-170722 return pxz;
2420 return 1 / pxz; // AZ-170722
2421}
2422// -------------------------------------------------------------------------
2423
2424set<int> BmnStsVectorFinder::KalmanWindow(candvec& cand, int hitIndx)
2425{
2426 // Find hits around extrapolated track
2427
2428 CbmStsTrack track;
2429 string hitcode;
2430 MakeStsTrack(cand, hitcode, track);
2431 fitter.DoFit(&track);
2432 FairTrackParam param;
2433 CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(hitIndx);
2434 Double_t z = hit->GetZ();
2435 fitter.Extrapolate(&track, z, &param);
2436 Double_t sigx = 5 * TMath::Sqrt(param.GetCovariance(0, 0));
2437 Double_t sigy = 5 * TMath::Sqrt(param.GetCovariance(1, 1));
2438 if (fVerbose > 0)
2439 cout << " Kalman: " << z << " " << hit->GetStationNr() << " " << sigx << " " << sigy << endl;
2440
2441 // Get hits on the downstream station
2442 int istanext = hit->GetStationNr() - 1;
2443 multimap<Double_t, Int_t>::iterator mityb = fmapY[istanext].lower_bound(param.GetY() - sigy);
2444 multimap<Double_t, Int_t>::iterator mitye = fmapY[istanext].upper_bound(param.GetY() + sigy);
2445 multimap<Double_t, Int_t>::iterator mitxb = fmapX[istanext].lower_bound(param.GetX() - sigx);
2446 multimap<Double_t, Int_t>::iterator mitxe = fmapX[istanext].upper_bound(param.GetX() + sigx);
2447 // Get hits from the acceptance window
2448 set<Int_t> setX, setY, intersect;
2449 for (multimap<Double_t, Int_t>::iterator mit = mitxb; mit != mitxe; ++mit)
2450 if (fmapHits[mit->second].used == 0)
2451 setX.insert(mit->second);
2452 for (multimap<Double_t, Int_t>::iterator mit = mityb; mit != mitye; ++mit)
2453 if (fmapHits[mit->second].used == 0)
2454 setY.insert(mit->second);
2455 set_intersection(setX.begin(), setX.end(), setY.begin(), setY.end(), std::inserter(intersect, intersect.begin()));
2456 return intersect;
2457}
2458
2459// -------------------------------------------------------------------------
2460
2461std::string BmnStsVectorFinder::MakeCode(candvec& cand)
2462{
2463 // Make track code
2464
2465 string code("-");
2466 for (map<int, int>::iterator it = cand.stahit.begin(); it != cand.stahit.end(); ++it)
2467 code += (to_string(it->second) + "-");
2468 return code;
2469}
2470
2471// -------------------------------------------------------------------------
2472
2473void BmnStsVectorFinder::PrintHits(candvec& cand)
2474{
2475 // Print hits
2476
2477 for (map<int, int>::iterator mit = cand.stahit.begin(); mit != cand.stahit.end(); ++mit) {
2478 CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(mit->second);
2479 cout << mit->first << "(" << hit->GetZ() << "):" << mit->second << " - ";
2480 }
2481 cout << endl;
2482}
2483
2484// -------------------------------------------------------------------------
2485
2486void BmnStsVectorFinder::InitTMVA3()
2487{
2488 // Initialize TMVA package for triplets - AZ-300922
2489
2490 // This loads the library
2491 // TMVA::Tools::Instance();
2492
2493 // Default MVA methods to be trained + tested
2494 std::map<std::string, int> Use;
2495
2496 Use["MLPBNN"] = 1; // Recommended ANN with BFGS training method and bayesian regulator
2497 Use["BDTD"] = 1; // decorrelation + Adaptive Boost
2498
2499 // --- Create the Reader object
2500
2501 TMVA::Reader* reader = new TMVA::Reader("!Color:!Silent");
2502 fReaderTMVA3 = reader;
2503
2504 // Create a set of variables and declare them to the reader
2505 // The variable names MUST corresponds in name and type to those given in the weight file(s) used
2506
2507 // FIXME
2508 Float_t* var = fVarTMVA;
2509 int iii = 0;
2510 reader->AddVariable("xb", &var[iii++]);
2511 reader->AddVariable("yb", &var[iii++]);
2512 reader->AddVariable("zb", &var[iii++]);
2513 reader->AddVariable("tx", &var[iii++]);
2514 reader->AddVariable("ty", &var[iii++]);
2515 // reader->AddVariable( "dy", &var[iii++] );
2516 // reader->AddVariable( "dz", &var[iii++] );
2517 // reader->AddVariable( "dphix", &var[iii++] );
2518 // reader->AddVariable( "dphixsc", &var[iii++] );
2519 // reader->AddVariable( "dtany", &var[iii++] );
2520 reader->AddVariable("dtanysc1-dtanysc2", &var[iii++]);
2521 reader->AddVariable("lenxz1", &var[iii++]);
2522 reader->AddVariable("lenxz2", &var[iii++]);
2523 reader->AddVariable("1/pxz1", &var[iii++]);
2524 reader->AddVariable("(1/pxz1-1/pxz2)/(abs(1/pxz1)+abs(1/pxz2))*sqrt(lenxz1+lenxz2)", &var[iii++]);
2525 reader->AddVariable("pass", &var[iii++]);
2526
2527 // --- Book the MVA methods
2528
2529 // TString dir = "dataset3diff/weights/";
2530 TString dir = "dataset4diff/weights/";
2531 TString prefix = "TMVAClassification";
2532
2533 // Book method(s)
2534 for (std::map<std::string, int>::iterator it = Use.begin(); it != Use.end(); it++) {
2535 if (it->second) {
2536 TString methodName = TString(it->first) + TString(" method");
2537 TString weightfile = dir + prefix + TString("_") + TString(it->first) + TString(".weights.xml");
2538 reader->BookMVA(methodName, weightfile);
2539 }
2540 }
2541}
2542
2543// -------------------------------------------------------------------------
2544
2545Float_t BmnStsVectorFinder::TMVAOutput3(candvec& aaa1, candvec& aaa2)
2546{
2547 // Compute TMVA output for triplets
2548
2549 Float_t* var = fVarTMVA /*, out = 0.0*/;
2550 map<Int_t, Int_t>& hitMap = aaa1.stahit;
2551 CbmStsHit* hit = nullptr;
2552 int iii = 0, indx = 0;
2553
2554 /// hitMap ---- aaa.stahit
2555 for (map<Int_t, Int_t>::iterator mit = hitMap.begin(); mit != hitMap.end(); ++mit) {
2556 hit = (CbmStsHit*)fHitArray->UncheckedAt(mit->second);
2557 if (indx == 0) {
2558 var[iii++] = hit->GetX();
2559 var[iii++] = hit->GetY();
2560 var[iii++] = hit->GetZ();
2561 } else if (indx == 1) {
2562 Double_t dx = hit->GetX() - var[0];
2563 Double_t dy = hit->GetY() - var[1];
2564 Double_t dz = hit->GetZ() - var[2];
2565 var[iii++] = dx / dz;
2566 var[iii++] = dy / dz;
2567 var[iii++] = aaa1.ty - aaa2.ty;
2568 var[iii++] = aaa1.lengxz;
2569 var[iii++] = aaa2.lengxz;
2570 var[iii++] = aaa1.momxz;
2571 var[iii++] =
2572 (aaa1.momxz - aaa2.momxz) / (abs(aaa1.momxz) + abs(aaa2.momxz)) * sqrt(aaa1.lengxz + aaa2.lengxz);
2573 var[iii++] = (fPass == 2) ? 1 : fPass; // fPass==2 has the same triplet cuts as fPass==1
2574 }
2575 ++indx;
2576 }
2577 return fReaderTMVA3->EvaluateMVA("MLPBNN method");
2578}
2579
2580// -------------------------------------------------------------------------
2581
2582void BmnStsVectorFinder::InitTMVA2()
2583{
2584 // Initialize TMVA package for doublets - AZ-221022
2585
2586 // This loads the library
2587 // TMVA::Tools::Instance();
2588
2589 // Default MVA methods to be trained + tested
2590 std::map<std::string, int> Use;
2591
2592 Use["MLPBNN"] = 1; // Recommended ANN with BFGS training method and bayesian regulator
2593 Use["BDTD"] = 1; // decorrelation + Adaptive Boost
2594
2595 // --- Create the Reader object
2596
2597 TMVA::Reader* reader = new TMVA::Reader("!Color:!Silent");
2598 fReaderTMVA2 = reader;
2599
2600 // Create a set of variables and declare them to the reader
2601 // The variable names MUST corresponds in name and type to those given in the weight file(s) used
2602
2603 // FIXME
2604 Float_t* var = fVarTMVA;
2605 int iii = 0;
2606 reader->AddVariable("xb", &var[iii++]);
2607 reader->AddVariable("yb", &var[iii++]);
2608 reader->AddVariable("zb", &var[iii++]);
2609 reader->AddVariable("tx", &var[iii++]);
2610 reader->AddVariable("ty", &var[iii++]);
2611 // reader->AddVariable( "dphixsc", &var[iii++] );
2612 reader->AddVariable("dtanysc", &var[iii++]);
2613 reader->AddVariable("lenxz", &var[iii++]);
2614 reader->AddVariable("1/pxz", &var[iii++]);
2615 reader->AddVariable("corrq[1]-corrq[0]", &var[iii++]);
2616 reader->AddVariable("corrn[1]-corrn[0]", &var[iii++]);
2617 reader->AddVariable("pass", &var[iii++]);
2618
2619 // --- Book the MVA methods
2620
2621 TString dir = "dataset4diff/weights2/";
2622 TString prefix = "TMVAClassification";
2623
2624 // Book method(s)
2625 for (std::map<std::string, int>::iterator it = Use.begin(); it != Use.end(); it++) {
2626 if (it->second) {
2627 TString methodName = TString(it->first) + TString(" method");
2628 TString weightfile = dir + prefix + TString("_") + TString(it->first) + TString(".weights.xml");
2629 reader->BookMVA(methodName, weightfile);
2630 }
2631 }
2632}
2633
2634// -------------------------------------------------------------------------
2635
2636Float_t BmnStsVectorFinder::TMVAOutput2(candvec& aaa1)
2637{
2638 // Compute TMVA output for doublets
2639
2640 Float_t* var = fVarTMVA /*, out = 0.0*/;
2641 map<Int_t, Int_t>& hitMap = aaa1.stahit;
2642 CbmStsHit *hit = nullptr, *hit1 = nullptr;
2643 int iii = 0, indx = 0;
2644
2645 /// hitMap ---- aaa.stahit
2646 for (map<Int_t, Int_t>::iterator mit = hitMap.begin(); mit != hitMap.end(); ++mit) {
2647 hit = (CbmStsHit*)fHitArray->UncheckedAt(mit->second);
2648 if (indx == 0) {
2649 var[iii++] = hit->GetX();
2650 var[iii++] = hit->GetY();
2651 var[iii++] = hit->GetZ();
2652 hit1 = hit;
2653 } else if (indx == 1) {
2654 Double_t dx = hit->GetX() - var[0];
2655 Double_t dy = hit->GetY() - var[1];
2656 Double_t dz = hit->GetZ() - var[2];
2657 var[iii++] = dx / dz;
2658 var[iii++] = dy / dz;
2659 var[iii++] = aaa1.ty;
2660 var[iii++] = aaa1.lengxz;
2661 var[iii++] = aaa1.momxz;
2662 var[iii++] = hit->GetSignalDiv() - hit1->GetSignalDiv();
2663 var[iii++] = hit->GetTimeStamp() - hit1->GetTimeStamp();
2664 var[iii++] = (fPass == 2) ? 1 : fPass; // fPass==2 has the same triplet cuts as fPass==1
2665 }
2666 ++indx;
2667 }
2668 return fReaderTMVA2->EvaluateMVA("MLPBNN method");
2669}
kGEM
@ kGEM
Definition BmnDetectorList.h:15
workTime
static Double_t workTime
Definition BmnStsVectorFinder.cxx:41
i
int i
Definition P4_F32vec4.h:22
Bool_t
bool Bool_t
Int_t
int Int_t
Float_t
float Float_t
kFALSE
constexpr Bool_t kFALSE
Double_t
double Double_t
kTRUE
constexpr Bool_t kTRUE
Option_t
const char Option_t
Fatal
void Fatal(const char *location, const char *msgfmt,...)
Cross
TGLVector3 Cross(const TGLVector3 &v1, const TGLVector3 &v2)
BmnStsVectorFinder::PrintHits
void PrintHits(candvec &cand)
Definition BmnStsVectorFinder.cxx:2473
BmnStsVectorFinder::Proxim
Double_t Proxim(Double_t ang1, Double_t ang2)
Definition BmnStsVectorFinder.cxx:2291
BmnStsVectorFinder::GetHitId
std::set< Int_t > GetHitId(CbmStsHit *hit, Int_t &idmaxP)
Definition BmnStsVectorFinder.cxx:480
BmnStsVectorFinder::GetHitIdBmn
std::set< Int_t > GetHitIdBmn(CbmStsHit *hit, Int_t &idmaxP)
Definition BmnStsVectorFinder.cxx:557
BmnStsVectorFinder::fmapPhx
std::multimap< Double_t, Int_t > fmapPhx[19]
Definition BmnStsVectorFinder.h:148
BmnStsVectorFinder::fmapTy
std::multimap< Double_t, Int_t > fmapTy[19]
Definition BmnStsVectorFinder.h:148
BmnStsVectorFinder::fMap3Out
std::map< std::string, int > fMap3Out
Definition BmnStsVectorFinder.h:158
BmnStsVectorFinder::fTracks
std::multimap< Double_t, CbmStsTrack > fTracks
Definition BmnStsVectorFinder.h:169
BmnStsVectorFinder::Exec
virtual void Exec(Option_t *opt)
Virtual method Exec.
Definition BmnStsVectorFinder.cxx:280
BmnStsVectorFinder::FitTrack
Double_t FitTrack(candvec &cand)
Definition BmnStsVectorFinder.cxx:1358
BmnStsVectorFinder::fmapX
std::multimap< Double_t, Int_t > fmapX[19]
Definition BmnStsVectorFinder.h:148
BmnStsVectorFinder::fDigiMatches
TClonesArray * fDigiMatches
Definition BmnStsVectorFinder.h:137
BmnStsVectorFinder::Parabola
TVector3 Parabola(TVector3 &pos0, TVector3 &pos1, TVector3 &pos2)
Definition BmnStsVectorFinder.cxx:1563
BmnStsVectorFinder::fMap2
std::unordered_multimap< int, int > fMap2[19]
Definition BmnStsVectorFinder.h:161
BmnStsVectorFinder::fReaderTMVA3
TMVA::Reader * fReaderTMVA3
Definition BmnStsVectorFinder.h:185
BmnStsVectorFinder::fHit2id
std::unordered_multimap< Int_t, Int_t > fHit2id
Definition BmnStsVectorFinder.h:173
BmnStsVectorFinder::Finish
virtual void Finish()
Virtual method Finish.
Definition BmnStsVectorFinder.cxx:1556
BmnStsVectorFinder::Init
virtual InitStatus Init()
Virtual method Init.
Definition BmnStsVectorFinder.cxx:81
BmnStsVectorFinder::fCandCodes
std::set< string > fCandCodes[19]
Definition BmnStsVectorFinder.h:151
BmnStsVectorFinder::fCandVec2
std::vector< candvec > fCandVec2[19]
Definition BmnStsVectorFinder.h:155
BmnStsVectorFinder::MakeCode
std::string MakeCode(candvec &cand)
Definition BmnStsVectorFinder.cxx:2461
BmnStsVectorFinder::fSeedVec
std::vector< candvec > fSeedVec[19]
Definition BmnStsVectorFinder.h:149
BmnStsVectorFinder::fCandSet
std::set< string > fCandSet[19]
Definition BmnStsVectorFinder.h:178
BmnStsVectorFinder::TMVAOutput2
Float_t TMVAOutput2(candvec &aaa1)
Definition BmnStsVectorFinder.cxx:2636
BmnStsVectorFinder::KalmanWindow
std::set< int > KalmanWindow(candvec &cand, int hitIndx)
Definition BmnStsVectorFinder.cxx:2424
BmnStsVectorFinder::fmapHits
std::unordered_map< Int_t, hitinfo > fmapHits
Definition BmnStsVectorFinder.h:171
BmnStsVectorFinder::fMap2OutPtr
std::map< std::string, int > * fMap2OutPtr
Definition BmnStsVectorFinder.h:157
BmnStsVectorFinder::fTrackArray
TClonesArray * fTrackArray
Definition BmnStsVectorFinder.h:137
BmnStsVectorFinder::FinishEvent
virtual void FinishEvent()
Virtual method FinishEvent.
Definition BmnStsVectorFinder.cxx:1546
BmnStsVectorFinder::BmnStsVectorFinder
BmnStsVectorFinder()
Default constructor.
Definition BmnStsVectorFinder.cxx:44
BmnStsVectorFinder::fCandVec
std::vector< candvec > fCandVec[19]
Definition BmnStsVectorFinder.h:150
BmnStsVectorFinder::fStsPoints
TClonesArray * fStsPoints
Definition BmnStsVectorFinder.h:137
BmnStsVectorFinder::fMap3OutPtr
std::map< std::string, int > * fMap3OutPtr
Definition BmnStsVectorFinder.h:159
BmnStsVectorFinder::fReaderTMVA2
TMVA::Reader * fReaderTMVA2
Definition BmnStsVectorFinder.h:186
BmnStsVectorFinder::fmapY
std::multimap< Double_t, Int_t > fmapY[19]
Definition BmnStsVectorFinder.h:148
BmnStsVectorFinder::fCandVec3
std::vector< candvec > fCandVec3[19]
Definition BmnStsVectorFinder.h:163
BmnStsVectorFinder::fClusMaps
std::unordered_multimap< Int_t, Int_t > fClusMaps[2]
Definition BmnStsVectorFinder.h:173
BmnStsVectorFinder::fClusArray
TClonesArray * fClusArray[4]
Input / output array of CbmStsHits.
Definition BmnStsVectorFinder.h:137
BmnStsVectorFinder::~BmnStsVectorFinder
~BmnStsVectorFinder()
Destructor.
Definition BmnStsVectorFinder.cxx:71
BmnStsVectorFinder::fSilPoints
TClonesArray * fSilPoints
Definition BmnStsVectorFinder.h:137
BmnStsVectorFinder::LinearFit
Double_t LinearFit(candvec &cand, candvec &cand2, CbmStsTrack *tr, Int_t newtr, Double_t &ty)
Definition BmnStsVectorFinder.cxx:2303
BmnStsVectorFinder::ExtendTrack
void ExtendTrack(candvec cand)
Definition BmnStsVectorFinder.cxx:1142
BmnStsVectorFinder::fMap3
std::unordered_multimap< int, int > fMap3[19]
Definition BmnStsVectorFinder.h:165
BmnStsVectorFinder::ExtendTracks
void ExtendTracks(Int_t ista)
Definition BmnStsVectorFinder.cxx:1530
BmnStsVectorFinder::fMap2Out
std::map< std::string, int > fMap2Out
Definition BmnStsVectorFinder.h:156
BmnStsVectorFinder::fitter
BmnStsKFTrackFitter fitter
Definition BmnStsVectorFinder.h:177
BmnStsVectorFinder::FilterHit
Double_t FilterHit(candvec &cand, CbmStsTrack &track)
Definition BmnStsVectorFinder.cxx:1415
BmnStsVectorFinder::fMapCode3
std::unordered_map< std::string, int > fMapCode3[19]
Definition BmnStsVectorFinder.h:167
BmnStsVectorFinder::ExcludeHits
Int_t ExcludeHits(Int_t minHits, TClonesArray *trArray)
Definition BmnStsVectorFinder.cxx:394
BmnStsVectorFinder::fVectorArray
TClonesArray * fVectorArray
Definition BmnStsVectorFinder.h:138
BmnStsVectorFinder::MakeStsTrack
void MakeStsTrack(candvec &cand, std::string &hitcode, CbmStsTrack &track)
Definition BmnStsVectorFinder.cxx:1485
BmnStsVectorFinder::AreTracksDoubles
Bool_t AreTracksDoubles(CbmStsTrack &tr1, CbmStsTrack &tr2)
Definition BmnStsVectorFinder.cxx:1903
BmnStsVectorFinder::TMVAOutput3
Float_t TMVAOutput3(candvec &aaa1, candvec &aaa2)
Definition BmnStsVectorFinder.cxx:2545
BmnStsVectorFinder::Curv3
Double_t Curv3(candvec &cand1, candvec &cand2, candvec &cand3, int newtr3)
Definition BmnStsVectorFinder.cxx:2366
BmnStsVectorFinder::CheckVarx
Bool_t CheckVarx(Double_t dx, Double_t dz, Double_t tx, Double_t distxz, Double_t &varx)
Definition BmnStsVectorFinder.cxx:886
BmnStsVectorFinder::DxVsMom
Double_t DxVsMom(Int_t ista, candvec &aaa)
Definition BmnStsVectorFinder.cxx:2264
CbmKFHit
Abstract base class for hit is being used by the CbmKF Kalman Filter.
Definition CbmKFHit.h:17
CbmKFHit::Filter
virtual Int_t Filter(CbmKFTrackInterface &track, Bool_t downstream, Double_t &QP0)=0
CbmKFTrack::SetStsTrack
void SetStsTrack(CbmStsTrack &track, bool first=1)
Definition CbmKFTrack.cxx:40
CbmKFTrack::GetHit
CbmKFHit * GetHit(Int_t i)
Number of hits.
Definition CbmKFTrack.h:81
CbmKFTrack::GetTrackParam
void GetTrackParam(FairTrackParam &track)
Definition CbmKFTrack.cxx:48
CbmKFTrack::GetRefChi2
Double_t & GetRefChi2()
array[15] of covariance matrix
Definition CbmKFTrack.h:76
CbmKF::Instance
static CbmKF * Instance()
Definition CbmKF.h:34
CbmKF::GetNStsStations
int GetNStsStations() const
Definition CbmKF.h:87
CbmStsCluster
CbmStsCluster.h.
Definition CbmStsCluster.h:32
CbmStsCluster::GetNDigis
Int_t GetNDigis() const
Definition CbmStsCluster.h:66
CbmStsCluster::GetDigi
Int_t GetDigi(Int_t inum)
Definition CbmStsCluster.h:67
CbmStsDigiMatch
CbmStsDigiMatch.h.
Definition CbmStsDigiMatch.h:27
CbmStsDigiMatch::GetRefIndex
Int_t GetRefIndex(Int_t i=0) const
Accessor to StsPoint indizes.
Definition CbmStsDigiMatch.cxx:57
CbmStsDigiScheme::Instance
static CbmStsDigiScheme * Instance(int version=1)
Definition CbmStsDigiScheme.cxx:65
CbmStsDigiScheme::GetStationByNr
CbmStsStation * GetStationByNr(Int_t stationNr)
Station by index.
Definition CbmStsDigiScheme.cxx:193
CbmStsHit::GetStationNr
virtual Int_t GetStationNr() const
Definition CbmStsHit.h:67
CbmStsHit::GetSignalDiv
Double_t GetSignalDiv() const
Definition CbmStsHit.h:78
CbmStsHit::GetSystemId
Int_t GetSystemId() const
Accessors.
Definition CbmStsHit.h:66
CbmStsHit::GetDigi
Int_t GetDigi(Int_t iSide) const
Definition CbmStsHit.cxx:120
CbmStsSector::GetDx
Double_t GetDx() const
Definition CbmStsSector.h:83
CbmStsStation::GetSector
CbmStsSector * GetSector(Int_t iSector)
Definition CbmStsStation.h:65
CbmStsTrack::SetNDF
void SetNDF(Int_t ndf)
Definition CbmStsTrack.h:87
CbmStsTrack::GetParamLast
FairTrackParam * GetParamLast()
Definition CbmStsTrack.h:70
CbmStsTrack::SetParamLast
void SetParamLast(FairTrackParam &par)
Definition CbmStsTrack.h:84
CbmStsTrack::GetStsHits
TArrayI * GetStsHits()
Definition CbmStsTrack.h:71
CbmStsTrack::SetParamFirst
void SetParamFirst(FairTrackParam &par)
Definition CbmStsTrack.h:83
CbmStsTrack::GetFlag
Int_t GetFlag() const
Definition CbmStsTrack.h:65
CbmStsTrack::SetChi2
void SetChi2(Double_t chi2)
Definition CbmStsTrack.h:86
CbmStsTrack::SetFlag
void SetFlag(Int_t flag)
Definition CbmStsTrack.h:85
CbmStsTrack::GetNStsHits
Int_t GetNStsHits() const
Accessors.
Definition CbmStsTrack.h:60
CbmStsTrack::GetChi2
Double_t GetChi2() const
Definition CbmStsTrack.h:66
CbmStsTrack::GetStsHitIndex
Int_t GetStsHitIndex(Int_t iHit) const
Definition CbmStsTrack.h:62
CbmStsTrack::GetParamFirst
FairTrackParam * GetParamFirst()
Definition CbmStsTrack.h:69
CbmStsTrack::GetNDF
Int_t GetNDF() const
Definition CbmStsTrack.h:67
TArrayI::Set
void Set(Int_t n) override
TArray::GetSize
Int_t GetSize() const
TClonesArray::Compress
void Compress() override
TClonesArray::Remove
TObject * Remove(TObject *obj) override
TMVA::Reader::BookMVA
IMethod * BookMVA(const TString &methodTag, const TString &weightfile)
TMVA::Reader::AddVariable
void AddVariable(const TString &expression, Float_t *)
TObjArray::GetEntriesFast
Int_t GetEntriesFast() const
TObjArray::UncheckedAt
TObject * UncheckedAt(Int_t i) const
TStopwatch::RealTime
Double_t RealTime()
TStopwatch::Start
void Start(Bool_t reset=kTRUE)
TStopwatch::Stop
void Stop()
TVector3::SetY
void SetY(Double_t)
TVector3::SetXYZ
void SetXYZ(Double_t x, Double_t y, Double_t z)
TVector3::Mag
Double_t Mag() const
abs
RVec< PromoteType< T > > abs(const RVec< T > &v)
sqrt
RVec< PromoteType< T > > sqrt(const RVec< T > &v)
TMath::Min
Double_t Min(Double_t a, Double_t b)
TMath::ACos
Double_t ACos(Double_t)
TMath::Sign
Double_t Sign(Double_t a, Double_t b)
TMath::Exp
Double_t Exp(Double_t x)
TMath::ATan2
Double_t ATan2(Double_t y, Double_t x)
TMath::Power
Double_t Power(Double_t x, Double_t y)
TMath::Sqrt
Double_t Sqrt(Double_t x)
TMath::Pi
constexpr Double_t Pi()
TMath::Abs
Double_t Abs(Double_t d)
TMath::Sin
Double_t Sin(Double_t)
TMath::Max
Double_t Max(Double_t a, Double_t b)
TMath::Log10
Double_t Log10(Double_t x)
BmnStsVectorFinder::candvec::stahit
std::map< Int_t, Int_t > stahit
Definition BmnStsVectorFinder.h:62