[u/mrichter/AliRoot.git] / PWGCF / Correlations / DPhi / DiHadronPID / AliAODTrackCutsDiHadronPID.h

#ifndef ALIAODTRACKCUTSDIHADRONPID_H
#define ALIAODTRACKCUTSDIHADRONPID_H

#include "TFormula.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TH3F.h"
#include "TList.h"

class AliAODTrackCutsDiHadronPID : public TNamed 

{

public:
	enum HistoClass {kAllCharged = 0, kPositive = 1, kNegative = 2,
		kAllPion = 3, kPosPion = 4, kNegPion = 5, 
		kAllKaon = 6, kPosKaon = 7, kNegKaon = 8, 
		kAllProton = 9, kPosProton = 10, kNegProton = 11};

public:
	AliAODTrackCutsDiHadronPID();					// Default Constructor
	AliAODTrackCutsDiHadronPID(const char* name);	// Named Constructor
	virtual ~AliAODTrackCutsDiHadronPID();			// Destructor
	virtual Long64_t Merge(TCollection* list);		// Merger

private:
	AliAODTrackCutsDiHadronPID(const AliAODTrackCutsDiHadronPID&);
	AliAODTrackCutsDiHadronPID& operator=(const AliAODTrackCutsDiHadronPID&);

// -------------------------------------------------------------------------
//  Interface, methods used to get information about the track cuts, and to
//  retrieve filled histograms. 
// -------------------------------------------------------------------------

public:
	void PrintCuts();								// Gives an overview of the cuts.

// Return the list of QA histos
	TList* GetListOfDataQAHistos() const{
		if (fDataTrackQAHistos) {return fDataTrackQAHistos;}
		else return 0x0;
	}
	TList* GetListOfPrimRecMCTrackQAHistos() const{
		if (fPrimRecMCTrackQAHistos) {return fPrimRecMCTrackQAHistos;}
		else return 0x0;
	}
	TList* GetListOfPrimGenMCTrackQAHistos() const {
		if (fPrimGenMCTrackQAHistos) {return fPrimGenMCTrackQAHistos;}
		else return 0x0;
	}
	TList* GetListOfSecRecMCTrackQAHistos() const {
		if (fSecRecMCTrackQAHistos) return fSecRecMCTrackQAHistos;
		else return 0x0;
	}
	TList* GetListOfSecGenMCTrackQAHistos() const {
		if (fSecGenMCTrackQAHistos) return fSecGenMCTrackQAHistos;
		else return 0x0;
	}

// Note that the PID histograms have in principle a different number of pT bins.
// FIXME: This is not very nice...
	Int_t GetNPtBins() const {return fNPtBins;}
	Int_t GetNPtBinsPID(Int_t ptclass = -1) const {
		// if class = -1, then return the sum.
		if (ptclass == -1) {
			Int_t nptbinspid = 0;
			for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
				nptbinspid += fNPtBinsPID[iPtClass];
			}
			return nptbinspid;
		} else if (ptclass >= 0 && ptclass < 5) {
			return fNPtBinsPID[ptclass];
		} else {return -999;}
	}	

// Returns the Pt axis for PID and other histograms.
	Double_t* GetPtAxis() {return fPtAxis;}
	Double_t* GetPtAxisPID() const {
		const Int_t nptbinspid = GetNPtBinsPID();
		Double_t* ptaxis = new Double_t[nptbinspid + 1];
		for (Int_t iPtBin = 0; iPtBin < nptbinspid; iPtBin++) {
			ptaxis[iPtBin] = GetPtMinPID(iPtBin + 1); 
		}
		ptaxis[nptbinspid] = GetPtMaxPID(nptbinspid);
		return ptaxis;
	}

// Return data histogram with a specific name. Since the histograms are not streamed, and only the
// TList containing them is, we have to retrieve them from the list.
	TObject* GetHistData(const char* name) const {return fDataTrackQAHistos->FindObject(name);} 
	TObject* GetHistPrimRecMC(const char* name) const {return fPrimRecMCTrackQAHistos->FindObject(name);}
	TObject* GetHistPrimGenMC(const char* name) const {return fPrimGenMCTrackQAHistos->FindObject(name);}
	TObject* GetHistSecRecMC(const char* name) const {return fSecRecMCTrackQAHistos->FindObject(name);}
	TObject* GetHistSecGenMC(const char* name) const {return fSecGenMCTrackQAHistos->FindObject(name);}

// Since we will often want to have TOF histograms, here are a few methods which return the 
// appropriate projections. The class does not own these projections, and the user must take care of them.
	TH1F* GetHistDataTOFProjection(Int_t charge, Int_t species, Int_t ptbin);
 	TH1F* GetHistDataTOFMismatch(Int_t charge, Int_t species, Int_t ptbin);
	Double_t GetPtMinPID(Int_t bin) const {
		Int_t ptclass = GetPtClass(bin);
		Int_t bininptclass = GetBinInPtClass(bin);
		Double_t minpt = fPtBoundaryPID[ptclass];
		Double_t maxpt = fPtBoundaryPID[ptclass+1];
		Double_t ptres = (maxpt - minpt)/((Double_t)fNPtBinsPID[ptclass]);
		return (minpt + ptres * ((Double_t)(bininptclass - 1)) );
	}
	Double_t GetPtMaxPID(Int_t bin) const {
		Int_t ptclass = GetPtClass(bin);
		Int_t bininptclass = GetBinInPtClass(bin);
		Double_t minpt = fPtBoundaryPID[ptclass];
		Double_t maxpt = fPtBoundaryPID[ptclass+1];
		Double_t ptres = (maxpt - minpt)/((Double_t)fNPtBinsPID[ptclass]);
		return (minpt + ptres * ((Double_t)(bininptclass)) );
	}
	Double_t GetPtClassMin(Int_t ptclass) const {
		if (ptclass >= 0 && ptclass < 5) {
			return fPtBoundaryPID[ptclass];
		} else {return -999;}
	}
	Double_t GetPtClassMax(Int_t ptclass) const {
		if (ptclass >= 0 && ptclass < 5) {
			return fPtBoundaryPID[ptclass+1];
		} else {return -999;}
	}

	Double_t GetPtBinWidthPID(Int_t bin) const {return (GetPtMaxPID(bin) - GetPtMinPID(bin)); }

