MUON/AliMUONTrackReconstructor.cxx

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