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