// BE CAREFUL! Following methods do not apply to the ptaxis of the PID histograms! For that, call
// the GetPtMinPID and GetPtMaxPID methods. 
	Double_t GetPtMin(Int_t bin) const {
		if ((bin < 1) || (bin > fNPtBins + 1)) {cout<<"Bin is out of range..."<<endl; return -999.;}
		else {return fPtAxis[bin - 1];}
	}
	Double_t GetPtMax(Int_t bin) const {
		if ((bin < 1) || (bin > fNPtBins + 1)) {cout<<"Bin is out of range..."<<endl; return -999.;}
		else {return fPtAxis[bin];}
	}
	Double_t GetPtBinWidth(Int_t bin) const {return (GetPtMax(bin) - GetPtMin(bin)); }

	Int_t GetNTOFbins(Int_t ptclass, Int_t species) const {return fTOFbins[ptclass][species];}
	Double_t GetTOFmin(Int_t ptclass, Int_t species) const {return fTOFLowerBound[ptclass][species];}
	Double_t GetTOFmax(Int_t ptclass, Int_t species) const {return fTOFUpperBound[ptclass][species];}

	Int_t GetNTPCbins(Int_t ptclass, Int_t species) const {return fTPCbins[ptclass][species];}
	Double_t GetTPCmin(Int_t ptclass, Int_t species) const {return fTPCLowerBound[ptclass][species];}
	Double_t GetTPCmax(Int_t ptclass, Int_t species) const {return fTPCUpperBound[ptclass][species];}

// Getters (Cuts)
	UInt_t GetFilterMask() const {return fFilterMask;}
	Double_t GetMaxEta() const {return fMaxEta;}
	ULong_t GetDemandedFlags() const {return fDemandedFlags;}
	TFormula* GetPtDeptDCACutFormula() const {return fPtDeptDCAxyCutFormula;}
	Double_t GetDCAzCut() const {return fDCAzCut;}
	UInt_t GetMinSPDHitsForPtDeptDCACut() const {return fMinSPDHitsForPtDeptDCAcut;}

// Getters (Settings)
	Bool_t GetIsMC() const {return fIsMC;}

	Int_t GetDebugLevel() const {return fDebug;}

// -------------------------------------------------------------------------
//  Methods used to configure the track cuts object, to be called at
//  initialization, i.e., before the object is added to an analysis task.  
// -------------------------------------------------------------------------

public:

// Request Certain QA histograms being filled.
	Bool_t RequestQAHistos(Int_t histoclass, Bool_t Enable3DSpectra = kFALSE, Bool_t EnablePIDHistos = kFALSE) {
		if ((histoclass > -1) && (histoclass < 12)) {
			fHistRequested[histoclass] = kTRUE;
			f3DSpectraEnabeled[histoclass] = Enable3DSpectra;
			fPIDHistosEnabeled[histoclass] = EnablePIDHistos;
			//cout<<"histoclass: "<<histoclass<<" requested: "<<fHistRequested[histoclass]<<endl;
			return kTRUE;
		} else { 
			return kFALSE;
		}
	}
	
// Setters (Cuts)
	void SetPtRange(Double_t minpt, Double_t maxpt) {
		fMinPt = minpt;
		fMaxPt = maxpt;
		fTestPt = kTRUE;
	}
	void SetFilterMask(UInt_t filtermask) {
		fFilterMask = filtermask;
		fTestFilterMask = kTRUE;
	}
	void SetMaxEta(Double_t maxeta) {
		fMaxEta = maxeta;
		fTestMaxEta = kTRUE;
	}
	void SetMaxRapidity(Double_t maxrapidity) {
		fMaxRapidity = maxrapidity;
		fTestMaxRapidity = kTRUE;
	}
	void SetDemandNoMismatch() {
		fTestTOFmismatch = kTRUE;
	}
	void SetDemandFlags(ULong_t demandedflags) {
		fDemandedFlags = demandedflags;
		fTestFlags = kTRUE;
	}
	void SetPtDeptDCACut(TFormula* DCAxyCutFormula, Double_t DCAzCut, UInt_t MinSPDHits = 1) {
		fPtDeptDCAxyCutFormula = DCAxyCutFormula;
		fDCAzCut = DCAzCut;
		fMinSPDHitsForPtDeptDCAcut = MinSPDHits;
		fTestPtDeptDCAcut = kTRUE;
	}

// Setters (Settings)
	void SetIsMC(Bool_t ismc = kTRUE) {fIsMC = ismc;}

	void SetDebugLevel(Int_t debuglevel) {fDebug = debuglevel;}

// -------------------------------------------------------------------------
//  Methods called by the analysis task. 
// -------------------------------------------------------------------------

public:

// These two functions signal the beginning and the end of a new event.
	void StartNewEvent();				// Some things are set to zero.
	void EventIsDone(Bool_t IsMC);		// Some final histograms are filled.
	void CreateHistos(); 				// Should be called by the UserCreateOutput() function of the analysis task.

// Is Selected, for different types of tracks.
	Bool_t IsSelectedData(AliTrackDiHadronPID* track, Double_t randomhittime = -1.e20);
	Bool_t IsSelectedGeneratedMC(AliAODMCParticle* particle);
	Bool_t IsSelectedReconstructedMC(AliTrackDiHadronPID* track);

