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