Flexible pt range for the efficiency histogramming
[u/mrichter/AliRoot.git] / MUON / MUONRefit.C
index e7bed8a..7b94be6 100644 (file)
 #include <TROOT.h>
 
 // STEER includes
-#include "AliMagFMaps.h"
-#include "AliTracker.h"
 #include "AliESDEvent.h"
 #include "AliESDMuonTrack.h"
-#include "AliRecoParam.h"
 #include "AliCDBManager.h"
 #include "AliGeomManager.h"
 
 // MUON includes
-#include "AliMpCDB.h"
+#include "AliMUONCDB.h"
 #include "AliMUONRecoParam.h"
 #include "AliMUONESDInterface.h"
 #include "AliMUONRefitter.h"
 #include "AliMUONVDigit.h"
-#include "AliMUONVCluster.h"
-#include "AliMUONVClusterStore.h"
 #include "AliMUONTrack.h"
 #include "AliMUONVTrackStore.h"
 #include "AliMUONTrackParam.h"
-#include "AliMUONTrackExtrap.h"
 #endif
 
 const Int_t printLevel = 1;
 
-void Prepare();
 TTree* GetESDTree(TFile *esdFile);
 
 //-----------------------------------------------------------------------
-void MUONRefit(Int_t nevents = -1, const char* esdFileNameIn = "AliESDs.root", const char* esdFileNameOut = "AliESDs_New.root")
+void MUONRefit(Int_t nevents = -1, const char* esdFileNameIn = "AliESDs.root", const char* esdFileNameOut = "AliESDs_New.root",
+              const char* geoFilename = "geometry.root", const char* ocdbPath = "local://$ALICE_ROOT/OCDB")
 {
+  /// Example of muon refitting:
   /// refit ESD tracks from ESD pads (i.e. re-clusterized the attached ESD clusters);
   /// reset the charge of the digit using their raw charge before refitting;
   /// compare results with original ESD tracks; 
   /// write results in a new ESD file
   
-  // prepare the refitting
-  gRandom->SetSeed(1);
-  Prepare();
-  AliMUONESDInterface esdInterface;
-  AliMUONRefitter refitter;
-  refitter.Connect(&esdInterface);
-  
   // open the ESD file and tree
   TFile* esdFile = TFile::Open(esdFileNameIn);
   TTree* esdTree = GetESDTree(esdFile);
   
-  // create the ESD output tree
-  gROOT->cd();
-  TTree* newESDTree = esdTree->CloneTree(0);
-  
   // connect ESD event to the ESD tree
   AliESDEvent* esd = new AliESDEvent();
   esd->ReadFromTree(esdTree);
-
+  
+  // create the ESD output file and tree
+  TFile* newESDFile = TFile::Open(esdFileNameOut, "RECREATE");
+  newESDFile->SetCompressionLevel(2);
+  
+  // connect ESD event to the ESD tree (recreate track branch for backward compatibility)
+  esdTree->SetBranchStatus("*MuonTracks*",0);
+  TTree* newESDTree = esdTree->CloneTree(0);
+  esdTree->SetBranchStatus("*MuonTracks*",1);
+  esd->WriteToTree(newESDTree);
+  
+  // get run number
+  if (esdTree->GetEvent(0) <= 0) {
+    Error("MUONRefit", "no ESD object found for event 0");
+    return;
+  }
+  Int_t runNumber = esd->GetRunNumber();
+  
+  // Import TGeo geometry
+  if (!gGeoManager) {
+    AliGeomManager::LoadGeometry(geoFilename);
+    if (!gGeoManager) {
+      Error("MUONRefit", "getting geometry from file %s failed", "generated/galice.root");
+      exit(-1);
+    }
+  }
+  
+  // load necessary data from OCDB
+  AliCDBManager::Instance()->SetDefaultStorage(ocdbPath);
+  AliCDBManager::Instance()->SetRun(runNumber);
+  if (!AliMUONCDB::LoadField()) return;
+  if (!AliMUONCDB::LoadMapping(kTRUE)) return;
+  
+  // reconstruction parameters for the refitting
+  AliMUONRecoParam* recoParam = AliMUONRecoParam::GetLowFluxParam();
+  Info("MUONRefit", "\n Reconstruction parameters for refitting:");
+  recoParam->Print("FULL");
+  
+  AliMUONESDInterface esdInterface;
+  AliMUONRefitter refitter(recoParam);
+  refitter.Connect(&esdInterface);
+  gRandom->SetSeed(1);
+  
   // timer start...
   TStopwatch timer;
   
@@ -101,16 +127,20 @@ void MUONRefit(Int_t nevents = -1, const char* esdFileNameIn = "AliESDs.root", c
     // get the ESD of current event
     esdTree->GetEvent(iEvent);
     if (!esd) {
-      Error("CheckESD", "no ESD object found for event %d", iEvent);
+      Error("MUONRefit", "no ESD object found for event %d", iEvent);
       return;
     }
     Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks();
     if (nTracks < 1) continue;
     
     // load the current event
-    esdInterface.LoadEvent(*esd);
+    esdInterface.LoadEvent(*esd, kFALSE);
+    
+    // remove digits and clusters from ESD
+    esd->FindListObject("MuonClusters")->Clear("C");
+    esd->FindListObject("MuonPads")->Clear("C");
     
-    // loop over digit to modify their charge
+    // loop over digits to modify their charge
     AliMUONVDigit *digit;
     TIter next(esdInterface.CreateDigitIterator());
     while ((digit = static_cast<AliMUONVDigit*>(next()))) {
@@ -119,6 +149,7 @@ void MUONRefit(Int_t nevents = -1, const char* esdFileNameIn = "AliESDs.root", c
     }
     
     // refit the tracks from digits
+    refitter.SetFirstClusterIndex(0);
     AliMUONVTrackStore* newTrackStore = refitter.ReconstructFromDigits();
     
     //----------------------------------------------//
@@ -141,12 +172,19 @@ void MUONRefit(Int_t nevents = -1, const char* esdFileNameIn = "AliESDs.root", c
       AliMUONTrack* newTrack = (AliMUONTrack*) newTrackStore->FindObject(esdTrack->GetUniqueID());
       
       // replace the content of the current ESD track or remove it if the refitting has failed
-      if (newTrack) {
+      if (newTrack && (!recoParam->ImproveTracks() || newTrack->IsImproved())) {
+       
+       // fill the track info
        Double_t vertex[3] = {esdTrack->GetNonBendingCoor(), esdTrack->GetBendingCoor(), esdTrack->GetZ()};
-       AliMUONESDInterface::MUONToESD(*newTrack, *esdTrack, vertex, esdInterface.GetDigits());
-      } else {
-       esdTracks->Remove(esdTrack);
-      }
+       AliMUONESDInterface::MUONToESD(*newTrack, *esdTrack, vertex);
+       
+       // add the clusters (and the digits if previously there)
+       for (Int_t i = 0; i < newTrack->GetNClusters(); i++) {
+         AliMUONTrackParam *param = static_cast<AliMUONTrackParam*>(newTrack->GetTrackParamAtCluster()->UncheckedAt(i));
+         AliMUONESDInterface::MUONToESD(*(param->GetClusterPtr()), *esd, esdInterface.GetDigits());
+       }
+       
+      } else esdTracks->Remove(esdTrack);
       
       // print initial and re-fitted track parameters at first cluster if any
       if (printLevel>0) {
@@ -179,57 +217,20 @@ void MUONRefit(Int_t nevents = -1, const char* esdFileNameIn = "AliESDs.root", c
   timer.Stop();
   
   // write output ESD tree
-  TFile* newESDFile = TFile::Open(esdFileNameOut, "RECREATE");
-  newESDFile->SetCompressionLevel(2);
+  newESDFile->cd();
   newESDTree->Write();
+  delete newESDTree;
   newESDFile->Close();
   
   // free memory
   esdFile->Close();
-  delete newESDTree;
   delete esd;
+  delete recoParam;
   
   cout<<endl<<"time to refit: R:"<<timer.RealTime()<<" C:"<<timer.CpuTime()<<endl<<endl;
 }
 
 //-----------------------------------------------------------------------
-void Prepare()
-{
-  /// Set the geometry, the magnetic field, the mapping and the reconstruction parameters
-  
-  // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
-  if (!gGeoManager) {
-    AliGeomManager::LoadGeometry("geometry.root");
-    if (!gGeoManager) {
-      Error("MUONRefit", "getting geometry from file %s failed", "generated/galice.root");
-      return;
-    }
-  }
-  
-  // set mag field
-  printf("Loading field map...\n");
-  AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 1, 1., 10., AliMagFMaps::k5kG);
-  AliTracker::SetFieldMap(field, kFALSE);
-  AliMUONTrackExtrap::SetField(AliTracker::GetFieldMap());
-  
-  // Load mapping
-  AliCDBManager* man = AliCDBManager::Instance();
-  man->SetDefaultStorage("local://$ALICE_ROOT");
-  man->SetRun(0);
-  if ( ! AliMpCDB::LoadDDLStore() ) {
-    Error("MUONRefit","Could not access mapping from OCDB !");
-    exit(-1);
-  }
-  
-  // set reconstruction parameters
-  AliMUONRecoParam *muonRecoParam = AliMUONRecoParam::GetLowFluxParam();
-  muonRecoParam->CombineClusterTrackReco(kFALSE);
-  muonRecoParam->Print("FULL");
-  AliRecoParam::Instance()->RegisterRecoParam(muonRecoParam);
-  
-}
-
-//-----------------------------------------------------------------------
 TTree* GetESDTree(TFile *esdFile)
 {
   /// Check that the file is properly open