BmnRoot
Loading...
Searching...
No Matches
BmnStsVectorFinderV9.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 "BmnStsVectorFinderV9.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"
20
21#include <TClonesArray.h>
22#include <TF1.h>
23#include <TGraph.h>
24#include <TMath.h>
25#include <TObjString.h>
26#include <iomanip>
27#include <iostream>
28#include <map>
29#include <set>
30#include <vector>
31
32using std::cout;
33using std::endl;
34using std::map;
35using std::multimap;
36using std::pair;
37using std::set;
38using std::vector;
39
40#include <TStopwatch.h>
41
42// std::ofstream outf("log.txt");
43
44static Double_t workTime = 0.0;
45
46// ----- Default constructor -------------------------------------------
47BmnStsVectorFinderV9::BmnStsVectorFinderV9()
48 : FairTask("STS Vector Finder V9")
49 , fHitArray(NULL)
50 , fTrackArray(NULL)
51 , fEvent(0)
52 , fPass(0)
53 , fNsta(0)
54 ,
55 // fNbranches(5),
56 fNbranches(7)
57 , // AZ-230625
58 // fNbranches(10), // AZ-270625
59 fExact(0)
60 , fExactSel(-9)
61 , discarded(0)
62 , fOut(0)
63 , // 0,2,3,6
64 fNhitsMin(nullptr)
65 , fMatBudgetFileName("")
66 , fUseTMVA(0)
67 , // 0,2,3 //AZ-300922
68 fTMVAtracks(1)
69 , // AZ-040525 - discriminate good from fake tracks
70 fNSi(4)
71 , // AZ-051022
72 fGemConfigFile("GemRun8.xml")
73 , fSilConfigFile("SiliconRun8.xml")
74{
75 for (int j = 0; j < 4; ++j)
76 fClusArray[j] = nullptr;
77 // fOutf = new ofstream("log.txt"); // AZ-290525 - for auxiliary output
78}
79// -------------------------------------------------------------------------
80
81// ----- Destructor ----------------------------------------------------
82BmnStsVectorFinderV9::~BmnStsVectorFinderV9() {}
83// -------------------------------------------------------------------------
84
85// ----- Public method Init --------------------------------------------
86InitStatus BmnStsVectorFinderV9::Init()
87{
88
89 // Get RootManager
90 FairRootManager* ioman = FairRootManager::Instance();
91 if (!ioman) {
92 cout << "-E- BmnStsVectorFinderV9::Init: " << "RootManager not instantiated!" << endl;
93 return kFATAL;
94 }
95
96 // Get input arrays
97 fClusArray[0] = (TClonesArray*)ioman->GetObject("StsCluster");
98 if (!fClusArray[0]) {
99 // cout << "-W- BmnStsVectorFinderV9::Init: "
100 // << "No StsCluster array!" << endl;
101 // return kERROR;
102 fClusArray[0] = (TClonesArray*)ioman->GetObject("BmnSiliconUpperCluster");
103 fClusArray[1] = (TClonesArray*)ioman->GetObject("BmnSiliconLowerCluster");
104 fClusArray[2] = (TClonesArray*)ioman->GetObject("BmnGemUpperCluster");
105 fClusArray[3] = (TClonesArray*)ioman->GetObject("BmnGemLowerCluster");
106 for (int j = 0; j < 4; ++j) {
107 // AZ-241222 if ( ! fClusArray[j] ) {
108 if (!fClusArray[0] && !fClusArray[2]) { // AZ-241222
109 cout << "-W- BmnStsVectorFinderV9::Init: " << "No BmnCluster array!" << j << endl;
110 return kERROR;
111 }
112 }
113 }
114
115 fHitArray = (TClonesArray*)ioman->GetObject("StsHit");
116 if (!fHitArray) {
117 // cout << "-W- BmnStsVectorFinderV9::Init: "
118 // << "No StsHit array!" << endl;
119 return kERROR;
120 }
121
122 fTrackArray = (TClonesArray*)ioman->GetObject("StsTrack");
123 if (!fTrackArray) {
124 // cout << "-W- BmnStsVectorFinderV9::Init: "
125 // << "No StsTrack array!" << endl;
126 // AZ-100722 return kERROR;
127 }
128
129 fDigiMatches = (TClonesArray*)ioman->GetObject("StsDigiMatch");
130 if (!fDigiMatches) {
131 // cout << "-W- BmnStsVectorFinderV9::Init: "
132 // << "No StsDigiMatch array!" << endl;
133 // return kERROR;
134 }
135
136 fStsPoints = (TClonesArray*)ioman->GetObject("StsPoint");
137 if (!fStsPoints) {
138 // cout << "-W- BmnStsVectorFinderV9::Init: "
139 // << "No StsPoint array!" << endl;
140 // AZ-181222 return kERROR;
141 }
142
143 fSilPoints = (TClonesArray*)ioman->GetObject("SiliconPoint");
144 if (!fSilPoints) {
145 // cout << "-E- BmnStsVectorFinderV9::Init: No SiliconPoint array!" << endl;
146 }
147
148 fVspPoints = (TClonesArray*)ioman->GetObject("VSPPoint");
149 if (!fVspPoints) {
150 // cout << "-E- BmnStsVectorFinder::Init: No SiliconPoint array!" << endl;
151 }
152
153 // Create and register output array
154 fVectorArray = new TClonesArray("CbmStsTrack");
155 ioman->Register("StsVector", "STS", fVectorArray, kTRUE);
156
157 if (fOut / 3 > 0) { // AZ-260922
158 // Store triplets
159 // AZ-240623 fMap3OutPtr = &fMap3Out;
160 // AZ-240623 ioman->RegisterAny("Triplets", fMap3OutPtr, kTRUE);
161 fStream3out.open("triplets.txt", std::ofstream::out);
162 }
163 if (fOut && fOut % 2 == 0) { // AZ-260922
164 // Store doublets
165 // AZ-240623 fMap2OutPtr = &fMap2Out;
166 // AZ-240623 ioman->RegisterAny("Doublets", fMap2OutPtr, kTRUE);
167 fStream2out.open("doublets.txt", std::ofstream::out);
168 }
169
170 // AZ-140226 fXyzv[0] = 0.4;
171 fXyzv[0] = 0.7; // AZ-140226
172 fXyzv[1] = 0.1;
173 fXyzv[2] = 0.0;
174
175 fitter.Init();
176
177 // Define logic
178 static Int_t nHitsMin[20] = {7, 5, 4, 4}; // min number of hits per track vs iteration - Run 7
179 // static Int_t nHitsMin[20] = {7, 5, 4, 3}; // min number of hits per track vs iteration - Run 7
180 // static Double_t dTanX[20] = {0.001, 0.003, 0.002, 0.008}; // window size in TanX vs iteration
181 static Double_t dTanX[20] = {0.002, 0.003, 0.002, 0.008}; // window size in TanX vs iteration
182 static Double_t dTanXTMVA[20] = {0.003, 0.0045, 0.003, 0.012}; // window size in TanX vs iteration - for TMVA
183 static Double_t dTanXB0[20] = {0.001, 0.002, 0.002, 0.002}; // window size in TanX vs iteration - for zero field
184 // static Double_t dTanY[20] = {0.01, 0.02, 0.01, 0.02}; // window size in TanY vs iteration
185 // AZ-230623 static Double_t dTanY[20] = {0.01, 0.02, 0.02, 0.03}; // window size in TanY vs iteration
186 static Double_t dTanY[20] = {0.02, 0.02, 0.02, 0.03}; // window size in TanY vs iteration
187 static Double_t dTanYTMVA[20] = {0.015, 0.03, 0.03, 0.045}; // window size in TanY vs iteration - for TMVA
188 static Double_t dTanYB0[20] = {0.01, 0.02, 0.02, 0.02}; // window size in TanY vs iteration - for zero field
189 // AZ-220623 static Double_t dTanY3[20] = {0.01, 0.02, 0.02, 0.03}; // window size in TanY vs iteration
190 static Double_t dTanY3[20] = {0.03, 0.02, 0.015, 0.015}; // AZ-220623 window size in TanY vs iteration
191 static Double_t dTanY3TMVA[20] = {0.015, 0.03, 0.03, 0.045}; // window size in TanY vs iteration - for TMVA
192 static Double_t dTanY3B0[20] = {0.01, 0.02, 0.02, 0.02}; // window size in TanY vs iteration - for zero field
193 static Double_t dTanXzero[20] = {0.01, 0.02, 0.02, 0.02}; // window size in TanX vs iteration - for zero field
194 // static Double_t dVarX[20] = {0.004, 0.01, 0.006, 0.02}; // window size in VarX vs iteration - coherent with
195 // ptCut AZ-120623 static Double_t ptCut[20] = {1.0, 0.3, 0.7, 0.1}; // min pxz
196 static Double_t ptCut[20] = {1.0, 0.7, 0.3, 0.1}; // min pxz - 120623
197 static Double_t ptCutB0[20] = {1., 1., 1., 1.}; // AZ-010123: min pxz - zero field
198 // static Double_t ptCutB0[20] = {0.5, 0.5, 0.5, 0.5 }; //AZ-010123: min pxz - zero field
199 static Double_t curvSta[20] = {
200 12.0, 12.0, 12.0, 7.5, 5.5, 4.5,
201 3.6, 3.0, 3.0, 3.0, 3.0}; // AZ-310722 - max curvature in stations (for tracks starting in Si)
202 static Double_t curvSta41[20] = {
203 12.0, 12.0, 12.0, 12.0, 7.5, 5.5, 4.5,
204 3.6, 3.0, 3.0, 3.0, 3.0}; // AZ-190625 - max curvature in stations (for tracks starting in Si) 4Si+VSP
205 // 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 -
206 // max curvature in stations (for tracks starting in Si) - 4 Si stations
207 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
208 // stations
209 // 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
210 // phi_XZ in stations - 4 Si
211 static Double_t phiXZ4[20] = {0.60, 0.63, 0.61, 0.62, 0.95, 0.74,
212 0.60, 0.51, 0.44, 0.39, 0.35}; // AZ-051122 - 4 Si 6_10_14_18
213 static Double_t phiXZ41[20] = {1.00, 0.60, 0.63, 0.61, 0.62, 0.95,
214 0.74, 0.60, 0.51, 0.44, 0.39, 0.35}; // AZ-190625 - 4 Si 6_10_14_18 + VSP
215 static Double_t tanXmax[20] = {0.95, 0.83, 0.83, 1.40, 0.91,
216 0.68, 0.56, 0.47, 0.41, 0.37}; // Max tan_X in stations
217 // 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
218 // tan_X in stations - 4 Si
219 static Double_t tanXmax4[20] = {0.68, 0.73, 0.69, 0.71, 1.40, 0.91,
220 0.68, 0.56, 0.47, 0.41, 0.37}; // AZ-051122 - 4 Si 6_10_14_18
221 static Double_t tanXmax41[20] = {1.50, 0.68, 0.73, 0.69, 0.71, 1.40,
222 0.91, 0.68, 0.56, 0.47, 0.41, 0.37}; // AZ-190625 - 4 Si 6_10_14_18 + VSP
223 static Double_t zMean[20] = {18.5, 27.2, 35.8, 60.9, 91.7,
224 123.4, 154.0, 185.8, 216.5, 248.0}; // Mean Z of stations
225 // 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
226 // 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,
227 // 185.8, 216.5, 248.0}; //AZ-051122 - 4Si 6_10_14_18
228 static Double_t zMean4[20] = {17.6, 26.4, 35.1, 43.7, 63.6, 94.2,
229 126.5, 156.8, 189.0, 219.2, 251.4}; // AZ-301222 - 4Si 6_10_14_18 - beam data
230 static Double_t zMean41[20] = {7.6, 17.6, 26.4, 35.1, 43.7, 63.6,
231 94.2, 126.5, 156.8, 189.0, 219.2, 251.4}; // AZ-190625 - 4Si 6_10_14_18 + VSP
232 //*/
233
234 FairField* magField = FairRunAna::Instance()->GetField();
235 if (magField == nullptr) {
236 LOG(error) << "BmnStsVectorFinderV9::Init: magField has null pointer";
237 return kERROR;
238 }
239 fBy = TMath::Abs(magField->GetBy(0.0, 0.0, 135.0)) / 10; // max. field in Tesla
240
241 fNhitsMin = nHitsMin;
242 fdTanX = (fUseTMVA) ? dTanXTMVA : dTanX;
243 fdTanY = (fUseTMVA) ? dTanYTMVA : dTanY;
244 fdTanY3 = (fUseTMVA) ? dTanY3TMVA : dTanY3;
245 fPTcut = ptCut;
246 // AZ-311222
247 if (fBy < 0.2) {
248 fdTanX = dTanXB0;
249 fdTanY = dTanYB0;
250 fdTanY3 = dTanY3B0;
251 fPTcut = ptCutB0;
252 } // AZ-311222
253 fdTanXB0 = dTanXzero; // AZ-010123
254 fPhiXZ = (fNSi == 4) ? phiXZ4 : phiXZ;
255 if (fNSi == 5)
256 fPhiXZ = phiXZ41; // AZ-190625 - 4Si+VSP
257 fTanXmax = (fNSi == 4) ? tanXmax4 : tanXmax;
258 if (fNSi == 5)
259 fTanXmax = tanXmax41; // AZ-190625 - 4Si+VSP
260 fZmean = (fNSi == 4) ? zMean4 : zMean;
261 if (fNSi == 5)
262 fZmean = zMean41; // AZ-190625 - 4Si+VSP
263 // fVarx = dVarX;
264 fCurvSta = curvSta;
265 if (fNSi == 5)
266 fCurvSta = curvSta41;
267
268 for (int j = 0; j < 20; ++j)
269 if (fPTcut[j] > 0.001)
270 fPTcut[j] = 1 / fPTcut[j]; // AZ-170722 - curvature
271
272 // fExact = 1; // for debug // removed DZ 13.10.2021
273 // fExactSel = 13; // for debug
274
275 if (fExact) {
276 // if (1) {
277 // Open windows
278 // Double_t scale = 3;
279 Double_t scale = 1;
280 for (Int_t j = 0; j < 20; ++j) {
281 dTanX[j] *= scale;
282 dTanY[j] *= scale;
283 // dX[j] *= scale;
284 }
285 }
286
287 // AZ-271222 for (int j = 0; j < 20; ++j) fCurvSta[j] *= (0.8 / fBy); //AZ-290822 - max curvature vs field
288 if (fBy > 0.1)
289 for (int j = 0; j < 20; ++j)
290 fCurvSta[j] *= (0.8 / fBy);
291 // AZ-311222 else fBy = 0.1; // zero mag. field
292
293 // Read system info
294 fTxBinw = 0.01;
295 TString fileGraphs = gSystem->ExpandPathName("$VMCWORKDIR");
296 fileGraphs += "/parameters/reco/";
297 if (fNSi == 4)
298 fileGraphs += "fitGraphs_8kG.txt"; // Run 8
299 else
300 fileGraphs += "fitGraphs_8kG_vsp.txt"; // AZ-220625 - Run 9
301 ReadTxInfo(fileGraphs);
302
303 if (fMatBudgetFileName != "")
304 fitter.ReadMatBudget(fMatBudgetFileName);
305
306 if (fUseTMVA == 2)
307 InitTMVA2(); // AZ-221022 - for doublets
308 if (fUseTMVA)
309 InitTMVA3(); // AZ-300922 - for triplets
310 if (fTMVAtracks)
311 InitTMVAtracks(); // AZ-040525 - for tracks
312
313 // Create tracking detector setups------------------------------------------------------
314 TString pathConfig = gSystem->Getenv("VMCWORKDIR");
315 fGemStationSet = new BmnGemStripStationSet(pathConfig + "/parameters/gem/XMLConfigs/" + fGemConfigFile);
316 fSilStationSet = new BmnSiliconStationSet(pathConfig + "/parameters/silicon/XMLConfigs/" + fSilConfigFile);
317
318 cout << "-I- BmnStsVectorFinderV9: Intialisation successfull" << endl;
319
320 return kSUCCESS;
321}
322// -------------------------------------------------------------------------
323
324void BmnStsVectorFinderV9::ReadTxInfo(TString fileName)
325{
326 // Read information about Tx windows
327
328 ifstream fin;
329 fin.open(fileName);
330 if (fin.fail()) {
331 cout << " !!! Error opening Tx file " << fileName << endl;
332 exit(1);
333 }
334 TString chline;
335 TObjArray* tokens = NULL;
336 int ista0 = -1, ista2;
337 vector<Float_t> params(10);
338
339 while (1) {
340 chline.ReadLine(fin);
341 if (fin.eof())
342 break;
343 // cout << chline << endl;
344 tokens = chline.Tokenize(" ");
345 int ntok = tokens->GetEntriesFast();
346 // if (chline.Contains("Station")) {
347 if (ntok == 4) {
348 // Station info
349 int j = 1;
350 ista0 = ((TObjString*)tokens->UncheckedAt(j++))->String().Atoi();
351 int ista1 = ((TObjString*)tokens->UncheckedAt(j++))->String().Atoi();
352 ista2 = 100 * ista0 + ista1;
353 if (fStatData.find(ista2) != fStatData.end())
354 continue;
355 // Float_t txmax = TMath::Abs (((TObjString*)tokens->UncheckedAt(j++))->String().Atof());
356 fTxStep = ((TObjString*)tokens->UncheckedAt(j))->String().Atof();
357 unordered_map<int, vector<Float_t>> mstat;
358 fStatData[ista2] = mstat;
359 // fTanXmax[ista0] = TMath::Abs(txmax); // average max Tx (for upper and lower parts)
360 fTanXmax[ista0] = -9;
361 } else {
362 // int ibin = ((TObjString*)tokens->UncheckedAt(0))->String().Atoi();
363 Float_t binc = ((TObjString*)tokens->UncheckedAt(0))->String().Atof();
364 int ilr = ((TObjString*)tokens->UncheckedAt(1))->String().Atoi();
365 Float_t binc1 = binc;
366 if (binc1 > 90.0)
367 binc1 -= 100.0;
368 if (fTanXmax[ista0] < -8) {
369 fTanXmax[ista0] = TMath::Abs(binc1 - fTxStep / 2);
370 // cout << ista0 << " " << fTanXmax[ista0] << endl;
371 }
372 int ibin = FindBin(ista0, fTxStep, binc1);
373 if (binc > 90.0)
374 ibin += 1000;
375 for (int j = 0; j < 5; ++j)
376 params[ilr * 5 + j] = ((TObjString*)tokens->UncheckedAt(j + 2))->String().Atof();
377 // if (ilr > 0) fStatData[ista2].insert(pair<int,vector<Float_t> >(ibin,params));
378 if (ilr > 0)
379 fStatData[ista2][ibin] = params;
380 }
381 if (tokens) {
382 tokens->Delete();
383 delete tokens;
384 tokens = NULL;
385 }
386 }
387 fin.close();
388 cout << "-I- BmnStsVectorFinderV9: Station info has been read from " << fileName << endl;
389}
390
391// ----- Public method Exec --------------------------------------------
392void BmnStsVectorFinderV9::Exec(Option_t* opt)
393{
394
395 TStopwatch sw;
396 sw.Start();
397
398 //*
399 // Reset output array
400 if (!fHitArray)
401 Fatal("Exec", "No StsHit array");
402
403 fVectorArray->Delete();
404 fMap2Out.clear(); // AZ-210922
405 fMap3Out.clear(); // AZ-260922
406
407 for (Int_t j = 0; j < 19; ++j) {
408 fCandSet[j].clear();
409 // fTripleCodes[j].clear();
410 fCandCodes[j].clear();
411 }
412
413 // AZ-100722 Int_t nTracks = fTrackArray->GetEntriesFast();
414 // AZ-100722 if (nTracks == 0) return;
415 Int_t nTracks = 0;
416 if (fTrackArray) {
417 nTracks = fTrackArray->GetEntriesFast();
418 if (nTracks == 0)
419 return;
420 }
421
422 const Int_t nPass = 4; // 4; //5; //8;
423 Int_t nsta = CbmKF::Instance()->GetNStsStations();
424 for (Int_t j = 0; j < nPass; ++j)
425 fNskips[j] = nsta - fNhitsMin[j];
426
427 Int_t nHits0 = fHitArray->GetEntriesFast(), idmaxP = 0;
428
429 // Fill cluster-to-hit maps
430 fHit2id.clear();
431 fClusMaps[0].clear();
432 fClusMaps[1].clear();
433 fmapHits.clear();
434
435 for (Int_t ihit = 0; ihit < nHits0; ++ihit) {
436 CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(ihit);
437 if (hit->GetZ() > 400.0)
438 continue; // AZ-290323 - exclude small GEM downstream
439 Int_t iclusF = hit->GetDigi(0);
440 Int_t iclusB = hit->GetDigi(1);
441 fClusMaps[0].insert(pair<Int_t, Int_t>(iclusF, ihit));
442 fClusMaps[1].insert(pair<Int_t, Int_t>(iclusB, ihit));
443 Int_t ista = hit->GetStationNr() - 0;
444 fNsta = TMath::Max(fNsta, ista);
445 // Code below for debugging
446 if (fVerbose > 0) {
447 set<Int_t> idset = GetHitId(hit, idmaxP);
448 for (set<Int_t>::iterator sit = idset.begin(); sit != idset.end(); ++sit)
449 fHit2id.insert(pair<Int_t, Int_t>(ihit, *sit));
450 }
451 }
452
453 Int_t nHitsOut = 0;
454 discarded = 0;
455 // AZ-080623 BuildTrackCand();
456 FillHitInfo(); // AZ-080623
457
458 for (Int_t ipass = 0; ipass < nPass; ++ipass) {
459 // for (Int_t ipass = nPass-1; ipass < nPass; ++ipass) {
460 fPass = ipass;
461 // Int_t minHits = 5; //4; // minimum number of hits on tracks to accept
462 // Int_t minHits = (ipass == 0) ? 15 : 5; // minimum number of hits on tracks to accept
463 Int_t minHits = (ipass == 0) ? 15 : -fNhitsMin[ipass - 1]; // minimum number of hits on tracks to accept
464 TClonesArray* trArray = (ipass == 0) ? fTrackArray : fVectorArray;
465 nHitsOut += ExcludeHits(minHits, trArray);
466
467 if (fVerbose > 0) {
468 cout << "-I- BmnStsVectorFinderV9: start - " << nHits0 << ", end - " << nHits0 - nHitsOut << endl;
469 }
470
471 Int_t ntr0 = fVectorArray->GetEntriesFast();
472 // BuildTrackCand();
473 /*
474 //ExtendTracks(fNsta-1);
475 FitTracks();
476 RemoveFakes();
477 RemoveDoubles();
478 cout << "\n ***** Pass " << ipass << ": Number of found tracks = "
479 << fVectorArray->GetEntriesFast() - ntr0 << " " << fVectorArray->GetEntriesFast() << "\n" << endl;
480 //if (fNhitsMin[ipass] <= 4) RemoveFakes();
481 */
482 BuildDoublets();
483 BuildTriplets();
484 fTracks.clear(); //
485 BuildTracks();
486 FitTracks();
487 RemoveDoubles();
488 if (fVerbose > 0) {
489 // if (1) {
490 cout << "\n ***** Pass " << ipass << ": Number of found tracks = " << fVectorArray->GetEntriesFast() - ntr0
491 << " " << fVectorArray->GetEntriesFast() << "\n"
492 << endl;
493 }
494 // GoToTarget(); //AZ-240625
495 }
496
497 if (fVerbose > 0)
498 cout << "discarded " << discarded << " track candidates" << endl;
499
500 // Post-processing - try to exclude fake tracks (with too many shared clusters)
501 ExcludeFakes();
502 if (fTMVAtracks)
503 TMVAOutputTra(); // AZ-040525
504
505 sw.Stop();
506 workTime += sw.RealTime();
507}
508
509// -------------------------------------------------------------------------
510
511Int_t BmnStsVectorFinderV9::ExcludeHits(Int_t minHits, TClonesArray* trArray)
512{
513 // Exclude hits used for tracking
514
515 // AZ-100722 Int_t nTracks = trArray->GetEntriesFast();
516 Int_t nTracks = 0; // AZ-100722
517 if (trArray)
518 nTracks = trArray->GetEntriesFast(); // AZ-100722
519
520 if (nTracks == 0)
521 return 0;
522 Int_t nHitsOut = 0;
523 // AZ-130722 multimap<Int_t,Int_t>::iterator mit;
524 // AZ-130722 pair<multimap<Int_t,Int_t>::iterator,multimap<Int_t,Int_t>::iterator> ret;
525 unordered_multimap<Int_t, Int_t>::iterator mit;
526 pair<unordered_multimap<Int_t, Int_t>::iterator, unordered_multimap<Int_t, Int_t>::iterator> ret; // AZ-130722
527
528 for (Int_t itra = 0; itra < nTracks; ++itra) {
529 CbmStsTrack* track = (CbmStsTrack*)trArray->UncheckedAt(itra);
530 if (track->GetNStsHits() < minHits) {
531 trArray->Remove(track);
532 continue;
533 }
534
535 Int_t nHitsTr = track->GetNStsHits();
536
537 for (Int_t ihit = 0; ihit < nHitsTr; ++ihit) {
538 int indx = track->GetStsHitIndex(ihit); // AZ-130722
539 // AZ-130722 CbmStsHit *hit = (CbmStsHit*) fHitArray->UncheckedAt(track->GetStsHitIndex(ihit));
540 CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(indx); // AZ-130722
541 if (!hit)
542 continue;
543
544 if (fNhitsMin[fPass] < -5) {
545 // if (fNhitsMin[fPass] > 6) { //AZ-180722
546 // Exclude all hits created from given clusters
547 Int_t iclusF = hit->GetDigi(0);
548 ret = fClusMaps[0].equal_range(iclusF);
549 for (mit = ret.first; mit != ret.second; ++mit) {
550 CbmStsHit* hit1 = (CbmStsHit*)fHitArray->UncheckedAt(mit->second);
551 if (hit1->GetUniqueID() == 0)
552 ++nHitsOut; // excluded hit - used in track
553 hit1->SetUniqueID(1); // exclude hit - used in track
554 // AZ-130722 fmapHits[ihit].used = 1;
555 fmapHits[mit->second].used = 1; // AZ-130722
556 }
557 Int_t iclusB = hit->GetDigi(1);
558 ret = fClusMaps[1].equal_range(iclusB);
559 for (mit = ret.first; mit != ret.second; ++mit) {
560 CbmStsHit* hit1 = (CbmStsHit*)fHitArray->UncheckedAt(mit->second);
561 if (hit1->GetUniqueID() == 0)
562 ++nHitsOut; // excluded hit - used in track
563 hit1->SetUniqueID(1); // exclude hit - used in track
564 // AZ-130722 fmapHits[ihit].used = 1;
565 fmapHits[mit->second].used = 1; // AZ-130722
566 }
567 }
568 if (hit->GetUniqueID() == 0)
569 ++nHitsOut;
570 hit->SetUniqueID(1);
571 // AZ-130722 fmapHits[ihit].used = 1;
572 fmapHits[indx].used = 1; // AZ-130722
573
574 // AZ-130722 Remove used doublets
575 //*
576 int ista = hit->GetStationNr() - 1;
577 ret = fMap2[ista].equal_range(indx);
578 fMap2[ista].erase(ret.first, ret.second);
579 // AZ-140722 Remove used triplets
580 ret = fMap3[ista].equal_range(indx);
581 for (mit = ret.first; mit != ret.second; ++mit) {
582 candvec cand = fCandVec3[ista][mit->second];
583 string code = MakeCode(cand);
584 fMapCode3[ista].erase(code);
585 }
586 fMap3[ista].erase(ret.first, ret.second);
587 //*/
588 }
589 }
590
591 trArray->Compress();
592 return nHitsOut;
593}
594
595// -------------------------------------------------------------------------
596
597set<Int_t> BmnStsVectorFinderV9::GetHitId(CbmStsHit* hit, Int_t& idmaxP)
598{
599 // Get IDs contributing to given hit (for debugging purposes)
600
601 // AZ-241222 if (fClusArray[1]) return GetHitIdBmn (hit, idmaxP); // BmnRoot data structures
602 if (1)
603 return GetHitIdBmn(hit, idmaxP); // BmnRoot data structures - AZ-241222
604
605 set<Int_t> ids;
606 // return ids; // FIXME
607
608 Int_t iclusF = hit->GetDigi(0);
609 Int_t iclusB = hit->GetDigi(1);
610 CbmStsCluster* clusF = (CbmStsCluster*)fClusArray[0]->UncheckedAt(iclusF);
611 CbmStsCluster* clusB = (CbmStsCluster*)fClusArray[0]->UncheckedAt(iclusB);
612
613 map<Int_t, Double_t> indF, indB;
614 // set<Int_t> ids;
615 Int_t nDigis = clusF->GetNDigis();
616 Double_t pmax = 0.0;
617 TVector3 mom3;
618
619 for (Int_t j = 0; j < nDigis; ++j) {
620 CbmStsDigiMatch* digiMatch = (CbmStsDigiMatch*)fDigiMatches->UncheckedAt(clusF->GetDigi(j));
621 for (Int_t j1 = 0; j1 < 3; ++j1) {
622 Int_t ip = digiMatch->GetRefIndex(j1);
623 if (ip < 0)
624 break;
625 CbmStsPoint* p = (CbmStsPoint*)fStsPoints->UncheckedAt(ip);
626 p->Momentum(mom3);
627 indF[p->GetTrackID()] = mom3.Mag();
628 }
629 }
630
631 nDigis = clusB->GetNDigis();
632
633 for (Int_t j = 0; j < nDigis; ++j) {
634 CbmStsDigiMatch* digiMatch = (CbmStsDigiMatch*)fDigiMatches->UncheckedAt(clusB->GetDigi(j));
635 for (Int_t j1 = 0; j1 < 3; ++j1) {
636 Int_t ip = digiMatch->GetRefIndex(j1);
637 if (ip < 0)
638 break;
639 CbmStsPoint* p = (CbmStsPoint*)fStsPoints->UncheckedAt(ip);
640 // Not needed p->Momentum(mom3);
641 indB[p->GetTrackID()] = 1.0; // mom3.Mag();
642 }
643 }
644
645 idmaxP = -1;
646
647 for (map<Int_t, Double_t>::iterator mit = indF.begin(); mit != indF.end(); ++mit) {
648 Int_t id = (indB.count(mit->first) == 0) ? -mit->first : mit->first;
649 ids.insert(id); // negative ID for ghost crossings
650 if (id >= 0 && mit->second > pmax) {
651 pmax = mit->second;
652 idmaxP = id;
653 }
654 }
655 // This is to be consistent with CBM hit matching procedure
656 // (hit is matched with the closest point)
657 if (hit->GetRefIndex() >= 0) {
658 CbmStsPoint* p = (CbmStsPoint*)fStsPoints->UncheckedAt(hit->GetRefIndex());
659 idmaxP = p->GetTrackID();
660 } else
661 idmaxP = -1;
662
663 for (map<Int_t, Double_t>::iterator mit = indB.begin(); mit != indB.end(); ++mit) {
664 if (ids.count(mit->first) > 0)
665 continue;
666 ids.insert(-mit->first); // negative ID for ghost crossings
667 }
668
669 return ids;
670}
671
672// -------------------------------------------------------------------------
673
674set<Int_t> BmnStsVectorFinderV9::GetHitIdBmn(CbmStsHit* hit, Int_t& idmaxP)
675{
676 // Get IDs contributing to given hit (for debugging purposes) - for BmnRoot data structures
677
678 set<Int_t> ids;
679 // int ind = 0; // Si hit
680 // if (hit->GetSystemId() == kGEM) ind = 2; // GEM hit
681
682 if (hit->GetRefIndex() == -1 || fStsPoints == nullptr || fSilPoints == nullptr)
683 ids.insert(-1); // fake hit
684 else {
685 // AZ-150722 CbmStsPoint *p = (CbmStsPoint*) fStsPoints->UncheckedAt(hit->GetRefIndex());
686 FairMCPoint* p = nullptr;
687 if (hit->GetSystemId() == kGEM)
688 p = (FairMCPoint*)fStsPoints->UncheckedAt(hit->GetRefIndex());
689 else if (hit->GetSystemId() == kVSP)
690 p = (FairMCPoint*)fVspPoints->UncheckedAt(hit->GetRefIndex());
691 else
692 p = (FairMCPoint*)fSilPoints->UncheckedAt(hit->GetRefIndex());
693 ids.insert(p->GetTrackID());
694 }
695 idmaxP = *ids.begin(); // AZ-290922
696 return ids;
697}
698
699// -------------------------------------------------------------------------
700
701// AZ-080623 void BmnStsVectorFinder::BuildTrackCand()
702void BmnStsVectorFinderV9::FillHitInfo() // AZ-080623
703{
704 // Fill hit information
705
706 for (Int_t ist = 0; ist < fNsta; ++ist) {
707 // AZ-080623 fmapPhx[ist].clear();
708 fmapTx[ist].clear();
709 fmapTy[ist].clear();
710 fmapX[ist].clear();
711 fmapY[ist].clear();
712 fSeedVec[ist].clear();
713 fCandVec[ist].clear();
714 // fCandMap2[ist].clear();
715 // fCandMap3[ist].clear();
716 fCandVec2[ist].clear();
717 fMap2[ist].clear();
718 fCandVec3[ist].clear();
719 fMap3[ist].clear();
720 fMapCode3[ist].clear();
721 fTxBins[ist].clear();
722 }
723 // fmapHits.clear();
724
725 // AZ-220725 fTxBinw = 0.01; // bin width for hits
726 Int_t nHits = fHitArray->GetEntriesFast(), idmaxP = 0;
727 if (fVerbose > 0)
728 cout << "nHits " << nHits << endl;
729 TVector3 pos;
730
731 for (Int_t ih = 0; ih < nHits; ++ih) {
732 CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(ih);
733 if (hit->GetUniqueID())
734 continue; // used hit
735 if (hit->GetZ() > 400.0)
736 continue; // AZ-290323 - exclude small GEM downstream
737 // if (hit->GetSystemId() == kGEM) continue; //!!!AZ-100123 - exclude GEMs
738 // if (hit->GetStationNr() > 6) continue; // AZ-220226 - keep only 1'st GEM
739 // if (hit->GetSystemId() != kGEM) continue; //!!!AZ-100123 - exclude Si
740 // else if (hit->GetSystemId() != kGEM && hit->GetStationNr() == 1) continue; //!!!AZ-110123 - exclude station
741 // 1 else if (hit->GetSystemId() != kGEM && hit->GetStationNr() == 4) continue; //!!!AZ-110123 - exclude
742 // station 1 if (hit->GetSystemId() == kGEM && hit->GetStationNr() > 7 && TMath::Abs(hit->GetX()) < 10.0 &&
743 // TMath::Abs(hit->GetY()) < 15.0) continue; //!!!AZ-130123 - exclude beam region
744 // if (hit->GetSystemId() == kGEM && hit->GetStationNr() > 10) continue; // AZ-190226 - exclude 2 last
745 // stations
746 /*if (hit->GetSystemId() == kGEM && hit->GetSectorNr() > 4) continue; //300123 - exclude lower GEMs
747 else if (hit->GetSystemId() != kGEM) {
748 //AZ-300123 - exclude lower Si
749 if (hit->GetStationNr() == 1 && hit->GetSectorNr() > 3) continue;
750 if (hit->GetStationNr() == 2 && hit->GetSectorNr() > 5) continue;
751 if (hit->GetStationNr() == 3 && hit->GetSectorNr() > 7) continue;
752 if (hit->GetStationNr() == 4 && hit->GetSectorNr() > 9) continue;
753 }*/
754 // AZ-260423 - apply cluster charge correlation cut
755 // if (hit->GetSystemId() == kGEM && TMath::Abs(hit->GetSignalDiv()+0.08) > 0.5) continue;
756 // else if (hit->GetSystemId() != kGEM && TMath::Abs(hit->GetSignalDiv()-0.11) > 0.5) continue;
757 if (fExactSel >= 0) {
758 // Accept only hits with given ID - for debug
759 set<Int_t> ids = GetHitId(hit, idmaxP);
760 if (ids.count(fExactSel) == 0)
761 continue;
762 }
763 if (fExact) {
764 // Accept only true hits
765 set<Int_t> ids = GetHitId(hit, idmaxP);
766 if (idmaxP < 0)
767 continue;
768 }
769 Int_t ista = hit->GetStationNr() - 1;
771 // if (ista > 2) continue;
772 // if (ista >= 4) continue; //AZ-201222 - Si only
774 hit->Position(pos);
775 Double_t dx = hit->GetX() - fXyzv[0];
776 Double_t dy = hit->GetY() - fXyzv[1];
777 Double_t dz = hit->GetZ() - fXyzv[2];
778 Double_t by = (fBy < 0.2) ? 0.8 : fBy; // AZ-271222
779 // AZ-250822 if (fPass >= 0) fmapHits[ih] = hitinfo(pos,dx/dz,dy/dz);
780 // 271222 if (fPass >= 0) fmapHits[ih] = hitinfo (pos, TMath::ATan2(dx,dz)/fPhiXZ[ista]/fZmean[ista]/fBy,
781 // dy/dz/fTanXmax[ista]); //AZ-280822 AZ-010123 if (fPass >= 0) fmapHits[ih] = hitinfo (pos,
782 // TMath::ATan2(dx,dz)/fPhiXZ[ista]/fZmean[ista]/by, dy/dz/fTanXmax[ista]); //AZ-271222
783 Double_t tanx = dx / dz;
784 // Double_t tanx = dx / (fZmean[ista] - fXyzv[2]); //AZ-090623 - w.r.t. nominal vertex (mean station Z!!!)
785 if (fPass >= 0)
786 fmapHits[ih] = hitinfo(pos, TMath::ATan2(dx, dz) / fPhiXZ[ista] / fZmean[ista] / by,
787 // AZ-070623 dx/dz/fTanXmax[ista], dy/dz/fTanXmax[ista]); //AZ-010123
788 // tanx/fTanXmax[ista], dy/dz/fTanXmax[ista]); //AZ-070623
789 // AZ-230623 tanx, dy/dz/fTanXmax[ista]); //AZ-080623
790 tanx, dy / dz); // AZ-230623
791 fmapX[ista].insert(pair<Double_t, Int_t>(pos[0], ih));
792 fmapY[ista].insert(pair<Double_t, Int_t>(pos[1], ih));
793 // AZ-080623 fmapPhx[ista].insert(pair<Double_t,Int_t>(fmapHits[ih].phx,ih));
794 fmapTx[ista].insert(pair<Double_t, Int_t>(fmapHits[ih].tx, ih));
795 fmapTy[ista].insert(pair<Double_t, Int_t>(fmapHits[ih].ty, ih));
796 // AZ-070623 - hit binning in Tx
797 int ista2 = ista * 100 + ista + 1;
798 int bin = FindBin(ista, fTxBinw, fmapHits[ih].tx);
799 if (bin < 0)
800 continue; // AZ-270623
801 int binStat = FindBin(ista, fTxStep, BinCenterTx(ista, bin, fTxBinw));
802 if (dy > 0)
803 binStat += 1000; // AZ-270623
804 if (fStatData[ista2].find(binStat) == fStatData[ista2].end())
805 continue;
806 // int bin = (tanx + fTanXmax[ista]) / fTxBinw;
807 if (dy > 0)
808 bin += 1000;
809 // int bin = (tanxMean + fTanXmax[ista]) / fTxBinw;
810 fTxBins[ista].insert(pair<int, int>(bin, ih));
811 } // for (Int_t ih = 0;
812 /*
813 //const Int_t stastop[9] = { fNsta-1, fNsta-2, fNsta-4}; // stop station vs iteration
814 //const Int_t stastop[9] = { fNsta-3, fNsta-2, fNsta-3}; // stop station vs iteration
815 Int_t stastop = fNsta - 1 - fNskips[fPass];
816 stastop = 0;
817
818 if (fVerbose > 0) cout << "fNskips fPass " << fNskips[fPass] << " " << fPass << " stastop " << stastop << " fNsta-1
819 " << fNsta-1 << endl;
820
821 for (Int_t ista = fNsta-1; ista >= stastop; --ista) {
822 if (fVerbose > 0) cout << "xmapsize " << fmapX[ista].size() << endl;
823 for (multimap<Double_t,Int_t>::iterator mit = fmapX[ista].begin(); mit != fmapX[ista].end(); ++mit) {
824 candvec aaa;
825 //AZ-280323 aaa.nskips = fNsta - 1 - ista;
826 aaa.nskips = 0; //AZ-280323
827 aaa.momxz = 0.0;
828 aaa.stahit[ista] = mit->second;
829 aaa.code = "-" + to_string(aaa.stahit[ista]) + "-"; // hit index coded
830 //if (fExact) {
831
832 set<Int_t> ids = GetHitId(mit->second, idmaxP);
833 aaa.idmaxP = idmaxP;
834 //}
835 //cout << "aaa.z " <<((CbmStsHit*) fHitArray->UncheckedAt(aaa.stahit[ista]))->GetZ() << " ";
836 fSeedVec[ista].push_back(aaa);
837 } //for (multimap<Double_t,Int_t>::iterator mit
838
839 if (fVerbose > 0) {
840 Int_t ncand = fSeedVec[ista].size();
841 cout << " Vector stat: " << ista << " " << ncand;// << endl;
842 //AZ-130722 pair<multimap<Int_t, Int_t>::iterator, multimap<Int_t, Int_t>::iterator> ret;
843 pair<unordered_multimap<Int_t, Int_t>::iterator, unordered_multimap<Int_t, Int_t>::iterator> ret; //AZ-130722
844
845 for (Int_t j = 0; j < ncand; ++j) {
846 Int_t ih = fSeedVec[ista][j].stahit[ista];
847 ret = fHit2id.equal_range(ih);
848 cout << " (" << ih << "*";
849 //AZ-130722 for (multimap<Int_t, Int_t>::iterator mit = ret.first; mit != ret.second; ++mit) {
850 for (unordered_multimap<Int_t, Int_t>::iterator mit = ret.first; mit != ret.second; ++mit) { //AZ-130722
851 if (mit != ret.first) cout << ":";
852 cout << mit->second;
853 }
854 cout << ")";
855 }
856 cout << "\n";
857 }
858 } // for (Int_t ista = fNsta-1;
859 */
860}
861
862// -------------------------------------------------------------------------
863
864void BmnStsVectorFinderV9::BuildDoublets()
865{
866 // Build doublets
867
868 int rebuild2 = 1; // AZ-140722
869
870 // Run8
871 const unordered_map<string, float> scaleMap4 = {
872 {"01", 1.9}, {"02", 2.8}, {"12", 1.0}, {"13", 1.5}, {"23", 0.65}, {"24", 1.4}, {"34", 0.95},
873 {"35", 1.5}, {"45", 1.0}, {"46", 1.5}, {"56", 0.7}, {"57", 1.0}, {"67", 0.6}, {"68", 0.7},
874 {"78", 0.35}, {"79", 0.5}, {"89", 0.3}, {"810", 0.35}, {"910", 0.25}};
875 // Run 9
876 const unordered_map<string, float> scaleMap5 = {
877 {"01", 11.5}, {"02", 15.0}, {"12", 1.9}, {"13", 2.8}, {"23", 1.0}, {"24", 1.5}, {"34", 0.65},
878 {"35", 1.4}, {"45", 0.95}, {"46", 1.5}, {"56", 1.0}, {"57", 1.5}, {"67", 0.7}, {"68", 1.0},
879 {"78", 0.6}, {"79", 0.7}, {"89", 0.35}, {"810", 0.5}, {"910", 0.3}, {"911", 0.35}, {"1011", 0.25}};
880
881 for (Int_t ist = 0; ist < fNsta; ++ist) {
882 fCandVec[ist].clear();
883 fCandVec2[ist].clear();
884 fMap2[ist].clear();
885 fCandVec3[ist].clear();
886 fMap3[ist].clear();
887 fMapCode3[ist].clear();
888 // AZ-140722 - Does not seem to speed up
889 /*
890 if ( (fPass == 0) || (TMath::Abs(fdTanX[fPass]-fdTanX[fPass-1]) > 0.001) || //AZ-140722
891 (TMath::Abs(fdTanY[fPass]-fdTanY[fPass-1]) > 0.001) ) {
892 fCandVec2[ist].clear();
893 fMap2[ist].clear();
894 rebuild2 = 1;
895 }
896 */
897 }
898 if (rebuild2 == 0)
899 return; // AZ-140722 - do not rebuild doublets
900 const unordered_map<string, float>& scaleMap = (fNSi == 4) ? scaleMap4 : scaleMap5;
901
902 Int_t idmaxP = 0;
903
904 for (Int_t ista = 0; ista < fNsta - 1; ++ista) {
905 int nHits = fTxBins[ista].size();
906 // cout << " nHits: " << ista << " " << nHits << endl;
907 if (nHits == 0)
908 continue;
909 string ssta0 = to_string(ista);
910
911 for (int istanext = ista + 1; istanext < ista + 3; ++istanext) {
912 // for (int istanext = ista + 1; istanext < ista + 2; ++istanext) { // no jumps over station
913 if (istanext > fNsta - 1)
914 continue;
915 if (fmapX[istanext].size() == 0)
916 continue;
917 Double_t dty = fdTanY[fPass]; //, dphx = fdTanX[fPass];
918 // dty *= 5; //AZ-040822
919 // dtx *= 2; //AZ-040822
920 int nskips = istanext - (ista + 1);
921 int ibin0 = -99, ista2 = ista * 100 + istanext;
922 set<Int_t> setTx;
923 Float_t tx = 0;
924 string ssta2 = ssta0 + to_string(istanext);
925 dty *= (scaleMap.find(ssta2)->second);
926
927 for (unordered_multimap<int, int>::iterator mit = fTxBins[ista].begin(); mit != fTxBins[ista].end(); ++mit)
928 {
929 // AZ-070623 - Loop over binned hits
930 int ih = mit->second;
931 if (fmapHits[ih].used)
932 continue; // used hit
933 candvec aaa;
934 aaa.nskips = 0;
935 aaa.momxz = 0.0;
936 aaa.stahit[ista] = mit->second;
937 aaa.code = "-" + to_string(aaa.stahit[ista]) + "-"; // hit index coded
938 set<Int_t> ids = GetHitId(mit->second, idmaxP);
939 aaa.idmaxP = idmaxP;
940
941 int ibin = mit->first;
942 if (ibin != ibin0) {
943 // Compute dTx constraints for this bin
944 ibin0 = ibin;
945 setTx.clear();
946 Float_t dtx[2] = {0};
947 tx = BinCenterTx(ista, ibin, fTxBinw);
948 if (fBy < 0.2) {
949 // AZ-140226 - for B = 0
950 dtx[0] = -fdTanXB0[fPass] * 2;
951 dtx[1] = fdTanXB0[fPass] * 2;
952 } else {
953 DTxWindow(ista2, tx, ibin / 1000, fPTcut[fPass], dtx);
954 // dtx[0] *= 2; dtx[1] *= 2; //
955 if (istanext < fNSi) {
956 dtx[0] *= 2;
957 dtx[1] *= 2;
958 } // AZ-120623
959 }
960 // else { dtx[0] *= 1.5; dtx[1] *= 1.5; }
961 // else { dtx[0] *= 1.5; dtx[1] *= 1.5; } // AZ-260625
962 // Get hits on the downstream station
963 // multimap<Double_t,Int_t>::iterator mitxb = fmapPhx[istanext].lower_bound(phx-dphx);
964 // multimap<Double_t,Int_t>::iterator mitxe = fmapPhx[istanext].upper_bound(phx+dphx);
965 multimap<Double_t, Int_t>::iterator mitxb = fmapTx[istanext].lower_bound(tx + dtx[0]);
966 multimap<Double_t, Int_t>::iterator mitxe = fmapTx[istanext].upper_bound(tx + dtx[1]);
967 for (multimap<Double_t, Int_t>::iterator mit1 = mitxb; mit1 != mitxe; ++mit1)
968 if (fmapHits[mit1->second].used == 0)
969 setTx.insert(mit1->second);
970 }
971
972 Double_t ty = fmapHits[ih].ty;
973 // AZ-230623 ty *= (fTanXmax[ista] / fTanXmax[istanext]); //AZ-280822
974 // Double_t tx = fmapHits[ih].tx;
975 // tx *= (fTanXmax[ista] / fTanXmax[istanext]); //AZ-010123
976 multimap<Double_t, candvec> mapDoublets; // AZ-251022 TMVA output -> doublets
977
978 // Get hits on the downstream station
979 multimap<Double_t, Int_t>::iterator mityb = fmapTy[istanext].lower_bound(ty - dty);
980 multimap<Double_t, Int_t>::iterator mitye = fmapTy[istanext].upper_bound(ty + dty);
981 // Get hits from the acceptance window
982 set<Int_t> setTy, intersect;
983 for (multimap<Double_t, Int_t>::iterator mit1 = mityb; mit1 != mitye; ++mit1)
984 if (fmapHits[mit1->second].used == 0)
985 setTy.insert(mit1->second);
986 set_intersection(setTx.begin(), setTx.end(), setTy.begin(), setTy.end(),
987 std::inserter(intersect, intersect.begin()));
988 // if (fVerbose > 0) cout << " Intersect: " << ista << " " << setTx.size() << " " << setTy.size() << " "
989 // << intersect.size() << endl;
990 Double_t xp = fmapHits[ih].xyz[0]; // AZ-280722
991 Double_t zp = fmapHits[ih].xyz[2]; // AZ-280722
992 // int n0 = fMap2[ista].size(); //AZ-241022 - for debug
993
994 for (set<Int_t>::iterator sit = intersect.begin(); sit != intersect.end(); ++sit) {
995 CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(*sit);
996 if (hit->GetUniqueID())
997 continue; // used hit
998 if (fmapHits[*sit].used)
999 continue;
1000 if (fBy < 0.2 && TMath::Abs(tx - fmapHits[*sit].tx) > fdTanXB0[fPass])
1001 continue; // AZ-010123 - zero mag. field
1002 if (fExact) {
1003 // Exact track ID match
1004 set<Int_t> ids1 = GetHitId(*sit, idmaxP);
1005 // if (ids.count(aaa.idmaxP) == 0) continue;
1006 if (idmaxP != aaa.idmaxP)
1007 continue;
1008 }
1009 candvec aaa1 = aaa;
1010 aaa1.stahit[istanext] = *sit; // second hit of the doublet
1011 Double_t dx = xp - fmapHits[*sit].xyz[0]; // AZ-280722
1012 Double_t dz = zp - fmapHits[*sit].xyz[2]; // AZ-280722
1013 aaa1.lengxz = TMath::Sqrt(dx * dx + dz * dz); // AZ-280722 - doublet length in XZ-plane
1014 aaa1.momxz =
1015 Curv3(aaa1, aaa1, aaa1, 0); // AZ-310722 - pxz-estimate with the origin at the primary vertex
1016 if (fBy > 0.2
1017 && (TMath::Abs(aaa1.momxz) > fPTcut[fPass] || TMath::Abs(aaa1.momxz) > fCurvSta[ista]))
1018 continue; // AZ-010123
1019 else if (fBy < 0.2 && TMath::Abs(aaa1.momxz) / fZmean[1] * fZmean[istanext] > fPTcut[fPass])
1020 continue; // AZ-010123 - zero mag. field
1021 // AZ-300922 aaa1.ty = fmapHits[*sit].ty; //AZ-290822
1022 // AZ-220623 aaa1.ty = ty - fmapHits[*sit].ty; //AZ-300922
1023 aaa1.ty = (fmapHits[ih].xyz[1] - fmapHits[*sit].xyz[1]) / aaa1.lengxz; // AZ-220623
1024 aaa1.tx = tx - fmapHits[*sit].tx; // AZ-010123 - for zero mag. field
1025 // AZ-250822 aaa1.varx = varx;
1026 // cout << " Doublet: " << aaa1.idmaxP << " " << aaa1.stahit.begin()->first << endl; // debug output
1027
1028 if (fOut && fOut % 2 == 0) { // AZ-210922 for writing out
1029 string code = to_string(ih) + "-" + to_string(*sit);
1030 if (fMap2Out.find(code) == fMap2Out.end())
1031 fMap2Out[code] = fPass;
1032 }
1033
1034 // AZ-221022
1035 Double_t tmvaOut = 9.9;
1036 if (fUseTMVA == 2)
1037 tmvaOut = TMVAOutput2(aaa1);
1038 // if (fUseTMVA) outf << tmvaOut << " " << fPass << endl;
1039 // if (tmvaOut < 0.02) continue;
1040 if (tmvaOut < 0.01)
1041 continue;
1042 aaa1.nskips = nskips;
1043
1044 mapDoublets.insert(pair<Double_t, candvec>(-tmvaOut, aaa1)); // AZ-251022
1045 } // for (set<Int_t>::iterator sit = intersect.begin();
1046
1047 // AZ-251022
1048 int n2tot = 0, n2max = (fUseTMVA == 2) ? fNbranches * 2 : 999;
1049
1050 for (auto ait = mapDoublets.begin(); ait != mapDoublets.end(); ++ait) {
1051 fCandVec2[ista].push_back(ait->second);
1052 fMap2[ista].insert(pair<int, int>(ih, fCandVec2[ista].size() - 1));
1053 ++n2tot;
1054 if (n2tot >= n2max)
1055 break;
1056 }
1057 // cout << " *** Doublets : " << ista << " " << ih << " " << fMap2[ista].size()-n0 << endl;
1058 } // for (unordered_multimap<int,int>::iterator mit = fTxBins[ista].begin();
1059 } // for (int istanext = ista + 1;
1060
1061 Int_t ncand = fCandVec2[ista].size();
1062 if (fVerbose >= 1) {
1063 cout << " Doublet stat: " << ista << " " << ncand << endl;
1064 for (int i2 = 0; i2 < ncand; ++i2)
1065 PrintHits(fCandVec2[ista][i2]);
1066 }
1067 } // for (Int_t ista = 0;
1068}
1069
1070// -------------------------------------------------------------------------
1071
1072int BmnStsVectorFinderV9::FindBin(int ista, Float_t binw, Float_t tx)
1073{
1074 // Find bin number
1075
1076 return (tx + fTanXmax[ista]) / binw;
1077}
1078
1079// -------------------------------------------------------------------------
1080
1081Float_t BmnStsVectorFinderV9::BinCenterTx(int ista, int ibin, Float_t binw)
1082{
1083 // Get Tx bin center
1084
1085 return -fTanXmax[ista] + binw * (ibin % 1000 + 0.5);
1086}
1087
1088// -------------------------------------------------------------------------
1089
1090void BmnStsVectorFinderV9::DTxWindow(int ista2, Float_t tx, int idu, Float_t curvCut, Float_t* dtx)
1091{
1092 // Compute dTx window
1093
1094 // void DTxWindow(Double_t curv)
1095 // Float_t params[2][4] = { {-0.00995846, -81.1826, -181.066, -94.9851},
1096 // {0.0032616, 79.4583, 163.907, -415.319} };
1097 // Float_t tMinMax[2] = {-0.21, 0.1};
1098 TF1 ff("ff", "pol3", -1, 1);
1099 TVector3 xyz[3];
1100
1101 int bintx = FindBin(ista2 / 100, fTxStep, tx) + 1000 * idu;
1102 // if (fStatData[ista2].find(bintx) == fStatData[ista2].end()) {cout << ista2 << " " << tx << " " << bintx << " " <<
1103 // fStatData[ista2].size() << " " << endl; exit(0);}
1104 vector<Float_t>& params = fStatData[ista2][bintx];
1105 Float_t tMinMax[2] = {params[4], params[9]};
1106
1107 for (int ilr = 0; ilr < 2; ++ilr) {
1108 int ip = ilr * 5;
1109 ff.SetParameters(params[ip], params[ip + 1], params[ip + 2], params[ip + 3]);
1110 Float_t cmax = ff.Eval(tMinMax[ilr]);
1111 // if (cmax <= curv) { cout << tMinMax[ilr] << endl; continue; }
1112 if (cmax <= curvCut) {
1113 dtx[ilr] = tMinMax[ilr];
1114 continue;
1115 }
1116 xyz[1][0] = tMinMax[ilr] * curvCut / cmax;
1117 xyz[1][1] = ff.Eval(xyz[1][0]);
1118 Float_t dx = tMinMax[ilr] / 10;
1119 for (int j = -1; j < 2; j += 2) {
1120 xyz[j + 1][0] = xyz[1][0] + j * dx;
1121 xyz[j + 1][1] = ff.Eval(xyz[j + 1][0]);
1122 }
1123 // xyz[0].Print();
1124 // xyz[1].Print();
1125 // xyz[2].Print();
1126 // Local parabolic approximation
1127 TVector3 cba = Parabola(xyz[0], xyz[1], xyz[2]);
1128 // cba.Print();
1129 Float_t d = cba[1] * cba[1] - 4 * cba[2] * (cba[0] - curvCut);
1130 // cout << d << endl;
1131 Float_t x1 = (-cba[1] + TMath::Sqrt(d)) / 2 / cba[2];
1132 Float_t x2 = (-cba[1] - TMath::Sqrt(d)) / 2 / cba[2];
1133 Float_t dthx = x1;
1134 if (x1 * x2 < 0) {
1135 if ((ilr == 0 && dthx > 0) || (ilr == 1 && dthx < 0))
1136 dthx = x2;
1137 } else
1138 dthx = TMath::Sign(TMath::Min(TMath::Abs(x1), TMath::Abs(x2)), x1);
1139 dtx[ilr] = dthx;
1140 if (ilr == 0)
1141 dtx[ilr] = TMath::Min(dtx[ilr], -0.01f);
1142 else
1143 dtx[ilr] = TMath::Max(dtx[ilr], 0.01f);
1144 // cout << x1 << " " << x2 << " " << dthx << endl;
1145 }
1146}
1147
1148// -------------------------------------------------------------------------
1149
1150Bool_t BmnStsVectorFinderV9::CheckVarx(Double_t dx, Double_t dz, Double_t tx, Double_t distxz, Double_t& varx)
1151{
1152 // AZ-310722 - Check Varx
1153
1154 varx = (dx / dz - tx) / distxz;
1155 // AZ-250822 if (TMath::Abs(varx) > fVarx[fPass]) return kFALSE;
1156 return kTRUE;
1157}
1158
1159// -------------------------------------------------------------------------
1160
1161void BmnStsVectorFinderV9::BuildTriplets()
1162{
1163 // Build triplets
1164
1165 // Int_t idmaxP = 0;
1166
1167 // AZ-130722 pair<multimap<Int_t,int>::iterator,multimap<Int_t,int>::iterator> ret;
1168 // AZ-130722 multimap<Int_t,int>::iterator mit, mit1;
1169 pair<unordered_multimap<Int_t, int>::iterator, unordered_multimap<Int_t, int>::iterator> ret; // AZ-130722
1170 unordered_multimap<Int_t, int>::iterator mit, mit1; // AZ-130722
1171 // AZ-140725 Double_t pxzCut = (fUseTMVA) ? 0.3 : 0.2; //AZ-041022 - pxz difference cut
1172 // pxzCut = 0.3; // AZ-220625
1173 // if (fBy < 0.2) pxzCut = 1.0; //AZ-261222 - zero magnetic field
1174 // const Float_t pxzScale4[19] = {5.0, 3.8, 3.0, 2.0, 1.2, 1.0, 0.9, 1.1, 1.0}; // AZ-260623 - Run 8
1175 // const Float_t pxzScale5[19] = {7.0, 7.0, 3.8, 3.0, 2.0, 1.2, 1.0, 0.9, 1.1, 1.0}; // AZ-220625 - Run 9
1176 // const Float_t dtymax9[19] = {0.014, 0.019, 0.020, 0.016, 0.009, 0.007, 0.007, 0.006, 0.006, 0.006, 0.006, 0.006};
1177 // // AZ-260625 - sigma
1178 const Float_t dtymax9[19] = {0.070, 0.095, 0.100, 0.080, 0.045, 0.035,
1179 0.035, 0.030, 0.030, 0.030, 0.03, 0.03}; // AZ-260625 - 5*sigma
1180 // const Float_t beamWidth9[19] = {0.19, 0.12, 0.10, 0.10, 0.11, 0.12, 0.14, 0.17, 0.19, 0.22}; // AZ-260624 - sigma
1181 // const Float_t beamWidth9[19] = {0.95, 0.60, 0.50, 0.50, 0.55, 0.60, 0.70, 0.85, 0.95, 1.10, 1.1, 1.1}; //
1182 // AZ-260624 - 5*sigma const Float_t beamWidth9[19] = {0.19, 0.12, 0.10, 0.09, 0.10, 0.10, 0.11, 0.13, 0.15, 0.19};
1183 // // AZ-260624 - sigma, fit with "w"
1184 const Float_t beamWidth9[19] = {0.95, 0.60, 0.50, 0.45, 0.50, 0.50,
1185 0.55, 0.65, 0.75, 0.95, 1.0, 1.0}; // AZ-260624 - 5*sigma
1186
1187 // const Float_t *pxzScale = (fNSi == 4) ? pxzScale4 : pxzScale5; // AZ-230625
1188 const Float_t *dtymax = dtymax9, *beamWidth = beamWidth9; // AZ-270625
1189 if (fNSi == 4) { // AZ-270625
1190 // pxzScale = pxzScale4;
1191 dtymax = &dtymax9[1];
1192 beamWidth = &beamWidth9[1];
1193 }
1194
1195 for (Int_t ista = 0; ista < fNsta - 2; ++ista) {
1196 discarded = 0;
1197 Int_t nTra = fCandVec2[ista].size();
1198 // AZ-280323 Int_t istanext = ista + 1;
1199 // AZ-280323 if (nTra == 0 || fCandVec2[istanext].size() == 0) continue;
1200 if (nTra == 0)
1201 continue;
1202
1203 for (mit = fMap2[ista].begin(); mit != fMap2[ista].end(); ++mit) {
1204 candvec& aaa = fCandVec2[ista][mit->second]; // first doublet
1205 Int_t istanext = ista + 1 + aaa.nskips; // AZ-280323
1206 if (fCandVec2[istanext].size() == 0)
1207 continue; // AZ-280323
1208 // AZ-300922 Double_t tany = fmapHits[aaa.stahit.begin()->second].ty; //AZ-2908822
1209 // AZ-300922 tany *= (fTanXmax[ista] / fTanXmax[istanext+1]); //AZ-290822
1210 // Double_t tany = aaa.ty; //AZ-300922
1211 Double_t tanx = aaa.tx; // AZ-010123
1212
1213 ret = fMap2[istanext].equal_range(aaa.stahit.rbegin()->second);
1214 multimap<Double_t, candvec> mapDoublets; // AZ-110722 second doublets, attached to one first doublet
1215 int nbr = 0, newtr = 0, newtr3 = 0; //, mult2 = 0;
1216 // Double_t c2max = 999.0;
1217
1218 // AZ-090525 - for TMVA
1219 /*
1220 if (fUseTMVA && fBy > 0.2) {
1221 for (mit1 = ret.first; mit1 != ret.second; ++mit1) {
1222 candvec &aaa2 = fCandVec2[istanext][mit1->second]; // second doublet
1223 if (TMath::Abs(tany - aaa2.ty) > fdTanY3[fPass]) continue;
1224 ++mult2; // doublet multiplicity in triplets
1225 }
1226 }
1227 */
1228
1229 for (mit1 = ret.first; mit1 != ret.second; ++mit1) {
1230 candvec aaa1;
1231 // AZ-150722 Int_t newtr = (mit1 == ret.first) ? 1 : 0;
1232
1233 // Int_t nhits = 3;
1234 aaa1 = aaa;
1235 candvec& aaa2 = fCandVec2[istanext][mit1->second]; // second doublet
1236 if (fExact) {
1237 // Exact track ID match
1238 if (aaa.idmaxP != aaa2.idmaxP)
1239 continue;
1240 }
1241 // AZ-310722 - Check doublet matching criteria
1242 // if (TMath::Abs(aaa.varx - aaa2.varx) > 0.002) continue; //AZ-050822
1243 // AZ-070822 if (TMath::Abs(aaa.ty - aaa2.ty) > 0.01) continue; //AZ-040822
1244 // if (TMath::Abs(aaa.ty - aaa2.ty) > fdTanY3[fPass]) continue; //AZ-250822
1245 // cout << TMath::Abs(tany - aaa2.ty) << " " << fdTanY3[fPass] << " "
1246 // << TMath::Abs(tanx - aaa2.tx) << " " << fdTanXB0[fPass] << endl;
1247 // AZ-230623 if (TMath::Abs(tany - aaa2.ty) > fdTanY3[fPass]) continue; //AZ-290822
1248 if (fBy < 0.2 && TMath::Abs(tanx - aaa2.tx) > fdTanXB0[fPass])
1249 continue; // AZ-010123 - for zero mag. field
1250 /* AZ-090525
1251 if (fOut / 3 > 0) { // AZ-260922 for writing out
1252 string code = to_string(aaa.stahit.begin()->second) + "-" +
1253 to_string(aaa.stahit.rbegin()->second)
1254 + "-" + to_string(aaa2.stahit.rbegin()->second);
1255 if (fMap3Out.find(code) == fMap3Out.end())
1256 fMap3Out[code] = fPass;
1257 }
1258 */
1259
1260 /* AZ-260625 - reldiff is not a good criterion for first stations and high momenta for a wide beam
1261 Float_t reldiff = TMath::Abs(aaa.momxz-aaa2.momxz) / (TMath::Abs(aaa.momxz)+TMath::Abs(aaa2.momxz));
1262 if (fBy > 0.2 && reldiff > 0.9) continue; //AZ-260623
1263 if (fBy > 0.2 && reldiff / pxzScale[ista] > pxzCut) continue; //AZ-260623
1264 //if (fBy < 0.2 && TMath::Abs(aaa.momxz-aaa2.momxz) / TMath::Sqrt(fZmean[ista]) > pxzCut) continue;
1265 //AZ-311222
1266 */
1267 // AZ-280625 if (TMath::Abs (aaa.ty - aaa2.ty) > dtymax[ista]) continue; // AZ-260625
1268 Double_t dtym = dtymax[ista]; // AZ-280625
1269 int dsta = aaa2.stahit.rbegin()->first - ista; // AZ-280625
1270 if (dsta > 2)
1271 dtym *= (2.0 / dsta); // AZ-280625
1272 if (TMath::Abs(aaa.ty - aaa2.ty) > dtym)
1273 continue; // AZ-260625
1274
1275 // AZ-260625 Use parabolic approximation to find X-position of triplet at Ztarg
1276 map<int, int>::iterator mitt = aaa.stahit.begin();
1277 TVector3 points[3];
1278 points[0] = fmapHits[mitt->second].xyz;
1279 ++mitt;
1280 points[1] = fmapHits[mitt->second].xyz;
1281 points[2] = fmapHits[aaa2.stahit.rbegin()->second].xyz;
1282 for (int j3 = 0; j3 < 3; ++j3)
1283 points[j3] -= fXyzv;
1284 TVector3 cba = Parabola(points[0], points[1], points[2], 2);
1285 Double_t dist = TMath::Abs(cba[0]) / TMath::Sqrt(points[0].Z());
1286 // cout << " ddd " << ista << " " << points[0].Z() << " " << points[1].Z() << " " << points[2].Z() << "
1287 // " << dist << endl;
1288 Double_t beamW = beamWidth[ista] * 2; // AZ-290625
1289 if (dsta > 2)
1290 beamW *= (2.0 / dsta); // AZ-290625
1291 if (dist > beamW)
1292 continue;
1293
1294 aaa1.stahit[aaa2.stahit.rbegin()->first] = aaa2.stahit.rbegin()->second; // third hit of the triplet
1295
1296 Double_t tmvaOut = 9.9; // AZ-300922
1297 /* AZ-200525
1298 if (fUseTMVA)
1299 tmvaOut = TMVAOutput3(aaa1, aaa2); // AZ-300922
1300 // if (fUseTMVA) outf << tmvaOut << " " << fPass << endl; //AZ-300922
1301 if (tmvaOut < 0.02)
1302 continue; // AZ-300922 - use TMVA MLPBNN method
1303 */
1304
1305 map<Int_t, Int_t>& hitMap = aaa1.stahit;
1306
1307 CbmStsHit* hit = nullptr;
1308 CbmStsTrack track;
1309 TArrayI& hits = *track.GetStsHits();
1310 hits.Set(3);
1311 Int_t indx = 0; //, iok = 1;
1312
1314 for (map<Int_t, Int_t>::iterator mit2 = hitMap.begin(); mit2 != hitMap.end(); ++mit2) {
1315 hits[indx++] = mit2->second;
1316 hit = (CbmStsHit*)fHitArray->UncheckedAt(mit2->second);
1317 // cout << " hit: " << mit2->first << " " << mit2->second << endl;
1318 // hit->SetDx(0.08/TMath::Sqrt(12.0));
1319 hit->SetDx(0.08 / TMath::Sqrt(12.0) * 1.2); // AZ-210125 - as in vers.1
1321 // hit->SetDx(0.015);
1322 hit->SetDy(0.1234); // AZ-210125 - as in vers.1
1323 // AZ-030123 if (hit->GetStationNr() <= 3) hit->SetDx(0.02/TMath::Sqrt(12.0));
1324 if (hit->GetStationNr() <= fNSi)
1325 hit->SetDx(0.02 / TMath::Sqrt(12.0)); // AZ-210125 - as in vers.1
1326 // if (hit->GetStationNr() <= fNSi) hit->SetDx(0.01); //AZ-040123 - early estimate
1328 // if (hit->GetStationNr() % 2 != 0) hit->SetCovXY(1.968e-3); // this is for 2-D hits
1329 // else hit->SetCovXY(-1.968e-3);
1330 // hit->SetDy(0.03); // for test
1331 if (indx == 3)
1332 break;
1333 }
1334
1335 // fitter.DoFit(&track);
1336 // float chi2 = track.GetChi2();
1337 // AZ-310722 Double_t ty = 0.0;
1338 Double_t ty = aaa.ty; // AZ-310722
1339 // AZ-280722 float chi2 = LinearFit (&track, newtr, ty) / hit->GetDy() / hit->GetDy();
1340 float chi2 = LinearFit(aaa, aaa2, &track, newtr, ty) / hit->GetDy() / hit->GetDy(); // AZ-280722
1341 ++newtr; // AZ-150722
1342 // float chi2 = LinearFit (&track, newtr, ty) / 0.09 / 0.09;
1343 // cout << "3hit track Chi2 " << track.GetChi2() << " " << nhits << endl;
1344 // if Chi2 > 10.0 exclude track cand
1345 // TODO DZ add index so I know how many candidates were discarded
1346 // if (aaa1.idmaxP == 78 || aaa1.idmaxP == 222 || aaa1.idmaxP == 127)
1347 // cout << " Triplet: " << aaa1.idmaxP << " " << aaa1.stahit.begin()->first << " " << chi2 << endl; //
1348 // debug output
1349
1350 if (chi2 > 10.0) {
1351 // if (chi2 > 20.0) { //AZ-010123
1352 discarded++;
1353 continue;
1354 }
1355 // Triplet curvature
1356 aaa1.momxz = Curv3(aaa, aaa2, aaa1, newtr3);
1357 newtr3 = 1;
1358 // AZ-150722 if (TMath::Abs(aaa1.momxz) < 0.2) {
1359 // Double_t ptcut = (fPass < 2) ? 0.3 : 0.2; //AZ-150722
1360 // if (TMath::Abs(aaa1.momxz) < ptcut) { //AZ-150722
1361 if (TMath::Abs(aaa1.momxz) > fPTcut[fPass]) { // AZ-170722 - curvature
1362 // if (TMath::Abs(aaa1.momxz) > fPTcut[fPass] || TMath::Abs(aaa1.momxz) > fCurvSta[ista]) {
1363 // //AZ-310722 - curvature
1364 discarded++;
1365 continue;
1366 }
1367
1368 if (fOut / 3 > 0) { // AZ-090525 for writing out
1369 string code = to_string(aaa.stahit.begin()->second) + "-" + to_string(aaa.stahit.rbegin()->second)
1370 + "-" + to_string(aaa2.stahit.rbegin()->second);
1371 if (fMap3Out.find(code) == fMap3Out.end())
1372 fMap3Out[code] = fPass;
1373 }
1374
1375 aaa1.ty = ty;
1376 /*AZ-110722
1377 fCandVec3[ista].push_back(aaa1);
1378 fMap3[ista].insert(pair<Int_t,int>(mit->first,fCandVec3[ista].size()-1));
1379 string code = "-";
1380 for (map<Int_t,Int_t>::iterator mitr = hitMap.begin(); mitr != hitMap.end(); ++mitr)
1381 code += (to_string(mitr->second) + "-");
1382 fMapCode3[ista][code] = fCandVec3[ista].size() - 1;
1383 */
1384 /*AZ-200525
1385 if (fUseTMVA) {
1386 // AZ-300922 Keep fNbranches with the highest NN response
1387 if (-tmvaOut > c2max && nbr >= fNbranches)
1388 continue;
1389 } else if (chi2 > c2max && nbr >= fNbranches)
1390 continue;
1391 */
1392 aaa1.chi2 = chi2; // AZ-140722
1393 Double_t qual = (fUseTMVA) ? -tmvaOut : chi2; // AZ-300922
1394 aaa1.nskips = aaa.nskips + aaa2.nskips; // AZ-280323
1395 // AZ-300922 mapDoublets.insert(pair<Double_t,candvec>(chi2,aaa1)); //AZ-110722
1396 mapDoublets.insert(pair<Double_t, candvec>(qual, aaa1)); // AZ-300922
1397 // mapDoublets.insert(pair<Double_t,candvec>(TMath::Abs(aaa1.momxz),aaa1)); //AZ-280722
1398 nbr = mapDoublets.size();
1399 // AZ-310722 if (nbr > fNbranches) mapDoublets.erase(mapDoublets.rbegin()->first);
1400 /*AZ-200525
1401 if (nbr > fNbranches) {
1402 multimap<Double_t, candvec>::iterator mit0 = mapDoublets.end();
1403 --mit0;
1404 mapDoublets.erase(mit0);
1405 }
1406 c2max = mapDoublets.rbegin()->first;
1407 */
1408 } // for (mit1 = ret.first;
1409
1410 // AZ-200525
1411 if (fUseTMVA && nbr > fNbranches) {
1412 multimap<Double_t, candvec> mapTmp;
1413 for (multimap<Double_t, candvec>::iterator it = mapDoublets.begin(); it != mapDoublets.end(); ++it) {
1414 auto aaa1 = it->second;
1415 // Double_t tmvaOut = TMVAOutput3 (aaa1, aaa2, mult2); // AZ-090525
1416 // AZ-280525 Double_t tmvaOut = TMVAOutput3 (aaa1, aaa, mult2); // AZ-090525
1417 Double_t tmvaOut = TMVAOutput3(aaa1, aaa, nbr); // AZ-280525
1418 if (tmvaOut < -0.2) {
1419 ++discarded;
1420 continue;
1421 } // AZ-090525 - use TMVA BDTD method
1422 mapTmp.insert(pair<Double_t, candvec>(-tmvaOut, aaa1));
1423 }
1424 mapDoublets.clear();
1425 mapDoublets = mapTmp;
1426 }
1427 //
1428
1429 // Take only first fNbranches with lowest chi2
1430 int nok = 0;
1431 for (multimap<Double_t, candvec>::iterator it = mapDoublets.begin(); it != mapDoublets.end(); ++it) {
1432 if (nok >= fNbranches)
1433 break;
1434 ++nok;
1435 fCandVec3[ista].push_back(it->second);
1436 fMap3[ista].insert(pair<Int_t, int>(mit->first, fCandVec3[ista].size() - 1));
1437 map<Int_t, Int_t>& hitMap = it->second.stahit;
1438 string code = "-";
1439 for (map<Int_t, Int_t>::iterator mitr = hitMap.begin(); mitr != hitMap.end(); ++mitr)
1440 code += (to_string(mitr->second) + "-");
1441 fMapCode3[ista][code] = fCandVec3[ista].size() - 1;
1442 // cout << ista << " " << code << endl;
1443 }
1444 } // for (mit = fCandMap2[ista].begin();
1445
1446 Int_t ncand = fCandVec3[ista].size();
1447 if (fVerbose > 0)
1448 cout << " Triplet stat: " << ista << " " << ncand << " " << discarded << " " << fPass << endl;
1449 } // for (Int_t ista = fNsta-1;
1450}
1451
1452// -------------------------------------------------------------------------
1453
1454void BmnStsVectorFinderV9::BuildTracks()
1455{
1456 // Build tracks starting from triplets
1457
1458 int istaEnd = fNsta - fNhitsMin[fPass] + 1;
1459 istaEnd = TMath::Min(istaEnd, fNsta - 3); // do not allow simple triplets on last 3 stations
1460 // std::multimap<Int_t,int> &candMap = fMap3[ista];
1461 // vector<candvec> &candVec = fCandVec3[ista];
1462 // std::multimap<Int_t,int> *candMap = &fMap2[ista];
1463 // vector<candvec> *candVec = &fCandVec2[ista];
1464
1465 for (Int_t ista = 0; ista < istaEnd; ++ista) {
1466 // std::multimap<Int_t,int> *candMap = &fMap2[ista];
1467 // vector<candvec> *candVec = &fCandVec2[ista];
1468 // AZ-130722 std::multimap<Int_t,int> *candMap = &fMap3[ista];
1469 std::unordered_multimap<Int_t, int>* candMap = &fMap3[ista]; // AZ-130722
1470 vector<candvec>* candVec = &fCandVec3[ista];
1471 Int_t nTra = candMap->size();
1472 // cout << " BuildTracks: " << ista << " " << nTra << " " << istaEnd << " " << fNsta << endl;
1473 if (nTra == 0)
1474 continue;
1475 // Int_t istaup = ista - 2;
1476
1477 // AZ-130722 for (multimap<Int_t,int>::iterator mit = candMap->begin(); mit != candMap->end(); ++mit) {
1478 for (unordered_multimap<Int_t, int>::iterator mit = candMap->begin(); mit != candMap->end(); ++mit) {
1479 candvec& cand = (*candVec)[mit->second];
1480 if (fPass == 0 && cand.nskips > 0)
1481 continue;
1482 map<int, int>::iterator it = cand.stahit.begin();
1483 int ista0 = it->first;
1484 ++it;
1485 int ista1 = it->first;
1486 if (ista1 - ista0 > 1)
1487 continue; // AZ-250423 - isolated first hit is not allowed
1488 ExtendTrack(cand);
1489 }
1490 }
1491}
1492
1493// -------------------------------------------------------------------------
1494
1495void BmnStsVectorFinderV9::ExtendTrack(candvec cand)
1496{
1497 // Extend the track by triplets or doublets
1498
1499 const Double_t c2ndfMax = 30.0; // cut on chi2/NDF
1500 Double_t pxzCut = 0.3; // AZ-261222 - pxz difference cut
1501 // if (fBy < 0.2) pxzCut = 1.0; //AZ-261222 - zero magnetic field
1502 // const Float_t pxzCuts9pos[19] = {0.043, 0.046, 0.049, 0.060, 0.053, 0.047, 0.050, 0.055, 0.064}; // AZ-290625 -
1503 // sigma
1504 const Float_t pxzCuts9pos[19] = {0.215, 0.230, 0.245, 0.300, 0.265,
1505 0.235, 0.250, 0.275, 0.320}; // AZ-270625 - 5*sigma
1506 // const Float_t pxzCuts9neg[19] = {0.047, 0.043, 0.037, 0.034, 0.032, 0.027, 0.025, 0.024, 0.029}; // AZ-290625 -
1507 // sigma
1508 const Float_t pxzCuts9neg[19] = {0.235, 0.215, 0.185, 0.170, 0.160,
1509 0.135, 0.125, 0.120, 0.145}; // AZ-290625 - 5*sigma
1510 const Float_t* pxzCutsPos = (fNSi == 5) ? pxzCuts9pos : &pxzCuts9pos[1]; // AZ-290625
1511 const Float_t* pxzCutsNeg = (fNSi == 5) ? pxzCuts9neg : &pxzCuts9neg[1]; // AZ-290625
1512
1513 Int_t nTra = 0, ista = cand.stahit.rbegin()->first;
1514 // cout << ista << " " << cand.stahit.size() << endl;
1515 /*
1517 if (ista == fNsta-2) {
1518 cand.stahit.erase(cand.stahit.rbegin()->first);
1519 ista = cand.stahit.rbegin()->first;
1520 }
1521 */
1522 // std::multimap<Int_t,int> &candMap = (ista == fNsta-2) ? fMap2[ista] : fMap3[ista]; // doublets or triplets
1523 // vector<candvec> &candVec = (ista == fNsta-2) ? fCandVec2[ista] : fCandVec3[ista]; // doublets or triplets
1524 // std::multimap<Int_t,int> *pCandMap = (ista == fNsta-2) ? &fMap2[ista] : &fMap2[ista]; // doublets or triplets
1525 // vector<candvec> *pCandVec = (ista == fNsta-2) ? &fCandVec2[ista] : &fCandVec2[ista]; // doublets or triplets
1526 // AZ-130722 std::multimap<Int_t,int> *pCandMap = (ista == fNsta-2) ? &fMap2[ista] : &fMap2[ista]; // doublets or
1527 // triplets
1528 std::unordered_multimap<Int_t, int>* pCandMap =
1529 (ista == fNsta - 2) ? &fMap2[ista] : &fMap2[ista]; // AZ-130722 doublets or triplets
1530 vector<candvec>* pCandVec = (ista == fNsta - 2) ? &fCandVec2[ista] : &fCandVec2[ista]; // doublets or triplets
1531 // AZ-130722 multimap<Int_t,int> &candMap = *pCandMap;
1532 unordered_multimap<Int_t, int>& candMap = *pCandMap; // AZ-130722
1533 vector<candvec>& candVec = *pCandVec;
1534
1535 if (ista < fNsta - 1) {
1536 /*
1537 if (ista == fNsta-2) {
1538 // Doublets
1539 candMap = fMap2[ista];
1540 candVec = fCandVec2[ista];
1541 }
1542 */
1543 nTra = candMap.count(cand.stahit.rbegin()->second);
1544 }
1545
1546 // cout << " ExtendTrack " << ista << " " << nTra << endl;
1547 if (nTra == 0) {
1548 // No extension found
1549 if (int(cand.stahit.size()) < fNhitsMin[fPass])
1550 return; // too short track
1551 // cout << " ExtendTrack - store " << cand.stahit.size() << " " << cand.idmaxP << endl;
1552 // FitTrack(cand);
1553 string code("-");
1554 for (map<int, int>::iterator it = cand.stahit.begin(); it != cand.stahit.end(); ++it)
1555 code += (to_string(it->second) + "-");
1556 // cout << cand.stahit.rbegin()->first << " " << ista << endl;
1557 if (fCandCodes[cand.stahit.size()].find(code) == fCandCodes[cand.stahit.size()].end()) {
1558 fCandCodes[cand.stahit.size()].insert(code);
1559 // cout << " 1 " << endl;
1560 // AZ-170122 if (cand.stahit.size() >= 4 && FitTrack(cand) > c2ndfMax) return;
1561 // AZ-110123 if (cand.stahit.size() >= 4 && cand.stahit.rbegin()->first != ista && FitTrack(cand) >
1562 // c2ndfMax) return;
1563 if (int(cand.stahit.size()) >= fNhitsMin[fPass] && cand.stahit.rbegin()->first != ista
1564 && FitTrack(cand) > c2ndfMax)
1565 return; // AZ-110123
1566 // cout << " ExtendTrack - store " << cand.stahit.size() << " " << cand.idmaxP << endl;
1567 fCandVec[cand.stahit.begin()->first].push_back(cand);
1568 }
1569 return;
1570 }
1571
1572 // AZ-130722 multimap<Int_t,int>::iterator mit;
1573 // AZ-130722 pair<multimap<Int_t,int>::iterator,multimap<Int_t,int>::iterator> ret;
1574 unordered_multimap<Int_t, int>::iterator mit; // AZ-130722
1575 pair<unordered_multimap<Int_t, int>::iterator, unordered_multimap<Int_t, int>::iterator> ret; // AZ-130722
1576 ret = candMap.equal_range(cand.stahit.rbegin()->second);
1577 string code = "-";
1578 int istaMid = 0;
1579
1580 // Use Kalman track constraints around extrapolated position
1581 /*
1582 set<int> kalmWind;
1583 if (nTra > 5)
1584 if (cand.stahit.size() >= 4) kalmWind = KalmanWindow(cand, candVec[ret.first->second].stahit.begin()->second);
1585 */
1586
1587 for (mit = ret.first; mit != ret.second; ++mit) {
1588 // candvec &aaa = candVec[mit->second];
1589 candvec aaa = candVec[mit->second];
1590 if (fExact && aaa.idmaxP != cand.idmaxP)
1591 continue;
1592 if (fPass == 0 && aaa.nskips > 0)
1593 continue;
1594
1595 // if (nTra > 5 && cand.stahit.size() >= 4 &&
1596 // kalmWind.find(aaa.stahit.begin()->second) == kalmWind.end()) continue; // outside Kalman window
1597
1598 // Check if there is a triplet started from the last but one point of the current one (cand)
1599 map<Int_t, Int_t>::iterator mitr = aaa.stahit.begin();
1600 ++mitr;
1601 if (code.size() == 1) {
1602 map<Int_t, Int_t>::reverse_iterator mit1 = cand.stahit.rbegin();
1603 ++mit1;
1604 istaMid = mit1->first;
1605 code += (to_string(mit1->second) + "-");
1606 --mit1;
1607 code += (to_string(mit1->second) + "-");
1608 code += (to_string(mitr->second) + "-");
1609 } else {
1610 // Replace last index
1611 int pos = code.rfind("-", code.size() - 2);
1612 code.erase(code.begin() + pos + 1, code.end());
1613 code += (to_string(mitr->second) + "-");
1614 }
1615 // cout << " Check code: " << code << " " << istaMid << endl;
1616 // Check extension possibility
1617 int ok = 1;
1618 if (fMapCode3[istaMid].find(code) == fMapCode3[istaMid].end())
1619 ok = 0; // no triplet found
1621 //*
1622 candvec* bbb = nullptr;
1623 if (ok) {
1624 bbb = &fCandVec3[istaMid][fMapCode3[istaMid][code]];
1625 // Double_t dty = TMath::Abs (cand.ty - bbb.ty);
1626 // if (dty * 2 / (TMath::Abs(cand.ty) + TMath::Abs(bbb.ty)) > 0.1) continue; // different slopes
1627 // if (dty > 0.005) continue; // different slopes
1628 // Curvatures
1629 // cout << cand.momxz << " " << bbb.momxz << endl;
1630 if ((fPass == 0 && bbb->nskips > 0) || (cand.nskips + bbb->nskips > 2))
1631 ok = 0;
1632 if (ok && cand.stahit.size() > 2) {
1633 Double_t dp = TMath::Abs(cand.momxz - bbb->momxz);
1634 // if (dp * 2 / (TMath::Abs (cand.momxz) + TMath::Abs(bbb->momxz)) > 0.3) {
1635 // AZ-060822 if (dp * 2 / (TMath::Abs (cand.momxz) + TMath::Abs(bbb->momxz)) > 0.2) {
1636 // AZ-261222 if (dp / (TMath::Abs (cand.momxz) + TMath::Abs(bbb->momxz)) > 0.3) { //AZ-060822
1637 // AZ-270625 if (TMath::Abs(fBy) > 0.2 && dp / (TMath::Abs (cand.momxz) + TMath::Abs(bbb->momxz)) >
1638 // pxzCut) { //AZ-311222
1639 int istabeg = cand.stahit.begin()->first; // AZ-290625
1640 pxzCut = (cand.momxz > 0) ? pxzCutsPos[istabeg] : pxzCutsNeg[istabeg]; // AZ-290625
1641 int dsta = bbb->stahit.rbegin()->first - istabeg; // AZ-290625
1642 if (dsta > 3)
1643 pxzCut *= (3.0 / dsta); // AZ-290625
1644 if (TMath::Abs(fBy) > 0.2
1645 && dp / (TMath::Abs(cand.momxz) + TMath::Abs(bbb->momxz)) / TMath::Sqrt(TMath::Abs(1 / cand.momxz))
1646 > pxzCut)
1647 { // AZ-270625
1648 // Different curvatures
1649 ok = 0;
1650 }
1651 }
1652 }
1653 if (ok) {
1654 if ((fPass == 0 && aaa.nskips > 0) || (cand.nskips + aaa.nskips > 2))
1655 ok = 0;
1656 if (aaa.stahit.size() == 3) {
1657 // Extension by a triplet
1658 if (ok) {
1659 Double_t dp = TMath::Abs(cand.momxz - aaa.momxz);
1660 // if (dp * 2 / (TMath::Abs (cand.momxz) + TMath::Abs(aaa.momxz)) > 0.3) {
1661 // AZ-070822 if (dp * 2 / (TMath::Abs (cand.momxz) + TMath::Abs(aaa.momxz)) > 0.2) {
1662 // AZ-261222 if (dp / (TMath::Abs (cand.momxz) + TMath::Abs(aaa.momxz)) > 0.3) { //AZ-070822
1663 // AZ-270625 if (TMath::Abs(fBy) > 0.2 && dp / (TMath::Abs (cand.momxz) + TMath::Abs(aaa.momxz)) >
1664 // pxzCut) { //AZ-020123
1665 int istabeg = cand.stahit.begin()->first; // AZ-290625
1666 pxzCut = (cand.momxz > 0) ? pxzCutsPos[istabeg] : pxzCutsNeg[istabeg]; // AZ-290625
1667 int dsta = aaa.stahit.rbegin()->first - istabeg; // AZ-290625
1668 if (dsta > 3)
1669 pxzCut *= (3.0 / dsta); // AZ-290625
1670 if (TMath::Abs(fBy) > 0.2
1671 && dp / (TMath::Abs(cand.momxz) + TMath::Abs(aaa.momxz))
1672 / TMath::Sqrt(TMath::Abs(1 / cand.momxz))
1673 > pxzCut)
1674 { // AZ-270625
1675 // different curvatures
1676 ok = 0;
1677 }
1678 }
1679 }
1680 }
1681 if (!ok) {
1682 // No extension found
1683 if (int(cand.stahit.size()) < fNhitsMin[fPass])
1684 continue; // too short track
1685 // cout << " ExtendTrack - store " << cand.stahit.size() << " " << cand.idmaxP << endl;
1686 string codeVec("-");
1687 for (map<int, int>::iterator it = cand.stahit.begin(); it != cand.stahit.end(); ++it)
1688 codeVec += (to_string(it->second) + "-");
1689 if (fCandCodes[cand.stahit.size()].find(codeVec) == fCandCodes[cand.stahit.size()].end()) {
1690 fCandCodes[cand.stahit.size()].insert(codeVec);
1691 if (cand.stahit.size() >= 4) {
1692 // cout << " 2 " << endl;
1693 if (cand.stahit.rbegin()->first == ista)
1694 continue; // AZ-171121
1695 Double_t chi2 = FitTrack(cand);
1696 // cout << " --- Chi2: " << chi2 << endl;
1697 if (chi2 > c2ndfMax)
1698 continue;
1699 }
1700 // cout << " ExtendTrack - store " << cand.stahit.size() << " " << cand.idmaxP << endl;
1701 fCandVec[cand.stahit.begin()->first].push_back(cand);
1702 }
1703 continue; // AZ-280323
1704 }
1705 //*/
1707 candvec candext = cand;
1708 if (cand.stahit.size() < 3)
1709 candext.momxz = bbb->momxz; // update momentum
1710 // candext.ty = bbb.ty;
1711 // Add triplet (or doublet)
1712 for (mitr = aaa.stahit.begin(); mitr != aaa.stahit.end(); ++mitr)
1713 candext.stahit[mitr->first] = mitr->second;
1714 candext.nskips += aaa.nskips;
1715 // cout << " ids: " << aaa.idmaxP << " " << cand.idmaxP << " " << candext.stahit.size() << " " << candext.nskips
1716 // << " " << fPass << endl; if (candext.nskips > 2) exit(0); AZ-120722 if (aaa.idmaxP != cand.idmaxP)
1717 // cand.idmaxP = -1;
1718 if (aaa.idmaxP != cand.idmaxP)
1719 candext.idmaxP = -1; // AZ-120722
1720 // Fit track to fight with combinatorics
1721 if (candext.stahit.size() >= 4) {
1722 Double_t chi2ndf = FitTrack(candext);
1723 // cout << " -3- " << chi2ndf << endl;
1724 if (chi2ndf > c2ndfMax)
1725 continue;
1726 /* This slows down
1727 if (chi2ndf > c2ndfMax) {
1728 string codext = MakeCode(candext);
1729 fCandCodes[candext.stahit.size()].insert(codext);
1730 continue;
1731 }
1732 */
1733 }
1734 ExtendTrack(candext);
1735 } // for (mit = ret.first;
1736}
1737
1738// -------------------------------------------------------------------------
1739
1740Double_t BmnStsVectorFinderV9::FitTrack(candvec& cand)
1741{
1742 // Fit track candidate
1743
1744 const Int_t gkChi2Cut = 5.0; // 10.0;
1745 // Int_t gkChi2Cut = 10.0; //AZ-010123
1746
1747 CbmStsTrack track;
1748 string hitcode;
1749
1750 MakeStsTrack(cand, hitcode, track);
1751 int nhits = track.GetNStsHits();
1752 if (fCandSet[nhits].find(hitcode) != fCandSet[nhits].end())
1753 return 0.0; // track has been saved already
1754
1755 // AZ-160122 fitter.DoFit(&track);
1756 if (nhits <= 4)
1757 fitter.DoFit(&track);
1758 else
1759 FilterHit(cand, track);
1760 /*
1761 else {
1762 // Debug
1763 CbmStsTrack track1(track);
1764 FilterHit (cand, track);
1765 fitter.DoFit(&track1);
1766 cout << " debug " << track.GetChi2() << " " << track1.GetChi2() << endl;
1767 }
1768 */
1769 // track.GetParamFirst()->Print();
1770 // cout << " aaaaaa " << track.GetChi2() << " " << track.GetNStsHits() << endl;
1771 int ndf = TMath::Max(track.GetNDF(), 1);
1772 Double_t chi2ndf = track.GetChi2() / ndf;
1773 cand.param = *(track.GetParamLast());
1774 cand.chi2 = track.GetChi2();
1775 if (chi2ndf < gkChi2Cut && nhits >= fNhitsMin[fPass]) {
1776 // Refit with material budget
1777 // fitter.DoFit(&track);
1778 fitter.Fit(&track); // AZ
1779 if (track.GetChi2() / ndf > gkChi2Cut) {
1780 // cout << " yyyyyy " << " " << track.GetChi2() / ndf << endl;
1781 // fitter.Fit(&track); //AZ - for debugging
1782 return chi2ndf;
1783 }
1784 // Save good tracks
1785 Double_t qual = nhits;
1786 qual += (100.0 - TMath::Min(TMath::Abs(track.GetChi2()), 100.0)) / 101.0;
1787 fTracks.insert(pair<Double_t, CbmStsTrack>(-qual, track));
1788 fCandSet[nhits].insert(hitcode);
1789 if (fVerbose > 0) {
1790 cout << " FitTrack: " << track.GetChi2() << " ";
1791 PrintHits(cand);
1792 }
1793 }
1794 // fCandSet[nhits].insert(hitcode); - don't use this !!!
1795 return chi2ndf;
1796}
1797
1798// -------------------------------------------------------------------------
1799//*
1800Double_t BmnStsVectorFinderV9::FilterHit(candvec& cand, CbmStsTrack& track)
1801{
1802 // Filter last hit on track candidate
1803
1804 // const Int_t gkChi2Cut = 10.0;
1805 // AZ const Int_t gkChi2Cut = 20.0; //AZ-010123
1806 /*
1807 CbmStsTrack track;
1808 string hitcode;
1809
1810 MakeStsTrack(cand, hitcode, track);
1811 */
1812 Int_t nhits = cand.stahit.size();
1813 // CbmStsHit *hit = (CbmStsHit*) fHitArray->UncheckedAt(track.GetStsHitIndex(nhits-1));
1814 /*cout << " !!! parameters: " << cand.stahit.size() << " " << cand.param.GetX() << " " << cand.param.GetY()
1815 << " " << cand.param.GetZ() << " "
1816 << hit->GetZ() << " " << cand.chi2 << " " << cand.code << endl;*/
1817 // AZ-050222 fitter.Extrapolate(&track,hit->GetZ(),track.GetParamLast());
1818 // cout << hit->GetZ() << " " << track.GetParamLast()->GetZ() << endl;
1819
1820 CbmKFTrack kftr;
1821 kftr.SetStsTrack(track, kFALSE); // last parameter
1822 fitter.SetKFHits(kftr, &track);
1823 Bool_t downstream = kTRUE;
1824 // Double_t qp0 = 0.0;
1825 Double_t qp0 = track.GetParamLast()->GetQp();
1826 // Double_t *cov = kftr.GetCovMatrix();
1827 // for (Int_t j = 0; j < 15; ++j) cov[j] = 0.0;
1828 // cov[0] = cov[2] = cov[5] = cov[9] = cov[14] = 10000.0;
1829 // kftr.GetRefChi2() = 0.0;
1830 // Double_t chi20 = 0.0, dchi2 = 0.0;
1831 // Double_t chi20 = kftr.GetRefChi2();
1832
1833 /*CbmKF *KF =*/CbmKF::Instance();
1834 CbmKFHit* h = NULL;
1835 // Int_t istold = 0, imax = 0;
1836 Bool_t err = kFALSE;
1837
1838 // for (Int_t i = 0; i < nhits; i++) {
1839 h = kftr.GetHit(nhits - 1);
1840 err = err || h->Filter(kftr, downstream, qp0);
1841 // Double_t dchi2 = kftr.GetRefChi2() - chi20;
1842 // cout << err << " " << dchi2 << " " << kftr.GetRefChi2() << " " << cand.chi2+dchi2 << endl;
1843 kftr.GetTrackParam(*track.GetParamLast());
1844 track.SetChi2(kftr.GetRefChi2() + err * 1000);
1845 track.SetNDF(2 * nhits - 5);
1846 /*
1847 if (fCandSet[track.GetNStsHits()].find(hitcode) != fCandSet[track.GetNStsHits()].end()) return 0.0; // track has
1848 been saved already
1849
1850 fitter.DoFit(&track);
1851 //track.GetParamFirst()->Print();
1852 //cout << " aaaaaa " << track.GetChi2() << " " << track.GetNStsHits() << endl;
1853 int ndf = TMath::Max (track.GetNDF(), 1);
1854 Double_t chi2ndf = track.GetChi2() / ndf;
1855 if (chi2ndf < gkChi2Cut && track.GetNStsHits() >= fNhitsMin[fPass]) {
1856 // Save good tracks
1857 Double_t qual = track.GetNStsHits();
1858 qual += (100.0 - TMath::Min(TMath::Abs(track.GetChi2()),100.0)) / 101.0;
1859 fTracks.insert(pair<Double_t,CbmStsTrack>(-qual,track));
1860 fCandSet[track.GetNStsHits()].insert(hitcode);
1861 //cout << " FitTrack: " << hitcode << " " << track.GetChi2() << endl;
1862 }
1863 return chi2ndf;
1864 */
1865 return 0.0;
1866}
1867//*/
1868// -------------------------------------------------------------------------
1869
1870void BmnStsVectorFinderV9::MakeStsTrack(candvec& cand, string& hitcode, CbmStsTrack& track)
1871{
1872 // Fill STS track information
1873
1874 map<Int_t, Int_t>& hitMap = cand.stahit;
1875 CbmStsHit* hit = nullptr;
1876 TArrayI& hits = *track.GetStsHits();
1877 hits.Set(hitMap.size());
1878 Int_t indx = 0; //, iok = 1;
1879 hitcode = "-";
1880
1882 for (map<Int_t, Int_t>::iterator mit2 = hitMap.begin(); mit2 != hitMap.end(); ++mit2) {
1883 hits[indx++] = mit2->second;
1884 hit = (CbmStsHit*)fHitArray->UncheckedAt(mit2->second);
1885 // cout << " hit: " << mit2->first << " " << mit2->second << endl;
1886 hit->SetDx(0.08 / TMath::Sqrt(12.0)); // AZ-210125 - as in vers.1
1888 // hit->SetDx(0.015);
1889 hit->SetDy(0.1234); // AZ-210125 - as in vers.1
1890 // AZ-030123 if (hit->GetStationNr() <= 3) hit->SetDx(0.02/TMath::Sqrt(12.0));
1891 if (hit->GetStationNr() <= fNSi)
1892 hit->SetDx(0.02 / TMath::Sqrt(12.0)); // AZ-210125 - as in vers.1
1893 // if (hit->GetStationNr() % 2 != 0) hit->SetCovXY(1.968e-3); // this is for 2-D hits
1894 // else hit->SetCovXY(-1.968e-3);
1895 // hit->SetDy(0.03); // for test
1896 CbmStsStation* stat = CbmStsDigiScheme::Instance()->GetStationByNr(hit->GetStationNr());
1897 Double_t dx = stat->GetSector(0)->GetDx();
1898 // if (dx < 0.02) dx *= 2; // scale up for Si
1899 // AZ-120722 if (dx > 0.01) dx *= 1.2; // for GEMs
1900 if (dx > 0.02)
1901 dx *= 1.2; // AZ-120722 for GEMs
1902 hit->SetDx(dx / TMath::Sqrt(12.0)); // AZ-210125 - as in vers.1
1903 // if (dx < 0.02) hit->SetDx(0.01); //AZ-040123 - early estimate for Si
1905 hitcode += (to_string(mit2->second) + "-");
1906 }
1907 track.SetParamFirst(cand.param);
1908 track.SetParamLast(cand.param);
1909 track.SetChi2(cand.chi2);
1910 cand.code = hitcode;
1911}
1912
1913// -------------------------------------------------------------------------
1914
1915void BmnStsVectorFinderV9::ExtendTracks(Int_t ista)
1916{
1917 // Extend tracks by going upstream
1918
1919 // const Double_t zsta[19] = {1.132, 32.750, 64.150, 96.750, 129.300, 160.800, 192.100}; // average hit Z-position -
1920 // 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}; //
1921 // average hit Z-position - Run 7 !!!
1922 /*
1923 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
1924 hit Z-position - Run 8 !!! const Double_t pars[19][4] = { {1.30358, 20.8953, 0.857369, 0.56339}, {1.64202, 7.18813,
1925 0.596002, 0.694601}, {3.10765, 3.21564, 0.291252, 1.08062}, {2.25484, 3.81810, 0.430724, 0.840519},
1926 {3.38309, 4.33504, 0.349704, 1.07882}, {3.99252, 3.35587, 0.302402, 1.65050},
1927 {10.3321, 4.24742, 0.168995, 1.71113} };
1928 */
1929}
1930
1931// -------------------------------------------------------------------------
1932
1933void BmnStsVectorFinderV9::FinishEvent()
1934{
1935 // Finish event
1936 /*
1937 FairRootManager* ioman = FairRootManager::Instance();
1938 */
1939 if (fStream2out.is_open()) {
1940 fStream2out << " Event: " << fEvent << " " << fMap2Out.size() << endl;
1941 for (unordered_map<string, int>::iterator it = fMap2Out.begin(); it != fMap2Out.end(); ++it)
1942 fStream2out << it->first << " " << it->second << "\n";
1943 }
1944 if (fStream3out.is_open()) {
1945 fStream3out << " Event: " << fEvent << " " << fMap3Out.size() << endl;
1946 for (unordered_map<string, int>::iterator it = fMap3Out.begin(); it != fMap3Out.end(); ++it)
1947 fStream3out << it->first << " " << it->second << "\n";
1948 }
1949 ++fEvent;
1950}
1951
1952// -------------------------------------------------------------------------
1953
1954void BmnStsVectorFinderV9::Finish()
1955{
1956 printf("Work time of BmnStsVectorFinderV9: %4.2f sec.\n", workTime);
1957 if (fStream2out.is_open())
1958 fStream2out.close();
1959 if (fStream3out.is_open())
1960 fStream3out.close();
1961}
1962
1963//__________________________________________________________________________
1964
1965// AZ-270625 TVector3 BmnStsVectorFinderV9::Parabola(TVector3& pos0, TVector3& pos1, TVector3& pos2)
1966TVector3 BmnStsVectorFinderV9::Parabola(TVector3& pos0, TVector3& pos1, TVector3& pos2, int ix) // AZ-270625
1967{
1968 // Get parabolic track approximation from 3 points (X vs Z for ix = 2 or Y vs X for ix = 0)
1969 // y = a*x^2 + bx + c
1970
1971 Double_t x2[3] = {pos0[2] - fXyzv[2], pos1[2] - fXyzv[2], pos2[2] - fXyzv[2]};
1972 Double_t y2[3] = {pos0[0] - fXyzv[0], pos1[0] - fXyzv[0], pos2[0] - fXyzv[0]};
1973 Double_t x0[3] = {pos0[0], pos1[0], pos2[0]};
1974 Double_t y0[3] = {pos0[1], pos1[1], pos2[1]};
1975 Double_t *x = x0, *y = y0;
1976 if (ix == 2) {
1977 x = x2;
1978 y = y2;
1979 }
1980
1981 Double_t denom = (x[0] - x[1]) * (x[0] - x[2]) * (x[1] - x[2]);
1982 Double_t dy10 = y[1] - y[0];
1983 Double_t dy02 = y[0] - y[2];
1984 Double_t dy21 = y[2] - y[1];
1985 // Double_t a = x[2] * (y[1] - y[0]) + x[1] * (y[0] - y[2]) + x[0] * (y[2] - y[1]);
1986 Double_t a = x[2] * dy10 + x[1] * dy02 + x[0] * dy21;
1987 a /= denom;
1988 // Double_t b = x[0]*x[0] * (y[1] - y[2]) + x[2]*x[2] * (y[0] - y[1]) +
1989 // x[1]*x[1] * (y[2] - y[0]);
1990 Double_t b = -x[0] * x[0] * dy21 - x[2] * x[2] * dy10 - x[1] * x[1] * dy02;
1991 b /= denom;
1992 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])
1993 + x[0] * x[2] * (x[2] - x[0]) * y[1];
1994 c /= denom;
1995 return TVector3(c, b, a);
1996 // 𝐴=𝑥3(𝑦2−𝑦1)+𝑥2(𝑦1−𝑦3)+𝑥1(𝑦3−𝑦2) / (𝑥1−𝑥2)(𝑥1−𝑥3)(𝑥2−𝑥3)
1997 // 𝐵=𝑥21(𝑦2−𝑦3)+𝑥23(𝑦1−𝑦2)+𝑥22(𝑦3−𝑦1) / (𝑥1−𝑥2)(𝑥1−𝑥3)(𝑥2−𝑥3)
1998 // 𝐶=𝑥22(𝑥3𝑦1−𝑥1𝑦3)+𝑥2(𝑥21𝑦3−𝑥23𝑦1)+𝑥1𝑥3(𝑥3−𝑥1)𝑦2 / (𝑥1−𝑥2)(𝑥1−𝑥3)(𝑥2−𝑥3)
1999}
2000
2001// -------------------------------------------------------------------------
2002
2003void BmnStsVectorFinderV9::FitTracks()
2004{
2005 // Fit tracks
2006
2007 const Double_t gkChi2Cut = 5.0; // 10.0;
2008 // const Double_t gkChi2Cut = 10.0; //AZ-010123
2009 // fTracks.clear();
2010
2011 // for (Int_t ista = 0; ista < fNsta-3; ++ista) {
2012 for (Int_t ista = 0; ista < fNsta; ++ista) {
2013 Int_t ntra = fCandVec[ista].size();
2014 if (fVerbose > 0)
2015 cout << " FitTracks: " << ista << " " << ntra << endl;
2016
2017 for (Int_t itra = 0; itra < ntra; ++itra) {
2018 candvec& cand = fCandVec[ista][itra];
2019 // if (cand.nextind >= 0) continue; // track has an extension already processed
2020
2021 map<Int_t, Int_t>& hitMap = cand.stahit;
2022 Int_t nhits = hitMap.size();
2023 // if (nhits < 3) continue; // too short track
2024 if (nhits < fNhitsMin[fPass])
2025 continue; // too short track
2026 CbmStsTrack track;
2027 TArrayI& hits = *track.GetStsHits();
2028 hits.Set(nhits);
2029 Int_t indx = 0, iok = 1;
2030
2031 // Check if the track candidate has been checked already
2032 string hitcode("-");
2033 for (map<Int_t, Int_t>::iterator mit = hitMap.begin(); mit != hitMap.end(); ++mit)
2034 hitcode += (to_string(mit->second) + "-");
2035 // cout << itra << " " << hitcode << endl;
2036
2037 /* AZ-140725
2038 Int_t same = 0; //, sta0 = hitMap.begin()->first;
2039
2040 //same = fCandSet[sta0].count(hitcode);
2041 same = fCandSet[nhits].count(hitcode);
2042 //if (same) cout << " same: " << sta0 << " " << same << " " << endl;
2043 */
2044
2046 // fCandSet[sta0].insert(hitcode);
2047 fCandSet[nhits].insert(hitcode);
2048
2049 for (map<Int_t, Int_t>::iterator mit = hitMap.begin(); mit != hitMap.end(); ++mit) {
2050 hits[indx++] = mit->second;
2051 // Debug
2052 CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(mit->second);
2053 // cout << " hit: " << mit->first << " " << hit->GetX() << " " << hit->GetY() << " " << hit->GetZ() << "
2054 // "
2055 // << hit->GetDx() << " " << hit->GetDy() << endl;
2056 // cout << " hit: " << mit->first << " " << mit->second << endl;
2057 // hit->SetDx(0.08/TMath::Sqrt(12.0));
2059 // hit->SetDx(0.015);
2061 // AZ-030123 if (hit->GetStationNr() <= 3) hit->SetDx(0.02/TMath::Sqrt(12.0));
2062 // if (hit->GetStationNr() <= fNSi) hit->SetDx(0.02/TMath::Sqrt(12.0)); //AZ-030123
2063 // if (hit->GetStationNr() <= 3) hit->SetDx(0.0025);
2064 // if (hit->GetStationNr() % 2 != 0) hit->SetCovXY(1.968e-3); // this is for 2-D hits (CbmFitter)
2065 // else hit->SetCovXY(-1.968e-3);
2066 // hit->SetDy(0.03); // for test
2067 CbmStsStation* stat = CbmStsDigiScheme::Instance()->GetStationByNr(hit->GetStationNr());
2068 Double_t dx = stat->GetSector(0)->GetDx();
2069 // if (dx < 0.02) dx *= 2; // scale up for Si
2070 // AZ-120722 if (dx > 0.01) dx *= 1.2; // scale up for GEMs
2071 if (dx > 0.02)
2072 dx *= 1.2; // AZ-120722 scale up for GEMs
2073 hit->SetDx(dx / TMath::Sqrt(12.0)); // AZ-210125 - as in vers.1
2074 // if (dx < 0.02) hit->SetDx(0.01); //AZ-040123 - early estimate for Si
2076 /*
2077 if (fMatBudgetFileName != "" && hit->GetDz() < 1.e-7) {
2078 // Store material budget for the hit (along the beam)
2079 Double_t zshift = (dx < 0.02) ? 0.5 : 1.5; // Si or GEM
2080 TProfile2D *prof = fMatHistos.lower_bound(hit->GetZ()-zshift)->second;
2081 int ix = prof->FindBin(hit->GetX());
2082 int iy = prof->FindBin(hit->GetY());
2083 Double_t radThick = prof->GetBinContent(ix,iy);
2084 hit->SetDz(radThick/100);
2085 }
2086 */
2087 }
2088
2089 // cout << " before " << track.GetParamFirst()->GetZ() << endl;
2090 // fitter.DoFit(&track);
2091 fitter.Fit(&track); // AZ
2092 // cout << " after " << track.GetParamFirst()->GetZ() << " " << track.GetNStsHits() << " " <<
2093 // track.GetChi2() << endl;
2094 /*
2095 CbmKFTrack kftr = CbmKFTrack(track);
2096 fitter.SetKFHits(kftr, &track);
2097 kftr.Fit(kTRUE); // downstream
2098 kftr.Fit(kFALSE); // upstream
2099 */
2100 if (fVerbose > 0) {
2101 cout << " aaaaaaa " << endl;
2102 // track.GetParamFirst()->Print();
2103 cout << track.GetChi2() << " " << nhits << " " << endl;
2104 }
2105 // AZ-111221 if (track.GetChi2() / track.GetNDF() > gkChi2Cut) {
2106 if (track.GetChi2() / track.GetNDF() > gkChi2Cut && track.GetNStsHits() > 3) {
2107 /*
2108 // Too high chi2/NDF - try to remove one outlier with the largest dChi2
2109 // Based on the code from CbmKFTrackInterface::Fit()
2110 CbmKFTrack kftr;
2111 //kftr.SetStsTrack(track, kFALSE); // last parameter
2112 //Bool_t downstream = kFALSE;
2113 kftr.SetStsTrack(track, kTRUE); // first parameter
2114 fitter.SetKFHits(kftr,&track);
2115 Bool_t downstream = kTRUE;
2116 Double_t qp0 = 0.0;
2117 Double_t *cov = kftr.GetCovMatrix();
2118 for (Int_t j = 0; j < 15; ++j) cov[j] = 0.0;
2119 cov[0] = cov[2] = cov[5] = cov[9] = cov[14] = 10000.0;
2120 kftr.GetRefChi2() = 0.0;
2121 Double_t chi20 = 0.0, dchi2 = 0.0, dchi2max = 0.0;
2122
2123 CbmKF *KF = CbmKF::Instance();
2124 CbmKFHit *h = NULL;
2125 Int_t istold = 0, imax = 0;
2126 Bool_t err = kFALSE;
2127 set<Int_t> indok;
2128
2129 //for (Int_t i = nhits - 1; i >= 0; i--) {
2130 // h = kftr.GetHit( i );
2131 // Int_t ist = h->MaterialIndex;
2132 // if (i < nhits - 1) {
2133 // for (Int_t j = istold - 1; j > ist; j--)
2134 // err = err || KF->vMaterial[j]->Pass( kftr, downstream, qp0 );
2135 // }
2136
2137 for (Int_t i = 0; i < nhits; i++) {
2138 h = kftr.GetHit( i );
2139 Int_t ist = h->MaterialIndex;
2140 if (i > 0) {
2141 for (Int_t j = istold + 1; j < ist; j++)
2142 err = err || KF->vMaterial[j]->Pass( kftr, downstream, qp0 );
2143 }
2144 err = err || h->Filter( kftr, downstream, qp0 );
2145 dchi2 = kftr.GetRefChi2() - chi20;
2146 cout << dchi2 << " " << kftr.GetRefChi2() << endl;
2147 //if (dchi2 <= gkChi2Cut) indok.insert(i);
2148 //if (dchi2 <= 2*gkChi2Cut) indok.insert(i);
2149 indok.insert(i);
2150 if (dchi2 > dchi2max) {
2151 dchi2max = dchi2;
2152 imax = i;
2153 }
2154 chi20 = kftr.GetRefChi2();
2155 istold = ist;
2156 }
2157 */
2158 // Remove last hit and refit
2159 set<Int_t> indok;
2160 for (Int_t i = 0; i < nhits; i++)
2161 indok.insert(i);
2162 int imax = nhits - 1;
2163
2164 if (indok.size() > 3) {
2165 // if (indok.size() >= fNhitsMin[fPass]) {
2166 // if (indok.size()-1 >= fNhitsMin[fPass]) {
2167 indok.erase(imax);
2168 // Good short track - refit
2169 indx = 0;
2170 hitcode = "-";
2171 for (set<Int_t>::iterator sit = indok.begin(); sit != indok.end(); ++sit) {
2172 hits[indx++] = hits[*sit];
2173 // hitcode += hits[*sit];
2174 hitcode += (to_string(hits[*sit]) + "-");
2175 }
2176 // if (fCandSet[sta0].count(hitcode)) continue;
2177 // fCandSet[sta0].insert(hitcode);
2178 // cout << " failed " << hitcode << endl;
2179 if (fCandSet[indok.size()].find(hitcode) != fCandSet[indok.size()].end())
2180 continue;
2181 fCandSet[indok.size()].insert(hitcode);
2182
2183 hits.Set(indok.size());
2184 // fitter.DoFit(&track);
2185 fitter.Fit(&track); // AZ
2186 int ndf = (track.GetFlag() == 0) ? track.GetNDF() : 1;
2187 // cout << " bbbbbb " << endl;
2188 // track.GetParamFirst()->Print();
2189 if (fVerbose > 0)
2190 cout << " bbbbbb " << track.GetChi2() << " " << hits.GetSize() << " " << track.GetChi2() / ndf
2191 << endl;
2192 if (track.GetChi2() / track.GetNDF() > gkChi2Cut)
2193 iok = 0;
2194 } else
2195 iok = 0;
2196
2197 } // if (track.GetChi2() / track.GetNDF() > gkChi2Cut)
2198
2199 // Good track - extra printout
2200 if (iok) {
2201 Double_t qual = track.GetNStsHits();
2202 // qual += (100.0 - TMath::Min(TMath::Abs(track.GetChi2()),100.0)) / 101.0;
2203 // if (track.GetNStsHits() <= 4) qual += TMath::Min(1/TMath::Abs(track.GetParamFirst()->GetQp()),10.0)
2204 // / 11.0;
2205 if (track.GetNStsHits() < 4) {
2206 // Vertex constraint
2207 /*
2208 TVector3 mom3, pos3;
2209 track.GetParamFirst()->Momentum(mom3);
2210 track.GetParamFirst()->Position(pos3);
2211 //Double_t ang = mom3.Angle(pos3);
2212 Double_t phi1 = TMath::ATan2 (track.GetParamFirst()->GetY(), track.GetParamFirst()->GetX());
2213 Double_t phi2 = TMath::ATan2 (track.GetParamFirst()->GetTy(), track.GetParamFirst()->GetTx());
2214 Double_t dphi = phi1 - Proxim(phi1,phi2);
2215 //qual += (2.1 - TMath::Min(TMath::Abs(dphi),2.1)) / 2.2;
2216 //Double_t www = TMath::Exp(-TMath::Abs(dphi)/2.1);
2217 Double_t www = TMath::Min (TMath::Exp(-TMath::Abs(dphi)/0.2), 0.999); // 0.2 - rms
2218 if (TMath::Abs(dphi) > 2.1) www = 0.0;
2219 */
2220 // Linear extrapolation to Zvert
2221 Double_t xextr; // = pos3.X() - pos3.Z() * track.GetParamFirst()->GetTx() + 0.04;
2222 Double_t yextr; // = pos3.Y() - pos3.Z() * track.GetParamFirst()->GetTy();
2223 // Extrapolation to Zvert
2224 FairTrackParam param = *track.GetParamFirst();
2225 fitter.Extrapolate(&param, fXyzv[2], &param);
2226 xextr = param.GetX() - fXyzv[0];
2227 yextr = param.GetY() - fXyzv[1];
2228 Double_t rad =
2229 TMath::Sqrt(xextr * xextr / 0.035 / 0.035 + yextr * yextr / 0.2 / 0.2); // 0.035, 0.2 - sigmas
2230 Double_t www = TMath::Min(TMath::Exp(-rad), 0.999);
2231 qual += TMath::Exp(-TMath::Abs(track.GetChi2())) * www;
2232 } else
2233 qual += (100.0 - TMath::Min(TMath::Abs(track.GetChi2()), 100.0)) / 101.0;
2234 fTracks.insert(pair<Double_t, CbmStsTrack>(-qual, track));
2235 // AZ-130722 pair<multimap<Int_t,Int_t>::iterator,multimap<Int_t,Int_t>::iterator> ret;
2236 pair<unordered_multimap<Int_t, Int_t>::iterator, unordered_multimap<Int_t, Int_t>::iterator>
2237 ret; // AZ-130722
2238 nhits = hits.GetSize();
2239 if (fVerbose > 0) {
2240 cout << " Good track: " << endl;
2241 // track.GetParamFirst()->Print();
2242 cout << track.GetChi2() << " " << nhits << " " << track.GetParamFirst()->GetZ() << endl;
2243
2244 for (Int_t j = nhits - 1; j >= 0; --j) {
2245 ret = fHit2id.equal_range(hits[j]);
2246 cout << " (";
2247 // AZ-130722 for (multimap<Int_t,Int_t>::iterator mit = ret.first; mit != ret.second; ++mit) {
2248 for (unordered_multimap<Int_t, Int_t>::iterator mit = ret.first; mit != ret.second; ++mit) {
2249 if (mit == ret.first)
2250 cout << mit->first << "*";
2251 else
2252 cout << ":";
2253 cout << mit->second;
2254 }
2255 cout << ")";
2256 }
2257 cout << "\n";
2258 }
2259 }
2260
2261 } // for (Int_t itra = 0;
2262 }
2263}
2264
2265// -------------------------------------------------------------------------
2266
2267void BmnStsVectorFinderV9::RemoveDoubles()
2268{
2269 // Remove double tracks (keep the ones with better quality)
2270
2271 multimap<Double_t, CbmStsTrack>::iterator mit = fTracks.begin(), mit1;
2272 Int_t nOK = fTracks.size(), nOut = fVectorArray->GetEntriesFast();
2273
2274 for (; mit != fTracks.end(); ++mit) {
2275 if (mit->second.GetFlag())
2276 continue; // discarded track
2277 mit1 = mit;
2278 ++mit1;
2279
2280 for (; mit1 != fTracks.end(); ++mit1) {
2281 if (mit1->second.GetFlag())
2282 continue; // discarded track
2283 if (!AreTracksDoubles(mit1->second, mit->second))
2284 continue;
2285 mit1->second.SetFlag(1); // double - discard it
2286 --nOK;
2287 /*
2288 // Debug printout
2289 TArrayI &hits = *mit1->second.GetStsHits();
2290 Int_t nhits = mit1->second.GetNStsHits();
2291 pair<multimap<Int_t,Int_t>::iterator,multimap<Int_t,Int_t>::iterator> ret;
2292 cout << " Double: " << mit->first << " " << mit1->first << endl;
2293 for (Int_t j = nhits-1; j >= 0; --j) {
2294 ret = fHit2id.equal_range(hits[j]);
2295 cout << " (";
2296 for (multimap<Int_t,Int_t>::iterator mit2 = ret.first; mit2 != ret.second; ++mit2) {
2297 if (mit2 == ret.first) cout << mit2->first << "*";
2298 else cout << ":";
2299 cout << mit2->second;
2300 }
2301 cout << ")";
2302 }
2303 cout << "\n";
2304 */
2305 }
2306
2307 new ((*fVectorArray)[nOut++]) CbmStsTrack(mit->second);
2308 // cout << " yyyyy " << mit->second.GetParamFirst()->GetZ() << " " << mit->second.GetNStsHits() << endl;
2309 }
2310 if (fVerbose > 0)
2311 cout << " RemoveDoubles: " << fTracks.size() << " " << nOK << endl;
2312}
2313
2314// -------------------------------------------------------------------------
2315
2316Bool_t BmnStsVectorFinderV9::AreTracksDoubles(CbmStsTrack& tr1, CbmStsTrack& tr2)
2317{
2320
2321 Int_t limCommonPoint = (tr1.GetNStsHits() + 1) / 2; // at least so many common hits should be found
2322 // limCommonPoint = TMath::Min (limCommonPoint, 2); // 2 common hits at max
2323 limCommonPoint = TMath::Min(limCommonPoint, 1); // 1 common hit at max
2324 TArrayI &hits1 = *tr1.GetStsHits(), &hits2 = *tr2.GetStsHits();
2325 Int_t nh1 = hits1.GetSize(), nh2 = hits2.GetSize(), nHitsCommon = 0, j = nh2 - 1;
2326 if (nh1 == 3)
2327 limCommonPoint = 0;
2328
2329 for (Int_t i = nh1 - 1; i >= 0; i--) {
2330 // for (Int_t i = 0; i <= nh1 - 1; ++i){
2331 CbmStsHit* hit1 = (CbmStsHit*)fHitArray->UncheckedAt(tr1.GetStsHitIndex(i));
2332
2333 for (; j >= 0; j--) {
2334 // for ( ; j <= nh2 - 1; ++j){
2335 CbmStsHit* hit2 = (CbmStsHit*)fHitArray->UncheckedAt(tr2.GetStsHitIndex(j));
2336
2337 // Is the hit common for two tracks?
2338 if (hit1 == hit2) {
2339 nHitsCommon++;
2340 // AZ-281121 if (nHitsCommon == limCommonPoint) return kTRUE; // already enough common hits
2341 if (nHitsCommon > limCommonPoint)
2342 return kTRUE; // already enough common hits
2343 break;
2344 }
2345
2346 if (hit2->GetStationNr() < hit1->GetStationNr())
2347 break; // already closer to target
2348 }
2349
2350 // AZ-281121 if ((nh1 - i - 1) + limCommonPoint - nHitsCommon > nh1) return kFALSE; // there'll be not enough
2351 // common hits already
2352 if (i + nHitsCommon <= limCommonPoint)
2353 return kFALSE; // there'll be not enough common hits already
2354 }
2355
2356 // AZ-281121 if (nHitsCommon < limCommonPoint) return kFALSE; // not too many common hits
2357 if (nHitsCommon <= limCommonPoint)
2358 return kFALSE; // not too many common hits
2359
2360 return kTRUE;
2361}
2362
2363// -------------------------------------------------------------------------
2364/*
2365void BmnStsVectorFinderV9::RemoveFakes()
2366{
2367 // Remove fake tracks (according to some empirical rules)
2368
2369 Int_t ncand = fVectorArray->GetEntriesFast();
2370
2371 for (Int_t j = 0; j < ncand; ++j) {
2372 CbmStsTrack *tr = (CbmStsTrack*) fVectorArray->UncheckedAt(j);
2373 Int_t nhits = tr->GetNStsHits();
2374 if (1/TMath::Abs(tr->GetParamFirst()->GetQp()) < 0.2) { fVectorArray->Remove(tr); continue; }
2375 if (nhits < 4) {
2376 CbmStsHit *hit = (CbmStsHit*) fHitArray->UncheckedAt(tr->GetStsHitIndex(0));
2377 if (hit->GetStationNr() <= 3) continue;
2378 fVectorArray->Remove(tr); // short track starts in GEM
2379 }
2380 }
2381 fVectorArray->Compress();
2382
2383 cout << " Remove fakes: " << ncand << " " << fVectorArray->GetEntriesFast() << endl;
2384}
2385*/
2386// -------------------------------------------------------------------------
2387void BmnStsVectorFinderV9::RemoveFakes()
2388{
2389 // Remove fake tracks (according to some empirical rules)
2390
2391 multimap<Double_t, CbmStsTrack>::iterator mit = fTracks.begin();
2392 Int_t nOK = fTracks.size(), nOut = nOK;
2393
2394 for (; mit != fTracks.end(); ++mit) {
2395 if (mit->second.GetFlag())
2396 continue; // discarded track
2397 CbmStsTrack* tr = &mit->second;
2398 Int_t nhits = tr->GetNStsHits();
2399 if (1 / TMath::Abs(tr->GetParamFirst()->GetQp()) < 0.2) {
2400 tr->SetFlag(1);
2401 --nOut;
2402 continue;
2403 }
2404 if (nhits < 4) {
2405 CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(tr->GetStsHitIndex(0));
2406 if (hit->GetStationNr() <= 3)
2407 continue;
2408 tr->SetFlag(1); // short track starts in GEM
2409 --nOut;
2410 }
2411 }
2412
2413 if (fVerbose > 0)
2414 cout << " Remove fakes: " << nOK << " " << nOut << endl;
2415}
2416
2417// -------------------------------------------------------------------------
2418
2419void BmnStsVectorFinderV9::ExcludeFakes()
2420{
2421 // Exclude fake tracks (with too many shared clusters)
2422
2423 Int_t ncand = fVectorArray->GetEntriesFast();
2424
2425 map<Int_t, set<Int_t>> mClusTr;
2426 map<Int_t, Double_t> mWeight;
2427
2428 for (Int_t j = 0; j < ncand; ++j) {
2429 CbmStsTrack* tr = (CbmStsTrack*)fVectorArray->UncheckedAt(j);
2430 Int_t nhits = tr->GetNStsHits();
2431
2432 for (Int_t jh = 0; jh < nhits; ++jh) {
2433 Int_t indx = tr->GetStsHitIndex(jh);
2434 CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(indx);
2435
2436 for (Int_t side = 0; side < 2; ++side) {
2437 Int_t iclus = hit->GetDigi(side);
2438 if (mClusTr.count(iclus) == 0) {
2439 set<Int_t> aaa;
2440 mClusTr[iclus] = aaa;
2441 }
2442 mClusTr[iclus].insert(j);
2443 }
2444 }
2445
2446 Double_t www = TMath::Exp(-tr->GetChi2() / tr->GetNDF());
2447 // Extrapolation to Zvert
2448 if (nhits == 3) {
2449 FairTrackParam param = *tr->GetParamFirst();
2450 fitter.Extrapolate(&param, fXyzv[2], &param);
2451 Double_t xextr = param.GetX() - fXyzv[0];
2452 Double_t yextr = param.GetY() - fXyzv[1];
2453 Double_t rad =
2454 TMath::Sqrt(xextr * xextr / 0.035 / 0.035 + yextr * yextr / 0.2 / 0.2); // 0.035, 0.2 - sigmas
2455 www *= TMath::Exp(-rad);
2456 }
2457 mWeight[j] = www;
2458 }
2459
2460 while (1) {
2461 multimap<Double_t, pair<Int_t, Int_t>> mapQual;
2462 set<Int_t> set2kill, set2kill1;
2463
2464 for (Int_t j = 0; j < ncand; ++j) {
2465 CbmStsTrack* tr = (CbmStsTrack*)fVectorArray->UncheckedAt(j);
2466 if (tr == NULL)
2467 continue; // removed track
2468 Int_t nhits = tr->GetNStsHits();
2469 Int_t nover = 0;
2470
2471 for (Int_t jh = 0; jh < nhits; ++jh) {
2472 Int_t indx = tr->GetStsHitIndex(jh);
2473 CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(indx);
2474 // if (mClusTr[hit->GetDigi(0)].size() > 1 || mClusTr[hit->GetDigi(1)].size() > 1) ++nover;
2475
2476 for (Int_t side = 0; side < 2; ++side) {
2477 Int_t iclus = hit->GetDigi(side);
2478 if (mClusTr[iclus].size() > 1)
2479 ++nover;
2480 }
2481 }
2482 // Double_t quality = nhits * 1000 + (1 - Double_t(nover) / nhits) * 100
2483 //+ (1 - TMath::Min(tr->GetChi2(),100.0) / 100);
2484 Double_t quality = nhits * 1000;
2485 Double_t www = TMath::Exp(-Double_t(nover) / nhits / 2) * mWeight[j];
2486 mapQual.insert(pair<Double_t, pair<Int_t, Int_t>>(quality + www, pair<Int_t, Int_t>(j, nover)));
2487 // if ((nover == nhits) || (nhits == 3 && nover == 2)) set2kill.insert(j);
2488 // if ((nhits <= 4 && nover == nhits) || (nhits == 3 && nover == 2)) set2kill.insert(j);
2489 // if (nhits <= 4 && Double_t(nover)/nhits/2 > 0.81) set2kill.insert(j);
2490 if (nhits <= 5 && Double_t(nover) / nhits / 2 > 0.6)
2491 set2kill.insert(j);
2492 else if (!fExact && nhits <= 3 && www < 1.e-8)
2493 set2kill1.insert(j); // 1e-8 - empirical value
2494 }
2495
2496 if (set2kill.size())
2497 set2kill1.clear(); // first kill in set2kill
2498 if (set2kill.size() || set2kill1.size()) {
2499 // Remove track and update cluster-to-track info
2500
2501 for (multimap<Double_t, pair<Int_t, Int_t>>::iterator mit = mapQual.begin(); mit != mapQual.end(); ++mit) {
2502 Int_t itr = mit->second.first;
2503 CbmStsTrack* tr = (CbmStsTrack*)fVectorArray->UncheckedAt(itr);
2504 Int_t nhits = tr->GetNStsHits(); //, nover = mit->second.second;
2505 if (set2kill.count(itr) == 0 && set2kill1.count(itr) == 0)
2506 continue;
2507
2508 for (Int_t jh = 0; jh < nhits; ++jh) {
2509 Int_t indx = tr->GetStsHitIndex(jh);
2510 CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(indx);
2511 if (mClusTr[hit->GetDigi(0)].count(itr) == 0 || mClusTr[hit->GetDigi(1)].count(itr) == 0)
2512 // AZ-300822 { cout << " !!! No track found: " << itr << endl; exit(0); }
2513 {
2514 if (fVerbose > 0)
2515 cout << " !!! No track found: " << itr << endl; /*exit(0);*/
2516 }
2517 mClusTr[hit->GetDigi(0)].erase(itr);
2518 mClusTr[hit->GetDigi(1)].erase(itr);
2519 }
2520 fVectorArray->Remove(tr);
2521 break;
2522 }
2523 } else
2524 break;
2525 } // while (1)
2526
2527 fVectorArray->Compress();
2528
2529 // AZ-280623 - kill short tracks (3 or 4 hits without consecutive triplet)
2530 int nKilled = 0;
2531 ncand = fVectorArray->GetEntriesFast();
2532 /* AZ-190125 - does not kill anything
2533 for (Int_t j = 0; j < ncand; ++j) {
2534 CbmStsTrack *tr = (CbmStsTrack*) fVectorArray->UncheckedAt(j);
2535 Int_t nhits = tr->GetNStsHits();
2536 // AZ-190125 if (nhits > 4) continue;
2537 if (nhits > 44) continue; // AZ-190125
2538 // AZ-190125 int sta0 = -1, ncons = 0;
2539 int sta0 = -1, ncons = 1; // AZ-190125
2540
2541 for (Int_t jh = 0; jh < nhits; ++jh) {
2542 Int_t indx = tr->GetStsHitIndex(jh);
2543 CbmStsHit *hit = (CbmStsHit*) fHitArray->UncheckedAt(indx);
2544 int sta = hit->GetStationNr();
2545 if (jh == 0) sta0 = hit->GetStationNr();
2546 if (sta - sta0 < 2) {
2547 //fVectorArray->RemoveAt(j);
2548 ++ncons;
2549 } else ncons = 1;
2550 if (ncons >= 3) break;
2551 sta0 = sta;
2552 }
2553 if (ncons < 3) { fVectorArray->RemoveAt(j); ++nKilled; }
2554 }
2555 fVectorArray->Compress();
2556 cout << " !!! Killed w/out triplet: " << nKilled << endl;
2557 */
2558 // AZ-030723 - check short tracks for missing hits (remove if nmissMax hits in the acceptance are missing)
2559 ncand = fVectorArray->GetEntriesFast();
2560 nKilled = 0;
2561
2562 for (Int_t j = 0; j < ncand; ++j) {
2563 CbmStsTrack* tr = (CbmStsTrack*)fVectorArray->UncheckedAt(j);
2564 Int_t nhits = tr->GetNStsHits();
2565 // if (nhits > 4) continue;
2566 if (nhits > -4)
2567 continue; // AZ-200125 - do not kill
2568 // if (nhits > 44) continue;
2569 set<int> setStat;
2570
2571 for (int ih = 0; ih < nhits; ++ih) {
2572 int indx = tr->GetStsHitIndex(ih);
2573 CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(indx);
2574 setStat.insert(hit->GetStationNr());
2575 }
2576
2577 CbmKFTrack kftr = CbmKFTrack(*tr);
2578 fitter.SetKFHits(kftr, tr);
2579 FairTrackParam param;
2580
2581 int firstSta = *setStat.begin(), nmisSi = 0, nmisGem = 0, nmiss = 0;
2582
2583 // Go upstream
2584
2585 for (int ista = firstSta - 2; ista >= 0; --ista) {
2586 nmiss = 0;
2587 if (ista >= fNSi) {
2588 BmnGemStripStation* stat = fGemStationSet->GetStation(ista - fNSi);
2589 kftr.Smooth(stat->GetZPosition());
2590 kftr.GetTrackParam(param);
2591 // fitter.Extrapolate (&param, stat->GetZPosition(), &param);
2592 // if (param.GetX() < stat->GetXMinStation() || param.GetX() > stat->GetXMaxStation()) continue;
2593 // if (param.GetY() < stat->GetYMinStation() || param.GetY() > stat->GetYMaxStation()) continue;
2594
2595 for (int imod = 0; imod < stat->GetNModules(); ++imod) {
2596 BmnGemStripModule* mod = stat->GetModule(imod);
2597 if (mod->IsPointInsideModule(param.GetX(), param.GetY())) {
2598 ++nmiss;
2599 ++nmisGem;
2600 // cout << " +++ " << param.GetX() << " " << param.GetY() << endl;
2601 break;
2602 } // else cout << " --- " << param.GetX() << " " << param.GetY() << endl;
2603 }
2604 } else {
2605 BmnSiliconStation* stat = fSilStationSet->GetStation(ista);
2606 kftr.Smooth(stat->GetZPosition());
2607 kftr.GetTrackParam(param);
2608 // fitter.Extrapolate (&param, stat->GetZPosition(), &param);
2609 // if (param.GetX() < stat->GetXMinStation() || param.GetX() > stat->GetXMaxStation()) continue;
2610 // if (param.GetY() < stat->GetYMinStation() || param.GetY() > stat->GetYMaxStation()) continue;
2611
2612 for (int imod = 0; imod < stat->GetNModules(); ++imod) {
2613 BmnSiliconModule* mod = stat->GetModule(imod);
2614 // AZ-180125 if (mod->IsPointInsideModule(param.GetX(),param.GetY())) {
2615 if (mod->IsPointInsideModule(param.GetX(), param.GetY(), kFALSE)) { // AZ-180125
2616 ++nmiss;
2617 ++nmisSi;
2618 // cout << " +++ " << param.GetX() << " " << param.GetY() << endl;
2619 break;
2620 } // else cout << " --- " << param.GetX() << " " << param.GetY() << endl;
2621 }
2622 }
2623 if (nmiss == 0)
2624 break; // no need for checking other stations
2625 } // for (int ista = firstSta - 2;
2626
2627 // Check downstream stations
2628 nmiss = 0;
2629
2630 for (int ista = firstSta; ista <= fNsta; ++ista) {
2631 if (setStat.find(ista + 1) != setStat.end())
2632 continue;
2633 nmiss = 0;
2634 if (ista >= fNSi) {
2635 BmnGemStripStation* stat = fGemStationSet->GetStation(ista - fNSi);
2636 kftr.Smooth(stat->GetZPosition());
2637 kftr.GetTrackParam(param);
2638 // fitter.Extrapolate (&param, stat->GetZPosition(), &param);
2639 // if (param.GetX() < stat->GetXMinStation() || param.GetX() > stat->GetXMaxStation()) continue;
2640 // if (param.GetY() < stat->GetYMinStation() || param.GetY() > stat->GetYMaxStation()) continue;
2641
2642 for (int imod = 0; imod < stat->GetNModules(); ++imod) {
2643 BmnGemStripModule* mod = stat->GetModule(imod);
2644 if (mod->IsPointInsideModule(param.GetX(), param.GetY())) {
2645 ++nmiss;
2646 ++nmisGem;
2647 // cout << " +++ " << param.GetX() << " " << param.GetY() << endl;
2648 break;
2649 } // else cout << " --- " << param.GetX() << " " << param.GetY() << endl;
2650 }
2651 } else {
2652 BmnSiliconStation* stat = fSilStationSet->GetStation(ista);
2653 kftr.Smooth(stat->GetZPosition());
2654 kftr.GetTrackParam(param);
2655 // fitter.Extrapolate (&param, stat->GetZPosition(), &param);
2656 // if (param.GetX() < stat->GetXMinStation() || param.GetX() > stat->GetXMaxStation()) continue;
2657 // if (param.GetY() < stat->GetYMinStation() || param.GetY() > stat->GetYMaxStation()) continue;
2658
2659 for (int imod = 0; imod < stat->GetNModules(); ++imod) {
2660 BmnSiliconModule* mod = stat->GetModule(imod);
2661 // AZ-180125 if (mod->IsPointInsideModule(param.GetX(),param.GetY())) {
2662 if (mod->IsPointInsideModule(param.GetX(), param.GetY(), kFALSE)) { // AZ-180125
2663 ++nmiss;
2664 ++nmisSi;
2665 // cout << " +++ " << param.GetX() << " " << param.GetY() << endl;
2666 break;
2667 } // else cout << " --- " << param.GetX() << " " << param.GetY() << endl;
2668 }
2669 }
2670 if (nmiss == 0)
2671 break; // no need for checking other stations
2672 } // for (int ista = firstSta; ista <= fNsta;
2673
2674 if (nmisSi >= 2) {
2675 fVectorArray->RemoveAt(j);
2676 ++nKilled;
2677 continue;
2678 }
2679 tr->SetFlag(nmisSi * 10 + nmisGem);
2680
2681 } // for (Int_t j = 0; j < ncand;
2682 // cout << " !!! Killed due to missing: " << nKilled << endl;
2683 fVectorArray->Compress();
2684
2685 if (fVerbose > 0)
2686 cout << " Exclude fakes: " << ncand << " " << fVectorArray->GetEntriesFast() << endl;
2687}
2688
2689// -------------------------------------------------------------------------
2690
2691Double_t BmnStsVectorFinderV9::DxVsMom(Int_t ista, candvec& vec)
2692{
2693 // Compute momentum correction for parabolic extrapolation
2694
2695 Double_t pars[19][4] = {{0, 0, 0, 0},
2696 {0, 0, 0, 0},
2697 {8.28698e-06, 13.3515, 0, 0},
2698 {0.031189, 4.56933, 0.031189, 4.56933},
2699 {0.395468, 3.77954, 0.190848, 3.41027},
2700 {0.128467, 4.26652, 0.128467, 4.26652},
2701 {0.160979, 4.95095, 0.160979, 4.95095}};
2702
2703 Double_t pxz = vec.momxz;
2704 Int_t endsta = vec.stahit.rbegin()->first;
2705 Double_t p0 = pars[ista][0];
2706 Double_t p1 = pars[ista][1];
2707 if (endsta > ista + 2) {
2708 p0 = pars[ista][2];
2709 p1 = pars[ista][3];
2710 }
2711
2712 Double_t dx = p0 / TMath::Power(TMath::Log10(10 * pxz), p1);
2713 return dx;
2714}
2715
2716// -------------------------------------------------------------------------
2717
2718Double_t BmnStsVectorFinderV9::Proxim(Double_t phi0, Double_t phi)
2719{
2721
2722 Double_t dPhi = phi0 - phi;
2723 if (TMath::Abs(dPhi) > TMath::Pi())
2724 phi += TMath::Pi() * 2 * TMath::Sign(1., dPhi);
2725 return phi;
2726}
2727
2728// -------------------------------------------------------------------------
2729// AZ-280722 Double_t BmnStsVectorFinderV9::LinearFit(CbmStsTrack *tr, Int_t newtr, Double_t& ty)
2730Double_t BmnStsVectorFinderV9::LinearFit(candvec& cand,
2731 candvec& cand2,
2732 CbmStsTrack* tr,
2733 Int_t newtr,
2734 Double_t& ty) // AZ-280722
2735{
2737
2738 // static Double_t x[3], y[3], z[3], xs[3], ys[3], xys[3], x2s[3];
2739 static Double_t /*x[3],*/ y[3], /*z[3],*/ l[3] = {0}, lens[3] = {0}, ys[3], lys[3], l2s[3];
2740
2741 // AZ-150722 Int_t i0 = (newtr == 1) ? 0 : 2;
2742 Int_t i0 = (newtr == 0) ? 0 : 2;
2743
2744 for (Int_t i = i0; i < 3; ++i) {
2745 CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(tr->GetStsHitIndex(i));
2746 // AZ-280722 x[i] = hit->GetX();
2747 // AZ-280722 z[i] = hit->GetZ();
2748 if (i) {
2749 /*AZ-280722
2750 Double_t dx = x[i] - x[i-1];
2751 Double_t dz = z[i] - z[i-1];
2752 l[i] = lens[i] = TMath::Sqrt (dx * dx + dz * dz) + l[i-1];
2753 */
2754 if (i == 1)
2755 l[i] = lens[i] = cand.lengxz; // AZ-280722
2756 else
2757 l[i] = lens[i] = cand2.lengxz + l[i - 1]; // AZ-280722
2758 }
2759 y[i] = ys[i] = hit->GetY();
2760 // cout << " ily " << i << " " << l[i] << " " << y[i] << endl;
2761 }
2762
2763 for (Int_t i = i0; i < 3; ++i) {
2764 lys[i] = l[i] * y[i];
2765 l2s[i] = l[i] * l[i];
2766 if (i) {
2767 lens[i] += lens[i - 1];
2768 ys[i] += ys[i - 1];
2769 lys[i] += lys[i - 1];
2770 l2s[i] += l2s[i - 1];
2771 }
2772 }
2773 Double_t b = 3 * lys[2] - lens[2] * ys[2];
2774 b /= (3 * l2s[2] - lens[2] * lens[2]);
2775 Double_t a = (ys[2] - b * lens[2]) / 3;
2776
2777 Double_t chi2 = 0.0;
2778
2779 for (int i = 0; i < 3; ++i) {
2780 Double_t dy = a + b * l[i] - y[i];
2781 chi2 += (dy * dy);
2782 }
2783 // Check
2784 // TGraph gr(3,l,y);
2785 // gr.Fit("pol1");
2786 // cout << " LS fit: " << a << " " << b << " " << chi2 << endl;
2787 ty = b;
2788 return chi2;
2789}
2790
2791// -------------------------------------------------------------------------
2792
2793Double_t BmnStsVectorFinderV9::Curv3(candvec& cand1, candvec& cand2, candvec& cand3, int newtr)
2794{
2795 // Compute triplet curvature in ZX plane
2796
2797 static TVector3 points3[3], midPoint;
2798 int indx = 1;
2799 static Double_t by = 0.0;
2800 // static int first = 99;
2801
2802 // if (first != fPass) {
2803 // first = fPass;
2804 // AZ-270625 if (first) {
2805 // if (cand3.stahit.size() <= 2 && first) { // AZ-270625 - for doublets
2806 if (cand3.stahit.size() <= 2) { // AZ-270625 - for doublets
2807 // first = 0;
2808 points3[0].SetXYZ(fXyzv[0], fXyzv[1], fXyzv[2]);
2809 }
2810
2811 map<int, int>::iterator mit = cand3.stahit.begin();
2812 if (0) { // fPass > 1) { //0 - assuming tracks coming from (0,0,0)
2813 // if (newtr == 0 && cand3.stahit.size() > 2) {
2814 points3[0] = fmapHits[mit->second].xyz; // AZ-280722
2815 points3[0].SetY(0.0); // AZ-280722
2816 }
2817 // AZ-310722 ++mit;
2818 // AZ-270625 if (cand3.stahit.size() > 2) ++mit; //AZ-310722
2819 if (cand3.stahit.size() > 2) {
2820 indx = 0;
2821 // first = 99;
2822 } // AZ-270625
2823
2824 for (; mit != cand3.stahit.end(); ++mit) {
2825 // if (newtr > 0 && indx < 2) { ++indx; continue; }
2826 points3[indx] = fmapHits[mit->second].xyz;
2827 if (indx == 1)
2828 midPoint = points3[indx];
2829 // AZ-270625 points3[indx] -= points3[0];
2830 if (indx)
2831 points3[indx] -= points3[0]; // AZ-270625
2832 points3[indx++].SetY(0.0);
2833 }
2834
2835 TVector3 vec21 = points3[2] - points3[1];
2836 // AZ-280722 Double_t cosAlpha = points3[2] * vec21 / points3[2].Mag() / vec21.Mag();
2837 // Double_t cosAlpha = points3[1] * vec21 / points3[1].Mag() / vec21.Mag(); //AZ-280722
2838 Double_t cosAlpha = points3[1] * vec21 / points3[1].Mag() / cand2.lengxz; // AZ-280722
2839 // AZ-280722 Double_t rad = points3[1].Mag() / 2. / TMath::Sin(TMath::ACos(cosAlpha));
2840 Double_t rad = points3[2].Mag() / 2. / TMath::Sin(TMath::ACos(cosAlpha)); // AZ-280722
2841 // Double_t bz = FairRunAna::Instance()->GetField()->GetBz(0.,0.,0.);
2842 // Double_t factor = 0.003 * bz / 10.; // 0.3 * 0.01 * 5kG / 10
2843 // AZ-310722 Double_t sign = TMath::Sign (1.0, points3[1].Cross(points3[2]).Y());
2844 Double_t sign = TMath::Sign(3e-4, points3[1].Cross(points3[2]).Y()); // AZ-310722
2845 if (newtr == 0) {
2846 FairField* magField = FairRunAna::Instance()->GetField();
2847 by = magField->GetBy(midPoint[0], midPoint[1], midPoint[2]);
2848 if (TMath::Abs(by) < 0.1)
2849 by = -10.0; // AZ-261222 - zero magnetic field
2850 }
2851 // AZ-270625 Double_t pxz = by * rad * sign;
2852 Double_t pxz = -by * rad * sign;
2853 // cout << " pt " << pxz << endl;
2854 // AZ-170722 return pxz;
2855 return 1 / pxz; // AZ-170722
2856}
2857// -------------------------------------------------------------------------
2858
2859set<int> BmnStsVectorFinderV9::KalmanWindow(candvec& cand, int hitIndx)
2860{
2861 // Find hits around extrapolated track
2862
2863 CbmStsTrack track;
2864 string hitcode;
2865 MakeStsTrack(cand, hitcode, track);
2866 fitter.DoFit(&track);
2867 FairTrackParam param;
2868 CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(hitIndx);
2869 Double_t z = hit->GetZ();
2870 fitter.Extrapolate(&track, z, &param);
2871 Double_t sigx = 5 * TMath::Sqrt(param.GetCovariance(0, 0));
2872 Double_t sigy = 5 * TMath::Sqrt(param.GetCovariance(1, 1));
2873 if (fVerbose > 0)
2874 cout << " Kalman: " << z << " " << hit->GetStationNr() << " " << sigx << " " << sigy << endl;
2875
2876 // Get hits on the downstream station
2877 int istanext = hit->GetStationNr() - 1;
2878 multimap<Double_t, Int_t>::iterator mityb = fmapY[istanext].lower_bound(param.GetY() - sigy);
2879 multimap<Double_t, Int_t>::iterator mitye = fmapY[istanext].upper_bound(param.GetY() + sigy);
2880 multimap<Double_t, Int_t>::iterator mitxb = fmapX[istanext].lower_bound(param.GetX() - sigx);
2881 multimap<Double_t, Int_t>::iterator mitxe = fmapX[istanext].upper_bound(param.GetX() + sigx);
2882 // Get hits from the acceptance window
2883 set<Int_t> setX, setY, intersect;
2884 for (multimap<Double_t, Int_t>::iterator mit = mitxb; mit != mitxe; ++mit)
2885 if (fmapHits[mit->second].used == 0)
2886 setX.insert(mit->second);
2887 for (multimap<Double_t, Int_t>::iterator mit = mityb; mit != mitye; ++mit)
2888 if (fmapHits[mit->second].used == 0)
2889 setY.insert(mit->second);
2890 set_intersection(setX.begin(), setX.end(), setY.begin(), setY.end(), std::inserter(intersect, intersect.begin()));
2891 return intersect;
2892}
2893
2894// -------------------------------------------------------------------------
2895
2896std::string BmnStsVectorFinderV9::MakeCode(candvec& cand)
2897{
2898 // Make track code
2899
2900 string code("-");
2901 for (map<int, int>::iterator it = cand.stahit.begin(); it != cand.stahit.end(); ++it)
2902 code += (to_string(it->second) + "-");
2903 return code;
2904}
2905
2906// -------------------------------------------------------------------------
2907
2908void BmnStsVectorFinderV9::PrintHits(candvec& cand)
2909{
2910 // Print hits
2911
2912 for (map<int, int>::iterator mit = cand.stahit.begin(); mit != cand.stahit.end(); ++mit) {
2913 int ista = mit->first, ih = mit->second;
2914 CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(ih);
2915 // cout << mit->first << "(" << hit->GetZ() << "):" << mit->second << " - ";
2916 printf("%d(%5.2f,%5.2f,%5.2f):%d", ista, hit->GetZ(), fmapHits[ih].tx, fmapHits[ih].ty, ih);
2917 auto ret = fHit2id.equal_range(ih);
2918 for (unordered_multimap<Int_t, Int_t>::iterator mit1 = ret.first; mit1 != ret.second; ++mit1) {
2919 if (mit1 == ret.first)
2920 cout << "*";
2921 else
2922 cout << ":";
2923 cout << mit1->second;
2924 }
2925 cout << " - ";
2926 }
2927 cout << endl;
2928}
2929
2930// -------------------------------------------------------------------------
2931
2932void BmnStsVectorFinderV9::InitTMVA3()
2933{
2934 // Initialize TMVA package for triplets - AZ-300922
2935
2936 // This loads the library
2937 // TMVA::Tools::Instance();
2938
2939 // Default MVA methods to be trained + tested
2940 std::map<std::string, int> Use;
2941
2942 // Use["MLPBNN"] = 1; // Recommended ANN with BFGS training method and bayesian regulator
2943 Use["BDTD"] = 1; // decorrelation + Adaptive Boost
2944 Float_t* var = fVarTMVA;
2945
2946 fReadersTMVA3.insert(pair<int, TMVA::Reader*>(-13, nullptr)); // pass3 for negative
2947 fReadersTMVA3.insert(pair<int, TMVA::Reader*>(-11, nullptr));
2948 fReadersTMVA3.insert(pair<int, TMVA::Reader*>(-10, nullptr));
2949 fReadersTMVA3.insert(pair<int, TMVA::Reader*>(10, nullptr));
2950 fReadersTMVA3.insert(pair<int, TMVA::Reader*>(11, nullptr));
2951 fReadersTMVA3.insert(pair<int, TMVA::Reader*>(13, nullptr));
2952
2953 // --- Create the Reader objects
2954
2955 for (auto mit = fReadersTMVA3.begin(); mit != fReadersTMVA3.end(); ++mit) {
2956 TMVA::Reader* reader = new TMVA::Reader("!Color:!Silent");
2957 fReadersTMVA3[mit->first] = reader;
2958
2959 // Create a set of variables and declare them to the reader
2960 // The variable names MUST corresponds in name and type to those given in the weight file(s) used
2961
2962 // FIXME
2963 int iii = 0;
2964
2965 reader->AddVariable("tanx", &var[iii++]);
2966 reader->AddVariable("tany", &var[iii++]);
2967 reader->AddVariable("zb", &var[iii++]);
2968 reader->AddVariable("dty", &var[iii++]);
2969 // reader->AddVariable( "1/pxz1", &var[iii++] );
2970 // reader->AddVariable( "log(abs(pxz1))", &var[iii++] );
2971 reader->AddVariable("log(abs(pxz2))", &var[iii++]);
2972 reader->AddVariable("dpxz", &var[iii++]);
2973 reader->AddVariable("corrq[0]", &var[iii++]);
2974 reader->AddVariable("corrq[1]", &var[iii++]);
2975 reader->AddVariable("corrq[2]", &var[iii++]);
2976 // reader->AddVariable( "log(mult+rndm)", &var[iii++] );
2977 reader->AddVariable("log(mult)", &var[iii++]);
2978
2979 // --- Book the MVA methods
2980
2981 TString dir = "datasetTMVA3/weights";
2982 dir += mit->first;
2983 dir += "/";
2984 TString prefix = "TMVAClassification";
2985
2986 // Book method(s)
2987 for (std::map<std::string, int>::iterator it = Use.begin(); it != Use.end(); it++) {
2988 if (it->second) {
2989 TString methodName = TString(it->first) + TString(" method");
2990 TString weightfile = dir + prefix + TString("_") + TString(it->first) + TString(".weights.xml");
2991 reader->BookMVA(methodName, weightfile);
2992 }
2993 }
2994 } // for (auto mit = fReadersTMVA3.begin();
2995}
2996
2997// -------------------------------------------------------------------------
2998
2999Float_t BmnStsVectorFinderV9::TMVAOutput3(candvec& aaa1, candvec& aaa2, int mult)
3000{
3001 // Compute TMVA output for triplets: aaa1 - triplet, aaa2 - first doublet
3002
3003 Float_t* var = fVarTMVA;
3004 map<Int_t, Int_t>& hitMap = aaa1.stahit;
3005 CbmStsHit* hit[3] = {nullptr};
3006 int iii = 0, indx = 0, inds[3];
3007 inds[0] = inds[1] = inds[2] = 0;
3008
3010 for (map<Int_t, Int_t>::iterator mit = hitMap.begin(); mit != hitMap.end(); ++mit) {
3011 inds[indx] = mit->second;
3012 hit[indx++] = (CbmStsHit*)fHitArray->UncheckedAt(mit->second);
3013 }
3014
3015 Double_t dx = hit[1]->GetX() - hit[0]->GetX();
3016 Double_t dz = hit[1]->GetZ() - hit[0]->GetZ();
3017 Double_t lenxz1 = TMath::Sqrt(dx * dx + dz * dz);
3018 Double_t dx1 = hit[2]->GetX() - hit[1]->GetX();
3019 Double_t dz1 = hit[2]->GetZ() - hit[1]->GetZ();
3020 Double_t lenxz2 = TMath::Sqrt(dx1 * dx1 + dz1 * dz1);
3021 Double_t dty = (hit[2]->GetY() - hit[1]->GetY()) / lenxz2 - (hit[1]->GetY() - hit[0]->GetY()) / lenxz1;
3022 Double_t pxz1 = aaa1.momxz; // triplet curvature
3023 Double_t dpxz = (aaa2.momxz - pxz1) / (TMath::Abs(aaa2.momxz) + TMath::Abs(pxz1));
3024
3025 var[iii++] = hit[0]->GetX() / hit[0]->GetZ();
3026 var[iii++] = hit[0]->GetY() / hit[0]->GetZ();
3027 var[iii++] = hit[0]->GetZ();
3028 var[iii++] = dty;
3029 // var[iii++] = pxz1;
3030 var[iii++] = TMath::Log(TMath::Abs(1 / pxz1));
3031 var[iii++] = dpxz;
3032 var[iii++] = hit[0]->GetSignalDiv();
3033 var[iii++] = hit[1]->GetSignalDiv();
3034 var[iii++] = hit[2]->GetSignalDiv();
3035 var[iii++] = TMath::Log(mult);
3036
3037 int iread = (fPass == 2) ? 1 : fPass; // fPass==2 has the same triplet cuts as fPass==1
3038 iread += 10;
3039 if (pxz1 < 0)
3040 iread = -iread;
3041 Float_t tmvaout = fReadersTMVA3[iread]->EvaluateMVA("BDTD method");
3042 /*
3043 if (fOutf) {
3044 for (int j = 0; j < iii; ++j) *fOutf << var[j] << " ";
3045 *fOutf << mult << " " << tmvaout << " " << iread << " " << inds[0] << " " << inds[1] << " "
3046 << inds[2] << " " << fEvNo << endl;
3047 }
3048 */
3049 return tmvaout;
3050}
3051
3052// -------------------------------------------------------------------------
3053
3054void BmnStsVectorFinderV9::InitTMVA2()
3055{
3056 // Initialize TMVA package for doublets - AZ-221022
3057
3058 // This loads the library
3059 // TMVA::Tools::Instance();
3060
3061 // Default MVA methods to be trained + tested
3062 std::map<std::string, int> Use;
3063
3064 Use["MLPBNN"] = 1; // Recommended ANN with BFGS training method and bayesian regulator
3065 Use["BDTD"] = 1; // decorrelation + Adaptive Boost
3066
3067 // --- Create the Reader object
3068
3069 TMVA::Reader* reader = new TMVA::Reader("!Color:!Silent");
3070 fReaderTMVA2 = reader;
3071
3072 // Create a set of variables and declare them to the reader
3073 // The variable names MUST corresponds in name and type to those given in the weight file(s) used
3074
3075 // FIXME
3076 Float_t* var = fVarTMVA;
3077 int iii = 0;
3078 reader->AddVariable("xb", &var[iii++]);
3079 reader->AddVariable("yb", &var[iii++]);
3080 reader->AddVariable("zb", &var[iii++]);
3081 reader->AddVariable("tx", &var[iii++]);
3082 reader->AddVariable("ty", &var[iii++]);
3083 // reader->AddVariable( "dphixsc", &var[iii++] );
3084 reader->AddVariable("dtanysc", &var[iii++]);
3085 reader->AddVariable("lenxz", &var[iii++]);
3086 reader->AddVariable("1/pxz", &var[iii++]);
3087 reader->AddVariable("corrq[1]-corrq[0]", &var[iii++]);
3088 reader->AddVariable("corrn[1]-corrn[0]", &var[iii++]);
3089 reader->AddVariable("pass", &var[iii++]);
3090
3091 // --- Book the MVA methods
3092
3093 TString dir = "dataset4diff/weights2/";
3094 TString prefix = "TMVAClassification";
3095
3096 // Book method(s)
3097 for (std::map<std::string, int>::iterator it = Use.begin(); it != Use.end(); it++) {
3098 if (it->second) {
3099 TString methodName = TString(it->first) + TString(" method");
3100 TString weightfile = dir + prefix + TString("_") + TString(it->first) + TString(".weights.xml");
3101 reader->BookMVA(methodName, weightfile);
3102 }
3103 }
3104}
3105
3106// -------------------------------------------------------------------------
3107
3108Float_t BmnStsVectorFinderV9::TMVAOutput2(candvec& aaa1)
3109{
3110 // Compute TMVA output for doublets
3111
3112 Float_t* var = fVarTMVA;
3113 map<Int_t, Int_t>& hitMap = aaa1.stahit;
3114 CbmStsHit *hit = nullptr, *hit1 = nullptr;
3115 int iii = 0, indx = 0;
3116
3118 for (map<Int_t, Int_t>::iterator mit = hitMap.begin(); mit != hitMap.end(); ++mit) {
3119 hit = (CbmStsHit*)fHitArray->UncheckedAt(mit->second);
3120 if (indx == 0) {
3121 var[iii++] = hit->GetX();
3122 var[iii++] = hit->GetY();
3123 var[iii++] = hit->GetZ();
3124 hit1 = hit;
3125 } else if (indx == 1) {
3126 Double_t dx = hit->GetX() - var[0];
3127 Double_t dy = hit->GetY() - var[1];
3128 Double_t dz = hit->GetZ() - var[2];
3129 var[iii++] = dx / dz;
3130 var[iii++] = dy / dz;
3131 var[iii++] = aaa1.ty;
3132 var[iii++] = aaa1.lengxz;
3133 var[iii++] = aaa1.momxz;
3134 var[iii++] = hit->GetSignalDiv() - hit1->GetSignalDiv();
3135 var[iii++] = hit->GetTimeStamp() - hit1->GetTimeStamp();
3136 var[iii++] = (fPass == 2) ? 1 : fPass; // fPass==2 has the same triplet cuts as fPass==1
3137 }
3138 ++indx;
3139 }
3140 return fReaderTMVA2->EvaluateMVA("MLPBNN method");
3141}
3142
3143// -------------------------------------------------------------------------
3144
3145void BmnStsVectorFinderV9::InitTMVAtracks()
3146{
3147 // Initialize TMVA package for rejecting fake tracks (4- and 5-hits)
3148
3149 // This loads the library
3150 // TMVA::Tools::Instance();
3151
3152 // Default MVA methods to be trained + tested
3153 std::map<std::string, int> Use;
3154
3155 Use["BDTD"] = 1; // decorrelation + Adaptive Boost
3156 Float_t* var = fVarTMVA;
3157 fTmvaReaders.insert(pair<int, TMVA::Reader*>(-5, nullptr));
3158 fTmvaReaders.insert(pair<int, TMVA::Reader*>(-4, nullptr));
3159 fTmvaReaders.insert(pair<int, TMVA::Reader*>(4, nullptr));
3160 fTmvaReaders.insert(pair<int, TMVA::Reader*>(5, nullptr));
3161
3162 // --- Create the Reader objects
3163
3164 for (auto mit = fTmvaReaders.begin(); mit != fTmvaReaders.end(); ++mit) {
3165 TMVA::Reader* reader = new TMVA::Reader("!Color:!Silent");
3166 int nhits = mit->first;
3167 fTmvaReaders[nhits] = reader;
3168
3169 // Create a set of variables and declare them to the reader
3170 // The variable names MUST corresponds in name and type to those given in the weight file(s) used
3171
3172 // FIXME
3173 int iii = 0;
3174 // reader->AddVariable( "sqrt(xb*xb+yb*yb)", &var[iii++] );
3175 reader->AddVariable("atan2(yb,xb)", &var[iii++]);
3176 // reader->AddVariable( "zb", &var[iii++] );
3177 reader->AddVariable("tx", &var[iii++]);
3178 // reader->AddVariable( "ty", &var[iii++] );
3179 reader->AddVariable("log(abs(momr))", &var[iii++]);
3180 reader->AddVariable("chi2", &var[iii++]);
3181 // reader->AddVariable( "log(sqrt(xextr*xextr+yextr*yextr))", &var[iii++] );
3182 // reader->AddVariable( "atan2(yextr,xextr)", &var[iii++] );
3183 reader->AddVariable("corrq[0]", &var[iii++]);
3184 reader->AddVariable("corrq[1]", &var[iii++]);
3185 reader->AddVariable("corrq[2]", &var[iii++]);
3186 reader->AddVariable("corrq[3]", &var[iii++]);
3187 if (TMath::Abs(nhits) == 5)
3188 reader->AddVariable("corrq[4]", &var[iii++]);
3189 reader->AddVariable("corrq3", &var[iii++]);
3190 // reader->AddVariable( "nmiss+rndm", &var[iii++] );
3191 reader->AddVariable("nseqmx", &var[iii++]);
3192 // reader->AddVariable( "leng", &var[iii++] );
3193 reader->AddVariable("corrn[0]", &var[iii++]);
3194 reader->AddVariable("corrn[1]", &var[iii++]);
3195 reader->AddVariable("corrn[2]", &var[iii++]);
3196 reader->AddVariable("corrn[3]", &var[iii++]);
3197 if (TMath::Abs(nhits) == 5)
3198 reader->AddVariable("corrn[4]", &var[iii++]);
3199 reader->AddVariable("nsigem", &var[iii++]);
3200
3201 // --- Book the MVA methods
3202
3203 TString dir = gSystem->Getenv("VMCWORKDIR");
3204 // TString dir = "datasetTMVA1/weights";
3205 dir += "/parameters/reco/datasetTMVAV8/weights";
3206 // dir += nhits;
3207 dir += (nhits < 10) ? TMath::Abs(nhits) : nhits / 10;
3208 if (nhits < 0)
3209 dir += "neg";
3210 else
3211 dir += "pos";
3212 dir += "/";
3213 TString prefix = "TMVAClassification";
3214
3215 // Book method(s)
3216 for (std::map<std::string, int>::iterator it = Use.begin(); it != Use.end(); it++) {
3217 if (it->second) {
3218 TString methodName = TString(it->first) + TString(" method");
3219 TString weightfile = dir + prefix + TString("_") + TString(it->first) + TString(".weights.xml");
3220 reader->BookMVA(methodName, weightfile);
3221 }
3222 }
3223 } // for (auto mit = fTmvaReaders.begin();
3224}
3225
3226// -------------------------------------------------------------------------
3227
3228void BmnStsVectorFinderV9::TMVAOutputTra()
3229{
3230 // Compute TMVA output for tracks
3231
3232 // Fill cluster overlap map
3233 map<Int_t, set<Int_t>> mClusTr;
3234 int nTracks = fVectorArray->GetEntriesFast();
3235
3236 for (int itr = 0; itr < nTracks; ++itr) {
3237 CbmStsTrack* tr = (CbmStsTrack*)fVectorArray->UncheckedAt(itr);
3238 int nhit = tr->GetNStsHits();
3239 if (nhit < 4)
3240 continue; // exclude 3-hit tracks
3241
3242 for (int ih = 0; ih < nhit; ++ih) {
3243 CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(tr->GetStsHitIndex(ih));
3244 // Check clusters
3245
3246 for (Int_t side = 0; side < 2; ++side) {
3247 Int_t iclus = hit->GetDigi(side);
3248 if (mClusTr.count(iclus) == 0) {
3249 set<Int_t> aaa;
3250 mClusTr[iclus] = aaa;
3251 }
3252 mClusTr[iclus].insert(itr);
3253 }
3254 }
3255 } // for (int itr = 0; itr < ntracks;
3256
3257 Float_t* var = fVarTMVA;
3258
3259 for (int itr = 0; itr < nTracks; ++itr) {
3260 CbmStsTrack* track = (CbmStsTrack*)fVectorArray->UncheckedAt(itr);
3261 int nhit = track->GetNStsHits();
3262 // if (nhit < 4 || nhit > 5) return 1.0;
3263 if (nhit < 4 || nhit > 5) {
3264 track->SetB(1.0);
3265 continue;
3266 }
3267 const FairTrackParam* param = track->GetParamFirst();
3268 Float_t momr = 1. / param->GetQp();
3269 int iq = (momr >= 0) ? 1 : -1;
3270 TVector3 posb, pose, mom3;
3271 Float_t corrq[19], corrn[19] = {0}, corrn3 = 0, corrq3 = 0;
3272 int nmiss = 0, ista0 = -1, nsigem = 0, nseqmx = 1, nseq = 1;
3273
3274 for (int ih = 0; ih < nhit; ++ih) {
3275 CbmStsHit* hit = (CbmStsHit*)fHitArray->UncheckedAt(track->GetStsHitIndex(ih));
3276 Int_t ista = hit->GetStationNr() - 1;
3277 if (ih == 0) {
3278 hit->Position(posb);
3279 ista0 = ista;
3280 } else {
3281 if (ista - ista0 > 1) {
3282 ++nmiss;
3283 nseqmx = TMath::Max(nseqmx, nseq);
3284 nseq = 1;
3285 } else
3286 ++nseq;
3287 if (ih == nhit - 1)
3288 hit->Position(pose);
3289 }
3290 corrq[ih] = hit->GetSignalDiv();
3291 corrq3 += TMath::Abs(corrq[ih]);
3292 ista0 = ista;
3293
3294 for (Int_t side = 0; side < 2; ++side) {
3295 Int_t iclus = hit->GetDigi(side);
3296 int over = mClusTr[iclus].size();
3297 corrn[ih] += (over - 1);
3298 }
3299 corrn3 += corrn[ih];
3300 if (ista < 4)
3301 nsigem += 10;
3302 else
3303 nsigem += 1;
3304 } // for (int ih = 0; ih < nhit;
3305
3306 nseqmx = TMath::Max(nseqmx, nseq);
3307 pose -= posb;
3308 // Float_t leng = pose.Mag() / (nhit-1);
3309
3310 // Extrapolate track to the target
3311 CbmKFTrack trKF = CbmKFTrack(*track);
3312 FairTrackParam par;
3313 // track.Extrapolate(xyzv[2]);
3314 trKF.Extrapolate(0.0);
3315 trKF.GetTrackParam(par);
3316 par.Momentum(mom3);
3317 momr = TMath::Sqrt(mom3.X() * mom3.X() + mom3.Z() * mom3.Z());
3318 par.Position(mom3);
3319
3320 int iii = 0;
3321 // var[iii++] = TMath::Sqrt (xb*xb+yb*yb);
3322 // var[iii++] = posb.Pt();
3323 // var[iii++] = TMath::ATan2 (yb, xb);
3324 var[iii++] = posb.Phi();
3325 // var[iii++] = zb;
3326 // var[iii++] = posb.Z();
3327 var[iii++] = param->GetTx();
3328 var[iii++] = TMath::Log(TMath::Abs(momr));
3329 var[iii++] = track->GetChi2() / track->GetNDF();
3330 // var[iii++] = TMath::Log (TMath::Sqrt (xextr*xextr+yextr*yextr));
3331 // var[iii++] = TMath::ATan2 (yextr, xextr);
3332 // var[iii++] = TMath::Log (mom3.Pt());
3333 // var[iii++] = mom3.Phi();
3334 var[iii++] = corrq[0];
3335 var[iii++] = corrq[1];
3336 var[iii++] = corrq[2];
3337 var[iii++] = corrq[3];
3338 if (nhit == 5)
3339 var[iii++] = corrq[4];
3340 var[iii++] = corrq3;
3341 // var[iii++] = nmiss + gRandom->Rndm();
3342 var[iii++] = nseqmx;
3343 // var[iii++] = leng;
3344 // var[iii++] = corrn3;
3345 var[iii++] = corrn[0];
3346 var[iii++] = corrn[1];
3347 var[iii++] = corrn[2];
3348 var[iii++] = corrn[3];
3349 if (nhit == 5)
3350 var[iii++] = corrn[4];
3351 var[iii++] = nsigem;
3352
3353 // return fTmvaReaders[nhit*iq]->EvaluateMVA("BDTD method");
3354 track->SetB(fTmvaReaders[nhit * iq]->EvaluateMVA("BDTD method"));
3355 } // for (int itr = 0; itr < nTracks;
3356}
3357
3358// -------------------------------------------------------------------------
3359
3360void BmnStsVectorFinderV9::GoToTarget()
3361{
3362 // Evaluate vertex coordinates
3363
3364 int ntracks = fVectorArray->GetEntriesFast();
3365 Double_t xv = 0, yv = 0, wwwxtot = 0, wwwytot = 0;
3366
3367 for (int itr = 0; itr < ntracks; ++itr) {
3368 CbmStsTrack* tr = (CbmStsTrack*)fVectorArray->UncheckedAt(itr);
3369 FairTrackParam param = *(tr->GetParamFirst());
3370 fitter.Extrapolate(&param, fXyzv[2], &param);
3371 // Double_t wwwx = 1 / TMath::Sqrt(param.GetCovariance(0,0));
3372 // Double_t wwwy = 1 / TMath::Sqrt(param.GetCovariance(1,1));
3373 Double_t wwwx = 1 / param.GetCovariance(0, 0);
3374 Double_t wwwy = 1 / param.GetCovariance(1, 1);
3375 xv += (param.GetX() * wwwx);
3376 yv += (param.GetY() * wwwy);
3377 wwwxtot += wwwx;
3378 wwwytot += wwwy;
3379 // cout << itr << " " << tr->GetNStsHits() << " " << param.GetX() << " " <<
3380 // TMath::Sqrt(param.GetCovariance(0,0)) << " "
3381 // << param.GetY() << " " << TMath::Sqrt(param.GetCovariance(1,1)) << " " << tr->GetParamFirst()->GetZ() <<
3382 // endl;
3383 }
3384 cout << " Vertex coords: " << fEvent << " " << fPass << " " << xv / wwwxtot << " " << yv / wwwytot << endl;
3385}
3386
3387// -------------------------------------------------------------------------
3388
3389ClassImp(BmnStsVectorFinderV9)
ClassImp
ClassImp(BmnH3LTripleFinder)
dist
vector< Double_t > dist(vector< Double_t > qp, Double_t mu)
Definition BmnMath.cxx:869
d
const Float_t d
Z-ccordinate of the first GEM-station.
Definition BmnMwpcHit.cxx:7
i
int i
Definition P4_F32vec4.h:22
kGEM
@ kGEM
Definition BmnDetectorList.h:15
kVSP
@ kVSP
Definition BmnDetectorList.h:21
BmnGemStripModule::IsPointInsideModule
Bool_t IsPointInsideModule(Double_t x, Double_t y)
Definition BmnGemStripModule.cxx:160
BmnGemStripStationSet::GetStation
BmnGemStripStation * GetStation(Int_t station_num)
Definition BmnGemStripStationSet.h:45
BmnGemStripStation::GetModule
BmnGemStripModule * GetModule(Int_t module_num)
Definition BmnGemStripStation.cxx:183
BmnSiliconModule::IsPointInsideModule
Bool_t IsPointInsideModule(Double_t x, Double_t y, Bool_t isLocal)
Definition BmnSiliconModule.cxx:251
BmnSiliconStationSet::GetStation
BmnSiliconStation * GetStation(Int_t station_num)
Definition BmnSiliconStationSet.h:41
BmnSiliconStation::GetModule
BmnSiliconModule * GetModule(Int_t module_num)
Definition BmnSiliconStation.cxx:95
BmnStsKFTrackFitter::Extrapolate
void Extrapolate(CbmStsTrack *track, Double_t z, FairTrackParam *e_track)
Definition BmnStsKFTrackFitter.cxx:480
BmnStsKFTrackFitter::ReadMatBudget
void ReadMatBudget(TString &matBudgetFileName)
Definition BmnStsKFTrackFitter.cxx:159
BmnStsKFTrackFitter::SetKFHits
void SetKFHits(CbmKFTrack &T, CbmStsTrack *track)
Definition BmnStsKFTrackFitter.cxx:83
BmnStsKFTrackFitter::DoFit
Int_t DoFit(CbmStsTrack *track, Int_t pidHypo=211)
Definition BmnStsKFTrackFitter.cxx:114
BmnStsKFTrackFitter::Fit
Int_t Fit(CbmStsTrack *track, Int_t pidHypo=211)
Definition BmnStsKFTrackFitter.cxx:194
BmnStsVectorFinderV9::Init
virtual InitStatus Init()
Definition BmnStsVectorFinderV9.cxx:86
BmnStsVectorFinderV9::Exec
virtual void Exec(Option_t *opt)
Definition BmnStsVectorFinderV9.cxx:392
CbmKFHit::Filter
virtual Int_t Filter(CbmKFTrackInterface &track, Bool_t downstream, Double_t &QP0)=0
CbmKFTrackInterface::Extrapolate
Int_t Extrapolate(Double_t z, Double_t *QP0=0, Bool_t line=false)
Access to i-th hit.
Definition CbmKFTrackInterface.cxx:37
CbmKFTrack::SetStsTrack
void SetStsTrack(CbmStsTrack &track, bool first=1)
Definition CbmKFTrack.cxx:33
CbmKFTrack::GetHit
CbmKFHit * GetHit(Int_t i)
Number of hits.
Definition CbmKFTrack.h:58
CbmKFTrack::GetTrackParam
void GetTrackParam(FairTrackParam &track)
Definition CbmKFTrack.cxx:41
CbmKFTrack::GetRefChi2
Double_t & GetRefChi2()
array[15] of covariance matrix
Definition CbmKFTrack.h:53
CbmKF::Instance
static CbmKF * Instance()
Definition CbmKF.h:35
CbmKF::GetNStsStations
int GetNStsStations() const
Definition CbmKF.h:88
CbmStsCluster::GetNDigis
Int_t GetNDigis() const
Definition CbmStsCluster.h:69
CbmStsCluster::GetDigi
Int_t GetDigi(Int_t inum)
Definition CbmStsCluster.h:70
CbmStsDigiMatch::GetRefIndex
Int_t GetRefIndex(Int_t i=0) const
Definition CbmStsDigiMatch.cxx:65
CbmStsDigiScheme::Instance
static CbmStsDigiScheme * Instance(int version=1)
Definition CbmStsDigiScheme.cxx:65
CbmStsDigiScheme::GetStationByNr
CbmStsStation * GetStationByNr(Int_t stationNr)
Definition CbmStsDigiScheme.cxx:193
CbmStsHit::GetStationNr
virtual Int_t GetStationNr() const
Definition CbmStsHit.h:66
CbmStsHit::GetSignalDiv
Double_t GetSignalDiv() const
Definition CbmStsHit.h:74
CbmStsHit::GetSystemId
Int_t GetSystemId() const
Definition CbmStsHit.h:64
CbmStsHit::GetDigi
Int_t GetDigi(Int_t iSide) const
Definition CbmStsHit.cxx:116
CbmStsSector::GetDx
Double_t GetDx() const
Definition CbmStsSector.h:95
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
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
CbmStsTrack::SetB
void SetB(Double_t b)
Definition CbmStsTrack.h:88
file_inspector.code
code
Definition file_inspector.py:30
nlohmann::to_string
NLOHMANN_BASIC_JSON_TPL_DECLARATION std::string to_string(const NLOHMANN_BASIC_JSON_TPL &j)
user-defined to_string function for JSON values
Definition json.hpp:21838