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   Bool_t *hitForRecCh1Used = new Bool_t[fNHitsForRecPerChamber[ch1]];
 419   for (Int_t hit1 = 0; hit1 < fNHitsForRecPerChamber[ch1]; hit1++) hitForRecCh1Used[hit1] = kFALSE;
 420   //
 421   //Extrapolate trackCandidate to chamber "ch2" to save computing time in the next steps
 422   AliMUONTrackParam *extrapTrackParamPtr = trackCandidate->GetExtrapTrackParam();
 423   *extrapTrackParamPtr = *((AliMUONTrackParam*)(trackCandidate->GetTrackParamAtHit()->First()));
 424   AliMUONTrackExtrap::ExtrapToZCov(extrapTrackParamPtr, zCh2);
 425   AliMUONTrackParam extrapTrackParamSave(*extrapTrackParamPtr);
 426   //
 427   // Printout for debuging
 428   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
 429     TMatrixD* paramCovForDebug = extrapTrackParamPtr->GetCovariances();
 430     cout<<endl<<"Track parameter covariances at first hit extrapolated to z = "<<zCh2<<":"<<endl;
 431     paramCovForDebug->Print();
 432   }
 433   //
 434   // look for candidates in chamber 2
 435   // Printout for debuging
 436   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 437     cout << "FollowTrackInStation: look for hits in chamber(1..): " << ch2+1 << endl;
 438   }
 439   for (Int_t hit2 = 0; hit2 < fNHitsForRecPerChamber[ch2]; hit2++) {
 440     hitForRecCh2 = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[ch2]+hit2);
 441     // extrapolate track parameters and covariances only once for this hit
 442     AliMUONTrackExtrap::ExtrapToZCov(extrapTrackParamPtr, hitForRecCh2->GetZ());
 443     chi2WithOneHitForRec = trackCandidate->TryOneHitForRec(hitForRecCh2);
 444     // if good chi2 then try to attach a hitForRec in the other chamber too
 445     if (chi2WithOneHitForRec < maxChi2WithOneHitForRec) {
 446       // Printout for debuging
 447       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 448         cout << "FollowTrackInStation: look for second hits in chamber(1..): " << ch1+1 << endl;
 449       }
 450       Bool_t foundSecondHit = kFALSE;
 451       for (Int_t hit1 = 0; hit1 < fNHitsForRecPerChamber[ch1]; hit1++) {
 452         hitForRecCh1 = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[ch1]+hit1);
 453         chi2WithTwoHitForRec = trackCandidate->TryTwoHitForRec(hitForRecCh2, hitForRecCh1); // order hits like that to save computing time
 454         // if good chi2 then create a new track by adding the 2 hitForRec to the "trackCandidate"
 455         if (chi2WithTwoHitForRec < maxChi2WithTwoHitForRec) {
 456           foundSecondHit = kTRUE;
 457           if (fgkTrackAllTracks) {
 458             // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new hitForRec's
 459             newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(*trackCandidate);
 460             fNRecTracks++;
 461             AliMUONTrackParam TrackParam1(extrapTrackParamSave);
 462             AliMUONTrackExtrap::ExtrapToZ(&TrackParam1, hitForRecCh1->GetZ());
 463             newTrack->AddTrackParamAtHit(&TrackParam1,hitForRecCh1);
 464             AliMUONTrackParam TrackParam2(extrapTrackParamSave);
 465             AliMUONTrackExtrap::ExtrapToZ(&TrackParam2, hitForRecCh2->GetZ());
 466             newTrack->AddTrackParamAtHit(&TrackParam2,hitForRecCh2);
 467             // Sort TrackParamAtHit according to increasing Z
 468             newTrack->GetTrackParamAtHit()->Sort();
 469             // Update the chi2 of the new track
 470             if (newTrack->GetFitFMin()<0) newTrack->SetFitFMin(chi2WithTwoHitForRec);
 471             else newTrack->SetFitFMin(newTrack->GetFitFMin() + chi2WithTwoHitForRec);
 472             // Tag hitForRecCh1 as used
 473             hitForRecCh1Used[hit1] = kTRUE;
 474             // Printout for debuging
 475             if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 476               cout << "FollowTrackInStation: added two hits in station(1..): " << nextStation+1
 477                    << " (Chi2 = " << chi2WithTwoHitForRec << ")" << endl;
 478               if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
 479             }
 480           } else if (chi2WithTwoHitForRec < bestChi2WithTwoHitForRec) {
 481             // keep track of the best couple of hits
 482             bestChi2WithTwoHitForRec = chi2WithTwoHitForRec;
 483             bestHitForRec1 = hitForRecCh1;
 484             bestHitForRec2 = hitForRecCh2;
 485           }
 486         }
 487       }
 488       // if no hitForRecCh1 found then consider to add hitForRecCh2 only
 489       if (!foundSecondHit) {
 490         if (fgkTrackAllTracks) {
 491           // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new hitForRec's
 492           newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(*trackCandidate);
 493           fNRecTracks++;
 494           AliMUONTrackParam TrackParam1(extrapTrackParamSave);
 495           AliMUONTrackExtrap::ExtrapToZ(&TrackParam1, hitForRecCh2->GetZ());
 496           newTrack->AddTrackParamAtHit(&TrackParam1,hitForRecCh2);
 497           // Sort TrackParamAtHit according to increasing Z
 498           newTrack->GetTrackParamAtHit()->Sort();
 499           // Update the chi2 of the new track
 500           if (newTrack->GetFitFMin()<0) newTrack->SetFitFMin(chi2WithOneHitForRec);
 501           else newTrack->SetFitFMin(newTrack->GetFitFMin() + chi2WithOneHitForRec);
 502           // Printout for debuging
 503           if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 504             cout << "FollowTrackInStation: added one hit in chamber(1..): " << ch2+1
 505                  << " (Chi2 = " << chi2WithOneHitForRec << ")" << endl;
 506             if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
 507           }
 508         } if (!bestHitForRec2 && chi2WithOneHitForRec < bestChi2WithOneHitForRec) {
 509           // keep track of the best single hitForRec except if a couple
 510           // of hits has already been found (i.e. bestHitForRec2!=0x0)
 511           bestChi2WithOneHitForRec = chi2WithOneHitForRec;
 512           bestHitForRec1 = hitForRecCh2;
 513         }
 514       }
 515     }
 516     // reset the extrapolated track parameter for next step
 517     trackCandidate->SetExtrapTrackParam(&extrapTrackParamSave);
 518   }
 519   //
 520   // look for candidates in chamber 1 not already attached to a track
 521   // if we want to keep all possible tracks or if no good couple of hitForRec has been found
 522   if (fgkTrackAllTracks || !bestHitForRec2) {
 523     if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 524       cout << "FollowTrackInStation: look for single hits in chamber(1..): " << ch1+1 << endl;
 525     }
 526     for (Int_t hit1 = 0; hit1 < fNHitsForRecPerChamber[ch1]; hit1++) {
 527       hitForRecCh1 = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[ch1]+hit1);
 528       if (hitForRecCh1Used[hit1]) continue; // Skip hitForRec already used
 529       chi2WithOneHitForRec = trackCandidate->TryOneHitForRec(hitForRecCh1);
 530       // if good chi2 then create a new track by adding the good hitForRec in "ch1" to the "trackCandidate"
 531       // We do not try to attach a hitForRec in the other chamber too since it has already been done above
 532       if (chi2WithOneHitForRec < maxChi2WithOneHitForRec) {
 533         if (fgkTrackAllTracks) {
 534           // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new hitForRec's
 535           newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(*trackCandidate);
 536           fNRecTracks++;
 537           AliMUONTrackParam TrackParam1(extrapTrackParamSave);
 538           AliMUONTrackExtrap::ExtrapToZ(&TrackParam1, hitForRecCh1->GetZ());
 539           newTrack->AddTrackParamAtHit(&TrackParam1,hitForRecCh1);
 540           // Sort TrackParamAtHit according to increasing Z
 541           newTrack->GetTrackParamAtHit()->Sort();
 542           // Update the chi2 of the new track
 543           if (newTrack->GetFitFMin()<0) newTrack->SetFitFMin(chi2WithOneHitForRec);
 544           else newTrack->SetFitFMin(newTrack->GetFitFMin() + chi2WithOneHitForRec);
 545           // Printout for debuging
 546           if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 547             cout << "FollowTrackInStation: added one hit in chamber(1..): " << ch1+1
 548                  << " (Chi2 = " << chi2WithOneHitForRec << ")" << endl;
 549             if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
 550           }
 551         } if (!bestHitForRec2 && chi2WithOneHitForRec < bestChi2WithOneHitForRec) {
 552           // keep track of the best single hitForRec except if a couple
 553           // of hits has already been found (i.e. bestHitForRec1!=0x0)
 554           bestChi2WithOneHitForRec = chi2WithOneHitForRec;
 555           bestHitForRec1 = hitForRecCh1;
 556         }
 557       }
 558     }
 559   }
 560   //
 561   // fill out the best track if required else clean up the fRecTracksPtr array
 562   if (!fgkTrackAllTracks && bestHitForRec1) {
 563     AliMUONTrackParam TrackParam1(extrapTrackParamSave);
 564     AliMUONTrackExtrap::ExtrapToZ(&TrackParam1, bestHitForRec1->GetZ());
 565     trackCandidate->AddTrackParamAtHit(&TrackParam1,bestHitForRec1);
 566     if (bestHitForRec2) {
 567       AliMUONTrackParam TrackParam2(extrapTrackParamSave);
 568       AliMUONTrackExtrap::ExtrapToZ(&TrackParam2, bestHitForRec2->GetZ());
 569       trackCandidate->AddTrackParamAtHit(&TrackParam2,bestHitForRec2);
 570       // Sort TrackParamAtHit according to increasing Z
 571       trackCandidate->GetTrackParamAtHit()->Sort();
 572       // Update the chi2 of the new track
 573       if (trackCandidate->GetFitFMin()<0) trackCandidate->SetFitFMin(bestChi2WithTwoHitForRec);
 574       else trackCandidate->SetFitFMin(trackCandidate->GetFitFMin() + bestChi2WithTwoHitForRec);
 575       // Printout for debuging
 576       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 577         cout << "FollowTrackInStation: added the two best hits in station(1..): " << nextStation+1
 578              << " (Chi2 = " << bestChi2WithTwoHitForRec << ")" << endl;
 579         if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
 580       }
 581     } else {
 582       // Sort TrackParamAtHit according to increasing Z
 583       trackCandidate->GetTrackParamAtHit()->Sort();
 584       // Update the chi2 of the new track
 585       if (trackCandidate->GetFitFMin()<0) trackCandidate->SetFitFMin(bestChi2WithOneHitForRec);
 586       else trackCandidate->SetFitFMin(trackCandidate->GetFitFMin() + bestChi2WithOneHitForRec);
 587       // Printout for debuging
 588       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 589         cout << "FollowTrackInStation: added the best hit in station(1..): " << nextStation+1
 590              << " (Chi2 = " << bestChi2WithOneHitForRec << ")" << endl;
 591         if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
 592       }
 593     }
 594   } else {
 595     fRecTracksPtr->Remove(trackCandidate); // obsolete track
 596     fNRecTracks--;
 597   }
 598
 599 }
 600
 601   //__________________________________________________________________________
 602 void AliMUONTrackReconstructor::SetVertexForFit(AliMUONTrack* trackCandidate)
 603 {
 604   /// Add the vertex as a measured hit to constrain the fit of the "trackCandidate"
 605   /// Compute the vertex resolution from natural vertex dispersion and
 606   /// multiple scattering effets according to trackCandidate path in absorber
 607   /// It is necessary to account for multiple scattering effects here instead of during the fit of
 608   /// the "trackCandidate" to do not influence the result by changing track resolution at vertex
 609   AliDebug(1,"Enter SetVertexForFit");
 610
 611   Double_t nonBendingReso2 = fNonBendingVertexDispersion * fNonBendingVertexDispersion;
 612   Double_t bendingReso2 = fBendingVertexDispersion * fBendingVertexDispersion;
 613   // add multiple scattering effets
 614   AliMUONTrackParam paramAtVertex(*((AliMUONTrackParam*)(trackCandidate->GetTrackParamAtHit()->First())));
 615   paramAtVertex.DeleteCovariances(); // to be sure to account only for multiple scattering
 616   AliMUONTrackExtrap::ExtrapToVertexUncorrected(&paramAtVertex,0.);
 617   TMatrixD* paramCov = paramAtVertex.GetCovariances();
 618   nonBendingReso2 += (*paramCov)(0,0);
 619   bendingReso2 += (*paramCov)(2,2);
 620   // Set the vertex
 621   AliMUONHitForRec vertex; // Coordinates set to (0.,0.,0.) by default
 622   vertex.SetNonBendingReso2(nonBendingReso2);
 623   vertex.SetBendingReso2(bendingReso2);
 624   trackCandidate->SetVertex(&vertex);
 625 }
 626
 627   //__________________________________________________________________________
 628 void AliMUONTrackReconstructor::Fit(AliMUONTrack *track, Bool_t includeMCS, Bool_t calcCov)
 629 {
 630   /// Fit the track "track" with or without multiple Coulomb scattering according to "includeMCS".
 631
 632   Double_t benC, errorParam, invBenP, nonBenC, x, y;
 633   AliMUONTrackParam *trackParam;
 634   Double_t arg[1], fedm, errdef, fitFMin;
 635   Int_t npari, nparx;
 636   Int_t status, covStatus;
 637
 638   // Clear MINUIT parameters
 639   gMinuit->mncler();
 640   // Give the fitted track to MINUIT
 641   gMinuit->SetObjectFit(track);
 642   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
 643     // Define print level
 644     arg[0] = 1;
 645     gMinuit->mnexcm("SET PRI", arg, 1, status);
 646     // Print covariance matrix
 647     gMinuit->mnexcm("SHO COV", arg, 0, status);
 648   } else {
 649     arg[0] = -1;
 650     gMinuit->mnexcm("SET PRI", arg, 1, status);
 651   }
 652   // No warnings
 653   gMinuit->mnexcm("SET NOW", arg, 0, status);
 654   // Define strategy
 655   //arg[0] = 2;
 656   //gMinuit->mnexcm("SET STR", arg, 1, status);
 657
 658   // Switch between available FCN according to "includeMCS"
 659   if (includeMCS) gMinuit->SetFCN(TrackChi2MCS);
 660   else gMinuit->SetFCN(TrackChi2);
 661
 662   // Set fitted parameters (!! The order is very important for the covariance matrix !!)
 663   trackParam = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First());
 664   // could be tried with no limits for the search (min=max=0) ????
 665   // mandatory limits in non Bending to avoid NaN values of parameters
 666   gMinuit->mnparm(0, "X", trackParam->GetNonBendingCoor(), 0.03, -500.0, 500.0, status);
 667   gMinuit->mnparm(1, "NonBenS", trackParam->GetNonBendingSlope(), 0.001, -0.5, 0.5, status);
 668   // mandatory limits in Bending to avoid NaN values of parameters
 669   gMinuit->mnparm(2, "Y", trackParam->GetBendingCoor(), 0.10, -500.0, 500.0, status);
 670   gMinuit->mnparm(3, "BenS", trackParam->GetBendingSlope(), 0.001, -0.5, 0.5, status);
 671   gMinuit->mnparm(4, "InvBenP", trackParam->GetInverseBendingMomentum(), 0.003, -0.4, 0.4, status);
 672
 673   // minimization
 674   gMinuit->mnexcm("MIGRAD", arg, 0, status);
 675
 676   // Calculate the covariance matrix more accurately if required
 677   if (calcCov) gMinuit->mnexcm("HESSE", arg, 0, status);
 678
 679   // get results into "invBenP", "benC", "nonBenC" ("x", "y")
 680   gMinuit->GetParameter(0, x, errorParam);
 681   trackParam->SetNonBendingCoor(x);
 682   gMinuit->GetParameter(1, nonBenC, errorParam);
 683   trackParam->SetNonBendingSlope(nonBenC);
 684   gMinuit->GetParameter(2, y, errorParam);
 685   trackParam->SetBendingCoor(y);
 686   gMinuit->GetParameter(3, benC, errorParam);
 687   trackParam->SetBendingSlope(benC);
 688   gMinuit->GetParameter(4, invBenP, errorParam);
 689   trackParam->SetInverseBendingMomentum(invBenP);
 690
 691   // global result of the fit
 692   gMinuit->mnstat(fitFMin, fedm, errdef, npari, nparx, covStatus);
 693   track->SetFitFMin(fitFMin);
 694
 695   // Get the covariance matrix if required
 696   if (calcCov) {
 697     // Covariance matrix according to HESSE status
 698     // If problem then keep only the diagonal terms (variances)
 699     Double_t matrix[5][5];
 700     gMinuit->mnemat(&matrix[0][0],5);
 701     if (covStatus == 3) trackParam->SetCovariances(matrix);
 702     else trackParam->SetVariances(matrix);
 703   } else *(trackParam->GetCovariances()) = 0;
 704
 705 }
 706
 707   //__________________________________________________________________________
 708 void TrackChi2(Int_t & /*NParam*/, Double_t * /*Gradient*/, Double_t &Chi2, Double_t *Param, Int_t /*Flag*/)
 709 {
 710   /// Return the "Chi2" to be minimized with Minuit for track fitting,
 711   /// with "NParam" parameters
 712   /// and their current values in array pointed to by "Param".
 713   /// Assumes that the track hits are sorted according to increasing Z.
 714   /// Track parameters at each TrackHit are updated accordingly.
 715   /// Multiple Coulomb scattering is not taken into account
 716
 717   AliMUONTrack *trackBeingFitted = (AliMUONTrack*) gMinuit->GetObjectFit();
 718 //  AliMUONTrack *trackBeingFitted = (AliMUONTrack*) AliMUONTrackReconstructor::Fitter()->GetObjectFit();
 719   AliMUONTrackParam param1;
 720   AliMUONTrackParam* trackParamAtHit;
 721   AliMUONHitForRec* hitForRec;
 722   Chi2 = 0.0; // initialize Chi2
 723   Double_t dX, dY;
 724
 725   // copy of track parameters to be fitted
 726   param1 = *((AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->First()));
 727   param1.SetNonBendingCoor(Param[0]);
 728   param1.SetNonBendingSlope(Param[1]);
 729   param1.SetBendingCoor(Param[2]);
 730   param1.SetBendingSlope(Param[3]);
 731   param1.SetInverseBendingMomentum(Param[4]);
 732
 733   // Take the vertex into account in the fit if required
 734   if (trackBeingFitted->GetFitWithVertex()) {
 735     AliMUONTrackParam paramAtVertex(param1);
 736     AliMUONTrackExtrap::ExtrapToZ(&paramAtVertex, 0.);
 737     AliMUONHitForRec *vertex = trackBeingFitted->GetVertex();
 738     if (!vertex) {
 739       cout<<"Error in TrackChi2: Want to use the vertex in tracking but it has not been created!!"<<endl;
 740       exit(-1);
 741     }
 742     dX = vertex->GetNonBendingCoor() - paramAtVertex.GetNonBendingCoor();
 743     dY = vertex->GetBendingCoor() - paramAtVertex.GetBendingCoor();
 744     Chi2 += dX * dX / vertex->GetNonBendingReso2() + dY * dY / vertex->GetBendingReso2();
 745   }
 746
 747   // Follow track through all planes of track hits
 748   trackParamAtHit = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->First());
 749   while (trackParamAtHit) {
 750     hitForRec = trackParamAtHit->GetHitForRecPtr();
 751     // extrapolation to the plane of the hitForRec attached to the current trackParamAtHit
 752     AliMUONTrackExtrap::ExtrapToZ(&param1, hitForRec->GetZ());
 753     // update track parameters of the current hit
 754     trackParamAtHit->SetTrackParam(param1);
 755     // Increment Chi2
 756     // done hit per hit, with hit resolution,
 757     // and not with point and angle like in "reco_muon.F" !!!!
 758     dX = hitForRec->GetNonBendingCoor() - param1.GetNonBendingCoor();
 759     dY = hitForRec->GetBendingCoor() - param1.GetBendingCoor();
 760     Chi2 = Chi2 + dX * dX / hitForRec->GetNonBendingReso2() + dY * dY / hitForRec->GetBendingReso2();
 761     trackParamAtHit = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->After(trackParamAtHit));
 762   }
 763 }
 764
 765   //__________________________________________________________________________
 766 void TrackChi2MCS(Int_t & /*NParam*/, Double_t * /*Gradient*/, Double_t &Chi2, Double_t *Param, Int_t /*Flag*/)
 767 {
 768   /// Return the "Chi2" to be minimized with Minuit for track fitting,
 769   /// with "NParam" parameters
 770   /// and their current values in array pointed to by "Param".
 771   /// Assumes that the track hits are sorted according to increasing Z.
 772   /// Track parameters at each TrackHit are updated accordingly.
 773   /// Multiple Coulomb scattering is taken into account with covariance matrix.
 774
 775   AliMUONTrack *trackBeingFitted = (AliMUONTrack*) gMinuit->GetObjectFit();
 776 //  AliMUONTrack *trackBeingFitted = (AliMUONTrack*) AliMUONTrackReconstructor::Fitter()->GetObjectFit();
 777   AliMUONTrackParam param1;
 778   AliMUONTrackParam* trackParamAtHit;
 779   AliMUONHitForRec* hitForRec;
 780   Chi2 = 0.0; // initialize Chi2
 781   Int_t chCurrent, chPrev = 0, hitNumber, hitNumber1, hitNumber2, hitNumber3;
 782   Double_t z1, z2, z3;
 783   AliMUONTrackParam *trackParamAtHit1, *trackParamAtHit2, *trackParamAtHit3;
 784   AliMUONHitForRec *hitForRec1, *hitForRec2;
 785   Double_t hbc1, hbc2, pbc1, pbc2;
 786   Double_t hnbc1, hnbc2, pnbc1, pnbc2;
 787   Int_t numberOfHit = trackBeingFitted->GetNTrackHits();
 788   TMatrixD *covBending = new TMatrixD(numberOfHit, numberOfHit);
 789   TMatrixD *covNonBending = new TMatrixD(numberOfHit, numberOfHit);
 790   Double_t *msa2 = new Double_t[numberOfHit];
 791
 792   // copy of track parameters to be fitted
 793   param1 = *((AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->First()));
 794   param1.SetNonBendingCoor(Param[0]);
 795   param1.SetNonBendingSlope(Param[1]);
 796   param1.SetBendingCoor(Param[2]);
 797   param1.SetBendingSlope(Param[3]);
 798   param1.SetInverseBendingMomentum(Param[4]);
 799
 800   // Take the vertex into account in the fit if required
 801   if (trackBeingFitted->GetFitWithVertex()) {
 802     AliMUONTrackParam paramAtVertex(param1);
 803     AliMUONTrackExtrap::ExtrapToZ(&paramAtVertex, 0.);
 804     AliMUONHitForRec *vertex = trackBeingFitted->GetVertex();
 805     if (!vertex) {
 806       cout<<"Error in TrackChi2MCS: Want to use the vertex in tracking but it has not been created!!"<<endl;
 807       exit(-1);
 808     }
 809     Double_t dX = vertex->GetNonBendingCoor() - paramAtVertex.GetNonBendingCoor();
 810     Double_t dY = vertex->GetBendingCoor() - paramAtVertex.GetBendingCoor();
 811     Chi2 += dX * dX / vertex->GetNonBendingReso2() + dY * dY / vertex->GetBendingReso2();
 812   }
 813
 814   // Predicted coordinates and multiple scattering angles are first calculated
 815   for (hitNumber = 0; hitNumber < numberOfHit; hitNumber++) {
 816     trackParamAtHit = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber));
 817     hitForRec = trackParamAtHit->GetHitForRecPtr();
 818     // extrapolation to the plane of the hitForRec attached to the current trackParamAtHit
 819     AliMUONTrackExtrap::ExtrapToZ(&param1, hitForRec->GetZ());
 820     // update track parameters of the current hit
 821     trackParamAtHit->SetTrackParam(param1);
 822     // square of multiple scattering angle at current hit, with one chamber
 823     msa2[hitNumber] = MultipleScatteringAngle2(&param1);
 824     // correction for eventual missing hits or multiple hits in a chamber,
 825     // according to the number of chambers
 826     // between the current hit and the previous one
 827     chCurrent = hitForRec->GetChamberNumber();
 828     if (hitNumber > 0) msa2[hitNumber] = msa2[hitNumber] * (chCurrent - chPrev);
 829     chPrev = chCurrent;
 830   }
 831
 832   // Calculates the covariance matrix
 833   for (hitNumber1 = 0; hitNumber1 < numberOfHit; hitNumber1++) {
 834     trackParamAtHit1 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber1));
 835     hitForRec1 = trackParamAtHit1->GetHitForRecPtr();
 836     z1 = hitForRec1->GetZ();
 837     for (hitNumber2 = hitNumber1; hitNumber2 < numberOfHit; hitNumber2++) {
 838       trackParamAtHit2 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber2));
 839       z2 = trackParamAtHit2->GetHitForRecPtr()->GetZ();
 840       // initialization to 0 (diagonal plus upper triangular part)
 841       (*covBending)(hitNumber2, hitNumber1) = 0.0;
 842       // contribution from multiple scattering in bending plane:
 843       // loop over upstream hits
 844       for (hitNumber3 = 0; hitNumber3 < hitNumber1; hitNumber3++) {
 845         trackParamAtHit3 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber3));
 846         z3 = trackParamAtHit3->GetHitForRecPtr()->GetZ();
 847         (*covBending)(hitNumber2, hitNumber1) = (*covBending)(hitNumber2, hitNumber1) + ((z1 - z3) * (z2 - z3) * msa2[hitNumber3]);
 848       }
 849       // equal contribution from multiple scattering in non bending plane
 850       (*covNonBending)(hitNumber2, hitNumber1) = (*covBending)(hitNumber2, hitNumber1);
 851       if (hitNumber1 == hitNumber2) {
 852         // Diagonal elements: add contribution from position measurements
 853         // in bending plane
 854         (*covBending)(hitNumber2, hitNumber1) = (*covBending)(hitNumber2, hitNumber1) + hitForRec1->GetBendingReso2();
 855         // and in non bending plane
 856         (*covNonBending)(hitNumber2, hitNumber1) = (*covNonBending)(hitNumber2, hitNumber1) + hitForRec1->GetNonBendingReso2();
 857       } else {
 858         // Non diagonal elements: symmetrization
 859         // for bending plane
 860         (*covBending)(hitNumber1, hitNumber2) = (*covBending)(hitNumber2, hitNumber1);
 861         // and non bending plane
 862         (*covNonBending)(hitNumber1, hitNumber2) = (*covNonBending)(hitNumber2, hitNumber1);
 863       }
 864     } // for (hitNumber2 = hitNumber1;...
 865   } // for (hitNumber1 = 0;...
 866
 867   // Inversion of covariance matrices
 868   Int_t ifailBending;
 869   gMinuit->mnvert(&((*covBending)(0,0)), numberOfHit, numberOfHit, numberOfHit, ifailBending);
 870   Int_t ifailNonBending;
 871   gMinuit->mnvert(&((*covNonBending)(0,0)), numberOfHit, numberOfHit, numberOfHit, ifailNonBending);
 872
 873   // It would be worth trying to calculate the inverse of the covariance matrix
 874   // only once per fit, since it cannot change much in principle,
 875   // and it would save a lot of computing time !!!!
 876
 877   // Calculates Chi2
 878   if ((ifailBending == 0) && (ifailNonBending == 0)) {
 879     // with Multiple Scattering if inversion correct
 880     for (hitNumber1 = 0; hitNumber1 < numberOfHit ; hitNumber1++) {
 881       trackParamAtHit1 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber1));
 882       hitForRec1 = trackParamAtHit1->GetHitForRecPtr();
 883       hbc1 = hitForRec1->GetBendingCoor();
 884       pbc1 = trackParamAtHit1->GetBendingCoor();
 885       hnbc1 = hitForRec1->GetNonBendingCoor();
 886       pnbc1 = trackParamAtHit1->GetNonBendingCoor();
 887       for (hitNumber2 = 0; hitNumber2 < numberOfHit; hitNumber2++) {
 888         trackParamAtHit2 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber2));
 889         hitForRec2 = trackParamAtHit2->GetHitForRecPtr();
 890         hbc2 = hitForRec2->GetBendingCoor();
 891         pbc2 = trackParamAtHit2->GetBendingCoor();
 892         hnbc2 = hitForRec2->GetNonBendingCoor();
 893         pnbc2 = trackParamAtHit2->GetNonBendingCoor();
 894         Chi2 += ((*covBending)(hitNumber2, hitNumber1) * (hbc1 - pbc1) * (hbc2 - pbc2)) +
 895                 ((*covNonBending)(hitNumber2, hitNumber1) * (hnbc1 - pnbc1) * (hnbc2 - pnbc2));
 896       }
 897     }
 898   } else {
 899     // without Multiple Scattering if inversion impossible
 900     for (hitNumber1 = 0; hitNumber1 < numberOfHit ; hitNumber1++) {
 901       trackParamAtHit1 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber1));
 902       hitForRec1 = trackParamAtHit1->GetHitForRecPtr();
 903       hbc1 = hitForRec1->GetBendingCoor();
 904       pbc1 = trackParamAtHit1->GetBendingCoor();
 905       hnbc1 = hitForRec1->GetNonBendingCoor();
 906       pnbc1 = trackParamAtHit1->GetNonBendingCoor();
 907       Chi2 += ((hbc1 - pbc1) * (hbc1 - pbc1) / hitForRec1->GetBendingReso2()) +
 908               ((hnbc1 - pnbc1) * (hnbc1 - pnbc1) / hitForRec1->GetNonBendingReso2());
 909     }
 910   }
 911
 912   delete covBending;
 913   delete covNonBending;
 914   delete [] msa2;
 915 }
 916
 917   //__________________________________________________________________________
 918 Double_t MultipleScatteringAngle2(AliMUONTrackParam *param)
 919 {
 920   /// Returns square of multiple Coulomb scattering angle
 921   /// from TrackParamAtHit pointed to by "param"
 922   Double_t slopeBending, slopeNonBending, radiationLength, inverseBendingMomentum2, inverseTotalMomentum2;
 923   Double_t varMultipleScatteringAngle;
 924   slopeBending = param->GetBendingSlope();
 925   slopeNonBending = param->GetNonBendingSlope();
 926   // thickness in radiation length for the current track,
 927   // taking local angle into account
 928   radiationLength = AliMUONConstants::ChamberThicknessInX0() *
 929                     TMath::Sqrt(1.0 + slopeBending*slopeBending + slopeNonBending*slopeNonBending);
 930   inverseBendingMomentum2 =  param->GetInverseBendingMomentum() * param->GetInverseBendingMomentum();
 931   inverseTotalMomentum2 = inverseBendingMomentum2 * (1.0 + slopeBending * slopeBending) /
 932                           (1.0 + slopeBending *slopeBending + slopeNonBending * slopeNonBending);
 933   varMultipleScatteringAngle = 0.0136 * (1.0 + 0.038 * TMath::Log(radiationLength));
 934   // The velocity is assumed to be 1 !!!!
 935   varMultipleScatteringAngle = inverseTotalMomentum2 * radiationLength * varMultipleScatteringAngle * varMultipleScatteringAngle;
 936   return varMultipleScatteringAngle;
 937 }
 938
 939   //__________________________________________________________________________
 940 void AliMUONTrackReconstructor::ExtrapTracksToVertex()
 941 {
 942   /// propagate tracks to the vertex through the absorber (Branson)
 943   AliMUONTrack *track;
 944   AliMUONTrackParam trackParamVertex;
 945   AliMUONHitForRec *vertex;
 946   Double_t vX, vY, vZ;
 947
 948   for (Int_t iRecTrack = 0; iRecTrack <fNRecTracks; iRecTrack++) {
 949     track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iRecTrack);
 950     trackParamVertex = *((AliMUONTrackParam*)(track->GetTrackParamAtHit()->First()));
 951     vertex = track->GetVertex();
 952     if (vertex) { // if track as a vertex defined, use it
 953       vX = vertex->GetNonBendingCoor();
 954       vY = vertex->GetBendingCoor();
 955       vZ = vertex->GetZ();
 956     } else { // else vertex = (0.,0.,0.)
 957       vX = 0.;
 958       vY = 0.;
 959       vZ = 0.;
 960     }
 961     AliMUONTrackExtrap::ExtrapToVertex(&trackParamVertex, vX, vY, vZ);
 962     track->SetTrackParamAtVertex(&trackParamVertex);
 963
 964     if (AliLog::GetGlobalDebugLevel() > 0) {
 965       cout << "FollowTracks: track candidate(0..): " << iRecTrack
 966            << " after extrapolation to vertex" << endl;
 967       track->RecursiveDump();
 968     }
 969   }
 970
 971 }
 972
 973   //__________________________________________________________________________
 974 void AliMUONTrackReconstructor::FillMUONTrack()
 975 {
 976   /// Fill fHitForRecAtHit of AliMUONTrack's
 977   AliMUONTrack *track;
 978   AliMUONTrackParam *trackParamAtHit;
 979   track = (AliMUONTrack*) fRecTracksPtr->First();
 980   while (track) {
 981     trackParamAtHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First());
 982     while (trackParamAtHit) {
 983       track->AddHitForRecAtHit(trackParamAtHit->GetHitForRecPtr());
 984       trackParamAtHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->After(trackParamAtHit));
 985     }
 986     track = (AliMUONTrack*) fRecTracksPtr->After(track);
 987   }
 988   return;
 989 }
 990
 991   //__________________________________________________________________________
 992 void AliMUONTrackReconstructor::EventDump(void)
 993 {
 994   /// Dump reconstructed event (track parameters at vertex and at first hit),
 995   /// and the particle parameters
 996   AliMUONTrack *track;
 997   AliMUONTrackParam *trackParam, *trackParam1;
 998   Double_t bendingSlope, nonBendingSlope, pYZ;
 999   Double_t pX, pY, pZ, x, y, z, c;
1000   Int_t trackIndex, nTrackHits;
1001
1002   AliDebug(1,"****** enter EventDump ******");
1003   AliDebug(1, Form("Number of Reconstructed tracks : %d", fNRecTracks));
1004
1005   fRecTracksPtr->Compress(); // for simple loop without "Next" since no hole
1006   // Loop over reconstructed tracks
1007   for (trackIndex = 0; trackIndex < fNRecTracks; trackIndex++) {
1008     AliDebug(1, Form("track number: %d", trackIndex));
1009     // function for each track for modularity ????
1010     track = (AliMUONTrack*) ((*fRecTracksPtr)[trackIndex]);
1011     nTrackHits = track->GetNTrackHits();
1012     AliDebug(1, Form("Number of track hits: %d ", nTrackHits));
1013     // track parameters at Vertex
1014     trackParam = track->GetTrackParamAtVertex();
1015     x = trackParam->GetNonBendingCoor();
1016     y = trackParam->GetBendingCoor();
1017     z = trackParam->GetZ();
1018     bendingSlope = trackParam->GetBendingSlope();
1019     nonBendingSlope = trackParam->GetNonBendingSlope();
1020     pYZ = 1/TMath::Abs(trackParam->GetInverseBendingMomentum());
1021     pZ = pYZ/TMath::Sqrt(1+bendingSlope*bendingSlope);
1022     pX = pZ * nonBendingSlope;
1023     pY = pZ * bendingSlope;
1024     c = TMath::Sign(1.0, trackParam->GetInverseBendingMomentum());
1025     AliDebug(1, Form("Track parameters at Vertex z= %f: X= %f Y= %f pX= %f pY= %f pZ= %f c= %f\n",
1026                      z, x, y, pX, pY, pZ, c));
1027
1028     // track parameters at first hit
1029     trackParam1 = (AliMUONTrackParam*) track->GetTrackParamAtHit()->First();
1030     x = trackParam1->GetNonBendingCoor();
1031     y = trackParam1->GetBendingCoor();
1032     z = trackParam1->GetZ();
1033     bendingSlope = trackParam1->GetBendingSlope();
1034     nonBendingSlope = trackParam1->GetNonBendingSlope();
1035     pYZ = 1/TMath::Abs(trackParam1->GetInverseBendingMomentum());
1036     pZ = pYZ/TMath::Sqrt(1.0 + bendingSlope * bendingSlope);
1037     pX = pZ * nonBendingSlope;
1038     pY = pZ * bendingSlope;
1039     c = TMath::Sign(1.0, trackParam1->GetInverseBendingMomentum());
1040     AliDebug(1, Form("track parameters at z= %f: X= %f Y= %f pX= %f pY= %f pZ= %f c= %f\n",
1041                      z, x, y, pX, pY, pZ, c));
1042   }
1043   // informations about generated particles NO !!!!!!!!
1044
1045 //    for (Int_t iPart = 0; iPart < np; iPart++) {
1046 //      p = gAlice->Particle(iPart);
1047 //      printf(" particle %d: type= %d px= %f py= %f pz= %f pdg= %d\n",
1048 //         iPart, p->GetPdgCode(), p->Px(), p->Py(), p->Pz(), p->GetPdgCode());
1049 //    }
1050   return;
1051 }
1052
1053