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