]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - MUON/AliMUONTrackReconstructor.cxx
Additional protection in case of negative indexes. More investigation is needed
[u/mrichter/AliRoot.git] / MUON / AliMUONTrackReconstructor.cxx
index c19000feda7abf334d4fe81c37b298a74a6c0533..2b894b52203ab772aa266a972ca24225ea52c2cc 100644 (file)
@@ -13,6 +13,8 @@
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
+/* $Id$ */
+
 ////////////////////////////////////
 //
 // MUON track reconstructor in ALICE (class renamed from AliMUONEventReconstructor)
@@ -34,7 +36,7 @@
 #include <Riostream.h>
 #include <TDirectory.h>
 #include <TFile.h>
-#include <TMatrixD.h> //AZ
+#include <TMatrixD.h>
 
 #include "AliMUONTrackReconstructor.h"
 #include "AliMUON.h"
@@ -52,7 +54,7 @@
 #include "AliRun.h" // for gAlice
 #include "AliRunLoader.h"
 #include "AliLoader.h"
-#include "AliMUONTrackK.h" //AZ
+#include "AliMUONTrackK.h" 
 #include "AliMC.h"
 #include "AliLog.h"
 #include "AliTrackReference.h"
@@ -79,8 +81,8 @@ const Double_t AliMUONTrackReconstructor::fgkDefaultEfficiency = 0.95;
 
 ClassImp(AliMUONTrackReconstructor) // Class implementation in ROOT context
 
-  //__________________________________________________________________________
-AliMUONTrackReconstructor::AliMUONTrackReconstructor(AliLoader* loader)
+//__________________________________________________________________________
+AliMUONTrackReconstructor::AliMUONTrackReconstructor(AliLoader* loader, AliMUONData* data)
   : TObject()
 {
   // Constructor for class AliMUONTrackReconstructor
@@ -102,7 +104,6 @@ AliMUONTrackReconstructor::AliMUONTrackReconstructor(AliLoader* loader)
   fNRecTracks = 0; // really needed or GetEntriesFast sufficient ????
   // Memory allocation for the TClonesArray of hits on reconstructed tracks
   // Is 100 the right size ????
-  //AZ fRecTrackHitsPtr = new TClonesArray("AliMUONTrack", 100);
   fRecTrackHitsPtr = new TClonesArray("AliMUONTrackHit", 100);
   fNRecTrackHits = 0; // really needed or GetEntriesFast sufficient ????
 
@@ -132,7 +133,8 @@ AliMUONTrackReconstructor::AliMUONTrackReconstructor(AliLoader* loader)
   fLoader = loader;
 
   // initialize container
-  fMUONData  = new AliMUONData(fLoader,"MUON","MUON");
+  // fMUONData  = new AliMUONData(fLoader,"MUON","MUON");
+  fMUONData  = data;
 
   return;
 }
