Adding code to patch St1 HV mapping
[u/mrichter/AliRoot.git] / MUON / MUONFakes.C
index 70ed78f..5827617 100644 (file)
@@ -1,30 +1,24 @@
 #if !defined(__CINT__) || defined(__MAKECINT__)
 // ROOT includes
-#include <TTree.h>
 #include <TFile.h>
 #include <TH1.h>
 #include <TCanvas.h>
 #include <Riostream.h>
 #include <TROOT.h>
-#include <TClonesArray.h>
-#include <TGeoGlobalMagField.h>
+#include <TObjArray.h>
 #include <TArrayI.h>
 
 // STEER includes
 #include "AliLog.h"
-#include "AliMagF.h"
 #include "AliESDEvent.h"
 #include "AliESDMuonTrack.h"
-#include "AliESDMuonCluster.h"
-#include "AliCDBPath.h"
-#include "AliCDBEntry.h"
 #include "AliCDBManager.h"
   
 // MUON includes
+#include "AliMUONCDB.h"
 #include "AliMUONTrack.h"
 #include "AliMUONVTrackStore.h"
 #include "AliMUONTrackParam.h"
-#include "AliMUONTrackExtrap.h"
 #include "AliMUONESDInterface.h"
 #include "AliMUONRecoCheck.h"
 #include "AliMUONVCluster.h"
 /// Results are saved in the root file Fakes.root
 /// Results are relevent provided that you use the same recoParams as for the reconstruction
 
-void Prepare(AliMUONRecoParam *&recoParam, Double_t &sigmaCut);
-TTree* GetESDTree(TFile *esdFile);
-Bool_t TrackMatched(AliMUONTrack &track, AliMUONTrack &trackRef, Float_t &fractionOfMatchCluster, Double_t sigmaCut);
-Bool_t IsRecontructible(AliMUONTrack &track, AliMUONRecoParam &recoParam);
-AliMUONTrack* MatchWithTrackRef(AliESDMuonTrack &muonTrack, AliMUONVTrackStore &trackRefStore,
-                               Float_t &fractionOfMatchCluster, Bool_t useLabel, Double_t sigmaCut);
-Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore, AliMUONRecoParam &recoParam,
-                          Bool_t useLabel, Double_t sigmaCut, TH1F &hFractionOfConnectedClusters);
+UInt_t requestedStationMask = 0;
+Bool_t request2ChInSameSt45 = kFALSE;
+Double_t sigmaCut = -1.;
+
+Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore,
+                          Bool_t useLabel, TH1F &hFractionOfConnectedClusters);
 
 //-----------------------------------------------------------------------
 void MUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent = -1,
