]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - MUON/AliMUONEventReconstructor.cxx
Code revision:
[u/mrichter/AliRoot.git] / MUON / AliMUONEventReconstructor.cxx
index 76121eddbc0fefc4e634115224f1fc74a9f96ff5..5a41a46fe8fc90a01e75c56cbb8be883fb5a95f0 100644 (file)
 //
 ////////////////////////////////////
 
-#include <Riostream.h> // for cout
-#include <stdlib.h> // for exit()
-
-#include <TTree.h>
+#include <Riostream.h>
+#include <TDirectory.h>
+#include <TFile.h>
+#include <TMatrixD.h> //AZ
 
-#include "AliMUON.h"
-//#include "AliMUONChamber.h"
 #include "AliMUONEventReconstructor.h"
+#include "AliMUON.h"
+#include "AliMUONHit.h"
 #include "AliMUONHitForRec.h"
 #include "AliMUONTriggerTrack.h"
-//#include "AliMUONTriggerConstants.h"
 #include "AliMUONTriggerCircuit.h"
 #include "AliMUONRawCluster.h"
 #include "AliMUONLocalTrigger.h"
 #include "AliMUONTrackHit.h"
 #include "AliMagF.h"
 #include "AliRun.h" // for gAlice
-#include "AliConfig.h"
 #include "AliRunLoader.h"
 #include "AliLoader.h"
 #include "AliMUONTrackK.h" //AZ
-#include <TMatrixD.h> //AZ
 #include "AliMC.h"
 
 //************* Defaults parameters for reconstruction
-static const Double_t kDefaultMinBendingMomentum = 3.0;
-static const Double_t kDefaultMaxBendingMomentum = 500.0;
-static const Double_t kDefaultMaxChi2 = 100.0;
-static const Double_t kDefaultMaxSigma2Distance = 16.0;
-static const Double_t kDefaultBendingResolution = 0.01;
-static const Double_t kDefaultNonBendingResolution = 0.144;
-static const Double_t kDefaultChamberThicknessInX0 = 0.03;
+const Double_t AliMUONEventReconstructor::fgkDefaultMinBendingMomentum = 3.0;
+const Double_t AliMUONEventReconstructor::fgkDefaultMaxBendingMomentum = 500.0;
+const Double_t AliMUONEventReconstructor::fgkDefaultMaxChi2 = 100.0;
+const Double_t AliMUONEventReconstructor::fgkDefaultMaxSigma2Distance = 16.0;
+const Double_t AliMUONEventReconstructor::fgkDefaultBendingResolution = 0.01;
+const Double_t AliMUONEventReconstructor::fgkDefaultNonBendingResolution = 0.144;
+const Double_t AliMUONEventReconstructor::fgkDefaultChamberThicknessInX0 = 0.03;
 // Simple magnetic field:
 // Value taken from macro MUONtracking.C: 0.7 T, hence 7 kG
 // Length and Position from reco_muon.F, with opposite sign:
 // Length = ZMAGEND-ZCOIL
 // Position = (ZMAGEND+ZCOIL)/2
 // to be ajusted differently from real magnetic field ????
-static const Double_t kDefaultSimpleBValue = 7.0;
-static const Double_t kDefaultSimpleBLength = 428.0;
-static const Double_t kDefaultSimpleBPosition = 1019.0;
-static const Int_t kDefaultRecGeantHits = 0;
-static const Double_t kDefaultEfficiency = 0.95;
+const Double_t AliMUONEventReconstructor::fgkDefaultSimpleBValue = 7.0;
+const Double_t AliMUONEventReconstructor::fgkDefaultSimpleBLength = 428.0;
+const Double_t AliMUONEventReconstructor::fgkDefaultSimpleBPosition = 1019.0;
+const Int_t    AliMUONEventReconstructor::fgkDefaultRecGeantHits = 0;
+const Double_t AliMUONEventReconstructor::fgkDefaultEfficiency = 0.95;
 
