Correlations/DPhi/AliAnalysisTaskLeadingTrackUE.cxx
Correlations/DPhi/AliAnalysisTaskMinijet.cxx
Correlations/DPhi/AliAnalysisTaskDiHadron.cxx
- Correlations/DPhi/AliAnalysisTaskDiHadronPID.cxx
+ Correlations/DPhi/DiHadronPID/AliTrackDiHadronPID.cxx
+ Correlations/DPhi/DiHadronPID/AliAODTrackCutsDiHadronPID.cxx
+ Correlations/DPhi/DiHadronPID/AliAODEventCutsDiHadronPID.cxx
+ Correlations/DPhi/DiHadronPID/AliHistToolsDiHadronPID.cxx
+ Correlations/DPhi/DiHadronPID/AliAnalysisTaskDiHadronPID.cxx
Correlations/DPhi/FourierDecomposition/AliPool.cxx
Correlations/DPhi/FourierDecomposition/AliDhcTask.cxx
Correlations/DPhi/FourierDecomposition/AliMuonEffMC.cxx
+++ /dev/null
-// ----------------------------------------------------------------------------
-// This class makes di-hadron correlations, with TOF and TPC signals for
-// the associated particles.
-//
-// ----------------------------------------------------------------------------
-// Author: Misha Veldhoen (misha.veldhoen@cern.ch)
-// Last edit: May 2nd 2012. (v 8.00)
-// ----------------------------------------------------------------------------
-
-#include <iostream>
-#include "TChain.h"
-#include "TTree.h"
-#include "TH1F.h"
-#include "TH2F.h"
-#include "TH3F.h"
-#include "TCanvas.h"
-#include "TFile.h"
-#include "TClonesArray.h"
-
-#include "AliAODTrack.h"
-#include "AliAODEvent.h"
-#include "AliAODInputHandler.h"
-#include "AliAODVertex.h"
-//#include "AliAODPid.h"
-
-#include "AliAnalysisManager.h"
-#include "AliInputEventHandler.h"
-#include "AliPIDResponse.h"
-#include "AliTPCPIDResponse.h"
-//#include "AliTOFPIDResponse.h"
-
-// Includes for a MC run.
-#include "AliAODMCParticle.h"
-
-using namespace std;
-
-#include "AliAnalysisTaskDiHadronPID.h"
-
-ClassImp(AliAnalysisTaskDiHadronPID);
-
-//_____________________________________________________________________________
-AliAnalysisTaskDiHadronPID::AliAnalysisTaskDiHadronPID():
- AliAnalysisTaskSE(),
- fPIDResponse(0x0),
- fAODEvent(0x0),
- fAODHeader(0x0),
- fAODVertex(0x0),
- fAODTrack(0x0),
- fGlobalTracks(0x0),
- fMCTracks(0x0),
- fCentrality(0x0),
- fVertexZ(0x0),
- fDCA(0x0),
- fDCAZoomed(0x0),
- fDCAZoomedTwice(0x0),
- fDCACut(0x0),
- fDCAZoomedCut(0x0),
- fDCAZoomedTwiceCut(0x0),
- fITSHits(0x0),
- fTrackCutsCount(0x0),
- fTrackCutsPt(0x0),
- fTrackCutsEta(0x0),
- fTrackCutsPhi(0x0),
- fEtaSpectrumTrig(0x0),
- fEtaSpectrumAssoc(0x0),
- fPhiSpectrumAssoc(0x0),
- fTPCnSigmaProton(0x0),
- fTOFnSigmaProton(0x0),
- fTPCnSigmaPion(0x0),
- fTOFnSigmaPion(0x0),
- fTPCnSigmaKaon(0x0),
- fTOFnSigmaKaon(0x0),
- fTPCSignal(0x0),
- fTOFSignal(0x0),
- fDiHadron(0x0),
- fMixedEvents(0x0),
- fHistoList(0x0),
- fCalculateMixedEvents(kFALSE),
- fBeamType("PbPb"),
- fMC(kFALSE),
- fMaxEta(0.8),
- fMaxPlotEta(0.8),
- fMaxRap(0.5),
- fMaxPt(10.),
- fNEtaBins(32),
- fNPhiBins(36),
- fVertexZMixedEvents(2.),
- fCentralityCutMax(0.),
- fCentralityCutMin(10.),
- fZoomed(kFALSE),
- fDoITSCut(kFALSE),
- fDoDCACut(kFALSE),
- fDemandNoMismatch(kFALSE),
- fVerbose(0),
- fPrintBufferSize(kFALSE),
- fTrigBufferIndex(0),
- fTrigBufferSize(0),
- fTrigBufferMaxSize(100)
-
-{
-
- //
- // Default Constructor.
- //
-
- // Trigger buffer.
- for(Int_t ii=0; ii<25000; ii++) {
- for(Int_t jj=0; jj<4; jj++) {
- fTrigBuffer[ii][jj]=0;
- }
- }
-
- // The identified di-hadron correlations.
- for (Int_t ii = 0; ii < 3; ii++) {
- //fMixedEventsTPCTOFCut[ii]=0x0;
- for (Int_t jj = 0; jj < 10; jj++) {
- fDiHadronTPCTOF[ii][jj]=0x0;
- fMixedEventsTPCTOF[ii][jj]=0x0;
- fInclusiveTPCTOF[ii][jj]=0x0;
- fInclusiveTPCTOFRap[ii][jj]=0x0;
- }
- }
-
- // Track cut labels.
- for (Int_t ii=0; ii<8; ii++) {
- fTrackCutLabelNumbers[ii]=ii+1;
- }
-
- // Monte Carlo Efficiency Histograms.
- for (Int_t iSpecies=0; iSpecies<6; iSpecies++) {
- fPtEtaDistrDataPrim[iSpecies]=0x0;
- fPtEtaDistrDataSec[iSpecies]=0x0;
- fPtRapDistrDataPrimRapCut[iSpecies]=0x0;
- fPtRapDistrDataSecRapCut[iSpecies]=0x0;
-
- fPtEtaDistrMCPrim[iSpecies]=0x0;
- fPtEtaDistrMCSec[iSpecies]=0x0;
- fPtRapDistrMCPrimRapCut[iSpecies]=0x0;
- fPtRapDistrMCSecRapCut[iSpecies]=0x0;
-
- fDiHadronMC[iSpecies]=0x0;
-
- }
-
-}
-
-//_____________________________________________________________________________
-AliAnalysisTaskDiHadronPID::AliAnalysisTaskDiHadronPID(const char *name):
- AliAnalysisTaskSE(name),
- fPIDResponse(0x0),
- fAODEvent(0x0),
- fAODHeader(0x0),
- fAODVertex(0x0),
- fAODTrack(0x0),
- fGlobalTracks(0x0),
- fMCTracks(0x0),
- fCentrality(0x0),
- fVertexZ(0x0),
- fDCA(0x0),
- fDCAZoomed(0x0),
- fDCAZoomedTwice(0x0),
- fDCACut(0x0),
- fDCAZoomedCut(0x0),
- fDCAZoomedTwiceCut(0x0),
- fITSHits(0x0),
- fTrackCutsCount(0x0),
- fTrackCutsPt(0x0),
- fTrackCutsEta(0x0),
- fTrackCutsPhi(0x0),
- fEtaSpectrumTrig(0x0),
- fEtaSpectrumAssoc(0x0),
- fPhiSpectrumAssoc(0x0),
- fTPCnSigmaProton(0x0),
- fTOFnSigmaProton(0x0),
- fTPCnSigmaPion(0x0),
- fTOFnSigmaPion(0x0),
- fTPCnSigmaKaon(0x0),
- fTOFnSigmaKaon(0x0),
- fTPCSignal(0x0),
- fTOFSignal(0x0),
- fDiHadron(0x0),
- fMixedEvents(0x0),
- fHistoList(0x0),
- fCalculateMixedEvents(kFALSE),
- fBeamType("PbPb"),
- fMC(kFALSE),
- fMaxEta(0.8),
- fMaxPlotEta(0.8),
- fMaxRap(0.5),
- fMaxPt(10.),
- fNEtaBins(32),
- fNPhiBins(36),
- fVertexZMixedEvents(2.),
- fCentralityCutMax(0.),
- fCentralityCutMin(10.),
- fZoomed(kFALSE),
- fDoITSCut(kFALSE),
- fDoDCACut(kFALSE),
- fDemandNoMismatch(kFALSE),
- fVerbose(0),
- fPrintBufferSize(kFALSE),
- fTrigBufferIndex(0),
- fTrigBufferSize(0),
- fTrigBufferMaxSize(100)
-
-{
-
- //
- // Named Constructor.
- //
-
- DefineInput(0, TChain::Class());
- DefineOutput(1, TList::Class());
-
- // Trigger buffer.
- for(Int_t ii=0; ii<25000; ii++) {
- for(Int_t jj=0; jj<4; jj++) {
- fTrigBuffer[ii][jj]=0;
- }
- }
-
- // The identified di-hadron correlations.
- for (Int_t ii = 0; ii < 3; ii++) {
- //fMixedEventsTPCTOFCut[ii]=0x0;
- for (Int_t jj = 0; jj < 10; jj++) {
- fDiHadronTPCTOF[ii][jj]=0x0;
- fMixedEventsTPCTOF[ii][jj]=0x0;
- fInclusiveTPCTOF[ii][jj]=0x0;
- fInclusiveTPCTOFRap[ii][jj]=0x0;
- }
- }
-
- // Track cut labels.
- for (Int_t ii=0; ii<8; ii++) {
- fTrackCutLabelNumbers[ii]=ii+1;
- }
-
- // Monte Carlo Efficiency Histograms.
- for (Int_t iSpecies=0; iSpecies<6; iSpecies++) {
- fPtEtaDistrDataPrim[iSpecies]=0x0;
- fPtEtaDistrDataSec[iSpecies]=0x0;
- fPtRapDistrDataPrimRapCut[iSpecies]=0x0;
- fPtRapDistrDataSecRapCut[iSpecies]=0x0;
-
- fPtEtaDistrMCPrim[iSpecies]=0x0;
- fPtEtaDistrMCSec[iSpecies]=0x0;
- fPtRapDistrMCPrimRapCut[iSpecies]=0x0;
- fPtRapDistrMCSecRapCut[iSpecies]=0x0;
-
- fDiHadronMC[iSpecies]=0x0;
-
- }
-
-}
-
-//_____________________________________________________________________________
-AliAnalysisTaskDiHadronPID::~AliAnalysisTaskDiHadronPID() {
-
- //
- // Destructor.
- //
-
- if(fGlobalTracks) {
- delete fGlobalTracks;
- fGlobalTracks=0x0;
- }
- /*
- if (fMCTracks) {
- delete fMCTracks;
- fMCTracks=0x0;
- }
-*/
-}
-
-//_____________________________________________________________________________
-void AliAnalysisTaskDiHadronPID::UserCreateOutputObjects()
-
-{
- //
- // Output objects.
- //
-
- // Print the settings of the analysis task.
- cout<<endl;
- cout<<"-----------------------------------------------"<<endl;
- cout<<" AliAnalysisTaskDiHadronPID Settings:"<<endl;
- cout<<"-----------------------------------------------"<<endl;
- cout<<"Verbose Level: "<<fVerbose<<endl;
- cout<<"Mixed Events Calculated: "<<fCalculateMixedEvents<<endl;
- cout<<"Beam Type: "<<fBeamType<<endl;
- cout<<"Run over MC: "<<fMC<<endl;
- cout<<Form("Max eta: %3.1f",fMaxEta)<<endl;
- cout<<Form("Max eta plotted: %3.1f",fMaxPlotEta)<<endl;
- cout<<Form("Max rapidity in spectra: %3.1f",fMaxRap)<<endl;
- cout<<Form("Max p_T for the triggers: %3.1f GeV/c.",fMaxPt)<<endl;
- cout<<"Nbins in eta and delta eta: "<<fNEtaBins<<endl;
- cout<<"Nbins in phi and delta phi: "<<fNPhiBins<<endl;
- cout<<Form("Events mixed if vertices differ %3.1f cm.",fVertexZMixedEvents)<<endl;
- cout<<Form("Centrality between %3.1f and %3.1f percent.",fCentralityCutMax,fCentralityCutMin)<<endl;
- cout<<"Tracks are cut if less than 2 SPD hits: "<<fDoITSCut<<endl;
- cout<<"Tracks cut if DCA is too big: "<<fDoDCACut<<endl;
- cout<<"Tracks cut if there is a TPC-TOF mismatch: "<<fDemandNoMismatch<<endl;
- cout<<"Maximum number of triggers stored: "<<fTrigBufferMaxSize<<endl;
- cout<<"-----------------------------------------------"<<endl;
- cout<<endl;
-
- // Obtain a pointer to the analysis manager.
- AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
- if (!manager) {
- if (fVerbose>0) cout<<"AliAnalysisTaskDiHadronPID::UserCreateOutputObjects -> ERROR: Analysis manager not found."<<endl;
- return;
- }
- if (fVerbose>1) {
- cout<<"AliAnalysisTaskDiHadronPID::UserCreateOutputObjects -> Analysis manager found."<<endl;
- }
-
- // Obtain a pointer to the input handler.
- AliInputEventHandler* inputHandler = dynamic_cast<AliInputEventHandler*> (manager->GetInputEventHandler());
- if (!inputHandler) {
- if (fVerbose>0) cout<<"AliAnalysisTaskDiHadronPID::UserCreateOutputObjects -> ERROR: Input handler not found."<<endl;
- return;
- }
- if (fVerbose>1) {
- cout<<"AliAnalysisTaskDiHadronPID::UserCreateOutputObjects -> Input handler found."<<endl;
- }
-
- // Obtain a pointer to the PID response object.
- fPIDResponse = inputHandler->GetPIDResponse();
- if (!fPIDResponse) {
- if (fVerbose>0) cout<<"AliAnalysisTaskDiHadronPID::UserCreateOutputObjects -> ERROR: PID response object not found."<<endl;
- return;
- }
- if (fVerbose>1) {
- cout<<"AliAnalysisTaskDiHadronPID::UserCreateOutputObjects -> PID response object found."<<endl;
- }
-
- // Create the output of the task.
- fHistoList = new TList();
- fHistoList->SetOwner(kTRUE);
-
- // Ranges in dPhi, dEta, and the number of bins for di-hadron correlations.
- Int_t binsHisto[4] = {fNPhiBins,fNEtaBins};
- Double_t minHisto[4] = {-TMath::Pi()/2.,-2*fMaxPlotEta};
- Double_t maxHisto[4] = {3.*TMath::Pi()/2,2*fMaxPlotEta};
-
- // --- EVENT QA PLOTS ---
-
- fCentrality = new TH1F("fCentrality","Centrality;Centrality;N_{evt}",100,0,100);
- fHistoList->Add(fCentrality);
- fVertexZ = new TH1F("fVertexZ","Vertex Z position;z (cm);N_{evt}",30,-15,15);
- fHistoList->Add(fVertexZ);
-
- // --- TRACK QA PLOTS ---
-
- fDCA = new TH2F("fDCA","DCA positions TPC only cuts;xy (cm);z (cm)",100,-5,5,100,-5,5);
- fHistoList->Add(fDCA);
- fDCAZoomed = new TH2F("fDCAZoomed","DCA positions TPC only cuts;xy (cm);z (cm)",100,-.5,.5,100,-.5,.5);
- fHistoList->Add(fDCAZoomed);
- fDCAZoomedTwice = new TH2F("fDCAZoomedTwice","DCA positions TPC only cuts;xy (cm);z (cm)",100,-.05,.05,100,-.05,.05);
- fHistoList->Add(fDCAZoomedTwice);
- fDCACut = new TH2F("fDCACut","DCA positions after DCA cut;xy (cm);z (cm)",100,-5,5,100,-5,5);
- fHistoList->Add(fDCACut);
- fDCAZoomedCut = new TH2F("fDCAZoomedCut","DCA positions after DCA cut;xy (cm);z (cm)",100,-.5,.5,100,-.5,.5);
- fHistoList->Add(fDCAZoomedCut);
- fDCAZoomedTwiceCut = new TH2F("fDCAZoomedTwiceCut","DCA positions after DCA cut;xy (cm);z (cm)",100,-.05,.05,100,-.05,.05);
- fHistoList->Add(fDCAZoomedTwiceCut);
-
- fITSHits = new TH1F("fITSHits","ITS hits",3,0,3);
- (fITSHits->GetXaxis())->SetBinLabel(1,"No SPD hit");
- (fITSHits->GetXaxis())->SetBinLabel(2,"One SPD hit");
- (fITSHits->GetXaxis())->SetBinLabel(3,"Two SPD hits");
- fHistoList->Add(fITSHits);
-
- // Determine the bins used in the track cut histograms.
- const Int_t ncuts = 5 + fDoITSCut + fDoDCACut + fDemandNoMismatch;
-
- if (!fDoITSCut) {
- fTrackCutLabelNumbers[3]--;
- fTrackCutLabelNumbers[4]--;
- fTrackCutLabelNumbers[5]--;
- fTrackCutLabelNumbers[6]--;
- fTrackCutLabelNumbers[7]--;
- }
-
- if (!fDoDCACut) {
- fTrackCutLabelNumbers[4]--;
- fTrackCutLabelNumbers[5]--;
- fTrackCutLabelNumbers[6]--;
- fTrackCutLabelNumbers[7]--;
- }
-
- fTrackCutsCount = new TH1F("fTrackCutsCount","Track Cuts;;Count",ncuts,0,ncuts);
- (fTrackCutsCount->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[0],"Std. TPC Only");
- (fTrackCutsCount->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[1],Form("#eta < %3.1f",fMaxEta));
- if (fDoITSCut) (fTrackCutsCount->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[2],"1-2 SPD hits");
- if (fDoDCACut) (fTrackCutsCount->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[3],"DCA cut");
- (fTrackCutsCount->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[4],"p_{T} < 5.0 GeV/c");
- (fTrackCutsCount->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[5],"TPCpid (p_{T} < 5.0 GeV/c)");
- (fTrackCutsCount->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[6],"TOFpid (p_{T} < 5.0 GeV/c)");
- if (fDemandNoMismatch) (fTrackCutsCount->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[7],"No Mismatch (p_{T} < 5.0 GeV/c)");
- fHistoList->Add(fTrackCutsCount);
-
- fTrackCutsPt = new TH2F("fTrackCutsPt","Track Cuts vs p_{T};;p_{T};Count",ncuts,0,ncuts,40,0,fMaxPt);
- (fTrackCutsPt->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[0],"Std. TPC Only");
- (fTrackCutsPt->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[1],Form("#eta < %3.1f",fMaxEta));
- if (fDoITSCut) (fTrackCutsPt->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[2],"1-2 SPD hits");
- if (fDoDCACut) (fTrackCutsPt->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[3],"DCA cut");
- (fTrackCutsPt->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[4],"p_{T} < 5.0 GeV/c");
- (fTrackCutsPt->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[5],"TPCpid (p_{T} < 5.0 GeV/c)");
- (fTrackCutsPt->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[6],"TOFpid (p_{T} < 5.0 GeV/c)");
- if (fDemandNoMismatch) (fTrackCutsPt->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[7],"No Mismatch (p_{T} < 5.0 GeV/c)");
- fHistoList->Add(fTrackCutsPt);
-
- fTrackCutsEta = new TH2F("fTrackCutsEta","Track Cuts vs #eta;;#eta;Count",ncuts,0,ncuts,4*binsHisto[1],-fMaxPlotEta,fMaxPlotEta);
- (fTrackCutsEta->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[0],"Std. TPC Only");
- (fTrackCutsEta->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[1],Form("#eta < %3.1f",fMaxEta));
- if (fDoITSCut) (fTrackCutsEta->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[2],"1-2 SPD hits");
- if (fDoDCACut) (fTrackCutsEta->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[3],"DCA cut");
- (fTrackCutsEta->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[4],"p_{T} < 5.0 GeV/c");
- (fTrackCutsEta->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[5],"TPCpid (p_{T} < 5.0 GeV/c)");
- (fTrackCutsEta->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[6],"TOFpid (p_{T} < 5.0 GeV/c)");
- if (fDemandNoMismatch) (fTrackCutsEta->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[7],"No Mismatch (p_{T} < 5.0 GeV/c)");
- fHistoList->Add(fTrackCutsEta);
-
- fTrackCutsPhi = new TH2F("fTrackCutsPhi","Track Cuts vs #phi;;#phi;Count",ncuts,0,ncuts,72,0,2.*TMath::Pi());
- (fTrackCutsPhi->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[0],"Std. TPC Only");
- (fTrackCutsPhi->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[1],Form("#eta < %3.1f",fMaxEta));
- if (fDoITSCut) (fTrackCutsPhi->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[2],"1-2 SPD hits");
- if (fDoDCACut) (fTrackCutsPhi->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[3],"DCA cut");
- (fTrackCutsPhi->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[4],"p_{T} < 5.0 GeV/c");
- (fTrackCutsPhi->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[5],"TPCpid (p_{T} < 5.0 GeV/c)");
- (fTrackCutsPhi->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[6],"TOFpid (p_{T} < 5.0 GeV/c)");
- if (fDemandNoMismatch) (fTrackCutsPhi->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[7],"No Mismatch (p_{T} < 5.0 GeV/c)");
- fHistoList->Add(fTrackCutsPhi);
-
- fEtaSpectrumTrig = new TH1F("fEtaSpectrumTrig","#eta Spectrum Triggers;#eta;Count",4*binsHisto[1],-fMaxPlotEta,fMaxPlotEta);
- fHistoList->Add(fEtaSpectrumTrig);
- fEtaSpectrumAssoc = new TH2F("fEtaSpectrumAssoc","#eta Spectrum Associateds;#eta;p_{T};Count",4*binsHisto[1],-fMaxPlotEta,fMaxPlotEta,10,0.,5.);
- fHistoList->Add(fEtaSpectrumAssoc);
- fPhiSpectrumAssoc = new TH2F("fPhiSpectrumAssoc","#phi Spectrum Associateds;#phi;p_{T};Count",72,0,2.*TMath::Pi(),10,0.,5.);
- fHistoList->Add(fPhiSpectrumAssoc);
-
- // --- QA PLOTS PID ---
-
- fTPCnSigmaProton = new TH2F("fTPCnSigmaProton","n#sigma plot for the TPC (Protons);p (GeV/c);n#sigma",100,0.,5.,100,-10.,10.);
- fHistoList->Add(fTPCnSigmaProton);
- fTOFnSigmaProton = new TH2F("fTOFnSigmaProton","n#sigma plot for the TOF (Protons);p (GeV/c);n#sigma",100,0.,5.,100,-10.,10.);
- fHistoList->Add(fTOFnSigmaProton);
- fTPCnSigmaPion = new TH2F("fTPCnSigmaPion","n#sigma plot for the TPC (Pions);p (GeV/c);n#sigma",100,0.,5.,100,-10.,10.);
- fHistoList->Add(fTPCnSigmaPion);
- fTOFnSigmaPion = new TH2F("fTOFnSigmaPion","n#sigma plot for the TOF (Pions);p (GeV/c);n#sigma",100,0.,5.,100,-10.,10.);
- fHistoList->Add(fTOFnSigmaPion);
- fTPCnSigmaKaon = new TH2F("fTPCnSigmaKaon","n#sigma plot for the TPC (Kaons);p (GeV/c);n#sigma",100,0.,5.,100,-10.,10.);
- fHistoList->Add(fTPCnSigmaKaon);
- fTOFnSigmaKaon = new TH2F("fTOFnSigmaKaon","n#sigma plot for the TOF (Kaons);p (GeV/c);n#sigma",100,0.,5.,100,-10.,10.);
- fHistoList->Add(fTOFnSigmaKaon);
- fTPCSignal = new TH3F("fTPCSignal","TPC Signal;#eta;p_{T} GeV/c;dE/dx",25,-fMaxEta,fMaxEta,100,0.,5.,150,0.,300.);
- fHistoList->Add(fTPCSignal);
- fTOFSignal = new TH3F("fTOFSignal","TOF Signal;#eta;p_{T} GeV/c;t",25,-fMaxEta,fMaxEta,100,0.,5.,150,10000.,25000.);
- fHistoList->Add(fTOFSignal);
-
- // --- DI-HADRON CORRELATIONS AND MIXED EVENTS WITH TPC AND TOF SIGNALS ---
-
- // Unzoomed pictures.
- Int_t binsTPCnoZoom[3][10] = {{100,100,100,100,100,100,100,100,100,100},
- {100,100,100,100,100,100,100,100,100,100},
- {100,100,100,100,100,100,100,100,100,100}};
-
- Double_t minTPCnoZoom[3][10] = {{ -50., -50.,-30.,-30.,-30.,-30.,-30.,-30.,-30.,-30.},
- { -100., -50.,-20.,-20.,-20.,-20.,-20.,-20.,-20.,-20.},
- { -200.,-150.,-50.,-30.,-20.,-20.,-20.,-20.,-20.,-20.}};
- Double_t maxTPCnoZoom[3][10] = {{ 150., 100., 50., 30., 25., 25., 25., 25., 25., 25.},
- { 50., 100., 40., 30., 30., 30., 30., 30., 30., 30.},
- { 100., 50., 50., 30., 30., 30., 30., 30., 30., 30.}};
-
- Int_t binsTOFnoZoom[3][10] = {{100,100,100,100,100,100,100,100,100,100},
- {100,100,100,100,100,100,100,100,100,100},
- {100,100,100,100,100,100,100,100,100,100}};
-
- Double_t minTOFnoZoom[3][10] = {{-2000.,-2000.,-1000., -500., -500., -500., -500., -500., -500., -500.},
- {-2000.,-2000.,-2000.,-1000.,-1000.,-1000.,-1000.,-1000.,-1000.,-1000.},
- {-2000.,-2000.,-6000.,-3000.,-2000.,-1500.,-1500.,-1500.,-1500.,-1500.}};
- Double_t maxTOFnoZoom[3][10] = {{ 2000., 2000., 5000., 3000., 2000., 1500., 1500., 1500., 1500., 1500.},
- { 2000., 2000., 5000., 2000., 1500., 1500., 1500., 1500., 1500., 1500.},
- { 2000., 2000., 2000., 1000., 1000., 1000., 1000., 1000., 1000., 1000.}};
-
- // Zoomed pictures.
- Int_t binsTPCZoom[3][10] = {{100,100,100,100,100,100,100,100,100,100},
- {100,100,100,100,100,100,100,100,100,100},
- {100,100,100,100,100,100,100,100,100,100}};
-
- Double_t minTPCZoom[3][10] = {{ -20., -20.,-20.,-20.,-20.,-20.,-20.,-20.,-20.,-20.},
- { -20., -20.,-20.,-20.,-20.,-20.,-20.,-20.,-20.,-20.},
- { -40., -20.,-20.,-20.,-20.,-20.,-20.,-20.,-20.,-20.}};
- Double_t maxTPCZoom[3][10] = {{ 20., 20., 20., 20., 20., 20., 20., 20., 20., 20.},
- { 20., 20., 20., 20., 20., 20., 20., 20., 20., 20.},
- { 40., 20., 20., 20., 20., 30., 30., 30., 30., 30.}};
-
- Int_t binsTOFZoom[3][10] = {{100,100,100,100,100,100,100,100,100,100},
- {100,100,100,100,100,100,100,100,100,100},
- {100,100,100,100,100,100,100,100,100,100}};
-
- Double_t minTOFZoom[3][10] = {{-1000.,-1000., -500., -500., -500., -400., -400., -400., -400., -400.},
- { -800., -800., -800., -800., -800., -600., -500., -500., -400., -400.},
- {-1000.,-1000.,-1000.,-1000., -800.,-1000.,-1000., -800., -700., -700.}};
- Double_t maxTOFZoom[3][10] = {{ 1000., 1000., 1000., 1000., 1000., 1000., 1000., 900., 800., 700.},
- { 1000., 1000., 500., 500., 500., 900., 700., 700., 600., 500.},
- { 1000., 1000., 1000., 500., 500., 500., 500., 500., 500., 500.}};
-
- const char* species[3] = {"Pion","Kaon","Proton"};
-
- for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
-
- for (Int_t iPtBin = 0; iPtBin < 10; iPtBin++) {
-
- if (fZoomed) {
- binsHisto[2] = binsTPCZoom[iSpecies][iPtBin];
- binsHisto[3] = binsTOFZoom[iSpecies][iPtBin];
- minHisto[2] = minTPCZoom[iSpecies][iPtBin];
- minHisto[3] = minTOFZoom[iSpecies][iPtBin];
- maxHisto[2] = maxTPCZoom[iSpecies][iPtBin];
- maxHisto[3] = maxTOFZoom[iSpecies][iPtBin];
- } else {
- binsHisto[2] = binsTPCnoZoom[iSpecies][iPtBin];
- binsHisto[3] = binsTOFnoZoom[iSpecies][iPtBin];
- minHisto[2] = minTPCnoZoom[iSpecies][iPtBin];
- minHisto[3] = minTOFnoZoom[iSpecies][iPtBin];
- maxHisto[2] = maxTPCnoZoom[iSpecies][iPtBin];
- maxHisto[3] = maxTOFnoZoom[iSpecies][iPtBin];
- }
-
- // Di-Hadron Correlations
- fDiHadronTPCTOF[iSpecies][iPtBin] = new THnSparseF(Form("fDiHadronTPCTOF_%s_%i",species[iSpecies],iPtBin),
- Form("Di-Hadron Correlation (%s) %3.1f < p_{T} < %3.1f;#Delta#phi;#Delta#eta;dE/dx;t (ms);",species[iSpecies],iPtBin*0.5,(iPtBin+1)*0.5),
- 4,binsHisto,minHisto,maxHisto);
- fHistoList->Add(fDiHadronTPCTOF[iSpecies][iPtBin]);
-
- // Mixed Events.
- fMixedEventsTPCTOF[iSpecies][iPtBin] = new THnSparseF(Form("fMixedEventsTPCTOF_%s_%i",species[iSpecies],iPtBin),
- Form("Mixed Events (%s) %3.1f < p_{T} < %3.1f;#Delta#phi;#Delta#eta;dE/dx;t (ms);",species[iSpecies],iPtBin*0.5,(iPtBin+1)*0.5),
- 4,binsHisto,minHisto,maxHisto);
- //fHistoList->Add(fMixedEventsTPCTOF[iSpecies][iPtBin]);
-
- // Inclusive spectra.
- fInclusiveTPCTOF[iSpecies][iPtBin] = new TH3F(Form("fInclusiveTPCTOF_%s_%i",species[iSpecies],iPtBin),
- Form("Inclusive Spectrum (%s) %3.1f < p_{T} < %3.1f;dE/dx;t (ms);#eta",species[iSpecies],iPtBin*0.5,(iPtBin+1)*0.5),
- binsHisto[2],minHisto[2],maxHisto[2],binsHisto[3],minHisto[3],maxHisto[3],fNEtaBins,-fMaxEta,fMaxEta);
- fHistoList->Add(fInclusiveTPCTOF[iSpecies][iPtBin]);
-
- // Inclusive spectra with additional rapidity cut.
- fInclusiveTPCTOFRap[iSpecies][iPtBin] = new TH3F(Form("fInclusiveTPCTOFRap_%s_%i",species[iSpecies],iPtBin),
- Form("Inclusive Spectrum (%s) %3.1f < p_{T} < %3.1f |Y| < %3.1f;dE/dx;t (ms);Y",species[iSpecies],iPtBin*0.5,(iPtBin+1)*0.5,fMaxRap),
- binsHisto[2],minHisto[2],maxHisto[2],binsHisto[3],minHisto[3],maxHisto[3],20,-fMaxRap,fMaxRap);
- fHistoList->Add(fInclusiveTPCTOFRap[iSpecies][iPtBin]);
-
-
- }
- }
-
- // Correlations without PID signals
- fDiHadron = new TH3F("fDiHadron","Di-Hadron Correlations;#Delta#phi;#Delta#eta;p_{T_assoc}",binsHisto[0],minHisto[0],maxHisto[0],binsHisto[1],minHisto[1],maxHisto[1],10,0.,5.);
- fDiHadron->Sumw2();
- fHistoList->Add(fDiHadron);
-
- fMixedEvents = new TH3F("fMixedEvents","Mixed Events;#Delta#phi;#Delta#eta;p_{T_assoc}",binsHisto[0],minHisto[0],maxHisto[0],binsHisto[1],minHisto[1],maxHisto[1],10,0.,5.);
- fMixedEvents->Sumw2();
- fHistoList->Add(fMixedEvents);
-
- const char* MCSpeciesName[6] = {"piplus","pimin","Kplus","Kmin","p","pbar"};
- const char* MCSpeciesTitle[6] = {"#pi^{+}","#pi^{-}","K^{+}","K^{-}","p","#bar{p}"};
- //const char* Cuts[3] = {"nocut","DCAcut","PIDandDCAcut"};
-
- // Efficiency Plots (Monte Carlo)
- if (fMC) {
-
- for (Int_t iSpecies=0; iSpecies<6; iSpecies++) {
-
- // Without rapidity cut.
- fPtEtaDistrDataPrim[iSpecies]=new TH2F(Form("fPtEtaDistrDataPrim_%s",MCSpeciesName[iSpecies]),
- Form("p_{T}-#eta distribution Data Primaries (%s);p_{T};#eta;Count",MCSpeciesTitle[iSpecies]),
- 20,0.0,5.0,fNEtaBins,-fMaxEta,fMaxEta);
- fPtEtaDistrDataPrim[iSpecies]->Sumw2();
- fHistoList->Add(fPtEtaDistrDataPrim[iSpecies]);
-
- fPtEtaDistrDataSec[iSpecies]=new TH2F(Form("fPtEtaDistrDataSec_%s",MCSpeciesName[iSpecies]),
- Form("p_{T}-#eta distribution Data Secondaries (%s);p_{T};#eta;Count",MCSpeciesTitle[iSpecies]),
- 20,0.0,5.0,fNEtaBins,-fMaxEta,fMaxEta);
- fPtEtaDistrDataSec[iSpecies]->Sumw2();
- fHistoList->Add(fPtEtaDistrDataSec[iSpecies]);
-
- fPtEtaDistrMCPrim[iSpecies]=new TH2F(Form("fPtEtaDistrMCPrim_%s",MCSpeciesName[iSpecies]),
- Form("p_{T}-#eta distribution MC Primaries (%s);p_{T};#eta;Count",MCSpeciesTitle[iSpecies]),
- 20,0.0,5.0,fNEtaBins,-fMaxEta,fMaxEta);
- fPtEtaDistrMCPrim[iSpecies]->Sumw2();
- fHistoList->Add(fPtEtaDistrMCPrim[iSpecies]);
-
- fPtEtaDistrMCSec[iSpecies]=new TH2F(Form("fPtEtaDistrMCSec_%s",MCSpeciesName[iSpecies]),
- Form("p_{T}-#eta distribution MC Secondaries (%s);p_{T};#eta;Count",MCSpeciesTitle[iSpecies]),
- 20,0.0,5.0,fNEtaBins,-fMaxEta,fMaxEta);
- fPtEtaDistrMCSec[iSpecies]->Sumw2();
- fHistoList->Add(fPtEtaDistrMCSec[iSpecies]);
-
- //fDiHadronMC[iSpecies]=new TH3F(Form("fDiHadronMC_%s",MCSpeciesName[iSpecies]),
- // Form("Di-Hadron Correlations MC (%s);#Delta#phi;#Delta#eta;p_{T_assoc}",MCSpeciesTitle[iSpecies]),
- // binsHisto[0],minHisto[0],maxHisto[0],binsHisto[1],minHisto[1],maxHisto[1],10,0.,5.);
- //fHistoList->Add(fDiHadronMC[iSpecies]);
-
- // With rapidity cut.
- fPtRapDistrDataPrimRapCut[iSpecies]=new TH2F(Form("fPtRapDistrDataPrimRapCut_%s",MCSpeciesName[iSpecies]),
- Form("p_{T}-Y distribution Data Primaries (%s) |Y| < %3.1f;p_{T};Y;Count",MCSpeciesTitle[iSpecies],fMaxRap),
- 20,0.0,5.0,20,-fMaxRap,fMaxRap);
- fPtRapDistrDataPrimRapCut[iSpecies]->Sumw2();
- fHistoList->Add(fPtRapDistrDataPrimRapCut[iSpecies]);
-
- fPtRapDistrDataSecRapCut[iSpecies]=new TH2F(Form("fPtRapDistrDataSecRapCut_%s",MCSpeciesName[iSpecies]),
- Form("p_{T}-Y distribution Data Secondaries (%s) |Y| < %3.1f;p_{T};Y;Count",MCSpeciesTitle[iSpecies],fMaxRap),
- 20,0.0,5.0,20,-fMaxRap,fMaxRap);
- fPtRapDistrDataSecRapCut[iSpecies]->Sumw2();
- fHistoList->Add(fPtRapDistrDataSecRapCut[iSpecies]);
-
- fPtRapDistrMCPrimRapCut[iSpecies]=new TH2F(Form("fPtRapDistrMCPrimRapCut_%s",MCSpeciesName[iSpecies]),
- Form("p_{T}-Y distribution MC Primaries (%s) |Y| < %3.1f;p_{T};Y;Count",MCSpeciesTitle[iSpecies],fMaxRap),
- 20,0.0,5.0,20,-fMaxRap,fMaxRap);
- fPtRapDistrMCPrimRapCut[iSpecies]->Sumw2();
- fHistoList->Add(fPtRapDistrMCPrimRapCut[iSpecies]);
-
- fPtRapDistrMCSecRapCut[iSpecies]=new TH2F(Form("fPtRapDistrMCSecRapCut_%s",MCSpeciesName[iSpecies]),
- Form("p_{T}-Y distribution MC Secondaries (%s) |Y| < %3.1f;p_{T};Y;Count",MCSpeciesTitle[iSpecies],fMaxRap),
- 20,0.0,5.0,20,-fMaxRap,fMaxRap);
- fPtRapDistrMCSecRapCut[iSpecies]->Sumw2();
- fHistoList->Add(fPtRapDistrMCSecRapCut[iSpecies]);
-
- }
- }
-
- PostData(1, fHistoList);
-
-}
-
-//_____________________________________________________________________________
-Int_t AliAnalysisTaskDiHadronPID::ClassifyTrack(AliAODTrack *track)
-
-{
- //
- // This function both classifies tracks, and fills the track QA histograms.
- //
- // Classifies the track in:
- // 0 -> Not Useful,
- // 1 -> Associated,
- // 2 -> Trigger.
- //
- // IDEA: later we can do this with filterbits.
- //
- // The following track cuts are applied for triggers:
- //
- // 1a) pT > 5.0 GeV/c, pT < fPtMax,
- // 2) StandardTPCOnlyTrackCuts (filterbit 7 for AOD's),
- // 3) eta < fMaxEta,
- // 4) ITS track cut, a track is only selected if it has at least one hit in the SPD,
- // that is, in one of the first two layers of the ITS. (can be switched on and off
- // using (Bool_t)fDoITSCut),
- // 5) DCA track cut, of all tracks with at least one SPD hit, the DCA is constrained as
- // a function of pT, in order to remove lots of the secondaries. (can be switched on and
- // off by using (Bool_t)fDoDCACut),
- //
- // For associateds all the same track cuts are applied, except 1a is
- // replaced by 2b. Moreover the following cuts are applied too:
- //
- // 1b) pT < 5.0 GeV/c,
- // 6) TPCpid,
- // 7) TOFpid,
- // 8) no TPCTOF mismatch.
- //
-
- // Check if a track is supplied.
- if (!track) {
- if (fVerbose>0) cout<<"AliAnalysisTaskDiHadronPID::SelectTrack -> ERROR: No track found."<<endl;
- return 0;
- }
-
- // Get basic track information.
- Int_t classification=0;
- Double_t pt = track->Pt();
- Double_t eta = track->Eta();
- Double_t phi = track->Phi();
-
- // 1) pT cut: First separation between triggers and associateds.
- if (pt<5.0) classification=1;
- if ((pt>5.0)&&(pt<fMaxPt)) classification = 2;
- if (!classification) return 0;
-
- // 2) StandardTPCOnlyTrackCuts.
- if (!(track->TestFilterMask(1<<7))) {
- //if (fVerbose>3) cout<<"Track Ignored: Did Not pass filterbit."<<endl;
- return 0;
- }
-
- if (fVerbose>3) {
- cout<<endl;
- cout<<"pt: "<<pt<<" eta: "<<eta<<" phi: "<<phi<<endl;
- }
-
- fTrackCutsCount->Fill(fTrackCutLabelNumbers[0]-.5);
- fTrackCutsPt->Fill(fTrackCutLabelNumbers[0]-.5,pt);
- fTrackCutsEta->Fill(fTrackCutLabelNumbers[0]-.5,eta);
- fTrackCutsPhi->Fill(fTrackCutLabelNumbers[0]-.5,phi);
-
- // 3) eta cut.
- if (TMath::Abs(eta)>fMaxEta) {
- if (fVerbose>3) cout<<"Track Ignored: Eta too large."<<endl;
- return 0;
- }
-
- fTrackCutsCount->Fill(fTrackCutLabelNumbers[1]-.5);
- fTrackCutsPt->Fill(fTrackCutLabelNumbers[1]-.5,pt);
- fTrackCutsEta->Fill(fTrackCutLabelNumbers[1]-.5,eta);
- fTrackCutsPhi->Fill(fTrackCutLabelNumbers[1]-.5,phi);
-
- // Obtaining ITS information.
- AliAODTrack* globaltrack = GetGlobalTrack(track);
- if (!globaltrack) {
- return 0;
- }
-
- Bool_t ITSLayerHit[6];
- for (Int_t iITSLayer=0; iITSLayer<6; iITSLayer++) {
- ITSLayerHit[iITSLayer] = globaltrack->HasPointOnITSLayer(iITSLayer);
- }
- Int_t SPDHits=ITSLayerHit[0]+ITSLayerHit[1];
-
- if (fVerbose>3) cout<<"SPD hits: "<<SPDHits<<endl;
-
- // Fill the ITS hist.
- fITSHits->Fill(SPDHits+0.5);
-
- // 4) ITS cut.
- if (fDoITSCut) {
- if (!SPDHits) {
- if (fVerbose>3) cout<<"Track Ignored: Not enough SPD hits."<<endl;
- return 0;
- }
- fTrackCutsCount->Fill(fTrackCutLabelNumbers[2]-.5);
- fTrackCutsPt->Fill(fTrackCutLabelNumbers[2]-.5,pt);
- fTrackCutsEta->Fill(fTrackCutLabelNumbers[2]-.5,eta);
- fTrackCutsPhi->Fill(fTrackCutLabelNumbers[2]-.5,phi);
- }
-
- // Propagate the global track to the DCA.
- Double_t PosAtDCA[2] = {-999,-999};
- Double_t covar[3] = {-999,-999,-999};
- globaltrack->PropagateToDCA(fAODVertex,fAODEvent->GetMagneticField(),100.,PosAtDCA,covar);
-
- // Fill the DCA hist (before cut)
- fDCA->Fill(PosAtDCA[0],PosAtDCA[1]);
- fDCAZoomed->Fill(PosAtDCA[0],PosAtDCA[1]);
- fDCAZoomedTwice->Fill(PosAtDCA[0],PosAtDCA[1]);
-
- // 5) DCA cut (See R_AA paper).
- Double_t DCAcutvalue[2];
- DCAcutvalue[0] = 0.018 + 0.035*TMath::Power(pt,-1.01);
- DCAcutvalue[1] = 2.;
- if (fDoDCACut) {
- if (SPDHits&&((TMath::Abs(PosAtDCA[0])>DCAcutvalue[0])||(TMath::Abs(PosAtDCA[1])>DCAcutvalue[1]))) {
- if (fVerbose>3) cout<<"Track Ignored: Enough SPD hits, but out of DCA range."<<endl;
- return 0;
- }
- fTrackCutsCount->Fill(fTrackCutLabelNumbers[3]-.5);
- fTrackCutsPt->Fill(fTrackCutLabelNumbers[3]-.5,pt);
- fTrackCutsEta->Fill(fTrackCutLabelNumbers[3]-.5,eta);
- fTrackCutsPhi->Fill(fTrackCutLabelNumbers[3]-.5,phi);
-
- // Fill the DCA hist (after cut)
- fDCACut->Fill(PosAtDCA[0],PosAtDCA[1]);
- fDCAZoomedCut->Fill(PosAtDCA[0],PosAtDCA[1]);
- fDCAZoomedTwiceCut->Fill(PosAtDCA[0],PosAtDCA[1]);
- }
-
- // Now all the common cuts have been performed. Tracks identified as a trigger will
- // be returned. Note that they will not appear in the histograms after track cut 5.
- if (classification==2) return 2;
-
- // Now we're left with only tracks with pT < 5.0 GeV/c.
- fTrackCutsCount->Fill(fTrackCutLabelNumbers[4]-.5);
- fTrackCutsPt->Fill(fTrackCutLabelNumbers[4]-.5,pt);
- fTrackCutsEta->Fill(fTrackCutLabelNumbers[4]-.5,eta);
- fTrackCutsPhi->Fill(fTrackCutLabelNumbers[4]-.5,phi);
-
- // Obtain the status of the track.
- ULong_t status = globaltrack->GetStatus();
-
- // 6) TPCpid
- if (!((status&AliAODTrack::kTPCpid)==AliAODTrack::kTPCpid)) {
- if (fVerbose>3) cout<<"Track Ignored: No TPC pid."<<endl;
- return 0;
- }
-
- fTrackCutsCount->Fill(fTrackCutLabelNumbers[5]-.5);
- fTrackCutsPt->Fill(fTrackCutLabelNumbers[5]-.5,pt);
- fTrackCutsEta->Fill(fTrackCutLabelNumbers[5]-.5,eta);
- fTrackCutsPhi->Fill(fTrackCutLabelNumbers[5]-.5,phi);
-
- // 7) TOFpid
- if (!((status&AliAODTrack::kTOFpid)==AliAODTrack::kTOFpid)) {
- if (fVerbose>3) cout<<"Track Ignored: No TOF pid."<<endl;
- return 0;
- }
-
- fTrackCutsCount->Fill(fTrackCutLabelNumbers[6]-.5);
- fTrackCutsPt->Fill(fTrackCutLabelNumbers[6]-.5,pt);
- fTrackCutsEta->Fill(fTrackCutLabelNumbers[6]-.5,eta);
- fTrackCutsPhi->Fill(fTrackCutLabelNumbers[6]-.5,phi);
-
- // 8) TPC, TOF mismatch.
- if (fDemandNoMismatch) {
- if ((status&AliAODTrack::kTOFmismatch)==AliAODTrack::kTOFmismatch) {
- if (fVerbose>3) cout<<"Track Ignored: TOF mismatch found."<<endl;
- return 0;
- }
- fTrackCutsCount->Fill(fTrackCutLabelNumbers[7]-.5);
- fTrackCutsPt->Fill(fTrackCutLabelNumbers[7]-.5,pt);
- fTrackCutsEta->Fill(fTrackCutLabelNumbers[7]-.5,eta);
- fTrackCutsPhi->Fill(fTrackCutLabelNumbers[7]-.5,phi);
- }
-
- // All tracks which made it up to here are classified as associateds.
-
- // Return associated.
- return 1;
-
-}
-
-//____________________________________________________________________________
-Bool_t AliAnalysisTaskDiHadronPID::SelectEvent(AliAODVertex* vertex)
-
-{
-
- //
- // Event Selection.
- //
-
- Double_t primVtx[3];
- vertex->GetXYZ(primVtx);
- if (TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2])>10.) {
- if (fVerbose>2) cout << "AliAnalysisTaskDiHadronPID::SelectEvent -> Vertex Out of Range." << endl;
- if (fVerbose>2) cout << "AliAnalysisTaskDiHadronPID::SelectEvent -> Event not selected." << endl;
- return kFALSE;
- }
- if (fVerbose>2) cout << "AliAnalysisTaskDiHadronPID::SelectEvent -> Vertex is OK." << endl;
-
- // We also wish to make centrality cut, but only if run on Pb+Pb.
- if (fBeamType=="PbPb") {
- Double_t cent = fAODHeader->GetCentrality();
- if (cent>fCentralityCutMin||cent<fCentralityCutMax) {
- if (fVerbose>2) cout<<"AliAnalysisTaskDiHadronPID::SelectEvent -> Event did not pass centaltity cut."<<endl;
- return kFALSE;
- }
- }
-
- if (fVerbose>2) cout << "AliAnalysisTaskDiHadronPID::SelectEvent -> Event selected." << endl;
- return kTRUE;
-
-}
-
-//_____________________________________________________________________________
-void AliAnalysisTaskDiHadronPID::FillGlobalTracksArray() {
-
- // Initialize the mapping for corresponding PID tracks. (see
- // GetGlobalTrack(AliAODTrack* track)).
- //
-
- if (!fAODEvent) {
- if (fVerbose>0) cout << "AliAnalysisTaskDiHadronPID::FillGlobalTracksArray -> ERROR: fAODEvent not set." << endl;
- return;
- }
-
- fGlobalTracks = new TObjArray();
- AliAODTrack* track = 0x0;
-
- for (Int_t iTrack = 0; iTrack < fAODEvent->GetNumberOfTracks(); iTrack++) {
-
- track = fAODEvent->GetTrack(iTrack);
-
- // I.e., if it does NOT pass the filtermask.
- if (!(track->TestFilterMask(1<<7))) {
- if (track->GetID()>-1) fGlobalTracks->AddAtAndExpand(track,track->GetID());
- //if (track->GetID()<1) cout<<"Track ID: "<<track->GetID()<<" Partner ID: "<<(-track->GetID()-1)<<endl;
- }
-
- }
-
-}
-
-//_____________________________________________________________________________
-AliAODTrack* AliAnalysisTaskDiHadronPID::GetGlobalTrack(AliAODTrack* track) {
-
- //
- // Returns the "parner track" of track, which contains the pid information.
- //
-
- AliAODTrack* partner = 0x0;
-
- partner = (AliAODTrack*)(fGlobalTracks->At(-track->GetID()-1));
-
- if (!partner&&(fVerbose>3)) cout<<"AliAnalysisTaskDiHadronPID::GetGlobalTrack -> No Global Track Found!"<<endl;
-
- return partner;
-
-}
-
-//_____________________________________________________________________________
-Double_t AliAnalysisTaskDiHadronPID::PhiRange(Double_t DPhi)
-
-{
- //
- // Puts the argument in the range [-pi/2,3 pi/2].
- //
-
- if (DPhi < -TMath::Pi()/2) DPhi += 2*TMath::Pi();
- if (DPhi > 3*TMath::Pi()/2) DPhi -= 2*TMath::Pi();
-
- return DPhi;
-
-}
-
-//_____________________________________________________________________________
-Int_t AliAnalysisTaskDiHadronPID::ConvertPdgCode(Int_t pdgcode) {
-
- //
- // Converts the pdg code to the bin number (0-5)
- //
-
- if (pdgcode==211) return 0; // Pi +
- if (pdgcode==-211) return 1; // Pi -
- if (pdgcode==321) return 2; // K +
- if (pdgcode==-321) return 3; // K -
- if (pdgcode==2212) return 4; // p +
- if (pdgcode==-2212) return 5; // p -
-
- return -999; // Any other particle.
-
-}
-
-//_____________________________________________________________________________
-void AliAnalysisTaskDiHadronPID::UserExec(Option_t *)
-
-{
- //
- // UserExec.
- //
-
- // Input the event.
- fAODEvent = dynamic_cast<AliAODEvent*>(InputEvent());
- if (!fAODEvent) {
- if (fVerbose>0) cout << "AliAnalysisTaskDiHadronPID::UserExec -> ERROR: No AliAODEvent pointer could be created." << endl;
- return;
- }
-
- // Get the event header.
- fAODHeader = fAODEvent->GetHeader();
- if (!fAODHeader) {
- if (fVerbose>0) cout << "AliAnalysisTaskDiHadronPID::UserExec -> ERROR: No AliAODHeader pointer could be created."<<endl;
- return;
- }
-
- // Get the event vertex.
- fAODVertex = fAODEvent->GetPrimaryVertex();
- if (!fAODVertex) {
- if (fVerbose>0) cout << "AliAnalysisTaskDiHadronPID::UserExec -> ERROR: No AliAODVertex pointer could be created." << endl;
- return;
- }
-
- // Get the MC tracks if ran on MC.
- if (fMC) {
- fMCTracks = dynamic_cast<TClonesArray*>(fAODEvent->FindListObject(AliAODMCParticle::StdBranchName()));
- if (!fMCTracks) {
- if (fVerbose>0) {
- cout<<"AliAnalysisTaskDiHadronPID::UserExec -> ERROR: No MC particles found."<<endl;
- return;
- }
- }
- }
-
- AliTPCPIDResponse& TPCPIDResponse = fPIDResponse->GetTPCResponse();
-
- // See if the event passes the event selection.
- if (!SelectEvent(fAODVertex)) return;
-
- // Display basic event information.
- if ((fVerbose>2)) cout << endl;
- if ((fVerbose>2)&&(fBeamType=="PbPb")) cout << "Event centrality: " << fAODHeader->GetCentrality() << endl;
- if ((fVerbose>2)) cout << "Event Vertex Z: " << fAODVertex->GetZ() << endl;
- if ((fVerbose>2)) cout << "Event tracks in AOD: " << fAODEvent->GetNumberOfTracks() << endl;
- if ((fVerbose>2)) cout << endl;
-
- // Filling Event QA plots.
- if (fBeamType=="PbPb") fCentrality->Fill(fAODHeader->GetCentrality());
- fVertexZ->Fill(fAODVertex->GetZ());
-
- // Fill the TObjArray which holds Global tracks.
- FillGlobalTracksArray();
-
- // Create object arrays for triggers and associateds.
- TObjArray *triggers = new TObjArray();
- TObjArray *associateds = new TObjArray();
-
- // Loop over MC tracks to fill the spectra if ran on MC.
- if (fMC) {
- for (Int_t iTrack=0; iTrack<fMCTracks->GetEntriesFast(); iTrack++) {
-
- AliAODMCParticle* mcparticle = (AliAODMCParticle*)fMCTracks->At(iTrack);
- Double_t mcPt = mcparticle->Pt();
- Double_t mcEta = mcparticle->Eta();
- //Double_t mcPhi = mcparticle->Phi();
- Double_t mcY = mcparticle->Y();
- Int_t mcPDG = mcparticle->PdgCode();
- Int_t mcSpeciesBin = ConvertPdgCode(mcPDG);
- if (mcSpeciesBin==-999) continue;
-
- if (mcparticle->IsPhysicalPrimary()) {
- fPtEtaDistrMCPrim[mcSpeciesBin]->Fill(mcPt,mcEta);
- if (TMath::Abs(mcY)<fMaxRap) fPtRapDistrMCPrimRapCut[mcSpeciesBin]->Fill(mcPt,mcY);
- } else {
- fPtEtaDistrMCSec[mcSpeciesBin]->Fill(mcPt,mcEta);
- if (TMath::Abs(mcY)<fMaxRap) fPtRapDistrMCSecRapCut[mcSpeciesBin]->Fill(mcPt,mcY);
- }
- }
- }
-
- // In this loop the triggers and associateds will be identified, track QA and PID QA histograms will be filled.
- for (Int_t iTrack = 0; iTrack < fAODEvent->GetNumberOfTracks(); iTrack++) {
-
- // Obtain a pointer to the track.
- fAODTrack = fAODEvent->GetTrack(iTrack);
- if (!fAODTrack&&(fVerbose>0)) {
- cout << "AliAnalysisTaskDiHadronPID::UserExec -> ERROR: Track object not found." << endl;
- continue;
- }
-
- //if (TMath::Abs(fAODTrack->Eta())>0.8&&fAODTrack->Pt()>0.5&&TMath::Abs(fAODTrack->Y(AliAODTrack::kProton))<0.5) cout<<"Eta: "<<fAODTrack->Eta()<<" Pt: "<<fAODTrack->Pt()<<" Ypion: "<<fAODTrack->Y(AliAODTrack::kPion)<<" Ykaon: "<<fAODTrack->Y(AliAODTrack::kKaon)<<" Yproton: "<<fAODTrack->Y(AliAODTrack::kProton)<<endl;
-
- // Find the track classification.
- Int_t tracktype = ClassifyTrack(fAODTrack);
-
- if (tracktype==0) {
- continue;
- }
-
- if (tracktype==1) {
- if (fVerbose>3) cout<<"Track added to associated buffer."<<endl;
- associateds->AddLast(fAODTrack);
- fEtaSpectrumAssoc->Fill(fAODTrack->Eta(),fAODTrack->Pt());
- fPhiSpectrumAssoc->Fill(fAODTrack->Phi(),fAODTrack->Pt());
- }
-
- if (tracktype==2) {
- if (fVerbose>3) cout<<"Track added to trigger buffer."<<endl;
- triggers->AddLast(fAODTrack);
- fEtaSpectrumTrig->Fill(fAODTrack->Eta());
-
- }
-
- if (tracktype==1) {
-
- // Fill the PID QA histograms for the associateds
- AliAODTrack* globaltrack = GetGlobalTrack(fAODTrack);
- Double_t mom, nSigma;
-
- Double_t pt = fAODTrack->Pt();
- Double_t eta = fAODTrack->Eta();
- const Int_t ptbin = (Int_t)(2*pt);
-
- //cout<<pt<<" "<<ptbin<<endl;
-
- mom = globaltrack->GetTPCmomentum();
- nSigma = fPIDResponse->NumberOfSigmasTPC(globaltrack,AliPID::kProton);
- fTPCnSigmaProton->Fill(mom,nSigma);
- nSigma = fPIDResponse->NumberOfSigmasTPC(globaltrack,AliPID::kPion);
- fTPCnSigmaPion->Fill(mom,nSigma);
- nSigma = fPIDResponse->NumberOfSigmasTPC(globaltrack,AliPID::kKaon);
- fTPCnSigmaKaon->Fill(mom,nSigma);
-
- fTPCSignal->Fill(eta,pt,globaltrack->GetTPCsignal());
-
- mom =globaltrack->P();
- nSigma = fPIDResponse->NumberOfSigmasTOF(globaltrack,AliPID::kProton);
- fTOFnSigmaProton->Fill(mom,nSigma);
- nSigma = fPIDResponse->NumberOfSigmasTOF(globaltrack,AliPID::kPion);
- fTOFnSigmaPion->Fill(mom,nSigma);
- nSigma = fPIDResponse->NumberOfSigmasTOF(globaltrack,AliPID::kKaon);
- fTOFnSigmaKaon->Fill(mom,nSigma);
-
- fTOFSignal->Fill(eta,pt,globaltrack->GetTOFsignal());
-
- // Fill the inclusives.
- Double_t TPCmom = globaltrack->GetTPCmomentum();
- Double_t TPCsignal = globaltrack->GetTPCsignal();
- Double_t expectedTPCsignalPion = TPCPIDResponse.GetExpectedSignal(TPCmom,AliPID::kPion);
- Double_t expectedTPCsignalKaon = TPCPIDResponse.GetExpectedSignal(TPCmom,AliPID::kKaon);
- Double_t expectedTPCsignalProton = TPCPIDResponse.GetExpectedSignal(TPCmom,AliPID::kProton);
-
- Double_t TOFsignal = globaltrack->GetTOFsignal();
- Double_t times[AliPID::kSPECIES];
- globaltrack->GetIntegratedTimes(times);
- Double_t expectedTOFsignalPion = times[AliPID::kPion];
- Double_t expectedTOFsignalKaon = times[AliPID::kKaon];
- Double_t expectedTOFsignalProton = times[AliPID::kProton];
-
- Double_t TPCsignalSubtracted = TPCsignal - expectedTPCsignalPion;
- Double_t TOFsignalSubtracted = TOFsignal - expectedTOFsignalPion;
- fInclusiveTPCTOF[0][ptbin]->Fill(TPCsignalSubtracted,TOFsignalSubtracted,eta);
- if (fAODTrack->Y(AliAODTrack::kPion)<fMaxRap) fInclusiveTPCTOFRap[0][ptbin]->Fill(TPCsignalSubtracted,TOFsignalSubtracted,fAODTrack->Y(AliAODTrack::kPion));
-
- TPCsignalSubtracted = TPCsignal - expectedTPCsignalKaon;
- TOFsignalSubtracted = TOFsignal - expectedTOFsignalKaon;
- fInclusiveTPCTOF[1][ptbin]->Fill(TPCsignalSubtracted,TOFsignalSubtracted,eta);
- if (fAODTrack->Y(AliAODTrack::kKaon)<fMaxRap) fInclusiveTPCTOFRap[1][ptbin]->Fill(TPCsignalSubtracted,TOFsignalSubtracted,fAODTrack->Y(AliAODTrack::kKaon));
-
- TPCsignalSubtracted = TPCsignal - expectedTPCsignalProton;
- TOFsignalSubtracted = TOFsignal - expectedTOFsignalProton;
- fInclusiveTPCTOF[2][ptbin]->Fill(TPCsignalSubtracted,TOFsignalSubtracted,eta);
- if (fAODTrack->Y(AliAODTrack::kProton)<fMaxRap) fInclusiveTPCTOFRap[2][ptbin]->Fill(TPCsignalSubtracted,TOFsignalSubtracted,fAODTrack->Y(AliAODTrack::kProton));
-
- // Fill the MC reconstructed spectra histograms.
- if (fMC) {
-
- Int_t aodlabel = TMath::Abs(fAODTrack->GetLabel());
- AliAODMCParticle* mcparticle = (AliAODMCParticle*)fMCTracks->At(aodlabel);
- Int_t mcPDG = mcparticle->PdgCode();
- Int_t mcSpeciesBin = ConvertPdgCode(mcPDG);
-
- // Only fill hisotos for pions, kaons and protons.
- if (mcSpeciesBin!=-999) {
- Double_t dataPt = fAODTrack->Pt();
- Double_t dataY=-999;
- if (mcSpeciesBin==0||mcSpeciesBin==1) {dataY = fAODTrack->Y(AliAODTrack::kPion);}
- if (mcSpeciesBin==2||mcSpeciesBin==3) {dataY = fAODTrack->Y(AliAODTrack::kKaon);}
- if (mcSpeciesBin==4||mcSpeciesBin==5) {dataY = fAODTrack->Y(AliAODTrack::kProton);}
-
- //if (TMath::Abs(fAODTrack->Eta())>0.8) cout<<"Eta: "<<fAODTrack->Eta()<<" Ypion: "<<fAODTrack->Y(AliAODTrack::kPion)<<" Ykaon: "<<fAODTrack->Y(AliAODTrack::kKaon)<<" Yproton: "<<fAODTrack->Y(AliAODTrack::kProton)<<endl;
-
- if (mcparticle->IsPhysicalPrimary()) {
- fPtEtaDistrDataPrim[mcSpeciesBin]->Fill(dataPt,eta);
- if (TMath::Abs(dataY)<fMaxRap) fPtRapDistrDataPrimRapCut[mcSpeciesBin]->Fill(dataPt,dataY);
- } else {
- fPtEtaDistrDataSec[mcSpeciesBin]->Fill(dataPt,eta);
- if (TMath::Abs(dataY)<fMaxRap) fPtRapDistrDataSecRapCut[mcSpeciesBin]->Fill(dataPt,dataY);
- }
- }
- }
-
-
-
- }
- }
-
- // In This Loop the di-hadron correlation will be made.
- Double_t histoFill[4]; // {Dphi, Deta, TPC signal, TOF signal}
- AliAODTrack* currentTrigger = 0x0;
- AliAODTrack* currentAssociated = 0x0;
- AliAODTrack* currentAssociatedGlobal = 0x0;
-
- for (Int_t iTrig = 0; iTrig < triggers->GetEntriesFast(); iTrig++){
-
- currentTrigger = (AliAODTrack*)(triggers->At(iTrig));
-
- for (Int_t iAssoc = 0; iAssoc < associateds->GetEntriesFast(); iAssoc++) {
-
- currentAssociated = (AliAODTrack*)(associateds->At(iAssoc));
- currentAssociatedGlobal = GetGlobalTrack(currentAssociated);
-
- Double_t pt = currentAssociated->Pt();
- histoFill[0] = PhiRange(currentTrigger->Phi() - currentAssociated->Phi());
- histoFill[1] = currentTrigger->Eta() - currentAssociated->Eta();
-
- // Is there a caveat here when Pt = 5.00000000?
- const Int_t ptbin = (Int_t)(2*pt);
- // cout<<"pt: "<<pt<<" ptbin: "<<ptbin<<endl; // Works OK!
-
- if (currentAssociatedGlobal) {
-
- // Get TPC (expected) signals.
- Double_t TPCmom = currentAssociatedGlobal->GetTPCmomentum();
- Double_t TPCsignal = currentAssociatedGlobal->GetTPCsignal();
- Double_t expectedTPCsignalPion = TPCPIDResponse.GetExpectedSignal(TPCmom,AliPID::kPion);
- Double_t expectedTPCsignalKaon = TPCPIDResponse.GetExpectedSignal(TPCmom,AliPID::kKaon);
- Double_t expectedTPCsignalProton = TPCPIDResponse.GetExpectedSignal(TPCmom,AliPID::kProton);
-
- // Get TOF (expected) signals.
- Double_t TOFsignal = currentAssociatedGlobal->GetTOFsignal();
- Double_t times[AliPID::kSPECIES];
- currentAssociatedGlobal->GetIntegratedTimes(times);
- Double_t expectedTOFsignalPion = times[AliPID::kPion];
- Double_t expectedTOFsignalKaon = times[AliPID::kKaon];
- Double_t expectedTOFsignalProton = times[AliPID::kProton];
-
- // Fill the histograms.
- histoFill[2] = TPCsignal - expectedTPCsignalPion;
- histoFill[3] = TOFsignal - expectedTOFsignalPion;
- fDiHadronTPCTOF[0][ptbin]->Fill(histoFill);
-
- histoFill[2] = TPCsignal - expectedTPCsignalKaon;
- histoFill[3] = TOFsignal - expectedTOFsignalKaon;
- fDiHadronTPCTOF[1][ptbin]->Fill(histoFill);
-
- histoFill[2] = TPCsignal - expectedTPCsignalProton;
- histoFill[3] = TOFsignal - expectedTOFsignalProton;
- fDiHadronTPCTOF[2][ptbin]->Fill(histoFill);
-
- fDiHadron->Fill(histoFill[0],histoFill[1],pt);
-
- }
- }
- }
-
- // In this loop we calculate the mixed events.
- if (fCalculateMixedEvents) {
-
- // Loop over the trigger buffer.
- if (fVerbose>3) cout << "AliAnalysisTaskDiHadronPID::UserExec -> Mixing the events with "<<fTrigBufferSize<<" triggers from the buffer." <<endl;
- if (fVerbose>3) cout << "AliAnalysisTaskDiHadronPID::UserExec -> Buffer size: "<<fTrigBufferIndex<<endl;
-
- for (Int_t iTrig=0;iTrig<fTrigBufferSize;iTrig++) {
-
- // Check if the trigger and the associated have a reconstructed
- // vertext no further than 2cm apart.
-
- // fTrigBuffer[i][0] = z
- // fTrigBuffer[i][1] = phi
- // fTrigBuffer[i][2] = eta
- // fTrigBuffer[i][3] = p_t
-
- if (TMath::Abs(fTrigBuffer[iTrig][0]-fAODVertex->GetZ())<fVertexZMixedEvents) {
-
- if (fVerbose>3) cout<<"AliAnalysisTaskDiHadronPID::UserExec -> Mixing with trigger Z: "<<fTrigBuffer[iTrig][0]<<", Pt: "<<fTrigBuffer[iTrig][3]<<endl;
-
- for (Int_t iAssoc = 0; iAssoc < associateds->GetEntriesFast(); iAssoc++) {
-
- currentAssociated = (AliAODTrack*)(associateds->At(iAssoc));
- currentAssociatedGlobal = GetGlobalTrack(currentAssociated);
-
- if (currentAssociatedGlobal) {
-
- Double_t pt = currentAssociated->Pt();
- histoFill[0] = PhiRange(fTrigBuffer[iTrig][1] - currentAssociated->Phi());
- histoFill[1] = fTrigBuffer[iTrig][2] - currentAssociated->Eta();
-
- //const Int_t ptbin = (Int_t)(2*pt);
-
- // Get TPC (expected) signals.
- Double_t TPCmom = currentAssociatedGlobal->GetTPCmomentum();
- Double_t TPCsignal = currentAssociatedGlobal->GetTPCsignal();
- Double_t expectedTPCsignalPion = TPCPIDResponse.GetExpectedSignal(TPCmom,AliPID::kPion);
- Double_t expectedTPCsignalKaon = TPCPIDResponse.GetExpectedSignal(TPCmom,AliPID::kKaon);
- Double_t expectedTPCsignalProton = TPCPIDResponse.GetExpectedSignal(TPCmom,AliPID::kProton);
-
- // Get TOF (expected) signals.
- Double_t TOFsignal = currentAssociatedGlobal->GetTOFsignal();
- Double_t times[AliPID::kSPECIES];
- currentAssociatedGlobal->GetIntegratedTimes(times);
- Double_t expectedTOFsignalPion = times[AliPID::kPion];
- Double_t expectedTOFsignalKaon = times[AliPID::kKaon];
- Double_t expectedTOFsignalProton = times[AliPID::kProton];
-
- // Fill the histograms.
- histoFill[2] = TPCsignal - expectedTPCsignalPion;
- histoFill[3] = TOFsignal - expectedTOFsignalPion;
- //if (ptbin==1) fMixedEventsTPCTOF[0][ptbin]->Fill(histoFill);
-
- histoFill[2] = TPCsignal - expectedTPCsignalKaon;
- histoFill[3] = TOFsignal - expectedTOFsignalKaon;
- //if (ptbin==1) fMixedEventsTPCTOF[1][ptbin]->Fill(histoFill);
-
- histoFill[2] = TPCsignal - expectedTPCsignalProton;
- histoFill[3] = TOFsignal - expectedTOFsignalProton;
- //if (ptbin==1) fMixedEventsTPCTOF[2][ptbin]->Fill(histoFill);
-
- fMixedEvents->Fill(histoFill[0],histoFill[1],pt);
-
- /*
- Double_t PionFillTPC = TPCsignal - expectedTPCsignalPion;
- Double_t PionFillTOF = TOFsignal - expectedTOFsignalPion;
-
- Int_t PionMinTPC = (fDiHadronTPCTOF[0][ptbin]->GetAxis(2))->GetXmin();
- Int_t PionMaxTPC = (fDiHadronTPCTOF[0][ptbin]->GetAxis(2))->GetXmax();
- Int_t PionMinTOF = (fDiHadronTPCTOF[0][ptbin]->GetAxis(3))->GetXmin();
- Int_t PionMaxTOF = (fDiHadronTPCTOF[0][ptbin]->GetAxis(3))->GetXmax();
-
- //cout<<PionMinTPC<<"<"<<PionFillTPC<<"<"<<PionMaxTPC<<"? "<<PionMinTOF<<"<"<<PionFillTOF<<"<"<<PionMaxTOF<<"?"<<endl;
- if (PionFillTPC<PionMaxTPC&&PionFillTPC>PionMinTPC&&PionFillTOF<PionMaxTOF&&PionFillTOF>PionMinTOF) {
- fMixedEventsTPCTOFCut[0]->Fill(DPhi,DEta,pt);
- //cout<<"Yes! Histogram will be filled."<<endl;
- }
-
- Double_t KaonFillTPC = TPCsignal - expectedTPCsignalKaon;
- Double_t KaonFillTOF = TOFsignal - expectedTOFsignalKaon;
-
- Int_t KaonMinTPC = (fDiHadronTPCTOF[1][ptbin]->GetAxis(2))->GetXmin();
- Int_t KaonMaxTPC = (fDiHadronTPCTOF[1][ptbin]->GetAxis(2))->GetXmax();
- Int_t KaonMinTOF = (fDiHadronTPCTOF[1][ptbin]->GetAxis(3))->GetXmin();
- Int_t KaonMaxTOF = (fDiHadronTPCTOF[1][ptbin]->GetAxis(3))->GetXmax();
- if (KaonFillTPC<KaonMaxTPC&&KaonFillTPC>KaonMinTPC&&KaonFillTOF<KaonMaxTOF&&KaonFillTOF>KaonMinTOF) {
- fMixedEventsTPCTOFCut[1]->Fill(DPhi,DEta,pt);
- }
-
- Double_t ProtonFillTPC = TPCsignal - expectedTPCsignalProton;
- Double_t ProtonFillTOF = TOFsignal - expectedTOFsignalProton;
- Int_t ProtonMinTPC = (fDiHadronTPCTOF[2][ptbin]->GetAxis(2))->GetXmin();
- Int_t ProtonMaxTPC = (fDiHadronTPCTOF[2][ptbin]->GetAxis(2))->GetXmax();
- Int_t ProtonMinTOF = (fDiHadronTPCTOF[2][ptbin]->GetAxis(3))->GetXmin();
- Int_t ProtonMaxTOF = (fDiHadronTPCTOF[2][ptbin]->GetAxis(3))->GetXmax();
- if (ProtonFillTPC<ProtonMaxTPC&&ProtonFillTPC>ProtonMinTPC&&ProtonFillTOF<ProtonMaxTOF&&ProtonFillTOF>ProtonMinTOF) {
- fMixedEventsTPCTOFCut[2]->Fill(DPhi,DEta,pt);
- }
- */
-
- }
- }
- }
- }
-
- // Copy the triggers from the current event into the buffer.
- if (fAODVertex->GetZ()<10.) {
-
- if (fVerbose>3) cout<<"AliAnalysisTaskDiHadronPID::UserExec -> Copying "<<triggers->GetEntriesFast()<<" triggers with vertex z = "<<fAODVertex->GetZ()<<" to the buffer."<<endl;
-
- for (Int_t iTrig = 0; iTrig<triggers->GetEntriesFast(); iTrig++) {
-
- currentTrigger = (AliAODTrack*)(triggers->At(iTrig));
- if (fVerbose>3) cout<<"AliAnalysisTaskDiHadronPID::UserExec -> Trigger pt = "<<currentTrigger->Pt()<<endl;
-
- fTrigBuffer[fTrigBufferIndex][0] = fAODVertex->GetZ();
- fTrigBuffer[fTrigBufferIndex][1] = currentTrigger->Phi();
- fTrigBuffer[fTrigBufferIndex][2] = currentTrigger->Eta();
- fTrigBuffer[fTrigBufferIndex][3] = currentTrigger->Pt();
- fTrigBufferIndex++;
- if (fTrigBufferSize<fTrigBufferMaxSize) {fTrigBufferSize++;} // 250 triggers should be enough to get 10 times more data in mixed events.
- if (fTrigBufferIndex==fTrigBufferMaxSize) {fTrigBufferIndex=0;}
- }
- }
- }
-
- if (fPrintBufferSize) cout<<"AliAnalysisTaskDiHadronPID::UserExec -> Trigger buffer index: "<<fTrigBufferIndex<<", and size: "<<fTrigBufferSize<<endl;
-
- delete triggers;
- delete associateds;
-
- PostData(1,fHistoList);
-
-}
-
-//_____________________________________________________________________________
-void AliAnalysisTaskDiHadronPID::Terminate(Option_t *)
-
-{
- //
- // Terminate.
- //
-
-}
-
+++ /dev/null
-// ----------------------------------------------------------------------------
-// AliAnalysisTaskDiHadronPID.h
-// ----------------------------------------------------------------------------
-// Author: Misha Veldhoen (misha.veldhoen@cern.ch)
-// Last Rev.: May 2nd 2012. (v 8.00)
-// ----------------------------------------------------------------------------
-
-#ifndef ALIANALYSISTASKDIHADRONPID_H
-#define ALIANALYSISTASKDIHADRONPID_H
-
-#include <iostream>
-#include "AliAnalysisTaskSE.h"
-#include "THnSparse.h"
-#include "TMath.h"
-
-using namespace std;
-
-class TH1F;
-class TH2F;
-class TH3F;
-class TList;
-class TObjArray;
-class TClonesArray;
-class TString;
-
-class AliAODTrack;
-class AliAODEvent;
-class AliAODVertex;
-
-class AliPIDResponse;
-
-class AliAnalysisTaskDiHadronPID: public AliAnalysisTaskSE {
-
-public:
- // Required functions.
- AliAnalysisTaskDiHadronPID();
- AliAnalysisTaskDiHadronPID(const char *name);
- virtual ~AliAnalysisTaskDiHadronPID();
-
- virtual void UserCreateOutputObjects();
- virtual void UserExec(Option_t *option);
- virtual void Terminate(Option_t *);
-
- // Setters
- void SetVerbose(Int_t verbose) {fVerbose=verbose;}
- void SetPrintBufferSize(Bool_t printbuffersize=kTRUE) {fPrintBufferSize=printbuffersize;}
- void SetCalculateMixedEvents(Bool_t mixedevents=kTRUE) {fCalculateMixedEvents=mixedevents;}
- void SetMC(Bool_t mc=kTRUE) {fMC = mc;}
- void SetBeamType(TString beamtype) {
- if ((beamtype!="pp")&&(beamtype!="PbPb")) {
- cout<<"SetBeamType -> Beamtype must be pp or PbPb"<<endl;
- return;
- }
- fBeamType=beamtype;
- }
-
- void SetMaxEta(Double_t maxeta) {
- if (TMath::Abs(maxeta)>0.9) {
- cout<<"SetMaxEta -> |eta| must be < 0.9"<<endl;
- return;
- }
- fMaxEta = maxeta;
- }
-
- void SetMaxRapidityInInclusiveSpectra(Double_t maxrap) {fMaxRap=maxrap;}
-
- void SetMaxPlotEta(Double_t maxploteta) {
- if (TMath::Abs(maxploteta)>1.0) {
- cout<<"SetMaxPlotEta -> |eta| must be < 1.0"<<endl;
- return;
- }
- fMaxPlotEta = maxploteta;
- }
-
- void SetMaxPt(Double_t maxpt) {
- if (maxpt<5.) {
- cout<<"SetMaxPt -> Maximum pT must be > 5.0 GeV/c."<<endl;
- return;
- }
- fMaxPt = maxpt;
- }
-
- void SetNEtaBins(Int_t netabins) {
- if (netabins<1||netabins>72) {
- cout<<"SetNEtaBins -> Number of bins must be between 1 and 72"<<endl;
- return;
- }
- fNEtaBins = netabins;
- }
-
- void SetNPhiBins(Int_t nphibins) {
- if (nphibins<1||nphibins>72) {
- cout<<"SetNPhiBins -> Number of bins must be between 1 and 72"<<endl;
- return;
- }
- fNPhiBins = nphibins;
-
- }
-
- void SetVertexZMixedEvents(Double_t vertexzmixedevents) {
- if (vertexzmixedevents<0.||vertexzmixedevents>10.) {
- cout<<"SetVertexZMixedEvents -> must be 0 < z < 10"<<endl;
- return;
- }
- fVertexZMixedEvents=vertexzmixedevents;
- }
-
- void SetZoomed(Bool_t zoomed=kTRUE) {fZoomed=zoomed;}
- void SetDoDCACut(Bool_t dodcacut=kTRUE) {fDoDCACut=dodcacut;}
- void SetDoITSCut(Bool_t doitscut=kTRUE) {fDoITSCut=doitscut;}
- void SetDemandNoMismatch(Bool_t demandnomismatch=kTRUE) {fDemandNoMismatch=demandnomismatch;}
- void SetTrigBufferMaxSize(Int_t trigbuffermaxsize) {
- if (trigbuffermaxsize<10||trigbuffermaxsize>25000) {
- cout<<"SetTrigBufferMaxSize -> Max buffer size must be between 10 and 25000."<<endl;
- return;
- }
- fTrigBufferMaxSize=trigbuffermaxsize;
- }
-
- void SetCentralityCut(Double_t centralitycutmax, Double_t centralitycutmin) {
- if (centralitycutmax<0.) {
- cout<<"SetCentralityCut -> Centrality cannot be lower than 0."<<endl;
- return;
- }
- if (centralitycutmin<centralitycutmax) {
- cout<<"SetCentralityCut -> Maximum centrality needs to be smaller than the minimum centrality. (It's confusing I know)"<<endl;
- return;
- }
- if (centralitycutmin>100.) {
- cout<<"SetCentralityCut -> Minimum centrality cannot exceed 100%."<<endl;
- return;
- }
- fCentralityCutMax=centralitycutmax;
- fCentralityCutMin=centralitycutmin;
-
- }
-
- // Getters
- Bool_t GetVerbose() {return fVerbose;}
- Bool_t GetPrintBufferSize() {return fPrintBufferSize;}
- Bool_t GetCalculateMixedEvents() {return fCalculateMixedEvents;}
- TString GetBeamType() {return fBeamType;}
- Double_t GetMaxEta() {return fMaxEta;}
- Double_t GetMaxPlotEta() {return fMaxPlotEta;}
- Double_t GetMaxPt() {return fMaxPt;}
- Int_t GetNEtaBins() {return fNEtaBins;}
- Int_t GetNPhiBins() {return fNPhiBins;}
- Double_t GetVertexZMixedEvents() {return fVertexZMixedEvents;}
- Bool_t GetZoomed() {return fZoomed;}
- Bool_t GetDoDCACut() {return fDoDCACut;}
- Bool_t GetDoITSCut() {return fDoITSCut;}
- Bool_t GetDemandNoMismatch() {return fDemandNoMismatch;}
- Int_t GetTrigBufferMaxSize() {return fTrigBufferMaxSize;}
- Double_t GetCentralityCutMax() {return fCentralityCutMax;}
- Double_t GetCentralityCutMin() {return fCentralityCutMin;}
-
-private:
- // Private Functions.
- AliAnalysisTaskDiHadronPID(const AliAnalysisTaskDiHadronPID&); // NOT IMPLEMENTED.
- AliAnalysisTaskDiHadronPID& operator=(const AliAnalysisTaskDiHadronPID&); // NOT IMPLEMENTED.
-
- void FillGlobalTracksArray();
- AliAODTrack* GetGlobalTrack(AliAODTrack* track);
-
- Bool_t SelectEvent(AliAODVertex *vertex);
- Int_t ClassifyTrack(AliAODTrack* track);
-
- Double_t PhiRange(Double_t DPhi);
- Int_t ConvertPdgCode(Int_t pdgcode);
-
-private:
- // PID object.
- AliPIDResponse *fPIDResponse; //! PID Response Handler.
-
- // Event and Track related objects.
- AliAODEvent *fAODEvent; //! The AOD Event.
- AliAODHeader *fAODHeader; //! The AOD Header.
- AliAODVertex *fAODVertex; //! The AOD Vertex.
-
- AliAODTrack *fAODTrack; //! Current AOD Track.
-
- TObjArray *fGlobalTracks; //! Partner Tracks.
- TClonesArray *fMCTracks; //! MC tracks, indexed by their track label.
-
- // HISTOGRAMS.
-
- // Event QA plots.
- TH1F *fCentrality; //! Centrality Histogram.
- TH1F *fVertexZ; //! Vertex Z position.
-
- // Track QA plots.
- TH2F *fDCA; //! DCA XY vs Z before DCA cut.
- TH2F *fDCAZoomed; //!
- TH2F *fDCAZoomedTwice; //!
- TH2F *fDCACut; //! DCA XY vs Z after DCA cut (if performed!).
- TH2F *fDCAZoomedCut; //!
- TH2F *fDCAZoomedTwiceCut; //!
-
- TH1F *fITSHits; //! 3 bins, [no hits in first 2 layers, 1 hit, 2 hits]
-
- TH1F *fTrackCutsCount; //! Counts of used tracks after cuts
- TH2F *fTrackCutsPt; //! pT spectrum after cuts.
- TH2F *fTrackCutsEta; //! eta spectrum after cuts.
- TH2F *fTrackCutsPhi; //! phi spectrum after cuts.
-
- TH1F *fEtaSpectrumTrig; //! eta spectrum of triggers (pT > 5.0 tracks, trigger track cuts.)
- TH2F *fEtaSpectrumAssoc; //! eta spectrum of associateds as a function of pT
- TH2F *fPhiSpectrumAssoc; //! phi spectrum of associateds as a function of pT
-
- // PID QA plots.
- TH2F *fTPCnSigmaProton; //! TPC nSigma plot for Protons.
- TH2F *fTOFnSigmaProton; //! TOF nSigma plot for Protons.
- TH2F *fTPCnSigmaPion; //! TPC nSigma plot for Pions.
- TH2F *fTOFnSigmaPion; //! TOF nSigma plot for Pions.
- TH2F *fTPCnSigmaKaon; //! TPC nSigma plot for Kaons.
- TH2F *fTOFnSigmaKaon; //! TOF nSigma plot for Kaons.
-
- TH3F *fTPCSignal; //! TPC signal (pt,eta).
- TH3F *fTOFSignal; //! TOF signal (pt,eta).
- TH3F *fInclusiveTPCTOF[3][10]; //! inclusive TPC-TOF histogram as a function of pT (and eta)
- TH3F *fInclusiveTPCTOFRap[3][10];//! inclusive TPC-TOF histogram as a function of pT and rapidity, with additional rapidity cut.
-
- // Efficiency Plots (Monte Carlo)
- TH2F *fPtEtaDistrDataPrim[6]; //! pT distribution of physical primaries [species: pi+,pi-,K+,K-,p,pbar].
- TH2F *fPtEtaDistrDataSec[6]; //!
-
- TH2F *fPtRapDistrDataPrimRapCut[6];//! pT distribution of physical primaries [species: pi+,pi-,K+,K-,p,pbar].
- TH2F *fPtRapDistrDataSecRapCut[6];//! with an additional rapidity cut.
-
- TH2F *fPtEtaDistrMCPrim[6]; //! pT distribution of MCParticles. [species: pi+,pi-,K+,K-,p,pbar].
- TH2F *fPtEtaDistrMCSec[6]; //!
-
- TH2F *fPtRapDistrMCPrimRapCut[6];//! pT distribution of MCParticles. [species: pi+,pi-,K+,K-,p,pbar].
- TH2F *fPtRapDistrMCSecRapCut[6]; //! with an additional rapidity cut.
-
- TH3F *fDiHadronMC[6]; //! DPhiDEta Plot per species [species]
-
- // Di-Hadron Correlations.
- TH3F *fDiHadron; //! regular di-hadron correlation, accepting all associateds.
- THnSparseF *fDiHadronTPCTOF[3][10]; //! Di-Hadron correlations with both TPC and TOF signal.
-
- // Mixed Events.
- TH3F *fMixedEvents; //! Mixed Events, associated track cuts.
- //TH3F *fMixedEventsTPCTOFCut[3]; //! For every species seperately we keep track of mixed events.
- THnSparseF *fMixedEventsTPCTOF[3][10]; //!
-
- // List of Histograms.
- TList *fHistoList; //! List of Histograms.
-
- // Analysis Task Configuration Variables.
- Bool_t fCalculateMixedEvents; //
- TString fBeamType; // pp or PbPb
- Bool_t fMC; // runs over MC.
- Double_t fMaxEta; // Q: Do we need to take extra care of the binning?
- Double_t fMaxPlotEta; //
- Double_t fMaxRap; // Max rapidity, applied to the inclusive spectra.
- Double_t fMaxPt; //
- Int_t fNEtaBins; // Number of bins in eta
- Int_t fNPhiBins; // Number of bins in phi
- Double_t fVertexZMixedEvents; // Events with a vertex z difference smaller than
- // this number (standard 2cm) will be mixed.
-
- Double_t fCentralityCutMax; // Maximum centrality (standard 0%)
- Double_t fCentralityCutMin; // Minimum centrality (standard 10%)
- Bool_t fZoomed; //
- Bool_t fDoITSCut; // Cut the tracks with not at least one SPD hit.
- Bool_t fDoDCACut; // Perform a DCA cut to get rid of secondaries.
- Bool_t fDemandNoMismatch; //
-
- Int_t fTrackCutLabelNumbers[8]; // Track Cut labels.
-
- // Level of verbal output.
- // 0 -> No output.
- // 1 -> Only error messages.
- // 2 -> Information about output creation (beginning of the job)
- // 3 -> Event information.
- // 4 -> Track information.
- Int_t fVerbose; //
- Bool_t fPrintBufferSize; //
-
- // Trigger buffer.
- Double_t fTrigBuffer[25000][4]; //!
- Int_t fTrigBufferIndex; //!
- Int_t fTrigBufferSize; //!
- Int_t fTrigBufferMaxSize; //!
-
-
- ClassDef(AliAnalysisTaskDiHadronPID,1);
-
-};
-
-#endif
-
--- /dev/null
+// -------------------------------------------------------------------------
+// Object managing event cuts, and holding QA histograms.
+// -------------------------------------------------------------------------
+// Author: Misha Veldhoen (misha.veldhoen@cern.ch)
+
+#include <iostream>
+#include "TList.h"
+#include "TMath.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TNamed.h"
+#include "TIterator.h"
+#include "AliAODEvent.h"
+#include "AliAODHeader.h"
+#include "AliAODVertex.h"
+#include "AliCentrality.h"
+#include "AliAnalysisManager.h"
+#include "AliInputEventHandler.h"
+#include "AliAODEventCutsDiHadronPID.h"
+
+using namespace std;
+
+ClassImp(AliAODEventCutsDiHadronPID);
+
+// -------------------------------------------------------------------------
+AliAODEventCutsDiHadronPID::AliAODEventCutsDiHadronPID():
+ TNamed(),
+ fIsPbPb(kTRUE),
+ fIsMC(kFALSE),
+ fTrigger(AliVEvent::kMB),
+ fMinCentrality(5.),
+ fMaxCentrality(0.),
+ fCentralityEstimator("V0M"),
+ fMaxVertexZ(10.),
+ fMinRefMult(0),
+ fTestTrigger(kFALSE),
+ fTestCentrality(kFALSE),
+ fTestVertexZ(kFALSE),
+ fTestMinRefMult(kFALSE),
+ fSelectedEventQAHistos(0x0),
+ fAllEventQAHistos(0x0),
+ fHistTrigger(0x0),
+ fHistRefMultiplicity(0x0),
+ fHistCentrality(0x0),
+ fHistCentralityQuality(0x0),
+ fHistVertexZ(0x0),
+ fDebug(0)
+
+{
+
+ //
+ // Default Constructor
+ //
+
+ cout<<"AliAODEventCutsDiHadronPID Default Constructor Called."<<endl;
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+}
+
+// -------------------------------------------------------------------------
+AliAODEventCutsDiHadronPID::AliAODEventCutsDiHadronPID(const char* name):
+ TNamed(name,"AOD Event Cuts"),
+ fIsPbPb(kTRUE),
+ fIsMC(kFALSE),
+ fTrigger(AliVEvent::kMB),
+ fMinCentrality(5.),
+ fMaxCentrality(0.),
+ fCentralityEstimator("V0M"),
+ fMaxVertexZ(10.),
+ fMinRefMult(0),
+ fTestTrigger(kFALSE),
+ fTestCentrality(kFALSE),
+ fTestVertexZ(kFALSE),
+ fTestMinRefMult(kFALSE),
+ fSelectedEventQAHistos(0x0),
+ fAllEventQAHistos(0x0),
+ fHistTrigger(0x0),
+ fHistRefMultiplicity(0x0),
+ fHistCentrality(0x0),
+ fHistCentralityQuality(0x0),
+ fHistVertexZ(0x0),
+ fDebug(0)
+
+{
+
+ //
+ // Named Constructor
+ //
+
+ cout<<"AliAODEventCutsDiHadronPID Named Constructor Called."<<endl;
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+}
+
+// -------------------------------------------------------------------------
+AliAODEventCutsDiHadronPID::AliAODEventCutsDiHadronPID(const AliAODEventCutsDiHadronPID &source):
+ TNamed(source),
+ fTrigger(source.fTrigger),
+ fMinCentrality(source.fMinCentrality),
+ fMaxCentrality(source.fMaxCentrality),
+ fCentralityEstimator(source.fCentralityEstimator),
+ fMaxVertexZ(source.fMaxVertexZ),
+ fTestTrigger(source.fTestTrigger),
+ fTestCentrality(source.fTestCentrality),
+ fTestVertexZ(source.fTestVertexZ),
+ fSelectedEventQAHistos(0x0),
+ fAllEventQAHistos(0x0),
+ fHistTrigger(0x0),
+ fHistRefMultiplicity(0x0),
+ fHistCentrality(0x0),
+ fHistVertexZ(0x0)
+
+{
+
+ //
+ // Copy Constructor
+ //
+
+ cout<<"AliAODEventCutsDiHadronPID Copy Constructor Called."<<endl;
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ //if (source.fSelectedEventQAHistos) fSelectedEventQAHistos = (TList*)source.fHistTrigger->Clone();
+ //if (source.fHistTrigger) fHistTrigger = (TH1F*)source.fHistTrigger->Clone();
+ //if (source.fHistRefMultiplicity) fHistRefMultiplicity = (TH1F*)source.fHistRefMultiplicity->Clone();
+ //if (source.fHistCentrality) fHistCentrality = (TH1F*)source.fHistCentrality->Clone();
+ //if (source.fHistVertexZ) fHistVertexZ = (TH1F*)source.fHistVertexZ->Clone();
+
+}
+
+// -------------------------------------------------------------------------
+AliAODEventCutsDiHadronPID::~AliAODEventCutsDiHadronPID() {
+
+ //
+ // Destructor
+ //
+
+ cout<<"AliAODEventCutsDiHadronPID Destructor Called."<<endl;
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ if (fSelectedEventQAHistos) delete fSelectedEventQAHistos;
+ fSelectedEventQAHistos = 0x0;
+ if (fAllEventQAHistos) delete fAllEventQAHistos;
+ fAllEventQAHistos = 0x0;
+
+}
+
+// -------------------------------------------------------------------------
+Long64_t AliAODEventCutsDiHadronPID::Merge(TCollection* list) {
+
+ //
+ // Merger.
+ //
+
+ cout<<"AliAODEventCutsDiHadronPID Merger Called."<<endl;
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ if (!list) return 0;
+ if (list->IsEmpty()) return 1;
+
+ if (!fSelectedEventQAHistos||!fAllEventQAHistos) {
+ cout<<"AliAODEventCutsDiHadronPID::Merge() - Warning, current object's histograms are missing... Generating."<<endl;
+ CreateHistos();
+ }
+
+ TIterator* iter = list->MakeIterator();
+ TObject* obj;
+
+ // List of collections
+ TList collection_fSelectedEventQAHistos;
+ TList collection_fAllEventQAHistos;
+
+ Int_t count = 0;
+
+ while ((obj = iter->Next())) {
+ AliAODEventCutsDiHadronPID* entry = dynamic_cast<AliAODEventCutsDiHadronPID*> (obj);
+ if (entry == 0) continue;
+
+ // Check if the object to be merged really has the same name! (FIXME!)
+
+ // Getting the lists from obj.
+ TList* list_fSelectedEventQAHistos = entry->GetListOfSelectedEventQAHistos();
+ TList* list_fAllEventQAHistos = entry->GetListOfAllEventQAHistos();
+
+ // Adding the retrieved lists to the collection.
+ if (list_fSelectedEventQAHistos) collection_fSelectedEventQAHistos.Add(list_fSelectedEventQAHistos);
+ if (list_fAllEventQAHistos) collection_fAllEventQAHistos.Add(list_fAllEventQAHistos);
+
+ count++;
+ }
+
+ // Merging. Note that we require the original list to exist.
+ // * Assume that if the collection happens to be empty, then nothing will happen.
+ // * All other variables are taken from the original object.
+ if (fSelectedEventQAHistos) fSelectedEventQAHistos->Merge(&collection_fSelectedEventQAHistos);
+ if (fAllEventQAHistos) fAllEventQAHistos->Merge(&collection_fAllEventQAHistos);
+
+ delete iter;
+
+ return count+1;
+
+}
+
+// -------------------------------------------------------------------------
+void AliAODEventCutsDiHadronPID::CreateHistos() {
+
+ cout<<"AliAODEventCutsDiHadronPID - Creating histograms"<<endl;
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ // Create list of Event related QA histograms (selected events).
+ fSelectedEventQAHistos = new TList();
+ fSelectedEventQAHistos->SetName("SelectedEventQAHistos");
+ fSelectedEventQAHistos->SetOwner(kTRUE);
+
+ // The same, but for all events.
+ fAllEventQAHistos = new TList();
+ fAllEventQAHistos->SetName("AllEventQAHistos");
+ fAllEventQAHistos->SetOwner(kTRUE);
+
+ // Creating arrays of pointers to the QA histos.
+ fHistTrigger = new TH1F*[2];
+ fHistRefMultiplicity = new TH1F*[2];
+ fHistCentrality = new TH1F*[2];
+ fHistCentralityQuality = new TH1F*[2];
+ fHistVertexZ = new TH1F*[2];
+
+ const char* HistType[2] = {"Selected","All"};
+
+ for (Int_t iHistType = 0; iHistType < 2; iHistType++) {
+
+ // Trigger Histogram.
+ fHistTrigger[iHistType] = new TH1F(Form("fHistTrigger%s",HistType[iHistType]),"Trigger;;Count",2,-0.5,1.5);
+ fHistTrigger[iHistType]->GetXaxis()->SetBinLabel(1,"kMB");
+ fHistTrigger[iHistType]->GetXaxis()->SetBinLabel(2,"Other");
+ if (iHistType == 0) fSelectedEventQAHistos->Add(fHistTrigger[iHistType]);
+ else fAllEventQAHistos->Add(fHistTrigger[iHistType]);
+
+ // Ref Multiplicity Histogram.
+ fHistRefMultiplicity[iHistType] = new TH1F(Form("fHistRefMultiplicity%s",HistType[iHistType]),"Reference Multiplicity;N_{tracks};Count",100,0.,10000.);
+ if (iHistType == 0) fSelectedEventQAHistos->Add(fHistRefMultiplicity[iHistType]);
+ else fAllEventQAHistos->Add(fHistRefMultiplicity[iHistType]);
+
+ // Centrality Histogram.
+ fHistCentrality[iHistType] = new TH1F(Form("fHistCentrality%s",HistType[iHistType]),"Centrality;Centrality;Count",20,0,100);
+ if (iHistType == 0) fSelectedEventQAHistos->Add(fHistCentrality[iHistType]);
+ else fAllEventQAHistos->Add(fHistCentrality[iHistType]);
+
+ // Centrality Quality.
+ fHistCentralityQuality[iHistType] = new TH1F(Form("fHistCentralityQuality%s",HistType[iHistType]),"Centrality Quality;Quality;Count",2,-0.5,1.5);
+ fHistCentralityQuality[iHistType]->GetXaxis()->SetBinLabel(1,"0");
+ fHistCentralityQuality[iHistType]->GetXaxis()->SetBinLabel(2,"Other");
+ if (iHistType == 0) fSelectedEventQAHistos->Add(fHistCentralityQuality[iHistType]);
+ else fAllEventQAHistos->Add(fHistCentralityQuality[iHistType]);
+
+ // VertexZ Histogram.
+ fHistVertexZ[iHistType] = new TH1F(Form("fHistVertexZ%s",HistType[iHistType]),"VertexZ;z (cm);Count",60,-15.,15.);
+ if (iHistType == 0) fSelectedEventQAHistos->Add(fHistVertexZ[iHistType]);
+ else fAllEventQAHistos->Add(fHistVertexZ[iHistType]);
+
+ }
+
+}
+
+// -------------------------------------------------------------------------
+Bool_t AliAODEventCutsDiHadronPID::IsSelected(AliAODEvent* event) {
+
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ if (!event) return kFALSE;
+
+ if (!fAllEventQAHistos||!fSelectedEventQAHistos) {cout<<"AliAODEventCutsDiHadronPID - Histograms were not created, you should have called CreateHistos()..."<<endl;}
+
+ // Input the event handler.
+ AliInputEventHandler* InputHandler = (AliInputEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
+ if (!InputHandler) return kFALSE;
+
+ Bool_t select = kTRUE;
+
+ // Test Trigger.
+ UInt_t trigger = InputHandler->IsEventSelected();
+ Int_t triggerselect = 0; // 0 = selected.
+ if (fTestTrigger) {
+ if (!(trigger & fTrigger)) {
+ select = kFALSE;
+ triggerselect = 1; // 1 = not selected.
+ }
+ }
+
+ AliCentrality* CurrentCentrality = 0x0;
+ Int_t CurrentCentralityQuality = -999;
+ Float_t percentile = -999.;
+
+ if (fIsPbPb) {
+
+ // Get the centrality object.
+ CurrentCentrality = event->GetCentrality();
+ if (!CurrentCentrality) select = kFALSE;
+
+ // Check the quality of the centrality estimation.
+ // If 0 then quality is OK, c.f. TOF/PbPb276/macros/TOFmatchEff.C
+ CurrentCentralityQuality = CurrentCentrality->GetQuality();
+ //cout<<"Centrality: "<<CurrentCentrality->GetCentralityPercentile(fCentralityEstimator.Data())<<" Quality: "<<CurrentCentrality->GetQuality()<<endl;
+ if (CurrentCentralityQuality) select = kFALSE;
+
+ // Test Centrality.
+ percentile = CurrentCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
+ if (fTestCentrality) {
+ if ((percentile < fMaxCentrality)||(percentile > fMinCentrality)) select = kFALSE;
+ }
+
+ }
+
+ // Get the primary vertex.
+ AliAODVertex* CurrentPrimaryVertex = event->GetPrimaryVertex();
+ if (!CurrentPrimaryVertex) select = kFALSE;
+
+ // Test Vertex Z.
+ Double_t vtxz = CurrentPrimaryVertex->GetZ();
+ if (fTestVertexZ) {
+ if (TMath::Abs(vtxz) > fMaxVertexZ) select = kFALSE;
+ }
+
+ // Get the event header.
+ AliAODHeader* CurrentHeader = event->GetHeader();
+
+ // Test minimum reference multiplicity.
+ Int_t CurrentRefMultiplicity = CurrentHeader->GetRefMultiplicity();
+ if (fTestMinRefMult) {
+ if (CurrentRefMultiplicity < fMinRefMult) select = kFALSE;
+ }
+
+ // Fill the histograms for selected events.
+ if (select) {
+ fHistTrigger[0]->Fill(triggerselect);
+ fHistRefMultiplicity[0]->Fill(CurrentHeader->GetRefMultiplicity());
+ if (fIsPbPb) fHistCentrality[0]->Fill(percentile);
+ if (fIsPbPb) fHistCentralityQuality[0]->Fill(CurrentCentralityQuality);
+ fHistVertexZ[0]->Fill(vtxz);
+ }
+
+ // Fill the histograms for all events.
+ fHistTrigger[1]->Fill(triggerselect);
+ fHistRefMultiplicity[1]->Fill(CurrentHeader->GetRefMultiplicity());
+ if (fIsPbPb) fHistCentrality[1]->Fill(percentile);
+ if (fIsPbPb) fHistCentralityQuality[1]->Fill(CurrentCentralityQuality);
+ fHistVertexZ[1]->Fill(vtxz);
+
+ cout<<"Event Selected: "<<select<<endl;
+
+ return select;
+
+}
+
+// -------------------------------------------------------------------------
+void AliAODEventCutsDiHadronPID::PrintCuts() {
+
+ // NOT IMPLEMENTED.
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ return;
+
+}
--- /dev/null
+#ifndef ALIAODEVENTCUTSDIHADRONPID_H
+#define ALIAODEVENTCUTSDIHADRONPID_H
+
+#include "TString.h"
+#include "TH1F.h"
+#include "TList.h"
+
+class AliAODEventCutsDiHadronPID : public TNamed
+
+{
+
+public:
+ AliAODEventCutsDiHadronPID(); // Default Constructor
+ AliAODEventCutsDiHadronPID(const char* name); // Named Constructor
+ AliAODEventCutsDiHadronPID(const AliAODEventCutsDiHadronPID &source); // Copy Constructor
+ virtual ~AliAODEventCutsDiHadronPID(); // Destructor
+ virtual Long64_t Merge(TCollection* list); // Merger
+
+ void CreateHistos(); // Create QA histograms
+
+public:
+ Bool_t IsSelected(AliAODEvent* event);
+ void PrintCuts();
+
+// Setters
+ void SetIsPbPb(Bool_t ispbpb = kTRUE) {fIsPbPb = ispbpb;}
+ void SetIsMC(Bool_t ismc = kTRUE) {fIsMC = ismc;}
+
+ void SetTrigger(UInt_t trigger) {
+ fTrigger = trigger;
+ fTestTrigger = kTRUE;
+ }
+
+ // Note that minCentrality is expected to be the biggest number.
+ void SetCentrality(Float_t maxCentrality, Float_t minCentrality) {
+ if (minCentrality > maxCentrality) {
+ fMinCentrality = minCentrality;
+ fMaxCentrality = maxCentrality;
+ } else {
+ fMinCentrality = maxCentrality;
+ fMaxCentrality = minCentrality;
+ }
+ fTestCentrality = kTRUE;
+ }
+
+ void SetCentralityEstimator(const char* centralityestimator) {
+ fCentralityEstimator = centralityestimator;
+ }
+ void SetMaxVertexZ(Float_t maxVertexZ) {
+ fMaxVertexZ = maxVertexZ;
+ fTestVertexZ = kTRUE;
+ }
+ void SetMinReferenceMultiplicity(Int_t minrefmult) {
+ fMinRefMult = minrefmult;
+ fTestMinRefMult = kTRUE;
+ }
+
+ void SetDebugLevel(Int_t debuglevel) {fDebug = debuglevel;}
+
+// Getters
+ Bool_t GetIsPbPb() const {return fIsPbPb;}
+ Bool_t GetIsMC() const {return fIsMC;}
+ UInt_t GetTrigger() const {return fTrigger;}
+ Float_t GetMinCentrality() const {return fMinCentrality;}
+ Float_t GetMaxCentrality() const {return fMaxCentrality;}
+ TString GetCentralityEstimator() const {return fCentralityEstimator;}
+ Float_t GetMaxVertexZ() const {return fMaxVertexZ;}
+
+ // Functions returning pointer data members aren't very safe.
+ TList* GetListOfSelectedEventQAHistos() {
+ if (fSelectedEventQAHistos) return fSelectedEventQAHistos;
+ else return 0x0;
+ }
+ TList* GetListOfAllEventQAHistos() {
+ if (fAllEventQAHistos) return fAllEventQAHistos;
+ else return 0x0;
+ }
+
+ TObject* GetHistSelectedEvents(const char* name) {return fSelectedEventQAHistos->FindObject(name);}
+ TObject* GetHistAllEvents(const char* name) {return fAllEventQAHistos->FindObject(name);}
+
+ // Cannot be made const because GetHistSelectedEvents() isn't safe.
+ Int_t GetNAcceptedEvents() {return ((TH1F*)GetHistSelectedEvents("fHistTriggerSelected"))->GetEntries();}
+
+ Int_t GetDebugLevel() const {return fDebug;}
+
+private:
+// Expected Event Details.
+ Bool_t fIsPbPb;
+ Bool_t fIsMC;
+
+// Event Cuts.
+ UInt_t fTrigger;
+ Float_t fMinCentrality;
+ Float_t fMaxCentrality;
+ TString fCentralityEstimator;
+ Float_t fMaxVertexZ;
+ Int_t fMinRefMult;
+
+// Which cuts to be checked.
+ Bool_t fTestTrigger;
+ Bool_t fTestCentrality;
+ Bool_t fTestVertexZ;
+ Bool_t fTestMinRefMult;
+
+// QA histograms (don't stream)
+ TList* fSelectedEventQAHistos;
+ TList* fAllEventQAHistos;
+ TH1F** fHistTrigger; //! Trigger
+ TH1F** fHistRefMultiplicity; //! Number of tracks
+ TH1F** fHistCentrality; //! Centrality
+ TH1F** fHistCentralityQuality; //! Centrality Quality
+ TH1F** fHistVertexZ; //! VertexZ
+
+ Int_t fDebug; // Debug flag.
+
+ ClassDef(AliAODEventCutsDiHadronPID,2);
+
+};
+
+#endif
--- /dev/null
+// -------------------------------------------------------------------------
+// Object to handle AOD Track cuts for the DiHadronPID analysis. Should
+// inherit from: ?TObject or TNamed?, so that it can be written to a file.
+//
+// Still open question: should this object hold multiple track cuts, or
+// should the analysis hold a TObjArray, with all the track cuts?
+// -------------------------------------------------------------------------
+// Author: Misha Veldhoen (misha.veldhoen@cern.ch)
+
+#include <iostream>
+
+#include "TMath.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TH3F.h"
+#include "TNamed.h"
+#include "TFormula.h"
+#include "TIterator.h"
+
+// AOD includes.
+#include "AliAODTrack.h"
+#include "AliAODEvent.h"
+#include "AliAODVertex.h"
+#include "AliAODMCParticle.h"
+
+// PID includes.
+#include "AliPID.h"
+#include "AliPIDResponse.h"
+#include "AliTPCPIDResponse.h"
+
+#include "AliTrackDiHadronPID.h"
+
+#include "AliAODTrackCutsDiHadronPID.h"
+
+using namespace std;
+
+ClassImp(AliAODTrackCutsDiHadronPID);
+
+// -------------------------------------------------------------------------
+AliAODTrackCutsDiHadronPID::AliAODTrackCutsDiHadronPID():
+ TNamed(),
+ fMinPt(0.),
+ fMaxPt(10.),
+ fFilterMask(0),
+ fMaxEta(-999.),
+ fMaxRapidity(-999.),
+ fDemandedFlags(0),
+ fMinSPDHitsForPtDeptDCAcut(0),
+ fPtDeptDCAxyCutFormula(0x0),
+ fDCAzCut(999.),
+ fIsMC(kFALSE),
+ fTestFilterMask(kFALSE),
+ fTestMaxEta(kFALSE),
+ fTestMaxRapidity(kFALSE),
+ fTestFlags(kFALSE),
+ fTestTOFmismatch(kFALSE),
+ fTestPtDeptDCAcut(kFALSE),
+ fDataTrackQAHistos(0x0),
+ fPrimRecMCTrackQAHistos(0x0),
+ fPrimGenMCTrackQAHistos(0x0),
+ fSecRecMCTrackQAHistos(0x0),
+ fSecGenMCTrackQAHistos(0x0),
+ fDebug(0)
+
+ {
+
+ //
+ // Default constructor
+ //
+
+ // Q: When an object is loaded from a TFile, then the default constructor is called.
+ // I am assuming that the valiables initialized here are overwritten by the values
+ // which were stored in the file, except for the data members which have an
+ // exclamation mark in the header file. -> Check this!
+
+ cout<<"AliAODTrackCutsDiHadronPID Default Constructor Called."<<endl;
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ for (Int_t iHistoName = 0; iHistoName < 12; iHistoName++) {
+
+ // Initialize the data histograms.
+ if (iHistoName < 3) {
+ fHistDataPt[iHistoName] = 0x0;
+ fHistDataNTracks[iHistoName] = 0x0;
+ fNTracks[iHistoName] = 0;
+
+ // DCA histograms.
+ fHistDataDCAxy[iHistoName] = 0x0;
+ fHistDataDCAz[iHistoName] = 0x0;
+
+ // Initialize PID histograms.
+ for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
+ for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
+ fHistDataPID[iHistoName][iSpecies][iPtClass] = 0x0;
+ fHistTOFMismatch[iHistoName][iSpecies][iPtClass] = 0x0;
+ }
+ }
+ }
+ // Initialize data DCA for one sigma.
+ fHistDataDCAxyOneSigma[iHistoName] = 0x0;
+
+ // Initialize MC histograms.
+ fHistPrimGenMCPt[iHistoName] = 0x0;
+ fHistPrimRecMCPt[iHistoName] = 0x0;
+ fHistPrimRecNTracks[iHistoName] = 0x0;
+
+ fHistSecGenMCPt[iHistoName] = 0x0;
+ fHistSecRecMCPt[iHistoName] = 0x0;
+
+ // Initialze MC DCA histograms
+ fHistPrimRecMCDCA[iHistoName] = 0x0;
+ fHistSecRecMCDCAMat[iHistoName] = 0x0;
+ fHistSecRecMCDCAWeak[iHistoName] = 0x0;
+
+ // Initialize other stuff.
+ fHistRequested[iHistoName] = kFALSE;
+ f3DSpectraEnabeled[iHistoName] = kFALSE;
+ fPIDHistosEnabeled[iHistoName] = kFALSE;
+ }
+
+ // TODO: It would be a bit neater to initialize all values to zero
+ // instead of the default values...
+ InitializeDefaultHistoNamesAndAxes();
+
+}
+
+// -------------------------------------------------------------------------
+AliAODTrackCutsDiHadronPID::AliAODTrackCutsDiHadronPID(const char* name):
+ TNamed(name,"AOD Track Cuts"),
+ fMinPt(0.),
+ fMaxPt(10.),
+ fFilterMask(0),
+ fMaxEta(-999.),
+ fMaxRapidity(-999.),
+ fDemandedFlags(0),
+ fMinSPDHitsForPtDeptDCAcut(0),
+ fPtDeptDCAxyCutFormula(0x0),
+ fDCAzCut(999.),
+ fIsMC(kFALSE),
+ fTestFilterMask(kFALSE),
+ fTestMaxEta(kFALSE),
+ fTestMaxRapidity(kFALSE),
+ fTestFlags(kFALSE),
+ fTestTOFmismatch(kFALSE),
+ fTestPtDeptDCAcut(kFALSE),
+ fDataTrackQAHistos(0x0),
+ fPrimRecMCTrackQAHistos(0x0),
+ fPrimGenMCTrackQAHistos(0x0),
+ fSecRecMCTrackQAHistos(0x0),
+ fSecGenMCTrackQAHistos(0x0),
+ fDebug(0)
+
+{
+
+ //
+ // Named constructor
+ //
+
+ cout<<"AliAODTrackCutsDiHadronPID Named Constructor Called."<<endl;
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ for (Int_t iHistoName = 0; iHistoName < 12; iHistoName++) {
+
+ // Initialize the data histograms.
+ if (iHistoName < 3) {
+ fHistDataPt[iHistoName] = 0x0;
+ fHistDataNTracks[iHistoName] = 0x0;
+ fNTracks[iHistoName] = 0;
+
+ // DCA histograms.
+ fHistDataDCAxy[iHistoName] = 0x0;
+ fHistDataDCAz[iHistoName] = 0x0;
+
+ // Initialize PID histograms.
+ for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
+ for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
+ fHistDataPID[iHistoName][iSpecies][iPtClass] = 0x0;
+ fHistTOFMismatch[iHistoName][iSpecies][iPtClass] = 0x0;
+ }
+ }
+ }
+ // Initialize data DCA for one sigma.
+ fHistDataDCAxyOneSigma[iHistoName] = 0x0;
+
+ // Initialize MC histograms.
+ fHistPrimGenMCPt[iHistoName] = 0x0;
+ fHistPrimRecMCPt[iHistoName] = 0x0;
+ fHistPrimRecNTracks[iHistoName] = 0x0;
+
+ fHistSecGenMCPt[iHistoName] = 0x0;
+ fHistSecRecMCPt[iHistoName] = 0x0;
+
+ // Initialze MC DCA histograms
+ fHistPrimRecMCDCA[iHistoName] = 0x0;
+ fHistSecRecMCDCAMat[iHistoName] = 0x0;
+ fHistSecRecMCDCAWeak[iHistoName] = 0x0;
+
+ // Initialize other stuff.
+ fHistRequested[iHistoName] = kFALSE;
+ f3DSpectraEnabeled[iHistoName] = kFALSE;
+ fPIDHistosEnabeled[iHistoName] = kFALSE;
+ }
+
+ InitializeDefaultHistoNamesAndAxes();
+
+}
+
+// -------------------------------------------------------------------------
+void AliAODTrackCutsDiHadronPID::InitializeDefaultHistoNamesAndAxes() {
+
+ // Initializes the histogram name conventions and the ranges of all the histograms
+ // by their default values. This method should only be called by the
+ // (default) constructor.
+
+ // TODO: User should be able to change these standard values at initialization.
+ // TODO: User should be able to retrieve all these variables with appropriate Getters.
+
+ cout<<"AliAODTrackCutsDiHadronPID - Initializing Default Histogram Names and axes..."<<endl;
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ // Setting the Pt axis for all histograms except the PID and Mismatch histograms.
+ Double_t ptaxis[47] = {0.20,0.25,0.30,0.35,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.90,0.95,1.00,
+ 1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3.0,
+ 3.2,3.4,3.6,3.8,4.0,4.2,4.4,4.6,4.8,5.0};
+ for (Int_t iPtBins = 0; iPtBins < 47; iPtBins++) {fPtAxis[iPtBins] = ptaxis[iPtBins];}
+ fNPtBins = 46;
+
+ // Setting the pt range of the five PID histogram pt classes.
+ Double_t ptboundarypid[6] = {0.5,0.7,1.0,1.7,3.0,5.0};
+ for (Int_t iPtBoundary = 0; iPtBoundary < 6; iPtBoundary++) {fPtBoundaryPID[iPtBoundary] = ptboundarypid[iPtBoundary];}
+
+ // Setting the number of bins in pt for the five PID histogram pt classes.
+ Int_t nptbinspid[5] = {8,6,14,13,10};
+ for (Int_t iPtBins = 0; iPtBins < 5; iPtBins++) {fNPtBinsPID[iPtBins] = nptbinspid[iPtBins];}
+
+ // Setting the TOF axes for the PID histograms.
+ Double_t toflowerbound[5][3] = {{-2000.,-6000.,-10000.},{-2000.,-4000.,-10000.},{-1000.,-2000.,-5000.},{-1000.,-1000.,-2500.},{-500.,-500.,-1000.}};
+ Double_t tofupperbound[5][3] = {{10000.,10000.,6000.},{10000.,8000.,6000.},{6000.,6000.,6000.},{6000.,6000.,6000.},{4000.,4000.,6000.}};
+ Int_t tofbins[5][3] = {{120,160,160},{120,120,160},{140,140,165},{140,140,170},{90,90,140}};
+ for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
+ for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
+ fTOFLowerBound[iPtClass][iSpecies] = toflowerbound[iPtClass][iSpecies];
+ fTOFUpperBound[iPtClass][iSpecies] = tofupperbound[iPtClass][iSpecies];
+ fTOFbins[iPtClass][iSpecies] = tofbins[iPtClass][iSpecies];
+ }
+ }
+
+ // Setting the TPC axes for the PID histograms.
+ Double_t tpclowerbound[5][3] = {{-20.,-50.,-100.},{-20.,-30.,-80.},{-25.,-25.,-45.},{-25.,-25.,-45.},{-25.,-20.,-20.}};
+ Double_t tpcupperbound[5][3] = {{60.,30.,20.},{60.,40.,20.},{50.,50.,25.},{45.,45.,25.},{25.,30.,30.}}; // Check highest pT bin boundaries for K,p
+ Int_t tpcbins[5][3] = {{80,80,120},{80,70,100},{75,75,70},{70,70,70},{50,50,50}};
+ for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
+ for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
+ fTPCLowerBound[iPtClass][iSpecies] = tpclowerbound[iPtClass][iSpecies];
+ fTPCUpperBound[iPtClass][iSpecies] = tpcupperbound[iPtClass][iSpecies];
+ fTPCbins[iPtClass][iSpecies] = tpcbins[iPtClass][iSpecies];
+ }
+ }
+
+ // Names for the 12 (species,charge) combinations considered.
+ const char* histoname[12] = {"AllCharged","Pos","Neg","AllPion","PosPion","NegPion","AllKaon","PosKaon","NegKaon","AllProton","PosProton","NegProton"};
+ const char* histolatex[12] = {"ch","ch^{+}","ch^{-}","#pi^{+/-}","#pi^{+}","#pi^{-}","K^{+/-}","K^{+}","K^{-}","p/#bar{p}","p","#bar{p}"};
+ for (Int_t iHistoName = 0; iHistoName < 12; iHistoName++) {
+ fHistoName[iHistoName] = histoname[iHistoName];
+ fHistoLatex[iHistoName] = histolatex[iHistoName];
+ }
+
+ // Names of the 3 species considered.
+ const char* particlename[3] = {"Pion","Kaon","Proton"};
+ for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {fParticleName[iSpecies] = particlename[iSpecies];}
+
+ // Names of the 5 pt classes considered for the PID histograms.
+ const char* ptclassname[5] = {"VeryLowPt","LowPt","MedPt","MedHighPt","HighPt"};
+ for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {fPtClassName[iPtClass] = ptclassname[iPtClass];}
+
+}
+
+// -------------------------------------------------------------------------
+AliAODTrackCutsDiHadronPID::~AliAODTrackCutsDiHadronPID() {
+
+ //
+ // Destructor
+ //
+
+ cout<<"AliAODTrackCutsDiHadronPID Destructor Called."<<endl;
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ if (fDataTrackQAHistos) delete fDataTrackQAHistos;
+ fDataTrackQAHistos = 0x0;
+ if (fPrimRecMCTrackQAHistos) delete fPrimRecMCTrackQAHistos;
+ fPrimRecMCTrackQAHistos = 0x0;
+ if (fPrimGenMCTrackQAHistos) delete fPrimGenMCTrackQAHistos;
+ fPrimGenMCTrackQAHistos = 0x0;
+ if (fSecRecMCTrackQAHistos) delete fSecRecMCTrackQAHistos;
+ fSecRecMCTrackQAHistos = 0x0;
+ if (fSecGenMCTrackQAHistos) delete fSecGenMCTrackQAHistos;
+ fSecGenMCTrackQAHistos = 0x0;
+
+}
+
+// -------------------------------------------------------------------------
+Long64_t AliAODTrackCutsDiHadronPID::Merge(TCollection* list) {
+
+ //
+ // Merger.
+ //
+
+ cout<<"AliAODTrackCutsDiHadronPID Merger Called."<<endl;
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ Bool_t HistosOK = kTRUE;
+
+ if (fIsMC) {
+ if (!fPrimGenMCTrackQAHistos||!fSecGenMCTrackQAHistos||!fPrimRecMCTrackQAHistos||!fSecRecMCTrackQAHistos) {HistosOK = kFALSE;}
+ } else {
+ if (!fDataTrackQAHistos) {HistosOK = kFALSE;}
+ }
+
+ if (!HistosOK) {
+ cout<<"AliAODTrackCutsDiHadronPID::Merge() - Warning, current object's histograms are missing... Generating."<<endl;
+ CreateHistos();
+ }
+
+ if (!list) {
+ //cout<<"No list available..."<<endl;
+ return 0;
+ }
+ if (list->IsEmpty()) {
+ //cout<<"List is empty..."<<endl;
+ return 1;
+ }
+
+ //Int_t NEntries = list->GetEntries();
+ //cout<<"Supplied list has "<<NEntries<<" entries."<<endl;
+
+ TIterator* iter = list->MakeIterator();
+ TObject* obj;
+
+ // List of collections
+ TList collection_fDataTrackQAHistos;
+ TList collection_fPrimRecMCTrackQAHistos;
+ TList collection_fPrimGenMCTrackQAHistos;
+ TList collection_fSecRecMCTrackQAHistos;
+ TList collection_fSecGenMCTrackQAHistos;
+
+ Int_t count = 0;
+
+ while ((obj = iter->Next())) {
+ AliAODTrackCutsDiHadronPID* entry = dynamic_cast<AliAODTrackCutsDiHadronPID*> (obj);
+ if (entry == 0) continue;
+
+ // Check if the object to be merged really has the same name! (FIXME!)
+
+ // Getting the lists from obj.
+ TList* list_fDataTrackQAHistos = entry->GetListOfDataQAHistos();
+ TList* list_fPrimRecMCTrackQAHistos = entry->GetListOfPrimRecMCTrackQAHistos();
+ TList* list_fPrimGenMCTrackQAHistos = entry->GetListOfPrimGenMCTrackQAHistos();
+ TList* list_fSecRecMCTrackQAHistos = entry->GetListOfSecRecMCTrackQAHistos();
+ TList* list_fSecGenMCTrackQAHistos = entry->GetListOfSecGenMCTrackQAHistos();
+
+ // Adding the retrieved lists to the collection.
+ if (list_fDataTrackQAHistos) collection_fDataTrackQAHistos.Add(list_fDataTrackQAHistos);
+ if (list_fPrimRecMCTrackQAHistos) collection_fPrimRecMCTrackQAHistos.Add(list_fPrimRecMCTrackQAHistos);
+ if (list_fPrimGenMCTrackQAHistos) collection_fPrimGenMCTrackQAHistos.Add(list_fPrimGenMCTrackQAHistos);
+ if (list_fSecRecMCTrackQAHistos) collection_fSecRecMCTrackQAHistos.Add(list_fSecRecMCTrackQAHistos);
+ if (list_fSecGenMCTrackQAHistos) collection_fSecGenMCTrackQAHistos.Add(list_fSecGenMCTrackQAHistos);
+
+ //cout<<"Entries intermediate list Data: "<<collection_fDataTrackQAHistos.GetEntries()<<endl;
+ //cout<<"Entries intermediate list RecMC: "<<collection_fPrimRecMCTrackQAHistos.GetEntries()<<endl;
+
+ count++;
+ }
+
+ // Merging. Note that we require the original list to exist.
+ // * Assume that if the collection happens to be empty, then nothing will happen.
+ // |-> This of course leads to some problems, since if the first file does not
+ // | have such a list, then the merged file will have no results...
+ // | IDEA: if fDataTrackQAHistos does not exist, then create a new one. (TO DO)
+ // * All other variables are taken from the original object.
+ if (fDataTrackQAHistos) fDataTrackQAHistos->Merge(&collection_fDataTrackQAHistos);
+ if (fPrimRecMCTrackQAHistos) fPrimRecMCTrackQAHistos->Merge(&collection_fPrimRecMCTrackQAHistos);
+ if (fPrimGenMCTrackQAHistos) fPrimGenMCTrackQAHistos->Merge(&collection_fPrimGenMCTrackQAHistos);
+ if (fSecRecMCTrackQAHistos) fSecRecMCTrackQAHistos->Merge(&collection_fSecRecMCTrackQAHistos);
+ if (fSecGenMCTrackQAHistos) fSecGenMCTrackQAHistos->Merge(&collection_fSecGenMCTrackQAHistos);
+
+ delete iter;
+
+ return count+1;
+
+}
+
+// -------------------------------------------------------------------------
+void AliAODTrackCutsDiHadronPID::PrintCuts() {}
+
+// -------------------------------------------------------------------------
+void AliAODTrackCutsDiHadronPID::CreateHistos() {
+
+ // Function should be called by the UserCreateOutput() function of the
+ // analysis task. This function then generates all the histograms that
+ // were requested locally on the workernode. Even if case that the
+ // histograms are not filled, it is imperative that they are still
+ // created, and in the same order, otherwise problems with merging will
+ // arise.
+
+ // TODO: In principle it should never happen that if this method is called,
+ // that the lists of QA histograms already exist, so this may become a fatal error
+ // instead of a warning.
+
+ cout<<"AliAODTrackCutsDiHadronPID - Creating histograms for cuts: "<<fName<<endl;
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ if (fIsMC) {
+
+ if (!fPrimGenMCTrackQAHistos) {
+ cout<<"AliAODTrackCutsDiHadronPID - Creating Prim. Gen. MC Track QA TList..."<<endl;
+ fPrimGenMCTrackQAHistos = new TList();
+ fPrimGenMCTrackQAHistos->SetName("PrimGenMCTrackQAHistos");
+ fPrimGenMCTrackQAHistos->SetOwner(kTRUE);
+ } else {cout<<"AliAODTrackCutsDiHadronPID - Warning, Prim. Gen. MC Track QA TList was already created..."<<endl;}
+
+ if (!fSecGenMCTrackQAHistos) {
+ cout<<"AliAODTrackCutsDiHadronPID - Creating Sec. Gen. MC Track QA TList..."<<endl;
+ fSecGenMCTrackQAHistos = new TList();
+ fSecGenMCTrackQAHistos->SetName("SecGenMCTrackQAHistos");
+ fSecGenMCTrackQAHistos->SetOwner(kTRUE);
+ for (Int_t iHistoClass = 0; iHistoClass < 12; iHistoClass++) {if (fHistRequested[iHistoClass]) {InitializeGenMCHistos(iHistoClass);}}
+ } else {cout<<"AliAODTrackCutsDiHadronPID - Warning, Sec. Gen. MC Track QA TList was already created..."<<endl;}
+
+ if (!fPrimRecMCTrackQAHistos) {
+ cout<<"AliAODTrackCutsDiHadronPID - Creating Prim. Rec. MC Track QA TList..."<<endl;
+ fPrimRecMCTrackQAHistos = new TList();
+ fPrimRecMCTrackQAHistos->SetName("PrimRecMCTrackQAHistos");
+ fPrimRecMCTrackQAHistos->SetOwner(kTRUE);
+ } else {cout<<"AliAODTrackCutsDiHadronPID - Warning, Prim. Rec. MC Track QA TList was already created..."<<endl;}
+
+ if (!fSecRecMCTrackQAHistos) {
+ cout<<"AliAODTrackCutsDiHadronPID - Creating Sec. Rec. MC Track QA TList..."<<endl;
+ fSecRecMCTrackQAHistos = new TList();
+ fSecRecMCTrackQAHistos->SetName("SecRecMCTrackQAHistos");
+ fSecRecMCTrackQAHistos->SetOwner(kTRUE);
+ for (Int_t iHistoClass = 0; iHistoClass < 12; iHistoClass++) {
+ if (fHistRequested[iHistoClass]) {
+ InitializeRecMCHistos(iHistoClass);
+ }
+ }
+ } else {cout<<"AliAODTrackCutsDiHadronPID - Warning, Sec. Rec. MC Track QA TList was already created..."<<endl;}
+
+ } else {
+
+ if (!fDataTrackQAHistos) {
+ cout<<"AliAODTrackCutsDiHadronPID - Creating Data Track QA TList..."<<endl;
+ fDataTrackQAHistos = new TList();
+ fDataTrackQAHistos->SetName("DataTrackQAHistos");
+ fDataTrackQAHistos->SetOwner(kTRUE);
+ for (Int_t iHistoClass = 0; iHistoClass < 12; iHistoClass++) {
+ fHistDataDCAxyOneSigma[iHistoClass] = InitializeDCASpectrum("fHistDataDCAxyOneSigma",iHistoClass);
+ fDataTrackQAHistos->Add(fHistDataDCAxyOneSigma[iHistoClass]);
+ if (fHistRequested[iHistoClass]) { if (iHistoClass < 3) {InitializeDataHistos(iHistoClass);}}
+ }
+ } else {cout<<"AliAODTrackCutsDiHadronPID - Warning, Data QA TList was already created..."<<endl;}
+
+ }
+}
+
+// -------------------------------------------------------------------------
+void AliAODTrackCutsDiHadronPID::StartNewEvent() {
+
+ // FIXME: This method is now only suited for Data (3 histo classes only.)
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ // Resetting the counters.
+ for (Int_t iHistoClass = 0; iHistoClass < 3; iHistoClass++) {
+ fNTracks[iHistoClass] = 0;
+ }
+}
+
+// -------------------------------------------------------------------------
+void AliAODTrackCutsDiHadronPID::EventIsDone(Bool_t IsMC) {
+
+ // FIXME: This method is now only suited for Data (3 histo classes only.)
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ // Fill NTracks histos.
+
+ for (Int_t iHistoClass = 0; iHistoClass < 3; iHistoClass++) {
+
+ if (fHistRequested[iHistoClass]) {
+
+ if (IsMC) {
+ // THE FOLLOWING SHOULD NEVER HAPPEN.
+ //if (!fHistPrimRecMCPt[iHistoClass]) InitializeRecMCHistos(iHistoClass);
+ fHistPrimRecNTracks[iHistoClass]->Fill(fNTracks[iHistoClass]);
+ } else {
+ // THE FOLLOWING SHOULD NEVER HAPPEN.
+ //if (!fHistDataPt[iHistoClass]) InitializeDataHistos(iHistoClass);
+ fHistDataNTracks[iHistoClass]->Fill(fNTracks[iHistoClass]);
+ }
+ }
+ }
+}
+
+// -------------------------------------------------------------------------
+Bool_t AliAODTrackCutsDiHadronPID::IsSelectedData(AliTrackDiHadronPID* track, Double_t randomhittime) {
+
+ //
+ // Checks performed on a data track.
+ //
+
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ if (!track) return kFALSE;
+
+ if (!fDataTrackQAHistos) {cout<<"AliAODTrackCutsDiHadronPID - Histograms were not created, you should have called CreateHistos()..."<<endl;}
+
+ // Check the track cuts (NB. check function return kTRUE if track is accepted.)
+ if (!CheckPt(track->Pt())) return kFALSE;
+ if (!CheckMaxEta(track->Eta())) return kFALSE;
+ if (!CheckFilterMask(track->GetFilterMap())) return kFALSE;
+ if (!CheckFlags(track->GetFlags())) return kFALSE;
+ if (!CheckTOFmismatch(track->IsTOFmismatch())) return kFALSE;
+
+ Int_t NSPDhits = 0;
+ if (track->HasPointOnITSLayer(0)) NSPDhits++;
+ if (track->HasPointOnITSLayer(1)) NSPDhits++;
+ if (!CheckPtDeptDCACut(track->GetZAtDCA(),track->GetXYAtDCA(),track->Pt(),NSPDhits)) return kFALSE;
+
+ // Track has passed the cuts, fill QA histograms.
+ for (Int_t iHistoClass = 0; iHistoClass < 3; iHistoClass++) {
+
+ if (fHistRequested[iHistoClass]) {
+
+ // Check the charge (could be neater).
+ if ((iHistoClass == 0) && (track->Charge() == 0)) continue;
+ if ((iHistoClass == 1) && (track->Charge() <= 0)) continue;
+ if ((iHistoClass == 2) && (track->Charge() >= 0)) continue;
+
+ FillDataHistos(iHistoClass, track);
+
+ // Ignore if random hit is < -1.e20.
+ if (randomhittime > -1.e20) FillTOFMismatchHistos(iHistoClass,track,randomhittime);
+
+ fNTracks[iHistoClass]++;
+
+ }
+ }
+
+ // Track Has Passed.
+ return kTRUE;
+
+}
+
+// -------------------------------------------------------------------------
+Bool_t AliAODTrackCutsDiHadronPID::IsSelectedGeneratedMC(AliAODMCParticle* particle) {
+
+ //
+ // Checks performed on a generated MC particle.
+ //
+
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ // PDG codes for particles:
+ // pi+ 211, K+ 321, p 2212
+
+ if (!particle) return kFALSE;
+
+ if (!fPrimGenMCTrackQAHistos||!fSecGenMCTrackQAHistos) {cout<<"AliAODTrackCutsDiHadronPID - Histograms were not created, you should have called CreateHistos()..."<<endl;}
+
+ // Check the track cuts.
+ if (!CheckPt(particle->Pt())) return kFALSE;
+ if (!CheckMaxEta(particle->Eta())) return kFALSE;
+ if (!CheckRapidity(particle->Y())) return kFALSE; // NEW: rapidity cut. (not done on data)
+
+ // NOT YET IMPLEMENTED! - Histoclasses for specific particles.
+ for (Int_t iHistoClass = 0; iHistoClass < 12; iHistoClass++) {
+
+ if (fHistRequested[iHistoClass]) {
+
+ // Check the charge (could be neater).
+ if ((iHistoClass == 0) && (particle->Charge() == 0)) continue;
+ if ((iHistoClass == 1) && (particle->Charge() <= 0)) continue;
+ if ((iHistoClass == 2) && (particle->Charge() >= 0)) continue;
+ if ((iHistoClass == 3) && (TMath::Abs(particle->GetPdgCode()) != 211)) continue;
+ if ((iHistoClass == 4) && (particle->GetPdgCode()) != 211) continue;
+ if ((iHistoClass == 5) && (particle->GetPdgCode()) != -211) continue;
+ if ((iHistoClass == 6) && (TMath::Abs(particle->GetPdgCode()) != 321)) continue;
+ if ((iHistoClass == 7) && (particle->GetPdgCode()) != 321) continue;
+ if ((iHistoClass == 8) && (particle->GetPdgCode()) != -321) continue;
+ if ((iHistoClass == 9) && (TMath::Abs(particle->GetPdgCode()) != 2212)) continue;
+ if ((iHistoClass == 10) && (particle->GetPdgCode()) != 2212) continue;
+ if ((iHistoClass == 11) && (particle->GetPdgCode()) != -2212) continue;
+
+ // Secondary specification not set.
+ //cout<<"Charge: "<<particle->Charge()<<" PDG code: "<<particle->GetPdgCode()<<endl;
+ //cout<<particle->IsPhysicalPrimary()<<" "<<particle->IsSecondaryFromWeakDecay()<<" "<<particle->IsSecondaryFromMaterial()<<endl;
+
+ // These two functions are not implemented...
+ if (particle->IsSecondaryFromWeakDecay()) cout<<"Secondary From Weak Decay!"<<endl;
+ if (particle->IsSecondaryFromMaterial()) cout<<"Secondary From Material!"<<endl;
+
+ FillGenMCHistos(iHistoClass, particle);
+
+ }
+ }
+
+ // Particle has Passed.
+ return kTRUE;
+
+}
+
+// -------------------------------------------------------------------------
+Bool_t AliAODTrackCutsDiHadronPID::IsSelectedReconstructedMC(AliTrackDiHadronPID* track) {
+
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ if (!track) return kFALSE;
+
+ if (!fPrimRecMCTrackQAHistos||!fSecRecMCTrackQAHistos) {cout<<"AliAODTrackCutsDiHadronPID - Histograms were not created, you should have called CreateHistos()..."<<endl;}
+
+ // Check the track cuts.
+ if (!CheckPt(track->MCPt())) return kFALSE; // Kinematic cuts are done on the MC particles
+ if (!CheckMaxEta(track->MCEta())) return kFALSE; // and not on the data.
+ if (!CheckRapidity(track->MCY())) return kFALSE; // NEW: rapidity cut.
+ if (!CheckFilterMask(track->GetFilterMap())) return kFALSE;
+ if (!CheckFlags(track->GetFlags())) return kFALSE;
+
+ Int_t NSPDhits = 0;
+ if (track->HasPointOnITSLayer(0)) NSPDhits++;
+ if (track->HasPointOnITSLayer(1)) NSPDhits++;
+ if (!CheckPtDeptDCACut(track->GetZAtDCA(),track->GetXYAtDCA(),track->Pt(),NSPDhits)) return kFALSE;
+
+ // Track has passed the cuts, fill QA histograms.
+ for (Int_t iHistoClass = 0; iHistoClass < 12; iHistoClass++) {
+
+ if (fHistRequested[iHistoClass]) {
+
+ // Check the charge (could be neater).
+ if ((iHistoClass == 0) && (track->Charge() == 0)) continue;
+ if ((iHistoClass == 1) && (track->Charge() <= 0)) continue;
+ if ((iHistoClass == 2) && (track->Charge() >= 0)) continue;
+ if ((iHistoClass == 3) && (TMath::Abs(track->GetPdgCode()) != 211)) continue;
+ if ((iHistoClass == 4) && (track->GetPdgCode()) != 211) continue;
+ if ((iHistoClass == 5) && (track->GetPdgCode()) != -211) continue;
+ if ((iHistoClass == 6) && (TMath::Abs(track->GetPdgCode()) != 321)) continue;
+ if ((iHistoClass == 7) && (track->GetPdgCode()) != 321) continue;
+ if ((iHistoClass == 8) && (track->GetPdgCode()) != -321) continue;
+ if ((iHistoClass == 9) && (TMath::Abs(track->GetPdgCode()) != 2212)) continue;
+ if ((iHistoClass == 10) && (track->GetPdgCode()) != 2212) continue;
+ if ((iHistoClass == 11) && (track->GetPdgCode()) != -2212) continue;
+
+ FillRecMCHistos(iHistoClass, track);
+
+ fNTracks[iHistoClass]++;
+ }
+ }
+
+ // Track Has Passed.
+ return kTRUE;
+
+}
+
+// -------------------------------------------------------------------------
+Bool_t AliAODTrackCutsDiHadronPID::FillDataHistos(Int_t histoclass, AliTrackDiHadronPID* track) {
+
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ // Fill the histograms.
+ fHistDataPt[histoclass]->Fill(track->Pt());
+
+ // Fill DCA histograms.
+ fHistDataDCAxy[histoclass]->Fill(track->GetXYAtDCA());
+ fHistDataDCAz[histoclass]->Fill(track->GetZAtDCA());
+
+ // Fill DCA_{xy} one sigma histograms.
+ Int_t checkSum = 0;
+ // histoclass = 0: all charges; histoclass = 1: positive; histoclass = 2: negative;
+ if (TMath::Sqrt(track->GetNumberOfSigmasTOF(0) * track->GetNumberOfSigmasTOF(0) +
+ track->GetNumberOfSigmasTPC(0) * track->GetNumberOfSigmasTPC(0)) < 1.) {
+ fHistDataDCAxyOneSigma[0 + histoclass]->Fill(track->Pt(),track->GetXYAtDCA()); // All species.
+ fHistDataDCAxyOneSigma[3 + histoclass]->Fill(track->Pt(),track->GetXYAtDCA()); // Pions.
+ // checkSum++; cout<<"Pion found: nSigTOF: "<<track->GetNumberOfSigmasTOF(0)<<"; nSigTPC: "<<track->GetNumberOfSigmasTPC(0)<<endl;
+ }
+ if (TMath::Sqrt(track->GetNumberOfSigmasTOF(1) * track->GetNumberOfSigmasTOF(1) +
+ track->GetNumberOfSigmasTPC(1) * track->GetNumberOfSigmasTPC(1)) < 1.) {
+ fHistDataDCAxyOneSigma[0 + histoclass]->Fill(track->Pt(),track->GetXYAtDCA()); // All species.
+ fHistDataDCAxyOneSigma[6 + histoclass]->Fill(track->Pt(),track->GetXYAtDCA()); // Kaons.
+ // checkSum++; cout<<"Kaon found: nSigTOF: "<<track->GetNumberOfSigmasTOF(1)<<"; nSigTPC: "<<track->GetNumberOfSigmasTPC(1)<<endl;
+ }
+ if (TMath::Sqrt(track->GetNumberOfSigmasTOF(2) * track->GetNumberOfSigmasTOF(2) +
+ track->GetNumberOfSigmasTPC(2) * track->GetNumberOfSigmasTPC(2)) < 1.) {
+ fHistDataDCAxyOneSigma[0 + histoclass]->Fill(track->Pt(),track->GetXYAtDCA()); // All species.
+ fHistDataDCAxyOneSigma[9 + histoclass]->Fill(track->Pt(),track->GetXYAtDCA()); // Protons.
+ // checkSum++; cout<<"Proton found: nSigTOF: "<<track->GetNumberOfSigmasTOF(2)<<"; nSigTPC: "<<track->GetNumberOfSigmasTPC(2)<<endl;
+ }
+
+ // check for double identification (or triple..)
+ if(checkSum > 1) {AliError("More than one particle identified for the same track!"); }
+
+ // Fill PID histos.
+ for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
+
+ // Note that a possible Y cut is only done on the PID cuts!
+ if (!CheckRapidity(track->Y(iSpecies))) continue;
+
+ for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
+ fHistDataPID[histoclass][iSpecies][iPtClass]->Fill(track->Pt(),track->GetTOFsignalMinusExpected(iSpecies),track->GetTPCsignalMinusExpected(iSpecies));
+ }
+ }
+
+ return kTRUE;
+
+}
+
+// -------------------------------------------------------------------------
+Bool_t AliAODTrackCutsDiHadronPID::FillTOFMismatchHistos(Int_t histoclass, AliTrackDiHadronPID* track, Double_t randomhittime) {
+
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ // Fill TOF mismatch.
+ for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
+
+ // Note that a possible Y cut is only done on the PID cuts!
+ if (!CheckRapidity(track->Y(iSpecies))) continue;
+
+ for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
+ fHistTOFMismatch[histoclass][iSpecies][iPtClass]->Fill(track->Pt(),randomhittime - track->GetTOFsignalExpected(iSpecies));
+ }
+ }
+
+ return kTRUE;
+
+}
+
+// -------------------------------------------------------------------------
+Bool_t AliAODTrackCutsDiHadronPID::FillGenMCHistos(Int_t histoclass, AliAODMCParticle* particle) {
+
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ // Fill the histograms.
+ if (particle->IsPhysicalPrimary()) {
+ fHistPrimGenMCPt[histoclass]->Fill(particle->Pt());
+ } else {
+ fHistSecGenMCPt[histoclass]->Fill(particle->Pt());
+ }
+
+ return kTRUE;
+
+}
+
+// -------------------------------------------------------------------------
+Bool_t AliAODTrackCutsDiHadronPID::FillRecMCHistos(Int_t histoclass, AliTrackDiHadronPID* track) {
+
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ // Fill the Pt histograms.
+ if (track->IsPhysicalPrimary()) {
+ fHistPrimRecMCPt[histoclass]->Fill(track->Pt());
+ } else {
+ fHistSecRecMCPt[histoclass]->Fill(track->Pt());
+ }
+
+ // Fill the DCA histograms.
+ if (track->IsPhysicalPrimary()) {
+ fHistPrimRecMCDCA[histoclass]->Fill(track->Pt(), track->GetXYAtDCA());
+ }
+ if (track->IsSecondaryFromMaterial()) {
+ fHistSecRecMCDCAMat[histoclass]->Fill(track->Pt(), track->GetXYAtDCA());
+ }
+ if (track->IsSecondaryFromWeakDecay()) {
+ fHistSecRecMCDCAWeak[histoclass]->Fill(track->Pt(), track->GetXYAtDCA());
+ }
+
+ return kTRUE;
+
+}
+
+// -------------------------------------------------------------------------
+Bool_t AliAODTrackCutsDiHadronPID::InitializeDataHistos(Int_t histoclass) {
+
+ cout<<"AliAODTrackCutsDiHadronPID - Creating Data Histograms of Class "<<histoclass<<endl;
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ // Add the Pt spectra.
+ fHistDataPt[histoclass] = InitializePtSpectrum("fHistDataPt",histoclass);
+ fDataTrackQAHistos->Add(fHistDataPt[histoclass]);
+
+ // Add the NTrack histograms.
+ fHistDataNTracks[histoclass] = InitializeNTracksHisto("fHistDataNTracks",histoclass);
+ fDataTrackQAHistos->Add(fHistDataNTracks[histoclass]);
+
+ // Add the DCA histograms.
+ fHistDataDCAxy[histoclass] = InitializeDCAxyHisto("fHistDCAxy",histoclass);
+ fDataTrackQAHistos->Add(fHistDataDCAxy[histoclass]);
+ fHistDataDCAz[histoclass] = InitializeDCAzHisto("fHistDCAz",histoclass);
+ fDataTrackQAHistos->Add(fHistDataDCAz[histoclass]);
+
+ // Add the PID histograms. (FIXME shoudl be able to turn the creation of these histos off.)
+ for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
+ for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
+ fHistDataPID[histoclass][iSpecies][iPtClass] = InitializePIDHisto("fHistDataPID",histoclass,iSpecies,iPtClass);
+ fDataTrackQAHistos->Add(fHistDataPID[histoclass][iSpecies][iPtClass]);
+ fHistTOFMismatch[histoclass][iSpecies][iPtClass] = InitializeTOFMismatchHisto("fHistTOFMismatch",histoclass,iSpecies,iPtClass);
+ fDataTrackQAHistos->Add(fHistTOFMismatch[histoclass][iSpecies][iPtClass]);
+ }
+ }
+
+ return kTRUE;
+
+}
+
+// -------------------------------------------------------------------------
+Bool_t AliAODTrackCutsDiHadronPID::InitializeGenMCHistos(Int_t histoclass) {
+
+ cout<<"AliAODTrackCutsDiHadronPID - Creating Generated MC Histograms of Class "<<histoclass<<endl;
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ // Primary Particles.
+ fHistPrimGenMCPt[histoclass] = InitializePtSpectrum("fHistPrimGenMCPt",histoclass);
+ fPrimGenMCTrackQAHistos->Add(fHistPrimGenMCPt[histoclass]);
+
+ // Secondary Particles.
+ fHistSecGenMCPt[histoclass] = InitializePtSpectrum("fHistSecGenMCPt",histoclass);
+ fSecGenMCTrackQAHistos->Add(fHistSecGenMCPt[histoclass]);
+
+ return kTRUE;
+
+}
+
+// -------------------------------------------------------------------------
+Bool_t AliAODTrackCutsDiHadronPID::InitializeRecMCHistos(Int_t histoclass) {
+
+ cout<<"AliAODTrackCutsDiHadronPID - Creating Reconstructed MC Histograms of Class "<<histoclass<<endl;
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ // Primary Particles.
+ fHistPrimRecMCPt[histoclass] = InitializePtSpectrum("fHistPrimRecMCPt",histoclass);
+ fPrimRecMCTrackQAHistos->Add(fHistPrimRecMCPt[histoclass]);
+
+ fHistPrimRecNTracks[histoclass] = InitializeNTracksHisto("fHistPrimRecNTracks",histoclass);
+ fPrimRecMCTrackQAHistos->Add(fHistPrimRecNTracks[histoclass]);
+
+ fHistPrimRecMCDCA[histoclass] = InitializeDCASpectrum("fHistPrimRecDCA",histoclass);
+ fPrimRecMCTrackQAHistos->Add(fHistPrimRecMCDCA[histoclass]);
+
+ // Secondary Particles.
+ fHistSecRecMCPt[histoclass] = InitializePtSpectrum("fHistSecRecMCPt",histoclass);
+ fSecRecMCTrackQAHistos->Add(fHistSecRecMCPt[histoclass]);
+
+ fHistSecRecMCDCAMat[histoclass] = InitializeDCASpectrum("fHistSecRecDCAMat",histoclass);
+ fSecRecMCTrackQAHistos->Add(fHistSecRecMCDCAMat[histoclass]);
+
+ fHistSecRecMCDCAWeak[histoclass] = InitializeDCASpectrum("fHistSecRecDCAWeak",histoclass);
+ fSecRecMCTrackQAHistos->Add(fHistSecRecMCDCAWeak[histoclass]);
+
+ return kTRUE;
+
+}
+
+// -------------------------------------------------------------------------
+TH1F* AliAODTrackCutsDiHadronPID::InitializePtSpectrum(const char* name, Int_t histoclass) {
+
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ TH1F* hout = new TH1F(Form("%s%s",name,fHistoName[histoclass].Data()),
+ Form("p_{T} Spectrum (%s);p_{T} (GeV/c);Count",fHistoName[histoclass].Data()),fNPtBins,fPtAxis);
+
+ hout->SetDirectory(0);
+
+ return hout;
+
+}
+
+// -------------------------------------------------------------------------
+TH2F* AliAODTrackCutsDiHadronPID::InitializeDCASpectrum(const char* name, Int_t histoclass){
+
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ TH2F* hout = new TH2F(Form("%s%s",name,fHistoName[histoclass].Data()),
+ Form("DCA_{xy} (%s); p_{T} (GeV/c); DCA_{xy} (cm)",fHistoName[histoclass].Data()),
+ fNPtBins,fPtAxis,
+ 300,-3.,3.);
+
+ hout->SetDirectory(0);
+
+ return hout;
+
+}
+
+// -------------------------------------------------------------------------
+TH1F* AliAODTrackCutsDiHadronPID::InitializeNTracksHisto(const char* name, Int_t histoclass) {
+
+ TH1F* hout = new TH1F(Form("%s%s",name,fHistoName[histoclass].Data()),
+ Form("Number of Accepted Tracks (%s);N%s;Count",fHistoName[histoclass].Data(),fHistoLatex[histoclass].Data()),
+ 100,0,4000);
+
+ hout->SetDirectory(0);
+
+ return hout;
+
+}
+
+// -------------------------------------------------------------------------
+TH1F* AliAODTrackCutsDiHadronPID::InitializeDCAxyHisto(const char* name, Int_t histoclass) {
+
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ TH1F* hout = new TH1F(Form("%s%s",name,fHistoName[histoclass].Data()),
+ Form("DCAxy (%s);DCAxy (cm);Count",fHistoName[histoclass].Data()),300,-15.,15.);
+
+ hout->SetDirectory(0);
+
+ return hout;
+
+}
+
+// -------------------------------------------------------------------------
+TH1F* AliAODTrackCutsDiHadronPID::InitializeDCAzHisto(const char* name, Int_t histoclass) {
+
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ TH1F* hout = new TH1F(Form("%s%s",name,fHistoName[histoclass].Data()),
+ Form("DCAz (%s);DCAz (cm);Count",fHistoName[histoclass].Data()),300,-15.,15.);
+
+ hout->SetDirectory(0);
+
+ return hout;
+
+}
+
+// -------------------------------------------------------------------------
+TH3F* AliAODTrackCutsDiHadronPID::InitializeAcceptanceHisto(const char* /*name*/, Int_t /*histoclass*/) {
+
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+ return 0x0;
+
+}
+
+// -------------------------------------------------------------------------
+TH3F* AliAODTrackCutsDiHadronPID::InitializePIDHisto(const char* name, Int_t histoclass, Int_t expspecies, Int_t ptclass) {
+
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ TH3F* hout = new TH3F(Form("%s%s%s%s",name,fHistoName[histoclass].Data(),fParticleName[expspecies].Data(),fPtClassName[ptclass].Data()),
+ Form("PID %s (Exp: %s);p_{T} (GeV/c);#Delta t (ms);dE/dx (a.u.)",fHistoName[histoclass].Data(),fParticleName[expspecies].Data()),
+ fNPtBinsPID[ptclass],fPtBoundaryPID[ptclass],fPtBoundaryPID[ptclass+1],
+ fTOFbins[ptclass][expspecies],fTOFLowerBound[ptclass][expspecies],fTOFUpperBound[ptclass][expspecies],
+ fTPCbins[ptclass][expspecies],fTPCLowerBound[ptclass][expspecies],fTPCUpperBound[ptclass][expspecies]);
+
+ hout->SetDirectory(0);
+
+ return hout;
+
+}
+
+// -------------------------------------------------------------------------
+TH2F* AliAODTrackCutsDiHadronPID::InitializeTOFMismatchHisto(const char* name, Int_t histoclass, Int_t expspecies, Int_t ptclass) {
+
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ TH2F* hout = new TH2F(Form("%s%s%s%s",name,fHistoName[histoclass].Data(),fParticleName[expspecies].Data(),fPtClassName[ptclass].Data()),
+ Form("TOF Mismatch %s (Exp: %s);p_{T} (GeV/c);#Delta t (ms)",fHistoName[histoclass].Data(),fParticleName[expspecies].Data()),
+ fNPtBinsPID[ptclass],fPtBoundaryPID[ptclass],fPtBoundaryPID[ptclass+1],
+ fTOFbins[ptclass][expspecies],fTOFLowerBound[ptclass][expspecies],fTOFUpperBound[ptclass][expspecies]);
+
+ hout->SetDirectory(0);
+
+ return hout;
+
+}
+
+// -------------------------------------------------------------------------
+TH1F* AliAODTrackCutsDiHadronPID::GetHistDataTOFProjection(Int_t charge, Int_t species, Int_t ptbin) {
+
+ // Returns a projection in TOF of the data histogram.
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ Int_t ptclass = GetPtClass(ptbin);
+ if (ptclass == -1) return 0x0;
+ Int_t bininptclass = GetBinInPtClass(ptbin);
+
+ // Retrieve original the 3D histogram.
+ TH3F* htmp = (TH3F*)GetHistData(Form("fHistDataPID%s%s%s",fHistoName[charge].Data(),fParticleName[species].Data(),fPtClassName[ptclass].Data()));
+
+ // Make the projection.
+ TH1F* htmp_proj = (TH1F*)htmp->ProjectionY(Form("TOFprojection_%i_%i_%i",charge,species,ptbin),bininptclass,bininptclass);
+
+ // Some settings of the output histogram.
+ htmp_proj->SetDirectory(0);
+ htmp_proj->SetTitle(Form("%5.3f < p_{T} < %5.3f",GetPtMinPID(ptbin),GetPtMaxPID(ptbin)));
+ htmp_proj->Sumw2();
+
+ return htmp_proj;
+
+}
+
+// -------------------------------------------------------------------------
+TH1F* AliAODTrackCutsDiHadronPID::GetHistDataTOFMismatch(Int_t charge, Int_t species, Int_t ptbin) {
+
+ // Returns a projection in TOF of the mismatch histogram.
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ Int_t ptclass = GetPtClass(ptbin);
+ if (ptclass == -1) return 0x0;
+ Int_t bininptclass = GetBinInPtClass(ptbin);
+
+ // Retrieve the original 2D histogram.
+ TH2F* htmp = (TH2F*)GetHistData(Form("fHistTOFMismatch%s%s%s",fHistoName[charge].Data(),fParticleName[species].Data(),fPtClassName[ptclass].Data()));
+
+ // Make the projection.
+ TH1F* htmp_proj = (TH1F*)htmp->ProjectionY(Form("TOFprojectionsMismatch_%i_%i_%i",charge,species,ptbin),bininptclass,bininptclass);
+
+ // Some settings of the output histogram.
+ htmp_proj->SetDirectory(0);
+ htmp_proj->SetTitle(Form("%5.3f < p_{T} < %5.3f",GetPtMinPID(ptbin),GetPtMaxPID(ptbin)));
+ htmp_proj->Sumw2();
+
+ return htmp_proj;
+
+}
--- /dev/null
+#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
+
+// -------------------------------------------------------------------------
+// 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;}
+ Float_t* GetPtAxisPID() const {
+ const Int_t nptbinspid = GetNPtBinsPID();
+ Float_t* ptaxis = new Float_t[nptbinspid];
+ 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);
+ 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.
+ 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[47]; // Pt axis used in all histograms, except PID and Mismatch histograms.
+ Int_t fNPtBins; // Number of bins in the pt-axis.
+
+ 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
--- /dev/null
+// -------------------------------------------------------------------------
+// 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;
+
+}
--- /dev/null
+#ifndef ALIANALYSYSTASKDIHADRONPID_H
+#define ALIANALYSYSTASKDIHADRONPID_H
+
+#include "AliAnalysisTaskSE.h"
+#include "AliEventPoolManager.h"
+#include "AliAODTrackCutsDiHadronPID.h"
+#include "AliAODEventCutsDiHadronPID.h"
+#include "TObjArray.h"
+#include "THn.h"
+
+class AliAnalysisTaskDiHadronPID : public AliAnalysisTaskSE {
+
+public:
+ // Constructors/Destructors.
+ AliAnalysisTaskDiHadronPID();
+ AliAnalysisTaskDiHadronPID(const char* name);
+ virtual ~AliAnalysisTaskDiHadronPID();
+
+ // Methods from AliAnalysisTaskSE.
+ void UserCreateOutputObjects();
+ void LocalInit();
+ void UserExec(Option_t*);
+ void Terminate(Option_t*);
+
+ // Are all cut objects provided to the task?
+ Bool_t ReadyToStart() const {return (fEventCuts && fTrackCutsTrigger && fTrackCutsAssociated);}
+
+ // Setters.
+ void SetDebugLevel(Int_t debuglevel) {fDebug = debuglevel;}
+ void SetEventCuts(AliAODEventCutsDiHadronPID* eventcuts) {fEventCuts = eventcuts;}
+ void SetTrackCutsTrigger(AliAODTrackCutsDiHadronPID* trackcuts) {fTrackCutsTrigger = trackcuts;}
+ void SetTrackCutsAssociated(AliAODTrackCutsDiHadronPID* trackcuts) {fTrackCutsAssociated = trackcuts;}
+
+ void SetNDEtaBins(Int_t nbins) {fNDEtaBins = nbins;}
+ void SetNDPhiBins(Int_t nbins) {fNDPhiBins = nbins;}
+ void SetMinEventsForMixing(Int_t nevents) {fMinNEventsForMixing = nevents;}
+ void SetPoolTrackDepth(Int_t trackdepth) {fPoolTrackDepth = trackdepth;}
+ void SetPoolSize(Int_t poolsize) {fPoolSize = poolsize;}
+
+ // Getters.
+ Int_t GetNDEtaBins() const {return fNDEtaBins;}
+ Int_t GetNDPhiBins() const {return fNDPhiBins;}
+ Int_t GetMinEventsForMixing() const {return fMinNEventsForMixing;}
+ Int_t GetPoolTrackDepth() const {return fPoolTrackDepth;}
+ Int_t GetPoolSize() const {return fPoolSize;}
+ Int_t GetDebugLevel() const {return fDebug;}
+
+private:
+ //void FillGlobalTracksArray();
+ Bool_t LoadExtMismatchHistos();
+ Double_t GenerateRandomHit(Double_t eta);
+
+private:
+
+ // PID Response Object.
+ AliPIDResponse* fPIDResponse; //! PID Response.
+
+ // Event Cuts Object.
+ AliAODEventCutsDiHadronPID* fEventCuts; //
+
+ // Track Cuts Object.
+ AliAODTrackCutsDiHadronPID* fTrackCutsTrigger; //
+ AliAODTrackCutsDiHadronPID* fTrackCutsAssociated; //
+
+ // Event Pool Manager.
+ AliEventPoolManager* fPoolMgr; //! Event pool manager.
+
+ // Track Arrays.
+ TObjArray* fTriggerTracks; //!
+ TObjArray* fAssociatedTracks; //!
+ // TObjArray* fGlobalTracksArray;
+
+ // Current Event.
+ AliAODEvent* fCurrentAODEvent; //! Current AOD Event.
+
+ // Output List.
+ TList* fOutputList; //! Output List.
+
+ // Histograms.
+ TH1F* fPtSpectrum; //! Pt Spectrum.
+ TH3F* fCorrelations; //! Correlations Histogram.
+ TH3F* fMixedEvents; //! Mixed Events Histogram.
+
+ TObjArray* fTOFhistos; //! Array holding all correlation functions.
+ THnF* fCorrelationsTOF[5][3]; //! Correlations with TOF info.
+ THnF* fCorrelationsTOFTPC[5][3]; //! Correlations with TPC and TOF info.
+
+ // Settings.
+ Int_t fNDEtaBins; //
+ Int_t fNDPhiBins; //
+ Int_t fMinNEventsForMixing; // Pool needs at least this many events for mixing.
+ Int_t fPoolTrackDepth; // For the pool.
+ Int_t fPoolSize; //
+ Bool_t fCalculateTOFmismatch; //
+
+ // TOF mismatch stuff.
+ TH1F* fT0Fill; //
+ TH2F* fLvsEta; //
+ TObjArray* fLvsEtaProjections; //
+
+ // Flags.
+ Int_t fDebug; // Debug flag.
+
+ ClassDef(AliAnalysisTaskDiHadronPID,2);
+
+};
+
+#endif
--- /dev/null
+// -------------------------------------------------------------------------
+// Tools for drawing/ manipulating histograms.
+// ! NEEDS CLEANUP !
+// -------------------------------------------------------------------------
+// Author: Misha Veldhoen (misha.veldhoen@cern.ch)
+
+#include <iostream>
+#include "TCanvas.h"
+#include "TLegend.h"
+#include "TMath.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TH3F.h"
+#include "THn.h"
+#include "TObjArray.h"
+#include "TStyle.h"
+#include "TArray.h"
+
+#include "AliHistToolsDiHadronPID.h"
+
+using namespace std;
+
+// -------------------------------------------------------------------------
+TObjArray* AliHistToolsDiHadronPID::CreateSpectraComparison(const char* name,
+ const char* title, TH1F* h1, TH1F* h2, Int_t markerstyle, Bool_t logy) {
+
+ // - Creates a window comparing two histograms h1, and h2.
+ // - Returns an array of pointers to the objects created
+ // in this function.
+
+ //Int_t optstat = gStyle->GetOptStat();
+ gStyle->SetOptStat(0);
+ //c->UseCurrentStyle();
+ //gStyle->SetOptStat(optstat);
+
+ TH1F* spectra[2];
+ TH1F* division;
+ Int_t spectracolor[2] = {1,2};
+ Int_t divisioncolor = 4;
+
+ // Cloning and configuring spectra.
+ spectra[0] = (TH1F*)h1->Clone();
+ //spectra[0]->SetDirectory(0);
+ spectra[1] = (TH1F*)h2->Clone();
+ //spactra[1]->SetDirectory(0);
+ for (Int_t iSpectra = 0; iSpectra < 2; iSpectra++) {
+ spectra[iSpectra]->Sumw2();
+ spectra[iSpectra]->SetMarkerStyle(markerstyle);
+ spectra[iSpectra]->SetLineColor(spectracolor[iSpectra]);
+ spectra[iSpectra]->SetMarkerColor(spectracolor[iSpectra]);
+ }
+
+ // Creating the division.
+ division = (TH1F*)spectra[0]->Clone("hDivision");
+ division->Divide(spectra[1]);
+ division->SetLineColor(divisioncolor);
+ division->SetMarkerColor(divisioncolor);
+ division->GetYaxis()->SetTitle("Ratio");
+
+ // Creating window
+ //TCanvas* c = new TCanvas(name,title,0,0,400,400);
+ TCanvas* c = TCanvas::MakeDefCanvas();
+ c->SetName(name);
+ c->SetTitle(title);
+ TPad* p1 = new TPad("p1","pad1",0,0.25,1,1);
+ p1->SetRightMargin(0.05);
+ p1->SetBottomMargin(0.);
+ TPad* p2 = new TPad("p2","pad2",0,0,1,0.25);
+ p2->SetTopMargin(0.);
+ p2->SetRightMargin(0.05);
+ p2->SetBottomMargin(0.3);
+ p2->SetFillStyle(4000);
+ p1->Draw();
+ p2->Draw();
+
+ // Determining plotting range.
+ Double_t max = TMath::Max(spectra[0]->GetMaximum(),spectra[1]->GetMaximum());
+ Double_t min = TMath::Min(spectra[0]->GetMinimum(),spectra[1]->GetMinimum());
+ Double_t range = max-min;
+ spectra[0]->SetMinimum(min-0.05*range);
+ spectra[0]->SetMaximum(max+0.05*range);
+
+ // Drawing
+ p1->cd();
+ if (logy) {p1->SetLogy();}
+ spectra[0]->Draw("e");
+ spectra[0]->GetXaxis()->SetLabelSize(0);
+ spectra[0]->GetYaxis()->SetLabelSize(0.05);
+ spectra[0]->GetYaxis()->SetTitleSize(0.05);
+ spectra[0]->GetYaxis()->SetTitleOffset(1.1);
+ //Double_t labelsize = spectra[0]->GetYaxis()->GetLabelSize();
+ spectra[1]->Draw("same e");
+ p2->cd();
+ division->SetTitle("");
+ division->GetXaxis()->SetLabelSize(0.1);
+ division->GetXaxis()->SetTitleSize(0.15);
+ division->GetXaxis()->SetTitleOffset(0.6);
+ division->GetYaxis()->SetLabelSize(0.1);
+ division->GetYaxis()->SetTitleSize(0.12);
+ division->GetYaxis()->SetTitleOffset(0.3);
+ division->Draw("e");
+
+ //c->UseCurrentStyle();
+
+ //gStyle->SetOptStat(optstat);
+
+ TLegend* legend = new TLegend(0.75,0.8,0.95,0.95);
+ legend->AddEntry(spectra[0],h1->GetTitle(),"lpe");
+ legend->AddEntry(spectra[1],h2->GetTitle(),"lpe");
+ p1->cd();
+ legend->Draw("same");
+
+ // returning the created objects.
+ TObjArray* returnobjects = new TObjArray(6);
+ returnobjects->AddLast(c);
+ returnobjects->AddLast(p1);
+ returnobjects->AddLast(p2);
+ returnobjects->AddLast(spectra[0]);
+ returnobjects->AddLast(spectra[1]);
+ returnobjects->AddLast(division);
+
+ return returnobjects;
+
+}
+
+// -------------------------------------------------------------------------
+TH1F* AliHistToolsDiHadronPID::RebinVariableBinning(
+ TH1F* histIn, Double_t* binsx, Int_t Nbinsx, Bool_t density) {
+
+ // Rebins a histogram (hin) with a variable binning to a histogram
+ // with another variable binning (binsx). If the "density" flag is set,
+ // then it is expected that the bins are per unit x-axis, otherwise an
+ // absolute count is assumed.
+ //
+ // TODO: determine over/under-flow bins.
+
+ // Gather info from the original histogram.
+ const TArrayD* binsxIn = histIn->GetXaxis()->GetXbins();
+ TArrayD* binsxOut = new TArrayD(Nbinsx + 1,binsx);
+ Int_t NbinsxIn = binsxIn->GetSize() - 1;
+ const char* nameIn = histIn->GetName();
+ const char* titleIn = histIn->GetTitle();
+ const char* xAxisTitleIn = histIn->GetXaxis()->GetTitle();
+ const char* yAxisTitleIn = histIn->GetYaxis()->GetTitle();
+
+ // Gather info from the output histogram and create it.
+ Int_t NbinsxOut = binsxOut->GetSize() - 1;
+ const char* nameOut = Form("%s (rebinned)",nameIn);
+ const char* titleOut = Form("%s_rebinned;%s;%s",titleIn,xAxisTitleIn,yAxisTitleIn);
+ TH1F* histOut = new TH1F(nameOut,titleOut,NbinsxOut,binsxOut->GetArray());
+ //histOut->Sumw2();
+
+ // Filling the bins.
+ for (Int_t iBinOut = 0; iBinOut < NbinsxOut; iBinOut++) {
+
+ Double_t binOutLowEdge = binsxOut->At(iBinOut);
+ Double_t binOutHighEdge = binsxOut->At(iBinOut+1);
+ Double_t binOutWidth = binOutHighEdge - binOutLowEdge;
+ //Double_t binOutCenter = binOutHighEdge - binOutWidth/2.;
+
+ // Setting all errors of the outgoing histogram to zero (needed later):
+ histOut->SetBinError(iBinOut+1,0.);
+
+ for (Int_t iBinIn = 0; iBinIn < NbinsxIn; iBinIn++) {
+
+ Double_t binInLowEdge = binsxIn->At(iBinIn);
+ Double_t binInHighEdge = binsxIn->At(iBinIn+1);
+ Double_t binInWidth = binInHighEdge - binInLowEdge;
+ //Double_t binInCenter = binInHighEdge - binInWidth/2.;
+ Double_t binInContent = histIn->GetBinContent(iBinIn+1);
+ Double_t binInError2 = histIn->GetBinError(iBinIn+1)*histIn->GetBinError(iBinIn+1);
+
+ /* -- Determining what percentage of the in-bin is included
+ in the current out-bin -- */
+ Double_t PercentageIncluded = 0.;
+
+ // Position of the low edge of the in-bin:
+ // In 1 2 3
+ // Out | |
+ Int_t binInLowEdgePos = 0;
+ if (binInLowEdge < binOutLowEdge) binInLowEdgePos = 1;
+ else if (binInLowEdge > binOutHighEdge) binInLowEdgePos = 3;
+ else binInLowEdgePos = 2;
+
+ // Same for the high edge of the in-bin:
+ Int_t binInHighEdgePos = 0;
+ if (binInHighEdge < binOutLowEdge) binInHighEdgePos = 1;
+ else if (binInHighEdge > binOutHighEdge) binInHighEdgePos = 3;
+ else binInHighEdgePos = 2;
+
+ // In-bin not included.
+ // In | | or | |
+ // Out | | | |
+ if ( binInLowEdgePos == 3 || binInHighEdgePos == 1 ) {
+ //cout<<"In-bin not included."<<endl;
+ PercentageIncluded = 0.;
+ }
+
+ // In-bin partially included (left side).
+ // In | |
+ // Out | |
+ else if ( binInLowEdgePos == 2 && binInHighEdgePos == 3 ) {
+ //cout<<"In-bin partially included (left side)."<<endl;
+ PercentageIncluded = (binOutHighEdge - binInLowEdge)/(binInHighEdge - binInLowEdge);
+ }
+
+ // In-bin partially included (middle).
+ // In | |
+ // Out | |
+ else if ( binInLowEdgePos == 1 && binInHighEdgePos == 3 ) {
+ //cout<<"In-bin partially included (middle)."<<endl;
+ PercentageIncluded = (binOutHighEdge - binOutLowEdge)/(binInHighEdge - binInLowEdge);
+ }
+
+ // In-bin partially included (right side).
+ // In | |
+ // Out | |
+ else if ( binInLowEdgePos == 1 && binInHighEdgePos == 2 ) {
+ //cout<<"In-bin partially included (right side)."<<endl;
+ PercentageIncluded = (binInHighEdge - binOutLowEdge)/(binInHighEdge - binInLowEdge);
+ }
+
+ // In-bin completely included.
+ // In | |
+ // Out | |
+ else if ( binInLowEdgePos == 2 && binInHighEdgePos == 2 ) {
+ //cout<<"In-bin completely included."<<endl;
+ PercentageIncluded = 1.;
+ }
+
+ // Give some status:
+ //cout<<"Bin Out: "<<iBinOut+1<<" ("<<binOutLowEdge<<","<<binOutHighEdge<<") Bin In: "<<iBinIn+1<<" ("<<binInLowEdge<<","<<binInHighEdge<<"), content: "<<binInContent<<" +/- "<<TMath::Sqrt(binInError2)<<". Bin in Pos: "<<binInLowEdgePos<<" "<<binInHighEdgePos<<" pct: "<<PercentageIncluded<<endl;
+
+ //else cout<<"Something weird just happened."<<endl;
+
+ Double_t binOutAddedError2 = 0.;
+ binOutAddedError2 = binInError2 * PercentageIncluded * PercentageIncluded;
+ Double_t binOutAddedContent = 0.;
+ binOutAddedContent = binInContent * PercentageIncluded;
+
+ if (density) {
+ binOutAddedContent*=binInWidth;
+ }
+
+ //histOut->Fill(binOutCenter,binOutAddedContent);
+ histOut->SetBinContent(iBinOut+1,histOut->GetBinContent(iBinOut+1)+binOutAddedContent);
+
+ //cout<<histOut->GetBinError(iBinOut+1)<<endl;
+ //cout<<binOutAddedError2<<endl;
+ histOut->SetBinError(iBinOut+1,histOut->GetBinError(iBinOut+1)+binOutAddedError2);
+ //cout<<histOut->GetBinError(iBinOut+1)<<endl;
+ }
+
+ // Once a certain binOut has been filled all the way we have to scale
+ // it if the "density" flag is set.
+
+ histOut->SetBinContent(iBinOut + 1,histOut->GetBinContent(iBinOut+1)/binOutWidth);
+ histOut->SetBinError(iBinOut + 1,TMath::Sqrt(histOut->GetBinError(iBinOut+1)));
+ //cout<<histOut->GetBinError(iBinOut+1)<<endl;
+ }
+
+ return histOut;
+
+}
+
+// -------------------------------------------------------------------------
+TH1F* AliHistToolsDiHadronPID::TrimHisto(TH1F* histo, Int_t firstbin, Int_t lastbin) {
+
+ const char* name = histo->GetName();
+ const char* title = histo->GetTitle();
+ if (firstbin < 0) return 0x0;
+ if (lastbin > histo->GetNbinsX()) return 0x0;
+ if (firstbin > lastbin) return 0x0;
+
+ Int_t Nbins = lastbin - firstbin + 1;
+ cout<<firstbin<<" "<<lastbin<<" "<<Nbins<<endl;
+ //const Int_t NbinsC = Nbins;
+ Double_t newaxis[41];
+
+ for (Int_t ii = firstbin; ii < lastbin+1; ii++) {
+ newaxis[ii - firstbin] = histo->GetBinLowEdge(ii);
+ //cout<<ii-firstbin<<" "<<newaxis[ii - firstbin]<<endl;
+ }
+ newaxis[Nbins] = histo->GetBinLowEdge(lastbin+1);
+
+ TH1F* hout = new TH1F("hout","hout",Nbins,newaxis);
+ hout->SetDirectory(0);
+ hout->SetTitle(title);
+ hout->SetName(name);
+
+ for (Int_t ii = 0; ii < Nbins+1; ii++) {
+ hout->SetBinContent(ii,histo->GetBinContent(firstbin+ii));
+ hout->SetBinError(ii,histo->GetBinError(firstbin+ii));
+ }
+
+ return hout;
+
+}
--- /dev/null
+#ifndef ALIHISTTOOLSDIHADRONPID_H
+#define ALIHISTTOOLSDIHADRONPID_H
+
+class AliHistToolsDiHadronPID {
+
+public:
+ AliHistToolsDiHadronPID() {};
+
+protected:
+ ~AliHistToolsDiHadronPID() {};
+
+public:
+
+ // Histogram Manipulation.
+ static TH1F* RebinVariableBinning(TH1F* histIn, Double_t* binsx, Int_t Nbinsx, Bool_t density = kTRUE);
+ static TH1F* RebinVariableBinning(TH1F* histIn, TH1F* histAxis, Bool_t density = kTRUE) {
+
+ // Rebins histogram histIn to the x-axis of histAxis
+ TAxis* xaxis = histAxis->GetXaxis();
+ Int_t nbinsx = xaxis->GetNbins();
+ const Double_t* binsx = (xaxis->GetXbins())->GetArray();
+ return RebinVariableBinning(histIn, const_cast<Double_t*>(binsx), nbinsx, density);
+
+ }
+ static TH1F* TrimHisto(TH1F* histo, Int_t firstbin, Int_t lastbin);
+ static void ConstMinusHist(TH1F* histo, const Float_t cc = 1) {
+
+ // h -> (c-h)
+ Int_t nbins = histo->GetNbinsX();
+ for (Int_t iBin = 0; iBin < (nbins + 1); iBin++) {
+ Float_t bincontent = histo->GetBinContent(iBin);
+ histo->SetBinContent(iBin,(cc - bincontent));
+ }
+
+ }
+ static TH3F* MakeHist3D(const char* name, const char* title,
+ Int_t nbinsX, Double_t minX, Double_t maxX,
+ Int_t nbinsY, Double_t minY, Double_t maxY,
+ Int_t nbinsZ, const Float_t* zaxis) {
+
+ const Float_t* xaxis = const_cast<Float_t*>(CreateAxis(nbinsX,minX,maxX));
+ const Float_t* yaxis = const_cast<Float_t*>(CreateAxis(nbinsY,minY,maxY));
+
+ TH3F* hout = new TH3F(name,title,nbinsX,xaxis,nbinsY,yaxis,nbinsZ,zaxis);
+
+ return hout;
+ }
+
+ // Histogram Visualization.
+ static TObjArray* CreateSpectraComparison(const char* name, const char* title, TH1F* h1, TH1F* h2, Int_t markerstyle = 8, Bool_t logy = kTRUE);
+
+private:
+ static Float_t* CreateAxis(const Int_t nbins, Double_t min, Double_t max) {
+ if (nbins <= 0) return 0x0;
+ if (max < min) return 0x0;
+
+ Float_t* axis = new Float_t[nbins + 1];
+ Float_t binsize = (max - min)/((Float_t)nbins);
+ for (Int_t iBin = 0; iBin < nbins + 1; iBin++) {
+ axis[iBin] = min + ((Float_t)iBin) * binsize;
+ }
+
+ return axis;
+ }
+
+
+};
+
+#endif
--- /dev/null
+// -------------------------------------------------------------------------
+// Copies all info that is needed for the DiHadronPID analysis.
+// Possible Extension: At this moment the object is protected for returning
+// pointers to the original tracks. It could at some point be beneficial to
+// be able to access this information.
+// -------------------------------------------------------------------------
+// Author: Misha Veldhoen (misha.veldhoen@cern.ch)
+
+#include <iostream>
+using namespace std;
+
+// AOD includes.
+#include "AliAODTrack.h"
+#include "AliAODEvent.h"
+#include "AliAODVertex.h"
+#include "AliAODMCParticle.h"
+
+// PID includes.
+#include "AliPID.h"
+#include "AliPIDResponse.h"
+#include "AliTPCPIDResponse.h"
+
+// Objects own include.
+#include "AliTrackDiHadronPID.h"
+
+ClassImp(AliTrackDiHadronPID);
+
+// -------------------------------------------------------------------------
+AliTrackDiHadronPID::AliTrackDiHadronPID():
+ TObject(),
+ fAODTrack(0x0),
+ fAODGlobalTrack(0x0),
+ fAODEvent(0x0),
+ fAODMCParticle(0x0),
+ fPIDResponse(0x0),
+ fBasicInfoAvailable(kFALSE),
+ fFlagsAvailable(kFALSE),
+ fDCAInfoAvailable(kFALSE),
+ fITSInfoAvailable(kFALSE),
+ fTPCInfoAvailable(kFALSE),
+ fTOFInfoAvailable(kFALSE),
+ fMCInfoAvailable(kFALSE),
+ fPt(-999.),
+ fEta(-999.),
+ fPhi(-999.),
+ fFlags(0),
+ fFilterMap(0),
+ fID(0),
+ fLabel(0),
+ fCharge(0),
+ fDCAz(-999.),
+ fDCAxy(-999.),
+ fTOFsignal(-999.),
+ fIsTOFmismatch(kFALSE),
+ fTPCsignal(-999.),
+ fTPCmomentum(-999.),
+ fMCPt(-999.),
+ fMCEta(-999.),
+ fMCPhi(-999.),
+ fMCY(-999.),
+ fPdgCode(0),
+ fIsPhysicalPrimary(kFALSE),
+ fIsSecondaryFromWeakDecay(kFALSE),
+ fIsSecondaryFromMaterial(kFALSE),
+ fDebug(0)
+
+{
+
+ //
+ // Default Constructor.
+ //
+
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
+ fTOFsignalMinusExpected[iSpecies] = -999.;
+ fTOFNsigma[iSpecies] = -999.;
+ fTPCsignalMinusExpected[iSpecies] = -999.;
+ fTPCNsigma[iSpecies] = -999.;
+ fY[iSpecies] = -999.;
+ }
+
+ for (Int_t iITSlayer = 0; iITSlayer < 6; iITSlayer++) {
+ fITSHits[iITSlayer] = kFALSE;
+ }
+
+}
+
+// -------------------------------------------------------------------------
+AliTrackDiHadronPID::AliTrackDiHadronPID(AliAODTrack* track, AliAODTrack* globaltrack, AliAODMCParticle* mcparticle, AliPIDResponse* pidresponse):
+ TObject(),
+ fAODTrack(0x0),
+ fAODGlobalTrack(0x0),
+ fAODEvent(0x0),
+ fAODMCParticle(0x0),
+ fPIDResponse(0x0),
+ fBasicInfoAvailable(kFALSE),
+ fFlagsAvailable(kFALSE),
+ fDCAInfoAvailable(kFALSE),
+ fITSInfoAvailable(kFALSE),
+ fTPCInfoAvailable(kFALSE),
+ fTOFInfoAvailable(kFALSE),
+ fMCInfoAvailable(kFALSE),
+ fPt(-999.),
+ fEta(-999.),
+ fPhi(-999.),
+ fFlags(0),
+ fFilterMap(0),
+ fID(0),
+ fLabel(0),
+ fCharge(0),
+ fDCAz(-999.),
+ fDCAxy(-999.),
+ fTOFsignal(-999.),
+ fIsTOFmismatch(kFALSE),
+ fTPCsignal(-999.),
+ fTPCmomentum(-999.),
+ fMCPt(-999.),
+ fMCEta(-999.),
+ fMCPhi(-999.),
+ fMCY(-999.),
+ fPdgCode(0),
+ fIsPhysicalPrimary(kFALSE),
+ fIsSecondaryFromWeakDecay(kFALSE),
+ fIsSecondaryFromMaterial(kFALSE),
+ fDebug(0)
+{
+
+ //
+ // Constructor.
+ //
+
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
+ fTOFsignalMinusExpected[iSpecies] = -999.;
+ fTOFNsigma[iSpecies] = -999.;
+ fTPCsignalMinusExpected[iSpecies] = -999.;
+ fTPCNsigma[iSpecies] = -999.;
+ fY[iSpecies] = -999.;
+ }
+
+ for (Int_t iITSlayer = 0; iITSlayer < 6; iITSlayer++) {
+ fITSHits[iITSlayer] = kFALSE;
+ }
+
+ if (track) {
+ fAODTrack = track;
+ fAODEvent = track->GetAODEvent();
+ }
+ if (globaltrack) fAODGlobalTrack = globaltrack;
+ if (mcparticle) fAODMCParticle = mcparticle;
+ if (pidresponse) fPIDResponse = pidresponse;
+
+ // Copy AOD Track info.
+ if (fAODTrack) {
+ CopyBasicTrackInfo();
+ } else {
+ AliError("No Track Supplied.");
+ }
+
+ // Find the Global Track.
+ if (fID >= 0) fAODGlobalTrack = fAODTrack;
+
+ // Copy DCA and PID info.
+ if (fAODGlobalTrack) {
+ CopyFlags();
+ if (fAODEvent) CopyDCAInfo();
+ else AliError("Couln't find AOD Event.");
+ CopyITSInfo();
+ if (fPIDResponse) CopyTPCInfo();
+ CopyTOFInfo();
+ } else {
+ AliError("Couldn't find Global Track.");
+ }
+
+ // Copy MC info.
+ if (fAODMCParticle) {
+ CopyMCInfo();
+ }
+
+ // Test
+ /* Double_t sigmaTOFProton = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(fAODTrack, AliPID::kProton));
+ if ( sigmaTOFProton < 1.0) {cout<<"tofsigmabelowone: "<<sigmaTOFProton<<endl;}
+
+ Double_t sigmaTPCProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(fAODTrack, AliPID::kProton));
+ if ( sigmaTPCProton < 1.0) {cout<<"tpcsigmabelowone: "<<sigmaTPCProton<<endl;}*/
+}
+
+// -------------------------------------------------------------------------
+Bool_t AliTrackDiHadronPID::CopyBasicTrackInfo() {
+
+ //
+ // Copies everything available in every AOD track.
+ //
+
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ fPt = fAODTrack->Pt();
+ fEta = fAODTrack->Eta();
+ fPhi = fAODTrack->Phi();
+
+ fY[0] = fAODTrack->Y(AliAODTrack::kPion);
+ fY[1] = fAODTrack->Y(AliAODTrack::kKaon);
+ fY[2] = fAODTrack->Y(AliAODTrack::kProton);
+
+ //fFlags = fAODTrack->GetFlags(); // FLAGS MUST BE TAKEN FROM GLOBAL TRACKS.
+ fFilterMap = fAODTrack->GetFilterMap();
+
+ fID = fAODTrack->GetID();
+ fLabel = fAODTrack->GetLabel();
+
+ fCharge = fAODTrack->Charge();
+
+ fBasicInfoAvailable = kTRUE;
+ return fBasicInfoAvailable;
+
+}
+
+// -------------------------------------------------------------------------
+Bool_t AliTrackDiHadronPID::CopyFlags() {
+
+ //
+ // Copies Flags (properly stored in global track)
+ //
+
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ // Copy Flags
+ fFlags = fAODGlobalTrack->GetFlags();
+
+ // Is TOF mismatch?
+ if (AliAODTrack::kTOFmismatch&fFlags) {
+ fIsTOFmismatch = kTRUE;
+ //cout<<"Found TOF mismatch!"<<endl;
+ }
+ else fIsTOFmismatch = kFALSE;
+
+ fFlagsAvailable = kTRUE;
+ return fFlagsAvailable;
+
+}
+
+// -------------------------------------------------------------------------
+Bool_t AliTrackDiHadronPID::CopyDCAInfo() {
+
+ //
+ // Copies DCA info. (only stored in a global track)
+ //
+
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ // Propagate track to DCA.
+ Double_t PosAtDCA[2] = {-999,-999};
+ Double_t covar[3] = {-999,-999,-999};
+ //cout<<fAODTrack<<" "<<fAODGlobalTrack<<endl;
+ Bool_t propagate = fAODGlobalTrack->PropagateToDCA(fAODEvent->GetPrimaryVertex(),fAODEvent->GetMagneticField(),100.,PosAtDCA,covar);
+
+ if (propagate) {
+ fDCAxy = PosAtDCA[0];
+ fDCAz = PosAtDCA[1];
+ } else {
+ //AliError("Could not propagate track to DCA.");
+ }
+
+ if (propagate) fDCAInfoAvailable = kTRUE;
+ return fDCAInfoAvailable;
+
+}
+
+// -------------------------------------------------------------------------
+Bool_t AliTrackDiHadronPID::CopyITSInfo() {
+
+ //
+ // Copies ITS info.
+ //
+
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ // Get the ITS clustermap
+ fITSClusterMap = fAODGlobalTrack->GetITSClusterMap();
+
+ // Copy the ITS hits.
+ for (Int_t iITSlayer = 0; iITSlayer < 6; iITSlayer++) {
+ fITSHits[iITSlayer] = fAODGlobalTrack->HasPointOnITSLayer(iITSlayer);
+ }
+
+ fITSInfoAvailable = kTRUE;
+ return fITSInfoAvailable;
+
+}
+
+// -------------------------------------------------------------------------
+Bool_t AliTrackDiHadronPID::CopyTPCInfo() {
+
+ //
+ // Copies TPC info. (needs global track and pid response)
+ //
+
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ // Get TPC signal.
+ fTPCsignal = fAODGlobalTrack->GetTPCsignal();
+
+ // Compute expected TPC signal under pi/K/p mass assumption.
+ AliTPCPIDResponse& TPCPIDResponse = fPIDResponse->GetTPCResponse();
+ fTPCmomentum = fAODGlobalTrack->GetTPCmomentum();
+
+ fTPCsignalMinusExpected[0] = fTPCsignal - TPCPIDResponse.GetExpectedSignal(fTPCmomentum,AliPID::kPion);
+ fTPCsignalMinusExpected[1] = fTPCsignal - TPCPIDResponse.GetExpectedSignal(fTPCmomentum,AliPID::kKaon);
+ fTPCsignalMinusExpected[2] = fTPCsignal - TPCPIDResponse.GetExpectedSignal(fTPCmomentum,AliPID::kProton);
+
+ fTPCNsigma[0] = fPIDResponse->NumberOfSigmasTPC(fAODGlobalTrack, AliPID::kPion);
+ fTPCNsigma[1] = fPIDResponse->NumberOfSigmasTPC(fAODGlobalTrack, AliPID::kKaon);
+ fTPCNsigma[2] = fPIDResponse->NumberOfSigmasTPC(fAODGlobalTrack, AliPID::kProton);
+
+ fTPCInfoAvailable = kTRUE;
+ return fTPCInfoAvailable;
+
+}
+
+// -------------------------------------------------------------------------
+Bool_t AliTrackDiHadronPID::CopyTOFInfo() {
+
+ //
+ // Copies TOF info. (needs global track)
+ //
+
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ // Get TOF signal.
+ fTOFsignal = fAODGlobalTrack->GetTOFsignal();
+
+ // Compute expected TOF signal under pi/K/p mass assumption.
+ Double_t times[AliPID::kSPECIES];
+ fAODGlobalTrack->GetIntegratedTimes(times);
+ fTOFsignalMinusExpected[0] = fTOFsignal - times[AliPID::kPion];
+ fTOFsignalMinusExpected[1] = fTOFsignal - times[AliPID::kKaon];
+ fTOFsignalMinusExpected[2] = fTOFsignal - times[AliPID::kProton];
+
+ fTOFNsigma[0] = fPIDResponse->NumberOfSigmasTOF(fAODGlobalTrack, AliPID::kPion);
+ fTOFNsigma[1] = fPIDResponse->NumberOfSigmasTOF(fAODGlobalTrack, AliPID::kKaon);
+ fTOFNsigma[2] = fPIDResponse->NumberOfSigmasTOF(fAODGlobalTrack, AliPID::kProton);
+
+ fTOFInfoAvailable = kTRUE;
+ return fTOFInfoAvailable;
+
+}
+
+// -------------------------------------------------------------------------
+Bool_t AliTrackDiHadronPID::CopyMCInfo() {
+
+ // Copies MC info (needs an MC track with the same label)
+
+ // Check if the label of the current track matches the label of the
+ // generated particle. Note that the label of the AOD track can be
+ // negative. This means that the quality of this track is not awesome,
+ // but that it does correspond to the MC particle.
+
+ if (fDebug) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
+
+ /*
+ if (fAODMCParticle->Label() != TMath::Abs(fAODTrack->GetLabel())) {
+ cout<<"Label of supplied MC particle and reconstructed track do not match."<<endl;
+ return kFALSE;
+ }
+ */
+ // Note: It seems like the Label of the AOD track points to the INDEX of the
+ // MCPArticle, not to the label (See AliAnalysisTaskCompareAODTrackCuts.cxx)
+
+ fMCPt = fAODMCParticle->Pt();
+ fMCEta = fAODMCParticle->Eta();
+ fMCPhi = fAODMCParticle->Phi();
+ fMCY = fAODMCParticle->Y();
+ fPdgCode = fAODMCParticle->PdgCode();
+
+ TClonesArray* mcArray = 0x0;
+ mcArray = dynamic_cast<TClonesArray*>(fAODEvent->FindListObject(AliAODMCParticle::StdBranchName()));
+ if (!mcArray) {
+ AliFatal("No MC array found in the AOD.");
+ }
+
+ // Primary particle
+ if ( fAODMCParticle->IsPhysicalPrimary() ){
+ fIsPhysicalPrimary = kTRUE;
+ } else {
+ // Safety check for mother existence.
+ if (fAODMCParticle->GetMother() >= 0){
+
+ Int_t mcMotherPDG = -999;
+ Int_t firstInt = -999;
+
+ AliAODMCParticle* mcMother = (AliAODMCParticle*) mcArray->At(TMath::Abs(fAODMCParticle->GetMother()));
+ mcMotherPDG = TMath::Abs(mcMother->GetPdgCode());
+
+ // Need a way to get the first intiger, for now Marek's method:
+ firstInt = Int_t (mcMotherPDG/ TMath::Power(10, Int_t(TMath::Log10(mcMotherPDG))));
+ // cout<<"Mother PDG: "<<mcMotherPDG<<"; Firt integer: "<<firstInt<<endl;
+
+ // Weak decay
+ if( firstInt == 3){
+ fIsSecondaryFromWeakDecay = kTRUE;
+ // Material decay
+ } else {
+ fIsSecondaryFromMaterial = kTRUE;
+ }
+ }
+ }
+
+ fMCInfoAvailable = kTRUE;
+ return fMCInfoAvailable;
+
+}
--- /dev/null
+#ifndef ALITRACKDIHADRONPID_H
+#define ALITRACKDIHADRONPID_H
+
+#include <iostream>
+using namespace std;
+
+#include "AliAODTrack.h"
+#include "AliAODEvent.h"
+#include "AliAODMCParticle.h"
+#include "AliPIDResponse.h"
+
+class AliTrackDiHadronPID : public TObject {
+
+public:
+ AliTrackDiHadronPID();
+ AliTrackDiHadronPID(AliAODTrack* track, AliAODTrack* globaltrack = 0x0, AliAODMCParticle* mcparticle = 0x0, AliPIDResponse* pidresponse = 0x0);
+ //virtual ~AliTrackDiHadronPID();
+
+// Internal copy functions.
+private:
+ Bool_t CopyBasicTrackInfo();
+ Bool_t CopyFlags();
+ Bool_t CopyDCAInfo();
+ Bool_t CopyITSInfo();
+ Bool_t CopyTPCInfo();
+ Bool_t CopyTOFInfo();
+ Bool_t CopyMCInfo();
+
+// Check functions.
+public:
+ Bool_t IsBasicTrackInfoAvailable() const {
+ if (!fBasicInfoAvailable) cout<<"Basic Track Info not available."<<endl;
+ return fBasicInfoAvailable;
+ }
+ Bool_t IsFlagInfoAvailable() const {
+ if (!fFlagsAvailable) cout<<"Flag Info not available."<<endl;
+ return fFlagsAvailable;
+ }
+ Bool_t IsDCAInfoAvailable() const {
+ if (!fDCAInfoAvailable) cout<<"DCA Info not available."<<endl;
+ return fDCAInfoAvailable;
+ }
+ Bool_t IsITSInfoAvailable() const {
+ if (!fITSInfoAvailable) cout<<"ITS Info not available."<<endl;
+ return fITSInfoAvailable;
+ }
+ Bool_t IsTPCInfoAvailable() const {
+ if (!fTPCInfoAvailable) cout<<"TPC Info not available."<<endl;
+ return fTPCInfoAvailable;
+ }
+ Bool_t IsTOFInfoAvailable() const {
+ if (!fTOFInfoAvailable) cout<<"TOF Info not available."<<endl;
+ return fTOFInfoAvailable;
+ }
+ Bool_t IsMCInfoAvailable() const {
+ if (!fMCInfoAvailable) cout<<"MC Info not available."<<endl;
+ return fMCInfoAvailable;
+ }
+
+public:
+// Getting Track Parameters. Functionality is the same as AOD track,
+// unless stated otherwise.
+ Double_t Pt() const {return fPt;}
+ Double_t Eta() const {return fEta;}
+ Double_t Phi() const {return fPhi;}
+ Double_t Y(Int_t species) {
+ if (species >= 0 && species < 3) return fY[species];
+ else return -999.;
+ }
+
+ ULong_t GetFlags() const {return fFlags;}
+ ULong_t GetStatus() const {return GetFlags();}
+ UInt_t GetFilterMap() const {return fFilterMap;}
+ Bool_t TestFilterMask(UInt_t filterMask) const {return (Bool_t)((filterMask & fFilterMap) == filterMask);}
+
+ Int_t GetID() const {return fID;}
+ Int_t GetLabel() const {return fLabel;}
+ Short_t Charge() const {return fCharge;}
+
+ Double_t GetZAtDCA() const {return fDCAz;}
+ Double_t GetXYAtDCA() const {return fDCAxy;}
+
+ Double_t GetTOFsignal() const {return fTOFsignal;}
+ Double_t GetTOFsignalMinusExpected(Int_t species) const {
+ if (species < 0 || species > 3) {
+ cout<<"ERROR: Unknown species"<<endl;
+ return -10e10;
+ }
+ return fTOFsignalMinusExpected[species];
+ }
+ Double_t GetTOFsignalExpected(Int_t species) const {
+ if (species < 0 || species > 3) {
+ cout<<"ERROR: Unknown species"<<endl;
+ return -10e10;
+ }
+ return (fTOFsignal - fTOFsignalMinusExpected[species]);
+ }
+ Double_t GetNumberOfSigmasTOF(Int_t species) const {
+ if (species < 0 || species > 3) {
+ cout<<"ERROR: Unknown species"<<endl;
+ return -10e10;
+ }
+ return fTOFNsigma[species];
+ }
+ Double_t GetNumberOfSigmasTPC(Int_t species) const {
+ if (species < 0 || species > 3) {
+ cout<<"ERROR: Unknown species"<<endl;
+ return -10e10;
+ }
+ return fTPCNsigma[species];
+ }
+ Bool_t IsTOFmismatch() const {return fIsTOFmismatch;}
+
+ Double_t GetTPCsignal() const {return fTPCsignal;}
+ Double_t GetTPCsignalMinusExpected(Int_t species) const {return fTPCsignalMinusExpected[species];}
+ Double_t GetTPCmomentum() const {return fTPCmomentum;}
+
+ Char_t GetITSClusterMap() const {return fITSClusterMap;}
+ Bool_t HasPointOnITSLayer(Int_t layer) const {
+ if ((layer > -1) && (layer < 6)) return fITSHits[layer];
+ else {
+ cout<<"ITS has only 6 layers."<<endl;
+ return kFALSE;
+ }
+ }
+
+ Double_t MCPt() const {return fMCPt;}
+ Double_t MCEta() const {return fMCEta;}
+ Double_t MCPhi() const {return fMCPhi;}
+ Double_t MCY() const {return fMCY;}
+ Int_t GetMCSpecies() const {
+ Int_t abspdg = TMath::Abs(GetPdgCode());
+ if (abspdg == 211) return 0;
+ else if (abspdg == 321) return 1;
+ else if (abspdg == 2212) return 2;
+ else return -999;
+ }
+
+ Int_t GetPdgCode() const {return fPdgCode;}
+ Bool_t IsPhysicalPrimary() const {return fIsPhysicalPrimary;}
+ Bool_t IsSecondaryFromWeakDecay() const {return fIsSecondaryFromWeakDecay;}
+ Bool_t IsSecondaryFromMaterial() const {return fIsSecondaryFromMaterial;}
+
+ void ForgetAboutPointers() {
+ // AOD tracks are usually deleted, while the
+ // AliTrackDiHadronPID is not. This method ensures that
+ // the pointers to the track/event objects etc. don't
+ // point to deleted objects.
+ fAODTrack = 0x0;
+ fAODGlobalTrack = 0x0;
+ fAODEvent = 0x0;
+ fAODMCParticle = 0x0;
+ fPIDResponse = 0x0;
+ }
+
+ void SetDebugLevel(Int_t debuglevel) {fDebug = debuglevel;}
+ Int_t GetDebugLevel() const {return fDebug;}
+
+private:
+// Pointers to original Tracks etc (cannot be returned, will not be streamed/saved).
+ AliAODTrack* fAODTrack; //! Original AOD Track.
+ AliAODTrack* fAODGlobalTrack; //! Corresponding Global AOD Track.
+ const AliAODEvent* fAODEvent; //! Original AOD Event.
+ AliAODMCParticle* fAODMCParticle; //! Original MC Particle.
+ AliPIDResponse* fPIDResponse; //! Original PID Response.
+
+// Flags if a certain type of information is available for this track.
+ Bool_t fBasicInfoAvailable;// Basic Track Info.
+ Bool_t fFlagsAvailable;
+ Bool_t fDCAInfoAvailable; // DCA Info.
+ Bool_t fITSInfoAvailable; // ITS Info.
+ Bool_t fTPCInfoAvailable; // TPC Info.
+ Bool_t fTOFInfoAvailable; // TOF Info.
+ Bool_t fMCInfoAvailable; // MC Info.
+
+// Basic Track Info.
+ Double_t fPt; // Reconstructed Pt.
+ Double_t fEta; // Reconstructed Eta.
+ Double_t fY[3]; // Reconstructed Rapidity for pi,K,p.
+ Double_t fPhi; // Reconstructed Phi.
+
+ ULong_t fFlags; // Reconstruction Flags.
+ UInt_t fFilterMap; // FilterMap.
+
+ Short_t fID; // Track ID.
+ Int_t fLabel; // Track Label.
+
+ Short_t fCharge; // Charge (is a Char_t in AliAODTrack)
+
+// DCA Info.
+ Double_t fDCAz; // z at DCA.
+ Double_t fDCAxy; // xy at DCA.
+
+// PID Info.
+ Double_t fTOFsignal;
+ Double_t fTOFsignalMinusExpected[3];
+ Double_t fTOFNsigma[3];
+ Bool_t fIsTOFmismatch;
+ Double_t fTPCsignal;
+ Double_t fTPCsignalMinusExpected[3];
+ Double_t fTPCNsigma[3];
+ Double_t fTPCmomentum;
+ UChar_t fITSClusterMap;
+ Bool_t fITSHits[6];
+
+// MC Info.
+ Double_t fMCPt;
+ Double_t fMCEta;
+ Double_t fMCPhi;
+ Double_t fMCY;
+ Int_t fPdgCode;
+ Bool_t fIsPhysicalPrimary;
+ Bool_t fIsSecondaryFromWeakDecay;
+ Bool_t fIsSecondaryFromMaterial;
+
+// Debug.
+ Int_t fDebug; // Debug flag.
+
+ ClassDef(AliTrackDiHadronPID,1);
+
+};
+
+#endif
+++ /dev/null
-// AddTask Macro (v 8.00).
-// Updated: May 2nd 2012.
-// Author: Misha Veldhoen (m.veldhoen@cern.ch)
-
-AliAnalysisTaskDiHadronPID *AddTaskDiHadronPID(Int_t verbose = 1,
- Bool_t printbuffersize = kFALSE,
- Bool_t mixedevents = kTRUE,
- Bool_t setmc = kFALSE,
- TString beamtype = "PbPb",
- Double_t MaxEta = 0.8,
- Double_t maxrap = 0.5,
- Double_t MaxPlotEta = 0.8,
- Double_t MaxPt = 10.,
- Int_t NEtaBins = 32,
- Int_t NPhiBins = 36,
- Double_t VertexZMixedEvents = 2.,
- Bool_t zoomed = kFALSE,
- Bool_t DoITSCut = kFALSE,
- Bool_t DoDCACut = kTRUE,
- Bool_t DemandNoMismatch = kTRUE,
- Int_t trigbuffermaxsize=2500,
- Double_t centralitycutmax=0.,
- Double_t centralitycutmin=10.,
- const char* outputFileName = 0,
- const char* containerName = "DiHadronPID",
- const char* folderName = "PWGCF_DiHadronPID")
-
-
-
-{
-
- // Get the current analysis manager.
- AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
- if (!mgr) {
- Error("AddTaskDiHadronPID.C", "No analysis manager found.");
- return 0;
- }
-
- // Create an instance of the task.
- AliAnalysisTaskDiHadronPID *task = new AliAnalysisTaskDiHadronPID(containerName);
-
- // Configure the task.
- task->SetVerbose(verbose);
- task->SetPrintBufferSize(printbuffersize);
- task->SetCalculateMixedEvents(mixedevents);
- task->SetMC(setmc);
- task->SetBeamType(beamtype);
- task->SetMaxEta(MaxEta);
- task->SetMaxRapidityInInclusiveSpectra(maxrap);
- task->SetMaxPlotEta(MaxPlotEta);
- task->SetMaxPt(MaxPt);
- task->SetNEtaBins(NEtaBins);
- task->SetNPhiBins(NPhiBins);
- task->SetVertexZMixedEvents(VertexZMixedEvents);
- task->SetZoomed(zoomed);
- task->SetDoITSCut(DoITSCut);
- task->SetDoDCACut(DoDCACut);
- task->SetDemandNoMismatch(DemandNoMismatch);
- task->SetTrigBufferMaxSize(trigbuffermaxsize);
- task->SetCentralityCut(centralitycutmax,centralitycutmin);
-
- // Add the task.
- mgr->AddTask(task);
-
- // Data containers.
- AliAnalysisDataContainer *cinput = mgr->GetCommonInputContainer();
- mgr->ConnectInput(task, 0, cinput);
-
- if (!outputFileName)
- outputFileName = AliAnalysisManager::GetCommonFileName();
-
- AliAnalysisDataContainer *coutput1 =
- mgr->CreateContainer(containerName, TList::Class(),
- AliAnalysisManager::kOutputContainer,Form("%s:%s", outputFileName, folderName));
-
- mgr->ConnectOutput (task, 1, coutput1);
-
- return task;
-
-}
-
-
--- /dev/null
+AliAnalysisTaskDiHadronPID* AddTaskDiHadronPID(\r
+ Int_t NDEtaBins = 32,\r
+ Int_t NDPhiBins = 32,\r
+ Int_t MinEventsForMixing = 5,\r
+ Double_t MinCentrality = 5.,\r
+ Double_t MaxCentrality = 0.,\r
+ const char* CentralityEstimator = "V0M",\r
+ Double_t maxVertexZ = 7.,\r
+ Double_t maxEta = 0.8,\r
+ Double_t minAssociatedPt = 0.2,\r
+ Double_t maxAssociatedPt = 5.0,\r
+ Double_t minTriggerPt = 5.,\r
+ Double_t maxTriggerPt = 10.,\r
+ Bool_t requestAllSingleTrackHistos = kTRUE,\r
+ Int_t FilterMaskTrigger = 7,\r
+ Int_t FilterMaskAssociated = 5,\r
+ Bool_t isPbPb = kTRUE,\r
+ Bool_t isMC = kFALSE,\r
+ Int_t DebugLevel = 0,\r
+ const char* outputFileName = 0,\r
+ const char* containerName = "DiHadronPID",\r
+ const char* folderName = "PWGCF_DiHadronPID")\r
+\r
+{\r
+ // Get a pointer to the analysis manager.\r
+ AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();\r
+ if (!mgr) {\r
+ cout<<"AddTaskDiHadronPID.C -> No analysis manager found."<<endl;\r
+ return 0x0;\r
+ } \r
+\r
+ // Create an instance of the task.\r
+ AliAnalysisTaskDiHadronPID* DiHadronPIDTask = new AliAnalysisTaskDiHadronPID(containerName);\r
+\r
+ // Configure the task.\r
+ DiHadronPIDTask->SetNDEtaBins(NDEtaBins);\r
+ DiHadronPIDTask->SetNDPhiBins(NDPhiBins);\r
+ DiHadronPIDTask->SetMinEventsForMixing(MinEventsForMixing);\r
+ DiHadronPIDTask->SetDebugLevel(DebugLevel);\r
+\r
+ // Configure and add Event Cuts.\r
+ AliAODEventCutsDiHadronPID* eventcuts = new AliAODEventCutsDiHadronPID("EventCuts");\r
+ eventcuts->SetTrigger(AliVEvent::kMB);\r
+ eventcuts->SetCentrality(MaxCentrality, MinCentrality);\r
+ eventcuts->SetMaxVertexZ(maxVertexZ);\r
+ eventcuts->SetCentralityEstimator(CentralityEstimator);\r
+ eventcuts->SetIsPbPb(isPbPb);\r
+ eventcuts->SetDebugLevel(DebugLevel);\r
+ DiHadronPIDTask->SetEventCuts(eventcuts);\r
+\r
+ // Configure and add track cuts for trigger.\r
+ AliAODTrackCutsDiHadronPID* triggercuts = new AliAODTrackCutsDiHadronPID("TrackCutsTrigger");\r
+ triggercuts->SetIsMC(isMC);\r
+ triggercuts->SetFilterMask(1<<FilterMaskTrigger);\r
+ triggercuts->SetPtRange(minTriggerPt,maxTriggerPt);\r
+ triggercuts->SetMaxEta(maxEta);\r
+ triggercuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kAllCharged);\r
+ triggercuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kPositive);\r
+ triggercuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kNegative);\r
+ if (requestAllSingleTrackHistos) {\r
+ triggercuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kAllPion);\r
+ triggercuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kPosPion);\r
+ triggercuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kNegPion); \r
+ triggercuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kAllKaon);\r
+ triggercuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kPosKaon);\r
+ triggercuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kNegKaon); \r
+ triggercuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kAllProton);\r
+ triggercuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kPosProton);\r
+ triggercuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kNegProton);\r
+ } triggercuts->SetDebugLevel(DebugLevel);\r
+ DiHadronPIDTask->SetTrackCutsTrigger(triggercuts);\r
+\r
+ // Configure and add track cuts for associateds.\r
+ AliAODTrackCutsDiHadronPID* associatedscuts = new AliAODTrackCutsDiHadronPID("TrackCutsAssociated");\r
+ associatedscuts->SetIsMC(isMC);\r
+ associatedscuts->SetFilterMask(1<<FilterMaskAssociated);\r
+ associatedscuts->SetPtRange(minAssociatedPt,maxAssociatedPt);\r
+ associatedscuts->SetMaxEta(maxEta);\r
+ ULong_t associatedflags = (UInt_t)(AliAODTrack::kTOFout)|(UInt_t)(AliAODTrack::kTIME); \r
+ associatedscuts->SetDemandFlags(associatedflags);\r
+ associatedscuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kAllCharged);\r
+ associatedscuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kPositive);\r
+ associatedscuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kNegative);\r
+ if (requestAllSingleTrackHistos) {\r
+ associatedscuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kAllPion);\r
+ associatedscuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kPosPion);\r
+ associatedscuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kNegPion); \r
+ associatedscuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kAllKaon);\r
+ associatedscuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kPosKaon);\r
+ associatedscuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kNegKaon); \r
+ associatedscuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kAllProton);\r
+ associatedscuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kPosProton);\r
+ associatedscuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kNegProton);\r
+ }\r
+ associatedscuts->SetDebugLevel(DebugLevel);\r
+ DiHadronPIDTask->SetTrackCutsAssociated(associatedscuts);\r
+\r
+ // Add the task.\r
+ mgr->AddTask(DiHadronPIDTask);\r
+ \r
+ // Data containers.\r
+ AliAnalysisDataContainer* cinput = mgr->GetCommonInputContainer();\r
+ mgr->ConnectInput(DiHadronPIDTask, 0, cinput); \r
+ \r
+ if (!outputFileName) {outputFileName = AliAnalysisManager::GetCommonFileName();}\r
+ \r
+ AliAnalysisDataContainer* coutput1 = mgr->CreateContainer(containerName, TList::Class(),\r
+ AliAnalysisManager::kOutputContainer,Form("%s:%s", outputFileName, folderName));\r
+ \r
+ mgr->ConnectOutput(DiHadronPIDTask,1,coutput1);\r
+ \r
+ return DiHadronPIDTask;\r
+\r
+}\r
#pragma link C++ class AliAnalysisTaskLeadingTrackUE+;
#pragma link C++ class AliAnalysisTaskMinijet+;
#pragma link C++ class AliAnalysisTaskDiHadron+;
+#pragma link C++ class AliTrackDiHadronPID+;
+#pragma link C++ class AliAODTrackCutsDiHadronPID+;
+#pragma link C++ class AliAODEventCutsDiHadronPID+;
+#pragma link C++ class AliHistToolsDiHadronPID+;
#pragma link C++ class AliAnalysisTaskDiHadronPID+;
#pragma link C++ class AliEvtPoolManager+;
#pragma link C++ class AliEvtPool+;
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff