[u/mrichter/AliRoot.git] / MUON / DIMUONFakes.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>
#include <TLorentzVector.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 DIMUONFakes.C
///
/// \author Ph. Pillot, Subatech, March. 2009
///
/// Macro to study the effects of fake tracks on the dimuon spectra
/// Results are saved in the root file DiFakes.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);
AliMUONTrack* MatchWithTrackRef(AliESDMuonTrack &muonTrack, AliMUONVTrackStore &trackRefStore,
				Float_t &fractionOfMatchCluster, Bool_t useLabel, Double_t sigmaCut);

//-----------------------------------------------------------------------
void DIMUONFakes(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("DiFakes.root", "RECREATE");
  
  TH1F *hMass = new TH1F("hMass", "Muon mass distribution (GeV/c^{2})", 100, 0., 12.);
  TH1F *hMassM = new TH1F("hMassM", "matched track mass distribution (GeV/c^{2})", 100, 0., 12.);
  TH1F *hMassF = new TH1F("hMassF", "fake track mass distribution (GeV/c^{2})", 100, 0., 12.);
  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 *hY = new TH1F("hY"," Muon rapidity distribution",100,-10,0);
  TH1F *hYM = new TH1F("hYM"," matched track rapidity distribution",100,-10,0);
  TH1F *hYF = new TH1F("hYF"," fake track rapidity distribution",100,-10,0);
  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);
  
  TLorentzVector vMu1, vMu2, vDiMu;
  
  // 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 and flag them
    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;
      
      // 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) {
	
	// flag matched tracks
	muonTrack->SetLabel(matchedTrackRef->GetUniqueID());
	
	// remove already matched trackRefs
	trackRefStore->Remove(*matchedTrackRef);
	
      } else {
	
	// flag fake tracks
	muonTrack->SetLabel(-1);
	
      }
      
    }
    
    // double loop over ESD tracks, build pairs and fill histograms according to their label
    for (Int_t iTrack1 = 0; iTrack1 <  nTracks;  iTrack1++) {
      AliESDMuonTrack* muonTrack1 = esd->GetMuonTrack(iTrack1);
      
      // skip ghosts
      if (!muonTrack1->ContainTrackerData()) continue;
      
      // get track info
      Short_t charge1 = muonTrack1->Charge();
      Int_t label1 = muonTrack1->GetLabel();
      muonTrack1->LorentzP(vMu1);
      
      for (Int_t iTrack2 = iTrack1+1; iTrack2 <  nTracks;  iTrack2++) {
	AliESDMuonTrack* muonTrack2 = esd->GetMuonTrack(iTrack2);
	
	// skip ghosts
	if (!muonTrack2->ContainTrackerData()) continue;
	
	// keep only opposite sign pairs
	Short_t charge2 = muonTrack2->Charge();
	if (charge1*charge2 > 0) continue;
	
	// get track info
	Int_t label2 = muonTrack2->GetLabel();
	muonTrack2->LorentzP(vMu2);
	
	// compute kinematics of the pair
	vDiMu = vMu1 + vMu2;
	Float_t mass = vDiMu.M();
	Float_t p = vDiMu.P();
	Float_t pt = vDiMu.Pt();
	Float_t y = vDiMu.Rapidity();
	Float_t eta = vDiMu.Eta();
	Float_t phi = vMu1.Phi();
	if (phi < 0) phi += 2.*TMath::Pi();
	
	// fill global histograms
	hMass->Fill(mass);
	hP->Fill(p);
	hPt->Fill(pt);
	hY->Fill(y);
	hEta->Fill(eta);
	hPhi->Fill(phi);
	
	// fill histograms according to labels
	if (label1 >= 0 && label2 >= 0) {
	  
	  hMassM->Fill(mass);
	  hPM->Fill(p);
	  hPtM->Fill(pt);
	  hYM->Fill(y);
	  hEtaM->Fill(eta);
	  hPhiM->Fill(phi);
	  
	} else {
	  
	  hMassF->Fill(mass);
	  hPF->Fill(p);
	  hPtF->Fill(pt);
	  hYF->Fill(y);
	  hEtaF->Fill(eta);
	  hPhiF->Fill(phi);
	  
	}
	
      }
      
    }
      
  } // end of loop over events
  
  // plot results
  TCanvas cDiFakesSummary("cDiFakesSummary","cDiFakesSummary",900,600);
  cDiFakesSummary.Divide(3,2);
  cDiFakesSummary.cd(1);
  cDiFakesSummary.GetPad(1)->SetLogy();
  hMass->Draw();
  hMass->SetMinimum(0.5);
  hMassM->Draw("same");
  hMassM->SetLineColor(4);
  hMassF->Draw("same");
  hMassF->SetLineColor(2);
  hMassF->SetFillColor(2);
  hMassF->SetFillStyle(3017);
  cDiFakesSummary.cd(2);
  cDiFakesSummary.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);
  cDiFakesSummary.cd(3);
  cDiFakesSummary.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);
  cDiFakesSummary.cd(4);
  cDiFakesSummary.GetPad(2)->SetLogy();
  hY->Draw();
  hY->SetMinimum(0.5);
  hYM->Draw("same");
  hYM->SetLineColor(4);
  hYF->Draw("same");
  hYF->SetLineColor(2);
  hYF->SetFillColor(2);
  hYF->SetFillStyle(3017);
  cDiFakesSummary.cd(5);
  cDiFakesSummary.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);
  cDiFakesSummary.cd(6);
  cDiFakesSummary.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();
  cDiFakesSummary.Write();
  histoFile->Close();
  
  // 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;
  
}