-              const TString esdFileName = "AliESDs.root", const TString SimDir = "./generated/")
+              const TString esdFileName = "AliESDs.root", const TString SimDir = "./generated/",
+              const TString ocdbPath = "local://$ALICE_ROOT/OCDB")
 {
   
   //Reset ROOT and connect tree file
@@ -86,19 +79,21 @@ void MUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent =
   TH1F *hPhiM = new TH1F("hPhiM"," matched track phi distribution",100,-1.,9.);
   TH1F *hPhiF = new TH1F("hPhiF"," fake track phi distribution",100,-1.,9.);
   
-  // prepare for analysis
-  AliMUONRecoParam *recoParam = 0x0;
-  Double_t sigmaCut = -1;
-  Prepare(recoParam, sigmaCut);
+  // link to reconstructed and simulated tracks
+  AliMUONRecoCheck rc(esdFileName, SimDir);
   
-  // link to reconstructed tracks
-  TFile* esdFile = TFile::Open(esdFileName);
-  TTree* esdTree = GetESDTree(esdFile);
-  AliESDEvent* esd = new AliESDEvent();
-  esd->ReadFromTree(esdTree);
+  // load necessary data from OCDB
+  AliCDBManager::Instance()->SetDefaultStorage(ocdbPath);
+  AliCDBManager::Instance()->SetRun(rc.GetRunNumber());
+  AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
+  if (!recoParam) return;
   
-  // link to simulated tracks
-  AliMUONRecoCheck rc(esdFileName, SimDir);
+  // get sigma cut from recoParam to associate clusters with TrackRefs in case the label are not used
+  sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking();
+  // compute the mask of requested stations from recoParam
+  for (Int_t i = 0; i < 5; i++) if (recoParam->RequestStation(i)) requestedStationMask |= ( 1 << i );
+  // get from recoParam whether a track need 2 chambers hit in the same station (4 or 5) or not to be reconstructible
+  request2ChInSameSt45 = !recoParam->MakeMoreTrackCandidates();
   
   // initialize global counters
   Int_t nReconstructibleTracks = 0;
@@ -118,46 +113,45 @@ void MUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent =
   
   // Loop over ESD events
   FirstEvent = TMath::Max(0, FirstEvent);
-  LastEvent = (LastEvent>=0) ? TMath::Min((Int_t)esdTree->GetEntries() - 1, LastEvent) : (Int_t)esdTree->GetEntries() - 1;
+  LastEvent = (LastEvent>=0) ? TMath::Min(rc.NumberOfEvents() - 1, LastEvent) : rc.NumberOfEvents() - 1;
   for (Int_t iEvent = FirstEvent; iEvent <= LastEvent; iEvent++) {
     
-    // get the ESD of current event
-    esdTree->GetEvent(iEvent);
-    if (!esd) {
-      Error("CheckESD", "no ESD object found for event %d", iEvent);
-      return;
-    }
-    
-    // convert TrackRef to MUON tracks
+    // get reconstructed and simulated tracks
+    AliMUONVTrackStore* muonTrackStore = rc.ReconstructedTracks(iEvent, kFALSE);
     AliMUONVTrackStore* trackRefStore = rc.TrackRefs(iEvent);
+    if (!muonTrackStore || !trackRefStore) continue;
     
     // count the number of reconstructible tracks
     TIter next(trackRefStore->CreateIterator());
     AliMUONTrack* trackRef;
     while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
-      if (IsRecontructible(*trackRef,*recoParam)) nReconstructibleTracks++;
+      if (trackRef->IsValid(requestedStationMask, request2ChInSameSt45)) nReconstructibleTracks++;
     }
     
     // loop over ESD tracks
     Int_t nTrackerTracks = 0;
     trackReconstructedYet = kFALSE;
     AliMUONVTrackStore *fakeTrackStore = AliMUONESDInterface::NewTrackStore();
+    const AliESDEvent* esd = rc.GetESDEvent();
     Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
     for (Int_t iTrack = 0; iTrack <  nTracks;  iTrack++) {
       
-      AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
+      AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack);
       
       // skip ghosts
-      if (!muonTrack->ContainTrackerData()) continue;
+      if (!esdTrack->ContainTrackerData()) continue;
       nTrackerTracks++;
       
+      // find the corresponding MUON track
+      AliMUONTrack* muonTrack = (AliMUONTrack*) muonTrackStore->FindObject(esdTrack->GetUniqueID());
+      
       // get track info
-      Int_t nClusters = muonTrack->GetNClusters();
-      Double_t normalizedChi2 = muonTrack->GetChi2() / (2. * muonTrack->GetNHit() - 5);
-      Double_t p = muonTrack->P();
-      Double_t pT = muonTrack->Pt();
-      Double_t eta = muonTrack->Eta();
-      Double_t phi = muonTrack->Phi();
+      Int_t nClusters = esdTrack->GetNClusters();
+      Double_t normalizedChi2 = esdTrack->GetChi2() / (2. * esdTrack->GetNHit() - 5);
+      Double_t p = esdTrack->P();
+      Double_t pT = esdTrack->Pt();
+      Double_t eta = esdTrack->Eta();
+      Double_t phi = esdTrack->Phi();
       
       // fill global histograms
       hNumberOfClusters->Fill(nClusters);
@@ -168,21 +162,21 @@ void MUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent =
       hPhi->Fill(phi);
       
       // try to match the reconstructed track with a simulated one
-      Float_t fractionOfMatchCluster = 0.;
-      AliMUONTrack* matchedTrackRef = MatchWithTrackRef(*muonTrack, *trackRefStore, fractionOfMatchCluster, useLabel, sigmaCut);
+      Int_t nMatchClusters = 0;
+      AliMUONTrack* matchedTrackRef = rc.FindCompatibleTrack(*muonTrack, *trackRefStore, nMatchClusters, useLabel, sigmaCut);
       
       // take actions according to matching result
       if (matchedTrackRef) {
        
        // global counter
        nTotMatchedTracks++;
-       if (!IsRecontructible(*matchedTrackRef,*recoParam)) {
+       if (!matchedTrackRef->IsValid(requestedStationMask, request2ChInSameSt45)) {
          trackReconstructedYet = kTRUE;
          nTotTracksReconstructedYet++;
        }
        
        // fill histograms
-       hFractionOfMatchedClusters->Fill(fractionOfMatchCluster);
+       hFractionOfMatchedClusters->Fill(((Float_t) nMatchClusters) / ((Float_t) nClusters));
        hNumberOfClustersMC->Fill(matchedTrackRef->GetNClusters());
        hNumberOfClustersM->Fill(nClusters);
        hChi2PerDofM->Fill(normalizedChi2);
@@ -208,7 +202,7 @@ void MUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent =
        hPhiF->Fill(phi);
        
        // store fake tracks
-       AliMUONESDInterface::Add(*muonTrack, *fakeTrackStore);
+       fakeTrackStore->Add(*muonTrack);
        
       }
       
@@ -223,8 +217,7 @@ void MUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent =
     if (fakeTrackStore->GetSize() > 0) {
       
       // remove the most connected fake tracks
-      Int_t nFreeMissingTracks = RemoveConnectedFakes(*fakeTrackStore, *trackRefStore, *recoParam,
-                                                     useLabel, sigmaCut, *hFractionOfConnectedClusters);
+      Int_t nFreeMissingTracks = RemoveConnectedFakes(*fakeTrackStore, *trackRefStore,  useLabel, *hFractionOfConnectedClusters);
       
       // remove the remaining free reconstructible tracks
       Int_t nAdditionalTracks = fakeTrackStore->GetSize() - nFreeMissingTracks;
@@ -354,179 +347,11 @@ void MUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent =
   cout << "REMINDER: results are relevent provided that you use the same recoParams as for the reconstruction" << endl;
   cout << endl;
   
-  // clear memory
-  delete esd;
-  esdFile->Close();
-  delete recoParam;
-  
-}
-
-//-----------------------------------------------------------------------
-void Prepare(AliMUONRecoParam *&recoParam, Double_t &sigmaCut)
-{
-  /// Set the magnetic field and return recoParam and sigmaCut to associate cluster and trackRef
-  
-  // prepare OCDB access
-  AliCDBManager* man = AliCDBManager::Instance();
-  man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
-  man->SetRun(0);
-  
-  // set  mag field 
-  // waiting for mag field in CDB 
-  if (!TGeoGlobalMagField::Instance()->GetField()) {
-    printf("Loading field map...\n");
-    AliMagF* field = new AliMagF("Maps","Maps",2,1.,1., 10.,AliMagF::k5kG);
-    TGeoGlobalMagField::Instance()->SetField(field);
-  }
-  // set the magnetic field for track extrapolations
-  AliMUONTrackExtrap::SetField();
-  
-  // Load initial reconstruction parameters from OCDB
-  AliCDBPath path("MUON","Calib","RecoParam");
-  AliCDBEntry *entry=man->Get(path.GetPath());
-  if(entry) {
-    recoParam = dynamic_cast<AliMUONRecoParam*>(entry->GetObject());
-    entry->SetOwner(0);
-    AliCDBManager::Instance()->UnloadFromCache(path.GetPath());
-  }
-  if (!recoParam) {
-    printf("Couldn't find RecoParam object in OCDB: create default one");
-    recoParam = AliMUONRecoParam::GetLowFluxParam();
-  }
-  
-  Info("MUONFakes", "\n recontruction parameters:");
-  recoParam->Print("FULL");
-  AliMUONESDInterface::ResetTracker(recoParam);
-  
-  // sigma cut to associate clusters with TrackRefs in case the label are not used
-  sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking(); 
-  
-}
-
-//-----------------------------------------------------------------------
-TTree* GetESDTree(TFile *esdFile)
-{
-  /// Check that the file is properly open
-  /// Return pointer to the ESD Tree
-  
-  if (!esdFile || !esdFile->IsOpen()) {
-    Error("GetESDTree", "opening ESD file failed");
-    exit(-1);
-  }
-  
-  TTree* tree = (TTree*) esdFile->Get("esdTree");
-  if (!tree) {
-    Error("GetESDTree", "no ESD tree found");
-    exit(-1);
-  }
-  
-  return tree;
-  
-}
-
-//-----------------------------------------------------------------------
-Bool_t TrackMatched(AliMUONTrack &track, AliMUONTrack &trackRef, Float_t &fractionOfMatchCluster, Double_t sigmaCut)
-{
-  /// Try to match 2 tracks
-  
-  Bool_t compTrack[10];
-  Int_t nMatchClusters = track.CompatibleTrack(&trackRef, sigmaCut, compTrack);
-  fractionOfMatchCluster = ((Float_t)nMatchClusters) / ((Float_t)track.GetNClusters());
-  
-  if ((compTrack[0] || compTrack[1] || compTrack[2] || compTrack[3]) && // at least 1 cluster matched in st 1 & 2
-      (compTrack[6] || compTrack[7] || compTrack[8] || compTrack[9]) && // at least 1 cluster matched in st 4 & 5
-      fractionOfMatchCluster > 0.5) return kTRUE;                       // more than 50% of clusters matched
-  else return kFALSE;
-  
-}
-
-//-----------------------------------------------------------------------
-Bool_t IsRecontructible(AliMUONTrack &track, AliMUONRecoParam &recoParam)
-{
-  /// Check il the track is reconstructible
-  Int_t nMinChHitInSt45 = (recoParam.MakeMoreTrackCandidates()) ? 2 : 3;
-  Int_t currentCh, previousCh = -1, nChHitInSt45 = 0;
-  Bool_t clusterInSt[5];
-  for (Int_t iSt = 0; iSt < 5; iSt++) clusterInSt[iSt] = !recoParam.RequestStation(iSt);
-  
-  AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->First());
-  while (trackParam) {
-    
-    currentCh = trackParam->GetClusterPtr()->GetChamberId();
-    
-    clusterInSt[currentCh/2] = kTRUE;
-    
-    if (currentCh > 5 && currentCh != previousCh) {
-      nChHitInSt45++;
-      previousCh = currentCh;
-    }
-    
-    trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->After(trackParam));
-  }
-  
-  return (clusterInSt[0] && clusterInSt[1] && clusterInSt[2] &&
-          clusterInSt[3] && clusterInSt[4] && nChHitInSt45 >= nMinChHitInSt45);
-  
-}
-
-//-----------------------------------------------------------------------
-AliMUONTrack* MatchWithTrackRef(AliESDMuonTrack &muonTrack, AliMUONVTrackStore &trackRefStore,
-                               Float_t &fractionOfMatchCluster, Bool_t useLabel, Double_t sigmaCut)
-{
-  /// Return if the trackRef matched with the reconstructed track and the fraction of matched clusters
-  
-  AliMUONTrack *matchedTrackRef = 0x0;
-  fractionOfMatchCluster = 0.;
-  
-  if (useLabel) { // by using the MC label
-    
-    // get the corresponding simulated track if any
-    Int_t label = muonTrack.GetLabel();
-    matchedTrackRef = (AliMUONTrack*) trackRefStore.FindObject(label);
-    
-    // get the fraction of matched clusters
-    if (matchedTrackRef) {
-      Int_t nMatchClusters = 0;
-      if (muonTrack.ClustersStored()) {
-       AliESDMuonCluster* cluster = (AliESDMuonCluster*) muonTrack.GetClusters().First();
-       while (cluster) {
-         if (cluster->GetLabel() == label) nMatchClusters++;
-         cluster = (AliESDMuonCluster*) muonTrack.GetClusters().After(cluster);
-       }
-      }
-      fractionOfMatchCluster = ((Float_t)nMatchClusters) / ((Float_t)muonTrack.GetNClusters());
-    }
-    
-  } else { // by comparing cluster/TrackRef positions
-    
-    // convert ESD track to MUON track
-    AliMUONTrack track;
-    AliMUONESDInterface::ESDToMUON(muonTrack,track);
-    
-    // look for the corresponding simulated track if any
-    TIter next(trackRefStore.CreateIterator());
-    AliMUONTrack* trackRef;
-    while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
-      
-      // check compatibility
-      Float_t f = 0.;
-      if (TrackMatched(track, *trackRef, f, sigmaCut)) {
-       matchedTrackRef = trackRef;
-       fractionOfMatchCluster = f;
-       break;
-      }
-      
-    }
-    
-  }
-  
-  return matchedTrackRef;
-  
 }
 
 //-----------------------------------------------------------------------
-Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore, AliMUONRecoParam &recoParam,
-                         Bool_t useLabel, Double_t sigmaCut, TH1F &hFractionOfConnectedClusters)
+Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore,
+                          Bool_t useLabel, TH1F &hFractionOfConnectedClusters)
 {
   /// loop over reconstructible TrackRef not associated with reconstructed track:
   /// for each of them, find and remove the most connected the fake track, if any,
@@ -541,13 +366,13 @@ Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStor
   while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
     
     // skip not reconstructible trackRefs
-    if (!IsRecontructible(*trackRef,recoParam)) continue;
+    if (!trackRef->IsValid(requestedStationMask, request2ChInSameSt45)) continue;
     
     Int_t label = trackRef->GetUniqueID();
     
     // look for the most connected fake track
     AliMUONTrack *connectedFake = 0x0;
-    Float_t fractionOfConnectedClusters = 0.;
+    Double_t fractionOfConnectedClusters = 0.;
     TIter next2(fakeTrackStore.CreateIterator());
     AliMUONTrack* fakeTrack;
     while ( ( fakeTrack = static_cast<AliMUONTrack*>(next2()) ) ) {
@@ -560,14 +385,14 @@ Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStor
            nConnectedClusters++;
       } else { // by comparing cluster/TrackRef positions
        Bool_t compTrack[10];
-       nConnectedClusters = fakeTrack->CompatibleTrack(trackRef, sigmaCut, compTrack);
+       nConnectedClusters = fakeTrack->FindCompatibleClusters(*trackRef, sigmaCut, compTrack);
       }
       
       // skip non-connected fake tracks
       if (nConnectedClusters == 0) continue;
       
       // check if it is the most connected fake track
-      Float_t f = ((Float_t)nConnectedClusters) / ((Float_t)fakeTrack->GetNClusters());
+      Double_t f = ((Double_t)nConnectedClusters) / ((Double_t)fakeTrack->GetNClusters());
       if (f > fractionOfConnectedClusters) {
        connectedFake = fakeTrack;
        fractionOfConnectedClusters = f;