]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - MUON/AliMUONRecoCheck.cxx
Updated alignment macros and related changes (Raffaele)
[u/mrichter/AliRoot.git] / MUON / AliMUONRecoCheck.cxx
index 207a9c9fffd93726c13e97561e492c6b15bdfe48..0d7d8ab40dacbb9bc506bc4ed81b435c4a222be2 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
-// -----------------------
-// Utility class to check the muon 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.     
+//-----------------------------------------------------------------------------
+/// \class AliMUONRecoCheck
+/// 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.     
+//-----------------------------------------------------------------------------
 
-#include <TParticle.h>
-
-#include "AliRun.h" // for gAlice
-#include "AliLog.h" 
-#include "AliLoader.h" 
-#include "AliRunLoader.h" 
-#include "AliTrackReference.h"
-#include "AliHeader.h"
-#include "AliMC.h"
-#include "AliStack.h"
-#include "AliMUON.h"
 #include "AliMUONRecoCheck.h"
+#include "AliMUONHitForRec.h"
 #include "AliMUONTrack.h"
-#include "AliMUONData.h"
 #include "AliMUONConstants.h"
+#include "AliMUONMCDataInterface.h"
+#include "AliMUONDataInterface.h"
+#include "AliStack.h"
+#include "AliTrackReference.h"
+#include "AliLog.h" 
+#include "AliMUONTrackStoreV1.h"
+#include <Riostream.h>
+#include <TParticle.h>
+#include <TParticlePDG.h>
 
+/// \cond CLASSIMP
 ClassImp(AliMUONRecoCheck)
+/// \endcond
 
 //_____________________________________________________________________________
-AliMUONRecoCheck::AliMUONRecoCheck(Char_t *chLoader)
+AliMUONRecoCheck::AliMUONRecoCheck(Char_t *chLoader, Char_t *chLoaderSim)
+: TObject(),
+fMCDataInterface(new AliMUONMCDataInterface(chLoaderSim)),
+fDataInterface(new AliMUONDataInterface(chLoader))
 {
-/// Constructor
-  fMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
-
-  // open the run loader
-  fRunLoader = AliRunLoader::Open(chLoader);
-  if (!fRunLoader) {
-    AliError(Form("no run loader found in file %s","galice.root" ));
-    return;
-  }
-  // initialize loader's
-  AliLoader *loader = fRunLoader->GetLoader("MUONLoader");
-
-  // initialize container
-  fMUONData  = new AliMUONData(loader,"MUON","MUON");
-
-   // Loading AliRun master
-  if (fRunLoader->GetAliRun() == 0x0) fRunLoader->LoadgAlice();
-
-  fRunLoader->LoadKinematics("READ");
-  fRunLoader->LoadTrackRefs("READ");
-  loader->LoadTracks("READ");
-
-  fReconstructibleTracks = 0; 
-  fRecoTracks = 0;
+  /// Normal ctor
 }
 
-//____________________________________________________________________
-AliMUONRecoCheck::AliMUONRecoCheck(const AliMUONRecoCheck& rhs)
- : TObject(rhs)
+//_____________________________________________________________________________
+AliMUONRecoCheck::~AliMUONRecoCheck()
 {
-/// Protected copy constructor
-
-  AliFatal("Not implemented.");
+  /// Destructor
+  delete fMCDataInterface;
+  delete fDataInterface;
 }
 
 //_____________________________________________________________________________
-AliMUONRecoCheck::~AliMUONRecoCheck()
+AliMUONVTrackStore* 
+AliMUONRecoCheck::ReconstructedTracks(Int_t event)
 {
-/// Destructor
-
-  fRunLoader->UnloadKinematics();
-  fRunLoader->UnloadTrackRefs();
-  fRunLoader->UnloadTracks();
-  fMuonTrackRef->Delete();
-  delete fMuonTrackRef;
-  delete fMUONData;
+  /// Return the reconstructed track store for a given event
+  return fDataInterface->TrackStore(event);
 }
 
-//________________________________________________________________________
-AliMUONRecoCheck& AliMUONRecoCheck::operator = (const AliMUONRecoCheck& rhs)
+//_____________________________________________________________________________
+AliMUONVTrackStore* 
+AliMUONRecoCheck::TrackRefs(Int_t event)
 {
-/// Protected assignement operator
-
-  if (this == &rhs) return *this;
-
-  AliFatal("Not implemented.");
-    
-  return *this;  
+  /// Return a track store containing the track references (converted into 
+  /// MUONTrack objects) for a given event
+  return MakeTrackRefs(event);
 }
 
 //_____________________________________________________________________________
