[u/mrichter/AliRoot.git] / PWGCF / Correlations / DPhi / DiHadronPID / AliAnalysisTaskDiHadronPID.cxx

// -------------------------------------------------------------------------
//  INFO
// -------------------------------------------------------------------------

#include <iostream>

// Basic Includes
#include "TH1F.h"
#include "TH2F.h"
#include "TH3F.h"
#include "THn.h"
#include "TFile.h"
#include "TChain.h"
#include "TObject.h"
#include "TObjArray.h"

// Manager/ Handler
#include "AliAnalysisManager.h"
#include "AliInputEventHandler.h"

// Event pool includes.
#include "AliEventPoolManager.h"

// PID includes.
#include "AliPIDResponse.h"

// AOD includes.
#include "AliAODEvent.h"
#include "AliAODTrack.h"
#include "AliAODHandler.h"
#include "AliAODVertex.h"
#include "AliAODInputHandler.h"
#include "AliAODMCParticle.h"
#include "AliAODMCHeader.h"

// Additional includes.
#include "AliTrackDiHadronPID.h"
#include "AliAODTrackCutsDiHadronPID.h"
#include "AliAODEventCutsDiHadronPID.h"
#include "AliHistToolsDiHadronPID.h"

// AnalysisTask headers.
#include "AliAnalysisTaskSE.h"
#include "AliAnalysisTaskDiHadronPID.h"

using namespace std;

ClassImp(AliAnalysisTaskDiHadronPID);

// -------------------------------------------------------------------------
AliAnalysisTaskDiHadronPID::AliAnalysisTaskDiHadronPID():
	AliAnalysisTaskSE(),
	fPIDResponse(0x0),
	fEventCuts(0x0),
	fTrackCutsTrigger(0x0),
	fTrackCutsAssociated(0x0),
	fPoolMgr(0x0),
	fTriggerTracks(0x0),
	fAssociatedTracks(0x0),
	fCurrentAODEvent(0x0),
	fOutputList(0x0),
	fPtSpectrum(0x0),
	fCorrelations(0x0),
	fMixedEvents(0x0),
	fTOFhistos(0x0),
	fNDEtaBins(32),
	fNDPhiBins(32),	
	fMinNEventsForMixing(5),
	fPoolTrackDepth(2000),
	fPoolSize(1000),
	fCalculateTOFmismatch(kTRUE),
	fT0Fill(0x0),
	fLvsEta(0x0),
	fLvsEtaProjections(0x0),		
	fDebug(0)

{

	//
	// Default Constructor.
	//

	if (fDebug > 0) {AliInfo("AliAnalysisTaskDiHadronPID Default Constructor.");}		

	for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
			for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
			fCorrelationsTOF[iPtClass][iSpecies] = 0x0;
			fCorrelationsTOFTPC[iPtClass][iSpecies] = 0x0;
		}
	}


}

// -------------------------------------------------------------------------
AliAnalysisTaskDiHadronPID::AliAnalysisTaskDiHadronPID(const char* name):
	AliAnalysisTaskSE(name),
	fPIDResponse(0x0),
	fEventCuts(0x0),
	fTrackCutsTrigger(0x0),
	fTrackCutsAssociated(0x0),
	fPoolMgr(0x0),
	fTriggerTracks(0x0),
	fAssociatedTracks(0x0),	
	fCurrentAODEvent(0x0),
	fOutputList(0x0),
	fPtSpectrum(0x0),
	fCorrelations(0x0),
	fMixedEvents(0x0),
	fTOFhistos(0x0),	
	fNDEtaBins(32),
	fNDPhiBins(32),
	fMinNEventsForMixing(5),
	fPoolTrackDepth(2000),
	fPoolSize(1000),
	fCalculateTOFmismatch(kTRUE),
	fT0Fill(0x0),
	fLvsEta(0x0),
	fLvsEtaProjections(0x0),						
	fDebug(0) 

{

	//
	// Named Constructor.
	//

	if (fDebug > 0) {AliInfo("AliAnalysisTaskDiHadronPID Named Constructor.");}	
	if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}

	for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
			for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
			fCorrelationsTOF[iPtClass][iSpecies] = 0x0;
			fCorrelationsTOFTPC[iPtClass][iSpecies] = 0x0;
		}
	}

	DefineInput(0,TChain::Class());
	DefineOutput(1,TList::Class());

}

// -------------------------------------------------------------------------
AliAnalysisTaskDiHadronPID::~AliAnalysisTaskDiHadronPID() {;

	//
	// Destructor.
	//

	if (fDebug > 0) {AliInfo("AliAnalysisTaskDiHadronPID Destructor.");}
	if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}

}

