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

#if !defined(__CINT__) || defined(__MAKECINT__)
// ROOT includes
#include <TFile.h>
#include <TH1.h>
#include <TCanvas.h>
#include <Riostream.h>
#include <TROOT.h>
#include <TObjArray.h>
#include <TArrayI.h>

// STEER includes
#include "AliLog.h"
#include "AliESDEvent.h"
#include "AliESDMuonTrack.h"
#include "AliCDBManager.h"
  
// MUON includes
#include "AliMUONCDB.h"
#include "AliMUONTrack.h"
#include "AliMUONVTrackStore.h"
#include "AliMUONTrackParam.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

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 ocdbPath = "local://$ALICE_ROOT/OCDB")
{
  
  //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.);
  
  // link to reconstructed and simulated tracks
  AliMUONRecoCheck rc(esdFileName, SimDir);
  
  // load necessary data from OCDB
  AliCDBManager::Instance()->SetDefaultStorage(ocdbPath);
  AliCDBManager::Instance()->SetRun(rc.GetRunNumber());
  AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
  if (!recoParam) return;
  
  // 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;
  Int_t nReconstructedTracks = 0;
  Int_t nEventsWithTrackReconstructedYet = 0;
  Int_t nEventsWithFake = 0;
  Int_t nEventsWithAdditionalFake = 0;
  Int_t nTotMatchedTracks = 0;
  Int_t nTotTracksReconstructedYet = 0;
  Int_t nTotFakeTracks = 0;
  Int_t nTotConnectedTracks = 0;
  Int_t nTotAdditionalTracks = 0;
  Bool_t trackReconstructedYet;
  TArrayI eventsWithTrackReconstructedYet(10);
  TArrayI eventsWithFake(10);
  TArrayI eventsWithAdditionalFake(10);
  
  // Loop over ESD events
  FirstEvent = TMath::Max(0, FirstEvent);
  LastEvent = (LastEvent>=0) ? TMath::Min(rc.NumberOfEvents() - 1, LastEvent) : rc.NumberOfEvents() - 1;
  for (Int_t iEvent = FirstEvent; iEvent <= LastEvent; iEvent++) {
    
    // 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 (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* esdTrack = esd->GetMuonTrack(iTrack);
      
      // skip ghosts
      if (!esdTrack->ContainTrackerData()) continue;
      nTrackerTracks++;
      
      // find the corresponding MUON track
      AliMUONTrack* muonTrack = (AliMUONTrack*) muonTrackStore->FindObject(esdTrack->GetUniqueID());
      
      // get track info
      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);
      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
      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 (!matchedTrackRef->IsValid(requestedStationMask, request2ChInSameSt45)) {
	  trackReconstructedYet = kTRUE;
	  nTotTracksReconstructedYet++;
	}
	
	// fill histograms
	hFractionOfMatchedClusters->Fill(((Float_t) nMatchClusters) / ((Float_t) nClusters));
	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
	fakeTrackStore->Add(*muonTrack);
	
      }
      
    } // end of loop over ESD tracks
    
    // fill histograms
    hNumberOfTracks->Fill(nTrackerTracks);
    nReconstructedTracks += nTrackerTracks;
    if (trackReconstructedYet) eventsWithTrackReconstructedYet[nEventsWithTrackReconstructedYet++] = iEvent;
    
    // count the number the additional fake tracks
    if (fakeTrackStore->GetSize() > 0) {
      
      // remove the most connected fake tracks
      Int_t nFreeMissingTracks = RemoveConnectedFakes(*fakeTrackStore, *trackRefStore,  useLabel, *hFractionOfConnectedClusters);
      
      // remove the remaining free reconstructible tracks
      Int_t nAdditionalTracks = fakeTrackStore->GetSize() - nFreeMissingTracks;
      
      // fill histograms
      eventsWithFake[nEventsWithFake] = iEvent;
      nEventsWithFake++;
      if (nAdditionalTracks > 0) {
	eventsWithAdditionalFake[nEventsWithAdditionalFake] = iEvent;
	nEventsWithAdditionalFake++;
	nTotAdditionalTracks += nAdditionalTracks;
	hNumberOfAdditionalTracks->Fill(nAdditionalTracks);
      }
      
    }
    
    delete fakeTrackStore;
    
  } // end of loop over events

  // total number of connected tracks
  nTotConnectedTracks = hFractionOfConnectedClusters->GetEntries();
  
  // 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 reconstructible tracks: " << nReconstructibleTracks << endl;
  cout << "- Number of reconstructed tracks: " << nReconstructedTracks << endl;
  cout << "- Number of matched tracks: " << nTotMatchedTracks << endl;
  cout << "  (including " << nTotTracksReconstructedYet << " track(s) matched with a TrackRef that is not reconstructible";
  if (nTotTracksReconstructedYet > 0) {
    for(Int_t i=0; i<nEventsWithTrackReconstructedYet; i++){
      if (i==0) cout << " (eventID = " << eventsWithTrackReconstructedYet[i];
      else cout << ", " << eventsWithTrackReconstructedYet[i];
    }
    cout << "))" << endl;
  } else cout << ")" << endl;
  cout << "- Number of fake tracks: " << nTotFakeTracks << endl;
  cout << "  (including " << nTotConnectedTracks << " track(s) still connected to a reconstructible one)" << endl;
  cout << "  (including " << nTotAdditionalTracks << " additional track(s) (compared to the number of expected ones))" << endl;
  cout << "- Number of events with fake track(s): " << nEventsWithFake;
  if (nEventsWithFake > 0) {
    for(Int_t i=0; i<nEventsWithFake; i++){
      if (i==0) cout << " (eventID = " << eventsWithFake[i];
      else cout << ", " << eventsWithFake[i];
    }
    cout << ")" << endl;
  } else cout << endl;
  cout << "  (including " << nEventsWithAdditionalFake << " events with additional track(s)";
  if (nEventsWithAdditionalFake > 0) {
    for(Int_t i=0; i<nEventsWithAdditionalFake; i++){
      if (i==0) cout << " (eventID = " << eventsWithAdditionalFake[i];
      else cout << ", " << eventsWithAdditionalFake[i];
    }
    cout << "))" << endl;
  } else cout << ")" << endl;
  cout << endl;
  cout << "REMINDER: results are relevent provided that you use the same recoParams as for the reconstruction" << endl;
  cout << endl;
  
}

