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 ////////////////////////////////////
  19 //
  20 // MUON track reconstructor in ALICE (class renamed from AliMUONEventReconstructor)
  21 //
  22 // This class contains as data:
  23 // * the parameters for the track reconstruction
  24 // * a pointer to the array of hits to be reconstructed (the event)
  25 // * a pointer to the array of segments made with these hits inside each station
  26 // * a pointer to the array of reconstructed tracks
  27 //
  28 // It contains as methods, among others:
  29 // * MakeEventToBeReconstructed to build the array of hits to be reconstructed
  30 // * MakeSegments to build the segments
  31 // * MakeTracks to build the tracks
  32 //
  33 ////////////////////////////////////
  34
  35 #include <stdlib.h>
  36 #include <Riostream.h>
  37 #include <TDirectory.h>
  38 #include <TFile.h>
  39 #include <TMatrixD.h>
  40
  41 #include "AliMUONTrackReconstructor.h"
  42 #include "AliMUONData.h"
  43 #include "AliMUONConstants.h"
  44 #include "AliMUONHitForRec.h"
  45 #include "AliMUONTriggerTrack.h"
  46 #include "AliMUONTriggerCircuitNew.h"
  47 #include "AliMUONRawCluster.h"
  48 #include "AliMUONLocalTrigger.h"
  49 #include "AliMUONGlobalTrigger.h"
  50 #include "AliMUONSegment.h"
  51 #include "AliMUONTrack.h"
  52 #include "AliMagF.h"
  53 #include "AliMUONTrackK.h"
  54 #include "AliLog.h"
  55 #include "AliTracker.h"
  56 #include <TVirtualFitter.h>
  57
  58 //************* Defaults parameters for reconstruction
  59 const Double_t AliMUONTrackReconstructor::fgkDefaultMinBendingMomentum = 3.0;
  60 const Double_t AliMUONTrackReconstructor::fgkDefaultMaxBendingMomentum = 3000.0;
  61 const Double_t AliMUONTrackReconstructor::fgkDefaultMaxChi2 = 100.0;
  62 const Double_t AliMUONTrackReconstructor::fgkDefaultMaxSigma2Distance = 16.0;
  63 const Double_t AliMUONTrackReconstructor::fgkDefaultBendingResolution = 0.01;
  64 const Double_t AliMUONTrackReconstructor::fgkDefaultNonBendingResolution = 0.144;
  65 // Simple magnetic field:
  66 // Value taken from macro MUONtracking.C: 0.7 T, hence 7 kG
  67 // Length and Position from reco_muon.F, with opposite sign:
  68 // Length = ZMAGEND-ZCOIL
  69 // Position = (ZMAGEND+ZCOIL)/2
  70 // to be ajusted differently from real magnetic field ????
  71 const Double_t AliMUONTrackReconstructor::fgkDefaultSimpleBValue = 7.0;
  72 const Double_t AliMUONTrackReconstructor::fgkDefaultSimpleBLength = 428.0;
  73 const Double_t AliMUONTrackReconstructor::fgkDefaultSimpleBPosition = 1019.0;
  74 const Double_t AliMUONTrackReconstructor::fgkDefaultEfficiency = 0.95;
  75
  76 TVirtualFitter* AliMUONTrackReconstructor::fgFitter = NULL;
  77
  78 // Functions to be minimized with Minuit
  79 void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag);
  80 void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag);
  81
  82 void mnvertLocal(Double_t* a, Int_t l, Int_t m, Int_t n, Int_t& ifail);
  83
  84 Double_t MultipleScatteringAngle2(AliMUONTrackParam *param);
  85
  86 ClassImp(AliMUONTrackReconstructor) // Class implementation in ROOT context
  87
  88 //__________________________________________________________________________
  89 AliMUONTrackReconstructor::AliMUONTrackReconstructor(AliMUONData* data)
  90   : TObject(),
  91     fTrackMethod(1), //AZ - tracking method (1-default, 2-Kalman)
  92     fMinBendingMomentum(fgkDefaultMinBendingMomentum),
  93     fMaxBendingMomentum(fgkDefaultMaxBendingMomentum),
  94     fMaxChi2(fgkDefaultMaxChi2),
  95     fMaxSigma2Distance(fgkDefaultMaxSigma2Distance),
  96     fRMin(0x0),
  97     fRMax(0x0),
  98     fSegmentMaxDistBending(0x0),
  99     fSegmentMaxDistNonBending(0x0),
 100     fBendingResolution(fgkDefaultBendingResolution),
 101     fNonBendingResolution(fgkDefaultNonBendingResolution),
 102     fChamberThicknessInX0(AliMUONConstants::DefaultChamberThicknessInX0()),
 103     fSimpleBValue(fgkDefaultSimpleBValue),
 104     fSimpleBLength(fgkDefaultSimpleBLength),
 105     fSimpleBPosition(fgkDefaultSimpleBPosition),
 106     fEfficiency(fgkDefaultEfficiency),
 107     fHitsForRecPtr(0x0),
 108     fNHitsForRec(0),
 109     fNHitsForRecPerChamber(0x0),
 110     fIndexOfFirstHitForRecPerChamber(0x0),
 111     fSegmentsPtr(0x0),
 112     fNSegments(0x0),
 113     fRecTracksPtr(0x0),
 114     fNRecTracks(0),
 115     fMUONData(data),
 116     fMuons(0),
 117     fTriggerTrack(new AliMUONTriggerTrack()),
 118     fTriggerCircuit(0x0)
 119 {
 120
 121   // Memory allocation
 122   fRMin = new Double_t[AliMUONConstants::NTrackingCh()];
 123   fRMax = new Double_t[AliMUONConstants::NTrackingCh()];
 124   fSegmentMaxDistBending = new Double_t[AliMUONConstants::NTrackingSt()];
 125   fSegmentMaxDistNonBending = new Double_t[AliMUONConstants::NTrackingSt()];
 126   fNHitsForRecPerChamber = new Int_t[AliMUONConstants::NTrackingCh()];
 127   fIndexOfFirstHitForRecPerChamber = new Int_t[AliMUONConstants::NTrackingCh()];
 128
 129   // Constructor for class AliMUONTrackReconstructor
 130   SetReconstructionParametersToDefaults();
 131
 132   // Memory allocation for the TClonesArray of hits for reconstruction
 133   // Is 10000 the right size ????
 134   fHitsForRecPtr = new TClonesArray("AliMUONHitForRec", 10000);
 135
 136   // Memory allocation for the TClonesArray's of segments in stations
 137   // Is 2000 the right size ????
 138   fSegmentsPtr = new TClonesArray*[AliMUONConstants::NTrackingSt()];
 139   fNSegments = new Int_t[AliMUONConstants::NTrackingSt()];
 140   for (Int_t st = 0; st < AliMUONConstants::NTrackingSt(); st++) {
 141     fSegmentsPtr[st] = new TClonesArray("AliMUONSegment", 2000);
 142     fNSegments[st] = 0; // really needed or GetEntriesFast sufficient ????
 143   }
 144   // Memory allocation for the TClonesArray of reconstructed tracks
 145   // Is 10 the right size ????
 146   fRecTracksPtr = new TClonesArray("AliMUONTrack", 10);
 147
 148   const AliMagF* kField = AliTracker::GetFieldMap();
 149   if (!kField) AliFatal("No field available");
 150  // Sign of fSimpleBValue according to sign of Bx value at (50,50,-950).
 151   Float_t b[3], x[3];
 152   x[0] = 50.; x[1] = 50.; x[2] = -950.;
 153   kField->Field(x,b);
 154
 155   fSimpleBValue    = TMath::Sign(fSimpleBValue,(Double_t) b[0]);
 156   fSimpleBPosition = TMath::Sign(fSimpleBPosition,(Double_t) x[2]);
 157   // See how to get fSimple(BValue, BLength, BPosition)
 158   // automatically calculated from the actual magnetic field ????
 159
 160   return;
 161 }
 162
 163   //__________________________________________________________________________
 164 AliMUONTrackReconstructor::~AliMUONTrackReconstructor(void)
 165 {
 166   // Destructor for class AliMUONTrackReconstructor
 167   delete [] fRMin;
 168   delete [] fRMax;
 169   delete [] fSegmentMaxDistBending;
 170   delete [] fSegmentMaxDistNonBending;
 171   delete [] fNHitsForRecPerChamber;
 172   delete [] fIndexOfFirstHitForRecPerChamber;
 173   delete fTriggerTrack;
 174   delete fHitsForRecPtr;
 175   // Correct destruction of everything ????
 176   for (Int_t st = 0; st < AliMUONConstants::NTrackingSt(); st++) delete fSegmentsPtr[st];
 177   delete [] fSegmentsPtr;
 178   delete [] fNSegments;
 179   delete fRecTracksPtr;
 180 }
 181   //__________________________________________________________________________
 182 void AliMUONTrackReconstructor::SetReconstructionParametersToDefaults(void)
 183 {
 184   // Set reconstruction parameters to default values
 185   // Would be much more convenient with a structure (or class) ????
 186
 187   // ******** Parameters for making HitsForRec
 188   // minimum radius,
 189   // like in TRACKF_STAT:
 190   // 2 degrees for stations 1 and 2, or ch(0...) from 0 to 3;
 191   // 30 cm for stations 3 to 5, or ch(0...) from 4 to 9
 192   for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
 193     if (ch < 4) fRMin[ch] = TMath::Abs(AliMUONConstants::DefaultChamberZ(ch)) *
 194                   2.0 * TMath::Pi() / 180.0;
 195     else fRMin[ch] = 30.0;
 196     // maximum radius at 10 degrees and Z of chamber
 197     fRMax[ch] = TMath::Abs(AliMUONConstants::DefaultChamberZ(ch)) *
 198                   10.0 * TMath::Pi() / 180.0;
 199   }
 200
 201   // ******** Parameters for making segments
 202   // should be parametrized ????
 203   // according to interval between chambers in a station ????
 204   // Maximum distance in non bending plane
 205   // 5 * 0.22 just to remember the way it was made in TRACKF_STAT
 206   // SIGCUT*DYMAX(IZ)
 207   for (Int_t st = 0; st < AliMUONConstants::NTrackingCh()/2; st++)
 208     fSegmentMaxDistNonBending[st] = 5. * 0.22;
 209   // Maximum distance in bending plane:
 210   // values from TRACKF_STAT, corresponding to (J psi 20cm),
 211   // scaled to the real distance between chambers in a station
 212   fSegmentMaxDistBending[0] = TMath::Abs( 1.5 *
 213                                           (AliMUONConstants::DefaultChamberZ(1) - AliMUONConstants::DefaultChamberZ(0)) / 20.0);
 214   fSegmentMaxDistBending[1] = TMath::Abs( 1.5 *
 215                                           (AliMUONConstants::DefaultChamberZ(3) - AliMUONConstants::DefaultChamberZ(2)) / 20.0);
 216   fSegmentMaxDistBending[2] = TMath::Abs( 3.0 *
 217                                           (AliMUONConstants::DefaultChamberZ(5) - AliMUONConstants::DefaultChamberZ(4)) / 20.0);
 218   fSegmentMaxDistBending[3] = TMath::Abs( 6.0 *
 219                                           (AliMUONConstants::DefaultChamberZ(7) - AliMUONConstants::DefaultChamberZ(6)) / 20.0);
 220   fSegmentMaxDistBending[4] = TMath::Abs( 6.0 *
 221                                           (AliMUONConstants::DefaultChamberZ(9) - AliMUONConstants::DefaultChamberZ(8)) / 20.0);
 222
 223   return;
 224 }
 225
 226 //__________________________________________________________________________
 227 Double_t AliMUONTrackReconstructor::GetImpactParamFromBendingMomentum(Double_t BendingMomentum) const
 228 {
 229   // Returns impact parameter at vertex in bending plane (cm),
 230   // from the signed bending momentum "BendingMomentum" in bending plane (GeV/c),
 231   // using simple values for dipole magnetic field.
 232   // The sign of "BendingMomentum" is the sign of the charge.
 233   return (-0.0003 * fSimpleBValue * fSimpleBLength * fSimpleBPosition /
 234           BendingMomentum);
 235 }
 236
 237 //__________________________________________________________________________
 238 Double_t AliMUONTrackReconstructor::GetBendingMomentumFromImpactParam(Double_t ImpactParam) const
 239 {
 240   // Returns signed bending momentum in bending plane (GeV/c),
 241   // the sign being the sign of the charge for particles moving forward in Z,
 242   // from the impact parameter "ImpactParam" at vertex in bending plane (cm),
 243   // using simple values for dipole magnetic field.
 244   return (-0.0003 * fSimpleBValue * fSimpleBLength * fSimpleBPosition /
 245           ImpactParam);
 246 }
 247
 248 //__________________________________________________________________________
 249 void AliMUONTrackReconstructor::EventReconstruct(void)
 250 {
 251   // To reconstruct one event
 252   AliDebug(1,"Enter EventReconstruct");
 253   ResetTracks(); //AZ
 254   ResetSegments(); //AZ
 255   ResetHitsForRec(); //AZ
 256   MakeEventToBeReconstructed();
 257   MakeSegments();
 258   MakeTracks();
 259   if (fMUONData->IsTriggerTrackBranchesInTree())
 260     ValidateTracksWithTrigger();
 261
 262   // Add tracks to MUON data container
 263   for(Int_t i=0; i<fNRecTracks; i++) {
 264     AliMUONTrack * track = (AliMUONTrack*) fRecTracksPtr->At(i);
 265     fMUONData->AddRecTrack(*track);
 266   }
 267 }
 268
 269 //__________________________________________________________________________
 270 void AliMUONTrackReconstructor::EventReconstructTrigger(void)
 271 {
 272   // To reconstruct one event
 273   AliDebug(1,"Enter EventReconstructTrigger");
 274   MakeTriggerTracks();
 275 }
 276
 277   //__________________________________________________________________________
 278 void AliMUONTrackReconstructor::ResetHitsForRec(void)
 279 {
 280   // To reset the array and the number of HitsForRec,
 281   // and also the number of HitsForRec
 282   // and the index of the first HitForRec per chamber
 283   if (fHitsForRecPtr) fHitsForRecPtr->Delete();
 284   fNHitsForRec = 0;
 285   for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++)
 286     fNHitsForRecPerChamber[ch] = fIndexOfFirstHitForRecPerChamber[ch] = 0;
 287   return;
 288 }
 289
 290   //__________________________________________________________________________
 291 void AliMUONTrackReconstructor::ResetSegments(void)
 292 {
 293   // To reset the TClonesArray of segments and the number of Segments
 294   // for all stations
 295   for (Int_t st = 0; st < AliMUONConstants::NTrackingCh()/2; st++) {
 296     if (fSegmentsPtr[st]) fSegmentsPtr[st]->Clear();
 297     fNSegments[st] = 0;
 298   }
 299   return;
 300 }
 301
 302   //__________________________________________________________________________
 303 void AliMUONTrackReconstructor::ResetTracks(void)
 304 {
 305   // To reset the TClonesArray of reconstructed tracks
 306   if (fRecTracksPtr) fRecTracksPtr->Delete();
 307   // Delete in order that the Track destructors are called,
 308   // hence the space for the TClonesArray of pointers to TrackHit's is freed
 309   fNRecTracks = 0;
 310   return;
 311 }
 312
 313   //__________________________________________________________________________
 314 void AliMUONTrackReconstructor::MakeEventToBeReconstructed(void)
 315 {
 316   // To make the list of hits to be reconstructed,
 317   // either from the track ref. hits or from the raw clusters
 318   // according to the parameter set for the reconstructor
 319
 320   AliDebug(1,"Enter MakeEventToBeReconstructed");
 321   //AZ ResetHitsForRec();
 322
 323   // Reconstruction from raw clusters
 324   // AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
 325   // Security on MUON ????
 326   // TreeR assumed to be be "prepared" in calling function
 327   // by "MUON->GetTreeR(nev)" ????
 328   TTree *treeR = fMUONData->TreeR();
 329
 330   //AZ? fMUONData->SetTreeAddress("RC");
 331   AddHitsForRecFromRawClusters(treeR);
 332   // No sorting: it is done automatically in the previous function
 333
 334
 335   AliDebug(1,"End of MakeEventToBeReconstructed");
 336     if (AliLog::GetGlobalDebugLevel() > 0) {
 337       AliDebug(1, Form("NHitsForRec: %d",fNHitsForRec));
 338       for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
 339         AliDebug(1, Form("Chamber(0...): %d",ch));
 340         AliDebug(1, Form("NHitsForRec: %d", fNHitsForRecPerChamber[ch]));
 341         AliDebug(1, Form("Index(first HitForRec): %d", fIndexOfFirstHitForRecPerChamber[ch]));
 342         for (Int_t hit = fIndexOfFirstHitForRecPerChamber[ch];
 343              hit < fIndexOfFirstHitForRecPerChamber[ch] + fNHitsForRecPerChamber[ch];
 344              hit++) {
 345           AliDebug(1, Form("HitForRec index(0...): %d",hit));
 346           ((*fHitsForRecPtr)[hit])->Dump();
 347       }
 348     }
 349   }
 350   return;
 351 }
 352
 353   //__________________________________________________________________________
 354 void AliMUONTrackReconstructor::SortHitsForRecWithIncreasingChamber()
 355 {
 356   // Sort HitsForRec's in increasing order with respect to chamber number.
 357   // Uses the function "Compare".
 358   // Update the information for HitsForRec per chamber too.
 359   Int_t ch, nhits, prevch;
 360   fHitsForRecPtr->Sort();
 361   for (ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
 362     fNHitsForRecPerChamber[ch] = 0;
 363     fIndexOfFirstHitForRecPerChamber[ch] = 0;
 364   }
 365   prevch = 0; // previous chamber
 366   nhits = 0; // number of hits in current chamber
 367   // Loop over HitsForRec
 368   for (Int_t hit = 0; hit < fNHitsForRec; hit++) {
 369     // chamber number (0...)
 370     ch = ((AliMUONHitForRec*)  ((*fHitsForRecPtr)[hit]))->GetChamberNumber();
 371     // increment number of hits in current chamber
 372     (fNHitsForRecPerChamber[ch])++;
 373     // update index of first HitForRec in current chamber
 374     // if chamber number different from previous one
 375     if (ch != prevch) {
 376       fIndexOfFirstHitForRecPerChamber[ch] = hit;
 377       prevch = ch;
 378     }
 379   }
 380   return;
 381 }
 382
 383   //__________________________________________________________________________
 384 void AliMUONTrackReconstructor::AddHitsForRecFromRawClusters(TTree* TR)
 385 {
 386   // To add to the list of hits for reconstruction all the raw clusters
 387   // No condition added, like in Fortran TRACKF_STAT,
 388   // on the radius between RMin and RMax.
 389   AliMUONHitForRec *hitForRec;
 390   AliMUONRawCluster *clus;
 391   Int_t iclus, nclus, nTRentries;
 392   TClonesArray *rawclusters;
 393   AliDebug(1,"Enter AddHitsForRecFromRawClusters");
 394
 395   if (fTrackMethod != 3) { //AZ
 396     fMUONData->SetTreeAddress("RC"); //AZ
 397     nTRentries = Int_t(TR->GetEntries());
 398     if (nTRentries != 1) {
 399       AliError(Form("nTRentries = %d not equal to 1 ",nTRentries));
 400       exit(0);
 401     }
 402     fMUONData->GetRawClusters(); // only one entry
 403   }
 404
 405   // Loop over tracking chambers
 406   for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
 407     // number of HitsForRec to 0 for the chamber
 408     fNHitsForRecPerChamber[ch] = 0;
 409     // index of first HitForRec for the chamber
 410     if (ch == 0) fIndexOfFirstHitForRecPerChamber[ch] = 0;
 411     else fIndexOfFirstHitForRecPerChamber[ch] = fNHitsForRec;
 412     rawclusters =fMUONData->RawClusters(ch);
 413     nclus = (Int_t) (rawclusters->GetEntries());
 414     // Loop over (cathode correlated) raw clusters
 415     for (iclus = 0; iclus < nclus; iclus++) {
 416       clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(iclus);
 417       // new AliMUONHitForRec from raw cluster
 418       // and increment number of AliMUONHitForRec's (total and in chamber)
 419       hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(clus);
 420       fNHitsForRec++;
 421       (fNHitsForRecPerChamber[ch])++;
 422       // more information into HitForRec
 423       //  resolution: info should be already in raw cluster and taken from it ????
 424       //hitForRec->SetBendingReso2(fBendingResolution * fBendingResolution);
 425       //hitForRec->SetNonBendingReso2(fNonBendingResolution * fNonBendingResolution);
 426       hitForRec->SetBendingReso2(clus->GetErrY() * clus->GetErrY());
 427       hitForRec->SetNonBendingReso2(clus->GetErrX() * clus->GetErrX());
 428       //  original raw cluster
 429       hitForRec->SetChamberNumber(ch);
 430       hitForRec->SetHitNumber(iclus);
 431       // Z coordinate of the raw cluster (cm)
 432       hitForRec->SetZ(clus->GetZ(0));
 433
 434       StdoutToAliDebug(3,
 435                        cout << "Chamber " << ch <<
 436                        " raw cluster  " << iclus << " : " << endl;
 437                        clus->Print("full");
 438                        cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl;
 439                        hitForRec->Print("full");
 440                        );
 441     } // end of cluster loop
 442   } // end of chamber loop
 443   SortHitsForRecWithIncreasingChamber();
 444   return;
 445 }
 446
 447   //__________________________________________________________________________
 448 void AliMUONTrackReconstructor::MakeSegments(void)
 449 {
 450   // To make the list of segments in all stations,
 451   // from the list of hits to be reconstructed
 452   AliDebug(1,"Enter MakeSegments");
 453   //AZ ResetSegments();
 454   // Loop over stations
 455   Int_t nb = (fTrackMethod != 1) ? 3 : 0; //AZ
 456   for (Int_t st = nb; st < AliMUONConstants::NTrackingCh()/2; st++)
 457   {
 458     MakeSegmentsPerStation(st);
 459   }
 460
 461   StdoutToAliDebug(3,
 462     cout << "end of MakeSegments" << endl;
 463     for (Int_t st = 0; st < AliMUONConstants::NTrackingCh()/2; st++)
 464     {
 465       cout << "station " << st
 466             << "  has " << fNSegments[st] << " segments:"
 467             << endl;
 468       for (Int_t seg = 0; seg < fNSegments[st]; seg++)
 469       {
 470               ((*fSegmentsPtr[st])[seg])->Print();
 471       }
 472     }
 473                    );
 474 }
 475
 476   //__________________________________________________________________________
 477 void AliMUONTrackReconstructor::MakeSegmentsPerStation(Int_t Station)
 478 {
 479   // To make the list of segments in station number "Station" (0...)
 480   // from the list of hits to be reconstructed.
 481   // Updates "fNSegments"[Station].
 482   // Segments in stations 4 and 5 are sorted
 483   // according to increasing absolute value of "impact parameter"
 484   AliMUONHitForRec *hit1Ptr, *hit2Ptr;
 485   AliMUONSegment *segment;
 486   Bool_t last2st;
 487   Double_t bendingSlope, distBend, distNonBend, extBendCoor, extNonBendCoor,
 488       impactParam = 0., maxImpactParam = 0., minImpactParam = 0.; // =0 to avoid compilation warnings.
 489   AliDebug(1,Form("Enter MakeSegmentsPerStation (0...) %d",Station));
 490   // first and second chambers (0...) in the station
 491   Int_t ch1 = 2 * Station;
 492   Int_t ch2 = ch1 + 1;
 493   // variable true for stations downstream of the dipole:
 494   // Station(0..4) equal to 3 or 4
 495   if ((Station == 3) || (Station == 4)) {
 496     last2st = kTRUE;
 497     // maximum impact parameter (cm) according to fMinBendingMomentum...
 498     maxImpactParam =
 499       TMath::Abs(GetImpactParamFromBendingMomentum(fMinBendingMomentum));
 500     // minimum impact parameter (cm) according to fMaxBendingMomentum...
 501     minImpactParam =
 502       TMath::Abs(GetImpactParamFromBendingMomentum(fMaxBendingMomentum));
 503   }
 504   else last2st = kFALSE;
 505   // extrapolation factor from Z of first chamber to Z of second chamber
 506   // dZ to be changed to take into account fine structure of chambers ????
 507   Double_t extrapFact;
 508   // index for current segment
 509   Int_t segmentIndex = 0;
 510   // Loop over HitsForRec in the first chamber of the station
 511   for (Int_t hit1 = fIndexOfFirstHitForRecPerChamber[ch1];
 512        hit1 < fIndexOfFirstHitForRecPerChamber[ch1] + fNHitsForRecPerChamber[ch1];
 513        hit1++) {
 514     // pointer to the HitForRec
 515     hit1Ptr = (AliMUONHitForRec*) ((*fHitsForRecPtr)[hit1]);
 516     // extrapolation,
 517     // on the straight line joining the HitForRec to the vertex (0,0,0),
 518     // to the Z of the second chamber of the station
 519     // Loop over HitsForRec in the second chamber of the station
 520     for (Int_t hit2 = fIndexOfFirstHitForRecPerChamber[ch2];
 521          hit2 < fIndexOfFirstHitForRecPerChamber[ch2] + fNHitsForRecPerChamber[ch2];
 522          hit2++) {
 523       // pointer to the HitForRec
 524       hit2Ptr = (AliMUONHitForRec*) ((*fHitsForRecPtr)[hit2]);
 525       // absolute values of distances, in bending and non bending planes,
 526       // between the HitForRec in the second chamber
 527       // and the previous extrapolation
 528       extrapFact = hit2Ptr->GetZ()/ hit1Ptr->GetZ();
 529       extBendCoor = extrapFact * hit1Ptr->GetBendingCoor();
 530       extNonBendCoor = extrapFact * hit1Ptr->GetNonBendingCoor();
 531       distBend = TMath::Abs(hit2Ptr->GetBendingCoor() - extBendCoor);
 532       distNonBend = TMath::Abs(hit2Ptr->GetNonBendingCoor() - extNonBendCoor);
 533       if (last2st) {
 534         // bending slope
 535         if ( hit1Ptr->GetZ() - hit2Ptr->GetZ() != 0.0 ) {
 536           bendingSlope = (hit1Ptr->GetBendingCoor() - hit2Ptr->GetBendingCoor()) /
 537             (hit1Ptr->GetZ() - hit2Ptr->GetZ());
 538           // absolute value of impact parameter
 539           impactParam =
 540             TMath::Abs(hit1Ptr->GetBendingCoor() - hit1Ptr->GetZ() * bendingSlope);
 541          }
 542          else {
 543            AliWarning("hit1Ptr->GetZ() = hit2Ptr->GetZ(): impactParam set to maxImpactParam");
 544            impactParam = maxImpactParam;
 545          }
 546       }
 547       // check for distances not too large,
 548       // and impact parameter not too big if stations downstream of the dipole.
 549       // Conditions "distBend" and "impactParam" correlated for these stations ????
 550       if ((distBend < fSegmentMaxDistBending[Station]) &&
 551           (distNonBend < fSegmentMaxDistNonBending[Station]) &&
 552           (!last2st || (impactParam < maxImpactParam)) &&
 553           (!last2st || (impactParam > minImpactParam))) {
 554         // make new segment
 555         segment = new ((*fSegmentsPtr[Station])[segmentIndex])
 556           AliMUONSegment(hit1Ptr, hit2Ptr);
 557         // update "link" to this segment from the hit in the first chamber
 558         if (hit1Ptr->GetNSegments() == 0)
 559           hit1Ptr->SetIndexOfFirstSegment(segmentIndex);
 560         hit1Ptr->SetNSegments(hit1Ptr->GetNSegments() + 1);
 561         if (AliLog::GetGlobalDebugLevel() > 1) {
 562           cout << "segmentIndex(0...): " << segmentIndex
 563                << "  distBend: " << distBend
 564                << "  distNonBend: " << distNonBend
 565                << endl;
 566           segment->Dump();
 567           cout << "HitForRec in first chamber" << endl;
 568           hit1Ptr->Dump();
 569           cout << "HitForRec in second chamber" << endl;
 570           hit2Ptr->Dump();
 571         };
 572         // increment index for current segment
 573         segmentIndex++;
 574       }
 575     } //for (Int_t hit2
 576   } // for (Int_t hit1...
 577   fNSegments[Station] = segmentIndex;
 578   // Sorting according to "impact parameter" if station(1..5) 4 or 5,
 579   // i.e. Station(0..4) 3 or 4, using the function "Compare".
 580   // After this sorting, it is impossible to use
 581   // the "fNSegments" and "fIndexOfFirstSegment"
 582   // of the HitForRec in the first chamber to explore all segments formed with it.
 583   // Is this sorting really needed ????
 584   if ((Station == 3) || (Station == 4)) (fSegmentsPtr[Station])->Sort();
 585   AliDebug(1,Form("Station: %d  NSegments:  %d ", Station, fNSegments[Station]));
 586   return;
 587 }
 588
 589   //__________________________________________________________________________
 590 void AliMUONTrackReconstructor::MakeTracks(void)
 591 {
 592   // To make the tracks,
 593   // from the list of segments and points in all stations
 594   AliDebug(1,"Enter MakeTracks");
 595   // The order may be important for the following Reset's
 596   //AZ ResetTracks();
 597   if (fTrackMethod != 1) { //AZ - Kalman filter
 598     MakeTrackCandidatesK();
 599     if (fRecTracksPtr->GetEntriesFast() == 0) return;
 600     // Follow tracks in stations(1..) 3, 2 and 1
 601     FollowTracksK();
 602     // Remove double tracks
 603     RemoveDoubleTracksK();
 604     // Propagate tracks to the vertex thru absorber
 605     GoToVertex();
 606     // Fill AliMUONTrack data members
 607     FillMUONTrack();
 608   } else {
 609     // Look for candidates from at least 3 aligned points in stations(1..) 4 and 5
 610     MakeTrackCandidates();
 611     // Follow tracks in stations(1..) 3, 2 and 1
 612     FollowTracks();
 613     // Remove double tracks
 614     RemoveDoubleTracks();
 615     UpdateHitForRecAtHit();
 616   }
 617   return;
 618 }
 619
 620   //__________________________________________________________________________
 621 void AliMUONTrackReconstructor::ValidateTracksWithTrigger(void)
 622 {
 623   // Try to match track from tracking system with trigger track
 624   AliMUONTrack *track;
 625   AliMUONTrackParam trackParam;
 626   AliMUONTriggerTrack *triggerTrack;
 627
 628   fMUONData->SetTreeAddress("RL");
 629   fMUONData->GetRecTriggerTracks();
 630   TClonesArray *recTriggerTracks = fMUONData->RecTriggerTracks();
 631
 632   Bool_t MatchTrigger;
 633   Double_t distSigma[3]={1,1,0.02}; // sigma of distributions (trigger-track) X,Y,slopeY
 634   Double_t distTriggerTrack[3];
 635   Double_t Chi2MatchTrigger, xTrack, yTrack, ySlopeTrack, dTrigTrackMin2, dTrigTrack2 = 0.;
 636
 637   track = (AliMUONTrack*) fRecTracksPtr->First();
 638   while (track) {
 639     MatchTrigger = kFALSE;
 640     Chi2MatchTrigger = 0.;
 641
 642     trackParam = *((AliMUONTrackParam*) (track->GetTrackParamAtHit()->Last()));
 643     trackParam.ExtrapToZ(AliMUONConstants::DefaultChamberZ(10)); // extrap to 1st trigger chamber
 644
 645     xTrack = trackParam.GetNonBendingCoor();
 646     yTrack = trackParam.GetBendingCoor();
 647     ySlopeTrack = trackParam.GetBendingSlope();
 648     dTrigTrackMin2 = 999.;
 649
 650     triggerTrack = (AliMUONTriggerTrack*) recTriggerTracks->First();
 651     while(triggerTrack){
 652       distTriggerTrack[0] = (triggerTrack->GetX11()-xTrack)/distSigma[0];
 653       distTriggerTrack[1] = (triggerTrack->GetY11()-yTrack)/distSigma[1];
 654       distTriggerTrack[2] = (TMath::Tan(triggerTrack->GetThetay())-ySlopeTrack)/distSigma[2];
 655       dTrigTrack2 = 0.;
 656       for (Int_t iVar = 0; iVar < 3; iVar++) dTrigTrack2 += distTriggerTrack[iVar]*distTriggerTrack[iVar];
 657       if (dTrigTrack2 < dTrigTrackMin2 && dTrigTrack2 < GetMaxSigma2Distance()) {
 658         dTrigTrackMin2 = dTrigTrack2;
 659         MatchTrigger = kTRUE;
 660         Chi2MatchTrigger =  dTrigTrack2/3.; // Normalized Chi2, 3 variables (X,Y,slopeY)
 661       }
 662       triggerTrack = (AliMUONTriggerTrack*) recTriggerTracks->After(triggerTrack);
 663     }
 664
 665     track->SetMatchTrigger(MatchTrigger);
 666     track->SetChi2MatchTrigger(Chi2MatchTrigger);
 667
 668     track = (AliMUONTrack*) fRecTracksPtr->After(track);
 669   }
 670
 671 }
 672
 673   //__________________________________________________________________________
 674 Bool_t AliMUONTrackReconstructor::MakeTriggerTracks(void)
 675 {
 676     // To make the trigger tracks from Local Trigger
 677   AliDebug(1, "Enter MakeTriggerTracks");
 678
 679     Int_t nTRentries;
 680     Long_t gloTrigPat;
 681     TClonesArray *localTrigger;
 682     TClonesArray *globalTrigger;
 683     AliMUONLocalTrigger *locTrg;
 684     AliMUONGlobalTrigger *gloTrg;
 685
 686     TTree* treeR = fMUONData->TreeR();
 687
 688     nTRentries = Int_t(treeR->GetEntries());
 689
 690     treeR->GetEvent(0); // only one entry
 691
 692     if (!(fMUONData->IsTriggerBranchesInTree())) {
 693       AliWarning(Form("Trigger information is not avalaible, nTRentries = %d not equal to 1",nTRentries));
 694       return kFALSE;
 695     }
 696
 697     fMUONData->SetTreeAddress("TC");
 698     fMUONData->GetTrigger();
 699
 700     // global trigger for trigger pattern
 701     gloTrigPat = 0;
 702     globalTrigger = fMUONData->GlobalTrigger();
 703     gloTrg = (AliMUONGlobalTrigger*)globalTrigger->UncheckedAt(0);
 704
 705     if (gloTrg)
 706       gloTrigPat = gloTrg->GetGlobalPattern();
 707
 708
 709     // local trigger for tracking
 710     localTrigger = fMUONData->LocalTrigger();
 711     Int_t nlocals = (Int_t) (localTrigger->GetEntries());
 712
 713     Float_t z11 = AliMUONConstants::DefaultChamberZ(10);
 714     Float_t z21 = AliMUONConstants::DefaultChamberZ(12);
 715
 716     Float_t y11 = 0.;
 717     Int_t stripX21 = 0;
 718     Float_t y21 = 0.;
 719     Float_t x11 = 0.;
 720
 721     for (Int_t i=0; i<nlocals; i++) { // loop on Local Trigger
 722       locTrg = (AliMUONLocalTrigger*)localTrigger->UncheckedAt(i);
 723
 724       AliDebug(1, "AliMUONTrackReconstructor::MakeTriggerTrack using NEW trigger \n");
 725       AliMUONTriggerCircuitNew* circuit =
 726         (AliMUONTriggerCircuitNew*)fTriggerCircuit->At(locTrg->LoCircuit()-1); // -1 !!!
 727
 728       y11 = circuit->GetY11Pos(locTrg->LoStripX());
 729       stripX21 = locTrg->LoStripX()+locTrg->LoDev()+1;
 730       y21 = circuit->GetY21Pos(stripX21);
 731       x11 = circuit->GetX11Pos(locTrg->LoStripY());
 732
 733       AliDebug(1, Form(" MakeTriggerTrack %d %d %d %d %d %f %f %f \n",i,locTrg->LoCircuit(),
 734                        locTrg->LoStripX(),locTrg->LoStripX()+locTrg->LoDev()+1,locTrg->LoStripY(),y11, y21, x11));
 735
 736       Float_t thetax = TMath::ATan2( x11 , z11 );
 737       Float_t thetay = TMath::ATan2( (y21-y11) , (z21-z11) );
 738
 739       fTriggerTrack->SetX11(x11);
 740       fTriggerTrack->SetY11(y11);
 741       fTriggerTrack->SetThetax(thetax);
 742       fTriggerTrack->SetThetay(thetay);
 743       fTriggerTrack->SetGTPattern(gloTrigPat);
 744
 745       fMUONData->AddRecTriggerTrack(*fTriggerTrack);
 746     } // end of loop on Local Trigger
 747     return kTRUE;
 748 }
 749
 750   //__________________________________________________________________________
 751 void AliMUONTrackReconstructor::MakeTrackCandidates(void)
 752 {
 753   // To make track candidates
 754   // with at least 3 aligned points in stations(1..) 4 and 5
 755   // (two Segment's or one Segment and one HitForRec)
 756   Int_t begStation, iBegSegment, nbCan1Seg1Hit, nbCan2Seg;
 757   AliMUONSegment *begSegment;
 758   AliDebug(1,"Enter MakeTrackCandidates");
 759   // Loop over stations(1..) 5 and 4 for the beginning segment
 760   for (begStation = 4; begStation > 2; begStation--) {
 761     // Loop over segments in the beginning station
 762     for (iBegSegment = 0; iBegSegment < fNSegments[begStation]; iBegSegment++) {
 763       // pointer to segment
 764       begSegment = (AliMUONSegment*) ((*fSegmentsPtr[begStation])[iBegSegment]);
 765       AliDebug(2,Form("Look for TrackCandidate's with Segment %d  in Station(0..) %d", iBegSegment, begStation));
 766       // Look for track candidates with two segments,
 767       // "begSegment" and all compatible segments in other station.
 768       // Only for beginning station(1..) 5
 769       // because candidates with 2 segments have to looked for only once.
 770       if (begStation == 4)
 771         nbCan2Seg = MakeTrackCandidatesWithTwoSegments(begSegment);
 772       // Look for track candidates with one segment and one point,
 773       // "begSegment" and all compatible HitForRec's in other station.
 774       // Only if "begSegment" does not belong already to a track candidate.
 775       // Is that a too strong condition ????
 776       if (!(begSegment->GetInTrack()))
 777         nbCan1Seg1Hit = MakeTrackCandidatesWithOneSegmentAndOnePoint(begSegment);
 778     } // for (iBegSegment = 0;...
 779   } // for (begStation = 4;...
 780   return;
 781 }
 782
 783   //__________________________________________________________________________
 784 Int_t AliMUONTrackReconstructor::MakeTrackCandidatesWithTwoSegments(AliMUONSegment *BegSegment)
 785 {
 786   // To make track candidates with two segments in stations(1..) 4 and 5,
 787   // the first segment being pointed to by "BegSegment".
 788   // Returns the number of such track candidates.
 789   Int_t endStation, iEndSegment, nbCan2Seg;
 790   AliMUONSegment *endSegment;
 791   AliMUONSegment *extrapSegment = NULL;
 792   AliMUONTrack *recTrack;
 793   Double_t mcsFactor;
 794   AliDebug(1,"Enter MakeTrackCandidatesWithTwoSegments");
 795   // Station for the end segment
 796   endStation = 7 - (BegSegment->GetHitForRec1())->GetChamberNumber() / 2;
 797   // multiple scattering factor corresponding to one chamber
 798   mcsFactor = 0.0136 /
 799     GetBendingMomentumFromImpactParam(BegSegment->GetBendingImpact());
 800   mcsFactor     = fChamberThicknessInX0 * mcsFactor * mcsFactor;
 801   // linear extrapolation to end station
 802   // number of candidates with 2 segments to 0
 803   nbCan2Seg = 0;
 804   // Loop over segments in the end station
 805   for (iEndSegment = 0; iEndSegment < fNSegments[endStation]; iEndSegment++) {
 806     // pointer to segment
 807     endSegment = (AliMUONSegment*) ((*fSegmentsPtr[endStation])[iEndSegment]);
 808     // test compatibility between current segment and "extrapSegment"
 809     // 4 because 4 quantities in chi2
 810     extrapSegment =
 811       BegSegment->CreateSegmentFromLinearExtrapToStation(endSegment->GetZ(), mcsFactor);
 812     if ((endSegment->
 813          NormalizedChi2WithSegment(extrapSegment,
 814                                    fMaxSigma2Distance)) <= 4.0) {
 815       // both segments compatible:
 816       // make track candidate from "begSegment" and "endSegment"
 817       AliDebug(2,Form("TrackCandidate with Segment %d in Station(0..) %d", iEndSegment, endStation));
 818       // flag for both segments in one track:
 819       // to be done in track constructor ????
 820       BegSegment->SetInTrack(kTRUE);
 821       endSegment->SetInTrack(kTRUE);
 822       recTrack = new ((*fRecTracksPtr)[fNRecTracks]) AliMUONTrack(BegSegment, endSegment);
 823       // Set track parameters at vertex from last stations 4 & 5
 824       CalcTrackParamAtVertex(recTrack);
 825       fNRecTracks++;
 826       if (AliLog::GetGlobalDebugLevel() > 1) recTrack->RecursiveDump();
 827       // increment number of track candidates with 2 segments
 828       nbCan2Seg++;
 829     }
 830     delete extrapSegment; // should not delete HitForRec's it points to !!!!
 831   } // for (iEndSegment = 0;...
 832   return nbCan2Seg;
 833 }
 834
 835   //__________________________________________________________________________
 836 Int_t AliMUONTrackReconstructor::MakeTrackCandidatesWithOneSegmentAndOnePoint(AliMUONSegment *BegSegment)
 837 {
 838   // To make track candidates with one segment and one point
 839   // in stations(1..) 4 and 5,
 840   // the segment being pointed to by "BegSegment".
 841   Int_t ch, ch1, ch2, endStation, iHit, iHitMax, iHitMin, nbCan1Seg1Hit;
 842   AliMUONHitForRec *extrapHitForRec= NULL;
 843   AliMUONHitForRec *hit;
 844   AliMUONTrack *recTrack;
 845   Double_t mcsFactor;
 846   AliDebug(1,"Enter MakeTrackCandidatesWithOneSegmentAndOnePoint");
 847   // station for the end point
 848   endStation = 7 - (BegSegment->GetHitForRec1())->GetChamberNumber() / 2;
 849   // multiple scattering factor corresponding to one chamber
 850   mcsFactor = 0.0136 /
 851     GetBendingMomentumFromImpactParam(BegSegment->GetBendingImpact());
 852   mcsFactor     = fChamberThicknessInX0 * mcsFactor * mcsFactor;
 853   // first and second chambers(0..) in the end station
 854   ch1 = 2 * endStation;
 855   ch2 = ch1 + 1;
 856   // number of candidates to 0
 857   nbCan1Seg1Hit = 0;
 858   // Loop over chambers of the end station
 859   for (ch = ch2; ch >= ch1; ch--) {
 860     // limits for the hit index in the loop
 861     iHitMin = fIndexOfFirstHitForRecPerChamber[ch];
 862     iHitMax = iHitMin + fNHitsForRecPerChamber[ch];
 863     // Loop over HitForRec's in the chamber
 864     for (iHit = iHitMin; iHit < iHitMax; iHit++) {
 865       // pointer to HitForRec
 866       hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[iHit]);
 867       // test compatibility between current HitForRec and "extrapHitForRec"
 868       // 2 because 2 quantities in chi2
 869       // linear extrapolation to chamber
 870       extrapHitForRec =
 871         BegSegment->CreateHitForRecFromLinearExtrapToChamber( hit->GetZ(), mcsFactor);
 872       if ((hit->
 873            NormalizedChi2WithHitForRec(extrapHitForRec,
 874                                        fMaxSigma2Distance)) <= 2.0) {
 875         // both HitForRec's compatible:
 876         // make track candidate from begSegment and current HitForRec
 877         AliDebug(2, Form("TrackCandidate with HitForRec  %d in Chamber(0..) %d", iHit, ch));
 878         // flag for beginning segments in one track:
 879         // to be done in track constructor ????
 880         BegSegment->SetInTrack(kTRUE);
 881         recTrack = new ((*fRecTracksPtr)[fNRecTracks]) AliMUONTrack(BegSegment, hit);
 882         // Set track parameters at vertex from last stations 4 & 5
 883         CalcTrackParamAtVertex(recTrack);
 884         // the right place to eliminate "double counting" ???? how ????
 885         fNRecTracks++;
 886         if (AliLog::GetGlobalDebugLevel() > 1) recTrack->RecursiveDump();
 887         // increment number of track candidates
 888         nbCan1Seg1Hit++;
 889       }
 890       delete extrapHitForRec;
 891     } // for (iHit = iHitMin;...
 892   } // for (ch = ch2;...
 893   return nbCan1Seg1Hit;
 894 }
 895
 896   //__________________________________________________________________________
 897 void AliMUONTrackReconstructor::CalcTrackParamAtVertex(AliMUONTrack *Track)
 898 {
 899   // Set track parameters at vertex.
 900   // TrackHit's are assumed to be only in stations(1..) 4 and 5,
 901   // and sorted according to increasing Z..
 902   // Parameters are calculated from information in HitForRec's
 903   // of first and last TrackHit's.
 904   AliMUONTrackParam *trackParamVertex = Track->GetTrackParamAtVertex(); // pointer to track parameters at vertex
 905   // Pointer to HitForRec attached to first TrackParamAtHit
 906   AliMUONHitForRec *firstHit = ((AliMUONTrackParam*) (Track->GetTrackParamAtHit()->First()))->GetHitForRecPtr();
 907   // Pointer to HitForRec attached to last TrackParamAtHit
 908   AliMUONHitForRec *lastHit = ((AliMUONTrackParam*) (Track->GetTrackParamAtHit()->Last()))->GetHitForRecPtr();
 909   // Z difference between first and last hits
 910   Double_t deltaZ = firstHit->GetZ() - lastHit->GetZ();
 911   // bending slope in stations(1..) 4 and 5
 912   Double_t bendingSlope = (firstHit->GetBendingCoor() - lastHit->GetBendingCoor()) / deltaZ;
 913   trackParamVertex->SetBendingSlope(bendingSlope);
 914   // impact parameter
 915   Double_t impactParam = firstHit->GetBendingCoor() - bendingSlope * firstHit->GetZ();
 916   // signed bending momentum
 917   trackParamVertex->SetInverseBendingMomentum(1.0 / GetBendingMomentumFromImpactParam(impactParam));
 918   // bending slope at vertex
 919   trackParamVertex->SetBendingSlope(bendingSlope + impactParam / GetSimpleBPosition());
 920   // non bending slope
 921   trackParamVertex->SetNonBendingSlope((firstHit->GetNonBendingCoor() - lastHit->GetNonBendingCoor()) / deltaZ);
 922   // vertex coordinates at (0,0,0)
 923   trackParamVertex->SetZ(0.0);
 924   trackParamVertex->SetBendingCoor(0.0);
 925   trackParamVertex->SetNonBendingCoor(0.0);
 926 }
 927
 928   //__________________________________________________________________________
 929 void AliMUONTrackReconstructor::FollowTracks(void)
 930 {
 931   // Follow tracks in stations(1..) 3, 2 and 1
 932   // too long: should be made more modular !!!!
 933   AliMUONHitForRec *bestHit, *extrapHit, *hit;
 934   AliMUONSegment *bestSegment, *extrapSegment, *segment;
 935   AliMUONTrack *track, *nextTrack;
 936   AliMUONTrackParam *trackParam1, trackParam[2], trackParamVertex;
 937   // -1 to avoid compilation warnings
 938   Int_t ch = -1, chInStation, chBestHit = -1, iHit, iSegment, station, trackIndex;
 939   Double_t bestChi2, chi2, dZ1, dZ2, dZ3, maxSigma2Distance, mcsFactor;
 940   Double_t bendingMomentum, chi2Norm = 0.;
 941
 942
 943   // local maxSigma2Distance, for easy increase in testing
 944   maxSigma2Distance = fMaxSigma2Distance;
 945   AliDebug(2,"Enter FollowTracks");
 946   // Loop over track candidates
 947   track = (AliMUONTrack*) fRecTracksPtr->First();
 948   trackIndex = -1;
 949   while (track) {
 950     trackIndex++;
 951     nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track); // prepare next track
 952     AliDebug(2,Form("FollowTracks: track candidate(0..): %d", trackIndex));
 953     // Fit track candidate from parameters at vertex
 954     // -> with 3 parameters (X_vertex and Y_vertex are fixed)
 955     // without multiple Coulomb scattering
 956     Fit(track,0,0);
 957     if (AliLog::GetGlobalDebugLevel()> 2) {
 958       cout << "FollowTracks: track candidate(0..): " << trackIndex
 959            << " after fit in stations(0..) 3 and 4" << endl;
 960       track->RecursiveDump();
 961     }
 962     // Loop over stations(1..) 3, 2 and 1
 963     // something SPECIAL for stations 2 and 1 for majority 3 coincidence ????
 964     // otherwise: majority coincidence 2 !!!!
 965     for (station = 2; station >= 0; station--) {
 966       // Track parameters at first track hit (smallest Z)
 967       trackParam1 = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First());
 968       // extrapolation to station
 969       trackParam1->ExtrapToStation(station, trackParam);
 970       extrapSegment = new AliMUONSegment(); //  empty segment
 971       // multiple scattering factor corresponding to one chamber
 972       // and momentum in bending plane (not total)
 973       mcsFactor = 0.0136 * trackParam1->GetInverseBendingMomentum();
 974       mcsFactor = fChamberThicknessInX0 * mcsFactor * mcsFactor;
 975       // Z difference from previous station
 976       dZ1 = AliMUONConstants::DefaultChamberZ(2 * station) -
 977             AliMUONConstants::DefaultChamberZ(2 * station + 2);
 978       // Z difference between the two previous stations
 979       dZ2 = AliMUONConstants::DefaultChamberZ(2 * station + 2) -
 980             AliMUONConstants::DefaultChamberZ(2 * station + 4);
 981       // Z difference between the two chambers in the previous station
 982       dZ3 = AliMUONConstants::DefaultChamberZ(2 * station) -
 983             AliMUONConstants::DefaultChamberZ(2 * station + 1);
 984       extrapSegment->SetBendingCoorReso2(fBendingResolution * fBendingResolution);
 985       extrapSegment->SetNonBendingCoorReso2(fNonBendingResolution * fNonBendingResolution);
 986       extrapSegment->UpdateFromStationTrackParam(trackParam, mcsFactor, dZ1, dZ2, dZ3, station,
 987                                                  trackParam1->GetInverseBendingMomentum());
 988       bestChi2 = 5.0;
 989       bestSegment = NULL;
 990       if (AliLog::GetGlobalDebugLevel() > 2) {
 991         cout << "FollowTracks: track candidate(0..): " << trackIndex
 992              << " Look for segment in station(0..): " << station << endl;
 993       }
 994
 995       // Loop over segments in station
 996       for (iSegment = 0; iSegment < fNSegments[station]; iSegment++) {
 997         // Look for best compatible Segment in station
 998         // should consider all possibilities ????
 999         // multiple scattering ????
1000         // separation in 2 functions: Segment and HitForRec ????
1001         segment = (AliMUONSegment*) ((*fSegmentsPtr[station])[iSegment]);
1002         // correction of corrected segment (fBendingCoor and fNonBendingCoor)
1003         // according to real Z value of "segment" and slopes of "extrapSegment"
1004         trackParam[0].ExtrapToZ(segment->GetZ());
1005         trackParam[1].ExtrapToZ(segment->GetZ()); // now same as trackParam[0] !?!?!?!?!?!
1006         extrapSegment->SetBendingCoor((&(trackParam[0]))->GetBendingCoor());
1007         extrapSegment->SetNonBendingCoor((&(trackParam[0]))->GetNonBendingCoor());
1008         extrapSegment->SetBendingSlope((&(trackParam[0]))->GetBendingSlope());
1009         extrapSegment->SetNonBendingSlope((&(trackParam[0]))->GetNonBendingSlope());
1010         chi2 = segment->NormalizedChi2WithSegment(extrapSegment, maxSigma2Distance);
1011         if (chi2 < bestChi2) {
1012           // update best Chi2 and Segment if better found
1013           bestSegment = segment;
1014           bestChi2 = chi2;
1015         }
1016       }
1017       if (bestSegment) {
1018         // best segment found: add it to track candidate
1019         trackParam[0].ExtrapToZ(bestSegment->GetZ());
1020         track->AddTrackParamAtHit(&(trackParam[0]),bestSegment->GetHitForRec1());
1021         trackParam[1].ExtrapToZ(bestSegment->GetZ()); // now same as trackParam[0] !?!?!?!?!?!
1022         track->AddTrackParamAtHit(&(trackParam[1]),bestSegment->GetHitForRec2());
1023         AliDebug(3, Form("FollowTracks: track candidate(0..): %d  Added segment in station(0..): %d", trackIndex, station));
1024         if (AliLog::GetGlobalDebugLevel()>2) track->RecursiveDump();
1025       } else {
1026         // No best segment found:
1027         // Look for best compatible HitForRec in station:
1028         // should consider all possibilities ????
1029         // multiple scattering ???? do about like for extrapSegment !!!!
1030         extrapHit = new AliMUONHitForRec(); //  empty hit
1031         bestChi2 = 3.0;
1032         bestHit = NULL;
1033         AliDebug(3, Form("FollowTracks: track candidate(0..): %d Segment not found, look for hit in station(0..): %d ",
1034                          trackIndex, station));
1035
1036         // Loop over chambers of the station
1037         for (chInStation = 0; chInStation < 2; chInStation++) {
1038           ch = 2 * station + chInStation;
1039           for (iHit = fIndexOfFirstHitForRecPerChamber[ch]; iHit < fIndexOfFirstHitForRecPerChamber[ch]+fNHitsForRecPerChamber[ch]; iHit++) {
1040             hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[iHit]);
1041             // coordinates of extrapolated hit
1042             trackParam[chInStation].ExtrapToZ(hit->GetZ());
1043             extrapHit->SetBendingCoor((&(trackParam[chInStation]))->GetBendingCoor());
1044             extrapHit->SetNonBendingCoor((&(trackParam[chInStation]))->GetNonBendingCoor());
1045             // resolutions from "extrapSegment"
1046             extrapHit->SetBendingReso2(extrapSegment->GetBendingCoorReso2());
1047             extrapHit->SetNonBendingReso2(extrapSegment->GetNonBendingCoorReso2());
1048             // Loop over hits in the chamber
1049             // condition for hit not already in segment ????
1050             chi2 = hit->NormalizedChi2WithHitForRec(extrapHit, maxSigma2Distance);
1051             if (chi2 < bestChi2) {
1052               // update best Chi2 and HitForRec if better found
1053               bestHit = hit;
1054               bestChi2 = chi2;
1055               chBestHit = chInStation;
1056             }
1057           }
1058         }
1059         if (bestHit) {
1060           // best hit found: add it to track candidate
1061           trackParam[chBestHit].ExtrapToZ(bestHit->GetZ());
1062           track->AddTrackParamAtHit(&(trackParam[chBestHit]),bestHit);
1063           if (AliLog::GetGlobalDebugLevel() > 2) {
1064             cout << "FollowTracks: track candidate(0..): " << trackIndex
1065                  << " Added hit in station(0..): " << station << endl;
1066             track->RecursiveDump();
1067           }
1068         } else {
1069           // Remove current track candidate
1070           // and corresponding TrackHit's, ...
1071           fRecTracksPtr->Remove(track);
1072           fNRecTracks--;
1073           delete extrapSegment;
1074           delete extrapHit;
1075           break; // stop the search for this candidate:
1076           // exit from the loop over station
1077         }
1078         delete extrapHit;
1079       }
1080       delete extrapSegment;
1081       // Sort TrackParamAtHit according to increasing Z
1082       track->GetTrackParamAtHit()->Sort();
1083       // Update track parameters at first track hit (smallest Z)
1084       trackParam1 = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First());
1085       bendingMomentum = 0.;
1086       if (TMath::Abs(trackParam1->GetInverseBendingMomentum()) > 0.)
1087         bendingMomentum = TMath::Abs(1/(trackParam1->GetInverseBendingMomentum()));
1088       // Track removed if bendingMomentum not in window [min, max]
1089       if ((bendingMomentum < fMinBendingMomentum) || (bendingMomentum > fMaxBendingMomentum)) {
1090         fRecTracksPtr->Remove(track);
1091         fNRecTracks--;
1092         break; // stop the search for this candidate:
1093         // exit from the loop over station
1094       }
1095       // Track fit from parameters at first hit
1096       // -> with 5 parameters (momentum and position)
1097       // with multiple Coulomb scattering if all stations
1098       if (station == 0) Fit(track,1,1);
1099       // without multiple Coulomb scattering if not all stations
1100       else Fit(track,1,0);
1101       Double_t numberOfDegFree = (2.0 * track->GetNTrackHits() - 5);
1102       if (numberOfDegFree > 0) {
1103         chi2Norm =  track->GetFitFMin() / numberOfDegFree;
1104       } else {
1105         chi2Norm = 1.e10;
1106       }
1107       // Track removed if normalized chi2 too high
1108       if (chi2Norm > fMaxChi2) {
1109         fRecTracksPtr->Remove(track);
1110         fNRecTracks--;
1111         break; // stop the search for this candidate:
1112         // exit from the loop over station
1113       }
1114       if (AliLog::GetGlobalDebugLevel() > 2) {
1115         cout << "FollowTracks: track candidate(0..): " << trackIndex
1116              << " after fit from station(0..): " << station << " to 4" << endl;
1117         track->RecursiveDump();
1118       }
1119       // Track extrapolation to the vertex through the absorber (Branson)
1120       // after going through the first station
1121       if (station == 0) {
1122         trackParamVertex = *((AliMUONTrackParam*) (track->GetTrackParamAtHit()->First()));
1123         (&trackParamVertex)->ExtrapToVertex(0.,0.,0.);
1124         track->SetTrackParamAtVertex(&trackParamVertex);
1125         if (AliLog::GetGlobalDebugLevel() > 0) {
1126           cout << "FollowTracks: track candidate(0..): " << trackIndex
1127                << " after extrapolation to vertex" << endl;
1128           track->RecursiveDump();
1129         }
1130       }
1131     } // for (station = 2;...
1132     // go really to next track
1133     track = nextTrack;
1134   } // while (track)
1135   // Compression of track array (necessary after Remove)
1136   fRecTracksPtr->Compress();
1137   return;
1138 }
1139
1140   //__________________________________________________________________________
1141 void AliMUONTrackReconstructor::Fit(AliMUONTrack *Track, Int_t FitStart, Int_t FitMCS)
1142 {
1143   // Fit the track "Track",
1144   // with or without multiple Coulomb scattering according to "FitMCS",
1145   // starting, according to "FitStart",
1146   // with track parameters at vertex or at the first TrackHit.
1147
1148   if ((FitStart != 0) && (FitStart != 1)) {
1149     cout << "ERROR in AliMUONTrackReconstructor::Fit(...)" << endl;
1150     cout << "FitStart = " << FitStart << " is neither 0 nor 1" << endl;
1151     exit(0);
1152   }
1153   if ((FitMCS != 0) && (FitMCS != 1)) {
1154     cout << "ERROR in AliMUONTrackReconstructor::Fit(...)" << endl;
1155     cout << "FitMCS = " << FitMCS << " is neither 0 nor 1" << endl;
1156     exit(0);
1157   }
1158
1159   Double_t arg[1], benC, errorParam, invBenP, lower, nonBenC, upper, x, y;
1160   char parName[50];
1161   AliMUONTrackParam *trackParam;
1162   // Check if Minuit is initialized...
1163   fgFitter = TVirtualFitter::Fitter(Track,5);
1164   fgFitter->Clear(); // necessary ???? probably yes
1165   // how to reset the printout number at every fit ????
1166   // is there any risk to leave it like that ????
1167   // how to go faster ???? choice of Minuit parameters like EDM ????
1168   // choice of function to be minimized according to fFitMCS
1169   if (FitMCS == 0) fgFitter->SetFCN(TrackChi2);
1170   else fgFitter->SetFCN(TrackChi2MCS);
1171   // Switch off printout
1172   arg[0] = -1;
1173   fgFitter->ExecuteCommand("SET PRINT", arg, 1); // More printing !!!!
1174   // No warnings
1175   fgFitter->ExecuteCommand("SET NOW", arg, 0);
1176   // Parameters according to "fFitStart"
1177   // (should be a function to be used at every place where needed ????)
1178   if (FitStart == 0) trackParam = Track->GetTrackParamAtVertex();
1179   else trackParam = (AliMUONTrackParam*) (Track->GetTrackParamAtHit()->First());
1180   // set first 3 Minuit parameters
1181   // could be tried with no limits for the search (min=max=0) ????
1182   fgFitter->SetParameter(0, "InvBenP",
1183                          trackParam->GetInverseBendingMomentum(),
1184                          0.003, -0.4, 0.4);
1185   fgFitter->SetParameter(1, "BenS",
1186                          trackParam->GetBendingSlope(),
1187                          0.001, -0.5, 0.5);
1188   fgFitter->SetParameter(2, "NonBenS",
1189                          trackParam->GetNonBendingSlope(),
1190                          0.001, -0.5, 0.5);
1191   if (FitStart == 1) {
1192     // set last 2 Minuit parameters when we start from first track hit
1193     // mandatory limits in Bending to avoid NaN values of parameters
1194     fgFitter->SetParameter(3, "X",
1195                            trackParam->GetNonBendingCoor(),
1196                            0.03, -500.0, 500.0);
1197     // mandatory limits in non Bending to avoid NaN values of parameters
1198     fgFitter->SetParameter(4, "Y",
1199                            trackParam->GetBendingCoor(),
1200                            0.10, -500.0, 500.0);
1201   }
1202   // search without gradient calculation in the function
1203   fgFitter->ExecuteCommand("SET NOGRADIENT", arg, 0);
1204   // minimization
1205   fgFitter->ExecuteCommand("MINIMIZE", arg, 0);
1206   // exit from Minuit
1207   //  fgFitter->ExecuteCommand("EXIT", arg, 0); // necessary ????
1208   // get results into "invBenP", "benC", "nonBenC" ("x", "y")
1209   fgFitter->GetParameter(0, parName, invBenP, errorParam, lower, upper);
1210   trackParam->SetInverseBendingMomentum(invBenP);
1211   fgFitter->GetParameter(1, parName, benC, errorParam, lower, upper);
1212   trackParam->SetBendingSlope(benC);
1213   fgFitter->GetParameter(2, parName, nonBenC, errorParam, lower, upper);
1214   trackParam->SetNonBendingSlope(nonBenC);
1215   if (FitStart == 1) {
1216     fgFitter->GetParameter(3, parName, x, errorParam, lower, upper);
1217     trackParam->SetNonBendingCoor(x);
1218     fgFitter->GetParameter(4, parName, y, errorParam, lower, upper);
1219     trackParam->SetBendingCoor(y);
1220   }
1221   // global result of the fit
1222   Double_t fedm, errdef, FitFMin;
1223   Int_t npari, nparx;
1224   fgFitter->GetStats(FitFMin, fedm, errdef, npari, nparx);
1225   Track->SetFitFMin(FitFMin);
1226 }
1227
1228   //__________________________________________________________________________
1229 void TrackChi2(Int_t &NParam, Double_t * /*Gradient*/, Double_t &Chi2, Double_t *Param, Int_t /*Flag*/)
1230 {
1231   // Return the "Chi2" to be minimized with Minuit for track fitting,
1232   // with "NParam" parameters
1233   // and their current values in array pointed to by "Param".
1234   // Assumes that the track hits are sorted according to increasing Z.
1235   // Track parameters at each TrackHit are updated accordingly.
1236   // Multiple Coulomb scattering is not taken into account
1237   AliMUONTrack *trackBeingFitted;
1238   AliMUONTrackParam param1;
1239   AliMUONTrackParam* TrackParamAtHit;
1240   AliMUONHitForRec* HitForRec;
1241   Chi2 = 0.0; // initialize Chi2
1242   // copy of track parameters to be fitted
1243   trackBeingFitted = (AliMUONTrack*) AliMUONTrackReconstructor::Fitter()->GetObjectFit();
1244   // 3 parameters means fit track candidate from parameters at vertex (X_vertex and Y_vertex are fixed)
1245   if (NParam == 3) param1 = *(trackBeingFitted->GetTrackParamAtVertex());
1246   else param1 = *((AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->First()));
1247   // Minuit parameters to be fitted into this copy
1248   param1.SetInverseBendingMomentum(Param[0]);
1249   param1.SetBendingSlope(Param[1]);
1250   param1.SetNonBendingSlope(Param[2]);
1251   if (NParam == 5) {
1252     param1.SetNonBendingCoor(Param[3]);
1253     param1.SetBendingCoor(Param[4]);
1254   }
1255   // Follow track through all planes of track hits
1256   TrackParamAtHit = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->First());
1257   while (TrackParamAtHit) {
1258     HitForRec = TrackParamAtHit->GetHitForRecPtr();
1259     // extrapolation to the plane of the HitForRec attached to the current TrackParamAtHit
1260     param1.ExtrapToZ(HitForRec->GetZ());
1261     // update track parameters of the current hit
1262     TrackParamAtHit->SetTrackParam(param1);
1263     // Increment Chi2
1264     // done hit per hit, with hit resolution,
1265     // and not with point and angle like in "reco_muon.F" !!!!
1266     // Needs to add multiple scattering contribution ????
1267     Double_t dX = HitForRec->GetNonBendingCoor() - param1.GetNonBendingCoor();
1268     Double_t dY = HitForRec->GetBendingCoor() - param1.GetBendingCoor();
1269     Chi2 = Chi2 + dX * dX / HitForRec->GetNonBendingReso2() + dY * dY / HitForRec->GetBendingReso2();
1270     TrackParamAtHit = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->After(TrackParamAtHit));
1271   }
1272 }
1273
1274   //__________________________________________________________________________
1275 void TrackChi2MCS(Int_t &NParam, Double_t * /*Gradient*/, Double_t &Chi2, Double_t *Param, Int_t /*Flag*/)
1276 {
1277   // Return the "Chi2" to be minimized with Minuit for track fitting,
1278   // with "NParam" parameters
1279   // and their current values in array pointed to by "Param".
1280   // Assumes that the track hits are sorted according to increasing Z.
1281   // Track parameters at each TrackHit are updated accordingly.
1282   // Multiple Coulomb scattering is taken into account with covariance matrix.
1283   AliMUONTrack *trackBeingFitted;
1284   AliMUONTrackParam param1;
1285   AliMUONTrackParam* TrackParamAtHit;
1286   AliMUONHitForRec* HitForRec;
1287   Chi2 = 0.0; // initialize Chi2
1288   // copy of track parameters to be fitted
1289   trackBeingFitted = (AliMUONTrack*) AliMUONTrackReconstructor::Fitter()->GetObjectFit();
1290   // 3 parameters means fit track candidate from parameters at vertex (X_vertex and Y_vertex are fixed)
1291   if (NParam == 3) param1 = *(trackBeingFitted->GetTrackParamAtVertex());
1292   else param1 = *((AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->First()));
1293   // Minuit parameters to be fitted into this copy
1294   param1.SetInverseBendingMomentum(Param[0]);
1295   param1.SetBendingSlope(Param[1]);
1296   param1.SetNonBendingSlope(Param[2]);
1297   if (NParam == 5) {
1298     param1.SetNonBendingCoor(Param[3]);
1299     param1.SetBendingCoor(Param[4]);
1300   }
1301
1302   Int_t chCurrent, chPrev = 0, hitNumber, hitNumber1, hitNumber2, hitNumber3;
1303   Double_t z1, z2, z3;
1304   AliMUONTrackParam *TrackParamAtHit1, *TrackParamAtHit2, *TrackParamAtHit3;
1305   AliMUONHitForRec *HitForRec1, *HitForRec2;
1306   Double_t hbc1, hbc2, pbc1, pbc2;
1307   Double_t hnbc1, hnbc2, pnbc1, pnbc2;
1308   Int_t numberOfHit = trackBeingFitted->GetNTrackHits();
1309   TMatrixD *covBending = new TMatrixD(numberOfHit, numberOfHit);
1310   TMatrixD *covNonBending = new TMatrixD(numberOfHit, numberOfHit);
1311   Double_t *msa2 = new Double_t[numberOfHit];
1312
1313   // Predicted coordinates and  multiple scattering angles are first calculated
1314   for (hitNumber = 0; hitNumber < numberOfHit; hitNumber++) {
1315     TrackParamAtHit = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber));
1316     HitForRec = TrackParamAtHit->GetHitForRecPtr();
1317     // extrapolation to the plane of the HitForRec attached to the current TrackParamAtHit
1318     param1.ExtrapToZ(HitForRec->GetZ());
1319     // update track parameters of the current hit
1320     TrackParamAtHit->SetTrackParam(param1);
1321     // square of multiple scattering angle at current hit, with one chamber
1322     msa2[hitNumber] = MultipleScatteringAngle2(&param1);
1323     // correction for eventual missing hits or multiple hits in a chamber,
1324     // according to the number of chambers
1325     // between the current hit and the previous one
1326     chCurrent = HitForRec->GetChamberNumber();
1327     if (hitNumber > 0) msa2[hitNumber] = msa2[hitNumber] * (chCurrent - chPrev);
1328     chPrev = chCurrent;
1329   }
1330
1331   // Calculates the covariance matrix
1332   for (hitNumber1 = 0; hitNumber1 < numberOfHit; hitNumber1++) {
1333     TrackParamAtHit1 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber1));
1334     HitForRec1 = TrackParamAtHit1->GetHitForRecPtr();
1335     z1 = HitForRec1->GetZ();
1336     for (hitNumber2 = hitNumber1; hitNumber2 < numberOfHit; hitNumber2++) {
1337       TrackParamAtHit2 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber2));
1338       z2 = TrackParamAtHit2->GetHitForRecPtr()->GetZ();
1339       // initialization to 0 (diagonal plus upper triangular part)
1340       (*covBending)(hitNumber2, hitNumber1) = 0.0;
1341       // contribution from multiple scattering in bending plane:
1342       // loop over upstream hits
1343       for (hitNumber3 = 0; hitNumber3 < hitNumber1; hitNumber3++) {
1344         TrackParamAtHit3 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber3));
1345         z3 = TrackParamAtHit3->GetHitForRecPtr()->GetZ();
1346         (*covBending)(hitNumber2, hitNumber1) = (*covBending)(hitNumber2, hitNumber1) + ((z1 - z3) * (z2 - z3) * msa2[hitNumber3]);
1347       }
1348       // equal contribution from multiple scattering in non bending plane
1349       (*covNonBending)(hitNumber2, hitNumber1) = (*covBending)(hitNumber2, hitNumber1);
1350       if (hitNumber1 == hitNumber2) {
1351         // Diagonal elements: add contribution from position measurements
1352         // in bending plane
1353         (*covBending)(hitNumber2, hitNumber1) = (*covBending)(hitNumber2, hitNumber1) + HitForRec1->GetBendingReso2();
1354         // and in non bending plane
1355         (*covNonBending)(hitNumber2, hitNumber1) = (*covNonBending)(hitNumber2, hitNumber1) + HitForRec1->GetNonBendingReso2();
1356       } else {
1357         // Non diagonal elements: symmetrization
1358         // for bending plane
1359         (*covBending)(hitNumber1, hitNumber2) = (*covBending)(hitNumber2, hitNumber1);
1360         // and non bending plane
1361         (*covNonBending)(hitNumber1, hitNumber2) = (*covNonBending)(hitNumber2, hitNumber1);
1362       }
1363     } // for (hitNumber2 = hitNumber1;...
1364   } // for (hitNumber1 = 0;...
1365
1366   // Inversion of covariance matrices
1367   // with "mnvertLocal", local "mnvert" function of Minuit.
1368   // One cannot use directly "mnvert" since "TVirtualFitter" does not know it.
1369   // One will have to replace this local function by the right inversion function
1370   // from a specialized Root package for symmetric positive definite matrices,
1371   // when available!!!!
1372   Int_t ifailBending;
1373   mnvertLocal(&((*covBending)(0,0)), numberOfHit, numberOfHit, numberOfHit, ifailBending);
1374   Int_t ifailNonBending;
1375   mnvertLocal(&((*covNonBending)(0,0)), numberOfHit, numberOfHit, numberOfHit, ifailNonBending);
1376
1377   // It would be worth trying to calculate the inverse of the covariance matrix
1378   // only once per fit, since it cannot change much in principle,
1379   // and it would save a lot of computing time !!!!
1380
1381   // Calculates Chi2
1382   if ((ifailBending == 0) && (ifailNonBending == 0)) {
1383     // with Multiple Scattering if inversion correct
1384     for (hitNumber1 = 0; hitNumber1 < numberOfHit ; hitNumber1++) {
1385       TrackParamAtHit1 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber1));
1386       HitForRec1 = TrackParamAtHit1->GetHitForRecPtr();
1387       hbc1 = HitForRec1->GetBendingCoor();
1388       pbc1 = TrackParamAtHit1->GetBendingCoor();
1389       hnbc1 = HitForRec1->GetNonBendingCoor();
1390       pnbc1 = TrackParamAtHit1->GetNonBendingCoor();
1391       for (hitNumber2 = 0; hitNumber2 < numberOfHit; hitNumber2++) {
1392         TrackParamAtHit2 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber2));
1393         HitForRec2 = TrackParamAtHit2->GetHitForRecPtr();
1394         hbc2 = HitForRec2->GetBendingCoor();
1395         pbc2 = TrackParamAtHit2->GetBendingCoor();
1396         hnbc2 = HitForRec2->GetNonBendingCoor();
1397         pnbc2 = TrackParamAtHit2->GetNonBendingCoor();
1398         Chi2 += ((*covBending)(hitNumber2, hitNumber1) * (hbc1 - pbc1) * (hbc2 - pbc2)) +
1399                 ((*covNonBending)(hitNumber2, hitNumber1) * (hnbc1 - pnbc1) * (hnbc2 - pnbc2));
1400       }
1401     }
1402   } else {
1403     // without Multiple Scattering if inversion impossible
1404     for (hitNumber1 = 0; hitNumber1 < numberOfHit ; hitNumber1++) {
1405       TrackParamAtHit1 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber1));
1406       HitForRec1 = TrackParamAtHit1->GetHitForRecPtr();
1407       hbc1 = HitForRec1->GetBendingCoor();
1408       pbc1 = TrackParamAtHit1->GetBendingCoor();
1409       hnbc1 = HitForRec1->GetNonBendingCoor();
1410       pnbc1 = TrackParamAtHit1->GetNonBendingCoor();
1411       Chi2 += ((hbc1 - pbc1) * (hbc1 - pbc1) / HitForRec1->GetBendingReso2()) +
1412               ((hnbc1 - pnbc1) * (hnbc1 - pnbc1) / HitForRec1->GetNonBendingReso2());
1413     }
1414   }
1415
1416   delete covBending;
1417   delete covNonBending;
1418   delete [] msa2;
1419 }
1420
1421 Double_t MultipleScatteringAngle2(AliMUONTrackParam *param)
1422 {
1423   // Returns square of multiple Coulomb scattering angle
1424   // from TrackParamAtHit pointed to by "param"
1425   Double_t slopeBending, slopeNonBending, radiationLength, inverseBendingMomentum2, inverseTotalMomentum2;
1426   Double_t varMultipleScatteringAngle;
1427   // Better implementation in AliMUONTrack class ????
1428   slopeBending = param->GetBendingSlope();
1429   slopeNonBending = param->GetNonBendingSlope();
1430   // thickness in radiation length for the current track,
1431   // taking local angle into account
1432   radiationLength = AliMUONConstants::DefaultChamberThicknessInX0() *
1433                     TMath::Sqrt(1.0 + slopeBending*slopeBending + slopeNonBending*slopeNonBending);
1434   inverseBendingMomentum2 =  param->GetInverseBendingMomentum() * param->GetInverseBendingMomentum();
1435   inverseTotalMomentum2 = inverseBendingMomentum2 * (1.0 + slopeBending * slopeBending) /
1436                           (1.0 + slopeBending *slopeBending + slopeNonBending * slopeNonBending);
1437   varMultipleScatteringAngle = 0.0136 * (1.0 + 0.038 * TMath::Log(radiationLength));
1438   // The velocity is assumed to be 1 !!!!
1439   varMultipleScatteringAngle = inverseTotalMomentum2 * radiationLength * varMultipleScatteringAngle * varMultipleScatteringAngle;
1440   return varMultipleScatteringAngle;
1441 }
1442
1443 //______________________________________________________________________________
1444  void mnvertLocal(Double_t *a, Int_t l, Int_t, Int_t n, Int_t &ifail)
1445 {
1446 //*-*-*-*-*-*-*-*-*-*-*-*Inverts a symmetric matrix*-*-*-*-*-*-*-*-*-*-*-*-*
1447 //*-*                    ==========================
1448 //*-*        inverts a symmetric matrix.   matrix is first scaled to
1449 //*-*        have all ones on the diagonal (equivalent to change of units)
1450 //*-*        but no pivoting is done since matrix is positive-definite.
1451 //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
1452
1453   // taken from TMinuit package of Root (l>=n)
1454   // fVERTs, fVERTq and fVERTpp changed to localVERTs, localVERTq and localVERTpp
1455   //  Double_t localVERTs[n], localVERTq[n], localVERTpp[n];
1456   Double_t * localVERTs = new Double_t[n];
1457   Double_t * localVERTq = new Double_t[n];
1458   Double_t * localVERTpp = new Double_t[n];
1459   // fMaxint changed to localMaxint
1460   Int_t localMaxint = n;
1461
1462     /* System generated locals */
1463     Int_t aOffset;
1464
1465     /* Local variables */
1466     Double_t si;
1467     Int_t i, j, k, kp1, km1;
1468
1469     /* Parameter adjustments */
1470     aOffset = l + 1;
1471     a -= aOffset;
1472
1473     /* Function Body */
1474     ifail = 0;
1475     if (n < 1) goto L100;
1476     if (n > localMaxint) goto L100;
1477 //*-*-                  scale matrix by sqrt of diag elements
1478     for (i = 1; i <= n; ++i) {
1479         si = a[i + i*l];
1480         if (si <= 0) goto L100;
1481         localVERTs[i-1] = 1 / TMath::Sqrt(si);
1482     }
1483     for (i = 1; i <= n; ++i) {
1484         for (j = 1; j <= n; ++j) {
1485             a[i + j*l] = a[i + j*l]*localVERTs[i-1]*localVERTs[j-1];
1486         }
1487     }
1488 //*-*-                                       . . . start main loop . . . .
1489     for (i = 1; i <= n; ++i) {
1490         k = i;
1491 //*-*-                  preparation for elimination step1
1492         if (a[k + k*l] != 0) localVERTq[k-1] = 1 / a[k + k*l];
1493         else goto L100;
1494         localVERTpp[k-1] = 1;
1495         a[k + k*l] = 0;
1496         kp1 = k + 1;
1497         km1 = k - 1;
1498         if (km1 < 0) goto L100;
1499         else if (km1 == 0) goto L50;
1500         else               goto L40;
1501 L40:
1502         for (j = 1; j <= km1; ++j) {
1503             localVERTpp[j-1] = a[j + k*l];
1504             localVERTq[j-1]  = a[j + k*l]*localVERTq[k-1];
1505             a[j + k*l]   = 0;
1506         }
1507 L50:
1508         if (k - n < 0) goto L51;
1509         else if (k - n == 0) goto L60;
1510         else                goto L100;
1511 L51:
1512         for (j = kp1; j <= n; ++j) {
1513             localVERTpp[j-1] = a[k + j*l];
1514             localVERTq[j-1]  = -a[k + j*l]*localVERTq[k-1];
1515             a[k + j*l]   = 0;
1516         }
1517 //*-*-                  elimination proper
1518 L60:
1519         for (j = 1; j <= n; ++j) {
1520             for (k = j; k <= n; ++k) { a[j + k*l] += localVERTpp[j-1]*localVERTq[k-1]; }
1521         }
1522     }
1523 //*-*-                  elements of left diagonal and unscaling
1524     for (j = 1; j <= n; ++j) {
1525         for (k = 1; k <= j; ++k) {
1526             a[k + j*l] = a[k + j*l]*localVERTs[k-1]*localVERTs[j-1];
1527             a[j + k*l] = a[k + j*l];
1528         }
1529     }
1530     delete [] localVERTs;
1531     delete [] localVERTq;
1532     delete [] localVERTpp;
1533     return;
1534 //*-*-                  failure return
1535 L100:
1536     delete [] localVERTs;
1537     delete [] localVERTq;
1538     delete [] localVERTpp;
1539     ifail = 1;
1540 } /* mnvertLocal */
1541
1542   //__________________________________________________________________________
1543 void AliMUONTrackReconstructor::RemoveDoubleTracks(void)
1544 {
1545   // To remove double tracks.
1546   // Tracks are considered identical
1547   // if they have at least half of their hits in common.
1548   // Among two identical tracks, one keeps the track with the larger number of hits
1549   // or, if these numbers are equal, the track with the minimum Chi2.
1550   AliMUONTrack *track1, *track2, *trackToRemove;
1551   Int_t hitsInCommon, nHits1, nHits2;
1552   Bool_t removedTrack1;
1553   // Loop over first track of the pair
1554   track1 = (AliMUONTrack*) fRecTracksPtr->First();
1555   while (track1) {
1556     removedTrack1 = kFALSE;
1557     nHits1 = track1->GetNTrackHits();
1558     // Loop over second track of the pair
1559     track2 = (AliMUONTrack*) fRecTracksPtr->After(track1);
1560     while (track2) {
1561       nHits2 = track2->GetNTrackHits();
1562       // number of hits in common between two tracks
1563       hitsInCommon = track1->HitsInCommon(track2);
1564       // check for identical tracks
1565       if ((4 * hitsInCommon) >= (nHits1 + nHits2)) {
1566         // decide which track to remove
1567         if ((nHits1 > nHits2) || ((nHits1 == nHits2) && (track1->GetFitFMin() < track2->GetFitFMin()))) {
1568           // remove track2 and continue the second loop with the track next to track2
1569           trackToRemove = track2;
1570           track2 = (AliMUONTrack*) fRecTracksPtr->After(track2);
1571           fRecTracksPtr->Remove(trackToRemove);
1572           fNRecTracks--;
1573           fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards
1574         } else {
1575           // else remove track1 and continue the first loop with the track next to track1
1576           trackToRemove = track1;
1577           track1 = (AliMUONTrack*) fRecTracksPtr->After(track1);
1578           fRecTracksPtr->Remove(trackToRemove);
1579           fNRecTracks--;
1580           fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards
1581           removedTrack1 = kTRUE;
1582           break;
1583         }
1584       } else track2 = (AliMUONTrack*) fRecTracksPtr->After(track2);
1585     } // track2
1586     if (removedTrack1) continue;
1587     track1 = (AliMUONTrack*) fRecTracksPtr->After(track1);
1588   } // track1
1589   return;
1590 }
1591
1592   //__________________________________________________________________________
1593 void AliMUONTrackReconstructor::UpdateHitForRecAtHit()
1594 {
1595   // Set cluster parameters after track fitting. Fill fHitForRecAtHit of AliMUONTrack's
1596   AliMUONTrack *track;
1597   AliMUONTrackParam *trackParamAtHit;
1598   track = (AliMUONTrack*) fRecTracksPtr->First();
1599   while (track) {
1600     trackParamAtHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First());
1601     while (trackParamAtHit) {
1602       track->AddHitForRecAtHit(trackParamAtHit->GetHitForRecPtr());
1603       trackParamAtHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->After(trackParamAtHit));
1604     }
1605     track = (AliMUONTrack*) fRecTracksPtr->After(track);
1606   }
1607   return;
1608 }
1609
1610   //__________________________________________________________________________
1611 void AliMUONTrackReconstructor::FillMUONTrack()
1612 {
1613   // Set track parameters at hits for Kalman track. Fill fTrackParamAtHit of AliMUONTrack's
1614   AliMUONTrackK *track;
1615   track = (AliMUONTrackK*) fRecTracksPtr->First();
1616   while (track) {
1617     track->FillMUONTrack();
1618     track = (AliMUONTrackK*) fRecTracksPtr->After(track);
1619   }
1620   return;
1621 }
1622
1623   //__________________________________________________________________________
1624 void AliMUONTrackReconstructor::EventDump(void)
1625 {
1626   // Dump reconstructed event (track parameters at vertex and at first hit),
1627   // and the particle parameters
1628
1629   AliMUONTrack *track;
1630   AliMUONTrackParam *trackParam, *trackParam1;
1631   Double_t bendingSlope, nonBendingSlope, pYZ;
1632   Double_t pX, pY, pZ, x, y, z, c;
1633   Int_t trackIndex, nTrackHits;
1634
1635   AliDebug(1,"****** enter EventDump ******");
1636   AliDebug(1, Form("Number of Reconstructed tracks : %d", fNRecTracks));
1637
1638   fRecTracksPtr->Compress(); // for simple loop without "Next" since no hole
1639   // Loop over reconstructed tracks
1640   for (trackIndex = 0; trackIndex < fNRecTracks; trackIndex++) {
1641     if (fTrackMethod != 1) continue; //AZ - skip the rest for now
1642     AliDebug(1, Form("track number: %d", trackIndex));
1643     // function for each track for modularity ????
1644     track = (AliMUONTrack*) ((*fRecTracksPtr)[trackIndex]);
1645     nTrackHits = track->GetNTrackHits();
1646     AliDebug(1, Form("Number of track hits: %d ", nTrackHits));
1647     // track parameters at Vertex
1648     trackParam = track->GetTrackParamAtVertex();
1649     x = trackParam->GetNonBendingCoor();
1650     y = trackParam->GetBendingCoor();
1651     z = trackParam->GetZ();
1652     bendingSlope = trackParam->GetBendingSlope();
1653     nonBendingSlope = trackParam->GetNonBendingSlope();
1654     pYZ = 1/TMath::Abs(trackParam->GetInverseBendingMomentum());
1655     pZ = pYZ/TMath::Sqrt(1+bendingSlope*bendingSlope);
1656     pX = pZ * nonBendingSlope;
1657     pY = pZ * bendingSlope;
1658     c = TMath::Sign(1.0, trackParam->GetInverseBendingMomentum());
1659     AliDebug(1, Form("Track parameters at Vertex z= %f: X= %f Y= %f pX= %f pY= %f pZ= %f c= %f\n",
1660                      z, x, y, pX, pY, pZ, c));
1661
1662     // track parameters at first hit
1663     trackParam1 = (AliMUONTrackParam*) track->GetTrackParamAtHit()->First();
1664     x = trackParam1->GetNonBendingCoor();
1665     y = trackParam1->GetBendingCoor();
1666     z = trackParam1->GetZ();
1667     bendingSlope = trackParam1->GetBendingSlope();
1668     nonBendingSlope = trackParam1->GetNonBendingSlope();
1669     pYZ = 1/TMath::Abs(trackParam1->GetInverseBendingMomentum());
1670     pZ = pYZ/TMath::Sqrt(1.0 + bendingSlope * bendingSlope);
1671     pX = pZ * nonBendingSlope;
1672     pY = pZ * bendingSlope;
1673     c = TMath::Sign(1.0, trackParam1->GetInverseBendingMomentum());
1674     AliDebug(1, Form("track parameters at z= %f: X= %f Y= %f pX= %f pY= %f pZ= %f c= %f\n",
1675                      z, x, y, pX, pY, pZ, c));
1676   }
1677   // informations about generated particles NO !!!!!!!!
1678
1679 //    for (Int_t iPart = 0; iPart < np; iPart++) {
1680 //      p = gAlice->Particle(iPart);
1681 //      printf(" particle %d: type= %d px= %f py= %f pz= %f pdg= %d\n",
1682 //         iPart, p->GetPdgCode(), p->Px(), p->Py(), p->Pz(), p->GetPdgCode());
1683 //    }
1684   return;
1685 }
1686
1687
1688 //__________________________________________________________________________
1689 void AliMUONTrackReconstructor::EventDumpTrigger(void)
1690 {
1691   // Dump reconstructed trigger event
1692   // and the particle parameters
1693
1694   AliMUONTriggerTrack *triggertrack ;
1695   Int_t nTriggerTracks = fMUONData->RecTriggerTracks()->GetEntriesFast();
1696
1697   AliDebug(1, "****** enter EventDumpTrigger ******");
1698   AliDebug(1, Form("Number of Reconstructed tracks : %d ",  nTriggerTracks));
1699
1700   // Loop over reconstructed tracks
1701   for (Int_t trackIndex = 0; trackIndex < nTriggerTracks; trackIndex++) {
1702     triggertrack = (AliMUONTriggerTrack*)fMUONData->RecTriggerTracks()->At(trackIndex);
1703       printf(" trigger track number %i x11=%f y11=%f thetax=%f thetay=%f \n",
1704              trackIndex,
1705              triggertrack->GetX11(),triggertrack->GetY11(),
1706              triggertrack->GetThetax(),triggertrack->GetThetay());
1707   }
1708 }
1709
1710 //__________________________________________________________________________
1711 void AliMUONTrackReconstructor::MakeTrackCandidatesK(void)
1712 {
1713   // To make initial tracks for Kalman filter from the list of segments
1714   Int_t istat, iseg;
1715   AliMUONSegment *segment;
1716   AliMUONTrackK *trackK;
1717
1718   AliDebug(1,"Enter MakeTrackCandidatesK");
1719
1720   AliMUONTrackK a(this, fHitsForRecPtr);
1721   // Loop over stations(1...) 5 and 4
1722   for (istat=4; istat>=3; istat--) {
1723     // Loop over segments in the station
1724     for (iseg=0; iseg<fNSegments[istat]; iseg++) {
1725       // Transform segments to tracks and evaluate covariance matrix
1726       segment = (AliMUONSegment*) ((*fSegmentsPtr[istat])[iseg]);
1727       trackK = new ((*fRecTracksPtr)[fNRecTracks++]) AliMUONTrackK(segment);
1728     } // for (iseg=0;...)
1729   } // for (istat=4;...)
1730   return;
1731 }
1732
1733 //__________________________________________________________________________
1734 void AliMUONTrackReconstructor::FollowTracksK(void)
1735 {
1736   // Follow tracks using Kalman filter
1737   Bool_t ok;
1738   Int_t icand, ichamBeg = 0, ichamEnd, chamBits;
1739   Double_t zDipole1, zDipole2;
1740   AliMUONTrackK *trackK;
1741   AliMUONHitForRec *hit;
1742   AliMUONRawCluster *clus;
1743   TClonesArray *rawclusters;
1744   clus = 0; rawclusters = 0;
1745
1746   zDipole1 = GetSimpleBPosition() + GetSimpleBLength()/2;
1747   zDipole2 = zDipole1 - GetSimpleBLength();
1748
1749   // Print hits
1750   trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[0]);
1751
1752   if (trackK->DebugLevel() > 0) {
1753     for (Int_t i1=0; i1<fNHitsForRec; i1++) {
1754       hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[i1]);
1755       printf(" Hit # %d %10.4f %10.4f %10.4f",
1756              hit->GetChamberNumber(), hit->GetBendingCoor(),
1757              hit->GetNonBendingCoor(), hit->GetZ());
1758
1759       // from raw clusters
1760       rawclusters = fMUONData->RawClusters(hit->GetChamberNumber());
1761       clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->
1762                                                            GetHitNumber());
1763       printf(" %d", clus->GetTrack(1));
1764       if (clus->GetTrack(2) != -1) printf(" %d \n", clus->GetTrack(2));
1765       else printf("\n");
1766
1767     }
1768   } // if (trackK->DebugLevel() > 0)
1769
1770   icand = -1;
1771   Int_t nSeeds;
1772   nSeeds = fNRecTracks; // starting number of seeds
1773   // Loop over track candidates
1774   while (icand < fNRecTracks-1) {
1775     icand ++;
1776     if (trackK->DebugLevel()>0) cout << " *** Kalman track candidate No. " << icand << endl;
1777     trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[icand]);
1778     if (trackK->GetRecover() < 0) continue; // failed track
1779
1780     // Discard candidate which will produce the double track
1781     /*
1782     if (icand > 0) {
1783       ok = CheckCandidateK(icand,nSeeds);
1784       if (!ok) {
1785         trackK->SetRecover(-1); // mark candidate to be removed
1786         continue;
1787       }
1788     }
1789     */
1790
1791     ok = kTRUE;
1792     if (trackK->GetRecover() == 0)
1793       hit = (AliMUONHitForRec*) trackK->GetTrackHits()->Last(); // last hit
1794     else
1795       hit = trackK->GetHitLastOk(); // hit where track stopped
1796
1797     if (hit) ichamBeg = hit->GetChamberNumber();
1798     ichamEnd = 0;
1799     // Check propagation direction
1800     if (!hit) { ichamBeg = ichamEnd; AliFatal(" ??? "); }
1801     else if (trackK->GetTrackDir() < 0) {
1802       ichamEnd = 9; // forward propagation
1803       ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2);
1804       if (ok) {
1805         ichamBeg = ichamEnd;
1806         ichamEnd = 6; // backward propagation
1807         // Change weight matrix and zero fChi2 for backpropagation
1808         trackK->StartBack();
1809         trackK->SetTrackDir(1);
1810         ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kTRUE,zDipole1,zDipole2);
1811         ichamBeg = ichamEnd;
1812         ichamEnd = 0;
1813       }
1814     } else {
1815       if (trackK->GetBPFlag()) {
1816         // backpropagation
1817         ichamEnd = 6; // backward propagation
1818         // Change weight matrix and zero fChi2 for backpropagation
1819         trackK->StartBack();
1820         ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kTRUE,zDipole1,zDipole2);
1821         ichamBeg = ichamEnd;
1822         ichamEnd = 0;
1823       }
1824     }
1825
1826     if (ok) {
1827       trackK->SetTrackDir(1);
1828       trackK->SetBPFlag(kFALSE);
1829       ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2);
1830     }
1831     if (!ok) { trackK->SetRecover(-1); continue; } // mark candidate to be removed
1832
1833     // Apply smoother
1834     if (trackK->GetRecover() >= 0) {
1835       ok = trackK->Smooth();
1836       if (!ok) trackK->SetRecover(-1); // mark candidate to be removed
1837     }
1838
1839     // Majority 3 of 4 in first 2 stations
1840     if (!ok) continue;
1841     chamBits = 0;
1842     Double_t chi2max = 0;
1843     for (Int_t i=0; i<trackK->GetNTrackHits(); i++) {
1844       hit = (AliMUONHitForRec*) (*trackK->GetTrackHits())[i];
1845       chamBits |= BIT(hit->GetChamberNumber());
1846       if (trackK->GetChi2PerPoint(i) > chi2max) chi2max = trackK->GetChi2PerPoint(i);
1847     }
1848     if (!((chamBits&3)==3 || (chamBits>>2&3)==3) && chi2max > 25) {
1849       //trackK->Recover();
1850       trackK->SetRecover(-1); //mark candidate to be removed
1851       continue;
1852     }
1853     if (ok) trackK->SetTrackQuality(0); // compute "track quality"
1854   } // while
1855
1856   for (Int_t i=0; i<fNRecTracks; i++) {
1857     trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[i]);
1858     if (trackK->GetRecover() < 0) fRecTracksPtr->RemoveAt(i);
1859   }
1860
1861   // Compress TClonesArray
1862   fRecTracksPtr->Compress();
1863   fNRecTracks = fRecTracksPtr->GetEntriesFast();
1864   return;
1865 }
1866
1867 //__________________________________________________________________________
1868 Bool_t AliMUONTrackReconstructor::CheckCandidateK(Int_t icand, Int_t nSeeds) const
1869 {
1870   // Discards track candidate if it will produce the double track (having
1871   // the same seed segment hits as hits of a good track found before)
1872   AliMUONTrackK *track1, *track2;
1873   AliMUONHitForRec *hit1, *hit2, *hit;
1874
1875   track1 = (AliMUONTrackK*) ((*fRecTracksPtr)[icand]);
1876   hit1 = (AliMUONHitForRec*) (*track1->GetTrackHits())[0]; // 1'st hit
1877   hit2 = (AliMUONHitForRec*) (*track1->GetTrackHits())[1]; // 2'nd hit
1878
1879   for (Int_t i=0; i<icand; i++) {
1880     track2 = (AliMUONTrackK*) ((*fRecTracksPtr)[i]);
1881     //if (track2->GetRecover() < 0) continue;
1882     if (track2->GetRecover() < 0 && icand >= nSeeds) continue;
1883
1884     if (track1->GetStartSegment() == track2->GetStartSegment()) {
1885       return kFALSE;
1886     } else {
1887       Int_t nSame = 0;
1888       for (Int_t j=0; j<track2->GetNTrackHits(); j++) {
1889         hit = (AliMUONHitForRec*) (*track2->GetTrackHits())[j];
1890         if (hit == hit1 || hit == hit2) {
1891           nSame++;
1892           if (nSame == 2) return kFALSE;
1893         }
1894       } // for (Int_t j=0;
1895     }
1896   } // for (Int_t i=0;
1897   return kTRUE;
1898 }
1899
1900 //__________________________________________________________________________
1901 void AliMUONTrackReconstructor::RemoveDoubleTracksK(void)
1902 {
1903   // Removes double tracks (sharing more than half of their hits). Keeps
1904   // the track with higher quality
1905   AliMUONTrackK *track1, *track2, *trackToKill;
1906
1907   // Sort tracks according to their quality
1908   fRecTracksPtr->Sort();
1909
1910   // Loop over first track of the pair
1911   track1 = (AliMUONTrackK*) fRecTracksPtr->First();
1912   Int_t debug = track1->DebugLevel();
1913   while (track1) {
1914     // Loop over second track of the pair
1915     track2 = (AliMUONTrackK*) fRecTracksPtr->After(track1);
1916     while (track2) {
1917       // Check whether or not to keep track2
1918       if (!track2->KeepTrack(track1)) {
1919         if (debug >= 0) cout << " Killed track: " << 1/(*track2->GetTrackParameters())(4,0) <<
1920           " " << track2->GetTrackQuality() << endl;
1921         trackToKill = track2;
1922         track2 = (AliMUONTrackK*) fRecTracksPtr->After(track2);
1923         trackToKill->Kill();
1924         fRecTracksPtr->Compress();
1925       } else track2 = (AliMUONTrackK*) fRecTracksPtr->After(track2);
1926     } // track2
1927     track1 = (AliMUONTrackK*) fRecTracksPtr->After(track1);
1928   } // track1
1929
1930   fNRecTracks = fRecTracksPtr->GetEntriesFast();
1931   if (debug >= 0) cout << " Number of Kalman tracks: " << fNRecTracks << endl;
1932 }
1933
1934 //__________________________________________________________________________
1935 void AliMUONTrackReconstructor::GoToVertex(void)
1936 {
1937   // Propagates track to the vertex thru absorber
1938   // (using Branson correction for now)
1939
1940   Double_t zVertex;
1941   zVertex = 0;
1942   for (Int_t i=0; i<fNRecTracks; i++) {
1943     //((AliMUONTrackK*)(*fRecTracksPtr)[i])->Branson();
1944     ((AliMUONTrackK*)(*fRecTracksPtr)[i])->SetTrackQuality(1); // compute Chi2
1945     //((AliMUONTrackK*)(*fRecTracksPtr)[i])->GoToZ(zVertex); // w/out absorber
1946     ((AliMUONTrackK*)(*fRecTracksPtr)[i])->GoToVertex(1); // with absorber
1947   }
1948 }
1949
1950 //__________________________________________________________________________
1951 void AliMUONTrackReconstructor::SetTrackMethod(Int_t iTrackMethod)
1952 {
1953   // Set track method and recreate track container if necessary
1954
1955   fTrackMethod = TMath::Min (iTrackMethod, 3);
1956   fTrackMethod = TMath::Max (fTrackMethod, 1);
1957   if (fTrackMethod != 1) {
1958     if (fRecTracksPtr) delete fRecTracksPtr;
1959     fRecTracksPtr = new TClonesArray("AliMUONTrackK", 10);
1960     if (fTrackMethod == 2) cout << " *** Tracking with the Kalman filter *** " << endl;
1961     else cout << " *** Combined cluster / track finder ***" << endl;
1962   } else cout << " *** Original tracking *** " << endl;
1963
1964 }