// -------------------------------------------------------------------------
//  Internal methods.
// -------------------------------------------------------------------------

public:

// For PID histograms we have a certain number of bins in pT, spread out over five
// large histograms, i.e., one for the lowest pT, and the biggest range in TOF/TPC,
// one for the higher pT and smaller range in TOF/TPC, etc. The following methods
// are a mapping between the total pT bin (what the user uses), and the pt bin
// within one of the five histograms (what's used internally)
	Int_t GetPtClass(const Int_t ptbin) const {
		
		// Returns a number [0..4]
		Int_t currentptclass = 0;
		Int_t currentptbin = fNPtBinsPID[0];
		while (currentptbin < ptbin) {
			currentptclass++;
			if (currentptclass == 5) {break;}
			currentptbin += fNPtBinsPID[currentptclass];
		}
		if (currentptclass == 5) {cout<<"GetPtClass -> ptbin out of range!"<<endl; return -1;}
		return currentptclass;
	}
	Int_t GetBinInPtClass(const Int_t ptbin) const {

		// Returns a number [1..Nbins]
		Int_t ptclass = GetPtClass(ptbin);
		if (ptclass == -1) {return -1;}

		Int_t ptbinout = ptbin;
		for (Int_t iPtClass = 0; iPtClass < ptclass; iPtClass++) {ptbinout -= fNPtBinsPID[iPtClass];}

		return ptbinout;

	}

private:

// Checks, return kTRUE if track passes the cut.
	Bool_t CheckPt(Double_t pt) const {
		if (!fTestPt) return kTRUE;
		if ((pt > fMinPt) && (pt < fMaxPt)) return kTRUE;
		return kFALSE;
	}
	Bool_t CheckMaxEta(Double_t eta) const {
		if (!fTestMaxEta) return kTRUE;				// Accepted if there is no check on this parameter.
		if (TMath::Abs(eta) < fMaxEta) return kTRUE;
		return kFALSE;
	}
	Bool_t CheckRapidity(Double_t rap) const {
		if (!fTestMaxRapidity) return kTRUE;
		if (TMath::Abs(rap) < fMaxRapidity) return kTRUE;
		return kFALSE;
	}
	Bool_t CheckFilterMask(UInt_t filtermap) const {
		if (!fTestFilterMask) return kTRUE;
		if ((fFilterMask & filtermap) == fFilterMask) return kTRUE;
		return kFALSE;
	}
	Bool_t CheckFlags(ULong_t flags) const {
		if (!fTestFlags) return kTRUE;
		if ((flags & fDemandedFlags) == fDemandedFlags) return kTRUE;
		return kFALSE;
	}
	Bool_t CheckTOFmismatch(Bool_t ismismatch) const {
		if (!fTestTOFmismatch) return kTRUE; // if we're not cutting on mismatch, then it's accepted.
		if (!ismismatch) return kTRUE; 		// so if the track is not a mismatch, then it is accepted.
		return kFALSE; 						// if it is a mismatch, then it's not accepted.
	}
	Bool_t CheckPtDeptDCACut(Double_t dcaz, Double_t dcaxy, Double_t pt, UInt_t SPDhits) const {
		if (!fTestPtDeptDCAcut) return kTRUE;
		if (SPDhits < fMinSPDHitsForPtDeptDCAcut) return kTRUE; // If there are not enough SPD hits to do the cut.
		if ((dcaz < fDCAzCut) && (dcaxy < fPtDeptDCAxyCutFormula->Eval(pt))) return kTRUE;
		return kFALSE;
	}

// Filling QA histograms.
	Bool_t FillDataHistos(Int_t histoclass, AliTrackDiHadronPID* track);
	Bool_t FillTOFMismatchHistos(Int_t histoclass, AliTrackDiHadronPID* track, Double_t randomhittime);
	Bool_t FillGenMCHistos(Int_t histoclass, AliAODMCParticle* particle);
	Bool_t FillRecMCHistos(Int_t histoclass, AliTrackDiHadronPID* track);

// Initializing QA histograms.
	Bool_t InitializeDataHistos(Int_t histoclass);
	Bool_t InitializeGenMCHistos(Int_t histoclass);
	Bool_t InitializeRecMCHistos(Int_t histoclass);

	void InitializeDefaultHistoNamesAndAxes();

	TH1F* InitializePtSpectrum(const char* name, Int_t histoclass);
	TH3F* InitializePhiEtaPt(const char* name, Int_t histoclass);
	TH1F* InitializeNTracksHisto(const char* name, Int_t histoclass);
	TH1F* InitializeDCAxyHisto(const char* name, Int_t histoclass);
	TH1F* InitializeDCAzHisto(const char* name, Int_t histoclass);
	TH3F* InitializeAcceptanceHisto(const char* /*name*/, Int_t /*histoclass*/); // TO BE IMPLEMENTED.
	TH2F* InitializeDCASpectrum(const char* name, Int_t histoclass);


	TH3F* InitializePIDHisto(const char* name, Int_t histoclass, Int_t expspecies, Int_t ptclass);
	TH2F* InitializeTOFMismatchHisto(const char* name, Int_t histoclass, Int_t expspecies, Int_t ptclass);

// -------------------------------------------------------------------------
//  Data members.
// -------------------------------------------------------------------------

private:
// Track Cuts
	Double_t				fMinPt;
	Double_t				fMaxPt;
	UInt_t 					fFilterMask;					// FilterMask to-be-checked.
	Double_t 				fMaxEta;						// Max Eta of the track.
	Double_t				fMaxRapidity;					// Rapidity cut (only done for PID plots!!)
	ULong_t 				fDemandedFlags;					// Flags demanded on the track.
	UInt_t 					fMinSPDHitsForPtDeptDCAcut;		// Required number of SPD hits for performing Pt-Dept DCA cut.
	TFormula* 				fPtDeptDCAxyCutFormula;			// Formula for the Pt-Dept DCA cut.
	Double_t				fDCAzCut;						// Max z at DCA.