//-----------------------------------------------------------------------
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,
  /// 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 (!trackRef->IsValid(requestedStationMask, request2ChInSameSt45)) continue;
    
    Int_t label = trackRef->GetUniqueID();
    
    // look for the most connected fake track
    AliMUONTrack *connectedFake = 0x0;
    Double_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->FindCompatibleClusters(*trackRef, sigmaCut, compTrack);
      }
      
      // skip non-connected fake tracks
      if (nConnectedClusters == 0) continue;
      
      // check if it is the most connected fake track
      Double_t f = ((Double_t)nConnectedClusters) / ((Double_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__)
66085093	2	// ROOT includes
66085093	3	#include <TFile.h>
	4	#include <TH1.h>
	5	#include <TCanvas.h>
	6	#include <Riostream.h>
	7	#include <TROOT.h>
ce8cd162	8	#include <TObjArray.h>
23035434	9	#include <TArrayI.h>
66085093	10
	11	// STEER includes
	12	#include "AliLog.h"
66085093	13	#include "AliESDEvent.h"
66085093	14	#include "AliESDMuonTrack.h"
66085093	15	#include "AliCDBManager.h"
	16
	17	// MUON includes
a99c3449	18	#include "AliMUONCDB.h"
66085093	19	#include "AliMUONTrack.h"
	20	#include "AliMUONVTrackStore.h"
	21	#include "AliMUONTrackParam.h"
66085093	22	#include "AliMUONESDInterface.h"
	23	#include "AliMUONRecoCheck.h"
	24	#include "AliMUONVCluster.h"
	25	#include "AliMUONRecoParam.h"
	26	#endif
	27
	28	/// \ingroup macros
	29	/// \file MUONFakes.C
	30	///
	31	/// \author Ph. Pillot, Subatech, March. 2009
	32	///
	33	/// Macro to study fake tracks by comparing reconstructed tracks with TrackRefs
	34	/// Results are saved in the root file Fakes.root
	35	/// Results are relevent provided that you use the same recoParams as for the reconstruction
	36
f202486b	37	UInt_t requestedStationMask = 0;
	38	Bool_t request2ChInSameSt45 = kFALSE;
	39	Double_t sigmaCut = -1.;
	40
	41	Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore,
	42	Bool_t useLabel, TH1F &hFractionOfConnectedClusters);