@@ -347,6 +349,10 @@ void AliMUONTrackReconstructor::EventReconstruct(void)
 {
   // To reconstruct one event
   AliDebug(1,"Enter EventReconstruct");
+  ResetTracks(); //AZ
+  ResetTrackHits(); //AZ 
+  ResetSegments(); //AZ
+  ResetHitsForRec(); //AZ
   MakeEventToBeReconstructed();
   MakeSegments();
   MakeTracks();
@@ -442,7 +448,7 @@ void AliMUONTrackReconstructor::MakeEventToBeReconstructed(void)
   AliRunLoader *runLoader = fLoader->GetRunLoader();
 
   AliDebug(1,"Enter MakeEventToBeReconstructed");
-  ResetHitsForRec();
+  //AZ ResetHitsForRec();
   if (fRecTrackRefHits == 1) {
     // Reconstruction from track ref. hits
     // Back to the signal file
@@ -477,7 +483,7 @@ void AliMUONTrackReconstructor::MakeEventToBeReconstructed(void)
     // TreeR assumed to be be "prepared" in calling function
     // by "MUON->GetTreeR(nev)" ????
     TTree *treeR = fLoader->TreeR();
-    fMUONData->SetTreeAddress("RC");
+    //AZ? fMUONData->SetTreeAddress("RC");
     AddHitsForRecFromRawClusters(treeR);
     // No sorting: it is done automatically in the previous function
   }
@@ -711,12 +717,15 @@ void AliMUONTrackReconstructor::AddHitsForRecFromRawClusters(TTree* TR)
   TClonesArray *rawclusters;
   AliDebug(1,"Enter AddHitsForRecFromRawClusters");
 
-  nTRentries = Int_t(TR->GetEntries());
-  if (nTRentries != 1) {
-    AliError(Form("nTRentries = %d not equal to 1 ",nTRentries));
-    exit(0);
+  if (fTrackMethod != 3) { //AZ
+    fMUONData->SetTreeAddress("RC"); //AZ
+    nTRentries = Int_t(TR->GetEntries());
+    if (nTRentries != 1) {
+      AliError(Form("nTRentries = %d not equal to 1 ",nTRentries));
+      exit(0);
+    }
+    fLoader->TreeR()->GetEvent(0); // only one entry  
   }
-  fLoader->TreeR()->GetEvent(0); // only one entry  
 
   // Loop over tracking chambers
   for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
@@ -739,8 +748,10 @@ void AliMUONTrackReconstructor::AddHitsForRecFromRawClusters(TTree* TR)
       (fNHitsForRecPerChamber[ch])++;
       // more information into HitForRec
       //  resolution: info should be already in raw cluster and taken from it ????
-      hitForRec->SetBendingReso2(fBendingResolution * fBendingResolution);
-      hitForRec->SetNonBendingReso2(fNonBendingResolution * fNonBendingResolution);
+      //hitForRec->SetBendingReso2(fBendingResolution * fBendingResolution);
+      //hitForRec->SetNonBendingReso2(fNonBendingResolution * fNonBendingResolution);
+      hitForRec->SetBendingReso2(clus->GetErrY() * clus->GetErrY());
+      hitForRec->SetNonBendingReso2(clus->GetErrX() * clus->GetErrX());
       //  original raw cluster
       hitForRec->SetChamberNumber(ch);
       hitForRec->SetHitNumber(iclus);
@@ -754,7 +765,7 @@ void AliMUONTrackReconstructor::AddHitsForRecFromRawClusters(TTree* TR)
        hitForRec->Dump();}
     } // end of cluster loop
   } // end of chamber loop
-  SortHitsForRecWithIncreasingChamber(); //AZ 
+  SortHitsForRecWithIncreasingChamber(); 
   return;
 }
 
@@ -764,11 +775,10 @@ void AliMUONTrackReconstructor::MakeSegments(void)
   // To make the list of segments in all stations,
   // from the list of hits to be reconstructed
   AliDebug(1,"Enter MakeSegments");
-  ResetSegments();
+  //AZ ResetSegments();
   // Loop over stations