-static const Int_t kDefaultPrintLevel = -1;
+const Int_t    AliMUONEventReconstructor::fgkDefaultPrintLevel = -1;
 
 ClassImp(AliMUONEventReconstructor) // Class implementation in ROOT context
 
   //__________________________________________________________________________
 AliMUONEventReconstructor::AliMUONEventReconstructor(AliLoader* loader)
+  : TObject()
 {
   // Constructor for class AliMUONEventReconstructor
   SetReconstructionParametersToDefaults();
@@ -96,7 +94,7 @@ AliMUONEventReconstructor::AliMUONEventReconstructor(AliLoader* loader)
   fNHitsForRec = 0; // really needed or GetEntriesFast sufficient ????
   // Memory allocation for the TClonesArray's of segments in stations
   // Is 2000 the right size ????
-  for (Int_t st = 0; st < kMaxMuonTrackingStations; st++) {
+  for (Int_t st = 0; st < fgkMaxMuonTrackingStations; st++) {
     fSegmentsPtr[st] = new TClonesArray("AliMUONSegment", 2000);
     fNSegments[st] = 0; // really needed or GetEntriesFast sufficient ????
   }
@@ -155,15 +153,24 @@ AliMUONEventReconstructor::AliMUONEventReconstructor(AliLoader* loader)
   return;
 }
   //__________________________________________________________________________
-AliMUONEventReconstructor::AliMUONEventReconstructor (const AliMUONEventReconstructor& Reconstructor):TObject(Reconstructor)
+AliMUONEventReconstructor::AliMUONEventReconstructor (const AliMUONEventReconstructor& rhs)
+  : TObject(rhs)
 {
-  // Dummy copy constructor
+// Protected copy constructor
+
+  Fatal("AliMUONEventReconstructor", "Not implemented.");
 }
 
-AliMUONEventReconstructor & AliMUONEventReconstructor::operator=(const AliMUONEventReconstructor& /*Reconstructor*/)
+AliMUONEventReconstructor & 
+AliMUONEventReconstructor::operator=(const AliMUONEventReconstructor& rhs)
 {
-  // Dummy assignment operator
-    return *this;
+// Protected assignement operator
+
+  if (this == &rhs) return *this;
+
+  Fatal("operator=", "Not implemented.");
+    
+  return *this;  
 }
 
   //__________________________________________________________________________
@@ -177,7 +184,7 @@ AliMUONEventReconstructor::~AliMUONEventReconstructor(void)
 //  if (fEventTree) delete fEventTree;
   if (fRecoEvent) delete fRecoEvent;
   delete fHitsForRecPtr; // Correct destruction of everything ???? or delete [] ????
-  for (Int_t st = 0; st < kMaxMuonTrackingStations; st++)
+  for (Int_t st = 0; st < fgkMaxMuonTrackingStations; st++)
     delete fSegmentsPtr[st]; // Correct destruction of everything ????
   return;
 }