66085093	43
	44	//-----------------------------------------------------------------------
	45	void MUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent = -1,
a99c3449	46	const TString esdFileName = "AliESDs.root", const TString SimDir = "./generated/",
a99c3449	47	const TString ocdbPath = "local://$ALICE_ROOT/OCDB")
66085093	48	{
	49
	50	//Reset ROOT and connect tree file
	51	gROOT->Reset();
	52
	53	// File for histograms and histogram booking
	54	TFile *histoFile = new TFile("Fakes.root", "RECREATE");
	55
	56	TH1F *hNumberOfTracks = new TH1F("hNumberOfTracks","nb of tracks /evt",20,0.,20.);
	57	TH1F *hNumberOfAdditionalTracks = new TH1F("hNumberOfAdditionalTracks","nb of fake - nb of missing track",20,0.,20.);
	58
	59	TH1F *hNumberOfClusters = new TH1F("hNumberOfClusters","nb of clusters /track",20,0.,20.);
	60	TH1F *hNumberOfClustersM = new TH1F("hNumberOfClustersM","nb of clusters /matched track",20,0.,20.);
	61	TH1F *hNumberOfClustersF = new TH1F("hNumberOfClustersF","nb of clusters /fake track",20,0.,20.);
	62	TH1F *hNumberOfClustersMC = new TH1F("hNumberOfClustersMC","nb of clusters /MC track",20,0.,20.);
	63	TH1F *hFractionOfMatchedClusters = new TH1F("hFractionOfMatchedClusters","nb of matched clusters / nb of clusters",110,0.,1.1);
	64	TH1F *hFractionOfConnectedClusters = new TH1F("hFractionOfConnectedClusters","nb of connected clusters / nb of clusters in fake tracks",110,0.,1.1);
	65
	66	TH1F *hChi2PerDof = new TH1F("hChi2PerDof", "track chi2/d.o.f.", 100, 0., 20.);
	67	TH1F *hChi2PerDofM = new TH1F("hChi2PerDofM", "matched track chi2/d.o.f.", 100, 0., 20.);
	68	TH1F *hChi2PerDofF = new TH1F("hChi2PerDofF", "fake track chi2/d.o.f.", 100, 0., 20.);
	69	TH1F *hP = new TH1F("hP", "Muon P distribution (GeV/c)", 100, 0., 200.);
	70	TH1F *hPM = new TH1F("hPM", "matched track P distribution (GeV/c)", 100, 0., 200.);
	71	TH1F *hPF = new TH1F("hPF", "fake track P distribution (GeV/c)", 100, 0., 200.);
	72	TH1F *hPt = new TH1F("hPt", "Muon Pt distribution (GeV/c)", 100, 0., 20.);
	73	TH1F *hPtM = new TH1F("hPtM", "matched track Pt distribution (GeV/c)", 100, 0., 20.);
	74	TH1F *hPtF = new TH1F("hPtF", "fake track Pt distribution (GeV/c)", 100, 0., 20.);
	75	TH1F *hEta = new TH1F("hEta"," Muon pseudo-rapidity distribution",100,-10,0);
	76	TH1F *hEtaM = new TH1F("hEtaM"," matched track pseudo-rapidity distribution",100,-10,0);
	77	TH1F *hEtaF = new TH1F("hEtaF"," fake track pseudo-rapidity distribution",100,-10,0);
	78	TH1F *hPhi = new TH1F("hPhi"," Muon phi distribution",100,-1.,9.);
	79	TH1F *hPhiM = new TH1F("hPhiM"," matched track phi distribution",100,-1.,9.);
	80	TH1F *hPhiF = new TH1F("hPhiF"," fake track phi distribution",100,-1.,9.);
	81
a99c3449	82	// link to reconstructed and simulated tracks
a99c3449	83	AliMUONRecoCheck rc(esdFileName, SimDir);
66085093	84
a99c3449	85	// load necessary data from OCDB
	86	AliCDBManager::Instance()->SetDefaultStorage(ocdbPath);
	87	AliCDBManager::Instance()->SetRun(rc.GetRunNumber());
	88	AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
	89	if (!recoParam) return;
66085093	90
a99c3449	91	// get sigma cut from recoParam to associate clusters with TrackRefs in case the label are not used
f202486b	92	sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking();
	93	// compute the mask of requested stations from recoParam
	94	for (Int_t i = 0; i < 5; i++) if (recoParam->RequestStation(i)) requestedStationMask \|= ( 1 << i );
	95	// get from recoParam whether a track need 2 chambers hit in the same station (4 or 5) or not to be reconstructible
	96	request2ChInSameSt45 = !recoParam->MakeMoreTrackCandidates();
66085093	97
66085093	98	// initialize global counters
23035434	99	Int_t nReconstructibleTracks = 0;
	100	Int_t nReconstructedTracks = 0;
	101	Int_t nEventsWithTrackReconstructedYet = 0;
66085093	102	Int_t nEventsWithFake = 0;
	103	Int_t nEventsWithAdditionalFake = 0;
	104	Int_t nTotMatchedTracks = 0;
	105	Int_t nTotTracksReconstructedYet = 0;
	106	Int_t nTotFakeTracks = 0;
23035434	107	Int_t nTotConnectedTracks = 0;
66085093	108	Int_t nTotAdditionalTracks = 0;
23035434	109	Bool_t trackReconstructedYet;
	110	TArrayI eventsWithTrackReconstructedYet(10);
	111	TArrayI eventsWithFake(10);
	112	TArrayI eventsWithAdditionalFake(10);
66085093	113
	114	// Loop over ESD events
	115	FirstEvent = TMath::Max(0, FirstEvent);
a99c3449	116	LastEvent = (LastEvent>=0) ? TMath::Min(rc.NumberOfEvents() - 1, LastEvent) : rc.NumberOfEvents() - 1;
66085093	117	for (Int_t iEvent = FirstEvent; iEvent <= LastEvent; iEvent++) {
66085093	118
a99c3449	119	// get reconstructed and simulated tracks
a99c3449	120	AliMUONVTrackStore* muonTrackStore = rc.ReconstructedTracks(iEvent, kFALSE);
66085093	121	AliMUONVTrackStore* trackRefStore = rc.TrackRefs(iEvent);
a99c3449	122	if (!muonTrackStore \|\| !trackRefStore) continue;
66085093	123
23035434	124	// count the number of reconstructible tracks
	125	TIter next(trackRefStore->CreateIterator());
	126	AliMUONTrack* trackRef;
	127	while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
f202486b	128	if (trackRef->IsValid(requestedStationMask, request2ChInSameSt45)) nReconstructibleTracks++;
23035434	129	}
23035434	130
66085093	131	// loop over ESD tracks
66085093	132	Int_t nTrackerTracks = 0;
23035434	133	trackReconstructedYet = kFALSE;
66085093	134	AliMUONVTrackStore *fakeTrackStore = AliMUONESDInterface::NewTrackStore();
a99c3449	135	const AliESDEvent* esd = rc.GetESDEvent();
66085093	136	Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
	137	for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
	138
a99c3449	139	AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack);
66085093	140
66085093	141	// skip ghosts
a99c3449	142	if (!esdTrack->ContainTrackerData()) continue;
66085093	143	nTrackerTracks++;
66085093	144
a99c3449	145	// find the corresponding MUON track
	146	AliMUONTrack* muonTrack = (AliMUONTrack*) muonTrackStore->FindObject(esdTrack->GetUniqueID());
	147