-  Int_t nb = (fTrackMethod == 2) ? 3 : 0; //AZ
-  //AZ for (Int_t st = 0; st < fgkMaxMuonTrackingStations; st++)
-  for (Int_t st = nb; st < AliMUONConstants::NTrackingCh()/2; st++) //AZ
+  Int_t nb = (fTrackMethod != 1) ? 3 : 0; //AZ
+  for (Int_t st = nb; st < AliMUONConstants::NTrackingCh()/2; st++) 
     MakeSegmentsPerStation(st); 
   if (AliLog::GetGlobalDebugLevel() > 1) {
     cout << "end of MakeSegments" << endl;
@@ -901,9 +911,9 @@ void AliMUONTrackReconstructor::MakeTracks(void)
   // from the list of segments and points in all stations
   AliDebug(1,"Enter MakeTracks");
   // The order may be important for the following Reset's
-  ResetTracks();
-  ResetTrackHits();
-  if (fTrackMethod == 2) { //AZ - Kalman filter
+  //AZ ResetTracks();
+  //AZ ResetTrackHits();
+  if (fTrackMethod != 1) { //AZ - Kalman filter
     MakeTrackCandidatesK();
     if (fRecTracksPtr->GetEntriesFast() == 0) return;
     // Follow tracks in stations(1..) 3, 2 and 1
@@ -1188,7 +1198,7 @@ void AliMUONTrackReconstructor::FollowTracks(void)
   Int_t ch = -1, chInStation, chBestHit = -1, iHit, iSegment, station, trackIndex; 
   Double_t bestChi2, chi2, dZ1, dZ2, dZ3, maxSigma2Distance, mcsFactor;
   Double_t bendingMomentum, chi2Norm = 0.;
-  AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
+
   // local maxSigma2Distance, for easy increase in testing
   maxSigma2Distance = fMaxSigma2Distance;
   AliDebug(2,"Enter FollowTracks");
@@ -1225,14 +1235,14 @@ void AliMUONTrackReconstructor::FollowTracks(void)
       mcsFactor = 0.0136 * trackParam1->GetInverseBendingMomentum();
       mcsFactor        = fChamberThicknessInX0 * mcsFactor * mcsFactor;
       // Z difference from previous station
-      dZ1 = (&(pMUON->Chamber(2 * station)))->Z() -
-       (&(pMUON->Chamber(2 * station + 2)))->Z();
+      dZ1 = AliMUONConstants::DefaultChamberZ(2 * station) -
+           AliMUONConstants::DefaultChamberZ(2 * station + 2);
       // Z difference between the two previous stations
-      dZ2 = (&(pMUON->Chamber(2 * station + 2)))->Z() -
-       (&(pMUON->Chamber(2 * station + 4)))->Z();
+      dZ2 = AliMUONConstants::DefaultChamberZ(2 * station + 2) -
+           AliMUONConstants::DefaultChamberZ(2 * station + 4);
       // Z difference between the two chambers in the previous station
-      dZ3 = (&(pMUON->Chamber(2 * station)))->Z() -
-       (&(pMUON->Chamber(2 * station + 1)))->Z();
+      dZ3 = AliMUONConstants::DefaultChamberZ(2 * station) -
+           AliMUONConstants::DefaultChamberZ(2 * station + 1);
       extrapSegment->SetBendingCoorReso2(fBendingResolution * fBendingResolution);
       extrapSegment->
        SetNonBendingCoorReso2(fNonBendingResolution * fNonBendingResolution);
@@ -1599,21 +1609,15 @@ void AliMUONTrackReconstructor::MakeTrackCandidatesK(void)
   AliMUONTrackK *trackK;
 
   AliDebug(1,"Enter MakeTrackCandidatesK");
-  // Reset the TClonesArray of reconstructed tracks
-  if (fRecTracksPtr) fRecTracksPtr->Delete();
-  // Delete in order that the Track destructors are called,
-  // hence the space for the TClonesArray of pointers to TrackHit's is freed
-  fNRecTracks = 0;
 
-  AliMUONTrackK a(this, fHitsForRecPtr); // bad idea ???
+  AliMUONTrackK a(this, fHitsForRecPtr);
   // Loop over stations(1...) 5 and 4
   for (istat=4; istat>=3; istat--) {
     // Loop over segments in the station
     for (iseg=0; iseg<fNSegments[istat]; iseg++) {
       // Transform segments to tracks and evaluate covariance matrix
       segment = (AliMUONSegment*) ((*fSegmentsPtr[istat])[iseg]);
-      trackK = new ((*fRecTracksPtr)[fNRecTracks]) AliMUONTrackK(segment);
-      fNRecTracks++;
+      trackK = new ((*fRecTracksPtr)[fNRecTracks++]) AliMUONTrackK(segment);
     } // for (iseg=0;...)
   } // for (istat=4;...)
   return;
@@ -1632,8 +1636,6 @@ void AliMUONTrackReconstructor::FollowTracksK(void)
   TClonesArray *rawclusters;
   clus = 0; rawclusters = 0;
 
-  //AZ(z->-z) zDipole1 = GetSimpleBPosition() - GetSimpleBLength()/2;
-  //AZ(z->-z) zDipole2 = zDipole1 + GetSimpleBLength();
   zDipole1 = GetSimpleBPosition() + GetSimpleBLength()/2;
   zDipole2 = zDipole1 - GetSimpleBLength();
 
@@ -1654,8 +1656,8 @@ void AliMUONTrackReconstructor::FollowTracksK(void)
        rawclusters = fMUONData->RawClusters(hit->GetChamberNumber());
        clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->
                                                  GetHitNumber());
-       printf("%3d", clus->GetTrack(1)-1);
-       if (clus->GetTrack(2) != 0) printf("%3d \n", clus->GetTrack(2)-1);
+       printf(" %d", clus->GetTrack(1));
+       if (clus->GetTrack(2) != -1) printf(" %d \n", clus->GetTrack(2));
        else printf("\n");
       } // if (fRecTrackRefHits)
     }
@@ -1684,14 +1686,12 @@ void AliMUONTrackReconstructor::FollowTracksK(void)
 
     ok = kTRUE;
     if (trackK->GetRecover() == 0) hit = (AliMUONHitForRec*) 
-                                   trackK->GetHitOnTrack()->Last(); // last hit
-    //else hit = (AliMUONHitForRec*) (*trackK->GetHitOnTrack())[1]; // 2'nd hit
+                                   trackK->GetTrackHits()->Last(); // last hit
     else hit = trackK->GetHitLastOk(); // hit where track stopped
     if (hit) ichamBeg = hit->GetChamberNumber();
     ichamEnd = 0;
     // Check propagation direction
