MUON/AliMUONTrackReconstructor.cxx

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