From: miweber Date: Mon, 8 Apr 2013 11:34:40 +0000 (+0000) Subject: DiHadronPID task update (Misha.Veldhoen@cern.ch) X-Git-Url: http://git.uio.no/git/?a=commitdiff_plain;h=6788af99495abdee5e24c4a7018c1267309be286;p=u%2Fmrichter%2FAliRoot.git DiHadronPID task update (Misha.Veldhoen@cern.ch) --- diff --git a/PWGCF/CMakelibPWGCFCorrelationsDPhi.pkg b/PWGCF/CMakelibPWGCFCorrelationsDPhi.pkg index 47c0101a665..640fb9bf239 100644 --- a/PWGCF/CMakelibPWGCFCorrelationsDPhi.pkg +++ b/PWGCF/CMakelibPWGCFCorrelationsDPhi.pkg @@ -34,7 +34,11 @@ set ( SRCS 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 diff --git a/PWGCF/Correlations/DPhi/AliAnalysisTaskDiHadronPID.cxx b/PWGCF/Correlations/DPhi/AliAnalysisTaskDiHadronPID.cxx deleted file mode 100644 index e31c03096fa..00000000000 --- a/PWGCF/Correlations/DPhi/AliAnalysisTaskDiHadronPID.cxx +++ /dev/null @@ -1,1360 +0,0 @@ -// ---------------------------------------------------------------------------- -// 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 -#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<0) cout<<"AliAnalysisTaskDiHadronPID::UserCreateOutputObjects -> ERROR: Analysis manager not found."<1) { - cout<<"AliAnalysisTaskDiHadronPID::UserCreateOutputObjects -> Analysis manager found."< (manager->GetInputEventHandler()); - if (!inputHandler) { - if (fVerbose>0) cout<<"AliAnalysisTaskDiHadronPID::UserCreateOutputObjects -> ERROR: Input handler not found."<1) { - cout<<"AliAnalysisTaskDiHadronPID::UserCreateOutputObjects -> Input handler found."<GetPIDResponse(); - if (!fPIDResponse) { - if (fVerbose>0) cout<<"AliAnalysisTaskDiHadronPID::UserCreateOutputObjects -> ERROR: PID response object not found."<1) { - cout<<"AliAnalysisTaskDiHadronPID::UserCreateOutputObjects -> PID response object found."<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."<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)&&(ptTestFilterMask(1<<7))) { - //if (fVerbose>3) cout<<"Track Ignored: Did Not pass filterbit."<3) { - cout<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."<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: "<Fill(SPDHits+0.5); - - // 4) ITS cut. - if (fDoITSCut) { - if (!SPDHits) { - if (fVerbose>3) cout<<"Track Ignored: Not enough SPD hits."<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."<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."<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."<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."<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||cent2) cout<<"AliAnalysisTaskDiHadronPID::SelectEvent -> Event did not pass centaltity cut."<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: "<GetID()<<" Partner ID: "<<(-track->GetID()-1)<At(-track->GetID()-1)); - - if (!partner&&(fVerbose>3)) cout<<"AliAnalysisTaskDiHadronPID::GetGlobalTrack -> No Global Track Found!"< 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(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."<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(fAODEvent->FindListObject(AliAODMCParticle::StdBranchName())); - if (!fMCTracks) { - if (fVerbose>0) { - cout<<"AliAnalysisTaskDiHadronPID::UserExec -> ERROR: No MC particles found."<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; iTrackGetEntriesFast(); 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)Fill(mcPt,mcY); - } else { - fPtEtaDistrMCSec[mcSpeciesBin]->Fill(mcPt,mcEta); - if (TMath::Abs(mcY)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: "<Eta()<<" Pt: "<Pt()<<" Ypion: "<Y(AliAODTrack::kPion)<<" Ykaon: "<Y(AliAODTrack::kKaon)<<" Yproton: "<Y(AliAODTrack::kProton)<3) cout<<"Track added to associated buffer."<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."<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<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)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)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)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: "<Eta()<<" Ypion: "<Y(AliAODTrack::kPion)<<" Ykaon: "<Y(AliAODTrack::kKaon)<<" Yproton: "<Y(AliAODTrack::kProton)<IsPhysicalPrimary()) { - fPtEtaDistrDataPrim[mcSpeciesBin]->Fill(dataPt,eta); - if (TMath::Abs(dataY)Fill(dataPt,dataY); - } else { - fPtEtaDistrDataSec[mcSpeciesBin]->Fill(dataPt,eta); - if (TMath::Abs(dataY)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: "<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 "<3) cout << "AliAnalysisTaskDiHadronPID::UserExec -> Buffer size: "<GetZ())3) cout<<"AliAnalysisTaskDiHadronPID::UserExec -> Mixing with trigger Z: "<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&&PionFillTOFPionMinTOF) { - fMixedEventsTPCTOFCut[0]->Fill(DPhi,DEta,pt); - //cout<<"Yes! Histogram will be filled."<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 (KaonFillTPCKaonMinTPC&&KaonFillTOFKaonMinTOF) { - 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 (ProtonFillTPCProtonMinTPC&&ProtonFillTOFProtonMinTOF) { - 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 "<GetEntriesFast()<<" triggers with vertex z = "<GetZ()<<" to the buffer."<GetEntriesFast(); iTrig++) { - - currentTrigger = (AliAODTrack*)(triggers->At(iTrig)); - if (fVerbose>3) cout<<"AliAnalysisTaskDiHadronPID::UserExec -> Trigger pt = "<Pt()<GetZ(); - fTrigBuffer[fTrigBufferIndex][1] = currentTrigger->Phi(); - fTrigBuffer[fTrigBufferIndex][2] = currentTrigger->Eta(); - fTrigBuffer[fTrigBufferIndex][3] = currentTrigger->Pt(); - fTrigBufferIndex++; - if (fTrigBufferSize Trigger buffer index: "< -#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"<0.9) { - cout<<"SetMaxEta -> |eta| must be < 0.9"<1.0) { - cout<<"SetMaxPlotEta -> |eta| must be < 1.0"< Maximum pT must be > 5.0 GeV/c."<72) { - cout<<"SetNEtaBins -> Number of bins must be between 1 and 72"<72) { - cout<<"SetNPhiBins -> Number of bins must be between 1 and 72"<10.) { - cout<<"SetVertexZMixedEvents -> must be 0 < z < 10"<25000) { - cout<<"SetTrigBufferMaxSize -> Max buffer size must be between 10 and 25000."< Centrality cannot be lower than 0."< Maximum centrality needs to be smaller than the minimum centrality. (It's confusing I know)"<100.) { - cout<<"SetCentralityCut -> Minimum centrality cannot exceed 100%."< 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 - diff --git a/PWGCF/Correlations/DPhi/DiHadronPID/AliAODEventCutsDiHadronPID.cxx b/PWGCF/Correlations/DPhi/DiHadronPID/AliAODEventCutsDiHadronPID.cxx new file mode 100644 index 00000000000..68a7047d915 --- /dev/null +++ b/PWGCF/Correlations/DPhi/DiHadronPID/AliAODEventCutsDiHadronPID.cxx @@ -0,0 +1,361 @@ +// ------------------------------------------------------------------------- +// Object managing event cuts, and holding QA histograms. +// ------------------------------------------------------------------------- +// Author: Misha Veldhoen (misha.veldhoen@cern.ch) + +#include +#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."<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."<IsEmpty()) return 1; + + if (!fSelectedEventQAHistos||!fAllEventQAHistos) { + cout<<"AliAODEventCutsDiHadronPID::Merge() - Warning, current object's histograms are missing... Generating."<MakeIterator(); + TObject* obj; + + // List of collections + TList collection_fSelectedEventQAHistos; + TList collection_fAllEventQAHistos; + + Int_t count = 0; + + while ((obj = iter->Next())) { + AliAODEventCutsDiHadronPID* entry = dynamic_cast (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"<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()..."<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: "<GetCentralityPercentile(fCentralityEstimator.Data())<<" Quality: "<GetQuality()<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: "< 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 diff --git a/PWGCF/Correlations/DPhi/DiHadronPID/AliAODTrackCutsDiHadronPID.cxx b/PWGCF/Correlations/DPhi/DiHadronPID/AliAODTrackCutsDiHadronPID.cxx new file mode 100644 index 00000000000..7ab5d2a0e0f --- /dev/null +++ b/PWGCF/Correlations/DPhi/DiHadronPID/AliAODTrackCutsDiHadronPID.cxx @@ -0,0 +1,1018 @@ +// ------------------------------------------------------------------------- +// 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 + +#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."<IsEmpty()) { + //cout<<"List is empty..."<GetEntries(); + //cout<<"Supplied list has "<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 (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: "< 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: "<SetName("PrimGenMCTrackQAHistos"); + fPrimGenMCTrackQAHistos->SetOwner(kTRUE); + } else {cout<<"AliAODTrackCutsDiHadronPID - Warning, Prim. Gen. MC Track QA TList was already created..."<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..."<SetName("PrimRecMCTrackQAHistos"); + fPrimRecMCTrackQAHistos->SetOwner(kTRUE); + } else {cout<<"AliAODTrackCutsDiHadronPID - Warning, Prim. Rec. MC Track QA TList was already created..."<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..."<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..."<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()..."<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()..."<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: "<Charge()<<" PDG code: "<GetPdgCode()<IsPhysicalPrimary()<<" "<IsSecondaryFromWeakDecay()<<" "<IsSecondaryFromMaterial()<IsSecondaryFromWeakDecay()) cout<<"Secondary From Weak Decay!"<IsSecondaryFromMaterial()) cout<<"Secondary From Material!"<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: "<GetNumberOfSigmasTOF(0)<<"; nSigTPC: "<GetNumberOfSigmasTPC(0)<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: "<GetNumberOfSigmasTOF(1)<<"; nSigTPC: "<GetNumberOfSigmasTPC(1)<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: "<GetNumberOfSigmasTOF(2)<<"; nSigTPC: "<GetNumberOfSigmasTPC(2)< 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 "<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 "<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 "<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; + +} diff --git a/PWGCF/Correlations/DPhi/DiHadronPID/AliAODTrackCutsDiHadronPID.h b/PWGCF/Correlations/DPhi/DiHadronPID/AliAODTrackCutsDiHadronPID.h new file mode 100644 index 00000000000..29ef6dba834 --- /dev/null +++ b/PWGCF/Correlations/DPhi/DiHadronPID/AliAODTrackCutsDiHadronPID.h @@ -0,0 +1,424 @@ +#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..."< fNPtBins + 1)) {cout<<"Bin is out of range..."< -1) && (histoclass < 12)) { + fHistRequested[histoclass] = kTRUE; + f3DSpectraEnabeled[histoclass] = Enable3DSpectra; + fPIDHistosEnabeled[histoclass] = EnablePIDHistos; + //cout<<"histoclass: "< ptbin out of range!"< 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 diff --git a/PWGCF/Correlations/DPhi/DiHadronPID/AliAnalysisTaskDiHadronPID.cxx b/PWGCF/Correlations/DPhi/DiHadronPID/AliAnalysisTaskDiHadronPID.cxx new file mode 100644 index 00000000000..60d313ababa --- /dev/null +++ b/PWGCF/Correlations/DPhi/DiHadronPID/AliAnalysisTaskDiHadronPID.cxx @@ -0,0 +1,543 @@ +// ------------------------------------------------------------------------- +// INFO +// ------------------------------------------------------------------------- + +#include + +// 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 (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(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: "<GetEntriesFast()<<" Associateds: "<GetEntriesFast()<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 ..."<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: "<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; + +} diff --git a/PWGCF/Correlations/DPhi/DiHadronPID/AliAnalysisTaskDiHadronPID.h b/PWGCF/Correlations/DPhi/DiHadronPID/AliAnalysisTaskDiHadronPID.h new file mode 100644 index 00000000000..5df43c15adf --- /dev/null +++ b/PWGCF/Correlations/DPhi/DiHadronPID/AliAnalysisTaskDiHadronPID.h @@ -0,0 +1,108 @@ +#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 diff --git a/PWGCF/Correlations/DPhi/DiHadronPID/AliHistToolsDiHadronPID.cxx b/PWGCF/Correlations/DPhi/DiHadronPID/AliHistToolsDiHadronPID.cxx new file mode 100644 index 00000000000..a3114e63074 --- /dev/null +++ b/PWGCF/Correlations/DPhi/DiHadronPID/AliHistToolsDiHadronPID.cxx @@ -0,0 +1,298 @@ +// ------------------------------------------------------------------------- +// Tools for drawing/ manipulating histograms. +// ! NEEDS CLEANUP ! +// ------------------------------------------------------------------------- +// Author: Misha Veldhoen (misha.veldhoen@cern.ch) + +#include +#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."<Fill(binOutCenter,binOutAddedContent); + histOut->SetBinContent(iBinOut+1,histOut->GetBinContent(iBinOut+1)+binOutAddedContent); + + //cout<GetBinError(iBinOut+1)<SetBinError(iBinOut+1,histOut->GetBinError(iBinOut+1)+binOutAddedError2); + //cout<GetBinError(iBinOut+1)<SetBinContent(iBinOut + 1,histOut->GetBinContent(iBinOut+1)/binOutWidth); + histOut->SetBinError(iBinOut + 1,TMath::Sqrt(histOut->GetBinError(iBinOut+1))); + //cout<GetBinError(iBinOut+1)<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<GetBinLowEdge(ii); + //cout<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; + +} diff --git a/PWGCF/Correlations/DPhi/DiHadronPID/AliHistToolsDiHadronPID.h b/PWGCF/Correlations/DPhi/DiHadronPID/AliHistToolsDiHadronPID.h new file mode 100644 index 00000000000..cd2667c44ef --- /dev/null +++ b/PWGCF/Correlations/DPhi/DiHadronPID/AliHistToolsDiHadronPID.h @@ -0,0 +1,69 @@ +#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(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(CreateAxis(nbinsX,minX,maxX)); + const Float_t* yaxis = const_cast(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 diff --git a/PWGCF/Correlations/DPhi/DiHadronPID/AliTrackDiHadronPID.cxx b/PWGCF/Correlations/DPhi/DiHadronPID/AliTrackDiHadronPID.cxx new file mode 100644 index 00000000000..dc66a7443a3 --- /dev/null +++ b/PWGCF/Correlations/DPhi/DiHadronPID/AliTrackDiHadronPID.cxx @@ -0,0 +1,413 @@ +// ------------------------------------------------------------------------- +// 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 +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: "<NumberOfSigmasTPC(fAODTrack, AliPID::kProton)); + if ( sigmaTPCProton < 1.0) {cout<<"tpcsigmabelowone: "<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!"<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."<Pt(); + fMCEta = fAODMCParticle->Eta(); + fMCPhi = fAODMCParticle->Phi(); + fMCY = fAODMCParticle->Y(); + fPdgCode = fAODMCParticle->PdgCode(); + + TClonesArray* mcArray = 0x0; + mcArray = dynamic_cast(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: "< +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."<= 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"< 3) { + cout<<"ERROR: Unknown species"< 3) { + cout<<"ERROR: Unknown species"< 3) { + cout<<"ERROR: Unknown species"< -1) && (layer < 6)) return fITSHits[layer]; + else { + cout<<"ITS has only 6 layers."<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; - -} - - diff --git a/PWGCF/Correlations/macros/DiHadronPID/AddTaskDiHadronPID.C b/PWGCF/Correlations/macros/DiHadronPID/AddTaskDiHadronPID.C new file mode 100644 index 00000000000..a2a6cc22498 --- /dev/null +++ b/PWGCF/Correlations/macros/DiHadronPID/AddTaskDiHadronPID.C @@ -0,0 +1,114 @@ +AliAnalysisTaskDiHadronPID* AddTaskDiHadronPID( + Int_t NDEtaBins = 32, + Int_t NDPhiBins = 32, + Int_t MinEventsForMixing = 5, + Double_t MinCentrality = 5., + Double_t MaxCentrality = 0., + const char* CentralityEstimator = "V0M", + Double_t maxVertexZ = 7., + Double_t maxEta = 0.8, + Double_t minAssociatedPt = 0.2, + Double_t maxAssociatedPt = 5.0, + Double_t minTriggerPt = 5., + Double_t maxTriggerPt = 10., + Bool_t requestAllSingleTrackHistos = kTRUE, + Int_t FilterMaskTrigger = 7, + Int_t FilterMaskAssociated = 5, + Bool_t isPbPb = kTRUE, + Bool_t isMC = kFALSE, + Int_t DebugLevel = 0, + const char* outputFileName = 0, + const char* containerName = "DiHadronPID", + const char* folderName = "PWGCF_DiHadronPID") + +{ + // Get a pointer to the analysis manager. + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if (!mgr) { + cout<<"AddTaskDiHadronPID.C -> No analysis manager found."<SetNDEtaBins(NDEtaBins); + DiHadronPIDTask->SetNDPhiBins(NDPhiBins); + DiHadronPIDTask->SetMinEventsForMixing(MinEventsForMixing); + DiHadronPIDTask->SetDebugLevel(DebugLevel); + + // Configure and add Event Cuts. + AliAODEventCutsDiHadronPID* eventcuts = new AliAODEventCutsDiHadronPID("EventCuts"); + eventcuts->SetTrigger(AliVEvent::kMB); + eventcuts->SetCentrality(MaxCentrality, MinCentrality); + eventcuts->SetMaxVertexZ(maxVertexZ); + eventcuts->SetCentralityEstimator(CentralityEstimator); + eventcuts->SetIsPbPb(isPbPb); + eventcuts->SetDebugLevel(DebugLevel); + DiHadronPIDTask->SetEventCuts(eventcuts); + + // Configure and add track cuts for trigger. + AliAODTrackCutsDiHadronPID* triggercuts = new AliAODTrackCutsDiHadronPID("TrackCutsTrigger"); + triggercuts->SetIsMC(isMC); + triggercuts->SetFilterMask(1<SetPtRange(minTriggerPt,maxTriggerPt); + triggercuts->SetMaxEta(maxEta); + triggercuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kAllCharged); + triggercuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kPositive); + triggercuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kNegative); + if (requestAllSingleTrackHistos) { + triggercuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kAllPion); + triggercuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kPosPion); + triggercuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kNegPion); + triggercuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kAllKaon); + triggercuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kPosKaon); + triggercuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kNegKaon); + triggercuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kAllProton); + triggercuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kPosProton); + triggercuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kNegProton); + } triggercuts->SetDebugLevel(DebugLevel); + DiHadronPIDTask->SetTrackCutsTrigger(triggercuts); + + // Configure and add track cuts for associateds. + AliAODTrackCutsDiHadronPID* associatedscuts = new AliAODTrackCutsDiHadronPID("TrackCutsAssociated"); + associatedscuts->SetIsMC(isMC); + associatedscuts->SetFilterMask(1<SetPtRange(minAssociatedPt,maxAssociatedPt); + associatedscuts->SetMaxEta(maxEta); + ULong_t associatedflags = (UInt_t)(AliAODTrack::kTOFout)|(UInt_t)(AliAODTrack::kTIME); + associatedscuts->SetDemandFlags(associatedflags); + associatedscuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kAllCharged); + associatedscuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kPositive); + associatedscuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kNegative); + if (requestAllSingleTrackHistos) { + associatedscuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kAllPion); + associatedscuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kPosPion); + associatedscuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kNegPion); + associatedscuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kAllKaon); + associatedscuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kPosKaon); + associatedscuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kNegKaon); + associatedscuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kAllProton); + associatedscuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kPosProton); + associatedscuts->RequestQAHistos(AliAODTrackCutsDiHadronPID::kNegProton); + } + associatedscuts->SetDebugLevel(DebugLevel); + DiHadronPIDTask->SetTrackCutsAssociated(associatedscuts); + + // Add the task. + mgr->AddTask(DiHadronPIDTask); + + // Data containers. + AliAnalysisDataContainer* cinput = mgr->GetCommonInputContainer(); + mgr->ConnectInput(DiHadronPIDTask, 0, cinput); + + if (!outputFileName) {outputFileName = AliAnalysisManager::GetCommonFileName();} + + AliAnalysisDataContainer* coutput1 = mgr->CreateContainer(containerName, TList::Class(), + AliAnalysisManager::kOutputContainer,Form("%s:%s", outputFileName, folderName)); + + mgr->ConnectOutput(DiHadronPIDTask,1,coutput1); + + return DiHadronPIDTask; + +} diff --git a/PWGCF/PWGCFCorrelationsDPhiLinkDef.h b/PWGCF/PWGCFCorrelationsDPhiLinkDef.h index 446a8274242..e69b8b96053 100644 --- a/PWGCF/PWGCFCorrelationsDPhiLinkDef.h +++ b/PWGCF/PWGCFCorrelationsDPhiLinkDef.h @@ -12,6 +12,10 @@ #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+;