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