//-----------------------------------------------------------------------
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;
  
}
Commit	Line	Data
30b3fe0f	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	#include <TLorentzVector.h>
	12
	13	// STEER includes
	14	#include "AliLog.h"
	15	#include "AliMagF.h"
	16	#include "AliESDEvent.h"
	17	#include "AliESDMuonTrack.h"
	18	#include "AliESDMuonCluster.h"
	19	#include "AliCDBPath.h"
	20	#include "AliCDBEntry.h"
	21	#include "AliCDBManager.h"
	22
	23	// MUON includes
	24	#include "AliMUONTrack.h"
	25	#include "AliMUONVTrackStore.h"
	26	#include "AliMUONTrackParam.h"
	27	#include "AliMUONTrackExtrap.h"
	28	#include "AliMUONESDInterface.h"
	29	#include "AliMUONRecoCheck.h"
	30	#include "AliMUONVCluster.h"
	31	#include "AliMUONRecoParam.h"
	32	#endif
	33
	34	/// \ingroup macros
	35	/// \file DIMUONFakes.C
	36	///
	37	/// \author Ph. Pillot, Subatech, March. 2009
	38	///
	39	/// Macro to study the effects of fake tracks on the dimuon spectra
	40	/// Results are saved in the root file DiFakes.root
	41	/// Results are relevent provided that you use the same recoParams as for the reconstruction
	42
	43	void Prepare(AliMUONRecoParam *&recoParam, Double_t &sigmaCut);
	44	TTree* GetESDTree(TFile *esdFile);
	45	Bool_t TrackMatched(AliMUONTrack &track, AliMUONTrack &trackRef, Float_t &fractionOfMatchCluster, Double_t sigmaCut);
	46	AliMUONTrack* MatchWithTrackRef(AliESDMuonTrack &muonTrack, AliMUONVTrackStore &trackRefStore,
	47	Float_t &fractionOfMatchCluster, Bool_t useLabel, Double_t sigmaCut);
	48
	49	//-----------------------------------------------------------------------
	50	void DIMUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent = -1,
	51	const TString esdFileName = "AliESDs.root", const TString SimDir = "./generated/")
	52	{
	53
	54	//Reset ROOT and connect tree file
	55	gROOT->Reset();
	56
	57	// File for histograms and histogram booking
	58	TFile *histoFile = new TFile("DiFakes.root", "RECREATE");
	59
	60	TH1F *hMass = new TH1F("hMass", "Muon mass distribution (GeV/c^{2})", 100, 0., 12.);
	61	TH1F *hMassM = new TH1F("hMassM", "matched track mass distribution (GeV/c^{2})", 100, 0., 12.);
	62	TH1F *hMassF = new TH1F("hMassF", "fake track mass distribution (GeV/c^{2})", 100, 0., 12.);
	63	TH1F *hP = new TH1F("hP", "Muon P distribution (GeV/c)", 100, 0., 200.);
	64	TH1F *hPM = new TH1F("hPM", "matched track P distribution (GeV/c)", 100, 0., 200.);
65	TH1F *hPF = new TH1F("hPF", "fake track P distribution (GeV/c)", 100, 0., 200.);
66	TH1F *hPt = new TH1F("hPt", "Muon Pt distribution (GeV/c)", 100, 0., 20.);
67	TH1F *hPtM = new TH1F("hPtM", "matched track Pt distribution (GeV/c)", 100, 0., 20.);
68	TH1F *hPtF = new TH1F("hPtF", "fake track Pt distribution (GeV/c)", 100, 0., 20.);
69	TH1F *hY = new TH1F("hY"," Muon rapidity distribution",100,-10,0);
70	TH1F *hYM = new TH1F("hYM"," matched track rapidity distribution",100,-10,0);
71	TH1F *hYF = new TH1F("hYF"," fake track rapidity distribution",100,-10,0);
72	TH1F *hEta = new TH1F("hEta"," Muon pseudo-rapidity distribution",100,-10,0);
73	TH1F *hEtaM = new TH1F("hEtaM"," matched track pseudo-rapidity distribution",100,-10,0);
74	TH1F *hEtaF = new TH1F("hEtaF"," fake track pseudo-rapidity distribution",100,-10,0);
75	TH1F *hPhi = new TH1F("hPhi"," Muon phi distribution",100,-1.,9.);
76	TH1F *hPhiM = new TH1F("hPhiM"," matched track phi distribution",100,-1.,9.);
77	TH1F *hPhiF = new TH1F("hPhiF"," fake track phi distribution",100,-1.,9.);
78
79	// prepare for analysis
80	AliMUONRecoParam *recoParam = 0x0;
81	Double_t sigmaCut = -1;
82	Prepare(recoParam, sigmaCut);
83
84	// link to reconstructed tracks
85	TFile* esdFile = TFile::Open(esdFileName);
86	TTree* esdTree = GetESDTree(esdFile);
87	AliESDEvent* esd = new AliESDEvent();
88	esd->ReadFromTree(esdTree);
89
90	// link to simulated tracks
91	AliMUONRecoCheck rc(esdFileName, SimDir);
92
93	TLorentzVector vMu1, vMu2, vDiMu;
94
95	// Loop over ESD events
96	FirstEvent = TMath::Max(0, FirstEvent);
97	LastEvent = (LastEvent>=0) ? TMath::Min((Int_t)esdTree->GetEntries() - 1, LastEvent) : (Int_t)esdTree->GetEntries() - 1;
98	for (Int_t iEvent = FirstEvent; iEvent <= LastEvent; iEvent++) {
99
100	// get the ESD of current event
101	esdTree->GetEvent(iEvent);
102	if (!esd) {
103	Error("CheckESD", "no ESD object found for event %d", iEvent);
104	return;
105	}
106
107	// convert TrackRef to MUON tracks
108	AliMUONVTrackStore* trackRefStore = rc.TrackRefs(iEvent);
109
110	// loop over ESD tracks and flag them
111	Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
112	for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
113
114	AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
115
116	// skip ghosts
117	if (!muonTrack->ContainTrackerData()) continue;
118
119	// try to match the reconstructed track with a simulated one
120	Float_t fractionOfMatchCluster = 0.;
121	AliMUONTrack* matchedTrackRef = MatchWithTrackRef(muonTrack, trackRefStore, fractionOfMatchCluster, useLabel, sigmaCut);
122
123	// take actions according to matching result
124	if (matchedTrackRef) {
125
126	// flag matched tracks
127	muonTrack->SetLabel(matchedTrackRef->GetUniqueID());
128
129	// remove already matched trackRefs
130	trackRefStore->Remove(*matchedTrackRef);
131
132	} else {
133
134	// flag fake tracks
135	muonTrack->SetLabel(-1);
136
137	}
138
139	}
140
141	// double loop over ESD tracks, build pairs and fill histograms according to their label
142	for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
143	AliESDMuonTrack* muonTrack1 = esd->GetMuonTrack(iTrack1);
144
145	// skip ghosts
146	if (!muonTrack1->ContainTrackerData()) continue;
147
148	// get track info
149	Short_t charge1 = muonTrack1->Charge();
150	Int_t label1 = muonTrack1->GetLabel();
151	muonTrack1->LorentzP(vMu1);
152
153	for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
154	AliESDMuonTrack* muonTrack2 = esd->GetMuonTrack(iTrack2);
155
156	// skip ghosts
157	if (!muonTrack2->ContainTrackerData()) continue;
158
159	// keep only opposite sign pairs
160	Short_t charge2 = muonTrack2->Charge();
161	if (charge1*charge2 > 0) continue;
162
163	// get track info
164	Int_t label2 = muonTrack2->GetLabel();
165	muonTrack2->LorentzP(vMu2);
166
167	// compute kinematics of the pair
168	vDiMu = vMu1 + vMu2;
169	Float_t mass = vDiMu.M();
170	Float_t p = vDiMu.P();
171	Float_t pt = vDiMu.Pt();
172	Float_t y = vDiMu.Rapidity();
173	Float_t eta = vDiMu.Eta();
174	Float_t phi = vMu1.Phi();
175	if (phi < 0) phi += 2.*TMath::Pi();
176
177	// fill global histograms
178	hMass->Fill(mass);
179	hP->Fill(p);
180	hPt->Fill(pt);
181	hY->Fill(y);
182	hEta->Fill(eta);
183	hPhi->Fill(phi);
184
185	// fill histograms according to labels
186	if (label1 >= 0 && label2 >= 0) {
187
188	hMassM->Fill(mass);
189	hPM->Fill(p);
190	hPtM->Fill(pt);
191	hYM->Fill(y);
192	hEtaM->Fill(eta);
193	hPhiM->Fill(phi);
194
195	} else {
196
197	hMassF->Fill(mass);
198	hPF->Fill(p);
199	hPtF->Fill(pt);
200	hYF->Fill(y);
201	hEtaF->Fill(eta);
202	hPhiF->Fill(phi);
203
204	}
205
206	}
207
208	}
209
210	} // end of loop over events
211
212	// plot results
213	TCanvas cDiFakesSummary("cDiFakesSummary","cDiFakesSummary",900,600);
214	cDiFakesSummary.Divide(3,2);
215	cDiFakesSummary.cd(1);
216	cDiFakesSummary.GetPad(1)->SetLogy();
217	hMass->Draw();
218	hMass->SetMinimum(0.5);
219	hMassM->Draw("same");
220	hMassM->SetLineColor(4);
221	hMassF->Draw("same");
222	hMassF->SetLineColor(2);
223	hMassF->SetFillColor(2);
224	hMassF->SetFillStyle(3017);
225	cDiFakesSummary.cd(2);
226	cDiFakesSummary.GetPad(3)->SetLogy();
227	hP->Draw();
228	hP->SetMinimum(0.5);
229	hPM->Draw("same");
230	hPM->SetLineColor(4);
231	hPF->Draw("same");
232	hPF->SetLineColor(2);
233	hPF->SetFillColor(2);
234	hPF->SetFillStyle(3017);
235	cDiFakesSummary.cd(3);
236	cDiFakesSummary.GetPad(4)->SetLogy();
237	hPt->Draw();
238	hPt->SetMinimum(0.5);
239	hPtM->Draw("same");
240	hPtM->SetLineColor(4);
241	hPtF->Draw("same");
242	hPtF->SetLineColor(2);
243	hPtF->SetFillColor(2);
244	hPtF->SetFillStyle(3017);
245	cDiFakesSummary.cd(4);
246	cDiFakesSummary.GetPad(2)->SetLogy();
247	hY->Draw();
248	hY->SetMinimum(0.5);
249	hYM->Draw("same");
250	hYM->SetLineColor(4);
251	hYF->Draw("same");
252	hYF->SetLineColor(2);
253	hYF->SetFillColor(2);
254	hYF->SetFillStyle(3017);
255	cDiFakesSummary.cd(5);
256	cDiFakesSummary.GetPad(5)->SetLogy();
257	hEta->Draw();
258	hEta->SetMinimum(0.5);
259	hEtaM->Draw("same");
260	hEtaM->SetLineColor(4);
261	hEtaF->Draw("same");
262	hEtaF->SetLineColor(2);
263	hEtaF->SetFillColor(2);
264	hEtaF->SetFillStyle(3017);
265	cDiFakesSummary.cd(6);
266	cDiFakesSummary.GetPad(6)->SetLogy();
267	hPhi->Draw();
268	hPhi->SetMinimum(0.5);
269	hPhiM->Draw("same");
270	hPhiM->SetLineColor(4);
271	hPhiF->Draw("same");
272	hPhiF->SetLineColor(2);
273	hPhiF->SetFillColor(2);
274	hPhiF->SetFillStyle(3017);
275
276	// save results
277	histoFile->cd();
278	histoFile->Write();
279	cDiFakesSummary.Write();
280	histoFile->Close();
281
282	// clear memory
283	delete esd;
284	esdFile->Close();
285	delete recoParam;
286
287	}
288
289	//-----------------------------------------------------------------------
290	void Prepare(AliMUONRecoParam *&recoParam, Double_t &sigmaCut)
291	{
292	/// Set the magnetic field and return recoParam and sigmaCut to associate cluster and trackRef
293
294	// prepare OCDB access
295	AliCDBManager* man = AliCDBManager::Instance();
296	man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
297	man->SetRun(0);
298
299	// set mag field
300	// waiting for mag field in CDB
301	if (!TGeoGlobalMagField::Instance()->GetField()) {
302	printf("Loading field map...\n");
303	AliMagF* field = new AliMagF("Maps","Maps",2,1.,1., 10.,AliMagF::k5kG);
304	TGeoGlobalMagField::Instance()->SetField(field);
305	}
306	// set the magnetic field for track extrapolations
307	AliMUONTrackExtrap::SetField();
308
309	// Load initial reconstruction parameters from OCDB
310	AliCDBPath path("MUON","Calib","RecoParam");
311	AliCDBEntry *entry=man->Get(path.GetPath());
312	if(entry) {
313	recoParam = dynamic_cast<AliMUONRecoParam*>(entry->GetObject());
314	entry->SetOwner(0);
315	AliCDBManager::Instance()->UnloadFromCache(path.GetPath());
316	}
317	if (!recoParam) {
318	printf("Couldn't find RecoParam object in OCDB: create default one");
319	recoParam = AliMUONRecoParam::GetLowFluxParam();
320	}
321
322	Info("MUONFakes", "\n recontruction parameters:");
323	recoParam->Print("FULL");
324	AliMUONESDInterface::ResetTracker(recoParam);
325
326	// sigma cut to associate clusters with TrackRefs in case the label are not used
327	sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking();
328
329	}
330
331	//-----------------------------------------------------------------------
332	TTree* GetESDTree(TFile *esdFile)
333	{
334	/// Check that the file is properly open
335	/// Return pointer to the ESD Tree
336
337	if (!esdFile \|\| !esdFile->IsOpen()) {
338	Error("GetESDTree", "opening ESD file failed");
339	exit(-1);
340	}
341
342	TTree* tree = (TTree*) esdFile->Get("esdTree");
343	if (!tree) {
344	Error("GetESDTree", "no ESD tree found");
345	exit(-1);
346	}
347
348	return tree;
349
350	}
351
352	//-----------------------------------------------------------------------
353	Bool_t TrackMatched(AliMUONTrack &track, AliMUONTrack &trackRef, Float_t &fractionOfMatchCluster, Double_t sigmaCut)
354	{
355	/// Try to match 2 tracks
356
357	Bool_t compTrack[10];
358	Int_t nMatchClusters = track.CompatibleTrack(&trackRef, sigmaCut, compTrack);
359	fractionOfMatchCluster = ((Float_t)nMatchClusters) / ((Float_t)track.GetNClusters());
360
361	if ((compTrack[0] \|\| compTrack[1] \|\| compTrack[2] \|\| compTrack[3]) && // at least 1 cluster matched in st 1 & 2
362	(compTrack[6] \|\| compTrack[7] \|\| compTrack[8] \|\| compTrack[9]) && // at least 1 cluster matched in st 4 & 5
363	fractionOfMatchCluster > 0.5) return kTRUE; // more than 50% of clusters matched
364	else return kFALSE;
365
366	}
367
368	//-----------------------------------------------------------------------
369	AliMUONTrack* MatchWithTrackRef(AliESDMuonTrack &muonTrack, AliMUONVTrackStore &trackRefStore,
370	Float_t &fractionOfMatchCluster, Bool_t useLabel, Double_t sigmaCut)
371	{
372	/// Return if the trackRef matched with the reconstructed track and the fraction of matched clusters
373
374	AliMUONTrack *matchedTrackRef = 0x0;
375	fractionOfMatchCluster = 0.;
376
377	if (useLabel) { // by using the MC label
378
379	// get the corresponding simulated track if any
380	Int_t label = muonTrack.GetLabel();
381	matchedTrackRef = (AliMUONTrack*) trackRefStore.FindObject(label);
382
383	// get the fraction of matched clusters
384	if (matchedTrackRef) {
385	Int_t nMatchClusters = 0;
386	if (muonTrack.ClustersStored()) {
387	AliESDMuonCluster* cluster = (AliESDMuonCluster*) muonTrack.GetClusters().First();
388	while (cluster) {
389	if (cluster->GetLabel() == label) nMatchClusters++;
390	cluster = (AliESDMuonCluster*) muonTrack.GetClusters().After(cluster);
391	}
392	}
393	fractionOfMatchCluster = ((Float_t)nMatchClusters) / ((Float_t)muonTrack.GetNClusters());
394	}
395
396	} else { // by comparing cluster/TrackRef positions
397
398	// convert ESD track to MUON track
399	AliMUONTrack track;
400	AliMUONESDInterface::ESDToMUON(muonTrack,track);
401
402	// look for the corresponding simulated track if any
403	TIter next(trackRefStore.CreateIterator());
404	AliMUONTrack* trackRef;
405	while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
406
407	// check compatibility
408	Float_t f = 0.;
409	if (TrackMatched(track, *trackRef, f, sigmaCut)) {
410	matchedTrackRef = trackRef;
411	fractionOfMatchCluster = f;
412	break;
413	}
414
415	}
416
417	}
418
419	return matchedTrackRef;
420
421	}
422