-void AliMUONRecoCheck::MakeTrackRef()
+AliMUONVTrackStore* 
+AliMUONRecoCheck::ReconstructibleTracks(Int_t event)
 {
-/// Make reconstructible tracks
+  /// Return a track store containing the reconstructible tracks for a given event
+  AliMUONVTrackStore* tmp = MakeTrackRefs(event);
+  AliMUONVTrackStore* reconstructible = MakeReconstructibleTracks(*tmp);
+  delete tmp;
+  return reconstructible;
+}
 
-  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::MakeTrackRefs(Int_t event)
+{
+  /// Make reconstructible tracks
+    
+  AliMUONVTrackStore* trackRefStore = new AliMUONTrackStoreV1;
+  
+  Int_t nTrackRef = fMCDataInterface->NumberOfTrackRefs(event);
+  
+  Int_t trackSave(-999);
+  Int_t iHitMin(0);
   
-  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);  
-
-  Int_t nTrackRef = (Int_t)branch->GetEntries();
-  track = trackSave = -999;
   Bool_t isNewTrack;
-  Int_t iHitMin, iChamber;
-
-  trackParam = new AliMUONTrackParam();
-  hitForRec = new AliMUONHitForRec();
-  muonTrack = new AliMUONTrack();
-
-  for (Int_t iTrackRef  = 0; iTrackRef < nTrackRef; iTrackRef++) {
-    branch->GetEntry(iTrackRef);
+    
+  AliStack* stack = fMCDataInterface->Stack(event);
+  Int_t max = stack->GetNtrack();
+  
+  for (Int_t iTrackRef  = 0; iTrackRef < nTrackRef; ++iTrackRef) 
+  {
+    TClonesArray* trackRefs = fMCDataInterface->TrackRefs(event,iTrackRef);
     
     iHitMin = 0;
     isNewTrack = kTRUE;
     
     if (!trackRefs->GetEntries()) continue; 
-
+    
     while (isNewTrack) {
       
-      for (Int_t iHit = iHitMin; iHit < trackRefs->GetEntries(); iHit++) {
+      AliMUONTrack muonTrack;
       
-       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 != 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);  
-       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;
+      for (Int_t iHit = iHitMin; iHit < trackRefs->GetEntries(); ++iHit) 
+      {        
+        AliTrackReference* trackReference = static_cast<AliTrackReference*>(trackRefs->At(iHit));
+        
+        Float_t x = trackReference->X();
+        Float_t y = trackReference->Y();
+        Float_t z = trackReference->Z();
+        Float_t pX = trackReference->Px();
+        Float_t pY = trackReference->Py();
+        Float_t pZ = trackReference->Pz();
+        
+        Int_t 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;
+        }
+        
+        Float_t bendingSlope = 0;
+        Float_t nonBendingSlope = 0;
+        Float_t inverseBendingMomentum = 0;
+        
+        AliMUONTrackParam trackParam;
+        
+        // track parameters at hit
+        trackParam.SetBendingCoor(y);
+        trackParam.SetNonBendingCoor(x);
+        trackParam.SetZ(z);
+        
+        if (TMath::Abs(pZ) > 0) 
+        {
+          bendingSlope = pY/pZ;
+          nonBendingSlope = pX/pZ;
+        }
+        Float_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
+        if (pYZ >0) inverseBendingMomentum = 1/pYZ; 
+        
+        trackParam.SetBendingSlope(bendingSlope);
+        trackParam.SetNonBendingSlope(nonBendingSlope);
+        trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
+        
+        AliMUONHitForRec hitForRec;
+        
+        hitForRec.SetBendingCoor(y);
+        hitForRec.SetNonBendingCoor(x);
+        hitForRec.SetZ(z);
+        hitForRec.SetBendingReso2(0.0); 
+        hitForRec.SetNonBendingReso2(0.0);
+        Int_t detElemId = hitForRec.GetDetElemId();
+        Int_t iChamber;
+        if (detElemId) 
+        {
+          iChamber = detElemId / 100 - 1; 
+        }
+        else 
+        {
+          iChamber = AliMUONConstants::ChamberNumber(z);
+        }
+        hitForRec.SetChamberNumber(iChamber);
+        
+        muonTrack.AddTrackParamAtHit(&trackParam,0);
+        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);
+      TParticle* particle = stack->Particle(muonTrack.GetTrackID());
       
-       muonTrack->SetTrackParamAtVertex(trackParam);
+      if (particle) 
+      {        
+        Float_t x = particle->Vx();
+        Float_t y = particle->Vy();
+        Float_t z = particle->Vz();
+        Float_t pX = particle->Px();
+        Float_t pY = particle->Py();
+        Float_t pZ = particle->Pz();
+        
+        AliMUONTrackParam trackParam;
+        
+        trackParam.SetBendingCoor(y);
+        trackParam.SetNonBendingCoor(x);
+        trackParam.SetZ(z);
+        
+        Float_t bendingSlope = 0;
+        Float_t nonBendingSlope = 0;
+        Float_t inverseBendingMomentum = 0;
+                
+        if (TMath::Abs(pZ) > 0) 
+        {
+          bendingSlope = pY/pZ;
+          nonBendingSlope = pX/pZ;
+        }
+        
+        Float_t pYZ = TMath::Sqrt(pY*pY+pZ*pZ);
+        if (pYZ >0) inverseBendingMomentum = 1/pYZ;       
+        
+        TParticlePDG* ppdg = particle->GetPDG(1);
+        Int_t charge = (Int_t)(ppdg->Charge()/3.0);
+        inverseBendingMomentum *= charge;
+        
+        trackParam.SetBendingSlope(bendingSlope);
+        trackParam.SetNonBendingSlope(nonBendingSlope);
+        trackParam.SetInverseBendingMomentum(inverseBendingMomentum);
+        
+        muonTrack.SetTrackParamAtVertex(&trackParam);
       }
-      
-      AddMuonTrackReference(muonTrack);
-      muonTrack->ResetTrackParamAtHit();
-      muonTrack->ResetHitForRecAtHit();
-      
+
+      trackRefStore->Add(muonTrack);
     } // end while isNewTrack
-    
   }
   
