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