- AliMUONRecoParam.cxx:
[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 "AliMUONVClusterServer.h"
34 #include "AliMUONVClusterStore.h"
35 #include "AliMUONTrack.h"
36 #include "AliMUONTrackParam.h"
37 #include "AliMUONTrackExtrap.h"
38
39 #include "AliMpArea.h"
40
41 #include "AliLog.h"
42
43 #include <TMinuit.h>
44 #include <Riostream.h>
45 #include <TMath.h>
46 #include <TMatrixD.h>
47
48 // Functions to be minimized with Minuit
49 void TrackChi2(Int_t &nParam, Double_t *gradient, Double_t &chi2, Double_t *param, Int_t flag);
50
51 /// \cond CLASSIMP
52 ClassImp(AliMUONTrackReconstructor) // Class implementation in ROOT context
53 /// \endcond
54
55   //__________________________________________________________________________
56 AliMUONTrackReconstructor::AliMUONTrackReconstructor(const AliMUONRecoParam* recoParam, AliMUONVClusterServer* clusterServer)
57   : AliMUONVTrackReconstructor(recoParam,clusterServer)
58 {
59   /// Constructor
60 }
61
62   //__________________________________________________________________________
63 AliMUONTrackReconstructor::~AliMUONTrackReconstructor()
64 {
65 /// Destructor
66
67
68   //__________________________________________________________________________
69 void AliMUONTrackReconstructor::MakeTrackCandidates(AliMUONVClusterStore& clusterStore)
70 {
71   /// To make track candidates (assuming linear propagation if the flag fgkMakeTrackCandidatesFast is set to kTRUE):
72   /// Start with segments station(1..) 4 or 5 then follow track in station 5 or 4.
73   /// Good candidates are made of at least three clusters.
74   /// Keep only best candidates or all of them according to the flag fgkTrackAllTracks.
75   
76   TClonesArray *segments;
77   AliMUONTrack *track;
78   Int_t iCandidate = 0;
79   Bool_t clusterFound;
80
81   AliDebug(1,"Enter MakeTrackCandidates");
82
83   // Unless we're doing combined tracking, we'll clusterize all stations at once
84   Int_t firstChamber(0);
85   Int_t lastChamber(9);
86   
87   if (GetRecoParam()->CombineClusterTrackReco()) {
88     // ... Here's the exception : ask the clustering to reconstruct
89     // clusters *only* in station 4 and 5 for combined tracking
90     firstChamber = 6;
91   }
92
93   for (Int_t i = firstChamber; i <= lastChamber; ++i ) 
94   {
95     if (fClusterServer && GetRecoParam()->UseChamber(i)) fClusterServer->Clusterize(i, clusterStore, AliMpArea());
96   }
97   
98   // Loop over stations(1..) 5 and 4 and make track candidates
99   for (Int_t istat=4; istat>=3; istat--) {
100     
101     // Make segments in the station
102     segments = MakeSegmentsBetweenChambers(clusterStore, 2*istat, 2*istat+1);
103     
104     // Loop over segments
105     for (Int_t iseg=0; iseg<segments->GetEntriesFast(); iseg++) 
106     {
107       AliDebug(1,Form("Making primary candidate(1..) %d",++iCandidate));
108       
109       // Transform segments to tracks and put them at the end of fRecTracksPtr
110       track = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack((AliMUONObjectPair*)((*segments)[iseg]),GetRecoParam()->GetBendingVertexDispersion());
111       fNRecTracks++;
112       
113       // Look for compatible cluster(s) in the other station
114       if (GetRecoParam()->MakeTrackCandidatesFast())
115         clusterFound = FollowLinearTrackInStation(*track, clusterStore, 7-istat);
116       else clusterFound = FollowTrackInStation(*track, clusterStore, 7-istat);
117       
118       // Remove track if no cluster found
119       if (!clusterFound && GetRecoParam()->RequestStation(7-istat)) {
120         fRecTracksPtr->Remove(track);
121         fNRecTracks--;
122       }
123       
124     }
125     
126     // delete the array of segments
127     delete segments;
128   }
129   
130   fRecTracksPtr->Compress(); // this is essential before checking tracks
131   
132   // Keep all different tracks or only the best ones as required
133   if (GetRecoParam()->TrackAllTracks()) RemoveIdenticalTracks();
134   else RemoveDoubleTracks();
135   
136   AliDebug(1,Form("Number of good candidates = %d",fNRecTracks));
137   
138 }
139
140   //__________________________________________________________________________
141 void AliMUONTrackReconstructor::MakeMoreTrackCandidates(AliMUONVClusterStore& clusterStore)
142 {
143   /// To make extra track candidates (assuming linear propagation if the flag fgkMakeTrackCandidatesFast is set to kTRUE):
144   /// clustering is supposed to be already done
145   /// Start with segments made of 1 cluster in each of the stations 4 and 5 then follow track in remaining chambers.
146   /// Good candidates are made of at least three clusters if both station are requested (two otherwise).
147   /// Keep only best candidates or all of them according to the flag fgkTrackAllTracks.
148   
149   TClonesArray *segments;
150   AliMUONTrack *track;
151   Int_t iCandidate = 0, iCurrentTrack, nCurrentTracks;
152   Bool_t clusterFound;
153   
154   AliDebug(1,"Enter MakeMoreTrackCandidates");
155   
156   // Double loop over chambers in stations(1..) 4 and 5 to make track candidates
157   for (Int_t ich1 = 6; ich1 <= 7; ich1++) {
158     for (Int_t ich2 = 8; ich2 <= 9; ich2++) {
159       
160       // Make segments in the station
161       segments = MakeSegmentsBetweenChambers(clusterStore, ich1, ich2);
162       
163       // Loop over segments
164       for (Int_t iseg=0; iseg<segments->GetEntriesFast(); iseg++) 
165       {
166         AliDebug(1,Form("Making primary candidate(1..) %d",++iCandidate));
167         
168         // Transform segments to tracks and put them at the end of fRecTracksPtr
169         iCurrentTrack = fRecTracksPtr->GetLast()+1;
170         track = new ((*fRecTracksPtr)[iCurrentTrack]) AliMUONTrack((AliMUONObjectPair*)((*segments)[iseg]),GetRecoParam()->GetBendingVertexDispersion());
171         fNRecTracks++;
172         
173         // Look for compatible cluster(s) in the second chamber of station 5
174         if (GetRecoParam()->MakeTrackCandidatesFast())
175           clusterFound = FollowLinearTrackInChamber(*track, clusterStore, 17-ich2);
176         else clusterFound = FollowTrackInChamber(*track, clusterStore, 17-ich2);
177         
178         // skip the original track in case it has been removed
179         if (GetRecoParam()->TrackAllTracks() && clusterFound) iCurrentTrack++;
180         
181         // loop over every new tracks
182         nCurrentTracks = fRecTracksPtr->GetLast()+1;
183         while (iCurrentTrack < nCurrentTracks) {
184           track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iCurrentTrack);
185           
186           // Look for compatible cluster(s) in the second chamber of station 4
187           if (GetRecoParam()->MakeTrackCandidatesFast())
188             FollowLinearTrackInChamber(*track, clusterStore, 13-ich1);
189           else FollowTrackInChamber(*track, clusterStore, 13-ich1);
190           
191           iCurrentTrack++;
192         }
193         
194       }
195       
196       // delete the array of segments
197       delete segments;
198     }
199   }
200   
201   fRecTracksPtr->Compress(); // this is essential before checking tracks
202   
203   // Keep all different tracks or only the best ones as required
204   if (GetRecoParam()->TrackAllTracks()) RemoveIdenticalTracks();
205   else RemoveDoubleTracks();
206   
207   AliDebug(1,Form("Number of good candidates = %d",fNRecTracks));
208   
209 }
210
211   //__________________________________________________________________________
212 void AliMUONTrackReconstructor::FollowTracks(AliMUONVClusterStore& clusterStore)
213 {
214   /// Follow tracks in stations(1..) 3, 2 and 1
215   AliDebug(1,"Enter FollowTracks");
216   
217   AliMUONTrack *track, *nextTrack;
218   AliMUONTrackParam *trackParam, *nextTrackParam;
219   Int_t currentNRecTracks;
220   
221   Double_t sigmaCut2 = GetRecoParam()->GetSigmaCutForTracking() *
222                        GetRecoParam()->GetSigmaCutForTracking();
223   
224   for (Int_t station = 2; station >= 0; station--) {
225     
226     // Save the actual number of reconstructed track in case of
227     // tracks are added or suppressed during the tracking procedure
228     // !! Do not compress fRecTracksPtr until the end of the loop over tracks !!
229     currentNRecTracks = fNRecTracks;
230     
231     for (Int_t iRecTrack = 0; iRecTrack <currentNRecTracks; iRecTrack++) {
232       AliDebug(1,Form("FollowTracks: track candidate(1..) %d", iRecTrack+1));
233       
234       track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iRecTrack);
235       
236       // Fit the track:
237       // Do not take into account the multiple scattering to speed up the fit
238       // Calculate the track parameter covariance matrix
239       // If there is no cluster out of station 4 or 5 then use the vertex to better constrain the fit
240       if (((AliMUONTrackParam*) track->GetTrackParamAtCluster()->First())->GetClusterPtr()->GetChamberId() > 5)
241         Fit(*track, kFALSE, kTRUE, kTRUE);
242       else Fit(*track, kFALSE, kFALSE, kTRUE);
243       
244       // remove track with absolute bending momentum out of limits
245       // or if the normalized chi2 is too high
246       Double_t bendingMomentum = TMath::Abs(1. / ((AliMUONTrackParam*)track->GetTrackParamAtCluster()->First())->GetInverseBendingMomentum());
247       if (bendingMomentum < GetRecoParam()->GetMinBendingMomentum() ||
248           bendingMomentum > GetRecoParam()->GetMaxBendingMomentum() ||
249           track->GetNormalizedChi2() > sigmaCut2) {
250         fRecTracksPtr->Remove(track);
251         fNRecTracks--;
252         continue;
253       }
254       
255       // save parameters from fit into smoothed parameters to complete track afterward
256       if (GetRecoParam()->ComplementTracks()) {
257         
258         if (station==2) { // save track parameters on stations 4 and 5
259           
260           // extrapolate track parameters and covariances at each cluster
261           track->UpdateCovTrackParamAtCluster();
262           
263           // save them
264           trackParam = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->First();
265           while (trackParam) {
266             trackParam->SetSmoothParameters(trackParam->GetParameters());
267             trackParam->SetSmoothCovariances(trackParam->GetCovariances());
268             trackParam = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->After(trackParam);
269           }
270           
271         } else { // or save track parameters on last station only
272           
273           trackParam = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->First();
274           if (trackParam->GetClusterPtr()->GetChamberId() < 2*(station+2)) {
275             
276             // save parameters from fit
277             trackParam->SetSmoothParameters(trackParam->GetParameters());
278             trackParam->SetSmoothCovariances(trackParam->GetCovariances());
279             
280             // save parameters extrapolated to the second chamber of the same station if it has been hit
281             nextTrackParam = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->After(trackParam);
282             if (nextTrackParam->GetClusterPtr()->GetChamberId() < 2*(station+2)) {
283               
284               // reset parameters and covariances
285               nextTrackParam->SetParameters(trackParam->GetParameters());
286               nextTrackParam->SetZ(trackParam->GetZ());
287               nextTrackParam->SetCovariances(trackParam->GetCovariances());
288               
289               // extrapolate them to the z of the corresponding cluster
290               AliMUONTrackExtrap::ExtrapToZCov(nextTrackParam, nextTrackParam->GetClusterPtr()->GetZ());
291               
292               // save them
293               nextTrackParam->SetSmoothParameters(nextTrackParam->GetParameters());
294               nextTrackParam->SetSmoothCovariances(nextTrackParam->GetCovariances());
295               
296             }
297             
298           }
299           
300         }
301         
302       }
303       
304       // Look for compatible cluster(s) in station(0..) "station"
305       if (!FollowTrackInStation(*track, clusterStore, station)) {
306         
307         // Try to recover track if required
308         if (GetRecoParam()->RecoverTracks()) {
309           
310           // work on a copy of the track if this station is not required
311           // to keep the case where no cluster is reconstructed as a possible candidate
312           if (!GetRecoParam()->RequestStation(station)) {
313             track = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(*track);
314             fNRecTracks++;
315           }
316           
317           // try to recover
318           if (!RecoverTrack(*track, clusterStore, station)) {
319             // remove track if no cluster found
320             fRecTracksPtr->Remove(track);
321             fNRecTracks--;
322           }
323           
324         } else if (GetRecoParam()->RequestStation(station)) {
325           // remove track if no cluster found
326           fRecTracksPtr->Remove(track);
327           fNRecTracks--;
328         } 
329         
330       }
331       
332     }
333     
334     // Compress fRecTracksPtr for the next step
335     fRecTracksPtr->Compress();
336     
337     // Keep only the best tracks if required
338     if (!GetRecoParam()->TrackAllTracks()) RemoveDoubleTracks();
339     
340   }
341   
342   // Last fit of track candidates with all stations
343   // Take into account the multiple scattering and remove bad tracks
344   Int_t trackIndex = -1;
345   track = (AliMUONTrack*) fRecTracksPtr->First();
346   while (track) {
347     
348     trackIndex++;
349     nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track); // prepare next track
350     
351     Fit(*track, kTRUE, kFALSE, kTRUE);
352     
353     // Printout for debuging
354     if (AliLog::GetGlobalDebugLevel() >= 3) {
355       cout << "FollowTracks: track candidate(0..) " << trackIndex << " after final fit" << endl;
356       track->RecursiveDump();
357     } 
358     
359     // Remove the track if the normalized chi2 is too high
360     if (track->GetNormalizedChi2() > sigmaCut2) {
361       fRecTracksPtr->Remove(track);
362       fNRecTracks--;
363       track = nextTrack;
364       continue;
365     }
366     
367     // save parameters from fit into smoothed parameters to complete track afterward
368     if (GetRecoParam()->ComplementTracks()) {
369       
370       trackParam = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->First();
371       if (trackParam->GetClusterPtr()->GetChamberId() < 2) {
372         
373         // save parameters from fit
374         trackParam->SetSmoothParameters(trackParam->GetParameters());
375         trackParam->SetSmoothCovariances(trackParam->GetCovariances());
376         
377         // save parameters extrapolated to the second chamber of the same station if it has been hit
378         nextTrackParam = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->After(trackParam);
379         if (nextTrackParam->GetClusterPtr()->GetChamberId() < 2) {
380           
381           // reset parameters and covariances
382           nextTrackParam->SetParameters(trackParam->GetParameters());
383           nextTrackParam->SetZ(trackParam->GetZ());
384           nextTrackParam->SetCovariances(trackParam->GetCovariances());
385           
386           // extrapolate them to the z of the corresponding cluster
387           AliMUONTrackExtrap::ExtrapToZCov(nextTrackParam, nextTrackParam->GetClusterPtr()->GetZ());
388           
389           // save them
390           nextTrackParam->SetSmoothParameters(nextTrackParam->GetParameters());
391           nextTrackParam->SetSmoothCovariances(nextTrackParam->GetCovariances());
392           
393         }
394         
395       }
396       
397     }
398     
399     track = nextTrack;
400     
401   }
402   
403   fRecTracksPtr->Compress();
404   
405 }
406
407   //__________________________________________________________________________
408 Bool_t AliMUONTrackReconstructor::FollowTrackInChamber(AliMUONTrack &trackCandidate, AliMUONVClusterStore& clusterStore, Int_t nextChamber)
409 {
410   /// Follow trackCandidate in chamber(0..) nextChamber and search for compatible cluster(s)
411   /// Keep all possibilities or only the best one(s) according to the flag fgkTrackAllTracks:
412   /// kTRUE:  duplicate "trackCandidate" if there are several possibilities and add the new tracks at the end of
413   ///         fRecTracksPtr to avoid conficts with other track candidates at this current stage of the tracking procedure.
414   ///         Remove the obsolete "trackCandidate" at the end.
415   /// kFALSE: add only the best cluster(s) to the "trackCandidate". Try to add a couple of clusters in priority.
416   AliDebug(1,Form("Enter FollowTrackInChamber(1..) %d", nextChamber+1));
417   
418   Double_t chi2WithOneCluster = 1.e10;
419   Double_t maxChi2WithOneCluster = 2. * GetRecoParam()->GetSigmaCutForTracking() *
420                                         GetRecoParam()->GetSigmaCutForTracking(); // 2 because 2 quantities in chi2
421   Double_t bestChi2WithOneCluster = maxChi2WithOneCluster;
422   Bool_t foundOneCluster = kFALSE;
423   AliMUONTrack *newTrack = 0x0;
424   AliMUONVCluster *cluster;
425   AliMUONTrackParam extrapTrackParam;
426   AliMUONTrackParam extrapTrackParamAtCh;
427   AliMUONTrackParam extrapTrackParamAtCluster;
428   AliMUONTrackParam bestTrackParamAtCluster;
429   
430   // Get track parameters according to the propagation direction
431   if (nextChamber > 7) extrapTrackParamAtCh = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->Last();
432   else extrapTrackParamAtCh = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->First();
433   
434   // Printout for debuging
435   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
436     cout<<endl<<"Track parameters and covariances at first cluster:"<<endl;
437     extrapTrackParamAtCh.GetParameters().Print();
438     extrapTrackParamAtCh.GetCovariances().Print();
439   }
440   
441   // Add MCS effect
442   AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(),1.);
443   
444   // Add MCS in the missing chamber(s) if any
445   Int_t currentChamber = extrapTrackParamAtCh.GetClusterPtr()->GetChamberId();
446   while (currentChamber > nextChamber + 1) {
447     // extrapolation to the missing chamber
448     currentChamber--;
449     AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(currentChamber));
450     // add MCS effect
451     AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(),1.);
452   }
453   
454   //Extrapolate trackCandidate to chamber
455   AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(nextChamber));
456   
457   // Printout for debuging
458   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
459     cout<<endl<<"Track parameters and covariances at first cluster extrapolated to z = "<<AliMUONConstants::DefaultChamberZ(nextChamber)<<":"<<endl;
460     extrapTrackParamAtCh.GetParameters().Print();
461     extrapTrackParamAtCh.GetCovariances().Print();
462   }
463   
464   // Printout for debuging
465   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
466     cout << "FollowTrackInStation: look for clusters in chamber(1..): " << nextChamber+1 << endl;
467   }
468   
469   // Ask the clustering to reconstruct new clusters around the track position in the current chamber
470   // except for station 4 and 5 that are already entirely clusterized
471   if (GetRecoParam()->CombineClusterTrackReco()) {
472     if (nextChamber < 6) AskForNewClustersInChamber(extrapTrackParamAtCh, clusterStore, nextChamber);
473   }
474   
475   // Create iterators to loop over clusters in both chambers
476   TIter next(clusterStore.CreateChamberIterator(nextChamber,nextChamber));
477   
478   // look for cluster in chamber
479   while ( ( cluster = static_cast<AliMUONVCluster*>(next()) ) ) {
480     
481     // try to add the current cluster fast
482     if (!TryOneClusterFast(extrapTrackParamAtCh, cluster)) continue;
483     
484     // try to add the current cluster accuratly
485     chi2WithOneCluster = TryOneCluster(extrapTrackParamAtCh, cluster, extrapTrackParamAtCluster);
486     
487     // if good chi2 then consider to add cluster
488     if (chi2WithOneCluster < maxChi2WithOneCluster) {
489       foundOneCluster = kTRUE;
490       
491       // Printout for debuging
492       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
493         cout << "FollowTrackInStation: found one cluster in chamber(1..): " << nextChamber+1
494         << " (Chi2 = " << chi2WithOneCluster << ")" << endl;
495         cluster->Print();
496       }
497       
498       if (GetRecoParam()->TrackAllTracks()) {
499         // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new cluster
500         newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
501         UpdateTrack(*newTrack,extrapTrackParamAtCluster);
502         fNRecTracks++;
503         
504         // Printout for debuging
505         if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
506           cout << "FollowTrackInStation: added one cluster in chamber(1..): " << nextChamber+1 << endl;
507           if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
508         }
509         
510       } else if (chi2WithOneCluster < bestChi2WithOneCluster) {
511         // keep track of the best single cluster except if a couple of clusters has already been found
512         bestChi2WithOneCluster = chi2WithOneCluster;
513         bestTrackParamAtCluster = extrapTrackParamAtCluster;
514       }
515       
516     }
517     
518   }
519   
520   // fill out the best track if required else clean up the fRecTracksPtr array
521   if (!GetRecoParam()->TrackAllTracks()) {
522     if (foundOneCluster) {
523       UpdateTrack(trackCandidate,bestTrackParamAtCluster);
524       
525       // Printout for debuging
526       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
527         cout << "FollowTrackInStation: added the best cluster in chamber(1..): " << bestTrackParamAtCluster.GetClusterPtr()->GetChamberId()+1 << endl;
528         if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
529       }
530       
531     } else return kFALSE;
532     
533   } else if (foundOneCluster) {
534     
535     // remove obsolete track
536     fRecTracksPtr->Remove(&trackCandidate);
537     fNRecTracks--;
538     
539   } else return kFALSE;
540   
541   return kTRUE;
542   
543 }
544
545   //__________________________________________________________________________
546 Bool_t AliMUONTrackReconstructor::FollowTrackInStation(AliMUONTrack &trackCandidate, AliMUONVClusterStore& clusterStore, Int_t nextStation)
547 {
548   /// Follow trackCandidate in station(0..) nextStation and search for compatible cluster(s)
549   /// Keep all possibilities or only the best one(s) according to the flag fgkTrackAllTracks:
550   /// kTRUE:  duplicate "trackCandidate" if there are several possibilities and add the new tracks at the end of
551   ///         fRecTracksPtr to avoid conficts with other track candidates at this current stage of the tracking procedure.
552   ///         Remove the obsolete "trackCandidate" at the end.
553   /// kFALSE: add only the best cluster(s) to the "trackCandidate". Try to add a couple of clusters in priority.
554   AliDebug(1,Form("Enter FollowTrackInStation(1..) %d", nextStation+1));
555   
556   // Order the chamber according to the propagation direction (tracking starts with chamber 2):
557   // - nextStation == station(1...) 5 => forward propagation
558   // - nextStation < station(1...) 5 => backward propagation
559   Int_t ch1, ch2;
560   if (nextStation==4) {
561     ch1 = 2*nextStation+1;
562     ch2 = 2*nextStation;
563   } else {
564     ch1 = 2*nextStation;
565     ch2 = 2*nextStation+1;
566   }
567   
568   Double_t chi2WithOneCluster = 1.e10;
569   Double_t chi2WithTwoClusters = 1.e10;
570   Double_t maxChi2WithOneCluster = 2. * GetRecoParam()->GetSigmaCutForTracking() *
571                                         GetRecoParam()->GetSigmaCutForTracking(); // 2 because 2 quantities in chi2
572   Double_t maxChi2WithTwoClusters = 4. * GetRecoParam()->GetSigmaCutForTracking() *
573                                          GetRecoParam()->GetSigmaCutForTracking(); // 4 because 4 quantities in chi2
574   Double_t bestChi2WithOneCluster = maxChi2WithOneCluster;
575   Double_t bestChi2WithTwoClusters = maxChi2WithTwoClusters;
576   Bool_t foundOneCluster = kFALSE;
577   Bool_t foundTwoClusters = kFALSE;
578   AliMUONTrack *newTrack = 0x0;
579   AliMUONVCluster *clusterCh1, *clusterCh2;
580   AliMUONTrackParam extrapTrackParam;
581   AliMUONTrackParam extrapTrackParamAtCh;
582   AliMUONTrackParam extrapTrackParamAtCluster1;
583   AliMUONTrackParam extrapTrackParamAtCluster2;
584   AliMUONTrackParam bestTrackParamAtCluster1;
585   AliMUONTrackParam bestTrackParamAtCluster2;
586   
587   Int_t nClusters = clusterStore.GetSize();
588   Bool_t *clusterCh1Used = new Bool_t[nClusters];
589   for (Int_t i = 0; i < nClusters; i++) clusterCh1Used[i] = kFALSE;
590   Int_t iCluster1;
591   
592   // Get track parameters according to the propagation direction
593   if (nextStation==4) extrapTrackParamAtCh = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->Last();
594   else extrapTrackParamAtCh = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->First();
595   
596   // Printout for debuging
597   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
598     cout<<endl<<"Track parameters and covariances at first cluster:"<<endl;
599     extrapTrackParamAtCh.GetParameters().Print();
600     extrapTrackParamAtCh.GetCovariances().Print();
601   }
602   
603   // Add MCS effect
604   AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(),1.);
605   
606   // Add MCS in the missing chamber(s) if any
607   Int_t currentChamber = extrapTrackParamAtCh.GetClusterPtr()->GetChamberId();
608   while (ch1 < ch2 && currentChamber > ch2 + 1) {
609     // extrapolation to the missing chamber
610     currentChamber--;
611     AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(currentChamber));
612     // add MCS effect
613     AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(),1.);
614   }
615   
616   //Extrapolate trackCandidate to chamber "ch2"
617   AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(ch2));
618   
619   // Printout for debuging
620   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
621     cout<<endl<<"Track parameters and covariances at first cluster extrapolated to z = "<<AliMUONConstants::DefaultChamberZ(ch2)<<":"<<endl;
622     extrapTrackParamAtCh.GetParameters().Print();
623     extrapTrackParamAtCh.GetCovariances().Print();
624   }
625   
626   // Printout for debuging
627   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
628     cout << "FollowTrackInStation: look for clusters in chamber(1..): " << ch2+1 << endl;
629   }
630   
631   // Ask the clustering to reconstruct new clusters around the track position in the current station
632   // except for station 4 and 5 that are already entirely clusterized
633   if (GetRecoParam()->CombineClusterTrackReco()) {
634     if (nextStation < 3) AskForNewClustersInStation(extrapTrackParamAtCh, clusterStore, nextStation);
635   }
636   
637   // Create iterators to loop over clusters in both chambers
638   TIter nextInCh1(clusterStore.CreateChamberIterator(ch1,ch1));
639   TIter nextInCh2(clusterStore.CreateChamberIterator(ch2,ch2));
640   
641   // look for candidates in chamber 2
642   while ( ( clusterCh2 = static_cast<AliMUONVCluster*>(nextInCh2()) ) ) {
643     
644     // try to add the current cluster fast
645     if (!TryOneClusterFast(extrapTrackParamAtCh, clusterCh2)) continue;
646     
647     // try to add the current cluster accuratly
648     chi2WithOneCluster = TryOneCluster(extrapTrackParamAtCh, clusterCh2, extrapTrackParamAtCluster2);
649     
650     // if good chi2 then try to attach a cluster in the other chamber too
651     if (chi2WithOneCluster < maxChi2WithOneCluster) {
652       Bool_t foundSecondCluster = kFALSE;
653       
654       // Printout for debuging
655       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
656         cout << "FollowTrackInStation: found one cluster in chamber(1..): " << ch2+1
657              << " (Chi2 = " << chi2WithOneCluster << ")" << endl;
658         clusterCh2->Print();
659         cout << "                      look for second clusters in chamber(1..): " << ch1+1 << " ..." << endl;
660       }
661       
662       // add MCS effect for next step
663       AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCluster2,AliMUONConstants::ChamberThicknessInX0(),1.);
664       
665       // copy new track parameters for next step
666       extrapTrackParam = extrapTrackParamAtCluster2;
667       
668       //Extrapolate track parameters to chamber "ch1"
669       AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParam, AliMUONConstants::DefaultChamberZ(ch1));
670       
671       // reset cluster iterator of chamber 1
672       nextInCh1.Reset();
673       iCluster1 = -1;
674       
675       // look for second candidates in chamber 1
676       while ( ( clusterCh1 = static_cast<AliMUONVCluster*>(nextInCh1()) ) ) {
677         iCluster1++;
678         
679         // try to add the current cluster fast
680         if (!TryOneClusterFast(extrapTrackParam, clusterCh1)) continue;
681         
682         // try to add the current cluster accuratly
683         chi2WithTwoClusters = TryTwoClusters(extrapTrackParamAtCluster2, clusterCh1, extrapTrackParamAtCluster1);
684         
685         // if good chi2 then create a new track by adding the 2 clusters to the "trackCandidate"
686         if (chi2WithTwoClusters < maxChi2WithTwoClusters) {
687           foundSecondCluster = kTRUE;
688           foundTwoClusters = kTRUE;
689           
690           // Printout for debuging
691           if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
692             cout << "FollowTrackInStation: found second cluster in chamber(1..): " << ch1+1
693                  << " (Global Chi2 = " << chi2WithTwoClusters << ")" << endl;
694             clusterCh1->Print();
695           }
696           
697           if (GetRecoParam()->TrackAllTracks()) {
698             // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new clusters
699             newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
700             UpdateTrack(*newTrack,extrapTrackParamAtCluster1,extrapTrackParamAtCluster2);
701             fNRecTracks++;
702             
703             // Tag clusterCh1 as used
704             clusterCh1Used[iCluster1] = kTRUE;
705             
706             // Printout for debuging
707             if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
708               cout << "FollowTrackInStation: added two clusters in station(1..): " << nextStation+1 << endl;
709               if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
710             }
711             
712           } else if (chi2WithTwoClusters < bestChi2WithTwoClusters) {
713             // keep track of the best couple of clusters
714             bestChi2WithTwoClusters = chi2WithTwoClusters;
715             bestTrackParamAtCluster1 = extrapTrackParamAtCluster1;
716             bestTrackParamAtCluster2 = extrapTrackParamAtCluster2;
717           }
718           
719         }
720         
721       }
722       
723       // if no clusterCh1 found then consider to add clusterCh2 only
724       if (!foundSecondCluster) {
725         foundOneCluster = kTRUE;
726         
727         if (GetRecoParam()->TrackAllTracks()) {
728           // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new cluster
729           newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
730           UpdateTrack(*newTrack,extrapTrackParamAtCluster2);
731           fNRecTracks++;
732           
733           // Printout for debuging
734           if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
735             cout << "FollowTrackInStation: added one cluster in chamber(1..): " << ch2+1 << endl;
736             if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
737           }
738           
739         } else if (!foundTwoClusters && chi2WithOneCluster < bestChi2WithOneCluster) {
740           // keep track of the best single cluster except if a couple of clusters has already been found
741           bestChi2WithOneCluster = chi2WithOneCluster;
742           bestTrackParamAtCluster1 = extrapTrackParamAtCluster2;
743         }
744         
745       }
746       
747     }
748     
749   }
750   
751   // look for candidates in chamber 1 not already attached to a track
752   // if we want to keep all possible tracks or if no good couple of clusters has been found
753   if (GetRecoParam()->TrackAllTracks() || !foundTwoClusters) {
754     
755     // add MCS effect for next step
756     AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(),1.);
757     
758     //Extrapolate trackCandidate to chamber "ch1"
759     AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(ch1));
760     
761     // Printout for debuging
762     if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
763       cout<<endl<<"Track parameters and covariances at first cluster extrapolated to z = "<<AliMUONConstants::DefaultChamberZ(ch1)<<":"<<endl;
764       extrapTrackParamAtCh.GetParameters().Print();
765       extrapTrackParamAtCh.GetCovariances().Print();
766     }
767     
768     // Printout for debuging
769     if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
770       cout << "FollowTrackInStation: look for single clusters in chamber(1..): " << ch1+1 << endl;
771     }
772     
773     // reset cluster iterator of chamber 1
774     nextInCh1.Reset();
775     iCluster1 = -1;
776     
777     // look for second candidates in chamber 1
778     while ( ( clusterCh1 = static_cast<AliMUONVCluster*>(nextInCh1()) ) ) {
779       iCluster1++;
780       
781       if (clusterCh1Used[iCluster1]) continue; // Skip cluster already used
782       
783       // try to add the current cluster fast
784       if (!TryOneClusterFast(extrapTrackParamAtCh, clusterCh1)) continue;
785       
786       // try to add the current cluster accuratly
787       chi2WithOneCluster = TryOneCluster(extrapTrackParamAtCh, clusterCh1, extrapTrackParamAtCluster1);
788       
789       // if good chi2 then consider to add clusterCh1
790       // We do not try to attach a cluster in the other chamber too since it has already been done above
791       if (chi2WithOneCluster < maxChi2WithOneCluster) {
792         foundOneCluster = kTRUE;
793         
794         // Printout for debuging
795         if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
796           cout << "FollowTrackInStation: found one cluster in chamber(1..): " << ch1+1
797                << " (Chi2 = " << chi2WithOneCluster << ")" << endl;
798           clusterCh1->Print();
799         }
800         
801         if (GetRecoParam()->TrackAllTracks()) {
802           // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new cluster
803           newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
804           UpdateTrack(*newTrack,extrapTrackParamAtCluster1);
805           fNRecTracks++;
806           
807           // Printout for debuging
808           if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
809             cout << "FollowTrackInStation: added one cluster in chamber(1..): " << ch1+1 << endl;
810             if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
811           }
812           
813         } else if (chi2WithOneCluster < bestChi2WithOneCluster) {
814           // keep track of the best single cluster except if a couple of clusters has already been found
815           bestChi2WithOneCluster = chi2WithOneCluster;
816           bestTrackParamAtCluster1 = extrapTrackParamAtCluster1;
817         }
818         
819       }
820       
821     }
822     
823   }
824   
825   // fill out the best track if required else clean up the fRecTracksPtr array
826   if (!GetRecoParam()->TrackAllTracks()) {
827     if (foundTwoClusters) {
828       UpdateTrack(trackCandidate,bestTrackParamAtCluster1,bestTrackParamAtCluster2);
829       
830       // Printout for debuging
831       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
832         cout << "FollowTrackInStation: added the two best clusters in station(1..): " << nextStation+1 << endl;
833         if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
834       }
835       
836     } else if (foundOneCluster) {
837       UpdateTrack(trackCandidate,bestTrackParamAtCluster1);
838       
839       // Printout for debuging
840       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
841         cout << "FollowTrackInStation: added the best cluster in chamber(1..): " << bestTrackParamAtCluster1.GetClusterPtr()->GetChamberId()+1 << endl;
842         if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
843       }
844       
845     } else {
846       delete [] clusterCh1Used;
847       return kFALSE;
848     }
849     
850   } else if (foundOneCluster || foundTwoClusters) {
851     
852     // remove obsolete track
853     fRecTracksPtr->Remove(&trackCandidate);
854     fNRecTracks--;
855     
856   } else {
857     delete [] clusterCh1Used;
858     return kFALSE;
859   }
860   
861   delete [] clusterCh1Used;
862   return kTRUE;
863   
864 }
865
866   //__________________________________________________________________________
867 Double_t AliMUONTrackReconstructor::TryTwoClusters(const AliMUONTrackParam &trackParamAtCluster1, AliMUONVCluster* cluster2,
868                                                    AliMUONTrackParam &trackParamAtCluster2)
869 {
870 /// Test the compatibility between the track and the 2 clusters together (using trackParam's covariance matrix):
871 /// return the corresponding Chi2 accounting for covariances between the 2 clusters
872 /// return trackParamAtCluster1 & 2
873   
874   // extrapolate track parameters at the z position of the second cluster (no need to extrapolate the covariances)
875   trackParamAtCluster2.SetParameters(trackParamAtCluster1.GetParameters());
876   trackParamAtCluster2.SetZ(trackParamAtCluster1.GetZ());
877   AliMUONTrackExtrap::ExtrapToZ(&trackParamAtCluster2, cluster2->GetZ());
878   
879   // set pointer to cluster2 into trackParamAtCluster2
880   trackParamAtCluster2.SetClusterPtr(cluster2);
881   
882   // Set differences between track and the 2 clusters in the bending and non bending directions
883   AliMUONVCluster* cluster1 = trackParamAtCluster1.GetClusterPtr();
884   TMatrixD dPos(4,1);
885   dPos(0,0) = cluster1->GetX() - trackParamAtCluster1.GetNonBendingCoor();
886   dPos(1,0) = cluster1->GetY() - trackParamAtCluster1.GetBendingCoor();
887   dPos(2,0) = cluster2->GetX() - trackParamAtCluster2.GetNonBendingCoor();
888   dPos(3,0) = cluster2->GetY() - trackParamAtCluster2.GetBendingCoor();
889   
890   // Calculate the error matrix from the track parameter covariances at first cluster
891   TMatrixD error(4,4);
892   error.Zero();
893   if (trackParamAtCluster1.CovariancesExist()) {
894     // Save track parameters at first cluster
895     TMatrixD paramAtCluster1Save(trackParamAtCluster1.GetParameters());
896     
897     // Save track coordinates at second cluster
898     Double_t nonBendingCoor2 = trackParamAtCluster2.GetNonBendingCoor();
899     Double_t bendingCoor2    = trackParamAtCluster2.GetBendingCoor();
900     
901     // copy track parameters at first cluster for jacobian calculation
902     AliMUONTrackParam trackParam(trackParamAtCluster1);
903     
904     // add MCS effect to the covariance matrix at first cluster
905     AliMUONTrackExtrap::AddMCSEffect(&trackParam,AliMUONConstants::ChamberThicknessInX0(),1.);
906     
907     // Get the pointer to the parameter covariance matrix at first cluster
908     const TMatrixD& kParamCov = trackParam.GetCovariances();
909     
910     // Calculate the jacobian related to the transformation between track parameters
911     // at first cluster and track coordinates at the 2 cluster z-positions
912     TMatrixD jacob(4,5);
913     jacob.Zero();
914     // first derivative at the first cluster:
915     jacob(0,0) = 1.; // dx1/dx
916     jacob(1,2) = 1.; // dy1/dy
917     // first derivative at the second cluster:
918     TMatrixD dParam(5,1);
919     for (Int_t i=0; i<5; i++) {
920       // Skip jacobian calculation for parameters with no associated error
921       if (kParamCov(i,i) == 0.) continue;
922       // Small variation of parameter i only
923       for (Int_t j=0; j<5; j++) {
924         if (j==i) {
925           dParam(j,0) = TMath::Sqrt(kParamCov(i,i));
926           if (j == 4) dParam(j,0) *= TMath::Sign(1.,-paramAtCluster1Save(4,0)); // variation always in the same direction
927         } else dParam(j,0) = 0.;
928       }
929       
930       // Set new track parameters at first cluster
931       trackParam.SetParameters(paramAtCluster1Save);
932       trackParam.AddParameters(dParam);
933       trackParam.SetZ(cluster1->GetZ());
934       
935       // Extrapolate new track parameters to the z position of the second cluster
936       AliMUONTrackExtrap::ExtrapToZ(&trackParam,cluster2->GetZ());
937       
938       // Calculate the jacobian
939       jacob(2,i) = (trackParam.GetNonBendingCoor() - nonBendingCoor2) / dParam(i,0); // dx2/dParami
940       jacob(3,i) = (trackParam.GetBendingCoor()    - bendingCoor2   ) / dParam(i,0); // dy2/dParami
941     }
942     
943     // Calculate the error matrix
944     TMatrixD tmp(jacob,TMatrixD::kMult,kParamCov);
945     error = TMatrixD(tmp,TMatrixD::kMultTranspose,jacob);
946   }
947   
948   // Add cluster resolution to the error matrix
949   error(0,0) += cluster1->GetErrX2();
950   error(1,1) += cluster1->GetErrY2();
951   error(2,2) += cluster2->GetErrX2();
952   error(3,3) += cluster2->GetErrY2();
953   
954   // invert the error matrix for Chi2 calculation
955   if (error.Determinant() != 0) {
956     error.Invert();
957   } else {
958     AliWarning(" Determinant error=0");
959     return 1.e10;
960   }
961   
962   // Compute the Chi2 value
963   TMatrixD tmp2(dPos,TMatrixD::kTransposeMult,error);
964   TMatrixD result(tmp2,TMatrixD::kMult,dPos);
965   
966   return result(0,0);
967   
968 }
969
970   //__________________________________________________________________________
971 void AliMUONTrackReconstructor::UpdateTrack(AliMUONTrack &track, AliMUONTrackParam &trackParamAtCluster)
972 {
973   /// Add 1 cluster to the track candidate
974   /// Update chi2 of the track 
975   
976   // Compute local chi2
977   AliMUONVCluster* cluster = trackParamAtCluster.GetClusterPtr();
978   Double_t deltaX = trackParamAtCluster.GetNonBendingCoor() - cluster->GetX();
979   Double_t deltaY = trackParamAtCluster.GetBendingCoor() - cluster->GetY();
980   Double_t localChi2 = deltaX*deltaX / cluster->GetErrX2() +
981                        deltaY*deltaY / cluster->GetErrY2();
982   
983   // Flag cluster as being not removable
984   if (GetRecoParam()->RequestStation(cluster->GetChamberId()/2))
985     trackParamAtCluster.SetRemovable(kFALSE);
986   else trackParamAtCluster.SetRemovable(kTRUE);
987   trackParamAtCluster.SetLocalChi2(0.); // --> Local chi2 not used
988   
989   // Update the chi2 of the new track
990   track.SetGlobalChi2(track.GetGlobalChi2() + localChi2);
991   
992   // Update TrackParamAtCluster
993   track.AddTrackParamAtCluster(trackParamAtCluster,*cluster);
994   
995 }
996
997   //__________________________________________________________________________
998 void AliMUONTrackReconstructor::UpdateTrack(AliMUONTrack &track, AliMUONTrackParam &trackParamAtCluster1, AliMUONTrackParam &trackParamAtCluster2)
999 {
1000   /// Add 2 clusters to the track candidate
1001   /// Update track and local chi2
1002   
1003   // Update local chi2 at first cluster
1004   AliMUONVCluster* cluster1 = trackParamAtCluster1.GetClusterPtr();
1005   Double_t deltaX = trackParamAtCluster1.GetNonBendingCoor() - cluster1->GetX();
1006   Double_t deltaY = trackParamAtCluster1.GetBendingCoor() - cluster1->GetY();
1007   Double_t localChi2AtCluster1 = deltaX*deltaX / cluster1->GetErrX2() +
1008                                  deltaY*deltaY / cluster1->GetErrY2();
1009   trackParamAtCluster1.SetLocalChi2(localChi2AtCluster1);
1010   
1011   // Flag first cluster as being removable
1012   trackParamAtCluster1.SetRemovable(kTRUE);
1013   
1014   // Update local chi2 at second cluster
1015   AliMUONVCluster* cluster2 = trackParamAtCluster2.GetClusterPtr();
1016   deltaX = trackParamAtCluster2.GetNonBendingCoor() - cluster2->GetX();
1017   deltaY = trackParamAtCluster2.GetBendingCoor() - cluster2->GetY();
1018   Double_t localChi2AtCluster2 = deltaX*deltaX / cluster2->GetErrX2() +
1019                                  deltaY*deltaY / cluster2->GetErrY2();
1020   trackParamAtCluster2.SetLocalChi2(localChi2AtCluster2);
1021   
1022   // Flag first cluster as being removable
1023   trackParamAtCluster2.SetRemovable(kTRUE);
1024   
1025   // Update the chi2 of the new track
1026   track.SetGlobalChi2(track.GetGlobalChi2() + localChi2AtCluster1 + localChi2AtCluster2);
1027   
1028   // Update TrackParamAtCluster
1029   track.AddTrackParamAtCluster(trackParamAtCluster1,*cluster1);
1030   track.AddTrackParamAtCluster(trackParamAtCluster2,*cluster2);
1031   
1032 }
1033
1034   //__________________________________________________________________________
1035 Bool_t AliMUONTrackReconstructor::RecoverTrack(AliMUONTrack &trackCandidate, AliMUONVClusterStore& clusterStore, Int_t nextStation)
1036 {
1037   /// Try to recover the track candidate in the next station
1038   /// by removing the worst of the two clusters attached in the current station
1039   /// Return kTRUE if recovering succeeds
1040   AliDebug(1,"Enter RecoverTrack");
1041   
1042   // Do not try to recover track until we have attached cluster(s) on station(1..) 3
1043   if (nextStation > 1) return kFALSE;
1044   
1045   Int_t worstClusterNumber = -1;
1046   Double_t localChi2, worstLocalChi2 = -1.;
1047   
1048   // Look for the cluster to remove
1049   for (Int_t clusterNumber = 0; clusterNumber < 2; clusterNumber++) {
1050     AliMUONTrackParam *trackParamAtCluster = (AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->UncheckedAt(clusterNumber);
1051     
1052     // check if current cluster is in the previous station
1053     if (trackParamAtCluster->GetClusterPtr()->GetChamberId()/2 != nextStation+1) break;
1054     
1055     // check if current cluster is removable
1056     if (!trackParamAtCluster->IsRemovable()) return kFALSE;
1057     
1058     // reset the current cluster as beig not removable if it is on a required station
1059     if (GetRecoParam()->RequestStation(nextStation+1)) trackParamAtCluster->SetRemovable(kFALSE);
1060     
1061     // Pick up cluster with the worst chi2
1062     localChi2 = trackParamAtCluster->GetLocalChi2();
1063     if (localChi2 > worstLocalChi2) {
1064       worstLocalChi2 = localChi2;
1065       worstClusterNumber = clusterNumber;
1066     }
1067     
1068   }
1069   
1070   // check if worst cluster found
1071   if (worstClusterNumber < 0) return kFALSE;
1072   
1073   // Remove the worst cluster
1074   trackCandidate.RemoveTrackParamAtCluster((AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->UncheckedAt(worstClusterNumber));
1075   
1076   // Re-fit the track:
1077   // Do not take into account the multiple scattering to speed up the fit
1078   // Calculate the track parameter covariance matrix
1079   Fit(trackCandidate, kFALSE, kFALSE, kTRUE);
1080   
1081   // Look for new cluster(s) in next station
1082   return FollowTrackInStation(trackCandidate,clusterStore,nextStation);
1083   
1084 }
1085
1086   //__________________________________________________________________________
1087 void AliMUONTrackReconstructor::SetVertexErrXY2ForFit(AliMUONTrack &trackCandidate)
1088 {
1089   /// Compute the vertex resolution square from natural vertex dispersion and
1090   /// multiple scattering effets according to trackCandidate path in absorber
1091   /// It is necessary to account for multiple scattering effects here instead of during the fit of
1092   /// the "trackCandidate" to do not influence the result by changing track resolution at vertex
1093   AliDebug(1,"Enter SetVertexForFit");
1094   
1095   Double_t nonBendingReso2 = GetRecoParam()->GetNonBendingVertexDispersion() *
1096                              GetRecoParam()->GetNonBendingVertexDispersion();
1097   Double_t bendingReso2 = GetRecoParam()->GetBendingVertexDispersion() *
1098                           GetRecoParam()->GetBendingVertexDispersion();
1099   
1100   // add multiple scattering effets
1101   AliMUONTrackParam paramAtVertex(*((AliMUONTrackParam*)(trackCandidate.GetTrackParamAtCluster()->First())));
1102   paramAtVertex.DeleteCovariances(); // to be sure to account only for multiple scattering
1103   AliMUONTrackExtrap::ExtrapToVertexUncorrected(&paramAtVertex,0.);
1104   const TMatrixD& kParamCov = paramAtVertex.GetCovariances();
1105   nonBendingReso2 += kParamCov(0,0);
1106   bendingReso2 += kParamCov(2,2);
1107   
1108   // Set the vertex resolution square
1109   trackCandidate.SetVertexErrXY2(nonBendingReso2,bendingReso2);
1110 }
1111
1112   //__________________________________________________________________________
1113 void AliMUONTrackReconstructor::Fit(AliMUONTrack &track, Bool_t includeMCS, Bool_t fitWithVertex, Bool_t calcCov)
1114 {
1115   /// Fit the track
1116   /// w/wo multiple Coulomb scattering according to "includeMCS".
1117   /// w/wo constraining the vertex according to "fitWithVertex".
1118   /// calculating or not the covariance matrix according to "calcCov".
1119   AliDebug(1,"Enter Fit");
1120   
1121   Double_t benC, errorParam, invBenP, nonBenC, x, y;
1122   AliMUONTrackParam *trackParam;
1123   Double_t arg[1], fedm, errdef, globalChi2;
1124   Int_t npari, nparx;
1125   Int_t status, covStatus;
1126   
1127   // Instantiate gMinuit if not already done
1128   if (!gMinuit) gMinuit = new TMinuit(6);
1129   // Clear MINUIT parameters
1130   gMinuit->mncler();
1131   // Give the fitted track to MINUIT
1132   gMinuit->SetObjectFit(&track);
1133   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
1134     // Define print level
1135     arg[0] = 1;
1136     gMinuit->mnexcm("SET PRI", arg, 1, status);
1137     // Print covariance matrix
1138     gMinuit->mnexcm("SHO COV", arg, 0, status);
1139   } else {
1140     arg[0] = -1;
1141     gMinuit->mnexcm("SET PRI", arg, 1, status);
1142   }
1143   // No warnings
1144   gMinuit->mnexcm("SET NOW", arg, 0, status);
1145   // Define strategy
1146   //arg[0] = 2;
1147   //gMinuit->mnexcm("SET STR", arg, 1, status);
1148   
1149   // set flag w/wo multiple scattering according to "includeMCS"
1150   track.FitWithMCS(includeMCS);
1151   if (includeMCS) {
1152     // compute cluster weights only once
1153     if (!track.ComputeClusterWeights()) {
1154       AliWarning("cannot take into account the multiple scattering effects");
1155       track.FitWithMCS(kFALSE);
1156     }
1157   }
1158   
1159   track.FitWithVertex(fitWithVertex);
1160   if (fitWithVertex) SetVertexErrXY2ForFit(track);
1161   
1162   // Set fitting function
1163   gMinuit->SetFCN(TrackChi2);
1164   
1165   // Set fitted parameters (!! The order is very important for the covariance matrix !!)
1166   // Mandatory limits to avoid NaN values of parameters
1167   trackParam = (AliMUONTrackParam*) (track.GetTrackParamAtCluster()->First());
1168   Double_t maxIBM = 1. / GetRecoParam()->GetMinBendingMomentum();
1169   gMinuit->mnparm(0, "X", trackParam->GetNonBendingCoor(), 0.03, -500.0, 500.0, status);
1170   gMinuit->mnparm(1, "NonBenS", trackParam->GetNonBendingSlope(), 0.001, -1., 1., status);
1171   gMinuit->mnparm(2, "Y", trackParam->GetBendingCoor(), 0.10, -500.0, 500.0, status);
1172   gMinuit->mnparm(3, "BenS", trackParam->GetBendingSlope(), 0.001, -1.5, 1.5, status);
1173   gMinuit->mnparm(4, "InvBenP", trackParam->GetInverseBendingMomentum(), 0.003, -maxIBM, maxIBM, status);
1174   
1175   // minimization
1176   gMinuit->mnexcm("MIGRAD", arg, 0, status);
1177   
1178   // Calculate the covariance matrix more accurately if required
1179   if (calcCov) gMinuit->mnexcm("HESSE", arg, 0, status);
1180    
1181   // get results into "invBenP", "benC", "nonBenC" ("x", "y")
1182   gMinuit->GetParameter(0, x, errorParam);
1183   trackParam->SetNonBendingCoor(x);
1184   gMinuit->GetParameter(1, nonBenC, errorParam);
1185   trackParam->SetNonBendingSlope(nonBenC);
1186   gMinuit->GetParameter(2, y, errorParam);
1187   trackParam->SetBendingCoor(y);
1188   gMinuit->GetParameter(3, benC, errorParam);
1189   trackParam->SetBendingSlope(benC);
1190   gMinuit->GetParameter(4, invBenP, errorParam);
1191   trackParam->SetInverseBendingMomentum(invBenP);
1192   
1193   // global result of the fit
1194   gMinuit->mnstat(globalChi2, fedm, errdef, npari, nparx, covStatus);
1195   track.SetGlobalChi2(globalChi2);
1196   
1197   // Get the covariance matrix if required
1198   if (calcCov) {
1199     // Covariance matrix according to HESSE status
1200     // If problem then keep only the diagonal terms (variances)
1201     Double_t matrix[5][5];
1202     gMinuit->mnemat(&matrix[0][0],5);
1203     if (covStatus == 3) trackParam->SetCovariances(matrix);
1204     else trackParam->SetVariances(matrix);
1205   } else trackParam->DeleteCovariances();
1206   
1207 }
1208
1209   //__________________________________________________________________________
1210 void TrackChi2(Int_t & /*nParam*/, Double_t * /*gradient*/, Double_t &chi2, Double_t *param, Int_t /*flag*/)
1211 {
1212   /// Return the "Chi2" to be minimized with Minuit for track fitting.
1213   /// Assumes that the trackParamAtCluster are sorted according to increasing Z.
1214   /// Track parameters at each cluster are updated accordingly.
1215   /// Vertex is used according to the flag "trackBeingFitted->GetFitWithVertex()".
1216   /// Multiple Coulomb scattering is taken into account according to the flag "trackBeingFitted->GetFitWithMCS()".
1217   
1218   AliMUONTrack *trackBeingFitted = (AliMUONTrack*) gMinuit->GetObjectFit();
1219   AliMUONTrackParam* trackParamAtCluster = (AliMUONTrackParam*) trackBeingFitted->GetTrackParamAtCluster()->First();
1220   Double_t dX, dY;
1221   chi2 = 0.; // initialize chi2
1222   
1223   // update track parameters
1224   trackParamAtCluster->SetNonBendingCoor(param[0]);
1225   trackParamAtCluster->SetNonBendingSlope(param[1]);
1226   trackParamAtCluster->SetBendingCoor(param[2]);
1227   trackParamAtCluster->SetBendingSlope(param[3]);
1228   trackParamAtCluster->SetInverseBendingMomentum(param[4]);
1229   trackBeingFitted->UpdateTrackParamAtCluster();
1230   
1231   // Take the vertex into account in the fit if required
1232   if (trackBeingFitted->FitWithVertex()) {
1233     Double_t nonBendingReso2,bendingReso2;
1234     trackBeingFitted->GetVertexErrXY2(nonBendingReso2,bendingReso2);
1235     if (nonBendingReso2 == 0. || bendingReso2 == 0.) chi2 += 1.e10;
1236     else {
1237       AliMUONTrackParam paramAtVertex(*trackParamAtCluster);
1238       AliMUONTrackExtrap::ExtrapToZ(&paramAtVertex, 0.); // vextex position = (0,0,0)
1239       dX = paramAtVertex.GetNonBendingCoor();
1240       dY = paramAtVertex.GetBendingCoor();
1241       chi2 += dX * dX / nonBendingReso2 + dY * dY / bendingReso2;
1242     }
1243   }
1244   
1245   // compute chi2 w/wo multiple scattering
1246   chi2 += trackBeingFitted->ComputeGlobalChi2(trackBeingFitted->FitWithMCS());
1247   
1248 }
1249
1250   //__________________________________________________________________________
1251 void AliMUONTrackReconstructor::ComplementTracks(const AliMUONVClusterStore& clusterStore)
1252 {
1253   /// Complete tracks by adding missing clusters (if there is an overlap between
1254   /// two detection elements, the track may have two clusters in the same chamber)
1255   /// Re-fit track parameters and covariances
1256   AliDebug(1,"Enter ComplementTracks");
1257   
1258   Int_t chamberId, detElemId;
1259   Double_t chi2OfCluster, bestChi2OfCluster;
1260   Double_t sigmaCut2 = GetRecoParam()->GetSigmaCutForTracking() *
1261                        GetRecoParam()->GetSigmaCutForTracking();
1262   Bool_t foundOneCluster, trackModified;
1263   AliMUONVCluster* cluster;
1264   AliMUONTrackParam *trackParam, *nextTrackParam, copyOfTrackParam, trackParamAtCluster, bestTrackParamAtCluster;
1265   
1266   AliMUONTrack *track = (AliMUONTrack*) fRecTracksPtr->First();
1267   while (track) {
1268     trackModified = kFALSE;
1269     
1270     trackParam = (AliMUONTrackParam*)track->GetTrackParamAtCluster()->First();
1271     while (trackParam) {
1272       foundOneCluster = kFALSE;
1273       bestChi2OfCluster = 2. * sigmaCut2; // 2 because 2 quantities in chi2
1274       chamberId = trackParam->GetClusterPtr()->GetChamberId();
1275       detElemId = trackParam->GetClusterPtr()->GetDetElemId();
1276       
1277       // prepare nextTrackParam before adding new cluster because of the sorting
1278       nextTrackParam = (AliMUONTrackParam*)track->GetTrackParamAtCluster()->After(trackParam);
1279       
1280       // recover track parameters from local fit and put them into a copy of trackParam
1281       copyOfTrackParam.SetZ(trackParam->GetZ());
1282       copyOfTrackParam.SetParameters(trackParam->GetSmoothParameters());
1283       copyOfTrackParam.SetCovariances(trackParam->GetSmoothCovariances());
1284       
1285       // Create iterators to loop over clusters in current chamber
1286       TIter nextInCh(clusterStore.CreateChamberIterator(chamberId,chamberId));
1287       
1288       // look for one second candidate in the same chamber
1289       while ( ( cluster = static_cast<AliMUONVCluster*>(nextInCh()) ) ) {
1290         
1291         // look for a cluster in another detection element
1292         if (cluster->GetDetElemId() == detElemId) continue;
1293         
1294         // try to add the current cluster fast
1295         if (!TryOneClusterFast(copyOfTrackParam, cluster)) continue;
1296         
1297         // try to add the current cluster accurately
1298         chi2OfCluster = TryOneCluster(copyOfTrackParam, cluster, trackParamAtCluster);
1299         
1300         // if better chi2 then prepare to add this cluster to the track
1301         if (chi2OfCluster < bestChi2OfCluster) {
1302           bestChi2OfCluster = chi2OfCluster;
1303           bestTrackParamAtCluster = trackParamAtCluster;
1304           foundOneCluster = kTRUE;
1305         }
1306         
1307       }
1308       
1309       // add new cluster if any
1310       if (foundOneCluster) {
1311         
1312         // Printout for debuging
1313         if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
1314           cout << "ComplementTracks: found one cluster in chamber(1..): " << chamberId+1 << endl;
1315           bestTrackParamAtCluster.GetClusterPtr()->Print();
1316           cout<<endl<<"Track parameters and covariances at cluster:"<<endl;
1317           bestTrackParamAtCluster.GetParameters().Print();
1318           bestTrackParamAtCluster.GetCovariances().Print();
1319         }
1320         
1321         trackParam->SetRemovable(kTRUE);
1322         bestTrackParamAtCluster.SetRemovable(kTRUE);
1323         track->AddTrackParamAtCluster(bestTrackParamAtCluster,*(bestTrackParamAtCluster.GetClusterPtr()));
1324         trackModified = kTRUE;
1325       }
1326       
1327       trackParam = nextTrackParam;
1328     }
1329     
1330     // re-fit track parameters if needed
1331     if (trackModified) Fit(*track, kTRUE, kFALSE, kTRUE);
1332     
1333     track = (AliMUONTrack*) fRecTracksPtr->After(track);
1334   }
1335   
1336 }
1337
1338   //__________________________________________________________________________
1339 void AliMUONTrackReconstructor::ImproveTrack(AliMUONTrack &track)
1340 {
1341   /// Improve the given track by removing clusters with local chi2 highter than the defined cut
1342   /// Recompute track parameters and covariances at the remaining clusters
1343   AliDebug(1,"Enter ImproveTrack");
1344   
1345   Double_t localChi2, worstLocalChi2;
1346   AliMUONTrackParam *trackParamAtCluster, *worstTrackParamAtCluster;
1347   Double_t sigmaCut2 = GetRecoParam()->GetSigmaCutForImprovement() *
1348                        GetRecoParam()->GetSigmaCutForImprovement();
1349   
1350   while (!track.IsImproved()) {
1351     
1352     // identify removable clusters
1353     track.TagRemovableClusters(GetRecoParam()->RequestedStationMask());
1354     
1355     // Update track parameters and covariances
1356     track.UpdateCovTrackParamAtCluster();
1357     
1358     // Compute local chi2 of each clusters
1359     track.ComputeLocalChi2(kTRUE);
1360     
1361     // Look for the cluster to remove
1362     worstTrackParamAtCluster = NULL;
1363     worstLocalChi2 = 0.;
1364     trackParamAtCluster = (AliMUONTrackParam*) track.GetTrackParamAtCluster()->First();
1365     while (trackParamAtCluster) {
1366       
1367       // Pick up cluster with the worst chi2
1368       localChi2 = trackParamAtCluster->GetLocalChi2();
1369       if (localChi2 > worstLocalChi2) {
1370         worstLocalChi2 = localChi2;
1371         worstTrackParamAtCluster = trackParamAtCluster;
1372       }
1373       
1374     trackParamAtCluster = (AliMUONTrackParam*) track.GetTrackParamAtCluster()->After(trackParamAtCluster);
1375     }
1376     
1377     // Check if worst cluster found
1378     if (!worstTrackParamAtCluster) {
1379       AliWarning("Bad local chi2 values?");
1380       break;
1381     }
1382     
1383     // Check whether the worst chi2 is under requirement or not
1384     if (worstLocalChi2 < 2. * sigmaCut2) { // 2 because 2 quantities in chi2
1385       track.SetImproved(kTRUE);
1386       break;
1387     }
1388     
1389     // if the worst cluster is not removable then stop improvement
1390     if (!worstTrackParamAtCluster->IsRemovable()) break;
1391     
1392     // Remove the worst cluster
1393     track.RemoveTrackParamAtCluster(worstTrackParamAtCluster);
1394     
1395     // Re-fit the track:
1396     // Take into account the multiple scattering
1397     // Calculate the track parameter covariance matrix
1398     Fit(track, kTRUE, kFALSE, kTRUE);
1399     
1400     // Printout for debuging
1401     if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
1402       cout << "ImproveTracks: track " << fRecTracksPtr->IndexOf(&track)+1 << " improved " << endl;
1403     }
1404     
1405   }
1406   
1407 }
1408
1409   //__________________________________________________________________________
1410 void AliMUONTrackReconstructor::FinalizeTrack(AliMUONTrack &track)
1411 {
1412   /// Recompute track parameters and covariances at each attached cluster
1413   /// from those at the first one, if not already done
1414   AliDebug(1,"Enter FinalizeTrack");
1415   if (!track.IsImproved()) track.UpdateCovTrackParamAtCluster();
1416 }
1417
1418   //__________________________________________________________________________
1419 Bool_t AliMUONTrackReconstructor::RefitTrack(AliMUONTrack &track, Bool_t enableImprovement)
1420 {
1421   /// re-fit the given track
1422   
1423   // check validity of the track
1424   if (track.GetNClusters() < 3) {
1425     AliWarning("the track does not contain enough clusters --> unable to refit");
1426     return kFALSE;
1427   }
1428   
1429   // reset the seed (i.e. parameters at first cluster) before fitting
1430   AliMUONTrackParam* firstTrackParam = (AliMUONTrackParam*) track.GetTrackParamAtCluster()->First();
1431   if (firstTrackParam->GetInverseBendingMomentum() == 0.) {
1432     AliWarning("track parameters at first chamber are not initialized --> unable to refit");
1433     return kFALSE;
1434   }
1435   
1436   // compute track parameters at each cluster from parameters at the first one
1437   // necessary to compute multiple scattering effect during refitting
1438   track.UpdateTrackParamAtCluster();
1439   
1440   // Re-fit the track:
1441   // Take into account the multiple scattering
1442   // Calculate the track parameter covariance matrix
1443   Fit(track, kTRUE, kFALSE, kTRUE);
1444   
1445   // Improve the reconstructed tracks if required
1446   track.SetImproved(kFALSE);
1447   if (enableImprovement && GetRecoParam()->ImproveTracks()) ImproveTrack(track);
1448   
1449   // Fill AliMUONTrack data members
1450   FinalizeTrack(track);
1451   
1452   return kTRUE;
1453   
1454 }
1455