From 2b5564407dd402a89b7b870a1f0fc1a5d6dab63a Mon Sep 17 00:00:00 2001 From: hansena Date: Tue, 14 Feb 2012 12:03:58 +0000 Subject: [PATCH] Mega-commit by Alexander - added Event plane to analysis task, and other flow updates. Some clean-up. --- PWGLF/CMakelibPWGLFforward2.pkg | 5 +- PWGLF/FORWARD/analysis2/AddTaskForwardFlow.C | 48 +- .../analysis2/AddTaskMCParticleFilter.C | 20 + PWGLF/FORWARD/analysis2/AliAODForwardEP.cxx | 124 ++ PWGLF/FORWARD/analysis2/AliAODForwardEP.h | 165 +++ .../analysis2/AliCentralMultiplicityTask.cxx | 2 +- .../analysis2/AliFMDEventPlaneFinder.cxx | 555 ++++++++ .../analysis2/AliFMDEventPlaneFinder.h | 186 +++ .../analysis2/AliFMDMCTrackDensity.cxx | 145 +- .../FORWARD/analysis2/AliFMDMCTrackDensity.h | 18 +- .../analysis2/AliForwardFlowTaskQC.cxx | 1257 ++++++++++------- .../FORWARD/analysis2/AliForwardFlowTaskQC.h | 256 ++-- .../FORWARD/analysis2/AliForwardFlowUtil.cxx | 433 ------ PWGLF/FORWARD/analysis2/AliForwardFlowUtil.h | 176 --- .../analysis2/AliForwardMCFlowTaskQC.cxx | 441 ++++++ .../analysis2/AliForwardMCFlowTaskQC.h | 159 +++ .../AliForwardMCMultiplicityTask.cxx | 21 + .../analysis2/AliForwardMCMultiplicityTask.h | 16 + .../analysis2/AliForwardMultiplicityBase.cxx | 3 + .../analysis2/AliForwardMultiplicityBase.h | 14 + .../analysis2/AliForwardMultiplicityTask.cxx | 24 +- .../analysis2/AliForwardMultiplicityTask.h | 17 + .../analysis2/AliSPDMCTrackDensity.cxx | 134 +- .../FORWARD/analysis2/AliSPDMCTrackDensity.h | 18 +- PWGLF/FORWARD/analysis2/ForwardAODConfig.C | 3 + PWGLF/FORWARD/analysis2/MakeFlow.C | 92 +- PWGLF/FORWARD/analysis2/trains/MakeAODTrain.C | 47 +- .../FORWARD/analysis2/trains/MakeFlowTrain.C | 124 ++ PWGLF/PWGLFforward2LinkDef.h | 7 +- 29 files changed, 3216 insertions(+), 1294 deletions(-) create mode 100644 PWGLF/FORWARD/analysis2/AddTaskMCParticleFilter.C create mode 100644 PWGLF/FORWARD/analysis2/AliAODForwardEP.cxx create mode 100644 PWGLF/FORWARD/analysis2/AliAODForwardEP.h create mode 100644 PWGLF/FORWARD/analysis2/AliFMDEventPlaneFinder.cxx create mode 100644 PWGLF/FORWARD/analysis2/AliFMDEventPlaneFinder.h delete mode 100644 PWGLF/FORWARD/analysis2/AliForwardFlowUtil.cxx delete mode 100644 PWGLF/FORWARD/analysis2/AliForwardFlowUtil.h create mode 100644 PWGLF/FORWARD/analysis2/AliForwardMCFlowTaskQC.cxx create mode 100644 PWGLF/FORWARD/analysis2/AliForwardMCFlowTaskQC.h create mode 100644 PWGLF/FORWARD/analysis2/trains/MakeFlowTrain.C diff --git a/PWGLF/CMakelibPWGLFforward2.pkg b/PWGLF/CMakelibPWGLFforward2.pkg index 398357d9cb4..7589ea8a40d 100644 --- a/PWGLF/CMakelibPWGLFforward2.pkg +++ b/PWGLF/CMakelibPWGLFforward2.pkg @@ -39,6 +39,7 @@ set ( SRCS FORWARD/analysis2/AliFMDEnergyFitter.cxx FORWARD/analysis2/AliFMDEnergyFitterTask.cxx FORWARD/analysis2/AliFMDEventInspector.cxx + FORWARD/analysis2/AliFMDEventPlaneFinder.cxx FORWARD/analysis2/AliFMDHistCollector.cxx FORWARD/analysis2/AliFMDSharingFilter.cxx FORWARD/analysis2/AliFMDMCEventInspector.cxx @@ -57,13 +58,14 @@ set ( SRCS FORWARD/analysis2/AliCentralMCMultiplicityTask.cxx FORWARD/analysis2/AliCentralMCCorrectionsTask.cxx FORWARD/analysis2/AliAODCentralMult.cxx + FORWARD/analysis2/AliAODForwardEP.cxx FORWARD/analysis2/AliCentralCorrSecondaryMap.cxx FORWARD/analysis2/AliCentralCorrAcceptance.cxx FORWARD/analysis2/AliCentraldNdetaTask.cxx FORWARD/analysis2/AliBasedNdetaTask.cxx FORWARD/analysis2/AliMCTruthdNdetaTask.cxx - FORWARD/analysis2/AliForwardFlowUtil.cxx FORWARD/analysis2/AliForwardFlowTaskQC.cxx + FORWARD/analysis2/AliForwardMCFlowTaskQC.cxx FORWARD/analysis2/AliSPDMCTrackDensity.cxx FORWARD/analysis2/AliFMDMultCuts.cxx FORWARD/analysis2/AliPoissonCalculator.cxx @@ -159,6 +161,7 @@ install ( FILES FORWARD/analysis2/AddTaskCentraldNdeta.C FORWARD/analysis2/AddTaskForwardQA.C FORWARD/analysis2/AddTaskForwarddNdeta.C FORWARD/analysis2/AddTaskMCTruthdNdeta.C + FORWARD/analysis2/AddTaskMCParticleFilter.C FORWARD/analysis2/DrawdNdeta.C FORWARD/analysis2/ForwardAODConfig.C FORWARD/analysis2/CentralAODConfig.C diff --git a/PWGLF/FORWARD/analysis2/AddTaskForwardFlow.C b/PWGLF/FORWARD/analysis2/AddTaskForwardFlow.C index 5fc568f385e..9767dbb82a9 100644 --- a/PWGLF/FORWARD/analysis2/AddTaskForwardFlow.C +++ b/PWGLF/FORWARD/analysis2/AddTaskForwardFlow.C @@ -1,7 +1,7 @@ /** * @file AddTaskForwardFlow.C - * @author Christian Holm Christensen - * @date Wed Mar 23 12:14:17 2011 + * @author Alexander Hansen alexander.hansen@cern.ch + * @date Wed Sep 07 12:14:17 2011 * * @brief * @@ -25,7 +25,6 @@ * @ingroup pwglf_forward_flow */ void AddTaskForwardFlow(TString type = "", - Int_t etabins = 48, Bool_t mc = kFALSE, TString addFlow = "", Int_t addFType = 0, @@ -64,27 +63,36 @@ void AddTaskForwardFlow(TString type = "", if (!type.Contains("6")) v6 = kFALSE; } - // --- Create output containers and find input from fmd task --- // - - TString outputFile = AliAnalysisManager::GetCommonFileName(); - outputFile += ":FlowResults"; - - AliAnalysisDataContainer* qcout = mgr->CreateContainer("QCumulants", TList::Class(), AliAnalysisManager::kOutputContainer, outputFile); + // --- Create containers for output --- // + AliAnalysisDataContainer* sums = + mgr->CreateContainer("FlowQCSums", TList::Class(), + AliAnalysisManager::kOutputContainer, + AliAnalysisManager::GetCommonFileName()); + AliAnalysisDataContainer* output = + mgr->CreateContainer("FlowQCResults", TList::Class(), + AliAnalysisManager::kParamContainer, + AliAnalysisManager::GetCommonFileName()); // --- For the selected flow tasks the input and output is set --- // - AliForwardFlowTaskQC* qc = new AliForwardFlowTaskQC("QCumulants"); + AliForwardFlowTaskQC* task = 0; + if (mc) + task = new AliForwardMCFlowTaskQC("QCumulants"); + else + task = new AliForwardFlowTaskQC("QCumulants"); + mgr->AddTask(task); - qc->SetDoHarmonics(v1, v2, v3, v4, v5, v6); - qc->SetUseNEtaBins(etabins); - qc->SetMCinput(mc); - qc->AddFlow(addFlow); - qc->AddFlowType(addFType); - qc->AddFlowOrder(addFOrder); - qc->SetRefEtaBins(4); - - mgr->ConnectInput(qc, 0, mgr->GetCommonInputContainer()); - mgr->ConnectOutput(qc, 1, qcout); + task->SetDoHarmonics(v1, v2, v3, v4, v5, v6); + if (mc) { + AliForwardMCFlowTaskQC* mcTask = + static_casttask; + mcTask->AddFlow(addFlow); + mcTask->AddFlowType(addFType); + mcTask->AddFlowOrder(addFOrder); + } + mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer()); + mgr->ConnectOutput(task, 1, sums); + mgr->ConnectOutput(task, 2, output); return; } diff --git a/PWGLF/FORWARD/analysis2/AddTaskMCParticleFilter.C b/PWGLF/FORWARD/analysis2/AddTaskMCParticleFilter.C new file mode 100644 index 00000000000..104f72d29e7 --- /dev/null +++ b/PWGLF/FORWARD/analysis2/AddTaskMCParticleFilter.C @@ -0,0 +1,20 @@ +AliAnalysisTask* +AddTaskMCParticleFilter() +{ + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if (!mgr) { + Error("AddTaskFMD", "No analysis manager to connect to."); + return NULL; + } + + AliAnalysisTaskMCParticleFilter* mctask = new AliAnalysisTaskMCParticleFilter("mcfilter"); + mgr->AddTask(mctask); + AliAnalysisDataContainer* histOut = + mgr->CreateContainer("mcfilter", TList::Class(), + AliAnalysisManager::kOutputContainer, "mcParticles.root"); + mgr->ConnectInput(mctask, 0, mgr->GetCommonInputContainer()); + mgr->ConnectOutput(mctask, 0, mgr->GetCommonOutputContainer()); + mgr->ConnectOutput(mctask, 1, histOut); + + return mctask; +} diff --git a/PWGLF/FORWARD/analysis2/AliAODForwardEP.cxx b/PWGLF/FORWARD/analysis2/AliAODForwardEP.cxx new file mode 100644 index 00000000000..ab50c8331d7 --- /dev/null +++ b/PWGLF/FORWARD/analysis2/AliAODForwardEP.cxx @@ -0,0 +1,124 @@ +// +// Class that contains results from FMD eventplane calculations +// +#include "AliAODForwardEP.h" +#include +#include +#include +#include +#include +#include "AliLog.h" +ClassImp(AliAODForwardEP) +#ifdef DOXY_INPUT +; // For Emacs +#endif + +//____________________________________________________________________ +AliAODForwardEP::AliAODForwardEP() + : fIsMC(false), + fEpT(-1), + fEpA(-1), + fEpC(-1), + fEp1(-1), + fEp2(-1), + fHist() +{ + // + // Constructor + // +} + +//____________________________________________________________________ +AliAODForwardEP::AliAODForwardEP(Bool_t isMC) + : fIsMC(isMC), + fEpT(-1), + fEpA(-1), + fEpC(-1), + fEp1(-1), + fEp2(-1), + fHist() +{ + // + // Constructor + // + // Parameters: + // isMC If set to true this is for MC data (effects branch name) + // + fHist.SetXTitle("#eta"); + fHist.SetDirectory(0); + fHist.Sumw2(); +} + +//____________________________________________________________________ +void +AliAODForwardEP::Init(const TAxis& etaAxis) +{ + // Initialize the histogram with an eta axis + // + // Parameters: + // etaAxis Eta axis to use + // + fHist.SetBins(etaAxis.GetNbins()/10, etaAxis.GetXmin(), etaAxis.GetXmax()); +} + +//____________________________________________________________________ +void +AliAODForwardEP::Clear(Option_t* option) +{ + // Clear (or reset) internal values + // + // Parameters: + // option Passed to TH1::Reset + // + fHist.Reset(option); + fEpT = -1; + fEpA = -1; + fEpC = -1; + fEp1 = -1; + fEp2 = -1; +} +//____________________________________________________________________ +void +AliAODForwardEP::Browse(TBrowser* /*b*/) +{ + // Browse this object + // + // Parameters: + // b Browser to use + // TODO: Make nice +/* static TObjString ipz; + static TObjString trg; + static TObjString cnt; + static TObjString ncl; + ipz = Form("ip_z=%fcm", fIpZ); + trg = GetTriggerString(fTriggers); + cnt = Form("%+6.1f%%", fCentrality); + ncl = Form("%d clusters", fNClusters); + b->Add(&fHist); + b->Add(&ipz); + b->Add(&trg); + b->Add(&cnt); + b->Add(&ncl); +*/ +} + +//____________________________________________________________________ +void +AliAODForwardEP::Print(Option_t* option) const +{ + // Print this object + // + // Parameters: + // option Not used + fHist.Print(option); + std::cout << "Total FMD EP: \t" << fEpT < +#include +class TBrowser; +class TH1D; + +class AliAODForwardEP : public TObject +{ +public: + /** + * Default constructor + * + * Used by ROOT I/O sub-system - do not use + */ + AliAODForwardEP(); + /** + * Constructor + * + * @param isMC Whether this was from MC or not + */ + AliAODForwardEP(Bool_t isMC); + /** + * Destructor + */ + virtual ~AliAODForwardEP() {} // Destructor + /** + * Initialize + * + * @param etaAxis Pseudo-rapidity axis + */ + void Init(const TAxis& etaAxis); + /** + * Clear all data + * + * @param option Passed on to TH2::Reset verbatim + */ + void Clear(Option_t* option=""); + /** + * browse this object + * + * @param b Browser + */ + void Browse(TBrowser* b); + /** + * This is a folder + * + * @return Always true + */ + Bool_t IsFolder() const { return kTRUE; } // Always true + /** + * Print content + * + * @param option Passed verbatim to TH2::Print + */ + void Print(Option_t* option="") const; + /* + * Get the name of the AOD object + * + * @return for now ForwardEP + */ + const Char_t* GetName() const { return (fIsMC ? "ForwardEP" : "ForwardEP"); } + /* + * Set the overall FMD eventplane + * + * @param ep FMD EP + */ + void SetEventplane(Double_t ep) { fEpT = ep; } + /* + * Get the overall FMD eventplane + * + * @return FMD EP + */ + Double_t GetEventplane() { return fEpT; } + /* + * Set FMD eventplane from A side + * + * @param epA FMD A side EP + */ + void SetEventplaneA(Double_t epA) { fEpA = epA; } + /* + * Get FMD eventplane from A side + * + * @return FMD A side EP + */ + Double_t GetEventplaneA() { return fEpA; } + /* + * Set FMD eventplane from C side + * + * @param epC FMD C side EP + */ + void SetEventplaneC(Double_t epC) { fEpC = epC; } + /* + * Get FMD eventplane from C side + * + * @return FMD C side EP + */ + Double_t GetEventplaneC() { return fEpC; } + /* + * Set FMD eventplane 1 using random particles + * + * @param ep1 FMD EP 1 from random selection + */ + void SetEventplane1(Double_t ep1) { fEp1 = ep1; } + /* + * Get FMD eventplane 1 from random particles + * + * @return FMD EP 1 from random selection + */ + Double_t GetEventplane1() { return fEp1; } + /* + * Set FMD eventplane 2 using random particles + * + * @param ep2 FMD EP 2 from random selection + */ + void SetEventplane2(Double_t ep2) { fEp2 = ep2; } + /* + * Get FMD eventplane 2 from random particles + * + * @return FMD EP 2 from random selection + */ + Double_t GetEventplane2() { return fEp2; } + /* + * Get eta histogram from eta vs. Psi_2 calculatins + * + * @return Returns eta vs. Psi_2 histogram + */ + const TH1D& GetHistogram() const { return fHist; } // Get histogram + /* + * Get eta histogram from eta vs. Psi_2 calculatins + * + * @return Returns eta vs. Psi_2 histogram + */ + TH1D& GetHistogram() { return fHist; } // Get histogram + +protected: + Bool_t fIsMC; // Whether this is from MC + Double_t fEpT; // Total FMD event plane + Double_t fEpA; // FMD1+2 subevent plane + Double_t fEpC; // FMD3 subevent plane + Double_t fEp1; // Random FMD subevent plane 1 + Double_t fEp2; // Random FMD subevent plane 2 + TH1D fHist; // Histogram of subevent planes in eta for this event + + ClassDef(AliAODForwardEP,1); // AOD forward multiplicity +}; + +#endif +// Local Variables: +// mode: C++ +// End: diff --git a/PWGLF/FORWARD/analysis2/AliCentralMultiplicityTask.cxx b/PWGLF/FORWARD/analysis2/AliCentralMultiplicityTask.cxx index ba4d2058daa..38613aa3979 100644 --- a/PWGLF/FORWARD/analysis2/AliCentralMultiplicityTask.cxx +++ b/PWGLF/FORWARD/analysis2/AliCentralMultiplicityTask.cxx @@ -188,7 +188,7 @@ AliCentralMultiplicityTask::GetESDEvent() fNClusterTracklet = new TH2D("nClusterVsnTracklet", "Total number of cluster vs number of tracklets", - 100, 0, 100, 100, 0, 100); + 100, 0, 10000, 100, 0, 10000); fNClusterTracklet->SetDirectory(0); fNClusterTracklet->SetXTitle("# of free clusters"); fNClusterTracklet->SetYTitle("# of tracklets"); diff --git a/PWGLF/FORWARD/analysis2/AliFMDEventPlaneFinder.cxx b/PWGLF/FORWARD/analysis2/AliFMDEventPlaneFinder.cxx new file mode 100644 index 00000000000..31ce9cb3051 --- /dev/null +++ b/PWGLF/FORWARD/analysis2/AliFMDEventPlaneFinder.cxx @@ -0,0 +1,555 @@ +// This task finds the eventplane +// using the FMD +// +#include "AliFMDEventPlaneFinder.h" +#include +#include +#include "AliLog.h" +#include +#include +#include "TROOT.h" +#include +#include +#include +#include +#include "AliOADBContainer.h" +#include "AliAnalysisManager.h" +#include "AliVEvent.h" +#include "AliEventplane.h" +#include "AliCentrality.h" +#include "AliAODForwardEP.h" +#include "AliAODForwardMult.h" +#include "AliAODEvent.h" + +ClassImp(AliFMDEventPlaneFinder) +#if 0 +; // For Emacs +#endif + +//____________________________________________________________________ +AliFMDEventPlaneFinder::AliFMDEventPlaneFinder() + : TNamed(), + fList(0), + fEvent(0), + fQt(), + fQa(), + fQc(), + fQ1(), + fQ2(), + fQeta(), + fHistFMDEventplane(0), + fHistFMDEventplaneA(0), + fHistFMDEventplaneC(0), + fHistFMDEventplane1(0), + fHistFMDEventplane2(0), + fHistPhiDist(0), + fHistFMDmTPCep(0), + fHistFMDvsTPCep(0), + fHistFMDmVZEROep(0), + fHistFMDvsVZEROep(0), + fDebug(0), + fOADBFileName(0), + fOADBContainer(0), + fPhiDist(0), + fRunNumber(0), + fUsePhiWeights(0) +{ + // + // Constructor + // + +} + +//____________________________________________________________________ +AliFMDEventPlaneFinder::AliFMDEventPlaneFinder(const char* title) + : TNamed("fmdEventPlaneFinder", title), + fList(0), + fEvent(0), + fQt(), + fQa(), + fQc(), + fQ1(), + fQ2(), + fQeta(), + fHistFMDEventplane(0), + fHistFMDEventplaneA(0), + fHistFMDEventplaneC(0), + fHistFMDEventplane1(0), + fHistFMDEventplane2(0), + fHistPhiDist(0), + fHistFMDmTPCep(0), + fHistFMDvsTPCep(0), + fHistFMDmVZEROep(0), + fHistFMDvsVZEROep(0), + fDebug(0), + fOADBFileName(0), + fOADBContainer(0), + fPhiDist(0), + fRunNumber(0), + fUsePhiWeights(1) +{ + // + // Constructor + // + // Parameters: + // name Name of object + // +} + +//____________________________________________________________________ +AliFMDEventPlaneFinder::AliFMDEventPlaneFinder(const + AliFMDEventPlaneFinder& o) + : TNamed(o), + fList(o.fList), + fEvent(o.fEvent), + fQt(o.fQt), + fQa(o.fQa), + fQc(o.fQc), + fQ1(o.fQ1), + fQ2(o.fQ2), + fQeta(o.fQeta), + fHistFMDEventplane(o.fHistFMDEventplane), + fHistFMDEventplaneA(o.fHistFMDEventplaneA), + fHistFMDEventplaneC(o.fHistFMDEventplaneC), + fHistFMDEventplane1(o.fHistFMDEventplane1), + fHistFMDEventplane2(o.fHistFMDEventplane2), + fHistPhiDist(o.fHistPhiDist), + fHistFMDmTPCep(o.fHistFMDmTPCep), + fHistFMDvsTPCep(o.fHistFMDvsTPCep), + fHistFMDmVZEROep(o.fHistFMDmVZEROep), + fHistFMDvsVZEROep(o.fHistFMDvsVZEROep), + fDebug(o.fDebug), + fOADBFileName(o.fOADBFileName), + fOADBContainer(o.fOADBContainer), + fPhiDist(o.fPhiDist), + fRunNumber(o.fRunNumber), + fUsePhiWeights(o.fUsePhiWeights) +{ + // + // Copy constructor + // + // Parameters: + // o Object to copy from + // +} + +//____________________________________________________________________ +AliFMDEventPlaneFinder::~AliFMDEventPlaneFinder() +{ + // + // Destructor + // +} + +//____________________________________________________________________ +AliFMDEventPlaneFinder& +AliFMDEventPlaneFinder::operator=(const AliFMDEventPlaneFinder& o) +{ + // + // Assignement operator + // + // Parameters: + // o Object to assign from + // + // Return: + // Reference to this object + // + if (&o == this) return *this; + TNamed::operator=(o); + + fList = o.fList; + fEvent = o.fEvent; + fQt = o.fQt; + fQa = o.fQa; + fQc = o.fQc; + fQ1 = o.fQ1; + fQ2 = o.fQ2; + fQeta = o.fQeta; + fHistFMDEventplane = o.fHistFMDEventplane; + fHistFMDEventplaneA = o.fHistFMDEventplaneA; + fHistFMDEventplaneC = o.fHistFMDEventplaneC; + fHistFMDEventplane1 = o.fHistFMDEventplane1; + fHistFMDEventplane2 = o.fHistFMDEventplane2; + fHistPhiDist = o.fHistPhiDist; + fHistFMDmTPCep = o.fHistFMDmTPCep; + fHistFMDvsTPCep = o.fHistFMDvsTPCep; + fHistFMDmVZEROep = o.fHistFMDmVZEROep; + fHistFMDvsVZEROep = o.fHistFMDvsVZEROep; + fDebug = o.fDebug; + fOADBFileName = o.fOADBFileName; + fOADBContainer = o.fOADBContainer; + fPhiDist = o.fPhiDist; + fRunNumber = o.fRunNumber; + fUsePhiWeights = o.fUsePhiWeights; + + return *this; +} + +//____________________________________________________________________ +void +AliFMDEventPlaneFinder::Init(const TAxis& etaAxis) +{ + // Intialize this sub-algorithm + // + // Parameters: + // etaAxis fmd eta axis binning + // + fHistFMDEventplane = new TH1D("hFMDEventplane", "FMD eventplane", + 100, 0., TMath::Pi()); + fHistFMDEventplane->Sumw2(); + fHistFMDEventplane->SetDirectory(0); + fList->Add(fHistFMDEventplane); + + // + fHistFMDEventplaneA = new TH1D("hFMDEventplaneA", "FMD eventplaneA", + 100, 0., TMath::Pi()); + fHistFMDEventplaneA->Sumw2(); + fHistFMDEventplaneA->SetDirectory(0); + fList->Add(fHistFMDEventplaneA); + + // + fHistFMDEventplaneC = new TH1D("hFMDEventplaneC", "FMD eventplaneC", + 100, 0., TMath::Pi()); + fHistFMDEventplaneC->Sumw2(); + fHistFMDEventplaneC->SetDirectory(0); + fList->Add(fHistFMDEventplaneC); + + // + fHistFMDEventplane1 = new TH1D("hFMDEventplane1", "FMD eventplane1", + 100, 0., TMath::Pi()); + fHistFMDEventplane1->Sumw2(); + fHistFMDEventplane1->SetDirectory(0); + fList->Add(fHistFMDEventplane1); + + // + fHistFMDEventplane2 = new TH1D("hFMDEventplane2", "FMD eventplane2", + 100, 0., TMath::Pi()); + fHistFMDEventplane2->Sumw2(); + fHistFMDEventplane2->SetDirectory(0); + fList->Add(fHistFMDEventplane2); + + // + fHistPhiDist = new TH2D("hPhiDist", "Phi distribution in FMD", + etaAxis.GetNbins(), etaAxis.GetXmin(), etaAxis.GetXmax(), + 20, 0., TMath::TwoPi()); + fHistPhiDist->Sumw2(); + fHistPhiDist->SetDirectory(0); + fList->Add(fHistPhiDist); + + // + fHistFMDmTPCep = new TH1D("hFMDmTPCep", + "Eventplane from FMD - TPC", + 100, -TMath::Pi()/2., TMath::Pi()/2.); + fHistFMDmTPCep->Sumw2(); + fHistFMDmTPCep->SetDirectory(0); + fList->Add(fHistFMDmTPCep); + + // + fHistFMDvsTPCep = new TH2F("hFMDvsTPCep", + "Eventplane from FMD vs. TPC", + 100, 0., TMath::Pi(), + 100, 0., TMath::Pi()); + fHistFMDvsTPCep->Sumw2(); + fHistFMDvsTPCep->SetDirectory(0); + fList->Add(fHistFMDvsTPCep); + + // + fHistFMDmVZEROep = new TH1D("hFMDmVZEROep", + "Eventplane from FMD - VZERO", + 100, -TMath::Pi()/2., TMath::Pi()/2.); + fHistFMDmVZEROep->Sumw2(); + fHistFMDmVZEROep->SetDirectory(0); + fList->Add(fHistFMDmVZEROep); + + // + fHistFMDvsVZEROep = new TH2F("hFMDvsVZEROep", + "Eventplane from FMD vs. VZERO", + 100, 0., TMath::Pi(), + 100, 0., TMath::Pi()); + fHistFMDvsVZEROep->Sumw2(); + fHistFMDvsVZEROep->SetDirectory(0); + fList->Add(fHistFMDvsVZEROep); + + if (!fUsePhiWeights) return; + + if (fOADBFileName.Length()==0) + fOADBFileName = Form("%s/PWGLF/FORWARD/FMDEVENTPLANE/data/fmdEPoadb.root", + AliAnalysisManager::GetOADBPath()); + + TFile foadb(fOADBFileName); + if(!foadb.IsOpen()) { + AliError(Form("Cannot open OADB file %s", fOADBFileName.Data())); + return; + } + fOADBContainer = (AliOADBContainer*)foadb.Get("FMDphidist"); + if (!fOADBContainer) AliError(Form("No OADB container found in %s", fOADBFileName.Data())); + foadb.Close(); + + return; +} + +//_____________________________________________________________________ +void +AliFMDEventPlaneFinder::DefineOutput(TList* dir) +{ + // + // Output diagnostic histograms to directory + // + // Parameters: + // dir List to write in + // + fList = new TList(); + fList->SetOwner(); + fList->SetName(GetName()); + dir->Add(fList); + + return; +} + +//____________________________________________________________________ +Bool_t +AliFMDEventPlaneFinder::FindEventplane(AliVEvent* ev, + AliAODForwardEP& aodep, + TH2D* h, + AliForwardUtil::Histos* hists) +{ + // + // Do the calculations + // + // Parameters: + // hists Histogram cache + // ep calculated eventplane + // + // Return: + // true on successs + + fQt.Set(0., 0.); + fQa.Set(0., 0.); + fQc.Set(0., 0.); + fQ1.Set(0., 0.); + fQ2.Set(0., 0.); + + TH1D epEtaHist = aodep.GetHistogram(); + + fEvent = ev; + if (hists) { + for (UShort_t d=1; d<=3; d++) { + UShort_t nr = (d == 1 ? 1 : 2); + for (UShort_t q=0; qGet(d,r), &epEtaHist); + } // for q + } // for d + } + else if (h) { + CalcQVectors(h, &epEtaHist); + } + + aodep.SetEventplane(CalcEventplane(fQt)); + aodep.SetEventplaneA(CalcEventplane(fQa)); + aodep.SetEventplaneC(CalcEventplane(fQc)); + // aodep.SetEventplane1(CalcEventplane(fQ1)); + // aodep.SetEventplane2(CalcEventplane(fQ2)); + + FillHists(&aodep); + + return kTRUE; +} + +//_____________________________________________________________________ +void +AliFMDEventPlaneFinder::CalcQVectors(TH2D* h, TH1D* eHist) +{ + // + // Calculate the Q vectors + // + Double_t phi = 0, eta = 0, weight = 0; + for (Int_t e = 1; e <= h->GetNbinsX(); e++) { + Double_t qx = 0, qy = 0; + eta = h->GetXaxis()->GetBinCenter(e); + for (Int_t p = 1; p <= h->GetNbinsY(); p++) { + phi = h->GetYaxis()->GetBinCenter(p); + weight = h->GetBinContent(e, p); + if (fUsePhiWeights) weight *= GetPhiWeight(e, p); + // fix for FMD1 hole + if (e > 168 && p == 14) { + h->GetBinContent(e, 4); + if (fUsePhiWeights) weight *= GetPhiWeight(e, 4); + } + fHistPhiDist->Fill(eta, phi, weight); + // for underflowbin total Nch/eta + fHistPhiDist->Fill(eta, -1., weight); + + // increment Q vectors + qx += weight*TMath::Cos(2.*phi); + qy += weight*TMath::Sin(2.*phi); + } + TVector2 qVec = TVector2(qx, qy); + fQt += qVec; + if (eta < 0) fQa += qVec; + if (eta > 0) fQc += qVec; + // TODO: Add random ep increments + fQeta += qVec; + if (e % 10 == 0) { + eHist->Fill(eta, CalcEventplane(fQeta)); + fQeta.Set(0., 0.); + } + } + + return; +} + +//_____________________________________________________________________ +Double_t +AliFMDEventPlaneFinder::CalcEventplane(TVector2 v) const +{ + // + // Calculate the eventplane + // + Double_t ep = -1; + + if (v.Mod() == 0.) return ep; + + ep = v.Phi(); + ep /= 2.; + + if (fDebug > 0) + AliInfo(Form("Eventplane found to be: %f", ep)); + + return ep; +} + +//_____________________________________________________________________ +void +AliFMDEventPlaneFinder::FillHists(AliAODForwardEP* fmdEP) +{ + // + // Fill diagnostics histograms + // + Double_t fmdEPt = fmdEP->GetEventplane(); + Double_t fmdEPa = fmdEP->GetEventplaneA(); + Double_t fmdEPc = fmdEP->GetEventplaneC(); + // Double_t fmdEP1 = fmdEP->GetEventplane1(); + // Double_t fmdEP2 = fmdEP->GetEventplane2(); + + // FMD hists + fHistFMDEventplane->Fill(fmdEPt); + fHistFMDEventplaneA->Fill(fmdEPa); + fHistFMDEventplaneC->Fill(fmdEPc); +// fHistFMDEventplane1->Fill(fmdEP1); +// fHistFMDEventplane2->Fill(fmdEP2); + + // TPC comparison + AliEventplane* ep = fEvent->GetEventplane(); + Double_t tpcEP = -1; + if (ep) ep->GetEventplane("Q"); + + Double_t diffTPC = fmdEPt - tpcEP; + + if (diffTPC < TMath::Pi()/2.) diffTPC = TMath::Pi() - diffTPC; + if (diffTPC >= TMath::Pi()/2.) diffTPC = TMath::Pi() - diffTPC; + + fHistFMDmTPCep->Fill(diffTPC); + +// if (fmdEPt >= TMath::Pi()/2. && fmdEPt - tpcEP >= TMath::Pi()/2.) tpcEP = TMath::Pi()-tpcEP; +// if (fmdEPt < TMath::Pi()/2. && tpcEP - fmdEPt >= TMath::Pi()/2.) tpcEP = TMath::Pi()-tpcEP; + + fHistFMDvsTPCep->Fill(fmdEPt, tpcEP); + + // VZERO comparison + Double_t vzeroEP = -1; + if (ep) ep->GetEventplane("V0"); + + Double_t diffVZERO = fmdEPt - vzeroEP; + + if (diffVZERO < TMath::Pi()/2.) diffVZERO = TMath::Pi() - diffVZERO; + if (diffVZERO >= TMath::Pi()/2.) diffVZERO = TMath::Pi() - diffVZERO; + + fHistFMDmVZEROep->Fill(diffVZERO); + +// if (fmdEPt >= TMath::Pi()/2. && fmdEPt - vzeroEP >= TMath::Pi()/2.) vzeroEP = TMath::Pi()-vzeroEP; +// if (fmdEPt < TMath::Pi()/2. && vzeroEP - fmdEPt >= TMath::Pi()/2.) vzeroEP = TMath::Pi()-vzeroEP; + + fHistFMDvsVZEROep->Fill(fmdEPt, vzeroEP); + return; +} + +//_____________________________________________________________________ +Double_t +AliFMDEventPlaneFinder::GetPhiWeight(Int_t etaBin, Int_t phiBin) const +{ + // + // Get phi weight for flattening + // + Double_t phiWeight = 1; + + if (fPhiDist) { + Double_t nParticles = fPhiDist->GetBinContent(etaBin, 0); + Double_t nPhiBins = fPhiDist->GetNbinsY(); + Double_t phiDistValue = fPhiDist->GetBinContent(etaBin, phiBin); + + if (phiDistValue > 0) phiWeight = nParticles/nPhiBins/phiDistValue; + } + + return phiWeight; +} + +//____________________________________________________________________ +void +AliFMDEventPlaneFinder::SetRunNumber(Int_t run) +{ + // + // Set run number + // + if (fRunNumber == run) return; + + fRunNumber = run; + if (fUsePhiWeights) GetPhiDist(); + + return; +} + +//____________________________________________________________________ +void +AliFMDEventPlaneFinder::GetPhiDist() +{ + // + // Get phi dist from OADB + // + if (!fOADBContainer) return; + + fPhiDist = (TH2D*)fOADBContainer->GetObject(fRunNumber, "Default")->Clone(Form("fPhiDist_%d", fRunNumber)); + fList->Add(fPhiDist); + + return; +} + +//____________________________________________________________________ +void +AliFMDEventPlaneFinder::Print(Option_t* /*option*/) const +{ + // + // Print information + // + // Parameters: + // option Not used + // + char ind[gROOT->GetDirLevel()+3]; + for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' '; + ind[gROOT->GetDirLevel()] = '\0'; + std::cout << ind << ClassName() << ": " << GetName() << '\n' + << std::boolalpha + << ind << "Eventplane finder active!" << '\n' + << ind << "Loading OADB object for run number: " << fRunNumber << '\n' + << ind << "Using phi weights: " << (fUsePhiWeights ? "true" : "false") << '\n' + << std::noboolalpha + << std::flush; + std::cout << std::endl; +} +//____________________________________________________________________ +// +// EOF +// + + + diff --git a/PWGLF/FORWARD/analysis2/AliFMDEventPlaneFinder.h b/PWGLF/FORWARD/analysis2/AliFMDEventPlaneFinder.h new file mode 100644 index 00000000000..226f12f9c21 --- /dev/null +++ b/PWGLF/FORWARD/analysis2/AliFMDEventPlaneFinder.h @@ -0,0 +1,186 @@ +// This task finds the eventplane +// using the FMD +// +#ifndef ALIFMDEVENTPLANEFINDER_H +#define ALIFMDEVENTPLANEFINDER_H +/** + * @file AliFMDEventPlaneFinder.h + * @author Alexander Hansen + * @date Tue Feb 14 2012 + * + * @brief + * + * + * @ingroup pwglf_forward + */ +#include +#include +#include "AliForwardUtil.h" +class AliVEvent; +class TH1D; +class TH2F; +class TH2D; +class TString; +class AliOADBContainer; +class AliAODForwardEP; + +class AliFMDEventPlaneFinder : public TNamed +{ +public: + /** + * Constructor + */ + AliFMDEventPlaneFinder(); + /** + * Constructor + * + * @param name Name of object + */ + AliFMDEventPlaneFinder(const char* name); + /** + * Copy constructor + * + * @param o Object to copy from + */ + AliFMDEventPlaneFinder(const AliFMDEventPlaneFinder& o); + /** + * Destructor + */ + virtual ~AliFMDEventPlaneFinder(); + /** + * Assignement operator + * + * @param o Object to assign from + * + * @return Reference to this object + */ + AliFMDEventPlaneFinder& operator=(const AliFMDEventPlaneFinder& o); + /** + * Initialize this sub-algorithm + * + */ + virtual void Init(const TAxis& etaAxis); + /** + * Do the calculations + * + * @param hists Histogram cache + * + * @return true on successs + */ + Bool_t FindEventplane(AliVEvent* esd, + AliAODForwardEP& aodEp, + TH2D* h, + AliForwardUtil::Histos* hists); + /** + * Output diagnostic histograms to directory + * + * @param dir List to write in + */ + virtual void DefineOutput(TList* dir); + /** + * Print information + * + * @param option Print options + * - max Print max weights + */ + void Print(Option_t* option="") const; + /** + * Set the debug level. The higher the value the more output + * + * @param dbg Debug level + */ + void SetDebug(Int_t dbg=1) { fDebug = dbg; } + /* + * Calculate Q vectors + * + * @param h dN/detadphi histogram + * @param eHist histogram for ep vs. eta + */ + void CalcQVectors(TH2D* h, TH1D* eHist); + /* + * Calculate the eventplane from a vector + * + * @param v TVector2 of Q-vectors + * + * @return the eventplane as a double + */ + Double_t CalcEventplane(TVector2 v) const; + /* + * Set the run number, used for OADB object + * + * @param run Run number + */ + void SetRunNumber(Int_t run); + /* + * Get the run number + * + * @return returns the run number + */ + Int_t GetRunNumber() { return fRunNumber; } + /* + * Get the OADB phi distribution for flattening + */ + void GetPhiDist(); + /* + * Flag for setting the use of phi weights for flattening + * + * @param use true or false + */ + void SetUsePhiWeights(Bool_t use = kTRUE) { fUsePhiWeights = use; } + /* + * Fill diagnostics hists + * + * @param fmdEP Object containing results of FMD EP calculations + */ + void FillHists(AliAODForwardEP* fmdEP); + /* + * Set the OADB path, for using a custom OADB path and file + * + * @param fname Name of the custom OADB file, including path + */ + void SetOADBPath(Char_t* fname) { fOADBFileName = fname; } + +protected: + /* + * Get the phi weight from OADB histogram for the ep flattening + * + * @param etaBin which eta bin + * @param phiBin which phi bin + * + * @return phi weight for etaBin, phiBin as double + */ + Double_t GetPhiWeight(Int_t etaBin, Int_t phiBin) const; + + TList* fList; // List for diag. hists. + AliVEvent* fEvent; // Current event + TVector2 fQt; // Q vector for total ep + TVector2 fQa; // Q vector for sub-ep A + TVector2 fQc; // Q vector for sub-ep C + TVector2 fQ1; // Q vector for sub-ep 1 + TVector2 fQ2; // Q vector for sub-ep 2 + TVector2 fQeta; // Q vector for psi eta-dependence + TH1D* fHistFMDEventplane; // Diagnostics histogram + TH1D* fHistFMDEventplaneA;// Diagnostics histogram + TH1D* fHistFMDEventplaneC;// Diagnostics histogram + TH1D* fHistFMDEventplane1;// Diagnostics histogram + TH1D* fHistFMDEventplane2;// Diagnostics histogram + TH2D* fHistPhiDist; // Diagnostics histogram + TH1D* fHistFMDmTPCep; // Diagnostics histogram + TH2F* fHistFMDvsTPCep; // Diagnostics histogram + TH1D* fHistFMDmVZEROep; // Diagnostics histogram + TH2F* fHistFMDvsVZEROep; // Diagnostics histogram + Int_t fDebug; // Debug flag + TString fOADBFileName; // Path to OADB container + AliOADBContainer* fOADBContainer; // OADBContainer object + TH2D* fPhiDist; // Phi dist. for phi weights + Int_t fRunNumber; // Run number supplied + Bool_t fUsePhiWeights; // Flag for phi weights + + ClassDef(AliFMDEventPlaneFinder,1); // +}; + +#endif +// Local Variables: +// mode: C++ +// End: + diff --git a/PWGLF/FORWARD/analysis2/AliFMDMCTrackDensity.cxx b/PWGLF/FORWARD/analysis2/AliFMDMCTrackDensity.cxx index 63295dc558c..8972f271277 100644 --- a/PWGLF/FORWARD/analysis2/AliFMDMCTrackDensity.cxx +++ b/PWGLF/FORWARD/analysis2/AliFMDMCTrackDensity.cxx @@ -12,6 +12,9 @@ #include #include #include +#include "AliGenHijingEventHeader.h" +#include +#include //____________________________________________________________________ AliFMDMCTrackDensity::AliFMDMCTrackDensity() @@ -121,7 +124,8 @@ AliFMDMCTrackDensity::StoreParticle(AliMCParticle* particle, Double_t vz, UShort_t nC, UShort_t nT, - AliESDFMD& output) const + AliESDFMD& output, + Double_t w) const { // Store a particle. if (longest < 0) return; @@ -141,7 +145,7 @@ AliFMDMCTrackDensity::StoreParticle(AliMCParticle* particle, if (old == AliESDFMD::kInvalidMult) old = 0; // Increment count - output.SetMultiplicity(d,r,s,t,old+1); + output.SetMultiplicity(d,r,s,t,old+w); // Get track-reference stuff Double_t x = ref->X(); @@ -211,6 +215,7 @@ AliFMDMCTrackDensity::GetMother(Int_t iTr, return 0; } +// #define USE_FLOW_WEIGHTS 1 //____________________________________________________________________ Bool_t AliFMDMCTrackDensity::Calculate(const AliESDFMD& input, @@ -247,6 +252,14 @@ AliFMDMCTrackDensity::Calculate(const AliESDFMD& input, output.SetEta(ed, er, es, et, input.Eta(ed, er, es, et)); } } + +#ifdef USE_FLOW_WEIGHTS + AliGenHijingEventHeader* hd = dynamic_cast + (event.GenEventHeader()); + Double_t rp = (hd ? hd->ReactionPlaneAngle() : 0.); + Double_t b = (hd ? hd->ImpactParameter() : -1 ); +#endif + AliStack* stack = const_cast(event).Stack(); Int_t nTracks = stack->GetNtrack();//event.GetNumberOfTracks(); Int_t nPrim = stack->GetNprimary();//event.GetNumberOfPrimary(); @@ -332,7 +345,17 @@ AliFMDMCTrackDensity::Calculate(const AliESDFMD& input, fMaxConsequtiveStrips); Int_t nnT = TMath::Abs(oT - ooT) + 1; const AliMCParticle* mother = GetMother(iTr, event); - StoreParticle(particle, mother, longest, vz, nC, nnT, output); + +#ifdef USE_FLOW_WEIGHTS + phi = (mother ? mother->Phi() : particle->Phi()); + eta = (mother ? mother->Eta() : particle->Eta()); + Double_t pt = (mother ? mother->Pt() : particle->Pt()); + Int_t id = (mother ? mother->PdgCode() : 2212); + Double_t weight = CalculateWeight(eta, pt, b, phi, rp, id); +#else + Double_t weight = 1; // // +#endif + StoreParticle(particle, mother, longest, vz, nC, nnT, output, weight); longest = -1; angle = 0; nC = 1; // Reset track-ref counter - we have this->1 @@ -374,10 +397,124 @@ AliFMDMCTrackDensity::Calculate(const AliESDFMD& input, Info("Process", "I=%3d L=%3d nT=%3d (out of %3d)", iTr, longest, nT, fMaxConsequtiveStrips); const AliMCParticle* mother = GetMother(iTr, event); - StoreParticle(particle, mother, longest, vz, nC, nT, output); + +#ifdef USE_FLOW_WEIGHTS + phi = (mother ? mother->Phi() : particle->Phi()); + eta = (mother ? mother->Eta() : particle->Eta()); + Double_t pt = (mother ? mother->Pt() : particle->Pt()); + Int_t id = (mother ? mother->PdgCode() : 2212); + Double_t weight = CalculateWeight(eta, pt, b, phi, rp, id); +#else + Double_t weight = 1; // // +#endif + StoreParticle(particle, mother, longest, vz, nC, nT, output, weight); } // Loop over tracks return kTRUE; } + +//____________________________________________________________________ +Double_t +AliFMDMCTrackDensity::CalculateWeight(Double_t eta, Double_t pt, Double_t b, + Double_t phi, Double_t rp, Int_t id) const +{ + static TF1 gaus = TF1("gaus", "gaus", -6, 6); + gaus.SetParameters(0.1, 0., 9); + // gaus.SetParameters(0.1, 0., 3); + // gaus.SetParameters(0.1, 0., 15); + + const Double_t xCumulant2nd4050ALICE[] = {0.00, 0.25, 0.350, + 0.45, 0.55, 0.650, + 0.75, 0.85, 0.950, + 1.10, 1.30, 1.500, + 1.70, 1.90, 2.250, + 2.75, 3.25, 3.750, + 4.50}; + const Double_t yCumulant2nd4050ALICE[] = {0.00000, 0.043400, + 0.059911,0.073516, + 0.089756,0.105486, + 0.117391,0.128199, + 0.138013,0.158271, + 0.177726,0.196383, + 0.208277,0.216648, + 0.242954,0.249961, + 0.240131,0.269006, + 0.207796}; + const Int_t nPointsCumulant2nd4050ALICE = + sizeof(xCumulant2nd4050ALICE)/sizeof(Double_t); + static TGraph alicePointsPt2(nPointsCumulant2nd4050ALICE,xCumulant2nd4050ALICE,yCumulant2nd4050ALICE); +#if 0 + const Double_t xCumulant4th3040ALICE[] = {0.00,0.250,0.35, + 0.45,0.550,0.65, + 0.75,0.850,0.95, + 1.10,1.300,1.50, + 1.70,1.900,2.25, + 2.75,3.250,3.75, + 4.50,5.500,7.00, + 9.000000}; + const Double_t yCumulant4th3040ALICE[] = {0.000000,0.037071, + 0.048566,0.061083, + 0.070910,0.078831, + 0.091396,0.102026, + 0.109691,0.124449, + 0.139819,0.155561, + 0.165701,0.173678, + 0.191149,0.202015, + 0.204540,0.212560, + 0.195885,0.000000, + 0.000000,0.000000}; +#endif + const Double_t xCumulant4th4050ALICE[] = {0.00,0.25,0.350, + 0.45,0.55,0.650, + 0.75,0.85,0.950, + 1.10,1.30,1.500, + 1.70,1.90,2.250, + 2.75,3.25,3.750, + 4.50}; + const Double_t yCumulant4th4050ALICE[] = {0.000000,0.038646, + 0.049824,0.066662, + 0.075856,0.081583, + 0.099778,0.104674, + 0.118545,0.131874, + 0.152959,0.155348, + 0.169751,0.179052, + 0.178532,0.198851, + 0.185737,0.239901, + 0.186098}; + const Int_t nPointsCumulant4th4050ALICE = + sizeof(xCumulant4th4050ALICE)/sizeof(Double_t); + static TGraph alicePointsPt4(nPointsCumulant4th4050ALICE, + xCumulant4th4050ALICE, + yCumulant4th4050ALICE); + + const Double_t xCumulant4thTPCrefMultTPConlyAll[] = {1.75, + 4.225, + 5.965, + 7.765, + 9.215, + 10.46, + 11.565, + 12.575}; + const Double_t yCumulant4thTPCrefMultTPConlyAll[] = {0.017855,0.032440, + 0.055818,0.073137, + 0.083898,0.086690, + 0.082040,0.077777}; + const Int_t nPointsCumulant4thTPCrefMultTPConlyAll = + sizeof(xCumulant4thTPCrefMultTPConlyAll)/sizeof(Double_t); + TGraph aliceCent(nPointsCumulant4thTPCrefMultTPConlyAll, + xCumulant4thTPCrefMultTPConlyAll, + yCumulant4thTPCrefMultTPConlyAll); + + + Double_t weight = (20. * gaus.Eval(eta) * (alicePointsPt2.Eval(pt) * 0.5 + + alicePointsPt4.Eval(pt) * 0.5) + * (aliceCent.Eval(b) / aliceCent.Eval(10.46)) + * 2. * TMath::Cos(2. * (phi - rp))); + if (TMath::Abs(id) == 211) weight *= 1.3; //pion flow + else if (TMath::Abs(id) == 2212) weight *= 1.0; //proton flow + else weight *= 0.7; + + return weight; +} //____________________________________________________________________ void AliFMDMCTrackDensity::Print(Option_t* /*option*/) const diff --git a/PWGLF/FORWARD/analysis2/AliFMDMCTrackDensity.h b/PWGLF/FORWARD/analysis2/AliFMDMCTrackDensity.h index 25d28a4527a..08fa2d34e60 100644 --- a/PWGLF/FORWARD/analysis2/AliFMDMCTrackDensity.h +++ b/PWGLF/FORWARD/analysis2/AliFMDMCTrackDensity.h @@ -138,7 +138,8 @@ protected: Double_t vz, UShort_t nC, UShort_t nT, - AliESDFMD& output) const; + AliESDFMD& output, + Double_t w) const; /** * Get incident angle of this track reference * @@ -158,6 +159,21 @@ protected: * @return Pointer to mother or null */ const AliMCParticle* GetMother(Int_t iTr, const AliMCEvent& event) const; + /** + * Calculate flow weight + * + * @param eta Pseudo rapidity + * @param pt Transverse momemtum + * @param b Impact parameter + * @param phi Azimuthal angle + * @param rp Reaction plance angle + * @param id Particle PDG code + * + * @return Flow weight for the particle + */ + Double_t CalculateWeight(Double_t eta, Double_t pt, Double_t b, + Double_t phi, Double_t rp, Int_t id) const; + Bool_t fUseOnlyPrimary; // Only use primaries UShort_t fMaxConsequtiveStrips; // Max 'cluster' size TH1D* fNr; // Number of track-refs per cluster diff --git a/PWGLF/FORWARD/analysis2/AliForwardFlowTaskQC.cxx b/PWGLF/FORWARD/analysis2/AliForwardFlowTaskQC.cxx index e2813fdf41f..532b73a5f66 100644 --- a/PWGLF/FORWARD/analysis2/AliForwardFlowTaskQC.cxx +++ b/PWGLF/FORWARD/analysis2/AliForwardFlowTaskQC.cxx @@ -15,7 +15,6 @@ #include #include #include -#include "THnSparse.h" #include "TH3D.h" #include "TProfile2D.h" #include "AliLog.h" @@ -23,36 +22,27 @@ #include "AliAnalysisManager.h" #include "AliAODHandler.h" #include "AliAODInputHandler.h" -#include "AliAODMCParticle.h" #include "AliAODForwardMult.h" +#include "AliAODCentralMult.h" #include "AliAODEvent.h" -// -// Enumeration for adding and retrieving stuff from the histogram -// -enum { kW2Two = 1, kW2, kW4Four, kW4, kQnRe, kQnIm, kM, - kCosphi1phi2, kSinphi1phi2, kCosphi1phi2phi3m, kSinphi1phi2phi3m, kMm1m2, - kw2two, kw2, kw4four, kw4, kpnRe, kpnIm, kmp, - kCospsi1phi2, kSinpsi1phi2, kCospsi1phi2phi3m, kSinpsi1phi2phi3m, - kmpmq, kCospsi1phi2phi3p, kSinpsi1phi2phi3p }; - ClassImp(AliForwardFlowTaskQC) #if 0 ; // For emacs #endif AliForwardFlowTaskQC::AliForwardFlowTaskQC() - : fOutputList(0), // Output list - fFlowUtil(0), // AliForwardFlowUtil + : AliAnalysisTaskSE(), + fBinsFMD(), // List with FMD flow histos + fBinsSPD(), // List with SPD flow histos + fSumList(0), // Event sum list + fOutputList(0), // Result output list fAOD(0), // AOD input event - fMC(kFALSE), // MC flag - fEtaBins(48), // # of etaBin bins in histograms - fEtaRef(12), // # of Eta bins for reference flow - fAddFlow(0), // Add flow string - fAddType(0), // Add flow type # - fAddOrder(0), // Add flow order - fZvertex(1111), // Z vertex range - fCent(-1) // Centrality + fZvertex(1111), // Z vertex coordinate + fCent(-1), // Centrality + fHistCent(), // Histo for centrality + fHistVertexSel(), // Histo for selected vertices + fHistVertexAll() // Histo for all vertices { // // Default constructor @@ -60,19 +50,18 @@ AliForwardFlowTaskQC::AliForwardFlowTaskQC() for (Int_t n = 0; n <= 6; n++) fv[n] = kTRUE; } //_____________________________________________________________________ -AliForwardFlowTaskQC::AliForwardFlowTaskQC(const char* name) : - AliAnalysisTaskSE(name), - fOutputList(0), // Output list - fFlowUtil(0), // AliForwardFlowUtil - fAOD(0), // AOD input event - fMC(kFALSE), // MC flag - fEtaBins(48), // # of Eta bins - fEtaRef(12), // # of Eta bins for reference flow - fAddFlow(0), // Add flow string - fAddType(0), // Add flow type # - fAddOrder(0), // Add flow order - fZvertex(1111), // Z vertex range - fCent(-1) // Centrality +AliForwardFlowTaskQC::AliForwardFlowTaskQC(const char* name) + : AliAnalysisTaskSE(name), + fBinsFMD(), // List with FMD flow histos + fBinsSPD(), // List with SPD flow histos + fSumList(0), // Event sum list + fOutputList(0), // Result output list + fAOD(0), // AOD input event + fZvertex(1111), // Z vertex coordinate + fCent(-1), // Centrality + fHistCent(), // Histo for centrality + fHistVertexSel(), // Histo for selected vertices + fHistVertexAll() // Histo for all vertices { // // Constructor @@ -81,22 +70,27 @@ AliForwardFlowTaskQC::AliForwardFlowTaskQC(const char* name) : // name: Name of task // for (Int_t n = 0; n <= 6; n++) fv[n] = kTRUE; + + fHistCent = new TH1D("hCent", "Centralities", 100, 0, 100); + fHistVertexSel = new TH1D("hVertexSel", "Selected vertices", 40, -20, 20); + fHistVertexAll = new TH1D("hVertexAll", "All vertices", 40, -20, 20); + DefineOutput(1, TList::Class()); + DefineOutput(2, TList::Class()); } //_____________________________________________________________________ -AliForwardFlowTaskQC::AliForwardFlowTaskQC(const AliForwardFlowTaskQC& o) : - AliAnalysisTaskSE(o), - fOutputList(o.fOutputList), // Output list - fFlowUtil(o.fFlowUtil), // AliForwardFlowUtil - fAOD(o.fAOD), // AOD input event - fMC(o.fMC), // MC flag - fEtaBins(o.fEtaBins), // # of Eta bins - fEtaRef(o.fEtaRef), // # of Eta bins for reference flow - fAddFlow(o.fAddFlow), // Add flow string - fAddType(o.fAddType), // Add flow type # - fAddOrder(o.fAddOrder), // Add flow order - fZvertex(o.fZvertex), // Z vertex range - fCent(o.fCent) // Centrality +AliForwardFlowTaskQC::AliForwardFlowTaskQC(const AliForwardFlowTaskQC& o) + : AliAnalysisTaskSE(o), + fBinsFMD(), // List with FMD flow histos + fBinsSPD(), // List with SPD flow histos + fSumList(o.fSumList), // Event sum list + fOutputList(o.fOutputList), // Result output list + fAOD(o.fAOD), // AOD input event + fZvertex(o.fZvertex), // Z vertex coordinate + fCent(o.fCent), // Centrality + fHistCent(o.fHistCent), // Histo for centrality + fHistVertexSel(o.fHistVertexSel), // Histo for selected vertices + fHistVertexAll(o.fHistVertexAll) // Histo for all vertices { // // Copy constructor @@ -105,221 +99,557 @@ AliForwardFlowTaskQC::AliForwardFlowTaskQC(const AliForwardFlowTaskQC& o) : // o Object to copy from // for (Int_t n = 0; n <= 6; n++) fv[n] = o.fv[n]; - DefineOutput(1, TList::Class()); +} +//_____________________________________________________________________ +AliForwardFlowTaskQC& +AliForwardFlowTaskQC::operator=(const AliForwardFlowTaskQC& o) +{ + // + // Assignment operator + // + fSumList = o.fSumList; + fOutputList = o.fOutputList; + fAOD = o.fAOD; + fZvertex = o.fZvertex; + fCent = o.fCent; + fHistCent = o.fHistCent; + fHistVertexSel = o.fHistVertexSel; + fHistVertexAll = o.fHistVertexAll; + fHistVertexAll = o.fHistVertexAll; + + for (Int_t n = 0; n <= 6; n++) fv[n] = o.fv[n]; + return *this; } //_____________________________________________________________________ void AliForwardFlowTaskQC::UserCreateOutputObjects() { // // Create output objects - // - if (!fOutputList) - fOutputList = new TList(); - fOutputList->SetName("QCumulants"); - fOutputList->SetOwner(); - fFlowUtil = new AliForwardFlowUtil(fOutputList); - - Double_t x[101] = { 0. }; - // Double_t etaMin = -6; - // Double_t etaMax = 6; - - // First we have a number of options for eta binning, also if it is - // not really eta, but centrality or pt we want to do flow as a - // function of, then this is possible: - if (fEtaBins == 5) { - x[0] = 0.; - x[1] = 1.; - x[2] = 2.; - x[3] = 3.; - x[4] = 4.5; - x[5] = 6.0; - // etaMin = 0; - // etaMax = 6; - } + // + InitVertexBins(); + InitHists(); - else if (fEtaBins == 100) { - for (Int_t e = 0; e<= 100; e++) { - x[e] = e; - } - // etaMin = 0; - // etaMax = 100; - } - - else { - if (fEtaBins % 6) fEtaBins = 48; - for (Int_t e = 0; e <=fEtaBins; e++) { - x[e] = -6. + e*(12./(Double_t)fEtaBins); + PostData(1, fSumList); + PostData(2, fOutputList); + +} +//_____________________________________________________________________ +void AliForwardFlowTaskQC::InitVertexBins() +{ + // + // Init vertexbin objects for FMD and SPD, and add them to the lists + // + for(Int_t n = 1; n <= 6; n++) { + if (!fv[n]) continue; + for (Int_t v = -10; v < 10; v++) { + fBinsFMD.Add(new VertexBin(v, v+1, n, "FMD")); + fBinsSPD.Add(new VertexBin(v, v+1, n, "SPD")); } } - - // Reference flow binning - Double_t xR[101] = { 0. }; - for (Int_t r = 0; r <= fEtaRef; r++) { - xR[r] = -6 + r*(12./(Double_t)fEtaRef); - } - - // Phi binning - Double_t phi[21] = { 0. }; - for (Int_t p = 0; p <= 20; p++) { - phi[p] = p*2.*TMath::Pi() / 20.; + +} +//_____________________________________________________________________ +void AliForwardFlowTaskQC::InitHists() +{ + // + // Init histograms and add vertex bin histograms to the sum list + // + if (!fSumList) + fSumList = new TList(); + fSumList->SetName("Sums"); + fSumList->SetOwner(); + + TList* dList = new TList(); + dList->SetName("Diagnostics"); + dList->Add(fHistCent); + dList->Add(fHistVertexSel); + dList->Add(fHistVertexAll); + fSumList->Add(dList); + + TIter nextFMD(&fBinsFMD); + VertexBin* bin = 0; + while ((bin = static_cast(nextFMD()))) { + bin->AddOutput(fSumList); } - Double_t phiMC[201] = { 0. }; - for (Int_t p = 0; p <= 200; p++) { - phiMC[p] = p*2.*TMath::Pi() / 200.; + TIter nextSPD(&fBinsSPD); + while ((bin = static_cast(nextSPD()))) { + bin->AddOutput(fSumList); } - - // Histograms for cumulants analysis - - // We loop over flow histograms here to add different orders of harmonics - TString type = ""; - for (Int_t loop = 1; loop <= 5; loop++) { - - if (loop == 1) type = "FMD"; - if (loop == 2) type = "SPD"; - if (loop == 3) type = "MC"; - if (loop == 4) type = "FMDTR"; - if (loop == 5) type = "SPDTR"; - - TList* tList = new TList(); - tList->SetName(type.Data()); - fOutputList->Add(tList); - - for (Int_t n = 1; n <= 6; n++) { - if (!fv[n]) continue; - - TList* vList = new TList(); - vList->SetName(Form("v%d", n)); - tList->Add(vList); - - TString tag = TString(); - for (Int_t c = 0; c < 100; c++) { - // Output histograms - tag = Form("%d_%d", c, c+1); - - TH3D* hFlowHist = new TH3D(Form("hQ%dCumuHist%s_%s", n, type.Data(), tag.Data()), - Form("hQ%dCumuHist%s_%s", n, type.Data(), tag.Data()), fEtaBins, -6, 6, 10, -5, 5, 26, 0.5, 26.5); - TAxis* xAxis = hFlowHist->GetXaxis(); - xAxis->Set(fEtaBins, x); - hFlowHist->Sumw2(); - vList->Add(hFlowHist); - - TProfile* hCumulant2DiffFlow = new TProfile(Form("hQ%dCumulant2DiffFlow%s_%s", n, type.Data(), tag.Data()), - Form("hQ%dCumulant2DiffFlow%s_%s", n, type.Data(), tag.Data()), fEtaBins, x); - hCumulant2DiffFlow->Sumw2(); - vList->Add(hCumulant2DiffFlow); - - TProfile* hCumulant4DiffFlow = new TProfile(Form("hQ%dCumulant4DiffFlow%s_%s", n, type.Data(), tag.Data()), - Form("hQ%dCumulant4DiffFlow%s_%s", n, type.Data(), tag.Data()), fEtaBins, x); - hCumulant4DiffFlow->Sumw2(); - vList->Add(hCumulant4DiffFlow); - } // end of centrality loop - - } // end of v_{n} loop - - // Single Event histograms - TH2D* hdNdphiSE = new TH2D(Form("hdNdphiSE%s", type.Data()), - Form("hdNdphiSE%s", type.Data()), fEtaRef, xR, (type.Contains("MC") ? 200 : 20), (type.Contains("MC") ? phiMC : phi)); - hdNdphiSE->Sumw2(); - fOutputList->Add(hdNdphiSE); - - TH2D* hdNdetadphiSE = new TH2D(Form("hdNdetadphiSE%s", type.Data()), - Form("hdNdetadphiSE%s", type.Data()), fEtaBins, x, (type.Contains("MC") ? 200 : 20), (type.Contains("MC") ? phiMC : phi)); - hdNdetadphiSE->Sumw2(); - fOutputList->Add(hdNdetadphiSE); - } // end of type loop - - TProfile2D* pMCTruth = new TProfile2D("pMCTruth", "pMCTruth", 48, -6, 6, 100, 0, 100); - TAxis* xAxis = pMCTruth->GetXaxis(); - xAxis->Set(fEtaBins, x); - pMCTruth->Sumw2(); - fOutputList->Add(pMCTruth); - - // Monitoring plots - TH1D* cent = new TH1D("Centralities", "Centralities", 100, 0, 100); - fOutputList->Add(cent); - - TH1D* vertexSel = new TH1D("VertexSelected", "VertexSelected", 50, -10, 10); - fOutputList->Add(vertexSel); - TH1D* vertexAll = new TH1D("VertexAll", "VertexAll", 50, -10, 10); - fOutputList->Add(vertexAll); - - TH2D* vertex2D = new TH2D("CoverageVsVertex", "CoverageVsVertex", fEtaBins, -6, 6, 20, -10, 10); - fOutputList->Add(vertex2D); - - PostData(1, fOutputList); } //_____________________________________________________________________ void AliForwardFlowTaskQC::UserExec(Option_t */*option*/) { - // - // Process each event + // + // Calls the analyze function - called every event // // Parameters: // option: Not used // + Analyze(); + + PostData(1, fSumList); + +} +//_____________________________________________________________________ +Bool_t AliForwardFlowTaskQC::Analyze() +{ + // + // Load FMD and SPD objects from aod tree and call the cumulants + // calculation for the correct vertexbin + // + // Reset data members fCent = -1; fZvertex = 1111; // Get input event fAOD = dynamic_cast(InputEvent()); - if (!fAOD) return; + if (!fAOD) return kFALSE; - // Fill histograms - if (!fFlowUtil->LoopAODFMD(fAOD)) return; - if (!fFlowUtil->LoopAODSPD(fAOD)) return; + AliAODForwardMult* aodfmult = static_cast(fAOD->FindListObject("Forward")); + if (!aodfmult) return kFALSE; + if (!AODCheck(aodfmult)) return kFALSE; + TH2D fmddNdetadphi = aodfmult->GetHistogram(); - // Get centrality and vertex from flow utility and fill monitoring histograms - fCent = fFlowUtil->GetCentrality(); - fZvertex = fFlowUtil->GetVertex(); + AliAODCentralMult* aodcmult = static_cast(fAOD->FindListObject("CentralClusters")); + if (!aodcmult) return kFALSE; + TH2D spddNdetadphi = aodcmult->GetHistogram(); - TH1D* cent = (TH1D*)fOutputList->FindObject("Centralities"); - - cent->Fill(fCent); - TH1D* vertex = (TH1D*)fOutputList->FindObject("VertexSelected"); - vertex->Fill(fZvertex); + // TODO: remove me! +// fCent = 0.5; // Run analysis on FMD and SPD - for (Int_t n = 1; n <= 6; n++) { - if (fv[n]) { - CumulantsMethod("FMD", n); - CumulantsMethod("SPD", n); + TIter nextFMD(&fBinsFMD); + VertexBin* bin = 0; + while ((bin = static_cast(nextFMD()))) { + if (bin->CheckVertex(fZvertex)) { + if (!bin->FillHists(fmddNdetadphi)) return kFALSE; + bin->CumulantsAccumulate(fCent); } } - // And do analysis if there is. - if (fMC) { - ProcessPrimary(); + TIter nextSPD(&fBinsSPD); + while ((bin = static_cast(nextSPD()))) { + if (bin->CheckVertex(fZvertex)) { + if (!bin->FillHists(spddNdetadphi)) return kFALSE; + bin->CumulantsAccumulate(fCent); + } } - PostData(1, fOutputList); + return kTRUE; } //_____________________________________________________________________ -void AliForwardFlowTaskQC::CumulantsMethod(TString type = "", - Int_t harmonic = 2) +void AliForwardFlowTaskQC::Terminate(Option_t */*option*/) { - // - // Calculate the Q cumulant of order n + // + // Calls the finalize function, done at the end of the analysis // // Parameters: - // type: Determines which histograms should be used - // - "FMD" or "SPD" = data histograms - // - "FMDTR" or "SPDTR" = track reference histograms - // - "MC" = MC truth histograms - // harmonic: Which harmonic to calculate + // option: Not used // - Double_t n = harmonic; - - TList* tList = (TList*)fOutputList->FindObject(type.Data()); - TList* vList = (TList*)tList->FindObject(Form("v%d", harmonic)); - // We get the histograms - TH3D* flowHist = (TH3D*)vList->FindObject(Form("hQ%dCumuHist%s_%d_%d", harmonic, type.Data(), (Int_t)fCent, (Int_t)fCent+1)); - TH2D* dNdetadphi = (TH2D*)fOutputList->FindObject(Form("hdNdetadphiSE%s", type.Data())); - TH2D* dNdphi = (TH2D*)fOutputList->FindObject(Form("hdNdphiSE%s", /*"FMD"*/type.Data())); + + // Make sure pointers are set to the correct lists + fSumList = dynamic_cast (GetOutputData(1)); + if(!fSumList) { + AliError("Could not retrieve TList fSumList"); + return; + } + + if (!fOutputList) + fOutputList = new TList(); + fOutputList->SetName("Results"); + fOutputList->SetOwner(); + + // Run finalize on VertexBins + Finalize(); + + // Collect centralities + TProfile2D* hist2D = 0; + TList* centList = 0; + TH1D* hist1D = 0; + TIter nextProfile(fOutputList); + while ((hist2D = dynamic_cast(nextProfile()))) { + for (Int_t cBin = 1; cBin <= 100; ) { + Int_t cMin = cBin - 1; + Int_t cMax = (cMin < 80 ? (cMin < 20 ? cMin + 5 : cMin + 10) : cMin + 20); + TString name = Form("cent_%d-%d", cMin, cMax); + centList = (TList*)fOutputList->FindObject(name.Data()); + if (!centList) { + centList = new TList(); + centList->SetName(name.Data()); + fOutputList->Add(centList); + } + hist1D = hist2D->ProjectionX(Form("%s_%s", hist2D->GetName(), name.Data()), + cMin, cMax, "E"); + hist1D->SetTitle(hist1D->GetName()); + hist1D->Scale(1./(cMax-cMin)); + centList->Add(hist1D); + + cBin = cMax+1; + } + } + + PostData(2, fOutputList); + +} +//_____________________________________________________________________ +void AliForwardFlowTaskQC::Finalize() +{ + // + // Finalize command, called by Terminate() + // + + + // Reinitiate vertex bins if Terminate is called separately! + if (fBinsFMD.GetEntries() == 0) InitVertexBins(); + + // Iterate over all vertex bins objects and finalize cumulants + // calculations + TIter nextFMD(&fBinsFMD); + VertexBin* bin = 0; + while ((bin = static_cast(nextFMD()))) { + bin->CumulantsTerminate(fSumList, fOutputList); + } + TIter nextSPD(&fBinsSPD); + while ((bin = static_cast(nextSPD()))) { + bin->CumulantsTerminate(fSumList, fOutputList); + } + +} +// _____________________________________________________________________ +Bool_t AliForwardFlowTaskQC::AODCheck(const AliAODForwardMult* aodfm) +{ + // + // Function to check that and AOD event meets the cuts + // + // Parameters: + // AliAODForwardMult: forward mult object with trigger and vertex info + // + // Returns false if there is no trigger or if the centrality or vertex + // is out of range. Otherwise true. + // + + // First check for trigger + if (!aodfm->IsTriggerBits(AliAODForwardMult::kOffline)) return kFALSE; + + // Then check for centrality + fCent = (Double_t)aodfm->GetCentrality(); + if (0. >= fCent || fCent >= 100.) return kFALSE; + fHistCent->Fill(fCent); + + // And finally check for vertex + fZvertex = aodfm->GetIpZ(); + fHistVertexAll->Fill(fZvertex); + if (TMath::Abs(fZvertex) >= 10.) return kFALSE; + fHistVertexSel->Fill(fZvertex); + + return kTRUE; + +} +//_____________________________________________________________________ +AliForwardFlowTaskQC::VertexBin::VertexBin() + : TNamed(), + fMoment(0), // Flow moment for this vertexbin + fVzMin(0), // Vertex z-coordinate min + fVzMax(0), // Vertex z-coordinate max + fType(), // Data type (FMD/SPD/FMDTR/SPDTR/MC) + fCumuRef(), // Histogram for reference flow + fCumuDiff(), // Histogram for differential flow + fCumuHist(), // Sum histogram for cumulants + fHistTwoCorr(), // Diagnostics histogram for <<2>> + fHistW2(), // Diagnostics histogram for w_<<2>> + fHistFourCorr(), // Diagnostics histogram for <<4>> + fHistW4(), // Diagnostics histogram for w_<<4>> + fdNdedpAcc() // Diagnostics histogram to make acc. maps +{ + // + // Default constructor + // +} +//_____________________________________________________________________ +AliForwardFlowTaskQC::VertexBin::VertexBin(const Int_t vLow, const Int_t vHigh, + const Int_t moment, const Char_t* name) + : TNamed("", ""), + fMoment(moment), // Flow moment for this vertexbin + fVzMin(vLow), // Vertex z-coordinate min + fVzMax(vHigh), // Vertex z-coordinate max + fType(name), // Data type (FMD/SPD/FMDTR/SPDTR/MC) + fCumuRef(), // Histogram for reference flow + fCumuDiff(), // Histogram for differential flow + fCumuHist(), // Sum histogram for cumulants + fHistTwoCorr(), // Diagnostics histogram for <<2>> + fHistW2(), // Diagnostics histogram for w_<<2>> + fHistFourCorr(), // Diagnostics histogram for <<4>> + fHistW4(), // Diagnostics histogram for w_<<4>> + fdNdedpAcc() // Diagnostics histogram to make acc. maps +{ + // + // Constructor + // + // Parameters + // vLow: min z-coordinate + // vHigh: max z-coordinate + // moment: flow moment + // name: data type name (FMD/SPD/FMDTR/SPDTR/MC) + // + SetName(Form("%svertexBin%d_%d_%d", name, moment, vLow, vHigh)); + SetTitle(Form("%svertexBin%d_%d_%d", name, moment, vLow, vHigh)); + +} +//_____________________________________________________________________ +AliForwardFlowTaskQC::VertexBin& +AliForwardFlowTaskQC::VertexBin::operator=(const AliForwardFlowTaskQC::VertexBin& o) +{ + // + // Assignment operator + // + // Parameters + // o: AliForwardFlowTaskQC::VertexBin + // + fType = o.fType; + fCumuRef = o.fCumuRef; + fCumuDiff = o.fCumuDiff; + fCumuHist = o.fCumuHist; + fHistTwoCorr = o.fHistTwoCorr; + fHistW2 = o.fHistW2; + fHistFourCorr = o.fHistFourCorr; + fHistW4 = o.fHistW4; + fdNdedpAcc = o.fdNdedpAcc; + + return *this; +} +//_____________________________________________________________________ +void AliForwardFlowTaskQC::VertexBin::AddOutput(TList* outputlist) +{ + // + // Add histograms to outputlist + // + // Parameters + // outputlist: list of histograms + // + + // First we try to find an outputlist for this vertexbin + TList* list = (TList*)outputlist->FindObject(Form("%svertex_%d_%d", fType, fVzMin, fVzMax)); + + // If it doesn't exist we make one + if (!list) { + list = new TList(); + list->SetName(Form("%svertex_%d_%d", fType, fVzMin, fVzMax)); + outputlist->Add(list); + } + + // We initiate the reference histogram according to an acceptance correction map, + // so we don't shift the SPD coverage within a reference bin + fCumuRef = new TH2D(Form("%s_v%d_%d_%d_ref", fType, fMoment, fVzMin, fVzMax), + Form("%s_v%d_%d_%d_ref", fType, fMoment, fVzMin, fVzMax), + 24, -6., 6., 5, 0.5, 5.5); + TFile acc("$ALICE_ROOT/PWGLF/FORWARD/corrections/FlowCorrections/FlowAccMap.root", "READ"); + TH1D* accMap = (TH1D*)acc.Get(Form("%saccVertex_%d_%d", fType, fVzMin, fVzMax)); + if (accMap) { + Int_t nBins = accMap->GetNbinsX(); + Double_t eta[48] = { 0. }; + Int_t n = 0; + Double_t newOcc[48] = { 0. }; + Double_t prev = -1; + for (Int_t i = 0; i < nBins; i++) { + Double_t occ = accMap->GetBinContent(i+1); + if (prev != occ && (occ > 0.6 || occ == 0)) { + eta[n] = i*0.25-6.; + newOcc[n] = occ; + n++; +// printf("eta: %f \t occ: %f \t Vertex: %d \n", eta[n-1], occ, fVzMin); + } + prev = occ; + } + eta[n] = 6.; + fCumuRef->GetXaxis()->Set(n, eta); + } + acc.Close(); + + fCumuRef->Sumw2(); + list->Add(fCumuRef); + + // We initiate the differential histogram + fCumuDiff = new TH2D(Form("%s_v%d_%d_%d_diff", fType, fMoment, fVzMin, fVzMax), + Form("%s_v%d_%d_%d_diff", fType, fMoment, fVzMin, fVzMax), + 48, -6., 6., 5, 0.5, 5.5); + fCumuDiff->Sumw2(); + list->Add(fCumuDiff); + + // Initiate the cumulant sum histogram + fCumuHist = new TH3D(Form("%sv%d_vertex_%d_%d", fType, fMoment, fVzMin, fVzMax), + Form("%sv%d_vertex_%d_%d", fType, fMoment, fVzMin, fVzMax), + 48, -6., 6., 100, 0., 100., 26, 0.5, 26.5); + fCumuHist->Sumw2(); + + list->Add(fCumuHist); + + // We check for diagnostics histograms (only done per type and moment, not vertexbin) + // If they are not found we create them. + TList* dList = (TList*)outputlist->FindObject("Diagnostics"); + if (!dList) AliFatal("No diagnostics list found, what kind of game are you running here?!?!"); + + // Corr. hists are shared over all vertex bins... + fHistTwoCorr = (TH3D*)dList->FindObject(Form("hHistTwoCorr_%s_v%d", fType, fMoment)); + if (fHistTwoCorr) { + fHistW2 = (TH3D*)dList->FindObject(Form("hHistW2_%s_v%d", fType, fMoment)); + fHistFourCorr = (TH3D*)dList->FindObject(Form("hHistFourCorr_%s_v%d", fType, fMoment)); + fHistW4 = (TH3D*)dList->FindObject(Form("hHistW4_%s_v%d", fType, fMoment)); + } else { + fHistTwoCorr = new TH3D(Form("hHistTwoCorr_%s_v%d", fType, fMoment), + "Two particle correlator: w_{2}<<2>>", + 48, -6., 6., 100, 0., 150000, 100, 0., 100.); + fHistTwoCorr->Sumw2(); + fHistW2 = new TH3D(Form("hHistW2_%s_v%d", fType, fMoment), + "Two particle event weight: w_{2}", + 48, -6., 6., 100, 0., 2e+7, 100, 0., 100.); + fHistW2->Sumw2(); + fHistFourCorr = new TH3D(Form("hHistFourCorr_%s_v%d", fType, fMoment), + "Four particle correlator: w_{4}<<4>>", + 48, -6., 6., 100, 0., 1e+10, 100, 0., 100.); + fHistFourCorr->Sumw2(); + fHistW4 = new TH3D(Form("hHistW4_%s_v%d", fType, fMoment), + "Four particle event weight: w_{4}", + 48, -6., 6., 100, 0., 4e+14, 100, 0., 100.); + fHistW4->Sumw2(); + + dList->Add(fHistTwoCorr); + dList->Add(fHistW2); + dList->Add(fHistFourCorr); + dList->Add(fHistW4); + } - // We create the coordinate array used to fill the THnSparse, centrality and vertex is set from the beginning. - Double_t coord[4] = { 0., fCent, fZvertex, 0. }; + // Acceptance hists are shared over all moments + fdNdedpAcc = (TH2D*)dList->FindObject(Form("h%sdNdedpAcc_%d_%d", fType, fVzMin, fVzMax)); + if (!fdNdedpAcc) { + fdNdedpAcc = new TH2D(Form("h%sdNdedpAcc_%d_%d", fType, fVzMin, fVzMax), + Form("%s acceptance map for %d cm < v_{z} < %d cm", fType, fVzMin, fVzMax), + 48, -6, 6, 20, 0, TMath::TwoPi()); + fdNdedpAcc->Sumw2(); + dList->Add(fdNdedpAcc); + } + + +} +//_____________________________________________________________________ +Bool_t AliForwardFlowTaskQC::VertexBin::CheckVertex(Double_t vz) +{ + // + // We check if this is the correct bin for the current event's vertex + // + // Parameters: + // vZ: Current event vertex + // + // Returns false if out of range, true otherwise + // + if ((Double_t)fVzMin > vz) return kFALSE; + if ((Double_t)fVzMax <= vz) return kFALSE; + + return kTRUE; +} +//_____________________________________________________________________ +Bool_t AliForwardFlowTaskQC::VertexBin::FillHists(TH2D dNdetadphi) +{ + // + // Fill reference and differential eta-histograms + // + // Parameters: + // dNdetadphi: 2D histogram with input data + // + + if (!fCumuRef) AliFatal("You have not called AddOutput() - Terminating!"); + + // Fist we reset histograms + fCumuRef->Reset(); + fCumuDiff->Reset(); + + // Numbers to cut away bad events and acceptance. + Double_t runAvg = 0; + Double_t max = 0; + Int_t nInAvg = 0; + Int_t nBadBins = 0; + Int_t nBins = (dNdetadphi.GetNbinsX() * 6) / (fCumuDiff->GetNbinsX() * 5); + Int_t nInBin = 0; + Int_t nCurBin = 0, nPrevBin = 0; + + // Then we loop over the input and calculate sum cos(k*n*phi) + // and fill it in the reference and differential histograms + Double_t eta, phi, weight; + Double_t dQnRe = 0, dQ2nRe = 0, dQnIm = 0, dQ2nIm = 0; + for (Int_t etaBin = 1; etaBin <= dNdetadphi.GetNbinsX(); etaBin++) { + eta = dNdetadphi.GetXaxis()->GetBinCenter(etaBin); + nCurBin = fCumuDiff->GetXaxis()->FindBin(eta); + // If we have moved to a new bin in the flow hist, and less than half the eta + // region has been covered by it we cut it away. + if (nPrevBin && nCurBin != nPrevBin) { + if (nInBin <= nBins/2) { + for (Int_t pBin = 1; pBin <= fCumuDiff->GetNbinsY(); pBin++) { + fCumuDiff->SetBinContent(nPrevBin, pBin, 0); + fCumuDiff->SetBinError(nPrevBin, pBin, 0); + } + } + nInBin = 0; + nPrevBin = nCurBin; + runAvg = 0; + nInAvg = 0; + max = 0; + } + Bool_t data = kFALSE; + for (Int_t phiBin = 1; phiBin <= dNdetadphi.GetNbinsY(); phiBin++) { + phi = dNdetadphi.GetYaxis()->GetBinCenter(phiBin); + weight = dNdetadphi.GetBinContent(etaBin, phiBin); + if (!weight) continue; + if (!data) data = kTRUE; + // We calculate the running average Nch per. bin + runAvg *= nInAvg; + runAvg += weight; + nInAvg++; + runAvg /= nInAvg; + // if the current bin has more than write the avg we count a bad bin + if (weight > max) max = weight; + + dQnRe = weight*TMath::Cos(fMoment*phi); + dQnIm = weight*TMath::Sin(fMoment*phi); + dQ2nRe = weight*TMath::Cos(2.*fMoment*phi); + dQ2nIm = weight*TMath::Sin(2.*fMoment*phi); + + fCumuRef->Fill(eta, kHmult, weight); + fCumuRef->Fill(eta, kHQnRe, dQnRe); + fCumuRef->Fill(eta, kHQnIm, dQnIm); + fCumuRef->Fill(eta, kHQ2nRe, dQ2nRe); + fCumuRef->Fill(eta, kHQ2nIm, dQ2nIm); + + fCumuDiff->Fill(eta, kHmult, weight); + fCumuDiff->Fill(eta, kHQnRe, dQnRe); + fCumuDiff->Fill(eta, kHQnIm, dQnIm); + fCumuDiff->Fill(eta, kHQ2nRe, dQ2nRe); + fCumuDiff->Fill(eta, kHQ2nIm, dQ2nIm); + + // Fill acc. map + fdNdedpAcc->Fill(eta, phi, weight); + } + if (data) { + nInBin++; + if (max > 35) nBadBins++; + } + // If there are too many bad bins we throw the event away! + if (nBadBins > 2) return kFALSE; + } + return kTRUE; +} +//_____________________________________________________________________ +void AliForwardFlowTaskQC::VertexBin::CumulantsAccumulate(Double_t cent) +{ + // + // Calculate the Q cumulant of order fMoment + // + // Parameters: + // cent: Centrality of event + // + + if (!fCumuRef) AliFatal("You have not called AddOutput() - Terminating!"); + // We create the objects needed for the analysis Double_t dQnRe = 0, dQ2nRe = 0, dQnIm = 0, dQ2nIm = 0, mult = 0; Double_t pnRe = 0, p2nRe = 0, qnRe = 0, q2nRe = 0, pnIm = 0, p2nIm = 0, qnIm = 0, q2nIm = 0; @@ -328,16 +658,16 @@ void AliForwardFlowTaskQC::CumulantsMethod(TString type = "", Double_t sinPhi1Phi2 = 0, sinPhi1Phi2Phi3m = 0; Double_t cosPsi1Phi2 = 0, cosPsi1Phi2Phi3m = 0, cosPsi1Phi2Phi3p = 0; Double_t sinPsi1Phi2 = 0, sinPsi1Phi2Phi3m = 0, sinPsi1Phi2Phi3p = 0; - Double_t phi = 0, eta = 0; + Double_t eta = 0; Double_t multi = 0, multp = 0, mp = 0, mq = 0; Double_t w2 = 0, w4 = 0, w2p = 0, w4p = 0; Int_t refEtaBin = 0; - Bool_t kEventCount = kFALSE; - + Bool_t eventCount = kFALSE; + // We loop over the data 1 time! - for (Int_t etaBin = 1; etaBin <= dNdetadphi->GetNbinsX(); etaBin++) { - eta = dNdetadphi->GetXaxis()->GetBinCenter(etaBin); - refEtaBin = dNdphi->GetXaxis()->FindBin(eta); + for (Int_t etaBin = 1; etaBin <= fCumuDiff->GetNbinsX(); etaBin++) { + eta = fCumuDiff->GetXaxis()->GetBinCenter(etaBin); + refEtaBin = fCumuRef->GetXaxis()->FindBin(eta); // The values for each individual etaBin bins are reset mp = 0; pnRe = 0; @@ -351,53 +681,43 @@ void AliForwardFlowTaskQC::CumulantsMethod(TString type = "", dQ2nRe = 0; dQ2nIm = 0; - for (Int_t phiBin = 1; phiBin <= dNdphi->GetNbinsY(); phiBin++) { - phi = dNdphi->GetYaxis()->GetBinCenter(phiBin); - - // Reference flow - multi = dNdphi->GetBinContent(refEtaBin, phiBin); - dQnRe += multi*TMath::Cos(n*phi); - dQnIm += multi*TMath::Sin(n*phi); - dQ2nRe += multi*TMath::Cos(2.*n*phi); - dQ2nIm += multi*TMath::Sin(2.*n*phi); - mult += multi; - - // For each etaBin bin the necessary values for differential flow - // is calculated. Here is the loop over the phi's. - multp = dNdetadphi->GetBinContent(etaBin, phiBin); - mp += multp; - pnRe += multp*TMath::Cos(n*phi); - pnIm += multp*TMath::Sin(n*phi); - p2nRe += multp*TMath::Cos(2.*n*phi); - p2nIm += multp*TMath::Sin(2.*n*phi); - } - if (mult == 0) continue; + // Reference flow + multi = fCumuRef->GetBinContent(refEtaBin, kHmult); + dQnRe = fCumuRef->GetBinContent(refEtaBin, kHQnRe); + dQnIm = fCumuRef->GetBinContent(refEtaBin, kHQnIm); + dQ2nRe = fCumuRef->GetBinContent(refEtaBin, kHQ2nRe); + dQ2nIm = fCumuRef->GetBinContent(refEtaBin, kHQ2nIm); + mult += multi; + + // For each etaBin bin the necessary values for differential flow + // is calculated. Here is the loop over the phi's. + multp = fCumuDiff->GetBinContent(etaBin, kHmult); + pnRe = fCumuDiff->GetBinContent(etaBin, kHQnRe); + pnIm = fCumuDiff->GetBinContent(etaBin, kHQnIm); + p2nRe = fCumuDiff->GetBinContent(etaBin, kHQ2nRe); + p2nIm = fCumuDiff->GetBinContent(etaBin, kHQ2nIm); + mp += multp; + + if (mult <= 3) continue; - if (!kEventCount) { + if (!eventCount) { // Count number of events - coord[0] = -7; - coord[3] = -0.5; - flowHist->Fill(coord[0], coord[2], coord[3], 1.); - kEventCount = kTRUE; + fCumuHist->Fill(-7., cent, -0.5, 1.); + eventCount = kTRUE; } - + if (mp == 0) continue; // The reference flow is calculated - coord[0] = eta; // 2-particle w2 = mult * (mult - 1.); two = dQnRe*dQnRe + dQnIm*dQnIm - mult; - coord[3] = kW2Two; - flowHist->Fill(coord[0], coord[2], coord[3], two); - coord[3] = kW2; - flowHist->Fill(coord[0], coord[2], coord[3], w2); - - coord[3] = kQnRe; - flowHist->Fill(coord[0], coord[2], coord[3], dQnRe); - coord[3] = kQnIm; - flowHist->Fill(coord[0], coord[2], coord[3], dQnIm); - coord[3] = kM; - flowHist->Fill(coord[0], coord[2], coord[3], mult); + + fCumuHist->Fill(eta, cent, kW2Two, two); + fCumuHist->Fill(eta, cent, kW2, w2); + + fCumuHist->Fill(eta, cent, kQnRe, dQnRe); + fCumuHist->Fill(eta, cent, kQnIm, dQnIm); + fCumuHist->Fill(eta, cent, kM, mult); // 4-particle w4 = mult * (mult - 1.) * (mult - 2.) * (mult - 3.); @@ -407,10 +727,8 @@ void AliForwardFlowTaskQC::CumulantsMethod(TString type = "", -2.*(TMath::Power(dQnRe,2.)*dQ2nRe+2.*dQnRe*dQnIm*dQ2nIm-TMath::Power(dQnIm,2.)*dQ2nRe) +(TMath::Power(dQ2nRe,2.)+TMath::Power(dQ2nIm,2.)); - coord[3] = kW4Four; - flowHist->Fill(coord[0], coord[2], coord[3], four); - coord[3] = kW4; - flowHist->Fill(coord[0], coord[2], coord[3], w4); + fCumuHist->Fill(eta, cent, kW4Four, four); + fCumuHist->Fill(eta, cent, kW4, w4); cosPhi1Phi2 = dQnRe*dQnRe - dQnIm*dQnIm - dQ2nRe; sinPhi1Phi2 = 2.*dQnRe*dQnIm - dQ2nIm; @@ -419,18 +737,19 @@ void AliForwardFlowTaskQC::CumulantsMethod(TString type = "", sinPhi1Phi2Phi3m = -dQnIm*(TMath::Power(dQnRe,2)+TMath::Power(dQnIm,2))+dQnRe*dQ2nIm-dQnIm*dQ2nRe+2.*(mult-1)*dQnIm; - coord[3] = kCosphi1phi2; - flowHist->Fill(coord[0], coord[2], coord[3], cosPhi1Phi2); - coord[3] = kSinphi1phi2; - flowHist->Fill(coord[0], coord[2], coord[3], sinPhi1Phi2); - coord[3] = kCosphi1phi2phi3m; - flowHist->Fill(coord[0], coord[2], coord[3], cosPhi1Phi2Phi3m); - coord[3] = kSinphi1phi2phi3m; - flowHist->Fill(coord[0], coord[2], coord[3], sinPhi1Phi2Phi3m); - coord[3] = kMm1m2; - flowHist->Fill(coord[0], coord[2], coord[3], mult*(mult-1.)*(mult-2.)); - - // Differential flow calculations for each etaBin bin is done: + fCumuHist->Fill(eta, cent, kCosphi1phi2, cosPhi1Phi2); + fCumuHist->Fill(eta, cent, kSinphi1phi2, sinPhi1Phi2); + fCumuHist->Fill(eta, cent, kCosphi1phi2phi3m, cosPhi1Phi2Phi3m); + fCumuHist->Fill(eta, cent, kSinphi1phi2phi3m, sinPhi1Phi2Phi3m); + fCumuHist->Fill(eta, cent, kMm1m2, mult*(mult-1.)*(mult-2.)); + + // Diagnostics are filled + fHistTwoCorr->Fill(eta, two, cent); + fHistW2->Fill(eta, w2, cent); + fHistFourCorr->Fill(eta, four, cent); + fHistW4->Fill(eta, w4, cent); + + // Differential flow calculations for each eta bin bin is done: mq = mp; qnRe = pnRe; qnIm = pnIm; @@ -441,17 +760,12 @@ void AliForwardFlowTaskQC::CumulantsMethod(TString type = "", w2p = mp * mult - mq; twoPrime = pnRe*dQnRe + pnIm*dQnIm - mq; - coord[3] = kw2two; - flowHist->Fill(coord[0], coord[2], coord[3], twoPrime); - coord[3] = kw2; - flowHist->Fill(coord[0], coord[2], coord[3], w2p); - - coord[3] = kpnRe; - flowHist->Fill(coord[0], coord[2], coord[3], pnRe); - coord[3] = kpnIm; - flowHist->Fill(coord[0], coord[2], coord[3], pnIm); - coord[3] = kmp; - flowHist->Fill(coord[0], coord[2], coord[3], mp); + fCumuHist->Fill(eta, cent, kw2two, twoPrime); + fCumuHist->Fill(eta, cent, kw2, w2p); + + fCumuHist->Fill(eta, cent, kpnRe, pnRe); + fCumuHist->Fill(eta, cent, kpnIm, pnIm); + fCumuHist->Fill(eta, cent, kmp, mp); // 4-particle differential flow w4p = (mp * mult - 3.*mq)*(mult - 1.)*(mult - 2.); @@ -468,11 +782,9 @@ void AliForwardFlowTaskQC::CumulantsMethod(TString type = "", + 2.*(pnRe*dQnRe+pnIm*dQnIm) + 2.*mq*mult - 6.*mq; - - coord[3] = kw4four; - flowHist->Fill(coord[0], coord[2], coord[3], fourPrime); - coord[3] = kw4; - flowHist->Fill(coord[0], coord[2], coord[3], w4p); + + fCumuHist->Fill(eta, cent, kw4four, fourPrime); + fCumuHist->Fill(eta, cent, kw4, w4p); cosPsi1Phi2 = pnRe*dQnRe - pnIm*dQnIm - q2nRe; sinPsi1Phi2 = pnRe*dQnIm + pnIm*dQnRe - q2nIm; @@ -493,43 +805,51 @@ void AliForwardFlowTaskQC::CumulantsMethod(TString type = "", - 1.*(pnIm*dQ2nRe-pnRe*dQ2nIm) + 2.*mq*dQnIm-2.*qnIm; - coord[3] = kCospsi1phi2; - flowHist->Fill(coord[0], coord[2], coord[3], cosPsi1Phi2); - coord[3] = kSinpsi1phi2; - flowHist->Fill(coord[0], coord[2], coord[3], sinPsi1Phi2); - coord[3] = kCospsi1phi2phi3m; - flowHist->Fill(coord[0], coord[2], coord[3], cosPsi1Phi2Phi3m); - coord[3] = kSinpsi1phi2phi3m; - flowHist->Fill(coord[0], coord[2], coord[3], sinPsi1Phi2Phi3m); - coord[3] = kmpmq; - flowHist->Fill(coord[0], coord[2], coord[3], (mp*mult-2.*mq)*(mult-1.)); - coord[3] = kCospsi1phi2phi3p; - flowHist->Fill(coord[0], coord[2], coord[3], cosPsi1Phi2Phi3p); - coord[3] = kSinpsi1phi2phi3p; - flowHist->Fill(coord[0], coord[2], coord[3], sinPsi1Phi2Phi3p); - + fCumuHist->Fill(eta, cent, kCospsi1phi2, cosPsi1Phi2); + fCumuHist->Fill(eta, cent, kSinpsi1phi2, sinPsi1Phi2); + fCumuHist->Fill(eta, cent, kCospsi1phi2phi3m, cosPsi1Phi2Phi3m); + fCumuHist->Fill(eta, cent, kSinpsi1phi2phi3m, sinPsi1Phi2Phi3m); + fCumuHist->Fill(eta, cent, kmpmq, (mp*mult-2.*mq)*(mult-1.)); + fCumuHist->Fill(eta, cent, kCospsi1phi2phi3p, cosPsi1Phi2Phi3p); + fCumuHist->Fill(eta, cent, kSinpsi1phi2phi3p, sinPsi1Phi2Phi3p); } } //_____________________________________________________________________ -void AliForwardFlowTaskQC::Terminate(Option_t */*option*/) +void AliForwardFlowTaskQC::VertexBin::CumulantsTerminate(TList* inlist, TList* outlist) { // - // End of job // Finalizes the Q cumulant calculations // // Parameters: - // option Not used + // inlist: input sumlist + // outlist: output result list // - TH3D* cumulantsHist = 0; - TProfile* cumulant2diffHist = 0; - TProfile* cumulant4diffHist = 0; - TList* tList = 0; - TList* vList = 0; + // Re-find cumulants hist if Terminate is called separately + if (!fCumuHist) { + TList* list = (TList*)inlist->FindObject(Form("%svertex_%d_%d", fType, fVzMin, fVzMax)); + fCumuHist = (TH3D*)list->FindObject(Form("%sv%d_vertex_%d_%d", fType, fMoment, fVzMin, fVzMax)); + } + + // Create result profiles + TProfile2D* cumu2 = (TProfile2D*)outlist->FindObject(Form("%sQC2_v%d_unCorr", fType, fMoment)); + TProfile2D* cumu4 = (TProfile2D*)outlist->FindObject(Form("%sQC4_v%d_unCorr", fType, fMoment)); + if (!cumu2) { + cumu2 = new TProfile2D(Form("%sQC2_v%d_unCorr", fType, fMoment), + Form("%sQC2_v%d_unCorr", fType, fMoment), + 48, -6., 6., 100, 0., 100); + outlist->Add(cumu2); + } + if (!cumu4) { + cumu4 = new TProfile2D(Form("%sQC4_v%d_unCorr", fType, fMoment), + Form("%sQC4_v%d_unCorr", fType, fMoment), + 48, -6., 6., 100, 0., 100); + outlist->Add(cumu4); + } // For flow calculations - Double_t two = 0, qc2 = 0, /* vnTwo = 0, */ four = 0, qc4 = 0 /*, vnFour = 0*/; + Double_t two = 0, qc2 = 0, vnTwo = 0, four = 0, qc4 = 0, vnFour = 0; Double_t twoPrime = 0, qc2Prime = 0, vnTwoDiff = 0, fourPrime = 0, qc4Prime = 0, vnFourDiff = 0; Double_t w2 = 0, w4 = 0, w2p = 0, w4p = 0; Double_t w2Two = 0, w2pTwoPrime = 0, w4Four = 0, w4pFourPrime = 0; @@ -539,191 +859,126 @@ void AliForwardFlowTaskQC::Terminate(Option_t */*option*/) Double_t cosP1nPsi1M1nPhi2M1nPhi3 = 0, sinP1nPsi1M1nPhi2M1nPhi3 = 0, mpqMult = 0; Double_t cosP1nPsi1P1nPhi2M1nPhi3 = 0, sinP1nPsi1P1nPhi2M1nPhi3 = 0; - Int_t nLoops = (fMC ? 5 : 2); - TString type = ""; - - // Do a loop over the difference analysis types, calculating flow - // 2 loops for real data, 3 for MC data - // inside each is a nested loop over each harmonic (1, 2, 3 and 4 at the moment) - for (Int_t loop = 1; loop <= nLoops; loop++) { - - if (loop == 1) type = "FMD"; - if (loop == 2) type = "SPD"; - if (loop == 3) type = "MC"; - if (loop == 4) type = "FMDTR"; - if (loop == 5) type = "SPDTR"; - - for (Int_t n = 1; n <= 6; n++) { - if (!fv[n]) continue; - - tList = (TList*)fOutputList->FindObject(type.Data()); - vList = (TList*)tList->FindObject(Form("v%d", n)); - TString tag = ""; - - // Centrality loop - for (Int_t c = 0; c < 100; c++) { - - Double_t nEv = 0; - tag = Form("%d_%d", c, c+1); - cumulantsHist = (TH3D*)vList->FindObject(Form("hQ%dCumuHist%s_%s", n, type.Data(), tag.Data())); - - cumulant2diffHist = (TProfile*)vList->FindObject(Form("hQ%dCumulant2DiffFlow%s_%s", n, type.Data(), tag.Data())); - cumulant4diffHist = (TProfile*)vList->FindObject(Form("hQ%dCumulant4DiffFlow%s_%s", n, type.Data(), tag.Data())); + // Loop over cumulant histogram for final calculations + // Centrality loop + for (Int_t c = 1; c <= 80; c++) { + Double_t nEv = 0; + // Eta loop + for (Int_t etaBin = 1; etaBin <= fCumuHist->GetNbinsX(); etaBin++) { + Double_t eta = fCumuHist->GetXaxis()->GetBinCenter(etaBin); + // 2-particle reference flow + w2Two = fCumuHist->GetBinContent(etaBin, c, kW2Two); + w2 = fCumuHist->GetBinContent(etaBin, c, kW2); + mult = fCumuHist->GetBinContent(etaBin, c, kM); + if (!w2 || !mult) continue; + cosP1nPhi = fCumuHist->GetBinContent(etaBin, c, kQnRe); + sinP1nPhi = fCumuHist->GetBinContent(etaBin, c, kQnIm); - for (Int_t vertexBin = 1; vertexBin <= cumulantsHist->GetNbinsY(); vertexBin++) { - - for (Int_t etaBin = 1; etaBin <= cumulantsHist->GetNbinsX(); etaBin++) { - Double_t eta = cumulantsHist->GetXaxis()->GetBinCenter(etaBin); - // 2-particle reference flow - w2Two = cumulantsHist->GetBinContent(etaBin, vertexBin, kW2Two); - if (!w2Two) continue; - w2 = cumulantsHist->GetBinContent(etaBin, vertexBin, kW2); - cosP1nPhi = cumulantsHist->GetBinContent(etaBin, vertexBin, kQnRe); - sinP1nPhi = cumulantsHist->GetBinContent(etaBin, vertexBin, kQnIm); - mult = cumulantsHist->GetBinContent(etaBin, vertexBin, kM); - - cosP1nPhi /= mult; - sinP1nPhi /= mult; - two = w2Two / w2; - qc2 = two - TMath::Power(cosP1nPhi, 2) - TMath::Power(sinP1nPhi, 2); - if (qc2 <= 0) continue; - // vnTwo = TMath::Sqrt(qc2); - // if (!TMath::IsNaN(vnTwo*mult)) - // cumulant2diffHist->Fill(eta, vnTwo, cumulantsHist->GetBinContent(0,vertexBin,0)); - - // 4-particle reference flow - w4Four = cumulantsHist->GetBinContent(etaBin, vertexBin, kW4Four); - w4 = cumulantsHist->GetBinContent(etaBin, vertexBin, kW4); - cosP1nPhi1P1nPhi2 = cumulantsHist->GetBinContent(etaBin, vertexBin, kCosphi1phi2); - sinP1nPhi1P1nPhi2 = cumulantsHist->GetBinContent(etaBin, vertexBin, kSinphi1phi2); - cosP1nPhi1M1nPhi2M1nPhi3 = cumulantsHist->GetBinContent(etaBin, vertexBin, kCosphi1phi2phi3m); - sinP1nPhi1M1nPhi2M1nPhi3 = cumulantsHist->GetBinContent(etaBin, vertexBin, kSinphi1phi2phi3m); - multm1m2 = cumulantsHist->GetBinContent(etaBin, vertexBin, kMm1m2); - - cosP1nPhi1P1nPhi2 /= w2; - sinP1nPhi1P1nPhi2 /= w2; - cosP1nPhi1M1nPhi2M1nPhi3 /= multm1m2; - sinP1nPhi1M1nPhi2M1nPhi3 /= multm1m2; - four = w4Four / w4; - qc4 = four-2.*TMath::Power(two,2.) - - 4.*cosP1nPhi*cosP1nPhi1M1nPhi2M1nPhi3 - + 4.*sinP1nPhi*sinP1nPhi1M1nPhi2M1nPhi3-TMath::Power(cosP1nPhi1P1nPhi2,2.)-TMath::Power(sinP1nPhi1P1nPhi2,2.) - + 4.*cosP1nPhi1P1nPhi2*(TMath::Power(cosP1nPhi,2.)-TMath::Power(sinP1nPhi,2.)) - + 8.*sinP1nPhi1P1nPhi2*sinP1nPhi*cosP1nPhi - + 8.*two*(TMath::Power(cosP1nPhi,2.)+TMath::Power(sinP1nPhi,2.)) - - 6.*TMath::Power((TMath::Power(cosP1nPhi,2.)+TMath::Power(sinP1nPhi,2.)),2.); - - if (qc4 >= 0) continue; - // vnFour = TMath::Power(-qc4, 0.25); - // if (!TMath::IsNaN(vnFour*mult)) - // cumulant4diffHist->Fill(eta, vnFour, cumulantsHist->GetBinContent(0,vertexBin,0)); - - // 2-particle differential flow - w2pTwoPrime = cumulantsHist->GetBinContent(etaBin, vertexBin, kw2two); - if (!w2pTwoPrime) continue; - w2p = cumulantsHist->GetBinContent(etaBin, vertexBin, kw2); - cosP1nPsi = cumulantsHist->GetBinContent(etaBin, vertexBin, kpnRe); - sinP1nPsi = cumulantsHist->GetBinContent(etaBin, vertexBin, kpnIm); - mp = cumulantsHist->GetBinContent(etaBin, vertexBin, kmp); - - cosP1nPsi /= mp; - sinP1nPsi /= mp; - twoPrime = w2pTwoPrime / w2p; - qc2Prime = twoPrime - sinP1nPsi*sinP1nPhi - cosP1nPsi*cosP1nPhi; - - vnTwoDiff = qc2Prime / TMath::Sqrt(qc2); - if (!TMath::IsNaN(vnTwoDiff*mp) && vnTwoDiff > 0) cumulant2diffHist->Fill(eta, vnTwoDiff, cumulantsHist->GetBinContent(0,vertexBin,0)); - - // 4-particle differential flow - w4pFourPrime = cumulantsHist->GetBinContent(etaBin, vertexBin, kw4four); - w4p = cumulantsHist->GetBinContent(etaBin, vertexBin, kw4); - cosP1nPsi1P1nPhi2 = cumulantsHist->GetBinContent(etaBin, vertexBin, kCospsi1phi2); - sinP1nPsi1P1nPhi2 = cumulantsHist->GetBinContent(etaBin, vertexBin, kSinpsi1phi2); - cosP1nPsi1M1nPhi2M1nPhi3 = cumulantsHist->GetBinContent(etaBin, vertexBin, kCospsi1phi2phi3m); - sinP1nPsi1M1nPhi2M1nPhi3 = cumulantsHist->GetBinContent(etaBin, vertexBin, kSinpsi1phi2phi3m); - cosP1nPsi1P1nPhi2M1nPhi3 = cumulantsHist->GetBinContent(etaBin, vertexBin, kCospsi1phi2phi3p); - sinP1nPsi1P1nPhi2M1nPhi3 = cumulantsHist->GetBinContent(etaBin, vertexBin, kSinpsi1phi2phi3p); - mpqMult = cumulantsHist->GetBinContent(etaBin, vertexBin, kmpmq); - - cosP1nPsi1P1nPhi2 /= w2p; - sinP1nPsi1P1nPhi2 /= w2p; - cosP1nPsi1M1nPhi2M1nPhi3 /= mpqMult; - sinP1nPsi1M1nPhi2M1nPhi3 /= mpqMult; - cosP1nPsi1P1nPhi2M1nPhi3 /= mpqMult; - sinP1nPsi1P1nPhi2M1nPhi3 /= mpqMult; - fourPrime = w4pFourPrime / w4p; - - qc4Prime = fourPrime-2.*twoPrime*two - - cosP1nPsi*cosP1nPhi1M1nPhi2M1nPhi3 - + sinP1nPsi*sinP1nPhi1M1nPhi2M1nPhi3 - - cosP1nPhi*cosP1nPsi1M1nPhi2M1nPhi3 - + sinP1nPhi*sinP1nPsi1M1nPhi2M1nPhi3 - - 2.*cosP1nPhi*cosP1nPsi1P1nPhi2M1nPhi3 - - 2.*sinP1nPhi*sinP1nPsi1P1nPhi2M1nPhi3 - - cosP1nPsi1P1nPhi2*cosP1nPhi1P1nPhi2 - - sinP1nPsi1P1nPhi2*sinP1nPhi1P1nPhi2 - + 2.*cosP1nPhi1P1nPhi2*(cosP1nPsi*cosP1nPhi-sinP1nPsi*sinP1nPhi) - + 2.*sinP1nPhi1P1nPhi2*(cosP1nPsi*sinP1nPhi+sinP1nPsi*cosP1nPhi) - + 4.*two*(cosP1nPsi*cosP1nPhi+sinP1nPsi*sinP1nPhi) - + 2.*cosP1nPsi1P1nPhi2*(TMath::Power(cosP1nPhi,2.)-TMath::Power(sinP1nPhi,2.)) - + 4.*sinP1nPsi1P1nPhi2*cosP1nPhi*sinP1nPhi - + 4.*twoPrime*(TMath::Power(cosP1nPhi,2.)+TMath::Power(sinP1nPhi,2.)) - - 6.*(TMath::Power(cosP1nPhi,2.)-TMath::Power(sinP1nPhi,2.)) - * (cosP1nPsi*cosP1nPhi-sinP1nPsi*sinP1nPhi) - - 12.*cosP1nPhi*sinP1nPhi - * (sinP1nPsi*cosP1nPhi+cosP1nPsi*sinP1nPhi); - - vnFourDiff = - qc4Prime / TMath::Power(-qc4, 0.75); - if (!TMath::IsNaN(vnFourDiff*mp) && vnFourDiff > 0) cumulant4diffHist->Fill(eta, vnFourDiff, cumulantsHist->GetBinContent(0,vertexBin,0)); - } // End of etaBin loop - nEv += cumulantsHist->GetBinContent(0,vertexBin,0); - } // End of vertexBin loop - // Number of events: - cumulant2diffHist->Fill(7., nEv); - cumulant4diffHist->Fill(7., nEv); - - } // End of centrality loop - } // End of harmonics loop - } // End of type loop - - PostData(1, fOutputList); + cosP1nPhi /= mult; + sinP1nPhi /= mult; + two = w2Two / w2; + qc2 = two - TMath::Power(cosP1nPhi, 2) - TMath::Power(sinP1nPhi, 2); + if (qc2 <= 0) continue; + vnTwo = TMath::Sqrt(qc2); + // if (!TMath::IsNaN(vnTwo*mult)) + // cumu2->Fill(eta, vnTwo, fCumuHist->GetBinContent(0,c,0)); + + // 4-particle reference flow + w4Four = fCumuHist->GetBinContent(etaBin, c, kW4Four); + w4 = fCumuHist->GetBinContent(etaBin, c, kW4); + multm1m2 = fCumuHist->GetBinContent(etaBin, c, kMm1m2); + if (!w4 || !multm1m2) continue; + cosP1nPhi1P1nPhi2 = fCumuHist->GetBinContent(etaBin, c, kCosphi1phi2); + sinP1nPhi1P1nPhi2 = fCumuHist->GetBinContent(etaBin, c, kSinphi1phi2); + cosP1nPhi1M1nPhi2M1nPhi3 = fCumuHist->GetBinContent(etaBin, c, kCosphi1phi2phi3m); + sinP1nPhi1M1nPhi2M1nPhi3 = fCumuHist->GetBinContent(etaBin, c, kSinphi1phi2phi3m); + + cosP1nPhi1P1nPhi2 /= w2; + sinP1nPhi1P1nPhi2 /= w2; + cosP1nPhi1M1nPhi2M1nPhi3 /= multm1m2; + sinP1nPhi1M1nPhi2M1nPhi3 /= multm1m2; + four = w4Four / w4; + qc4 = four-2.*TMath::Power(two,2.) + - 4.*cosP1nPhi*cosP1nPhi1M1nPhi2M1nPhi3 + + 4.*sinP1nPhi*sinP1nPhi1M1nPhi2M1nPhi3-TMath::Power(cosP1nPhi1P1nPhi2,2.)-TMath::Power(sinP1nPhi1P1nPhi2,2.) + + 4.*cosP1nPhi1P1nPhi2*(TMath::Power(cosP1nPhi,2.)-TMath::Power(sinP1nPhi,2.)) + + 8.*sinP1nPhi1P1nPhi2*sinP1nPhi*cosP1nPhi + + 8.*two*(TMath::Power(cosP1nPhi,2.)+TMath::Power(sinP1nPhi,2.)) + - 6.*TMath::Power((TMath::Power(cosP1nPhi,2.)+TMath::Power(sinP1nPhi,2.)),2.); + + if (qc4 >= 0) continue; + vnFour = TMath::Power(-qc4, 0.25); + // if (!TMath::IsNaN(vnFour*mult)) + // cumu4->Fill(eta, vnFour, fCumuHist->GetBinContent(0,c,0)); + + // 2-particle differential flow + w2pTwoPrime = fCumuHist->GetBinContent(etaBin, c, kw2two); + w2p = fCumuHist->GetBinContent(etaBin, c, kw2); + mp = fCumuHist->GetBinContent(etaBin, c, kmp); + if (!w2p || !mp) continue; + cosP1nPsi = fCumuHist->GetBinContent(etaBin, c, kpnRe); + sinP1nPsi = fCumuHist->GetBinContent(etaBin, c, kpnIm); + + cosP1nPsi /= mp; + sinP1nPsi /= mp; + twoPrime = w2pTwoPrime / w2p; + qc2Prime = twoPrime - sinP1nPsi*sinP1nPhi - cosP1nPsi*cosP1nPhi; + + vnTwoDiff = qc2Prime / TMath::Sqrt(qc2); + if (!TMath::IsNaN(vnTwoDiff*mp)) cumu2->Fill(eta, (Double_t)c-1., vnTwoDiff, fCumuHist->GetBinContent(0,c,0)); + + // 4-particle differential flow + w4pFourPrime = fCumuHist->GetBinContent(etaBin, c, kw4four); + w4p = fCumuHist->GetBinContent(etaBin, c, kw4); + mpqMult = fCumuHist->GetBinContent(etaBin, c, kmpmq); + if (!w4p || !mpqMult) continue; + cosP1nPsi1P1nPhi2 = fCumuHist->GetBinContent(etaBin, c, kCospsi1phi2); + sinP1nPsi1P1nPhi2 = fCumuHist->GetBinContent(etaBin, c, kSinpsi1phi2); + cosP1nPsi1M1nPhi2M1nPhi3 = fCumuHist->GetBinContent(etaBin, c, kCospsi1phi2phi3m); + sinP1nPsi1M1nPhi2M1nPhi3 = fCumuHist->GetBinContent(etaBin, c, kSinpsi1phi2phi3m); + cosP1nPsi1P1nPhi2M1nPhi3 = fCumuHist->GetBinContent(etaBin, c, kCospsi1phi2phi3p); + sinP1nPsi1P1nPhi2M1nPhi3 = fCumuHist->GetBinContent(etaBin, c, kSinpsi1phi2phi3p); + + cosP1nPsi1P1nPhi2 /= w2p; + sinP1nPsi1P1nPhi2 /= w2p; + cosP1nPsi1M1nPhi2M1nPhi3 /= mpqMult; + sinP1nPsi1M1nPhi2M1nPhi3 /= mpqMult; + cosP1nPsi1P1nPhi2M1nPhi3 /= mpqMult; + sinP1nPsi1P1nPhi2M1nPhi3 /= mpqMult; + + fourPrime = w4pFourPrime / w4p; + + qc4Prime = fourPrime-2.*twoPrime*two + - cosP1nPsi*cosP1nPhi1M1nPhi2M1nPhi3 + + sinP1nPsi*sinP1nPhi1M1nPhi2M1nPhi3 + - cosP1nPhi*cosP1nPsi1M1nPhi2M1nPhi3 + + sinP1nPhi*sinP1nPsi1M1nPhi2M1nPhi3 + - 2.*cosP1nPhi*cosP1nPsi1P1nPhi2M1nPhi3 + - 2.*sinP1nPhi*sinP1nPsi1P1nPhi2M1nPhi3 + - cosP1nPsi1P1nPhi2*cosP1nPhi1P1nPhi2 + - sinP1nPsi1P1nPhi2*sinP1nPhi1P1nPhi2 + + 2.*cosP1nPhi1P1nPhi2*(cosP1nPsi*cosP1nPhi-sinP1nPsi*sinP1nPhi) + + 2.*sinP1nPhi1P1nPhi2*(cosP1nPsi*sinP1nPhi+sinP1nPsi*cosP1nPhi) + + 4.*two*(cosP1nPsi*cosP1nPhi+sinP1nPsi*sinP1nPhi) + + 2.*cosP1nPsi1P1nPhi2*(TMath::Power(cosP1nPhi,2.)-TMath::Power(sinP1nPhi,2.)) + + 4.*sinP1nPsi1P1nPhi2*cosP1nPhi*sinP1nPhi + + 4.*twoPrime*(TMath::Power(cosP1nPhi,2.)+TMath::Power(sinP1nPhi,2.)) + - 6.*(TMath::Power(cosP1nPhi,2.)-TMath::Power(sinP1nPhi,2.)) + * (cosP1nPsi*cosP1nPhi-sinP1nPsi*sinP1nPhi) + - 12.*cosP1nPhi*sinP1nPhi + * (sinP1nPsi*cosP1nPhi+cosP1nPsi*sinP1nPhi); + + vnFourDiff = - qc4Prime / TMath::Power(-qc4, 0.75); + if (!TMath::IsNaN(vnFourDiff*mp)) cumu4->Fill(eta, (Double_t)c-1., vnFourDiff, fCumuHist->GetBinContent(0,c,0)); + } // End of eta loop + // Number of events: + nEv += fCumuHist->GetBinContent(0,c,0); + cumu2->Fill(7., (Double_t)c-1., nEv); + cumu4->Fill(7., (Double_t)c-1., nEv); + } // End of centrality loop return; } -// _____________________________________________________________________ -void AliForwardFlowTaskQC::ProcessPrimary() -{ - // - // If fMC == kTRUE this function takes care of organizing the input - // Monte Carlo data and histograms so AliForwardFlowTaskQC::QCumulants - // can be run on it. - // - - if (fFlowUtil->LoopAODFMDtrrefHits(fAOD)) { - // Run analysis on TrackRefs - for (Int_t n = 1; n <= 6; n++) { - if (fv[n]) - CumulantsMethod("FMDTR", n); - } - } - - if (fFlowUtil->LoopAODSPDtrrefHits(fAOD)) { - // Run analysis on SPD TrackRefs - for (Int_t n = 1; n <= 6; n++) { - if (fv[n]) - CumulantsMethod("SPDTR", n); - } - } - - if (fFlowUtil->LoopAODmc(fAOD, fAddFlow, fAddType, fAddOrder)) { - // Run analysis on MC truth - for (Int_t n = 1; n <= 6; n++) { - if (fv[n]) - CumulantsMethod("MC", n); - } - } - -} //_____________________________________________________________________ // // diff --git a/PWGLF/FORWARD/analysis2/AliForwardFlowTaskQC.h b/PWGLF/FORWARD/analysis2/AliForwardFlowTaskQC.h index b734aa087c3..09d652c0a07 100644 --- a/PWGLF/FORWARD/analysis2/AliForwardFlowTaskQC.h +++ b/PWGLF/FORWARD/analysis2/AliForwardFlowTaskQC.h @@ -4,18 +4,20 @@ #ifndef ALIFORWARDFLOWTASKQC_H #define ALIFORWARDFLOWTASKQC_H /** - * @file AliForwardFlowTaskQC.h - * @author Alexander Hansen alexander.hansen@cern.ch - * @date Fri Mar 25 13:53:00 2011 + * @file AliForwardFlowTaskQC.h + * @author Alexander Hansen + * @date Tue Feb 14 2012 * - * @brief + * @brief * * * @ingroup pwglf_forward_flow */ #include "AliAnalysisTaskSE.h" -#include "AliForwardFlowUtil.h" -class AliAODEvent; +class AliAODForwardMult; +class TH1D; +class TH2D; +class TH3D; /** * @defgroup pwglf_forward_tasks_flow Flow tasks @@ -33,7 +35,6 @@ class AliAODEvent; * @ingroup pwglf_forward_tasks_flow * @ingroup pwglf_forward_flow * - * @todo Add centrality stuff * */ class AliForwardFlowTaskQC : public AliAnalysisTaskSE @@ -59,14 +60,12 @@ public: */ /** * Create output objects - * */ virtual void UserCreateOutputObjects(); /** * Initialize the task - * */ - virtual void Init() {} + virtual void Init() {} /** * Process each event * @@ -80,60 +79,152 @@ public: */ virtual void Terminate(Option_t *option); /* @} */ - /* + /** * Returns the outputlist - * - */ + * + * @return TList* + */ TList* GetOutputList() { return fOutputList; } - /* - * Set Number of @f$ \eta@f$ bins to be used in flow analysis - * - */ - void SetUseNEtaBins(Int_t nbins) { fEtaBins = nbins; } - /* + /** + * Check AODForwardMult object for trigger, vertex and centrality + * returns true if event is OK + * + * @param const aodfm + * + * @return Bool_t + */ + Bool_t AODCheck(const AliAODForwardMult* aodfm); + /* * Set which harmonics to calculate. @f$ v_{1}@f$ to @f$ v_{4}@f$ is * available and calculated as default - * - * @param v1 Do @f$ v_{1}@f$ - * @param v2 Do @f$ v_{2}@f$ - * @param v3 Do @f$ v_{3}@f$ - * @param v4 Do @f$ v_{4}@f$ + * + * @param v1 Do @f$ v_{1}$f$ + * @param v2 Do @f$ v_{2}$f$ + * @param v3 Do @f$ v_{3}$f$ + * @param v4 Do @f$ v_{4}$f$ + * @param v5 Do @f$ v_{5}$f$ + * @param v6 Do @f$ v_{6}$f$ + * + * @return void */ void SetDoHarmonics(Bool_t v1 = kTRUE, Bool_t v2 = kTRUE, Bool_t v3 = kTRUE, Bool_t v4 = kTRUE, Bool_t v5 = kTRUE, Bool_t v6 = kTRUE) { fv[1] = v1; fv[2] = v2; fv[3] = v3; fv[4] = v4; fv[5] = v5; fv[6] = v6;} - /* - * Set string to add flow to MC truth particles - * - * @param type String - */ - void AddFlow(TString type = "") { fAddFlow = type; } - /* - * Set which function fAddFlow should use - * - * @param type of AddFlow - */ - void AddFlowType(Int_t number = 0) { fAddType = number; } - /* - * Set which order of flow to add - * - * @param order Flow order - */ - void AddFlowOrder(Int_t order = 2) { fAddOrder = order; } - /* - * - * Set MC input flag - * - * @param mc MC input - */ - void SetMCinput(Bool_t mc = kTRUE) { fMC = mc; } - /* - * Set number of eta bins to be used in reference flow - * - * @param bins Ref Eta Bins - */ - void SetRefEtaBins(Int_t bins) { fEtaRef = bins; } + /** + * Nested class to handle cumulant calculations in vertex bins + */ + class VertexBin : public TNamed + { + public: + /* + * Constructor + */ + VertexBin(); + /** + * Constructor + * + * @param vLow Min vertex z-coordinate + * @param vHigh Max vertex z-coordinate + * @param moment Flow moment + * @param type Data type (FMD/SPD/FMDTR/SPDTR/MC) + */ + VertexBin(const Int_t vLow, const Int_t vHigh, + const Int_t moment, const Char_t* type); + /** + * Copy constructor + * + * @param o Object to copy from + * + * @return VertexBin + */ + VertexBin(const VertexBin& o); + /** + * Assignment operator + * + * @param VertexBin& + * + * @return VertexBin& + */ + VertexBin& operator=(const VertexBin&); + /** + * Destructor + */ + ~VertexBin(){} + /** + * Add vertex bin output to list + * + * @param list Histograms are added to this list + * + * @return void + */ + virtual void AddOutput(TList* list); + /** + * Check if vertex vZ is in range for this bin + * + * @param vZ z-coordinate of vertex to check + * + * @return Bool_t + */ + Bool_t CheckVertex(Double_t vZ); + /** + * Fill reference and differential flow histograms for analysis + * + * @param dNdetadphi 2D data histogram + * + * @return false if bad event (det. hotspot) + */ + Bool_t FillHists(TH2D dNdetadphi); + /** + * Do cumulants calculations for current event with + * centrality cent + * + * @param cent Event centrality + * + * @return void + */ + void CumulantsAccumulate(Double_t cent); + /** + * Finish cumulants calculations. Takes input and + * output lists in case Terminate is called separately + * + * @param inlist List with input histograms + * @param outlist List with output histograms + * + * @return void + */ + void CumulantsTerminate(TList* inlist, TList* outlist); + + protected: + /* + * Enumeration for ref/diff histograms + */ + enum { kHmult = 1, kHQnRe, kHQnIm, kHQ2nRe, kHQ2nIm }; + /* + * Enumeration for cumulant histogram + */ + enum { kW2Two = 1, kW2, kW4Four, kW4, kQnRe, kQnIm, kM, + kCosphi1phi2, kSinphi1phi2, kCosphi1phi2phi3m, kSinphi1phi2phi3m, kMm1m2, + kw2two, kw2, kw4four, kw4, kpnRe, kpnIm, kmp, + kCospsi1phi2, kSinpsi1phi2, kCospsi1phi2phi3m, kSinpsi1phi2phi3m, + kmpmq, kCospsi1phi2phi3p, kSinpsi1phi2phi3p }; + + const Int_t fMoment; // flow moment + const Int_t fVzMin; // z-vertex min + const Int_t fVzMax; // z-vertex max + const Char_t* fType; // data type + TH2D* fCumuRef; // histogram for reference flow + TH2D* fCumuDiff; // histogram for differential flow + TH3D* fCumuHist; // histogram for cumulants calculations + TH3D* fHistTwoCorr; // Diagnostics histogram for <2> + TH3D* fHistW2; // Diagnostics histogram for w_<2> + TH3D* fHistFourCorr; // Diagnostics histogram for <4> + TH3D* fHistW4; // Diagnostics histogram for w_<4> + TH2D* fdNdedpAcc; // Diagnostics histogram to make acc. maps + + ClassDef(VertexBin, 1); // object for cumulants ananlysis in FMD + }; + protected: /** * Copy constructor @@ -146,40 +237,37 @@ protected: * * @return Reference to this object */ - AliForwardFlowTaskQC& operator=(const AliForwardFlowTaskQC&) { return *this; } + AliForwardFlowTaskQC& operator=(const AliForwardFlowTaskQC&); /* - * if MC information is available do analysis on Monte Carlo truth - * and track references - * + * Initiate vertex bin objects */ - void ProcessPrimary(); - /** - * Calculate Q cumulant - * - * @param type Determines which histograms should be used - * - "FMD" = FMD data histograms - * - "SPD" = SPD data histograms - * - "TR" = track reference histograms - * - "MC" = MC truth histograms - * @param harmonic Which harmonic to calculate + virtual void InitVertexBins(); + /* + * Initiate diagnostics histograms + */ + virtual void InitHists(); + /* + * Analyze event + */ + virtual Bool_t Analyze(); + /* + * Finalize analysis */ - void CumulantsMethod(TString type, Int_t harmonic); + virtual void Finalize(); - TList* fOutputList; // Output list - AliForwardFlowUtil* fFlowUtil;// AliForwardFlowUtil - AliAODEvent* fAOD; // AOD event - Bool_t fMC; // Is MC flags - Int_t fEtaBins; // Number of eta bins in histograms - Int_t fEtaRef; // Number of eta bins for reference flow - Bool_t fv[7]; // Calculate v_{n} flag - TString fAddFlow; // Add flow string - Int_t fAddType; // Add flow type # - Int_t fAddOrder; // Add flow order - Float_t fZvertex; // Z vertex bin - Double_t fCent; // Centrality - + TList fBinsFMD; // list with FMD VertexBin objects + TList fBinsSPD; // list with SPD VertexBin objects + TList* fSumList; // sum list + TList* fOutputList; // Output list + AliAODEvent* fAOD; // AOD event + Bool_t fv[7]; // Calculate v_{n} flag + Float_t fZvertex; // Z vertex bin + Double_t fCent; // Centrality + TH1D* fHistCent; // Diagnostics hist for centrality + TH1D* fHistVertexSel; // Diagnostics hist for selected vertices + TH1D* fHistVertexAll; // Diagnostics hist for all vertices - ClassDef(AliForwardFlowTaskQC, 2); // Analysis task for FMD analysis + ClassDef(AliForwardFlowTaskQC, 1); // Analysis task for FMD analysis }; #endif diff --git a/PWGLF/FORWARD/analysis2/AliForwardFlowUtil.cxx b/PWGLF/FORWARD/analysis2/AliForwardFlowUtil.cxx deleted file mode 100644 index c991f033abd..00000000000 --- a/PWGLF/FORWARD/analysis2/AliForwardFlowUtil.cxx +++ /dev/null @@ -1,433 +0,0 @@ -// -// Class used to handle the input from AODs and put it into histograms -// the Forward Flow tasks can run on. -// It can also add flow to AliAODMCParticles. -// -#include -#include "AliForwardFlowUtil.h" -#include "AliAODCentralMult.h" -#include "AliAODMCParticle.h" -#include "AliAODMCHeader.h" -#include "AliLog.h" -#include "AliAODForwardMult.h" -#include "AliAODEvent.h" -#include "TProfile2D.h" -#include "TGraph.h" -#include "TF1.h" - -ClassImp(AliForwardFlowUtil) - -//_____________________________________________________________________ -AliForwardFlowUtil::AliForwardFlowUtil() : - fList(0), - fCent(-1), - fVertex(1111), - fAliceCent4th(), - fAlicePt2nd4050(), - fAlicePt4th3040(), - fAlicePt4th4050(), - fImpactParToCent() - {} - // - // Default Constructor - // -//_____________________________________________________________________ -AliForwardFlowUtil::AliForwardFlowUtil(TList* outputList) : - fList(0), - fCent(-1), - fVertex(1111), - fAliceCent4th(), - fAlicePt2nd4050(), - fAlicePt4th3040(), - fAlicePt4th4050(), - fImpactParToCent() -{ - // - // Constructor - // - // Parameters: - // TList: list containing histograms for flow analysis - // - fList = outputList; - - Double_t xCumulant4thTPCrefMultTPConlyAll[] = {2.5,7.5,15,25,35,45,55,65}; - Double_t yCumulant4thTPCrefMultTPConlyAll[] = {0.017855,0.032440,0.055818,0.073137,0.083898,0.086690,0.082040,0.077777}; - Int_t nPointsCumulant4thTPCrefMultTPConlyAll = sizeof(xCumulant4thTPCrefMultTPConlyAll)/sizeof(Double_t); - fAliceCent4th = new TGraph(nPointsCumulant4thTPCrefMultTPConlyAll,xCumulant4thTPCrefMultTPConlyAll, - yCumulant4thTPCrefMultTPConlyAll); - - Double_t xCumulant2nd4050ALICE[] = {0.000000,0.250000,0.350000,0.450000,0.550000,0.650000,0.750000,0.850000,0.950000, - 1.100000,1.300000,1.500000,1.700000,1.900000,2.250000,2.750000,3.250000,3.750000,4.500000}; - Double_t yCumulant2nd4050ALICE[] = {0.000000,0.043400,0.059911,0.073516,0.089756,0.105486,0.117391,0.128199,0.138013, - 0.158271,0.177726,0.196383,0.208277,0.216648,0.242954,0.249961,0.240131,0.269006,0.207796}; - Int_t nPointsCumulant2nd4050ALICE = sizeof(xCumulant2nd4050ALICE)/sizeof(Double_t); - fAlicePt2nd4050 = new TGraph(nPointsCumulant2nd4050ALICE,xCumulant2nd4050ALICE,yCumulant2nd4050ALICE); - - Double_t xCumulant4th3040ALICE[] = {0.00000,0.250000,0.350000,0.450000,0.550000,0.650000,0.750000,0.850000,0.950000, - 1.100000,1.300000,1.500000,1.700000,1.900000,2.250000,2.750000,3.250000,3.750000,4.500000, - 5.500000,7.000000,9.000000}; - Double_t yCumulant4th3040ALICE[] = {0.000000,0.037071,0.048566,0.061083,0.070910,0.078831,0.091396,0.102026,0.109691, - 0.124449,0.139819,0.155561,0.165701,0.173678,0.191149,0.202015,0.204540,0.212560,0.195885, - 0.000000,0.000000,0.000000}; - Int_t nPointsCumulant4th3040ALICE = sizeof(xCumulant4th3040ALICE)/sizeof(Double_t); - fAlicePt4th3040 = new TGraph(nPointsCumulant4th3040ALICE,xCumulant4th3040ALICE,yCumulant4th3040ALICE); - - Double_t xCumulant4th4050ALICE[] = {0.000000,0.250000,0.350000,0.450000,0.550000,0.650000,0.750000,0.850000,0.950000, - 1.100000,1.300000,1.500000,1.700000,1.900000,2.250000,2.750000,3.250000,3.750000,4.500000}; - Double_t yCumulant4th4050ALICE[] = {0.000000,0.038646,0.049824,0.066662,0.075856,0.081583,0.099778,0.104674,0.118545, - 0.131874,0.152959,0.155348,0.169751,0.179052,0.178532,0.198851,0.185737,0.239901,0.186098}; - Int_t nPointsCumulant4th4050ALICE = sizeof(xCumulant4th4050ALICE)/sizeof(Double_t); - fAlicePt4th4050 = new TGraph(nPointsCumulant4th4050ALICE, xCumulant4th4050ALICE, yCumulant4th4050ALICE); - - Double_t impactParam[] = {0.,1.75,4.225,5.965,7.765,9.215,10.46,11.565,12.575,13.515,16.679}; - Double_t centrality[] = {0.,2.5,7.5,15,25,35,45,55,65,75,90}; - - Int_t nPoints = sizeof(impactParam)/sizeof(Double_t); - fImpactParToCent = new TGraph(nPoints, impactParam, centrality); - -} -//_____________________________________________________________________ -Bool_t AliForwardFlowUtil::AODCheck(const AliAODForwardMult* aodfm) -{ - // - // Function to check that and AOD event meets the cuts - // - // Parameters: - // AliAODForwardMult: forward mult object with trigger and vertex info - // - fCent = -1; - fVertex = 1111; - - if (!aodfm->IsTriggerBits(AliAODForwardMult::kOffline)) return kFALSE; - fCent = (Double_t)aodfm->GetCentrality(); - if (0. >= fCent || fCent >= 100.) return kFALSE; - TH1D* vertex = (TH1D*)fList->FindObject("VertexAll"); - fVertex = aodfm->GetIpZ(); - vertex->Fill(fVertex); - if (TMath::Abs(fVertex) >= 5.) return kFALSE; - return kTRUE; -} -//_____________________________________________________________________ -Bool_t AliForwardFlowUtil::LoopAODFMD(const AliAODEvent* aodevent) -{ - // - // Loop over AliAODFowardMult object and fill histograms provided - // by flow task. - // - - AliAODForwardMult* aodfmult = static_cast(aodevent->FindListObject("Forward")); - if (!aodfmult) return kFALSE; - if (!AODCheck(aodfmult)) return kFALSE; - - TH2D hdNdetadphi = aodfmult->GetHistogram(); - TH2D* vertex = (TH2D*)fList->FindObject("CoverageVsVertex"); - TH2D* dNdphi = (TH2D*)fList->FindObject("hdNdphiSEFMD"); - TH2D* dNdetadphi = (TH2D*)fList->FindObject("hdNdetadphiSEFMD"); - dNdphi->Reset(); - dNdetadphi->Reset(); - - Double_t mult = 0; - Double_t eta = 0; - Double_t phi = 0; - Double_t weight = 0; - for (Int_t etaBin = 1; etaBin<=hdNdetadphi.GetNbinsX(); etaBin++) { - eta = hdNdetadphi.GetXaxis()->GetBinCenter(etaBin); - if (TMath::Abs(eta) < 1.75) continue; - for (Int_t phiBin = 0; phiBin<=hdNdetadphi.GetNbinsY(); phiBin++) { - phi = hdNdetadphi.GetYaxis()->GetBinCenter(phiBin); - weight = hdNdetadphi.GetBinContent(etaBin, phiBin); - if (phiBin == 0) { - vertex->Fill(eta, fVertex, weight); - continue; - } - dNdetadphi->Fill(eta, phi, weight); - dNdphi->Fill(eta, phi, weight); - dNdphi->Fill(-1.*eta, phi, weight); - mult += weight; - } - } -// fCent = GetCentFromMC(aodevent); -// fCent = 0.5; - - return kTRUE; -} -//_____________________________________________________________________ -Bool_t AliForwardFlowUtil::LoopAODSPD(const AliAODEvent* aodevent) const -{ - // - // Loop over AliAODCentralMult object and fill histograms - // provided by flow task - // - AliAODCentralMult* aodcmult = static_cast(aodevent->FindListObject("CentralClusters")); - if (!aodcmult) return kFALSE; - - TH2D hdNdetadphi = aodcmult->GetHistogram(); - TH2D* vertex = (TH2D*)fList->FindObject("CoverageVsVertex"); - TH2D* dNdetadphi = (TH2D*)fList->FindObject("hdNdetadphiSESPD"); - TH2D* dNdphi = (TH2D*)fList->FindObject("hdNdphiSESPD"); - dNdphi->Reset(); - dNdetadphi->Reset(); - - Double_t eta = 0; - Double_t phi = 0; - Double_t weight = 0; - Double_t mult = 0; - for (Int_t etaBin = 1; etaBin<=hdNdetadphi.GetNbinsX(); etaBin++) { - eta = hdNdetadphi.GetXaxis()->GetBinCenter(etaBin); - if (TMath::Abs(eta) > 1.75) continue; - for (Int_t phiBin = 0; phiBin<=hdNdetadphi.GetNbinsY(); phiBin++) { - phi = hdNdetadphi.GetYaxis()->GetBinCenter(phiBin); - weight = hdNdetadphi.GetBinContent(etaBin, phiBin); - if (phiBin == 0) { - vertex->Fill(eta, fVertex, weight); - continue; - } - dNdetadphi->Fill(eta, phi, weight); - if (TMath::Abs(eta) < 0.5) continue; - dNdphi->Fill(eta, phi, weight); -// dNdphi->Fill(-1.*eta, phi, weight); - mult += weight; - } - } - - return kTRUE; -} -//_____________________________________________________________________ -Bool_t AliForwardFlowUtil::LoopAODFMDtrrefHits(const AliAODEvent* aodevent) const -{ - // - // Loop over AliAODForwardMult object, get MC track ref information - // and fill flow histograms - // - - AliAODForwardMult* aodfmult = static_cast(aodevent->FindListObject("ForwardMC")); - if (!aodfmult) return kFALSE; - - TH2D hdNdetadphi = aodfmult->GetHistogram(); - TH2D* dNdphi = (TH2D*)fList->FindObject("hdNdphiSEFMDTR"); - TH2D* dNdetadphi = (TH2D*)fList->FindObject("hdNdetadphiSEFMDTR"); - dNdphi->Reset(); - dNdetadphi->Reset(); - - Double_t mult = 0; - Double_t eta = 0; - Double_t phi = 0; - Double_t weight = 0; - for(Int_t etaBin = 1; etaBin<=hdNdetadphi.GetNbinsX(); etaBin++) { - eta = hdNdetadphi.GetXaxis()->GetBinCenter(etaBin); - if (TMath::Abs(eta) < 1.75) continue; - for(Int_t phiBin = 1; phiBin<=hdNdetadphi.GetNbinsY(); phiBin++) { - phi = hdNdetadphi.GetYaxis()->GetBinCenter(phiBin); - weight = hdNdetadphi.GetBinContent(etaBin, phiBin); - dNdetadphi->Fill(eta, phi, weight); - dNdphi->Fill(eta, phi, weight); - dNdphi->Fill(-1.*eta, phi, weight); - mult += weight; - } - } - - return kTRUE; -} -//_____________________________________________________________________ -Bool_t AliForwardFlowUtil::LoopAODSPDtrrefHits(const AliAODEvent* aodevent) const -{ - // - // Loop over AliAODCentralMult object and fill histograms - // provided by flow task - // - AliAODCentralMult* aodcmult = static_cast(aodevent->FindListObject("CentralClustersMC")); - if (!aodcmult) return kFALSE; - - TH2D hdNdetadphi = aodcmult->GetHistogram(); - TH2D* dNdphi = (TH2D*)fList->FindObject("hdNdphiSESPDTR"); - TH2D* dNdetadphi = (TH2D*)fList->FindObject("hdNdetadphiSESPDTR"); - dNdphi->Reset(); - dNdetadphi->Reset(); - - Double_t eta = 0; - Double_t phi = 0; - Double_t weight = 0; - Double_t mult = 0; - for (Int_t etaBin = 1; etaBin<=hdNdetadphi.GetNbinsX(); etaBin++) { - eta = hdNdetadphi.GetXaxis()->GetBinCenter(etaBin); - if (TMath::Abs(eta) > 1.75) continue; - for (Int_t phiBin = 1; phiBin<=hdNdetadphi.GetNbinsY(); phiBin++) { - phi = hdNdetadphi.GetYaxis()->GetBinCenter(phiBin); - weight = hdNdetadphi.GetBinContent(etaBin, phiBin); - dNdetadphi->Fill(eta, phi, weight); - dNdphi->Fill(eta, phi, weight); - dNdphi->Fill(-1.*eta, phi, weight); - mult += weight; - } - } - - return kTRUE; -} -//_____________________________________________________________________ -Bool_t AliForwardFlowUtil::LoopAODmc(const AliAODEvent* aodevent, - TString addFlow = "", - Int_t type = 0, - Int_t order = 2) const -{ - // - // Loop over AliAODParticle branch and fill flow histograms - // - - //retreive MC particles from event - TClonesArray* mcArray = (TClonesArray*)aodevent->FindListObject(AliAODMCParticle::StdBranchName()); - if(!mcArray){ -// AliWarning("No MC array found in AOD. Try making it again."); - return kFALSE; - } - - Double_t rp = 0; - AliAODMCHeader* header = dynamic_cast(aodevent->FindListObject(AliAODMCHeader::StdBranchName())); - if (!header) { - AliWarning("No header file found."); - } - else { - rp = header->GetReactionPlaneAngle(); - } - - TH2D* dNdetadphi = (TH2D*)fList->FindObject("hdNdetadphiSEMC"); - TH2D* dNdphi = (TH2D*)fList->FindObject("hdNdphiSEMC"); - TProfile2D* mcTruth = (TProfile2D*)fList->FindObject("pMCTruth"); - dNdphi->Reset(); - dNdetadphi->Reset(); - - Int_t ntracks = mcArray->GetEntriesFast(); - - // Track loop: chek how many particles will be accepted - Double_t mult = 0; - Double_t weight = 0; - for (Int_t it = 0; it < ntracks; it++) { - weight = 1; - AliAODMCParticle* particle = (AliAODMCParticle*)mcArray->At(it); - if (!particle) { - AliError(Form("Could not receive track %d", it)); - continue; - } - if (!particle->IsPhysicalPrimary()) continue; - if (particle->Charge() == 0) continue; - Double_t pT = particle->Pt(); -// if (pT <= 0.2 || pT >= 5.) continue; - Double_t eta = particle->Eta(); - Double_t phi = particle->Phi(); - if (TMath::Abs(eta) < 6.) { - // Add flow if it is in the argument - if (addFlow.Length() > 1) { - if (addFlow.Contains("pt")) - weight *= AddptFlow(pT, type); -// weight *= AddptFlow(pT, 1); - if (addFlow.Contains("pid")) - weight *= AddpidFlow(particle->PdgCode(), 1); - if (addFlow.Contains("eta")) - weight *= AddetaFlow(eta, type); -// weight *= AddetaFlow(eta, 2); - if (addFlow.Contains("cent")) - weight *= fAliceCent4th->Eval(fCent)/fAliceCent4th->Eval(45); - - weight *= 20*2.*TMath::Cos((Double_t)order*(phi-rp)); - weight += 1; - } - dNdphi->Fill(eta, phi, weight); - dNdphi->Fill(-1.*eta, phi, weight); - dNdetadphi->Fill(eta, phi, weight); - mcTruth->Fill(eta, fCent, weight*TMath::Cos(2.*(phi-rp))); - mult += weight; - } - } - - return kTRUE; -} -//_____________________________________________________________________ -Double_t AliForwardFlowUtil::AddptFlow(Double_t pt = 0, Int_t type = 0) const -{ - // - // Add pt dependent flow factor - // - Double_t weight = 0; - - switch(type) - { - case 1: weight = fAlicePt2nd4050->Eval(pt)*0.5+fAlicePt4th4050->Eval(pt)*0.5; - break; - case 2: weight = fAlicePt2nd4050->Eval(pt); - break; - case 3: weight = fAlicePt4th3040->Eval(pt); - break; - case 4: weight = fAlicePt4th4050->Eval(pt); - break; - } - - return weight; -} -//_____________________________________________________________________ -Double_t AliForwardFlowUtil::AddpidFlow(Int_t id = 0, Int_t type = 0) const -{ - // - // Add pid dependent flow factor - // - Double_t weight = 0; - - switch(type) - { - case 1: if (TMath::Abs(id) == 211) // pion flow - weight = 1.; - else if (TMath::Abs(id) == 2212) // proton flow - weight = 1.3; - else - weight = 0.7; - break; - case 2: weight = 1.207; - break; - } - - return weight; -} -//_____________________________________________________________________ -Double_t AliForwardFlowUtil::AddetaFlow(Double_t eta = 0, Int_t type = 0) const -{ - // - // Add eta dependent flow factor - // - Double_t weight = 0; - - TF1 gaus = TF1("gaus", "gaus", -6, 6); - - switch(type) - { - case 1: gaus.SetParameters(0.1, 0., 9); - break; - case 2: gaus.SetParameters(0.1, 0., 3); - break; - case 3: gaus.SetParameters(0.1, 0., 15); - break; - } - - weight = gaus.Eval(eta); - - return weight; -} -//_____________________________________________________________________ -Double_t AliForwardFlowUtil::GetCentFromMC(const AliAODEvent* aodevent) const -{ - // - // Get centrality from MC impact parameter. - // Values taken from: https://twiki.cern.ch/twiki/bin/viewauth/ALICE/CentStudies - // - Double_t cent = -1.; - Double_t b = -1.; - AliAODMCHeader* header = (AliAODMCHeader*)aodevent->FindListObject(AliAODMCHeader::StdBranchName()); - if (!header) return cent; - b = header->GetImpactParameter(); - - cent = fImpactParToCent->Eval(b); - - return cent; -} -//_____________________________________________________________________ -// -// -// EOF - diff --git a/PWGLF/FORWARD/analysis2/AliForwardFlowUtil.h b/PWGLF/FORWARD/analysis2/AliForwardFlowUtil.h deleted file mode 100644 index f69c26a051d..00000000000 --- a/PWGLF/FORWARD/analysis2/AliForwardFlowUtil.h +++ /dev/null @@ -1,176 +0,0 @@ -// -// Class used to handle the input from AODs and put it into histograms -// the Forward Flow tasks can run on -// -#ifndef ALIFORWARDFLOWUTIL_H -#define ALIFORWARDFLOWUTIL_H -/** - * @file AliForwardFlowUtil.h - * @author Alexander Hansen alexander.hansen@cern.ch - * @date Fri Mar 25 13:15:40 2011 - * - * @brief - * - * - * @ingroup pwglf_forward_flow - */ -#include "TNamed.h" -class AliAODForwardMult; -class AliAODEvent; -class TList; -class TGraph; - -/** - * - * Class used to handle the input from AODs and put it into histograms - * the Forward Flow tasks can run on. - * - * @ingroup pwglf_forward_tasks_flow - * @ingroup pwglf_forward_flow - */ -class AliForwardFlowUtil : public TNamed -{ -public: - /** - * Constructor - */ - AliForwardFlowUtil(); - /* - * Constructor - * - * @param l list of histograms for flow analysis - */ - AliForwardFlowUtil(TList* l); - /** - * Check that AOD event meet trigger requirements - * - * @param aodfm Forward multplicity AOD event structure - * - * @return true on success - */ - Bool_t AODCheck(const AliAODForwardMult* aodfm); - /** - * Loop over AliAODForwardMult object and fill flow histograms - * - * @param AODevent AOD event structure - * - * @return true on success - */ - Bool_t LoopAODFMD(const AliAODEvent* AODevent); - /* - * Loop over AliAODCentralMult object and fill flow histograms - * - * @param AODevent AOD event structure - * - * @return true on success - */ - Bool_t LoopAODSPD(const AliAODEvent* AODevent) const; - /** - * Loop over AliAODForwardMult object and fill flow histograms from - * track refs - * - * @param AODevent AOD event structure - * - * @return true on success - */ - Bool_t LoopAODFMDtrrefHits(const AliAODEvent* AODevent) const; - /** - * Loop over AliAODCentralMult object and fill flow histograms from - * track refs - * - * @param AODevent AOD event structure - * - * @return true on success - */ - Bool_t LoopAODSPDtrrefHits(const AliAODEvent* AODevent) const; - /** - * Loop over AliAODMCParticle branch object and fill flow histograms - * add flow if arguments are set - * - * @param AODevent AOD event structure - * @param addFlow What to add flow to - * @param type Type of flow - * @param order Order of added flow - * - * @return true on success - */ - Bool_t LoopAODmc(const AliAODEvent* AODevent, TString addFlow, - Int_t type, Int_t order) const; - /** - * Get centrality from newest processed event - */ - Double_t GetCentrality() const { return fCent; } - /** - * Get z vertex coordinate from newest processed event - */ - Float_t GetVertex() const { return fVertex; } - /** - * Parametrize ALICE data points - */ - -protected: - /* - * Copy constructor - * - * @param o Object to copy from - */ - AliForwardFlowUtil(const AliForwardFlowUtil& o) : TNamed(), - fList(o.fList), - fCent(o.fCent), - fVertex(o.fVertex), - fAliceCent4th(o.fAliceCent4th), - fAlicePt2nd4050(o.fAlicePt2nd4050), - fAlicePt4th3040(o.fAlicePt4th3040), - fAlicePt4th4050(o.fAlicePt4th4050), - fImpactParToCent(o.fImpactParToCent) - {} - /** - * Assignment operator - * - * @return Reference to this object - */ - AliForwardFlowUtil& operator=(const AliForwardFlowUtil&) { return *this; } - /** - * Add pt dependent flow factor - * - * @param Pt @f$ p_T@f$ - * @param type Type of flow - */ - Double_t AddptFlow(Double_t Pt, Int_t type) const; - /** - * Add pid dependent flow factor - * - * @param ID Particle ID - * @param type Type of flow - */ - Double_t AddpidFlow(Int_t ID, Int_t type) const; - /** - * Add eta dependent flow factor - * - * @param Eta @f$\eta@f$ - * @param type Type of flow - */ - Double_t AddetaFlow(Double_t Eta, Int_t type) const; - /** - * Get centrality form MC impact parameter - * - * @param AODevent AOD event structure - */ - Double_t GetCentFromMC(const AliAODEvent* AODevent) const; - - TList* fList; // List of flow histograms - Double_t fCent; // centrality - Float_t fVertex; // z vertex coordinate - TGraph* fAliceCent4th; - TGraph* fAlicePt2nd4050; - TGraph* fAlicePt4th3040; - TGraph* fAlicePt4th4050; - TGraph* fImpactParToCent; - - ClassDef(AliForwardFlowUtil, 1); -}; - -#endif -// Local Variables: -// mode: C++ -// End: diff --git a/PWGLF/FORWARD/analysis2/AliForwardMCFlowTaskQC.cxx b/PWGLF/FORWARD/analysis2/AliForwardMCFlowTaskQC.cxx new file mode 100644 index 00000000000..bfc407e4927 --- /dev/null +++ b/PWGLF/FORWARD/analysis2/AliForwardMCFlowTaskQC.cxx @@ -0,0 +1,441 @@ +// +// Calculate flow on MC data in the forward and central regions using the Q cumulants method. +// +// Inputs: +// - AliAODEvent +// +// Outputs: +// - AnalysisResults.root +// +#include "AliForwardMCFlowTaskQC.h" +#include "AliAODMCParticle.h" +#include "AliAODMCHeader.h" +#include "TGraph.h" +#include "TF1.h" +#include "TProfile2D.h" +#include "AliAODEvent.h" +#include "AliAODForwardMult.h" +#include "AliAODCentralMult.h" + +ClassImp(AliForwardMCFlowTaskQC) + +//_____________________________________________________________________ +AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC() : + AliForwardFlowTaskQC(), + fBinsFMDTR(), // List of FMDTR analysis objects + fBinsSPDTR(), // List of SPDTR analysis objects + fBinsMC(), // List of MC truth analysis objects + fdNdedpMC(), // MC truth d^2N/detadphi histogram + fAliceCent4th(), // Alice QC4 vs. centrality data points + fAlicePt2nd4050(), // Alice QC2 vs. pT data points + fAlicePt4th3040(), // Alice QC4 vs. pT data points + fAlicePt4th4050(), // Alice QC4 vs. pT data points + fImpactParToCent(), // Impact parameter to centrality graph + fAddFlow(0), // Add flow to MC truth + fAddType(0), // Add type of flow to MC truth + fAddOrder(0) // Add order of flow to MC truth + {} + // + // Default Constructor + // +//_____________________________________________________________________ +AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC(const char* name) : + AliForwardFlowTaskQC(name), + fBinsFMDTR(), // List of FMDTR analysis objects + fBinsSPDTR(), // List of SPDTR analysis objects + fBinsMC(), // List of MC truth analysis objects + fdNdedpMC(), // MC truth d^2N/detadphi histogram + fAliceCent4th(), // Alice QC4 vs. centrality data points + fAlicePt2nd4050(), // Alice QC2 vs. pT data points + fAlicePt4th3040(), // Alice QC4 vs. pT data points + fAlicePt4th4050(), // Alice QC4 vs. pT data points + fImpactParToCent(), // Impact parameter to centrality graph + fAddFlow(0), // Add flow to MC truth + fAddType(0), // Add type of flow to MC truth + fAddOrder(0) // Add order of flow to MC truth +{ + // + // Constructor + // Parameters: + // name: Name of task + // + fdNdedpMC = TH2D("fdNdedpMC", "fdNdedpMC", 48, -6., 6., 200, 0., 2.*TMath::Pi()); + fdNdedpMC.Sumw2(); + + // Add parametrizations of ALICE data + Double_t xCumulant4thTPCrefMultTPConlyAll[] = {2.5,7.5,15,25,35,45,55,65}; + Double_t yCumulant4thTPCrefMultTPConlyAll[] = {0.017855,0.032440,0.055818,0.073137,0.083898,0.086690,0.082040,0.077777}; + Int_t nPointsCumulant4thTPCrefMultTPConlyAll = sizeof(xCumulant4thTPCrefMultTPConlyAll)/sizeof(Double_t); + fAliceCent4th = new TGraph(nPointsCumulant4thTPCrefMultTPConlyAll,xCumulant4thTPCrefMultTPConlyAll, + yCumulant4thTPCrefMultTPConlyAll); + + Double_t xCumulant2nd4050ALICE[] = {0.000000,0.250000,0.350000,0.450000,0.550000,0.650000,0.750000,0.850000,0.950000, + 1.100000,1.300000,1.500000,1.700000,1.900000,2.250000,2.750000,3.250000,3.750000,4.500000}; + Double_t yCumulant2nd4050ALICE[] = {0.000000,0.043400,0.059911,0.073516,0.089756,0.105486,0.117391,0.128199,0.138013, + 0.158271,0.177726,0.196383,0.208277,0.216648,0.242954,0.249961,0.240131,0.269006,0.207796}; + Int_t nPointsCumulant2nd4050ALICE = sizeof(xCumulant2nd4050ALICE)/sizeof(Double_t); + fAlicePt2nd4050 = new TGraph(nPointsCumulant2nd4050ALICE,xCumulant2nd4050ALICE,yCumulant2nd4050ALICE); + + Double_t xCumulant4th3040ALICE[] = {0.00000,0.250000,0.350000,0.450000,0.550000,0.650000,0.750000,0.850000,0.950000, + 1.100000,1.300000,1.500000,1.700000,1.900000,2.250000,2.750000,3.250000,3.750000,4.500000, + 5.500000,7.000000,9.000000}; + Double_t yCumulant4th3040ALICE[] = {0.000000,0.037071,0.048566,0.061083,0.070910,0.078831,0.091396,0.102026,0.109691, + 0.124449,0.139819,0.155561,0.165701,0.173678,0.191149,0.202015,0.204540,0.212560,0.195885, + 0.000000,0.000000,0.000000}; + Int_t nPointsCumulant4th3040ALICE = sizeof(xCumulant4th3040ALICE)/sizeof(Double_t); + fAlicePt4th3040 = new TGraph(nPointsCumulant4th3040ALICE,xCumulant4th3040ALICE,yCumulant4th3040ALICE); + + Double_t xCumulant4th4050ALICE[] = {0.000000,0.250000,0.350000,0.450000,0.550000,0.650000,0.750000,0.850000,0.950000, + 1.100000,1.300000,1.500000,1.700000,1.900000,2.250000,2.750000,3.250000,3.750000,4.500000}; + Double_t yCumulant4th4050ALICE[] = {0.000000,0.038646,0.049824,0.066662,0.075856,0.081583,0.099778,0.104674,0.118545, + 0.131874,0.152959,0.155348,0.169751,0.179052,0.178532,0.198851,0.185737,0.239901,0.186098}; + Int_t nPointsCumulant4th4050ALICE = sizeof(xCumulant4th4050ALICE)/sizeof(Double_t); + fAlicePt4th4050 = new TGraph(nPointsCumulant4th4050ALICE, xCumulant4th4050ALICE, yCumulant4th4050ALICE); + + Double_t impactParam[] = {0.,1.75,4.225,5.965,7.765,9.215,10.46,11.565,12.575,13.515,16.679}; + Double_t centrality[] = {0.,2.5,7.5,15,25,35,45,55,65,75,90}; + + Int_t nPoints = sizeof(impactParam)/sizeof(Double_t); + fImpactParToCent = new TGraph(nPoints, impactParam, centrality); + +} +//_____________________________________________________________________ +AliForwardMCFlowTaskQC::AliForwardMCFlowTaskQC(const AliForwardMCFlowTaskQC& o) : + AliForwardFlowTaskQC(o), + fBinsFMDTR(), // List of FMDTR analysis objects + fBinsSPDTR(), // List of SPDTR analysis objects + fBinsMC(), // List of MC truth analysis objects + fdNdedpMC(o.fdNdedpMC), // MC truth d^2N/detadphi histogram + fAliceCent4th(o.fAliceCent4th), // Alice QC4 vs. centrality data points + fAlicePt2nd4050(o.fAlicePt2nd4050), // Alice QC2 vs. pT data points + fAlicePt4th3040(o.fAlicePt4th3040), // Alice QC4 vs. pT data points + fAlicePt4th4050(o.fAlicePt4th4050), // Alice QC4 vs. pT data points + fImpactParToCent(o.fImpactParToCent), // Impact parameter to centrality graph + fAddFlow(o.fAddFlow), // Add flow to MC truth + fAddType(o.fAddType), // Add type of flow to MC truth + fAddOrder(o.fAddOrder) // Add order of flow to MC truth + {} + // + // Copy Constructor + // +//_____________________________________________________________________ +AliForwardMCFlowTaskQC& +AliForwardMCFlowTaskQC::operator=(const AliForwardMCFlowTaskQC& o) +{ + // + // Assignment operator + // Parameters: + // o Object to copy from + // + fdNdedpMC = o.fdNdedpMC; + fAliceCent4th = o.fAliceCent4th; + fAlicePt2nd4050 = o.fAlicePt2nd4050; + fAlicePt4th3040 = o.fAlicePt4th3040; + fAlicePt4th4050 = o.fAlicePt4th4050; + fImpactParToCent = o.fImpactParToCent; + fAddFlow = o.fAddFlow; + fAddType = o.fAddType; + fAddOrder = o.fAddOrder; + return *this; +} +//_____________________________________________________________________ +void AliForwardMCFlowTaskQC::InitVertexBins() +{ + // + // Initiate VertexBin lists + // + AliForwardFlowTaskQC::InitVertexBins(); + + for(Int_t n = 1; n <= 6; n++) { + if (!fv[n]) continue; + for (Int_t v = -10; v < 10; v++) { + fBinsFMDTR.Add(new VertexBin(v, v+1, n, "FMDTR")); + fBinsSPDTR.Add(new VertexBin(v, v+1, n, "SPDTR")); + fBinsMC.Add(new VertexBin(v, v+1, n, "MC")); + } + } + +} +//_____________________________________________________________________ +void AliForwardMCFlowTaskQC::InitHists() +{ + // + // Initiate diagnostics hists and add to outputlist + // + AliForwardFlowTaskQC::InitHists(); + + TIter nextFMDTR(&fBinsFMDTR); + VertexBin* bin = 0; + while ((bin = static_cast(nextFMDTR()))) { + bin->AddOutput(fSumList); + } + TIter nextSPDTR(&fBinsSPDTR); + while ((bin = static_cast(nextSPDTR()))) { + bin->AddOutput(fSumList); + } + TIter nextMC(&fBinsMC); + while ((bin = static_cast(nextMC()))) { + bin->AddOutput(fSumList); + } + +} +//_____________________________________________________________________ +Bool_t AliForwardMCFlowTaskQC::Analyze() +{ + // + // Load FMD and SPD MC objects from aod tree and call the cumulants + // calculation for the correct vertexbin + // + if (!AliForwardFlowTaskQC::Analyze()) return kFALSE; + + // Run analysis on trackrefs from FMD and SPD + AliAODForwardMult* aodfmult = static_cast(fAOD->FindListObject("ForwardMC")); + if (!aodfmult) return kFALSE; + TH2D fmdTRdNdetadphi = aodfmult->GetHistogram(); + + AliAODCentralMult* aodcmult = static_cast(fAOD->FindListObject("CentralClustersMC")); + if (!aodcmult) return kFALSE; + TH2D spdTRdNdetadphi = aodcmult->GetHistogram(); + + TIter nextFMDTR(&fBinsFMDTR); + VertexBin* bin = 0; + while ((bin = static_cast(nextFMDTR()))) { + if (bin->CheckVertex(fZvertex)) { + if (!bin->FillHists(fmdTRdNdetadphi)) return kFALSE; + bin->CumulantsAccumulate(fCent); + } + } + + TIter nextSPDTR(&fBinsSPDTR); + while ((bin = static_cast(nextSPDTR()))) { + if (bin->CheckVertex(fZvertex)) { + if (!bin->FillHists(spdTRdNdetadphi)) return kFALSE; + bin->CumulantsAccumulate(fCent); + } + } + + // Run analysis on MC branch + if (!LoopAODMC()) return kFALSE; + + TIter nextMC(&fBinsMC); + while ((bin = static_cast(nextMC()))) { + if (bin->CheckVertex(fZvertex)) { + if (!bin->FillHists(fdNdedpMC)) return kFALSE; + bin->CumulantsAccumulate(fCent); + } + } + + return kTRUE; +} +//_____________________________________________________________________ +void AliForwardMCFlowTaskQC::Finalize() +{ + // + // Finalize command, called by Terminate() + // + AliForwardFlowTaskQC::Finalize(); + + TIter nextFMDTR(&fBinsFMDTR); + VertexBin* bin = 0; + while ((bin = static_cast(nextFMDTR()))) { + bin->CumulantsTerminate(fSumList, fOutputList); + } + TIter nextSPDTR(&fBinsSPDTR); + while ((bin = static_cast(nextSPDTR()))) { + bin->CumulantsTerminate(fSumList, fOutputList); + } + TIter nextMC(&fBinsMC); + while ((bin = static_cast(nextMC()))) { + bin->CumulantsTerminate(fSumList, fOutputList); + } + + TProfile2D* fmdHist = 0; + TProfile2D* spdHist = 0; + TProfile2D* mcHist = 0; + TList* correctionList = new TList(); + correctionList->SetName("CorrectionList"); +// fOutputList->Add(correctionList); + + for (Int_t i = 2; i <= 4; i += 2) { + for (Int_t n = 1; n <= 6; n++) { + if (!fv[n]) continue; + fmdHist = (TProfile2D*)fOutputList->FindObject(Form("FMDQC%d_v%d_unCorr", i, n)) + ->Clone(Form("FMDQC%d_v%d_Correction", i, n)); + spdHist = (TProfile2D*)fOutputList->FindObject(Form("SPDQC%d_v%d_unCorr", i, n)) + ->Clone(Form("SPDQC%d_v%d_Correction", i, n)); + mcHist = (TProfile2D*)fOutputList->FindObject(Form("MCQC%d_v%d_unCorr", i, n)); + + if (!fmdHist || !spdHist || !mcHist) { + AliError(Form("Histogram missing, correction object not created for v%d", n)); + continue; + } + + fmdHist->Divide(mcHist); + spdHist->Divide(mcHist); + fmdHist->SetTitle(Form("FMD QC{%d} v_{%d} Correction Object", i, n)); + fmdHist->SetTitle(Form("SPD QC{%d} v_{%d} Correction Object", i, n)); + +// correctionList->Add(fmdHist); +// correctionList->Add(spdHist); + } + } + +} +//_____________________________________________________________________ +Bool_t AliForwardMCFlowTaskQC::LoopAODMC() +{ + // + // Loop over AliAODParticle branch and fill d^2N/detadphi-histograms. + // Add flow if set to do so in AddTask function + fdNdedpMC.Reset(); + + //retreive MC particles from event + TClonesArray* mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName()); + if(!mcArray){ +// AliWarning("No MC array found in AOD. Try making it again."); + return kFALSE; + } + + Double_t rp = 0; + AliAODMCHeader* header = dynamic_cast(fAOD->FindListObject(AliAODMCHeader::StdBranchName())); + if (!header) { + AliWarning("No header file found."); + } + else { + rp = header->GetReactionPlaneAngle(); + } + + Int_t ntracks = mcArray->GetEntriesFast(); + + // Track loop: chek how many particles will be accepted + Double_t weight = 0; + for (Int_t it = 0; it < ntracks; it++) { + weight = 1; + AliAODMCParticle* particle = (AliAODMCParticle*)mcArray->At(it); + if (!particle) { + AliError(Form("Could not receive track %d", it)); + continue; + } + if (!particle->IsPhysicalPrimary()) continue; + if (particle->Charge() == 0) continue; + Double_t pT = particle->Pt(); + Double_t eta = particle->Eta(); + Double_t phi = particle->Phi(); + if (TMath::Abs(eta) < 6.) { + // Add flow if it is in the argument + if (fAddFlow.Length() > 1) { + if (fAddFlow.Contains("pt")) + weight *= AddptFlow(pT); + if (fAddFlow.Contains("pid")) + weight *= AddpidFlow(particle->PdgCode()); + if (fAddFlow.Contains("eta")) + weight *= AddetaFlow(eta); + if (fAddFlow.Contains("cent")) + weight *= fAliceCent4th->Eval(fCent)/fAliceCent4th->Eval(45); + + weight *= 20*2.*TMath::Cos((Double_t)fAddOrder*(phi-rp)); + weight += 1; + } + fdNdedpMC.Fill(eta, phi, weight); + } + } + + return kTRUE; +} +//_____________________________________________________________________ +Double_t AliForwardMCFlowTaskQC::AddptFlow(Double_t pt = 0) const +{ + // + // Add pt dependent flow factor + // Parameters: + // pt: pT parametrization to use + // + Double_t weight = 0; + + switch(fAddType) + { + case 1: weight = fAlicePt2nd4050->Eval(pt)*0.5+fAlicePt4th4050->Eval(pt)*0.5; + break; + case 2: weight = fAlicePt2nd4050->Eval(pt); + break; + case 3: weight = fAlicePt4th3040->Eval(pt); + break; + case 4: weight = fAlicePt4th4050->Eval(pt); + break; + } + + return weight; +} +//_____________________________________________________________________ +Double_t AliForwardMCFlowTaskQC::AddpidFlow(Int_t id = 0) const +{ + // + // Add pid dependent flow factor + // Parameters: + // id: choose PID dependent setup + // + Double_t weight = 0; + + switch(fAddType) + { + case 1: if (TMath::Abs(id) == 211) // pion flow + weight = 1.3; + else if (TMath::Abs(id) == 2212) // proton flow + weight = 1.; + else + weight = 0.7; + break; + case 2: weight = 1.207; + break; + } + + return weight; +} +//_____________________________________________________________________ +Double_t AliForwardMCFlowTaskQC::AddetaFlow(Double_t eta = 0) const +{ + // + // Add eta dependent flow factor + // Parameters: + // eta: choose v_n(eta) shape + // + Double_t weight = 0; + + TF1 gaus = TF1("gaus", "gaus", -6, 6); + + switch(fAddType) + { + case 1: gaus.SetParameters(0.1, 0., 9); + break; + case 2: gaus.SetParameters(0.1, 0., 3); + break; + case 3: gaus.SetParameters(0.1, 0., 15); + break; + } + + weight = gaus.Eval(eta); + + return weight; +} +//_____________________________________________________________________ +Double_t AliForwardMCFlowTaskQC::GetCentFromB() const +{ + // + // Get centrality from MC impact parameter. + // Values taken from: https://twiki.cern.ch/twiki/bin/viewauth/ALICE/CentStudies + // + Double_t cent = -1.; + Double_t b = -1.; + AliAODMCHeader* header = (AliAODMCHeader*)fAOD->FindListObject(AliAODMCHeader::StdBranchName()); + if (!header) return cent; + b = header->GetImpactParameter(); + + cent = fImpactParToCent->Eval(b); + + return cent; +} +//_____________________________________________________________________ +// +// +// EOF + diff --git a/PWGLF/FORWARD/analysis2/AliForwardMCFlowTaskQC.h b/PWGLF/FORWARD/analysis2/AliForwardMCFlowTaskQC.h new file mode 100644 index 00000000000..84cf9496f8b --- /dev/null +++ b/PWGLF/FORWARD/analysis2/AliForwardMCFlowTaskQC.h @@ -0,0 +1,159 @@ +// +// Calculate the flow in the forward regions using the Q cumulants method +// +#ifndef ALIFORWARDMCFLOWTASKQC_H +#define ALIFORWARDMCFLOWTASKQC_H +/** + * @file AliForwardMCFlowTaskQC.h + * @author Alexander Hansen alexander.hansen@cern.ch + * @date Tue Feb 14 2012 + * + * @brief + * + * + * @ingroup pwglf_forward_flow + */ +#include "AliForwardFlowTaskQC.h" +#include +class TGraph; + + /** + * @defgroup pwg2_forward_tasks_flow Flow tasks + * @ingroup pwg2_forward_tasks + */ +/** + * Calculate the flow in the forward regions using the Q cumulants method + * + * @par Inputs: + * - AliAODEvent + * + * Outputs: + * - AnalysisResults.root + * + * @ingroup pwg2_forward_tasks_flow + * @ingroup pwg2_forward_flow + * + * + */ +class AliForwardMCFlowTaskQC : public AliForwardFlowTaskQC +{ +public: + /** + * Constructor + */ + AliForwardMCFlowTaskQC(); + /* + * Constructor + * + * @param name Name of task + */ + AliForwardMCFlowTaskQC(const char* name); + /** + * Destructor + */ + virtual ~AliForwardMCFlowTaskQC() {} + /** + * @{ + * @name Task interface methods + */ + /** + * Loop over AliAODMCParticle branch object and fill d^2N/detadphi histograms + * add flow if arguments are set + * + * @return true on success + */ + Bool_t LoopAODMC(); + /* + * Set string to add flow to MC truth particles + * + * @param type String + */ + void AddFlow(TString type = "") { fAddFlow = type; } + /* + * Set which function fAddFlow should use + * + * @param type of AddFlow + */ + void AddFlowType(Int_t number = 0) { fAddType = number; } + /* + * Set which order of flow to add + * + * @param order Flow order + */ + void AddFlowOrder(Int_t order = 2) { fAddOrder = order; } + +protected: + /* + * Copy constructor + * + * @param o Object to copy from + */ + AliForwardMCFlowTaskQC(const AliForwardMCFlowTaskQC& o); + /** + * Assignment operator + * + * @return Reference to this object + */ + AliForwardMCFlowTaskQC& operator=(const AliForwardMCFlowTaskQC& o); + /* + * Initiate vertex bin objects + */ + void InitVertexBins(); + /* + * Initiate diagnostics histograms + */ + void InitHists(); + /* + * Analyze event + */ + Bool_t Analyze(); + /* + * Finalize analysis + */ + void Finalize(); + /** + * Add pt dependent flow factor + * + * @param Pt @f$ p_T@f$ + * @param type Type of flow + */ + Double_t AddptFlow(Double_t pt) const; + /** + * Add pid dependent flow factor + * + * @param ID Particle ID + * @param type Type of flow + */ + Double_t AddpidFlow(Int_t id) const; + /** + * Add eta dependent flow factor + * + * @param Eta @f$\eta@f$ + * @param type Type of flow + */ + Double_t AddetaFlow(Double_t eta) const; + /** + * Get centrality form MC impact parameter + */ + Double_t GetCentFromB() const; + + TList fBinsFMDTR; // List with FMDTR VertexBin objects + TList fBinsSPDTR; // List with SPDTR VertexBin objects + TList fBinsMC; // List with MC VertexBin objects + TH2D fdNdedpMC; // d^2N/detadphi MC truth histogram + TGraph* fAliceCent4th; // Parametrization of ALICE QC4 vs. cent. data + TGraph* fAlicePt2nd4050; // Parametrization of ALICE QC2 vs. pT data + TGraph* fAlicePt4th3040; // Parametrization of ALICE QC4 vs. pT data + TGraph* fAlicePt4th4050; // Parametrization of ALICE QC4 vs. pT data + TGraph* fImpactParToCent; // Parametrization of b to centrality datapoints + TString fAddFlow; // Add flow string + Int_t fAddType; // Add flow type # + Int_t fAddOrder; // Add flow order + + ClassDef(AliForwardMCFlowTaskQC, 1); // FMD MC analysis task +}; + +#endif +// Local Variables: +// mode: C++ +// End: diff --git a/PWGLF/FORWARD/analysis2/AliForwardMCMultiplicityTask.cxx b/PWGLF/FORWARD/analysis2/AliForwardMCMultiplicityTask.cxx index 1fcad32d8de..7b1e6bc0d74 100644 --- a/PWGLF/FORWARD/analysis2/AliForwardMCMultiplicityTask.cxx +++ b/PWGLF/FORWARD/analysis2/AliForwardMCMultiplicityTask.cxx @@ -35,6 +35,7 @@ AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask() fESDFMD(), fHistos(), fAODFMD(), + fAODEP(), fMCESDFMD(), fMCHistos(), fMCAODFMD(), @@ -46,6 +47,7 @@ AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask() fDensityCalculator(), fCorrections(), fHistCollector(), + fEventPlaneFinder(), fList(0) { // @@ -60,6 +62,7 @@ AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const char* name) fESDFMD(), fHistos(), fAODFMD(kFALSE), + fAODEP(kFALSE), fMCESDFMD(), fMCHistos(), fMCAODFMD(kTRUE), @@ -71,6 +74,7 @@ AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const char* name) fDensityCalculator("density"), fCorrections("corrections"), fHistCollector("collector"), + fEventPlaneFinder("eventplane"), fList(0) { // @@ -89,6 +93,7 @@ AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const AliForwardMCMul fESDFMD(o.fESDFMD), fHistos(o.fHistos), fAODFMD(o.fAODFMD), + fAODEP(o.fAODEP), fMCESDFMD(o.fMCESDFMD), fMCHistos(o.fMCHistos), fMCAODFMD(o.fMCAODFMD), @@ -100,6 +105,7 @@ AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const AliForwardMCMul fDensityCalculator(o.fDensityCalculator), fCorrections(o.fCorrections), fHistCollector(o.fHistCollector), + fEventPlaneFinder(o.fEventPlaneFinder), fList(o.fList) { // @@ -133,8 +139,10 @@ AliForwardMCMultiplicityTask::operator=(const AliForwardMCMultiplicityTask& o) fDensityCalculator = o.fDensityCalculator; fCorrections = o.fCorrections; fHistCollector = o.fHistCollector; + fEventPlaneFinder = o.fEventPlaneFinder; fHistos = o.fHistos; fAODFMD = o.fAODFMD; + fAODEP = o.fAODEP; fMCHistos = o.fMCHistos; fMCAODFMD = o.fMCAODFMD; fRingSums = o.fRingSums; @@ -160,6 +168,7 @@ AliForwardMCMultiplicityTask::SetDebug(Int_t dbg) fDensityCalculator.SetDebug(dbg); fCorrections.SetDebug(dbg); fHistCollector.SetDebug(dbg); + fEventPlaneFinder.SetDebug(dbg); } //____________________________________________________________________ void @@ -184,6 +193,7 @@ AliForwardMCMultiplicityTask::InitializeSubs() fHistos.Init(*pe); fAODFMD.Init(*pe); + fAODEP.Init(*pe); fMCHistos.Init(*pe); fMCAODFMD.Init(*pe); fRingSums.Init(*pe); @@ -234,6 +244,7 @@ AliForwardMCMultiplicityTask::InitializeSubs() fDensityCalculator.Init(*pe); fCorrections.Init(*pe); fHistCollector.Init(*pv,*pe); + fEventPlaneFinder.Init(*pe); this->Print(); } @@ -261,6 +272,9 @@ AliForwardMCMultiplicityTask::UserCreateOutputObjects() TObject* mcobj = &fMCAODFMD; ah->AddBranch("AliAODForwardMult", &mcobj); + TObject* epobj = &fAODFMD; + ah->AddBranch("AliAODForwardEP", &epobj); + fPrimary = new TH2D("primary", "MC Primaries", 200, -4, 6, 20, 0, 2*TMath::Pi()); fPrimary->SetXTitle("#eta"); @@ -276,6 +290,7 @@ AliForwardMCMultiplicityTask::UserCreateOutputObjects() fDensityCalculator.DefineOutput(fList); fCorrections.DefineOutput(fList); fHistCollector.DefineOutput(fList); + fEventPlaneFinder.DefineOutput(fList); PostData(1, fList); } @@ -302,6 +317,7 @@ AliForwardMCMultiplicityTask::UserExec(Option_t*) fHistos.Clear(); fESDFMD.Clear(); fAODFMD.Clear(); + fAODEP.Clear(); fMCHistos.Clear(); fMCESDFMD.Clear(); fMCAODFMD.Clear(); @@ -396,6 +412,11 @@ AliForwardMCMultiplicityTask::UserExec(Option_t*) } fDensityCalculator.CompareResults(fHistos, fMCHistos); + if (fEventInspector.GetCollisionSystem() == AliFMDEventInspector::kPbPb) { + if (!fEventPlaneFinder.FindEventplane(esd, fAODEP, &(fAODFMD.GetHistogram()) , &fHistos)) + AliWarning("Eventplane finder failed!"); + } + // Do the secondary and other corrections. if (!fCorrections.Correct(fHistos, ivz)) { AliWarning("Corrections failed"); diff --git a/PWGLF/FORWARD/analysis2/AliForwardMCMultiplicityTask.h b/PWGLF/FORWARD/analysis2/AliForwardMCMultiplicityTask.h index 65f21aa7f86..ca2d20a22eb 100644 --- a/PWGLF/FORWARD/analysis2/AliForwardMCMultiplicityTask.h +++ b/PWGLF/FORWARD/analysis2/AliForwardMCMultiplicityTask.h @@ -21,7 +21,9 @@ #include "AliFMDMCCorrector.h" #include "AliFMDHistCollector.h" #include "AliAODForwardMult.h" +#include "AliAODForwardEP.h" #include "AliFMDEnergyFitter.h" +#include "AliFMDEventPlaneFinder.h" #include class AliESDEvent; class TH2D; @@ -169,6 +171,18 @@ public: * @return Reference to AliFMDHistCollector object */ const AliFMDHistCollector& GetHistCollector() const { return fHistCollector; } + /** + * Get reference to the EventPlaneFinder algorithm + * + * @return Reference to AliFMDEventPlaneFinder object + */ + AliFMDEventPlaneFinder& GetEventPlaneFinder() { return fEventPlaneFinder; } + /** + * Get reference to the EventPlaneFinder algorithm + * + * @return Reference to AliFMDEventPlaneFinder object + */ + const AliFMDEventPlaneFinder& GetEventPlaneFinder() const { return fEventPlaneFinder; } /** * @} */ @@ -189,6 +203,7 @@ protected: AliESDFMD fESDFMD; // Sharing corrected ESD object AliForwardUtil::Histos fHistos; // Cache histograms AliAODForwardMult fAODFMD; // Output object + AliAODForwardEP fAODEP; // Output object AliESDFMD fMCESDFMD; // MC 'Sharing corrected' ESD object AliForwardUtil::Histos fMCHistos; // MC Cache histograms AliAODForwardMult fMCAODFMD; // MC Output object @@ -201,6 +216,7 @@ protected: AliFMDMCDensityCalculator fDensityCalculator; // Algorithm AliFMDMCCorrector fCorrections; // Algorithm AliFMDHistCollector fHistCollector; // Algorithm + AliFMDEventPlaneFinder fEventPlaneFinder; // Algorithm TList* fList; // Output list diff --git a/PWGLF/FORWARD/analysis2/AliForwardMultiplicityBase.cxx b/PWGLF/FORWARD/analysis2/AliForwardMultiplicityBase.cxx index 34218555c65..3f2cdf7b10c 100644 --- a/PWGLF/FORWARD/analysis2/AliForwardMultiplicityBase.cxx +++ b/PWGLF/FORWARD/analysis2/AliForwardMultiplicityBase.cxx @@ -24,6 +24,7 @@ #include "AliFMDDensityCalculator.h" #include "AliFMDCorrector.h" #include "AliFMDHistCollector.h" +#include "AliFMDEventPlaneFinder.h" #include "AliESDEvent.h" #include #include @@ -191,6 +192,7 @@ AliForwardMultiplicityBase::GetESDEvent() fFirstEvent = false; + GetEventPlaneFinder().SetRunNumber(esd->GetRunNumber()); InitializeSubs(); } return esd; @@ -315,6 +317,7 @@ AliForwardMultiplicityBase::Print(Option_t* option) const GetDensityCalculator().Print(option); GetCorrections() .Print(option); GetHistCollector() .Print(option); + GetEventPlaneFinder() .Print(option); gROOT->DecreaseDirLevel(); } diff --git a/PWGLF/FORWARD/analysis2/AliForwardMultiplicityBase.h b/PWGLF/FORWARD/analysis2/AliForwardMultiplicityBase.h index e78a615164c..4a25f5ea4f4 100644 --- a/PWGLF/FORWARD/analysis2/AliForwardMultiplicityBase.h +++ b/PWGLF/FORWARD/analysis2/AliForwardMultiplicityBase.h @@ -22,6 +22,7 @@ class AliFMDDensityCalculator; class AliFMDCorrector; class AliFMDHistCollector; class AliForwardCorrectionManager; +class AliFMDEventPlaneFinder; class AliESDEvent; class TH2D; class TList; @@ -234,6 +235,19 @@ public: * @return Reference to AliFMDHistCollector object */ virtual const AliFMDHistCollector& GetHistCollector() const = 0; + /** + * Get reference to the EventPlaneFinder algorithm + * + * @return Reference to AliFMDEventPlaneFinder object + */ + virtual AliFMDEventPlaneFinder& GetEventPlaneFinder() = 0; + /** + * Get reference to the EventPlaneFinder algorithm + * + * @return Reference to AliFMDEventPlaneFinder object + */ + virtual const AliFMDEventPlaneFinder& GetEventPlaneFinder() const = 0; + /** * @} */ diff --git a/PWGLF/FORWARD/analysis2/AliForwardMultiplicityTask.cxx b/PWGLF/FORWARD/analysis2/AliForwardMultiplicityTask.cxx index b7508e60aa9..20b3052299a 100644 --- a/PWGLF/FORWARD/analysis2/AliForwardMultiplicityTask.cxx +++ b/PWGLF/FORWARD/analysis2/AliForwardMultiplicityTask.cxx @@ -34,12 +34,14 @@ AliForwardMultiplicityTask::AliForwardMultiplicityTask() fESDFMD(), fHistos(), fAODFMD(), + fAODEP(), fRingSums(), fEventInspector(), fSharingFilter(), fDensityCalculator(), fCorrections(), fHistCollector(), + fEventPlaneFinder(), fList(0) { // @@ -54,12 +56,14 @@ AliForwardMultiplicityTask::AliForwardMultiplicityTask(const char* name) fESDFMD(), fHistos(), fAODFMD(false), + fAODEP(false), fRingSums(), fEventInspector("event"), fSharingFilter("sharing"), fDensityCalculator("density"), fCorrections("corrections"), fHistCollector("collector"), + fEventPlaneFinder("eventplane"), fList(0) { // @@ -78,12 +82,14 @@ AliForwardMultiplicityTask::AliForwardMultiplicityTask(const AliForwardMultiplic fESDFMD(o.fESDFMD), fHistos(o.fHistos), fAODFMD(o.fAODFMD), + fAODEP(o.fAODEP), fRingSums(o.fRingSums), fEventInspector(o.fEventInspector), fSharingFilter(o.fSharingFilter), fDensityCalculator(o.fDensityCalculator), fCorrections(o.fCorrections), fHistCollector(o.fHistCollector), + fEventPlaneFinder(o.fEventPlaneFinder), fList(o.fList) { // @@ -117,8 +123,10 @@ AliForwardMultiplicityTask::operator=(const AliForwardMultiplicityTask& o) fDensityCalculator = o.fDensityCalculator; fCorrections = o.fCorrections; fHistCollector = o.fHistCollector; + fEventPlaneFinder = o.fEventPlaneFinder; fHistos = o.fHistos; fAODFMD = o.fAODFMD; + fAODEP = o.fAODEP; fRingSums = o.fRingSums; fList = o.fList; @@ -140,6 +148,7 @@ AliForwardMultiplicityTask::SetDebug(Int_t dbg) fDensityCalculator.SetDebug(dbg); fCorrections.SetDebug(dbg); fHistCollector.SetDebug(dbg); + fEventPlaneFinder.SetDebug(dbg); } //____________________________________________________________________ @@ -157,6 +166,7 @@ AliForwardMultiplicityTask::InitializeSubs() fHistos.Init(*pe); fAODFMD.Init(*pe); + fAODEP.Init(*pe); fRingSums.Init(*pe); fHData = static_cast(fAODFMD.GetHistogram().Clone("d2Ndetadphi")); @@ -185,6 +195,7 @@ AliForwardMultiplicityTask::InitializeSubs() fDensityCalculator.Init(*pe); fCorrections.Init(*pe); fHistCollector.Init(*pv,*pe); + fEventPlaneFinder.Init(*pe); this->Print(); } @@ -208,12 +219,15 @@ AliForwardMultiplicityTask::UserCreateOutputObjects() TObject* obj = &fAODFMD; ah->AddBranch("AliAODForwardMult", &obj); + TObject* epobj = &fAODEP; + ah->AddBranch("AliAODForwardEP", &epobj); fEventInspector.DefineOutput(fList); fSharingFilter.DefineOutput(fList); fDensityCalculator.DefineOutput(fList); fCorrections.DefineOutput(fList); fHistCollector.DefineOutput(fList); + fEventPlaneFinder.DefineOutput(fList); PostData(1, fList); } @@ -237,6 +251,7 @@ AliForwardMultiplicityTask::UserExec(Option_t*) fHistos.Clear(); fESDFMD.Clear(); fAODFMD.Clear(); + fAODEP.Clear(); Bool_t lowFlux = kFALSE; UInt_t triggers = 0; @@ -257,7 +272,7 @@ AliForwardMultiplicityTask::UserExec(Option_t*) fAODFMD.SetCentrality(cent); fAODFMD.SetNClusters(nClusters); MarkEventForStore(); - + if (found & AliFMDEventInspector::kNoSPD) return; if (found & AliFMDEventInspector::kNoFMD) return; if (found & AliFMDEventInspector::kNoVertex) return; @@ -278,13 +293,18 @@ AliForwardMultiplicityTask::UserExec(Option_t*) AliWarning("Sharing filter failed!"); return; } - + // Calculate the inclusive charged particle density if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux, cent)) { // if (!fDensityCalculator.Calculate(*esdFMD, fHistos, ivz, lowFlux)) { AliWarning("Density calculator failed!"); return; } + + if (fEventInspector.GetCollisionSystem() == AliFMDEventInspector::kPbPb) { + if (!fEventPlaneFinder.FindEventplane(esd, fAODEP, &(fAODFMD.GetHistogram()), &fHistos)) + AliWarning("Eventplane finder failed!"); + } // Do the secondary and other corrections. if (!fCorrections.Correct(fHistos, ivz)) { diff --git a/PWGLF/FORWARD/analysis2/AliForwardMultiplicityTask.h b/PWGLF/FORWARD/analysis2/AliForwardMultiplicityTask.h index 4792d2f6878..bcef9606f5e 100644 --- a/PWGLF/FORWARD/analysis2/AliForwardMultiplicityTask.h +++ b/PWGLF/FORWARD/analysis2/AliForwardMultiplicityTask.h @@ -21,7 +21,9 @@ #include "AliFMDCorrector.h" #include "AliFMDHistCollector.h" #include "AliAODForwardMult.h" +#include "AliAODForwardEP.h" #include "AliFMDEnergyFitter.h" +#include "AliFMDEventPlaneFinder.h" #include class AliESDEvent; class TH2D; @@ -159,9 +161,22 @@ public: * @return Reference to AliFMDHistCollector object */ const AliFMDHistCollector& GetHistCollector() const { return fHistCollector; } + /** + * Get reference to the EventPlaneFinder algorithm + * + * @return Reference to AliFMDEventPlaneFinder object + */ + AliFMDEventPlaneFinder& GetEventPlaneFinder() { return fEventPlaneFinder; } + /** + * Get reference to the EventPlaneFinder algorithm + * + * @return Reference to AliFMDEventPlaneFinder object + */ + const AliFMDEventPlaneFinder& GetEventPlaneFinder() const { return fEventPlaneFinder; } /** * @} */ + /** * Set debug level * @@ -179,6 +194,7 @@ protected: AliESDFMD fESDFMD; // Sharing corrected ESD object AliForwardUtil::Histos fHistos; // Cache histograms AliAODForwardMult fAODFMD; // Output object + AliAODForwardEP fAODEP; // Output object AliForwardUtil::Histos fRingSums; // Cache histograms AliFMDEventInspector fEventInspector; // Algorithm @@ -186,6 +202,7 @@ protected: AliFMDDensityCalculator fDensityCalculator; // Algorithm AliFMDCorrector fCorrections; // Algorithm AliFMDHistCollector fHistCollector; // Algorithm + AliFMDEventPlaneFinder fEventPlaneFinder; // Algorithm TList* fList; // Output list diff --git a/PWGLF/FORWARD/analysis2/AliSPDMCTrackDensity.cxx b/PWGLF/FORWARD/analysis2/AliSPDMCTrackDensity.cxx index ad072af5532..e398497b8b6 100644 --- a/PWGLF/FORWARD/analysis2/AliSPDMCTrackDensity.cxx +++ b/PWGLF/FORWARD/analysis2/AliSPDMCTrackDensity.cxx @@ -10,6 +10,9 @@ #include #include #include +#include "AliGenHijingEventHeader.h" +#include +#include //____________________________________________________________________ AliSPDMCTrackDensity::AliSPDMCTrackDensity() @@ -117,7 +120,8 @@ AliSPDMCTrackDensity::StoreParticle(AliMCParticle* particle, const AliMCParticle* mother, Int_t refNo, Double_t vz, - TH2D& output) const + TH2D& output, + Double_t w) const { // Store a particle. if (refNo < 0) return; @@ -136,14 +140,14 @@ AliSPDMCTrackDensity::StoreParticle(AliMCParticle* particle, Double_t et = -TMath::Log(TMath::Tan(th/2)); Double_t ph = TMath::ATan2(y,x); if (ph < 0) ph += 2*TMath::Pi(); - output.Fill(et,ph); + output.Fill(et,ph,w); const AliMCParticle* mp = (mother ? mother : particle); Double_t dEta = mp->Eta() - et; Double_t dPhi = (mp->Phi() - ph) * 180 / TMath::Pi(); if (dPhi > 180) dPhi -= 360; if (dPhi < -180) dPhi += 360; - fBinFlow->Fill(dEta, dPhi); + fBinFlow->Fill(dEta, dPhi,w); } @@ -166,6 +170,7 @@ AliSPDMCTrackDensity::GetMother(Int_t iTr, return 0; } +// #define USE_FLOW_WEIGHTS //____________________________________________________________________ Bool_t AliSPDMCTrackDensity::Calculate(const AliMCEvent& event, @@ -188,6 +193,13 @@ AliSPDMCTrackDensity::Calculate(const AliMCEvent& event, // True on succes, false otherwise // +#ifdef USE_FLOW_WEIGHTS + AliGenHijingEventHeader* hd = dynamic_cast + (event.GenEventHeader()); + Double_t rp = (hd ? hd->ReactionPlaneAngle() : 0.); + Double_t b = (hd ? hd->ImpactParameter() : -1 ); +#endif + AliStack* stack = const_cast(event).Stack(); Int_t nTracks = stack->GetNtrack(); Int_t nPrim = stack->GetNtrack(); @@ -243,13 +255,127 @@ AliSPDMCTrackDensity::Calculate(const AliMCEvent& event, // Only fill first reference if (nRef == 1) { const AliMCParticle* mother = GetMother(iTr, event); - StoreParticle(particle, mother, iTrRef, vz, output); +#ifdef USE_FLOW_WEIGHTS + Double_t phi = (mother ? mother->Phi() : particle->Phi()); + Double_t eta = (mother ? mother->Eta() : particle->Eta()); + Double_t pt = (mother ? mother->Pt() : particle->Pt()); + Int_t id = (mother ? mother->PdgCode() : 2212); + Double_t weight = CalculateWeight(eta, pt, b, phi, rp, id); +#else + Double_t weight = 1.; +#endif + StoreParticle(particle, mother, iTrRef, vz, output, weight); } } fNRefs->Fill(nRef); } return kTRUE; } + +//____________________________________________________________________ +Double_t +AliSPDMCTrackDensity::CalculateWeight(Double_t eta, Double_t pt, Double_t b, + Double_t phi, Double_t rp, Int_t id) const +{ + static TF1 gaus = TF1("gaus", "gaus", -6, 6); + gaus.SetParameters(0.1, 0., 9); + // gaus.SetParameters(0.1, 0., 3); + // gaus.SetParameters(0.1, 0., 15); + + const Double_t xCumulant2nd4050ALICE[] = {0.00, 0.25, 0.350, + 0.45, 0.55, 0.650, + 0.75, 0.85, 0.950, + 1.10, 1.30, 1.500, + 1.70, 1.90, 2.250, + 2.75, 3.25, 3.750, + 4.50}; + const Double_t yCumulant2nd4050ALICE[] = {0.00000, 0.043400, + 0.059911,0.073516, + 0.089756,0.105486, + 0.117391,0.128199, + 0.138013,0.158271, + 0.177726,0.196383, + 0.208277,0.216648, + 0.242954,0.249961, + 0.240131,0.269006, + 0.207796}; + const Int_t nPointsCumulant2nd4050ALICE = + sizeof(xCumulant2nd4050ALICE)/sizeof(Double_t); + static TGraph alicePointsPt2(nPointsCumulant2nd4050ALICE,xCumulant2nd4050ALICE,yCumulant2nd4050ALICE); +#if 0 + const Double_t xCumulant4th3040ALICE[] = {0.00,0.250,0.35, + 0.45,0.550,0.65, + 0.75,0.850,0.95, + 1.10,1.300,1.50, + 1.70,1.900,2.25, + 2.75,3.250,3.75, + 4.50,5.500,7.00, + 9.000000}; + const Double_t yCumulant4th3040ALICE[] = {0.000000,0.037071, + 0.048566,0.061083, + 0.070910,0.078831, + 0.091396,0.102026, + 0.109691,0.124449, + 0.139819,0.155561, + 0.165701,0.173678, + 0.191149,0.202015, + 0.204540,0.212560, + 0.195885,0.000000, + 0.000000,0.000000}; +#endif + const Double_t xCumulant4th4050ALICE[] = {0.00,0.25,0.350, + 0.45,0.55,0.650, + 0.75,0.85,0.950, + 1.10,1.30,1.500, + 1.70,1.90,2.250, + 2.75,3.25,3.750, + 4.50}; + const Double_t yCumulant4th4050ALICE[] = {0.000000,0.038646, + 0.049824,0.066662, + 0.075856,0.081583, + 0.099778,0.104674, + 0.118545,0.131874, + 0.152959,0.155348, + 0.169751,0.179052, + 0.178532,0.198851, + 0.185737,0.239901, + 0.186098}; + const Int_t nPointsCumulant4th4050ALICE = + sizeof(xCumulant4th4050ALICE)/sizeof(Double_t); + static TGraph alicePointsPt4(nPointsCumulant4th4050ALICE, + xCumulant4th4050ALICE, + yCumulant4th4050ALICE); + + const Double_t xCumulant4thTPCrefMultTPConlyAll[] = {1.75, + 4.225, + 5.965, + 7.765, + 9.215, + 10.46, + 11.565, + 12.575}; + const Double_t yCumulant4thTPCrefMultTPConlyAll[] = {0.017855,0.032440, + 0.055818,0.073137, + 0.083898,0.086690, + 0.082040,0.077777}; + const Int_t nPointsCumulant4thTPCrefMultTPConlyAll = + sizeof(xCumulant4thTPCrefMultTPConlyAll)/sizeof(Double_t); + TGraph aliceCent(nPointsCumulant4thTPCrefMultTPConlyAll, + xCumulant4thTPCrefMultTPConlyAll, + yCumulant4thTPCrefMultTPConlyAll); + + + Double_t weight = (20. * gaus.Eval(eta) * (alicePointsPt2.Eval(pt) * 0.5 + + alicePointsPt4.Eval(pt) * 0.5) + * (aliceCent.Eval(b) / aliceCent.Eval(10.46)) + * 2. * TMath::Cos(2. * (phi - rp))); + if (TMath::Abs(id) == 211) weight *= 1.3; //pion flow + else if (TMath::Abs(id) == 2212) weight *= 1.0; //proton flow + else weight *= 0.7; + + return weight; +} + //____________________________________________________________________ void AliSPDMCTrackDensity::Print(Option_t* /*option*/) const diff --git a/PWGLF/FORWARD/analysis2/AliSPDMCTrackDensity.h b/PWGLF/FORWARD/analysis2/AliSPDMCTrackDensity.h index a0e5f64da4a..36af38438a5 100644 --- a/PWGLF/FORWARD/analysis2/AliSPDMCTrackDensity.h +++ b/PWGLF/FORWARD/analysis2/AliSPDMCTrackDensity.h @@ -117,7 +117,8 @@ protected: const AliMCParticle* mother, Int_t refNo, Double_t vz, - TH2D& output) const; + TH2D& output, + Double_t w) const; /** * Get ultimate mother of a track * @@ -137,6 +138,21 @@ protected: */ Double_t GetTrackRefTheta(const AliTrackReference* ref, Double_t vz) const; + /** + * Calculate flow weight + * + * @param eta Pseudo rapidity + * @param pt Transverse momemtum + * @param b Impact parameter + * @param phi Azimuthal angle + * @param rp Reaction plance angle + * @param id Particle PDG code + * + * @return Flow weight for the particle + */ + Double_t CalculateWeight(Double_t eta, Double_t pt, Double_t b, + Double_t phi, Double_t rp, Int_t id) const; + Bool_t fUseOnlyPrimary; // Only use primaries Double_t fMinR; // Min radius Double_t fMaxR; // Max radius diff --git a/PWGLF/FORWARD/analysis2/ForwardAODConfig.C b/PWGLF/FORWARD/analysis2/ForwardAODConfig.C index aed6bf972c2..6d2edd3431d 100644 --- a/PWGLF/FORWARD/analysis2/ForwardAODConfig.C +++ b/PWGLF/FORWARD/analysis2/ForwardAODConfig.C @@ -151,6 +151,9 @@ ForwardAODConfig(AliForwardMultiplicityBase* task) // Set the debug level of a single algorithm task->GetSharingFilter().SetDebug(0); + // --- Eventplane Finder ------------------------------------------- + task->GetEventPlaneFinder().SetUsePhiWeights(false); + // --- Set limits on fits the energy ------------------------------- // Maximum relative error on parameters AliFMDCorrELossFit::ELossFit::fgMaxRelError = .12; diff --git a/PWGLF/FORWARD/analysis2/MakeFlow.C b/PWGLF/FORWARD/analysis2/MakeFlow.C index ba7d7bc8bcf..424847aea84 100644 --- a/PWGLF/FORWARD/analysis2/MakeFlow.C +++ b/PWGLF/FORWARD/analysis2/MakeFlow.C @@ -5,7 +5,7 @@ * * @brief * - * @ingroup pwglf_forward_scripts_makers + * @ingroup pwg2_forward_scripts_makers * */ /** @@ -21,30 +21,35 @@ * @par Outputs: * - * - * @ingroup pwglf_forward_flow + * @ingroup pwg2_forward_flow */ void MakeFlow(TString data = "", - Int_t nevents = 0, + Int_t nEvents = 0, TString type = "", - Int_t etabins = 48, Bool_t mc = kFALSE, + const char* name = 0, + Int_t proof = 0, TString addFlow = "", Int_t addFType = 0, Int_t addFOrder = 0, - Bool_t proof = kFALSE) + Bool_t gdb = kFALSE) { - Bool_t proof = kFALSE; - // --- Load libs --------------------------------------------------- gROOT->Macro("$ALICE_ROOT/PWGLF/FORWARD/analysis2/scripts/LoadLibs.C"); - // --- Check for proof mode, and possibly upload pars -------------- - if (proof> 0) { - gROOT->LoadMacro("$ALICE_ROOT/PWGLF/FORWARD/analysis2/scripts/LoadPars.C"); - if (!LoadPars(proof)) { - AliError("MakeFlow", "Failed to load PARs"); - return; - } + // --- Possibly use plug-in for this ------------------------------- + if ((name && name[0] != '\0') && gSystem->Load("libRAliEn") >= 0) { + + gROOT->LoadMacro("$ALICE_ROOT/PWGLF/FORWARD/analysis2/trains/TrainSetup.C+"); + gROOT->LoadMacro("$ALICE_ROOT/PWGLF/FORWARD/analysis2/trains/MakeFlowTrain.C+"); + + MakeFlowTrain t(name, type.Data(), mc, addFlow.Data(), addFType, addFOrder, false); + t.SetDataDir(data.Data()); + t.SetDataSet(""); + t.SetProofServer(Form("workers=%d", proof)); + t.SetUseGDB(gdb); + t.Run(proof > 0 ? "proof" : "local", "full", nEvents, proof > 0); + return; } // --- Set the macro path ------------------------------------------ @@ -57,19 +62,10 @@ void MakeFlow(TString data = "", AliError("You didn't add a data file"); return; } - TChain* chain = new TChain("aodTree"); - - if (data.Contains(".txt")) MakeChain(data, chain); - - if (data.Contains(".root")) { - TFile* test = TFile::Open(data.Data()); - if (!test) { - AliError(Form("AOD file %s not found", data.Data())); - return; - } - test->Close(); // Remember to close! - chain->Add(data.Data()); - } + gROOT->LoadMacro("$ALICE_ROOT/PWGLF/FORWARD/analysis2/scripts/MakeChain.C"); + TChain* chain = MakeChain("AOD", data.Data(), true); + // If 0 or less events is select, choose all + if (nEvents <= 0) nEvents = chain->GetEntries(); // --- Initiate the event handlers -------------------------------- AliAnalysisManager *mgr = new AliAnalysisManager("Forward Flow", @@ -81,7 +77,7 @@ void MakeFlow(TString data = "", // --- Add the tasks --------------------------------------------- gROOT->LoadMacro("AddTaskForwardFlow.C"); - AddTaskForwardFlow(type, etabins, mc, addFlow, addFType, addFOrder); + AddTaskForwardFlow(type, mc, addFlow, addFType, addFOrder); // --- Run the analysis -------------------------------------------- TStopwatch t; @@ -91,51 +87,21 @@ void MakeFlow(TString data = "", } mgr->PrintStatus(); Printf("****************************************"); - Printf("Doing flow analysis on %d Events", nevents == 0 ? chain->GetEntries() : nevents); + Printf("Doing flow analysis on %d Events", nEvents); Printf("****************************************"); // - if (proof) mgr->SetDebugLevel(3); - if (mgr->GetDebugLevel() < 1 && !proof) - mgr->SetUseProgressBar(kTRUE,chain->GetEntries() < 10000 ? 100 : 1000); + mgr->SetDebugLevel(0); + if (mgr->GetDebugLevel() < 1) + mgr->SetUseProgressBar(kTRUE, nEvents < 10000 ? 100 : 1000); // mgr->SetSkipTerminate(true); t.Start(); - if (nevents == 0) mgr->StartAnalysis("local", chain); - if (nevents != 0) mgr->StartAnalysis("local", chain, nevents); + mgr->StartAnalysis("local", chain, nEvents); t.Stop(); t.Print(); } //---------------------------------------------------------------- -void MakeChain(TString data = "", TChain* chain = 0) -{ - // creates chain of files in a given directory or file containing a list. - - // Open the input stream - ifstream in; - in.open(data.Data()); - - // Read the input list of files and add them to the chain - TString line; - TFile* file; - while(in.good()) - { - in >> line; - - if (line.Length() == 0) - continue; - if (!(file = TFile::Open(line))) - gROOT->ProcessLine(Form(".!rm %s", line.Data())); - else { - chain->Add(line); - file->Close(); - } - } - - in.close(); - - return; -} // // EOF // diff --git a/PWGLF/FORWARD/analysis2/trains/MakeAODTrain.C b/PWGLF/FORWARD/analysis2/trains/MakeAODTrain.C index 3ac1c707a75..9f0c7d035f1 100644 --- a/PWGLF/FORWARD/analysis2/trains/MakeAODTrain.C +++ b/PWGLF/FORWARD/analysis2/trains/MakeAODTrain.C @@ -1,3 +1,5 @@ +#include "TrainSetup.C" + /** * @file MakeAODTrain.C * @author Christian Holm Christensen @@ -8,7 +10,6 @@ * @ingroup pwglf_forward_trains */ //==================================================================== -#include "TrainSetup.C" /** * Analysis train to make Forward and Central multiplicity @@ -117,10 +118,15 @@ protected: // --- Set load path --------------------------------------------- gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWGLF/FORWARD/analysis2", gROOT->GetMacroPath())); + gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/ANALYSIS/macros", + gROOT->GetMacroPath())); // --- Check if this is MC --------------------------------------- Bool_t mc = mgr->GetMCtruthEventHandler() != 0; + // --- Add TPC eventplane task + gROOT->Macro("AddTaskEventplane.C"); + // --- Task to copy header information --------------------------- gROOT->Macro("AddTaskCopyHeader.C"); @@ -133,6 +139,10 @@ protected: gROOT->Macro(Form("AddTaskCentralMult.C(%d,%d,%d,%d)", mc, fSys, fSNN, fField)); AddExtraFile(gSystem->Which(gROOT->GetMacroPath(), "CentralAODConfig.C")); + + // --- Add MC particle task -------------------------------------- + if (mc) gROOT->Macro("AddTaskMCParticleFilter.C"); + } //__________________________________________________________________ /** @@ -161,6 +171,41 @@ protected: // --- Ignore trigger class when selecting events. This means --- // --- that we get offline+(A,C,E) events too -------------------- // ps->SetSkipTriggerClassSelection(true); +/* + if (mc) { + AliOADBPhysicsSelection * oadbDefaultPbPb = new AliOADBPhysicsSelection("oadbDefaultPbPb"); + oadbDefaultPbPb->AddCollisionTriggerClass ( AliVEvent::kHighMult,"+C0SMH-B-NOPF-ALLNOTRD","B",0); + oadbDefaultPbPb->AddBGTriggerClass ( AliVEvent::kHighMult,"+C0SMH-A-NOPF-ALLNOTRD","A",0); + oadbDefaultPbPb->AddBGTriggerClass ( AliVEvent::kHighMult,"+C0SMH-C-NOPF-ALLNOTRD","C",0); + oadbDefaultPbPb->AddBGTriggerClass ( AliVEvent::kHighMult,"+C0SMH-E-NOPF-ALLNOTRD","E",0); + oadbDefaultPbPb->SetHardwareTrigger ( 0,"SPDGFO >= 100"); + oadbDefaultPbPb->SetOfflineTrigger ( 0,"SPDGFO >= 100 && !V0ABG && !V0CBG"); + + oadbDefaultPbPb->AddCollisionTriggerClass ( AliVEvent::kMB,"+CMBACS2-B-NOPF-ALLNOTRD","B",1); + oadbDefaultPbPb->AddBGTriggerClass ( AliVEvent::kMB,"+CMBACS2-A-NOPF-ALLNOTRD","A",1); + oadbDefaultPbPb->AddBGTriggerClass ( AliVEvent::kMB,"+CMBACS2-C-NOPF-ALLNOTRD","C",1); + oadbDefaultPbPb->AddBGTriggerClass ( AliVEvent::kMB,"+CMBACS2-E-NOPF-ALLNOTRD","E",1); + oadbDefaultPbPb->SetHardwareTrigger ( 1,"(V0A && V0C && SPDGFO > 1)"); + oadbDefaultPbPb->SetOfflineTrigger ( 1,"(V0A && V0C && SPDGFO > 1) && !V0ABG && !V0CBG"); + // LHC10h8 + AliOADBPhysicsSelection * oadbLHC10h8 = new AliOADBPhysicsSelection("oadbLHC10h8"); + oadbLHC10h8->AddCollisionTriggerClass ( AliVEvent::kHighMult,"+C0SMH-B-NOPF-ALL","B",0); + oadbLHC10h8->AddBGTriggerClass ( AliVEvent::kHighMult,"+C0SMH-A-NOPF-ALL","A",0); + oadbLHC10h8->AddBGTriggerClass ( AliVEvent::kHighMult,"+C0SMH-C-NOPF-ALL","C",0); + oadbLHC10h8->AddBGTriggerClass ( AliVEvent::kHighMult,"+C0SMH-E-NOPF-ALL","E",0); + oadbLHC10h8->SetHardwareTrigger ( 0,"SPDGFO >= 100"); + oadbLHC10h8->SetOfflineTrigger ( 0,"SPDGFO >= 100 && !V0ABG && !V0CBG"); + + oadbLHC10h8->AddCollisionTriggerClass ( AliVEvent::kMB,"+CMBACS2-B-NOPF-ALL","B",1); + oadbLHC10h8->AddBGTriggerClass ( AliVEvent::kMB,"+CMBACS2-A-NOPF-ALL","A",1); + oadbLHC10h8->AddBGTriggerClass ( AliVEvent::kMB,"+CMBACS2-C-NOPF-ALL","C",1); + oadbLHC10h8->AddBGTriggerClass ( AliVEvent::kMB,"+CMBACS2-E-NOPF-ALL","E",1); + oadbLHC10h8->SetHardwareTrigger ( 1,"(V0A && V0C && SPDGFOL1 > 1)"); + oadbLHC10h8->SetOfflineTrigger ( 1,"(V0A && V0C && SPDGFOL1 > 1) && !V0ABG && !V0CBG"); + + ps->SetCustomOADBObjects(oadbDefaultPbPb, 0, 0); + ps->SetCustomOADBObjects(oadbLHC10h8, 0, 0); + }*/ } //__________________________________________________________________ /** diff --git a/PWGLF/FORWARD/analysis2/trains/MakeFlowTrain.C b/PWGLF/FORWARD/analysis2/trains/MakeFlowTrain.C new file mode 100644 index 00000000000..2f921427972 --- /dev/null +++ b/PWGLF/FORWARD/analysis2/trains/MakeFlowTrain.C @@ -0,0 +1,124 @@ +#include "TrainSetup.C" + +//==================================================================== +/** + * Analysis train to make @f$ flow@f$ + * + * To run, do + * @code + * gROOT->LoadMacro("TrainSetup.C"); + * // Make train + * MakeFlowTrain t("My Analysis"); + * // Set variaous parameters on the train + * t.SetDataDir("/home/of/data"); + * t.AddRun(118506) + * // Run it + * t.Run("LOCAL", "FULL", -1, false, false); + * @endcode + * + * @ingroup pwg2_forward_flow + * @ingroup pwg2_forward_trains + */ +class MakeFlowTrain : public TrainSetup +{ +public: + /** + * Constructor. Date and time must be specified when running this + * in Termiante mode on Grid + * + * @param name Name of train (free form) + * @param dateTime Append date and time to name + * @param year Year - if not specified, current year + * @param month Month - if not specified, current month + * @param day Day - if not specified, current day + * @param hour Hour - if not specified, current hour + * @param min Minutes - if not specified, current minutes + */ + MakeFlowTrain(const char* name, + char* type="", + Bool_t mc = false, + char* addFlow = "", + Int_t addFType = 0, + Int_t addFOrder = 0, + Bool_t dateTime=false, + UShort_t year = 0, + UShort_t month = 0, + UShort_t day = 0, + UShort_t hour = 0, + UShort_t min = 0) + : TrainSetup(name, dateTime, year, month, day, hour, min), + fType(type), + fMC(mc), + fAddFlow(addFlow), + fAddFType(addFType), + fAddFOrder(addFOrder) + { + } + /** + * Run this analysis + * + * @param mode Mode - see TrainSetup::EMode + * @param oper Operation - see TrainSetup::EOperation + * @param nEvents Number of events (negative means all) + * @param usePar If true, use PARs + */ + void Run(const char* mode, const char* oper, + Int_t nEvents=-1, Bool_t usePar=false) + { + Exec("AOD", mode, oper, nEvents, false, usePar); + } + /** + * Run this analysis + * + * @param mode Mode - see TrainSetup::EMode + * @param oper Operation - see TrainSetup::EOperation + * @param nEvents Number of events (negative means all) + * @param usePar If true, use PARs + */ + void Run(EMode mode, EOper oper, Int_t nEvents=-1, + Bool_t usePar=false) + { + Exec(kAOD, mode, oper, nEvents, false, usePar); + } +protected: + /** + * Create the tasks + * + * @param mode Processing mode + * @param par Whether to use par files + */ + void CreateTasks(EMode mode, Bool_t par, AliAnalysisManager*) + { + // --- Output file name ------------------------------------------ + AliAnalysisManager::SetCommonFileName("AnalysisResults.root"); + + // --- Load libraries/pars --------------------------------------- + LoadLibrary("PWGLFforward2", mode, par, true); + + // --- Set load path --------------------------------------------- + gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWGLF/FORWARD/analysis2", + gROOT->GetMacroPath())); + + // --- Add the task ---------------------------------------------- + gROOT->Macro(Form("AddTaskForwardFlow.C(\"%s\",%d,\"%s\",%d,%d)", + fType, fMC, fAddFlow, fAddFType, fAddFOrder)); + } + /** + * Do not the centrality selection + */ + void CreateCentralitySelection(Bool_t, AliAnalysisManager*) {} + /** + * Crete output handler - we don't want one here. + * + * @return 0 + */ + AliVEventHandler* CreateOutputHandler(EType) { return 0; } + char* fType; + Bool_t fMC; + char* fAddFlow; + Int_t fAddFType; + Int_t fAddFOrder; +}; +// +// EOF +// diff --git a/PWGLF/PWGLFforward2LinkDef.h b/PWGLF/PWGLFforward2LinkDef.h index 99514c092b6..9e30e89fa3e 100644 --- a/PWGLF/PWGLFforward2LinkDef.h +++ b/PWGLF/PWGLFforward2LinkDef.h @@ -34,6 +34,7 @@ #pragma link C++ class AliForwardUtil::RingHistos+; #pragma link C++ class AliFMDEventInspector+; #pragma link C++ class AliFMDMCEventInspector+; +#pragma link C++ class AliFMDEventPlaneFinder+; #pragma link C++ class AliFMDSharingFilter+; #pragma link C++ class AliFMDSharingFilter::RingHistos+; #pragma link C++ class AliFMDMCSharingFilter+; @@ -57,7 +58,8 @@ #pragma link C++ class AliFMDCorrDoubleHit+; #pragma link C++ class AliFMDCorrVertexBias+; #pragma link C++ class AliFMDCorrMergingEfficiency+; -#pragma link C++ class AliAODForwardMult+; +#pragma link C++ class AliAODForwardMult+; +#pragma link C++ class AliAODForwardEP+; #pragma link C++ class AliForwardMultiplicityBase+; #pragma link C++ class AliForwardMultiplicityTask+; #pragma link C++ class AliForwardMCMultiplicityTask+; @@ -81,8 +83,9 @@ #pragma link C++ class AliCentralCorrSecondaryMap+; #pragma link C++ class AliCentralCorrAcceptance+; #pragma link C++ class AliCentraldNdetaTask+; -#pragma link C++ class AliForwardFlowUtil+; #pragma link C++ class AliForwardFlowTaskQC+; +#pragma link C++ class AliForwardFlowTaskQC::VertexBin+; +#pragma link C++ class AliForwardMCFlowTaskQC+; #pragma link C++ class AliSPDMCTrackDensity+; #pragma link C++ class AliFMDMultCuts+; #pragma link C++ class AliPoissonCalculator+; -- 2.43.0