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

#if !defined(__CINT__) || defined(__MAKECINT__)
// ROOT includes
#include <TFile.h>
#include <TH1.h>
#include <TCanvas.h>
#include <Riostream.h>
#include <TROOT.h>
#include <TClonesArray.h>
#include <TLorentzVector.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 "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

Double_t sigmaCut = -1.;

//-----------------------------------------------------------------------
void DIMUONFakes(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();
  
  AliLog::SetClassDebugLevel("AliMCEvent",-1);
  
  // File for histograms and histogram booking
  TFile *histoFile = new TFile("DiFakes.root", "RECREATE");
  
  TH1F *hMass = new TH1F("hMass", "Dimuon 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", "Dimuon 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", "Dimuon 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"," Dimuon 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"," Dimuon 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"," Dimuon 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()->SetSpecificStorage("GRP/GRP/Data","local://.");
  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
  UInt_t requestedStationMask = 0;
  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
  Bool_t request2ChInSameSt45 = !recoParam->MakeMoreTrackCandidates();
  
  TLorentzVector vMu1, vMu2, vDiMu;
  Int_t nReconstructiblePairs = 0;
  Int_t nReconstructedPairs = 0;
  Int_t nMatchedPairs = 0;
  
  // 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 pairs
    Int_t nMuPlus = 0, nMuMinus = 0;
    TIter next(trackRefStore->CreateIterator());
    AliMUONTrack* trackRef;
    while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
      if (trackRef->IsValid(requestedStationMask, request2ChInSameSt45)) {
	if (trackRef->GetTrackParamAtVertex()->GetCharge() > 0) nMuPlus++;
	else nMuMinus++;
      }
    }
    nReconstructiblePairs += nMuPlus*nMuMinus;
    
    // loop over ESD tracks and flag them
    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;
      
      // find the corresponding MUON track
      AliMUONTrack* muonTrack = (AliMUONTrack*) muonTrackStore->FindObject(esdTrack->GetUniqueID());
      
      // 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) {
	
	// flag matched tracks
	esdTrack->SetLabel(matchedTrackRef->GetUniqueID());
	
	// remove already matched trackRefs
	trackRefStore->Remove(*matchedTrackRef);
	
      } else {
	
	// flag fake tracks
	esdTrack->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;
	
	nReconstructedPairs++;
	
	// 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 = vDiMu.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) {
	  
	  nMatchedPairs++;
	  
	  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();
  
  // print results
  printf("\n");
  printf("nb of reconstructible OS pairs: %d \n", nReconstructiblePairs);
  printf("nb of reconstructed OS pairs: %d \n", nReconstructedPairs);
  printf("nb of reconstructed OS pairs matched with trackRefs: %d \n", nMatchedPairs);
  printf("\n");
  
}
Commit	Line	Data
30b3fe0f	1	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
30b3fe0f	2	// ROOT includes
30b3fe0f	3	#include <TFile.h>
	4	#include <TH1.h>
	5	#include <TCanvas.h>
	6	#include <Riostream.h>
	7	#include <TROOT.h>
	8	#include <TClonesArray.h>
30b3fe0f	9	#include <TLorentzVector.h>
	10
	11	// STEER includes
	12	#include "AliLog.h"
30b3fe0f	13	#include "AliESDEvent.h"
30b3fe0f	14	#include "AliESDMuonTrack.h"
30b3fe0f	15	#include "AliCDBManager.h"
	16
	17	// MUON includes
a99c3449	18	#include "AliMUONCDB.h"
30b3fe0f	19	#include "AliMUONTrack.h"
	20	#include "AliMUONVTrackStore.h"
	21	#include "AliMUONTrackParam.h"
30b3fe0f	22	#include "AliMUONRecoCheck.h"
	23	#include "AliMUONVCluster.h"
	24	#include "AliMUONRecoParam.h"
	25	#endif
	26
	27	/// \ingroup macros
	28	/// \file DIMUONFakes.C
	29	///
	30	/// \author Ph. Pillot, Subatech, March. 2009
	31	///
	32	/// Macro to study the effects of fake tracks on the dimuon spectra
	33	/// Results are saved in the root file DiFakes.root
	34	/// Results are relevent provided that you use the same recoParams as for the reconstruction
	35
f202486b	36	Double_t sigmaCut = -1.;
30b3fe0f	37
	38	//-----------------------------------------------------------------------
	39	void DIMUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent = -1,
