CbmRoot
LitTrackFinderNN.cxx
Go to the documentation of this file.
1 
7 #include "LitTrackFinderNN.h"
8 #include "LitAddMaterial.h"
9 #include "LitExtrapolation.h"
10 #include "LitFiltration.h"
11 #include "LitMath.h"
12 #include "LitTrackSelection.h"
13 
14 #include <algorithm>
15 #include <iostream>
16 #include <limits>
17 #include <map>
18 using std::cout;
19 using std::endl;
20 using std::for_each;
21 using std::map;
22 using std::numeric_limits;
23 
25 
27  : fTracks()
28  , fHitData()
29  , fUsedHitsSet()
30  , fUsedSeedsSet()
31  , fLayout()
32  , fNofIterations(0)
33  , fIteration(0)
34  , fMaxNofMissingHits()
35  , fPDG()
36  , fChiSqStripHitCut()
37  , fChiSqPixelHitCut()
38  , fSigmaCoef() {}
39 
41 
43  const vector<lit::parallel::LitScalPixelHit*>& hits,
44  const vector<lit::parallel::LitScalTrack*>& trackSeeds,
45  vector<lit::parallel::LitScalTrack*>& tracks) {
46  fTracks.clear();
47  fUsedSeedsSet.clear();
48  fUsedHitsSet.clear();
49  fHitData.SetNofStations(fLayout.GetNofStations());
50 
51  for (fIteration = 0; fIteration < fNofIterations; fIteration++) {
52  // cout << "LitTrackFinderNN::DoFind: iteration=" << fIteration << endl;
53  ArrangeHits(hits);
54  // cout << fHitData.ToString();
55 
56  InitTrackSeeds(trackSeeds);
57  FollowTracks();
58  SelectTracks();
59  RemoveHits();
60  CopyToOutput(tracks);
61 
62  for_each(fTracks.begin(), fTracks.end(), DeleteObject());
63  fTracks.clear();
64  fHitData.Clear();
65  }
66 
67  static int eventNo = 0;
68  cout << "LitTrackFinderNN::DoFind: eventNo=" << eventNo++ << endl;
69 }
70 
72  const vector<lit::parallel::LitScalPixelHit*>& hits) {
73  unsigned int nofHits = hits.size();
74  for (unsigned int iHit = 0; iHit < nofHits; iHit++) {
75  LitScalPixelHit* hit = hits[iHit];
76  if (fUsedHitsSet.find(hit->refId) != fUsedHitsSet.end()) { continue; }
77  fHitData.AddHit(hit);
78  }
79  fHitData.Arrange();
80 }
81 
83  const vector<lit::parallel::LitScalTrack*>& trackSeeds) {
84  unsigned int nofSeeds = trackSeeds.size();
85  for (unsigned int iTrack = 0; iTrack < nofSeeds; iTrack++) {
86  LitScalTrack* seed = trackSeeds[iTrack];
87  // Apply cuts here
88  if (std::fabs(seed->GetParamFirst().Qp) > 10)
89  continue; // momentum cut 0.1 GeV
90 
91  //if (fUsedSeedsSet.find(seed->GetPreviousTrackId()) != fUsedSeedsSet.end()) { continue; }
92  LitScalTrack* track = new LitScalTrack();
93  // FIXME if several iterations this procedure will be repeated!!
94  LitTrackParamScal par = seed->GetParamFirst();
95  PropagateVirtualStations(par);
96  track->SetParamFirst(par);
97  track->SetParamLast(par);
98  track->SetPreviousTrackId(seed->GetPreviousTrackId());
99  fTracks.push_back(track);
100  }
101 }
102 
104  LitTrackParamScal& par) {
105  unsigned char nofVirtualStations = fLayout.GetNofVirtualStations();
106  unsigned char nofSteps = (nofVirtualStations - 1) / 2;
107  for (unsigned char iStep = 0; iStep < nofSteps; iStep++) {
108  const LitVirtualStationScal& vsFront =
109  fLayout.GetVirtualStation(2 * iStep + 0);
110  const LitVirtualStationScal& vsMiddle =
111  fLayout.GetVirtualStation(2 * iStep + 1);
112  const LitVirtualStationScal& vsBack =
113  fLayout.GetVirtualStation(2 * iStep + 2);
114  if (vsFront.GetField().IsEmpty() || vsMiddle.GetField().IsEmpty()
115  || vsBack.GetField().IsEmpty()) {
116  LitLineExtrapolation(par, vsBack.GetZ());
117  } else {
119  vsBack.GetZ(),
120  vsFront.GetField(),
121  vsMiddle.GetField(),
122  vsBack.GetField());
123  }
124 
125  if (!vsFront.GetMaterial().IsEmpty()) {
126  fscal thickness = vsFront.GetMaterial().GetMaterial(par.X, par.Y);
127  if (thickness > 0) LitAddMaterial<fscal>(par, thickness);
128  }
129 
130  if (!vsMiddle.GetMaterial().IsEmpty()) {
131  fscal thickness = vsMiddle.GetMaterial().GetMaterial(par.X, par.Y);
132  if (thickness > 0) LitAddMaterial<fscal>(par, thickness);
133  }
134 
135  if (!vsBack.GetMaterial().IsEmpty()) {
136  fscal thickness = vsBack.GetMaterial().GetMaterial(par.X, par.Y);
137  if (thickness > 0) LitAddMaterial<fscal>(par, thickness);
138  }
139  }
140 }
141 
143  unsigned char stationId,
144  LitTrackParamScal& par) {
145  const LitStationScal& station = fLayout.GetStation(stationId);
146  unsigned char nofVirtualStations = station.GetNofVirtualStations();
147  if (nofVirtualStations == 1) {
148  const LitVirtualStationScal& vs = station.GetVirtualStation(0);
149  fscal z = vs.GetZ();
150  LitLineExtrapolation(par, z);
151  if (!vs.GetMaterial().IsEmpty()) {
152  fscal thickness = vs.GetMaterial().GetMaterial(par.X, par.Y);
153  if (thickness > 0) LitAddMaterial<fscal>(par, thickness);
154  }
155  } else {
156  unsigned char nofSteps = (nofVirtualStations - 1) / 2;
157  for (unsigned char iStep = 0; iStep < nofSteps; iStep++) {
158  const LitVirtualStationScal& vsFront =
159  station.GetVirtualStation(2 * iStep + 0);
160  const LitVirtualStationScal& vsMiddle =
161  station.GetVirtualStation(2 * iStep + 1);
162  const LitVirtualStationScal& vsBack =
163  station.GetVirtualStation(2 * iStep + 2);
164  if (vsFront.GetField().IsEmpty() || vsMiddle.GetField().IsEmpty()
165  || vsBack.GetField().IsEmpty()) {
166  LitLineExtrapolation(par, vsBack.GetZ());
167  } else {
169  vsBack.GetZ(),
170  vsFront.GetField(),
171  vsMiddle.GetField(),
172  vsBack.GetField());
173  }
174 
175  if (!vsFront.GetMaterial().IsEmpty()) {
176  fscal thickness = vsFront.GetMaterial().GetMaterial(par.X, par.Y);
177  if (thickness > 0) LitAddMaterial<fscal>(par, thickness);
178  }
179 
180  if (!vsMiddle.GetMaterial().IsEmpty()) {
181  fscal thickness = vsMiddle.GetMaterial().GetMaterial(par.X, par.Y);
182  if (thickness > 0) LitAddMaterial<fscal>(par, thickness);
183  }
184 
185  if (!vsBack.GetMaterial().IsEmpty()) {
186  fscal thickness = vsBack.GetMaterial().GetMaterial(par.X, par.Y);
187  if (thickness > 0) LitAddMaterial<fscal>(par, thickness);
188  }
189  }
190  }
191 }
192 
194  unsigned int nofTracks = fTracks.size();
195  for (unsigned int iTrack = 0; iTrack < nofTracks; iTrack++) {
196  LitScalTrack* track = fTracks[iTrack];
197 
198  unsigned char nofStations = fLayout.GetNofStations();
199  for (int iStation = 0; iStation < nofStations; iStation++) {
200  LitTrackParamScal par(track->GetParamLast());
201  // const LitStationScal& station = fLayout.GetStation(iStation);
202 
203  fscal zMin = fHitData.GetMinZPos(iStation);
204  LitLineExtrapolation(par, zMin);
205 
206  const vector<int>& bins = fHitData.GetZPosBins(iStation);
207  map<int, LitTrackParamScal> binParamMap;
208  vector<int>::const_iterator itBins;
209  for (itBins = bins.begin(); itBins != bins.end(); itBins++) {
210  binParamMap[*itBins] = LitTrackParamScal();
211  }
212 
213  // Extrapolate track parameters to each Z position in the map.
214  // This is done to improve calculation speed.
215  // In case of planar station only 1 track extrapolation is required,
216  // since all hits located at the same Z.
217  map<int, LitTrackParamScal>::iterator itMap;
218  for (itMap = binParamMap.begin(); itMap != binParamMap.end(); itMap++) {
219  (*itMap).second = par;
220  // fscal z = fHitData.GetZPosByBin(iStation, (*itMap).first);
221  LitTrackParamScal& tpar = (*itMap).second;
222  PropagateToStation(iStation, tpar);
223  // LitLineExtrapolation(tpar, z);
224  // if (!station.GetVirtualStation().GetMaterial().IsEmpty()) {
225  // fscal thickness = station.GetVirtualStation().GetMaterial().GetMaterial(tpar.X, tpar.Y);
226  // if (thickness > 0) LitAddMaterial<fscal>(tpar, thickness);
227  // }
228  }
229 
230  // Loop over hits
231  fscal minChiSq =
232  numeric_limits<fscal>::max(); // minimum chi-square of hit
233  const LitScalPixelHit* minHit =
234  NULL; // Pointer to hit with minimum chi-square
235  LitTrackParamScal minPar; // Track parameters for closest hit
236  const vector<LitScalPixelHit*>& hits = fHitData.GetHits(iStation);
237  unsigned int nofHits = hits.size();
238  for (unsigned int iHit = 0; iHit < nofHits; iHit++) {
239  const LitScalPixelHit* hit = hits[iHit];
240  int bin = fHitData.GetBinByZPos(iStation, hit->Z);
241  assert(binParamMap.find(bin) != binParamMap.end());
242  LitTrackParamScal tpar(binParamMap[bin]);
243 
244  // Check preliminary if hit is in the validation gate.
245  // This is done in order to speed up the algorithm.
246  // Based on the predicted track position (w/o KF update step)
247  // and maximum hit position error.
248  fscal maxErrX = fHitData.GetMaxErrX(iStation);
249  fscal devX =
250  fSigmaCoef[fIteration] * (sqrt(tpar.C0 + maxErrX * maxErrX));
251  fscal maxErrY = fHitData.GetMaxErrY(iStation);
252  fscal devY =
253  fSigmaCoef[fIteration] * (sqrt(tpar.C5 + maxErrY * maxErrY));
254  bool hitInside =
255  (hit->X < (tpar.X + devX)) && (hit->X > (tpar.X - devX))
256  && (hit->Y < (tpar.Y + devY)) && (hit->Y > (tpar.Y - devY));
257  if (!hitInside) continue;
258 
260  LitFiltration(tpar, *hit, chi);
261  bool hitInValidationGate = chi < fChiSqStripHitCut[fIteration];
262  if (
263  hitInValidationGate
264  && chi
265  < minChiSq) { // Check if hit is inside validation gate and closer to the track.
266  minChiSq = chi;
267  minHit = hit;
268  minPar = tpar;
269  }
270  }
271 
272  if (minHit != NULL) { // Check if hit was added
273  if (track->GetNofHits() == 0) track->SetParamFirst(minPar);
274  track->AddHit(minHit);
275  track->SetParamLast(minPar);
276  track->IncChiSq(minChiSq);
277  track->SetNDF(NDF(*track));
278  track->SetLastStationId(iStation);
279  } else { // Missing hit
280  track->SetNofMissingHits(track->GetNofMissingHits() + 1);
281  if (track->GetNofMissingHits() > fMaxNofMissingHits[fIteration]) {
282  break;
283  }
284  }
285  }
286  }
287 }
288 
290  unsigned int nofTracks = fTracks.size();
291  for (unsigned int iTrack = 0; iTrack < nofTracks; iTrack++) {
292  LitScalTrack* track = fTracks[iTrack];
293  track->IsGood(track->GetNofHits() > 0);
294  }
295  DoSelectSharedHits(fTracks);
296 }
297 
299  unsigned int nofTracks = fTracks.size();
300  for (unsigned int iTrack = 0; iTrack < nofTracks; iTrack++) {
301  LitScalTrack* track = fTracks[iTrack];
302  if (!track->IsGood()) { continue; }
303  for (unsigned int iHit = 0; iHit < track->GetNofHits(); iHit++) {
304  fUsedHitsSet.insert(track->GetHit(iHit)->refId);
305  }
306  }
307 }
308 
310  vector<lit::parallel::LitScalTrack*>& tracks) {
311  unsigned int nofTracks = fTracks.size();
312  for (unsigned int iTrack = 0; iTrack < nofTracks; iTrack++) {
313  LitScalTrack* track = fTracks[iTrack];
314  if (!track->IsGood()) { continue; }
315  fUsedSeedsSet.insert(track->GetPreviousTrackId());
316  tracks.push_back(new LitScalTrack(*track));
317  }
318 }
lit::parallel::LitTrackParamScal
LitTrackParam< fscal > LitTrackParamScal
Scalar version of LitTrackParam.
Definition: LitTrackParam.h:113
LitMath.h
Useful math functions.
fscal
float fscal
Definition: L1/vectors/P4_F32vec4.h:250
lit::parallel::LitScalPixelHit
Base class for scalar pixel hits.
Definition: LitScalPixelHit.h:31
LitExtrapolation.h
Functions for track parameters extrapolation.
lit::parallel::LitTrackFinderNN::PropagateToStation
void PropagateToStation(unsigned char stationId, LitTrackParamScal &par)
Definition: LitTrackFinderNN.cxx:142
lit::parallel::LitFiltration
void LitFiltration(LitTrackParam< T > &par, const LitPixelHit< T > &hit, T &chiSq)
Function implements Kalman filter update step for pixel hit.
Definition: LitFiltration.h:34
lit::parallel::LitTrackParam::X
T X
Definition: LitTrackParam.h:92
lit::parallel::LitTrackParam::Y
T Y
Definition: LitTrackParam.h:93
lit::parallel::LitTrackFinderNN::LitTrackFinderNN
LitTrackFinderNN()
Constructor.
Definition: LitTrackFinderNN.cxx:26
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
lit::parallel::LitScalTrack::GetHit
const LitScalPixelHit * GetHit(unsigned short index) const
Returns pointer to the hit.
Definition: LitScalTrack.h:82
lit::parallel::LitTrackFinderNN::SelectTracks
void SelectTracks()
Definition: LitTrackFinderNN.cxx:289
lit::parallel::LitTrackFinderNN::FollowTracks
void FollowTracks()
Follow tracks through detector.
Definition: LitTrackFinderNN.cxx:193
lit::parallel::LitMaterialGrid::GetMaterial
fscal GetMaterial(fscal x, fscal y) const
Return material thickness in silicon equivalent for (X, Y) position (scalar version).
Definition: LitMaterialGrid.h:91
LitAddMaterial.h
Functions for calculation of the material effects.
lit::parallel::LitTrackParam::C0
T C0
Definition: LitTrackParam.h:100
lit::parallel::LitScalPixelHit::Z
fscal Z
Definition: LitScalPixelHit.h:78
lit::parallel::LitScalTrack::SetLastStationId
void SetLastStationId(unsigned short lastStationId)
Set last station ID.
Definition: LitScalTrack.h:166
lit::parallel::LitStation
Detector station.
Definition: LitStation.h:31
lit::parallel::LitTrackFinderNN::PropagateVirtualStations
void PropagateVirtualStations(LitTrackParamScal &par)
Definition: LitTrackFinderNN.cxx:103
lit::parallel::LitVirtualStation
Virtual detector station which stores information needed for track propagation.
Definition: LitVirtualStation.h:31
lit::parallel::LitScalTrack
Scalar track data class.
Definition: LitScalTrack.h:33
lit::parallel::LitScalTrack::GetNofHits
unsigned short GetNofHits() const
Returns number of hits in track.
Definition: LitScalTrack.h:66
lit::parallel::DeleteObject
Functor class for convenient memory release.
Definition: LitUtils.h:51
lit::parallel::LitScalTrack::SetPreviousTrackId
void SetPreviousTrackId(unsigned short previousTrackId)
Sets previous trackId.
Definition: LitScalTrack.h:188
lit::parallel::LitFieldGrid::IsEmpty
bool IsEmpty() const
Check if field was set.
Definition: LitFieldGrid.h:193
lit::parallel::LitScalTrack::AddHit
void AddHit(const LitScalPixelHit *hit)
Adds hit to track.
Definition: LitScalTrack.h:60
lit::parallel::LitTrackParam< fscal >
lit::parallel::DoSelectSharedHits
void DoSelectSharedHits(vector< LitScalTrack * > &tracks)
Definition: LitTrackSelection.h:86
lit::parallel::LitScalTrack::IncChiSq
void IncChiSq(fscal dChiSq)
Increases chi square by dChiSq.
Definition: LitScalTrack.h:128
lit::parallel::LitTrackFinderNN::CopyToOutput
void CopyToOutput(vector< lit::parallel::LitScalTrack * > &tracks)
Copy tracks to output array.
Definition: LitTrackFinderNN.cxx:309
lit::parallel::LitScalTrack::GetNofMissingHits
unsigned short GetNofMissingHits() const
Returns number of missing hits.
Definition: LitScalTrack.h:146
LitTrackSelection.h
Implementation of the track selection algorithms.
lit::parallel::LitScalTrack::IsGood
bool IsGood() const
Returns true if track is good.
Definition: LitScalTrack.h:196
lit::parallel::LitScalTrack::SetParamFirst
void SetParamFirst(const LitTrackParamScal &param)
Sets first track parameter.
Definition: LitScalTrack.h:96
tracks
TClonesArray * tracks
Definition: Analyze_matching.h:17
lit::parallel::LitMaterialGrid::IsEmpty
bool IsEmpty() const
Check if material was set.
Definition: LitMaterialGrid.h:161
lit::parallel::LitRK4Extrapolation
void LitRK4Extrapolation(LitTrackParam< T > &par, T zOut, const LitFieldGrid &field1, const LitFieldGrid &field2, const LitFieldGrid &field3)
Definition: LitExtrapolation.h:73
lit::parallel::LitScalPixelHit::Y
fscal Y
Definition: LitScalPixelHit.h:73
lit::parallel::LitScalTrack::SetNofMissingHits
void SetNofMissingHits(unsigned short nofMissingHits)
Sets number of missing hits.
Definition: LitScalTrack.h:152
lit::parallel::LitScalPixelHit::X
fscal X
Definition: LitScalPixelHit.h:73
lit::parallel::LitScalTrack::GetPreviousTrackId
unsigned short GetPreviousTrackId() const
Return Previous track index.
Definition: LitScalTrack.h:182
lit::parallel::LitTrackFinderNN::DoFind
void DoFind(const vector< lit::parallel::LitScalPixelHit * > &hits, const vector< lit::parallel::LitScalTrack * > &trackSeeds, vector< lit::parallel::LitScalTrack * > &tracks)
Main function for track reconstruction.
Definition: LitTrackFinderNN.cxx:42
LitTrackFinderNN.h
Parallel implementation of the nearest neighbor tracking algorithm.
lit::parallel::LitVirtualStation::GetMaterial
const LitMaterialGrid & GetMaterial() const
Definition: LitVirtualStation.h:51
lit::parallel::LitTrackFinderNN::ArrangeHits
void ArrangeHits(const vector< lit::parallel::LitScalPixelHit * > &hits)
Definition: LitTrackFinderNN.cxx:71
LitFiltration.h
lit::parallel::LitTrackParam::C5
T C5
Definition: LitTrackParam.h:100
lit::parallel::LitVirtualStation::GetZ
T GetZ() const
Definition: LitVirtualStation.h:53
lit::parallel::LitScalTrack::GetParamFirst
const LitTrackParamScal & GetParamFirst() const
Returns first parameter of the track.
Definition: LitScalTrack.h:90
lit::parallel::LitStation::GetVirtualStation
const LitVirtualStation< T > & GetVirtualStation(unsigned char virtualStation) const
Return virtual station by index.
Definition: LitStation.h:65
lit::parallel::LitScalPixelHit::refId
unsigned short refId
Definition: LitScalPixelHit.h:77
lit::parallel::LitScalTrack::SetParamLast
void SetParamLast(const LitTrackParamScal &param)
Sets last track parameter.
Definition: LitScalTrack.h:110
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:60
lit::parallel::LitScalTrack::GetParamLast
const LitTrackParamScal & GetParamLast() const
Returns last parameter of the track.
Definition: LitScalTrack.h:104
hits
static vector< vector< QAHit > > hits
Definition: CbmTofHitFinderTBQA.cxx:114
lit::parallel::LitVirtualStation::GetField
const LitFieldGrid & GetField() const
Definition: LitVirtualStation.h:52
lit::parallel::LitTrackParam::Qp
T Qp
Definition: LitTrackParam.h:97
lit::parallel::NDF
unsigned short NDF(const LitScalTrack &track)
Returns number of degrees of freedom for the track.
Definition: LitMath.h:50
lit::parallel::LitTrackFinderNN::InitTrackSeeds
void InitTrackSeeds(const vector< lit::parallel::LitScalTrack * > &trackSeeds)
Initialize track seeds and copy to local array.
Definition: LitTrackFinderNN.cxx:82
lit::parallel::LitHitData::EPSILON
static const fscal EPSILON
Definition: LitHitData.h:227
lit::parallel::LitStation::GetNofVirtualStations
unsigned char GetNofVirtualStations() const
Return number of virtual stations.
Definition: LitStation.h:55
lit::parallel::LitLineExtrapolation
void LitLineExtrapolation(LitTrackParam< T > &par, T zOut)
Line track extrapolation for the field free regions.
Definition: LitExtrapolation.h:37
max
friend F32vec4 max(const F32vec4 &a, const F32vec4 &b)
Definition: L1/vectors/P4_F32vec4.h:36
lit::parallel::LitTrackFinderNN::~LitTrackFinderNN
virtual ~LitTrackFinderNN()
Destructor.
Definition: LitTrackFinderNN.cxx:40
lit::parallel::LitTrackFinderNN::RemoveHits
void RemoveHits()
Write already used hits to a used hits set.
Definition: LitTrackFinderNN.cxx:298
lit::parallel::LitScalTrack::SetNDF
void SetNDF(unsigned short NDF)
Sets number of degrees of freedom.
Definition: LitScalTrack.h:140