[u/mrichter/AliRoot.git] / MUON / MUONFakes.C

#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>

// 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 "AliMUONTrack.h"
#include "AliMUONVTrackStore.h"
#include "AliMUONTrackParam.h"
#include "AliMUONTrackExtrap.h"
#include "AliMUONESDInterface.h"
#include "AliMUONRecoCheck.h"
#include "AliMUONVCluster.h"
#include "AliMUONRecoParam.h"
#endif

/// \ingroup macros
/// \file MUONFakes.C
///
/// \author Ph. Pillot, Subatech, March. 2009
///
/// Macro to study fake tracks by comparing reconstructed tracks with TrackRefs
/// 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);

//-----------------------------------------------------------------------
void MUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent = -1,
	       const TString esdFileName = "AliESDs.root", const TString SimDir = "./generated/")
{
  
  //Reset ROOT and connect tree file
  gROOT->Reset();
  
  // File for histograms and histogram booking
  TFile *histoFile = new TFile("Fakes.root", "RECREATE");
  
  TH1F *hNumberOfTracks = new TH1F("hNumberOfTracks","nb of tracks /evt",20,0.,20.);
  TH1F *hNumberOfAdditionalTracks = new TH1F("hNumberOfAdditionalTracks","nb of fake - nb of missing track",20,0.,20.);
  
  TH1F *hNumberOfClusters = new TH1F("hNumberOfClusters","nb of clusters /track",20,0.,20.);
  TH1F *hNumberOfClustersM = new TH1F("hNumberOfClustersM","nb of clusters /matched track",20,0.,20.);
  TH1F *hNumberOfClustersF = new TH1F("hNumberOfClustersF","nb of clusters /fake track",20,0.,20.);
  TH1F *hNumberOfClustersMC = new TH1F("hNumberOfClustersMC","nb of clusters /MC track",20,0.,20.);
  TH1F *hFractionOfMatchedClusters = new TH1F("hFractionOfMatchedClusters","nb of matched clusters / nb of clusters",110,0.,1.1);
  TH1F *hFractionOfConnectedClusters = new TH1F("hFractionOfConnectedClusters","nb of connected clusters / nb of clusters in fake tracks",110,0.,1.1);
  
  TH1F *hChi2PerDof = new TH1F("hChi2PerDof", "track chi2/d.o.f.", 100, 0., 20.);
  TH1F *hChi2PerDofM = new TH1F("hChi2PerDofM", "matched track chi2/d.o.f.", 100, 0., 20.);
  TH1F *hChi2PerDofF = new TH1F("hChi2PerDofF", "fake track chi2/d.o.f.", 100, 0., 20.);
  TH1F *hP = new TH1F("hP", "Muon P distribution (GeV/c)", 100, 0., 200.);
  TH1F *hPM = new TH1F("hPM", "matched track P distribution (GeV/c)", 100, 0., 200.);
  TH1F *hPF = new TH1F("hPF", "fake track P distribution (GeV/c)", 100, 0., 200.);
  TH1F *hPt = new TH1F("hPt", "Muon Pt distribution (GeV/c)", 100, 0., 20.);
  TH1F *hPtM = new TH1F("hPtM", "matched track Pt distribution (GeV/c)", 100, 0., 20.);
  TH1F *hPtF = new TH1F("hPtF", "fake track Pt distribution (GeV/c)", 100, 0., 20.);
  TH1F *hEta = new TH1F("hEta"," Muon pseudo-rapidity distribution",100,-10,0);
  TH1F *hEtaM = new TH1F("hEtaM"," matched track pseudo-rapidity distribution",100,-10,0);
  TH1F *hEtaF = new TH1F("hEtaF"," fake track pseudo-rapidity distribution",100,-10,0);
  TH1F *hPhi = new TH1F("hPhi"," Muon phi distribution",100,-1.,9.);
  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 tracks
  TFile* esdFile = TFile::Open(esdFileName);
  TTree* esdTree = GetESDTree(esdFile);
  AliESDEvent* esd = new AliESDEvent();
  esd->ReadFromTree(esdTree);
  
  // link to simulated tracks
  AliMUONRecoCheck rc(esdFileName, SimDir);
  
  // initialize global counters
  Int_t nEventsWithFake = 0;
  Int_t nEventsWithAdditionalFake = 0;
  Int_t nTotMatchedTracks = 0;
  Int_t nTotTracksReconstructedYet = 0;
  Int_t nTotFakeTracks = 0;
  Int_t nTotAdditionalTracks = 0;
  
  // 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;
  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
    AliMUONVTrackStore* trackRefStore = rc.TrackRefs(iEvent);
    
    // loop over ESD tracks
    Int_t nTrackerTracks = 0;
    AliMUONVTrackStore *fakeTrackStore = AliMUONESDInterface::NewTrackStore();
    Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
    for (Int_t iTrack = 0; iTrack <  nTracks;  iTrack++) {
      
      AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
      
      // skip ghosts
      if (!muonTrack->ContainTrackerData()) continue;
      nTrackerTracks++;
      
      // 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();
      
      // fill global histograms
      hNumberOfClusters->Fill(nClusters);
      hChi2PerDof->Fill(normalizedChi2);
      hP->Fill(p);
      hPt->Fill(pT);
      hEta->Fill(eta);
      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);
      
      // take actions according to matching result
      if (matchedTrackRef) {
	
	// global counter
	nTotMatchedTracks++;
	if (!IsRecontructible(*matchedTrackRef,*recoParam)) nTotTracksReconstructedYet++;
	
	// fill histograms
	hFractionOfMatchedClusters->Fill(fractionOfMatchCluster);
	hNumberOfClustersMC->Fill(matchedTrackRef->GetNClusters());
	hNumberOfClustersM->Fill(nClusters);
	hChi2PerDofM->Fill(normalizedChi2);
	hPM->Fill(p);
	hPtM->Fill(pT);
	hEtaM->Fill(eta);
	hPhiM->Fill(phi);
	
	// remove already matched trackRefs
	trackRefStore->Remove(*matchedTrackRef);
	
      } else {
	
	// global counter
	nTotFakeTracks++;
	
	// fill histograms
	hNumberOfClustersF->Fill(nClusters);
	hChi2PerDofF->Fill(normalizedChi2);
	hPF->Fill(p);
	hPtF->Fill(pT);
	hEtaF->Fill(eta);
	hPhiF->Fill(phi);
	
	// store fake tracks
	AliMUONESDInterface::Add(*muonTrack, *fakeTrackStore);
	
      }
      
    } // end of loop over ESD tracks
    
    // fill histograms
    hNumberOfTracks->Fill(nTrackerTracks);
    
    // count the number the additional fake tracks
    if (fakeTrackStore->GetSize() > 0) {
      
      // remove the most connected fake tracks
      Int_t nFreeMissingTracks = RemoveConnectedFakes(*fakeTrackStore, *trackRefStore, *recoParam,
						      useLabel, sigmaCut, *hFractionOfConnectedClusters);
      
      // remove the remaining free reconstructible tracks
      Int_t nAdditionalTracks = fakeTrackStore->GetSize() - nFreeMissingTracks;
      
      // fill histograms
      nEventsWithFake++;
      if (nAdditionalTracks > 0) {
	nEventsWithAdditionalFake++;
	nTotAdditionalTracks += nAdditionalTracks;
	hNumberOfAdditionalTracks->Fill(nAdditionalTracks);
      }
      
    }
    
    delete fakeTrackStore;
    
  } // end of loop over events
  
  // plot results
  TCanvas cFakesSummary("cFakesSummary","cFakesSummary",900,600);
  cFakesSummary.Divide(3,2);
  cFakesSummary.cd(1);
  cFakesSummary.GetPad(1)->SetLogy();
  hNumberOfClusters->Draw();
  hNumberOfClusters->SetMinimum(0.5);
  hNumberOfClustersM->Draw("same");
  hNumberOfClustersM->SetLineColor(4);
  hNumberOfClustersF->Draw("same");
  hNumberOfClustersF->SetLineColor(2);
  hNumberOfClustersF->SetFillColor(2);
  hNumberOfClustersF->SetFillStyle(3017);
  cFakesSummary.cd(2);
  cFakesSummary.GetPad(2)->SetLogy();
  hChi2PerDof->Draw();
  hChi2PerDof->SetMinimum(0.5);
  hChi2PerDofM->Draw("same");
  hChi2PerDofM->SetLineColor(4);
  hChi2PerDofF->Draw("same");
  hChi2PerDofF->SetLineColor(2);
  hChi2PerDofF->SetFillColor(2);
  hChi2PerDofF->SetFillStyle(3017);
  cFakesSummary.cd(3);
  cFakesSummary.GetPad(3)->SetLogy();
  hP->Draw();
  hP->SetMinimum(0.5);
  hPM->Draw("same");
  hPM->SetLineColor(4);
  hPF->Draw("same");
  hPF->SetLineColor(2);
  hPF->SetFillColor(2);
  hPF->SetFillStyle(3017);
  cFakesSummary.cd(4);
  cFakesSummary.GetPad(4)->SetLogy();
  hPt->Draw();
  hPt->SetMinimum(0.5);
  hPtM->Draw("same");
  hPtM->SetLineColor(4);
  hPtF->Draw("same");
  hPtF->SetLineColor(2);
  hPtF->SetFillColor(2);
  hPtF->SetFillStyle(3017);
  cFakesSummary.cd(5);
  cFakesSummary.GetPad(5)->SetLogy();
  hEta->Draw();
  hEta->SetMinimum(0.5);
  hEtaM->Draw("same");
  hEtaM->SetLineColor(4);
  hEtaF->Draw("same");
  hEtaF->SetLineColor(2);
  hEtaF->SetFillColor(2);
  hEtaF->SetFillStyle(3017);
  cFakesSummary.cd(6);
  cFakesSummary.GetPad(6)->SetLogy();
  hPhi->Draw();
  hPhi->SetMinimum(0.5);
  hPhiM->Draw("same");
  hPhiM->SetLineColor(4);
  hPhiF->Draw("same");
  hPhiF->SetLineColor(2);
  hPhiF->SetFillColor(2);
  hPhiF->SetFillStyle(3017);
  
  // save results
  histoFile->cd();
  histoFile->Write();
  cFakesSummary.Write();
  histoFile->Close();
  
  // print results
  cout << endl;
  cout << "- Number of matched tracks: " << nTotMatchedTracks << endl;
  cout << "  (including " << nTotTracksReconstructedYet << " tracks matched with a TrackRef that is not reconstructible)" << endl;
  cout << "- Number of fake tracks: " << nTotFakeTracks << endl;
  cout << "  (including " << nTotAdditionalTracks << " additional tracks (compared to the number of expected ones))" << endl;
  cout << "- Number of events with fake track(s): " << nEventsWithFake << endl;
  cout << "  (including " << nEventsWithAdditionalFake << " events with additional tracks)" << endl;
  cout << endl;
  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)
{
  /// loop over reconstructible TrackRef not associated with reconstructed track:
  /// for each of them, find and remove the most connected the fake track, if any,
  /// and fill the histograms with the fraction of connected clusters.
  /// Return the number of reconstructible track not connected to any fake
  
  Int_t nFreeMissingTracks = 0;
  
  // loop over trackRefs
  TIter next(trackRefStore.CreateIterator());
  AliMUONTrack* trackRef;
  while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
    
    // skip not reconstructible trackRefs
    if (!IsRecontructible(*trackRef,recoParam)) continue;
    
    Int_t label = trackRef->GetUniqueID();
    
    // look for the most connected fake track
    AliMUONTrack *connectedFake = 0x0;
    Float_t fractionOfConnectedClusters = 0.;
    TIter next2(fakeTrackStore.CreateIterator());
    AliMUONTrack* fakeTrack;
    while ( ( fakeTrack = static_cast<AliMUONTrack*>(next2()) ) ) {
      
      // get the number of connected clusters
      Int_t nConnectedClusters = 0;
      if (useLabel) { // by using the MC label
	for (Int_t iCl = 0; iCl < fakeTrack->GetNClusters(); iCl++)
	  if (((AliMUONTrackParam*) fakeTrack->GetTrackParamAtCluster()->UncheckedAt(iCl))->GetClusterPtr()->GetMCLabel() == label)
	    nConnectedClusters++;
      } else { // by comparing cluster/TrackRef positions
	Bool_t compTrack[10];
	nConnectedClusters = fakeTrack->CompatibleTrack(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());
      if (f > fractionOfConnectedClusters) {
	connectedFake = fakeTrack;
	fractionOfConnectedClusters = f;
      }
      
    }
    
    // remove the most connected fake track
    if (connectedFake) {
      hFractionOfConnectedClusters.Fill(fractionOfConnectedClusters);
      fakeTrackStore.Remove(*connectedFake);
    } else nFreeMissingTracks++;
    
  }
  
  return nFreeMissingTracks;
  
}
Commit	Line	Data
66085093	1	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	2	// ROOT includes
	3	#include <TTree.h>
	4	#include <TFile.h>
	5	#include <TH1.h>
	6	#include <TCanvas.h>
	7	#include <Riostream.h>
	8	#include <TROOT.h>
	9	#include <TClonesArray.h>
	10	#include <TGeoGlobalMagField.h>
	11
	12	// STEER includes
	13	#include "AliLog.h"
	14	#include "AliMagF.h"
	15	#include "AliESDEvent.h"
	16	#include "AliESDMuonTrack.h"
	17	#include "AliESDMuonCluster.h"
	18	#include "AliCDBPath.h"
	19	#include "AliCDBEntry.h"
	20	#include "AliCDBManager.h"
	21
	22	// MUON includes
	23	#include "AliMUONTrack.h"
	24	#include "AliMUONVTrackStore.h"
	25	#include "AliMUONTrackParam.h"
	26	#include "AliMUONTrackExtrap.h"
	27	#include "AliMUONESDInterface.h"
	28	#include "AliMUONRecoCheck.h"
	29	#include "AliMUONVCluster.h"
	30	#include "AliMUONRecoParam.h"
	31	#endif
	32
	33	/// \ingroup macros
	34	/// \file MUONFakes.C
	35	///
	36	/// \author Ph. Pillot, Subatech, March. 2009
	37	///
	38	/// Macro to study fake tracks by comparing reconstructed tracks with TrackRefs
	39	/// Results are saved in the root file Fakes.root
	40	/// Results are relevent provided that you use the same recoParams as for the reconstruction
	41
	42	void Prepare(AliMUONRecoParam *&recoParam, Double_t &sigmaCut);
	43	TTree* GetESDTree(TFile *esdFile);
	44	Bool_t TrackMatched(AliMUONTrack &track, AliMUONTrack &trackRef, Float_t &fractionOfMatchCluster, Double_t sigmaCut);
	45	Bool_t IsRecontructible(AliMUONTrack &track, AliMUONRecoParam &recoParam);
	46	AliMUONTrack* MatchWithTrackRef(AliESDMuonTrack &muonTrack, AliMUONVTrackStore &trackRefStore,
	47	Float_t &fractionOfMatchCluster, Bool_t useLabel, Double_t sigmaCut);
	48	Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore, AliMUONRecoParam &recoParam,
	49	Bool_t useLabel, Double_t sigmaCut, TH1F &hFractionOfConnectedClusters);
	50
	51	//-----------------------------------------------------------------------
	52	void MUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent = -1,
	53	const TString esdFileName = "AliESDs.root", const TString SimDir = "./generated/")
	54	{
	55
	56	//Reset ROOT and connect tree file
	57	gROOT->Reset();
	58
	59	// File for histograms and histogram booking
	60	TFile *histoFile = new TFile("Fakes.root", "RECREATE");
	61
	62	TH1F *hNumberOfTracks = new TH1F("hNumberOfTracks","nb of tracks /evt",20,0.,20.);
	63	TH1F *hNumberOfAdditionalTracks = new TH1F("hNumberOfAdditionalTracks","nb of fake - nb of missing track",20,0.,20.);
	64
65	TH1F *hNumberOfClusters = new TH1F("hNumberOfClusters","nb of clusters /track",20,0.,20.);
66	TH1F *hNumberOfClustersM = new TH1F("hNumberOfClustersM","nb of clusters /matched track",20,0.,20.);
67	TH1F *hNumberOfClustersF = new TH1F("hNumberOfClustersF","nb of clusters /fake track",20,0.,20.);
68	TH1F *hNumberOfClustersMC = new TH1F("hNumberOfClustersMC","nb of clusters /MC track",20,0.,20.);
69	TH1F *hFractionOfMatchedClusters = new TH1F("hFractionOfMatchedClusters","nb of matched clusters / nb of clusters",110,0.,1.1);
70	TH1F *hFractionOfConnectedClusters = new TH1F("hFractionOfConnectedClusters","nb of connected clusters / nb of clusters in fake tracks",110,0.,1.1);
71
72	TH1F *hChi2PerDof = new TH1F("hChi2PerDof", "track chi2/d.o.f.", 100, 0., 20.);
73	TH1F *hChi2PerDofM = new TH1F("hChi2PerDofM", "matched track chi2/d.o.f.", 100, 0., 20.);
74	TH1F *hChi2PerDofF = new TH1F("hChi2PerDofF", "fake track chi2/d.o.f.", 100, 0., 20.);
75	TH1F *hP = new TH1F("hP", "Muon P distribution (GeV/c)", 100, 0., 200.);
76	TH1F *hPM = new TH1F("hPM", "matched track P distribution (GeV/c)", 100, 0., 200.);
77	TH1F *hPF = new TH1F("hPF", "fake track P distribution (GeV/c)", 100, 0., 200.);
78	TH1F *hPt = new TH1F("hPt", "Muon Pt distribution (GeV/c)", 100, 0., 20.);
79	TH1F *hPtM = new TH1F("hPtM", "matched track Pt distribution (GeV/c)", 100, 0., 20.);
80	TH1F *hPtF = new TH1F("hPtF", "fake track Pt distribution (GeV/c)", 100, 0., 20.);
81	TH1F *hEta = new TH1F("hEta"," Muon pseudo-rapidity distribution",100,-10,0);
82	TH1F *hEtaM = new TH1F("hEtaM"," matched track pseudo-rapidity distribution",100,-10,0);
83	TH1F *hEtaF = new TH1F("hEtaF"," fake track pseudo-rapidity distribution",100,-10,0);
84	TH1F *hPhi = new TH1F("hPhi"," Muon phi distribution",100,-1.,9.);
85	TH1F *hPhiM = new TH1F("hPhiM"," matched track phi distribution",100,-1.,9.);
86	TH1F *hPhiF = new TH1F("hPhiF"," fake track phi distribution",100,-1.,9.);
87
88	// prepare for analysis
89	AliMUONRecoParam *recoParam = 0x0;
90	Double_t sigmaCut = -1;
91	Prepare(recoParam, sigmaCut);
92
93	// link to reconstructed tracks
94	TFile* esdFile = TFile::Open(esdFileName);
95	TTree* esdTree = GetESDTree(esdFile);
96	AliESDEvent* esd = new AliESDEvent();
97	esd->ReadFromTree(esdTree);
98
99	// link to simulated tracks
100	AliMUONRecoCheck rc(esdFileName, SimDir);
101
102	// initialize global counters
103	Int_t nEventsWithFake = 0;
104	Int_t nEventsWithAdditionalFake = 0;
105	Int_t nTotMatchedTracks = 0;
106	Int_t nTotTracksReconstructedYet = 0;
107	Int_t nTotFakeTracks = 0;
108	Int_t nTotAdditionalTracks = 0;
109
110	// Loop over ESD events
111	FirstEvent = TMath::Max(0, FirstEvent);
112	LastEvent = (LastEvent>=0) ? TMath::Min((Int_t)esdTree->GetEntries() - 1, LastEvent) : (Int_t)esdTree->GetEntries() - 1;
113	for (Int_t iEvent = FirstEvent; iEvent <= LastEvent; iEvent++) {
114
115	// get the ESD of current event
116	esdTree->GetEvent(iEvent);
117	if (!esd) {
118	Error("CheckESD", "no ESD object found for event %d", iEvent);
119	return;
120	}
121
122	// convert TrackRef to MUON tracks
123	AliMUONVTrackStore* trackRefStore = rc.TrackRefs(iEvent);
124
125	// loop over ESD tracks
126	Int_t nTrackerTracks = 0;
127	AliMUONVTrackStore *fakeTrackStore = AliMUONESDInterface::NewTrackStore();
128	Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
129	for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
130
131	AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
132
133	// skip ghosts
134	if (!muonTrack->ContainTrackerData()) continue;
135	nTrackerTracks++;
136
137	// get track info
138	Int_t nClusters = muonTrack->GetNClusters();
139	Double_t normalizedChi2 = muonTrack->GetChi2() / (2. * muonTrack->GetNHit() - 5);
140	Double_t p = muonTrack->P();
141	Double_t pT = muonTrack->Pt();
142	Double_t eta = muonTrack->Eta();
143	Double_t phi = muonTrack->Phi();
144
145	// fill global histograms
146	hNumberOfClusters->Fill(nClusters);
147	hChi2PerDof->Fill(normalizedChi2);
148	hP->Fill(p);
149	hPt->Fill(pT);
150	hEta->Fill(eta);
151	hPhi->Fill(phi);
152
153	// try to match the reconstructed track with a simulated one
154	Float_t fractionOfMatchCluster = 0.;
155	AliMUONTrack* matchedTrackRef = MatchWithTrackRef(muonTrack, trackRefStore, fractionOfMatchCluster, useLabel, sigmaCut);
156
157	// take actions according to matching result
158	if (matchedTrackRef) {
159
160	// global counter
161	nTotMatchedTracks++;
162	if (!IsRecontructible(matchedTrackRef,recoParam)) nTotTracksReconstructedYet++;
163
164	// fill histograms
165	hFractionOfMatchedClusters->Fill(fractionOfMatchCluster);
166	hNumberOfClustersMC->Fill(matchedTrackRef->GetNClusters());
167	hNumberOfClustersM->Fill(nClusters);
168	hChi2PerDofM->Fill(normalizedChi2);
169	hPM->Fill(p);
170	hPtM->Fill(pT);
171	hEtaM->Fill(eta);
172	hPhiM->Fill(phi);
173
174	// remove already matched trackRefs
175	trackRefStore->Remove(*matchedTrackRef);
176
177	} else {
178
179	// global counter
180	nTotFakeTracks++;
181
182	// fill histograms
183	hNumberOfClustersF->Fill(nClusters);
184	hChi2PerDofF->Fill(normalizedChi2);
185	hPF->Fill(p);
186	hPtF->Fill(pT);
187	hEtaF->Fill(eta);
188	hPhiF->Fill(phi);
189
190	// store fake tracks
191	AliMUONESDInterface::Add(muonTrack, fakeTrackStore);
192
193	}
194
195	} // end of loop over ESD tracks
196
197	// fill histograms
198	hNumberOfTracks->Fill(nTrackerTracks);
199
200	// count the number the additional fake tracks
201	if (fakeTrackStore->GetSize() > 0) {
202
203	// remove the most connected fake tracks
204	Int_t nFreeMissingTracks = RemoveConnectedFakes(fakeTrackStore, trackRefStore, *recoParam,
205	useLabel, sigmaCut, *hFractionOfConnectedClusters);
206
207	// remove the remaining free reconstructible tracks
208	Int_t nAdditionalTracks = fakeTrackStore->GetSize() - nFreeMissingTracks;
209
210	// fill histograms
211	nEventsWithFake++;
212	if (nAdditionalTracks > 0) {
213	nEventsWithAdditionalFake++;
214	nTotAdditionalTracks += nAdditionalTracks;
215	hNumberOfAdditionalTracks->Fill(nAdditionalTracks);
216	}
217
218	}
219
220	delete fakeTrackStore;
221
222	} // end of loop over events
223
224	// plot results
225	TCanvas cFakesSummary("cFakesSummary","cFakesSummary",900,600);
226	cFakesSummary.Divide(3,2);
227	cFakesSummary.cd(1);
228	cFakesSummary.GetPad(1)->SetLogy();
229	hNumberOfClusters->Draw();
230	hNumberOfClusters->SetMinimum(0.5);
231	hNumberOfClustersM->Draw("same");
232	hNumberOfClustersM->SetLineColor(4);
233	hNumberOfClustersF->Draw("same");
234	hNumberOfClustersF->SetLineColor(2);
235	hNumberOfClustersF->SetFillColor(2);
236	hNumberOfClustersF->SetFillStyle(3017);
237	cFakesSummary.cd(2);
238	cFakesSummary.GetPad(2)->SetLogy();
239	hChi2PerDof->Draw();
240	hChi2PerDof->SetMinimum(0.5);
241	hChi2PerDofM->Draw("same");
242	hChi2PerDofM->SetLineColor(4);
243	hChi2PerDofF->Draw("same");
244	hChi2PerDofF->SetLineColor(2);
245	hChi2PerDofF->SetFillColor(2);
246	hChi2PerDofF->SetFillStyle(3017);
247	cFakesSummary.cd(3);
248	cFakesSummary.GetPad(3)->SetLogy();
249	hP->Draw();
250	hP->SetMinimum(0.5);
251	hPM->Draw("same");
252	hPM->SetLineColor(4);
253	hPF->Draw("same");
254	hPF->SetLineColor(2);
255	hPF->SetFillColor(2);
256	hPF->SetFillStyle(3017);
257	cFakesSummary.cd(4);
258	cFakesSummary.GetPad(4)->SetLogy();
259	hPt->Draw();
260	hPt->SetMinimum(0.5);
261	hPtM->Draw("same");
262	hPtM->SetLineColor(4);
263	hPtF->Draw("same");
264	hPtF->SetLineColor(2);
265	hPtF->SetFillColor(2);
266	hPtF->SetFillStyle(3017);
267	cFakesSummary.cd(5);
268	cFakesSummary.GetPad(5)->SetLogy();
269	hEta->Draw();
270	hEta->SetMinimum(0.5);
271	hEtaM->Draw("same");
272	hEtaM->SetLineColor(4);
273	hEtaF->Draw("same");
274	hEtaF->SetLineColor(2);
275	hEtaF->SetFillColor(2);
276	hEtaF->SetFillStyle(3017);
277	cFakesSummary.cd(6);
278	cFakesSummary.GetPad(6)->SetLogy();
279	hPhi->Draw();
280	hPhi->SetMinimum(0.5);
281	hPhiM->Draw("same");
282	hPhiM->SetLineColor(4);
283	hPhiF->Draw("same");
284	hPhiF->SetLineColor(2);
285	hPhiF->SetFillColor(2);
286	hPhiF->SetFillStyle(3017);
287
288	// save results
289	histoFile->cd();
290	histoFile->Write();
291	cFakesSummary.Write();
292	histoFile->Close();
293
294	// print results
295	cout << endl;
296	cout << "- Number of matched tracks: " << nTotMatchedTracks << endl;
297	cout << " (including " << nTotTracksReconstructedYet << " tracks matched with a TrackRef that is not reconstructible)" << endl;
298	cout << "- Number of fake tracks: " << nTotFakeTracks << endl;
299	cout << " (including " << nTotAdditionalTracks << " additional tracks (compared to the number of expected ones))" << endl;
300	cout << "- Number of events with fake track(s): " << nEventsWithFake << endl;
301	cout << " (including " << nEventsWithAdditionalFake << " events with additional tracks)" << endl;
302	cout << endl;
303	cout << "REMINDER: results are relevent provided that you use the same recoParams as for the reconstruction" << endl;
304	cout << endl;
305
306	// clear memory
307	delete esd;
308	esdFile->Close();
309	delete recoParam;
310
311	}
312
313	//-----------------------------------------------------------------------
314	void Prepare(AliMUONRecoParam *&recoParam, Double_t &sigmaCut)
315	{
316	/// Set the magnetic field and return recoParam and sigmaCut to associate cluster and trackRef
317
318	// prepare OCDB access
319	AliCDBManager* man = AliCDBManager::Instance();
320	man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
321	man->SetRun(0);
322
323	// set mag field
324	// waiting for mag field in CDB
325	if (!TGeoGlobalMagField::Instance()->GetField()) {
326	printf("Loading field map...\n");
327	AliMagF* field = new AliMagF("Maps","Maps",2,1.,1., 10.,AliMagF::k5kG);
328	TGeoGlobalMagField::Instance()->SetField(field);
329	}
330	// set the magnetic field for track extrapolations
331	AliMUONTrackExtrap::SetField();
332
333	// Load initial reconstruction parameters from OCDB
334	AliCDBPath path("MUON","Calib","RecoParam");
335	AliCDBEntry *entry=man->Get(path.GetPath());
336	if(entry) {
337	recoParam = dynamic_cast<AliMUONRecoParam*>(entry->GetObject());
338	entry->SetOwner(0);
339	AliCDBManager::Instance()->UnloadFromCache(path.GetPath());
340	}
341	if (!recoParam) {
342	printf("Couldn't find RecoParam object in OCDB: create default one");
343	recoParam = AliMUONRecoParam::GetLowFluxParam();
344	}
345
346	Info("MUONFakes", "\n recontruction parameters:");
347	recoParam->Print("FULL");
348	AliMUONESDInterface::ResetTracker(recoParam);
349
350	// sigma cut to associate clusters with TrackRefs in case the label are not used
351	sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking();
352
353	}
354
355	//-----------------------------------------------------------------------
356	TTree* GetESDTree(TFile *esdFile)
357	{
358	/// Check that the file is properly open
359	/// Return pointer to the ESD Tree
360
361	if (!esdFile \|\| !esdFile->IsOpen()) {
362	Error("GetESDTree", "opening ESD file failed");
363	exit(-1);
364	}
365
366	TTree* tree = (TTree*) esdFile->Get("esdTree");
367	if (!tree) {
368	Error("GetESDTree", "no ESD tree found");
369	exit(-1);
370	}
371
372	return tree;
373
374	}
375
376	//-----------------------------------------------------------------------
377	Bool_t TrackMatched(AliMUONTrack &track, AliMUONTrack &trackRef, Float_t &fractionOfMatchCluster, Double_t sigmaCut)
378	{
379	/// Try to match 2 tracks
380
381	Bool_t compTrack[10];
382	Int_t nMatchClusters = track.CompatibleTrack(&trackRef, sigmaCut, compTrack);
383	fractionOfMatchCluster = ((Float_t)nMatchClusters) / ((Float_t)track.GetNClusters());
384
385	if ((compTrack[0] \|\| compTrack[1] \|\| compTrack[2] \|\| compTrack[3]) && // at least 1 cluster matched in st 1 & 2
386	(compTrack[6] \|\| compTrack[7] \|\| compTrack[8] \|\| compTrack[9]) && // at least 1 cluster matched in st 4 & 5
387	fractionOfMatchCluster > 0.5) return kTRUE; // more than 50% of clusters matched
388	else return kFALSE;
389
390	}
391
392	//-----------------------------------------------------------------------
393	Bool_t IsRecontructible(AliMUONTrack &track, AliMUONRecoParam &recoParam)
394	{
395	/// Check il the track is reconstructible
396	Int_t nMinChHitInSt45 = (recoParam.MakeMoreTrackCandidates()) ? 2 : 3;
397	Int_t currentCh, previousCh = -1, nChHitInSt45 = 0;
398	Bool_t clusterInSt[5];
399	for (Int_t iSt = 0; iSt < 5; iSt++) clusterInSt[iSt] = !recoParam.RequestStation(iSt);
400
401	AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->First());
402	while (trackParam) {
403
404	currentCh = trackParam->GetClusterPtr()->GetChamberId();
405
406	clusterInSt[currentCh/2] = kTRUE;
407
408	if (currentCh > 5 && currentCh != previousCh) {
409	nChHitInSt45++;
410	previousCh = currentCh;
411	}
412
413	trackParam = static_cast<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->After(trackParam));
414	}
415
416	return (clusterInSt[0] && clusterInSt[1] && clusterInSt[2] &&
417	clusterInSt[3] && clusterInSt[4] && nChHitInSt45 >= nMinChHitInSt45);
418
419	}
420
421	//-----------------------------------------------------------------------
422	AliMUONTrack* MatchWithTrackRef(AliESDMuonTrack &muonTrack, AliMUONVTrackStore &trackRefStore,
423	Float_t &fractionOfMatchCluster, Bool_t useLabel, Double_t sigmaCut)
424	{
425	/// Return if the trackRef matched with the reconstructed track and the fraction of matched clusters
426
427	AliMUONTrack *matchedTrackRef = 0x0;
428	fractionOfMatchCluster = 0.;
429
430	if (useLabel) { // by using the MC label
431
432	// get the corresponding simulated track if any
433	Int_t label = muonTrack.GetLabel();
434	matchedTrackRef = (AliMUONTrack*) trackRefStore.FindObject(label);
435
436	// get the fraction of matched clusters
437	if (matchedTrackRef) {
438	Int_t nMatchClusters = 0;
439	if (muonTrack.ClustersStored()) {
440	AliESDMuonCluster* cluster = (AliESDMuonCluster*) muonTrack.GetClusters().First();
441	while (cluster) {
442	if (cluster->GetLabel() == label) nMatchClusters++;
443	cluster = (AliESDMuonCluster*) muonTrack.GetClusters().After(cluster);
444	}
445	}
446	fractionOfMatchCluster = ((Float_t)nMatchClusters) / ((Float_t)muonTrack.GetNClusters());
447	}
448
449	} else { // by comparing cluster/TrackRef positions
450
451	// convert ESD track to MUON track
452	AliMUONTrack track;
453	AliMUONESDInterface::ESDToMUON(muonTrack,track);
454
455	// look for the corresponding simulated track if any
456	TIter next(trackRefStore.CreateIterator());
457	AliMUONTrack* trackRef;
458	while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
459
460	// check compatibility
461	Float_t f = 0.;
462	if (TrackMatched(track, *trackRef, f, sigmaCut)) {
463	matchedTrackRef = trackRef;
464	fractionOfMatchCluster = f;
465	break;
466	}
467
468	}
469
470	}
471
472	return matchedTrackRef;
473
474	}
475
476	//-----------------------------------------------------------------------
477	Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore, AliMUONRecoParam &recoParam,
478	Bool_t useLabel, Double_t sigmaCut, TH1F &hFractionOfConnectedClusters)
479	{
480	/// loop over reconstructible TrackRef not associated with reconstructed track:
481	/// for each of them, find and remove the most connected the fake track, if any,
482	/// and fill the histograms with the fraction of connected clusters.
483	/// Return the number of reconstructible track not connected to any fake
484
485	Int_t nFreeMissingTracks = 0;
486
487	// loop over trackRefs
488	TIter next(trackRefStore.CreateIterator());
489	AliMUONTrack* trackRef;
490	while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
491
492	// skip not reconstructible trackRefs
493	if (!IsRecontructible(*trackRef,recoParam)) continue;
494
495	Int_t label = trackRef->GetUniqueID();
496
497	// look for the most connected fake track
498	AliMUONTrack *connectedFake = 0x0;
499	Float_t fractionOfConnectedClusters = 0.;
500	TIter next2(fakeTrackStore.CreateIterator());
501	AliMUONTrack* fakeTrack;
502	while ( ( fakeTrack = static_cast<AliMUONTrack*>(next2()) ) ) {
503
504	// get the number of connected clusters
505	Int_t nConnectedClusters = 0;
506	if (useLabel) { // by using the MC label
507	for (Int_t iCl = 0; iCl < fakeTrack->GetNClusters(); iCl++)
508	if (((AliMUONTrackParam*) fakeTrack->GetTrackParamAtCluster()->UncheckedAt(iCl))->GetClusterPtr()->GetMCLabel() == label)
509	nConnectedClusters++;
510	} else { // by comparing cluster/TrackRef positions
511	Bool_t compTrack[10];
512	nConnectedClusters = fakeTrack->CompatibleTrack(trackRef, sigmaCut, compTrack);
513	}
514
515	// skip non-connected fake tracks
516	if (nConnectedClusters == 0) continue;
517
518	// check if it is the most connected fake track
519	Float_t f = ((Float_t)nConnectedClusters) / ((Float_t)fakeTrack->GetNClusters());
520	if (f > fractionOfConnectedClusters) {
521	connectedFake = fakeTrack;
522	fractionOfConnectedClusters = f;
523	}
524
525	}
526
527	// remove the most connected fake track
528	if (connectedFake) {
529	hFractionOfConnectedClusters.Fill(fractionOfConnectedClusters);
530	fakeTrackStore.Remove(*connectedFake);
531	} else nFreeMissingTracks++;
532
533	}
534
535	return nFreeMissingTracks;
536
537	}
538