]>
Commit | Line | Data |
---|---|---|
8cc77df0 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | /* $Id$ */ | |
17 | ||
06ca6d7b | 18 | /// \class AliMUONTrackReconstructorK |
19 | ///////////////////////////////////// | |
20 | /// | |
21 | /// MUON track reconstructor using the kalman method | |
22 | /// | |
23 | /// This class contains as data: | |
24 | /// - the parameters for the track reconstruction | |
25 | /// | |
26 | /// It contains as methods, among others: | |
27 | /// - MakeTracks to build the tracks | |
28 | /// | |
8cc77df0 | 29 | //////////////////////////////////// |
30 | ||
8cc77df0 | 31 | #include "AliMUONTrackReconstructorK.h" |
8cc77df0 | 32 | #include "AliMUONConstants.h" |
33 | #include "AliMUONHitForRec.h" | |
ea94c18b | 34 | #include "AliMUONTrack.h" |
35 | #include "AliMUONTrackParam.h" | |
36 | #include "AliMUONTrackExtrap.h" | |
8cde4af5 | 37 | |
8cc77df0 | 38 | #include "AliLog.h" |
39 | ||
8cde4af5 | 40 | #include <Riostream.h> |
ea94c18b | 41 | #include <TMath.h> |
42 | #include <TMatrixD.h> | |
8cde4af5 | 43 | |
44 | /// \cond CLASSIMP | |
8cc77df0 | 45 | ClassImp(AliMUONTrackReconstructorK) // Class implementation in ROOT context |
8cde4af5 | 46 | /// \endcond |
8cc77df0 | 47 | |
ea94c18b | 48 | //************* Defaults parameters for reconstruction |
49 | const Bool_t AliMUONTrackReconstructorK::fgkRunSmoother = kTRUE; | |
8cc77df0 | 50 | |
8cc77df0 | 51 | |
7ec3b9cf | 52 | //__________________________________________________________________________ |
ea94c18b | 53 | AliMUONTrackReconstructorK::AliMUONTrackReconstructorK() |
54 | : AliMUONVTrackReconstructor() | |
8cc77df0 | 55 | { |
ea94c18b | 56 | /// Constructor for class AliMUONTrackReconstructorK |
57 | AliInfo("*** Tracking with Kalman Filter ***"); | |
8cc77df0 | 58 | } |
59 | ||
7ec3b9cf | 60 | //__________________________________________________________________________ |
ea94c18b | 61 | AliMUONTrackReconstructorK::~AliMUONTrackReconstructorK() |
8cc77df0 | 62 | { |
ea94c18b | 63 | /// Destructor |
64 | } | |
8cc77df0 | 65 | |
66 | //__________________________________________________________________________ | |
ea94c18b | 67 | void AliMUONTrackReconstructorK::MakeTrackCandidates() |
8cc77df0 | 68 | { |
ea94c18b | 69 | /// To make track candidates: |
70 | /// Start with segments station(1..) 4 or 5 then follow track in station 5 or 4. | |
71 | /// Good candidates are made of at least three hitForRec's. | |
72 | /// Keep only best candidates or all of them according to the flag fgkTrackAllTracks. | |
208f139e | 73 | TClonesArray *segments; |
ea94c18b | 74 | AliMUONTrack *track; |
75 | Int_t iCandidate = 0; | |
8cc77df0 | 76 | |
ea94c18b | 77 | AliDebug(1,"Enter MakeTrackCandidates"); |
8cc77df0 | 78 | |
ea94c18b | 79 | // Loop over stations(1..) 5 and 4 and make track candidates |
80 | for (Int_t istat=4; istat>=3; istat--) { | |
81 | ||
208f139e | 82 | // Make segments in the station |
83 | segments = MakeSegmentsInStation(istat); | |
ea94c18b | 84 | |
85 | // Loop over segments | |
86 | for (Int_t iseg=0; iseg<segments->GetEntriesFast(); iseg++) { | |
87 | AliDebug(1,Form("Making primary candidate(1..) %d",++iCandidate)); | |
88 | ||
89 | // Transform segments to tracks and put them at the end of fRecTracksPtr | |
90 | track = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack((AliMUONObjectPair*)((*segments)[iseg])); | |
91 | fNRecTracks++; | |
92 | ||
93 | // Recompute track parameters and covariances on station(1..) 5 using Kalman filter | |
94 | // (to make sure all tracks are treated in the same way) | |
95 | if (istat == 4) RetraceTrack(*track,kFALSE); | |
96 | ||
97 | // Printout for debuging | |
98 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) { | |
99 | cout<<endl<<"Track parameter covariances at first hit:"<<endl; | |
100 | ((AliMUONTrackParam*) (track->GetTrackParamAtHit()->First()))->GetCovariances().Print(); | |
101 | } | |
102 | ||
103 | // Look for compatible hitForRec(s) in the other station | |
104 | // Remove track if no hit found | |
105 | if (!FollowTrackInStation(*track,7-istat)) { | |
106 | fRecTracksPtr->Remove(track); | |
107 | fNRecTracks--; | |
108 | } | |
109 | ||
110 | } | |
111 | // delete the array of segments | |
7ec3b9cf | 112 | delete segments; |
ea94c18b | 113 | } |
114 | ||
115 | fRecTracksPtr->Compress(); // this is essential before checking tracks | |
116 | ||
117 | // Remove bad tracks | |
118 | Int_t currentNRecTracks = fNRecTracks; | |
119 | for (Int_t iRecTrack = 0; iRecTrack < currentNRecTracks; iRecTrack++) { | |
120 | ||
121 | track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iRecTrack); | |
122 | ||
123 | // Remove the track if the normalized chi2 is too high | |
124 | if (track->GetNormalizedChi2() > fgkSigmaToCutForTracking * fgkSigmaToCutForTracking) { | |
125 | fRecTracksPtr->Remove(track); | |
126 | fNRecTracks--; | |
127 | } | |
128 | ||
129 | } | |
130 | ||
131 | fRecTracksPtr->Compress(); // this is essential before checking tracks | |
132 | ||
133 | // Keep all different tracks or only the best ones as required | |
134 | if (fgkTrackAllTracks) RemoveIdenticalTracks(); | |
135 | else RemoveDoubleTracks(); | |
136 | ||
137 | AliDebug(1,Form("Number of good candidates = %d",fNRecTracks)); | |
138 | ||
8cc77df0 | 139 | } |
140 | ||
ea94c18b | 141 | //__________________________________________________________________________ |
142 | void AliMUONTrackReconstructorK::RetraceTrack(AliMUONTrack &trackCandidate, Bool_t resetSeed) | |
8cc77df0 | 143 | { |
ea94c18b | 144 | /// Re-run the kalman filter from the most downstream hit to the most uptream one |
145 | AliDebug(1,"Enter RetraceTrack"); | |
146 | ||
147 | AliMUONTrackParam* startingTrackParam = (AliMUONTrackParam*) trackCandidate.GetTrackParamAtHit()->Last(); | |
148 | ||
149 | // Reset the "seed" (= track parameters and their covariances at last hit) if required | |
150 | if (resetSeed) { | |
151 | // => Shift track parameters at the position of the last hit | |
152 | AliMUONHitForRec* hitForRecAtHit = startingTrackParam->GetHitForRecPtr(); | |
153 | startingTrackParam->SetNonBendingCoor(hitForRecAtHit->GetNonBendingCoor()); | |
154 | startingTrackParam->SetBendingCoor(hitForRecAtHit->GetBendingCoor()); | |
155 | ||
156 | // => Re-compute and reset track parameters covariances at last hit (as if the other hits did not exist) | |
157 | const TMatrixD& kParamCov = startingTrackParam->GetCovariances(); | |
158 | TMatrixD newParamCov(5,5); | |
159 | newParamCov.Zero(); | |
160 | // Non bending plane | |
161 | newParamCov(0,0) = hitForRecAtHit->GetNonBendingReso2(); | |
162 | newParamCov(1,1) = 100.*kParamCov(1,1); | |
163 | // Bending plane | |
164 | newParamCov(2,2) = hitForRecAtHit->GetBendingReso2(); | |
165 | newParamCov(3,3) = 100.*kParamCov(3,3); | |
166 | // Inverse bending momentum | |
167 | newParamCov(4,4) = 0.5*startingTrackParam->GetInverseBendingMomentum() * 0.5*startingTrackParam->GetInverseBendingMomentum(); | |
168 | startingTrackParam->SetCovariances(newParamCov); | |
169 | ||
170 | // Reset the track chi2 | |
171 | startingTrackParam->SetTrackChi2(0.); | |
172 | ||
173 | } | |
174 | ||
175 | // Redo the tracking | |
176 | RetracePartialTrack(trackCandidate, startingTrackParam); | |
177 | ||
178 | } | |
8cc77df0 | 179 | |
ea94c18b | 180 | //__________________________________________________________________________ |
181 | void AliMUONTrackReconstructorK::RetracePartialTrack(AliMUONTrack &trackCandidate, const AliMUONTrackParam* startingTrackParam) | |
182 | { | |
183 | /// Re-run the kalman filter from the hit attached to startingTrackParam to the most uptream hit | |
184 | AliDebug(1,"Enter RetracePartialTrack"); | |
185 | ||
186 | // Printout for debuging | |
187 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { | |
188 | cout << "RetracePartialTrack: track chi2 before re-tracking: " << trackCandidate.GetFitFMin() << endl; | |
189 | } | |
190 | ||
191 | // Reset the track chi2 | |
192 | trackCandidate.SetFitFMin(startingTrackParam->GetTrackChi2()); | |
193 | ||
194 | // loop over attached hits until the first one and recompute track parameters and covariances using kalman filter | |
195 | Int_t expectedChamber = startingTrackParam->GetHitForRecPtr()->GetChamberNumber() - 1; | |
196 | Int_t currentChamber; | |
197 | Double_t addChi2TrackAtHit; | |
198 | AliMUONTrackParam* trackParamAtHit = (AliMUONTrackParam*) trackCandidate.GetTrackParamAtHit()->Before(startingTrackParam); | |
199 | while (trackParamAtHit) { | |
200 | ||
201 | // reset track parameters and their covariances | |
202 | trackParamAtHit->SetParameters(startingTrackParam->GetParameters()); | |
203 | trackParamAtHit->SetZ(startingTrackParam->GetZ()); | |
204 | trackParamAtHit->SetCovariances(startingTrackParam->GetCovariances()); | |
205 | ||
206 | // add MCS effect | |
207 | AliMUONTrackExtrap::AddMCSEffect(trackParamAtHit,AliMUONConstants::ChamberThicknessInX0(),1.); | |
208 | ||
209 | // reset propagator for smoother | |
210 | if (fgkRunSmoother) trackParamAtHit->ResetPropagator(); | |
211 | ||
212 | // add MCS in missing chambers if any (at most 2 chambers can be missing according to tracking criteria) | |
213 | currentChamber = trackParamAtHit->GetHitForRecPtr()->GetChamberNumber(); | |
214 | while (currentChamber < expectedChamber) { | |
215 | // extrapolation to the missing chamber (update the propagator) | |
216 | AliMUONTrackExtrap::ExtrapToZCov(trackParamAtHit, AliMUONConstants::DefaultChamberZ(expectedChamber), fgkRunSmoother); | |
217 | // add MCS effect | |
218 | AliMUONTrackExtrap::AddMCSEffect(trackParamAtHit,AliMUONConstants::ChamberThicknessInX0(),1.); | |
219 | expectedChamber--; | |
220 | } | |
221 | ||
222 | // extrapolation to the plane of the hitForRec attached to the current trackParamAtHit (update the propagator) | |
223 | AliMUONTrackExtrap::ExtrapToZCov(trackParamAtHit, trackParamAtHit->GetHitForRecPtr()->GetZ(), fgkRunSmoother); | |
224 | ||
225 | if (fgkRunSmoother) { | |
226 | // save extrapolated parameters for smoother | |
227 | trackParamAtHit->SetExtrapParameters(trackParamAtHit->GetParameters()); | |
228 | ||
229 | // save extrapolated covariance matrix for smoother | |
230 | trackParamAtHit->SetExtrapCovariances(trackParamAtHit->GetCovariances()); | |
231 | } | |
232 | ||
233 | // Compute new track parameters including "hitForRecCh2" using kalman filter | |
234 | addChi2TrackAtHit = RunKalmanFilter(*trackParamAtHit); | |
235 | ||
236 | // Update the track chi2 | |
237 | trackCandidate.SetFitFMin(trackCandidate.GetFitFMin() + addChi2TrackAtHit); | |
238 | trackParamAtHit->SetTrackChi2(trackCandidate.GetFitFMin()); | |
239 | ||
240 | // prepare next step, add MCS effects in parameter covariances | |
241 | expectedChamber--; | |
242 | startingTrackParam = trackParamAtHit; | |
243 | trackParamAtHit = (AliMUONTrackParam*) (trackCandidate.GetTrackParamAtHit()->Before(startingTrackParam)); | |
244 | } | |
245 | ||
246 | // Printout for debuging | |
247 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { | |
248 | cout << "RetracePartialTrack: track chi2 after re-tracking: " << trackCandidate.GetFitFMin() << endl; | |
249 | } | |
250 | ||
251 | } | |
8cc77df0 | 252 | |
ea94c18b | 253 | //__________________________________________________________________________ |
254 | void AliMUONTrackReconstructorK::FollowTracks() | |
255 | { | |
256 | /// Follow tracks in stations(1..) 3, 2 and 1 | |
257 | AliDebug(1,"Enter FollowTracks"); | |
258 | ||
259 | AliMUONTrack *track; | |
260 | Int_t currentNRecTracks; | |
261 | Bool_t hitFound; | |
262 | ||
263 | for (Int_t station = 2; station >= 0; station--) { | |
264 | ||
265 | // Save the actual number of reconstructed track in case of | |
266 | // tracks are added or suppressed during the tracking procedure | |
267 | // !! Do not compress fRecTracksPtr until the end of the loop over tracks !! | |
268 | currentNRecTracks = fNRecTracks; | |
269 | ||
270 | for (Int_t iRecTrack = 0; iRecTrack <currentNRecTracks; iRecTrack++) { | |
271 | AliDebug(1,Form("FollowTracks: track candidate(1..) %d", iRecTrack+1)); | |
272 | ||
273 | track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iRecTrack); | |
274 | ||
275 | // Printout for debuging | |
276 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) { | |
277 | cout<<endl<<"Track parameter covariances at first hit:"<<endl; | |
278 | ((AliMUONTrackParam*) (track->GetTrackParamAtHit()->First()))->GetCovariances().Print(); | |
279 | } | |
280 | ||
281 | // Look for compatible hitForRec in station(0..) "station" | |
282 | hitFound = FollowTrackInStation(*track,station); | |
283 | ||
284 | // Try to recover track if required | |
285 | if (!hitFound && fgkRecoverTracks) hitFound = RecoverTrack(*track,station); | |
286 | ||
287 | // remove track if no hit found | |
288 | if (!hitFound) { | |
289 | fRecTracksPtr->Remove(track); | |
290 | fNRecTracks--; | |
8cc77df0 | 291 | } |
ea94c18b | 292 | |
8cc77df0 | 293 | } |
ea94c18b | 294 | |
295 | fRecTracksPtr->Compress(); // this is essential before checking tracks | |
296 | ||
297 | // Remove bad tracks | |
298 | Int_t currentNRecTracks = fNRecTracks; | |
299 | for (Int_t iRecTrack = 0; iRecTrack < currentNRecTracks; iRecTrack++) { | |
300 | track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iRecTrack); | |
301 | ||
302 | // Remove the track if the normalized chi2 is too high | |
303 | if (track->GetNormalizedChi2() > fgkSigmaToCutForTracking * fgkSigmaToCutForTracking) { | |
304 | fRecTracksPtr->Remove(track); | |
305 | fNRecTracks--; | |
306 | } | |
307 | ||
308 | } | |
309 | ||
310 | fRecTracksPtr->Compress(); // this is essential before checking tracks | |
311 | ||
312 | // Keep only the best tracks if required | |
313 | if (!fgkTrackAllTracks) RemoveDoubleTracks(); | |
314 | ||
315 | } | |
316 | ||
317 | } | |
b0e4c598 | 318 | |
ea94c18b | 319 | //__________________________________________________________________________ |
320 | Bool_t AliMUONTrackReconstructorK::FollowTrackInStation(AliMUONTrack &trackCandidate, Int_t nextStation) | |
321 | { | |
322 | /// Follow trackCandidate in station(0..) nextStation and search for compatible HitForRec(s) | |
323 | /// Keep all possibilities or only the best one(s) according to the flag fgkTrackAllTracks: | |
324 | /// kTRUE: duplicate "trackCandidate" if there are several possibilities and add the new tracks at the end of | |
325 | /// fRecTracksPtr to avoid conficts with other track candidates at this current stage of the tracking procedure. | |
326 | /// Remove the obsolete "trackCandidate" at the end. | |
327 | /// kFALSE: add only the best hit(s) to the "trackCandidate". Try to add a couple of hits in priority. | |
328 | /// return kTRUE if new hits have been found (otherwise return kFALSE) | |
329 | AliDebug(1,Form("Enter FollowTrackInStation(1..) %d", nextStation+1)); | |
330 | ||
331 | // Order the chamber according to the propagation direction (tracking starts with chamber 2): | |
332 | // - nextStation == station(1...) 5 => forward propagation | |
333 | // - nextStation < station(1...) 5 => backward propagation | |
334 | Int_t ch1, ch2; | |
335 | if (nextStation==4) { | |
336 | ch1 = 2*nextStation+1; | |
337 | ch2 = 2*nextStation; | |
338 | } else { | |
339 | ch1 = 2*nextStation; | |
340 | ch2 = 2*nextStation+1; | |
341 | } | |
342 | ||
343 | Double_t zCh2 = AliMUONConstants::DefaultChamberZ(ch2); | |
344 | Double_t chi2OfHitForRec; | |
345 | Double_t maxChi2OfHitForRec = 2. * fgkSigmaToCutForTracking * fgkSigmaToCutForTracking; // 2 because 2 quantities in chi2 | |
346 | Double_t addChi2TrackAtHit1; | |
347 | Double_t addChi2TrackAtHit2; | |
348 | Double_t bestAddChi2TrackAtHit1 = 1.e10; | |
349 | Double_t bestAddChi2TrackAtHit2 = 1.e10; | |
350 | Bool_t foundOneHit = kFALSE; | |
351 | Bool_t foundTwoHits = kFALSE; | |
352 | AliMUONTrack *newTrack = 0x0; | |
353 | AliMUONHitForRec *hitForRecCh1, *hitForRecCh2; | |
354 | AliMUONTrackParam extrapTrackParam; | |
355 | AliMUONTrackParam extrapTrackParamAtHit1; | |
356 | AliMUONTrackParam extrapTrackParamAtHit2; | |
357 | AliMUONTrackParam bestTrackParamAtHit1; | |
358 | AliMUONTrackParam bestTrackParamAtHit2; | |
359 | Bool_t *hitForRecCh1Used = new Bool_t[fNHitsForRecPerChamber[ch1]]; | |
360 | for (Int_t hit1 = 0; hit1 < fNHitsForRecPerChamber[ch1]; hit1++) hitForRecCh1Used[hit1] = kFALSE; | |
361 | ||
362 | // Get track parameters and extrapolate them to chamber "ch2" to save computing time in the next steps | |
363 | AliMUONTrackParam extrapTrackParamAtCh2(*(AliMUONTrackParam*)trackCandidate.GetTrackParamAtHit()->First()); | |
364 | ||
365 | // Add MCS effect | |
366 | AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh2,AliMUONConstants::ChamberThicknessInX0(),1.); | |
367 | ||
368 | // reset propagator for smoother | |
369 | if (fgkRunSmoother) extrapTrackParamAtCh2.ResetPropagator(); | |
370 | ||
371 | // Add MCS in the missing chamber if any (only 1 chamber can be missing according to tracking criteria) | |
372 | if (ch1 < ch2 && extrapTrackParamAtCh2.GetHitForRecPtr()->GetChamberNumber() > ch2 + 1) { | |
373 | // extrapolation to the missing chamber | |
374 | AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh2, AliMUONConstants::DefaultChamberZ(ch2 + 1), fgkRunSmoother); | |
375 | // add MCS effect | |
376 | AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh2,AliMUONConstants::ChamberThicknessInX0(),1.); | |
377 | } | |
378 | ||
379 | //Extrapolate trackCandidate to chamber "ch2" | |
380 | AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh2, zCh2, fgkRunSmoother); | |
381 | ||
382 | // Printout for debuging | |
383 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) { | |
384 | cout<<endl<<"Track parameter covariances at first hit extrapolated to z = "<<zCh2<<":"<<endl; | |
385 | extrapTrackParamAtCh2.GetCovariances().Print(); | |
386 | } | |
387 | ||
388 | // Printout for debuging | |
389 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { | |
390 | cout << "FollowTrackInStation: look for hits in chamber(1..): " << ch2+1 << endl; | |
391 | } | |
392 | ||
393 | // look for candidates in chamber 2 | |
394 | for (Int_t hit2 = 0; hit2 < fNHitsForRecPerChamber[ch2]; hit2++) { | |
395 | ||
396 | hitForRecCh2 = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[ch2]+hit2); | |
397 | ||
398 | // try to add the current hit | |
399 | chi2OfHitForRec = TryOneHitForRec(extrapTrackParamAtCh2, hitForRecCh2, extrapTrackParamAtHit2, fgkRunSmoother); | |
400 | ||
401 | // if good chi2 then try to attach a hitForRec in the other chamber too | |
402 | if (chi2OfHitForRec < maxChi2OfHitForRec) { | |
403 | Bool_t foundSecondHit = kFALSE; | |
404 | ||
405 | // Printout for debuging | |
406 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { | |
407 | cout << "FollowTrackInStation: found one hit in chamber(1..): " << ch2+1 | |
408 | << " (Chi2 = " << chi2OfHitForRec << ")" << endl; | |
409 | cout << " look for second hits in chamber(1..): " << ch1+1 << " ..." << endl; | |
8cc77df0 | 410 | } |
ea94c18b | 411 | |
412 | if (fgkRunSmoother) { | |
413 | // save extrapolated parameters for smoother | |
414 | extrapTrackParamAtHit2.SetExtrapParameters(extrapTrackParamAtHit2.GetParameters()); | |
415 | ||
416 | // save extrapolated covariance matrix for smoother | |
417 | extrapTrackParamAtHit2.SetExtrapCovariances(extrapTrackParamAtHit2.GetCovariances()); | |
418 | } | |
419 | ||
420 | // Compute new track parameters including "hitForRecCh2" using kalman filter | |
421 | addChi2TrackAtHit2 = RunKalmanFilter(extrapTrackParamAtHit2); | |
422 | ||
423 | // copy new track parameters for next step | |
424 | extrapTrackParam = extrapTrackParamAtHit2; | |
425 | ||
426 | // add MCS effect | |
427 | AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParam,AliMUONConstants::ChamberThicknessInX0(),1.); | |
428 | ||
429 | // reset propagator for smoother | |
430 | if (fgkRunSmoother) extrapTrackParam.ResetPropagator(); | |
431 | ||
432 | for (Int_t hit1 = 0; hit1 < fNHitsForRecPerChamber[ch1]; hit1++) { | |
433 | ||
434 | hitForRecCh1 = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[ch1]+hit1); | |
435 | ||
436 | // try to add the current hit | |
437 | chi2OfHitForRec = TryOneHitForRec(extrapTrackParam, hitForRecCh1, extrapTrackParamAtHit1, fgkRunSmoother); | |
438 | ||
439 | // if good chi2 then consider to add the 2 hitForRec to the "trackCandidate" | |
440 | if (chi2OfHitForRec < maxChi2OfHitForRec) { | |
441 | foundSecondHit = kTRUE; | |
442 | foundTwoHits = kTRUE; | |
443 | ||
444 | // Printout for debuging | |
445 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { | |
446 | cout << "FollowTrackInStation: found one hit in chamber(1..): " << ch1+1 | |
447 | << " (Chi2 = " << chi2OfHitForRec << ")" << endl; | |
448 | } | |
449 | ||
450 | if (fgkRunSmoother) { | |
451 | // save extrapolated parameters for smoother | |
452 | extrapTrackParamAtHit1.SetExtrapParameters(extrapTrackParamAtHit1.GetParameters()); | |
453 | ||
454 | // save extrapolated covariance matrix for smoother | |
455 | extrapTrackParamAtHit1.SetExtrapCovariances(extrapTrackParamAtHit1.GetCovariances()); | |
456 | } | |
457 | ||
458 | // Compute new track parameters including "hitForRecCh1" using kalman filter | |
459 | addChi2TrackAtHit1 = RunKalmanFilter(extrapTrackParamAtHit1); | |
460 | ||
461 | if (fgkTrackAllTracks) { | |
462 | // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new hitForRec's | |
463 | newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate); | |
464 | UpdateTrack(*newTrack,extrapTrackParamAtHit1,extrapTrackParamAtHit2,addChi2TrackAtHit1,addChi2TrackAtHit2); | |
465 | fNRecTracks++; | |
466 | ||
467 | // if we are arrived on station(1..) 5, recompute track parameters and covariances starting from this station | |
468 | // (going in the right direction) | |
469 | if (nextStation == 4) RetraceTrack(*newTrack,kTRUE); | |
470 | ||
471 | // Tag hitForRecCh1 as used | |
472 | hitForRecCh1Used[hit1] = kTRUE; | |
473 | ||
474 | // Printout for debuging | |
475 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { | |
476 | cout << "FollowTrackInStation: added two hits in station(1..): " << nextStation+1 << endl; | |
477 | if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump(); | |
478 | } | |
479 | ||
480 | } else if (addChi2TrackAtHit1+addChi2TrackAtHit2 < bestAddChi2TrackAtHit1+bestAddChi2TrackAtHit2) { | |
481 | // keep track of the best couple of hits | |
482 | bestAddChi2TrackAtHit1 = addChi2TrackAtHit1; | |
483 | bestAddChi2TrackAtHit2 = addChi2TrackAtHit2; | |
484 | bestTrackParamAtHit1 = extrapTrackParamAtHit1; | |
485 | bestTrackParamAtHit2 = extrapTrackParamAtHit2; | |
486 | } | |
487 | ||
b0e4c598 | 488 | } |
ea94c18b | 489 | |
8cc77df0 | 490 | } |
ea94c18b | 491 | |
492 | // if no hitForRecCh1 found then consider to add hitForRecCh2 only | |
493 | if (!foundSecondHit) { | |
494 | foundOneHit = kTRUE; | |
495 | ||
496 | if (fgkTrackAllTracks) { | |
497 | // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new hitForRec's | |
498 | newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate); | |
499 | UpdateTrack(*newTrack,extrapTrackParamAtHit2,addChi2TrackAtHit2); | |
500 | fNRecTracks++; | |
501 | ||
502 | // if we are arrived on station(1..) 5, recompute track parameters and covariances starting from this station | |
503 | // (going in the right direction) | |
504 | if (nextStation == 4) RetraceTrack(*newTrack,kTRUE); | |
505 | ||
506 | // Printout for debuging | |
507 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { | |
508 | cout << "FollowTrackInStation: added one hit in chamber(1..): " << ch2+1 << endl; | |
509 | if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump(); | |
510 | } | |
511 | ||
512 | } else if (!foundTwoHits && addChi2TrackAtHit2 < bestAddChi2TrackAtHit1) { | |
513 | // keep track of the best single hitForRec except if a couple | |
514 | // of hits has already been found (i.e. bestHitForRec2!=0x0) | |
515 | bestAddChi2TrackAtHit1 = addChi2TrackAtHit2; | |
516 | bestTrackParamAtHit1 = extrapTrackParamAtHit2; | |
517 | } | |
518 | ||
519 | } | |
520 | ||
8cc77df0 | 521 | } |
ea94c18b | 522 | |
523 | } | |
524 | ||
525 | // look for candidates in chamber 1 not already attached to a track | |
526 | // if we want to keep all possible tracks or if no good couple of hitForRec has been found | |
527 | if (fgkTrackAllTracks || !foundTwoHits) { | |
528 | ||
529 | // Printout for debuging | |
530 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { | |
531 | cout << "FollowTrackInStation: look for single hits in chamber(1..): " << ch1+1 << endl; | |
8cc77df0 | 532 | } |
ea94c18b | 533 | |
534 | // add MCS effect for next step | |
535 | AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh2,AliMUONConstants::ChamberThicknessInX0(),1.); | |
536 | ||
537 | for (Int_t hit1 = 0; hit1 < fNHitsForRecPerChamber[ch1]; hit1++) { | |
538 | ||
539 | hitForRecCh1 = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[ch1]+hit1); | |
540 | ||
541 | if (hitForRecCh1Used[hit1]) continue; // Skip hitForRec already used | |
542 | ||
543 | // try to add the current hit | |
544 | chi2OfHitForRec = TryOneHitForRec(extrapTrackParamAtCh2, hitForRecCh1, extrapTrackParamAtHit1, fgkRunSmoother); | |
545 | ||
546 | // if good chi2 then consider to add hitForRecCh1 | |
547 | // We do not try to attach a hitForRec in the other chamber too since it has already been done above | |
548 | if (chi2OfHitForRec < maxChi2OfHitForRec) { | |
549 | foundOneHit = kTRUE; | |
550 | ||
551 | // Printout for debuging | |
552 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { | |
553 | cout << "FollowTrackInStation: found one hit in chamber(1..): " << ch1+1 | |
554 | << " (Chi2 = " << chi2OfHitForRec << ")" << endl; | |
555 | } | |
556 | ||
557 | if (fgkRunSmoother) { | |
558 | // save extrapolated parameters for smoother | |
559 | extrapTrackParamAtHit1.SetExtrapParameters(extrapTrackParamAtHit1.GetParameters()); | |
560 | ||
561 | // save extrapolated covariance matrix for smoother | |
562 | extrapTrackParamAtHit1.SetExtrapCovariances(extrapTrackParamAtHit1.GetCovariances()); | |
563 | } | |
564 | ||
565 | // Compute new track parameters including "hitForRecCh1" using kalman filter | |
566 | addChi2TrackAtHit1 = RunKalmanFilter(extrapTrackParamAtHit1); | |
567 | ||
568 | if (fgkTrackAllTracks) { | |
569 | // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new hitForRec's | |
570 | newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate); | |
571 | UpdateTrack(*newTrack,extrapTrackParamAtHit1,addChi2TrackAtHit1); | |
572 | fNRecTracks++; | |
573 | ||
574 | // if we are arrived on station(1..) 5, recompute track parameters and covariances starting from this station | |
575 | // (going in the right direction) | |
576 | if (nextStation == 4) RetraceTrack(*newTrack,kTRUE); | |
577 | ||
578 | // Printout for debuging | |
579 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { | |
580 | cout << "FollowTrackInStation: added one hit in chamber(1..): " << ch1+1 << endl; | |
581 | if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump(); | |
582 | } | |
583 | ||
584 | } else if (addChi2TrackAtHit1 < bestAddChi2TrackAtHit1) { | |
585 | // keep track of the best single hitForRec except if a couple | |
586 | // of hits has already been found (i.e. bestHitForRec1!=0x0) | |
587 | bestAddChi2TrackAtHit1 = addChi2TrackAtHit1; | |
588 | bestTrackParamAtHit1 = extrapTrackParamAtHit1; | |
589 | } | |
590 | ||
591 | } | |
592 | ||
8cc77df0 | 593 | } |
ea94c18b | 594 | |
595 | } | |
596 | ||
597 | // fill out the best track if required else clean up the fRecTracksPtr array | |
598 | if (!fgkTrackAllTracks) { | |
599 | if (foundTwoHits) { | |
600 | UpdateTrack(trackCandidate,bestTrackParamAtHit1,bestTrackParamAtHit2,bestAddChi2TrackAtHit1,bestAddChi2TrackAtHit2); | |
601 | ||
602 | // if we are arrived on station(1..) 5, recompute track parameters and covariances starting from this station | |
603 | // (going in the right direction) | |
604 | if (nextStation == 4) RetraceTrack(trackCandidate,kTRUE); | |
8cc77df0 | 605 | |
ea94c18b | 606 | // Printout for debuging |
607 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { | |
608 | cout << "FollowTrackInStation: added the two best hits in station(1..): " << nextStation+1 << endl; | |
609 | if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump(); | |
610 | } | |
611 | ||
612 | } else if (foundOneHit) { | |
613 | UpdateTrack(trackCandidate,bestTrackParamAtHit1,bestAddChi2TrackAtHit1); | |
614 | ||
615 | // if we are arrived on station(1..) 5, recompute track parameters and covariances starting from this station | |
616 | // (going in the right direction) | |
617 | if (nextStation == 4) RetraceTrack(trackCandidate,kTRUE); | |
618 | ||
619 | // Printout for debuging | |
620 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { | |
621 | cout << "FollowTrackInStation: added the best hit in chamber(1..): " << bestTrackParamAtHit1.GetHitForRecPtr()->GetChamberNumber()+1 << endl; | |
622 | if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump(); | |
623 | } | |
624 | ||
625 | } else return kFALSE; | |
626 | ||
627 | } else if (foundOneHit || foundTwoHits) { | |
628 | ||
629 | // remove obsolete track | |
630 | fRecTracksPtr->Remove(&trackCandidate); | |
631 | fNRecTracks--; | |
632 | ||
633 | } else return kFALSE; | |
634 | ||
635 | return kTRUE; | |
636 | ||
637 | } | |
8cc77df0 | 638 | |
ea94c18b | 639 | //__________________________________________________________________________ |
640 | Double_t AliMUONTrackReconstructorK::RunKalmanFilter(AliMUONTrackParam &trackParamAtHit) | |
641 | { | |
642 | /// Compute new track parameters and their covariances including new hit using kalman filter | |
643 | /// return the additional track chi2 | |
644 | AliDebug(1,"Enter RunKalmanFilter"); | |
645 | ||
646 | // Get actual track parameters (p) | |
647 | TMatrixD param(trackParamAtHit.GetParameters()); | |
648 | ||
649 | // Get new hit parameters (m) | |
650 | AliMUONHitForRec *hitForRecAtHit = trackParamAtHit.GetHitForRecPtr(); | |
651 | TMatrixD hit(5,1); | |
652 | hit.Zero(); | |
653 | hit(0,0) = hitForRecAtHit->GetNonBendingCoor(); | |
654 | hit(2,0) = hitForRecAtHit->GetBendingCoor(); | |
655 | ||
656 | // Compute the actual parameter weight (W) | |
657 | TMatrixD paramWeight(trackParamAtHit.GetCovariances()); | |
658 | if (paramWeight.Determinant() != 0) { | |
659 | paramWeight.Invert(); | |
660 | } else { | |
661 | AliWarning(" Determinant = 0"); | |
662 | return 1.e10; | |
8cc77df0 | 663 | } |
ea94c18b | 664 | |
665 | // Compute the new hit weight (U) | |
666 | TMatrixD hitWeight(5,5); | |
667 | hitWeight.Zero(); | |
668 | hitWeight(0,0) = 1. / hitForRecAtHit->GetNonBendingReso2(); | |
669 | hitWeight(2,2) = 1. / hitForRecAtHit->GetBendingReso2(); | |
8cc77df0 | 670 | |
ea94c18b | 671 | // Compute the new parameters covariance matrix ( (W+U)^-1 ) |
672 | TMatrixD newParamCov(paramWeight,TMatrixD::kPlus,hitWeight); | |
673 | if (newParamCov.Determinant() != 0) { | |
674 | newParamCov.Invert(); | |
675 | } else { | |
676 | AliWarning(" Determinant = 0"); | |
677 | return 1.e10; | |
678 | } | |
679 | ||
680 | // Save the new parameters covariance matrix | |
681 | trackParamAtHit.SetCovariances(newParamCov); | |
682 | ||
683 | // Compute the new parameters (p' = ((W+U)^-1)U(m-p) + p) | |
684 | TMatrixD tmp(hit,TMatrixD::kMinus,param); | |
685 | TMatrixD tmp2(hitWeight,TMatrixD::kMult,tmp); // U(m-p) | |
686 | TMatrixD newParam(newParamCov,TMatrixD::kMult,tmp2); // ((W+U)^-1)U(m-p) | |
687 | newParam += param; // ((W+U)^-1)U(m-p) + p | |
688 | ||
689 | // Save the new parameters | |
690 | trackParamAtHit.SetParameters(newParam); | |
691 | ||
692 | // Compute the additional chi2 (= ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m)) | |
693 | tmp = newParam; // p' | |
694 | tmp -= param; // (p'-p) | |
695 | TMatrixD tmp3(paramWeight,TMatrixD::kMult,tmp); // W(p'-p) | |
696 | TMatrixD addChi2Track(tmp,TMatrixD::kTransposeMult,tmp3); // ((p'-p)^-1)W(p'-p) | |
697 | tmp = newParam; // p' | |
698 | tmp -= hit; // (p'-m) | |
699 | TMatrixD tmp4(hitWeight,TMatrixD::kMult,tmp); // U(p'-m) | |
700 | addChi2Track += TMatrixD(tmp,TMatrixD::kTransposeMult,tmp4); // ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m) | |
701 | ||
702 | return addChi2Track(0,0); | |
703 | ||
8cc77df0 | 704 | } |
705 | ||
ea94c18b | 706 | //__________________________________________________________________________ |
707 | void AliMUONTrackReconstructorK::UpdateTrack(AliMUONTrack &track, AliMUONTrackParam &trackParamAtHit, Double_t addChi2) | |
8cc77df0 | 708 | { |
ea94c18b | 709 | /// Add 1 hit to the track candidate |
710 | /// Update chi2 of the track | |
711 | ||
712 | // Update the track chi2 into TrackParamAtHit | |
713 | trackParamAtHit.SetTrackChi2(track.GetFitFMin() + addChi2); | |
714 | ||
715 | // Update the chi2 of the new track | |
716 | track.SetFitFMin(trackParamAtHit.GetTrackChi2()); | |
717 | ||
718 | // Update array of TrackParamAtHit | |
719 | track.AddTrackParamAtHit(&trackParamAtHit,trackParamAtHit.GetHitForRecPtr()); | |
720 | track.GetTrackParamAtHit()->Sort(); | |
721 | ||
722 | } | |
8cc77df0 | 723 | |
ea94c18b | 724 | //__________________________________________________________________________ |
725 | void AliMUONTrackReconstructorK::UpdateTrack(AliMUONTrack &track, AliMUONTrackParam &trackParamAtHit1, AliMUONTrackParam &trackParamAtHit2, | |
726 | Double_t addChi2AtHit1, Double_t addChi2AtHit2) | |
727 | { | |
728 | /// Add 2 hits to the track candidate | |
729 | /// Update track and local chi2 | |
730 | ||
731 | // Update local chi2 at first hit | |
732 | AliMUONHitForRec* hit1 = trackParamAtHit1.GetHitForRecPtr(); | |
733 | Double_t deltaX = trackParamAtHit1.GetNonBendingCoor() - hit1->GetNonBendingCoor(); | |
734 | Double_t deltaY = trackParamAtHit1.GetBendingCoor() - hit1->GetBendingCoor(); | |
735 | Double_t localChi2AtHit1 = deltaX*deltaX / hit1->GetNonBendingReso2() + | |
736 | deltaY*deltaY / hit1->GetBendingReso2(); | |
737 | trackParamAtHit1.SetLocalChi2(localChi2AtHit1); | |
738 | ||
739 | // Flag first hit as being removable | |
740 | trackParamAtHit1.SetRemovable(kTRUE); | |
741 | ||
742 | // Update local chi2 at second hit | |
743 | AliMUONHitForRec* hit2 = trackParamAtHit2.GetHitForRecPtr(); | |
744 | AliMUONTrackParam extrapTrackParamAtHit2(trackParamAtHit1); | |
745 | AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParamAtHit2, trackParamAtHit2.GetZ()); | |
746 | deltaX = extrapTrackParamAtHit2.GetNonBendingCoor() - hit2->GetNonBendingCoor(); | |
747 | deltaY = extrapTrackParamAtHit2.GetBendingCoor() - hit2->GetBendingCoor(); | |
748 | Double_t localChi2AtHit2 = deltaX*deltaX / hit2->GetNonBendingReso2() + | |
749 | deltaY*deltaY / hit2->GetBendingReso2(); | |
750 | trackParamAtHit2.SetLocalChi2(localChi2AtHit2); | |
751 | ||
752 | // Flag second hit as being removable | |
753 | trackParamAtHit2.SetRemovable(kTRUE); | |
754 | ||
755 | // Update the track chi2 into TrackParamAtHit1 | |
756 | trackParamAtHit1.SetTrackChi2(track.GetFitFMin() + addChi2AtHit1); | |
757 | ||
758 | // Update the track chi2 into TrackParamAtHit2 | |
759 | trackParamAtHit2.SetTrackChi2(trackParamAtHit1.GetTrackChi2() + addChi2AtHit2); | |
760 | ||
761 | // Update the chi2 of the new track | |
762 | track.SetFitFMin(trackParamAtHit2.GetTrackChi2()); | |
763 | ||
764 | // Update array of TrackParamAtHit | |
765 | track.AddTrackParamAtHit(&trackParamAtHit1,trackParamAtHit1.GetHitForRecPtr()); | |
766 | track.AddTrackParamAtHit(&trackParamAtHit2,trackParamAtHit2.GetHitForRecPtr()); | |
767 | track.GetTrackParamAtHit()->Sort(); | |
768 | ||
769 | } | |
8cc77df0 | 770 | |
ea94c18b | 771 | //__________________________________________________________________________ |
772 | Bool_t AliMUONTrackReconstructorK::RecoverTrack(AliMUONTrack &trackCandidate, Int_t nextStation) | |
773 | { | |
774 | /// Try to recover the track candidate in the next station | |
775 | /// by removing the worst of the two hits attached in the current station | |
776 | /// Return kTRUE if recovering succeeds | |
777 | AliDebug(1,"Enter RecoverTrack"); | |
778 | ||
779 | // Do not try to recover track until we have attached hit(s) on station(1..) 3 | |
780 | if (nextStation > 1) return kFALSE; | |
781 | ||
782 | Int_t worstHitNumber = -1; | |
783 | Double_t localChi2, worstChi2 = 0.; | |
784 | ||
785 | // Look for the hit to remove | |
786 | for (Int_t hitNumber = 0; hitNumber < 2; hitNumber++) { | |
787 | AliMUONTrackParam *trackParamAtHit = (AliMUONTrackParam*)trackCandidate.GetTrackParamAtHit()->UncheckedAt(hitNumber); | |
788 | ||
789 | // check if current hit is removable | |
790 | if (!trackParamAtHit->IsRemovable()) return kFALSE; | |
791 | ||
792 | // Pick up hit with the worst chi2 | |
793 | localChi2 = trackParamAtHit->GetLocalChi2(); | |
794 | if (localChi2 > worstChi2) { | |
795 | worstChi2 = localChi2; | |
796 | worstHitNumber = hitNumber; | |
8cc77df0 | 797 | } |
ea94c18b | 798 | } |
799 | ||
800 | // Reset best hit as being NOT removable | |
801 | ((AliMUONTrackParam*)trackCandidate.GetTrackParamAtHit()->UncheckedAt((worstHitNumber+1)%2))->SetRemovable(kFALSE); | |
802 | ||
803 | // Remove the worst hit | |
804 | trackCandidate.RemoveTrackParamAtHit((AliMUONTrackParam*)trackCandidate.GetTrackParamAtHit()->UncheckedAt(worstHitNumber)); | |
805 | ||
806 | // Re-calculate track parameters at the (new) first hit | |
807 | RetracePartialTrack(trackCandidate,(AliMUONTrackParam*)trackCandidate.GetTrackParamAtHit()->UncheckedAt(1)); | |
808 | ||
809 | // Look for new hit(s) in next station | |
810 | return FollowTrackInStation(trackCandidate,nextStation); | |
811 | ||
8cc77df0 | 812 | } |
813 | ||
814 | //__________________________________________________________________________ | |
ea94c18b | 815 | Bool_t AliMUONTrackReconstructorK::RunSmoother(AliMUONTrack &track) |
8cc77df0 | 816 | { |
ea94c18b | 817 | /// Compute new track parameters and their covariances using smoother |
818 | AliDebug(1,"Enter RunSmoother"); | |
819 | ||
820 | AliMUONTrackParam *previousTrackParam = (AliMUONTrackParam*) track.GetTrackParamAtHit()->First(); | |
821 | ||
822 | // Smoothed parameters and covariances at first hit = filtered parameters and covariances | |
823 | previousTrackParam->SetSmoothParameters(previousTrackParam->GetParameters()); | |
824 | previousTrackParam->SetSmoothCovariances(previousTrackParam->GetCovariances()); | |
825 | ||
826 | // Compute local chi2 at first hit | |
827 | AliMUONHitForRec *hitForRecAtHit = previousTrackParam->GetHitForRecPtr(); | |
828 | Double_t dX = hitForRecAtHit->GetNonBendingCoor() - previousTrackParam->GetNonBendingCoor(); | |
829 | Double_t dY = hitForRecAtHit->GetBendingCoor() - previousTrackParam->GetBendingCoor(); | |
830 | Double_t localChi2 = dX * dX / hitForRecAtHit->GetNonBendingReso2() + dY * dY / hitForRecAtHit->GetBendingReso2(); | |
831 | ||
832 | // Save local chi2 at first hit | |
833 | previousTrackParam->SetLocalChi2(localChi2); | |
834 | ||
835 | AliMUONTrackParam *currentTrackParam = (AliMUONTrackParam*) track.GetTrackParamAtHit()->After(previousTrackParam); | |
836 | while (currentTrackParam) { | |
837 | ||
838 | // Get variables | |
839 | const TMatrixD &extrapParameters = previousTrackParam->GetExtrapParameters(); // X(k+1 k) | |
840 | const TMatrixD &filteredParameters = currentTrackParam->GetParameters(); // X(k k) | |
841 | const TMatrixD &previousSmoothParameters = previousTrackParam->GetSmoothParameters(); // X(k+1 n) | |
842 | const TMatrixD &propagator = previousTrackParam->GetPropagator(); // F(k) | |
843 | const TMatrixD &extrapCovariances = previousTrackParam->GetExtrapCovariances(); // C(k+1 k) | |
844 | const TMatrixD &filteredCovariances = currentTrackParam->GetCovariances(); // C(k k) | |
845 | const TMatrixD &previousSmoothCovariances = previousTrackParam->GetSmoothCovariances(); // C(k+1 n) | |
846 | ||
847 | // Compute smoother gain: A(k) = C(kk) * F(f)^t * (C(k+1 k))^-1 | |
848 | TMatrixD extrapWeight(extrapCovariances); | |
849 | if (extrapWeight.Determinant() != 0) { | |
850 | extrapWeight.Invert(); // (C(k+1 k))^-1 | |
851 | } else { | |
852 | AliWarning(" Determinant = 0"); | |
853 | return kFALSE; | |
854 | } | |
855 | TMatrixD smootherGain(filteredCovariances,TMatrixD::kMultTranspose,propagator); // C(kk) * F(f)^t | |
856 | smootherGain *= extrapWeight; // C(kk) * F(f)^t * (C(k+1 k))^-1 | |
857 | ||
858 | // Compute smoothed parameters: X(k n) = X(k k) + A(k) * (X(k+1 n) - X(k+1 k)) | |
859 | TMatrixD tmpParam(previousSmoothParameters,TMatrixD::kMinus,extrapParameters); // X(k+1 n) - X(k+1 k) | |
860 | TMatrixD smoothParameters(smootherGain,TMatrixD::kMult,tmpParam); // A(k) * (X(k+1 n) - X(k+1 k)) | |
861 | smoothParameters += filteredParameters; // X(k k) + A(k) * (X(k+1 n) - X(k+1 k)) | |
862 | ||
863 | // Save smoothed parameters | |
864 | currentTrackParam->SetSmoothParameters(smoothParameters); | |
865 | ||
866 | // Compute smoothed covariances: C(k n) = C(k k) + A(k) * (C(k+1 n) - C(k+1 k)) * (A(k))^t | |
867 | TMatrixD tmpCov(previousSmoothCovariances,TMatrixD::kMinus,extrapCovariances); // C(k+1 n) - C(k+1 k) | |
868 | TMatrixD tmpCov2(tmpCov,TMatrixD::kMultTranspose,smootherGain); // (C(k+1 n) - C(k+1 k)) * (A(k))^t | |
869 | TMatrixD smoothCovariances(smootherGain,TMatrixD::kMult,tmpCov2); // A(k) * (C(k+1 n) - C(k+1 k)) * (A(k))^t | |
870 | smoothCovariances += filteredCovariances; // C(k k) + A(k) * (C(k+1 n) - C(k+1 k)) * (A(k))^t | |
871 | ||
872 | // Save smoothed covariances | |
873 | currentTrackParam->SetSmoothCovariances(smoothCovariances); | |
874 | ||
875 | // Compute smoothed residual: r(k n) = hit - X(k n) | |
876 | hitForRecAtHit = currentTrackParam->GetHitForRecPtr(); | |
877 | TMatrixD smoothResidual(2,1); | |
878 | smoothResidual.Zero(); | |
879 | smoothResidual(0,0) = hitForRecAtHit->GetNonBendingCoor() - smoothParameters(0,0); | |
880 | smoothResidual(1,0) = hitForRecAtHit->GetBendingCoor() - smoothParameters(2,0); | |
881 | ||
882 | // Compute weight of smoothed residual: W(k n) = (hitCov - C(k n))^-1 | |
883 | TMatrixD smoothResidualWeight(2,2); | |
884 | smoothResidualWeight(0,0) = hitForRecAtHit->GetNonBendingReso2() - smoothCovariances(0,0); | |
885 | smoothResidualWeight(0,1) = - smoothCovariances(0,2); | |
886 | smoothResidualWeight(1,0) = - smoothCovariances(2,0); | |
887 | smoothResidualWeight(1,1) = hitForRecAtHit->GetBendingReso2() - smoothCovariances(2,2); | |
888 | if (smoothResidualWeight.Determinant() != 0) { | |
889 | smoothResidualWeight.Invert(); | |
890 | } else { | |
891 | AliWarning(" Determinant = 0"); | |
892 | return kFALSE; | |
893 | } | |
894 | ||
895 | // Compute local chi2 = (r(k n))^t * W(k n) * r(k n) | |
896 | TMatrixD tmpChi2(smoothResidual,TMatrixD::kTransposeMult,smoothResidualWeight); // (r(k n))^t * W(k n) | |
897 | TMatrixD localChi2(tmpChi2,TMatrixD::kMult,smoothResidual); // (r(k n))^t * W(k n) * r(k n) | |
898 | ||
899 | // Save local chi2 | |
900 | currentTrackParam->SetLocalChi2(localChi2(0,0)); | |
901 | ||
902 | previousTrackParam = currentTrackParam; | |
903 | currentTrackParam = (AliMUONTrackParam*) track.GetTrackParamAtHit()->After(previousTrackParam); | |
904 | } | |
905 | ||
906 | return kTRUE; | |
907 | ||
8cc77df0 | 908 | } |
909 | ||
8cc77df0 | 910 | //__________________________________________________________________________ |
ea94c18b | 911 | void AliMUONTrackReconstructorK::ImproveTracks() |
8cc77df0 | 912 | { |
ea94c18b | 913 | /// Improve tracks by removing clusters with local chi2 highter than the defined cut |
914 | /// Recompute track parameters and covariances at the remaining clusters | |
915 | AliDebug(1,"Enter ImproveTracks"); | |
916 | ||
917 | Double_t localChi2, worstLocalChi2; | |
918 | Int_t worstChamber; | |
919 | AliMUONTrackParam *trackParamAtHit, *worstTrackParamAtHit; | |
920 | Bool_t smoothed; | |
921 | ||
922 | // Remove double track to improve only "good" tracks | |
923 | RemoveDoubleTracks(); | |
924 | ||
925 | AliMUONTrack *track = (AliMUONTrack*) fRecTracksPtr->First(); | |
8cc77df0 | 926 | while (track) { |
ea94c18b | 927 | |
928 | while (!track->IsImproved()) { | |
929 | ||
930 | // Run smoother if required | |
931 | smoothed = kFALSE; | |
932 | if (fgkRunSmoother) smoothed = RunSmoother(*track); | |
933 | ||
934 | // Use standard procedure to compute local chi2 if smoother not required or not working | |
935 | if (!smoothed) { | |
936 | ||
937 | // Update track parameters and covariances | |
938 | track->UpdateCovTrackParamAtHit(); | |
939 | ||
940 | // Compute local chi2 of each hits | |
941 | track->ComputeLocalChi2(kTRUE); | |
942 | } | |
943 | ||
944 | // Look for the hit to remove | |
945 | worstTrackParamAtHit = NULL; | |
946 | worstLocalChi2 = 0.; | |
947 | trackParamAtHit = (AliMUONTrackParam*)track->GetTrackParamAtHit()->First(); | |
948 | while (trackParamAtHit) { | |
949 | ||
950 | // Pick up hit with the worst chi2 | |
951 | localChi2 = trackParamAtHit->GetLocalChi2(); | |
952 | if (localChi2 > worstLocalChi2) { | |
953 | worstLocalChi2 = localChi2; | |
954 | worstTrackParamAtHit = trackParamAtHit; | |
955 | } | |
956 | ||
957 | trackParamAtHit = (AliMUONTrackParam*)track->GetTrackParamAtHit()->After(trackParamAtHit); | |
958 | } | |
959 | ||
960 | // Check if bad removable hit found | |
961 | if (!worstTrackParamAtHit) { | |
962 | track->SetImproved(kTRUE); | |
963 | break; | |
964 | } | |
965 | ||
966 | // check whether the worst hit is removable or not | |
967 | if (!worstTrackParamAtHit->IsRemovable()) { | |
968 | track->SetImproved(kTRUE); | |
969 | break; | |
970 | } | |
971 | ||
972 | // Check whether the worst chi2 is under requirement or not | |
973 | if (worstLocalChi2 < 2. * fgkSigmaToCutForImprovement * fgkSigmaToCutForImprovement) { // 2 because 2 quantities in chi2 | |
974 | track->SetImproved(kTRUE); | |
975 | break; | |
976 | } | |
977 | ||
978 | // Reset the second hit in the same station as the bad one as being NOT removable | |
979 | worstChamber = worstTrackParamAtHit->GetHitForRecPtr()->GetChamberNumber(); | |
980 | if (worstChamber%2 == 0) ((AliMUONTrackParam*)track->GetTrackParamAtHit()->After(worstTrackParamAtHit))->SetRemovable(kFALSE); | |
981 | else ((AliMUONTrackParam*)track->GetTrackParamAtHit()->Before(worstTrackParamAtHit))->SetRemovable(kFALSE); | |
982 | ||
983 | // Save pointer to the trackParamAtHit next to the one to be removed | |
984 | trackParamAtHit = (AliMUONTrackParam*)track->GetTrackParamAtHit()->After(worstTrackParamAtHit); | |
985 | ||
986 | // Remove the worst hit | |
987 | track->RemoveTrackParamAtHit(worstTrackParamAtHit); | |
988 | ||
989 | // Re-calculate track parameters | |
990 | // - from the hit immediately downstream the one suppressed | |
991 | // - or from the begining - if parameters have been re-computed using the standard method (kalman parameters have been lost) | |
992 | // - or if the removed hit was the last one | |
993 | if (smoothed && trackParamAtHit) RetracePartialTrack(*track,trackParamAtHit); | |
994 | else RetraceTrack(*track,kTRUE); | |
995 | ||
996 | // Printout for debuging | |
997 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { | |
998 | cout << "ImproveTracks: track " << fRecTracksPtr->IndexOf(track)+1 << " improved " << endl; | |
999 | } | |
1000 | ||
1001 | } | |
1002 | ||
1003 | track = (AliMUONTrack*) fRecTracksPtr->After(track); | |
1004 | } | |
1005 | ||
8cc77df0 | 1006 | } |
1007 | ||
b0e4c598 | 1008 | //__________________________________________________________________________ |
ea94c18b | 1009 | void AliMUONTrackReconstructorK::Finalize() |
b0e4c598 | 1010 | { |
ea94c18b | 1011 | /// Fill AliMUONTrack's fHitForRecAtHit array |
1012 | AliMUONTrack *track; | |
1013 | AliMUONTrackParam *trackParamAtHit; | |
1014 | Bool_t smoothed = kFALSE; | |
1015 | ||
1016 | track = (AliMUONTrack*) fRecTracksPtr->First(); | |
1017 | while (track) { | |
1018 | ||
1019 | // update track parameters (using smoother if required) if not already done | |
1020 | if (!track->IsImproved()) { | |
1021 | smoothed = kFALSE; | |
1022 | if (fgkRunSmoother) smoothed = RunSmoother(*track); | |
1023 | if (!smoothed) track->UpdateCovTrackParamAtHit(); | |
1024 | } | |
1025 | ||
1026 | trackParamAtHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First()); | |
1027 | while (trackParamAtHit) { | |
1028 | ||
1029 | // copy smoothed parameters and covariances if any | |
1030 | if (smoothed) { | |
1031 | trackParamAtHit->SetParameters(trackParamAtHit->GetSmoothParameters()); | |
1032 | trackParamAtHit->SetCovariances(trackParamAtHit->GetSmoothCovariances()); | |
1033 | } | |
1034 | ||
1035 | // update array of track hits | |
1036 | track->AddHitForRecAtHit(trackParamAtHit->GetHitForRecPtr()); | |
1037 | ||
1038 | trackParamAtHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->After(trackParamAtHit)); | |
1039 | } | |
1040 | ||
1041 | track = (AliMUONTrack*) fRecTracksPtr->After(track); | |
1042 | ||
b0e4c598 | 1043 | } |
ea94c18b | 1044 | |
b0e4c598 | 1045 | } |
ea94c18b | 1046 |