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