New class added to check the reconstruction of muon tracks using reference tracks
authorcussonno <cussonno@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 20 Oct 2004 05:35:39 +0000 (05:35 +0000)
committercussonno <cussonno@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 20 Oct 2004 05:35:39 +0000 (05:35 +0000)
14 files changed:
MUON/AliMUONEventReconstructor.cxx
MUON/AliMUONEventReconstructor.h
MUON/AliMUONHitForRec.cxx
MUON/AliMUONHitForRec.h
MUON/AliMUONRecoCheck.cxx [new file with mode: 0644]
MUON/AliMUONRecoCheck.h [new file with mode: 0644]
MUON/AliMUONTrack.cxx
MUON/AliMUONTrack.h
MUON/AliMUONTrackParam.cxx
MUON/AliMUONTrackParam.h
MUON/MUONLinkDef.h
MUON/MUONRecoCheck.C [new file with mode: 0644]
MUON/README
MUON/libMUON.pkg

index 2c21a553838bda807413b407eef0240eda17e8fa..098cab1ec5fa20861dbe1e4c4b399ecd1cfcae35 100644 (file)
@@ -987,6 +987,7 @@ void AliMUONEventReconstructor::MakeTracks(void)
     // Remove double tracks
     RemoveDoubleTracks();
     UpdateTrackParamAtHit();
+    UpdateHitForRecAtHit();
   }
   return;
 }
@@ -1547,6 +1548,26 @@ void AliMUONEventReconstructor::UpdateTrackParamAtHit()
   return;
 }
 
