96aac39a8acce958887d0c73aa04124ef01f1c08
[u/mrichter/AliRoot.git] / MUON / AliMUONTrackReconstructor.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 AliMUONTrackReconstructor
20 /// MUON track reconstructor using the original method
21 ///
22 /// This class contains as data:
23 /// - the parameters for the track reconstruction
24 ///
25 /// It contains as methods, among others:
26 /// - MakeTracks to build the tracks
27 //-----------------------------------------------------------------------------
28
29 #include "AliMUONTrackReconstructor.h"
30
31 #include "AliMUONConstants.h"
32 #include "AliMUONHitForRec.h"
33 #include "AliMUONTrack.h"
34 #include "AliMUONTrackParam.h"
35 #include "AliMUONTrackExtrap.h"
36
37 #include "AliLog.h"
38
39 #include <TMinuit.h>
40 #include <Riostream.h>
41 #include <TMath.h>
42 #include <TMatrixD.h>
43
44 // Functions to be minimized with Minuit
45 void TrackChi2(Int_t &nParam, Double_t *gradient, Double_t &chi2, Double_t *param, Int_t flag);
46
47 /// \cond CLASSIMP
48 ClassImp(AliMUONTrackReconstructor) // Class implementation in ROOT context
49 /// \endcond
50
51   //__________________________________________________________________________
52 AliMUONTrackReconstructor::AliMUONTrackReconstructor()
53   : AliMUONVTrackReconstructor()
54 {
55   /// Constructor for class AliMUONTrackReconstructor
56   AliInfo("*** Original tracking ***");
57 }
58
59   //__________________________________________________________________________
60 AliMUONTrackReconstructor::~AliMUONTrackReconstructor()
61 {
62 /// Destructor
63
64
65   //__________________________________________________________________________
66 void AliMUONTrackReconstructor::MakeTrackCandidates()
67 {
68   /// To make track candidates (assuming linear propagation if the flag fgkMakeTrackCandidatesFast is set to kTRUE):
69   /// Start with segments station(1..) 4 or 5 then follow track in station 5 or 4.
70   /// Good candidates are made of at least three hitForRec's.
71   /// Keep only best candidates or all of them according to the flag fgkTrackAllTracks.
72   
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     {
89       AliDebug(1,Form("Making primary candidate(1..) %d",++iCandidate));
90       
91       // Transform segments to tracks and put them at the end of fRecTracksPtr
92       track = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack((AliMUONObjectPair*)((*segments)[iseg]));
93       fNRecTracks++;
94       
95       // Printout for debuging
96       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2))
97       {
98         cout<<endl<<"Track parameter covariances at first hit with multiple Coulomb scattering effects:"<<endl;
99         ((AliMUONTrackParam*) track->GetTrackParamAtHit()->First())->GetCovariances().Print();
100       }
101       
102       // Look for compatible hitForRec(s) in the other station
103       if (AliMUONReconstructor::GetRecoParam()->MakeTrackCandidatesFast()) hitFound = FollowLinearTrackInStation(*track,7-istat);
104       else hitFound = FollowTrackInStation(*track,7-istat);
105       
106       // Remove track if no hit found
107       if (!hitFound) {
108         fRecTracksPtr->Remove(track);
109         fNRecTracks--;
110       }
111       
112     }
113     
114     // delete the array of segments
115     delete segments;
116   }
117   
118   fRecTracksPtr->Compress(); // this is essential before checking tracks
119   
120   // Keep all different tracks or only the best ones as required
121   if (AliMUONReconstructor::GetRecoParam()->TrackAllTracks()) RemoveIdenticalTracks();
122   else RemoveDoubleTracks();
123   
124   AliDebug(1,Form("Number of good candidates = %d",fNRecTracks));
125   
126 }
127
128   //__________________________________________________________________________
129 void AliMUONTrackReconstructor::FollowTracks()
130 {
131   /// Follow tracks in stations(1..) 3, 2 and 1
132   AliDebug(1,"Enter FollowTracks");
133   
134   AliMUONTrack *track, *nextTrack;
135   AliMUONTrackParam *trackParam, *nextTrackParam;
136   Int_t currentNRecTracks;
137   Bool_t hitFound;
138   
139   Double_t sigmaCut2 = AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking() *
140                        AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking();
141   
142   for (Int_t station = 2; station >= 0; station--) {
143     
144     // Save the actual number of reconstructed track in case of
145     // tracks are added or suppressed during the tracking procedure
146     // !! Do not compress fRecTracksPtr until the end of the loop over tracks !!
147     currentNRecTracks = fNRecTracks;
148     
149     for (Int_t iRecTrack = 0; iRecTrack <currentNRecTracks; iRecTrack++) {
150       AliDebug(1,Form("FollowTracks: track candidate(1..) %d", iRecTrack+1));
151       
152       track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iRecTrack);
153       
154       // Fit the track:
155       // Do not take into account the multiple scattering to speed up the fit
156       // Calculate the track parameter covariance matrix
157       // If "station" is station(1..) 3 then use the vertex to better constrain the fit
158       if (station==2) {
159         SetVertexForFit(*track);
160         track->SetFitWithVertex(kTRUE);
161       } else track->SetFitWithVertex(kFALSE);
162       Fit(*track, kFALSE, kTRUE);
163       
164       // Remove the track if the normalized chi2 is too high
165       if (track->GetNormalizedChi2() > sigmaCut2) {
166         fRecTracksPtr->Remove(track);
167         fNRecTracks--;
168         continue;
169       }
170       
171       // save parameters from fit into smoothed parameters to complete track afterward
172       if (AliMUONReconstructor::GetRecoParam()->ComplementTracks()) {
173         
174         if (station==2) { // save track parameters on stations 4 and 5
175           
176           // extrapolate track parameters and covariances at each cluster
177           track->UpdateCovTrackParamAtHit();
178           
179           // save them
180           trackParam = (AliMUONTrackParam*) track->GetTrackParamAtHit()->First();
181           while (trackParam) {
182             trackParam->SetSmoothParameters(trackParam->GetParameters());
183             trackParam->SetSmoothCovariances(trackParam->GetCovariances());
184             trackParam = (AliMUONTrackParam*) track->GetTrackParamAtHit()->After(trackParam);
185           }
186           
187         } else { // or save track parameters on last station only
188           
189           // save parameters from fit
190           trackParam = (AliMUONTrackParam*) track->GetTrackParamAtHit()->First();
191           trackParam->SetSmoothParameters(trackParam->GetParameters());
192           trackParam->SetSmoothCovariances(trackParam->GetCovariances());
193           
194           // save parameters extrapolated to the second chamber of the same station if it has been hit
195           nextTrackParam = (AliMUONTrackParam*) track->GetTrackParamAtHit()->After(trackParam);
196           if (nextTrackParam->GetHitForRecPtr()->GetChamberNumber() < 2*(station+2)) {
197             
198             // reset parameters and covariances
199             nextTrackParam->SetParameters(trackParam->GetParameters());
200             nextTrackParam->SetZ(trackParam->GetZ());
201             nextTrackParam->SetCovariances(trackParam->GetCovariances());
202             
203             // extrapolate them to the z of the corresponding cluster
204             AliMUONTrackExtrap::ExtrapToZCov(nextTrackParam, nextTrackParam->GetHitForRecPtr()->GetZ());
205             
206             // save them
207             nextTrackParam->SetSmoothParameters(nextTrackParam->GetParameters());
208             nextTrackParam->SetSmoothCovariances(nextTrackParam->GetCovariances());
209             
210           }
211           
212         }
213         
214       }
215       
216       // Printout for debuging
217       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
218         cout<<endl<<"Track parameter covariances at first hit with multiple Coulomb scattering effects:"<<endl;
219         ((AliMUONTrackParam*) track->GetTrackParamAtHit()->First())->GetCovariances().Print();
220       }
221       
222       // Look for compatible hitForRec in station(0..) "station"
223       hitFound = FollowTrackInStation(*track,station);
224       
225       // Try to recover track if required
226       if (!hitFound && AliMUONReconstructor::GetRecoParam()->RecoverTracks()) hitFound = RecoverTrack(*track,station);
227       
228       // remove track if no hit found
229       if (!hitFound) {
230         fRecTracksPtr->Remove(track);
231         fNRecTracks--;
232       }
233       
234     }
235     
236     // Compress fRecTracksPtr for the next step
237     fRecTracksPtr->Compress();
238     
239     // Keep only the best tracks if required
240     if (!AliMUONReconstructor::GetRecoParam()->TrackAllTracks()) RemoveDoubleTracks();
241     
242   }
243   
244   // Last fit of track candidates with all station
245   // Take into account the multiple scattering and remove bad tracks
246   Int_t trackIndex = -1;
247   track = (AliMUONTrack*) fRecTracksPtr->First();
248   while (track) {
249     
250     trackIndex++;
251     nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track); // prepare next track
252     
253     track->SetFitWithVertex(kFALSE); // just to be sure
254     Fit(*track, kTRUE, kTRUE);
255     
256     // Printout for debuging
257     if (AliLog::GetGlobalDebugLevel() >= 3) {
258       cout << "FollowTracks: track candidate(0..) " << trackIndex << " after final fit" << endl;
259       track->RecursiveDump();
260     } 
261     
262     // Remove the track if the normalized chi2 is too high
263     if (track->GetNormalizedChi2() > sigmaCut2) {
264       fRecTracksPtr->Remove(track);
265       fNRecTracks--;
266     }
267     
268     // save parameters from fit into smoothed parameters to complete track afterward
269     if (AliMUONReconstructor::GetRecoParam()->ComplementTracks()) {
270       
271       // save parameters from fit
272       trackParam = (AliMUONTrackParam*) track->GetTrackParamAtHit()->First();
273       trackParam->SetSmoothParameters(trackParam->GetParameters());
274       trackParam->SetSmoothCovariances(trackParam->GetCovariances());
275       
276       // save parameters extrapolated to the second chamber of the same station if it has been hit
277       nextTrackParam = (AliMUONTrackParam*) track->GetTrackParamAtHit()->After(trackParam);
278       if (nextTrackParam->GetHitForRecPtr()->GetChamberNumber() < 2) {
279         
280         // reset parameters and covariances
281         nextTrackParam->SetParameters(trackParam->GetParameters());
282         nextTrackParam->SetZ(trackParam->GetZ());
283         nextTrackParam->SetCovariances(trackParam->GetCovariances());
284         
285         // extrapolate them to the z of the corresponding cluster
286         AliMUONTrackExtrap::ExtrapToZCov(nextTrackParam, nextTrackParam->GetHitForRecPtr()->GetZ());
287         
288         // save them
289         nextTrackParam->SetSmoothParameters(nextTrackParam->GetParameters());
290         nextTrackParam->SetSmoothCovariances(nextTrackParam->GetCovariances());
291         
292       }
293       
294     }
295     
296     track = nextTrack;
297     
298   }
299   
300   fRecTracksPtr->Compress();
301   
302 }
303
304   //__________________________________________________________________________
305 Bool_t AliMUONTrackReconstructor::FollowTrackInStation(AliMUONTrack &trackCandidate, Int_t nextStation)
306 {
307   /// Follow trackCandidate in station(0..) nextStation and search for compatible HitForRec(s)
308   /// Keep all possibilities or only the best one(s) according to the flag fgkTrackAllTracks:
309   /// kTRUE:  duplicate "trackCandidate" if there are several possibilities and add the new tracks at the end of
310   ///         fRecTracksPtr to avoid conficts with other track candidates at this current stage of the tracking procedure.
311   ///         Remove the obsolete "trackCandidate" at the end.
312   /// kFALSE: add only the best hit(s) to the "trackCandidate". Try to add a couple of hits in priority.
313   AliDebug(1,Form("Enter FollowTrackInStation(1..) %d", nextStation+1));
314   
315   // Order the chamber according to the propagation direction (tracking starts with chamber 2):
316   // - nextStation == station(1...) 5 => forward propagation
317   // - nextStation < station(1...) 5 => backward propagation
318   Int_t ch1, ch2;
319   if (nextStation==4) {
320     ch1 = 2*nextStation+1;
321     ch2 = 2*nextStation;
322   } else {
323     ch1 = 2*nextStation;
324     ch2 = 2*nextStation+1;
325   }
326   
327   Double_t chi2WithOneHitForRec = 1.e10;
328   Double_t chi2WithTwoHitForRec = 1.e10;
329   Double_t maxChi2WithOneHitForRec = 2. * AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking() *
330                                           AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking(); // 2 because 2 quantities in chi2
331   Double_t maxChi2WithTwoHitForRec = 4. * AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking() *
332                                           AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking(); // 4 because 4 quantities in chi2
333   Double_t bestChi2WithOneHitForRec = maxChi2WithOneHitForRec;
334   Double_t bestChi2WithTwoHitForRec = maxChi2WithTwoHitForRec;
335   Bool_t foundOneHit = kFALSE;
336   Bool_t foundTwoHits = kFALSE;
337   AliMUONTrack *newTrack = 0x0;
338   AliMUONHitForRec *hitForRecCh1, *hitForRecCh2;
339   AliMUONTrackParam extrapTrackParam;
340   AliMUONTrackParam extrapTrackParamAtHit1;
341   AliMUONTrackParam extrapTrackParamAtHit2;
342   AliMUONTrackParam bestTrackParamAtHit1;
343   AliMUONTrackParam bestTrackParamAtHit2;
344   Bool_t *hitForRecCh1Used = new Bool_t[fNHitsForRecPerChamber[ch1]];
345   for (Int_t hit1 = 0; hit1 < fNHitsForRecPerChamber[ch1]; hit1++) hitForRecCh1Used[hit1] = kFALSE;
346   
347   // Get track parameters
348   AliMUONTrackParam extrapTrackParamAtCh(*(AliMUONTrackParam*)trackCandidate.GetTrackParamAtHit()->First());
349   
350   // Add MCS effect
351   AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(),1.);
352   
353   // Add MCS in the missing chamber if any (only 1 chamber can be missing according to tracking criteria)
354   if (ch1 < ch2 && extrapTrackParamAtCh.GetHitForRecPtr()->GetChamberNumber() > ch2 + 1) {
355     // extrapolation to the missing chamber
356     AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(ch2 + 1));
357     // add MCS effect
358     AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(),1.);
359   }
360   
361   //Extrapolate trackCandidate to chamber "ch2"
362   AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(ch2));
363   
364   // Printout for debuging
365   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
366     cout<<endl<<"Track parameter covariances at first hit extrapolated to z = "<<AliMUONConstants::DefaultChamberZ(ch2)<<":"<<endl;
367     extrapTrackParamAtCh.GetCovariances().Print();
368   }
369   
370   // Printout for debuging
371   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
372     cout << "FollowTrackInStation: look for hits in chamber(1..): " << ch2+1 << endl;
373   }
374   
375   // look for candidates in chamber 2 
376   for (Int_t hit2 = 0; hit2 < fNHitsForRecPerChamber[ch2]; hit2++) {
377     
378     hitForRecCh2 = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[ch2]+hit2);
379     
380     // try to add the current hit fast
381     if (!TryOneHitForRecFast(extrapTrackParamAtCh, hitForRecCh2)) continue;
382     
383     // try to add the current hit accuratly
384     chi2WithOneHitForRec = TryOneHitForRec(extrapTrackParamAtCh, hitForRecCh2, extrapTrackParamAtHit2);
385     
386     // if good chi2 then try to attach a hitForRec in the other chamber too
387     if (chi2WithOneHitForRec < maxChi2WithOneHitForRec) {
388       Bool_t foundSecondHit = kFALSE;
389       
390       // Printout for debuging
391       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
392         cout << "FollowTrackInStation: found one hit in chamber(1..): " << ch2+1
393              << " (Chi2 = " << chi2WithOneHitForRec << ")" << endl;
394         cout << "                      look for second hits in chamber(1..): " << ch1+1 << " ..." << endl;
395       }
396       
397       // add MCS effect for next step
398       AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtHit2,AliMUONConstants::ChamberThicknessInX0(),1.);
399       
400       // copy new track parameters for next step
401       extrapTrackParam = extrapTrackParamAtHit2;
402       
403       //Extrapolate track parameters to chamber "ch1"
404       AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParam, AliMUONConstants::DefaultChamberZ(ch1));
405       
406       for (Int_t hit1 = 0; hit1 < fNHitsForRecPerChamber[ch1]; hit1++) {
407         
408         hitForRecCh1 = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[ch1]+hit1);
409         
410         // try to add the current hit fast
411         if (!TryOneHitForRecFast(extrapTrackParam, hitForRecCh1)) continue;
412         
413         // try to add the current hit accuratly
414         chi2WithTwoHitForRec = TryTwoHitForRec(extrapTrackParamAtHit2, hitForRecCh1, extrapTrackParamAtHit1);
415         
416         // if good chi2 then create a new track by adding the 2 hitForRec to the "trackCandidate"
417         if (chi2WithTwoHitForRec < maxChi2WithTwoHitForRec) {
418           foundSecondHit = kTRUE;
419           foundTwoHits = kTRUE;
420           
421           // Printout for debuging
422           if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
423             cout << "FollowTrackInStation: found second hit in chamber(1..): " << ch1+1
424                  << " (Global Chi2 = " << chi2WithTwoHitForRec << ")" << endl;
425           }
426           
427           if (AliMUONReconstructor::GetRecoParam()->TrackAllTracks()) {
428             // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new hitForRec's
429             newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
430             UpdateTrack(*newTrack,extrapTrackParamAtHit1,extrapTrackParamAtHit2);
431             fNRecTracks++;
432             
433             // Tag hitForRecCh1 as used
434             hitForRecCh1Used[hit1] = kTRUE;
435             
436             // Printout for debuging
437             if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
438               cout << "FollowTrackInStation: added two hits in station(1..): " << nextStation+1 << endl;
439               if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
440             }
441             
442           } else if (chi2WithTwoHitForRec < bestChi2WithTwoHitForRec) {
443             // keep track of the best couple of hits
444             bestChi2WithTwoHitForRec = chi2WithTwoHitForRec;
445             bestTrackParamAtHit1 = extrapTrackParamAtHit1;
446             bestTrackParamAtHit2 = extrapTrackParamAtHit2;
447           }
448           
449         }
450         
451       }
452       
453       // if no hitForRecCh1 found then consider to add hitForRecCh2 only
454       if (!foundSecondHit) {
455         foundOneHit = kTRUE;
456         
457         if (AliMUONReconstructor::GetRecoParam()->TrackAllTracks()) {
458           // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new hitForRec's
459           newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
460           UpdateTrack(*newTrack,extrapTrackParamAtHit2);
461           fNRecTracks++;
462           
463           // Printout for debuging
464           if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
465             cout << "FollowTrackInStation: added one hit in chamber(1..): " << ch2+1 << endl;
466             if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
467           }
468           
469         } else if (!foundTwoHits && chi2WithOneHitForRec < bestChi2WithOneHitForRec) {
470           // keep track of the best single hitForRec except if a couple of hits has already been found
471           bestChi2WithOneHitForRec = chi2WithOneHitForRec;
472           bestTrackParamAtHit1 = extrapTrackParamAtHit2;
473         }
474         
475       }
476       
477     }
478     
479   }
480   
481   // look for candidates in chamber 1 not already attached to a track
482   // if we want to keep all possible tracks or if no good couple of hitForRec has been found
483   if (AliMUONReconstructor::GetRecoParam()->TrackAllTracks() || !foundTwoHits) {
484     
485     // Printout for debuging
486     if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
487       cout << "FollowTrackInStation: look for single hits in chamber(1..): " << ch1+1 << endl;
488     }
489     
490     // add MCS effect for next step
491     AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(),1.);
492     
493     //Extrapolate trackCandidate to chamber "ch1"
494     AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(ch1));
495     
496     for (Int_t hit1 = 0; hit1 < fNHitsForRecPerChamber[ch1]; hit1++) {
497       
498       hitForRecCh1 = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[ch1]+hit1);
499       
500       if (hitForRecCh1Used[hit1]) continue; // Skip hitForRec already used
501       
502       // try to add the current hit fast
503       if (!TryOneHitForRecFast(extrapTrackParamAtCh, hitForRecCh1)) continue;
504       
505       // try to add the current hit accuratly
506       chi2WithOneHitForRec = TryOneHitForRec(extrapTrackParamAtCh, hitForRecCh1, extrapTrackParamAtHit1);
507     
508       // if good chi2 then consider to add hitForRecCh1
509       // We do not try to attach a hitForRec in the other chamber too since it has already been done above
510       if (chi2WithOneHitForRec < maxChi2WithOneHitForRec) {
511         foundOneHit = kTRUE;
512           
513         // Printout for debuging
514         if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
515           cout << "FollowTrackInStation: found one hit in chamber(1..): " << ch1+1
516                << " (Chi2 = " << chi2WithOneHitForRec << ")" << endl;
517         }
518         
519         if (AliMUONReconstructor::GetRecoParam()->TrackAllTracks()) {
520           // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new hitForRec's
521           newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
522           UpdateTrack(*newTrack,extrapTrackParamAtHit1);
523           fNRecTracks++;
524           
525           // Printout for debuging
526           if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
527             cout << "FollowTrackInStation: added one hit in chamber(1..): " << ch1+1 << endl;
528             if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
529           }
530           
531         } else if (chi2WithOneHitForRec < bestChi2WithOneHitForRec) {
532           // keep track of the best single hitForRec except if a couple of hits has already been found
533           bestChi2WithOneHitForRec = chi2WithOneHitForRec;
534           bestTrackParamAtHit1 = extrapTrackParamAtHit1;
535         }
536         
537       }
538       
539     }
540     
541   }
542   
543   // fill out the best track if required else clean up the fRecTracksPtr array
544   if (!AliMUONReconstructor::GetRecoParam()->TrackAllTracks()) {
545     if (foundTwoHits) {
546       UpdateTrack(trackCandidate,bestTrackParamAtHit1,bestTrackParamAtHit2);
547       
548       // Printout for debuging
549       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
550         cout << "FollowTrackInStation: added the two best hits in station(1..): " << nextStation+1 << endl;
551         if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
552       }
553       
554     } else if (foundOneHit) {
555       UpdateTrack(trackCandidate,bestTrackParamAtHit1);
556       
557       // Printout for debuging
558       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
559         cout << "FollowTrackInStation: added the best hit in chamber(1..): " << bestTrackParamAtHit1.GetHitForRecPtr()->GetChamberNumber()+1 << endl;
560         if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
561       }
562       
563     } else {
564       delete [] hitForRecCh1Used;
565       return kFALSE;
566     }
567     
568   } else if (foundOneHit || foundTwoHits) {
569     
570     // remove obsolete track
571     fRecTracksPtr->Remove(&trackCandidate);
572     fNRecTracks--;
573     
574   } else {
575     delete [] hitForRecCh1Used;
576     return kFALSE;
577   }
578   
579   delete [] hitForRecCh1Used;
580   return kTRUE;
581   
582 }
583
584   //__________________________________________________________________________
585 Double_t AliMUONTrackReconstructor::TryTwoHitForRec(const AliMUONTrackParam &trackParamAtHit1, AliMUONHitForRec* hitForRec2,
586                                                     AliMUONTrackParam &trackParamAtHit2)
587 {
588 /// Test the compatibility between the track and the 2 hitForRec together (using trackParam's covariance matrix):
589 /// return the corresponding Chi2 accounting for covariances between the 2 hitForRec
590 /// return trackParamAtHit1 & 2
591   
592   // extrapolate track parameters at the z position of the second hit
593   trackParamAtHit2.SetParameters(trackParamAtHit1.GetParameters());
594   trackParamAtHit2.SetZ(trackParamAtHit1.GetZ());
595   AliMUONTrackExtrap::ExtrapToZ(&trackParamAtHit2, hitForRec2->GetZ());
596   
597   // set pointer to hit2 into trackParamAtHit2
598   trackParamAtHit2.SetHitForRecPtr(hitForRec2);
599   
600   // Set differences between track and the 2 hitForRec in the bending and non bending directions
601   AliMUONHitForRec* hitForRec1 = trackParamAtHit1.GetHitForRecPtr();
602   TMatrixD dPos(4,1);
603   dPos(0,0) = hitForRec1->GetNonBendingCoor() - trackParamAtHit1.GetNonBendingCoor();
604   dPos(1,0) = hitForRec1->GetBendingCoor() - trackParamAtHit1.GetBendingCoor();
605   dPos(2,0) = hitForRec2->GetNonBendingCoor() - trackParamAtHit2.GetNonBendingCoor();
606   dPos(3,0) = hitForRec2->GetBendingCoor() - trackParamAtHit2.GetBendingCoor();
607   
608   // Calculate the error matrix from the track parameter covariances at first hitForRec
609   TMatrixD error(4,4);
610   error.Zero();
611   if (trackParamAtHit1.CovariancesExist()) {
612     // Save track parameters at first hitForRec
613     AliMUONTrackParam trackParamAtHit1Save(trackParamAtHit1);
614     TMatrixD paramAtHit1Save(trackParamAtHit1Save.GetParameters());
615     Double_t z1 = trackParamAtHit1Save.GetZ();
616     
617     // Save track coordinates at second hitForRec
618     Double_t nonBendingCoor2         = trackParamAtHit2.GetNonBendingCoor();
619     Double_t bendingCoor2            = trackParamAtHit2.GetBendingCoor();
620     
621     // add MCS effect at first hitForRec
622     AliMUONTrackExtrap::AddMCSEffect(&trackParamAtHit1Save,AliMUONConstants::ChamberThicknessInX0(),1.);
623     
624     // Get the pointer to the parameter covariance matrix at first hitForRec
625     const TMatrixD& kParamCov = trackParamAtHit1Save.GetCovariances();
626     
627     // Calculate the jacobian related to the transformation between track parameters
628     // at first hitForRec and track coordinates at the 2 hitForRec z-position
629     TMatrixD jacob(4,5);
630     jacob.Zero();
631     // first derivative at the first hitForRec:
632     jacob(0,0) = 1.; // dx1/dx
633     jacob(1,2) = 1.; // dy1/dy
634     // first derivative at the second hitForRec:
635     TMatrixD dParam(5,1);
636     for (Int_t i=0; i<5; i++) {
637       // Skip jacobian calculation for parameters with no associated error
638       if (kParamCov(i,i) == 0.) continue;
639       // Small variation of parameter i only
640       for (Int_t j=0; j<5; j++) {
641         if (j==i) {
642           dParam(j,0) = TMath::Sqrt(kParamCov(i,i));
643           if (j == 4) dParam(j,0) *= TMath::Sign(1.,-paramAtHit1Save(4,0)); // variation always in the same direction
644         } else dParam(j,0) = 0.;
645       }
646       
647       // Set new track parameters at first hitForRec
648       trackParamAtHit1Save.SetParameters(paramAtHit1Save);
649       trackParamAtHit1Save.AddParameters(dParam);
650       trackParamAtHit1Save.SetZ(z1);
651       
652       // Extrapolate new track parameters to the z position of the second hitForRec
653       AliMUONTrackExtrap::ExtrapToZ(&trackParamAtHit1Save,hitForRec2->GetZ());
654       
655       // Calculate the jacobian
656       jacob(2,i) = (trackParamAtHit1Save.GetNonBendingCoor()  - nonBendingCoor2) / dParam(i,0); // dx2/dParami
657       jacob(3,i) = (trackParamAtHit1Save.GetBendingCoor()     - bendingCoor2   ) / dParam(i,0); // dy2/dParami
658     }
659     
660     // Calculate the error matrix
661     TMatrixD tmp(jacob,TMatrixD::kMult,kParamCov);
662     error = TMatrixD(tmp,TMatrixD::kMultTranspose,jacob);
663   }
664   
665   // Add hitForRec resolution to the error matrix
666   error(0,0) += hitForRec1->GetNonBendingReso2();
667   error(1,1) += hitForRec1->GetBendingReso2();
668   error(2,2) += hitForRec2->GetNonBendingReso2();
669   error(3,3) += hitForRec2->GetBendingReso2();
670   
671   // invert the error matrix for Chi2 calculation
672   if (error.Determinant() != 0) {
673     error.Invert();
674   } else {
675     AliWarning(" Determinant error=0");
676     return 1.e10;
677   }
678   
679   // Compute the Chi2 value
680   TMatrixD tmp2(dPos,TMatrixD::kTransposeMult,error);
681   TMatrixD result(tmp2,TMatrixD::kMult,dPos);
682   
683   return result(0,0);
684   
685 }
686
687   //__________________________________________________________________________
688 void AliMUONTrackReconstructor::UpdateTrack(AliMUONTrack &track, AliMUONTrackParam &trackParamAtHit)
689 {
690   /// Add 1 hit to the track candidate
691   /// Update chi2 of the track 
692   
693   // Compute local chi2
694   AliMUONHitForRec* hit = trackParamAtHit.GetHitForRecPtr();
695   Double_t deltaX = trackParamAtHit.GetNonBendingCoor() - hit->GetNonBendingCoor();
696   Double_t deltaY = trackParamAtHit.GetBendingCoor() - hit->GetBendingCoor();
697   Double_t localChi2 = deltaX*deltaX / hit->GetNonBendingReso2() +
698                        deltaY*deltaY / hit->GetBendingReso2();
699   
700   // Flag hit as being not removable
701   trackParamAtHit.SetRemovable(kFALSE);
702   trackParamAtHit.SetLocalChi2(0.); // --> Local chi2 not used
703   
704   // Update the chi2 of the new track
705   track.SetFitFMin(track.GetFitFMin() + localChi2);
706   
707   // Update TrackParamAtHit
708   track.AddTrackParamAtHit(&trackParamAtHit,trackParamAtHit.GetHitForRecPtr());
709   track.GetTrackParamAtHit()->Sort();
710   
711 }
712
713   //__________________________________________________________________________
714 void AliMUONTrackReconstructor::UpdateTrack(AliMUONTrack &track, AliMUONTrackParam &trackParamAtHit1, AliMUONTrackParam &trackParamAtHit2)
715 {
716   /// Add 2 hits to the track candidate
717   /// Update track and local chi2
718   
719   // Update local chi2 at first hit
720   AliMUONHitForRec* hit1 = trackParamAtHit1.GetHitForRecPtr();
721   Double_t deltaX = trackParamAtHit1.GetNonBendingCoor() - hit1->GetNonBendingCoor();
722   Double_t deltaY = trackParamAtHit1.GetBendingCoor() - hit1->GetBendingCoor();
723   Double_t localChi2AtHit1 = deltaX*deltaX / hit1->GetNonBendingReso2() +
724                              deltaY*deltaY / hit1->GetBendingReso2();
725   trackParamAtHit1.SetLocalChi2(localChi2AtHit1);
726   
727   // Flag first hit as being removable
728   trackParamAtHit1.SetRemovable(kTRUE);
729   
730   // Update local chi2 at second hit
731   AliMUONHitForRec* hit2 = trackParamAtHit2.GetHitForRecPtr();
732   deltaX = trackParamAtHit2.GetNonBendingCoor() - hit2->GetNonBendingCoor();
733   deltaY = trackParamAtHit2.GetBendingCoor() - hit2->GetBendingCoor();
734   Double_t localChi2AtHit2 = deltaX*deltaX / hit2->GetNonBendingReso2() +
735                              deltaY*deltaY / hit2->GetBendingReso2();
736   trackParamAtHit2.SetLocalChi2(localChi2AtHit2);
737   
738   // Flag first hit as being removable
739   trackParamAtHit2.SetRemovable(kTRUE);
740   
741   // Update the chi2 of the new track
742   track.SetFitFMin(track.GetFitFMin() + localChi2AtHit1 + localChi2AtHit2);
743   
744   // Update TrackParamAtHit
745   track.AddTrackParamAtHit(&trackParamAtHit1,trackParamAtHit1.GetHitForRecPtr());
746   track.AddTrackParamAtHit(&trackParamAtHit2,trackParamAtHit2.GetHitForRecPtr());
747   track.GetTrackParamAtHit()->Sort();
748   
749 }
750
751   //__________________________________________________________________________
752 Bool_t AliMUONTrackReconstructor::RecoverTrack(AliMUONTrack &trackCandidate, Int_t nextStation)
753 {
754   /// Try to recover the track candidate in the next station
755   /// by removing the worst of the two hits attached in the current station
756   /// Return kTRUE if recovering succeeds
757   AliDebug(1,"Enter RecoverTrack");
758   
759   // Do not try to recover track until we have attached hit(s) on station(1..) 3
760   if (nextStation > 1) return kFALSE;
761   
762   Int_t worstHitNumber = -1;
763   Double_t localChi2, worstLocalChi2 = 0.;
764   
765   // Look for the hit to remove
766   for (Int_t hitNumber = 0; hitNumber < 2; hitNumber++) {
767     AliMUONTrackParam *trackParamAtHit = (AliMUONTrackParam*)trackCandidate.GetTrackParamAtHit()->UncheckedAt(hitNumber);
768     
769     // check if current hit is removable
770     if (!trackParamAtHit->IsRemovable()) return kFALSE;
771     
772     // Pick up hit with the worst chi2
773     localChi2 = trackParamAtHit->GetLocalChi2();
774     if (localChi2 > worstLocalChi2) {
775       worstLocalChi2 = localChi2;
776       worstHitNumber = hitNumber;
777     }
778   }
779   
780   // Reset best hit as being NOT removable
781   ((AliMUONTrackParam*)trackCandidate.GetTrackParamAtHit()->UncheckedAt((worstHitNumber+1)%2))->SetRemovable(kFALSE);
782   
783   // Remove the worst hit
784   trackCandidate.RemoveTrackParamAtHit((AliMUONTrackParam*)trackCandidate.GetTrackParamAtHit()->UncheckedAt(worstHitNumber));
785   
786   // Re-fit the track:
787   // Do not take into account the multiple scattering to speed up the fit
788   // Calculate the track parameter covariance matrix
789   trackCandidate.SetFitWithVertex(kFALSE); // To be sure
790   Fit(trackCandidate, kFALSE, kTRUE);
791   
792   // Look for new hit(s) in next station
793   return FollowTrackInStation(trackCandidate,nextStation);
794   
795 }
796
797   //__________________________________________________________________________
798 void AliMUONTrackReconstructor::SetVertexForFit(AliMUONTrack &trackCandidate)
799 {
800   /// Add the vertex as a measured hit to constrain the fit of the "trackCandidate"
801   /// Compute the vertex resolution from natural vertex dispersion and
802   /// multiple scattering effets according to trackCandidate path in absorber
803   /// It is necessary to account for multiple scattering effects here instead of during the fit of
804   /// the "trackCandidate" to do not influence the result by changing track resolution at vertex
805   AliDebug(1,"Enter SetVertexForFit");
806   
807   Double_t nonBendingReso2 = AliMUONReconstructor::GetRecoParam()->GetNonBendingVertexDispersion() *
808                              AliMUONReconstructor::GetRecoParam()->GetNonBendingVertexDispersion();
809   Double_t bendingReso2 = AliMUONReconstructor::GetRecoParam()->GetBendingVertexDispersion() *
810                           AliMUONReconstructor::GetRecoParam()->GetBendingVertexDispersion();
811   // add multiple scattering effets
812   AliMUONTrackParam paramAtVertex(*((AliMUONTrackParam*)(trackCandidate.GetTrackParamAtHit()->First())));
813   paramAtVertex.DeleteCovariances(); // to be sure to account only for multiple scattering
814   AliMUONTrackExtrap::ExtrapToVertexUncorrected(&paramAtVertex,0.);
815   const TMatrixD& kParamCov = paramAtVertex.GetCovariances();
816   nonBendingReso2 += kParamCov(0,0);
817   bendingReso2 += kParamCov(2,2);
818   // Set the vertex
819   AliMUONHitForRec vertex; // Coordinates set to (0.,0.,0.) by default
820   vertex.SetNonBendingReso2(nonBendingReso2);
821   vertex.SetBendingReso2(bendingReso2);
822   trackCandidate.SetVertex(&vertex);
823 }
824
825   //__________________________________________________________________________
826 void AliMUONTrackReconstructor::Fit(AliMUONTrack &track, Bool_t includeMCS, Bool_t calcCov)
827 {
828   /// Fit the track "track" w/wo multiple Coulomb scattering according to "includeMCS".
829   
830   Double_t benC, errorParam, invBenP, nonBenC, x, y;
831   AliMUONTrackParam *trackParam;
832   Double_t arg[1], fedm, errdef, fitFMin;
833   Int_t npari, nparx;
834   Int_t status, covStatus;
835   
836   // Instantiate gMinuit if not already done
837   if (!gMinuit) gMinuit = new TMinuit(6);
838   // Clear MINUIT parameters
839   gMinuit->mncler();
840   // Give the fitted track to MINUIT
841   gMinuit->SetObjectFit(&track);
842   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
843     // Define print level
844     arg[0] = 1;
845     gMinuit->mnexcm("SET PRI", arg, 1, status);
846     // Print covariance matrix
847     gMinuit->mnexcm("SHO COV", arg, 0, status);
848   } else {
849     arg[0] = -1;
850     gMinuit->mnexcm("SET PRI", arg, 1, status);
851   }
852   // No warnings
853   gMinuit->mnexcm("SET NOW", arg, 0, status);
854   // Define strategy
855   //arg[0] = 2;
856   //gMinuit->mnexcm("SET STR", arg, 1, status);
857   
858   // set flag w/wo multiple scattering according to "includeMCS"
859   track.SetFitWithMCS(includeMCS);
860   if (includeMCS) {
861     // compute hit weights only once
862     if (!track.ComputeHitWeights()) {
863       AliWarning("cannot take into account the multiple scattering effects");
864       track.SetFitWithMCS(kFALSE);
865     }
866   }
867   
868   // Set fitting function
869   gMinuit->SetFCN(TrackChi2);
870   
871   // Set fitted parameters (!! The order is very important for the covariance matrix !!)
872   trackParam = (AliMUONTrackParam*) (track.GetTrackParamAtHit()->First());
873   // could be tried with no limits for the search (min=max=0) ????
874   // mandatory limits in non Bending to avoid NaN values of parameters
875   gMinuit->mnparm(0, "X", trackParam->GetNonBendingCoor(), 0.03, -500.0, 500.0, status);
876   gMinuit->mnparm(1, "NonBenS", trackParam->GetNonBendingSlope(), 0.001, -0.5, 0.5, status);
877   // mandatory limits in Bending to avoid NaN values of parameters
878   gMinuit->mnparm(2, "Y", trackParam->GetBendingCoor(), 0.10, -500.0, 500.0, status);
879   gMinuit->mnparm(3, "BenS", trackParam->GetBendingSlope(), 0.001, -0.5, 0.5, status);
880   gMinuit->mnparm(4, "InvBenP", trackParam->GetInverseBendingMomentum(), 0.003, -0.4, 0.4, status);
881   
882   // minimization
883   gMinuit->mnexcm("MIGRAD", arg, 0, status);
884   
885   // Calculate the covariance matrix more accurately if required
886   if (calcCov) gMinuit->mnexcm("HESSE", arg, 0, status);
887    
888   // get results into "invBenP", "benC", "nonBenC" ("x", "y")
889   gMinuit->GetParameter(0, x, errorParam);
890   trackParam->SetNonBendingCoor(x);
891   gMinuit->GetParameter(1, nonBenC, errorParam);
892   trackParam->SetNonBendingSlope(nonBenC);
893   gMinuit->GetParameter(2, y, errorParam);
894   trackParam->SetBendingCoor(y);
895   gMinuit->GetParameter(3, benC, errorParam);
896   trackParam->SetBendingSlope(benC);
897   gMinuit->GetParameter(4, invBenP, errorParam);
898   trackParam->SetInverseBendingMomentum(invBenP);
899   
900   // global result of the fit
901   gMinuit->mnstat(fitFMin, fedm, errdef, npari, nparx, covStatus);
902   track.SetFitFMin(fitFMin);
903   
904   // Get the covariance matrix if required
905   if (calcCov) {
906     // Covariance matrix according to HESSE status
907     // If problem then keep only the diagonal terms (variances)
908     Double_t matrix[5][5];
909     gMinuit->mnemat(&matrix[0][0],5);
910     if (covStatus == 3) trackParam->SetCovariances(matrix);
911     else trackParam->SetVariances(matrix);
912   } else trackParam->DeleteCovariances();
913   
914 }
915
916   //__________________________________________________________________________
917 void TrackChi2(Int_t & /*nParam*/, Double_t * /*gradient*/, Double_t &chi2, Double_t *param, Int_t /*flag*/)
918 {
919   /// Return the "Chi2" to be minimized with Minuit for track fitting.
920   /// Assumes that the track hits are sorted according to increasing Z.
921   /// Track parameters at each TrackHit are updated accordingly.
922   /// Vertex is used according to the flag "trackBeingFitted->GetFitWithVertex()".
923   /// Multiple Coulomb scattering is taken into account according to the flag "trackBeingFitted->GetFitWithMCS()".
924   
925   AliMUONTrack *trackBeingFitted = (AliMUONTrack*) gMinuit->GetObjectFit();
926   AliMUONTrackParam* trackParamAtHit = (AliMUONTrackParam*) trackBeingFitted->GetTrackParamAtHit()->First();
927   Double_t dX, dY;
928   chi2 = 0.; // initialize chi2
929   
930   // update track parameters
931   trackParamAtHit->SetNonBendingCoor(param[0]);
932   trackParamAtHit->SetNonBendingSlope(param[1]);
933   trackParamAtHit->SetBendingCoor(param[2]);
934   trackParamAtHit->SetBendingSlope(param[3]);
935   trackParamAtHit->SetInverseBendingMomentum(param[4]);
936   trackBeingFitted->UpdateTrackParamAtHit();
937   
938   // Take the vertex into account in the fit if required
939   if (trackBeingFitted->GetFitWithVertex()) {
940     AliMUONTrackParam paramAtVertex(*trackParamAtHit);
941     AliMUONTrackExtrap::ExtrapToZ(&paramAtVertex, 0.);
942     AliMUONHitForRec *vertex = trackBeingFitted->GetVertex();
943     if (!vertex) {
944       cout<<"Error in TrackChi2: Want to use the vertex in tracking but it has not been created!!"<<endl;
945       exit(-1);
946     }
947     dX = vertex->GetNonBendingCoor() - paramAtVertex.GetNonBendingCoor();
948     dY = vertex->GetBendingCoor() - paramAtVertex.GetBendingCoor();
949     chi2 += dX * dX / vertex->GetNonBendingReso2() + dY * dY / vertex->GetBendingReso2();
950   }
951   
952   // compute chi2 w/wo multiple scattering
953   chi2 += trackBeingFitted->ComputeGlobalChi2(trackBeingFitted->GetFitWithMCS());
954   
955 }
956
957   //__________________________________________________________________________
958 void AliMUONTrackReconstructor::ComplementTracks()
959 {
960   /// Complete tracks by adding missing clusters (if there is an overlap between
961   /// two detection elements, the track may have two clusters in the same chamber)
962   /// Re-fit track parameters and covariances
963   AliDebug(1,"Enter ComplementTracks");
964   
965   Int_t chamberId, detElemId;
966   Double_t chi2OfHitForRec, bestChi2OfHitForRec;
967   Bool_t foundOneHit, trackModified;
968   AliMUONHitForRec* hitForRec;
969   AliMUONTrackParam* trackParam, *nextTrackParam;
970   AliMUONTrackParam copyOfTrackParam;
971   AliMUONTrackParam trackParamAtHit;
972   AliMUONTrackParam bestTrackParamAtHit;
973   Double_t sigmaCut2 = AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking() *
974                        AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking();
975   
976   // Remove double track to complete only "good" tracks
977   RemoveDoubleTracks();
978   
979   AliMUONTrack *track = (AliMUONTrack*) fRecTracksPtr->First();
980   while (track) {
981     trackModified = kFALSE;
982     
983     trackParam = (AliMUONTrackParam*)track->GetTrackParamAtHit()->First();
984     while (trackParam) {
985       foundOneHit = kFALSE;
986       bestChi2OfHitForRec = 2. * sigmaCut2; // 2 because 2 quantities in chi2
987       
988       // prepare nextTrackParam before adding new cluster because of the sorting
989       nextTrackParam = (AliMUONTrackParam*)track->GetTrackParamAtHit()->After(trackParam);
990       
991       chamberId = trackParam->GetHitForRecPtr()->GetChamberNumber();
992       detElemId = trackParam->GetHitForRecPtr()->GetDetElemId();
993       
994       // recover track parameters from local fit and put them into a copy of trackParam
995       copyOfTrackParam.SetZ(trackParam->GetZ());
996       copyOfTrackParam.SetParameters(trackParam->GetSmoothParameters());
997       copyOfTrackParam.SetCovariances(trackParam->GetSmoothCovariances());
998       
999       // look for one second candidate in the same chamber
1000       for (Int_t hit = 0; hit < fNHitsForRecPerChamber[chamberId]; hit++) {
1001         
1002         hitForRec = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[chamberId]+hit);
1003         
1004         // look for a cluster in another detection element
1005         if (hitForRec->GetDetElemId() == detElemId) continue;
1006         
1007         // try to add the current hit fast
1008         if (!TryOneHitForRecFast(copyOfTrackParam, hitForRec)) continue;
1009         
1010         // try to add the current hit accurately
1011         chi2OfHitForRec = TryOneHitForRec(copyOfTrackParam, hitForRec, trackParamAtHit);
1012         
1013         // if better chi2 then prepare to add this cluster to the track
1014         if (chi2OfHitForRec < bestChi2OfHitForRec) {
1015           bestChi2OfHitForRec = chi2OfHitForRec;
1016           bestTrackParamAtHit = trackParamAtHit;
1017           foundOneHit = kTRUE;
1018         }
1019         
1020       }
1021       
1022       // add new cluster if any
1023       if (foundOneHit) {
1024         UpdateTrack(*track,bestTrackParamAtHit);
1025         bestTrackParamAtHit.SetAloneInChamber(kFALSE);
1026         trackParam->SetAloneInChamber(kFALSE);
1027         trackModified = kTRUE;
1028       }
1029       
1030       trackParam = nextTrackParam;
1031     }
1032     
1033     // re-fit track parameters if needed
1034     if (trackModified) Fit(*track, kTRUE, kTRUE);
1035     
1036     track = (AliMUONTrack*) fRecTracksPtr->After(track);
1037   }
1038   
1039 }
1040
1041   //__________________________________________________________________________
1042 void AliMUONTrackReconstructor::ImproveTracks()
1043 {
1044   /// Improve tracks by removing clusters with local chi2 highter than the defined cut
1045   /// Recompute track parameters and covariances at the remaining clusters
1046   AliDebug(1,"Enter ImproveTracks");
1047   
1048   Double_t localChi2, worstLocalChi2;
1049   Int_t worstChamber, previousChamber;
1050   AliMUONTrack *track, *nextTrack;
1051   AliMUONTrackParam *trackParamAtHit, *worstTrackParamAtHit, *previousTrackParam, *nextTrackParam;
1052   Double_t sigmaCut2 = AliMUONReconstructor::GetRecoParam()->GetSigmaCutForImprovement() *
1053                        AliMUONReconstructor::GetRecoParam()->GetSigmaCutForImprovement();
1054   
1055   // Remove double track to improve only "good" tracks
1056   RemoveDoubleTracks();
1057   
1058   track = (AliMUONTrack*) fRecTracksPtr->First();
1059   while (track) {
1060     
1061     // prepare next track in case the actual track is suppressed
1062     nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track);
1063     
1064     while (!track->IsImproved()) {
1065       
1066       // Update track parameters and covariances
1067       track->UpdateCovTrackParamAtHit();
1068       
1069       // Compute local chi2 of each hits
1070       track->ComputeLocalChi2(kTRUE);
1071       
1072       // Look for the hit to remove
1073       worstTrackParamAtHit = NULL;
1074       worstLocalChi2 = 0.;
1075       trackParamAtHit = (AliMUONTrackParam*) track->GetTrackParamAtHit()->First();
1076       while (trackParamAtHit) {
1077         
1078         // Pick up hit with the worst chi2
1079         localChi2 = trackParamAtHit->GetLocalChi2();
1080         if (localChi2 > worstLocalChi2) {
1081           worstLocalChi2 = localChi2;
1082           worstTrackParamAtHit = trackParamAtHit;
1083         }
1084         
1085       trackParamAtHit = (AliMUONTrackParam*) track->GetTrackParamAtHit()->After(trackParamAtHit);
1086       }
1087       
1088       // Check if bad hit found
1089       if (!worstTrackParamAtHit) {
1090         track->SetImproved(kTRUE);
1091         break;
1092       }
1093       
1094       // Check whether the worst chi2 is under requirement or not
1095       if (worstLocalChi2 < 2. * sigmaCut2) { // 2 because 2 quantities in chi2
1096         track->SetImproved(kTRUE);
1097         break;
1098       }
1099       
1100       // if the worst hit is not removable then remove the entire track
1101       if (!worstTrackParamAtHit->IsRemovable() && worstTrackParamAtHit->IsAloneInChamber()) {
1102         fRecTracksPtr->Remove(track);
1103         fNRecTracks--;
1104         break;
1105       }
1106       
1107       // Reset the second hit in the same station as being not removable
1108       // or reset the second hit in the same chamber as being alone
1109       worstChamber = worstTrackParamAtHit->GetHitForRecPtr()->GetChamberNumber();
1110       previousTrackParam = (AliMUONTrackParam*) track->GetTrackParamAtHit()->Before(worstTrackParamAtHit);
1111       nextTrackParam = (AliMUONTrackParam*) track->GetTrackParamAtHit()->After(worstTrackParamAtHit);
1112       if (worstTrackParamAtHit->IsAloneInChamber()) { // Worst hit removable and alone in chamber
1113         
1114         if (worstChamber%2 == 0) { // Modify flags in next chamber
1115           
1116           nextTrackParam->SetRemovable(kFALSE);
1117           if (!nextTrackParam->IsAloneInChamber()) // Make sure both hits in second chamber are not removable anymore
1118             ((AliMUONTrackParam*) track->GetTrackParamAtHit()->After(nextTrackParam))->SetRemovable(kFALSE);
1119           
1120         } else { // Modify flags in previous chamber
1121           
1122           previousTrackParam->SetRemovable(kFALSE);
1123           if (!previousTrackParam->IsAloneInChamber()) // Make sure both hits in second chamber are not removable anymore
1124             ((AliMUONTrackParam*) track->GetTrackParamAtHit()->Before(previousTrackParam))->SetRemovable(kFALSE);
1125           
1126         }
1127         
1128       } else { // Worst hit not alone in its chamber
1129         
1130         if (previousTrackParam) previousChamber = previousTrackParam->GetHitForRecPtr()->GetChamberNumber();
1131         else previousChamber = -1;
1132         
1133         if (previousChamber == worstChamber) { // the second hit on the same chamber is the previous one
1134           
1135           previousTrackParam->SetAloneInChamber(kTRUE);
1136           // transfert the removability to the second hit
1137           if (worstTrackParamAtHit->IsRemovable()) previousTrackParam->SetRemovable(kTRUE);
1138           
1139         } else { // the second hit on the same chamber is the next one
1140           
1141           nextTrackParam->SetAloneInChamber(kTRUE);
1142           // transfert the removability to the second hit
1143           if (worstTrackParamAtHit->IsRemovable()) nextTrackParam->SetRemovable(kTRUE);
1144           
1145         }
1146         
1147       }
1148       
1149       // Remove the worst hit
1150       track->RemoveTrackParamAtHit(worstTrackParamAtHit);
1151       
1152       // Re-fit the track:
1153       // Take into account the multiple scattering
1154       // Calculate the track parameter covariance matrix
1155       Fit(*track, kTRUE, kTRUE);
1156       
1157       // Printout for debuging
1158       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
1159         cout << "ImproveTracks: track " << fRecTracksPtr->IndexOf(track)+1 << " improved " << endl;
1160       }
1161       
1162     }
1163     
1164     track = nextTrack;
1165   }
1166   
1167   // compress the array in case of some tracks have been removed
1168   fRecTracksPtr->Compress();
1169   
1170 }
1171
1172   //__________________________________________________________________________
1173 void AliMUONTrackReconstructor::Finalize()
1174 {
1175   /// Fill AliMUONTrack's fHitForRecAtHit array
1176   /// Recompute track parameters and covariances at each attached cluster from those at the first one
1177   
1178   AliMUONTrack *track;
1179   AliMUONTrackParam *trackParamAtHit;
1180   
1181   track = (AliMUONTrack*) fRecTracksPtr->First();
1182   while (track) {
1183     
1184     // update track parameters if not already done
1185     if (!track->IsImproved()) track->UpdateCovTrackParamAtHit();
1186     
1187     trackParamAtHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First());
1188     while (trackParamAtHit) {
1189       
1190       // update array of track hit
1191       track->AddHitForRecAtHit(trackParamAtHit->GetHitForRecPtr());
1192       
1193       trackParamAtHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->After(trackParamAtHit));
1194     }
1195     
1196     track = (AliMUONTrack*) fRecTracksPtr->After(track);
1197     
1198   }
1199     
1200 }
1201
1202