]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - MUON/AliMUONRecoCheck.cxx
Create the magnetic field map in the reconstruction macros (Yuri)
[u/mrichter/AliRoot.git] / MUON / AliMUONRecoCheck.cxx
index a6a3ca63885ace6e0bda424e35999c453f7f91d6..4cfa4b784c4f080bbdeadb8c52c73df84b81a3be 100644 (file)
 /**************************************************************************
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- *                                                                        *
- * Author: The ALICE Off-line Project.                                    *
- * Contributors are mentioned in the code where appropriate.              *
- *                                                                        *
- * Permission to use, copy, modify and distribute this software and its   *
- * documentation strictly for non-commercial purposes is hereby granted   *
- * without fee, provided that the above copyright notice appears in all   *
- * copies and that both the copyright notice and this permission notice   *
- * appear in the supporting documentation. The authors make no claims     *
- * about the suitability of this software for any purpose. It is          *
- * provided "as is" without express or implied warranty.                  *
- **************************************************************************/
+* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+*                                                                        *
+* Author: The ALICE Off-line Project.                                    *
+* Contributors are mentioned in the code where appropriate.              *
+*                                                                        *
+* Permission to use, copy, modify and distribute this software and its   *
+* documentation strictly for non-commercial purposes is hereby granted   *
+* without fee, provided that the above copyright notice appears in all   *
+* copies and that both the copyright notice and this permission notice   *
+* appear in the supporting documentation. The authors make no claims     *
+* about the suitability of this software for any purpose. It is          *
+* provided "as is" without express or implied warranty.                  *
+**************************************************************************/
 
 /* $Id$ */
 
+//-----------------------------------------------------------------------------
 /// \class AliMUONRecoCheck
-/// \brief Utility class to check reconstruction
+/// Utility class to check reconstruction
 /// Reconstructed tracks are compared to reference tracks. 
 /// The reference tracks are built from AliTrackReference for the
 /// hit in chamber (0..9) and from kinematics for the vertex parameters.     
+//-----------------------------------------------------------------------------
+
+// 13 Nov 2007:
+// Added a method to create a list of reconstructed AliMUONTrack objects from
+// ESD data. This is necessary since the track objects that are actually created
+// during offline reconstruction are no longer stored to disk.
+//  - Artur Szostak <artursz@iafrica.com>
+// 25 Jan 2008:
+// Use the new ESDInterface to create MUON objects from ESD data
+// - Philippe Pillot
 
-#include "AliMUON.h"
 #include "AliMUONRecoCheck.h"
 #include "AliMUONTrack.h"
-#include "AliMUONData.h"
+#include "AliMUONVTrackStore.h"
+#include "AliMUONVCluster.h"
+#include "AliMUONVClusterStore.h"
 #include "AliMUONConstants.h"
+#include "AliMUONESDInterface.h"
 
-#include "AliRun.h" // for gAlice
-#include "AliLoader.h" 
-#include "AliRunLoader.h" 
-#include "AliHeader.h"
+#include "AliMCEventHandler.h"
+#include "AliMCEvent.h"
 #include "AliStack.h"
 #include "AliTrackReference.h"
-#include "AliLog.h" 
+#include "AliLog.h"
+#include "AliESDEvent.h"
+#include "AliESDMuonTrack.h"
 
+#include <TFile.h>
+#include <TTree.h>
 #include <TParticle.h>
+#include <TParticlePDG.h>
+#include <Riostream.h>
 
 /// \cond CLASSIMP
 ClassImp(AliMUONRecoCheck)
 /// \endcond
 
 //_____________________________________________________________________________