+  //__________________________________________________________________________
+void AliMUONEventReconstructor::UpdateHitForRecAtHit()
+{
+  // Set cluster parameterss after track fitting. Fill fHitForRecAtHit of AliMUONTrack's
+  AliMUONTrack *track;
+  AliMUONTrackHit *trackHit;
+  AliMUONHitForRec *hitForRec;
+  track = (AliMUONTrack*) fRecTracksPtr->First();
+  while (track) {
+    trackHit = (AliMUONTrackHit*) (track->GetTrackHitsPtr())->First();
+    while (trackHit) {
+      hitForRec = trackHit->GetHitForRecPtr();
+      track->AddHitForRecAtHit(hitForRec);
+      trackHit = (AliMUONTrackHit*) (track->GetTrackHitsPtr())->After(trackHit); 
+    } // trackHit    
+    track = (AliMUONTrack*) fRecTracksPtr->After(track);
+  } // track
+  return;
+}
+
   //__________________________________________________________________________
 void AliMUONEventReconstructor::EventDump(void)
 {
index 6b2b9c7a9460f616fc035eefdbcd2a842060c9f0..56f1efaadea7401b22df7abc3c062d50bdcd9164 100644 (file)
@@ -219,6 +219,7 @@ class AliMUONEventReconstructor : public TObject {
   void FollowTracks(void);
   void RemoveDoubleTracks(void);
   void UpdateTrackParamAtHit(void);
+  void UpdateHitForRecAtHit(void);
   void ValidateTracksWithTrigger(void);
 
 
index 67f23a1a8f610ed71ef5c0f5fa747c3acdf977a5..528ad4a2c66ac99d20a5c9231ac299b314d10590 100644 (file)
@@ -98,24 +98,45 @@ AliMUONHitForRec::AliMUONHitForRec(AliMUONRawCluster* theRawCluster)
 }
 
   //__________________________________________________________________________
-AliMUONHitForRec::AliMUONHitForRec (const AliMUONHitForRec& rhs)
-  : TObject(rhs)
+AliMUONHitForRec::AliMUONHitForRec (const AliMUONHitForRec& theMUONHitForRec)
+  :  TObject(theMUONHitForRec)
 {
-// Protected copy constructor
+  fBendingCoor = theMUONHitForRec.fBendingCoor;
+  fNonBendingCoor = theMUONHitForRec.fNonBendingCoor;
+  fZ = theMUONHitForRec.fZ;
+  fBendingReso2 = theMUONHitForRec.fBendingReso2;
+  fNonBendingReso2 = theMUONHitForRec.fNonBendingReso2;
+  fChamberNumber = theMUONHitForRec.fChamberNumber;
+  fHitNumber = theMUONHitForRec.fHitNumber;
+  fTHTrack = theMUONHitForRec.fTHTrack;
+  fGeantSignal = theMUONHitForRec.fGeantSignal;
+  fIndexOfFirstSegment = theMUONHitForRec.fIndexOfFirstSegment;
+  fNSegments = theMUONHitForRec.fNSegments;
+  fFirstTrackHitPtr = theMUONHitForRec.fFirstTrackHitPtr;
+  fLastTrackHitPtr = theMUONHitForRec.fLastTrackHitPtr;
+  fNTrackHits = theMUONHitForRec.fNTrackHits;
 
-  AliFatal( "Not implemented.");
 }
 
   //__________________________________________________________________________
-AliMUONHitForRec & AliMUONHitForRec::operator=(const AliMUONHitForRec& rhs)
+AliMUONHitForRec & AliMUONHitForRec::operator=(const AliMUONHitForRec& theMUONHitForRec)
 {
-// Protected assignement operator
+  fBendingCoor = theMUONHitForRec.fBendingCoor;
+  fNonBendingCoor = theMUONHitForRec.fNonBendingCoor;
+  fZ = theMUONHitForRec.fZ;
+  fBendingReso2 = theMUONHitForRec.fBendingReso2;
+  fNonBendingReso2 = theMUONHitForRec.fNonBendingReso2;
+  fChamberNumber = theMUONHitForRec.fChamberNumber;
+  fHitNumber = theMUONHitForRec.fHitNumber;
+  fTHTrack = theMUONHitForRec.fTHTrack;
+  fGeantSignal = theMUONHitForRec.fGeantSignal;
+  fIndexOfFirstSegment = theMUONHitForRec.fIndexOfFirstSegment;
+  fNSegments = theMUONHitForRec.fNSegments;
+  fFirstTrackHitPtr = theMUONHitForRec.fFirstTrackHitPtr;
+  fLastTrackHitPtr = theMUONHitForRec.fLastTrackHitPtr;
+  fNTrackHits = theMUONHitForRec.fNTrackHits;
 
-  if (this == &rhs) return *this;
-
-  AliFatal( "Not implemented.");
-    
-  return *this;  
+  return *this;
 }
   //__________________________________________________________________________
 /*AZ
index 73759cac7b74cc394f8b92c3b01431a802b491d6..e412fa790c6fb9dae80ea01477241e312938c63c 100644 (file)
@@ -17,6 +17,8 @@ class AliMUONHitForRec : public TObject {
  public:
   AliMUONHitForRec(); // Constructor
   virtual ~AliMUONHitForRec(){} // Destructor
+  AliMUONHitForRec (const AliMUONHitForRec& AliMUONHitForRec); // copy constructor
+  AliMUONHitForRec& operator=(const AliMUONHitForRec& AliMUONHitForRec); // assignment operator
   AliMUONHitForRec(AliMUONHit* mHit); // Constructor from GEANT hit
   AliMUONHitForRec(AliMUONRawCluster* theRawCluster); // Constructor from raw cluster
 
@@ -58,9 +60,6 @@ class AliMUONHitForRec : public TObject {
   Bool_t IsSortable() const { return kTRUE; }
   Int_t Compare(const TObject* HitForRec) const; // "Compare" function for sorting
 
- protected:
-  AliMUONHitForRec (const AliMUONHitForRec& AliMUONHitForRec); // copy constructor
-  AliMUONHitForRec& operator=(const AliMUONHitForRec& AliMUONHitForRec); // assignment operator
 
  private:
   Double_t fBendingCoor; // coordinate (cm) in bending plane
@@ -78,13 +77,13 @@ class AliMUONHitForRec : public TObject {
   Int_t fGeantSignal; // Geant signal (1) or background (0)
 
   // links forward to the segment(s) if HitForRec in first chamber of a station
-  Int_t fIndexOfFirstSegment; // index of first Segment
-  Int_t fNSegments; // number of Segments
+  Int_t fIndexOfFirstSegment; //! index of first Segment
+  Int_t fNSegments; //! number of Segments
 
   // links forward to reconstructed track hits
-  AliMUONTrackHit *fFirstTrackHitPtr ; // pointer to first TrackHit made with HitForRec
-  AliMUONTrackHit *fLastTrackHitPtr ; // pointer to last TrackHit made with HitForRec
-  Int_t fNTrackHits; // number of TrackHit's made with HitForRec
+  AliMUONTrackHit *fFirstTrackHitPtr ; //! pointer to first TrackHit made with HitForRec
+  AliMUONTrackHit *fLastTrackHitPtr ; //! pointer to last TrackHit made with HitForRec
+  Int_t fNTrackHits; //! number of TrackHit's made with HitForRec
   
   ClassDef(AliMUONHitForRec, 1) // Hit for reconstruction in ALICE dimuon spectrometer
     };
diff --git a/MUON/AliMUONRecoCheck.cxx b/MUON/AliMUONRecoCheck.cxx
new file mode 100644 (file)
index 0000000..3ac9ff5
--- /dev/null
@@ -0,0 +1,452 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+//////////////////////////////////////////////////////////////////////////////////////
+//                                                                     
+// 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.     
+//                                                                      
+//////////////////////////////////////////////////////////////////////////////////////
+
+
+#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 "AliMUONTrack.h"
+#include "AliMUONData.h"
+#include "AliMUONChamber.h"
+#include "AliMUONConstants.h"
+
+ClassImp(AliMUONRecoCheck)
+
+//_____________________________________________________________________________
+AliMUONRecoCheck::AliMUONRecoCheck(Char_t *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;
+}
+
+//_____________________________________________________________________________
+AliMUONRecoCheck::~AliMUONRecoCheck()
+{
+  fRunLoader->UnloadKinematics();
+  fRunLoader->UnloadTrackRefs();
+  fRunLoader->UnloadTracks();
+  delete fMuonTrackRef;
+  delete fMUONData;
+}
+
+//_____________________________________________________________________________
+void AliMUONRecoCheck::MakeTrackRef()
+{
+  // Make reconstructible tracks
+
+  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;
+  
+  TTree* treeTR = fRunLoader->TreeTR();
+  if (treeTR == NULL) return;
+
+  TBranch* branch = treeTR->GetBranch("MUON");
+  if (branch == NULL) return;
+
+  TClonesArray* trackRefs = new TClonesArray("AliTrackReference", 10);
+  branch->SetAddress(&trackRefs);
+
+  Int_t nTrackRef = (Int_t)treeTR->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++) {
+    treeTR->GetEntry(iTrackRef);
+    
+    iHitMin = 0;
+    isNewTrack = kTRUE;
+
+    while (isNewTrack) {
+
+      for (Int_t iHit = iHitMin; iHit < trackRefs->GetEntries(); 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 != 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 = 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);
+      
+       muonTrack->SetTrackParamAtVertex(trackParam);
+      }
+      
+      AddMuonTrackReference(muonTrack);
+      muonTrack->ResetTrackParamAtHit();
+      muonTrack->ResetHitForRecAtHit();
+      
+    } // end while isNewTrack
+    
+  }
+  
+  CleanMuonTrackRef();
+  
+  ReconstructibleTracks();
+
+  delete muonTrack;
+  delete trackParam;
+  delete hitForRec;
+  delete trackRefs;
+
+}
+
+//____________________________________________________________________________
+TClonesArray* AliMUONRecoCheck::GetTrackReco()
+{
+  // Return TClonesArray of reconstructed tracks
+
+  GetMUONData()->SetTreeAddress("RT");
+  fTrackReco = GetMUONData()->RecTracks(); 
+  GetMUONData()->GetRecTracks();
+  fRecoTracks = fTrackReco->GetEntriesFast();
+  return fTrackReco;
+}
+//_____________________________________________________________________________
+void AliMUONRecoCheck::PrintEvent() 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);
+    }
+  }
+}
+
+//_____________________________________________________________________________
+void AliMUONRecoCheck::ResetTracks() const
+{
+  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;
+  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;
+    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 = 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 = 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;
+      zRec /= (Float_t)nRec;
+      bendingSlope /= (Float_t)nRec;
+      nonBendingSlope /= (Float_t)nRec;
+      bendingMomentum /= (Float_t)nRec;
+      
+      hitForRec->SetNonBendingCoor(xRec);
+      hitForRec->SetBendingCoor(yRec);
+      hitForRec->SetZ(zRec);
+      iChamber = 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);
+      trackParam->SetInverseBendingMomentum(1./bendingMomentum);
+
+      trackNew->AddHitForRecAtHit(hitForRec);
+      trackNew->AddTrackParamAtHit(trackParam);
+      
+      iHit1++;
+    } // end iHit1
+
+    trackNew->SetTrackID(track->GetTrackID());
+    trackNew->SetTrackParamAtVertex(trackParamAtVertex);
+    {new ((*newMuonTrackRef)[newMuonTrackRef->GetEntriesFast()]) AliMUONTrack(*trackNew);}    
+    trackNew->ResetHitForRecAtHit();
+    trackNew->ResetTrackParamAtHit();
+    
+  } // 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;
+  delete newMuonTrackRef;
+  
+}
+
+//_____________________________________________________________________________
+void AliMUONRecoCheck::ReconstructibleTracks()
+{
+  // calculate the number of reconstructible tracks
+  
+  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;
+      
+    }
+    // track is reconstructible if the particle is crossing every tracking chambers
+    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
+  }
+  fMuonTrackRef->Compress();
+}
+
+
+//_____________________________________________________________________________
+Int_t AliMUONRecoCheck::ChamberNumber(Float_t z) const
+{
+  // return chamber number according z position of hit. Should be taken from geometry ?
+  Float_t dMaxChamber =  AliMUONConstants::DzSlat() + AliMUONConstants::DzCh() + 0.25; // cm st 3 &4 & 5
+  if ( z >  (AliMUONConstants::DefaultChamberZ(4)+50.)) dMaxChamber = 7.; // cm stations 1 & 2
+  Int_t iChamber;
+
+  for (iChamber = 0; iChamber < 10; iChamber++) {
+    
+    if (TMath::Abs(z-AliMUONConstants::DefaultChamberZ(iChamber)) < dMaxChamber) {
+      return iChamber;
+    }
+  }
+  return -1;
+}
diff --git a/MUON/AliMUONRecoCheck.h b/MUON/AliMUONRecoCheck.h
new file mode 100644 (file)
index 0000000..48b6d1f
--- /dev/null
@@ -0,0 +1,61 @@
+#ifndef ALIMUONRECOCHECK_H
+#define ALIMUONRECOCHECK_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//////////////////////////////////////////////////////////////////////////
+//                                                                      //
+// AliMUONRecoCheck                                                     //
+//                                                                      //
+//////////////////////////////////////////////////////////////////////////
+#include <TObject.h>
+#include "AliMUONTrack.h"
+
+class TClonesArray;
+class AliMUONData;
+class AliRunLoader;
+
+
+class AliMUONRecoCheck : public TObject 
+{
+public:
+  AliMUONRecoCheck(Char_t *chLoader);
+  virtual          ~AliMUONRecoCheck();
+
+  AliMUONData*  GetMUONData() {return fMUONData;}
+  void MakeTrackRef();
+  void AddMuonTrackReference(const AliMUONTrack *muonTrack) 
+    {new ((*fMuonTrackRef)[fMuonTrackRef->GetEntriesFast()]) AliMUONTrack(*muonTrack);}
+  void PrintEvent() const;
+  void ResetTracks() const;
+  AliRunLoader* GetRunLoader() {return fRunLoader;}
+  void CleanMuonTrackRef();
+  void ReconstructibleTracks();
+  Int_t GetNumberOfReconstuctibleTracks() {return fReconstructibleTracks;}
+  Int_t ChamberNumber(Float_t z) const;
+  Int_t GetNumberOfRecoTracks() {return fRecoTracks;}
+  TClonesArray *GetTrackReco();
+  TClonesArray *GetMuonTrackRef() {return fMuonTrackRef;}
+
+private:
+  
+  AliRunLoader* fRunLoader;     // alice run loader 
+  AliMUONData*  fMUONData;      // Data container for MUON subsystem 
+  TClonesArray* fMuonTrackRef;  // reference muon tracks
+  TClonesArray* fTrackReco;     // reconstructed muon tracks
+  Int_t fReconstructibleTracks; // number of reconstructible tracks 
+  Int_t fRecoTracks; // number of reconstructed tracks 
+
+  ClassDef(AliMUONRecoCheck, 0)   //Utility class to check reconstruction
+};
+
+#endif
+
+
+
+
+
+
+
+
index 03375183a518926484d898c084756ebb22b3e0b3..1c61adf31dce189f6de5d642f4d862b0ff2ab987 100644 (file)
@@ -62,6 +62,8 @@ AliMUONTrack::AliMUONTrack()
   fEventReconstructor = 0;
   fTrackHitsPtr = new TObjArray(10);
   fTrackParamAtHit = new TClonesArray("AliMUONTrackParam",10);  
+  fHitForRecAtHit = new TClonesArray("AliMUONHitForRec",10); 
+  fTrackID = 0;
 }
 
   //__________________________________________________________________________
@@ -78,6 +80,7 @@ AliMUONTrack::AliMUONTrack(AliMUONSegment* BegSegment, AliMUONSegment* EndSegmen
   fTrackHitsPtr->Sort(); // sort TrackHits according to increasing Z
   SetTrackParamAtVertex(); // set track parameters at vertex
   fTrackParamAtHit = new TClonesArray("AliMUONTrackParam",10);
+  fHitForRecAtHit = new TClonesArray("AliMUONHitForRec",10);
   // set fit conditions...
   fFitMCS = 0;
   fFitNParam = 3;
@@ -85,6 +88,7 @@ AliMUONTrack::AliMUONTrack(AliMUONSegment* BegSegment, AliMUONSegment* EndSegmen
   fFitFMin = -1.0;
   fMatchTrigger = kFALSE;
   fChi2MatchTrigger = 0;
+  fTrackID = 0;
   return;
 }
 
@@ -102,6 +106,7 @@ AliMUONTrack::AliMUONTrack(AliMUONSegment* Segment, AliMUONHitForRec* HitForRec,
   fTrackHitsPtr->Sort(); // sort TrackHits according to increasing Z
   SetTrackParamAtVertex(); // set track parameters at vertex
   fTrackParamAtHit = new TClonesArray("AliMUONTrackParam",10);
+  fHitForRecAtHit = new TClonesArray("AliMUONHitForRec",10);
   // set fit conditions...
   fFitMCS = 0;
   fFitNParam = 3;
@@ -109,6 +114,7 @@ AliMUONTrack::AliMUONTrack(AliMUONSegment* Segment, AliMUONHitForRec* HitForRec,
   fFitFMin = -1.0;
   fMatchTrigger = kFALSE;
   fChi2MatchTrigger = 0;
+  fTrackID = 0;
   return;
 }
 
@@ -126,6 +132,12 @@ AliMUONTrack::~AliMUONTrack()
     delete fTrackParamAtHit;
     fTrackParamAtHit = NULL;
   }
+
+  if (fHitForRecAtHit) {
+    // delete the TClonesArray of pointers to HitForRec
+    delete fHitForRecAtHit;
+    fHitForRecAtHit = NULL;
+  }
 }
 
   //__________________________________________________________________________
@@ -153,6 +165,13 @@ AliMUONTrack::AliMUONTrack (const AliMUONTrack& theMUONTrack)
        AliMUONTrackParam(*(AliMUONTrackParam*)(theMUONTrack.fTrackParamAtHit)->At(index));}
   }
 
+  // necessary to make a copy of the objects and not only the pointers in TClonesArray.
+  fHitForRecAtHit  =  new TClonesArray("AliMUONHitForRec",10);
+  for (Int_t index = 0; index < (theMUONTrack.fHitForRecAtHit)->GetEntriesFast(); index++) {
+    {new ((*fHitForRecAtHit)[fHitForRecAtHit->GetEntriesFast()]) 
+       AliMUONHitForRec(*(AliMUONHitForRec*)(theMUONTrack.fHitForRecAtHit)->At(index));}
+  }
+
   fNTrackHits       =  theMUONTrack.fNTrackHits;
   fFitMCS           =  theMUONTrack.fFitMCS;
   fFitNParam        =  theMUONTrack.fFitNParam;
@@ -160,6 +179,7 @@ AliMUONTrack::AliMUONTrack (const AliMUONTrack& theMUONTrack)
   fFitStart         =  theMUONTrack.fFitStart;
   fMatchTrigger     =  theMUONTrack.fMatchTrigger;
   fChi2MatchTrigger =  theMUONTrack.fChi2MatchTrigger;
+  fTrackID          =  theMUONTrack.fTrackID;
 }
 
   //__________________________________________________________________________
@@ -192,6 +212,13 @@ AliMUONTrack & AliMUONTrack::operator=(const AliMUONTrack& theMUONTrack)
        AliMUONTrackParam(*(AliMUONTrackParam*)(theMUONTrack.fTrackParamAtHit)->At(index));}
   }
 
+  // necessary to make a copy of the objects and not only the pointers in TClonesArray.
+  fHitForRecAtHit  =  new TClonesArray("AliMUONHitForRec",10);
+  for (Int_t index = 0; index < (theMUONTrack.fHitForRecAtHit)->GetEntriesFast(); index++) {
+    {new ((*fHitForRecAtHit)[fHitForRecAtHit->GetEntriesFast()]) 
+       AliMUONHitForRec(*(AliMUONHitForRec*)(theMUONTrack.fHitForRecAtHit)->At(index));}
+  }
+
   fNTrackHits         =  theMUONTrack.fNTrackHits;
   fFitMCS             =  theMUONTrack.fFitMCS;
   fFitNParam          =  theMUONTrack.fFitNParam;
@@ -199,6 +226,7 @@ AliMUONTrack & AliMUONTrack::operator=(const AliMUONTrack& theMUONTrack)
   fFitStart           =  theMUONTrack.fFitStart;
   fMatchTrigger       =  theMUONTrack.fMatchTrigger;
   fChi2MatchTrigger   =  theMUONTrack.fChi2MatchTrigger;
+  fTrackID            =  theMUONTrack.fTrackID;
 
   return *this;
 }
@@ -306,6 +334,45 @@ void AliMUONTrack::RecursiveDump(void) const
   }
   return;
 }
+  
+  //__________________________________________________________________________
+Bool_t* AliMUONTrack::CompatibleTrack(AliMUONTrack * Track, Double_t Sigma2Cut) const
+{
+  // Return kTRUE/kFALSE for each chamber if hit is compatible or not 
+  TClonesArray *hitArray, *thisHitArray;
+  AliMUONHitForRec *hit, *thisHit;
+  Int_t chamberNumber;
+  Float_t deltaZ;
+  Float_t deltaZMax = 1.; // 1 cm
+  Float_t chi2 = 0;
+  Bool_t *nCompHit = new Bool_t[AliMUONConstants::NTrackingCh()]; 
+
+  for ( Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
+    nCompHit[ch] = kFALSE;
+  }
+
+  thisHitArray = this->GetHitForRecAtHit();
+
+  hitArray =  Track->GetHitForRecAtHit();
+
+  for (Int_t iHthis = 0; iHthis < thisHitArray->GetEntriesFast(); iHthis++) {
+    thisHit = (AliMUONHitForRec*) thisHitArray->At(iHthis);
+    chamberNumber = thisHit->GetChamberNumber();
+    if (chamberNumber < 0 || chamberNumber > AliMUONConstants::NTrackingCh()) continue; 
+    nCompHit[chamberNumber] = kFALSE;
+    for (Int_t iH = 0; iH < hitArray->GetEntriesFast(); iH++) {
+      hit = (AliMUONHitForRec*) hitArray->At(iH);
+      deltaZ = TMath::Abs(thisHit->GetZ() - hit->GetZ());
+      chi2 = thisHit->NormalizedChi2WithHitForRec(hit,Sigma2Cut); // set cut to 4 sigmas
+      if (chi2 < 3. && deltaZ < deltaZMax) {
+       nCompHit[chamberNumber] = kTRUE;
+       break;
+      }
+    }  
+  }
+  
+  return nCompHit;
+}
 
   //__________________________________________________________________________
 Int_t AliMUONTrack::HitsInCommon(AliMUONTrack* Track) const
index 57a89e9f4b2d476ee930ec59e654693352a7ba5e..0be63dc47bc3905720e41fe45939101b0410b01b 100644 (file)
@@ -14,6 +14,7 @@
 #include <TClonesArray.h>
 
 #include "AliMUONTrackParam.h" // object belongs to the class
+#include "AliMUONHitForRec.h"  // object belongs to the class
 
 //const Int_t kMaxTrackingChamber=10;
         // not used
@@ -41,9 +42,13 @@ class AliMUONTrack : public TObject
   void SetTrackParamAtVertex(void); // Set track parameters at vertex from last stations 4 & 5
   void SetTrackParamAtVertex(AliMUONTrackParam* TrackParam) {fTrackParamAtVertex = *TrackParam;}
   TClonesArray *GetTrackParamAtHit(void) const {return fTrackParamAtHit;}
+  TClonesArray *GetHitForRecAtHit(void) const {return fHitForRecAtHit;}
   void ResetTrackParamAtHit(void) { fTrackParamAtHit->Delete(); }
+  void ResetHitForRecAtHit(void) { fHitForRecAtHit->Delete(); }
   void AddTrackParamAtHit(const AliMUONTrackParam *trackParam) 
     {new ((*fTrackParamAtHit)[fTrackParamAtHit->GetEntriesFast()]) AliMUONTrackParam(*trackParam);}
+  void AddHitForRecAtHit(const AliMUONHitForRec *hitForRec) 
+    {new ((*fHitForRecAtHit)[fHitForRecAtHit->GetEntriesFast()]) AliMUONHitForRec(*hitForRec);}
 
   TObjArray* GetTrackHitsPtr(void) const {return fTrackHitsPtr;}
   Int_t GetNTrackHits(void) const {return fNTrackHits;}
@@ -66,6 +71,10 @@ class AliMUONTrack : public TObject
   void SetTrackParamAtHit(Int_t indexHit, AliMUONTrackParam *TrackParam) const;
   Int_t HitsInCommon(AliMUONTrack* Track) const;
   void MatchTriggerTrack(TClonesArray* TriggerTrackArray);
+  Bool_t* CompatibleTrack(AliMUONTrack* Track, Double_t Sigma2Cut) const; // return array of compatible chamber
+  
+  Int_t GetTrackID() const {return fTrackID;}
+  void SetTrackID(Int_t trackID) {fTrackID = trackID;}
 
   static TVirtualFitter* Fitter(void) {return fgFitter;}
 
@@ -75,6 +84,7 @@ class AliMUONTrack : public TObject
   AliMUONEventReconstructor* fEventReconstructor; //!   Pointer to EventReconstructor
   AliMUONTrackParam fTrackParamAtVertex; // Track parameters at vertex
   TClonesArray *fTrackParamAtHit; // Track parameters at hit
+  TClonesArray *fHitForRecAtHit; // Cluster parameters at hit
   TObjArray *fTrackHitsPtr; //!  Pointer to array of pointers to TrackHit's
   Int_t fNTrackHits; // Number of TrackHit's
   Int_t fFitMCS; // 0(1) for fit without(with) multiple Coulomb scattering
@@ -84,6 +94,8 @@ class AliMUONTrack : public TObject
   Bool_t fMatchTrigger; // 1 if track matches with trigger track, 0 if not
   Double_t fChi2MatchTrigger; // chi2 of trigger/track matching 
  
+  Int_t fTrackID; // track ID = track number in TrackRefs
+
   ClassDef(AliMUONTrack, 2) // Reconstructed track in ALICE dimuon spectrometer
     };
        
index f0019f7313a33102234e82de3c3ca220094e587f..922c8f081671554ee19eedc396121711145d0ead 100644 (file)
@@ -403,10 +403,10 @@ void AliMUONTrackParam::BransonCorrection()
 
   pYZ = TMath::Abs(1.0 / fInverseBendingMomentum);
   sign = 1;      
-  if (fInverseBendingMomentum < 0) sign = -1;     
-  pZ = -pYZ / (TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope)); // spectro (z<0)
-  pX = pZ * fNonBendingSlope
-  pY = pZ * fBendingSlope
+  if (fInverseBendingMomentum < 0) sign = -1;  
+  pZ = Pz();
+  pX = Px()
+  pY = Py()
   pTotal = TMath::Sqrt(pYZ *pYZ + pX * pX);
   xEndAbsorber = fNonBendingCoor; 
   yEndAbsorber = fBendingCoor; 
@@ -494,9 +494,9 @@ void AliMUONTrackParam::FieldCorrection(Double_t Z)
   pYZ = TMath::Abs(1.0 / fInverseBendingMomentum);
   c = TMath::Sign(1.0,fInverseBendingMomentum); // particle charge 
  
-  pZ = -pYZ / (TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope));  // spectro. (z<0)
-  pX = pZ * fNonBendingSlope
-  pY = pZ * fBendingSlope;
+  pZ = Pz();
+  pX = Px()
+  pY = Py();
   pT = TMath::Sqrt(pX*pX+pY*pY);
 
   if (TMath::Abs(pZ) <= 0) return;
@@ -520,4 +520,53 @@ void AliMUONTrackParam::FieldCorrection(Double_t Z)
   
   fInverseBendingMomentum = c / TMath::Sqrt(pYNew*pYNew+pZ*pZ);
  
+}
+  //__________________________________________________________________________
+Double_t AliMUONTrackParam::Px()
+{
+  // return px from track paramaters
+  Double_t pYZ, pZ, pX;
+  pYZ = 0;
+  if (  TMath::Abs(fInverseBendingMomentum) > 0 )
+    pYZ = TMath::Abs(1.0 / fInverseBendingMomentum);
+  pZ = -pYZ / (TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope));  // spectro. (z<0)
+  pX = pZ * fNonBendingSlope; 
+  return pX;
+}
+  //__________________________________________________________________________
+Double_t AliMUONTrackParam::Py()
+{
+  // return px from track paramaters
+  Double_t pYZ, pZ, pY;
+  pYZ = 0;
+  if (  TMath::Abs(fInverseBendingMomentum) > 0 )
+    pYZ = TMath::Abs(1.0 / fInverseBendingMomentum);
+  pZ = -pYZ / (TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope));  // spectro. (z<0)
+  pY = pZ * fBendingSlope; 
+  return pY;
+}
+  //__________________________________________________________________________
+Double_t AliMUONTrackParam::Pz()
+{
+  // return px from track paramaters
+  Double_t pYZ, pZ;
+  pYZ = 0;
+  if (  TMath::Abs(fInverseBendingMomentum) > 0 )
+    pYZ = TMath::Abs(1.0 / fInverseBendingMomentum);
+  pZ = -pYZ / (TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope));  // spectro. (z<0)
+  return pZ;
+}
+  //__________________________________________________________________________
+Double_t AliMUONTrackParam::P()
+{
+  // return p from track paramaters
+  Double_t  pYZ, pZ, p;
+  pYZ = 0;
+  if (  TMath::Abs(fInverseBendingMomentum) > 0 )
+    pYZ = TMath::Abs(1.0 / fInverseBendingMomentum);
+  pZ = -pYZ / (TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope));  // spectro. (z<0)
+  p = TMath::Abs(pZ) * 
+    TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope + fNonBendingSlope * fNonBendingSlope);
+  return p;
+  
 }
index ce09432017befe585e7ef3e6ae33c9d11e1026d5..1d010e1b9fe7ce3d7e0b24cf16127313e332f022 100644 (file)
@@ -33,6 +33,10 @@ class AliMUONTrackParam : public TObject
   void SetBendingCoor(Double_t BendingCoor) {fBendingCoor = BendingCoor;}
   Double_t GetNonBendingCoor(void) const {return fNonBendingCoor;}
   void SetNonBendingCoor(Double_t NonBendingCoor) {fNonBendingCoor = NonBendingCoor;}
+  Double_t Px();  // return px
+  Double_t Py();  // return py
+  Double_t Pz();  // return pz
+  Double_t P();   // return total momentum
 
   void ExtrapToZ(Double_t Z);
   void ExtrapToStation(Int_t Station, AliMUONTrackParam *TrackParam);
index 05b354b04494005ca0763f54f16af9226b8ac7d1..6fcc08cea449c90f3642ba527c5b2290ccd4a773 100644 (file)
@@ -98,6 +98,7 @@
 #pragma link C++ class  AliMUONSt2GeometryBuilder+;
 #pragma link C++ class  AliMUONSlatGeometryBuilder+;
 #pragma link C++ class  AliMUONTriggerGeometryBuilder+;
+#pragma link C++ class  AliMUONRecoCheck+;
 
 #endif
 
diff --git a/MUON/MUONRecoCheck.C b/MUON/MUONRecoCheck.C
new file mode 100644 (file)
index 0000000..9b6867b
--- /dev/null
@@ -0,0 +1,232 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+// ROOT includes
+#include "TClonesArray.h"
+#include "TH1.h"
+#include "TParticle.h"
+#include "TFile.h"
+
+// STEER includes
+#include "AliRun.h"
+#include "AliHeader.h"
+#include "AliMC.h"
+#include "AliStack.h"
+#include "AliRunLoader.h"
+
+// MUON includes
+#include "AliMUON.h"
+#include "AliMUONConstants.h"
+#include "AliMUONTrack.h"
+#include "AliMUONRecoCheck.h"
+#include "AliMUONTrackParam.h"
+
+Int_t TrackCheck( Bool_t *compTrack);
+
+void MUONRecoCheck (Int_t nEvent = 1, char * filename="galice.root"){
+
+  // Utility macro 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 (TreeK) for the vertex parameters.     
+  
+  Int_t nTrackReco, nTrackRef;
+  AliMUONTrack *trackReco, *trackRef;
+  Bool_t *compTrack;
+  Bool_t compTrackOK[10];
+  Int_t nHitOK = 0;
+  Int_t testTrack = 0; 
+  Int_t iTrack = 0;
+  Int_t indexOK = 0;
+  Int_t trackID = 0;
+  Double_t sigma2Cut = 16;  // 4 sigmas cut, sigma2Cut = 4*4
+  AliMUONTrackParam *trackParam;
+  TClonesArray *trackParamAtHit;
+  Double_t x1,y1,z1,pX1,pY1,pZ1,p1;
+  Double_t x2,y2,z2,pX2,pY2,pZ2,p2;
+  TParticle* particle = new TParticle();
+
+  TClonesArray *trackRecoArray = NULL;
+  TClonesArray *trackRefArray = NULL;
+
+
+  // File for histograms and histogram booking
+  TFile *histoFile = new TFile("MUONRecoCheck.root", "RECREATE");
+
+  TH1F *hReconstructible = new TH1F("hReconstructible"," Nb of reconstructible tracks ",15,-0.5,14.5);
+  TH1F *hReco = new TH1F("hReco"," Nb of reconstructed tracks / evt",15,-0.5,14.5);
+  TH1F *hNHitComp = new TH1F("hNHitComp"," Nb of compatible hits / track ",15,-0.5,14.5);
+  TH1F *hTestTrack = new TH1F("hTestTrack"," Reconstruction requirement / track",15,-0.5,14.5);
+  TH1F *hTrackRefID = new TH1F("hTrackRefID"," track reference ID ",100,-0.5,99.5);
+  
+  TH1F *hResMomVertex = new TH1F("hMomVertex"," delta P vertex (GeV/c)",100,-10.,10);
+  TH1F *hResMomFirstHit = new TH1F("hMomFirstHit"," delta P first hit (GeV/c)",100,-10.,10);
+
+
+  
+  AliMUONRecoCheck rc(filename);
+  AliRunLoader *runLoader = rc.GetRunLoader();
+
+    
+  Int_t nevents = runLoader->GetNumberOfEvents();
+  
+  if (nevents < nEvent) nEvent = nevents;
+  
+  Int_t ievent;
+  Int_t nReconstructibleTracks = 0;
+  Int_t nReconstructedTracks = 0;
+  Int_t nReconstructibleTracksCheck = 0;
+
+  for (ievent=0; ievent<nEvent; ievent++) {
+    if (!(ievent%10)) printf(" **** event # %d  \n",ievent);
+    runLoader->GetEvent(ievent);
+    rc.ResetTracks();
+    rc.MakeTrackRef(); // make reconstructible tracks
+//     rc.PrintEvent();
+
+    
+    trackRecoArray = rc.GetTrackReco();
+    trackRefArray = rc.GetMuonTrackRef();
+    
+    nTrackRef = trackRefArray->GetEntriesFast();
+    nTrackReco = trackRecoArray->GetEntriesFast();
+
+    hReconstructible->Fill(rc.GetNumberOfReconstuctibleTracks());
+    hReco->Fill(rc.GetNumberOfRecoTracks());
+
+    nReconstructibleTracks += rc.GetNumberOfReconstuctibleTracks();
+    nReconstructedTracks += rc.GetNumberOfRecoTracks();
+
+    for (Int_t index1 = 0; index1 < nTrackRef; index1++) {  
+      trackRef = (AliMUONTrack *)trackRefArray->At(index1);
+
+      testTrack = 0;
+      indexOK = 0;
+      for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) 
+       compTrackOK[ch] = kFALSE;
+      for (Int_t index2 = 0; index2 < nTrackReco; index2++) {
+       trackReco = (AliMUONTrack *)trackRecoArray->At(index2);
+
+       // check if trackRef is compatible with trackReco
+       compTrack = trackRef->CompatibleTrack(trackReco,sigma2Cut);
+
+       iTrack = TrackCheck(compTrack);
+       
+       if (iTrack > testTrack) {
+         nHitOK = 0;
+         for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
+           if (compTrack[ch]) nHitOK++;
+           compTrackOK[ch] = compTrack[ch];
+         }
+         testTrack = iTrack;
+         indexOK = index2;
+       }
+      }
+
+      hTestTrack->Fill(testTrack);
+      trackID = trackRef->GetTrackID();
+      hTrackRefID->Fill(trackID);
+
+      if (testTrack == 4) {     // tracking requirements verified, track is found
+       nReconstructibleTracksCheck++;
+       hNHitComp->Fill(nHitOK);
+       particle = runLoader->GetHeader()->Stack()->Particle(trackID);
+//     printf(" trackID: %d , PDG code: %d \n",trackID,particle->GetPdgCode());
+       trackParam = trackRef->GetTrackParamAtVertex();
+       x1 = trackParam->GetNonBendingCoor();
+       y1 = trackParam->GetBendingCoor();
+       z1 = trackParam->GetZ();
+       pX1 = trackParam->Px();
+       pY1 = trackParam->Py();
+       pZ1 = trackParam->Pz();
+       p1  = trackParam->P();
+       
+//     printf(" Ref. track at vertex: x,y,z: %f %f %f px,py,pz,p: %f %f %f %f \n",x1,y1,z1,pX1,pY1,pZ1,p1);
+       
+       trackParam = ((AliMUONTrack *)trackRecoArray->At(indexOK))->GetTrackParamAtVertex();
+       x2 = trackParam->GetNonBendingCoor();
+       y2 = trackParam->GetBendingCoor();
+       z2 = trackParam->GetZ();
+       pX2 = trackParam->Px();
+       pY2 = trackParam->Py();
+       pZ2 = trackParam->Pz();
+       p2  = trackParam->P();
+//     printf(" Reconst. track at vertex: x,y,z: %f %f %f px,py,pz: %f %f %f %f \n",x2,y2,z2,pX2,pY2,pZ2,p2);
+       
+       hResMomVertex->Fill(p2-p1);
+
+       trackParamAtHit =  trackRef->GetTrackParamAtHit();
+       trackParam = (AliMUONTrackParam*) trackParamAtHit->First();
+       x1 = trackParam->GetNonBendingCoor();
+       y1 = trackParam->GetBendingCoor();
+       z1 = trackParam->GetZ();
+       pX1 = trackParam->Px();
+       pY1 = trackParam->Py();
+       pZ1 = trackParam->Pz();
+       p1  = trackParam->P();
+//     printf(" Ref. track at 1st hit: x,y,z: %f %f %f px,py,pz: %f %f %f \n",x1,y1,z1,pX1,pY1,pZ1);
+       trackParamAtHit =  ((AliMUONTrack *) trackRecoArray->At(indexOK))->GetTrackParamAtHit();
+       trackParam = (AliMUONTrackParam*) trackParamAtHit->First();
+       x2 = trackParam->GetNonBendingCoor();
+       y2 = trackParam->GetBendingCoor();
+       z2 = trackParam->GetZ();
+       pX2 = trackParam->Px();
+       pY2 = trackParam->Py();
+       pZ2 = trackParam->Pz();
+       p2  = trackParam->P();
+//     printf(" Reconst. track at 1st hit: x,y,z: %f %f %f px,py,pz: %f %f %f \n",x2,y2,z2,pX2,pY2,pZ2);
+
+       hResMomFirstHit->Fill(p2-p1);
+              
+      }
+    } // end loop track ref.
+
+  } // end loop on event  
+
+  printf(" nb of reconstructible tracks: %d \n", nReconstructibleTracks);
+  printf(" nb of reconstructed tracks: %d \n", nReconstructedTracks);
+  printf(" nb of reconstructible tracks which are reconstructed: %d \n", nReconstructibleTracksCheck);
+  
+  histoFile->Write();
+  histoFile->Close();
+}
+
+
+Int_t TrackCheck( Bool_t *compTrack)
+{
+  // Apply reconstruction requirements
+  // Return number of validated conditions 
+  // If all the tests are verified then TrackCheck = 4 (good track)
+  Int_t iTrack = 0;
+  Int_t hitsInLastStations = 0;
+  
+  // apply reconstruction requirements
+  if (compTrack[0] || compTrack[1]) iTrack++; // at least one hit in st. 0
+  if (compTrack[2] || compTrack[3]) iTrack++; // at least one hit in st. 1
+  if (compTrack[4] || compTrack[5]) iTrack++; // at least one hit in st. 2
+  for (Int_t ch = 6; ch < AliMUONConstants::NTrackingCh(); ch++) {
+    if (compTrack[ch]) hitsInLastStations++; 
+  }
+  if (hitsInLastStations > 2) iTrack++; // at least 3 hits in st. 3 & 4
+  
+  return iTrack;
+
+}
+
+
+
+
+
+
+
index d542d3b29cc71eb219067a26d2ee99f38c010d1e..6959e509ea4ed476a1ceb9dc24d213889acbf6c9 100644 (file)
@@ -179,6 +179,25 @@ MUONmassPlot("galice.root",0,99999);
 .q
 EOF
 
+
+============================================================
+ How to run MUONRecoCheck macro
+============================================================
+To check the muon reconstruction by comparing the reconstructed tracks
+with the reference tracks made of "AliTrackReference" for the hits and
+kinematic informations (TParticle) for the vertex.
+This macro can be used to check the track reconstruction e.g. efficiency,
+momentum resolution ... but also to make physics analysis whenever
+track identification is needed.   
+
+To compile MUONRecoCheck.C
+.includepath $ALICE_ROOT/STEER
+.includepath $ALICE_ROOT/MUON
+.L $ALICE_ROOT/MUON/MUONRecoCheck.C+
+
+// To run MUONRecoCheck
+MUONRecoCheck(nEvent,"galice.root"); // nEvent = nb of events
+
 ===========================================================
  Still working ..............
 ===========================================================
index 89417ef94048e63b85348c65b3d7b5d2807a01f9..f5c257a2ec16c7f6a0d5d8748117e4bf20ff63aa 100644 (file)
@@ -42,7 +42,8 @@ SRCS         = AliMUONChamber.cxx AliMUONChamberTrigger.cxx \
                AliMUONTrackK.cxx AliMUONClusterFinderAZ.cxx AliMUONPixel.cxx \
                AliMUONLoader.cxx AliMUONData.cxx  AliMUONDataInterface.cxx \
                AliMUONDDLTrigger.cxx AliMUONSubEventTrigger.cxx AliMUONDDLTracker.cxx \
-               AliMUONRawData.cxx AliMUONSubEventTracker.cxx AliMUONRawStream.cxx
+               AliMUONRawData.cxx AliMUONSubEventTracker.cxx AliMUONRawStream.cxx \
+               AliMUONRecoCheck.cxx 
 
 
 SRCS     += AliMUONSt1Segmentation.cxx AliMUONSt1Response.cxx \