// Settings
	Bool_t					fIsMC;							// Is the current event MC or not.

// Requested Histograms;
	Bool_t					fHistRequested[12];				//
	Bool_t					f3DSpectraEnabeled[12];			//
	Bool_t					fPIDHistosEnabeled[12];			//

// Which Track Cuts will be tested.
	Bool_t					fTestPt;						//
	Bool_t 					fTestFilterMask;				//
	Bool_t 					fTestMaxEta;					//
	Bool_t					fTestMaxRapidity;				//
	Bool_t 					fTestFlags;						//
	Bool_t 					fTestTOFmismatch;				//
	Bool_t 					fTestPtDeptDCAcut;				//

// QA histograms for Data.
	TList*					fDataTrackQAHistos;				// 
	TH1F*					fHistDataPt[3];					//! Pt distribution of tracks passing this cut.
	TH3F*					fHistDataPhiEtaPt[3];			//! Pt, Eta, Phi distribution.
	TH1F*					fHistDataNTracks[3];			//! Number of tracks passing the cut per event (filling by EventIsDone()).
	TH1F*					fHistDataDCAxy[3];				//! DCA_{xy} distribution.
	TH1F*					fHistDataDCAz[3];				//! DCA_{z} distribution
	TH2F* 					fHistDataDCAxyOneSigma[12];		//! DCA_{xy} distribution of particles as identified by 1 sigma method.
	Int_t					fNTracks[3];					//! Number of tracks

	TH3F* 					fHistDataPID[3][3][5];			//! TPC/TOF v.s. pT, [charge][mass assumption][ptclass]
	TH2F*					fHistTOFMismatch[3][3][5];		//! TOF Mismatch histograms, [charge][mass assumption][ptclass]

// QA histograms for Primary Reconstructed MC tracks.
	TList*					fPrimRecMCTrackQAHistos;		//
	TH1F*					fHistPrimRecMCPt[12];			//! Pt distribution of reconstructed MC track passing this cut.
	TH1F*					fHistPrimRecNTracks[12];		//!
	TH2F*					fHistPrimRecMCDCA[12];			//! DCA_xy distribution of reconstructed MC track passing this cut.

// QA histograms for Primary Generated MC particles.
	TList*					fPrimGenMCTrackQAHistos;		//
	TH1F*					fHistPrimGenMCPt[12];			//! Pt distribution of generated MC particles passing this cut.

// QA histograms for Secondary Reconstructed MC tracks.
	TList*					fSecRecMCTrackQAHistos;			//
	TH1F*					fHistSecRecMCPt[12];			//! Pt distribution of reconstructed MC track passing this cut.
	TH2F*					fHistSecRecMCDCAMat[12];		//! DCA_xy distribution of material decay particles.
	TH2F*					fHistSecRecMCDCAWeak[12];		//! DCA_xy distribution of weak decay.

// QA histograms for Secondary Generated MC particles.
	TList*					fSecGenMCTrackQAHistos;			//
	TH1F*					fHistSecGenMCPt[12];			//! Pt distribution of generated MC particles passing this cut.

// Binning of all the histograms.
	Double_t				fPtAxis[57];					// Pt axis used in all histograms, except PID and Mismatch histograms.
	Int_t					fNPtBins;						// Number of bins in the pt-axis.
	Int_t					fNEtaBins;						// 
	Int_t					fNPhiBins;						//

	Double_t				fPtBoundaryPID[6];				// There are five different PID histo's. This array gives the pT range of these histograms.
	Int_t					fNPtBinsPID[5];					// This array gives the number of pT bins for each of these histograms.

	Double_t				fTOFLowerBound[5][3];			// These arrays give the lower and upper bound of the TOF axes,
	Double_t				fTOFUpperBound[5][3];			// for each species, as well as the number of bins. The numbers
	Int_t					fTOFbins[5][3];					// size of the array is [ptrange][species].

	Double_t				fTPCLowerBound[5][3];			// The same, but now for TPC.
	Double_t				fTPCUpperBound[5][3];
	Int_t					fTPCbins[5][3];

// Naming conventions of the histograms.
	TString					fHistoName[12];					// Names of the histogram classes.
	TString					fHistoLatex[12];				// Names of the histogram classes in LaTeX.
	TString					fParticleName[3];				// Names of the particles (Pion, Kaon, Proton)
	TString					fPtClassName[5];				// Names of the ptclasses (should only be for internal use)

	Int_t 					fDebug;							// Debug flag.

	ClassDef(AliAODTrackCutsDiHadronPID,4);

};