66085093	148	// get track info
a99c3449	149	Int_t nClusters = esdTrack->GetNClusters();
	150	Double_t normalizedChi2 = esdTrack->GetChi2() / (2. * esdTrack->GetNHit() - 5);
	151	Double_t p = esdTrack->P();
	152	Double_t pT = esdTrack->Pt();
	153	Double_t eta = esdTrack->Eta();
	154	Double_t phi = esdTrack->Phi();
66085093	155
	156	// fill global histograms
	157	hNumberOfClusters->Fill(nClusters);
	158	hChi2PerDof->Fill(normalizedChi2);
	159	hP->Fill(p);
	160	hPt->Fill(pT);
	161	hEta->Fill(eta);
	162	hPhi->Fill(phi);
	163
	164	// try to match the reconstructed track with a simulated one
f202486b	165	Int_t nMatchClusters = 0;
f202486b	166	AliMUONTrack* matchedTrackRef = rc.FindCompatibleTrack(muonTrack, trackRefStore, nMatchClusters, useLabel, sigmaCut);
66085093	167
	168	// take actions according to matching result
	169	if (matchedTrackRef) {
	170
	171	// global counter
	172	nTotMatchedTracks++;
f202486b	173	if (!matchedTrackRef->IsValid(requestedStationMask, request2ChInSameSt45)) {
23035434	174	trackReconstructedYet = kTRUE;
	175	nTotTracksReconstructedYet++;
	176	}
66085093	177
66085093	178	// fill histograms
f202486b	179	hFractionOfMatchedClusters->Fill(((Float_t) nMatchClusters) / ((Float_t) nClusters));
66085093	180	hNumberOfClustersMC->Fill(matchedTrackRef->GetNClusters());
	181	hNumberOfClustersM->Fill(nClusters);
	182	hChi2PerDofM->Fill(normalizedChi2);
	183	hPM->Fill(p);
	184	hPtM->Fill(pT);
	185	hEtaM->Fill(eta);
	186	hPhiM->Fill(phi);
	187
	188	// remove already matched trackRefs
	189	trackRefStore->Remove(*matchedTrackRef);
	190
	191	} else {
	192
	193	// global counter
	194	nTotFakeTracks++;
	195
	196	// fill histograms
	197	hNumberOfClustersF->Fill(nClusters);
	198	hChi2PerDofF->Fill(normalizedChi2);
	199	hPF->Fill(p);
	200	hPtF->Fill(pT);
	201	hEtaF->Fill(eta);
	202	hPhiF->Fill(phi);
	203
	204	// store fake tracks
a99c3449	205	fakeTrackStore->Add(*muonTrack);
66085093	206
	207	}
	208
	209	} // end of loop over ESD tracks
	210
	211	// fill histograms
	212	hNumberOfTracks->Fill(nTrackerTracks);