// -------------------------------------------------------------------------
void AliAnalysisTaskDiHadronPID::UserCreateOutputObjects() {

	//
	// Create Output objects.
	//

	if (fDebug > 0) {AliInfo("UserCreateOutputObjects()");}
	if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}

	// --- BEGIN: Initialization on the worker nodes ---
	AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
	AliInputEventHandler* inputHandler = dynamic_cast<AliInputEventHandler*> (manager->GetInputEventHandler());

	// Getting the pointer to the PID response object.
	fPIDResponse = inputHandler->GetPIDResponse();	

	// Not very neat - only set up for 0-5% analysis.
	Int_t nCentralityBins  = 5;
	Double_t centralityBins[] = {0.,1.,2.,3.,4.,5.};

	Int_t nZvtxBins  = 7;
	Double_t vertexBins[] = {-7.,-5.,-3.,-1.,1.,3.,5.,7.};

	fPoolMgr = new AliEventPoolManager(fPoolSize, fPoolTrackDepth, nCentralityBins, (Double_t*) centralityBins, nZvtxBins, (Double_t*) vertexBins);
    // --- END ---

	// Create the output list.
	fOutputList = new TList();
	fOutputList->SetOwner(kTRUE);

	// Creating all requested histograms locally.
	fEventCuts->CreateHistos();
	fOutputList->Add(fEventCuts);

	fTrackCutsTrigger->CreateHistos();
	fOutputList->Add(fTrackCutsTrigger);

	fTrackCutsAssociated->CreateHistos();
	fOutputList->Add(fTrackCutsAssociated);

	// Get the pT axis for the PID histograms.
	Float_t* ptaxis = fTrackCutsAssociated->GetPtAxisPID();
	Int_t nptbins = fTrackCutsAssociated->GetNPtBinsPID();

	// Create Pt spectrum histogram.
	fPtSpectrum = new TH1F("fPtSpectrum","p_{T} Spectrum;p_{T} (GeV/c);Count",nptbins,ptaxis);
	fOutputList->Add(fPtSpectrum);

	// Create unidentified correlations histogram.
	fCorrelations = AliHistToolsDiHadronPID::MakeHist3D("fCorrelations","Correlations;#Delta#phi;#Delta#eta;p_{T} (GeV/c)",
		fNDPhiBins,-TMath::Pi()/2.,3.*TMath::Pi()/2.,
		fNDEtaBins,-1.6,1.6,
		nptbins, ptaxis);
	fOutputList->Add(fCorrelations);

	// Create unidentified mixed events histogram.
	fMixedEvents = AliHistToolsDiHadronPID::MakeHist3D("fMixedEvents","Mixed Events;#Delta#phi;#Delta#eta;p_{T} (GeV/c)",
		fNDPhiBins,-TMath::Pi()/2.,3.*TMath::Pi()/2.,
		fNDEtaBins,-1.6,1.6,
		nptbins, ptaxis);
	fOutputList->Add(fMixedEvents);

	// Create TOF correlations histograms (DPhi,DEta,pt,TOF).
	fTOFhistos = new TObjArray(15);
	fTOFhistos->SetName("CorrelationsTOF");

	Int_t nbins[4] = {fNDPhiBins,fNDEtaBins,0.,0.};
	Double_t min[4] = {-TMath::Pi()/2.,-1.6,0.,0.};
	Double_t max[4] = {3.*TMath::Pi()/2.,1.6,0.,0.};

	for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
		
		nbins[2] = fTrackCutsAssociated->GetNPtBinsPID(iPtClass);
		min[2] = fTrackCutsAssociated->GetPtClassMin(iPtClass);
		max[2] = fTrackCutsAssociated->GetPtClassMax(iPtClass);

		for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {

			nbins[3] = fTrackCutsAssociated->GetNTOFbins(iPtClass,iSpecies);
			min[3] = fTrackCutsAssociated->GetTOFmin(iPtClass,iSpecies);
			max[3] = fTrackCutsAssociated->GetTOFmax(iPtClass,iSpecies);

			fCorrelationsTOF[iPtClass][iSpecies] = new THnF(
				Form("fCorrelationsTOF_%i_%i",iPtClass,iSpecies),
				Form("CorrelationsTOF_%i_%i",iPtClass,iSpecies),
				4,nbins,min,max);

			fTOFhistos->Add(fCorrelationsTOF[iPtClass][iSpecies]);

		}
	}

	fOutputList->Add(fTOFhistos);

	// TODO: Create TOF/TPC correlations histogram.

	// Load external TOF histograms if flag is set.
	if (fCalculateTOFmismatch) {LoadExtMismatchHistos();}

	PostData(1,fOutputList);

}

// -------------------------------------------------------------------------
void AliAnalysisTaskDiHadronPID::LocalInit() {

	//
	// Initialize on the client. (or on my computer?? - I think so...)
	//

	if (fDebug > 0) {AliInfo("LocalInit()");}
	if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 	
}

