+ Modifications of the standard tracking algorithm:
[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 /// \class AliMUONTrackReconstructor
19 /// MUON track reconstructor using the original method
20 ///
21 /// This class contains as data:
22 /// - the parameters for the track reconstruction
23 ///
24 /// It contains as methods, among others:
25 /// - MakeTracks to build the tracks
26 ///
27
28 #include <stdlib.h>
29 #include <Riostream.h>
30 #include <TMatrixD.h>
31
32 #include "AliMUONVTrackReconstructor.h"
33 #include "AliMUONTrackReconstructor.h"
34 #include "AliMUONData.h"
35 #include "AliMUONConstants.h"
36 #include "AliMUONRawCluster.h"
37 #include "AliMUONHitForRec.h"
38 #include "AliMUONObjectPair.h"
39 #include "AliMUONTrack.h"
40 #include "AliMUONTrackParam.h"
41 #include "AliMUONTrackExtrap.h"
42 #include "AliLog.h"
43 #include <TMinuit.h>
44
45 // Functions to be minimized with Minuit
46 void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag);
47 void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag);
48
49 Double_t MultipleScatteringAngle2(AliMUONTrackParam *param);
50
51 /// \cond CLASSIMP
52 ClassImp(AliMUONTrackReconstructor) // Class implementation in ROOT context
53 /// \endcond
54
55 //************* Defaults parameters for reconstruction
56 const Double_t AliMUONTrackReconstructor::fgkMaxNormChi2 = 100.0;
57 const Bool_t AliMUONTrackReconstructor::fgkTrackAllTracks = kFALSE;
58
59 //__________________________________________________________________________
60 AliMUONTrackReconstructor::AliMUONTrackReconstructor(AliMUONData* data)
61   : AliMUONVTrackReconstructor(data)
62 {
63   /// Constructor for class AliMUONTrackReconstructor
64   
65   // Memory allocation for the TClonesArray of reconstructed tracks
66   fRecTracksPtr = new TClonesArray("AliMUONTrack", 10);
67 }
68
69   //__________________________________________________________________________
70 AliMUONTrackReconstructor::~AliMUONTrackReconstructor(void)
71 {
72   /// Destructor for class AliMUONTrackReconstructor
73   delete fRecTracksPtr;
74 }
75
76   //__________________________________________________________________________
77 void AliMUONTrackReconstructor::AddHitsForRecFromRawClusters()
78 {
79   /// To add to the list of hits for reconstruction all the raw clusters
80   TTree *treeR;
81   AliMUONHitForRec *hitForRec;
82   AliMUONRawCluster *clus;
83   Int_t iclus, nclus;
84   TClonesArray *rawclusters;
85   AliDebug(1,"Enter AddHitsForRecFromRawClusters");
86   
87   treeR = fMUONData->TreeR();
88   if (!treeR) {
89     AliError("TreeR must be loaded");
90     exit(0);
91   }
92   
93   fMUONData->SetTreeAddress("RC");
94   fMUONData->GetRawClusters(); // only one entry  
95   
96   // Loop over tracking chambers
97   for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
98     rawclusters =fMUONData->RawClusters(ch);
99     nclus = (Int_t) (rawclusters->GetEntries());
100     // Loop over (cathode correlated) raw clusters
101     for (iclus = 0; iclus < nclus; iclus++) {
102       clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(iclus);
103       // new AliMUONHitForRec from raw cluster
104       // and increment number of AliMUONHitForRec's (total and in chamber)
105       hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(clus);
106       fNHitsForRec++;
107       // more information into HitForRec
108       hitForRec->SetBendingReso2(clus->GetErrY() * clus->GetErrY());
109       hitForRec->SetNonBendingReso2(clus->GetErrX() * clus->GetErrX());
110       //  original raw cluster
111       hitForRec->SetChamberNumber(ch);
112       hitForRec->SetHitNumber(iclus);
113       // Z coordinate of the raw cluster (cm)
114       hitForRec->SetZ(clus->GetZ(0));
115       StdoutToAliDebug(3,
116                        cout << "Chamber " << ch <<
117                        " raw cluster  " << iclus << " : " << endl;
118                        clus->Print("full");
119                        cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl;
120                        hitForRec->Print("full");
121                        );
122     } // end of cluster loop
123   } // end of chamber loop
124   SortHitsForRecWithIncreasingChamber(); 
125   
126   AliDebug(1,"End of AddHitsForRecFromRawClusters");
127     if (AliLog::GetGlobalDebugLevel() > 0) {
128       AliDebug(1, Form("NHitsForRec: %d",fNHitsForRec));
129       for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
130         AliDebug(1, Form("Chamber(0...): %d",ch));
131         AliDebug(1, Form("NHitsForRec: %d", fNHitsForRecPerChamber[ch]));
132         AliDebug(1, Form("Index(first HitForRec): %d", fIndexOfFirstHitForRecPerChamber[ch]));
133         for (Int_t hit = fIndexOfFirstHitForRecPerChamber[ch];
134              hit < fIndexOfFirstHitForRecPerChamber[ch] + fNHitsForRecPerChamber[ch];
135              hit++) {
136           AliDebug(1, Form("HitForRec index(0...): %d",hit));
137           ((*fHitsForRecPtr)[hit])->Dump();
138       }
139     }
140   }
141   
142   return;
143 }
144
145   //__________________________________________________________________________
146 void AliMUONTrackReconstructor::MakeTracks(void)
147 {
148   /// To make the tracks from the list of segments and points in all stations
149   AliDebug(1,"Enter MakeTracks");
150   // Look for candidates from at least 3 aligned points in stations(1..) 4 and 5
151   MakeTrackCandidates();
152   // Follow tracks in stations(1..) 3, 2 and 1
153   FollowTracks();
154   // Remove double tracks
155   RemoveDoubleTracks();
156   // Propagate tracks to the vertex through absorber
157   ExtrapTracksToVertex();
158   // Fill out the AliMUONTrack's
159   FillMUONTrack();
160 }
161
162   //__________________________________________________________________________
163 void AliMUONTrackReconstructor::MakeTrackCandidates(void)
164 {
165   /// To make track candidates:
166   /// Start with segments station(1..) 4 or 5 then follow track in station 5 or 4.
167   /// Good candidates are made of at least three hitForRec's.
168   /// Keep only best candidates or all of them according to the flag fgkTrackAllTracks.
169   TClonesArray *segments;
170   AliMUONObjectPair *segment;
171   AliMUONHitForRec *hitForRec1, *hitForRec2;
172   AliMUONTrack *track;
173   AliMUONTrackParam *trackParamAtFirstHit;
174   Int_t iCandidate = 0;
175
176   AliDebug(1,"Enter MakeTrackCandidates");
177
178   // Loop over stations(1..) 5 and 4 and make track candidates
179   for (Int_t istat=4; istat>=3; istat--) {
180     // Make segments in the station
181     segments = MakeSegmentsInStation(istat);
182     // Loop over segments
183     for (Int_t iseg=0; iseg<segments->GetEntriesFast(); iseg++) {
184       AliDebug(1,Form("Making primary candidate(1..) %d",++iCandidate));
185       // Transform segments to tracks and put them at the end of fRecTracksPtr
186       segment = (AliMUONObjectPair*) ((*segments)[iseg]);
187       hitForRec1 = (AliMUONHitForRec*) segment->First();
188       hitForRec2 = (AliMUONHitForRec*) segment->Second();
189       track = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(hitForRec1, hitForRec2);
190       fNRecTracks++;
191       // Add MCS effects in parameter covariances
192       trackParamAtFirstHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First());
193       AliMUONTrackExtrap::AddMCSEffectInTrackParamCov(trackParamAtFirstHit,AliMUONConstants::ChamberThicknessInX0(),1.);
194       // Printout for debuging
195       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
196         cout<<endl<<"Track parameter covariances at first hit with multiple Coulomb scattering effects:"<<endl;
197         trackParamAtFirstHit->GetCovariances()->Print();
198       }
199       // Look for compatible hitForRec(s) in the other station
200       FollowTrackInStation(track,7-istat);
201     }
202     // delete the array of segments
203     delete segments;
204   }
205   fRecTracksPtr->Compress(); // this is essential before checking tracks
206   
207   // Keep all different tracks or only the best ones as required
208   if (fgkTrackAllTracks) RemoveIdenticalTracks();
209   else RemoveDoubleTracks();
210   
211   AliDebug(1,Form("Number of good candidates = %d",fNRecTracks));
212   
213 }
214
215   //__________________________________________________________________________
216 void AliMUONTrackReconstructor::RemoveIdenticalTracks(void)
217 {
218   /// To remove identical tracks:
219   /// Tracks are considered identical if they have all their hits in common.
220   /// One keeps the track with the larger number of hits if need be
221   AliMUONTrack *track1, *track2, *trackToRemove;
222   Int_t hitsInCommon, nHits1, nHits2;
223   Bool_t removedTrack1;
224   // Loop over first track of the pair
225   track1 = (AliMUONTrack*) fRecTracksPtr->First();
226   while (track1) {
227     removedTrack1 = kFALSE;
228     nHits1 = track1->GetNTrackHits();
229     // Loop over second track of the pair
230     track2 = (AliMUONTrack*) fRecTracksPtr->After(track1);
231     while (track2) {
232       nHits2 = track2->GetNTrackHits();
233       // number of hits in common between two tracks
234       hitsInCommon = track1->HitsInCommon(track2);
235       // check for identical tracks
236       if ((hitsInCommon == nHits1) || (hitsInCommon == nHits2)) {
237         // decide which track to remove
238         if (nHits2 > nHits1) {
239           // remove track1 and continue the first loop with the track next to track1
240           trackToRemove = track1;
241           track1 = (AliMUONTrack*) fRecTracksPtr->After(track1);
242           fRecTracksPtr->Remove(trackToRemove);
243           fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards
244           fNRecTracks--;
245           removedTrack1 = kTRUE;
246           break;
247         } else {
248           // remove track2 and continue the second loop with the track next to track2
249           trackToRemove = track2;
250           track2 = (AliMUONTrack*) fRecTracksPtr->After(track2);
251           fRecTracksPtr->Remove(trackToRemove);
252           fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards
253           fNRecTracks--;
254         }
255       } else track2 = (AliMUONTrack*) fRecTracksPtr->After(track2);
256     } // track2
257     if (removedTrack1) continue;
258     track1 = (AliMUONTrack*) fRecTracksPtr->After(track1);
259   } // track1
260   return;
261 }
262
263   //__________________________________________________________________________
264 void AliMUONTrackReconstructor::RemoveDoubleTracks(void)
265 {
266   /// To remove double tracks:
267   /// Tracks are considered identical if more than half of the hits of the track
268   /// which has the smaller number of hits are in common with the other track.
269   /// Among two identical tracks, one keeps the track with the larger number of hits
270   /// or, if these numbers are equal, the track with the minimum chi2.
271   AliMUONTrack *track1, *track2, *trackToRemove;
272   Int_t hitsInCommon, nHits1, nHits2;
273   Bool_t removedTrack1;
274   // Loop over first track of the pair
275   track1 = (AliMUONTrack*) fRecTracksPtr->First();
276   while (track1) {
277     removedTrack1 = kFALSE;
278     nHits1 = track1->GetNTrackHits();
279     // Loop over second track of the pair
280     track2 = (AliMUONTrack*) fRecTracksPtr->After(track1);
281     while (track2) {
282       nHits2 = track2->GetNTrackHits();
283       // number of hits in common between two tracks
284       hitsInCommon = track1->HitsInCommon(track2);
285       // check for identical tracks
286       if (((nHits1 < nHits2) && (2 * hitsInCommon > nHits1)) || (2 * hitsInCommon > nHits2)) {
287         // decide which track to remove
288         if ((nHits1 > nHits2) || ((nHits1 == nHits2) && (track1->GetFitFMin() <= track2->GetFitFMin()))) {
289           // remove track2 and continue the second loop with the track next to track2
290           trackToRemove = track2;
291           track2 = (AliMUONTrack*) fRecTracksPtr->After(track2);
292           fRecTracksPtr->Remove(trackToRemove);
293           fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards
294           fNRecTracks--;
295         } else {
296           // else remove track1 and continue the first loop with the track next to track1
297           trackToRemove = track1;
298           track1 = (AliMUONTrack*) fRecTracksPtr->After(track1);
299           fRecTracksPtr->Remove(trackToRemove);
300           fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards
301           fNRecTracks--;
302           removedTrack1 = kTRUE;
303           break;
304         }
305       } else track2 = (AliMUONTrack*) fRecTracksPtr->After(track2);
306     } // track2
307     if (removedTrack1) continue;
308     track1 = (AliMUONTrack*) fRecTracksPtr->After(track1);
309   } // track1
310   return;
311 }
312
313   //__________________________________________________________________________
314 void AliMUONTrackReconstructor::FollowTracks(void)
315 {
316   /// Follow tracks in stations(1..) 3, 2 and 1
317   AliDebug(1,"Enter FollowTracks");
318   
319   AliMUONTrack *track, *nextTrack;
320   AliMUONTrackParam *trackParamAtFirstHit;
321   Double_t numberOfDegFree, chi2Norm;
322   Int_t CurrentNRecTracks;
323   
324   for (Int_t station = 2; station >= 0; station--) {
325     // Save the actual number of reconstructed track in case of
326     // tracks are added or suppressed during the tracking procedure
327     // !! Do not compress fRecTracksPtr until the end of the loop over tracks !!
328     CurrentNRecTracks = fNRecTracks;
329     for (Int_t iRecTrack = 0; iRecTrack <CurrentNRecTracks; iRecTrack++) {
330       AliDebug(1,Form("FollowTracks: track candidate(1..) %d", iRecTrack+1));
331       track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iRecTrack);
332       // Fit the track:
333       // Do not take into account the multiple scattering to speed up the fit
334       // Calculate the track parameter covariance matrix
335       // If "station" is station(1..) 3 then use the vertex to better constrain the fit
336       if (station==2) {
337         SetVertexForFit(track);
338         track->SetFitWithVertex(kTRUE);
339       } else track->SetFitWithVertex(kFALSE);
340       Fit(track,kFALSE, kTRUE);
341       // Remove the track if the normalized chi2 is too high
342       numberOfDegFree = (2. * track->GetNTrackHits() - 5.);
343       if (numberOfDegFree > 0) chi2Norm = track->GetFitFMin() / numberOfDegFree;
344       else chi2Norm = 1.e10;
345       if (chi2Norm > fgkMaxNormChi2) {
346         fRecTracksPtr->Remove(track);
347         fNRecTracks--;
348         continue;
349       }
350       // Add MCS effects in parameter covariances
351       trackParamAtFirstHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First());
352       AliMUONTrackExtrap::AddMCSEffectInTrackParamCov(trackParamAtFirstHit,AliMUONConstants::ChamberThicknessInX0(),1.);
353       // Printout for debuging
354       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
355         cout<<endl<<"Track parameter covariances at first hit with multiple Coulomb scattering effects:"<<endl;
356         trackParamAtFirstHit->GetCovariances()->Print();
357       }
358       // Look for compatible hitForRec in station(0..) "station"
359       FollowTrackInStation(track,station);
360     }
361     // Compress fRecTracksPtr for the next step
362     fRecTracksPtr->Compress();
363     // Keep only the best tracks if required
364     if (!fgkTrackAllTracks) RemoveDoubleTracks();
365   }
366   
367   // Last fit of track candidates with all station
368   // Take into account the multiple scattering and remove bad tracks
369   Int_t trackIndex = -1;
370   track = (AliMUONTrack*) fRecTracksPtr->First();
371   while (track) {
372     trackIndex++;
373     nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track); // prepare next track
374     track->SetFitWithVertex(kFALSE); // just to be sure
375     Fit(track,kTRUE, kFALSE);
376     // Printout for debuging
377     if (AliLog::GetGlobalDebugLevel() >= 3) {
378       cout << "FollowTracks: track candidate(0..) " << trackIndex << " after final fit" << endl;
379       track->RecursiveDump();
380     } 
381     // Remove the track if the normalized chi2 is too high
382     numberOfDegFree = (2.0 * track->GetNTrackHits() - 5);
383     if (numberOfDegFree > 0) chi2Norm = track->GetFitFMin() / numberOfDegFree;
384     else chi2Norm = 1.e10;
385     if (chi2Norm > fgkMaxNormChi2) {
386       fRecTracksPtr->Remove(track);
387       fNRecTracks--;
388     }
389     track = nextTrack;
390   }
391   fRecTracksPtr->Compress();
392   
393 }
394
395   //__________________________________________________________________________
396 void AliMUONTrackReconstructor::FollowTrackInStation(AliMUONTrack* trackCandidate, Int_t nextStation)
397 {
398   /// Follow trackCandidate in station(0..) nextStation and search for compatible HitForRec(s)
399   /// Keep all possibilities or only the best one(s) according to the flag fgkTrackAllTracks:
400   /// kTRUE:  duplicate "trackCandidate" if there are several possibilities and add the new tracks at the end of
401   ///         fRecTracksPtr to avoid conficts with other track candidates at this current stage of the tracking procedure.
402   ///         Remove the obsolete "trackCandidate" at the end.
403   /// kFALSE: add only the best hit(s) to the "trackCandidate". Try to add a couple of hits in priority.
404   AliDebug(1,Form("Enter FollowTrackInStation(1..) %d", nextStation+1));
405   
406   Int_t ch1 = 2*nextStation;
407   Int_t ch2 = 2*nextStation+1;
408   Double_t zCh2 = AliMUONConstants::DefaultChamberZ(ch2);
409   Double_t chi2WithOneHitForRec = 1.e10;
410   Double_t chi2WithTwoHitForRec = 1.e10;
411   Double_t maxChi2WithOneHitForRec = 2.*fgkMaxNormChi2; // 2 because 2 quantities in chi2
412   Double_t maxChi2WithTwoHitForRec = 4.*fgkMaxNormChi2; // 4 because 4 quantities in chi2
413   Double_t bestChi2WithOneHitForRec = maxChi2WithOneHitForRec;
414   Double_t bestChi2WithTwoHitForRec = maxChi2WithTwoHitForRec;
415   AliMUONTrack *newTrack = 0x0;
416   AliMUONHitForRec *hitForRecCh1, *hitForRecCh2;
417   AliMUONHitForRec *bestHitForRec1 = 0x0, *bestHitForRec2 = 0x0;
418   //
419   //Extrapolate trackCandidate to chamber "ch2" to save computing time in the next steps
420   AliMUONTrackParam *extrapTrackParamPtr = trackCandidate->GetExtrapTrackParam();
421   *extrapTrackParamPtr = *((AliMUONTrackParam*)(trackCandidate->GetTrackParamAtHit()->First()));
422   AliMUONTrackExtrap::ExtrapToZCov(extrapTrackParamPtr, zCh2);
423   AliMUONTrackParam extrapTrackParamSave(*extrapTrackParamPtr);
424   //
425   // Printout for debuging
426   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
427     TMatrixD* paramCovForDebug = extrapTrackParamPtr->GetCovariances();
428     cout<<endl<<"Track parameter covariances at first hit extrapolated to z = "<<zCh2<<":"<<endl;
429     paramCovForDebug->Print();
430   }
431   //
432   // look for candidates in chamber 2 
433   // Printout for debuging
434   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
435     cout << "FollowTrackInStation: look for hits in chamber(1..): " << ch2+1 << endl;
436   }
437   for (Int_t hit2 = 0; hit2 < fNHitsForRecPerChamber[ch2]; hit2++) {
438     hitForRecCh2 = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[ch2]+hit2);
439     // extrapolate track parameters and covariances only once for this hit
440     AliMUONTrackExtrap::ExtrapToZCov(extrapTrackParamPtr, hitForRecCh2->GetZ());
441     chi2WithOneHitForRec = trackCandidate->TryOneHitForRec(hitForRecCh2);
442     // if good chi2 then try to attach a hitForRec in the other chamber too
443     if (chi2WithOneHitForRec < maxChi2WithOneHitForRec) {
444       // Printout for debuging
445       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
446         cout << "FollowTrackInStation: look for second hits in chamber(1..): " << ch1+1 << endl;
447       }
448       Bool_t foundSecondHit = kFALSE;
449       for (Int_t hit1 = 0; hit1 < fNHitsForRecPerChamber[ch1]; hit1++) {
450         hitForRecCh1 = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[ch1]+hit1);
451         chi2WithTwoHitForRec = trackCandidate->TryTwoHitForRec(hitForRecCh2, hitForRecCh1); // order hits like that to save computing time
452         // if good chi2 then create a new track by adding the 2 hitForRec to the "trackCandidate"
453         if (chi2WithTwoHitForRec < maxChi2WithTwoHitForRec) {
454           foundSecondHit = kTRUE;
455           if (fgkTrackAllTracks) {
456             // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new hitForRec's
457             newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(*trackCandidate);
458             fNRecTracks++;
459             AliMUONTrackParam TrackParam1(extrapTrackParamSave);
460             AliMUONTrackExtrap::ExtrapToZ(&TrackParam1, hitForRecCh1->GetZ());
461             newTrack->AddTrackParamAtHit(&TrackParam1,hitForRecCh1);
462             AliMUONTrackParam TrackParam2(extrapTrackParamSave);
463             AliMUONTrackExtrap::ExtrapToZ(&TrackParam2, hitForRecCh2->GetZ());
464             newTrack->AddTrackParamAtHit(&TrackParam2,hitForRecCh2);
465             // Sort TrackParamAtHit according to increasing Z
466             newTrack->GetTrackParamAtHit()->Sort();
467             // Update the chi2 of the new track
468             if (newTrack->GetFitFMin()<0) newTrack->SetFitFMin(chi2WithTwoHitForRec);
469             else newTrack->SetFitFMin(newTrack->GetFitFMin() + chi2WithTwoHitForRec);
470             // Printout for debuging
471             if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
472               cout << "FollowTrackInStation: added two hits in station(1..): " << nextStation+1
473                    << " (Chi2 = " << chi2WithTwoHitForRec << ")" << endl;
474               if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
475             }
476           } else if (chi2WithTwoHitForRec < bestChi2WithTwoHitForRec) {
477             // keep track of the best couple of hits
478             bestChi2WithTwoHitForRec = chi2WithTwoHitForRec;
479             bestHitForRec1 = hitForRecCh1;
480             bestHitForRec2 = hitForRecCh2;
481           }
482         }
483       }
484       // if no hitForRecCh1 found then consider to add hitForRecCh2 only
485       if (!foundSecondHit) {
486         if (fgkTrackAllTracks) {
487           // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new hitForRec's
488           newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(*trackCandidate);
489           fNRecTracks++;
490           AliMUONTrackParam TrackParam1(extrapTrackParamSave);
491           AliMUONTrackExtrap::ExtrapToZ(&TrackParam1, hitForRecCh2->GetZ());
492           newTrack->AddTrackParamAtHit(&TrackParam1,hitForRecCh2);
493           // Sort TrackParamAtHit according to increasing Z
494           newTrack->GetTrackParamAtHit()->Sort();
495           // Update the chi2 of the new track
496           if (newTrack->GetFitFMin()<0) newTrack->SetFitFMin(chi2WithOneHitForRec);
497           else newTrack->SetFitFMin(newTrack->GetFitFMin() + chi2WithOneHitForRec);
498           // Printout for debuging
499           if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
500             cout << "FollowTrackInStation: added one hit in chamber(1..): " << ch2+1
501                  << " (Chi2 = " << chi2WithOneHitForRec << ")" << endl;
502             if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
503           }
504         } if (!bestHitForRec2 && chi2WithOneHitForRec < bestChi2WithOneHitForRec) {
505           // keep track of the best single hitForRec except if a couple
506           // of hits has already been found (i.e. bestHitForRec2!=0x0)
507           bestChi2WithOneHitForRec = chi2WithOneHitForRec;
508           bestHitForRec1 = hitForRecCh2;
509         }
510       }
511     }
512     // reset the extrapolated track parameter for next step
513     trackCandidate->SetExtrapTrackParam(&extrapTrackParamSave);
514   }
515   //
516   // look for candidates in chamber 1 not already attached to a track
517   // if we want to keep all possible tracks or if no good couple of hitForRec has been found
518   if (fgkTrackAllTracks || !bestHitForRec2) {
519     if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
520       cout << "FollowTrackInStation: look for single hits in chamber(1..): " << ch1+1 << endl;
521     }
522     for (Int_t hit1 = 0; hit1 < fNHitsForRecPerChamber[ch1]; hit1++) {
523       hitForRecCh1 = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[ch1]+hit1);
524       if (hitForRecCh1->GetNTrackHits() >= 1) continue; // Skip hitForRec already used
525       chi2WithOneHitForRec = trackCandidate->TryOneHitForRec(hitForRecCh1);
526       // if good chi2 then create a new track by adding the good hitForRec in "ch1" to the "trackCandidate"
527       // We do not try to attach a hitForRec in the other chamber too since it has already been done above
528       if (chi2WithOneHitForRec < maxChi2WithOneHitForRec) {
529         if (fgkTrackAllTracks) {
530           // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new hitForRec's
531           newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(*trackCandidate);
532           fNRecTracks++;
533           AliMUONTrackParam TrackParam1(extrapTrackParamSave);
534           AliMUONTrackExtrap::ExtrapToZ(&TrackParam1, hitForRecCh1->GetZ());
535           newTrack->AddTrackParamAtHit(&TrackParam1,hitForRecCh1);
536           // Sort TrackParamAtHit according to increasing Z
537           newTrack->GetTrackParamAtHit()->Sort();
538           // Update the chi2 of the new track
539           if (newTrack->GetFitFMin()<0) newTrack->SetFitFMin(chi2WithOneHitForRec);
540           else newTrack->SetFitFMin(newTrack->GetFitFMin() + chi2WithOneHitForRec);
541           // Printout for debuging
542           if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
543             cout << "FollowTrackInStation: added one hit in chamber(1..): " << ch1+1
544                  << " (Chi2 = " << chi2WithOneHitForRec << ")" << endl;
545             if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
546           }
547         } if (!bestHitForRec2 && chi2WithOneHitForRec < bestChi2WithOneHitForRec) {
548           // keep track of the best single hitForRec except if a couple
549           // of hits has already been found (i.e. bestHitForRec1!=0x0)
550           bestChi2WithOneHitForRec = chi2WithOneHitForRec;
551           bestHitForRec1 = hitForRecCh1;
552         }
553       }
554     }
555   }
556   //
557   // fill out the best track if required else clean up the fRecTracksPtr array
558   if (!fgkTrackAllTracks && bestHitForRec1) {
559     AliMUONTrackParam TrackParam1(extrapTrackParamSave);
560     AliMUONTrackExtrap::ExtrapToZ(&TrackParam1, bestHitForRec1->GetZ());
561     trackCandidate->AddTrackParamAtHit(&TrackParam1,bestHitForRec1);
562     if (bestHitForRec2) {
563       AliMUONTrackParam TrackParam2(extrapTrackParamSave);
564       AliMUONTrackExtrap::ExtrapToZ(&TrackParam2, bestHitForRec2->GetZ());
565       trackCandidate->AddTrackParamAtHit(&TrackParam2,bestHitForRec2);
566       // Sort TrackParamAtHit according to increasing Z
567       trackCandidate->GetTrackParamAtHit()->Sort();
568       // Update the chi2 of the new track
569       if (trackCandidate->GetFitFMin()<0) trackCandidate->SetFitFMin(bestChi2WithTwoHitForRec);
570       else trackCandidate->SetFitFMin(trackCandidate->GetFitFMin() + bestChi2WithTwoHitForRec);
571       // Printout for debuging
572       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
573         cout << "FollowTrackInStation: added the two best hits in station(1..): " << nextStation+1
574              << " (Chi2 = " << bestChi2WithTwoHitForRec << ")" << endl;
575         if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
576       }
577     } else {
578       // Sort TrackParamAtHit according to increasing Z
579       trackCandidate->GetTrackParamAtHit()->Sort();
580       // Update the chi2 of the new track
581       if (trackCandidate->GetFitFMin()<0) trackCandidate->SetFitFMin(bestChi2WithOneHitForRec);
582       else trackCandidate->SetFitFMin(trackCandidate->GetFitFMin() + bestChi2WithOneHitForRec);
583       // Printout for debuging
584       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
585         cout << "FollowTrackInStation: added the best hit in station(1..): " << nextStation+1
586              << " (Chi2 = " << bestChi2WithOneHitForRec << ")" << endl;
587         if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
588       }
589     }
590   } else {
591     fRecTracksPtr->Remove(trackCandidate); // obsolete track
592     fNRecTracks--;
593   }
594   
595 }
596
597   //__________________________________________________________________________
598 void AliMUONTrackReconstructor::SetVertexForFit(AliMUONTrack* trackCandidate)
599 {
600   /// Add the vertex as a measured hit to constrain the fit of the "trackCandidate"
601   /// Compute the vertex resolution from natural vertex dispersion and
602   /// multiple scattering effets according to trackCandidate path in absorber
603   /// It is necessary to account for multiple scattering effects here instead of during the fit of
604   /// the "trackCandidate" to do not influence the result by changing track resolution at vertex
605   AliDebug(1,"Enter SetVertexForFit");
606   
607   Double_t nonBendingReso2 = fNonBendingVertexDispersion * fNonBendingVertexDispersion;
608   Double_t bendingReso2 = fBendingVertexDispersion * fBendingVertexDispersion;
609   // add multiple scattering effets
610   AliMUONTrackParam paramAtVertex(*((AliMUONTrackParam*)(trackCandidate->GetTrackParamAtHit()->First())));
611   paramAtVertex.DeleteCovariances(); // to be sure to account only for multiple scattering
612   AliMUONTrackExtrap::ExtrapToVertexUncorrected(&paramAtVertex,0.);
613   TMatrixD* paramCov = paramAtVertex.GetCovariances();
614   nonBendingReso2 += (*paramCov)(0,0);
615   bendingReso2 += (*paramCov)(2,2);
616   // Set the vertex
617   AliMUONHitForRec vertex; // Coordinates set to (0.,0.,0.) by default
618   vertex.SetNonBendingReso2(nonBendingReso2);
619   vertex.SetBendingReso2(bendingReso2);
620   trackCandidate->SetVertex(&vertex);
621 }
622
623   //__________________________________________________________________________
624 void AliMUONTrackReconstructor::Fit(AliMUONTrack *track, Bool_t includeMCS, Bool_t calcCov)
625 {
626   /// Fit the track "track" with or without multiple Coulomb scattering according to "includeMCS".
627   
628   Double_t benC, errorParam, invBenP, nonBenC, x, y;
629   AliMUONTrackParam *trackParam;
630   Double_t arg[1], fedm, errdef, fitFMin;
631   Int_t npari, nparx;
632   Int_t status, covStatus;
633   
634   // Clear MINUIT parameters
635   gMinuit->mncler();
636   // Give the fitted track to MINUIT
637   gMinuit->SetObjectFit(track);
638   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
639     // Define print level
640     arg[0] = 1;
641     gMinuit->mnexcm("SET PRI", arg, 1, status);
642     // Print covariance matrix
643     gMinuit->mnexcm("SHO COV", arg, 0, status);
644   } else {
645     arg[0] = -1;
646     gMinuit->mnexcm("SET PRI", arg, 1, status);
647   }
648   // No warnings
649   gMinuit->mnexcm("SET NOW", arg, 0, status);
650   // Define strategy
651   //arg[0] = 2;
652   //gMinuit->mnexcm("SET STR", arg, 1, status);
653   
654   // Switch between available FCN according to "includeMCS"
655   if (includeMCS) gMinuit->SetFCN(TrackChi2MCS);
656   else gMinuit->SetFCN(TrackChi2);
657   
658   // Set fitted parameters (!! The order is very important for the covariance matrix !!)
659   trackParam = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First());
660   // could be tried with no limits for the search (min=max=0) ????
661   // mandatory limits in non Bending to avoid NaN values of parameters
662   gMinuit->mnparm(0, "X", trackParam->GetNonBendingCoor(), 0.03, -500.0, 500.0, status);
663   gMinuit->mnparm(1, "NonBenS", trackParam->GetNonBendingSlope(), 0.001, -0.5, 0.5, status);
664   // mandatory limits in Bending to avoid NaN values of parameters
665   gMinuit->mnparm(2, "Y", trackParam->GetBendingCoor(), 0.10, -500.0, 500.0, status);
666   gMinuit->mnparm(3, "BenS", trackParam->GetBendingSlope(), 0.001, -0.5, 0.5, status);
667   gMinuit->mnparm(4, "InvBenP", trackParam->GetInverseBendingMomentum(), 0.003, -0.4, 0.4, status);
668   
669   // minimization
670   gMinuit->mnexcm("MIGRAD", arg, 0, status);
671   
672   // Calculate the covariance matrix more accurately if required
673   if (calcCov) gMinuit->mnexcm("HESSE", arg, 0, status);
674   
675   // get results into "invBenP", "benC", "nonBenC" ("x", "y")
676   gMinuit->GetParameter(0, x, errorParam);
677   trackParam->SetNonBendingCoor(x);
678   gMinuit->GetParameter(1, nonBenC, errorParam);
679   trackParam->SetNonBendingSlope(nonBenC);
680   gMinuit->GetParameter(2, y, errorParam);
681   trackParam->SetBendingCoor(y);
682   gMinuit->GetParameter(3, benC, errorParam);
683   trackParam->SetBendingSlope(benC);
684   gMinuit->GetParameter(4, invBenP, errorParam);
685   trackParam->SetInverseBendingMomentum(invBenP);
686   
687   // global result of the fit
688   gMinuit->mnstat(fitFMin, fedm, errdef, npari, nparx, covStatus);
689   track->SetFitFMin(fitFMin);
690   
691   // Get the covariance matrix if required
692   if (calcCov) {
693     // Covariance matrix according to HESSE status
694     // If problem then keep only the diagonal terms (variances)
695     Double_t matrix[5][5];
696     gMinuit->mnemat(&matrix[0][0],5);
697     if (covStatus == 3) trackParam->SetCovariances(matrix);
698     else trackParam->SetVariances(matrix);
699   } else *(trackParam->GetCovariances()) = 0;
700   
701 }
702
703   //__________________________________________________________________________
704 void TrackChi2(Int_t & /*NParam*/, Double_t * /*Gradient*/, Double_t &Chi2, Double_t *Param, Int_t /*Flag*/)
705 {
706   /// Return the "Chi2" to be minimized with Minuit for track fitting,
707   /// with "NParam" parameters
708   /// and their current values in array pointed to by "Param".
709   /// Assumes that the track hits are sorted according to increasing Z.
710   /// Track parameters at each TrackHit are updated accordingly.
711   /// Multiple Coulomb scattering is not taken into account
712   
713   AliMUONTrack *trackBeingFitted = (AliMUONTrack*) gMinuit->GetObjectFit();
714 //  AliMUONTrack *trackBeingFitted = (AliMUONTrack*) AliMUONTrackReconstructor::Fitter()->GetObjectFit();
715   AliMUONTrackParam param1;
716   AliMUONTrackParam* trackParamAtHit;
717   AliMUONHitForRec* hitForRec;
718   Chi2 = 0.0; // initialize Chi2
719   Double_t dX, dY;
720   
721   // copy of track parameters to be fitted
722   param1 = *((AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->First()));
723   param1.SetNonBendingCoor(Param[0]);
724   param1.SetNonBendingSlope(Param[1]);
725   param1.SetBendingCoor(Param[2]);
726   param1.SetBendingSlope(Param[3]);
727   param1.SetInverseBendingMomentum(Param[4]);
728   
729   // Take the vertex into account in the fit if required
730   if (trackBeingFitted->GetFitWithVertex()) {
731     AliMUONTrackParam paramAtVertex(param1);
732     AliMUONTrackExtrap::ExtrapToZ(&paramAtVertex, 0.);
733     AliMUONHitForRec *vertex = trackBeingFitted->GetVertex();
734     if (!vertex) {
735       cout<<"Error in TrackChi2: Want to use the vertex in tracking but it has not been created!!"<<endl;
736       exit(-1);
737     }
738     dX = vertex->GetNonBendingCoor() - paramAtVertex.GetNonBendingCoor();
739     dY = vertex->GetBendingCoor() - paramAtVertex.GetBendingCoor();
740     Chi2 += dX * dX / vertex->GetNonBendingReso2() + dY * dY / vertex->GetBendingReso2();
741   }
742   
743   // Follow track through all planes of track hits
744   trackParamAtHit = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->First());
745   while (trackParamAtHit) {
746     hitForRec = trackParamAtHit->GetHitForRecPtr();
747     // extrapolation to the plane of the hitForRec attached to the current trackParamAtHit
748     AliMUONTrackExtrap::ExtrapToZ(&param1, hitForRec->GetZ());
749     // update track parameters of the current hit
750     trackParamAtHit->SetTrackParam(param1);
751     // Increment Chi2
752     // done hit per hit, with hit resolution,
753     // and not with point and angle like in "reco_muon.F" !!!!
754     dX = hitForRec->GetNonBendingCoor() - param1.GetNonBendingCoor();
755     dY = hitForRec->GetBendingCoor() - param1.GetBendingCoor();
756     Chi2 = Chi2 + dX * dX / hitForRec->GetNonBendingReso2() + dY * dY / hitForRec->GetBendingReso2();
757     trackParamAtHit = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->After(trackParamAtHit));
758   }
759 }
760
761   //__________________________________________________________________________
762 void TrackChi2MCS(Int_t & /*NParam*/, Double_t * /*Gradient*/, Double_t &Chi2, Double_t *Param, Int_t /*Flag*/)
763 {
764   /// Return the "Chi2" to be minimized with Minuit for track fitting,
765   /// with "NParam" parameters
766   /// and their current values in array pointed to by "Param".
767   /// Assumes that the track hits are sorted according to increasing Z.
768   /// Track parameters at each TrackHit are updated accordingly.
769   /// Multiple Coulomb scattering is taken into account with covariance matrix.
770   
771   AliMUONTrack *trackBeingFitted = (AliMUONTrack*) gMinuit->GetObjectFit();
772 //  AliMUONTrack *trackBeingFitted = (AliMUONTrack*) AliMUONTrackReconstructor::Fitter()->GetObjectFit();
773   AliMUONTrackParam param1;
774   AliMUONTrackParam* trackParamAtHit;
775   AliMUONHitForRec* hitForRec;
776   Chi2 = 0.0; // initialize Chi2
777   Int_t chCurrent, chPrev = 0, hitNumber, hitNumber1, hitNumber2, hitNumber3;
778   Double_t z1, z2, z3;
779   AliMUONTrackParam *trackParamAtHit1, *trackParamAtHit2, *trackParamAtHit3;
780   AliMUONHitForRec *hitForRec1, *hitForRec2;
781   Double_t hbc1, hbc2, pbc1, pbc2;
782   Double_t hnbc1, hnbc2, pnbc1, pnbc2;
783   Int_t numberOfHit = trackBeingFitted->GetNTrackHits();
784   TMatrixD *covBending = new TMatrixD(numberOfHit, numberOfHit);
785   TMatrixD *covNonBending = new TMatrixD(numberOfHit, numberOfHit);
786   Double_t *msa2 = new Double_t[numberOfHit];
787   
788   // copy of track parameters to be fitted
789   param1 = *((AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->First()));
790   param1.SetNonBendingCoor(Param[0]);
791   param1.SetNonBendingSlope(Param[1]);
792   param1.SetBendingCoor(Param[2]);
793   param1.SetBendingSlope(Param[3]);
794   param1.SetInverseBendingMomentum(Param[4]);
795
796   // Take the vertex into account in the fit if required
797   if (trackBeingFitted->GetFitWithVertex()) {
798     AliMUONTrackParam paramAtVertex(param1);
799     AliMUONTrackExtrap::ExtrapToZ(&paramAtVertex, 0.);
800     AliMUONHitForRec *vertex = trackBeingFitted->GetVertex();
801     if (!vertex) {
802       cout<<"Error in TrackChi2MCS: Want to use the vertex in tracking but it has not been created!!"<<endl;
803       exit(-1);
804     }
805     Double_t dX = vertex->GetNonBendingCoor() - paramAtVertex.GetNonBendingCoor();
806     Double_t dY = vertex->GetBendingCoor() - paramAtVertex.GetBendingCoor();
807     Chi2 += dX * dX / vertex->GetNonBendingReso2() + dY * dY / vertex->GetBendingReso2();
808   }
809   
810   // Predicted coordinates and multiple scattering angles are first calculated
811   for (hitNumber = 0; hitNumber < numberOfHit; hitNumber++) {
812     trackParamAtHit = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber));
813     hitForRec = trackParamAtHit->GetHitForRecPtr();
814     // extrapolation to the plane of the hitForRec attached to the current trackParamAtHit
815     AliMUONTrackExtrap::ExtrapToZ(&param1, hitForRec->GetZ());
816     // update track parameters of the current hit
817     trackParamAtHit->SetTrackParam(param1);
818     // square of multiple scattering angle at current hit, with one chamber
819     msa2[hitNumber] = MultipleScatteringAngle2(&param1);
820     // correction for eventual missing hits or multiple hits in a chamber,
821     // according to the number of chambers
822     // between the current hit and the previous one
823     chCurrent = hitForRec->GetChamberNumber();
824     if (hitNumber > 0) msa2[hitNumber] = msa2[hitNumber] * (chCurrent - chPrev);
825     chPrev = chCurrent;
826   }
827
828   // Calculates the covariance matrix
829   for (hitNumber1 = 0; hitNumber1 < numberOfHit; hitNumber1++) { 
830     trackParamAtHit1 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber1));
831     hitForRec1 = trackParamAtHit1->GetHitForRecPtr();
832     z1 = hitForRec1->GetZ();
833     for (hitNumber2 = hitNumber1; hitNumber2 < numberOfHit; hitNumber2++) {
834       trackParamAtHit2 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber2));
835       z2 = trackParamAtHit2->GetHitForRecPtr()->GetZ();
836       // initialization to 0 (diagonal plus upper triangular part)
837       (*covBending)(hitNumber2, hitNumber1) = 0.0;
838       // contribution from multiple scattering in bending plane:
839       // loop over upstream hits
840       for (hitNumber3 = 0; hitNumber3 < hitNumber1; hitNumber3++) {     
841         trackParamAtHit3 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber3));
842         z3 = trackParamAtHit3->GetHitForRecPtr()->GetZ();
843         (*covBending)(hitNumber2, hitNumber1) = (*covBending)(hitNumber2, hitNumber1) + ((z1 - z3) * (z2 - z3) * msa2[hitNumber3]); 
844       }
845       // equal contribution from multiple scattering in non bending plane
846       (*covNonBending)(hitNumber2, hitNumber1) = (*covBending)(hitNumber2, hitNumber1);
847       if (hitNumber1 == hitNumber2) {
848         // Diagonal elements: add contribution from position measurements
849         // in bending plane
850         (*covBending)(hitNumber2, hitNumber1) = (*covBending)(hitNumber2, hitNumber1) + hitForRec1->GetBendingReso2();
851         // and in non bending plane
852         (*covNonBending)(hitNumber2, hitNumber1) = (*covNonBending)(hitNumber2, hitNumber1) + hitForRec1->GetNonBendingReso2();
853       } else {
854         // Non diagonal elements: symmetrization
855         // for bending plane
856         (*covBending)(hitNumber1, hitNumber2) = (*covBending)(hitNumber2, hitNumber1);
857         // and non bending plane
858         (*covNonBending)(hitNumber1, hitNumber2) = (*covNonBending)(hitNumber2, hitNumber1);
859       }
860     } // for (hitNumber2 = hitNumber1;...
861   } // for (hitNumber1 = 0;...
862     
863   // Inversion of covariance matrices
864   Int_t ifailBending;
865   gMinuit->mnvert(&((*covBending)(0,0)), numberOfHit, numberOfHit, numberOfHit, ifailBending);
866   Int_t ifailNonBending;
867   gMinuit->mnvert(&((*covNonBending)(0,0)), numberOfHit, numberOfHit, numberOfHit, ifailNonBending);
868
869   // It would be worth trying to calculate the inverse of the covariance matrix
870   // only once per fit, since it cannot change much in principle,
871   // and it would save a lot of computing time !!!!
872   
873   // Calculates Chi2
874   if ((ifailBending == 0) && (ifailNonBending == 0)) {
875     // with Multiple Scattering if inversion correct
876     for (hitNumber1 = 0; hitNumber1 < numberOfHit ; hitNumber1++) { 
877       trackParamAtHit1 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber1));
878       hitForRec1 = trackParamAtHit1->GetHitForRecPtr();
879       hbc1 = hitForRec1->GetBendingCoor();
880       pbc1 = trackParamAtHit1->GetBendingCoor();
881       hnbc1 = hitForRec1->GetNonBendingCoor();
882       pnbc1 = trackParamAtHit1->GetNonBendingCoor();
883       for (hitNumber2 = 0; hitNumber2 < numberOfHit; hitNumber2++) {
884         trackParamAtHit2 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber2));
885         hitForRec2 = trackParamAtHit2->GetHitForRecPtr();
886         hbc2 = hitForRec2->GetBendingCoor();
887         pbc2 = trackParamAtHit2->GetBendingCoor();
888         hnbc2 = hitForRec2->GetNonBendingCoor();
889         pnbc2 = trackParamAtHit2->GetNonBendingCoor();
890         Chi2 += ((*covBending)(hitNumber2, hitNumber1) * (hbc1 - pbc1) * (hbc2 - pbc2)) +
891                 ((*covNonBending)(hitNumber2, hitNumber1) * (hnbc1 - pnbc1) * (hnbc2 - pnbc2));
892       }
893     }
894   } else {
895     // without Multiple Scattering if inversion impossible
896     for (hitNumber1 = 0; hitNumber1 < numberOfHit ; hitNumber1++) { 
897       trackParamAtHit1 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber1));
898       hitForRec1 = trackParamAtHit1->GetHitForRecPtr();
899       hbc1 = hitForRec1->GetBendingCoor();
900       pbc1 = trackParamAtHit1->GetBendingCoor();
901       hnbc1 = hitForRec1->GetNonBendingCoor();
902       pnbc1 = trackParamAtHit1->GetNonBendingCoor();
903       Chi2 += ((hbc1 - pbc1) * (hbc1 - pbc1) / hitForRec1->GetBendingReso2()) +
904               ((hnbc1 - pnbc1) * (hnbc1 - pnbc1) / hitForRec1->GetNonBendingReso2());
905     }
906   }
907   
908   delete covBending;
909   delete covNonBending;
910   delete [] msa2;
911 }
912
913   //__________________________________________________________________________
914 Double_t MultipleScatteringAngle2(AliMUONTrackParam *param)
915 {
916   /// Returns square of multiple Coulomb scattering angle
917   /// from TrackParamAtHit pointed to by "param"
918   Double_t slopeBending, slopeNonBending, radiationLength, inverseBendingMomentum2, inverseTotalMomentum2;
919   Double_t varMultipleScatteringAngle;
920   slopeBending = param->GetBendingSlope();
921   slopeNonBending = param->GetNonBendingSlope();
922   // thickness in radiation length for the current track,
923   // taking local angle into account
924   radiationLength = AliMUONConstants::ChamberThicknessInX0() *
925                     TMath::Sqrt(1.0 + slopeBending*slopeBending + slopeNonBending*slopeNonBending);
926   inverseBendingMomentum2 =  param->GetInverseBendingMomentum() * param->GetInverseBendingMomentum();
927   inverseTotalMomentum2 = inverseBendingMomentum2 * (1.0 + slopeBending * slopeBending) /
928                           (1.0 + slopeBending *slopeBending + slopeNonBending * slopeNonBending); 
929   varMultipleScatteringAngle = 0.0136 * (1.0 + 0.038 * TMath::Log(radiationLength));
930   // The velocity is assumed to be 1 !!!!
931   varMultipleScatteringAngle = inverseTotalMomentum2 * radiationLength * varMultipleScatteringAngle * varMultipleScatteringAngle;
932   return varMultipleScatteringAngle;
933 }
934
935   //__________________________________________________________________________
936 void AliMUONTrackReconstructor::ExtrapTracksToVertex()
937 {
938   /// propagate tracks to the vertex through the absorber (Branson)
939   AliMUONTrack *track;
940   AliMUONTrackParam trackParamVertex;
941   AliMUONHitForRec *vertex;
942   Double_t vX, vY, vZ;
943   
944   for (Int_t iRecTrack = 0; iRecTrack <fNRecTracks; iRecTrack++) {
945     track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iRecTrack);
946     trackParamVertex = *((AliMUONTrackParam*)(track->GetTrackParamAtHit()->First()));
947     vertex = track->GetVertex();
948     if (vertex) { // if track as a vertex defined, use it
949       vX = vertex->GetNonBendingCoor();
950       vY = vertex->GetBendingCoor();
951       vZ = vertex->GetZ();
952     } else { // else vertex = (0.,0.,0.)
953       vX = 0.;
954       vY = 0.;
955       vZ = 0.;
956     }
957     AliMUONTrackExtrap::ExtrapToVertex(&trackParamVertex, vX, vY, vZ);
958     track->SetTrackParamAtVertex(&trackParamVertex);
959     
960     if (AliLog::GetGlobalDebugLevel() > 0) {
961       cout << "FollowTracks: track candidate(0..): " << iRecTrack
962            << " after extrapolation to vertex" << endl;
963       track->RecursiveDump();
964     }
965   }
966   
967 }
968
969   //__________________________________________________________________________
970 void AliMUONTrackReconstructor::FillMUONTrack()
971 {
972   /// Fill fHitForRecAtHit of AliMUONTrack's
973   AliMUONTrack *track;
974   AliMUONTrackParam *trackParamAtHit;
975   track = (AliMUONTrack*) fRecTracksPtr->First();
976   while (track) {
977     trackParamAtHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First());
978     while (trackParamAtHit) {
979       track->AddHitForRecAtHit(trackParamAtHit->GetHitForRecPtr());
980       trackParamAtHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->After(trackParamAtHit)); 
981     }
982     track = (AliMUONTrack*) fRecTracksPtr->After(track);
983   }
984   return;
985 }
986
987   //__________________________________________________________________________
988 void AliMUONTrackReconstructor::EventDump(void)
989 {
990   /// Dump reconstructed event (track parameters at vertex and at first hit),
991   /// and the particle parameters
992   AliMUONTrack *track;
993   AliMUONTrackParam *trackParam, *trackParam1;
994   Double_t bendingSlope, nonBendingSlope, pYZ;
995   Double_t pX, pY, pZ, x, y, z, c;
996   Int_t trackIndex, nTrackHits;
997  
998   AliDebug(1,"****** enter EventDump ******");
999   AliDebug(1, Form("Number of Reconstructed tracks : %d", fNRecTracks)); 
1000   
1001   fRecTracksPtr->Compress(); // for simple loop without "Next" since no hole
1002   // Loop over reconstructed tracks
1003   for (trackIndex = 0; trackIndex < fNRecTracks; trackIndex++) {
1004     AliDebug(1, Form("track number: %d", trackIndex));
1005     // function for each track for modularity ????
1006     track = (AliMUONTrack*) ((*fRecTracksPtr)[trackIndex]);
1007     nTrackHits = track->GetNTrackHits();
1008     AliDebug(1, Form("Number of track hits: %d ", nTrackHits));
1009     // track parameters at Vertex
1010     trackParam = track->GetTrackParamAtVertex();
1011     x = trackParam->GetNonBendingCoor();
1012     y = trackParam->GetBendingCoor();
1013     z = trackParam->GetZ();
1014     bendingSlope = trackParam->GetBendingSlope();
1015     nonBendingSlope = trackParam->GetNonBendingSlope();
1016     pYZ = 1/TMath::Abs(trackParam->GetInverseBendingMomentum());
1017     pZ = pYZ/TMath::Sqrt(1+bendingSlope*bendingSlope);
1018     pX = pZ * nonBendingSlope;
1019     pY = pZ * bendingSlope;
1020     c = TMath::Sign(1.0, trackParam->GetInverseBendingMomentum());
1021     AliDebug(1, Form("Track parameters at Vertex z= %f: X= %f Y= %f pX= %f pY= %f pZ= %f c= %f\n",
1022                      z, x, y, pX, pY, pZ, c));
1023
1024     // track parameters at first hit
1025     trackParam1 = (AliMUONTrackParam*) track->GetTrackParamAtHit()->First();
1026     x = trackParam1->GetNonBendingCoor();
1027     y = trackParam1->GetBendingCoor();
1028     z = trackParam1->GetZ();
1029     bendingSlope = trackParam1->GetBendingSlope();
1030     nonBendingSlope = trackParam1->GetNonBendingSlope();
1031     pYZ = 1/TMath::Abs(trackParam1->GetInverseBendingMomentum());
1032     pZ = pYZ/TMath::Sqrt(1.0 + bendingSlope * bendingSlope);
1033     pX = pZ * nonBendingSlope;
1034     pY = pZ * bendingSlope;
1035     c = TMath::Sign(1.0, trackParam1->GetInverseBendingMomentum());
1036     AliDebug(1, Form("track parameters at z= %f: X= %f Y= %f pX= %f pY= %f pZ= %f c= %f\n",
1037                      z, x, y, pX, pY, pZ, c));
1038   }
1039   // informations about generated particles NO !!!!!!!!
1040   
1041 //    for (Int_t iPart = 0; iPart < np; iPart++) {
1042 //      p = gAlice->Particle(iPart);
1043 //      printf(" particle %d: type= %d px= %f py= %f pz= %f pdg= %d\n",
1044 //         iPart, p->GetPdgCode(), p->Px(), p->Py(), p->Pz(), p->GetPdgCode());    
1045 //    }
1046   return;
1047 }
1048
1049