a99c3449	40	const TString esdFileName = "AliESDs.root", const TString SimDir = "./generated/",
a99c3449	41	const TString ocdbPath = "local://$ALICE_ROOT/OCDB")
30b3fe0f	42	{
	43
	44	//Reset ROOT and connect tree file
	45	gROOT->Reset();
	46
590a335d	47	AliLog::SetClassDebugLevel("AliMCEvent",-1);
590a335d	48
30b3fe0f	49	// File for histograms and histogram booking
	50	TFile *histoFile = new TFile("DiFakes.root", "RECREATE");
	51
23035434	52	TH1F *hMass = new TH1F("hMass", "Dimuon mass distribution (GeV/c^{2})", 100, 0., 12.);
30b3fe0f	53	TH1F *hMassM = new TH1F("hMassM", "matched track mass distribution (GeV/c^{2})", 100, 0., 12.);
30b3fe0f	54	TH1F *hMassF = new TH1F("hMassF", "fake track mass distribution (GeV/c^{2})", 100, 0., 12.);
23035434	55	TH1F *hP = new TH1F("hP", "Dimuon P distribution (GeV/c)", 100, 0., 200.);
30b3fe0f	56	TH1F *hPM = new TH1F("hPM", "matched track P distribution (GeV/c)", 100, 0., 200.);
30b3fe0f	57	TH1F *hPF = new TH1F("hPF", "fake track P distribution (GeV/c)", 100, 0., 200.);
23035434	58	TH1F *hPt = new TH1F("hPt", "Dimuon Pt distribution (GeV/c)", 100, 0., 20.);
30b3fe0f	59	TH1F *hPtM = new TH1F("hPtM", "matched track Pt distribution (GeV/c)", 100, 0., 20.);
30b3fe0f	60	TH1F *hPtF = new TH1F("hPtF", "fake track Pt distribution (GeV/c)", 100, 0., 20.);
23035434	61	TH1F *hY = new TH1F("hY"," Dimuon rapidity distribution",100,-10,0);
30b3fe0f	62	TH1F *hYM = new TH1F("hYM"," matched track rapidity distribution",100,-10,0);
30b3fe0f	63	TH1F *hYF = new TH1F("hYF"," fake track rapidity distribution",100,-10,0);
23035434	64	TH1F *hEta = new TH1F("hEta"," Dimuon pseudo-rapidity distribution",100,-10,0);
30b3fe0f	65	TH1F *hEtaM = new TH1F("hEtaM"," matched track pseudo-rapidity distribution",100,-10,0);
30b3fe0f	66	TH1F *hEtaF = new TH1F("hEtaF"," fake track pseudo-rapidity distribution",100,-10,0);
23035434	67	TH1F *hPhi = new TH1F("hPhi"," Dimuon phi distribution",100,-1.,9.);
30b3fe0f	68	TH1F *hPhiM = new TH1F("hPhiM"," matched track phi distribution",100,-1.,9.);
	69	TH1F *hPhiF = new TH1F("hPhiF"," fake track phi distribution",100,-1.,9.);
	70
a99c3449	71	// link to reconstructed and simulated tracks
a99c3449	72	AliMUONRecoCheck rc(esdFileName, SimDir);
30b3fe0f	73
a99c3449	74	// load necessary data from OCDB
a99c3449	75	AliCDBManager::Instance()->SetDefaultStorage(ocdbPath);
590a335d	76	AliCDBManager::Instance()->SetSpecificStorage("GRP/GRP/Data","local://.");
a99c3449	77	AliCDBManager::Instance()->SetRun(rc.GetRunNumber());
	78	AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
	79	if (!recoParam) return;
30b3fe0f	80
a99c3449	81	// get sigma cut from recoParam to associate clusters with TrackRefs in case the label are not used
f202486b	82	sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking();
590a335d	83	// compute the mask of requested stations from recoParam
	84	UInt_t requestedStationMask = 0;
	85	for (Int_t i = 0; i < 5; i++) if (recoParam->RequestStation(i)) requestedStationMask \|= ( 1 << i );
	86	// get from recoParam whether a track need 2 chambers hit in the same station (4 or 5) or not to be reconstructible
	87	Bool_t request2ChInSameSt45 = !recoParam->MakeMoreTrackCandidates();
30b3fe0f	88
30b3fe0f	89	TLorentzVector vMu1, vMu2, vDiMu;
590a335d	90	Int_t nReconstructiblePairs = 0;
	91	Int_t nReconstructedPairs = 0;
	92	Int_t nMatchedPairs = 0;
30b3fe0f	93
	94	// Loop over ESD events
	95	FirstEvent = TMath::Max(0, FirstEvent);
a99c3449	96	LastEvent = (LastEvent>=0) ? TMath::Min(rc.NumberOfEvents() - 1, LastEvent) : rc.NumberOfEvents() - 1;
30b3fe0f	97	for (Int_t iEvent = FirstEvent; iEvent <= LastEvent; iEvent++) {
30b3fe0f	98
a99c3449	99	// get reconstructed and simulated tracks
a99c3449	100	AliMUONVTrackStore* muonTrackStore = rc.ReconstructedTracks(iEvent, kFALSE);
30b3fe0f	101	AliMUONVTrackStore* trackRefStore = rc.TrackRefs(iEvent);
a99c3449	102	if (!muonTrackStore \|\| !trackRefStore) continue;
30b3fe0f	103
590a335d	104	// count the number of reconstructible pairs
	105	Int_t nMuPlus = 0, nMuMinus = 0;
	106	TIter next(trackRefStore->CreateIterator());
	107	AliMUONTrack* trackRef;
	108	while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
	109	if (trackRef->IsValid(requestedStationMask, request2ChInSameSt45)) {
	110	if (trackRef->GetTrackParamAtVertex()->GetCharge() > 0) nMuPlus++;
	111	else nMuMinus++;
	112	}
	113	}
	114	nReconstructiblePairs += nMuPlus*nMuMinus;
	115
30b3fe0f	116	// loop over ESD tracks and flag them
a99c3449	117	const AliESDEvent* esd = rc.GetESDEvent();
30b3fe0f	118	Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
	119	for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
	120
a99c3449	121	AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack);
30b3fe0f	122
30b3fe0f	123	// skip ghosts
a99c3449	124	if (!esdTrack->ContainTrackerData()) continue;
	125
	126	// find the corresponding MUON track
	127	AliMUONTrack* muonTrack = (AliMUONTrack*) muonTrackStore->FindObject(esdTrack->GetUniqueID());
