CbmRoot
CbmTofTrackFinderNN.cxx
Go to the documentation of this file.
1 // ROOT Classes and includes
2 #include "TClonesArray.h"
3 #include "TDirectory.h"
4 #include "TGeoManager.h"
5 #include "TROOT.h"
6 #include <TCanvas.h>
7 #include <TF2.h>
8 #include <TGraph2D.h>
9 #include <TGraph2DErrors.h>
10 #include <TH1.h>
11 #include <TMath.h>
12 #include <TRandom2.h>
13 #include <TStyle.h>
14 
15 // FAIR classes and includes
16 #include "FairEventManager.h" // for FairEventManager
17 #include "FairLogger.h"
18 #include "FairRootManager.h"
19 #include "FairRunAna.h"
20 #include "FairRuntimeDb.h"
21 #include "TEveManager.h" // for TEveManager, gEve
22 
23 // CBMroot classes and includes
24 #include "CbmMatch.h"
25 #include "CbmTofAddress.h" // in cbmdata/tof
26 #include "CbmTofCell.h" // in tof/TofData
27 #include "CbmTofDetectorId_v12b.h" // in cbmdata/tof
28 #include "CbmTofDetectorId_v14a.h" // in cbmdata/tof
29 #include "CbmTofDigiBdfPar.h" // in tof/TofParam
30 #include "CbmTofDigiPar.h" // in tof/TofParam
31 #include "CbmTofFindTracks.h"
32 #include "CbmTofGeoHandler.h" // in tof/TofTools
33 #include "CbmTofHit.h" // in cbmdata/tof
34 #include "CbmTofTrackFinderNN.h"
35 #include "CbmTofTrackFitter.h"
36 #include "CbmTofTracklet.h"
37 #include "CbmTofTrackletParam.h"
38 #include "LKFMinuit.h"
39 #include <CbmEvDisTracks.h> // in eventdisplay/tof
40 
41 // C++ includes
42 #include <map>
43 #include <vector>
44 using std::cout;
45 using std::endl;
46 using std::map;
47 
48 const Int_t DetMask = 0x3FFFFF; // check for consistency with geometry
49 LKFMinuit CbmTofTrackFinderNN::fMinuit;
50 
51 CbmTofTrackFinderNN::CbmTofTrackFinderNN()
52  : fHits(NULL)
53  , fOutTracks(NULL)
54  , fiNtrks(0)
55  , fFitter(NULL)
56  , fFindTracks(NULL)
57  , fDigiPar(NULL)
58  , fMaxTofTimeDifference(0.)
59  , fTxLIM(0.)
60  , fTyLIM(0.)
61  , fTxMean(0.)
62  , fTyMean(0.)
63  , fSIGLIM(4.)
64  , fChiMaxAccept(3.)
65  , fPosYMaxScal(0.55)
66  , fTracks()
67  , fvTrkVec() {}
68 
69 CbmTofTrackFinderNN::~CbmTofTrackFinderNN() {}
70 
71 //Copy constructor
72 CbmTofTrackFinderNN::CbmTofTrackFinderNN(const CbmTofTrackFinderNN& finder)
73  : fHits(NULL)
74  , fOutTracks(NULL)
75  , fiNtrks(0)
76  , fFitter(NULL)
77  , fFindTracks(NULL)
78  , fDigiPar(NULL)
79  , fMaxTofTimeDifference(0.)
80  , fTxLIM(0.)
81  , fTyLIM(0.)
82  , fTxMean(0.)
83  , fTyMean(0.)
84  , fSIGLIM(4.)
85  , fChiMaxAccept(3.)
86  , fPosYMaxScal(0.55)
87  , fTracks()
88  , fvTrkVec() {
89  // action
90  fHits = finder.fHits;
91  fTracks = finder.fTracks;
92  fiNtrks = finder.fiNtrks;
93 }
94 
95 // assignement operator
96 CbmTofTrackFinderNN&
97 CbmTofTrackFinderNN::operator=(const CbmTofTrackFinderNN& /*fSource*/) {
98  // do copy
99  // ... (too lazy) ...
100  // return the existing object
101  return *this;
102 }
103 
104 void CbmTofTrackFinderNN::Init() {
105  /*
106  CbmTofFindTracks* fFindTracks = CbmTofFindTracks::Instance();
107  fFitter = fFindTracks->GetFitter();
108  */
109  FairRunAna* ana = FairRunAna::Instance();
110  FairRuntimeDb* rtdb = ana->GetRuntimeDb();
111  fDigiPar = (CbmTofDigiPar*) (rtdb->getContainer("CbmTofDigiPar"));
112  if (NULL == fDigiPar) {
113  LOG(error)
114  << "CbmTofTrackFinderNN::Init => Could not obtain the CbmTofDigiPar ";
115  }
116 
117  LOG(info) << Form(" CbmTofTrackFinderNN::Init : Fitter at 0x%p", fFitter);
118 
119  fFindTracks = CbmTofFindTracks::Instance();
120  if (NULL == fFindTracks)
121  LOG(fatal) << Form(
122  " CbmTofTrackFinderNN::Init : no FindTracks instance found");
123 
124  fMinuit.Initialize();
125 
126  LOG(info) << "MaxTofTimeDifference = " << fMaxTofTimeDifference;
127 }
128 
129 Int_t CbmTofTrackFinderNN::DoFind(TClonesArray* fTofHits,
130  TClonesArray* fTofTracks) {
131  fiNtrks = 0; // initialize
132  fHits = fTofHits;
133  fOutTracks = fTofTracks; //new TClonesArray("CbmTofTracklet");
134  //fTracks = new TClonesArray("CbmTofTracklet");
135  //if (0 == fFindTracks->GetStationType(0)){ // Generate Pseudo TofHit at origin
136  if (
137  0) { // == fFindTracks->GetAddrOfStation(0)) { // generate new track seed, disabled
138  fFindTracks->SetStation(0, 0, 0, 0);
139  const Int_t iDetId = CbmTofAddress::GetUniqueAddress(0, 0, 0, 0, 0);
140  const TVector3 hitPos(0., 0., 0.);
141  const TVector3 hitPosErr(0.5, 0.5, 0.5); // initialize fake hit error
142  const Double_t dTime0 = 0.; // FIXME
143 
144  Int_t iNbHits = fHits->GetEntries();
145  /*CbmTofHit *pHit0 =*/new ((*fHits)[iNbHits]) CbmTofHit(
146  iDetId,
147  hitPos,
148  hitPosErr, // local detector coordinates
149  iNbHits, // this number is used as reference!!
150  dTime0, // Time of hit
151  0, //vPtsRef.size(), // flag = number of TofPoints generating the cluster
152  0);
153  LOG(debug1)
154  << "CbmTofTrackFinderNN::DoFind: Fake Hit at origin added at position "
155  << iNbHits << Form(", DetId 0x%08x", iDetId);
156  }
157 
158  // fvTrkMap.resize(fHits->GetEntries());
159  fvTrkVec.resize(fHits->GetEntries());
160  LOG(debug2) << "<I> TrkMap/Vec resized for " << fHits->GetEntries()
161  << " entries ";
162  // for (Int_t iHit=0; iHit<fHits->GetEntries(); iHit++) { fvTrkMap[iHit].clear();}
163  for (Int_t iHit = 0; iHit < fHits->GetEntries(); iHit++) {
164  fvTrkVec[iHit].clear();
165  }
166 
167  LOG(debug) << "MaxTofTimeDifference = " << fMaxTofTimeDifference;
168  Int_t iNTrks = 0;
169  Int_t iSt0 = -1;
170  Int_t iSt1 = 0;
171  while (iSt0
172  < fFindTracks->GetNofStations()
173  - fFindTracks
174  ->GetMinNofHits()) { // seed loop, all combinations as seeds
175  iSt0++;
176  iSt1 = iSt0;
177  while (iSt1
178  < fFindTracks->GetNofStations() - fFindTracks->GetMinNofHits() + 1) {
179  iSt1++;
180  for (Int_t iHit = 0; iHit < fHits->GetEntries();
181  iHit++) { // loop over Hits
182  CbmTofHit* pHit = (CbmTofHit*) fHits->At(iHit);
183  Int_t iAddr = (pHit->GetAddress() & DetMask);
184  // Int_t iSmType = CbmTofAddress::GetSmType( iAddr ); (VF) not used
185  if (HitUsed(iHit) == 1 && iAddr != fFindTracks->GetBeamCounter())
186  continue; // skip used Hits except for BeamCounter
187  LOG(debug1) << Form(
188  "<I> TofTracklet Chkseed St0 %2d, St1 %2d, Mul %2d, Hit %2d, addr = "
189  "0x%08x - X %6.2f, Y %6.2f Z %6.2f R %6.2f T %6.2f TM %lu",
190  iSt0,
191  iSt1,
192  fiNtrks,
193  iHit,
194  pHit->GetAddress(),
195  pHit->GetX(),
196  pHit->GetY(),
197  pHit->GetZ(),
198  pHit->GetR(),
199  pHit->GetTime(),
200  fvTrkVec[iHit].size());
201  if (iAddr
202  == fFindTracks->GetAddrOfStation(
203  iSt0)) { // generate new track seed
204  LOG(debug) << Form(
205  "<I> TofTracklet seed St0 %2d, St1 %2d, Mul %2d, Hit %2d, addr = "
206  "0x%08x - X %6.2f, Y %6.2f Z %6.2f T %6.2f TM %lu",
207  iSt0,
208  iSt1,
209  fiNtrks,
210  iHit,
211  pHit->GetAddress(),
212  pHit->GetX(),
213  pHit->GetY(),
214  pHit->GetZ(),
215  pHit->GetTime(),
216  fvTrkVec[iHit].size());
217 
218  Int_t iChId = pHit->GetAddress();
219  CbmTofCell* fChannelInfo = fDigiPar->GetCell(iChId);
220  Int_t iCh = CbmTofAddress::GetChannelId(iChId);
221  Double_t hitpos[3] = {3 * 0.};
222  Double_t hitpos_local[3] = {3 * 0.};
223  Double_t dSizey = 1.;
224 
225  if (1) { // iSmType>0) { // prevent geometry inspection for FAKE hits
226  if (NULL == fChannelInfo) {
227  LOG(fatal) << "<D> CbmTofTrackFinderNN::DoFind0: Invalid Channel "
228  "Pointer for ChId "
229  << Form(" 0x%08x ", iChId) << ", Ch " << iCh;
230  // continue;
231  } else {
232  /*TGeoNode *fNode= */ // prepare global->local trafo
233  gGeoManager->FindNode(fChannelInfo->GetX(),
234  fChannelInfo->GetY(),
235  fChannelInfo->GetZ());
236  if (fFindTracks->GetAddrOfStation(iSt0)
237  != fFindTracks->GetBeamCounter()) {
238  hitpos[0] = pHit->GetX();
239  hitpos[1] = pHit->GetY();
240  hitpos[2] = pHit->GetZ();
241  /*TGeoNode* cNode=*/gGeoManager->GetCurrentNode();
242  gGeoManager->MasterToLocal(hitpos, hitpos_local);
243  }
244  dSizey = fChannelInfo->GetSizey();
245  LOG(debug) << Form(
246  "<D> TofTracklet start %d, Hit %d - yloc %6.2f, dy %6.2f, Scal "
247  "%6.2f -> stations 0x%08x, 0x%08x,",
248  fiNtrks,
249  iHit,
250  hitpos_local[1],
251  dSizey,
252  fPosYMaxScal,
253  fFindTracks->GetAddrOfStation(iSt0),
254  fFindTracks->GetAddrOfStation(iSt1));
255  }
256  }
257 
258  if (TMath::Abs(hitpos_local[1]) < dSizey * fPosYMaxScal)
259  for (Int_t iHit1 = 0; iHit1 < fHits->GetEntries();
260  iHit1++) // loop over all Hits (order unknown)
261  {
262  if (HitUsed(iHit1) == 1) continue; // skip used Hits
263  CbmTofHit* pHit1 = (CbmTofHit*) fHits->At(iHit1);
264  // Int_t iSmType1 = CbmTofAddress::GetSmType( pHit1->GetAddress() & DetMask ); (VF) not used
265  Int_t iAddr1 = (pHit1->GetAddress() & DetMask);
266  //if (iSmType1 == fFindTracks->GetStationType(1)) { // generate new track seed
267  if (iAddr1
268  == fFindTracks->GetAddrOfStation(
269  iSt1)) { // generate new track seed
270  Int_t iChId1 = pHit1->GetAddress();
271  CbmTofCell* fChannelInfo1 = fDigiPar->GetCell(iChId1);
272  Int_t iCh1 = CbmTofAddress::GetChannelId(iChId1);
273  Double_t hitpos1[3] = {3 * 0.};
274  Double_t hitpos1_local[3] = {3 * 0.};
275  Double_t dSizey1 = 1.;
276  if (NULL == fChannelInfo1) {
277  LOG(debug) << "CbmTofTrackFinderNN::DoFindi: Invalid Channel "
278  "Pointer for ChId "
279  << Form(" 0x%08x ", iChId1) << ", Ch " << iCh1;
280  // continue;
281  } else {
282  /*TGeoNode *fNode1=*/ // prepare global->local trafo
283  gGeoManager->FindNode(fChannelInfo1->GetX(),
284  fChannelInfo1->GetY(),
285  fChannelInfo1->GetZ());
286  hitpos1[0] = pHit1->GetX();
287  hitpos1[1] = pHit1->GetY();
288  hitpos1[2] = pHit1->GetZ();
289  /*TGeoNode* cNode1=*/gGeoManager->GetCurrentNode();
290  gGeoManager->MasterToLocal(hitpos1, hitpos1_local);
291  dSizey1 = fChannelInfo1->GetSizey();
292  }
293  Double_t dDT = pHit1->GetTime() - pHit->GetTime();
294  //if(dDT<0.) continue; // request forward propagation in time
295 
296  Double_t dLz = pHit1->GetZ() - pHit->GetZ();
297  Double_t dTx = (pHit1->GetX() - pHit->GetX()) / dLz;
298  Double_t dTy = (pHit1->GetY() - pHit->GetY()) / dLz;
299  Double_t dDist = TMath::Sqrt(
300  TMath::Power((pHit->GetX() - pHit1->GetX()), 2)
301  + TMath::Power((pHit->GetY() - pHit1->GetY()), 2)
302  + TMath::Power((pHit->GetZ() - pHit1->GetZ()), 2));
303  LOG(debug1) << Form(
304  "<I> TofTracklet %d, Hits %d, %d, add = 0x%08x,0x%08x - DT "
305  "%6.2f, Tx %6.2f Ty %6.2f Tt %6.2f pos %6.2f %6.2f %6.2f ",
306  fiNtrks,
307  iHit,
308  iHit1,
309  pHit->GetAddress(),
310  pHit1->GetAddress(),
311  dDT,
312  dTx,
313  dTy,
314  dDT / dLz,
315  hitpos1_local[0],
316  hitpos1_local[1],
317  hitpos1_local[2]);
318  LOG(debug3)
319  << Form(" selection: y %6.2f < %6.2f, T %6.2f < %6.2f, "
320  "dTpos %6.2f < %6.2f, Abs(%6.2f - %6.2f) < %6.2f",
321  hitpos1_local[1],
322  dSizey1 * fPosYMaxScal,
323  dDT / dLz,
324  fMaxTofTimeDifference,
325  dTx,
326  fTxLIM,
327  dTy,
328  fTyMean,
329  fTyLIM);
330 
331  if (TMath::Abs(hitpos1_local[1]) < dSizey1 * fPosYMaxScal)
332  if (TMath::Abs(dDT / dLz) < fMaxTofTimeDifference
333  && TMath::Abs(dDT / dDist)
334  > 0.005 // FIXME: numeric constant in code
335  && TMath::Abs(dTx - fTxMean) < fTxLIM
336  && TMath::Abs(dTy - fTyMean) < fTyLIM) {
337  CbmTofTracklet* pTrk = new CbmTofTracklet();
338  fTracks.push_back(pTrk);
339  pTrk->SetTofHitIndex(
340  iHit, iAddr, pHit); // store 1. Hit index
341  pTrk->AddTofHitIndex(
342  iHit1, iAddr1, pHit1); // store 2. Hit index
343  fiNtrks = fTracks.size();
344 
345  fvTrkVec[iHit].push_back(pTrk);
346  fvTrkVec[iHit1].push_back(pTrk);
347 
348  pTrk->SetTime(
349  pHit->GetTime()); // define reference time from 1. plane
350  //Double_t dR = pHit1->GetR() - pHit->GetR();
351  Double_t dR = pTrk->Dist3D(pHit1, pHit);
352  Double_t dTt =
353  fFindTracks
354  ->GetTtTarg(); // assume calibration target value
355  if (0) { //== iSmType) { // disabled
356  Double_t T0Fake = pHit->GetTime();
357  Double_t w = fvTrkVec[iHit].size();
358  T0Fake =
359  (T0Fake * (w - 1.) + (pHit1->GetTime() - dTt * dR)) / w;
360  LOG(debug1) << Form(
361  "<I> TofTracklet %d, Fake T0, old %8.0f -> new %8.0f",
362  fiNtrks,
363  pHit->GetTime(),
364  T0Fake);
365  pHit->SetTime(T0Fake);
366  }
367  Double_t dSign = 1.;
368  if (pHit1->GetZ() < pHit->GetZ()) dSign = -1.;
369  dTt = dSign * dDT / dR;
370  pTrk->SetTt(dTt); // store inverse velocity
371  pTrk->UpdateT0();
372  CbmTofTrackletParam* tPar = pTrk->GetTrackParameter();
373  tPar->SetX(pHit->GetX()); // fill TrackParam
374  tPar->SetY(pHit->GetY()); // fill TrackParam
375  tPar->SetZ(pHit->GetZ()); // fill TrackParam
376  tPar->SetLz(dLz);
377  tPar->SetTx(dTx);
378  tPar->SetTy(dTy);
379 
380  LOG(debug)
381  << Form("<I> TofTracklet %d, Hits %d, %d initialized, "
382  "add 0x%08x,0x%08x, DT %6.3f, Sgn %2.0f, DR "
383  "%6.3f, T0 %6.2f, Tt %6.4f ",
384  fiNtrks,
385  iHit,
386  iHit1,
387  pHit->GetAddress(),
388  pHit1->GetAddress(),
389  dDT,
390  dSign,
391  dR,
392  pTrk->GetT0(),
393  pTrk->GetTt())
394  // << tPar->ToString()
395  ;
396  }
397  }
398  }
399  } // iSt0 condition end
400  } // Loop on Hits end
401 
402  if (iNTrks >= 0 && static_cast<size_t>(iNTrks) == fTracks.size())
403  continue; // nothing new
404  iNTrks = fTracks.size();
405  PrintStatus((char*) Form(
406  "after seeds of St0 %d, St1 %d, Mul %d", iSt0, iSt1, iNTrks));
407 
408  const Int_t MAXNCAND =
409  1000; // Max number of hits matchable to current tracklets
410  // Propagate track seeds to remaining detectors
411  for (Int_t iDet = iSt1 + 1; iDet < fFindTracks->GetNStations(); iDet++) {
412  Int_t iNCand = 0;
413  Int_t iHitInd[MAXNCAND];
414  CbmTofTracklet* pTrkInd[MAXNCAND];
415  Double_t dChi2[MAXNCAND];
416  for (size_t iTrk = 0; iTrk < fTracks.size();
417  iTrk++) { // loop over Trackseeds
418  CbmTofTracklet* pTrk = (CbmTofTracklet*) fTracks[iTrk];
419  LOG(debug3) << " Propagate Loop " << iTrk << " pTrk " << pTrk
420  << Form(" to station %d, addr: 0x%08x ",
421  iDet,
422  fFindTracks->GetAddrOfStation(iDet));
423  if (NULL == pTrk) continue;
424 
425  for (Int_t iHit = 0; iHit < fHits->GetEntries();
426  iHit++) { // loop over Hits
427  if (HitUsed(iHit) == 1) continue; // skip used Hits
428  CbmTofHit* pHit = (CbmTofHit*) fHits->At(iHit);
429  Int_t iAddr = (pHit->GetAddress() & DetMask);
430  if (iAddr != fFindTracks->GetAddrOfStation(iDet)) continue;
431 
432  // Int_t iSmType = CbmTofAddress::GetSmType( pHit->GetAddress() & DetMask ); (VF) not used
433  Int_t iChId = pHit->GetAddress();
434  CbmTofCell* fChannelInfo = fDigiPar->GetCell(iChId);
435  Int_t iCh = CbmTofAddress::GetChannelId(iChId);
436  Double_t hitpos[3] = {3 * 0.};
437  Double_t hitpos_local[3] = {3 * 0.};
438  Double_t dSizey = 1.;
439 
440  if (NULL == fChannelInfo) {
441  LOG(debug) << "CbmTofTrackFinderNN::DoFind: Invalid Channel "
442  "Pointer from Hit "
443  << iHit << " for ChId " << Form(" 0x%08x ", iChId)
444  << ", Ch " << iCh;
445  // continue;
446  } else {
447  /*TGeoNode *fNode=*/ // prepare global->local trafo
448  gGeoManager->FindNode(fChannelInfo->GetX(),
449  fChannelInfo->GetY(),
450  fChannelInfo->GetZ());
451  hitpos[0] = pHit->GetX();
452  hitpos[1] = pHit->GetY();
453  hitpos[2] = pHit->GetZ();
454  /*TGeoNode* cNode=*/gGeoManager->GetCurrentNode();
455  gGeoManager->MasterToLocal(hitpos, hitpos_local);
456  dSizey = fChannelInfo->GetSizey();
457  }
458  if (TMath::Abs(hitpos_local[1])
459  < dSizey
460  * fPosYMaxScal) { // extrapolate Tracklet to this station
461  if (pTrk->GetStationHitIndex(iAddr) > -1)
462  continue; // Station already part of this tracklet
463  CbmTofTrackletParam* tPar = pTrk->GetTrackParameter();
464  Int_t iHit0 = pTrk->GetTofHitIndex(0);
465  Int_t iHit1 = pTrk->GetTofHitIndex(1);
466  // CbmTofHit* pHit0 = (CbmTofHit*) fHits->At( iHit0 ); (VF) not used
467  // CbmTofHit* pHit1 = (CbmTofHit*) fHits->At( iHit1 ); (VF) not used
468  if (iHit0 < 0 || iHit0 >= fHits->GetEntries())
469  LOG(fatal) << "CbmTofTrackFinderNN::DoFind Invalid Hit Index "
470  << iHit0 << " for Track " << iTrk << "("
471  << fTracks.size() << ")";
472  Double_t dDz = pHit->GetZ() - tPar->GetZ();
473  Double_t dXex =
474  tPar->GetX()
475  + tPar->GetTx() * dDz; // pTrk->GetFitX(pHit->GetZ());
476  Double_t dYex =
477  tPar->GetY()
478  + tPar->GetTy() * dDz; // pTrk->GetFitY(pHit->GetZ());
479  /*
480  Double_t dDr = pHit->GetR() - pHit0->GetR();
481  Double_t dTex = pHit0->GetTime() + (pHit1->GetTime()-pHit0->GetTime())/(pHit1->GetR()-pHit0->GetR())*dDr;
482  */
483  Double_t dTex = pTrk->GetTex(pHit);
484  // pTrk->GetFitT(pHit->GetR());
485 
486  Double_t dChi =
487  TMath::Sqrt((TMath::Power(TMath::Abs(dTex - pHit->GetTime())
488  / fFindTracks->GetSigT(iAddr),
489  2)
490  + TMath::Power(TMath::Abs(dXex - pHit->GetX())
491  / fFindTracks->GetSigX(iAddr),
492  2)
493  + TMath::Power(TMath::Abs(dYex - pHit->GetY())
494  / fFindTracks->GetSigY(iAddr),
495  2))
496  / 3);
497 
498  LOG(debug1) << Form(
499  "<IP> TofTracklet %lu, HMul %d, Hits %d, %d check %d, Station "
500  "0x%08x: DT %f, DX %f, DY %f, Chi %f",
501  iTrk,
502  pTrk->GetNofHits(),
503  iHit0,
504  iHit1,
505  iHit,
506  iAddr,
507  (dTex - pHit->GetTime()) / fFindTracks->GetSigT(iAddr),
508  (dXex - pHit->GetX()) / fFindTracks->GetSigX(iAddr),
509  (dYex - pHit->GetY()) / fFindTracks->GetSigY(iAddr),
510  dChi);
511 
512  if (
513  dChi
514  < fSIGLIM) // FIXME: should scale limit with material budget between hit and track reference
515  { // extend and update tracklet
516  LOG(debug) << Form("<IP> TofTracklet %lu, HMul %d, Hits %d, %d "
517  "mark for extension by %d, add = 0x%08x, DT "
518  "%6.2f, DX %6.2f, DY=%6.2f ",
519  iTrk,
520  pTrk->GetNofHits(),
521  iHit0,
522  iHit1,
523  iHit,
524  pHit->GetAddress(),
525  dTex - pHit->GetTime(),
526  dXex - pHit->GetX(),
527  dYex - pHit->GetY())
528  << tPar->ToString();
529 
530  if (iNCand > 0) {
531  LOG(debug) << Form("CbmTofTrackFinderNN::DoFind new match %d "
532  "of Hit %d, Trk %lu, chi2 = %f",
533  iNCand,
534  iHit,
535  iTrk,
536  dChi);
537  iNCand++;
538  if (iNCand == MAXNCAND) iNCand--; // Limit index to maximum
539 
540  for (Int_t iCand = 0; iCand < iNCand; iCand++) {
541  if (dChi < dChi2[iCand]) {
542  for (Int_t iCC = iNCand; iCC > iCand; iCC--) {
543  pTrkInd[iCC] = pTrkInd[iCC - 1];
544  iHitInd[iCC] = iHitInd[iCC - 1];
545  dChi2[iCC] = dChi2[iCC - 1];
546  }
547  pTrkInd[iCand] = pTrk;
548  iHitInd[iCand] = iHit;
549  dChi2[iCand] = dChi;
550  dChi2[iNCand] = 1.E8;
551  LOG(debug1)
552  << Form(" <D> candidate inserted at pos %d", iCand);
553  break;
554  }
555  }
556  } else {
557  LOG(debug) << Form("CbmTofTrackFinderNN::DoFind first match "
558  "%d of Hit %d, Trk %p, chi2 = %f",
559  iNCand,
560  iHit,
561  pTrk,
562  dChi);
563  pTrkInd[iNCand] = pTrk;
564  iHitInd[iNCand] = iHit;
565  dChi2[iNCand] = dChi; // relative quality measure
566  iNCand++;
567  dChi2[iNCand] = 1.E8;
568  }
569  }
570 
571  } // hit y position check end
572  } // hit loop end
573  } // loop over tracklets end
574 
575  while (iNCand > 0) { // at least one matching hit - trk pair found
576  CbmTofTracklet* pTrk = pTrkInd[0];
577  if (NULL == pTrk) continue;
578  LOG(debug) << Form(
579  "%d hit match candidates in station %d to %lu TofTracklets",
580  iNCand,
581  iDet,
582  fTracks.size());
583  for (Int_t iM = 0; iM < iNCand; iM++) {
584  pTrk = (CbmTofTracklet*) pTrkInd[iM];
585  if (NULL == pTrk) break;
586  std::vector<CbmTofTracklet*>::iterator it =
587  std::find(fTracks.begin(), fTracks.end(), pTrk);
588  if (it == fTracks.end()) break; // track candidate not existing
589 
590  LOG(debug1) << "\t"
591  << Form("Hit %d, Trk %p with chi2 %f (%f)",
592  iHitInd[iM],
593  pTrkInd[iM],
594  dChi2[iM],
595  pTrk->GetMatChi2(
596  fFindTracks->GetAddrOfStation(iDet)));
597  }
598  PrintStatus((char*) "starting NCand");
599 
600  // check if best pTrk still active
601  pTrk = (CbmTofTracklet*) pTrkInd[0];
602  size_t iTr = 0;
603  for (; iTr < fTracks.size(); iTr++) {
604  if (fTracks[iTr] == pTrk) {
605  LOG(debug) << "Track " << pTrk << " active at pos " << iTr;
606  break;
607  }
608  }
609  if (iTr == fTracks.size()) {
610  iNCand--;
611  for (Int_t iCand = 0; iCand < iNCand; iCand++) {
612  pTrkInd[iCand] = pTrkInd[iCand + 1];
613  iHitInd[iCand] = iHitInd[iCand + 1];
614  dChi2[iCand] = dChi2[iCand + 1];
615  }
616  continue;
617  }
618 
619  // CbmTofTrackletParam *tPar = pTrk->GetTrackParameter(); (VF) not used
620  Int_t iHit0 = pTrk->GetTofHitIndex(0);
621  Int_t iHit1 = pTrk->GetTofHitIndex(1);
622  if (pTrk->GetNofHits() > fFindTracks->GetNStations()
623  || pTrk->GetNofHits() <= 0)
624  LOG(fatal) << " No or more Tracklet hits than stations ! Stop ";
625  // check if tracklet already contains a hit of this layer
626  Int_t iHit = iHitInd[0];
627  CbmTofHit* pHit = (CbmTofHit*) fHits->At(iHit);
628  Int_t iAddr = (pHit->GetAddress() & DetMask);
629  if (Double_t dLastChi2 =
630  pTrk->GetMatChi2(fFindTracks->GetAddrOfStation(iDet)) == -1.) {
631  LOG(debug1) << Form(
632  " -D- Add hit %d at %p, Addr 0x%08x, Chi2 %6.2f to size %u",
633  iHit,
634  pHit,
635  iAddr,
636  dChi2[0],
637  pTrk->GetNofHits());
638  pTrk->AddTofHitIndex(
639  iHit,
640  iAddr,
641  pHit,
642  dChi2[0]); // store next Hit index with matching chi2
643  // pTrk->PrintInfo();
644  fvTrkVec[iHit].push_back(pTrk);
645  Line3Dfit(pTrk); // full MINUIT fit overwrites ParamLast!
646  Bool_t bkeep = kFALSE;
647  if (pTrk->GetChiSq() > fChiMaxAccept) {
648  LOG(debug) << Form("Add hit %d invalidates tracklet with Chi "
649  "%6.2f > %6.2f -> undo ",
650  iHit,
651  pTrk->GetChiSq(),
652  fChiMaxAccept);
653  fvTrkVec[iHit].pop_back();
654  pTrk->RemoveTofHitIndex(iHit, iAddr, pHit, dChi2[0]);
655  Line3Dfit(pTrk); //restore old status
656  bkeep = kTRUE;
657  }
658 
659  if (bkeep) {
660  iNCand--;
661  for (Int_t iCand = 0; iCand < iNCand; iCand++) {
662  pTrkInd[iCand] = pTrkInd[iCand + 1];
663  iHitInd[iCand] = iHitInd[iCand + 1];
664  dChi2[iCand] = dChi2[iCand + 1];
665  }
666  } else { // update chi2 array
667  PrintStatus((char*) "before UpdateTrackList");
668  UpdateTrackList(pTrk);
669  Int_t iNCandNew = 0;
670  for (Int_t iCand = 0; iCand < iNCand; iCand++) {
671  if (pTrk != pTrkInd[iCand] && iHit != iHitInd[iCand]) {
672  pTrkInd[iNCandNew] = pTrkInd[iCand];
673  iHitInd[iNCandNew] = iHitInd[iCand];
674  dChi2[iNCandNew] = dChi2[iCand];
675  iNCandNew++;
676  }
677  }
678  iNCand = iNCandNew;
679  }
680  } else {
681  if (dChi2[0] < dLastChi2) { // replace hit index
682  LOG(fatal) << Form(
683  "-D- Replace %d, Addr 0x%08x, at %p, Chi2 %6.2f",
684  iHit,
685  iAddr,
686  pHit,
687  dChi2[0]);
688  //cout << " -D- Replace " << endl;
689  pTrk->ReplaceTofHitIndex(iHit, iAddr, pHit, dChi2[0]);
690  // TODO remove tracklet assigment of old hit! FIXME
691  } else {
692  LOG(debug) << Form(
693  " -D- Ignore %d, Det %d, Addr 0x%08x, at 0x%p, Chi2 %6.2f",
694  iHit,
695  iDet,
696  iAddr,
697  pHit,
698  dChi2[0]);
699  // Form new seeds
700  //if (iDet<(fFindTracks->GetNStations()-1)) TrklSeed(fHits,fTracks,iHit);
701  break;
702  }
703  }
704 
705  // pTrk->SetParamLast(tPar); // Initialize FairTrackParam for KF
706  //fFitter->DoFit(pTrk); //whatever that means ... KF - Fitting
707  //pTrk->GetFairTrackParamLast(); // transfer fit result to CbmTofTracklet
708  //pTrk->SetTime(pHit->GetTime()); // update reference time
709 
710  //Line3Dfit(pTrk); // full MINUIT fit for debugging overwrites ParamLast!
711 
712  pTrk->SetTime(
713  pTrk->UpdateT0()); // update reference time (and fake hit time)
714 
715  //FairTrackParam paramExtr;
716  //fFitter->Extrapolate(pTrk->GetParamLast(),0.,&paramExtr);
717  //pTrk->GetParamFirst()->Print();
718  //pTrk->GetParamLast()->Print();
719  //paramExtr.Print();
720 
721  // check with ROOT fitting method
722  //TLinearFitter *lf=new TLinearFitter(3);
723  //lf->SetFormula("hyp3");
724 
725  // update inverse velocity
726  Double_t dTt = pTrk->GetTt();
727 
728  LOG(debug) << Form("<Res> TofTracklet %p, HMul %d, Hits %d, %d, %d, "
729  "NDF %d, Chi2 %6.2f, T0 %6.2f, Tt %6.4f ",
730  pTrk,
731  pTrk->GetNofHits(),
732  iHit0,
733  iHit1,
734  iHit,
735  pTrk->GetNDF(),
736  pTrk->GetChiSq(),
737  pTrk->GetTime(),
738  dTt);
739  PrintStatus((char*) "<Res> ");
740  LOG(debug1) << " Match loop status: NCand " << iNCand << ", iDet "
741  << iDet;
742 
743  /* live display insert
744  if(gLogger->IsLogNeeded(debug3)) // update event display, if initialized
745  {
746  Int_t ii;
747  CbmEvDisTracks* fDis = CbmEvDisTracks::Instance();
748  if(NULL != fDis) {
749  //gEve->Redraw3D(kTRUE);
750  //gEve->FullRedraw3D();
751  fiNtrks=0;
752  for(Int_t iTr=0; iTr<fTracks.size(); iTr++){
753  if(fTracks[iTr]==NULL) continue;
754  CbmTofTracklet* pTrkDis = new((*fTofTracks)[fiNtrks++]) CbmTofTracklet (*fTracks[iTr]);
755  }
756  fDis->Exec("");
757  }
758  cout << " fDis "<<fDis<<" with "<<fiNtrks<<" tracks, to continue type 0 ! "<<endl;
759  scanf("%d",&ii);
760  } // end of live display
761  */
762  } // end of while(iNCand>0)
763  } // detector loop (propagate) end
764  } // iSt1 while condition end
765  } // iSt0 while condition end
766  //fTracks->Compress();
767  //fTofTracks = fTracks;
768 
769  // copy fTracks -> fTofTracks / fOutTracks
770 
771 
772  fiNtrks = 0;
773  for (size_t iTr = 0; iTr < fTracks.size(); iTr++) {
774  if (fTracks[iTr] == NULL) continue;
775  if (fTracks[iTr]->GetNofHits() < 3)
776  continue; // request minimum number of hits (3)
777  if (fTracks[iTr]->GetChiSq() > fChiMaxAccept)
778  continue; // request minimum ChiSq (3)
779  CbmTofTracklet* pTrk =
780  new ((*fTofTracks)[fiNtrks++]) CbmTofTracklet(*fTracks[iTr]);
781 
782  if (gLogger->IsLogNeeded(fair::Severity::debug)) {
783  LOG(info) << "Found Trkl " << iTr << ", ";
784  pTrk->PrintInfo();
785  }
786  for (Int_t iHit = 0; iHit < pTrk->GetNofHits(); iHit++) { // mark used Hit
787  CbmTofHit* pHit = (CbmTofHit*) fHits->At(pTrk->GetHitIndex(iHit));
788  pHit->SetFlag(pHit->GetFlag() + 100.);
789  LOG(debug) << Form(" hit %d at %d flagged to %d ",
790  iHit,
791  pTrk->GetHitIndex(iHit),
792  pHit->GetFlag());
793  }
794  }
795  PrintStatus((char*) "<D> Final result");
796 
797  for (size_t iTr = 0; iTr < fTracks.size(); iTr++) {
798  if (fTracks[iTr] == NULL) continue;
799  fTracks[iTr]->Delete();
800  //delete fTracks[iTr];
801  LOG(debug) << Form("<I> TofTracklet %lu, %p deleted", iTr, fTracks[iTr]);
802  }
803  fTracks.resize(0); //cleanup
804  // fFindTracks->PrintSetup();
805  return 0;
806 } // DoFind end
807 
808 void CbmTofTrackFinderNN::TrklSeed(Int_t iHit) {
809  CbmTofHit* pHit = (CbmTofHit*) fHits->At(iHit);
810  // Int_t iSmType = CbmTofAddress::GetSmType( pHit->GetAddress() & DetMask ); (VF) not used
811  Int_t iAddr = (pHit->GetAddress() & DetMask);
812  //Int_t iDet = fFindTracks->GetTypeStation(iSmType);
813  Int_t iDet = fFindTracks->GetStationOfAddr(iAddr);
814  if (iDet == fFindTracks->GetNofStations())
815  return; // hit not in tracking setup
816  for (Int_t iDet1 = 0; iDet1 < iDet; iDet1++) { // build new seeds
817  for (Int_t iHit1 = 0; iHit1 < fHits->GetEntries();
818  iHit1++) { // loop over previous Hits
819  CbmTofHit* pHit1 = (CbmTofHit*) fHits->At(iHit1);
820  // Int_t iSmType1 = CbmTofAddress::GetSmType( pHit1->GetAddress() & DetMask ); (VF) not used
821  Int_t iAddr1 = (pHit1->GetAddress() & DetMask);
822  if (iAddr1 == iAddr) continue;
823  //if (iSmType1 == fFindTracks->GetStationType(iDet1)) { // generate candidate for new track seed
824  if (iAddr1
825  == fFindTracks->GetAddrOfStation(
826  iDet1)) { // generate candidate for new track seed
827  Int_t iChId1 = pHit1->GetAddress();
828  CbmTofCell* fChannelInfo1 = fDigiPar->GetCell(iChId1);
829  Int_t iCh1 = CbmTofAddress::GetChannelId(iChId1);
830  Double_t hitpos1[3] = {3 * 0.};
831  Double_t hitpos1_local[3] = {3 * 0.};
832  Double_t dSizey1 = 1.;
833 
834  if (NULL == fChannelInfo1) {
835  LOG(debug)
836  << "CbmTofTrackFinderNN::DoFindp: Invalid Channel Pointer for ChId "
837  << Form(" 0x%08x ", iChId1) << ", Ch " << iCh1;
838  // continue;
839  } else {
840  /*TGeoNode *fNode1=*/ // prepare global->local trafo
841  gGeoManager->FindNode(fChannelInfo1->GetX(),
842  fChannelInfo1->GetY(),
843  fChannelInfo1->GetZ());
844  hitpos1[0] = pHit1->GetX();
845  hitpos1[1] = pHit1->GetY();
846  hitpos1[2] = pHit1->GetZ();
847  /*TGeoNode* cNode1=*/gGeoManager->GetCurrentNode();
848  gGeoManager->MasterToLocal(hitpos1, hitpos1_local);
849  dSizey1 = fChannelInfo1->GetSizey();
850  }
851  Double_t dDT = pHit->GetTime() - pHit1->GetTime();
852  Double_t dLz = pHit->GetZ() - pHit1->GetZ();
853  Double_t dTx = (pHit->GetX() - pHit1->GetX()) / dLz;
854  Double_t dTy = (pHit->GetY() - pHit1->GetY()) / dLz;
855  Int_t iUsed = HitUsed(iHit1);
856  LOG(debug1) << Form(
857  "<ISeed> TofTracklet %d, Hits %d, %d, used %d check, add = "
858  "0x%08x,0x%08x - DT %6.2f, Tx %6.2f Ty %6.2f ",
859  fiNtrks,
860  iHit,
861  iHit1,
862  iUsed,
863  pHit->GetAddress(),
864  pHit1->GetAddress(),
865  dDT,
866  dTx,
867  dTy);
868  if (TMath::Abs(hitpos1_local[1]) < dSizey1 * fPosYMaxScal)
869  if (TMath::Abs(dDT / dLz) < fMaxTofTimeDifference
870  && TMath::Abs(dTx - fTxMean) < fTxLIM
871  && TMath::Abs(dTy - fTyMean) < fTyLIM && iUsed == 0) {
872  // CbmTofTracklet* pTrk = new((*fTracks)[++fiNtrks]) CbmTofTracklet(); // generate new track seed
873  CbmTofTracklet* pTrk = new CbmTofTracklet();
874  ++fiNtrks;
875  fTracks.push_back(pTrk);
876  pTrk->SetTofHitIndex(iHit1, iAddr1, pHit1); // store Hit index
877  //Int_t NTrks1=fvTrkMap[iHit1].size()+1;
878  //fvTrkMap[iHit1].insert(std::pair<CbmTofTracklet*,Int_t>(pTrk,NTrks1));
879  fvTrkVec[iHit1].push_back(pTrk);
880 
881  pTrk->AddTofHitIndex(iHit, iAddr, pHit); // store 2. Hit index
882  //Int_t NTrks=fvTrkMap[iHit].size()+1;
883  //fvTrkMap[iHit].insert(std::pair<CbmTofTracklet*,Int_t>(pTrk,NTrks));
884  fvTrkVec[iHit].push_back(pTrk);
885 
886  pTrk->SetTime(
887  pHit->GetTime()); // define reference time from 2. plane
888  Double_t dR = pHit->GetR() - pHit1->GetR();
889  Double_t dTt = 1. / 30.; // assume speed of light: 1 / 30 cm/ns
890  // if( 0 == iSmType) pHit1->SetTime(pHit->GetTime() - dTt * dR);
891  dTt = (pHit->GetTime() - pHit1->GetTime()) / dR;
892  pTrk->SetTt(dTt);
893  pTrk->UpdateT0();
894 
895  CbmTofTrackletParam* tPar = pTrk->GetTrackParameter();
896  tPar->SetX(pHit->GetX()); // fill TrackParam
897  tPar->SetY(pHit->GetY()); // fill TrackParam
898  tPar->SetZ(pHit->GetZ()); // fill TrackParam
899  tPar->SetLz(dLz);
900  tPar->SetTx(dTx);
901  tPar->SetTy(dTy);
902  LOG(debug) << Form("<DSeed> TofTracklet %d, Hits %d, %d add "
903  "initialized, add = 0x%08x,0x%08x ",
904  fiNtrks,
905  iHit,
906  iHit1,
907  pHit->GetAddress(),
908  pHit1->GetAddress());
909  PrintStatus((char*) "after DSeed");
910  }
911  }
912  } // hit loop end
913  } // Station loop end
914 }
915 
916 Int_t CbmTofTrackFinderNN::HitUsed(Int_t iHit) {
917  // CbmTofHit* pHit = (CbmTofHit*) fHits->At( iHit ); (VF) not used
918  //Int_t iSmType = CbmTofAddress::GetSmType( pHit->GetAddress() & DetMask );
919  //Int_t iDet = fFindTracks->GetTypeStation(iSmType);
920  Int_t iUsed = 0;
921 
922  // LOG(debug1)<<"CbmTofTrackFinderNN::HitUsed of Hind "<<iHit<<", TrkMap.size "<<fvTrkMap[iHit].size()
923  LOG(DEBUG4) << "CbmTofTrackFinderNN::HitUsed of Hind " << iHit
924  << ", TrkVec.size " << fvTrkVec[iHit].size();
925  /*
926  for ( std::map<CbmTofTracklet*,Int_t>::iterator it=fvTrkMap[iHit].begin(); it != fvTrkMap[iHit].end(); it++){
927  if(it->first->GetNofHits() > 2) return iUsed=1;
928  }
929  */
930  if (fvTrkVec[iHit].size() > 0) {
931  if (fvTrkVec[iHit][0]->GetNofHits() > 2) iUsed = 1;
932  }
933  return iUsed;
934 }
935 
936 void CbmTofTrackFinderNN::UpdateTrackList(Int_t iTrk) {
937  CbmTofTracklet* pTrk = (CbmTofTracklet*) fTracks[iTrk];
938  UpdateTrackList(pTrk);
939 }
940 
941 void CbmTofTrackFinderNN::UpdateTrackList(CbmTofTracklet* pTrk) {
942  for (Int_t iHit = 0; iHit < pTrk->GetNofHits();
943  iHit++) { // loop over Tracklet Hits
944  Int_t iHitInd = pTrk->GetHitIndex(iHit); // Hit index in fHits
945  //Int_t NTrks=fvTrkMap[iHitInd].size(); // Number of tracks containing this hit
946  Int_t NTrks =
947  fvTrkVec[iHitInd].size(); // Number of tracks containing this hit
948  Int_t iAddr = (pTrk->GetTofHitPointer(iHit)->GetAddress() & DetMask);
949  if (iAddr == fFindTracks->GetBeamCounter())
950  continue; // keep all tracklets from common beam reference counter
951 
952  // Int_t iSmType = CbmTofAddress::GetSmType( iAddr ); (VF) not used
953  //if(iSmType==0) continue; // keep all tracklets with common target faked hit
954 
955  if (NTrks == 0)
956  LOG(fatal) << "UpdateTrackList NTrks=0 for event "
957  << fFindTracks->GetEventNumber() << ", pTrk " << pTrk
958  << ", iHit " << iHit;
959  if (NTrks > 0) {
960  //PrintStatus((char*)"UpdateTrackList::cleanup1");
961  Int_t iterClean = 1;
962  while (iterClean > 0) {
963  LOG(debug2) << " <D1> UpdateTrackList for Trk " << pTrk
964  << Form(
965  ", %d.Hit(%d) at ind %d with %d(%d) registered tracks",
966  iHit,
967  pTrk->GetNofHits(),
968  iHitInd,
969  (int) fvTrkVec[iHitInd].size(),
970  NTrks);
971  //if(fvTrkVec[iHitInd].size()==1) break;
972  for (std::vector<CbmTofTracklet*>::iterator iT =
973  fvTrkVec[iHitInd].begin();
974  iT != fvTrkVec[iHitInd].end();
975  iT++) {
976  iterClean = 0;
977  if (!Active(*iT)) break; // check whether tracklet is still active
978  LOG(debug2) << " <D2> process Trk " << *iT << " with "
979  << (*iT)->GetNofHits() << " hits";
980 
981  for (Int_t iH = 0; iH < (*iT)->GetNofHits(); iH++) {
982  if (!Active(*iT)) break; // check whether tracklet is still active
983  Int_t iHi = (*iT)->GetTofHitIndex(iH);
984  LOG(debug2) << " <D3> process Hit " << iH << " at index " << iHi;
985  Int_t iAddri =
986  ((*iT)->GetTofHitPointer(iH)->GetAddress() & DetMask);
987  LOG(debug2)
988  << " --- iHitInd " << iHitInd << "(" << fvTrkVec.size()
989  << "), size " << fvTrkVec[iHitInd].size() << " - iH " << iH << "("
990  << (*iT)->GetNofHits() << "), iHi " << iHi << " Hi vec size "
991  << fvTrkVec[iHi].size()
992  << Form(" poi %p, iTpoi %p, SmAddr 0x%08x, 0x%08x, 0x%08x ",
993  pTrk,
994  *iT,
995  (*iT)->GetTofHitPointer(iH)->GetAddress(),
996  iAddri,
997  fFindTracks->GetBeamCounter());
998 
999  if (iAddri == fFindTracks->GetBeamCounter()) {
1000  LOG(debug2) << " Hit in beam counter, continue ...";
1001  continue;
1002  }
1003  if (fvTrkVec[iHi].size() == 0) {
1004  LOG(fatal) << "CbmTofTrackFinderNN::UpdateTrackList no track "
1005  << " for hit " << iH << ", Hind " << iHi << ", size "
1006  << fvTrkVec[iHi].size();
1007  break;
1008  } else { // loop over tracks referenced by hit iHi
1009  for (std::vector<CbmTofTracklet*>::iterator it =
1010  fvTrkVec[iHi].begin();
1011  it != fvTrkVec[iHi].end();
1012  it++) {
1013  LOG(debug2)
1014  << " UpdateTrackList for pTrk " << pTrk << " <-> " << *iT
1015  << " <-> " << *it << ", clean " << iterClean << ", hit "
1016  << iHi << ", size " << fvTrkVec[iHi].size();
1017  if (*it != pTrk) {
1018  size_t iTr = 0;
1019  ;
1020  for (iTr = 0; iTr < fTracks.size(); iTr++) {
1021  if (fTracks[iTr] == *it) {
1022  LOG(debug2) << Form(" found track entry %p(%d) at %lu "
1023  "of iHi %d, pTrk %p",
1024  *it,
1025  (int) fvTrkVec[iHi].size(),
1026  iTr,
1027  iHi,
1028  pTrk);
1029  break;
1030  }
1031  }
1032 
1033  if (iTr == fTracks.size()) {
1034  LOG(fatal) << "CbmTofTrackFinderNN::UpdateTrackList: "
1035  "Invalid iTr for pTrk "
1036  << pTrk << ", iTr " << iTr << ", size "
1037  << fvTrkVec[iHi].size();
1038  break;
1039  }
1040 
1041  LOG(debug2) << Form("<D4> number of registered hits %3d at "
1042  "%p while keeping iHi = %d, pTrk at %p",
1043  (*it)->GetNofHits(),
1044  (*it),
1045  iHi,
1046  pTrk);
1047 
1048  CbmTofTracklet* pKill = *it;
1049  // remove track link registered for each associated hit to the track that is going to be removed
1050  for (Int_t iht = 0; iht < pKill->GetNofHits(); iht++) {
1051  Int_t iHI = pKill->GetHitIndex(iht);
1052  LOG(debug2) << Form("<D5> remove track link %p for hit iHi "
1053  "= %d, loop index %d: iHI = %3d ",
1054  pKill,
1055  iHi,
1056  iht,
1057  iHI);
1058 
1059  for (std::vector<CbmTofTracklet*>::iterator itt =
1060  fvTrkVec[iHI].begin();
1061  itt != fvTrkVec[iHI].end();
1062  itt++) {
1063  if ((*itt) == pTrk) continue;
1064  if ((*itt) == pKill) {
1065  LOG(debug2) << Form("<D6> remove track link %p for hit "
1066  "iHi = %d, iHI = %3d, #Trks %3d",
1067  pKill,
1068  iHi,
1069  iHI,
1070  (int) fvTrkVec[iHI].size());
1071  if (fvTrkVec[iHI].size() == 1) {
1072  LOG(debug2)
1073  << "<D6a> clear vector fvTrkVec for " << iHI;
1074  fvTrkVec[iHI].clear();
1075  // it =fvTrkVec[iHi].begin();
1076  break;
1077  } else {
1078  itt = fvTrkVec[iHI].erase(itt); // costly operation
1079  LOG(debug2) << "<D6b> reduce fvTrkVec size of " << iHI
1080  << " to " << fvTrkVec[iHI].size();
1081  break;
1082  }
1083  }
1084  }
1085  LOG(debug2) << Form("<D7> removed track link %p for hit "
1086  "iHi = %d, loop %d: iHI = %3d ",
1087  pKill,
1088  iHi,
1089  iht,
1090  iHI);
1091 
1092  // PrintStatus((char*)"UpdateTrackList::Remove1");
1093  } // loop on associated hits end
1094  pKill->Delete(); //delete tracklet *it;
1095  LOG(debug2) << "<D8> remove tracklet at pos " << iTr;
1096  fTracks.erase(fTracks.begin() + iTr);
1097  fiNtrks--;
1098 
1099  LOG(debug2)
1100  << "Erase1 for pTrk " << pTrk << ", at " << iTr << ", hit "
1101  << iHi << ", size " << fvTrkVec[iHi].size();
1102 
1103  PrintStatus((char*) "UpdateTrackList::Erase1");
1104 
1105  /*
1106  if(fvTrkVec[iHi].size() == 1) {
1107  fvTrkVec[iHi].clear();
1108  LOG(debug2) << " clear1 for pTrk "<<pTrk<<", hit "<<iHi<<", size "<<fvTrkVec[iHi].size()
1109  ;
1110  goto loopclean;
1111  }else{
1112  it=fvTrkVec[iHi].erase(it);
1113  //NTrks--;
1114  LOG(debug2) << " erase3 for "<<iTrk<<" at "<<pTrk<<", hit "<<iHi
1115  <<", size "<<fvTrkVec[iHi].size()<<", "<<NTrks
1116  ;
1117  }
1118 
1119  */
1120  //if(iHi == iHitInd) NTrks--;
1121  //PrintStatus((char*)"UpdateTrackList::cleanup2");
1122  iterClean = 2;
1123  break;
1124  } else { // *it==pTrk
1125  if (fvTrkVec[iHi].size() < 2) break;
1126  // if(pTrk == *iT) goto loopclean; //
1127  }
1128  } // end of loop over tracks referenced by hit iHi
1129  }
1130  }
1131  }
1132  }
1133  }
1134  }
1135 }
1136 
1137 void CbmTofTrackFinderNN::PrintStatus(char* cComment) {
1138  LOG(debug) << Form(
1139  "<PS %s> for fiNtrks = %d tracks out of %d fTracks.size() ",
1140  cComment,
1141  fiNtrks,
1142  (int) fTracks.size());
1143 
1144  for (size_t it = 0; it < fTracks.size(); it++) {
1145  CbmTofTracklet* pTrk = (CbmTofTracklet*) fTracks[it];
1146  if (NULL == pTrk) continue;
1147  if (pTrk->GetNofHits() < 2) { LOG(fatal) << "Invalid track found"; }
1148 
1149  TString sTrk = "";
1150  sTrk += Form(" Track %lu at %p, Hits: ", it, pTrk);
1151  for (Int_t ih = 0; ih < pTrk->GetNofHits(); ih++) {
1152  sTrk += Form(" %3d ", pTrk->GetHitIndex(ih));
1153  }
1154  sTrk += Form(", ChiSq %7.1f", pTrk->GetChiSq());
1155  sTrk += Form(", Tt %6.4f", pTrk->GetTt());
1156  LOG(debug) << sTrk;
1157  }
1158 
1159  for (size_t ih = 0; ih < fvTrkVec.size(); ih++) {
1160  CbmTofHit* pHit = (CbmTofHit*) fHits->At(ih);
1161  Int_t iAddr = (pHit->GetAddress() & DetMask);
1162  Int_t iSt = fFindTracks->GetStationOfAddr(iAddr);
1163  TString sTrk = "";
1164  sTrk += Form(" Hit %lu, A 0x%08x, St %d, T %6.2f, Tracks(%d): ",
1165  ih,
1166  pHit->GetAddress(),
1167  iSt,
1168  pHit->GetTime(),
1169  (int) fvTrkVec[ih].size());
1170  if (iSt < fFindTracks->GetNStations()) {
1171  for (size_t it = 0; it < fvTrkVec[ih].size(); it++) {
1172  CbmTofTracklet* pTrk = fvTrkVec[ih][it];
1173  sTrk += Form(" %p ", pTrk);
1174  }
1175  LOG(debug) << sTrk;
1176  }
1177  }
1178 }
1179 
1180 Bool_t CbmTofTrackFinderNN::Active(CbmTofTracklet* pCheck) {
1181  for (size_t it = 0; it < fTracks.size(); it++) {
1182  CbmTofTracklet* pTrk = (CbmTofTracklet*) fTracks[it];
1183  if (NULL == pTrk) continue;
1184  if (pCheck == pTrk) return kTRUE;
1185  }
1186  return kFALSE;
1187 }
1188 
1189 void CbmTofTrackFinderNN::Line3Dfit(CbmTofTracklet* pTrk) {
1190  TGraph2DErrors* gr = new TGraph2DErrors();
1191 
1192  // Fill the 2D graph
1193  // generate graph with the 3d points
1194  for (Int_t N = 0; N < pTrk->GetNofHits(); N++) {
1195  double x, y, z = 0;
1196  x = (pTrk->GetTofHitPointer(N))->GetX();
1197  y = (pTrk->GetTofHitPointer(N))->GetY();
1198  z = (pTrk->GetTofHitPointer(N))->GetZ();
1199  gr->SetPoint(N, x, y, z);
1200  double dx, dy, dz = 0.;
1201  dx = (pTrk->GetTofHitPointer(N))->GetDx();
1202  dy = (pTrk->GetTofHitPointer(N))->GetDy();
1203  dz = (pTrk->GetTofHitPointer(N))->GetDz(); //FIXME
1204  gr->SetPointError(N, dx, dy, dz);
1205  LOG(debug) << "Line3Dfit add N = " << N << ",\t" << pTrk->GetTofHitIndex(N)
1206  << ",\t" << x << ",\t" << y << ",\t" << z << ",\t" << dx << ",\t"
1207  << dy << ",\t" << dz;
1208  }
1209  // fit the graph now
1210  Double_t pStart[4] = {0., 0., 0., 0.};
1211  pStart[0] = pTrk->GetFitX(0.);
1212  pStart[1] = (pTrk->GetTrackParameter())->GetTx();
1213  pStart[2] = pTrk->GetFitY(0.);
1214  pStart[3] = (pTrk->GetTrackParameter())->GetTy();
1215  LOG(debug) << "Line3Dfit init: X0 " << pStart[0] << ", TX " << pStart[1]
1216  << ", Y0 " << pStart[2] << ", TY " << pStart[3];
1217 
1218  fMinuit.DoFit(gr, pStart);
1219  //gr->Draw("err p0");
1220  gr->Delete();
1221  Double_t* dRes;
1222  dRes = fMinuit.GetParFit();
1223  LOG(debug) << "Line3Dfit result: " << gMinuit->fCstatu << " : X0 " << dRes[0]
1224  << ", TX " << dRes[1] << ", Y0 " << dRes[2] << ", TY " << dRes[3]
1225  << ", Chi2DoF: " << fMinuit.GetChi2DoF();
1226 
1227  (pTrk->GetTrackParameter())->SetX(dRes[0]);
1228  (pTrk->GetTrackParameter())->SetY(dRes[2]);
1229  (pTrk->GetTrackParameter())->SetZ(0.);
1230  (pTrk->GetTrackParameter())->SetTx(dRes[1]);
1231  (pTrk->GetTrackParameter())->SetTy(dRes[3]);
1232  (pTrk->GetTrackParameter())->SetQp(1.); // FIXME
1233  // pTrk->SetChiSq(fMinuit.GetChi2DoF());
1234  pTrk->SetChiSq(
1235  fMinuit.GetChi2DoF()
1236  / pTrk
1237  ->GetNofHits()); // empirical to equilibrate bias on hit multiplicity!!!
1238 }
1239 
1240 ClassImp(CbmTofTrackFinderNN)
CbmTofDigiBdfPar.h
CbmTofCell::GetZ
Double_t GetZ() const
Definition: CbmTofCell.h:38
CbmHit::GetZ
Double_t GetZ() const
Definition: CbmHit.h:70
CbmTofFindTracks::GetEventNumber
Int_t GetEventNumber() const
Definition: CbmTofFindTracks.h:125
CbmTofFindTracks::GetMinNofHits
Int_t GetMinNofHits() const
Definition: CbmTofFindTracks.h:107
LKFMinuit::GetParFit
double * GetParFit()
Definition: LKFMinuit.h:37
CbmTofTrackFinderNN::fMaxTofTimeDifference
Double_t fMaxTofTimeDifference
Definition: CbmTofTrackFinderNN.h:88
CbmTofTracklet::GetTrackParameter
CbmTofTrackletParam * GetTrackParameter()
Definition: CbmTofTracklet.h:96
CbmTofTrackletParam::SetLz
void SetLz(Double_t lz)
Definition: CbmTofTrackletParam.h:62
CbmPixelHit::GetX
Double_t GetX() const
Definition: CbmPixelHit.h:83
CbmTofTracklet::GetHitIndex
Int_t GetHitIndex(Int_t ind) const
Definition: CbmTofTracklet.h:51
CbmTofTrackFinderNN::UpdateTrackList
void UpdateTrackList(Int_t iTrk)
Definition: CbmTofTrackFinderNN.cxx:936
CbmTofTracklet::PrintInfo
virtual void PrintInfo()
Definition: CbmTofTracklet.cxx:541
CbmTofTrackFinderNN::Active
Bool_t Active(CbmTofTracklet *pTrk)
Definition: CbmTofTrackFinderNN.cxx:1180
CbmTofTracklet::GetNofHits
Int_t GetNofHits() const
Definition: CbmTofTracklet.h:48
CbmTofTrackletParam::SetTy
void SetTy(Double_t ty)
Definition: CbmTofTrackletParam.h:64
CbmTofTracklet::SetTofHitIndex
void SetTofHitIndex(Int_t tofHitIndex, Int_t iDet, CbmTofHit *pHit)
Definition: CbmTofTracklet.h:138
CbmTofTrackletParam::GetTy
Double_t GetTy() const
Definition: CbmTofTrackletParam.h:53
CbmPixelHit::GetY
Double_t GetY() const
Definition: CbmPixelHit.h:84
CbmTofTrackletParam::GetTx
Double_t GetTx() const
Definition: CbmTofTrackletParam.h:52
CbmTofTracklet::SetChiSq
void SetChiSq(Double_t chiSq)
Definition: CbmTofTracklet.h:247
CbmTofTrackFinderNN::fFitter
CbmTofTrackFitter * fFitter
Definition: CbmTofTrackFinderNN.h:85
CbmTofTracklet::GetStationHitIndex
Int_t GetStationHitIndex(Int_t iSm) const
Definition: CbmTofTracklet.h:65
LKFMinuit.h
CbmTofTrackFinderNN::~CbmTofTrackFinderNN
virtual ~CbmTofTrackFinderNN()
Destructor.
Definition: CbmTofTrackFinderNN.cxx:69
CbmTofTrackFinderNN::fTracks
std::vector< CbmTofTracklet * > fTracks
Definition: CbmTofTrackFinderNN.h:99
CbmTofTracklet::RemoveTofHitIndex
void RemoveTofHitIndex(Int_t, Int_t iDet, CbmTofHit *, Double_t)
Definition: CbmTofTracklet.h:205
CbmTofTracklet::GetTofHitIndex
Int_t GetTofHitIndex(Int_t ind) const
Definition: CbmTofTracklet.h:71
CbmTofTracklet
Provides information on attaching a TofHit to a TofTrack.
Definition: CbmTofTracklet.h:25
gr
Char_t * gr
Definition: CbmEvDisTracks.cxx:289
CbmTofTrackFinderNN::fHits
TClonesArray * fHits
Definition: CbmTofTrackFinderNN.h:82
CbmTofFindTracks::GetSigT
Double_t GetSigT() const
Definition: CbmTofFindTracks.h:128
CbmTofTrackFinderNN::fFindTracks
CbmTofFindTracks * fFindTracks
Definition: CbmTofTrackFinderNN.h:86
CbmTofTrackFinderNN.h
CbmTofTrackletParam
Definition: CbmTofTrackletParam.h:27
CbmTofTrackletParam::SetY
void SetY(Double_t y)
Definition: CbmTofTrackletParam.h:60
LKFMinuit
Definition: LKFMinuit.h:19
CbmTofFindTracks::GetNofStations
Int_t GetNofStations()
Definition: CbmTofFindTracks.h:97
CbmTofTracklet::GetTex
Double_t GetTex(CbmTofHit *pHit)
Definition: CbmTofTracklet.cxx:155
CbmTofDigiPar.h
CbmMatch.h
CbmTofFindTracks::GetSigX
Double_t GetSigX() const
Definition: CbmTofFindTracks.h:129
CbmTofTracklet::Dist3D
virtual Double_t Dist3D(CbmTofHit *pHit0, CbmTofHit *pHit1)
Definition: CbmTofTracklet.cxx:557
CbmTofTrackletParam::ToString
std::string ToString() const
Return string representation of class.
Definition: CbmTofTrackletParam.h:118
CbmTofFindTracks::GetTtTarg
Double_t GetTtTarg() const
Definition: CbmTofFindTracks.h:126
CbmTofDigiPar::GetCell
CbmTofCell * GetCell(Int_t i)
Definition: CbmTofDigiPar.h:48
CbmTofTracklet.h
CbmTofTrackFinderNN::fTxLIM
Double_t fTxLIM
Definition: CbmTofTrackFinderNN.h:89
fTofTracks
TClonesArray * fTofTracks
Definition: CbmHadronAnalysis.cxx:51
CbmTofGeoHandler.h
CbmTofTrackletParam::GetX
Double_t GetX() const
Definition: CbmTofTrackletParam.h:48
CbmTofTracklet::GetMatChi2
virtual Double_t GetMatChi2(Int_t iSm)
Definition: CbmTofTracklet.cxx:116
CbmTofFindTracks::Instance
static CbmTofFindTracks * Instance()
Definition: CbmTofFindTracks.h:65
CbmTofDetectorId_v14a.h
CbmTofFindTracks::GetBeamCounter
Int_t GetBeamCounter() const
Definition: CbmTofFindTracks.h:124
CbmTofTrackFinderNN::Init
void Init()
Inherited from CbmTofTrackFinder.
Definition: CbmTofTrackFinderNN.cxx:104
CbmTofTracklet::GetNDF
Int_t GetNDF() const
Definition: CbmTofTracklet.h:128
CbmTofTrackFinderNN::fOutTracks
TClonesArray * fOutTracks
Definition: CbmTofTrackFinderNN.h:83
CbmTofTrackFinderNN
Definition: CbmTofTrackFinderNN.h:20
CbmTofTrackFinderNN::fiNtrks
Int_t fiNtrks
Definition: CbmTofTrackFinderNN.h:84
CbmTofTrackFinderNN::HitUsed
Int_t HitUsed(Int_t iHit)
Definition: CbmTofTrackFinderNN.cxx:916
CbmTofTrackFinderNN::TrklSeed
void TrklSeed(Int_t iHit)
Definition: CbmTofTrackFinderNN.cxx:808
CbmTofTrackletParam::SetTx
void SetTx(Double_t tx)
Definition: CbmTofTrackletParam.h:63
CbmTofTracklet::SetTime
void SetTime(Double_t val)
Definition: CbmTofTracklet.h:240
CbmTofCell::GetX
Double_t GetX() const
Definition: CbmTofCell.h:36
CbmTofDetectorId_v12b.h
CbmTofTrackFinderNN::operator=
CbmTofTrackFinderNN & operator=(const CbmTofTrackFinderNN &fSource)
Definition: CbmTofTrackFinderNN.cxx:97
CbmHit::GetTime
Double_t GetTime() const
Definition: CbmHit.h:75
CbmTofCell
Definition: CbmTofCell.h:8
CbmTofAddress::GetChannelId
static Int_t GetChannelId(UInt_t address)
Definition: CbmTofAddress.h:96
CbmHit::GetAddress
Int_t GetAddress() const
Definition: CbmHit.h:73
CbmTofTrackFinderNN::fPosYMaxScal
Double_t fPosYMaxScal
Definition: CbmTofTrackFinderNN.h:95
CbmTofTracklet::SetTt
void SetTt(Double_t val)
Definition: CbmTofTracklet.h:241
CbmTofTrackFinderNN::fTyLIM
Double_t fTyLIM
Definition: CbmTofTrackFinderNN.h:90
CbmTofTracklet::GetTt
Double_t GetTt() const
Definition: CbmTofTracklet.h:54
CbmTofFindTracks.h
CbmTofAddress.h
CbmTofTracklet::GetT0
Double_t GetT0() const
Definition: CbmTofTracklet.h:53
CbmTofTrackFinderNN::CbmTofTrackFinderNN
CbmTofTrackFinderNN()
Constructor.
Definition: CbmTofTrackFinderNN.cxx:51
CbmTofHit::GetR
Double_t GetR() const
Definition: core/data/tof/CbmTofHit.h:94
CbmTofTrackFinderNN::fvTrkVec
std::vector< std::vector< CbmTofTracklet * > > fvTrkVec
Definition: CbmTofTrackFinderNN.h:102
ClassImp
ClassImp(CbmConverterManager) InitStatus CbmConverterManager
Definition: CbmConverterManager.cxx:12
CbmTofTracklet::ReplaceTofHitIndex
void ReplaceTofHitIndex(Int_t tofHitIndex, Int_t iDet, CbmTofHit *pHit, Double_t chi2)
Definition: CbmTofTracklet.h:191
CbmTofTrackletParam::GetZ
Double_t GetZ() const
Definition: CbmTofTrackletParam.h:50
CbmTofTrackFinderNN::fMinuit
static LKFMinuit fMinuit
Definition: CbmTofTrackFinderNN.h:96
CbmTofTrackFinderNN::DoFind
Int_t DoFind(TClonesArray *fTofHits, TClonesArray *fTofTracks)
Definition: CbmTofTrackFinderNN.cxx:129
CbmTofTrackFinderNN::Line3Dfit
static void Line3Dfit(CbmTofTracklet *pTrk)
Definition: CbmTofTrackFinderNN.cxx:1189
CbmTofCell.h
CbmTofTrackFinderNN::PrintStatus
void PrintStatus(char *cComm)
Definition: CbmTofTrackFinderNN.cxx:1137
CbmTofTracklet::AddTofHitIndex
void AddTofHitIndex(Int_t tofHitIndex, Int_t iDet, CbmTofHit *pHit)
Definition: CbmTofTracklet.h:167
CbmTofHit::GetFlag
Int_t GetFlag() const
Definition: core/data/tof/CbmTofHit.h:91
LKFMinuit::DoFit
int DoFit(TGraph2DErrors *gr, double pStart[])
Definition: LKFMinuit.cxx:28
CbmTofCell::GetSizey
Double_t GetSizey() const
Definition: CbmTofCell.h:41
CbmTofTrackFinderNN::fChiMaxAccept
Double_t fChiMaxAccept
Definition: CbmTofTrackFinderNN.h:94
CbmTofFindTracks::GetSigY
Double_t GetSigY() const
Definition: CbmTofFindTracks.h:130
CbmTofTracklet::UpdateT0
Double_t UpdateT0()
Definition: CbmTofTracklet.cxx:186
CbmEvDisTracks.h
CbmTofTracklet::GetTime
Double_t GetTime() const
Definition: CbmTofTracklet.h:112
CbmHit::SetTime
void SetTime(Double_t time)
Definition: CbmHit.h:84
CbmTofTrackletParam::SetX
void SetX(Double_t x)
Definition: CbmTofTrackletParam.h:59
LKFMinuit::GetChi2DoF
double GetChi2DoF()
Definition: LKFMinuit.h:39
CbmTofTrackFinderNN::fTyMean
Double_t fTyMean
Definition: CbmTofTrackFinderNN.h:92
CbmTofTracklet::GetFitX
Double_t GetFitX(Double_t Z)
Definition: CbmTofTracklet.cxx:516
CbmTofTrackFinderNN::fTxMean
Double_t fTxMean
Definition: CbmTofTrackFinderNN.h:91
CbmTofTrackFinderNN::fDigiPar
CbmTofDigiPar * fDigiPar
Definition: CbmTofTrackFinderNN.h:87
CbmTofTrackletParam.h
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmTofDigiPar
Definition: CbmTofDigiPar.h:18
CbmTofTracklet::GetFitY
Double_t GetFitY(Double_t Z)
Definition: CbmTofTracklet.cxx:520
CbmTofTracklet::GetChiSq
Double_t GetChiSq() const
Definition: CbmTofTracklet.h:127
y
Double_t y
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmTofFindTracks::GetAddrOfStation
Int_t GetAddrOfStation(Int_t iVal)
Definition: CbmTofFindTracks.h:116
LKFMinuit::Initialize
int Initialize()
Definition: LKFMinuit.cxx:16
DetMask
const Int_t DetMask
Definition: CbmTofTrackFinderNN.cxx:48
CbmTofTrackletParam::GetY
Double_t GetY() const
Definition: CbmTofTrackletParam.h:49
CbmTofHit
Definition: core/data/tof/CbmTofHit.h:26
CbmTofCell::GetY
Double_t GetY() const
Definition: CbmTofCell.h:37
CbmTofTrackletParam::SetZ
void SetZ(Double_t z)
Definition: CbmTofTrackletParam.h:61
CbmTofAddress::GetUniqueAddress
static UInt_t GetUniqueAddress(UInt_t Sm, UInt_t Rpc, UInt_t Channel, UInt_t Side=0, UInt_t SmType=0)
Definition: CbmTofAddress.h:124
CbmTofFindTracks::GetNStations
Int_t GetNStations() const
Definition: CbmTofFindTracks.h:108
CbmTofTrackFitter.h
CbmTofTracklet::GetTofHitPointer
CbmTofHit * GetTofHitPointer(Int_t ind)
Definition: CbmTofTracklet.h:72
CbmTofFindTracks::SetStation
void SetStation(Int_t iVal, Int_t iModType, Int_t iModId, Int_t iRpcId)
Definition: CbmTofFindTracks.cxx:2338
CbmTofTrackFinderNN::fSIGLIM
Double_t fSIGLIM
Definition: CbmTofTrackFinderNN.h:93
CbmTofFindTracks::GetStationOfAddr
Int_t GetStationOfAddr(Int_t iAddr)
Definition: CbmTofFindTracks.cxx:2358
CbmTofHit::SetFlag
void SetFlag(Int_t flag)
Definition: core/data/tof/CbmTofHit.h:108