-  AliMUONRecoCheck::AliMUONRecoCheck(AliRunLoader *runloader, AliMUONData *muondata)
-  : TObject(),
-  fRunLoader(0x0),
-  fMUONData(0x0),
-  fMuonTrackRef(0x0),
-  fTrackReco(0x0),
-  fReconstructibleTracks(0),
-  fRecoTracks(0)
+AliMUONRecoCheck::AliMUONRecoCheck(Char_t *esdFileName, Char_t *pathSim)
+: TObject(),
+fMCEventHandler(new AliMCEventHandler()),
+fESDEvent(new AliESDEvent()),
+fESDTree (0x0),
+fESDFile (0x0),
+fCurrentEvent(0),
+fTrackRefStore(0x0),
+fRecoTrackRefStore(0x0),
+fRecoTrackStore(0x0)
 {
-/// Constructor
-
-  fMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
-
-  // run loader
-  fRunLoader = runloader;
-  if (!fRunLoader) {
-    AliError(Form("no run loader found " ));
+  /// Normal ctor
+  
+  // TrackRefs and Particules
+  fMCEventHandler->SetInputPath(pathSim);
+  fMCEventHandler->InitIO("");
+  
+  // ESD MUON Tracks
+  fESDFile = TFile::Open(esdFileName); // open the file
+  if (!fESDFile || !fESDFile->IsOpen()) {
+    AliError(Form("opening ESD file %s failed", esdFileName));
+    fESDFile = 0x0;
     return;
   }
-
-  // container
-  fMUONData  = muondata;
-  if (!fMUONData) {
-    AliError(Form("no MUONData found " ));
+  fESDTree = (TTree*) fESDFile->Get("esdTree"); // get the tree
+  if (!fESDTree) {
+    AliError("no ESD tree found");
+    fESDFile->Close();
+    fESDFile = 0x0;
     return;
   }
-
+  fESDEvent->ReadFromTree(fESDTree); // link fESDEvent to the tree
 }
 
-
 //_____________________________________________________________________________
 AliMUONRecoCheck::~AliMUONRecoCheck()
 {
-/// Destructor
+  /// Destructor
+  delete fMCEventHandler;
+  delete fESDEvent;
+  if (fESDFile) fESDFile->Close();
+  ResetStores();
+}
 
-  fMuonTrackRef->Delete();
-  delete fMuonTrackRef;
+//_____________________________________________________________________________
+void AliMUONRecoCheck::ResetStores()
+{
+  /// Deletes all the store objects that have been created and resets the pointers to 0x0
+  delete fTrackRefStore;      fTrackRefStore = 0x0;
+  delete fRecoTrackRefStore;  fRecoTrackRefStore = 0x0;
+  delete fRecoTrackStore;     fRecoTrackStore = 0x0;
 }
 
 //_____________________________________________________________________________
-void AliMUONRecoCheck::MakeTrackRef()
+Int_t AliMUONRecoCheck::NumberOfEvents() const
 {
-/// Make reconstructible tracks
+  /// Return the number of events
+  if (fESDTree) return fESDTree->GetEntries();
+  return 0;
+}
 
-  AliTrackReference *trackReference;
-  AliMUONTrack *muonTrack;
-  AliMUONTrackParam *trackParam;
-  AliMUONHitForRec *hitForRec;
-  Float_t x, y, z, pX, pY, pZ, pYZ;
-  Int_t track, trackSave;
-  TParticle *particle;   
-  Float_t bendingSlope = 0;
-  Float_t  nonBendingSlope = 0;
-  Float_t inverseBendingMomentum = 0;
+//_____________________________________________________________________________
+AliMUONVTrackStore* AliMUONRecoCheck::ReconstructedTracks(Int_t event)
+{
+  /// Return a track store containing the reconstructed tracks (converted into 
+  /// MUONTrack objects) for a given event
+  if (event != fCurrentEvent) {
+    ResetStores();
+    fCurrentEvent = event;
+  }
   
-  TTree* treeTR = fRunLoader->TreeTR();
-  if (treeTR == NULL) return;
-
-  TBranch* branch = treeTR->GetBranch("MUON");
-  if (branch == NULL) return;
-
-  TClonesArray* trackRefs = 0;
-  branch->SetAddress(&trackRefs);
-  branch->SetAutoDelete(kTRUE);  
+  if (fRecoTrackStore != 0x0) return fRecoTrackStore;
+  else {
+    if (!fESDTree) return 0x0;
+    if (fESDTree->GetEvent(event) <= 0) {
+      AliError(Form("fails to read ESD object for event %d", event));
+      return 0x0;
+    }
+    MakeReconstructedTracks();
+    return fRecoTrackStore;
+  }
+}
 
-  Int_t nTrackRef = (Int_t)branch->GetEntries();
-  track = trackSave = -999;
-  Bool_t isNewTrack;
-  Int_t iHitMin, iChamber, detElemId;
+//_____________________________________________________________________________
+AliMUONVTrackStore* AliMUONRecoCheck::TrackRefs(Int_t event)
+{
+  /// Return a track store containing the track references (converted into 
+  /// MUONTrack objects) for a given event
+  if (event != fCurrentEvent) {
+    ResetStores();
+    fCurrentEvent = event;
+  }
+  
+  if (fTrackRefStore != 0x0) return fTrackRefStore;
+  else {
+    if (!fMCEventHandler->GetEvent(event)) {
+      AliError(Form("fails to read MC objects for event %d", event));
+      return 0x0;
+    }
+    MakeTrackRefs();
+    return fTrackRefStore;
+  }
+}
 
-  trackParam = new AliMUONTrackParam();
-  hitForRec = new AliMUONHitForRec();
-  muonTrack = new AliMUONTrack();
+//_____________________________________________________________________________
+AliMUONVTrackStore* AliMUONRecoCheck::ReconstructibleTracks(Int_t event)
+{
+  /// Return a track store containing the reconstructible tracks for a given event
+  if (event != fCurrentEvent) {
+    ResetStores();
+    fCurrentEvent = event;
+  }
+  
+  if (fRecoTrackRefStore != 0x0) return fRecoTrackRefStore;
+  else {
+    if (TrackRefs(event) == 0x0) return 0x0;
+    MakeReconstructibleTracks();
+    return fRecoTrackRefStore;
+  }
+}
 
-  Int_t max = fRunLoader->GetHeader()->Stack()->GetNtrack();
-  for (Int_t iTrackRef  = 0; iTrackRef < nTrackRef; iTrackRef++) {
+//_____________________________________________________________________________
+void AliMUONRecoCheck::MakeReconstructedTracks()
+{
+  /// Make reconstructed tracks
+  if (!(fRecoTrackStore = AliMUONESDInterface::NewTrackStore())) return;
+  
+  // loop over all reconstructed tracks and add them to the store
+  Int_t nTracks = (Int_t) fESDEvent->GetNumberOfMuonTracks();
+  for (Int_t iTrack = 0; iTrack < nTracks; iTrack++)
+    AliMUONESDInterface::Add(*(fESDEvent->GetMuonTrack(iTrack)), *fRecoTrackStore);
+  
+}
 
-    branch->GetEntry(iTrackRef);
+//_____________________________________________________________________________
+void AliMUONRecoCheck::MakeTrackRefs()
+{
+  /// Make reconstructible tracks
+  AliMUONVTrackStore *tmpTrackRefStore = AliMUONESDInterface::NewTrackStore();
+  if (!tmpTrackRefStore) return;
+  
+  Double_t x, y, z, pX, pY, pZ, bendingSlope, nonBendingSlope, inverseBendingMomentum;
+  TParticle* particle;
+  TClonesArray* trackRefs;
+  Int_t nTrackRef = fMCEventHandler->MCEvent()->GetNumberOfTracks();
+  AliMUONVClusterStore* cStore = AliMUONESDInterface::NewClusterStore();
+  if (!cStore) return;
+  AliMUONVCluster* hit = cStore->CreateCluster(0,0,0);
+  
+  // loop over simulated tracks
+  for (Int_t iTrackRef  = 0; iTrackRef < nTrackRef; ++iTrackRef) {
+    Int_t nHits = fMCEventHandler->GetParticleAndTR(iTrackRef, particle, trackRefs);
     
-    iHitMin = 0;
-    isNewTrack = kTRUE;
+    // skip empty trackRefs
+    if (nHits < 1) continue;
     
-    if (!trackRefs->GetEntries()) continue; 
-
-    while (isNewTrack) {
-      for (Int_t iHit = iHitMin; iHit < trackRefs->GetEntries(); iHit++) {
+    // get the particle charge for further calculation
+    TParticlePDG* ppdg = particle->GetPDG();
+    Int_t charge = (Int_t)(ppdg->Charge()/3.0);
+    
+    AliMUONTrack track;
+    
+    // loop over simulated track hits
+    for (Int_t iHit = 0; iHit < nHits; ++iHit) {        
+      AliTrackReference* trackReference = static_cast<AliTrackReference*>(trackRefs->UncheckedAt(iHit));
       
-       trackReference = (AliTrackReference*)trackRefs->At(iHit);
-       x = trackReference->X();
-       y = trackReference->Y();
-       z = trackReference->Z();
-       pX = trackReference->Px();
-       pY = trackReference->Py();
-       pZ = trackReference->Pz();
-
-       track = trackReference->GetTrack();
-
-       if(track >= max){
-         AliWarningStream()
-           << "Track ID " << track 
-           << " larger than max number of particles " << max << endl;
-         isNewTrack = kFALSE;
-         break;
-       }
-       if (track != trackSave && iHit != 0) {
-         iHitMin = iHit;
-         trackSave = track;
-         break;
-       }
-
-       // track parameters at hit
-       trackParam->SetBendingCoor(y);
-       trackParam->SetNonBendingCoor(x);
-       trackParam->SetZ(z);
-
-       if (TMath::Abs(pZ) > 0) {
-         bendingSlope = pY/pZ;
-         nonBendingSlope = pX/pZ;
-       }
-       pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
-       if (pYZ >0) inverseBendingMomentum = 1/pYZ; 
-
-       trackParam->SetBendingSlope(bendingSlope);
-       trackParam->SetNonBendingSlope(nonBendingSlope);
-       trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
-
-       hitForRec->SetBendingCoor(y);
-       hitForRec->SetNonBendingCoor(x);
-       hitForRec->SetZ(z);
-       hitForRec->SetBendingReso2(0.0); 
-       hitForRec->SetNonBendingReso2(0.0);
-       detElemId = hitForRec->GetDetElemId();
-       if (detElemId) iChamber = detElemId / 100 - 1; 
-       else iChamber = AliMUONConstants::ChamberNumber(z);
-       hitForRec->SetChamberNumber(iChamber);
-
-       muonTrack->AddTrackParamAtHit(trackParam);
-       muonTrack->AddHitForRecAtHit(hitForRec);
-       muonTrack->SetTrackID(track);
-
-       trackSave = track;
-       if (iHit == trackRefs->GetEntries()-1) isNewTrack = kFALSE;
-      }
-
-      // track parameters at vertex 
-      particle = fRunLoader->GetHeader()->Stack()->Particle(muonTrack->GetTrackID());
-
-      if (particle) {
-
-       x = particle->Vx();
-       y = particle->Vy();
-       z = particle->Vz();
-       pX = particle->Px();
-       pY = particle->Py();
-       pZ = particle->Pz();
-
-       trackParam->SetBendingCoor(y);
-       trackParam->SetNonBendingCoor(x);
-       trackParam->SetZ(z);
-       if (TMath::Abs(pZ) > 0) {
-         bendingSlope = pY/pZ;
-         nonBendingSlope = pX/pZ;
-       }
-       pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
-       if (pYZ >0) inverseBendingMomentum = 1/pYZ;       
-       trackParam->SetBendingSlope(bendingSlope);
-       trackParam->SetNonBendingSlope(nonBendingSlope);
-       trackParam->SetInverseBendingMomentum(inverseBendingMomentum);
+      // skip trackRefs not in MUON
+      if (trackReference->DetectorId() != AliTrackReference::kMUON) continue;
+    
+      // Get track parameters of current hit
+      x = trackReference->X();
+      y = trackReference->Y();
+      z = trackReference->Z();
+      pX = trackReference->Px();
+      pY = trackReference->Py();
+      pZ = trackReference->Pz();
+      
+      // check chamberId of current trackReference
+      Int_t chamberId = AliMUONConstants::ChamberNumber(z);
+      if (chamberId < 0 || chamberId >= AliMUONConstants::NTrackingCh()) continue;
       
-       muonTrack->SetTrackParamAtVertex(trackParam);
+      // set hit parameters
+      hit->SetUniqueID(AliMUONVCluster::BuildUniqueID(chamberId, 0, 0));
+      hit->SetXYZ(x,y,z);
+      hit->SetErrXY(0.,0.);
+      
+      // compute track parameters at hit
+      Double_t bendingSlope = 0;
+      Double_t nonBendingSlope = 0;
+      Double_t inverseBendingMomentum = 0;
+      if (TMath::Abs(pZ) > 0) {
+       bendingSlope = pY/pZ;
+       nonBendingSlope = pX/pZ;
       }
-
-      AddMuonTrackReference(muonTrack);
-      muonTrack->ResetTrackParamAtHit();
-      muonTrack->ResetHitForRecAtHit();
+      Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
+      if (pYZ >0) inverseBendingMomentum = 1/pYZ; 
+      inverseBendingMomentum *= charge;
+      
+      // set track parameters at hit
+      AliMUONTrackParam trackParam;
+      trackParam.SetNonBendingCoor(x);
+      trackParam.SetBendingCoor(y);
+      trackParam.SetZ(z);
+      trackParam.SetBendingSlope(bendingSlope);
+      trackParam.SetNonBendingSlope(nonBendingSlope);
+      trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
       
-    } // end while isNewTrack
+      // add track parameters at current hit to the track
+      track.AddTrackParamAtCluster(trackParam, *hit, kTRUE);
+    }
+    
+    // if none of the track hits was in MUON, goto the next track
+    if (track.GetNClusters() < 1) continue;
+    
+    // get track parameters at particle's vertex
+    x = particle->Vx();
+    y = particle->Vy();
+    z = particle->Vz();
+    pX = particle->Px();
+    pY = particle->Py();
+    pZ = particle->Pz();
+    
+    // compute rest of track parameters at particle's vertex
+    bendingSlope = 0;
+    nonBendingSlope = 0;
+    inverseBendingMomentum = 0;
+    if (TMath::Abs(pZ) > 0) {
+      bendingSlope = pY/pZ;
+      nonBendingSlope = pX/pZ;
+    }
+    Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
+    if (pYZ >0) inverseBendingMomentum = 1/pYZ;
+    inverseBendingMomentum *= charge;
+    
+    // set track parameters at particle's vertex
+    AliMUONTrackParam trackParamAtVertex;
+    trackParamAtVertex.SetNonBendingCoor(x);
+    trackParamAtVertex.SetBendingCoor(y);
+    trackParamAtVertex.SetZ(z);
+    trackParamAtVertex.SetBendingSlope(bendingSlope);
+    trackParamAtVertex.SetNonBendingSlope(nonBendingSlope);
+    trackParamAtVertex.SetInverseBendingMomentum(inverseBendingMomentum);
+    
+    // add track parameters at vertex
+    track.SetTrackParamAtVertex(&trackParamAtVertex);
+    
+    // store the track
+    track.SetTrackID(iTrackRef);
+    tmpTrackRefStore->Add(track);
   }
   
-  CleanMuonTrackRef();
+  CleanMuonTrackRef(tmpTrackRefStore);
   
-  ReconstructibleTracks();
-
-  delete muonTrack;
-  delete trackParam;
-  delete hitForRec;
-
+  delete hit;
+  delete cStore;
+  delete tmpTrackRefStore;
 }
 
-//____________________________________________________________________________
-TClonesArray* AliMUONRecoCheck::GetTrackReco()
-{
-/// Return TClonesArray of reconstructed tracks
-
-  fMUONData->ResetRecTracks();
-  fMUONData->SetTreeAddress("RT");
-  fTrackReco = fMUONData->RecTracks(); 
-  fMUONData->GetRecTracks();
-  fRecoTracks = fTrackReco->GetEntriesFast();
-  return fTrackReco;
-}
 //_____________________________________________________________________________
-void AliMUONRecoCheck::PrintEvent() const
+void AliMUONRecoCheck::CleanMuonTrackRef(const AliMUONVTrackStore *tmpTrackRefStore)
 {
-/// Debug facility
-
-  AliMUONTrack *track;
-  AliMUONHitForRec *hitForRec;
-  TClonesArray *  hitForRecAtHit = 0;
-  Float_t xRec,yRec,zRec;
-
-  Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
+  /// Re-calculate hits parameters because two AliTrackReferences are recorded for
+  /// each chamber (one when particle is entering + one when particle is leaving 
+  /// the sensitive volume) 
+  if (!(fTrackRefStore = AliMUONESDInterface::NewTrackStore())) return;
   
-  printf(" ******************************* \n");
-  printf(" nb of tracks %d \n",nTrackRef);
-
-  for (Int_t index = 0; index < nTrackRef; index++) {
-    track = (AliMUONTrack*)fMuonTrackRef->At(index);
-    hitForRecAtHit = track->GetHitForRecAtHit();
-    Int_t nTrackHits = hitForRecAtHit->GetEntriesFast();
-    printf(" track number %d \n",index);
-    for (Int_t iHit = 0; iHit < nTrackHits; iHit++){
-      hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit); 
-      xRec  = hitForRec->GetNonBendingCoor();
-      yRec  = hitForRec->GetBendingCoor();
-      zRec  = hitForRec->GetZ();
-      printf("  x,y,z: %f , %f , %f \n",xRec,yRec,zRec);
-    }
-  }
-}
-
-//_____________________________________________________________________________
-void AliMUONRecoCheck::ResetTracks() const
-{
-/// Reset tracks
-
-  if (fMuonTrackRef) fMuonTrackRef->Clear();
-}
-//_____________________________________________________________________________
-void AliMUONRecoCheck::CleanMuonTrackRef()
-{
-/// Re-calculate hits parameters because two AliTrackReferences are recorded for
-/// each chamber (one when particle is entering + one when particle is leaving 
-/// the sensitive volume) 
-
-  Float_t maxGasGap = 1.; // cm 
-  AliMUONTrack *track, *trackNew;
-  AliMUONHitForRec *hitForRec, *hitForRec1, *hitForRec2;
-  AliMUONTrackParam *trackParam, *trackParam1, *trackParam2, *trackParamAtVertex;
-  TClonesArray *  hitForRecAtHit = 0;
-  TClonesArray *  trackParamAtHit = 0;
-  Float_t xRec,yRec,zRec;
-  Float_t xRec1,yRec1,zRec1;
-  Float_t xRec2,yRec2,zRec2;
-  Float_t bendingSlope,nonBendingSlope,bendingMomentum;
-  Float_t bendingSlope1,nonBendingSlope1,bendingMomentum1;
-  Float_t bendingSlope2,nonBendingSlope2,bendingMomentum2;
-  TClonesArray *newMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
-  Int_t iHit1;
-  Int_t iChamber = 0, detElemId = 0;
-  Int_t nRec = 0;
-  Int_t nTrackHits = 0;
-
-  hitForRec = new AliMUONHitForRec();
-  trackParam = new AliMUONTrackParam();
-  trackNew = new AliMUONTrack();
-
-  Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
-  for (Int_t index = 0; index < nTrackRef; index++) {
-    track = (AliMUONTrack*)fMuonTrackRef->At(index);
-    hitForRecAtHit = track->GetHitForRecAtHit();
-    trackParamAtHit = track->GetTrackParamAtHit();
-    trackParamAtVertex = track->GetTrackParamAtVertex();
-    nTrackHits = hitForRecAtHit->GetEntriesFast();
-    iHit1 = 0;
+  Double_t maxGasGap = 1.; // cm 
+  Double_t x, y, z, pX, pY, pZ, x1, y1, z1, pX1, pY1, pZ1, z2;
+  Double_t bendingSlope,nonBendingSlope,inverseBendingMomentum;
+  AliMUONVClusterStore* cStore = AliMUONESDInterface::NewClusterStore();
+  if (!cStore) return;
+  AliMUONVCluster* hit = cStore->CreateCluster(0,0,0);
+  
+  // create iterator
+  TIter next(tmpTrackRefStore->CreateIterator());
+  
+  // loop over tmpTrackRef
+  AliMUONTrack* track;
+  while ( ( track = static_cast<AliMUONTrack*>(next()) ) ) {
+    
+    AliMUONTrack newTrack;
+    
+    // loop over tmpTrackRef's hits
+    Int_t iHit1 = 0;
+    Int_t nTrackHits = track->GetNClusters();
     while (iHit1 < nTrackHits) {
-      hitForRec1 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit1); 
-      trackParam1 = (AliMUONTrackParam*) trackParamAtHit->At(iHit1); 
-      xRec1  = hitForRec1->GetNonBendingCoor();
-      yRec1  = hitForRec1->GetBendingCoor();
-      zRec1  = hitForRec1->GetZ();     
-      xRec   = xRec1;
-      yRec   = yRec1;
-      zRec   = zRec1;
-      bendingSlope1 = trackParam1->GetBendingSlope();
-      nonBendingSlope1 = trackParam1->GetNonBendingSlope();
-      bendingMomentum1 = 0;
-      if (TMath::Abs(trackParam1->GetInverseBendingMomentum()) > 0)
-       bendingMomentum1 = 1./trackParam1->GetInverseBendingMomentum();
-      bendingSlope = bendingSlope1;
-      nonBendingSlope = nonBendingSlope1;
-      bendingMomentum = bendingMomentum1;
-      nRec = 1;  
+      AliMUONTrackParam *trackParam1 = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit1);
+      
+      // get track parameters at hit1
+      x1  = trackParam1->GetNonBendingCoor();
+      y1  = trackParam1->GetBendingCoor();
+      z1  = trackParam1->GetZ();
+      pX1 = trackParam1->Px();
+      pY1 = trackParam1->Py();
+      pZ1 = trackParam1->Pz();
+      
+      // prepare new track parameters
+      x  = x1;
+      y  = y1;
+      z  = z1;
+      pX = pX1;
+      pY = pY1;
+      pZ = pZ1;
+      
+      // loop over next tmpTrackRef's hits
+      Int_t nCombinedHits = 1;
       for (Int_t iHit2 = iHit1+1; iHit2 < nTrackHits; iHit2++) {
-       hitForRec2 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit2); 
-       trackParam2 = (AliMUONTrackParam*) trackParamAtHit->At(iHit2); 
-       xRec2  = hitForRec2->GetNonBendingCoor();
-       yRec2  = hitForRec2->GetBendingCoor();
-       zRec2  = hitForRec2->GetZ();      
-       bendingSlope2 = trackParam2->GetBendingSlope();
-       nonBendingSlope2 = trackParam2->GetNonBendingSlope();
-       bendingMomentum2 = 0;
-       if (TMath::Abs(trackParam2->GetInverseBendingMomentum()) > 0)
-         bendingMomentum2 = 1./trackParam2->GetInverseBendingMomentum();
-       
-       if ( TMath::Abs(zRec2-zRec1) < maxGasGap ) {
-
-         nRec++;
-         xRec += xRec2;
-         yRec += yRec2;
-         zRec += zRec2;
-         bendingSlope += bendingSlope2;
-         nonBendingSlope += nonBendingSlope2;
-         bendingMomentum += bendingMomentum2;
-         iHit1 = iHit2;
-       }
+       AliMUONTrackParam *trackParam2 = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit2);
+        
+       // get z position of hit2
+       z2 = trackParam2->GetZ();
        
-      } // end iHit2
-      xRec /= (Float_t)nRec;
-      yRec /= (Float_t)nRec;
-      zRec /= (Float_t)nRec;
-      bendingSlope /= (Float_t)nRec;
-      nonBendingSlope /= (Float_t)nRec;
-      bendingMomentum /= (Float_t)nRec;
+       // complete new track parameters if hit2 is on the same detection element
+        if ( TMath::Abs(z2-z1) < maxGasGap ) {
+         x  += trackParam2->GetNonBendingCoor();
+         y  += trackParam2->GetBendingCoor();
+         z  += z2;
+         pX += trackParam2->Px();
+         pY += trackParam2->Py();
+         pZ += trackParam2->Pz();
+         nCombinedHits++;
+          iHit1 = iHit2;
+        }
+        
+      }
       
-      hitForRec->SetNonBendingCoor(xRec);
-      hitForRec->SetBendingCoor(yRec);
-      hitForRec->SetZ(zRec);
-      detElemId = hitForRec->GetDetElemId();
-      if (detElemId) iChamber = detElemId / 100 - 1;
-      else iChamber = AliMUONConstants::ChamberNumber(zRec);
-      hitForRec->SetChamberNumber(iChamber);
-      hitForRec->SetBendingReso2(0.0); 
-      hitForRec->SetNonBendingReso2(0.0); 
-      trackParam->SetNonBendingCoor(xRec);
-      trackParam->SetBendingCoor(yRec);
-      trackParam->SetZ(zRec);
-      trackParam->SetNonBendingSlope(nonBendingSlope);
-      trackParam->SetBendingSlope(bendingSlope);
-      if (TMath::Abs(bendingMomentum) > 0)
-       trackParam->SetInverseBendingMomentum(1./bendingMomentum);
-
-      trackNew->AddHitForRecAtHit(hitForRec);
-      trackNew->AddTrackParamAtHit(trackParam);
+      // finalize new track parameters
+      x  /= (Double_t)nCombinedHits;
+      y  /= (Double_t)nCombinedHits;
+      z  /= (Double_t)nCombinedHits;
+      pX /= (Double_t)nCombinedHits;
+      pY /= (Double_t)nCombinedHits;
+      pZ /= (Double_t)nCombinedHits;
+      bendingSlope = 0;
+      nonBendingSlope = 0;
+      inverseBendingMomentum = 0;
+      if (TMath::Abs(pZ) > 0) {
+       bendingSlope = pY/pZ;
+       nonBendingSlope = pX/pZ;
+      }
+      Double_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
+      if (pYZ >0) inverseBendingMomentum = 1/pYZ; 
+      inverseBendingMomentum *= trackParam1->GetCharge();
+      
+      // set hit parameters
+      hit->SetUniqueID(AliMUONVCluster::BuildUniqueID(trackParam1->GetClusterPtr()->GetChamberId(), 0, 0));
+      hit->SetXYZ(x,y,z);
+      hit->SetErrXY(0.,0.);
+      
+      // set new track parameters at new hit
+      AliMUONTrackParam trackParam;
+      trackParam.SetNonBendingCoor(x);
+      trackParam.SetBendingCoor(y);
+      trackParam.SetZ(z);
+      trackParam.SetBendingSlope(bendingSlope);
+      trackParam.SetNonBendingSlope(nonBendingSlope);
+      trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
+      
+      // add track parameters at current hit to the track
+      newTrack.AddTrackParamAtCluster(trackParam, *hit, kTRUE);
       
       iHit1++;
-    } // end iHit1
-
-    trackNew->SetTrackID(track->GetTrackID());
-    trackNew->SetTrackParamAtVertex(trackParamAtVertex);
-    {new ((*newMuonTrackRef)[newMuonTrackRef->GetEntriesFast()]) AliMUONTrack(*trackNew);}    
-    trackNew->ResetHitForRecAtHit();
-    trackNew->ResetTrackParamAtHit();
+    }
+    
+    newTrack.SetTrackID(track->GetTrackID());
+    newTrack.SetTrackParamAtVertex(track->GetTrackParamAtVertex());
+    fTrackRefStore->Add(newTrack);
     