23035434	213	nReconstructedTracks += nTrackerTracks;
23035434	214	if (trackReconstructedYet) eventsWithTrackReconstructedYet[nEventsWithTrackReconstructedYet++] = iEvent;
66085093	215
	216	// count the number the additional fake tracks
	217	if (fakeTrackStore->GetSize() > 0) {
	218
	219	// remove the most connected fake tracks
f202486b	220	Int_t nFreeMissingTracks = RemoveConnectedFakes(fakeTrackStore, trackRefStore, useLabel, *hFractionOfConnectedClusters);
66085093	221
	222	// remove the remaining free reconstructible tracks
	223	Int_t nAdditionalTracks = fakeTrackStore->GetSize() - nFreeMissingTracks;
	224
	225	// fill histograms
23035434	226	eventsWithFake[nEventsWithFake] = iEvent;
66085093	227	nEventsWithFake++;
66085093	228	if (nAdditionalTracks > 0) {
23035434	229	eventsWithAdditionalFake[nEventsWithAdditionalFake] = iEvent;
66085093	230	nEventsWithAdditionalFake++;
	231	nTotAdditionalTracks += nAdditionalTracks;
	232	hNumberOfAdditionalTracks->Fill(nAdditionalTracks);
	233	}
	234
	235	}
	236
	237	delete fakeTrackStore;
	238
	239	} // end of loop over events
23035434	240
	241	// total number of connected tracks
	242	nTotConnectedTracks = hFractionOfConnectedClusters->GetEntries();