-    if (hit == NULL) ichamBeg = ichamEnd;
-    //AZ(z->-z) else if (trackK->GetTrackDir() > 0) {
+    if (!hit) { ichamBeg = ichamEnd; AliFatal(" ??? "); }
     else if (trackK->GetTrackDir() < 0) {
       ichamEnd = 9; // forward propagation
       ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2);
@@ -1700,7 +1700,6 @@ void AliMUONTrackReconstructor::FollowTracksK(void)
         ichamEnd = 6; // backward propagation
        // Change weight matrix and zero fChi2 for backpropagation
         trackK->StartBack();
-       //AZ trackK->SetTrackDir(-1);
        trackK->SetTrackDir(1);
         ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kTRUE,zDipole1,zDipole2);
         ichamBeg = ichamEnd;
@@ -1719,7 +1718,6 @@ void AliMUONTrackReconstructor::FollowTracksK(void)
     }
 
     if (ok) {
-      //AZ trackK->SetTrackDir(-1);
       trackK->SetTrackDir(1);
       trackK->SetBPFlag(kFALSE);
       ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2);
@@ -1737,7 +1735,7 @@ void AliMUONTrackReconstructor::FollowTracksK(void)
     chamBits = 0;
     Double_t chi2max = 0;
     for (Int_t i=0; i<trackK->GetNTrackHits(); i++) {
-      hit = (AliMUONHitForRec*) (*trackK->GetHitOnTrack())[i];
+      hit = (AliMUONHitForRec*) (*trackK->GetTrackHits())[i];
       chamBits |= BIT(hit->GetChamberNumber());
       if (trackK->GetChi2PerPoint(i) > chi2max) chi2max = trackK->GetChi2PerPoint(i);
     }
@@ -1769,8 +1767,8 @@ Bool_t AliMUONTrackReconstructor::CheckCandidateK(Int_t icand, Int_t nSeeds) con
   AliMUONHitForRec *hit1, *hit2, *hit;
 
   track1 = (AliMUONTrackK*) ((*fRecTracksPtr)[icand]);
-  hit1 = (AliMUONHitForRec*) (*track1->GetHitOnTrack())[0]; // 1'st hit
-  hit2 = (AliMUONHitForRec*) (*track1->GetHitOnTrack())[1]; // 2'nd hit
+  hit1 = (AliMUONHitForRec*) (*track1->GetTrackHits())[0]; // 1'st hit
+  hit2 = (AliMUONHitForRec*) (*track1->GetTrackHits())[1]; // 2'nd hit
 
   for (Int_t i=0; i<icand; i++) {
     track2 = (AliMUONTrackK*) ((*fRecTracksPtr)[i]);
@@ -1782,7 +1780,7 @@ Bool_t AliMUONTrackReconstructor::CheckCandidateK(Int_t icand, Int_t nSeeds) con
     } else {
       Int_t nSame = 0;
       for (Int_t j=0; j<track2->GetNTrackHits(); j++) {
-        hit = (AliMUONHitForRec*) (*track2->GetHitOnTrack())[j];
+        hit = (AliMUONHitForRec*) (*track2->GetTrackHits())[j];
         if (hit == hit1 || hit == hit2) {
           nSame++;
           if (nSame == 2) return kFALSE;
@@ -1848,11 +1846,13 @@ void AliMUONTrackReconstructor::SetTrackMethod(Int_t iTrackMethod)
 {
   // Set track method and recreate track container if necessary
   
-  fTrackMethod = iTrackMethod == 2 ? 2 : 1;
-  if (fTrackMethod == 2) {
+  fTrackMethod = TMath::Min (iTrackMethod, 3);
+  fTrackMethod = TMath::Max (fTrackMethod, 1);
+  if (fTrackMethod != 1) {
     if (fRecTracksPtr) delete fRecTracksPtr;
     fRecTracksPtr = new TClonesArray("AliMUONTrackK", 10);
-    cout << " *** Tracking with the Kalman filter *** " << endl;
-  } else cout << " *** Traditional tracking *** " << endl;
+    if (fTrackMethod == 2) cout << " *** Tracking with the Kalman filter *** " << endl;
+    else cout << " *** Combined cluster / track finder ***" << endl;
+  } else cout << " *** Original tracking *** " << endl;
 
 }