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