@@ -187,10 +194,10 @@ void AliMUONEventReconstructor::SetReconstructionParametersToDefaults(void)
 {
   // Set reconstruction parameters to default values
   // Would be much more convenient with a structure (or class) ????
-  fMinBendingMomentum = kDefaultMinBendingMomentum;
-  fMaxBendingMomentum = kDefaultMaxBendingMomentum;
-  fMaxChi2 = kDefaultMaxChi2;
-  fMaxSigma2Distance = kDefaultMaxSigma2Distance;
+  fMinBendingMomentum = fgkDefaultMinBendingMomentum;
+  fMaxBendingMomentum = fgkDefaultMaxBendingMomentum;
+  fMaxChi2 = fgkDefaultMaxChi2;
+  fMaxSigma2Distance = fgkDefaultMaxSigma2Distance;
 
   AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON");
   // ******** Parameters for making HitsForRec
@@ -198,7 +205,7 @@ void AliMUONEventReconstructor::SetReconstructionParametersToDefaults(void)
   // like in TRACKF_STAT:
   // 2 degrees for stations 1 and 2, or ch(0...) from 0 to 3;
   // 30 cm for stations 3 to 5, or ch(0...) from 4 to 9
-  for (Int_t ch = 0; ch < kMaxMuonTrackingChambers; ch++) {
+  for (Int_t ch = 0; ch < fgkMaxMuonTrackingChambers; ch++) {
     if (ch < 4) fRMin[ch] = TMath::Abs((&(pMUON->Chamber(ch)))->Z()) *
                  2.0 * TMath::Pi() / 180.0;
     else fRMin[ch] = 30.0;
@@ -213,7 +220,7 @@ void AliMUONEventReconstructor::SetReconstructionParametersToDefaults(void)
   // Maximum distance in non bending plane
   // 5 * 0.22 just to remember the way it was made in TRACKF_STAT
   // SIGCUT*DYMAX(IZ)
-  for (Int_t st = 0; st < kMaxMuonTrackingStations; st++)
+  for (Int_t st = 0; st < fgkMaxMuonTrackingStations; st++)
     fSegmentMaxDistNonBending[st] = 5. * 0.22;
   // Maximum distance in bending plane:
   // values from TRACKF_STAT, corresponding to (J psi 20cm),
@@ -229,15 +236,15 @@ void AliMUONEventReconstructor::SetReconstructionParametersToDefaults(void)
   fSegmentMaxDistBending[4] =  TMath::Abs( 6.0 *
     ((&(pMUON->Chamber(9)))->Z() - (&(pMUON->Chamber(8)))->Z()) / 20.0);
   
-  fBendingResolution = kDefaultBendingResolution;
-  fNonBendingResolution = kDefaultNonBendingResolution;
-  fChamberThicknessInX0 = kDefaultChamberThicknessInX0;
-  fSimpleBValue = kDefaultSimpleBValue;
-  fSimpleBLength = kDefaultSimpleBLength;
-  fSimpleBPosition = kDefaultSimpleBPosition;
-  fRecGeantHits = kDefaultRecGeantHits;
-  fEfficiency = kDefaultEfficiency;
-  fPrintLevel = kDefaultPrintLevel;
+  fBendingResolution = fgkDefaultBendingResolution;
+  fNonBendingResolution = fgkDefaultNonBendingResolution;
+  fChamberThicknessInX0 = fgkDefaultChamberThicknessInX0;
+  fSimpleBValue = fgkDefaultSimpleBValue;
+  fSimpleBLength = fgkDefaultSimpleBLength;
+  fSimpleBPosition = fgkDefaultSimpleBPosition;
+  fRecGeantHits = fgkDefaultRecGeantHits;
+  fEfficiency = fgkDefaultEfficiency;
+  fPrintLevel = fgkDefaultPrintLevel;
   return;
 }
 
@@ -400,7 +407,7 @@ void AliMUONEventReconstructor::ResetHitsForRec(void)
   // and the index of the first HitForRec per chamber
   if (fHitsForRecPtr) fHitsForRecPtr->Clear();
   fNHitsForRec = 0;
-  for (Int_t ch = 0; ch < kMaxMuonTrackingChambers; ch++)
+  for (Int_t ch = 0; ch < fgkMaxMuonTrackingChambers; ch++)
     fNHitsForRecPerChamber[ch] = fIndexOfFirstHitForRecPerChamber[ch] = 0;
   return;
 }
@@ -410,7 +417,7 @@ void AliMUONEventReconstructor::ResetSegments(void)
 {
   // To reset the TClonesArray of segments and the number of Segments
   // for all stations
-  for (Int_t st = 0; st < kMaxMuonTrackingStations; st++) {
+  for (Int_t st = 0; st < fgkMaxMuonTrackingStations; st++) {
     if (fSegmentsPtr[st]) fSegmentsPtr[st]->Clear();
     fNSegments[st] = 0;
   }
@@ -514,7 +521,7 @@ void AliMUONEventReconstructor::MakeEventToBeReconstructed(void)
   if (fPrintLevel >= 10) {
     cout << "end of MakeEventToBeReconstructed" << endl;
     cout << "NHitsForRec: " << fNHitsForRec << endl;
-    for (Int_t ch = 0; ch < kMaxMuonTrackingChambers; ch++) {
+    for (Int_t ch = 0; ch < fgkMaxMuonTrackingChambers; ch++) {
       cout << "chamber(0...): " << ch
           << "  NHitsForRec: " << fNHitsForRecPerChamber[ch]
           << "  index(first HitForRec): " << fIndexOfFirstHitForRecPerChamber[ch]
@@ -622,7 +629,7 @@ AliMUONHitForRec* AliMUONEventReconstructor::NewHitForRecFromGEANT(AliMUONHit* H
   Double_t bendCoor, nonBendCoor, radius;
   Int_t chamber = Hit->Chamber() - 1; // chamber(0...)
   // only in tracking chambers (fChamber starts at 1)
-  if (chamber >= kMaxMuonTrackingChambers) return NULL;
+  if (chamber >= fgkMaxMuonTrackingChambers) return NULL;
   // only if hit is efficient (keep track for checking ????)
   if (gRandom->Rndm() > fEfficiency) return NULL;
   // only if radius between RMin and RMax
@@ -665,7 +672,7 @@ void AliMUONEventReconstructor::SortHitsForRecWithIncreasingChamber()
   // Update the information for HitsForRec per chamber too.
   Int_t ch, nhits, prevch;
   fHitsForRecPtr->Sort();
-  for (ch = 0; ch < kMaxMuonTrackingChambers; ch++) {
+  for (ch = 0; ch < fgkMaxMuonTrackingChambers; ch++) {
     fNHitsForRecPerChamber[ch] = 0;
     fIndexOfFirstHitForRecPerChamber[ch] = 0;
   }
@@ -700,7 +707,7 @@ void AliMUONEventReconstructor::SortHitsForRecWithIncreasingChamber()
 //   AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
 //   // Security on MUON ????
 //   // Loop over tracking chambers
-//   for (Int_t ch = 0; ch < kMaxMuonTrackingChambers; ch++) {
+//   for (Int_t ch = 0; ch < fgkMaxMuonTrackingChambers; ch++) {
 //     // number of HitsForRec to 0 for the chamber
 //     fNHitsForRecPerChamber[ch] = 0;
 //     // index of first HitForRec for the chamber
@@ -780,7 +787,7 @@ void AliMUONEventReconstructor::AddHitsForRecFromRawClusters(TTree* TR)
   fLoader->TreeR()->GetEvent(0); // only one entry  
 
   // Loop over tracking chambers
-  for (Int_t ch = 0; ch < kMaxMuonTrackingChambers; ch++) {
+  for (Int_t ch = 0; ch < fgkMaxMuonTrackingChambers; ch++) {
     // number of HitsForRec to 0 for the chamber
     fNHitsForRecPerChamber[ch] = 0;
     // index of first HitForRec for the chamber
@@ -828,12 +835,12 @@ void AliMUONEventReconstructor::MakeSegments(void)
   ResetSegments();
   // Loop over stations
   Int_t nb = (fTrackMethod == 2) ? 3 : 0; //AZ
-  //AZ for (Int_t st = 0; st < kMaxMuonTrackingStations; st++)
-  for (Int_t st = nb; st < kMaxMuonTrackingStations; st++) //AZ
+  //AZ for (Int_t st = 0; st < fgkMaxMuonTrackingStations; st++)
+  for (Int_t st = nb; st < fgkMaxMuonTrackingStations; st++) //AZ
     MakeSegmentsPerStation(st); 
   if (fPrintLevel >= 10) {
     cout << "end of MakeSegments" << endl;
-    for (Int_t st = 0; st < kMaxMuonTrackingStations; st++) {
+    for (Int_t st = 0; st < fgkMaxMuonTrackingStations; st++) {
       cout << "station(0...): " << st
           << "  Segments: " << fNSegments[st]
           << endl;
@@ -1022,14 +1029,14 @@ Bool_t AliMUONEventReconstructor::MakeTriggerTracks(void)
     AliMUONTriggerCircuit *circuit;
     AliMUONTriggerTrack *recTriggerTrack = 0;
 
-    TTree* TR = fLoader->TreeR();
+    TTree* treeR = fLoader->TreeR();
     
     // Loading MUON subsystem
     AliMUON * pMUON = (AliMUON *) gAlice->GetDetector("MUON");
     
-    nTRentries = Int_t(TR->GetEntries());
+    nTRentries = Int_t(treeR->GetEntries());
      
-    TR->GetEvent(0); // only one entry  
+    treeR->GetEvent(0); // only one entry  
 
     if (!(fMUONData->IsTriggerBranchesInTree())) {
       cout << "Warning in AliMUONEventReconstructor::MakeTriggerTracks"
@@ -1733,7 +1740,7 @@ void AliMUONEventReconstructor::MakeTrackCandidatesK(void)
 void AliMUONEventReconstructor::FollowTracksK(void)
 {
   // Follow tracks using Kalman filter
-  Bool_t Ok;
+  Bool_t ok;
   Int_t icand, ichamBeg, ichamEnd, chamBits;
   Double_t zDipole1, zDipole2;
   AliMUONTrackK *trackK;
@@ -1787,14 +1794,14 @@ void AliMUONEventReconstructor::FollowTracksK(void)
 
     // Discard candidate which will produce the double track
     if (icand > 0) {
-      Ok = CheckCandidateK(icand,nSeeds);
-      if (!Ok) {
+      ok = CheckCandidateK(icand,nSeeds);
+      if (!ok) {
         //trackK->SetRecover(-1); // mark candidate to be removed
         //continue;
       }
     }
 
-    Ok = kTRUE;
+    ok = kTRUE;
     if (trackK->GetRecover() == 0) hit = (AliMUONHitForRec*) 
                                    trackK->GetHitOnTrack()->Last(); // last hit
     else hit = (AliMUONHitForRec*) (*trackK->GetHitOnTrack())[1]; // 2'nd hit
@@ -1803,13 +1810,13 @@ void AliMUONEventReconstructor::FollowTracksK(void)
     // Check propagation direction
     if (trackK->GetTrackDir() > 0) {
       ichamEnd = 9; // forward propagation
-      Ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2);
-      if (Ok) {
+      ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2);
+      if (ok) {
         ichamBeg = ichamEnd;
         ichamEnd = 6; // backward propagation
        // Change weight matrix and zero fChi2 for backpropagation
         trackK->StartBack();
-        Ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kTRUE,zDipole1,zDipole2);
+        ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kTRUE,zDipole1,zDipole2);
         ichamBeg = ichamEnd;
         ichamEnd = 0;
       }
@@ -1819,18 +1826,18 @@ void AliMUONEventReconstructor::FollowTracksK(void)
         ichamEnd = 6; // backward propagation
        // Change weight matrix and zero fChi2 for backpropagation
         trackK->StartBack();
-        Ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kTRUE,zDipole1,zDipole2);
+        ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kTRUE,zDipole1,zDipole2);
         ichamBeg = ichamEnd;
         ichamEnd = 0;
       }
     }
 
-    if (Ok) {
+    if (ok) {
       trackK->SetTrackDir(-1);
       trackK->SetBPFlag(kFALSE);
-      Ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2);
+      ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2);
     }
-    if (Ok) trackK->SetTrackQuality(0); // compute "track quality"
+    if (ok) trackK->SetTrackQuality(0); // compute "track quality"
     else trackK->SetRecover(-1); // mark candidate to be removed
 
     // Majority 3 of 4 in first 2 stations
@@ -1855,7 +1862,7 @@ void AliMUONEventReconstructor::FollowTracksK(void)
 }
 
 //__________________________________________________________________________
-Bool_t AliMUONEventReconstructor::CheckCandidateK(Int_t icand, Int_t nSeeds)
+Bool_t AliMUONEventReconstructor::CheckCandidateK(Int_t icand, Int_t nSeeds) const
 {
   // Discards track candidate if it will produce the double track (having
   // the same seed segment hits as hits of a good track found before)