66085093	243
	244	// plot results
	245	TCanvas cFakesSummary("cFakesSummary","cFakesSummary",900,600);
	246	cFakesSummary.Divide(3,2);
	247	cFakesSummary.cd(1);
	248	cFakesSummary.GetPad(1)->SetLogy();
	249	hNumberOfClusters->Draw();
	250	hNumberOfClusters->SetMinimum(0.5);
	251	hNumberOfClustersM->Draw("same");
	252	hNumberOfClustersM->SetLineColor(4);
	253	hNumberOfClustersF->Draw("same");
	254	hNumberOfClustersF->SetLineColor(2);
	255	hNumberOfClustersF->SetFillColor(2);
	256	hNumberOfClustersF->SetFillStyle(3017);
	257	cFakesSummary.cd(2);
	258	cFakesSummary.GetPad(2)->SetLogy();
	259	hChi2PerDof->Draw();
	260	hChi2PerDof->SetMinimum(0.5);
	261	hChi2PerDofM->Draw("same");
	262	hChi2PerDofM->SetLineColor(4);
	263	hChi2PerDofF->Draw("same");
	264	hChi2PerDofF->SetLineColor(2);
	265	hChi2PerDofF->SetFillColor(2);
	266	hChi2PerDofF->SetFillStyle(3017);
	267	cFakesSummary.cd(3);
	268	cFakesSummary.GetPad(3)->SetLogy();
	269	hP->Draw();
	270	hP->SetMinimum(0.5);
	271	hPM->Draw("same");
	272	hPM->SetLineColor(4);
	273	hPF->Draw("same");
	274	hPF->SetLineColor(2);
	275	hPF->SetFillColor(2);
	276	hPF->SetFillStyle(3017);
	277	cFakesSummary.cd(4);
	278	cFakesSummary.GetPad(4)->SetLogy();
	279	hPt->Draw();
	280	hPt->SetMinimum(0.5);
	281	hPtM->Draw("same");
	282	hPtM->SetLineColor(4);
	283	hPtF->Draw("same");
	284	hPtF->SetLineColor(2);
	285	hPtF->SetFillColor(2);
	286	hPtF->SetFillStyle(3017);
	287	cFakesSummary.cd(5);
	288	cFakesSummary.GetPad(5)->SetLogy();
	289	hEta->Draw();
	290	hEta->SetMinimum(0.5);
	291	hEtaM->Draw("same");
	292	hEtaM->SetLineColor(4);
	293	hEtaF->Draw("same");
	294	hEtaF->SetLineColor(2);
	295	hEtaF->SetFillColor(2);
	296	hEtaF->SetFillStyle(3017);
	297	cFakesSummary.cd(6);
	298	cFakesSummary.GetPad(6)->SetLogy();
	299	hPhi->Draw();
	300	hPhi->SetMinimum(0.5);
	301	hPhiM->Draw("same");
	302	hPhiM->SetLineColor(4);
	303	hPhiF->Draw("same");
	304	hPhiF->SetLineColor(2);
	305	hPhiF->SetFillColor(2);
	306	hPhiF->SetFillStyle(3017);
307
308	// save results
309	histoFile->cd();
310	histoFile->Write();
311	cFakesSummary.Write();
312	histoFile->Close();
313
314	// print results
315	cout << endl;
23035434	316	cout << "- Number of reconstructible tracks: " << nReconstructibleTracks << endl;
23035434	317	cout << "- Number of reconstructed tracks: " << nReconstructedTracks << endl;
66085093	318	cout << "- Number of matched tracks: " << nTotMatchedTracks << endl;
23035434	319	cout << " (including " << nTotTracksReconstructedYet << " track(s) matched with a TrackRef that is not reconstructible";
	320	if (nTotTracksReconstructedYet > 0) {
	321	for(Int_t i=0; i<nEventsWithTrackReconstructedYet; i++){
	322	if (i==0) cout << " (eventID = " << eventsWithTrackReconstructedYet[i];
	323	else cout << ", " << eventsWithTrackReconstructedYet[i];
	324	}
	325	cout << "))" << endl;
	326	} else cout << ")" << endl;