30b3fe0f	128
30b3fe0f	129	// try to match the reconstructed track with a simulated one
f202486b	130	Int_t nMatchClusters = 0;
f202486b	131	AliMUONTrack* matchedTrackRef = rc.FindCompatibleTrack(muonTrack, trackRefStore, nMatchClusters, useLabel, sigmaCut);
30b3fe0f	132
	133	// take actions according to matching result
	134	if (matchedTrackRef) {
	135
	136	// flag matched tracks
a99c3449	137	esdTrack->SetLabel(matchedTrackRef->GetUniqueID());
30b3fe0f	138
	139	// remove already matched trackRefs
	140	trackRefStore->Remove(*matchedTrackRef);
	141
	142	} else {
	143
	144	// flag fake tracks
a99c3449	145	esdTrack->SetLabel(-1);
30b3fe0f	146
	147	}
	148
	149	}
	150
	151	// double loop over ESD tracks, build pairs and fill histograms according to their label
	152	for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
	153	AliESDMuonTrack* muonTrack1 = esd->GetMuonTrack(iTrack1);
	154
	155	// skip ghosts
	156	if (!muonTrack1->ContainTrackerData()) continue;
	157
	158	// get track info
	159	Short_t charge1 = muonTrack1->Charge();
	160	Int_t label1 = muonTrack1->GetLabel();
	161	muonTrack1->LorentzP(vMu1);
	162
	163	for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
	164	AliESDMuonTrack* muonTrack2 = esd->GetMuonTrack(iTrack2);
	165
	166	// skip ghosts
	167	if (!muonTrack2->ContainTrackerData()) continue;
	168
	169	// keep only opposite sign pairs
	170	Short_t charge2 = muonTrack2->Charge();
	171	if (charge1*charge2 > 0) continue;
	172
590a335d	173	nReconstructedPairs++;
590a335d	174
30b3fe0f	175	// get track info
	176	Int_t label2 = muonTrack2->GetLabel();
	177	muonTrack2->LorentzP(vMu2);
	178
	179	// compute kinematics of the pair
	180	vDiMu = vMu1 + vMu2;
	181	Float_t mass = vDiMu.M();
	182	Float_t p = vDiMu.P();
	183	Float_t pt = vDiMu.Pt();
	184	Float_t y = vDiMu.Rapidity();
	185	Float_t eta = vDiMu.Eta();
