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