66085093	327	cout << "- Number of fake tracks: " << nTotFakeTracks << endl;
23035434	328	cout << " (including " << nTotConnectedTracks << " track(s) still connected to a reconstructible one)" << endl;
	329	cout << " (including " << nTotAdditionalTracks << " additional track(s) (compared to the number of expected ones))" << endl;
	330	cout << "- Number of events with fake track(s): " << nEventsWithFake;
	331	if (nEventsWithFake > 0) {
	332	for(Int_t i=0; i<nEventsWithFake; i++){
	333	if (i==0) cout << " (eventID = " << eventsWithFake[i];
	334	else cout << ", " << eventsWithFake[i];
	335	}
	336	cout << ")" << endl;
	337	} else cout << endl;
	338	cout << " (including " << nEventsWithAdditionalFake << " events with additional track(s)";
	339	if (nEventsWithAdditionalFake > 0) {
	340	for(Int_t i=0; i<nEventsWithAdditionalFake; i++){
	341	if (i==0) cout << " (eventID = " << eventsWithAdditionalFake[i];
	342	else cout << ", " << eventsWithAdditionalFake[i];
	343	}
	344	cout << "))" << endl;
	345	} else cout << ")" << endl;
66085093	346	cout << endl;
	347	cout << "REMINDER: results are relevent provided that you use the same recoParams as for the reconstruction" << endl;
	348	cout << endl;
	349
66085093	350	}
	351
	352	//-----------------------------------------------------------------------
f202486b	353	Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore,
f202486b	354	Bool_t useLabel, TH1F &hFractionOfConnectedClusters)
66085093	355	{
	356	/// loop over reconstructible TrackRef not associated with reconstructed track:
	357	/// for each of them, find and remove the most connected the fake track, if any,
	358	/// and fill the histograms with the fraction of connected clusters.
	359	/// Return the number of reconstructible track not connected to any fake
	360
	361	Int_t nFreeMissingTracks = 0;
	362
	363	// loop over trackRefs
	364	TIter next(trackRefStore.CreateIterator());
	365	AliMUONTrack* trackRef;
	366	while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
	367
	368	// skip not reconstructible trackRefs
f202486b	369	if (!trackRef->IsValid(requestedStationMask, request2ChInSameSt45)) continue;
66085093	370
	371	Int_t label = trackRef->GetUniqueID();
	372
	373	// look for the most connected fake track
	374	AliMUONTrack *connectedFake = 0x0;
a99c3449	375	Double_t fractionOfConnectedClusters = 0.;
66085093	376	TIter next2(fakeTrackStore.CreateIterator());
	377	AliMUONTrack* fakeTrack;
	378	while ( ( fakeTrack = static_cast<AliMUONTrack*>(next2()) ) ) {
	379
	380	// get the number of connected clusters
	381	Int_t nConnectedClusters = 0;
	382	if (useLabel) { // by using the MC label
	383	for (Int_t iCl = 0; iCl < fakeTrack->GetNClusters(); iCl++)
	384	if (((AliMUONTrackParam*) fakeTrack->GetTrackParamAtCluster()->UncheckedAt(iCl))->GetClusterPtr()->GetMCLabel() == label)
	385	nConnectedClusters++;
	386	} else { // by comparing cluster/TrackRef positions
	387	Bool_t compTrack[10];
f202486b	388	nConnectedClusters = fakeTrack->FindCompatibleClusters(*trackRef, sigmaCut, compTrack);
66085093	389	}
	390
	391	// skip non-connected fake tracks
	392	if (nConnectedClusters == 0) continue;
	393
	394	// check if it is the most connected fake track
a99c3449	395	Double_t f = ((Double_t)nConnectedClusters) / ((Double_t)fakeTrack->GetNClusters());
66085093	396	if (f > fractionOfConnectedClusters) {
	397	connectedFake = fakeTrack;
	398	fractionOfConnectedClusters = f;
	399	}
	400
	401	}
	402
	403	// remove the most connected fake track
	404	if (connectedFake) {
	405	hFractionOfConnectedClusters.Fill(fractionOfConnectedClusters);
	406	fakeTrackStore.Remove(*connectedFake);
	407	} else nFreeMissingTracks++;
	408
	409	}
	410
	411	return nFreeMissingTracks;
	412
	413	}
	414