// -------------------------------------------------------------------------
void AliAnalysisTaskDiHadronPID::UserExec(Option_t*) {

	//
	// Main Loop.
	//

	if (fDebug > 0) {AliInfo("UserExec()");}
	if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}

	// Input Current Event.
	fCurrentAODEvent = dynamic_cast<AliAODEvent*>(InputEvent());
	if (!fCurrentAODEvent) AliFatal("No Event Found!");

	if (!fEventCuts->IsSelected(fCurrentAODEvent)) {return;}

	// Fill the global tracks array. - NOT NEEDED I THINK, since we're not using
	// bit 1<<7 for the associated tracks!

	// Let the track cut objects know that a new event will start.
	fTrackCutsTrigger->StartNewEvent();
	fTrackCutsAssociated->StartNewEvent();

	// Create arrays for trigger/associated particles.
	fTriggerTracks = new TObjArray(350);
	fTriggerTracks->SetOwner(kTRUE);

	fAssociatedTracks = new TObjArray(3500);
	fAssociatedTracks->SetOwner(kTRUE);

	for (Int_t iTrack = 0; iTrack < fCurrentAODEvent->GetNumberOfTracks(); iTrack++) {

		AliAODTrack* track = (AliAODTrack*)fCurrentAODEvent->GetTrack(iTrack);
		AliTrackDiHadronPID* pidtrack = new AliTrackDiHadronPID(track,0x0,0x0,fPIDResponse);
		pidtrack->ForgetAboutPointers();
		pidtrack->SetDebugLevel(fDebug);

		Double_t rndhittime = -1.e21;
		if (fCalculateTOFmismatch) rndhittime = GenerateRandomHit(pidtrack->Eta());

		// Fill p_T spectrum.
		fPtSpectrum->Fill(pidtrack->Pt());

		// Fill the trigger/associated tracks array.
		if (fTrackCutsTrigger->IsSelectedData(pidtrack,rndhittime)) {fTriggerTracks->AddLast(pidtrack);}
		else if (fTrackCutsAssociated->IsSelectedData(pidtrack,rndhittime)) {fAssociatedTracks->AddLast(pidtrack);} 
		else {delete pidtrack;}

	}

	// Fill Correlation histograms.
	for (Int_t iTrigger = 0; iTrigger < fTriggerTracks->GetEntriesFast(); iTrigger++) {
		AliTrackDiHadronPID* triggertrack = (AliTrackDiHadronPID*)fTriggerTracks->At(iTrigger);

		for (Int_t iAssociated = 0; iAssociated < fAssociatedTracks->GetEntriesFast(); iAssociated++) {
			AliTrackDiHadronPID* associatedtrack = (AliTrackDiHadronPID*)fAssociatedTracks->At(iAssociated);

			Double_t DPhi = triggertrack->Phi() - associatedtrack->Phi();
			if (DPhi < -TMath::Pi()/2.) {DPhi += 2.*TMath::Pi();}
			if (DPhi > 3.*TMath::Pi()/2.) {DPhi -= 2.*TMath::Pi();}

			Double_t DEta = triggertrack->Eta() - associatedtrack->Eta();
			fCorrelations->Fill(DPhi,DEta,associatedtrack->Pt());

			Double_t tofhistfill[4] = {DPhi,DEta,associatedtrack->Pt(),-999.};

			for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
				tofhistfill[3] = associatedtrack->GetTOFsignalMinusExpected(iSpecies);

				for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {

					// prevent under/over-flow bins to be filled.
					Double_t ptmin = fTrackCutsAssociated->GetPtClassMin(iPtClass);
					Double_t ptmax = fTrackCutsAssociated->GetPtClassMax(iPtClass);
					Double_t apt = associatedtrack->Pt();

					if ( (apt > ptmin) && (apt < ptmax) ) {
						fCorrelationsTOF[iPtClass][iSpecies]->Fill(tofhistfill);
					}

				}
			}
		}
	}

	cout<<"Triggers: "<<fTriggerTracks->GetEntriesFast()<<" Associateds: "<<fAssociatedTracks->GetEntriesFast()<<endl;	

	// Determine centrality of current event.
	TString centralityestimator = fEventCuts->GetCentralityEstimator();
	AliCentrality* currentcentrality = fCurrentAODEvent->GetCentrality();
	Float_t percentile = currentcentrality->GetCentralityPercentile(centralityestimator.Data());

	// Determine vtxz of current event.
	AliAODVertex* currentprimaryvertex = fCurrentAODEvent->GetPrimaryVertex();
	Double_t vtxz = currentprimaryvertex->GetZ();

	// Obtain event pool for current centrality/ vtxz.
	AliEventPool* poolin = fPoolMgr->GetEventPool(percentile, vtxz); 
	if (!poolin) {AliFatal(Form("No pool found for centrality = %f, vtxz = %f", percentile, vtxz));}
	// TObjArray* fGlobalTracksArray; 

	// Mix events if there are enough events in the pool. (TODO: should be a data member.)
	if (poolin->GetCurrentNEvents() >= fMinNEventsForMixing) {
		if (fDebug) {AliInfo("Mixing Events.");}

		// Loop over all triggers in this event.
		for (Int_t iTrigger = 0; iTrigger < fTriggerTracks->GetEntriesFast(); iTrigger++) {
			AliTrackDiHadronPID* triggertrack = (AliTrackDiHadronPID*)fTriggerTracks->At(iTrigger);

			// Loop over all events in the event pool.
			for (Int_t iMixEvent = 0; iMixEvent < poolin->GetCurrentNEvents(); iMixEvent++) {
	    		TObjArray* mixtracks = poolin->GetEvent(iMixEvent);

	    		// Loop over all mixed tracks.
	    		for (Int_t iMixTrack = 0; iMixTrack < mixtracks->GetEntriesFast(); iMixTrack++) {
	    			AliTrackDiHadronPID* mixtrack = (AliTrackDiHadronPID*)mixtracks->At(iMixTrack);
						
					Double_t DPhi = triggertrack->Phi() - mixtrack->Phi();
					if (DPhi < -TMath::Pi()/2.) {DPhi += 2.*TMath::Pi();}
					if (DPhi > 3.*TMath::Pi()/2.) {DPhi -= 2.*TMath::Pi();}

					Double_t DEta = triggertrack->Eta() - mixtrack->Eta();
					fMixedEvents->Fill(DPhi,DEta,mixtrack->Pt());

	    		}
	   		}
	   	}
	}	

	// Update the event pool.
	AliEventPool* poolout = fPoolMgr->GetEventPool(percentile, vtxz); // Get the buffer associated with the current centrality and z-vtx
	if (!poolout) AliFatal(Form("No pool found for centrality = %f, vtx_z = %f", percentile, vtxz));

	// Q: is it a problem that the fAssociatedTracks array can be bigger than the number of tracks inside?
	poolout->UpdatePool(fAssociatedTracks);

	// Delete trigger array and its contents. Set pointer to zero.
	// Don't delete the associated array, since ownership has been transferred to the pool manager. Set pointer to zero.
	fTriggerTracks->Delete();
	delete fTriggerTracks;
	fTriggerTracks = 0x0;
	fAssociatedTracks = 0x0;

	// Tell the track cut object that the event is done.
	fTrackCutsTrigger->EventIsDone(0);
	fTrackCutsAssociated->EventIsDone(0);

	PostData(1,fOutputList);

}