#endif
Commit	Line	Data
6788af99	1	#ifndef ALIAODTRACKCUTSDIHADRONPID_H
	2	#define ALIAODTRACKCUTSDIHADRONPID_H
	3
	4	#include "TFormula.h"
	5	#include "TH1F.h"
	6	#include "TH2F.h"
	7	#include "TH3F.h"
	8	#include "TList.h"
	9
	10	class AliAODTrackCutsDiHadronPID : public TNamed
	11
	12	{
	13
	14	public:
	15	enum HistoClass {kAllCharged = 0, kPositive = 1, kNegative = 2,
	16	kAllPion = 3, kPosPion = 4, kNegPion = 5,
	17	kAllKaon = 6, kPosKaon = 7, kNegKaon = 8,
	18	kAllProton = 9, kPosProton = 10, kNegProton = 11};
	19
	20	public:
	21	AliAODTrackCutsDiHadronPID(); // Default Constructor
	22	AliAODTrackCutsDiHadronPID(const char* name); // Named Constructor
	23	virtual ~AliAODTrackCutsDiHadronPID(); // Destructor
	24	virtual Long64_t Merge(TCollection* list); // Merger
	25
f054df96	26	private:
	27	AliAODTrackCutsDiHadronPID(const AliAODTrackCutsDiHadronPID&);
	28	AliAODTrackCutsDiHadronPID& operator=(const AliAODTrackCutsDiHadronPID&);
	29
6788af99	30	// -------------------------------------------------------------------------
	31	// Interface, methods used to get information about the track cuts, and to
	32	// retrieve filled histograms.
	33	// -------------------------------------------------------------------------
	34
	35	public:
	36	void PrintCuts(); // Gives an overview of the cuts.
	37
	38	// Return the list of QA histos
	39	TList* GetListOfDataQAHistos() const{
	40	if (fDataTrackQAHistos) {return fDataTrackQAHistos;}
	41	else return 0x0;
	42	}
	43	TList* GetListOfPrimRecMCTrackQAHistos() const{
	44	if (fPrimRecMCTrackQAHistos) {return fPrimRecMCTrackQAHistos;}
	45	else return 0x0;
	46	}
	47	TList* GetListOfPrimGenMCTrackQAHistos() const {
	48	if (fPrimGenMCTrackQAHistos) {return fPrimGenMCTrackQAHistos;}
	49	else return 0x0;
	50	}
	51	TList* GetListOfSecRecMCTrackQAHistos() const {
	52	if (fSecRecMCTrackQAHistos) return fSecRecMCTrackQAHistos;
	53	else return 0x0;
	54	}
	55	TList* GetListOfSecGenMCTrackQAHistos() const {
	56	if (fSecGenMCTrackQAHistos) return fSecGenMCTrackQAHistos;
	57	else return 0x0;
	58	}
	59
	60	// Note that the PID histograms have in principle a different number of pT bins.
	61	// FIXME: This is not very nice...
	62	Int_t GetNPtBins() const {return fNPtBins;}
	63	Int_t GetNPtBinsPID(Int_t ptclass = -1) const {
	64	// if class = -1, then return the sum.
	65	if (ptclass == -1) {
	66	Int_t nptbinspid = 0;
	67	for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
	68	nptbinspid += fNPtBinsPID[iPtClass];
	69	}
	70	return nptbinspid;
	71	} else if (ptclass >= 0 && ptclass < 5) {
	72	return fNPtBinsPID[ptclass];
	73	} else {return -999;}
	74	}
	75
	76	// Returns the Pt axis for PID and other histograms.
	77	Double_t* GetPtAxis() {return fPtAxis;}
50dfda71	78	Double_t* GetPtAxisPID() const {
6788af99	79	const Int_t nptbinspid = GetNPtBinsPID();
50dfda71	80	Double_t* ptaxis = new Double_t[nptbinspid + 1];
6788af99	81	for (Int_t iPtBin = 0; iPtBin < nptbinspid; iPtBin++) {
	82	ptaxis[iPtBin] = GetPtMinPID(iPtBin + 1);
	83	}
	84	ptaxis[nptbinspid] = GetPtMaxPID(nptbinspid);
	85	return ptaxis;
	86	}
	87
	88	// Return data histogram with a specific name. Since the histograms are not streamed, and only the
	89	// TList containing them is, we have to retrieve them from the list.
	90	TObject* GetHistData(const char* name) const {return fDataTrackQAHistos->FindObject(name);}
	91	TObject* GetHistPrimRecMC(const char* name) const {return fPrimRecMCTrackQAHistos->FindObject(name);}
	92	TObject* GetHistPrimGenMC(const char* name) const {return fPrimGenMCTrackQAHistos->FindObject(name);}
	93	TObject* GetHistSecRecMC(const char* name) const {return fSecRecMCTrackQAHistos->FindObject(name);}
	94	TObject* GetHistSecGenMC(const char* name) const {return fSecGenMCTrackQAHistos->FindObject(name);}
	95
	96	// Since we will often want to have TOF histograms, here are a few methods which return the
	97	// appropriate projections. The class does not own these projections, and the user must take care of them.
	98	TH1F* GetHistDataTOFProjection(Int_t charge, Int_t species, Int_t ptbin);
	99	TH1F* GetHistDataTOFMismatch(Int_t charge, Int_t species, Int_t ptbin);
	100	Double_t GetPtMinPID(Int_t bin) const {
	101	Int_t ptclass = GetPtClass(bin);
	102	Int_t bininptclass = GetBinInPtClass(bin);
	103	Double_t minpt = fPtBoundaryPID[ptclass];
	104	Double_t maxpt = fPtBoundaryPID[ptclass+1];
	105	Double_t ptres = (maxpt - minpt)/((Double_t)fNPtBinsPID[ptclass]);
	106	return (minpt + ptres * ((Double_t)(bininptclass - 1)) );
	107	}
	108	Double_t GetPtMaxPID(Int_t bin) const {
	109	Int_t ptclass = GetPtClass(bin);
	110	Int_t bininptclass = GetBinInPtClass(bin);
	111	Double_t minpt = fPtBoundaryPID[ptclass];
	112	Double_t maxpt = fPtBoundaryPID[ptclass+1];
	113	Double_t ptres = (maxpt - minpt)/((Double_t)fNPtBinsPID[ptclass]);
	114	return (minpt + ptres * ((Double_t)(bininptclass)) );
	115	}
	116	Double_t GetPtClassMin(Int_t ptclass) const {
	117	if (ptclass >= 0 && ptclass < 5) {
	118	return fPtBoundaryPID[ptclass];
	119	} else {return -999;}
	120	}
	121	Double_t GetPtClassMax(Int_t ptclass) const {
	122	if (ptclass >= 0 && ptclass < 5) {
	123	return fPtBoundaryPID[ptclass+1];
	124	} else {return -999;}
	125	}
	126
	127	Double_t GetPtBinWidthPID(Int_t bin) const {return (GetPtMaxPID(bin) - GetPtMinPID(bin)); }
	128
	129	// BE CAREFUL! Following methods do not apply to the ptaxis of the PID histograms! For that, call
	130	// the GetPtMinPID and GetPtMaxPID methods.
	131	Double_t GetPtMin(Int_t bin) const {
	132	if ((bin < 1) \|\| (bin > fNPtBins + 1)) {cout<<"Bin is out of range..."<<endl; return -999.;}
	133	else {return fPtAxis[bin - 1];}
	134	}
	135	Double_t GetPtMax(Int_t bin) const {
	136	if ((bin < 1) \|\| (bin > fNPtBins + 1)) {cout<<"Bin is out of range..."<<endl; return -999.;}
	137	else {return fPtAxis[bin];}
	138	}
	139	Double_t GetPtBinWidth(Int_t bin) const {return (GetPtMax(bin) - GetPtMin(bin)); }
	140
	141	Int_t GetNTOFbins(Int_t ptclass, Int_t species) const {return fTOFbins[ptclass][species];}
	142	Double_t GetTOFmin(Int_t ptclass, Int_t species) const {return fTOFLowerBound[ptclass][species];}
	143	Double_t GetTOFmax(Int_t ptclass, Int_t species) const {return fTOFUpperBound[ptclass][species];}
	144