-  CleanMuonTrackRef();
+  AliMUONVTrackStore* rv = CleanMuonTrackRef(*trackRefStore);
   
-  ReconstructibleTracks();
-
-  delete muonTrack;
-  delete trackParam;
-  delete hitForRec;
+  delete trackRefStore;
 
+  return rv;
 }
 
-//____________________________________________________________________________
-TClonesArray* AliMUONRecoCheck::GetTrackReco()
-{
-/// Return TClonesArray of reconstructed tracks
-
-  GetMUONData()->ResetRecTracks();
-  GetMUONData()->SetTreeAddress("RT");
-  fTrackReco = GetMUONData()->RecTracks(); 
-  GetMUONData()->GetRecTracks();
-  fRecoTracks = fTrackReco->GetEntriesFast();
-  return fTrackReco;
-}
 //_____________________________________________________________________________
-void AliMUONRecoCheck::PrintEvent() const
+Int_t
+AliMUONRecoCheck::NumberOfEvents() const
 {
-/// Debug facility
-
-  AliMUONTrack *track;
-  AliMUONHitForRec *hitForRec;
-  TClonesArray *  hitForRecAtHit = 0;
-  Float_t xRec,yRec,zRec;
-
-  Int_t nTrackRef = fMuonTrackRef->GetEntriesFast();
-  
-  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);
-    }
+  /// Return the number of events
+  if ( fDataInterface ) 
+  {
+    return fDataInterface->NumberOfEvents();
   }
+  return 0;
 }
 
 //_____________________________________________________________________________
-void AliMUONRecoCheck::ResetTracks() const
+AliMUONVTrackStore*
+AliMUONRecoCheck::CleanMuonTrackRef(const AliMUONVTrackStore& trackRefs)
 {
-/// 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) 
-
+  /// 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;
@@ -311,26 +278,29 @@ void AliMUONRecoCheck::CleanMuonTrackRef()
   Float_t bendingSlope,nonBendingSlope,bendingMomentum;
   Float_t bendingSlope1,nonBendingSlope1,bendingMomentum1;
   Float_t bendingSlope2,nonBendingSlope2,bendingMomentum2;
-  TClonesArray *newMuonTrackRef = new TClonesArray("AliMUONTrack", 10);
+  
+  AliMUONVTrackStore* newMuonTrackRef = static_cast<AliMUONVTrackStore*>(trackRefs.Create());
   Int_t iHit1;
-  Int_t iChamber = 0;
+  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);
+  TIter next(trackRefs.CreateIterator());
+  AliMUONTrack* track;
+  
+  while ( ( track = static_cast<AliMUONTrack*>(next())) ) 
+  {
     hitForRecAtHit = track->GetHitForRecAtHit();
     trackParamAtHit = track->GetTrackParamAtHit();
     trackParamAtVertex = track->GetTrackParamAtVertex();
     nTrackHits = hitForRecAtHit->GetEntriesFast();
+    AliMUONTrack trackNew;
     iHit1 = 0;
-    while (iHit1 < nTrackHits) {
+    while (iHit1 < nTrackHits) 
+    {
       hitForRec1 = (AliMUONHitForRec*) hitForRecAtHit->At(iHit1); 
       trackParam1 = (AliMUONTrackParam*) trackParamAtHit->At(iHit1); 
       xRec1  = hitForRec1->GetNonBendingCoor();
@@ -343,34 +313,36 @@ void AliMUONRecoCheck::CleanMuonTrackRef()
       nonBendingSlope1 = trackParam1->GetNonBendingSlope();
       bendingMomentum1 = 0;
       if (TMath::Abs(trackParam1->GetInverseBendingMomentum()) > 0)
-       bendingMomentum1 = 1./trackParam1->GetInverseBendingMomentum();
+        bendingMomentum1 = 1./trackParam1->GetInverseBendingMomentum();
       bendingSlope = bendingSlope1;
       nonBendingSlope = nonBendingSlope1;
       bendingMomentum = bendingMomentum1;
       nRec = 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;
-       }
-       
+      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;
+        }
+        
       } // end iHit2
       xRec /= (Float_t)nRec;
       yRec /= (Float_t)nRec;