// -------------------------------------------------------------------------
Bool_t AliAnalysisTaskDiHadronPID::LoadExtMismatchHistos() {

	//
	// Attempting to load a root file containing information needed
	// to generate random TOF hits.
 	//

	if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}

	// Opening external TOF file.
	if (fDebug > 0) cout<<"Trying to open TOFmismatchHistos.root ..."<<endl;
	TFile* fin = 0x0;
	fin = TFile::Open("alien:///alice/cern.ch/user/m/mveldhoe/rootfiles/TOFmismatchHistos.root");
	if (!fin) {
		AliWarning("Couln't open TOFmismatchHistos, will not calculate mismatches...");
		fCalculateTOFmismatch = kFALSE;
		return kFALSE;
	}

	// Check if the required histograms are present.
	TH1F* tmp1 = (TH1F*)fin->Get("hNewT0Fill");
	if (!tmp1) {
		AliWarning("Couln't find hNewT0Fill, will not calculate mismatches...");
		fCalculateTOFmismatch = kFALSE;
		return kFALSE;	
	}
	TH2F* tmp2 = (TH2F*)fin->Get("hLvsEta");
	if (!tmp2) {
		AliWarning("Couln't find hLvsEta, will not calculate mismatches...");
		fCalculateTOFmismatch = kFALSE;
		return kFALSE;	
	}	

	// Make a deep copy of the files in the histogram.
	fT0Fill = (TH1F*)tmp1->Clone("fT0Fill");
	fLvsEta = (TH2F*)tmp2->Clone("fLvsEta");

	// Close the external file.
	AliInfo("Closing external file.");
	fin->Close();

	// Creating a TObjArray for LvsEta projections.
	const Int_t nbinseta = fLvsEta->GetNbinsX();
	fLvsEtaProjections = new TObjArray(nbinseta);
	fLvsEtaProjections->SetOwner(kTRUE);

	// Making the projections needed (excluding underflow/ overflow).
	for (Int_t iEtaBin = 1; iEtaBin < (nbinseta + 1); iEtaBin++) {
		TH1F* tmp = (TH1F*)fLvsEta->ProjectionY(Form("LvsEtaProjection_%i",iEtaBin),iEtaBin,iEtaBin);
		tmp->SetDirectory(0);
		fLvsEtaProjections->AddAt(tmp,iEtaBin);
	}

	return kTRUE;

}

// -------------------------------------------------------------------------
Double_t AliAnalysisTaskDiHadronPID::GenerateRandomHit(Double_t eta) {

	//
	// Returns a random TOF time.
	//

	if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}

	// Default (error) value:
	Double_t rndhittime = -1.e21;

	// TOF mismatch flag is not turned on.
	if (!fCalculateTOFmismatch) {
		AliFatal("Called GenerateRandomHit() method, but flag fCalculateTOFmismatch not set.");
		return rndhittime;
	}

	// TOF doesn't extend much further than 0.8.
	if (TMath::Abs(eta) > 0.8) {
		if (fDebug) {AliInfo("Tried to get a random hit for a track with eta > 0.8.");}
		return rndhittime;
	}

	// Finding the bin of the eta.
	TAxis* etaAxis = fLvsEta->GetXaxis();
	Int_t etaBin = etaAxis->FindBin(eta);
	//cout<<"Eta: "<<eta<<" bin: "<<etaBin<<endl;
	const TH1F* lengthDistribution = (const TH1F*)fLvsEtaProjections->At(etaBin);

	if (!lengthDistribution) {
		AliFatal("length Distribution not found.");
		return rndhittime;
	}

	Double_t currentRndLength = lengthDistribution->GetRandom(); // in cm.

	// Similar to Roberto's code.
	Double_t currentRndTime = currentRndLength / (TMath::C() * 1.e2 / 1.e12);
	Double_t t0fill = -1.26416e+04;
	rndhittime = fT0Fill->GetRandom() - t0fill + currentRndTime;

	return rndhittime;

}