145	Int_t GetNTPCbins(Int_t ptclass, Int_t species) const {return fTPCbins[ptclass][species];}
146	Double_t GetTPCmin(Int_t ptclass, Int_t species) const {return fTPCLowerBound[ptclass][species];}
147	Double_t GetTPCmax(Int_t ptclass, Int_t species) const {return fTPCUpperBound[ptclass][species];}
148
149	// Getters (Cuts)
150	UInt_t GetFilterMask() const {return fFilterMask;}
151	Double_t GetMaxEta() const {return fMaxEta;}
152	ULong_t GetDemandedFlags() const {return fDemandedFlags;}
153	TFormula* GetPtDeptDCACutFormula() const {return fPtDeptDCAxyCutFormula;}
154	Double_t GetDCAzCut() const {return fDCAzCut;}
155	UInt_t GetMinSPDHitsForPtDeptDCACut() const {return fMinSPDHitsForPtDeptDCAcut;}
156
157	// Getters (Settings)
158	Bool_t GetIsMC() const {return fIsMC;}
159
160	Int_t GetDebugLevel() const {return fDebug;}
161
162	// -------------------------------------------------------------------------
163	// Methods used to configure the track cuts object, to be called at
164	// initialization, i.e., before the object is added to an analysis task.
165	// -------------------------------------------------------------------------
166
167	public:
168
169	// Request Certain QA histograms being filled.
170	Bool_t RequestQAHistos(Int_t histoclass, Bool_t Enable3DSpectra = kFALSE, Bool_t EnablePIDHistos = kFALSE) {
171	if ((histoclass > -1) && (histoclass < 12)) {
172	fHistRequested[histoclass] = kTRUE;
173	f3DSpectraEnabeled[histoclass] = Enable3DSpectra;
174	fPIDHistosEnabeled[histoclass] = EnablePIDHistos;
175	//cout<<"histoclass: "<<histoclass<<" requested: "<<fHistRequested[histoclass]<<endl;
176	return kTRUE;
177	} else {
178	return kFALSE;
179	}
180	}
181
182	// Setters (Cuts)
183	void SetPtRange(Double_t minpt, Double_t maxpt) {
184	fMinPt = minpt;
185	fMaxPt = maxpt;
186	fTestPt = kTRUE;
187	}
188	void SetFilterMask(UInt_t filtermask) {
189	fFilterMask = filtermask;
190	fTestFilterMask = kTRUE;
191	}
192	void SetMaxEta(Double_t maxeta) {
193	fMaxEta = maxeta;
194	fTestMaxEta = kTRUE;
195	}
196	void SetMaxRapidity(Double_t maxrapidity) {
197	fMaxRapidity = maxrapidity;
198	fTestMaxRapidity = kTRUE;
199	}
200	void SetDemandNoMismatch() {
201	fTestTOFmismatch = kTRUE;
202	}
203	void SetDemandFlags(ULong_t demandedflags) {
204	fDemandedFlags = demandedflags;
205	fTestFlags = kTRUE;
206	}
207	void SetPtDeptDCACut(TFormula* DCAxyCutFormula, Double_t DCAzCut, UInt_t MinSPDHits = 1) {
208	fPtDeptDCAxyCutFormula = DCAxyCutFormula;
209	fDCAzCut = DCAzCut;
210	fMinSPDHitsForPtDeptDCAcut = MinSPDHits;
211	fTestPtDeptDCAcut = kTRUE;
212	}
213
214	// Setters (Settings)
215	void SetIsMC(Bool_t ismc = kTRUE) {fIsMC = ismc;}
216
217	void SetDebugLevel(Int_t debuglevel) {fDebug = debuglevel;}
218
219	// -------------------------------------------------------------------------
220	// Methods called by the analysis task.
221	// -------------------------------------------------------------------------
222
223	public:
224
225	// These two functions signal the beginning and the end of a new event.
226	void StartNewEvent(); // Some things are set to zero.
227	void EventIsDone(Bool_t IsMC); // Some final histograms are filled.
228	void CreateHistos(); // Should be called by the UserCreateOutput() function of the analysis task.
229
230	// Is Selected, for different types of tracks.
231	Bool_t IsSelectedData(AliTrackDiHadronPID* track, Double_t randomhittime = -1.e20);
232	Bool_t IsSelectedGeneratedMC(AliAODMCParticle* particle);
233	Bool_t IsSelectedReconstructedMC(AliTrackDiHadronPID* track);
234
235	// -------------------------------------------------------------------------
236	// Internal methods.
237	// -------------------------------------------------------------------------
238
239	public:
240
241	// For PID histograms we have a certain number of bins in pT, spread out over five
242	// large histograms, i.e., one for the lowest pT, and the biggest range in TOF/TPC,
243	// one for the higher pT and smaller range in TOF/TPC, etc. The following methods
244	// are a mapping between the total pT bin (what the user uses), and the pt bin
245	// within one of the five histograms (what's used internally)
246	Int_t GetPtClass(const Int_t ptbin) const {
247
248	// Returns a number [0..4]
249	Int_t currentptclass = 0;
250	Int_t currentptbin = fNPtBinsPID[0];
251	while (currentptbin < ptbin) {
252	currentptclass++;
253	if (currentptclass == 5) {break;}
254	currentptbin += fNPtBinsPID[currentptclass];
255	}
256	if (currentptclass == 5) {cout<<"GetPtClass -> ptbin out of range!"<<endl; return -1;}
257	return currentptclass;
258	}
259	Int_t GetBinInPtClass(const Int_t ptbin) const {
260
261	// Returns a number [1..Nbins]
262	Int_t ptclass = GetPtClass(ptbin);
263	if (ptclass == -1) {return -1;}
264
265	Int_t ptbinout = ptbin;
266	for (Int_t iPtClass = 0; iPtClass < ptclass; iPtClass++) {ptbinout -= fNPtBinsPID[iPtClass];}
267
268	return ptbinout;
269
270	}
271
272	private:
273
274	// Checks, return kTRUE if track passes the cut.
275	Bool_t CheckPt(Double_t pt) const {
276	if (!fTestPt) return kTRUE;
277	if ((pt > fMinPt) && (pt < fMaxPt)) return kTRUE;
278	return kFALSE;
279	}
280	Bool_t CheckMaxEta(Double_t eta) const {
281	if (!fTestMaxEta) return kTRUE; // Accepted if there is no check on this parameter.
282	if (TMath::Abs(eta) < fMaxEta) return kTRUE;
283	return kFALSE;
284	}
285	Bool_t CheckRapidity(Double_t rap) const {
286	if (!fTestMaxRapidity) return kTRUE;
287	if (TMath::Abs(rap) < fMaxRapidity) return kTRUE;
288	return kFALSE;
289	}
290	Bool_t CheckFilterMask(UInt_t filtermap) const {
291	if (!fTestFilterMask) return kTRUE;
292	if ((fFilterMask & filtermap) == fFilterMask) return kTRUE;
293	return kFALSE;
294	}
295	Bool_t CheckFlags(ULong_t flags) const {
296	if (!fTestFlags) return kTRUE;
297	if ((flags & fDemandedFlags) == fDemandedFlags) return kTRUE;
298	return kFALSE;
299	}
300	Bool_t CheckTOFmismatch(Bool_t ismismatch) const {
301	if (!fTestTOFmismatch) return kTRUE; // if we're not cutting on mismatch, then it's accepted.
302	if (!ismismatch) return kTRUE; // so if the track is not a mismatch, then it is accepted.
303	return kFALSE; // if it is a mismatch, then it's not accepted.
304	}
305	Bool_t CheckPtDeptDCACut(Double_t dcaz, Double_t dcaxy, Double_t pt, UInt_t SPDhits) const {
306	if (!fTestPtDeptDCAcut) return kTRUE;
307	if (SPDhits < fMinSPDHitsForPtDeptDCAcut) return kTRUE; // If there are not enough SPD hits to do the cut.
308	if ((dcaz < fDCAzCut) && (dcaxy < fPtDeptDCAxyCutFormula->Eval(pt))) return kTRUE;
309	return kFALSE;
310	}
311
312	// Filling QA histograms.
313	Bool_t FillDataHistos(Int_t histoclass, AliTrackDiHadronPID* track);
314	Bool_t FillTOFMismatchHistos(Int_t histoclass, AliTrackDiHadronPID* track, Double_t randomhittime);
315	Bool_t FillGenMCHistos(Int_t histoclass, AliAODMCParticle* particle);
316	Bool_t FillRecMCHistos(Int_t histoclass, AliTrackDiHadronPID* track);
317
318	// Initializing QA histograms.
319	Bool_t InitializeDataHistos(Int_t histoclass);
320	Bool_t InitializeGenMCHistos(Int_t histoclass);
321	Bool_t InitializeRecMCHistos(Int_t histoclass);
322
323	void InitializeDefaultHistoNamesAndAxes();
324
325	TH1F* InitializePtSpectrum(const char* name, Int_t histoclass);
50dfda71	326	TH3F* InitializePhiEtaPt(const char* name, Int_t histoclass);
6788af99	327	TH1F* InitializeNTracksHisto(const char* name, Int_t histoclass);
	328	TH1F* InitializeDCAxyHisto(const char* name, Int_t histoclass);
	329	TH1F* InitializeDCAzHisto(const char* name, Int_t histoclass);
	330	TH3F* InitializeAcceptanceHisto(const char* /name/, Int_t /histoclass/); // TO BE IMPLEMENTED.
	331	TH2F* InitializeDCASpectrum(const char* name, Int_t histoclass);
	332
	333
	334	TH3F* InitializePIDHisto(const char* name, Int_t histoclass, Int_t expspecies, Int_t ptclass);
	335	TH2F* InitializeTOFMismatchHisto(const char* name, Int_t histoclass, Int_t expspecies, Int_t ptclass);
	336
	337	// -------------------------------------------------------------------------
	338	// Data members.
	339	// -------------------------------------------------------------------------
	340
	341	private:
	342	// Track Cuts
	343	Double_t fMinPt;
	344	Double_t fMaxPt;
	345	UInt_t fFilterMask; // FilterMask to-be-checked.
	346	Double_t fMaxEta; // Max Eta of the track.
	347	Double_t fMaxRapidity; // Rapidity cut (only done for PID plots!!)
	348	ULong_t fDemandedFlags; // Flags demanded on the track.
	349	UInt_t fMinSPDHitsForPtDeptDCAcut; // Required number of SPD hits for performing Pt-Dept DCA cut.
	350	TFormula* fPtDeptDCAxyCutFormula; // Formula for the Pt-Dept DCA cut.
	351	Double_t fDCAzCut; // Max z at DCA.
	352
	353	// Settings
	354	Bool_t fIsMC; // Is the current event MC or not.
	355
	356	// Requested Histograms;
	357	Bool_t fHistRequested[12]; //
	358	Bool_t f3DSpectraEnabeled[12]; //
	359	Bool_t fPIDHistosEnabeled[12]; //
	360
	361	// Which Track Cuts will be tested.
	362	Bool_t fTestPt; //
	363	Bool_t fTestFilterMask; //
	364	Bool_t fTestMaxEta; //
	365	Bool_t fTestMaxRapidity; //
	366	Bool_t fTestFlags; //
	367	Bool_t fTestTOFmismatch; //
	368	Bool_t fTestPtDeptDCAcut; //
	369
	370	// QA histograms for Data.
	371	TList* fDataTrackQAHistos; //
	372	TH1F* fHistDataPt[3]; //! Pt distribution of tracks passing this cut.
