]> git.uio.no Git - u/mrichter/AliRoot.git/blob - MUON/AliMUONTrackReconstructorK.cxx
08e452b29dafbceea5814fe45c26d6ae0aeda8bb
[u/mrichter/AliRoot.git] / MUON / AliMUONTrackReconstructorK.cxx
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
18 //-----------------------------------------------------------------------------
19 /// \class AliMUONTrackReconstructorK
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 ///
29 //-----------------------------------------------------------------------------
30
31 #include "AliMUONTrackReconstructorK.h"
32 #include "AliMUONConstants.h"
33 #include "AliMUONHitForRec.h"
34 #include "AliMUONTrack.h"
35 #include "AliMUONTrackParam.h"
36 #include "AliMUONTrackExtrap.h"
37
38 #include "AliLog.h"
39
40 #include <Riostream.h>
41 #include <TMath.h>
42 #include <TMatrixD.h>
43
44 /// \cond CLASSIMP
45 ClassImp(AliMUONTrackReconstructorK) // Class implementation in ROOT context
46 /// \endcond
47
48 //************* Defaults parameters for reconstruction
49 const Bool_t AliMUONTrackReconstructorK::fgkRunSmoother = kTRUE;
50
51
52   //__________________________________________________________________________
53 AliMUONTrackReconstructorK::AliMUONTrackReconstructorK()
54   : AliMUONVTrackReconstructor()
55 {
56   /// Constructor for class AliMUONTrackReconstructorK
57   AliInfo("*** Tracking with Kalman Filter ***");
58 }
59
60   //__________________________________________________________________________
61 AliMUONTrackReconstructorK::~AliMUONTrackReconstructorK()
62 {
63 /// Destructor
64
65
66   //__________________________________________________________________________
67 void AliMUONTrackReconstructorK::MakeTrackCandidates()
68 {
69   /// To make track candidates (assuming linear propagation if the flag fgkMakeTrackCandidatesFast is set to kTRUE):
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.
73   TClonesArray *segments;
74   AliMUONTrack *track;
75   Int_t iCandidate = 0;
76   Bool_t hitFound;
77
78   AliDebug(1,"Enter MakeTrackCandidates");
79
80   // Loop over stations(1..) 5 and 4 and make track candidates
81   for (Int_t istat=4; istat>=3; istat--) {
82     
83     // Make segments in the station
84     segments = MakeSegmentsInStation(istat);
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       
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
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       
109       // Remove track if no hit found
110       if (!hitFound) {
111         fRecTracksPtr->Remove(track);
112         fNRecTracks--;
113       }
114       
115     }
116     // delete the array of segments
117     delete segments;
118   }
119   
120   
121   // Retrace tracks using Kalman filter and remove bad ones
122   if (fgkMakeTrackCandidatesFast) {
123     fRecTracksPtr->Compress(); // this is essential before checking tracks
124     
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       
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   
150 }
151
152   //__________________________________________________________________________
153 void AliMUONTrackReconstructorK::RetraceTrack(AliMUONTrack &trackCandidate, Bool_t resetSeed)
154 {
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 }
190
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 }
263
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--;
302       }
303       
304     }
305     
306     fRecTracksPtr->Compress(); // this is essential before checking tracks
307     
308     // Keep only the best tracks if required
309     if (!fgkTrackAllTracks) RemoveDoubleTracks();
310     
311   }
312   
313 }
314
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   
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   
357   // Get track parameters
358   AliMUONTrackParam extrapTrackParamAtCh(*(AliMUONTrackParam*)trackCandidate.GetTrackParamAtHit()->First());
359   
360   // Add MCS effect
361   AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(),1.);
362   
363   // reset propagator for smoother
364   if (fgkRunSmoother) extrapTrackParamAtCh.ResetPropagator();
365   
366   // Add MCS in the missing chamber if any (only 1 chamber can be missing according to tracking criteria)
367   if (ch1 < ch2 && extrapTrackParamAtCh.GetHitForRecPtr()->GetChamberNumber() > ch2 + 1) {
368     // extrapolation to the missing chamber
369     AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(ch2 + 1), fgkRunSmoother);
370     // add MCS effect
371     AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(),1.);
372   }
373   
374   //Extrapolate trackCandidate to chamber "ch2"
375   AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(ch2), fgkRunSmoother);
376   
377   // Printout for debuging
378   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
379     cout<<endl<<"Track parameter covariances at first hit extrapolated to z = "<<AliMUONConstants::DefaultChamberZ(ch2)<<":"<<endl;
380     extrapTrackParamAtCh.GetCovariances().Print();
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     
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);
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;
408       }
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       
430       //Extrapolate track parameters to chamber "ch1"
431       AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParam, AliMUONConstants::DefaultChamberZ(ch1), fgkRunSmoother);
432       
433       for (Int_t hit1 = 0; hit1 < fNHitsForRecPerChamber[ch1]; hit1++) {
434         
435         hitForRecCh1 = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[ch1]+hit1);
436         
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         
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           }
453           
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           
492         }
493         
494       }
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) {
517           // keep track of the best single hitForRec except if a couple of hits has already been found
518           bestAddChi2TrackAtHit1 = addChi2TrackAtHit2;
519           bestTrackParamAtHit1 = extrapTrackParamAtHit2;
520         }
521         
522       }
523       
524     }
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;
535     }
536     
537     // add MCS effect for next step
538     AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(),1.);
539     
540     //Extrapolate trackCandidate to chamber "ch1"
541     AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(ch1), fgkRunSmoother);
542     
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       
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       
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         }
565         
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) {
594           // keep track of the best single hitForRec except if a couple of hits has already been found
595           bestAddChi2TrackAtHit1 = addChi2TrackAtHit1;
596           bestTrackParamAtHit1 = extrapTrackParamAtHit1;
597         }
598         
599       }
600       
601     }
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);
613
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 }
646
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;
671   }
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();
678
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   
712 }
713
714   //__________________________________________________________________________
715 void AliMUONTrackReconstructorK::UpdateTrack(AliMUONTrack &track, AliMUONTrackParam &trackParamAtHit, Double_t addChi2)
716 {
717   /// Add 1 hit to the track candidate
718   /// Update chi2 of the track 
719   
720   // Flag hit as being not removable
721   trackParamAtHit.SetRemovable(kFALSE);
722   trackParamAtHit.SetLocalChi2(0.); // --> Local chi2 not used
723   
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 }
735
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 }
782
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;
809     }
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   
824 }
825
826   //__________________________________________________________________________
827 Bool_t AliMUONTrackReconstructorK::RunSmoother(AliMUONTrack &track)
828 {
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   
920 }
921
922   //__________________________________________________________________________
923 void AliMUONTrackReconstructorK::ImproveTracks()
924 {
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();
938   while (track) {
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   
1018 }
1019
1020   //__________________________________________________________________________
1021 void AliMUONTrackReconstructorK::Finalize()
1022 {
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     
1055   }
1056     
1057 }
1058