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