// -------------------------------------------------------------------------
void AliAnalysisTaskDiHadronPID::Terminate(Option_t*) {;

	//
	// Called when task is done.
	//

	if (fDebug > 0) {AliInfo("Terminate()");}
	if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}

	delete fT0Fill;
	fT0Fill = 0x0;
	delete fLvsEta;
	fLvsEta = 0x0;
	delete fLvsEtaProjections;
	fLvsEtaProjections = 0x0;

}
Commit	Line	Data
6788af99	1	// -------------------------------------------------------------------------
	2	// INFO
	3	// -------------------------------------------------------------------------
	4
	5	#include <iostream>
	6
	7	// Basic Includes
	8	#include "TH1F.h"
	9	#include "TH2F.h"
	10	#include "TH3F.h"
	11	#include "THn.h"
	12	#include "TFile.h"
	13	#include "TChain.h"
	14	#include "TObject.h"
	15	#include "TObjArray.h"
	16
	17	// Manager/ Handler
	18	#include "AliAnalysisManager.h"
	19	#include "AliInputEventHandler.h"
	20
	21	// Event pool includes.
	22	#include "AliEventPoolManager.h"
	23
	24	// PID includes.
	25	#include "AliPIDResponse.h"
	26
	27	// AOD includes.
	28	#include "AliAODEvent.h"
	29	#include "AliAODTrack.h"
	30	#include "AliAODHandler.h"
	31	#include "AliAODVertex.h"
	32	#include "AliAODInputHandler.h"
	33	#include "AliAODMCParticle.h"
	34	#include "AliAODMCHeader.h"
	35
	36	// Additional includes.
	37	#include "AliTrackDiHadronPID.h"
	38	#include "AliAODTrackCutsDiHadronPID.h"
	39	#include "AliAODEventCutsDiHadronPID.h"
	40	#include "AliHistToolsDiHadronPID.h"
	41
	42	// AnalysisTask headers.
	43	#include "AliAnalysisTaskSE.h"
	44	#include "AliAnalysisTaskDiHadronPID.h"
	45
	46	using namespace std;
	47
	48	ClassImp(AliAnalysisTaskDiHadronPID);
	49
	50	// -------------------------------------------------------------------------
	51	AliAnalysisTaskDiHadronPID::AliAnalysisTaskDiHadronPID():
	52	AliAnalysisTaskSE(),
	53	fPIDResponse(0x0),
	54	fEventCuts(0x0),
	55	fTrackCutsTrigger(0x0),
	56	fTrackCutsAssociated(0x0),
	57	fPoolMgr(0x0),
	58	fTriggerTracks(0x0),
	59	fAssociatedTracks(0x0),
	60	fCurrentAODEvent(0x0),
	61	fOutputList(0x0),
	62	fPtSpectrum(0x0),
	63	fCorrelations(0x0),
	64	fMixedEvents(0x0),
