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