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 "AliMUONTrackHit.h"
  53 #include "AliMagF.h"
  54 #include "AliMUONTrackK.h"
  55 #include "AliLog.h"
  56 #include "AliTracker.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 const Double_t AliMUONTrackReconstructor::fgkDefaultChamberThicknessInX0 = 0.03;
  66 // Simple magnetic field:
  67 // Value taken from macro MUONtracking.C: 0.7 T, hence 7 kG
  68 // Length and Position from reco_muon.F, with opposite sign:
  69 // Length = ZMAGEND-ZCOIL
  70 // Position = (ZMAGEND+ZCOIL)/2
  71 // to be ajusted differently from real magnetic field ????
  72 const Double_t AliMUONTrackReconstructor::fgkDefaultSimpleBValue = 7.0;
  73 const Double_t AliMUONTrackReconstructor::fgkDefaultSimpleBLength = 428.0;
  74 const Double_t AliMUONTrackReconstructor::fgkDefaultSimpleBPosition = 1019.0;
  75 const Double_t AliMUONTrackReconstructor::fgkDefaultEfficiency = 0.95;
  76
  77 ClassImp(AliMUONTrackReconstructor) // Class implementation in ROOT context
  78
  79 //__________________________________________________________________________
  80 AliMUONTrackReconstructor::AliMUONTrackReconstructor(AliMUONData* data)
  81   : TObject(),
  82     fTrackMethod(1), //AZ - tracking method (1-default, 2-Kalman)
  83     fMinBendingMomentum(fgkDefaultMinBendingMomentum),
  84     fMaxBendingMomentum(fgkDefaultMaxBendingMomentum),
  85     fMaxChi2(fgkDefaultMaxChi2),
  86     fMaxSigma2Distance(fgkDefaultMaxSigma2Distance),
  87     fBendingResolution(fgkDefaultBendingResolution),
  88     fNonBendingResolution(fgkDefaultNonBendingResolution),
  89     fChamberThicknessInX0(fgkDefaultChamberThicknessInX0),
  90     fSimpleBValue(fgkDefaultSimpleBValue),
  91     fSimpleBLength(fgkDefaultSimpleBLength),
  92     fSimpleBPosition(fgkDefaultSimpleBPosition),
  93     fEfficiency(fgkDefaultEfficiency),
  94     fHitsForRecPtr(0x0),
  95     fNHitsForRec(0),
  96     fRecTracksPtr(0x0),
  97     fNRecTracks(0),
  98     fRecTrackHitsPtr(0x0),
  99     fNRecTrackHits(0),
 100     fMUONData(data),
 101     fMuons(0),
 102     fTriggerTrack(new AliMUONTriggerTrack()),
 103     fTriggerCircuit(0x0)
 104 {
 105   // Constructor for class AliMUONTrackReconstructor
 106   SetReconstructionParametersToDefaults();
 107
 108   // Memory allocation for the TClonesArray of hits for reconstruction
 109   // Is 10000 the right size ????
 110   fHitsForRecPtr = new TClonesArray("AliMUONHitForRec", 10000);
 111
 112   // Memory allocation for the TClonesArray's of segments in stations
 113   // Is 2000 the right size ????
 114   for (Int_t st = 0; st < AliMUONConstants::NTrackingCh()/2; st++) {
 115     fSegmentsPtr[st] = new TClonesArray("AliMUONSegment", 2000);
 116     fNSegments[st] = 0; // really needed or GetEntriesFast sufficient ????
 117   }
 118   // Memory allocation for the TClonesArray of reconstructed tracks
 119   // Is 10 the right size ????
 120   fRecTracksPtr = new TClonesArray("AliMUONTrack", 10);
 121
 122   // Memory allocation for the TClonesArray of hits on reconstructed tracks
 123   // Is 100 the right size ????
 124   fRecTrackHitsPtr = new TClonesArray("AliMUONTrackHit", 100);
 125
 126   const AliMagF* kField = AliTracker::GetFieldMap();
 127   if (!kField) AliFatal("No field available");
 128  // Sign of fSimpleBValue according to sign of Bx value at (50,50,-950).
 129   Float_t b[3], x[3];
 130   x[0] = 50.; x[1] = 50.; x[2] = -950.;
 131   kField->Field(x,b);
 132
 133   fSimpleBValue    = TMath::Sign(fSimpleBValue,(Double_t) b[0]);
 134   fSimpleBPosition = TMath::Sign(fSimpleBPosition,(Double_t) x[2]);
 135   // See how to get fSimple(BValue, BLength, BPosition)
 136   // automatically calculated from the actual magnetic field ????
 137
 138   return;
 139 }
 140
 141   //__________________________________________________________________________
 142 AliMUONTrackReconstructor::~AliMUONTrackReconstructor(void)
 143 {
 144   // Destructor for class AliMUONTrackReconstructor
 145   delete fHitsForRecPtr; // Correct destruction of everything ???? or delete [] ????
 146   for (Int_t st = 0; st < AliMUONConstants::NTrackingCh()/2; st++)
 147     delete fSegmentsPtr[st]; // Correct destruction of everything ????
 148
 149   delete fTriggerTrack;
 150 }
 151   //__________________________________________________________________________
 152 void AliMUONTrackReconstructor::SetReconstructionParametersToDefaults(void)
 153 {
 154   // Set reconstruction parameters to default values
 155   // Would be much more convenient with a structure (or class) ????
 156
 157   // ******** Parameters for making HitsForRec
 158   // minimum radius,
 159   // like in TRACKF_STAT:
 160   // 2 degrees for stations 1 and 2, or ch(0...) from 0 to 3;
 161   // 30 cm for stations 3 to 5, or ch(0...) from 4 to 9
 162   for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
 163     if (ch < 4) fRMin[ch] = TMath::Abs(AliMUONConstants::DefaultChamberZ(ch)) *
 164                   2.0 * TMath::Pi() / 180.0;
 165     else fRMin[ch] = 30.0;
 166     // maximum radius at 10 degrees and Z of chamber
 167     fRMax[ch] = TMath::Abs(AliMUONConstants::DefaultChamberZ(ch)) *
 168                   10.0 * TMath::Pi() / 180.0;
 169   }
 170
 171   // ******** Parameters for making segments
 172   // should be parametrized ????
 173   // according to interval between chambers in a station ????
 174   // Maximum distance in non bending plane
 175   // 5 * 0.22 just to remember the way it was made in TRACKF_STAT
 176   // SIGCUT*DYMAX(IZ)
 177   for (Int_t st = 0; st < AliMUONConstants::NTrackingCh()/2; st++)
 178     fSegmentMaxDistNonBending[st] = 5. * 0.22;
 179   // Maximum distance in bending plane:
 180   // values from TRACKF_STAT, corresponding to (J psi 20cm),
 181   // scaled to the real distance between chambers in a station
 182   fSegmentMaxDistBending[0] = TMath::Abs( 1.5 *
 183                                           (AliMUONConstants::DefaultChamberZ(1) - AliMUONConstants::DefaultChamberZ(0)) / 20.0);
 184   fSegmentMaxDistBending[1] = TMath::Abs( 1.5 *
 185                                           (AliMUONConstants::DefaultChamberZ(3) - AliMUONConstants::DefaultChamberZ(2)) / 20.0);
 186   fSegmentMaxDistBending[2] = TMath::Abs( 3.0 *
 187                                           (AliMUONConstants::DefaultChamberZ(5) - AliMUONConstants::DefaultChamberZ(4)) / 20.0);
 188   fSegmentMaxDistBending[3] = TMath::Abs( 6.0 *
 189                                           (AliMUONConstants::DefaultChamberZ(7) - AliMUONConstants::DefaultChamberZ(6)) / 20.0);
 190   fSegmentMaxDistBending[4] = TMath::Abs( 6.0 *
 191                                           (AliMUONConstants::DefaultChamberZ(9) - AliMUONConstants::DefaultChamberZ(8)) / 20.0);
 192
 193   return;
 194 }
 195
 196 //__________________________________________________________________________
 197 Double_t AliMUONTrackReconstructor::GetImpactParamFromBendingMomentum(Double_t BendingMomentum) const
 198 {
 199   // Returns impact parameter at vertex in bending plane (cm),
 200   // from the signed bending momentum "BendingMomentum" in bending plane (GeV/c),
 201   // using simple values for dipole magnetic field.
 202   // The sign of "BendingMomentum" is the sign of the charge.
 203   return (-0.0003 * fSimpleBValue * fSimpleBLength * fSimpleBPosition /
 204           BendingMomentum);
 205 }
 206
 207 //__________________________________________________________________________
 208 Double_t AliMUONTrackReconstructor::GetBendingMomentumFromImpactParam(Double_t ImpactParam) const
 209 {
 210   // Returns signed bending momentum in bending plane (GeV/c),
 211   // the sign being the sign of the charge for particles moving forward in Z,
 212   // from the impact parameter "ImpactParam" at vertex in bending plane (cm),
 213   // using simple values for dipole magnetic field.
 214   return (-0.0003 * fSimpleBValue * fSimpleBLength * fSimpleBPosition /
 215           ImpactParam);
 216 }
 217
 218 //__________________________________________________________________________
 219 void AliMUONTrackReconstructor::EventReconstruct(void)
 220 {
 221   // To reconstruct one event
 222   AliDebug(1,"Enter EventReconstruct");
 223   ResetTracks(); //AZ
 224   ResetTrackHits(); //AZ
 225   ResetSegments(); //AZ
 226   ResetHitsForRec(); //AZ
 227   MakeEventToBeReconstructed();
 228   MakeSegments();
 229   MakeTracks();
 230   if (fMUONData->IsTriggerTrackBranchesInTree())
 231     ValidateTracksWithTrigger();
 232
 233   // Add tracks to MUON data container
 234   for(Int_t i=0; i<GetNRecTracks(); i++) {
 235     AliMUONTrack * track = (AliMUONTrack*) GetRecTracksPtr()->At(i);
 236     fMUONData->AddRecTrack(*track);
 237   }
 238
 239   return;
 240 }
 241
 242 //__________________________________________________________________________
 243 void AliMUONTrackReconstructor::EventReconstructTrigger(void)
 244 {
 245   // To reconstruct one event
 246   AliDebug(1,"Enter EventReconstructTrigger");
 247   MakeTriggerTracks();
 248   return;
 249 }
 250
 251   //__________________________________________________________________________
 252 void AliMUONTrackReconstructor::ResetHitsForRec(void)
 253 {
 254   // To reset the array and the number of HitsForRec,
 255   // and also the number of HitsForRec
 256   // and the index of the first HitForRec per chamber
 257   if (fHitsForRecPtr) fHitsForRecPtr->Delete();
 258   fNHitsForRec = 0;
 259   for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++)
 260     fNHitsForRecPerChamber[ch] = fIndexOfFirstHitForRecPerChamber[ch] = 0;
 261   return;
 262 }
 263
 264   //__________________________________________________________________________
 265 void AliMUONTrackReconstructor::ResetSegments(void)
 266 {
 267   // To reset the TClonesArray of segments and the number of Segments
 268   // for all stations
 269   for (Int_t st = 0; st < AliMUONConstants::NTrackingCh()/2; st++) {
 270     if (fSegmentsPtr[st]) fSegmentsPtr[st]->Clear();
 271     fNSegments[st] = 0;
 272   }
 273   return;
 274 }
 275
 276   //__________________________________________________________________________
 277 void AliMUONTrackReconstructor::ResetTracks(void)
 278 {
 279   // To reset the TClonesArray of reconstructed tracks
 280   if (fRecTracksPtr) fRecTracksPtr->Delete();
 281   // Delete in order that the Track destructors are called,
 282   // hence the space for the TClonesArray of pointers to TrackHit's is freed
 283   fNRecTracks = 0;
 284   return;
 285 }
 286
 287
 288   //__________________________________________________________________________
 289 void AliMUONTrackReconstructor::ResetTrackHits(void)
 290 {
 291   // To reset the TClonesArray of hits on reconstructed tracks
 292   if (fRecTrackHitsPtr) fRecTrackHitsPtr->Clear();
 293   fNRecTrackHits = 0;
 294   return;
 295 }
 296
 297   //__________________________________________________________________________
 298 void AliMUONTrackReconstructor::MakeEventToBeReconstructed(void)
 299 {
 300   // To make the list of hits to be reconstructed,
 301   // either from the track ref. hits or from the raw clusters
 302   // according to the parameter set for the reconstructor
 303
 304   AliDebug(1,"Enter MakeEventToBeReconstructed");
 305   //AZ ResetHitsForRec();
 306
 307   // Reconstruction from raw clusters
 308   // AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
 309   // Security on MUON ????
 310   // TreeR assumed to be be "prepared" in calling function
 311   // by "MUON->GetTreeR(nev)" ????
 312   TTree *treeR = fMUONData->TreeR();
 313
 314   //AZ? fMUONData->SetTreeAddress("RC");
 315   AddHitsForRecFromRawClusters(treeR);
 316   // No sorting: it is done automatically in the previous function
 317
 318
 319   AliDebug(1,"End of MakeEventToBeReconstructed");
 320     if (AliLog::GetGlobalDebugLevel() > 0) {
 321       AliDebug(1, Form("NHitsForRec: %d",fNHitsForRec));
 322       for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
 323         AliDebug(1, Form("Chamber(0...): %d",ch));
 324         AliDebug(1, Form("NHitsForRec: %d", fNHitsForRecPerChamber[ch]));
 325         AliDebug(1, Form("Index(first HitForRec): %d", fIndexOfFirstHitForRecPerChamber[ch]));
 326         for (Int_t hit = fIndexOfFirstHitForRecPerChamber[ch];
 327              hit < fIndexOfFirstHitForRecPerChamber[ch] + fNHitsForRecPerChamber[ch];
 328              hit++) {
 329           AliDebug(1, Form("HitForRec index(0...): %d",hit));
 330           ((*fHitsForRecPtr)[hit])->Dump();
 331       }
 332     }
 333   }
 334   return;
 335 }
 336
 337   //__________________________________________________________________________
 338 void AliMUONTrackReconstructor::SortHitsForRecWithIncreasingChamber()
 339 {
 340   // Sort HitsForRec's in increasing order with respect to chamber number.
 341   // Uses the function "Compare".
 342   // Update the information for HitsForRec per chamber too.
 343   Int_t ch, nhits, prevch;
 344   fHitsForRecPtr->Sort();
 345   for (ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
 346     fNHitsForRecPerChamber[ch] = 0;
 347     fIndexOfFirstHitForRecPerChamber[ch] = 0;
 348   }
 349   prevch = 0; // previous chamber
 350   nhits = 0; // number of hits in current chamber
 351   // Loop over HitsForRec
 352   for (Int_t hit = 0; hit < fNHitsForRec; hit++) {
 353     // chamber number (0...)
 354     ch = ((AliMUONHitForRec*)  ((*fHitsForRecPtr)[hit]))->GetChamberNumber();
 355     // increment number of hits in current chamber
 356     (fNHitsForRecPerChamber[ch])++;
 357     // update index of first HitForRec in current chamber
 358     // if chamber number different from previous one
 359     if (ch != prevch) {
 360       fIndexOfFirstHitForRecPerChamber[ch] = hit;
 361       prevch = ch;
 362     }
 363   }
 364   return;
 365 }
 366
 367   //__________________________________________________________________________
 368 void AliMUONTrackReconstructor::AddHitsForRecFromRawClusters(TTree* TR)
 369 {
 370   // To add to the list of hits for reconstruction all the raw clusters
 371   // No condition added, like in Fortran TRACKF_STAT,
 372   // on the radius between RMin and RMax.
 373   AliMUONHitForRec *hitForRec;
 374   AliMUONRawCluster *clus;
 375   Int_t iclus, nclus, nTRentries;
 376   TClonesArray *rawclusters;
 377   AliDebug(1,"Enter AddHitsForRecFromRawClusters");
 378
 379   if (fTrackMethod != 3) { //AZ
 380     fMUONData->SetTreeAddress("RC"); //AZ
 381     nTRentries = Int_t(TR->GetEntries());
 382     if (nTRentries != 1) {
 383       AliError(Form("nTRentries = %d not equal to 1 ",nTRentries));
 384       exit(0);
 385     }
 386     fMUONData->GetRawClusters(); // only one entry
 387   }
 388
 389   // Loop over tracking chambers
 390   for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
 391     // number of HitsForRec to 0 for the chamber
 392     fNHitsForRecPerChamber[ch] = 0;
 393     // index of first HitForRec for the chamber
 394     if (ch == 0) fIndexOfFirstHitForRecPerChamber[ch] = 0;
 395     else fIndexOfFirstHitForRecPerChamber[ch] = fNHitsForRec;
 396     rawclusters =fMUONData->RawClusters(ch);
 397     nclus = (Int_t) (rawclusters->GetEntries());
 398     // Loop over (cathode correlated) raw clusters
 399     for (iclus = 0; iclus < nclus; iclus++) {
 400       clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(iclus);
 401       // new AliMUONHitForRec from raw cluster
 402       // and increment number of AliMUONHitForRec's (total and in chamber)
 403       hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(clus);
 404       fNHitsForRec++;
 405       (fNHitsForRecPerChamber[ch])++;
 406       // more information into HitForRec
 407       //  resolution: info should be already in raw cluster and taken from it ????
 408       //hitForRec->SetBendingReso2(fBendingResolution * fBendingResolution);
 409       //hitForRec->SetNonBendingReso2(fNonBendingResolution * fNonBendingResolution);
 410       hitForRec->SetBendingReso2(clus->GetErrY() * clus->GetErrY());
 411       hitForRec->SetNonBendingReso2(clus->GetErrX() * clus->GetErrX());
 412       //  original raw cluster
 413       hitForRec->SetChamberNumber(ch);
 414       hitForRec->SetHitNumber(iclus);
 415       // Z coordinate of the raw cluster (cm)
 416       hitForRec->SetZ(clus->GetZ(0));
 417       if (AliLog::GetGlobalDebugLevel() > 1) {
 418         cout << "chamber (0...): " << ch <<
 419           " raw cluster (0...): " << iclus << endl;
 420         clus->Dump();
 421         cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl;
 422         hitForRec->Dump();}
 423     } // end of cluster loop
 424   } // end of chamber loop
 425   SortHitsForRecWithIncreasingChamber();
 426   return;
 427 }
 428
 429   //__________________________________________________________________________
 430 void AliMUONTrackReconstructor::MakeSegments(void)
 431 {
 432   // To make the list of segments in all stations,
 433   // from the list of hits to be reconstructed
 434   AliDebug(1,"Enter MakeSegments");
 435   //AZ ResetSegments();
 436   // Loop over stations
 437   Int_t nb = (fTrackMethod != 1) ? 3 : 0; //AZ
 438   for (Int_t st = nb; st < AliMUONConstants::NTrackingCh()/2; st++)
 439     MakeSegmentsPerStation(st);
 440   if (AliLog::GetGlobalDebugLevel() > 1) {
 441     cout << "end of MakeSegments" << endl;
 442     for (Int_t st = 0; st < AliMUONConstants::NTrackingCh()/2; st++) {
 443       cout << "station(0...): " << st
 444            << "  Segments: " << fNSegments[st]
 445            << endl;
 446       for (Int_t seg = 0;
 447            seg < fNSegments[st];
 448            seg++) {
 449         cout << "Segment index(0...): " << seg << endl;
 450         ((*fSegmentsPtr[st])[seg])->Dump();
 451       }
 452     }
 453   }
 454   return;
 455 }
 456
 457   //__________________________________________________________________________
 458 void AliMUONTrackReconstructor::MakeSegmentsPerStation(Int_t Station)
 459 {
 460   // To make the list of segments in station number "Station" (0...)
 461   // from the list of hits to be reconstructed.
 462   // Updates "fNSegments"[Station].
 463   // Segments in stations 4 and 5 are sorted
 464   // according to increasing absolute value of "impact parameter"
 465   AliMUONHitForRec *hit1Ptr, *hit2Ptr;
 466   AliMUONSegment *segment;
 467   Bool_t last2st;
 468   Double_t bendingSlope, distBend, distNonBend, extBendCoor, extNonBendCoor,
 469       impactParam = 0., maxImpactParam = 0., minImpactParam = 0.; // =0 to avoid compilation warnings.
 470   AliDebug(1,Form("Enter MakeSegmentsPerStation (0...) %d",Station));
 471   // first and second chambers (0...) in the station
 472   Int_t ch1 = 2 * Station;
 473   Int_t ch2 = ch1 + 1;
 474   // variable true for stations downstream of the dipole:
 475   // Station(0..4) equal to 3 or 4
 476   if ((Station == 3) || (Station == 4)) {
 477     last2st = kTRUE;
 478     // maximum impact parameter (cm) according to fMinBendingMomentum...
 479     maxImpactParam =
 480       TMath::Abs(GetImpactParamFromBendingMomentum(fMinBendingMomentum));
 481     // minimum impact parameter (cm) according to fMaxBendingMomentum...
 482     minImpactParam =
 483       TMath::Abs(GetImpactParamFromBendingMomentum(fMaxBendingMomentum));
 484   }
 485   else last2st = kFALSE;
 486   // extrapolation factor from Z of first chamber to Z of second chamber
 487   // dZ to be changed to take into account fine structure of chambers ????
 488   Double_t extrapFact;
 489   // index for current segment
 490   Int_t segmentIndex = 0;
 491   // Loop over HitsForRec in the first chamber of the station
 492   for (Int_t hit1 = fIndexOfFirstHitForRecPerChamber[ch1];
 493        hit1 < fIndexOfFirstHitForRecPerChamber[ch1] + fNHitsForRecPerChamber[ch1];
 494        hit1++) {
 495     // pointer to the HitForRec
 496     hit1Ptr = (AliMUONHitForRec*) ((*fHitsForRecPtr)[hit1]);
 497     // extrapolation,
 498     // on the straight line joining the HitForRec to the vertex (0,0,0),
 499     // to the Z of the second chamber of the station
 500     // Loop over HitsForRec in the second chamber of the station
 501     for (Int_t hit2 = fIndexOfFirstHitForRecPerChamber[ch2];
 502          hit2 < fIndexOfFirstHitForRecPerChamber[ch2] + fNHitsForRecPerChamber[ch2];
 503          hit2++) {
 504       // pointer to the HitForRec
 505       hit2Ptr = (AliMUONHitForRec*) ((*fHitsForRecPtr)[hit2]);
 506       // absolute values of distances, in bending and non bending planes,
 507       // between the HitForRec in the second chamber
 508       // and the previous extrapolation
 509       extrapFact = hit2Ptr->GetZ()/ hit1Ptr->GetZ();
 510       extBendCoor = extrapFact * hit1Ptr->GetBendingCoor();
 511       extNonBendCoor = extrapFact * hit1Ptr->GetNonBendingCoor();
 512       distBend = TMath::Abs(hit2Ptr->GetBendingCoor() - extBendCoor);
 513       distNonBend = TMath::Abs(hit2Ptr->GetNonBendingCoor() - extNonBendCoor);
 514       if (last2st) {
 515         // bending slope
 516         if ( hit1Ptr->GetZ() - hit2Ptr->GetZ() != 0.0 ) {
 517           bendingSlope = (hit1Ptr->GetBendingCoor() - hit2Ptr->GetBendingCoor()) /
 518             (hit1Ptr->GetZ() - hit2Ptr->GetZ());
 519           // absolute value of impact parameter
 520           impactParam =
 521             TMath::Abs(hit1Ptr->GetBendingCoor() - hit1Ptr->GetZ() * bendingSlope);
 522          }
 523          else {
 524            AliWarning("hit1Ptr->GetZ() = hit2Ptr->GetZ(): impactParam set to maxImpactParam");
 525            impactParam = maxImpactParam;
 526          }
 527       }
 528       // check for distances not too large,
 529       // and impact parameter not too big if stations downstream of the dipole.
 530       // Conditions "distBend" and "impactParam" correlated for these stations ????
 531       if ((distBend < fSegmentMaxDistBending[Station]) &&
 532           (distNonBend < fSegmentMaxDistNonBending[Station]) &&
 533           (!last2st || (impactParam < maxImpactParam)) &&
 534           (!last2st || (impactParam > minImpactParam))) {
 535         // make new segment
 536         segment = new ((*fSegmentsPtr[Station])[segmentIndex])
 537           AliMUONSegment(hit1Ptr, hit2Ptr);
 538         // update "link" to this segment from the hit in the first chamber
 539         if (hit1Ptr->GetNSegments() == 0)
 540           hit1Ptr->SetIndexOfFirstSegment(segmentIndex);
 541         hit1Ptr->SetNSegments(hit1Ptr->GetNSegments() + 1);
 542         if (AliLog::GetGlobalDebugLevel() > 1) {
 543           cout << "segmentIndex(0...): " << segmentIndex
 544                << "  distBend: " << distBend
 545                << "  distNonBend: " << distNonBend
 546                << endl;
 547           segment->Dump();
 548           cout << "HitForRec in first chamber" << endl;
 549           hit1Ptr->Dump();
 550           cout << "HitForRec in second chamber" << endl;
 551           hit2Ptr->Dump();
 552         };
 553         // increment index for current segment
 554         segmentIndex++;
 555       }
 556     } //for (Int_t hit2
 557   } // for (Int_t hit1...
 558   fNSegments[Station] = segmentIndex;
 559   // Sorting according to "impact parameter" if station(1..5) 4 or 5,
 560   // i.e. Station(0..4) 3 or 4, using the function "Compare".
 561   // After this sorting, it is impossible to use
 562   // the "fNSegments" and "fIndexOfFirstSegment"
 563   // of the HitForRec in the first chamber to explore all segments formed with it.
 564   // Is this sorting really needed ????
 565   if ((Station == 3) || (Station == 4)) (fSegmentsPtr[Station])->Sort();
 566   AliDebug(1,Form("Station: %d  NSegments:  %d ", Station, fNSegments[Station]));
 567   return;
 568 }
 569
 570   //__________________________________________________________________________
 571 void AliMUONTrackReconstructor::MakeTracks(void)
 572 {
 573   // To make the tracks,
 574   // from the list of segments and points in all stations
 575   AliDebug(1,"Enter MakeTracks");
 576   // The order may be important for the following Reset's
 577   //AZ ResetTracks();
 578   //AZ ResetTrackHits();
 579   if (fTrackMethod != 1) { //AZ - Kalman filter
 580     MakeTrackCandidatesK();
 581     if (fRecTracksPtr->GetEntriesFast() == 0) return;
 582     // Follow tracks in stations(1..) 3, 2 and 1
 583     FollowTracksK();
 584     // Remove double tracks
 585     RemoveDoubleTracksK();
 586     // Propagate tracks to the vertex thru absorber
 587     GoToVertex();
 588     // Fill AliMUONTrack data members
 589     FillMUONTrack();
 590   } else {
 591     // Look for candidates from at least 3 aligned points in stations(1..) 4 and 5
 592     MakeTrackCandidates();
 593     // Follow tracks in stations(1..) 3, 2 and 1
 594     FollowTracks();
 595     // Remove double tracks
 596     RemoveDoubleTracks();
 597     UpdateTrackParamAtHit();
 598     UpdateHitForRecAtHit();
 599   }
 600   return;
 601 }
 602
 603   //__________________________________________________________________________
 604 void AliMUONTrackReconstructor::ValidateTracksWithTrigger(void)
 605 {
 606   // Try to match track from tracking system with trigger track
 607   AliMUONTrack *track;
 608   TClonesArray *recTriggerTracks;
 609
 610   fMUONData->SetTreeAddress("RL");
 611   fMUONData->GetRecTriggerTracks();
 612   recTriggerTracks = fMUONData->RecTriggerTracks();
 613
 614   track = (AliMUONTrack*) fRecTracksPtr->First();
 615   while (track) {
 616     track->MatchTriggerTrack(recTriggerTracks);
 617     track = (AliMUONTrack*) fRecTracksPtr->After(track);
 618   }
 619
 620 }
 621
 622   //__________________________________________________________________________
 623 Bool_t AliMUONTrackReconstructor::MakeTriggerTracks(void)
 624 {
 625     // To make the trigger tracks from Local Trigger
 626   AliDebug(1, "Enter MakeTriggerTracks");
 627
 628     Int_t nTRentries;
 629     Long_t gloTrigPat;
 630     TClonesArray *localTrigger;
 631     TClonesArray *globalTrigger;
 632     AliMUONLocalTrigger *locTrg;
 633     AliMUONGlobalTrigger *gloTrg;
 634
 635     TTree* treeR = fMUONData->TreeR();
 636
 637     nTRentries = Int_t(treeR->GetEntries());
 638
 639     treeR->GetEvent(0); // only one entry
 640
 641     if (!(fMUONData->IsTriggerBranchesInTree())) {
 642       AliWarning(Form("Trigger information is not avalaible, nTRentries = %d not equal to 1",nTRentries));
 643       return kFALSE;
 644     }
 645
 646     fMUONData->SetTreeAddress("TC");
 647     fMUONData->GetTrigger();
 648
 649     // global trigger for trigger pattern
 650     gloTrigPat = 0;
 651     globalTrigger = fMUONData->GlobalTrigger();
 652     gloTrg = (AliMUONGlobalTrigger*)globalTrigger->UncheckedAt(0);
 653
 654     if (gloTrg)
 655       gloTrigPat = gloTrg->GetGlobalPattern();
 656
 657
 658     // local trigger for tracking
 659     localTrigger = fMUONData->LocalTrigger();
 660     Int_t nlocals = (Int_t) (localTrigger->GetEntries());
 661
 662     Float_t z11 = AliMUONConstants::DefaultChamberZ(10);
 663     Float_t z21 = AliMUONConstants::DefaultChamberZ(12);
 664
 665     Float_t y11 = 0.;
 666     Int_t stripX21 = 0;
 667     Float_t y21 = 0.;
 668     Float_t x11 = 0.;
 669
 670     for (Int_t i=0; i<nlocals; i++) { // loop on Local Trigger
 671       locTrg = (AliMUONLocalTrigger*)localTrigger->UncheckedAt(i);
 672
 673       AliDebug(1, "AliMUONTrackReconstructor::MakeTriggerTrack using NEW trigger \n");
 674       AliMUONTriggerCircuitNew* circuit =
 675         (AliMUONTriggerCircuitNew*)fTriggerCircuit->At(locTrg->LoCircuit()-1); // -1 !!!
 676
 677       y11 = circuit->GetY11Pos(locTrg->LoStripX());
 678       stripX21 = locTrg->LoStripX()+locTrg->LoDev()+1;
 679       y21 = circuit->GetY21Pos(stripX21);
 680       x11 = circuit->GetX11Pos(locTrg->LoStripY());
 681
 682       AliDebug(1, Form(" MakeTriggerTrack %d %d %d %d %d %f %f %f \n",i,locTrg->LoCircuit(),
 683                        locTrg->LoStripX(),locTrg->LoStripX()+locTrg->LoDev()+1,locTrg->LoStripY(),y11, y21, x11));
 684
 685       Float_t thetax = TMath::ATan2( x11 , z11 );
 686       Float_t thetay = TMath::ATan2( (y21-y11) , (z21-z11) );
 687
 688       fTriggerTrack->SetX11(x11);
 689       fTriggerTrack->SetY11(y11);
 690       fTriggerTrack->SetThetax(thetax);
 691       fTriggerTrack->SetThetay(thetay);
 692       fTriggerTrack->SetGTPattern(gloTrigPat);
 693
 694       fMUONData->AddRecTriggerTrack(*fTriggerTrack);
 695     } // end of loop on Local Trigger
 696     return kTRUE;
 697 }
 698
 699   //__________________________________________________________________________
 700 Int_t AliMUONTrackReconstructor::MakeTrackCandidatesWithTwoSegments(AliMUONSegment *BegSegment)
 701 {
 702   // To make track candidates with two segments in stations(1..) 4 and 5,
 703   // the first segment being pointed to by "BegSegment".
 704   // Returns the number of such track candidates.
 705   Int_t endStation, iEndSegment, nbCan2Seg;
 706   AliMUONSegment *endSegment;
 707   AliMUONSegment *extrapSegment = NULL;
 708   AliMUONTrack *recTrack;
 709   Double_t mcsFactor;
 710   AliDebug(1,"Enter MakeTrackCandidatesWithTwoSegments");
 711   // Station for the end segment
 712   endStation = 7 - (BegSegment->GetHitForRec1())->GetChamberNumber() / 2;
 713   // multiple scattering factor corresponding to one chamber
 714   mcsFactor = 0.0136 /
 715     GetBendingMomentumFromImpactParam(BegSegment->GetBendingImpact());
 716   mcsFactor     = fChamberThicknessInX0 * mcsFactor * mcsFactor;
 717   // linear extrapolation to end station
 718   // number of candidates with 2 segments to 0
 719   nbCan2Seg = 0;
 720   // Loop over segments in the end station
 721   for (iEndSegment = 0; iEndSegment < fNSegments[endStation]; iEndSegment++) {
 722     // pointer to segment
 723     endSegment = (AliMUONSegment*) ((*fSegmentsPtr[endStation])[iEndSegment]);
 724     // test compatibility between current segment and "extrapSegment"
 725     // 4 because 4 quantities in chi2
 726     extrapSegment =
 727       BegSegment->CreateSegmentFromLinearExtrapToStation(endSegment->GetZ(), mcsFactor);
 728     if ((endSegment->
 729          NormalizedChi2WithSegment(extrapSegment,
 730                                    fMaxSigma2Distance)) <= 4.0) {
 731       // both segments compatible:
 732       // make track candidate from "begSegment" and "endSegment"
 733       AliDebug(2,Form("TrackCandidate with Segment %d in Station(0..) %d", iEndSegment, endStation));
 734       // flag for both segments in one track:
 735       // to be done in track constructor ????
 736       BegSegment->SetInTrack(kTRUE);
 737       endSegment->SetInTrack(kTRUE);
 738       recTrack = new ((*fRecTracksPtr)[fNRecTracks])
 739         AliMUONTrack(BegSegment, endSegment, this);
 740       fNRecTracks++;
 741       if (AliLog::GetGlobalDebugLevel() > 1) recTrack->RecursiveDump();
 742       // increment number of track candidates with 2 segments
 743       nbCan2Seg++;
 744     }
 745     delete extrapSegment; // should not delete HitForRec's it points to !!!!
 746   } // for (iEndSegment = 0;...
 747   return nbCan2Seg;
 748 }
 749
 750   //__________________________________________________________________________
 751 Int_t AliMUONTrackReconstructor::MakeTrackCandidatesWithOneSegmentAndOnePoint(AliMUONSegment *BegSegment)
 752 {
 753   // To make track candidates with one segment and one point
 754   // in stations(1..) 4 and 5,
 755   // the segment being pointed to by "BegSegment".
 756   Int_t ch, ch1, ch2, endStation, iHit, iHitMax, iHitMin, nbCan1Seg1Hit;
 757   AliMUONHitForRec *extrapHitForRec= NULL;
 758   AliMUONHitForRec *hit;
 759   AliMUONTrack *recTrack;
 760   Double_t mcsFactor;
 761   AliDebug(1,"Enter MakeTrackCandidatesWithOneSegmentAndOnePoint");
 762   // station for the end point
 763   endStation = 7 - (BegSegment->GetHitForRec1())->GetChamberNumber() / 2;
 764   // multiple scattering factor corresponding to one chamber
 765   mcsFactor = 0.0136 /
 766     GetBendingMomentumFromImpactParam(BegSegment->GetBendingImpact());
 767   mcsFactor     = fChamberThicknessInX0 * mcsFactor * mcsFactor;
 768   // first and second chambers(0..) in the end station
 769   ch1 = 2 * endStation;
 770   ch2 = ch1 + 1;
 771   // number of candidates to 0
 772   nbCan1Seg1Hit = 0;
 773   // Loop over chambers of the end station
 774   for (ch = ch2; ch >= ch1; ch--) {
 775     // limits for the hit index in the loop
 776     iHitMin = fIndexOfFirstHitForRecPerChamber[ch];
 777     iHitMax = iHitMin + fNHitsForRecPerChamber[ch];
 778     // Loop over HitForRec's in the chamber
 779     for (iHit = iHitMin; iHit < iHitMax; iHit++) {
 780       // pointer to HitForRec
 781       hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[iHit]);
 782       // test compatibility between current HitForRec and "extrapHitForRec"
 783       // 2 because 2 quantities in chi2
 784       // linear extrapolation to chamber
 785       extrapHitForRec =
 786         BegSegment->CreateHitForRecFromLinearExtrapToChamber( hit->GetZ(), mcsFactor);
 787       if ((hit->
 788            NormalizedChi2WithHitForRec(extrapHitForRec,
 789                                        fMaxSigma2Distance)) <= 2.0) {
 790         // both HitForRec's compatible:
 791         // make track candidate from begSegment and current HitForRec
 792         AliDebug(2, Form("TrackCandidate with HitForRec  %d in Chamber(0..) %d", iHit, ch));
 793         // flag for beginning segments in one track:
 794         // to be done in track constructor ????
 795         BegSegment->SetInTrack(kTRUE);
 796         recTrack = new ((*fRecTracksPtr)[fNRecTracks])
 797           AliMUONTrack(BegSegment, hit, this);
 798         // the right place to eliminate "double counting" ???? how ????
 799         fNRecTracks++;
 800         if (AliLog::GetGlobalDebugLevel() > 1) recTrack->RecursiveDump();
 801         // increment number of track candidates
 802         nbCan1Seg1Hit++;
 803       }
 804       delete extrapHitForRec;
 805     } // for (iHit = iHitMin;...
 806   } // for (ch = ch2;...
 807   return nbCan1Seg1Hit;
 808 }
 809
 810   //__________________________________________________________________________
 811 void AliMUONTrackReconstructor::MakeTrackCandidates(void)
 812 {
 813   // To make track candidates
 814   // with at least 3 aligned points in stations(1..) 4 and 5
 815   // (two Segment's or one Segment and one HitForRec)
 816   Int_t begStation, iBegSegment, nbCan1Seg1Hit, nbCan2Seg;
 817   AliMUONSegment *begSegment;
 818   AliDebug(1,"Enter MakeTrackCandidates");
 819   // Loop over stations(1..) 5 and 4 for the beginning segment
 820   for (begStation = 4; begStation > 2; begStation--) {
 821     // Loop over segments in the beginning station
 822     for (iBegSegment = 0; iBegSegment < fNSegments[begStation]; iBegSegment++) {
 823       // pointer to segment
 824       begSegment = (AliMUONSegment*) ((*fSegmentsPtr[begStation])[iBegSegment]);
 825       AliDebug(2,Form("Look for TrackCandidate's with Segment %d  in Station(0..) %d", iBegSegment, begStation));
 826       // Look for track candidates with two segments,
 827       // "begSegment" and all compatible segments in other station.
 828       // Only for beginning station(1..) 5
 829       // because candidates with 2 segments have to looked for only once.
 830       if (begStation == 4)
 831         nbCan2Seg = MakeTrackCandidatesWithTwoSegments(begSegment);
 832       // Look for track candidates with one segment and one point,
 833       // "begSegment" and all compatible HitForRec's in other station.
 834       // Only if "begSegment" does not belong already to a track candidate.
 835       // Is that a too strong condition ????
 836       if (!(begSegment->GetInTrack()))
 837         nbCan1Seg1Hit = MakeTrackCandidatesWithOneSegmentAndOnePoint(begSegment);
 838     } // for (iBegSegment = 0;...
 839   } // for (begStation = 4;...
 840   return;
 841 }
 842
 843   //__________________________________________________________________________
 844 void AliMUONTrackReconstructor::FollowTracks(void)
 845 {
 846   // Follow tracks in stations(1..) 3, 2 and 1
 847   // too long: should be made more modular !!!!
 848   AliMUONHitForRec *bestHit, *extrapHit, *hit;
 849   AliMUONSegment *bestSegment, *extrapSegment, *segment;
 850   AliMUONTrack *track, *nextTrack;
 851   AliMUONTrackParam *trackParam1, trackParam[2], trackParamVertex;
 852   // -1 to avoid compilation warnings
 853   Int_t ch = -1, chInStation, chBestHit = -1, iHit, iSegment, station, trackIndex;
 854   Double_t bestChi2, chi2, dZ1, dZ2, dZ3, maxSigma2Distance, mcsFactor;
 855   Double_t bendingMomentum, chi2Norm = 0.;
 856
 857
 858   // local maxSigma2Distance, for easy increase in testing
 859   maxSigma2Distance = fMaxSigma2Distance;
 860   AliDebug(2,"Enter FollowTracks");
 861   // Loop over track candidates
 862   track = (AliMUONTrack*) fRecTracksPtr->First();
 863   trackIndex = -1;
 864   while (track) {
 865     // Follow function for each track candidate ????
 866     trackIndex++;
 867     nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track); // prepare next track
 868     AliDebug(2,Form("FollowTracks: track candidate(0..): %d", trackIndex));
 869     // Fit track candidate
 870     track->SetFitMCS(0); // without multiple Coulomb scattering
 871     track->SetFitNParam(3); // with 3 parameters (X = Y = 0)
 872     track->SetFitStart(0); // from parameters at vertex
 873     track->Fit();
 874     if (AliLog::GetGlobalDebugLevel()> 2) {
 875       cout << "FollowTracks: track candidate(0..): " << trackIndex
 876            << " after fit in stations(0..) 3 and 4" << endl;
 877       track->RecursiveDump();
 878     }
 879     // Loop over stations(1..) 3, 2 and 1
 880     // something SPECIAL for stations 2 and 1 for majority 3 coincidence ????
 881     // otherwise: majority coincidence 2 !!!!
 882     for (station = 2; station >= 0; station--) {
 883       // Track parameters at first track hit (smallest Z)
 884       trackParam1 = ((AliMUONTrackHit*)
 885                      (track->GetTrackHitsPtr()->First()))->GetTrackParam();
 886       // extrapolation to station
 887       trackParam1->ExtrapToStation(station, trackParam);
 888       extrapSegment = new AliMUONSegment(); //  empty segment
 889       // multiple scattering factor corresponding to one chamber
 890       // and momentum in bending plane (not total)
 891       mcsFactor = 0.0136 * trackParam1->GetInverseBendingMomentum();
 892       mcsFactor = fChamberThicknessInX0 * mcsFactor * mcsFactor;
 893       // Z difference from previous station
 894       dZ1 = AliMUONConstants::DefaultChamberZ(2 * station) -
 895             AliMUONConstants::DefaultChamberZ(2 * station + 2);
 896       // Z difference between the two previous stations
 897       dZ2 = AliMUONConstants::DefaultChamberZ(2 * station + 2) -
 898             AliMUONConstants::DefaultChamberZ(2 * station + 4);
 899       // Z difference between the two chambers in the previous station
 900       dZ3 = AliMUONConstants::DefaultChamberZ(2 * station) -
 901             AliMUONConstants::DefaultChamberZ(2 * station + 1);
 902       extrapSegment->SetBendingCoorReso2(fBendingResolution * fBendingResolution);
 903       extrapSegment->
 904         SetNonBendingCoorReso2(fNonBendingResolution * fNonBendingResolution);
 905       extrapSegment->UpdateFromStationTrackParam
 906         (trackParam, mcsFactor, dZ1, dZ2, dZ3, station,
 907          trackParam1->GetInverseBendingMomentum());
 908       bestChi2 = 5.0;
 909       bestSegment = NULL;
 910       if (AliLog::GetGlobalDebugLevel() > 2) {
 911         cout << "FollowTracks: track candidate(0..): " << trackIndex
 912              << " Look for segment in station(0..): " << station << endl;
 913       }
 914
 915       // Loop over segments in station
 916       for (iSegment = 0; iSegment < fNSegments[station]; iSegment++) {
 917         // Look for best compatible Segment in station
 918         // should consider all possibilities ????
 919         // multiple scattering ????
 920         // separation in 2 functions: Segment and HitForRec ????
 921         segment = (AliMUONSegment*) ((*fSegmentsPtr[station])[iSegment]);
 922         // correction of corrected segment (fBendingCoor and fNonBendingCoor)
 923         // according to real Z value of "segment" and slopes of "extrapSegment"
 924         (&(trackParam[0]))->ExtrapToZ(segment->GetZ());
 925         (&(trackParam[1]))->ExtrapToZ(segment->GetZ());
 926         extrapSegment->SetBendingCoor((&(trackParam[0]))->GetBendingCoor());
 927         extrapSegment->SetNonBendingCoor((&(trackParam[0]))->GetNonBendingCoor());
 928         extrapSegment->SetBendingSlope((&(trackParam[0]))->GetBendingSlope());
 929         extrapSegment->SetNonBendingSlope((&(trackParam[0]))->GetNonBendingSlope());
 930         chi2 = segment->
 931           NormalizedChi2WithSegment(extrapSegment, maxSigma2Distance);
 932         if (chi2 < bestChi2) {
 933           // update best Chi2 and Segment if better found
 934           bestSegment = segment;
 935           bestChi2 = chi2;
 936         }
 937       }
 938       if (bestSegment) {
 939         // best segment found: add it to track candidate
 940         (&(trackParam[0]))->ExtrapToZ(bestSegment->GetZ());
 941         (&(trackParam[1]))->ExtrapToZ(bestSegment->GetZ());
 942         track->AddSegment(bestSegment);
 943         // set track parameters at these two TrakHit's
 944         track->SetTrackParamAtHit(track->GetNTrackHits() - 2, &(trackParam[0]));
 945         track->SetTrackParamAtHit(track->GetNTrackHits() - 1, &(trackParam[1]));
 946         AliDebug(3, Form("FollowTracks: track candidate(0..): %d  Added segment in station(0..): %d", trackIndex, station));
 947         if (AliLog::GetGlobalDebugLevel()>2) track->RecursiveDump();
 948       }
 949       else {
 950         // No best segment found:
 951         // Look for best compatible HitForRec in station:
 952         // should consider all possibilities ????
 953         // multiple scattering ???? do about like for extrapSegment !!!!
 954         extrapHit = new AliMUONHitForRec(); //  empty hit
 955         bestChi2 = 3.0;
 956         bestHit = NULL;
 957         AliDebug(3, Form("FollowTracks: track candidate(0..): %d Segment not found, look for hit in station(0..): %d ",
 958                          trackIndex, station));
 959
 960         // Loop over chambers of the station
 961         for (chInStation = 0; chInStation < 2; chInStation++) {
 962           ch = 2 * station + chInStation;
 963           for (iHit = fIndexOfFirstHitForRecPerChamber[ch];
 964                iHit < fIndexOfFirstHitForRecPerChamber[ch] +
 965                  fNHitsForRecPerChamber[ch];
 966                iHit++) {
 967             hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[iHit]);
 968             // coordinates of extrapolated hit
 969             (&(trackParam[chInStation]))->ExtrapToZ(hit->GetZ());
 970             extrapHit->
 971               SetBendingCoor((&(trackParam[chInStation]))->GetBendingCoor());
 972             extrapHit->
 973               SetNonBendingCoor((&(trackParam[chInStation]))->GetNonBendingCoor());
 974             // resolutions from "extrapSegment"
 975             extrapHit->SetBendingReso2(extrapSegment->GetBendingCoorReso2());
 976             extrapHit->SetNonBendingReso2(extrapSegment->GetNonBendingCoorReso2());
 977             // Loop over hits in the chamber
 978             // condition for hit not already in segment ????
 979             chi2 = hit->NormalizedChi2WithHitForRec(extrapHit, maxSigma2Distance);
 980             if (chi2 < bestChi2) {
 981               // update best Chi2 and HitForRec if better found
 982               bestHit = hit;
 983               bestChi2 = chi2;
 984               chBestHit = chInStation;
 985             }
 986           }
 987         }
 988         if (bestHit) {
 989           // best hit found: add it to track candidate
 990           (&(trackParam[chBestHit]))->ExtrapToZ(bestHit->GetZ());
 991           track->AddHitForRec(bestHit);
 992           // set track parameters at this TrackHit
 993           track->SetTrackParamAtHit(track->GetNTrackHits() - 1,
 994                                     &(trackParam[chBestHit]));
 995           if (AliLog::GetGlobalDebugLevel() > 2) {
 996             cout << "FollowTracks: track candidate(0..): " << trackIndex
 997                  << " Added hit in station(0..): " << station << endl;
 998             track->RecursiveDump();
 999           }
1000         }
1001         else {
1002           // Remove current track candidate
1003           // and corresponding TrackHit's, ...
1004           track->Remove();
1005           delete extrapSegment;
1006           delete extrapHit;
1007           break; // stop the search for this candidate:
1008           // exit from the loop over station
1009         }
1010         delete extrapHit;
1011       }
1012       delete extrapSegment;
1013       // Sort track hits according to increasing Z
1014       track->GetTrackHitsPtr()->Sort();
1015       // Update track parameters at first track hit (smallest Z)
1016       trackParam1 = ((AliMUONTrackHit*)
1017                      (track->GetTrackHitsPtr()->First()))->GetTrackParam();
1018       bendingMomentum = 0.;
1019       if (TMath::Abs(trackParam1->GetInverseBendingMomentum()) > 0.)
1020         bendingMomentum = TMath::Abs(1/(trackParam1->GetInverseBendingMomentum()));
1021       // Track removed if bendingMomentum not in window [min, max]
1022       if ((bendingMomentum < fMinBendingMomentum) || (bendingMomentum > fMaxBendingMomentum)) {
1023         track->Remove();
1024         break; // stop the search for this candidate:
1025         // exit from the loop over station
1026       }
1027       // Track fit
1028       // with multiple Coulomb scattering if all stations
1029       if (station == 0) track->SetFitMCS(1);
1030       // without multiple Coulomb scattering if not all stations
1031       else track->SetFitMCS(0);
1032       track->SetFitNParam(5);  // with 5 parameters (momentum and position)
1033       track->SetFitStart(1);  // from parameters at first hit
1034       track->Fit();
1035       Double_t numberOfDegFree = (2.0 * track->GetNTrackHits() - 5);
1036       if (numberOfDegFree > 0) {
1037         chi2Norm =  track->GetFitFMin() / numberOfDegFree;
1038       } else {
1039         chi2Norm = 1.e10;
1040       }
1041       // Track removed if normalized chi2 too high
1042       if (chi2Norm > fMaxChi2) {
1043         track->Remove();
1044         break; // stop the search for this candidate:
1045         // exit from the loop over station
1046       }
1047       if (AliLog::GetGlobalDebugLevel() > 2) {
1048         cout << "FollowTracks: track candidate(0..): " << trackIndex
1049              << " after fit from station(0..): " << station << " to 4" << endl;
1050         track->RecursiveDump();
1051       }
1052       // Track extrapolation to the vertex through the absorber (Branson)
1053       // after going through the first station
1054       if (station == 0) {
1055         trackParamVertex = *trackParam1;
1056         (&trackParamVertex)->ExtrapToVertex(0.,0.,0.);
1057         track->SetTrackParamAtVertex(&trackParamVertex);
1058         if (AliLog::GetGlobalDebugLevel() > 0) {
1059           cout << "FollowTracks: track candidate(0..): " << trackIndex
1060                << " after extrapolation to vertex" << endl;
1061           track->RecursiveDump();
1062         }
1063       }
1064     } // for (station = 2;...
1065     // go really to next track
1066     track = nextTrack;
1067   } // while (track)
1068   // Compression of track array (necessary after Remove ????)
1069   fRecTracksPtr->Compress();
1070   return;
1071 }
1072
1073   //__________________________________________________________________________
1074 void AliMUONTrackReconstructor::RemoveDoubleTracks(void)
1075 {
1076   // To remove double tracks.
1077   // Tracks are considered identical
1078   // if they have at least half of their hits in common.
1079   // Among two identical tracks, one keeps the track with the larger number of hits
1080   // or, if these numbers are equal, the track with the minimum Chi2.
1081   AliMUONTrack *track1, *track2, *trackToRemove;
1082   Bool_t identicalTracks;
1083   Int_t hitsInCommon, nHits1, nHits2;
1084   identicalTracks = kTRUE;
1085   while (identicalTracks) {
1086     identicalTracks = kFALSE;
1087     // Loop over first track of the pair
1088     track1 = (AliMUONTrack*) fRecTracksPtr->First();
1089     while (track1 && (!identicalTracks)) {
1090       nHits1 = track1->GetNTrackHits();
1091       // Loop over second track of the pair
1092       track2 = (AliMUONTrack*) fRecTracksPtr->After(track1);
1093       while (track2 && (!identicalTracks)) {
1094         nHits2 = track2->GetNTrackHits();
1095         // number of hits in common between two tracks
1096         hitsInCommon = track1->HitsInCommon(track2);
1097         // check for identical tracks
1098         if ((4 * hitsInCommon) >= (nHits1 + nHits2)) {
1099           identicalTracks = kTRUE;
1100           // decide which track to remove
1101           if (nHits1 > nHits2) trackToRemove = track2;
1102           else if (nHits1 < nHits2) trackToRemove = track1;
1103           else if ((track1->GetFitFMin()) < (track2->GetFitFMin()))
1104             trackToRemove = track2;
1105           else trackToRemove = track1;
1106           // remove it
1107           trackToRemove->Remove();
1108         }
1109         track2 = (AliMUONTrack*) fRecTracksPtr->After(track2);
1110       } // track2
1111       track1 = (AliMUONTrack*) fRecTracksPtr->After(track1);
1112     } // track1
1113   }
1114   return;
1115 }
1116
1117   //__________________________________________________________________________
1118 void AliMUONTrackReconstructor::UpdateTrackParamAtHit()
1119 {
1120   // Set track parameters after track fitting. Fill fTrackParamAtHit of AliMUONTrack's
1121   AliMUONTrack *track;
1122   AliMUONTrackHit *trackHit;
1123   AliMUONTrackParam *trackParam;
1124   track = (AliMUONTrack*) fRecTracksPtr->First();
1125   while (track) {
1126     trackHit = (AliMUONTrackHit*) (track->GetTrackHitsPtr())->First();
1127     while (trackHit) {
1128       trackParam = trackHit->GetTrackParam();
1129       track->AddTrackParamAtHit(trackParam);
1130       trackHit = (AliMUONTrackHit*) (track->GetTrackHitsPtr())->After(trackHit);
1131     } // trackHit
1132     track = (AliMUONTrack*) fRecTracksPtr->After(track);
1133   } // track
1134   return;
1135 }
1136
1137   //__________________________________________________________________________
1138 void AliMUONTrackReconstructor::UpdateHitForRecAtHit()
1139 {
1140   // Set cluster parameterss after track fitting. Fill fHitForRecAtHit of AliMUONTrack's
1141   AliMUONTrack *track;
1142   AliMUONTrackHit *trackHit;
1143   AliMUONHitForRec *hitForRec;
1144   track = (AliMUONTrack*) fRecTracksPtr->First();
1145   while (track) {
1146     trackHit = (AliMUONTrackHit*) (track->GetTrackHitsPtr())->First();
1147     while (trackHit) {
1148       hitForRec = trackHit->GetHitForRecPtr();
1149       track->AddHitForRecAtHit(hitForRec);
1150       trackHit = (AliMUONTrackHit*) (track->GetTrackHitsPtr())->After(trackHit);
1151     } // trackHit
1152     track = (AliMUONTrack*) fRecTracksPtr->After(track);
1153   } // track
1154   return;
1155 }
1156
1157   //__________________________________________________________________________
1158 void AliMUONTrackReconstructor::FillMUONTrack()
1159 {
1160   // Set track parameters at hits for Kalman track. Fill fTrackParamAtHit of AliMUONTrack's
1161   AliMUONTrackK *track;
1162   track = (AliMUONTrackK*) fRecTracksPtr->First();
1163   while (track) {
1164     track->FillMUONTrack();
1165     track = (AliMUONTrackK*) fRecTracksPtr->After(track);
1166   }
1167   return;
1168 }
1169
1170   //__________________________________________________________________________
1171 void AliMUONTrackReconstructor::EventDump(void)
1172 {
1173   // Dump reconstructed event (track parameters at vertex and at first hit),
1174   // and the particle parameters
1175
1176   AliMUONTrack *track;
1177   AliMUONTrackParam *trackParam, *trackParam1;
1178   Double_t bendingSlope, nonBendingSlope, pYZ;
1179   Double_t pX, pY, pZ, x, y, z, c;
1180   Int_t trackIndex, nTrackHits;
1181
1182   AliDebug(1,"****** enter EventDump ******");
1183   AliDebug(1, Form("Number of Reconstructed tracks : %d", fNRecTracks));
1184
1185   fRecTracksPtr->Compress(); // for simple loop without "Next" since no hole
1186   // Loop over reconstructed tracks
1187   for (trackIndex = 0; trackIndex < fNRecTracks; trackIndex++) {
1188     if (fTrackMethod != 1) continue; //AZ - skip the rest for now
1189     AliDebug(1, Form("track number: %d", trackIndex));
1190     // function for each track for modularity ????
1191     track = (AliMUONTrack*) ((*fRecTracksPtr)[trackIndex]);
1192     nTrackHits = track->GetNTrackHits();
1193     AliDebug(1, Form("Number of track hits: %d ", nTrackHits));
1194     // track parameters at Vertex
1195     trackParam = track->GetTrackParamAtVertex();
1196     x = trackParam->GetNonBendingCoor();
1197     y = trackParam->GetBendingCoor();
1198     z = trackParam->GetZ();
1199     bendingSlope = trackParam->GetBendingSlope();
1200     nonBendingSlope = trackParam->GetNonBendingSlope();
1201     pYZ = 1/TMath::Abs(trackParam->GetInverseBendingMomentum());
1202     pZ = pYZ/TMath::Sqrt(1+bendingSlope*bendingSlope);
1203     pX = pZ * nonBendingSlope;
1204     pY = pZ * bendingSlope;
1205     c = TMath::Sign(1.0, trackParam->GetInverseBendingMomentum());
1206     AliDebug(1, Form("Track parameters at Vertex z= %f: X= %f Y= %f pX= %f pY= %f pZ= %f c= %f\n",
1207                      z, x, y, pX, pY, pZ, c));
1208
1209     // track parameters at first hit
1210     trackParam1 = ((AliMUONTrackHit*)
1211                    (track->GetTrackHitsPtr()->First()))->GetTrackParam();
1212     x = trackParam1->GetNonBendingCoor();
1213     y = trackParam1->GetBendingCoor();
1214     z = trackParam1->GetZ();
1215     bendingSlope = trackParam1->GetBendingSlope();
1216     nonBendingSlope = trackParam1->GetNonBendingSlope();
1217     pYZ = 1/TMath::Abs(trackParam1->GetInverseBendingMomentum());
1218     pZ = pYZ/TMath::Sqrt(1.0 + bendingSlope * bendingSlope);
1219     pX = pZ * nonBendingSlope;
1220     pY = pZ * bendingSlope;
1221     c = TMath::Sign(1.0, trackParam1->GetInverseBendingMomentum());
1222     AliDebug(1, Form("track parameters at z= %f: X= %f Y= %f pX= %f pY= %f pZ= %f c= %f\n",
1223                      z, x, y, pX, pY, pZ, c));
1224   }
1225   // informations about generated particles NO !!!!!!!!
1226
1227 //    for (Int_t iPart = 0; iPart < np; iPart++) {
1228 //      p = gAlice->Particle(iPart);
1229 //      printf(" particle %d: type= %d px= %f py= %f pz= %f pdg= %d\n",
1230 //         iPart, p->GetPdgCode(), p->Px(), p->Py(), p->Pz(), p->GetPdgCode());
1231 //    }
1232   return;
1233 }
1234
1235
1236 //__________________________________________________________________________
1237 void AliMUONTrackReconstructor::EventDumpTrigger(void)
1238 {
1239   // Dump reconstructed trigger event
1240   // and the particle parameters
1241
1242   AliMUONTriggerTrack *triggertrack ;
1243   Int_t nTriggerTracks = fMUONData->RecTriggerTracks()->GetEntriesFast();
1244
1245   AliDebug(1, "****** enter EventDumpTrigger ******");
1246   AliDebug(1, Form("Number of Reconstructed tracks : %d ",  nTriggerTracks));
1247
1248   // Loop over reconstructed tracks
1249   for (Int_t trackIndex = 0; trackIndex < nTriggerTracks; trackIndex++) {
1250     triggertrack = (AliMUONTriggerTrack*)fMUONData->RecTriggerTracks()->At(trackIndex);
1251       printf(" trigger track number %i x11=%f y11=%f thetax=%f thetay=%f \n",
1252              trackIndex,
1253              triggertrack->GetX11(),triggertrack->GetY11(),
1254              triggertrack->GetThetax(),triggertrack->GetThetay());
1255   }
1256 }
1257
1258 //__________________________________________________________________________
1259 void AliMUONTrackReconstructor::MakeTrackCandidatesK(void)
1260 {
1261   // To make initial tracks for Kalman filter from the list of segments
1262   Int_t istat, iseg;
1263   AliMUONSegment *segment;
1264   AliMUONTrackK *trackK;
1265
1266   AliDebug(1,"Enter MakeTrackCandidatesK");
1267
1268   AliMUONTrackK a(this, fHitsForRecPtr);
1269   // Loop over stations(1...) 5 and 4
1270   for (istat=4; istat>=3; istat--) {
1271     // Loop over segments in the station
1272     for (iseg=0; iseg<fNSegments[istat]; iseg++) {
1273       // Transform segments to tracks and evaluate covariance matrix
1274       segment = (AliMUONSegment*) ((*fSegmentsPtr[istat])[iseg]);
1275       trackK = new ((*fRecTracksPtr)[fNRecTracks++]) AliMUONTrackK(segment);
1276     } // for (iseg=0;...)
1277   } // for (istat=4;...)
1278   return;
1279 }
1280
1281 //__________________________________________________________________________
1282 void AliMUONTrackReconstructor::FollowTracksK(void)
1283 {
1284   // Follow tracks using Kalman filter
1285   Bool_t ok;
1286   Int_t icand, ichamBeg = 0, ichamEnd, chamBits;
1287   Double_t zDipole1, zDipole2;
1288   AliMUONTrackK *trackK;
1289   AliMUONHitForRec *hit;
1290   AliMUONRawCluster *clus;
1291   TClonesArray *rawclusters;
1292   clus = 0; rawclusters = 0;
1293
1294   zDipole1 = GetSimpleBPosition() + GetSimpleBLength()/2;
1295   zDipole2 = zDipole1 - GetSimpleBLength();
1296
1297   // Print hits
1298   trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[0]);
1299
1300   if (trackK->DebugLevel() > 0) {
1301     for (Int_t i1=0; i1<fNHitsForRec; i1++) {
1302       hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[i1]);
1303       printf(" Hit # %d %10.4f %10.4f %10.4f",
1304              hit->GetChamberNumber(), hit->GetBendingCoor(),
1305              hit->GetNonBendingCoor(), hit->GetZ());
1306
1307       // from raw clusters
1308       rawclusters = fMUONData->RawClusters(hit->GetChamberNumber());
1309       clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->
1310                                                            GetHitNumber());
1311       printf(" %d", clus->GetTrack(1));
1312       if (clus->GetTrack(2) != -1) printf(" %d \n", clus->GetTrack(2));
1313       else printf("\n");
1314
1315     }
1316   } // if (trackK->DebugLevel() > 0)
1317
1318   icand = -1;
1319   Int_t nSeeds;
1320   nSeeds = fNRecTracks; // starting number of seeds
1321   // Loop over track candidates
1322   while (icand < fNRecTracks-1) {
1323     icand ++;
1324     if (trackK->DebugLevel()>0) cout << " *** Kalman track candidate No. " << icand << endl;
1325     trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[icand]);
1326     if (trackK->GetRecover() < 0) continue; // failed track
1327
1328     // Discard candidate which will produce the double track
1329     /*
1330     if (icand > 0) {
1331       ok = CheckCandidateK(icand,nSeeds);
1332       if (!ok) {
1333         trackK->SetRecover(-1); // mark candidate to be removed
1334         continue;
1335       }
1336     }
1337     */
1338
1339     ok = kTRUE;
1340     if (trackK->GetRecover() == 0)
1341       hit = (AliMUONHitForRec*) trackK->GetTrackHits()->Last(); // last hit
1342     else
1343       hit = trackK->GetHitLastOk(); // hit where track stopped
1344
1345     if (hit) ichamBeg = hit->GetChamberNumber();
1346     ichamEnd = 0;
1347     // Check propagation direction
1348     if (!hit) { ichamBeg = ichamEnd; AliFatal(" ??? "); }
1349     else if (trackK->GetTrackDir() < 0) {
1350       ichamEnd = 9; // forward propagation
1351       ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2);
1352       if (ok) {
1353         ichamBeg = ichamEnd;
1354         ichamEnd = 6; // backward propagation
1355         // Change weight matrix and zero fChi2 for backpropagation
1356         trackK->StartBack();
1357         trackK->SetTrackDir(1);
1358         ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kTRUE,zDipole1,zDipole2);
1359         ichamBeg = ichamEnd;
1360         ichamEnd = 0;
1361       }
1362     } else {
1363       if (trackK->GetBPFlag()) {
1364         // backpropagation
1365         ichamEnd = 6; // backward propagation
1366         // Change weight matrix and zero fChi2 for backpropagation
1367         trackK->StartBack();
1368         ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kTRUE,zDipole1,zDipole2);
1369         ichamBeg = ichamEnd;
1370         ichamEnd = 0;
1371       }
1372     }
1373
1374     if (ok) {
1375       trackK->SetTrackDir(1);
1376       trackK->SetBPFlag(kFALSE);
1377       ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2);
1378     }
1379     if (!ok) { trackK->SetRecover(-1); continue; } // mark candidate to be removed
1380
1381     // Apply smoother
1382     if (trackK->GetRecover() >= 0) {
1383       ok = trackK->Smooth();
1384       if (!ok) trackK->SetRecover(-1); // mark candidate to be removed
1385     }
1386
1387     // Majority 3 of 4 in first 2 stations
1388     if (!ok) continue;
1389     chamBits = 0;
1390     Double_t chi2max = 0;
1391     for (Int_t i=0; i<trackK->GetNTrackHits(); i++) {
1392       hit = (AliMUONHitForRec*) (*trackK->GetTrackHits())[i];
1393       chamBits |= BIT(hit->GetChamberNumber());
1394       if (trackK->GetChi2PerPoint(i) > chi2max) chi2max = trackK->GetChi2PerPoint(i);
1395     }
1396     if (!((chamBits&3)==3 || (chamBits>>2&3)==3) && chi2max > 25) {
1397       //trackK->Recover();
1398       trackK->SetRecover(-1); //mark candidate to be removed
1399       continue;
1400     }
1401     if (ok) trackK->SetTrackQuality(0); // compute "track quality"
1402   } // while
1403
1404   for (Int_t i=0; i<fNRecTracks; i++) {
1405     trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[i]);
1406     if (trackK->GetRecover() < 0) fRecTracksPtr->RemoveAt(i);
1407   }
1408
1409   // Compress TClonesArray
1410   fRecTracksPtr->Compress();
1411   fNRecTracks = fRecTracksPtr->GetEntriesFast();
1412   return;
1413 }
1414
1415 //__________________________________________________________________________
1416 Bool_t AliMUONTrackReconstructor::CheckCandidateK(Int_t icand, Int_t nSeeds) const
1417 {
1418   // Discards track candidate if it will produce the double track (having
1419   // the same seed segment hits as hits of a good track found before)
1420   AliMUONTrackK *track1, *track2;
1421   AliMUONHitForRec *hit1, *hit2, *hit;
1422
1423   track1 = (AliMUONTrackK*) ((*fRecTracksPtr)[icand]);
1424   hit1 = (AliMUONHitForRec*) (*track1->GetTrackHits())[0]; // 1'st hit
1425   hit2 = (AliMUONHitForRec*) (*track1->GetTrackHits())[1]; // 2'nd hit
1426
1427   for (Int_t i=0; i<icand; i++) {
1428     track2 = (AliMUONTrackK*) ((*fRecTracksPtr)[i]);
1429     //if (track2->GetRecover() < 0) continue;
1430     if (track2->GetRecover() < 0 && icand >= nSeeds) continue;
1431
1432     if (track1->GetStartSegment() == track2->GetStartSegment()) {
1433       return kFALSE;
1434     } else {
1435       Int_t nSame = 0;
1436       for (Int_t j=0; j<track2->GetNTrackHits(); j++) {
1437         hit = (AliMUONHitForRec*) (*track2->GetTrackHits())[j];
1438         if (hit == hit1 || hit == hit2) {
1439           nSame++;
1440           if (nSame == 2) return kFALSE;
1441         }
1442       } // for (Int_t j=0;
1443     }
1444   } // for (Int_t i=0;
1445   return kTRUE;
1446 }
1447
1448 //__________________________________________________________________________
1449 void AliMUONTrackReconstructor::RemoveDoubleTracksK(void)
1450 {
1451   // Removes double tracks (sharing more than half of their hits). Keeps
1452   // the track with higher quality
1453   AliMUONTrackK *track1, *track2, *trackToKill;
1454
1455   // Sort tracks according to their quality
1456   fRecTracksPtr->Sort();
1457
1458   // Loop over first track of the pair
1459   track1 = (AliMUONTrackK*) fRecTracksPtr->First();
1460   Int_t debug = track1->DebugLevel();
1461   while (track1) {
1462     // Loop over second track of the pair
1463     track2 = (AliMUONTrackK*) fRecTracksPtr->After(track1);
1464     while (track2) {
1465       // Check whether or not to keep track2
1466       if (!track2->KeepTrack(track1)) {
1467         if (debug >= 0) cout << " Killed track: " << 1/(*track2->GetTrackParameters())(4,0) <<
1468           " " << track2->GetTrackQuality() << endl;
1469         trackToKill = track2;
1470         track2 = (AliMUONTrackK*) fRecTracksPtr->After(track2);
1471         trackToKill->Kill();
1472         fRecTracksPtr->Compress();
1473       } else track2 = (AliMUONTrackK*) fRecTracksPtr->After(track2);
1474     } // track2
1475     track1 = (AliMUONTrackK*) fRecTracksPtr->After(track1);
1476   } // track1
1477
1478   fNRecTracks = fRecTracksPtr->GetEntriesFast();
1479   if (debug >= 0) cout << " Number of Kalman tracks: " << fNRecTracks << endl;
1480 }
1481
1482 //__________________________________________________________________________
1483 void AliMUONTrackReconstructor::GoToVertex(void)
1484 {
1485   // Propagates track to the vertex thru absorber
1486   // (using Branson correction for now)
1487
1488   Double_t zVertex;
1489   zVertex = 0;
1490   for (Int_t i=0; i<fNRecTracks; i++) {
1491     //((AliMUONTrackK*)(*fRecTracksPtr)[i])->Branson();
1492     ((AliMUONTrackK*)(*fRecTracksPtr)[i])->SetTrackQuality(1); // compute Chi2
1493     //((AliMUONTrackK*)(*fRecTracksPtr)[i])->GoToZ(zVertex); // w/out absorber
1494     ((AliMUONTrackK*)(*fRecTracksPtr)[i])->GoToVertex(1); // with absorber
1495   }
1496 }
1497
1498 //__________________________________________________________________________
1499 void AliMUONTrackReconstructor::SetTrackMethod(Int_t iTrackMethod)
1500 {
1501   // Set track method and recreate track container if necessary
1502
1503   fTrackMethod = TMath::Min (iTrackMethod, 3);
1504   fTrackMethod = TMath::Max (fTrackMethod, 1);
1505   if (fTrackMethod != 1) {
1506     if (fRecTracksPtr) delete fRecTracksPtr;
1507     fRecTracksPtr = new TClonesArray("AliMUONTrackK", 10);
1508     if (fTrackMethod == 2) cout << " *** Tracking with the Kalman filter *** " << endl;
1509     else cout << " *** Combined cluster / track finder ***" << endl;
1510   } else cout << " *** Original tracking *** " << endl;
1511
1512 }