65	fTOFhistos(0x0),
66	fNDEtaBins(32),
67	fNDPhiBins(32),
68	fMinNEventsForMixing(5),
69	fPoolTrackDepth(2000),
70	fPoolSize(1000),
71	fCalculateTOFmismatch(kTRUE),
72	fT0Fill(0x0),
73	fLvsEta(0x0),
74	fLvsEtaProjections(0x0),
75	fDebug(0)
76
77	{
78
79	//
80	// Default Constructor.
81	//
82
83	if (fDebug > 0) {AliInfo("AliAnalysisTaskDiHadronPID Default Constructor.");}
84
85	for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
86	for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
87	fCorrelationsTOF[iPtClass][iSpecies] = 0x0;
88	fCorrelationsTOFTPC[iPtClass][iSpecies] = 0x0;
89	}
90	}
91
92
93	}
94
95	// -------------------------------------------------------------------------
96	AliAnalysisTaskDiHadronPID::AliAnalysisTaskDiHadronPID(const char* name):
97	AliAnalysisTaskSE(name),
98	fPIDResponse(0x0),
99	fEventCuts(0x0),
100	fTrackCutsTrigger(0x0),
101	fTrackCutsAssociated(0x0),
102	fPoolMgr(0x0),
103	fTriggerTracks(0x0),
104	fAssociatedTracks(0x0),
105	fCurrentAODEvent(0x0),
106	fOutputList(0x0),
107	fPtSpectrum(0x0),
108	fCorrelations(0x0),
109	fMixedEvents(0x0),
110	fTOFhistos(0x0),
111	fNDEtaBins(32),
112	fNDPhiBins(32),
113	fMinNEventsForMixing(5),
114	fPoolTrackDepth(2000),
115	fPoolSize(1000),
116	fCalculateTOFmismatch(kTRUE),
117	fT0Fill(0x0),
118	fLvsEta(0x0),
119	fLvsEtaProjections(0x0),
120	fDebug(0)
121
122	{
123
124	//
125	// Named Constructor.
126	//
127
128	if (fDebug > 0) {AliInfo("AliAnalysisTaskDiHadronPID Named Constructor.");}
129	if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
130
131	for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
132	for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
133	fCorrelationsTOF[iPtClass][iSpecies] = 0x0;
134	fCorrelationsTOFTPC[iPtClass][iSpecies] = 0x0;
135	}
136	}
137
138	DefineInput(0,TChain::Class());
139	DefineOutput(1,TList::Class());
140
141	}
142
143	// -------------------------------------------------------------------------
144	AliAnalysisTaskDiHadronPID::~AliAnalysisTaskDiHadronPID() {;
145
146	//
147	// Destructor.
148	//
149
150	if (fDebug > 0) {AliInfo("AliAnalysisTaskDiHadronPID Destructor.");}
151	if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
152
153	}
154
155	// -------------------------------------------------------------------------
156	void AliAnalysisTaskDiHadronPID::UserCreateOutputObjects() {
157
158	//
159	// Create Output objects.
160	//
161
162	if (fDebug > 0) {AliInfo("UserCreateOutputObjects()");}
163	if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
164
165	// --- BEGIN: Initialization on the worker nodes ---
166	AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
167	AliInputEventHandler* inputHandler = dynamic_cast<AliInputEventHandler*> (manager->GetInputEventHandler());
168
169	// Getting the pointer to the PID response object.
170	fPIDResponse = inputHandler->GetPIDResponse();
171
172	// Not very neat - only set up for 0-5% analysis.
173	Int_t nCentralityBins = 5;
174	Double_t centralityBins[] = {0.,1.,2.,3.,4.,5.};
175
176	Int_t nZvtxBins = 7;
177	Double_t vertexBins[] = {-7.,-5.,-3.,-1.,1.,3.,5.,7.};
178
179	fPoolMgr = new AliEventPoolManager(fPoolSize, fPoolTrackDepth, nCentralityBins, (Double_t) centralityBins, nZvtxBins, (Double_t) vertexBins);
180	// --- END ---
181
182	// Create the output list.
183	fOutputList = new TList();
184	fOutputList->SetOwner(kTRUE);
185
186	// Creating all requested histograms locally.
187	fEventCuts->CreateHistos();
188	fOutputList->Add(fEventCuts);
189
190	fTrackCutsTrigger->CreateHistos();
191	fOutputList->Add(fTrackCutsTrigger);
192
193	fTrackCutsAssociated->CreateHistos();
194	fOutputList->Add(fTrackCutsAssociated);
195
196	// Get the pT axis for the PID histograms.
197	Float_t* ptaxis = fTrackCutsAssociated->GetPtAxisPID();
198	Int_t nptbins = fTrackCutsAssociated->GetNPtBinsPID();
199
200	// Create Pt spectrum histogram.
201	fPtSpectrum = new TH1F("fPtSpectrum","p_{T} Spectrum;p_{T} (GeV/c);Count",nptbins,ptaxis);
202	fOutputList->Add(fPtSpectrum);
203
204	// Create unidentified correlations histogram.
205	fCorrelations = AliHistToolsDiHadronPID::MakeHist3D("fCorrelations","Correlations;#Delta#phi;#Delta#eta;p_{T} (GeV/c)",
206	fNDPhiBins,-TMath::Pi()/2.,3.*TMath::Pi()/2.,
207	fNDEtaBins,-1.6,1.6,
208	nptbins, ptaxis);
209	fOutputList->Add(fCorrelations);
210
211	// Create unidentified mixed events histogram.
212	fMixedEvents = AliHistToolsDiHadronPID::MakeHist3D("fMixedEvents","Mixed Events;#Delta#phi;#Delta#eta;p_{T} (GeV/c)",
213	fNDPhiBins,-TMath::Pi()/2.,3.*TMath::Pi()/2.,
214	fNDEtaBins,-1.6,1.6,
215	nptbins, ptaxis);
216	fOutputList->Add(fMixedEvents);
217
218	// Create TOF correlations histograms (DPhi,DEta,pt,TOF).
219	fTOFhistos = new TObjArray(15);
220	fTOFhistos->SetName("CorrelationsTOF");
221
222	Int_t nbins[4] = {fNDPhiBins,fNDEtaBins,0.,0.};
223	Double_t min[4] = {-TMath::Pi()/2.,-1.6,0.,0.};
224	Double_t max[4] = {3.*TMath::Pi()/2.,1.6,0.,0.};
225
226	for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
227
228	nbins[2] = fTrackCutsAssociated->GetNPtBinsPID(iPtClass);
229	min[2] = fTrackCutsAssociated->GetPtClassMin(iPtClass);
230	max[2] = fTrackCutsAssociated->GetPtClassMax(iPtClass);
231
232	for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
233
234	nbins[3] = fTrackCutsAssociated->GetNTOFbins(iPtClass,iSpecies);
235	min[3] = fTrackCutsAssociated->GetTOFmin(iPtClass,iSpecies);
236	max[3] = fTrackCutsAssociated->GetTOFmax(iPtClass,iSpecies);
237
238	fCorrelationsTOF[iPtClass][iSpecies] = new THnF(
239	Form("fCorrelationsTOF_%i_%i",iPtClass,iSpecies),
240	Form("CorrelationsTOF_%i_%i",iPtClass,iSpecies),
241	4,nbins,min,max);
242
243	fTOFhistos->Add(fCorrelationsTOF[iPtClass][iSpecies]);
244
245	}
246	}
247
248	fOutputList->Add(fTOFhistos);
249
250	// TODO: Create TOF/TPC correlations histogram.
251
252	// Load external TOF histograms if flag is set.
253	if (fCalculateTOFmismatch) {LoadExtMismatchHistos();}
254
255	PostData(1,fOutputList);
256
257	}
258
259	// -------------------------------------------------------------------------
260	void AliAnalysisTaskDiHadronPID::LocalInit() {
261
262	//
263	// Initialize on the client. (or on my computer?? - I think so...)
264	//
265
266	if (fDebug > 0) {AliInfo("LocalInit()");}
267	if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
268
269	}
270
271	// -------------------------------------------------------------------------
272	void AliAnalysisTaskDiHadronPID::UserExec(Option_t*) {
273
274	//
275	// Main Loop.
276	//
277
278	if (fDebug > 0) {AliInfo("UserExec()");}
279	if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
280
281	// Input Current Event.
282	fCurrentAODEvent = dynamic_cast<AliAODEvent*>(InputEvent());
283	if (!fCurrentAODEvent) AliFatal("No Event Found!");
284
285	if (!fEventCuts->IsSelected(fCurrentAODEvent)) {return;}
286
287	// Fill the global tracks array. - NOT NEEDED I THINK, since we're not using
288	// bit 1<<7 for the associated tracks!
289
290	// Let the track cut objects know that a new event will start.
291	fTrackCutsTrigger->StartNewEvent();
292	fTrackCutsAssociated->StartNewEvent();
293
294	// Create arrays for trigger/associated particles.
295	fTriggerTracks = new TObjArray(350);
296	fTriggerTracks->SetOwner(kTRUE);
297
298	fAssociatedTracks = new TObjArray(3500);
299	fAssociatedTracks->SetOwner(kTRUE);
300
301	for (Int_t iTrack = 0; iTrack < fCurrentAODEvent->GetNumberOfTracks(); iTrack++) {
302
303	AliAODTrack* track = (AliAODTrack*)fCurrentAODEvent->GetTrack(iTrack);
304	AliTrackDiHadronPID* pidtrack = new AliTrackDiHadronPID(track,0x0,0x0,fPIDResponse);
305	pidtrack->ForgetAboutPointers();
306	pidtrack->SetDebugLevel(fDebug);
307
308	Double_t rndhittime = -1.e21;
309	if (fCalculateTOFmismatch) rndhittime = GenerateRandomHit(pidtrack->Eta());
310
311	// Fill p_T spectrum.
312	fPtSpectrum->Fill(pidtrack->Pt());
313
314	// Fill the trigger/associated tracks array.
315	if (fTrackCutsTrigger->IsSelectedData(pidtrack,rndhittime)) {fTriggerTracks->AddLast(pidtrack);}
316	else if (fTrackCutsAssociated->IsSelectedData(pidtrack,rndhittime)) {fAssociatedTracks->AddLast(pidtrack);}
317	else {delete pidtrack;}
318
319	}
320
321	// Fill Correlation histograms.
322	for (Int_t iTrigger = 0; iTrigger < fTriggerTracks->GetEntriesFast(); iTrigger++) {
323	AliTrackDiHadronPID* triggertrack = (AliTrackDiHadronPID*)fTriggerTracks->At(iTrigger);
324
325	for (Int_t iAssociated = 0; iAssociated < fAssociatedTracks->GetEntriesFast(); iAssociated++) {
326	AliTrackDiHadronPID* associatedtrack = (AliTrackDiHadronPID*)fAssociatedTracks->At(iAssociated);
327
328	Double_t DPhi = triggertrack->Phi() - associatedtrack->Phi();
329	if (DPhi < -TMath::Pi()/2.) {DPhi += 2.*TMath::Pi();}
330	if (DPhi > 3.TMath::Pi()/2.) {DPhi -= 2.TMath::Pi();}
331
332	Double_t DEta = triggertrack->Eta() - associatedtrack->Eta();
333	fCorrelations->Fill(DPhi,DEta,associatedtrack->Pt());
334
335	Double_t tofhistfill[4] = {DPhi,DEta,associatedtrack->Pt(),-999.};
336
337	for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
338	tofhistfill[3] = associatedtrack->GetTOFsignalMinusExpected(iSpecies);
339
340	for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
341
342	// prevent under/over-flow bins to be filled.
343	Double_t ptmin = fTrackCutsAssociated->GetPtClassMin(iPtClass);
344	Double_t ptmax = fTrackCutsAssociated->GetPtClassMax(iPtClass);
345	Double_t apt = associatedtrack->Pt();
346
347	if ( (apt > ptmin) && (apt < ptmax) ) {
348	fCorrelationsTOF[iPtClass][iSpecies]->Fill(tofhistfill);
349	}
350
351	}
352	}
353	}
354	}
355
356	cout<<"Triggers: "<<fTriggerTracks->GetEntriesFast()<<" Associateds: "<<fAssociatedTracks->GetEntriesFast()<<endl;
357
358	// Determine centrality of current event.
359	TString centralityestimator = fEventCuts->GetCentralityEstimator();
360	AliCentrality* currentcentrality = fCurrentAODEvent->GetCentrality();
361	Float_t percentile = currentcentrality->GetCentralityPercentile(centralityestimator.Data());
362
363	// Determine vtxz of current event.
364	AliAODVertex* currentprimaryvertex = fCurrentAODEvent->GetPrimaryVertex();
365	Double_t vtxz = currentprimaryvertex->GetZ();
366
367	// Obtain event pool for current centrality/ vtxz.
368	AliEventPool* poolin = fPoolMgr->GetEventPool(percentile, vtxz);
369	if (!poolin) {AliFatal(Form("No pool found for centrality = %f, vtxz = %f", percentile, vtxz));}
370	// TObjArray* fGlobalTracksArray;
371
372	// Mix events if there are enough events in the pool. (TODO: should be a data member.)
373	if (poolin->GetCurrentNEvents() >= fMinNEventsForMixing) {
374	if (fDebug) {AliInfo("Mixing Events.");}
375
376	// Loop over all triggers in this event.
377	for (Int_t iTrigger = 0; iTrigger < fTriggerTracks->GetEntriesFast(); iTrigger++) {
378	AliTrackDiHadronPID* triggertrack = (AliTrackDiHadronPID*)fTriggerTracks->At(iTrigger);
379
380	// Loop over all events in the event pool.
381	for (Int_t iMixEvent = 0; iMixEvent < poolin->GetCurrentNEvents(); iMixEvent++) {
382	TObjArray* mixtracks = poolin->GetEvent(iMixEvent);
383
384	// Loop over all mixed tracks.
385	for (Int_t iMixTrack = 0; iMixTrack < mixtracks->GetEntriesFast(); iMixTrack++) {
386	AliTrackDiHadronPID* mixtrack = (AliTrackDiHadronPID*)mixtracks->At(iMixTrack);
387
388	Double_t DPhi = triggertrack->Phi() - mixtrack->Phi();
389	if (DPhi < -TMath::Pi()/2.) {DPhi += 2.*TMath::Pi();}
390	if (DPhi > 3.TMath::Pi()/2.) {DPhi -= 2.TMath::Pi();}
391
392	Double_t DEta = triggertrack->Eta() - mixtrack->Eta();
393	fMixedEvents->Fill(DPhi,DEta,mixtrack->Pt());
394
395	}
396	}
397	}
398	}
399
400	// Update the event pool.
401	AliEventPool* poolout = fPoolMgr->GetEventPool(percentile, vtxz); // Get the buffer associated with the current centrality and z-vtx
402	if (!poolout) AliFatal(Form("No pool found for centrality = %f, vtx_z = %f", percentile, vtxz));
403
404	// Q: is it a problem that the fAssociatedTracks array can be bigger than the number of tracks inside?
405	poolout->UpdatePool(fAssociatedTracks);
406
407	// Delete trigger array and its contents. Set pointer to zero.
408	// Don't delete the associated array, since ownership has been transferred to the pool manager. Set pointer to zero.
409	fTriggerTracks->Delete();
410	delete fTriggerTracks;
411	fTriggerTracks = 0x0;
412	fAssociatedTracks = 0x0;
413
414	// Tell the track cut object that the event is done.
415	fTrackCutsTrigger->EventIsDone(0);
416	fTrackCutsAssociated->EventIsDone(0);
417
418	PostData(1,fOutputList);
419
420	}
421
422	// -------------------------------------------------------------------------
423	Bool_t AliAnalysisTaskDiHadronPID::LoadExtMismatchHistos() {
424
425	//
426	// Attempting to load a root file containing information needed
427	// to generate random TOF hits.
428	//
429
430	if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
431
432	// Opening external TOF file.
433	if (fDebug > 0) cout<<"Trying to open TOFmismatchHistos.root ..."<<endl;
434	TFile* fin = 0x0;
435	fin = TFile::Open("alien:///alice/cern.ch/user/m/mveldhoe/rootfiles/TOFmismatchHistos.root");
436	if (!fin) {
437	AliWarning("Couln't open TOFmismatchHistos, will not calculate mismatches...");
438	fCalculateTOFmismatch = kFALSE;
439	return kFALSE;
440	}
441
442	// Check if the required histograms are present.
443	TH1F* tmp1 = (TH1F*)fin->Get("hNewT0Fill");
444	if (!tmp1) {
445	AliWarning("Couln't find hNewT0Fill, will not calculate mismatches...");
446	fCalculateTOFmismatch = kFALSE;
447	return kFALSE;
448	}
449	TH2F* tmp2 = (TH2F*)fin->Get("hLvsEta");
450	if (!tmp2) {
451	AliWarning("Couln't find hLvsEta, will not calculate mismatches...");
452	fCalculateTOFmismatch = kFALSE;
453	return kFALSE;
454	}
455
456	// Make a deep copy of the files in the histogram.
457	fT0Fill = (TH1F*)tmp1->Clone("fT0Fill");
458	fLvsEta = (TH2F*)tmp2->Clone("fLvsEta");
459
460	// Close the external file.
461	AliInfo("Closing external file.");
462	fin->Close();
463
464	// Creating a TObjArray for LvsEta projections.
465	const Int_t nbinseta = fLvsEta->GetNbinsX();
466	fLvsEtaProjections = new TObjArray(nbinseta);
467	fLvsEtaProjections->SetOwner(kTRUE);
468
469	// Making the projections needed (excluding underflow/ overflow).
470	for (Int_t iEtaBin = 1; iEtaBin < (nbinseta + 1); iEtaBin++) {
471	TH1F* tmp = (TH1F*)fLvsEta->ProjectionY(Form("LvsEtaProjection_%i",iEtaBin),iEtaBin,iEtaBin);
472	tmp->SetDirectory(0);
473	fLvsEtaProjections->AddAt(tmp,iEtaBin);
474	}
475
476	return kTRUE;
477
478	}
479
480	// -------------------------------------------------------------------------
481	Double_t AliAnalysisTaskDiHadronPID::GenerateRandomHit(Double_t eta) {
482
483	//
484	// Returns a random TOF time.
485	//
486
487	if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
488
489	// Default (error) value:
490	Double_t rndhittime = -1.e21;
491
492	// TOF mismatch flag is not turned on.
493	if (!fCalculateTOFmismatch) {
494	AliFatal("Called GenerateRandomHit() method, but flag fCalculateTOFmismatch not set.");
495	return rndhittime;
496	}
497
498	// TOF doesn't extend much further than 0.8.
499	if (TMath::Abs(eta) > 0.8) {
500	if (fDebug) {AliInfo("Tried to get a random hit for a track with eta > 0.8.");}
501	return rndhittime;
502	}
503
504	// Finding the bin of the eta.
505	TAxis* etaAxis = fLvsEta->GetXaxis();
506	Int_t etaBin = etaAxis->FindBin(eta);
507	//cout<<"Eta: "<<eta<<" bin: "<<etaBin<<endl;
508	const TH1F* lengthDistribution = (const TH1F*)fLvsEtaProjections->At(etaBin);
509
510	if (!lengthDistribution) {
511	AliFatal("length Distribution not found.");
512	return rndhittime;
513	}
514
515	Double_t currentRndLength = lengthDistribution->GetRandom(); // in cm.
516
517	// Similar to Roberto's code.
518	Double_t currentRndTime = currentRndLength / (TMath::C() * 1.e2 / 1.e12);
519	Double_t t0fill = -1.26416e+04;
520	rndhittime = fT0Fill->GetRandom() - t0fill + currentRndTime;
521
522	return rndhittime;
523
524	}
525
526	// -------------------------------------------------------------------------
527	void AliAnalysisTaskDiHadronPID::Terminate(Option_t*) {;
528
529	//
530	// Called when task is done.
531	//
532
533	if (fDebug > 0) {AliInfo("Terminate()");}
534	if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
535
536	delete fT0Fill;
537	fT0Fill = 0x0;
538	delete fLvsEta;
539	fLvsEta = 0x0;
540	delete fLvsEtaProjections;
541	fLvsEtaProjections = 0x0;
542
543	}