50dfda71	373	TH3F* fHistDataPhiEtaPt[3]; //! Pt, Eta, Phi distribution.
6788af99	374	TH1F* fHistDataNTracks[3]; //! Number of tracks passing the cut per event (filling by EventIsDone()).
	375	TH1F* fHistDataDCAxy[3]; //! DCA_{xy} distribution.
	376	TH1F* fHistDataDCAz[3]; //! DCA_{z} distribution
	377	TH2F* fHistDataDCAxyOneSigma[12]; //! DCA_{xy} distribution of particles as identified by 1 sigma method.
	378	Int_t fNTracks[3]; //! Number of tracks
	379
	380	TH3F* fHistDataPID[3][3][5]; //! TPC/TOF v.s. pT, [charge][mass assumption][ptclass]
	381	TH2F* fHistTOFMismatch[3][3][5]; //! TOF Mismatch histograms, [charge][mass assumption][ptclass]
	382
	383	// QA histograms for Primary Reconstructed MC tracks.
	384	TList* fPrimRecMCTrackQAHistos; //
	385	TH1F* fHistPrimRecMCPt[12]; //! Pt distribution of reconstructed MC track passing this cut.
	386	TH1F* fHistPrimRecNTracks[12]; //!
	387	TH2F* fHistPrimRecMCDCA[12]; //! DCA_xy distribution of reconstructed MC track passing this cut.
	388
	389	// QA histograms for Primary Generated MC particles.
	390	TList* fPrimGenMCTrackQAHistos; //
	391	TH1F* fHistPrimGenMCPt[12]; //! Pt distribution of generated MC particles passing this cut.
	392
	393	// QA histograms for Secondary Reconstructed MC tracks.
	394	TList* fSecRecMCTrackQAHistos; //
	395	TH1F* fHistSecRecMCPt[12]; //! Pt distribution of reconstructed MC track passing this cut.
	396	TH2F* fHistSecRecMCDCAMat[12]; //! DCA_xy distribution of material decay particles.
	397	TH2F* fHistSecRecMCDCAWeak[12]; //! DCA_xy distribution of weak decay.
	398
	399	// QA histograms for Secondary Generated MC particles.
	400	TList* fSecGenMCTrackQAHistos; //
	401	TH1F* fHistSecGenMCPt[12]; //! Pt distribution of generated MC particles passing this cut.
	402
	403	// Binning of all the histograms.
