Adding comment lines to class description needed for Root documentation,
[u/mrichter/AliRoot.git] / MUON / AliMUONTrackReconstructorK.cxx
CommitLineData
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 45ClassImp(AliMUONTrackReconstructorK) // Class implementation in ROOT context
8cde4af5 46/// \endcond
8cc77df0 47
ea94c18b 48//************* Defaults parameters for reconstruction
49const Bool_t AliMUONTrackReconstructorK::fgkRunSmoother = kTRUE;
8cc77df0 50
8cc77df0 51
019df241 52 //__________________________________________________________________________
ea94c18b 53AliMUONTrackReconstructorK::AliMUONTrackReconstructorK()
54 : AliMUONVTrackReconstructor()
8cc77df0 55{
ea94c18b 56 /// Constructor for class AliMUONTrackReconstructorK
57 AliInfo("*** Tracking with Kalman Filter ***");
8cc77df0 58}
59
019df241 60 //__________________________________________________________________________
ea94c18b 61AliMUONTrackReconstructorK::~AliMUONTrackReconstructorK()
8cc77df0 62{
ea94c18b 63/// Destructor
64}
8cc77df0 65
66 //__________________________________________________________________________
ea94c18b 67void 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 //__________________________________________________________________________
153void 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 //__________________________________________________________________________
192void 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 //__________________________________________________________________________
265void 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 //__________________________________________________________________________
316Bool_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 //__________________________________________________________________________
648Double_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 //__________________________________________________________________________
715void 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 //__________________________________________________________________________
737void 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 //__________________________________________________________________________
784Bool_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 827Bool_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
922 //__________________________________________________________________________
ea94c18b 923void 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 1021void 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