+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+// AliAnalysisTwoParticleResonanceFlowTask:
+// author: Redmer Alexander Bertens (rbertens@cern.ch)
+// You Zhou (you.zhou@cern.ch)
+// analyis task for
+// AliResonanceFlowHelperTrack provides a lightweight helper track for reconstruction
+// This analysis DOES NOT support ESDs
+// Version: v2012.9.7
+
+
+#include "TChain.h"
+#include "TTree.h"
+#include "TH1.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TCanvas.h"
+#include "TMath.h"
+#include "TObjArray.h"
+#include "AliAnalysisTask.h"
+#include "AliAnalysisTaskSE.h"
+#include "AliAnalysisManager.h"
+#include "AliAODEvent.h"
+#include "AliAODInputHandler.h"
+#include "AliCentrality.h"
+#include "AliVEvent.h"
+#include "AliAnalysisTwoParticleResonanceFlowTask.h"
+#include "AliFlowBayesianPID.h"
+#include "AliStack.h"
+#include "AliMCEvent.h"
+#include "TProfile.h"
+#include "AliFlowCandidateTrack.h"
+#include "AliFlowTrackCuts.h"
+#include "AliFlowEventSimple.h"
+#include "AliFlowTrackSimple.h"
+#include "AliFlowCommonConstants.h"
+#include "AliFlowEvent.h"
+#include "TVector3.h"
+#include "AliAODVZERO.h"
+#include "AliPID.h"
+#include "AliPIDResponse.h"
+#include "AliAODMCParticle.h"
+#include "AliAnalysisTaskVnV0.h"
+#include "AliEventPoolManager.h"
+
+class AliFlowTrackCuts;
+
+using std::cout;
+using std::endl;
+
+ClassImp(AliAnalysisTwoParticleResonanceFlowTask)
+
+AliAnalysisTwoParticleResonanceFlowTask::AliAnalysisTwoParticleResonanceFlowTask() : AliAnalysisTaskSE(),
+ fSpeciesA(0), fSpeciesB(0), fChargeA(0), fChargeB(0), fMassA(0), fMassB(0), fDebug(0), fIsMC(0), fEventMixing(0), fPhiMinusPsiMethod(0), fQA(0), fV0(0), fMassBins(1), fMinMass(-1.), fMaxMass(0.), fCutsRP(NULL), fNullCuts(0), fPIDResponse(0), fFlowEvent(0), fBayesianResponse(0), fCandidates(0), fCandidateEtaPtCut(0), fCandidateMinEta(0), fCandidateMaxEta(0), fCandidateMinPt(0), fCandidateMaxPt(0), fNPtBins(18), fNdPhiBins(18), fCentrality(999), fVertex(999), fAOD(0), fPoolManager(0), fOutputList(0), fEventStats(0), fCentralityPass(0), fCentralityNoPass(0), fNOPID(0), fPIDk(0), fPIDp(0), fNOPIDTOF(0), fPIDTOF(0), fPtP(0), fPtN(0), fPtKP(0), fPtKN(0), fCentralityMin(0), fCentralityMax(100), fkCentralityMethod(0), fPOICuts(0), fVertexRange(0), fPhi(0), fPt(0), fEta(0), fVZEROA(0), fVZEROC(0), fTPCM(0), fDeltaDipAngle(0), fDeltaDipPt(0), fApplyDeltaDipCut(0), fDCAAll(0), fDCAXYQA(0), fDCAZQA(0), fDCAPrim(0), fDCASecondaryWeak(0), fDCAMaterial(0), fSubEventDPhiv2(0), fAnalysisSummary(0)
+{
+ // Default constructor
+ for(Int_t i(0); i < 7; i++) fPIDConfig[i] = 0.;
+ for(Int_t i(0); i < 5; i++) fDCAConfig[i] = 0.;
+ for(Int_t i(0); i < 20; i++) {
+ fVertexMixingBins[i] = 0;
+ fCentralityMixingBins[i] = 0;
+ }
+ fMixingParameters[0] = 1000; fMixingParameters[1] = 50000; fMixingParameters[2] = 5;
+ for(Int_t i(0); i < 18; i++) {
+ for(Int_t j(0); j < 2; j++) {
+ fV0Data[i][j] = 0;
+ for(Int_t k(0); k < 15; k++) fPhiMinusPsiDataContainer[k][i][j] = 0x0;
+ for(Int_t k(0); k < 15; k++) fPhiMinusPsiBackgroundContainer[k][i][j] = 0x0;
+ }
+ fResonanceSignal[i] = 0; fResonanceBackground[i] = 0; fPtSpectra[i] = 0; fPtBins[i] = 0.; fdPhiBins[i] = 0.;
+ }
+}
+//_____________________________________________________________________________
+AliAnalysisTwoParticleResonanceFlowTask::AliAnalysisTwoParticleResonanceFlowTask(const char *name) : AliAnalysisTaskSE(name), fSpeciesA(0), fSpeciesB(0), fChargeA(0), fChargeB(0), fMassA(0), fMassB(0), fDebug(0), fIsMC(0), fEventMixing(0), fPhiMinusPsiMethod(0), fQA(0), fV0(0), fMassBins(1), fMinMass(-1.), fMaxMass(0.), fCutsRP(NULL), fNullCuts(0), fPIDResponse(0), fFlowEvent(0), fBayesianResponse(0), fCandidates(0), fCandidateEtaPtCut(0), fCandidateMinEta(0), fCandidateMaxEta(0), fCandidateMinPt(0), fCandidateMaxPt(0), fNPtBins(18), fNdPhiBins(18), fCentrality(999), fVertex(999), fAOD(0), fPoolManager(0), fOutputList(0), fEventStats(0), fCentralityPass(0), fCentralityNoPass(0), fNOPID(0), fPIDk(0), fPIDp(0), fNOPIDTOF(0), fPIDTOF(0), fPtP(0), fPtN(0), fPtKP(0), fPtKN(0), fCentralityMin(0), fCentralityMax(100), fkCentralityMethod(0), fPOICuts(0), fVertexRange(0), fPhi(0), fPt(0), fEta(0), fVZEROA(0), fVZEROC(0), fTPCM(0), fDeltaDipAngle(0), fDeltaDipPt(0), fApplyDeltaDipCut(0), fDCAAll(0), fDCAXYQA(0), fDCAZQA(0), fDCAPrim(0), fDCASecondaryWeak(0), fDCAMaterial(0), fSubEventDPhiv2(0), fAnalysisSummary(0)
+{
+ // Constructor
+ for(Int_t i(0); i < 7; i++) fPIDConfig[i] = 0.;
+ for(Int_t i(0); i < 5; i++) fDCAConfig[i] = 0.;
+ for(Int_t i(0); i < 20; i++) {
+ fVertexMixingBins[i] = 0;
+ fCentralityMixingBins[i] = 0;
+ }
+ fMixingParameters[0] = 1000; fMixingParameters[1] = 50000; fMixingParameters[2] = 5;
+ for(Int_t i(0); i < 18; i++) {
+ for(Int_t j(0); j < 2; j++) {
+ fV0Data[i][j] = 0;
+ for(Int_t k(0); k < 15; k++) fPhiMinusPsiDataContainer[k][i][j] = 0x0;
+ for(Int_t k(0); k < 15; k++) fPhiMinusPsiBackgroundContainer[k][i][j] = 0x0;
+ }
+ fResonanceSignal[i] = 0; fResonanceBackground[i] = 0; fPtSpectra[i] = 0; fPtBins[i] = 0.; fdPhiBins[i] = 0.;
+ }
+ DefineInput(0, TChain::Class());
+ DefineOutput(1, TList::Class());
+ DefineOutput(2, AliFlowEventSimple::Class());
+ if(fDebug > 0) cout << " === Created instance of AliAnaysisTaskGenericResonanceFlowAnalysis === " << endl;
+}
+//_____________________________________________________________________________
+AliAnalysisTwoParticleResonanceFlowTask::~AliAnalysisTwoParticleResonanceFlowTask()
+{
+ // Destructor
+ if (fNullCuts) delete fNullCuts;
+ if (fOutputList) delete fOutputList;
+ if (fCandidates) delete fCandidates;
+ if (fFlowEvent) delete fFlowEvent;
+ if (fBayesianResponse) delete fBayesianResponse;
+ if (fEventMixing) delete fPoolManager;
+ if (fDebug > 0) cout << " === Deleted instance of AliAnalysisTwoParticleResonanceFlowTask === " << endl;
+}
+//_____________________________________________________________________________
+void AliAnalysisTwoParticleResonanceFlowTask::ForceExit(Int_t type, const char* message)
+{
+ // force exit
+ if(type==0) cout << " --> Illegal options in AddTask <-- " << endl;
+ if(type==1) cout << " --> Unknown error <-- " << endl;
+ AliFatal(message);
+}
+//_____________________________________________________________________________
+TH1F* AliAnalysisTwoParticleResonanceFlowTask::BookHistogram(const char* name)
+{
+ // Return a pointer to a TH1 with predefined binning
+ if(fDebug > 0) cout << " *** BookHistogram() *** " << name << endl;
+ TH1F *hist = new TH1F(name, Form("M_{INV} (%s)", name), 2*fMassBins, fMinMass, fMaxMass);
+ hist->GetXaxis()->SetTitle("M_{INV} (GeV / c^{2})");
+ hist->GetYaxis()->SetTitle("No. of pairs");
+ hist->SetMarkerStyle(kFullCircle);
+ hist->Sumw2();
+ fOutputList->Add(hist);
+ return hist;
+}
+//_____________________________________________________________________________
+TH2F* AliAnalysisTwoParticleResonanceFlowTask::BookPIDHistogram(const char* name, Bool_t TPC)
+{
+ // Return a pointer to a TH2 with predefined binning
+ if(fDebug > 0) cout << " *** BookPIDHisotgram() *** " << endl;
+ TH2F *hist = 0x0;
+ if(TPC) {
+ hist = new TH2F(name, Form("PID (%s)", name), 100, 0, 5, 100, 0, 1000);
+ hist->GetYaxis()->SetTitle("dE/dx (a.u.)");
+ }
+ if(!TPC) {
+ hist = new TH2F(name, Form("PID (%s)", name), 100, 0, 5, 100, 0, 1.5);
+ hist->GetYaxis()->SetTitle("#beta");
+ }
+ hist->GetXaxis()->SetTitle("P (GeV / c)");
+ fOutputList->Add(hist);
+ return hist;
+}
+//_____________________________________________________________________________
+TH1F* AliAnalysisTwoParticleResonanceFlowTask::InitPtSpectraHistograms(Float_t nmin, Float_t nmax)
+{
+ // intialize p_t histograms for each p_t bin
+ if(fDebug > 0) cout << " *** InitPtSpectraHistograms() *** " << endl;
+ TH1F* hist = new TH1F(Form("%4.2f p_{t} %4.2f", nmin, nmax), Form("%f p_{t} %f", nmin, nmax), 2*fMassBins, nmin, nmax);
+ hist->GetXaxis()->SetTitle("p_{T} GeV / c");
+ fOutputList->Add(hist);
+ return hist;
+}
+//_____________________________________________________________________________
+TH1F* AliAnalysisTwoParticleResonanceFlowTask::BookPtHistogram(const char* name)
+{
+ // Return a pointer to a p_T spectrum histogram
+ if(fDebug > 0) cout << " *** BookPtHistogram() *** " << endl;
+ TH1F* ratio = new TH1F(name, name, 100, 0, 7);
+ ratio->GetXaxis()->SetTitle("p_{T} ( GeV / c^{2} )");
+ ratio->GetYaxis()->SetTitle("No. of events");
+ ratio->Sumw2();
+ fOutputList->Add(ratio);
+ return ratio;
+}
+//_____________________________________________________________________________
+Bool_t AliAnalysisTwoParticleResonanceFlowTask::InitializeAnalysis()
+{
+ // set up the analysis
+ fNullCuts = new AliFlowTrackCuts("null_cuts");
+ fBayesianResponse = new AliFlowBayesianPID();
+ Float_t t(0);
+ for(Int_t i = 0; i < 7; i++) t+=TMath::Abs(fPIDConfig[i]);
+ if(t < 0.1) AliFatal("No valid PID procedure recognized -- terminating analysis !!!");
+ if(fNPtBins > 18) AliFatal("Invalid number of pt bins initialied ( > 18 ) -- terminating analysis !!!");
+ AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
+ cc->SetNbinsQ(500); cc->SetNbinsPhi(180); cc->SetNbinsMult(10000);
+ cc->SetQMin(0.0); cc->SetPhiMin(0.0); cc->SetMultMin(0);
+ cc->SetQMax(3.0); cc->SetPhiMax(TMath::TwoPi()); cc->SetMultMax(10000);
+ cc->SetNbinsMass(fMassBins); cc->SetNbinsEta(200); (fMassBins == 1) ? cc->SetNbinsPt(15) : cc->SetNbinsPt(100); // high pt
+ cc->SetMassMin(fMinMass); cc->SetEtaMin(-5.0); cc->SetPtMin(0);
+ cc->SetMassMax(fMaxMass); cc->SetEtaMax(+5.0); (fMassBins == 1) ? cc->SetPtMax(15) : cc->SetPtMax(10); // high pt
+ AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
+ if (man) {
+ AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
+ if (inputHandler) fPIDResponse = inputHandler->GetPIDResponse();
+ }
+ fMassA = fMassA*fMassA; //we always use mass squared in calculations, so doing it once in the init reduces cpu load
+ fMassB = fMassB*fMassB;
+ return kTRUE;
+}
+//_____________________________________________________________________________
+void AliAnalysisTwoParticleResonanceFlowTask::UserCreateOutputObjects()
+{
+ // Create user defined output objects
+ if(!InitializeAnalysis()) ForceExit(0, " > Couldn't initialize the analysis !!! <");
+ if(fDebug > 0) cout << " *** UserCreateOutputObjects() *** " << endl;
+ // Create all output objects and store them to a list
+ fOutputList = new TList();
+ fOutputList->SetOwner(kTRUE);
+ if(fQA) {
+ fEventStats = new TH1F("fHistStats", "Event Statistics", 18, -.25, 4.25);
+ fEventStats->GetXaxis()->SetTitle("No. of events");
+ fEventStats->GetYaxis()->SetTitle("Statistics");
+ fOutputList->Add(fEventStats);
+ }
+ fCentralityPass = new TH1F("fCentralityPass", "Centrality Pass", 101, -1, 100);
+ fOutputList->Add(fCentralityPass);
+ if(fQA) {
+ fCentralityNoPass = new TH1F("fCentralityNoPass", "Centrality No Pass", 101, -1, 100);
+ fOutputList->Add(fCentralityNoPass);
+ fNOPID = BookPIDHistogram("TPC signal, all particles", kTRUE);
+ fPIDk = BookPIDHistogram("TPC signal, kaons", kTRUE);
+ fPIDp = BookPIDHistogram("TPC signal, pions", kTRUE);
+ fNOPIDTOF = BookPIDHistogram("TOF signal, all particles", kFALSE);
+ fPIDTOF = BookPIDHistogram("TOF signal, kaons", kFALSE);
+ }
+ for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
+ if(!fPhiMinusPsiMethod) {
+ fResonanceSignal[ptbin] = BookHistogram(Form("NP, %4.2f < p_{T} < %4.2f GeV", fPtBins[ptbin], fPtBins[ptbin+1]));
+ fResonanceBackground[ptbin] = BookHistogram(Form("PP, %4.2f < p_{T} < %4.2f GeV", fPtBins[ptbin], fPtBins[ptbin+1]));
+ }
+ else if(fPhiMinusPsiBackgroundContainer) {
+ for(Int_t i(0); i<fNdPhiBins; i++) {// for all delta phi bins
+ fPhiMinusPsiDataContainer[i][ptbin][0] = BookHistogram(Form("%4.2f < p_{T} < %4.2f GeV, %4.2f < (#phi - #Psi_{A}) < %4.2f", fPtBins[ptbin], fPtBins[ptbin+1], fdPhiBins[i], fdPhiBins[i+1]));
+ fPhiMinusPsiDataContainer[i][ptbin][1] = BookHistogram(Form("%4.2f < p_{T} < %4.2f GeV, %4.2f < (#phi - #Psi_{C}) < %4.2f", fPtBins[ptbin], fPtBins[ptbin+1], fdPhiBins[i], fdPhiBins[i+1]));
+ fPhiMinusPsiBackgroundContainer[i][ptbin][0] = BookHistogram(Form("BG, %4.2f < p_{T} < %4.2f GeV, %4.2f < (#phi - #Psi_{A}) < %4.2f", fPtBins[ptbin], fPtBins[ptbin+1], fdPhiBins[i], fdPhiBins[i+1]));
+ fPhiMinusPsiBackgroundContainer[i][ptbin][1] = BookHistogram(Form("BG, %4.2f < p_{T} < %4.2f GeV, %4.2f < (#phi - #Psi_{C}) < %4.2f", fPtBins[ptbin], fPtBins[ptbin+1], fdPhiBins[i], fdPhiBins[i+1]));
+ }
+ }
+ }
+ if(fQA) {
+ for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) fPtSpectra[ptbin] = InitPtSpectraHistograms(fPtBins[ptbin], fPtBins[ptbin+1]);
+ fPtP = BookPtHistogram("i^{+}");
+ fPtN = BookPtHistogram("i^{-}");
+ fPtKP = BookPtHistogram("K^{+}");
+ fPtKN = BookPtHistogram("K^{-}");
+ fPhi = new TH1F("fPhi", "#phi distribution", 100, -.5, 7);
+ fOutputList->Add(fPhi);
+ fPt = new TH1F("fPt", "p_{T}", 100, 0, 5.5);
+ fOutputList->Add(fPt);
+ fEta = new TH1F("fEta", "#eta distribution", 100, -1.1, 1.1);
+ fOutputList->Add(fEta);
+ fVZEROA = new TH1F("fVZEROA", "VZERO A Multiplicity", 1000, 0, 10000);
+ fOutputList->Add(fVZEROA);
+ fVZEROC = new TH1F("fVZEROC", "VZERO C Multiplicity", 1000, 0, 10000);
+ fOutputList->Add(fVZEROC);
+ fTPCM = new TH1F("fTPCM", "TPC multiplicity", 1000, 0, 10000);
+ fOutputList->Add(fTPCM);
+ fDCAXYQA = new TH1F("fDCAXYQA", "fDCAXYQA", 1000, -5, 5);
+ fOutputList->Add(fDCAXYQA);
+ fDCAZQA = new TH1F("fDCAZQA", "fDCAZQA", 1000, -5, 5);
+ fOutputList->Add(fDCAZQA);
+ }
+ if(fIsMC || fQA) {
+ fDCAAll = new TH2F("fDCAAll", "fDCAAll", 1000, 0, 10, 1000, -5, 5);
+ fOutputList->Add(fDCAAll);
+ fDCAPrim = new TH2F("fDCAprim","fDCAprim", 1000, 0, 10, 1000, -5, 5);
+ fOutputList->Add(fDCAPrim);
+ fDCASecondaryWeak = new TH2F("fDCASecondaryWeak","fDCASecondaryWeak", 1000, 0, 10, 1000, -5, 5);
+ fOutputList->Add(fDCASecondaryWeak);
+ fDCAMaterial = new TH2F("fDCAMaterial","fDCAMaterial", 1000, 0, 10, 1000, -5, 5);
+ fOutputList->Add(fDCAMaterial);
+ }
+ if(fV0||fPhiMinusPsiMethod) {
+ fSubEventDPhiv2 = new TProfile("fSubEventDPhiv2", "fSubEventDPhiv2", 3, 0, 3);
+ fSubEventDPhiv2->GetXaxis()->SetBinLabel(1, "<cos(2(#Psi_{a} - #Psi_{b}))>");
+ fSubEventDPhiv2->GetXaxis()->SetBinLabel(2, "<cos(2(#Psi_{a} - #Psi_{c}>))");
+ fSubEventDPhiv2->GetXaxis()->SetBinLabel(3, "<cos(2(#Psi_{b} - #Psi_{c}>))");
+ fOutputList->Add(fSubEventDPhiv2);
+ const char* V0[] = {"V0A", "V0C"};
+ for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++)
+ for(Int_t iV0(0); iV0 < 2; iV0++) {
+ fV0Data[ptbin][iV0] = new TProfile(Form("%s v2 %4.2f < p_{T} < %4.2f GeV", V0[iV0], fPtBins[ptbin], fPtBins[ptbin+1]), Form("%s v2 %4.2f < p_{T} < %4.2f GeV", V0[iV0], fPtBins[ptbin], fPtBins[ptbin+1]), fMassBins, fMinMass, fMaxMass);
+ fOutputList->Add(fV0Data[ptbin][iV0]);
+ }
+ }
+ fAnalysisSummary = new TH1F("Analysis summary","Analysis summary", 26, 0, 26);
+ fAnalysisSummary->GetXaxis()->SetBinLabel(1, "fIsMC");
+ fAnalysisSummary->GetXaxis()->SetBinLabel(2, "fEventMixing");
+ fAnalysisSummary->GetXaxis()->SetBinLabel(3, "fAODAnalysis");
+ fAnalysisSummary->GetXaxis()->SetBinLabel(4, "bayesian purity");
+ fAnalysisSummary->GetXaxis()->SetBinLabel(5, "dca mode");
+ fAnalysisSummary->GetXaxis()->SetBinLabel(6, "fVertex");
+ fAnalysisSummary->GetXaxis()->SetBinLabel(7, "fCentralityMin");
+ fAnalysisSummary->GetXaxis()->SetBinLabel(8, "fCentralityMax");
+ fAnalysisSummary->GetXaxis()->SetBinLabel(9, "V0 centrality");
+ fAnalysisSummary->GetXaxis()->SetBinLabel(10, "TPC centrality");
+ fAnalysisSummary->GetXaxis()->SetBinLabel(11, "fApplyDeltaDipCut");
+ fAnalysisSummary->GetXaxis()->SetBinLabel(12, "POI filterbit");
+ fAnalysisSummary->GetXaxis()->SetBinLabel(13, "RP filterbit");
+ fAnalysisSummary->GetXaxis()->SetBinLabel(14, "PhiMinusPsiMethod");
+ fAnalysisSummary->GetXaxis()->SetBinLabel(15, "SP");
+ fAnalysisSummary->GetXaxis()->SetBinLabel(16, "SPSUB");
+ fAnalysisSummary->GetXaxis()->SetBinLabel(17, "QC");
+ fAnalysisSummary->GetXaxis()->SetBinLabel(18, "EP");
+ fAnalysisSummary->GetXaxis()->SetBinLabel(19, "EP3sub");
+ fAnalysisSummary->GetXaxis()->SetBinLabel(20, "V0_SP");
+ fAnalysisSummary->GetXaxis()->SetBinLabel(21, "eta gap");
+ fAnalysisSummary->GetXaxis()->SetBinLabel(22, "harm");
+ fAnalysisSummary->GetXaxis()->SetBinLabel(23, "highPtMode");
+ fAnalysisSummary->GetXaxis()->SetBinLabel(24, "deltaDipAngle");
+ fAnalysisSummary->GetXaxis()->SetBinLabel(25, "POImaxPt");
+ fAnalysisSummary->GetXaxis()->SetBinLabel(26, "iterator");
+ fOutputList->Add(fAnalysisSummary);
+ fAnalysisSummary->SetBinContent(1, (Float_t)fIsMC);
+ fAnalysisSummary->SetBinContent(2, (Float_t)fEventMixing);
+ fAnalysisSummary->SetBinContent(3, (Float_t)1);
+ fAnalysisSummary->SetBinContent(4, (Float_t)fPIDConfig[6]);
+ fAnalysisSummary->SetBinContent(5, (Float_t)fDCAConfig[0]);
+ fAnalysisSummary->SetBinContent(6, (Float_t)fVertexRange);
+ fAnalysisSummary->SetBinContent(7, (Float_t)fCentralityMin);
+ fAnalysisSummary->SetBinContent(8, (Float_t)fCentralityMax);
+ const char* v = "V0M"; const char* t = "TRK"; // dont compare string literals
+ if(fkCentralityMethod==v) fAnalysisSummary->SetBinContent(9, (Float_t)1);
+ if(fkCentralityMethod==t) fAnalysisSummary->SetBinContent(10, (Float_t)1);
+ fAnalysisSummary->SetBinContent(11, (Float_t)fApplyDeltaDipCut);
+ fAnalysisSummary->SetBinContent(12, (Float_t)fPOICuts->GetAODFilterBit());
+ fAnalysisSummary->SetBinContent(13, (Float_t)fCutsRP->GetAODFilterBit());
+ for(Int_t i(0); i<12;i++) fAnalysisSummary->SetBinContent(14+i, fAddTaskMacroSummary[i]);
+ fAnalysisSummary->SetBinContent(26, (Float_t)1);
+ PostData(1, fOutputList);
+ // create candidate array
+ fCandidates = new TObjArray(1000);
+ fCandidates->SetOwner(kTRUE);
+ // create and post flowevent
+ fFlowEvent = new AliFlowEvent(10000);
+ PostData(2, fFlowEvent);
+ if(fEventMixing) fPoolManager = InitializeEventMixing();
+}
+//_____________________________________________________________________________
+AliEventPoolManager* AliAnalysisTwoParticleResonanceFlowTask::InitializeEventMixing()
+{
+ // initialize event mixing
+ if(fDebug > 0) cout << " *** InitializeEventMixing() *** " << endl;
+ Int_t _c(0), _v(0);
+ for(Int_t i(0); i < 19; i++) {
+ if (fCentralityMixingBins[i+1] < fCentralityMixingBins[i]) { _c = i; break; }
+ else _c = 19;
+ }
+ for(Int_t i(0); i < 19; i++) {
+ if (fVertexMixingBins[i+1] < fVertexMixingBins[i]) { _v = i; break; }
+ else _v = 19;
+ }
+ if(fDebug > 0 ) cout << Form(" --> found %d centrality bins for mixing, %d vertex bins for mixing", _c, _v) << endl;
+ Float_t centralityBins[_c];
+ Float_t vertexBins[_v];
+ for(Int_t i(0); i < _c + 1; i++) centralityBins[i] = fCentralityMixingBins[i];
+ for(Int_t i(0); i < _v + 1; i++) vertexBins[i] = fVertexMixingBins[i];
+ return new AliEventPoolManager(fMixingParameters[0], fMixingParameters[1], _c, (Double_t*)centralityBins, _v, (Double_t*)vertexBins);
+}
+//_____________________________________________________________________________
+template <typename T> Float_t AliAnalysisTwoParticleResonanceFlowTask::InvariantMass(const T* track1, const T* track2) const
+{
+ // Return the invariant mass of two tracks
+ if(fDebug > 1) cout << " *** InvariantMass() *** " << endl;
+ if ((!track2) || (!track1)) return 0.;
+ Float_t pxs = TMath::Power((track1->Px() + track2->Px()), 2);
+ Float_t pys = TMath::Power((track1->Py() + track2->Py()), 2);
+ Float_t pzs = TMath::Power((track1->Pz() + track2->Pz()), 2);
+ Float_t e1 = TMath::Sqrt(track1->P() * track1->P() + track1->Mass());
+ Float_t e2 = TMath::Sqrt(track2->P() * track2->P() + track2->Mass());
+ Float_t es = TMath::Power((e1 + e2), 2);
+ if ((es - (pxs + pys + pzs)) < 0) return 0.;
+ return TMath::Sqrt((es - (pxs + pys + pzs)));
+}
+//_____________________________________________________________________________
+template <typename T> Float_t AliAnalysisTwoParticleResonanceFlowTask::DeltaDipAngle(const T* track1, const T* track2) const
+{
+ // Calculate the delta dip angle between two particles
+ if(fDebug > 1) cout << " *** DeltaDipAngle() *** " << endl;
+ if (track1->P()*track2->P() == 0) return 999;
+ return TMath::ACos(((track1->Pt() * track2->Pt()) + (track1->Pz() * track2->Pz())) / (track1->P() * track2->P()));
+}
+//_____________________________________________________________________________
+template <typename T> Bool_t AliAnalysisTwoParticleResonanceFlowTask::CheckDeltaDipAngle(const T* track1, const T* track2) const
+{
+ // Check if pair passes delta dip angle cut within 0 < p_t < fDeltaDipPt
+ if(fDebug > 1) cout << " *** CheckDeltaDipAngle() *** " << endl;
+ if ((TMath::Abs(DeltaDipAngle(track1, track2)) < fDeltaDipAngle) && (PairPt(track1, track2) < fDeltaDipPt)) return kFALSE;
+ return kTRUE;
+}
+//_____________________________________________________________________________
+template <typename T> Bool_t AliAnalysisTwoParticleResonanceFlowTask::CheckCandidateEtaPtCut(const T* track1, const T* track2) const
+{
+ // Check if pair passes eta and pt cut
+ if(fDebug > 1) cout << " *** CheckCandidateEtaPtCut() *** " << endl;
+ if (fCandidateMinPt > PairPt(track1, track2) || fCandidateMaxPt < PairPt(track1, track2)) return kFALSE;
+ TVector3 a(track1->Px(), track1->Py(), track1->Pz());
+ TVector3 b(track2->Px(), track2->Py(), track2->Pz());
+ TVector3 c = a + b;
+ if (fCandidateMinEta > c.Eta() || fCandidateMaxEta < c.Eta()) return kFALSE;
+ return kTRUE;
+}
+//_____________________________________________________________________________
+template <typename T> Bool_t AliAnalysisTwoParticleResonanceFlowTask::EventCut(T* event)
+{
+ // Impose event cuts
+ if(fDebug > 0) cout << " *** EventCut() *** " << endl;
+ if (!event) return kFALSE;
+ if (!CheckVertex(event)) return kFALSE;
+ if (!CheckCentrality(event)) return kFALSE;
+ if(fQA) PlotMultiplcities(event);
+ return kTRUE;
+}
+//_____________________________________________________________________________
+template <typename T> void AliAnalysisTwoParticleResonanceFlowTask::PlotMultiplcities(const T* event) const
+{
+ // QA multiplicity plots
+ if(fDebug > 1) cout << " *** PlotMultiplcities() *** " << endl;
+ fVZEROA->Fill(event->GetVZEROData()->GetMTotV0A());
+ fVZEROC->Fill(event->GetVZEROData()->GetMTotV0C());
+ fTPCM->Fill(event->GetNumberOfTracks());
+}
+//_____________________________________________________________________________
+template <typename T> Bool_t AliAnalysisTwoParticleResonanceFlowTask::CheckVertex(const T* event)
+{
+ // Check if event vertex is within given range
+ if(fDebug > 0) cout << " *** CheckVertex() *** " << endl;
+ if (!event->GetPrimaryVertex()) return 0x0;
+ fVertex = event->GetPrimaryVertex()->GetZ();
+ if (TMath::Abs(fVertex) > fVertexRange) return 0x0;
+ return kTRUE;
+}
+//_____________________________________________________________________________
+template <typename T> Bool_t AliAnalysisTwoParticleResonanceFlowTask::CheckCentrality(T* event)
+{
+ // Check if event is within the set centrality range. Falls back to V0 centrality determination if no method is set
+ if(fDebug > 0) cout << " *** CheckCentrality() *** " << endl;
+ if (!fkCentralityMethod) AliFatal("No centrality method set! FATAL ERROR!");
+ fCentrality = event->GetCentrality()->GetCentralityPercentile(fkCentralityMethod);
+ if ((fCentrality <= fCentralityMin) || (fCentrality > fCentralityMax)) {
+ if(fQA) fCentralityNoPass->Fill(fCentrality) ;
+ return kFALSE;
+ }
+ fCentralityPass->Fill(fCentrality);
+ return kTRUE;
+}
+//_____________________________________________________________________________
+void AliAnalysisTwoParticleResonanceFlowTask::InitializeBayesianPID(AliAODEvent* event)
+{
+ // Initialize the Bayesian PID object for AOD
+ if(fDebug > 0) cout << " *** InitializeBayesianPID() *** " << endl;
+ fBayesianResponse->SetDetResponse(event, fCentrality);
+}
+//_____________________________________________________________________________
+template <typename T> Bool_t AliAnalysisTwoParticleResonanceFlowTask::PassesTPCbayesianCut(T* track, Int_t species) const
+{
+ // Check if the particle passes the TPC TOF bayesian cut.
+ if(fDebug > 1) cout << " *** PassesTPCbayesianCut() *** " << endl;
+ fBayesianResponse->ComputeProb(track);
+ if (!fBayesianResponse->GetCurrentMask(0)) return kFALSE; // return false if TPC has no response
+ Int_t contamination(0); //type of particle we don't want in the sample
+ if(species==3) contamination = 2; //we want to reject pions when selection kaons
+ if(species==2) contamination = 3; // we want to reject kaons when selecting pions
+ Float_t *probabilities = fBayesianResponse->GetProb();
+ if (probabilities[species] > fPIDConfig[6] && (probabilities[contamination] < probabilities[species])) { // 3 is kaon, 2 is pion
+ if(fQA) {fPhi->Fill(track->Phi()); fPt->Fill(track->Pt()); fEta->Fill(track->Eta());}
+ return kTRUE;
+ }
+ return kFALSE;
+}
+//_____________________________________________________________________________
+Bool_t AliAnalysisTwoParticleResonanceFlowTask::PassesDCACut(AliAODTrack* track) const
+{
+ // check if track passes dca cut according to dca routine
+ // setup the routine as follows:
+ // fDCAConfig[0] < -1 no pt dependence
+ // fDCAConfig[0] = 0 do nothing
+ // fDCAConfig[0] > 1 pt dependent dca cut
+ if(fDebug > 1) cout << " *** PassesDCACut() *** " << endl;
+ if(fIsMC) return kTRUE;
+ if( (fDCAConfig[0] < 0.1) && (fDCAConfig[0] > -0.1) ) {
+ if(fQA) {
+ Double_t b[2] = { -99., -99.};
+ Double_t bCov[3] = { -99., -99., -99.};
+ track->PropagateToDCA(fAOD->GetPrimaryVertex(), fAOD->GetMagneticField(), 100., b, bCov);
+ fDCAXYQA->Fill(b[0]);
+ fDCAZQA->Fill(b[1]);
+ fDCAPrim->Fill(track->Pt(), b[0]);
+ }
+ return kTRUE;
+ }
+ Double_t b[2] = { -99., -99.};
+ Double_t bCov[3] = { -99., -99., -99.};
+ track->PropagateToDCA(fAOD->GetPrimaryVertex(), fAOD->GetMagneticField(), 100., b, bCov);
+ if((!fIsMC)&&fQA) fDCAMaterial->Fill(track->Pt(), b[0]);
+ if( (fDCAConfig[0] < -.9) && ( (TMath::Abs(b[0]) > fDCAConfig[1]) || (TMath::Abs(b[1]) > fDCAConfig[2])) ) return kFALSE;
+ if(fDCAConfig[0] > .9) {
+ if(fDCAConfig[4] < TMath::Abs(b[1])) return kFALSE;
+ Float_t denom = TMath::Power(track->Pt(), fDCAConfig[3]);
+ if( denom < 0.0000001 ) return kFALSE; // avoid division by zero
+ if( (fDCAConfig[1] + fDCAConfig[2] / denom) < TMath::Abs(b[0]) ) return kFALSE;
+ }
+ if(fQA) {
+ fDCAXYQA->Fill(b[0]);
+ fDCAZQA->Fill(b[1]);
+ fDCAPrim->Fill(track->Pt(), b[0]);
+ }
+ return kTRUE;
+}
+//_____________________________________________________________________________
+Bool_t AliAnalysisTwoParticleResonanceFlowTask::AcceptTrack(AliAODTrack* track, Int_t species) const
+{
+ if(fDebug > 1) cout << " *** IsTrack() *** " << endl;
+ if (!PassesDCACut(track)) return kFALSE;
+ if (PassesTPCbayesianCut(track, species)) {
+ if(fQA) if(species==3) fPIDk->Fill(track->P(), track->GetTPCsignal());
+ if(fQA) if(species==2) fPIDp->Fill(track->P(), track->GetTPCsignal());
+ return kTRUE;
+ }
+ return kFALSE;
+}
+//_____________________________________________________________________________
+template <typename T> Float_t AliAnalysisTwoParticleResonanceFlowTask::PairPt(const T* track1, const T* track2, Bool_t phi) const
+{
+ // return p_t of track pair, select phi to return the azimuthal angle instead
+ TVector3 a(track1->Px(), track1->Py(), track1->Pz());
+ TVector3 b(track2->Px(), track2->Py(), track2->Pz());
+ TVector3 c = a + b;
+ if(phi) return c.Phi();
+ return c.Pt();
+}
+//_____________________________________________________________________________
+template <typename T> Float_t AliAnalysisTwoParticleResonanceFlowTask::PtSelector(Int_t tracktype, const T* track1, const T* track2, Float_t mass) const
+{
+ // plot m_inv spectra of like- and unlike sign pairs
+ Float_t pt = PairPt(track1, track2);
+ if (tracktype == 0) { // signal histograms
+ for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
+ if ((fPtBins[ptbin] <= pt) && (fPtBins[ptbin+1] > pt )) {
+ fResonanceSignal[ptbin]->Fill(mass);
+ if(fQA) fPtSpectra[ptbin]->Fill(pt);
+ }
+ }
+ }
+ if (tracktype == 1) for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) if ((fPtBins[ptbin] <= pt) && (fPtBins[ptbin+1] > pt )) fResonanceBackground[ptbin]->Fill(mass);
+ return pt;
+}
+//_____________________________________________________________________________
+template <typename T> Bool_t AliAnalysisTwoParticleResonanceFlowTask::QualityCheck(T* track) const
+{
+ // check if track meets quality cuts
+ if(!track) return kFALSE;
+ return fPOICuts->IsSelected(track);
+}
+//_____________________________________________________________________________
+template <typename T> void AliAnalysisTwoParticleResonanceFlowTask::SetNullCuts(T* event)
+{
+ // Set null cuts
+ if (fDebug > 0) cout << " *** SetNullCuts() *** " << fCutsRP << endl;
+ fCutsRP->SetEvent(event, MCEvent());
+ fNullCuts->SetParamType(AliFlowTrackCuts::kGlobal);
+ fNullCuts->SetPtRange(+1, -1); // select nothing QUICK
+ fNullCuts->SetEtaRange(+1, -1); // select nothing VZERO
+ fNullCuts->SetEvent(event, MCEvent());
+}
+//_____________________________________________________________________________
+void AliAnalysisTwoParticleResonanceFlowTask::PrepareFlowEvent(Int_t iMulti)
+{
+ // Prepare flow events
+ if (fDebug > 0 ) cout << " *** PrepareFlowEvent() *** " << endl;
+ fFlowEvent->ClearFast();
+ fFlowEvent->Fill(fCutsRP, fNullCuts);
+ fFlowEvent->SetReferenceMultiplicity(iMulti);
+ fFlowEvent->DefineDeadZone(0, 0, 0, 0);
+}
+//_____________________________________________________________________________
+void AliAnalysisTwoParticleResonanceFlowTask::PhiMinusPsiMethod(TObjArray* MixingCandidates)
+{
+ // extract v2 using the phi minus psi method
+ if(fDebug > 0) cout << " ** PhiMinusPsiMethod() *** " << endl;
+ if(!AliAnalysisTaskVnV0::IsPsiComputed()) { // AliAnalysisTaskVnV0 must be added to analysis que before this task !!!
+ if(fDebug > 0 ) cout << "Fatal error: unable to retrieve VZERO task output !!! Exiting ..." << endl;
+ return;
+ }
+ // retrieve data
+ Float_t abcPsi2[] = {AliAnalysisTaskVnV0::GetPsi2V0A(), AliAnalysisTaskVnV0::GetPsi2TPC(), AliAnalysisTaskVnV0::GetPsi2V0C()};
+ fSubEventDPhiv2->Fill(0.5, TMath::Cos(2.*(abcPsi2[0]-abcPsi2[1]))); // vzeroa - tpc
+ fSubEventDPhiv2->Fill(1.5, TMath::Cos(2.*(abcPsi2[0]-abcPsi2[2]))); // vzeroa - vzeroc
+ fSubEventDPhiv2->Fill(2.5, TMath::Cos(2.*(abcPsi2[1]-abcPsi2[2]))); // tpc - vzeroc
+// if(abcPsi2[0] < 0) abcPsi2[0] += TMath::Pi();
+// if(abcPsi2[2] < 0) abcPsi2[2] += TMath::Pi();
+ // if event plane stuff went ok, fill the histograms necessary for flow analysis with phi - psi method
+ if(fDebug > 0) cout << " *** ReconstructionWithEventMixing() *** " << endl;
+ AliEventPool* pool = fPoolManager->GetEventPool(fCentrality, fVertex);
+ if(!pool) AliFatal(Form("No pool found for centrality = %f, zVtx = %f", fCentrality, fVertex));
+ if(pool->IsReady() || pool->NTracksInPool() > fMixingParameters[1] / 10 || pool->GetCurrentNEvents() >= fMixingParameters[2]) {
+ Int_t nEvents = pool->GetCurrentNEvents();
+ if(fDebug > 0) cout << " --> " << nEvents << " events in mixing buffer ... " << endl;
+ for (Int_t iEvent(0); iEvent < nEvents; iEvent++) {
+ TObjArray* mixed_candidates = pool->GetEvent(iEvent);
+ if(!mixed_candidates) continue; // this should NEVER happen
+ Int_t bufferTracks = mixed_candidates->GetEntriesFast(); // buffered candidates
+ Int_t candidates = MixingCandidates->GetEntriesFast(); // mixing candidates
+ if(fDebug > 0) cout << Form(" - mixing %d tracks with %d buffered tracks from event %d ... ", candidates, bufferTracks, iEvent) << endl;
+ TObjArray* SpeciesA = new TObjArray();
+ SpeciesA->SetOwner(kTRUE);
+ TObjArray* SpeciesB = new TObjArray();
+ SpeciesB->SetOwner(kTRUE);
+ TObjArray* SpeciesAFromBuffer = new TObjArray();
+ SpeciesAFromBuffer->SetOwner(kTRUE);
+ TObjArray* SpeciesBFromBuffer = new TObjArray();
+ SpeciesBFromBuffer->SetOwner(kTRUE);
+ for (Int_t iTracks = 0; iTracks < candidates; iTracks++) {
+ AliResonanceFlowHelperTrack* track = (AliResonanceFlowHelperTrack*)MixingCandidates->At(iTracks);
+ if (!track) continue;
+ if (track->Charge() == fChargeA && track->Species() == fSpeciesA) SpeciesA->Add(track);
+ else if (track->Charge() == fChargeB && track->Species() == fSpeciesB) SpeciesB->Add(track);
+ }
+ for (Int_t iTracks = 0; iTracks < bufferTracks; iTracks++) {
+ AliResonanceFlowHelperTrack* track = (AliResonanceFlowHelperTrack*)mixed_candidates->At(iTracks);
+ if (!track) continue;
+ if (track->Charge() == fChargeA && track->Species() == fSpeciesA ) SpeciesAFromBuffer->Add(track);
+ else if (track->Charge() == fChargeB && track->Species() == fSpeciesB) SpeciesBFromBuffer->Add(track);
+ }
+ PhiMinusPsiMethodWriteData(kTRUE, SpeciesA, SpeciesB, abcPsi2); // write signal information to histograms
+ PhiMinusPsiMethodWriteData(kFALSE, SpeciesA, SpeciesBFromBuffer, abcPsi2);
+ PhiMinusPsiMethodWriteData(kFALSE, SpeciesAFromBuffer, SpeciesB, abcPsi2);
+ } // end of mixed events loop
+ } // end of checking to see whether pool is filled correctly
+ pool->UpdatePool(MixingCandidates); // update pool with current mixing candidates (for next event)
+ if(fDebug > 0 ) pool->PrintInfo();
+}
+//_____________________________________________________________________________
+void AliAnalysisTwoParticleResonanceFlowTask::PhiMinusPsiMethodWriteData(Bool_t signal, TObjArray* SpeciesA, TObjArray* SpeciesB, Float_t* abcPsi2)
+{
+ // write data for phi minus psi method into correct histograms
+ if(fDebug > 0 ) cout << " *** PhiMinusPsiMethodWriteData *** " << endl;
+ for (Int_t i = 0; i < SpeciesA->GetEntriesFast(); i++) { // create pairs
+ for(Int_t j = 0; j < SpeciesB->GetEntriesFast(); j++) {
+ AliResonanceFlowHelperTrack* trackA = (AliResonanceFlowHelperTrack*)(SpeciesA->At(i));
+ AliResonanceFlowHelperTrack* trackB = (AliResonanceFlowHelperTrack*)(SpeciesB->At(j));
+ if(!(trackA||trackB)) continue; // shouldn't happen
+ if(trackA->ID() == trackB->ID()) continue; // remove autocorrelations
+ if(fApplyDeltaDipCut && (!CheckDeltaDipAngle(trackA, trackB))) continue;
+ if(fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(trackA, trackB))) continue;
+ Float_t minv = InvariantMass(trackA, trackB);
+ Float_t pt = PtSelector(2, trackA, trackB, minv);
+ TVector3 a(trackA->Px(), trackA->Py(), trackA->Pz());
+ TVector3 b(trackB->Px(), trackB->Py(), trackB->Pz());
+ TVector3 c = a + b;
+ Float_t phi = c.Phi();
+// if(phi < 0) phi += TMath::Pi();
+ Float_t dPhi_a = phi - abcPsi2[0];
+ Float_t dPhi_c = phi - abcPsi2[2];
+// if (dPhi_a < 0) dPhi_a += TMath::Pi();
+// if (dPhi_c < 0) dPhi_c += TMath::Pi();
+ // now all necessary numbers are here, so we can put them into histograms
+ for(Int_t k(0); k < fNdPhiBins; k++) { // dPhi bin loop
+ if(dPhi_a >= fdPhiBins[k] && dPhi_a < fdPhiBins[k+1]) // see if track is in desired dPhi_a range
+ for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) { // loop over pt bins
+ if(pt <= fPtBins[ptbin] || pt > fPtBins[ptbin+1]) continue;
+ signal ? fPhiMinusPsiDataContainer[k][ptbin][0]->Fill(minv) : fPhiMinusPsiBackgroundContainer[k][ptbin][0]->Fill(minv);
+ }
+ if(dPhi_c >= fdPhiBins[k] && dPhi_c < fdPhiBins[k+1]) // see if track is in desired dPhi_c range
+ for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) { // loop over pt bins
+ if(pt <= fPtBins[ptbin] || pt > fPtBins[ptbin+1]) continue;
+ signal ? fPhiMinusPsiDataContainer[k][ptbin][1]->Fill(minv) : fPhiMinusPsiBackgroundContainer[k][ptbin][1]->Fill(minv);
+ }
+ }
+ } // end of loop on i- tracks
+ } // end of loop on i+ tracks
+}
+//_____________________________________________________________________________
+void AliAnalysisTwoParticleResonanceFlowTask::VZEROSubEventAnalysis()
+{
+ // vzero event plane analysis using three subevents
+ if(fDebug > 0) cout << " ** VZEROSubEventAnalysis() *** " << endl;
+ if(!AliAnalysisTaskVnV0::IsPsiComputed()) { // AliAnalysisTaskVnV0 must be added to analysis que before this task !!!
+ if(fDebug > 0 ) cout << "Fatal error: unable to retrieve VZERO task output !!! Exiting ..." << endl;
+ return;
+ }
+ // retrieve data
+ Float_t abcPsi2[] = {AliAnalysisTaskVnV0::GetPsi2V0A(), AliAnalysisTaskVnV0::GetPsi2TPC(), AliAnalysisTaskVnV0::GetPsi2V0C()};
+ for(Int_t i(0); i < 3; i++) if(abcPsi2[i] == 0) {
+ if(fDebug > 0) cout << " Warning: I've encountered 0 in one of the symmetry panes (TPC, VZERO) - skipping event !" << endl;
+ return;
+ }
+ // save info for resolution calculations
+ fSubEventDPhiv2->Fill(0.5, TMath::Cos(2.*(abcPsi2[0]-abcPsi2[1]))); // vzeroa - tpc
+ fSubEventDPhiv2->Fill(1.5, TMath::Cos(2.*(abcPsi2[0]-abcPsi2[2]))); // vzeroa - vzeroc
+ fSubEventDPhiv2->Fill(2.5, TMath::Cos(2.*(abcPsi2[1]-abcPsi2[2]))); // tpc - vzeroc
+ Float_t minv[fMassBins+1];
+ Float_t _dphi[fMassBins][fNPtBins][2]; //minv, pt, v0a-c
+ Int_t occurence[fMassBins][fNPtBins]; //minv, pt
+ Float_t _inc = (fMaxMass - fMinMass) / (Float_t)fMassBins;
+ if(_inc==0) return; // avoid division by zero later on
+ for(Int_t mb(0); mb < fMassBins+1; mb++) minv[mb] = fMinMass + mb * _inc;
+ for(Int_t i(0); i < fMassBins; i++) for (Int_t j(0); j < fNPtBins; j++) for(Int_t k(0); k < 2; k ++) {
+ _dphi[i][j][k] = 0;
+ if(k==0) occurence[i][j] = 0;
+ }
+ for(Int_t iCand(0); iCand < fCandidates->GetEntriesFast(); iCand++) { // loop over unlike sign tracks
+ AliFlowTrackSimple *track = dynamic_cast<AliFlowTrackSimple*>(fCandidates->At(iCand));
+ if(!track) {
+ if(fDebug > 1) cout << " --> dynamic_cast returned null-pointer, skipping candidate <-- " << endl;
+ continue;
+ }
+ for(Int_t mb(0); mb < fMassBins; mb++) { // loop over massbands
+ if(track->Mass() <= minv[mb] || track->Mass() > minv[mb+1]) continue;
+ for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) { // loop over pt bins
+ if(track->Pt() <= fPtBins[ptbin] || track->Pt() > fPtBins[ptbin+1]) continue;
+ _dphi[mb][ptbin][0]+=TMath::Cos(2.*(track->Phi() - abcPsi2[0]));
+ _dphi[mb][ptbin][1]+=TMath::Cos(2.*(track->Phi() - abcPsi2[2]));
+ occurence[mb][ptbin]+=1;
+ }
+ }
+ for(Int_t mb(0); mb < fMassBins; mb++) // write vn values to tprofiles
+ for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
+ if(occurence[mb][ptbin]==0) continue;
+ fV0Data[ptbin][0]->Fill(mb*_inc+fMinMass+0.5*_inc, _dphi[mb][ptbin][0]/(Float_t)occurence[mb][ptbin]);
+ fV0Data[ptbin][1]->Fill(mb*_inc+fMinMass+0.5*_inc, _dphi[mb][ptbin][1]/(Float_t)occurence[mb][ptbin]);
+ }
+ }
+}
+//_____________________________________________________________________________
+void AliAnalysisTwoParticleResonanceFlowTask::UserExec(Option_t *)
+{
+ // UserExec: called for each event. Commented where necessary
+ if(fDebug > 0 ) cout << " *** UserExec() *** " << endl;
+ TObjArray* MixingCandidates = 0x0; // init null pointer for event mixing
+ if(fEventMixing) {
+ MixingCandidates = new TObjArray();
+ MixingCandidates->SetOwner(kTRUE); // mixing candidates has ownership of objects in array
+ }
+
+ TObjArray* SpeciesA = new TObjArray(); // create arrays for the helper tracks
+ SpeciesA->SetOwner(kTRUE);
+ TObjArray* ocSpeciesA = new TObjArray(); // opposite charge particles
+ ocSpeciesA->SetOwner(kTRUE);
+ TObjArray* SpeciesB = new TObjArray();
+ SpeciesB->SetOwner(kTRUE);
+ TObjArray* ocSpeciesB = new TObjArray();
+ ocSpeciesB->SetOwner(kTRUE);
+
+ if (!fPIDResponse) { // kill the event if there isn't a pid response
+ if(fDebug > 0 ) cout << " Could not get PID response " << endl;
+ return;
+ }
+
+ fAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // check for aod data type
+ if (fAOD) {
+ if (!EventCut(fAOD)) return; // check for event cuts
+ InitializeBayesianPID(fAOD); // init event objects
+ SetNullCuts(fAOD);
+ Int_t iTracks = fAOD->GetNumberOfTracks();
+ PrepareFlowEvent(iTracks); // does number of tracks correspond to the set filterbit ??? FIXME !!!
+ fCandidates->SetLast(-1);
+ if(fIsMC) IsMC(); // launch mc stuff FIXME
+ if(fQA) fEventStats->Fill(0);
+ for (Int_t i = 0; i < iTracks; i++) { // select analysis candidates
+ AliAODTrack* track = fAOD->GetTrack(i);
+ if (!QualityCheck(track)) continue; // reject poor quality tracks
+ if (fQA) { // fill qa histo's
+ if(track->Charge() > 0) {fEventStats->Fill(1); fPtP->Fill(track->Pt());}
+ if(track->Charge() < 0) {fEventStats->Fill(2); fPtN->Fill(track->Pt());}
+ }
+ if (track->Charge() == fChargeA && AcceptTrack(track, fSpeciesA)) { // store species a
+ SpeciesA->Add(new AliResonanceFlowHelperTrack(track->Eta(), track->Phi(), track->P(), track->Px(), track->Py(), track->Pz(),
+ track->Pt(), track->Charge(), fMassA, track->GetID(), fSpeciesA));
+ if(fEventMixing) MixingCandidates->Add(new AliResonanceFlowHelperTrack(track->Eta(), track->Phi(), track->P(), track->Px(),
+ track->Py(), track->Pz(), track->Pt(), track->Charge(), fMassA,
+ track->GetID(), fSpeciesA));
+ }
+ if (track->Charge() == -1*fChargeA && AcceptTrack(track, fSpeciesA)) { // store opposite charge species a
+ ocSpeciesA->Add(new AliResonanceFlowHelperTrack(track->Eta(), track->Phi(), track->P(), track->Px(), track->Py(), track->Pz(),
+ track->Pt(), track->Charge(), fMassA, track->GetID(), fSpeciesA));
+ if(fEventMixing) MixingCandidates->Add(new AliResonanceFlowHelperTrack(track->Eta(), track->Phi(), track->P(), track->Px(),
+ track->Py(), track->Pz(), track->Pt(), track->Charge(), fMassA,
+ track->GetID(), fSpeciesA));
+ }
+ if (track->Charge() == fChargeB && AcceptTrack(track, fSpeciesB)) { // store species b
+ SpeciesB->Add(new AliResonanceFlowHelperTrack(track->Eta(), track->Phi(), track->P(), track->Px(), track->Py(), track->Pz(),
+ track->Pt(), track->Charge(), fMassB, track->GetID(), fSpeciesB));
+ if(fEventMixing) MixingCandidates->Add(new AliResonanceFlowHelperTrack(track->Eta(), track->Phi(), track->P(), track->Px(),
+ track->Py(), track->Pz(), track->Pt(), track->Charge(), fMassB,
+ track->GetID(), fSpeciesB));
+ }
+ if (track->Charge() == -1*fChargeB && AcceptTrack(track, fSpeciesB)) { // store opposite charge species b
+ ocSpeciesB->Add(new AliResonanceFlowHelperTrack(track->Eta(), track->Phi(), track->P(), track->Px(), track->Py(), track->Pz(),
+ track->Pt(), track->Charge(), fMassB, track->GetID(), fSpeciesB));
+ if(fEventMixing) MixingCandidates->Add(new AliResonanceFlowHelperTrack(track->Eta(), track->Phi(), track->P(), track->Px(),
+ track->Py(), track->Pz(), track->Pt(), track->Charge(), fMassB,
+ track->GetID(), fSpeciesB));
+ }
+ }
+// cout << SpeciesA->GetEntriesFast() << endl;
+// cout << SpeciesB->GetEntriesFast() << endl;
+// cout << ocSpeciesA->GetEntriesFast() << endl;
+// cout << ocSpeciesB->GetEntriesFast() << endl;
+ // do the phi minus psi method if that's specified FIXME currenlty only supports event mixing
+ if (fPhiMinusPsiMethod) { // for the phi minus psi method, no need for flow package etc
+ PhiMinusPsiMethod(MixingCandidates); // call the method
+ PostData(1, fOutputList); // save the output data
+ return; // return to retrieve next event
+ }
+ // start the resonance reconstruction, this fills the fCandidates array with flow tracks
+ ResonanceSignal(SpeciesA, SpeciesB);
+ //
+ if(fV0) VZEROSubEventAnalysis();
+ if (fDebug > 0) printf("I received %d candidates\n", fCandidates->GetEntriesFast()); // insert candidates into flow events
+ for (int iCand = 0; iCand != fCandidates->GetEntriesFast(); ++iCand) {
+ AliFlowCandidateTrack *cand = dynamic_cast<AliFlowCandidateTrack*>(fCandidates->At(iCand));
+ if (!cand) continue;
+ if (fDebug > 1) printf(" --> Checking at candidate %d with %d daughters: mass %f\n", iCand, cand->GetNDaughters(), cand->Mass());
+ for (int iDau = 0; iDau != cand->GetNDaughters(); ++iDau) {
+ if (fDebug>1) printf(" *** Daughter %d with fID %d ***", iDau, cand->GetIDDaughter(iDau));
+ for (int iRPs = 0; iRPs != fFlowEvent->NumberOfTracks(); ++iRPs) {
+ AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack(iRPs));
+ if (!iRP) continue;
+ if (!iRP->InRPSelection()) continue;
+ if (cand->GetIDDaughter(iDau) == iRP->GetID()) {
+ if (fDebug > 1) printf(" was in RP set");
+ iRP->SetForRPSelection(kFALSE);
+ fFlowEvent->SetNumberOfRPs(fFlowEvent->GetNumberOfRPs() - 1);
+ }
+ }
+ if (fDebug > 1) printf("\n");
+ }
+ cand->SetForPOISelection(kTRUE);
+ fFlowEvent->InsertTrack(((AliFlowTrack*) cand));
+ }
+ if(!fEventMixing) { // do the combinatorial background
+ ResonanceBackground(ocSpeciesA, SpeciesB); // mix opposite charge species A with species B
+ ResonanceBackground(SpeciesA, ocSpeciesB); // mix species A with opposite charge species B
+ }
+ delete SpeciesA;
+ delete SpeciesB;
+ delete ocSpeciesA;
+ delete ocSpeciesB; // clean heap memory
+ if(fEventMixing) ReconstructionWithEventMixing(MixingCandidates);
+ PostData(1, fOutputList);
+ PostData(2, fFlowEvent);
+ }
+}
+//_____________________________________________________________________________
+void AliAnalysisTwoParticleResonanceFlowTask::ResonanceSignal(TObjArray* SpeciesA, TObjArray* SpeciesB) const
+{
+ // fill signal histograms
+ for (Int_t i(0); i < SpeciesA->GetEntriesFast(); i++) { //track loop over species A
+ for (Int_t j(0); j < SpeciesB->GetEntriesFast(); j++) { //track loop over species B
+ AliResonanceFlowHelperTrack* trackA = (AliResonanceFlowHelperTrack*)SpeciesA->At(i);
+ AliResonanceFlowHelperTrack* trackB = (AliResonanceFlowHelperTrack*)SpeciesB->At(j);
+ if(!(trackA||trackB)) continue; // shouldn't happen
+ if(trackA->ID()==trackB->ID()) continue; // beware of autocorrelations
+ if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(trackA, trackB))) continue;
+ if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(trackA, trackB))) continue;
+ Float_t mass = InvariantMass(trackA, trackB);
+ Float_t pt = PtSelector(0, trackA, trackB, mass);
+ TVector3 a(trackA->Px(), trackA->Py(), trackA->Pz());
+ TVector3 b(trackB->Px(), trackB->Py(), trackB->Pz());
+ TVector3 c = a + b;
+ Float_t phi = c.Phi();
+ Float_t eta = c.Eta();
+ Int_t nIDs[2];
+ nIDs[0] = trackA->ID();
+ nIDs[1] = trackB->ID();
+ MakeTrack(mass, pt, phi, eta, 2, nIDs);
+ }
+ }
+}
+//_____________________________________________________________________________
+void AliAnalysisTwoParticleResonanceFlowTask::ResonanceBackground(TObjArray* SpeciesA, TObjArray* SpeciesB, Bool_t checkAutoCorrelations) const
+{
+ // fill background histograms
+ for (Int_t i(0); i < SpeciesA->GetEntriesFast(); i++) { //track loop over species A
+ for (Int_t j(0); j < SpeciesB->GetEntriesFast(); j++) { //track loop over species B
+ AliResonanceFlowHelperTrack* trackA = (AliResonanceFlowHelperTrack*)SpeciesA->At(i);
+ AliResonanceFlowHelperTrack* trackB = (AliResonanceFlowHelperTrack*)SpeciesB->At(j);
+ if(!(trackA||trackB)) continue; // shouldn't happen
+ if((trackA->ID() == trackB->ID()) && checkAutoCorrelations) continue; // remove autocorrelations
+ if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(trackA, trackB))) continue;
+ if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(trackA, trackB))) continue;
+ PtSelector(1, trackA, trackB, InvariantMass(trackA, trackB));
+ }
+ }
+}
+//_____________________________________________________________________________
+void AliAnalysisTwoParticleResonanceFlowTask::ReconstructionWithEventMixing(TObjArray* MixingCandidates) const
+{
+ // perform reconstruction with event mixing
+ if(fDebug > 0) cout << " *** ReconstructionWithEventMixing() *** " << endl;
+ AliEventPool* pool = fPoolManager->GetEventPool(fCentrality, fVertex);
+ if(!pool) AliFatal(Form("No pool found for centrality = %f, zVtx = %f", fCentrality, fVertex));
+ if(pool->IsReady() || pool->NTracksInPool() > fMixingParameters[1] / 10 || pool->GetCurrentNEvents() >= fMixingParameters[2]) {
+ Int_t nEvents = pool->GetCurrentNEvents();
+ if(fDebug > 0) cout << " --> " << nEvents << " events in mixing buffer ... " << endl;
+ for (Int_t iEvent(0); iEvent < nEvents; iEvent++) {
+ TObjArray* mixed_candidates = pool->GetEvent(iEvent);
+ if(!mixed_candidates) continue; // this should NEVER happen
+ Int_t bufferTracks = mixed_candidates->GetEntriesFast(); // buffered candidates
+ Int_t candidates = MixingCandidates->GetEntriesFast(); // mixing candidates
+ if(fDebug > 0) cout << Form(" - mixing %d tracks with %d buffered tracks from event %d ... ", candidates, bufferTracks, iEvent) << endl;
+ TObjArray* SpeciesA = new TObjArray();
+ SpeciesA->SetOwner(kTRUE);
+ TObjArray* SpeciesB = new TObjArray();
+ SpeciesB->SetOwner(kTRUE);
+ TObjArray* SpeciesAFromBuffer = new TObjArray();
+ SpeciesAFromBuffer->SetOwner(kTRUE);
+ TObjArray* SpeciesBFromBuffer = new TObjArray();
+ SpeciesBFromBuffer->SetOwner(kTRUE);
+ for (Int_t iTracks = 0; iTracks < candidates; iTracks++) {
+ AliResonanceFlowHelperTrack* track = (AliResonanceFlowHelperTrack*)MixingCandidates->At(iTracks);
+ if (!track) continue;
+ if (track->Charge() == fChargeA && track->Species() == fSpeciesA) SpeciesA->Add(track);
+ else if (track->Charge() == fChargeB && track->Species() == fSpeciesB) SpeciesB->Add(track);
+ }
+ for (Int_t iTracks = 0; iTracks < bufferTracks; iTracks++) {
+ AliResonanceFlowHelperTrack* track = (AliResonanceFlowHelperTrack*)mixed_candidates->At(iTracks);
+ if (!track) continue;
+ if (track->Charge() == fChargeA && track->Species() == fSpeciesA ) SpeciesAFromBuffer->Add(track);
+ else if (track->Charge() == fChargeB && track->Species() == fSpeciesB) SpeciesBFromBuffer->Add(track);
+ }
+ ResonanceBackground(SpeciesA, SpeciesBFromBuffer, kFALSE);
+ ResonanceBackground(SpeciesB, SpeciesAFromBuffer, kFALSE);
+ } // end of mixed events loop
+ } // end of checking to see whether pool is filled correctly
+ pool->UpdatePool(MixingCandidates); // update pool with current mixing candidates (for next event)
+ if(fDebug > 0 ) pool->PrintInfo();
+}
+//_____________________________________________________________________________
+void AliAnalysisTwoParticleResonanceFlowTask::Terminate(Option_t *)
+{
+ // terminate
+ if(fDebug > 0) cout << " *** Terminate() *** " << endl;
+}
+//______________________________________________________________________________
+void AliAnalysisTwoParticleResonanceFlowTask::MakeTrack(Float_t mass, Float_t pt, Float_t phi, Float_t eta, Int_t nDau, Int_t iID[]) const
+{
+ // Construct Flow Candidate Track from two selected candidates
+ if(fDebug > 1 ) cout << " *** MakeTrack() *** " << endl;
+ Bool_t overwrite = kTRUE;
+ AliFlowCandidateTrack *sTrack = static_cast<AliFlowCandidateTrack*>(fCandidates->At(fCandidates->GetLast() + 1));
+ if (!sTrack) {
+ sTrack = new AliFlowCandidateTrack(); //deleted by fCandidates
+ overwrite = kFALSE;
+ }
+ else sTrack->ClearMe();
+ sTrack->SetMass(mass);
+ sTrack->SetPt(pt);
+ sTrack->SetPhi(phi);
+ sTrack->SetEta(eta);
+ for (Int_t iDau = 0; iDau != nDau; ++iDau) sTrack->AddDaughter(iID[iDau]);
+ sTrack->SetForPOISelection(kTRUE);
+ sTrack->SetForRPSelection(kFALSE);
+ if (overwrite) fCandidates->SetLast(fCandidates->GetLast() + 1);
+ else fCandidates->AddLast(sTrack);
+ return;
+}
+//_____________________________________________________________________________
+void AliAnalysisTwoParticleResonanceFlowTask::IsMC()
+{
+ ForceExit(1,"No MC mode available at this stage of developement ... ");
+}
+//_____________________________________________________________________________