23035434	186	Float_t phi = vDiMu.Phi();
30b3fe0f	187	if (phi < 0) phi += 2.*TMath::Pi();
	188
	189	// fill global histograms
	190	hMass->Fill(mass);
	191	hP->Fill(p);
	192	hPt->Fill(pt);
	193	hY->Fill(y);
	194	hEta->Fill(eta);
	195	hPhi->Fill(phi);
	196
	197	// fill histograms according to labels
	198	if (label1 >= 0 && label2 >= 0) {
	199
590a335d	200	nMatchedPairs++;
590a335d	201
30b3fe0f	202	hMassM->Fill(mass);
	203	hPM->Fill(p);
	204	hPtM->Fill(pt);
	205	hYM->Fill(y);
	206	hEtaM->Fill(eta);
	207	hPhiM->Fill(phi);
	208
	209	} else {
	210
	211	hMassF->Fill(mass);
	212	hPF->Fill(p);
	213	hPtF->Fill(pt);
	214	hYF->Fill(y);
	215	hEtaF->Fill(eta);
	216	hPhiF->Fill(phi);
	217
	218	}
	219
	220	}
	221
	222	}
	223
	224	} // end of loop over events
	225
	226	// plot results
	227	TCanvas cDiFakesSummary("cDiFakesSummary","cDiFakesSummary",900,600);
	228	cDiFakesSummary.Divide(3,2);
	229	cDiFakesSummary.cd(1);
	230	cDiFakesSummary.GetPad(1)->SetLogy();
	231	hMass->Draw();
	232	hMass->SetMinimum(0.5);
	233	hMassM->Draw("same");
	234	hMassM->SetLineColor(4);
	235	hMassF->Draw("same");
	236	hMassF->SetLineColor(2);
	237	hMassF->SetFillColor(2);
	238	hMassF->SetFillStyle(3017);
	239	cDiFakesSummary.cd(2);
	240	cDiFakesSummary.GetPad(3)->SetLogy();
	241	hP->Draw();
	242	hP->SetMinimum(0.5);
	243	hPM->Draw("same");
	244	hPM->SetLineColor(4);
	245	hPF->Draw("same");
	246	hPF->SetLineColor(2);
	247	hPF->SetFillColor(2);
	248	hPF->SetFillStyle(3017);
	249	cDiFakesSummary.cd(3);
	250	cDiFakesSummary.GetPad(4)->SetLogy();
	251	hPt->Draw();
	252	hPt->SetMinimum(0.5);
	253	hPtM->Draw("same");
	254	hPtM->SetLineColor(4);
	255	hPtF->Draw("same");
	256	hPtF->SetLineColor(2);
	257	hPtF->SetFillColor(2);
	258	hPtF->SetFillStyle(3017);
	259	cDiFakesSummary.cd(4);
	260	cDiFakesSummary.GetPad(2)->SetLogy();
	261	hY->Draw();
	262	hY->SetMinimum(0.5);
	263	hYM->Draw("same");
	264	hYM->SetLineColor(4);
	265	hYF->Draw("same");
266	hYF->SetLineColor(2);
267	hYF->SetFillColor(2);
268	hYF->SetFillStyle(3017);
269	cDiFakesSummary.cd(5);
270	cDiFakesSummary.GetPad(5)->SetLogy();
271	hEta->Draw();
272	hEta->SetMinimum(0.5);
273	hEtaM->Draw("same");
274	hEtaM->SetLineColor(4);
275	hEtaF->Draw("same");
276	hEtaF->SetLineColor(2);
277	hEtaF->SetFillColor(2);
278	hEtaF->SetFillStyle(3017);
279	cDiFakesSummary.cd(6);
280	cDiFakesSummary.GetPad(6)->SetLogy();
281	hPhi->Draw();
282	hPhi->SetMinimum(0.5);
283	hPhiM->Draw("same");
284	hPhiM->SetLineColor(4);
285	hPhiF->Draw("same");
286	hPhiF->SetLineColor(2);
287	hPhiF->SetFillColor(2);
288	hPhiF->SetFillStyle(3017);
289
290	// save results
291	histoFile->cd();
292	histoFile->Write();
293	cDiFakesSummary.Write();
294	histoFile->Close();
295
590a335d	296	// print results
	297	printf("\n");
	298	printf("nb of reconstructible OS pairs: %d \n", nReconstructiblePairs);
	299	printf("nb of reconstructed OS pairs: %d \n", nReconstructedPairs);
	300	printf("nb of reconstructed OS pairs matched with trackRefs: %d \n", nMatchedPairs);
	301	printf("\n");
	302
30b3fe0f	303	}
30b3fe0f	304