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