@@ -382,7 +354,9 @@ void AliMUONRecoCheck::CleanMuonTrackRef()
       hitForRec->SetNonBendingCoor(xRec);
       hitForRec->SetBendingCoor(yRec);
       hitForRec->SetZ(zRec);
-      iChamber = AliMUONConstants::ChamberNumber(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); 
@@ -392,77 +366,75 @@ void AliMUONRecoCheck::CleanMuonTrackRef()
       trackParam->SetNonBendingSlope(nonBendingSlope);
       trackParam->SetBendingSlope(bendingSlope);
       if (TMath::Abs(bendingMomentum) > 0)
-       trackParam->SetInverseBendingMomentum(1./bendingMomentum);
-
-      trackNew->AddHitForRecAtHit(hitForRec);
-      trackNew->AddTrackParamAtHit(trackParam);
+        trackParam->SetInverseBendingMomentum(1./bendingMomentum);
+      
+      trackNew.AddHitForRecAtHit(hitForRec);
+      trackNew.AddTrackParamAtHit(trackParam,0);
       
       iHit1++;
     } // end iHit1
-
-    trackNew->SetTrackID(track->GetTrackID());
-    trackNew->SetTrackParamAtVertex(trackParamAtVertex);
-    {new ((*newMuonTrackRef)[newMuonTrackRef->GetEntriesFast()]) AliMUONTrack(*trackNew);}    
-    trackNew->ResetHitForRecAtHit();
-    trackNew->ResetTrackParamAtHit();
+    
+    trackNew.SetTrackID(track->GetTrackID());
+    trackNew.SetTrackParamAtVertex(trackParamAtVertex);
+    newMuonTrackRef->Add(trackNew);
     
   } // 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;
-  
+  return newMuonTrackRef;  
 }
 
 //_____________________________________________________________________________
-void AliMUONRecoCheck::ReconstructibleTracks()
+AliMUONVTrackStore*
+AliMUONRecoCheck::MakeReconstructibleTracks(const AliMUONVTrackStore& trackRefs)
 {
-/// Calculate the number of reconstructible tracks
+  /// Calculate the number of reconstructible tracks
+  
+  AliMUONVTrackStore* reconstructibleStore = static_cast<AliMUONVTrackStore*>(trackRefs.Create());
   
-  AliMUONTrack* track;
   TClonesArray* hitForRecAtHit = NULL;
   AliMUONHitForRec* hitForRec;
   Float_t zRec;
-  Int_t nTrackRef, nTrackHits;
+  Int_t nTrackHits;
   Int_t isChamberInTrack[10];
   Int_t iChamber = 0;
   Bool_t isTrackOK = kTRUE;
+    
+  TIter next(trackRefs.CreateIterator());
+  AliMUONTrack* track;
 
-  fReconstructibleTracks = 0;
-
-  nTrackRef = fMuonTrackRef->GetEntriesFast();
-  for (Int_t index = 0; index < nTrackRef; index++) {
-    track = (AliMUONTrack*)fMuonTrackRef->At(index);
+  while ( ( track = static_cast<AliMUONTrack*>(next()) ) )
+  {
     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();
+      iChamber = hitForRec->GetChamberNumber(); 
       if (iChamber < 0 || iChamber > 10) continue;
       isChamberInTrack[iChamber] = 1;
-      
-    }
-    // track is reconstructible if the particle is crossing every tracking chambers
+    } 
+    // track is reconstructible if the particle is depositing a hit
+    // in the following chamber combinations:
+    
     isTrackOK = kTRUE;
-    for (Int_t ch = 0; ch < 10; ch++) {
-      if (!isChamberInTrack[ch]) isTrackOK = kFALSE;
-    }    if (isTrackOK) fReconstructibleTracks++;
-    if (!isTrackOK) fMuonTrackRef->Remove(track); // remove non reconstructible tracks
+    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) 
+    {
+      reconstructibleStore->Add(*track);
+    }
   }
-  fMuonTrackRef->Compress();
-}
 
+  return reconstructibleStore;
+}