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