From: snelling Date: Fri, 12 Oct 2012 13:59:20 +0000 (+0000) Subject: test task to replace phi and K* analysis X-Git-Url: http://git.uio.no/git/?a=commitdiff_plain;h=fbd0da9cd2bd5b44e7ffdc78075b1f9ff6bcd960;p=u%2Fmrichter%2FAliRoot.git test task to replace phi and K* analysis --- diff --git a/PWG/CMakelibPWGflowTasks.pkg b/PWG/CMakelibPWGflowTasks.pkg index ecb8a7636fb..d3601ddc240 100644 --- a/PWG/CMakelibPWGflowTasks.pkg +++ b/PWG/CMakelibPWGflowTasks.pkg @@ -50,6 +50,7 @@ set ( SRCS FLOW/Tasks/AliAnalysisTaskQAPmdflow.cxx FLOW/Tasks/AliFlowBayesianPID.cxx FLOW/Tasks/AliAnalysisTaskPhiFlow.cxx + FLOW/Tasks/AliAnalysisTwoParticleResonanceFlowTask.cxx FLOW/Tasks/AliAnalysisTaskFilterFE.cxx FLOW/Tasks/AliAnalysisTaskVnV0.cxx FLOW/Tasks/AliFlowVZEROResults.cxx diff --git a/PWG/FLOW/Tasks/AliAnalysisTwoParticleResonanceFlowTask.cxx b/PWG/FLOW/Tasks/AliAnalysisTwoParticleResonanceFlowTask.cxx new file mode 100644 index 00000000000..6494a5ff1ed --- /dev/null +++ b/PWG/FLOW/Tasks/AliAnalysisTwoParticleResonanceFlowTask.cxx @@ -0,0 +1,972 @@ +/************************************************************************** + * 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); iAdd(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, ""); + fSubEventDPhiv2->GetXaxis()->SetBinLabel(2, "))"); + fSubEventDPhiv2->GetXaxis()->SetBinLabel(3, "))"); + 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 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 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 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 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 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 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 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 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 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 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 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 Bool_t AliAnalysisTwoParticleResonanceFlowTask::QualityCheck(T* track) const +{ + // check if track meets quality cuts + if(!track) return kFALSE; + return fPOICuts->IsSelected(track); +} +//_____________________________________________________________________________ +template 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(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(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(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(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(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 ... "); +} +//_____________________________________________________________________________ diff --git a/PWG/FLOW/Tasks/AliAnalysisTwoParticleResonanceFlowTask.h b/PWG/FLOW/Tasks/AliAnalysisTwoParticleResonanceFlowTask.h new file mode 100644 index 00000000000..c4f78da004d --- /dev/null +++ b/PWG/FLOW/Tasks/AliAnalysisTwoParticleResonanceFlowTask.h @@ -0,0 +1,243 @@ +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. */ +/* See cxx source for full Copyright notice */ +/* $Id$ */ + +// AliAnalysisTwoParticleResonanceFlowTask: +// origin: Redmer Alexander Bertens (rbertens@nikhef.nl) +// analyis task for Resonance-meson reconstruction and estimation of v_n + +#ifndef ALIANALYSISTWOPARTICLERESONANCEFLOWTASK_H +#define ALIANALYSISTWOPARTICLERESONANCEFLOWTASK_H + +class TH1F; +class TH2F; +class TProfile; +class AliFlowTrackCuts; +class AliFlowEvent; +class AliFlowCandidateTrack; +class AliFlowBayesianPID; +class AliEventPoolManager; +class AliResonanceFlowHelperTrack; + +#include "AliAnalysisTaskSE.h" +class AliResonanceFlowHelperTrack : public TObject +{ +public: + AliResonanceFlowHelperTrack(Double_t eta, Double_t phi, Double_t p, Double_t px, Double_t py, Double_t pz, Double_t pt, Int_t charge, Double_t mass, Int_t id, Int_t species) : fEta(eta), fPhi(phi), fp(p), fpX(px), fpY(py), fpZ(pz), fpT(pt), fCharge(charge), fMass(mass), fID(id), fSpecies(species) { } + ~AliResonanceFlowHelperTrack() {} + virtual Double_t P() const { return fp; } + virtual Double_t Px() const { return fpX; } + virtual Double_t Py() const { return fpY; } + virtual Double_t Pz() const { return fpZ; } + virtual Double_t Pt() const { return fpT; } + virtual Double_t Phi() const { return fPhi; } + virtual Double_t Eta() const { return fEta; } + virtual Int_t Charge() const { return fCharge; } + virtual Double_t Mass() const { return fMass; } + virtual Short_t ID() const { return fID; } + virtual Int_t Species() const { return fSpecies; } + void InitializeHelperTrack(Double_t eta, Double_t phi, Double_t p, Double_t px, Double_t py, Double_t pz, Double_t pt, Int_t charge, Double_t mass, Int_t id, Int_t species) { fEta = eta; fPhi = phi; fp = p; fpX = px; fpY = py; fpZ = pz; fpT = pt; fCharge = charge; fMass = mass; fID = id; fSpecies = species; } +private: + Double_t fEta; // eta + Double_t fPhi; // phi + Double_t fp; // p + Double_t fpX; // pX + Double_t fpY; // pY + Double_t fpZ; // pZ + Double_t fpT; // pT + Int_t fCharge; // charge + Double_t fMass; // mass + Short_t fID; // id + Int_t fSpecies; // species + ClassDef(AliResonanceFlowHelperTrack, 1); // lightweight helper track for phi reconstruction +}; + +class AliAnalysisTwoParticleResonanceFlowTask : public AliAnalysisTaskSE +{ +public: + AliAnalysisTwoParticleResonanceFlowTask(); + AliAnalysisTwoParticleResonanceFlowTask(const char *name); + virtual ~AliAnalysisTwoParticleResonanceFlowTask(); + // technical aspects of the analysis + void ForceExit(Int_t type, const char* message); + Int_t SetDebugLevelResonanceTask(Int_t debug) {fDebug = debug; return fDebug; } + Bool_t SetIsMC(Bool_t ismc) {fIsMC = ismc; return fIsMC; } + void IsMC(); + Bool_t UseEventMixing(Bool_t mix) { fEventMixing = mix; return mix; } + Bool_t UsePhiMinusPsiMethod(Bool_t p) {fPhiMinusPsiMethod = p; return p;} + Bool_t SetVZEROSubEvents(Bool_t v0) { fV0 = v0; return v0; } + // configure the output of the analysis + TH1F* BookHistogram(const char * name); + TH2F* BookPIDHistogram(const char * name, Bool_t TPC); + TH1F* InitPtSpectraHistograms(Float_t nmin, Float_t nmax); + TH1F* BookPtHistogram(const char* name); + Bool_t InitializeAnalysis(); + void AddResonanceIdentificationOutputObjects(); + virtual void UserCreateOutputObjects(); + // setters + void SetPtBins(Float_t bin[19], Int_t n) { for(Int_t i = 0; i < n+1; i++) fPtBins[i] = bin[i]; fNPtBins = n; } + void SetdPhiBins(Float_t bin[19], Int_t n) { for(Int_t i = 0; i < n+1; i++) fdPhiBins[i] = bin[i]; fNdPhiBins = n;} + void SetCentralityParameters(Float_t CentralityMin, Float_t CentralityMax, const char* CentralityMethod) { + fCentralityMin = CentralityMin; + fCentralityMax = CentralityMax; + fkCentralityMethod = CentralityMethod; } + void SetPOICuts(AliFlowTrackCuts *cutsPOI) { fPOICuts = cutsPOI; } + void SetRPCuts(AliFlowTrackCuts *cutsRP) { fCutsRP = cutsRP; } + void SetPIDConfiguration(Float_t prob[7]) { for(Int_t i = 0; i < 7; i++) fPIDConfig[i] = prob[i]; } + Bool_t SetQA(Bool_t qa) {fQA = qa; return fQA;} + void SetAddTaskMacroSummary(Float_t m[12]) {for(Int_t i(0); i < 12; i++) fAddTaskMacroSummary[i] = m[i];} + void SetPOIDCAXYZ(Float_t dca[5]) { for(Int_t i = 0; i < 5; i++) fDCAConfig[i] = dca[i]; } + void SetMixingParameters(Int_t p[3]) { for(Int_t i = 0; i < 3; i++) fMixingParameters[i] = p[i]; } + void SetupSpeciesA(Int_t species, Int_t charge, Float_t mass) {fSpeciesA = species; fChargeA = charge; fMassA = mass;} + void SetupSpeciesB(Int_t species, Int_t charge, Float_t mass) {fSpeciesB = species; fChargeB = charge; fMassB = mass;} + void SetCandidateEtaAndPt(Float_t mineta, Float_t maxeta, Float_t minpt, Float_t maxpt) { fCandidateMinEta = mineta; + fCandidateMaxEta = maxeta; + fCandidateMinPt = minpt; + fCandidateMaxPt = maxpt; + fCandidateEtaPtCut = kTRUE;} + void SetCommonConstants(Int_t massBins, Float_t minMass, Float_t maxMass) { fMassBins = massBins; + fMinMass = minMass; + fMaxMass = maxMass; } + void SetVertexZ(Float_t z) { fVertexRange = z; } + void SetMaxDeltaDipAngleAndPt(Float_t a, Float_t pt) { fDeltaDipAngle = a; + fDeltaDipPt = pt; + fApplyDeltaDipCut = kTRUE; }; + void SetMixingBins(Int_t c[20], Int_t v[20]) {for(Int_t i = 0; i < 20; i++) { fCentralityMixingBins[i] = c[i]; + fVertexMixingBins[i] = v[i]; } } + //getters + void GetMixingParameters(Int_t p[3]) const { for(Int_t i = 0; i < 3; i++) p[i] = fMixingParameters[i]; } + Float_t GetCenMin() const {return fCentralityMin; } + Float_t GetCenMax() const {return fCentralityMax; } + const char* GetCentralityMethod() const {return fkCentralityMethod; } + Float_t GetVertexZ() const { return fVertexRange; } + Float_t GetDeltaDipAngle() const {return fDeltaDipAngle; } + Float_t GetDeltaDipPt() const {return fDeltaDipPt; } + void GetPIDConfiguration(Float_t prob[7]) const {for(Int_t i = 0; i < 7; i++) prob[i] = fPIDConfig[i]; } + void GetPOIDCZXYZ(Float_t dca[5]) const { for(Int_t i = 0; i < 5; i++) dca[i] = fDCAConfig[i]; } + void GetCandidateEtaAndPt(Float_t etapt[4]) const { etapt[0] = fCandidateMinEta; + etapt[1] = fCandidateMaxEta; + etapt[2] = fCandidateMinPt; + etapt[3] = fCandidateMaxPt; } + // the analysis itself + AliEventPoolManager* InitializeEventMixing(); + template Float_t InvariantMass(const T* track1, const T* track2) const; + template Float_t DeltaDipAngle(const T* track1, const T* track2) const; + template Bool_t CheckDeltaDipAngle(const T* track1, const T* track2) const; + template Bool_t CheckCandidateEtaPtCut(const T* track1, const T* track2) const; + template Bool_t EventCut(T* event); + template void PlotMultiplcities(const T* event) const; + template Bool_t CheckVertex(const T* event); + template Bool_t CheckCentrality(T* event); + void InitializeBayesianPID(AliAODEvent* event); + template Bool_t PassesTPCbayesianCut(T* track, Int_t species) const; + Bool_t PassesDCACut(AliAODTrack* track) const; + Bool_t AcceptTrack(AliAODTrack* track, Int_t species) const; + template Float_t PairPt(const T* track_1, const T* track_2, Bool_t phi = kFALSE) const; + template Float_t PtSelector(Int_t _track_type, const T* track_1, const T* track_2, Float_t mass) const; + template Bool_t QualityCheck(T* track) const; + template void SetNullCuts(T* esd); + void PrepareFlowEvent(Int_t iMulti); + void PhiMinusPsiMethod(TObjArray* MixingCandidates); + void PhiMinusPsiMethodWriteData(Bool_t signal, TObjArray* SpeciesA, TObjArray* SpeciesB, Float_t* abcPsi2); + void VZEROSubEventAnalysis(); + virtual void UserExec(Option_t *option); + void ResonanceSignal(TObjArray* SpeciesA, TObjArray* SpeciesB) const; + void ResonanceBackground(TObjArray* SpeciesA, TObjArray* SpeciesB, Bool_t checkAutoCorrelations = kTRUE) const; + void ReconstructionWithEventMixing(TObjArray* MixingCandidates) const; + virtual void Terminate(Option_t *); + +private: + + Int_t fSpeciesA; // particle species a + Int_t fSpeciesB; // species b + Int_t fChargeA; // charge for species a + Int_t fChargeB; // charge for species b + Float_t fMassA; // mass for species a + Float_t fMassB; // mass for species b + Int_t fDebug; // debug level (0 none, 1 fcn calls, 2 verbose) + Bool_t fIsMC; // use mc mode + Bool_t fEventMixing; // use event mixing + Bool_t fPhiMinusPsiMethod; // use phi minus psi method (default is invariant mass fit method) + Bool_t fQA; // make qa plots + Bool_t fV0; // use three subevents including vzero + Int_t fMassBins; // mass bins + Float_t fMinMass; // mass range + Float_t fMaxMass; // mass range + AliFlowTrackCuts *fCutsRP; // track cuts for reference particles + AliFlowTrackCuts *fNullCuts; // dummy cuts for flow event tracks + AliPIDResponse *fPIDResponse; //! pid response object + AliFlowEvent *fFlowEvent; //! flow events (one for each inv mass band) + AliFlowBayesianPID *fBayesianResponse; //!PID response object + TObjArray *fCandidates; // candidate array + Bool_t fCandidateEtaPtCut; // set eta and pt cut for candidate tracks and combinatorial background + Float_t fCandidateMinEta; // minimum eta for candidates + Float_t fCandidateMaxEta; // maximum eta for candidates + Float_t fCandidateMinPt; // minimum pt for candidates + Float_t fCandidateMaxPt; // maximum pt for candidates + Float_t fPIDConfig[7]; // configure pid routine + Float_t fDCAConfig[5]; // configure dca routine + Int_t fMixingParameters[3]; // mixing: poolsize, mixing tracks, pool buffer + Int_t fCentralityMixingBins[20]; // configure centrality bins for event mixing + Int_t fVertexMixingBins[20]; // configure vertex bins for event mixing + Float_t fPtBins[19]; // pt bin borders + Float_t fdPhiBins[19]; // dPhi bin borders + Int_t fNPtBins; // no of pt bins + 1 + Int_t fNdPhiBins; // no of dphi bins + 1 + Float_t fCentrality; // event centrality + Float_t fVertex; // event vertex z + AliAODEvent *fAOD; //! AOD oject + AliEventPoolManager *fPoolManager; //! event pool manager + TList *fOutputList; // ! Output list + TH1F *fEventStats; // ! Histogram for event statistics + TH1F *fCentralityPass; // ! QA histogram of events that pass centrality cut + TH1F *fCentralityNoPass; //! QA histogram of events that do not pass centrality cut + TH2F *fNOPID;//! QA histogram of TPC response of all charged particles + TH2F *fPIDk;//! QA histogram of TPC response of kaons + TH2F *fPIDp; //! QA histogram of TPC response of pions + TH2F *fNOPIDTOF; //! QA histo of TOF repsonse charged particles + TH2F *fPIDTOF; //! QA histo of TOF response kaons + TH1F *fResonanceSignal[18]; //! signal histograms + TH1F *fResonanceBackground[18]; //! like-sign kaon pairs + TH1F *fPtSpectra[18]; //! pt spectra + TH1F *fPtP; //! QA histogram of p_t distribution of positive particles + TH1F *fPtN; //! QA histogram of p_t distribution of negative particles + TH1F *fPtKP; //! QA histogram of p_t distribution of positive kaons + TH1F *fPtKN; //! QA histogram of p_t distribution of negative kaons + Float_t fCentralityMin; // lower bound of cenrality bin + Float_t fCentralityMax; // upper bound of centrality bin + const char *fkCentralityMethod; // method used to determine centrality (V0 by default) + AliFlowTrackCuts *fPOICuts; // cuts for particles of interest (flow package) + Float_t fVertexRange; // absolute value of maximum distance of vertex along the z-axis + TH1F *fPhi; //! QA plot of azimuthal distribution of tracks used for event plane estimation + TH1F *fPt; //! QA plot of p_t sectrum of tracks used for event plane estimation + TH1F *fEta; //! QA plot of eta distribution of tracks used for event plane estimation + TH1F *fVZEROA; //! QA plot vzeroa multiplicity (all tracks in event) + TH1F *fVZEROC; //! QA plot vzeroc multiplicity (all tracks in event) + TH1F *fTPCM; //! QA plot TPC multiplicity (tracks used for event plane estimation) + Float_t fDeltaDipAngle; // absolute value of delta dip angle to be excluded + Float_t fDeltaDipPt; // upper value of pt range in which delta dip angle must be applied + Bool_t fApplyDeltaDipCut; // enforce delta dip cut + TH2F *fDCAAll;//! qa dca of all charged particles + TH1F *fDCAXYQA; //! qa plot of dca xz + TH1F *fDCAZQA; //!qa plot of dca z + TH2F *fDCAPrim; //!dca of primaries (mc) or kaons (data) + TH2F *fDCASecondaryWeak; //! dca of weak (mc) + TH2F *fDCAMaterial; //!dca material (mc) all (data) + TProfile *fSubEventDPhiv2; //! subevent resolution info for v2 + TProfile *fV0Data[18][2]; //! profiles for vzero vn(minv) + TH1F *fPhiMinusPsiDataContainer[18][18][2]; //! histograms for phi minus psi results + TH1F *fPhiMinusPsiBackgroundContainer[18][18][2]; //! histograms for phi minus psi background + TH1F *fAnalysisSummary; //! plot analysis flags + Float_t fAddTaskMacroSummary[12]; // add task macro summary + + AliAnalysisTwoParticleResonanceFlowTask(const AliAnalysisTwoParticleResonanceFlowTask&); // Not implemented + AliAnalysisTwoParticleResonanceFlowTask& operator=(const AliAnalysisTwoParticleResonanceFlowTask&); // Not implemented + void MakeTrack(Float_t, Float_t, Float_t, Float_t, Int_t , Int_t[]) const; + + ClassDef(AliAnalysisTwoParticleResonanceFlowTask, 1); +}; + +#endif + + + diff --git a/PWG/PWGflowTasksLinkDef.h b/PWG/PWGflowTasksLinkDef.h index bd38edb0747..e5c7de23ca0 100644 --- a/PWG/PWGflowTasksLinkDef.h +++ b/PWG/PWGflowTasksLinkDef.h @@ -29,6 +29,8 @@ #pragma link C++ class AliFlowBayesianPID+; #pragma link C++ class AliAnalysisTaskPhiFlow+; #pragma link C++ class AliPhiMesonHelperTrack+; +#pragma link C++ class AliAnalysisTwoParticleResonanceFlowTask+; +#pragma link C++ class AliResonanceFlowHelperTrack+; #pragma link C++ class AliAnalysisTaskFilterFE+; #pragma link C++ class AliAnalysisTaskVnV0+; #pragma link C++ class AliFlowVZEROResults+; diff --git a/PWGCF/FLOW/macros/AddTwoParticleResonanceFlowTask.C b/PWGCF/FLOW/macros/AddTwoParticleResonanceFlowTask.C new file mode 100644 index 00000000000..ab0f12d124e --- /dev/null +++ b/PWGCF/FLOW/macros/AddTwoParticleResonanceFlowTask.C @@ -0,0 +1,561 @@ +///////////////////////////////////////////////////////////////////////////////////////////// +// +// AddTaskPhiFlow macro +// Author: Redmer A. Bertens, Utrecht University, 2012 +// rbertens@cern.ch, r.a.bertens@uu.nl +// Commented where necessary +///////////////////////////////////////////////////////////////////////////////////////////// + +class AliAnalysisDataContainer; +class AliFlowTrackCuts; +class AliFlowTrackSimpleCuts; +class AliFlowEventCuts; +class AliFlowEventSimpleCuts; +class AliAnalysisDataContainer; + +AliAnalysisTwoParticleResonanceFlowTask* AddTwoParticleResonanceFlowTask(Bool_t SP = 1, // select flow analysis methods + Bool_t SPSUB = 0, + Bool_t QC = 0, + Bool_t EP = 0, + Bool_t EP3sub = 0, + Bool_t VZERO_SP = 0, // use vzero sp method + Float_t centrMin = 20., // centrality selection + Float_t centrMax = 30., + Float_t ITSsigma = 0., // pid mode (see task implementation) + Float_t ITSrange = 0., + Float_t TPCcontrol = 1., + Float_t TPCsigma = 3., + Float_t TPCrange = 0., + Float_t ITScontrol = -1., + Float_t Bpurity = 0.3, + TString suffixName = "UniqueID", // unique id for output objects + Bool_t bCentralTrigger = kTRUE, // trigger selection + Float_t EtaGap = 0., // eta gap for SPSUB + TString DCA = "pt", // dca mode (see task implementation) + Int_t harm = 2, // harmonic vn + UInt_t poi_filter = 1, // aod filterbits + UInt_t rp_filter = 1, + Float_t minMass = .71, // min mass for flow analysis + Float_t maxMaxx = 1.13, // max mass for flow analysis + Int_t speciesA = 3, // species a + Int_t speciesB = 2, // species b + Int_t chargeA = -1, // charge a + Int_t chargeB = 1, // charge b + Float_t massA = 4.93676999999999977e-01, // mass species a + Float_t massB = 1.39570e-01, // mass species b + Bool_t PhiMinusPsiMethod = 0, // use phi minus psi method + Bool_t event_mixing = kTRUE, // use event mixing + Bool_t highPtMode = kFALSE, // use with caution !!! disables invariant mass fit method + Float_t deltaPhiMass = 0.0003, // dM in which to look for phi + Float_t POIPtMax = 5., // max pt of daughterp particles + Bool_t shrinkSP = kTRUE, // shrink output + Bool_t debug = kTRUE) // macro debug mode, for task's debug mode see header +{ + /* some defaults that are used frequently: + * for phi: + * - minMass = 0.99 + * - maxMass = 1.092 + * for kstar + * - minMass = 0.71 + * - maxMass = 1.13 + * kaon mass = 4.93676999999999977e-01 + * pion mass = 1.39570e-01 + */ + // some defaults that have been removed as function arguments (august 30 2012) + Float_t POIPtMin = 0.2; // pt of daughters + Float_t deltaDip = 0.; + Float_t deltaDipMaxPt = 0.; + Bool_t TPCStandAloneTracks = kFALSE; // deprecated + Bool_t useGlobalRPCuts = kTRUE; // deprecated + Float_t vertexZ = 10.; + Float_t POIEtaMin = -0.8; + Float_t POIEtaMax = 0.8; + // start of macro + Float_t PIDconfig[] = {ITSsigma, ITSrange, TPCcontrol, TPCsigma, TPCrange, ITScontrol, Bpurity}; + // main function, create and add tasks + if(debug) cout << " === Adding Task PhiFlow === " << endl; + // set up main output container's name + TString fileName = AliAnalysisManager::GetCommonFileName(); + fileName += ":PhiReconstruction"; + suffixName += Form("%.0f", centrMin); + suffixName += Form("%.0f", centrMax); + fileName+=suffixName; + if(debug) cout << " --> Reconstruction data container: " << fileName << endl; + // check validity of arguments + if((!SPSUB)&&(EtaGap > 0.)) { + if(debug) cout << " --> Fatal error: Eta gap is introduced but method SPSUB is not enabled !!! <-- " << endl; + return 0x0; + } + else if ((QC||SP||EP||EP3sub)&&(EtaGap > 0.)) { + if(debug) cout << " --> Fatal error: one or more of the flow analyses is not compatible with the Eta Gap !!! <--" << endl; + return 0x0; + } + else if(EP3sub) { + if(highPtMode) { + if(debug) cout << " --> Can't launch 3 subevent method for high pt analysis, exiting ... <--" << endl; + return 0x0; + } + if(debug) cout << " --> Starting 3 subevent plane method - try at your own risk !!! <-- " << endl; + if(!(harm!=2||harm!=3)) { + if(debug) cout << " --> Fatal error: can only return v2 and v3 with 3 subevent method! " << endl; + return 0x0; + } + } + if((EP3sub||PhiMinusPsiMethod)&&debug) { + gROOT->LoadMacro("$ALICE_ROOT/PWGCF/FLOW/macros/AddTaskVZERO.C"); + AddTaskVZERO(0,0,0,0); + } + // get the manager and input event handler + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if (!mgr) { + if(debug) cout << " Fatal error: no analysis manager found! " << endl; + return 0x0; + } + if (!mgr->GetInputEventHandler()) { + if(debug) cout << " Fatal error: no imput event handler found!" << endl; + return 0x0; + } + if(EP3sub&&debug) { // don't use this on train ! this is why it's only enabled in macro debug mode + gROOT->LoadMacro("$ALICE_ROOT/PWGCF/FLOW/macros/AddTaskVZERO.C"); + AddTaskVZERO(0,0,0,0); + } + // create the main task + AliAnalysisTwoParticleResonanceFlowTask *task = new AliAnalysisTwoParticleResonanceFlowTask("TaskPhiFlow"); + if(debug) cout << " === AliAnalysisTaskPhiFlow === " << task << endl; + if(!task) { + if(debug) cout << " --> Unexpected error occurred: NO TASK WAS CREATED! (could be a library problem!) " << endl; + return 0x0; + } + if(task->UsePhiMinusPsiMethod(PhiMinusPsiMethod)) { + cout << " --> Using phi - Psi method <-- ... " << endl; + Float_t dphibins[] = {0., 0.63, 1.26, 1.89, 2.52, 3.15}; + task->SetdPhiBins(dphibins, (Int_t)(sizeof(dphibins)/sizeof(dphibins[1])-1)); + event_mixing = kTRUE; + } + task->SetupSpeciesA(speciesA, chargeA, massA); + task->SetupSpeciesB(speciesB, chargeB, massB); + + if(task->SetVZEROSubEvents(EP3sub)) cout << " --> Setting up VZERO subevents method ... " << endl; + if(event_mixing) { + if(debug) cout << " --> Enabeling event mixing for phi reconstruction - try at your own risk !!! <-- " << endl; + // set vertex and mixing bins - arrays MUST have length 20! + Int_t c[] = {0, 2, 4, 6, 8, 10, 20, 30, 40, 50, 60, 70, 80, 90, 101, 0, 0, 0, 0, 0}; + Int_t v[] = {-10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0}; + if(((Int_t)(sizeof(c)/sizeof(c[1]))!=20)||((Int_t)(sizeof(v)/sizeof(v[1]))!=20)) { + cout << " --> Fatal error: check mixing parameters ! <-- " << endl; + return 0x0; + } + else task->SetMixingBins(c, v); + } + // set triggers + if (bCentralTrigger) task->SelectCollisionCandidates(AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral); + else task->SelectCollisionCandidates(AliVEvent::kMB); + if(debug) cout << " --> Set trigger selection to "; + if(debug&&bCentralTrigger) cout << " kMB, kCentral, kSemiCentral " << endl; + if(debug&&(!bCentralTrigger)) cout << " kMB " << endl; + //set RP cuts for flow package analysis + AliFlowTrackCuts* cutsRP = new AliFlowTrackCuts("RFPcuts"); + if(!cutsRP) { + if(debug) cout << " Fatal error: no RP cuts found, could be a library problem! " << endl; + return 0x0; + } + if(useGlobalRPCuts&&(!VZERO_SP)) { + AliFlowTrackCuts::trackParameterType rptype = AliFlowTrackCuts::kGlobal; + cutsRP->SetParamType(rptype); + cutsRP->SetPtRange(0.2, 5.0); + cutsRP->SetEtaRange(-0.8, 0.8); + cutsRP->SetMinNClustersTPC(70); + cutsRP->SetMinChi2PerClusterTPC(0.1); + cutsRP->SetMaxChi2PerClusterTPC(4.0); + cutsRP->SetRequireTPCRefit(kTRUE); + cutsRP->SetMaxDCAToVertexXY(0.3); + cutsRP->SetMaxDCAToVertexZ(0.3); + cutsRP->SetAcceptKinkDaughters(kFALSE); + cutsRP->SetMinimalTPCdedx(10.); + if(rp_filter < 9999 ) { + if(debug) cout << " > set RP filterbit " << rp_filter << endl; + cutsRP->SetAODfilterBit(rp_filter); + } + if(debug) cout << " --> kGlobal RP's " << cutsRP << endl; + } + if(VZERO_SP) { // use vzero sub analysis + cutsRP = cutsRP->GetStandardVZEROOnlyTrackCuts(); // select vzero tracks + SP = kFALSE; // disable other methods + SPSUB = kTRUE; // calculate sp_qa and sp_qb + QC = kFALSE; + EP = kFALSE; + EP3sub = kFALSE; + EtaGap = 0.; // no eta gap, full tpc poi's + if(debug) cout << " --> VZERO RP's " << cutsRP << endl; + } + task->SetRPCuts(cutsRP); + // set POI cuts for kaon selection + AliFlowTrackCuts* cutsPOI = new AliFlowTrackCuts("GlobalPOI"); + if(!cutsPOI) { + if(debug) cout << " Fatal error: no POI cuts (could be a library problem)!" << endl; + return 0x0; + } + cutsPOI = cutsPOI->GetStandardGlobalTrackCuts2010(); + cutsPOI->SetPtRange(POIPtMin, POIPtMax); // pt range of DAUGHTERS !!! + cutsPOI->SetMaxDCAToVertexXY(0.3); // FIXME not implemented in aod086 aod095 see PassesDCACut() in implementation + cutsPOI->SetMaxDCAToVertexZ(0.3); + if(poi_filter < 9999 ) { + if(debug) cout << " > set POI filterbit " << poi_filter << endl; + cutsPOI->SetAODfilterBit(poi_filter); + } + if(debug) cout << " --> cutsPOI " << cutsPOI << endl; + task->SetPOICuts(cutsPOI); + //set POI cuts for aods XY Z - 3 distinct cases + Float_t POIDCA[] = {0., 0., 0., 0., 0.}; + if(DCA == "none" ) { // 1 --- do nothing + if (debug) cout << " --> No DCA cut on POI's <-- " << endl; + for (Int_t i = 0; i < 5; i++) POIDCA[i] = 0.; + } + if(DCA == "fix" ) { // 2 --- use fixed values for xy z + if (debug) cout << " --> Fixed DCA cut on POI's <-- " << endl; + POIDCA[0] = -1.; POIDCA[1] = 0.3; POIDCA[2] = 0.3; POIDCA[3] = 0.; POIDCA[4] = 0.; + } + if(DCA == "pt" ) { // 3 --- use pt dependent cut + if (debug) cout << " --> Pt dependent DCA cut on POI's <-- " << endl; + POIDCA[0] = 1.; POIDCA[1] = 0.0105; POIDCA[2] = 0.0350; POIDCA[3] = 1.1; POIDCA[4] = 2.; + } + task->SetPOIDCAXYZ(POIDCA); + Float_t AddTaskMacroSummary[] = {(Float_t)0, (Float_t)SP, (Float_t)SPSUB, (Float_t)QC, (Float_t)EP, (Float_t)EP3sub, (Float_t)VZERO_SP, (Float_t)EtaGap, (Float_t)harm, (Float_t)highPtMode, (Float_t)deltaDip, (Float_t)POIPtMax}; + task->SetAddTaskMacroSummary(AddTaskMacroSummary); + + if(highPtMode) { // high pt loop - loop will end macro + Float_t _pt[] = {0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 12., 15.}; + task->SetPtBins(_pt, (Int_t)(sizeof(_pt)/sizeof(_pt[1]))-1); + // general approach: use kinematic filters which select kaon pairs with a certain mass window + AliFlowTrackSimpleCuts* HighPtSubEventFilterA = new AliFlowTrackSimpleCuts("HighPtSubEventFilterA"); + HighPtSubEventFilterA->SetEtaMin(-0.8); + HighPtSubEventFilterA->SetEtaMax(0.0); + HighPtSubEventFilterA->SetMassMin(1.019445 - deltaPhiMass); + HighPtSubEventFilterA->SetMassMax(1.019445 + deltaPhiMass); + AliFlowTrackSimpleCuts* HighPtSubEventFilterB = new AliFlowTrackSimpleCuts("HighPtSubEventFilterB"); + HighPtSubEventFilterB->SetEtaMin(0.0); + HighPtSubEventFilterB->SetEtaMax(+0.8); + HighPtSubEventFilterB->SetMassMin(1.019445 - deltaPhiMass); + HighPtSubEventFilterB->SetMassMax(1.019445 + deltaPhiMass); + AliFlowTrackSimpleCuts* HighPtGenericFilter = new AliFlowTrackSimpleCuts("HighPtGenericFilter"); + HighPtGenericFilter->SetEtaMin(-0.8); + HighPtGenericFilter->SetEtaMax(+0.8); + HighPtGenericFilter->SetMassMin(1.019445 - deltaPhiMass); + HighPtGenericFilter->SetMassMax(1.019445 + deltaPhiMass); + if(debug) cout << " --> Created poi filters " << endl; + // set pair and event cuts + if((deltaDip>0.005)&&(deltaDipMaxPt>0.005)) task->SetMaxDeltaDipAngleAndPt(deltaDip, deltaDipMaxPt); + else cout << " --> Disabled Delta-Dip exclusion. <-- " << endl; + task->SetCandidateEtaAndPt(POIEtaMin, POIEtaMax, 0., 15.); + task->SetCentralityParameters(centrMin, centrMax, "TRK"); + task->SetVertexZ(vertexZ); + if(debug) cout << " --> Set pair cuts and event cuts" << endl; + // specify the PID procedure which will be used + task->SetPIDConfiguration(PIDconfig); + AliAnalysisDataContainer *cinput = mgr->GetCommonInputContainer(); + AliAnalysisDataContainer *coutHist = mgr->CreateContainer(Form("PhiV2_%s", OutputName(centrMin, centrMax,PIDconfig, suffixName, bCentralTrigger, EtaGap, POIEtaMin, POIEtaMax, 0., 15., deltaDip, deltaDipMaxPt, DCA, harm, TPCStandAloneTracks, vertexZ, debug, useGlobalRPCuts).Data()), TList::Class(), AliAnalysisManager::kOutputContainer, fileName); + if(debug) cout << " --> Created IO containers " << cinput << ", " << coutHist << endl; + mgr->AddTask(task); + if(debug) cout << " === ADDING MAIN TASK == " << endl; + mgr->ConnectInput(task, 0, cinput); + mgr->ConnectOutput(task, 1, coutHist); + if(debug) cout << " --> Connected IO containers " << endl; + if (SP || EP || QC || SPSUB) // if flow analysis should be done after reconstruction + { + Int_t mb(999); + if(debug) cout << " === RECEIVED REQUEST FOR FLOW ANALYSIS === " << endl; + AliAnalysisDataContainer *flowEvent = mgr->CreateContainer(Form("FC%s", suffixName.Data()), AliFlowEventSimple::Class(), AliAnalysisManager::kExchangeContainer); + mgr->ConnectOutput(task, 2, flowEvent); + if(debug) cout << " --> Created IO containers " << flowEvent << endl; + if(debug) cout << " --> suffixName " << suffixName << endl; + if (QC) { // add qc tasks + AddQCmethod(Form("QCTPCMB_%d_%s", mb, suffixName.Data()), harm, flowEvent, HighPtGenericFilter, debug, 0x0, suffixName.Data()); + if(debug) cout << " --> Hanging QC task ... " << mb << " succes! "<< endl; + } + if (SPSUB) { // add sp subevent tasks + AddSPmethod(Form("SPTPCMB_%d_%s", mb, suffixName.Data()), -0.8, -.5*EtaGap, .5*EtaGap, +0.8, "Qa", harm, flowEvent, HighPtSubEventFilterB, NULL, false, shrinkSP, debug, true, suffixName.Data()); + if(debug) cout << " --> Hanging SP Qa task " << mb << " succes!" << endl; + AddSPmethod(Form("SPTPCMB_%d_%s", mb, suffixName.Data()), -0.8, -.5*EtaGap, .5*EtaGap, +0.8, "Qb", harm, flowEvent, HighPtSubEventFilterA, NULL, false, shrinkSP, debug, true, suffixName.Data()); + if(debug) cout << " --> Hanging SP Qb task ..." << mb << " succes!"<< endl; + } + if (SP) { // add sp tasks + AddSPmethod(Form("SPTPCMB_%d_%s", mb, suffixName.Data()), -0.8, -0.0, +0.0, +0.8, "QaQb", harm, flowEvent, HighPtGenericFilter, NULL, false, shrinkSP, debug, 0x0, suffixName.Data()); + if(debug) cout << " --> Hanging SP task ... " << mb << " succes!" << endl; + } + if (EP) { // add ep tasks + AddSPmethod(Form("EPTPCMB_%d_%s", mb, suffixName.Data()), -0.8, -0.0, +0.0, +0.8, "QaQb", harm, flowEvent, HighPtGenericFilter, NULL, true, shrinkSP, debug, 0x0, suffixName.Data()); + if(debug) cout << " --> Hanging EP task ... " << mb << " succes!" << endl; + } + } + // print summary to screen + cout << endl << endl << " ==== AddTaskPhiFlow launched ===== " << endl; + cout << " ************ Configured PID routine ************ " << endl; + cout << " 0 < " << PIDconfig[1] << " p_t, ITS || TPC with s < " << PIDconfig[0] << endl; + if(PIDconfig[2] < 0.) cout << " --> TPC control disabled " << endl; + if(PIDconfig[2] > 0.) cout << " --> TPC control enabled " << endl; + cout << " " << PIDconfig[1] << " < " << PIDconfig[4] << " p_t, TPC || ITS with s < " << PIDconfig[3] << endl; + if(PIDconfig[5] < 0.) cout << " --> ITS control disabled " << endl; + if(PIDconfig[5] > 0.) cout << " --> ITS control enabled " << endl; + cout << " " << PIDconfig[4] << " < 7 p_t, TPC / TOF Bayesian with p < " << PIDconfig[6] << endl; + cout << " ************ Configured DCA setup ************** " << endl; + cout << " DCA type: " << DCA; + if (DCA == "") cout << "default"; + cout << endl << " ************* Task statisti:q!cs ****************** " << endl; + cout << " -> Launched PHI reconstruction " << endl; + if(SP) cout << " --> Launched QaQb SP filters and corresponding SP task " << endl; + if(EP) cout << " --> Launched QaQb QC filters and corresponding EP task " << endl; + if(SPSUB) cout << " --> Launched Qa&&Qb SP filters and corresponding SP tasks " << endl; + if(QC) cout << " --> Launched QaQb QC filters and corresponding QC task " << endl; + if(EP3sub) cout << " --> Launched VZERO subevent analysis alongside reconstruction - USE WITH CAUTION!" << endl; + cout << " ************************************************ " << endl; + TString condit = ""; + (task->SetQA(kFALSE)) ? condit+= " --> Enabled QA plots <-- " : condit+= " --> Disabled QA plots <-- "; + (task->SetIsMC(kFALSE)) ? condit+= " --> MC mode <-- " : condit+= " --> DATA mode <-- "; + (task->UseEventMixing(event_mixing)) ? condit+= " --> Using EVENT MIXING <--" : condit+= "--> Combinatorial background <--"; + cout << condit << endl; + cout << " --> Now go for a coffee! <-- " << endl; + cout << " ************************************************ " << endl; + return task; + } // end of high pt loop - high-pt task is ready to go at this point + Float_t _pt[] = {0.0, 0.3, 0.6, 0.9, 1.2, 1.5, 1.8, 2.1, 2.4, 2.7, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0}; + task->SetPtBins(_pt, (Int_t)(sizeof(_pt)/sizeof(_pt[1]))-1); + // POI filter cuts, will filter invm mass bands and subevents + AliFlowTrackSimpleCuts* POIfilterQC[30]; + AliFlowTrackSimpleCuts* POIfilterSP[30][2]; + Double_t flowBands[2][30]; + Double_t _inc = (maxMaxx-minMass)/30.; + for (Int_t mb = 0; mb < 30; mb++) { + flowBands[0][mb] = minMass + mb * _inc; + flowBands[1][mb] = minMass + (mb + 1) * _inc; + POIfilterSP[mb][0] = new AliFlowTrackSimpleCuts(Form("FilterPOISP_MB%d_ETANEG", mb)); + if(!POIfilterSP[mb][0]) { + if(debug) cout << " FAILED to initialize POI filters, could be a library problem!" << endl; + return 0x0; + } + POIfilterSP[mb][0]->SetEtaMin(-0.8); + POIfilterSP[mb][0]->SetEtaMax(0.0); + POIfilterSP[mb][0]->SetMassMin(flowBands[0][mb]); + POIfilterSP[mb][0]->SetMassMax(flowBands[1][mb]); + POIfilterSP[mb][1] = new AliFlowTrackSimpleCuts(Form("FilterPOISP_MB%d_ETAPOS", mb)); + POIfilterSP[mb][1]->SetEtaMin(0.0); + POIfilterSP[mb][1]->SetEtaMax(+0.8); + POIfilterSP[mb][1]->SetMassMin(flowBands[0][mb]); + POIfilterSP[mb][1]->SetMassMax(flowBands[1][mb]); + POIfilterQC[mb] = new AliFlowTrackSimpleCuts(Form("FilterPOIQC_MB%d", mb)); + POIfilterQC[mb]->SetEtaMin(-0.8); + POIfilterQC[mb]->SetEtaMax(+0.8); + POIfilterQC[mb]->SetMassMin(flowBands[0][mb]); + POIfilterQC[mb]->SetMassMax(flowBands[1][mb]); + } + if(debug) cout << " --> Created poi filters " << endl; + task->SetCommonConstants(30, flowBands[0][0], flowBands[1][29]); + if(debug) cout << " --> Set common constants " << endl; + // set pair and event cuts + if((deltaDip>0.005)&&(deltaDipMaxPt>0.005)) task->SetMaxDeltaDipAngleAndPt(deltaDip, deltaDipMaxPt); + else cout << " --> Disabled Delta-Dip exclusion. <-- " << endl; + task->SetCandidateEtaAndPt(POIEtaMin, POIEtaMax, 0., 10.); + task->SetCentralityParameters(centrMin, centrMax, "TRK"); + task->SetVertexZ(vertexZ); + if(debug) cout << " --> Set pair cuts and event cuts" << endl; + // set the kaon cuts, and specify the PID procedure which will be used + task->SetPIDConfiguration(PIDconfig); + if(debug) { + cout << " ************ Configured PID routine ************ " << endl; + cout << " 0 < " << PIDconfig[1] << " p_t, ITS || TPC with s < " << PIDconfig[0] << endl; + if(PIDconfig[2] < 0.) cout << " --> TPC control disabled " << endl; + if(PIDconfig[2] > 0.) cout << " --> TPC control enabled " << endl; + cout << " " << PIDconfig[1] << " < " << PIDconfig[4] << " p_t, TPC || ITS with s < " << PIDconfig[3] << endl; + if(PIDconfig[5] < 0.) cout << " --> ITS control disabled " << endl; + if(PIDconfig[5] > 0.) cout << " --> ITS control enabled " << endl; + cout << " " << PIDconfig[4] << " < 7 p_t, TPC / TOF Bayesian with p < " << PIDconfig[6] << endl; + cout << " ************************************************ " << endl; + } + if (TPCStandAloneTracks) + { // switch to tpc standalone tracks for POI selection + if(debug) cout << " --> Switching to TPC standalone analsis " << endl; + task->SetRequireStrictKaonCuts(); + task->SetRequireTPCStandAloneKaons(); + } + // create and configure IO containers + AliAnalysisDataContainer *cinput = mgr->GetCommonInputContainer(); + AliAnalysisDataContainer *coutHist = mgr->CreateContainer(Form("PhiV2_%s", OutputName(centrMin, centrMax,PIDconfig, suffixName, bCentralTrigger, EtaGap, POIEtaMin, POIEtaMax, 0., 10., deltaDip, deltaDipMaxPt, DCA, harm, TPCStandAloneTracks, vertexZ, debug, useGlobalRPCuts).Data()), TList::Class(), AliAnalysisManager::kOutputContainer, fileName); + if(debug) cout << " --> Created IO containers " << cinput << ", " << coutHist << endl; + mgr->AddTask(task); + if(debug) cout << " === ADDING MAIN TASK == " << endl; + mgr->ConnectInput(task, 0, cinput); + mgr->ConnectOutput(task, 1, coutHist); + if(debug) cout << " --> Connected IO containers " << endl; + if (SP || EP || QC || SPSUB) // if flow analysis should be done after reconstruction + { + if(debug) cout << " === RECEIVED REQUEST FOR FLOW ANALYSIS === " << endl; + AliAnalysisDataContainer *flowEvent = mgr->CreateContainer(Form("FC%s", suffixName.Data()), AliFlowEventSimple::Class(), AliAnalysisManager::kExchangeContainer); + mgr->ConnectOutput(task, 2, flowEvent); + if(debug) cout << " --> Created IO containers " << flowEvent << endl; + if(debug) cout << " --> suffixName " << suffixName << endl; + for (int mb = 0; mb != 30; ++mb) { + if (QC) { // add qc tasks + AddQCmethod(Form("QCTPCMB_%d_%s", mb, suffixName.Data()), harm, flowEvent, POIfilterQC[mb], debug, 0x0, suffixName.Data()); + if(debug) cout << " --> Hanging QC task ... " << mb << " succes! "<< endl; + } + if (SPSUB) { // add sp subevent tasks + AddSPmethod(Form("SPTPCMB_%d_%s", mb, suffixName.Data()), -0.8, -.5*EtaGap, .5*EtaGap, +0.8, "Qa", harm, flowEvent, POIfilterSP[mb][1], NULL, false, shrinkSP, debug, true, suffixName.Data(), VZERO_SP); + if(debug) cout << " --> Hanging SP Qa task " << mb << " succes!" << endl; + AddSPmethod(Form("SPTPCMB_%d_%s", mb, suffixName.Data()), -0.8, -.5*EtaGap, .5*EtaGap, +0.8, "Qb", harm, flowEvent, POIfilterSP[mb][0], NULL, false, shrinkSP, debug, true, suffixName.Data(), VZERO_SP); + if(debug) cout << " --> Hanging SP Qb task ..." << mb << " succes!"<< endl; + } + if (SP) { // add sp tasks + AddSPmethod(Form("SPTPCMB_%d_%s", mb, suffixName.Data()), -0.8, -0.0, +0.0, +0.8, "QaQb", harm, flowEvent, POIfilterQC[mb], NULL, false, shrinkSP, debug, 0x0, suffixName.Data()); + if(debug) cout << " --> Hanging SP task ... " << mb << " succes!" << endl; + } + if (EP) { // add ep tasks + AddSPmethod(Form("EPTPCMB_%d_%s", mb, suffixName.Data()), -0.8, -0.0, +0.0, +0.8, "QaQb", harm, flowEvent, POIfilterQC[mb], NULL, true, shrinkSP, debug, 0x0, suffixName.Data()); + if(debug) cout << " --> Hanging EP task ... " << mb << " succes!" << endl; + } + } + } + // print summary to screen + cout << endl << endl << " ==== AddTaskPhiFlow launched ===== " << endl; + cout << " ************ Configured PID routine ************ " << endl; + cout << " 0 < " << PIDconfig[1] << " p_t, ITS || TPC with s < " << PIDconfig[0] << endl; + if(PIDconfig[2] < 0.) cout << " --> TPC control disabled " << endl; + if(PIDconfig[2] > 0.) cout << " --> TPC control enabled " << endl; + cout << " " << PIDconfig[1] << " < " << PIDconfig[4] << " p_t, TPC || ITS with s < " << PIDconfig[3] << endl; + if(PIDconfig[5] < 0.) cout << " --> ITS control disabled " << endl; + if(PIDconfig[5] > 0.) cout << " --> ITS control enabled " << endl; + cout << " " << PIDconfig[4] << " < 7 p_t, TPC / TOF Bayesian with p < " << PIDconfig[6] << endl; + cout << " ************ Configured DCA setup ************** " << endl; + cout << " DCA type: " << DCA; + if (DCA == "") cout << "default"; + cout << endl << " ************* Task statisti:q!cs ****************** " << endl; + cout << " -> Launched PHI reconstruction " << endl; + if(SP) cout << " --> Launched 30 QaQb SP filters and corresponding 30 SP tasks " << endl; + if(EP) cout << " --> Launched 30 QaQb QC filters and corresponding 30 EP tasks " << endl; + if(SPSUB) cout << " --> Launched 30+30 Qa&&Qb SP filters and corresponding 30+30 SP tasks " << endl; + if(QC) cout << " --> Launched 30 QaQb QC filters and corresponding 30 QC tasks " << endl; + if(EP3sub) cout << " --> Launched VZERO subevent analysis alongside reconstruction - USE WITH CAUTION!" << endl; + cout << " ************************************************ " << endl; + TString condit = ""; + (task->SetQA(kFALSE)) ? condit+= " --> Enabled QA plots <-- " : condit+= " --> Disabled QA plots <-- "; + (task->SetIsMC(kFALSE)) ? condit+= " --> MC mode <-- " : condit+= " --> DATA mode <-- "; + (task->UseEventMixing(event_mixing)) ? condit+= " --> Using EVENT MIXING <--" : condit+= "--> Combinatorial background <--"; + cout << condit << endl; + cout << " --> Now go for a coffee! <-- " << endl; + cout << " ************************************************ " << endl; + return task; +} +//_____________________________________________________________________________ +void AddSPmethod(char *name, double minEtaA, double maxEtaA, double minEtaB, double maxEtaB, char *Qvector, int harmonic, AliAnalysisDataContainer *flowEvent, AliFlowTrackSimpleCuts *cutsPOI = NULL, AliFlowTrackSimpleCuts *cutsRFP = NULL, bool bEP, bool shrink = false, bool debug, bool etagap = false, char *suffixName, Bool_t VZERO_SP = kFALSE) +{ + // add sp task and invm filter tasks + if(debug) (bEP) ? cout << " ****** Reveived request for EP task ****** " << endl : cout << " ******* Switching to SP task ******* " << endl; + TString fileName = AliAnalysisManager::GetCommonFileName(); + (bEP) ? fileName+=":EP" : fileName+=":SP"; + fileName+=suffixName; + if(etagap) { + fileName+="_SUBEVENTS"; + if(debug) cout << " --> Setting up subevent analysis <-- " << endl; + } + if(debug) cout << " --> fileName " << fileName << endl; + TString myNameSP; + (bEP) ? myNameSP = Form("%sEPv%d%s", name, harmonic, Qvector): myNameSP = Form("%sSPv%d%s", name, harmonic, Qvector); + if(debug) cout << " Task and filter name: " << myNameSP << endl; + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + AliAnalysisDataContainer *flowEvent2 = mgr->CreateContainer(Form("Filter_%s", myNameSP.Data()), AliFlowEventSimple::Class(), AliAnalysisManager::kExchangeContainer); + AliAnalysisTaskFilterFE *tskFilter = new AliAnalysisTaskFilterFE(Form("TaskFilter_%s", myNameSP.Data()), cutsRFP, cutsPOI); + tskFilter->SetSubeventEtaRange(minEtaA, maxEtaA, minEtaB, maxEtaB); + if(VZERO_SP) tskFilter->SetSubeventEtaRange(-10, 0, 0, 10); + mgr->AddTask(tskFilter); + mgr->ConnectInput(tskFilter, 0, flowEvent); + mgr->ConnectOutput(tskFilter, 1, flowEvent2); + AliAnalysisDataContainer *outSP = mgr->CreateContainer(myNameSP.Data(), TList::Class(), AliAnalysisManager::kOutputContainer, fileName); + AliAnalysisTaskScalarProduct *tskSP = new AliAnalysisTaskScalarProduct(Form("TaskScalarProduct_%s", myNameSP.Data()), kFALSE); + tskSP->SetApplyCorrectionForNUA(kTRUE); + tskSP->SetHarmonic(harmonic); + tskSP->SetTotalQvector(Qvector); + if (bEP) tskSP->SetBehaveAsEP(); + if (shrink) tskSP->SetBookOnlyBasicCCH(kTRUE); + mgr->AddTask(tskSP); + mgr->ConnectInput(tskSP, 0, flowEvent2); + mgr->ConnectOutput(tskSP, 1, outSP); +} +//_____________________________________________________________________________ +void AddQCmethod(char *name, int harmonic, AliAnalysisDataContainer *flowEvent, AliFlowTrackSimpleCuts *cutsPOI = NULL, Bool_t debug, AliFlowTrackSimpleCuts *cutsRFP = NULL, char *suffixName) +{ + // add qc task and invm filter tasks + if(debug) cout << " ****** Received request for QC v" << harmonic << " task " << name << ", POI " << cutsPOI << ", IO ****** " << flowEvent << endl; + TString fileName = AliAnalysisManager::GetCommonFileName(); + fileName+=":QC"; + fileName+=suffixName; + if(debug) cout << " --> Common filename: " << fileName << endl; + TString myName = Form("%s", name); + if(debug) cout << " --> myName: " << myName << endl; + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + AliAnalysisDataContainer *flowEvent2 = mgr->CreateContainer(Form("Filter_%s", myName.Data()), AliFlowEventSimple::Class(),AliAnalysisManager::kExchangeContainer); + AliAnalysisTaskFilterFE *tskFilter = new AliAnalysisTaskFilterFE(Form("TaskFilter_%s", myName.Data()), cutsRFP, cutsPOI); + mgr->AddTask(tskFilter); + mgr->ConnectInput(tskFilter, 0, flowEvent); + mgr->ConnectOutput(tskFilter, 1, flowEvent2); + AliAnalysisDataContainer *outQC = mgr->CreateContainer(myName.Data(), TList::Class(), AliAnalysisManager::kOutputContainer, fileName); + AliAnalysisTaskQCumulants *tskQC = new AliAnalysisTaskQCumulants(Form("TaskQCumulants_%s", myName.Data()), kFALSE); + tskQC->SetApplyCorrectionForNUA(kTRUE); + tskQC->SetHarmonic(harmonic); + tskQC->SetBookOnlyBasicCCH(kTRUE); + mgr->AddTask(tskQC); + mgr->ConnectInput(tskQC, 0, flowEvent2); + mgr->ConnectOutput(tskQC, 1, outQC); +} +//_____________________________________________________________________________ +TString OutputName( Float_t centrMin, + Float_t centrMax, + Float_t PIDconfig[7], + TString suffixName, + Bool_t bCentralTrigger, + Float_t EtaGap, + Float_t POIEtaMin, + Float_t POIEtaMax, + Float_t POIPtMin, + Float_t POIPtMax, + Float_t deltaDip, + Float_t deltaDipMaxPt, + TString DCA, + Int_t harm, + Bool_t TPCStandAloneTracks, + Float_t vertexZ, + Bool_t debug, + Bool_t useGlobalRPCuts) +{ + // generate output name + TString centralityName = ""; + centralityName += suffixName; + centralityName += "_DCA"; + centralityName += DCA; + centralityName += Form("_vZ%.f", vertexZ); + centralityName += "_"; + for(Int_t i = 0; i < 7; i++) centralityName += Form("%.1f_", PIDconfig[i]); + centralityName += "_POIEta"; + centralityName += Form("%.1f", POIEtaMin); + centralityName += Form("%.1f", POIEtaMax); + centralityName += "_gap"; + centralityName += Form("%.1f", -0.5*EtaGap); + centralityName += Form("%.1f", 0.5*EtaGap); + centralityName += "_"; + centralityName += Form("dDip%.2f", deltaDip); + centralityName += "-"; + centralityName += Form("dDipPt%.2f", deltaDipMaxPt); + if (TPCStandAloneTracks) { + centralityName += "-"; + centralityName += "TPCStandAloneTracks"; + } + if (bCentralTrigger) { + centralityName += "-"; + centralityName += "kMBkCkSC"; + } + if (!useGlobalRPCuts) { + centralityName += "-"; + centralityName += "TPCRP"; + } + if(debug) cout << " --> centralityName " << centralityName << endl; + return centralityName; +} +//_____________________________________________________________________________ +