-  } // end trackRef
-
-  fMuonTrackRef->Clear();
-  nTrackRef = newMuonTrackRef->GetEntriesFast();
-  for (Int_t index = 0; index < nTrackRef; index++) {
-    track = (AliMUONTrack*)newMuonTrackRef->At(index);
-    AddMuonTrackReference(track);
   }
   
-
-  delete trackNew;
-  delete hitForRec;
-  delete trackParam;
-  newMuonTrackRef->Delete();
-  delete newMuonTrackRef;
-  
+  delete hit;
+  delete cStore;
 }
 
 //_____________________________________________________________________________
-void AliMUONRecoCheck::ReconstructibleTracks()
+void AliMUONRecoCheck::MakeReconstructibleTracks()
 {
-/// Calculate the number of reconstructible tracks
+  /// Isolate the reconstructible tracks
+  if (!(fRecoTrackRefStore = AliMUONESDInterface::NewTrackStore())) return;
   
+  // create iterator on trackRef
+  TIter next(fTrackRefStore->CreateIterator());
+  
+  // loop over trackRef
   AliMUONTrack* track;
-  TClonesArray* hitForRecAtHit = NULL;
-  AliMUONHitForRec* hitForRec;
-  Float_t zRec;
-  Int_t nTrackRef, nTrackHits;
-  Int_t isChamberInTrack[10];
-  Int_t iChamber = 0;
-  Bool_t isTrackOK = kTRUE;
-
-  fReconstructibleTracks = 0;
-
-  nTrackRef = fMuonTrackRef->GetEntriesFast();
-  for (Int_t index = 0; index < nTrackRef; index++) {
-    track = (AliMUONTrack*)fMuonTrackRef->At(index);
-    hitForRecAtHit = track->GetHitForRecAtHit();
-    nTrackHits = hitForRecAtHit->GetEntriesFast();
-    for (Int_t ch = 0; ch < 10; ch++) isChamberInTrack[ch] = 0;
-    for ( Int_t iHit = 0; iHit < nTrackHits; iHit++) {
-      hitForRec = (AliMUONHitForRec*) hitForRecAtHit->At(iHit); 
-      zRec  = hitForRec->GetZ();
-      iChamber = hitForRec->GetChamberNumber();
-      if (iChamber < 0 || iChamber > 10) continue;
-      isChamberInTrack[iChamber] = 1;
+  while ( ( track = static_cast<AliMUONTrack*>(next()) ) ) {
+    
+    Bool_t* chamberInTrack = new Bool_t(AliMUONConstants::NTrackingCh());
+    for (Int_t iCh = 0; iCh < AliMUONConstants::NTrackingCh(); iCh++) chamberInTrack[iCh] = kFALSE;
+    
+    // loop over trackRef's hits to get hit chambers
+    Int_t nTrackHits = track->GetNClusters();
+    for (Int_t iHit = 0; iHit < nTrackHits; iHit++) {
+      AliMUONVCluster* hit = ((AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit))->GetClusterPtr(); 
+      chamberInTrack[hit->GetChamberId()] = kTRUE;
     } 
+    
     // track is reconstructible if the particle is depositing a hit
     // in the following chamber combinations:
-
-    isTrackOK = kTRUE;
-    if (!isChamberInTrack[0] && !isChamberInTrack[1]) isTrackOK = kFALSE;
-    if (!isChamberInTrack[2] && !isChamberInTrack[3]) isTrackOK = kFALSE;
-    if (!isChamberInTrack[4] && !isChamberInTrack[5]) isTrackOK = kFALSE;
-    Int_t nHitsInLastStations=0;
-    for (Int_t ch = 6; ch < AliMUONConstants::NTrackingCh(); ch++)
-      if (isChamberInTrack[ch]) nHitsInLastStations++; 
-    if(nHitsInLastStations < 3) isTrackOK = kFALSE;
-
-    if (isTrackOK) fReconstructibleTracks++;
-    if (!isTrackOK) fMuonTrackRef->Remove(track); // remove non reconstructible tracks
+    Bool_t trackOK = kTRUE;
+    if (!chamberInTrack[0] && !chamberInTrack[1]) trackOK = kFALSE;
+    if (!chamberInTrack[2] && !chamberInTrack[3]) trackOK = kFALSE;
+    if (!chamberInTrack[4] && !chamberInTrack[5]) trackOK = kFALSE;
+    Int_t nHitsInLastStations = 0;
+    for (Int_t iCh = 6; iCh < AliMUONConstants::NTrackingCh(); iCh++)
+      if (chamberInTrack[iCh]) nHitsInLastStations++; 
+    if(nHitsInLastStations < 3) trackOK = kFALSE;
+    
+    // Add reconstructible tracks to fRecoTrackRefStore
+    if (trackOK) fRecoTrackRefStore->Add(*track);
+    
+    delete [] chamberInTrack;
   }
-  fMuonTrackRef->Compress();
+
 }