50dfda71	404	Double_t fPtAxis[57]; // Pt axis used in all histograms, except PID and Mismatch histograms.
6788af99	405	Int_t fNPtBins; // Number of bins in the pt-axis.
50dfda71	406	Int_t fNEtaBins; //
50dfda71	407	Int_t fNPhiBins; //
6788af99	408
	409	Double_t fPtBoundaryPID[6]; // There are five different PID histo's. This array gives the pT range of these histograms.
	410	Int_t fNPtBinsPID[5]; // This array gives the number of pT bins for each of these histograms.
	411
	412	Double_t fTOFLowerBound[5][3]; // These arrays give the lower and upper bound of the TOF axes,
	413	Double_t fTOFUpperBound[5][3]; // for each species, as well as the number of bins. The numbers
	414	Int_t fTOFbins[5][3]; // size of the array is [ptrange][species].
	415
	416	Double_t fTPCLowerBound[5][3]; // The same, but now for TPC.
	417	Double_t fTPCUpperBound[5][3];
	418	Int_t fTPCbins[5][3];
	419
	420	// Naming conventions of the histograms.
	421	TString fHistoName[12]; // Names of the histogram classes.
	422	TString fHistoLatex[12]; // Names of the histogram classes in LaTeX.
	423	TString fParticleName[3]; // Names of the particles (Pion, Kaon, Proton)
	424	TString fPtClassName[5]; // Names of the ptclasses (should only be for internal use)
	425
	426	Int_t fDebug; // Debug flag.
	427
	428	ClassDef(AliAODTrackCutsDiHadronPID,4);
	429
	430	};
	431
	432	#endif