From 58f5fae2f4dbad0a7830d1616b20d5646fa0d08e Mon Sep 17 00:00:00 2001 From: hdalsgaa Date: Mon, 14 Mar 2011 10:05:52 +0000 Subject: [PATCH] Upgrades and fixes for coding violations from FC --- PWG2/CMakelibPWG2forward2.pkg | 6 +- PWG2/FORWARD/analysis2/AddTaskForwardFlow.C | 13 +- .../analysis2/AliForwardFlowTaskQC.cxx | 1166 +++++++++-------- PWG2/FORWARD/analysis2/AliForwardFlowTaskQC.h | 34 +- PWG2/FORWARD/analysis2/AliForwardFlowUtil.cxx | 363 +++++ PWG2/FORWARD/analysis2/AliForwardFlowUtil.h | 102 ++ PWG2/FORWARD/analysis2/MakeFlow.C | 81 +- PWG2/PWG2forward2LinkDef.h | 2 +- 8 files changed, 1219 insertions(+), 548 deletions(-) create mode 100644 PWG2/FORWARD/analysis2/AliForwardFlowUtil.cxx create mode 100644 PWG2/FORWARD/analysis2/AliForwardFlowUtil.h diff --git a/PWG2/CMakelibPWG2forward2.pkg b/PWG2/CMakelibPWG2forward2.pkg index 5023a72f3f4..22d7cf7a749 100644 --- a/PWG2/CMakelibPWG2forward2.pkg +++ b/PWG2/CMakelibPWG2forward2.pkg @@ -57,18 +57,18 @@ set ( SRCS FORWARD/analysis2/AliCentralCorrAcceptance.cxx FORWARD/analysis2/AliCentraldNdetaTask.cxx FORWARD/analysis2/AliBasedNdetaTask.cxx - FORWARD/analysis2/AliForwardFlowBase.cxx + FORWARD/analysis2/AliForwardFlowUtil.cxx FORWARD/analysis2/AliForwardFlowTaskQC.cxx ) string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" ) set ( HDRS ${HDRS} FORWARD/analysis2/AliFMDStripIndex.h ) -set ( EINCLUDE ANALYSIS PWG2/FORWARD/analysis2) +set ( EINCLUDE ANALYSIS PWG2/FORWARD/analysis2 STEER) set ( EXPORT FORWARD/analysis2/AliAODForwardMult.h FORWARD/analysis2/AliAODCentralMult.h - FORWARD/analysis2/AliForwardUtil.h ) + FORWARD/analysis2/AliForwardUtil.h ) set ( DHDR PWG2forward2LinkDef.h) diff --git a/PWG2/FORWARD/analysis2/AddTaskForwardFlow.C b/PWG2/FORWARD/analysis2/AddTaskForwardFlow.C index 227e6e005bf..fd3df1c2c14 100644 --- a/PWG2/FORWARD/analysis2/AddTaskForwardFlow.C +++ b/PWG2/FORWARD/analysis2/AddTaskForwardFlow.C @@ -1,4 +1,9 @@ -void AddTaskForwardFlow(TString type = "", Int_t etabins = 20) +void AddTaskForwardFlow(TString type = "", + Int_t etabins = 40, + Int_t zVertex = 2, + TString addFlow = "", + Int_t addFType = 0, + Int_t addFOrder = 0) { AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); if (!mgr) { @@ -42,8 +47,12 @@ void AddTaskForwardFlow(TString type = "", Int_t etabins = 20) qc->SetDoHarmonics(v1, v2, v3, v4); qc->SetUseNEtaBins(etabins); + qc->AddFlow(addFlow); + qc->AddFlowType(addFType); + qc->AddFlowOrder(addFOrder); + qc->SetVertexRange(zVertex); - mgr->ConnectInput(qc ,0, mgr->GetCommonInputContainer()); + mgr->ConnectInput(qc, 0, mgr->GetCommonInputContainer()); mgr->ConnectOutput(qc, 0, mgr->GetCommonOutputContainer()); mgr->ConnectOutput(qc, 1, qcout); diff --git a/PWG2/FORWARD/analysis2/AliForwardFlowTaskQC.cxx b/PWG2/FORWARD/analysis2/AliForwardFlowTaskQC.cxx index 770a529fb34..c6f617e16bd 100644 --- a/PWG2/FORWARD/analysis2/AliForwardFlowTaskQC.cxx +++ b/PWG2/FORWARD/analysis2/AliForwardFlowTaskQC.cxx @@ -7,10 +7,6 @@ // Outputs: // - AnalysisResults.root // -// TODO! -// - Add centrality stuff -// END OF TODO! -// #include #include #include @@ -19,14 +15,27 @@ #include #include #include +#include "TH3D.h" #include "AliLog.h" #include "AliForwardFlowTaskQC.h" #include "AliAnalysisManager.h" #include "AliAODHandler.h" #include "AliAODInputHandler.h" #include "AliAODMCParticle.h" -#include "AliForwardFlowBase.h" #include "AliAODForwardMult.h" +#include "AliAODEvent.h" + +// +// Enumeration for adding and retrieving stuff from the histogram +// +enum { kW2avg2 = 1, kW2, kW2avg2sq, kW2w2pavg2avg2p, kW2w2p, kW4avg4, + kW4, kW4avg4sq, kW2w4, kW2w4avg2avg4, kW2w4p, kW2w4pavg2avg4p, + kW4w2p, kW4w2pavg4avg2p, kW4w4p, kW4w4pavg4avg4p, kQnRe, kQnIm, + kM, kCosphi1phi2, kSinphi1phi2, kCosphi1phi2phi3m, kSinphi1phi2phi3m, + kMm1m2, kCosphi1phi2phi3p, kSinphi1phi2phi3p, + kRW2avg2, kRW2, kRW2avg2sq, kRM, kRW4avg4, kRW4, kRW4avg4sq, + kRW2w4, kRW2w4avg2avg4, kRCosphi1phi2, kRSinphi1phi2, + kRCosphi1phi2phi3m, kRSinphi1phi2phi3m, kRMm1m2 }; ClassImp(AliForwardFlowTaskQC) #if 0 @@ -34,11 +43,16 @@ ClassImp(AliForwardFlowTaskQC) #endif AliForwardFlowTaskQC::AliForwardFlowTaskQC() - : fDebug(0), // Debug flag - fOutputList(0), // Output list + : fOutputList(0), // Output list + fFlowUtil(0), // AliForwardFlowUtil fAOD(0), // AOD input event - fMC(kFALSE), // MC flag - fEtaBins(20) // # of eta bins in histograms + fMC(kFALSE), // MC flag + fEtaBins(20), // # of etaBin bins in histograms + fAddFlow(0), // Add flow string + fAddType(0), // Add flow type # + fAddOrder(0), // Add flow order + fZvertex(0), // Z vertex range + fCent(0) // Centrality { // // Default constructor @@ -47,11 +61,16 @@ AliForwardFlowTaskQC::AliForwardFlowTaskQC() //_____________________________________________________________________ AliForwardFlowTaskQC::AliForwardFlowTaskQC(const char* name) : AliAnalysisTaskSE(name), - fDebug(0), // Debug flag fOutputList(0), // Output list + fFlowUtil(0), // AliForwardFlowUtil fAOD(0), // AOD input event fMC(kFALSE), // MC flag - fEtaBins(20) // # of Eta bins + fEtaBins(20), // # of Eta bins + fAddFlow(0), // Add flow string + fAddType(0), // Add flow type # + fAddOrder(0), // Add flow order + fZvertex(0), // Z vertex range + fCent(0) // Centrality { // // Constructor @@ -65,11 +84,16 @@ AliForwardFlowTaskQC::AliForwardFlowTaskQC(const char* name) : //_____________________________________________________________________ AliForwardFlowTaskQC::AliForwardFlowTaskQC(const AliForwardFlowTaskQC& o) : AliAnalysisTaskSE(o), - fDebug(o.fDebug), // Debug flag 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 + fEtaBins(o.fEtaBins), // # of Eta bins + 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(0) // Centrality { // // Copy constructor @@ -89,114 +113,126 @@ void AliForwardFlowTaskQC::CreateOutputObjects() if (!fOutputList) fOutputList = new TList(); fOutputList->SetName("QCumulants"); + fOutputList->SetOwner(); + fFlowUtil = new AliForwardFlowUtil(fOutputList); + fFlowUtil->SetVertexRange(fZvertex); if (fEtaBins % 20) fEtaBins = 20; // Histograms for cumulants analysis // We loop over flow histograms here to add different orders of harmonics - for (Int_t n = 1; n <= 4; n++) { - if (!fv[n]) continue; - // Only one flow histogram is needed for each type of data; - // x-axis is eta-bins with differential flow, integrated is in underflowbin - // y-axis bin 1: (w_<2> * <2>).Re() - // y-axis bin 2: (w_<2> * <2>).Im() - // y-axis bin 3: w_<2> = M(M-1) - // y-axis bin 4: (w_<2> * <2> * <2>).Re() - // y-axis bin 5: (w_<2> * <2> * <2>).Im() - // y-axis bin 6: (w_<2> * w_<2'> * <2> * <2'>).Re() - // y-axis bin 7: (w_<2> * w_<2'> * <2> * <2'>).Im() - // y-axis bin 8: w_<2> * w_<2'> - // y-axis bin 9: (w_<4> * <4>).Re() - // y-axis bin 10: (w_<4> * <4>).Im() - // y-axis bin 11: w_<4> - // y-axis bin 12: (w_<4> * <4> * <4>).Re() - // y-axis bin 13: (w_<4> * <4> * <4>).Im() - // y-axis bin 14: w_<2> * w_<4> - // y-axis bin 15: (w_<2> * w_<4> * <2> * <4>).Re() - // y-axis bin 16: (w_<2> * w_<4> * <2> * <4>).Im() - // y-axis bin 17: w_<2> * w_<4'> - // y-axis bin 18: (w_<2> * w_<4'> * <2> * <4'>).Re() - // y-axis bin 19: (w_<2> * w_<4'> * <2> * <4'>).Im() - // y-axis bin 20: w_<4> * w_<2'> - // y-axis bin 21: (w_<4> * w_<2'> * <4> * <2'>).Re() - // y-axis bin 22: (w_<4> * w_<2'> * <4> * <2'>).Im() - // y-axis bin 23: w_<4> * w_<4'> - // y-axis bin 24: (w_<4> * w_<4'> * <4> * <4'>).Re() - // y-axis bin 25: (w_<4> * w_<4'> * <4> * <4'>).Im() - // y-axis bin 26: Qn or pn.Re() = <> - // y-axis bin 27: Qn or pn.Im() = <> - // y-axis bin 28: M or mp - // y-axis bin 29: (Qn*Qn-Q2n).Re() = <> - // y-axis bin 30: (Qn*Qn-Q2n).Im() = <> - // y-axis bin 31: <> - // y-axis bin 32: <> - // y-axis bin 33: M*(M-1)*(M-2) or similar for diff - // y-axis bin 34: <> - // y-axis bin 35: <> - TH2D* hFlowHist = new TH2D(Form("hQ%dCumuHist", n), Form("hQ%dCumuHist", n), fEtaBins, -4, 6, 35, 0.5, 35.5); - hFlowHist->Sumw2(); - fOutputList->Add(hFlowHist); - - TH2D* hFlowHistMC = new TH2D(Form("hQ%dCumuHistMC", n), Form("hQ%dCumuHistMC", n), fEtaBins, -4, 6, 35, 0.5, 35.5); - hFlowHistMC->Sumw2(); - fOutputList->Add(hFlowHistMC); - - TH2D* hFlowHistTrRef = new TH2D(Form("hQ%dCumuHistTrRef", n), Form("hQ%dCumuHistTrRef", n), fEtaBins, -4, 6, 35, 0.5, 35.5); - hFlowHistTrRef->Sumw2(); - fOutputList->Add(hFlowHistTrRef); - - // Output histograms - TH1D* hCumulant2Flow = new TH1D(Form("hQ%dCumulant2Flow", n), Form("hQ%dCumulant2Flow", n), fEtaBins, -4, 6); - hCumulant2Flow->Sumw2(); - fOutputList->Add(hCumulant2Flow); - - TH1D* hCumulant2FlowMC = new TH1D(Form("hQ%dCumulant2FlowMC", n),Form("hQ%dCumulant2FlowMC", n), fEtaBins, -4, 6); - hCumulant2FlowMC->Sumw2(); - fOutputList->Add(hCumulant2FlowMC); - - TH1D* hCumulant2FlowTrRef = new TH1D(Form("hQ%dCumulant2FlowTrRef", n), Form("hQ%dCumulant2FlowTrRef", n), fEtaBins, -4, 6); - hCumulant2FlowTrRef->Sumw2(); - fOutputList->Add(hCumulant2FlowTrRef); - - - TH1D* hCumulant4Flow = new TH1D(Form("hQ%dCumulant4Flow", n), Form("hQ%dCumulant4Flow", n), fEtaBins, -4, 6); - hCumulant4Flow->Sumw2(); - fOutputList->Add(hCumulant4Flow); - - TH1D* hCumulant4FlowMC = new TH1D(Form("hQ%dCumulant4FlowMC", n), Form("hQ%dCumulant4FlowMC", n), fEtaBins, -4, 6); - hCumulant4FlowMC->Sumw2(); - fOutputList->Add(hCumulant4FlowMC); - - TH1D* hCumulant4FlowTrRef = new TH1D(Form("hQ%dCumulant4FlowTrRef", n), Form("hQ%dCumulant4FlowTrRef", n), fEtaBins, -4, 6); - hCumulant4FlowTrRef->Sumw2(); - fOutputList->Add(hCumulant4FlowTrRef); - } - - // Single Event histograms - TH1D* hdNdphiSE = new TH1D("hdNdphiSE","hdNdphiSE", 20, 0, 2*TMath::Pi()); - hdNdphiSE->Sumw2(); - fOutputList->Add(hdNdphiSE); - - TH2D* hdNdetadphiSE = new TH2D("hdNdetadphiSE", "hdNdetadphiSE", fEtaBins, -4, 6, 20, 0, 2*TMath::Pi()); - hdNdetadphiSE->Sumw2(); - fOutputList->Add(hdNdetadphiSE); - - TH1D* hdNdphiSEMC = new TH1D("hdNdphiSEMC","hdNdphiSEMC", 20, 0, 2*TMath::Pi()); - hdNdphiSEMC->Sumw2(); - fOutputList->Add(hdNdphiSEMC); + TString type = ""; + for (Int_t loop = 1; loop <= 4; loop++) { + + if (loop == 1) type = ""; + if (loop == 2) type = "SPD"; + if (loop == 3) type = "MC"; + if (loop == 4) type = "TrRef"; - TH2D* hdNdetadphiSEMC = new TH2D("hdNdetadphiSEMC", "hdNdetadphiSEMC", fEtaBins, -4, 6, 20, 0, 2*TMath::Pi()); - hdNdetadphiSEMC->Sumw2(); - fOutputList->Add(hdNdetadphiSEMC); + for (Int_t n = 1; n <= 4; n++) { + if (!fv[n]) continue; + // Only one flow histogram is needed for each type of data; + // x-axis is etaBin-bins with differential flow, integrated is in underflowbin + // z-axis bin 1: (w_<2> * <2>).Re() + // z-axis bin 2: w_<2> = mqM-mp + // z-axis bin 3: (w_<2> * <2> * <2>).Re() + // z-axis bin 4: (w_<2> * w_<2'> * <2> * <2'>).Re() + // z-axis bin 5: w_<2> * w_<2'> + // z-axis bin 6: (w_<4> * <4>).Re() + // z-axis bin 7: w_<4> + // z-axis bin 8: (w_<4> * <4> * <4>).Re() + // z-axis bin 9: w_<2> * w_<4> + // z-axis bin 10: (w_<2> * w_<4> * <2> * <4>).Re() + // z-axis bin 11: w_<2> * w_<4'> + // z-axis bin 12: (w_<2> * w_<4'> * <2> * <4'>).Re() + // z-axis bin 13: w_<4> * w_<2'> + // z-axis bin 14: (w_<4> * w_<2'> * <4> * <2'>).Re() + // z-axis bin 15: w_<4> * w_<4'> + // z-axis bin 16: (w_<4> * w_<4'> * <4> * <4'>).Re() + // z-axis bin 17: Qn or pn.Re() = <> + // z-axis bin 18: Qn or pn.Im() = <> + // z-axis bin 19: M or mp + // z-axis bin 20: (Qn*Qn-Q2n).Re() = <> + // z-axis bin 21: (Qn*Qn-Q2n).Im() = <> + // z-axis bin 22: <> + // z-axis bin 23: <> + // z-axis bin 24: M*(M-1)*(M-2) or similar for diff + // z-axis bin 25: <> + // z-axis bin 26: <> + // z-axis bin 27: ref: W_<2> * <2> + // z-axis bin 28: ref: W_<2> + // z-axis bin 29: ref: W_<2> * <2> * <2> + // z-axis bin 30: ref: Mult + // z-axis bin 31: ref: W_<4> * <4> + // z-axis bin 32: ref: W_<4> + // z-axis bin 33: ref: W_<4> * <4> * <4> + // z-axis bin 34: ref: W_<2> * W_<4> + // z-axis bin 35: ref: W_<2> * W_<4> * <2> * <4> + // z-axis bin 36: ref: + // z-axis bin 37: ref: + // z-axis bin 38: ref: + // z-axis bin 39: ref: + // z-axis bin 40: ref: M(M-1)(M-2) + Double_t x[41] = { 0. }; + for (Int_t e = 0; e <=fEtaBins; e++) { + x[e] = -4. + e*(10./(Double_t)fEtaBins); + } +// Double_t x[6] = { 0.0, 1.0, 2.0, 3.0, 4.5, 6.0 }; + Double_t y[11] = { 0., 5., 10., 20., 30., 40., 50., 60., 70., 80., 100. }; + Double_t z[41] = { 0. }; + for (Int_t k = 0; k <= 40; k++) { + z[k] = 0.5 + k*1.; + } + TH3D* hFlowHist = new TH3D(Form("hQ%dCumuHist%s", n, type.Data()), + Form("hQ%dCumuHist%s", n, type.Data()), fEtaBins, x, 10, y, 40, z); +// hFlowHist->RebinAxis(40/fEtaBins, hFlowHist->GetXaxis()); + hFlowHist->Sumw2(); + fOutputList->Add(hFlowHist); + TString tag = TString(); + for (Int_t c = -2; c < 10; c++) { + // Output histograms + if (c == -2) tag = "mb"; + else if (c == -1) tag = "0_40"; + else tag = Form("%d_%d", (Int_t)y[c], (Int_t)y[c+1]); + + TH1D* hCumulant2RefFlow = new TH1D(Form("hQ%dCumulant2RefFlow%s_%s", n, type.Data(), tag.Data()), + Form("hQ%dCumulant2RefFlow%s_%s", n, type.Data(), tag.Data()), fEtaBins, x); + hCumulant2RefFlow->Sumw2(); + fOutputList->Add(hCumulant2RefFlow); + + TH1D* hCumulant4RefFlow = new TH1D(Form("hQ%dCumulant4RefFlow%s_%s", n, type.Data(), tag.Data()), + Form("hQ%dCumulant4RefFlow%s_%s", n, type.Data(), tag.Data()), fEtaBins, x); + hCumulant4RefFlow->Sumw2(); + fOutputList->Add(hCumulant4RefFlow); + + TH1D* hCumulant2DiffFlow = new TH1D(Form("hQ%dCumulant2DiffFlow%s_%s", n, type.Data(), tag.Data()), + Form("hQ%dCumulant2DiffFlow%s_%s", n, type.Data(), tag.Data()), fEtaBins, x); + hCumulant2DiffFlow->Sumw2(); + fOutputList->Add(hCumulant2DiffFlow); + + TH1D* hCumulant4DiffFlow = new TH1D(Form("hQ%dCumulant4DiffFlow%s_%s", n, type.Data(), tag.Data()), + Form("hQ%dCumulant4DiffFlow%s_%s", n, type.Data(), tag.Data()), fEtaBins, x); + hCumulant4DiffFlow->Sumw2(); + fOutputList->Add(hCumulant4DiffFlow); + } // end of centrality loop + + } // end of v_{n} loop + + // Single Event histograms + TH1D* hdNdphiSE = new TH1D(Form("hdNdphiSE%s", type.Data()), + Form("hdNdphiSE%s", type.Data()), 20, 0, 2*TMath::Pi()); + hdNdphiSE->Sumw2(); + fOutputList->Add(hdNdphiSE); - TH1D* hdNdphiSETrRef = new TH1D("hdNdphiSETrRef","hdNdphiSETrRef", 20, 0, 2*TMath::Pi()); - hdNdphiSETrRef->Sumw2(); - fOutputList->Add(hdNdphiSETrRef); + TH2D* hdNdetaBindphiSE = new TH2D(Form("hdNdetaBindphiSE%s", type.Data()), + Form("hdNdetaBindphiSE%s", type.Data()), fEtaBins, -4, 6, 20, 0, 2*TMath::Pi()); + hdNdetaBindphiSE->Sumw2(); + fOutputList->Add(hdNdetaBindphiSE); + } // end of type loop - TH2D* hdNdetadphiSETrRef = new TH2D("hdNdetadphiSETrRef", "hdNdetadphiSETrRef", fEtaBins, -4, 6, 20, 0, 2*TMath::Pi()); - hdNdetadphiSETrRef->Sumw2(); - fOutputList->Add(hdNdetadphiSETrRef); + TH1D* cent = new TH1D("cent", "cent", 100, 0, 100); + fOutputList->Add(cent); PostData(1, fOutputList); } @@ -214,26 +250,16 @@ void AliForwardFlowTaskQC::UserExec(Option_t */*option*/) fAOD = dynamic_cast(InputEvent()); if (!fAOD) return; - // Load histograms and reset from last event - TH1D* dNdphi = (TH1D*)fOutputList->FindObject("hdNdphiSE"); - TH2D* dNdetadphi = (TH2D*)fOutputList->FindObject("hdNdetadphiSE"); - - dNdphi->Reset(); - dNdetadphi->Reset(); - // Initiate FlowCommon and fill histograms - AliForwardFlowBase* common = new AliForwardFlowBase(fOutputList); + // fill histograms + if (!fFlowUtil->LoopAODFMD(fAOD)) return; - if (!common->LoopAODFMD(fAOD)) return; - // else if (!common->LoopAODSPD(fAOD)) return; - // if (!common->LoopAODFMDandSPD(fAOD)) return; - - // Run analysis + // Run analysis on FMD for (Int_t n = 1; n <= 4; n++) { if (fv[n]) CumulantsMethod("", n); } - + // Find out if there's any MC data present if (!fMC) { TClonesArray* mcArray = 0; @@ -243,6 +269,12 @@ void AliForwardFlowTaskQC::UserExec(Option_t */*option*/) if (fMC) ProcessPrimary(); + if (!fFlowUtil->LoopAODSPD(fAOD)) return; + for (Int_t n = 1; n <= 4; n++) { + if (fv[n]) + CumulantsMethod("SPD", n); + } + } //_____________________________________________________________________ void AliForwardFlowTaskQC::CumulantsMethod(TString type = "", @@ -261,189 +293,252 @@ void AliForwardFlowTaskQC::CumulantsMethod(TString type = "", Double_t n = harmonic; // We get histograms depending on if it's real data or MC truth data - TH2D* flowHist = (TH2D*)fOutputList->FindObject(Form("hQ%dCumuHist%s", harmonic, type.Data())); - TH1D* dNdphi = (TH1D*)fOutputList->FindObject(Form("hdNdphiSE%s", type.Data())); - TH2D* dNdetadphi = (TH2D*)fOutputList->FindObject(Form("hdNdetadphiSE%s", type.Data())); + TH3D* flowHist = (TH3D*)fOutputList->FindObject(Form("hQ%dCumuHist%s", harmonic, type.Data())); + TH2D* dNdetaBindphi = (TH2D*)fOutputList->FindObject(Form("hdNdetaBindphiSE%s", type.Data())); + TH1D* dNdphi = 0; + if (!type.Contains("SPD")) dNdphi = (TH1D*)fOutputList->FindObject(Form("hdNdphiSE%s", type.Data())); + if ( type.Contains("SPD")) dNdphi = (TH1D*)fOutputList->FindObject("hdNdphiSE"); // We create the objects needed for the analysis - Double_t Mult = dNdphi->GetBinContent(0); + Double_t mult = dNdphi->GetBinContent(0); + if (type.Length() <= 1) fCent = dNdphi->GetBinContent(dNdphi->GetNbinsX()+1); + + TH1D* cent = (TH1D*)fOutputList->FindObject("cent"); + if (type.Length() <= 1) cent->Fill(fCent); - Double_t QnRe = 0, Q2nRe = 0, QnIm = 0, Q2nIm = 0; - Double_t pnRe = 0, p2nRe = 0, qnRe = 0, qnnRe = 0, pnIm = 0, p2nIm = 0, qnIm = 0, qnnIm = 0; + Double_t dQnRe = 0, dQ2nRe = 0, dQnIm = 0, dQ2nIm = 0; + Double_t pnRe = 0, p2nRe = 0, qnRe = 0, q2nRe = 0, pnIm = 0, p2nIm = 0, qnIm = 0, q2nIm = 0; Double_t avg2 = 0, avg4 = 0, avg2p = 0, avg4p = 0; Double_t w2avg2sq = 0, w2pavg2psq = 0, w4avg4sq = 0, w4pavg4psq = 0; Double_t w2w2pavg2avg2p = 0, w2w4avg2avg4 = 0, w2pw4pavg2pavg4p = 0; Double_t w2w4pavg2avg4p = 0, w4w2pavg4avg2p = 0, w4w4pavg4avg4p = 0; - Double_t CosPhi1Phi2 = 0, CosPhi1Phi2Phi3m = 0, CosPhi1Phi2Phi3p = 0; - Double_t SinPhi1Phi2 = 0, SinPhi1Phi2Phi3m = 0, SinPhi1Phi2Phi3p = 0; - Double_t Phii = 0; + Double_t cosPhi1Phi2 = 0, cosPhi1Phi2Phi3m = 0, cosPhi1Phi2Phi3p = 0; + Double_t sinPhi1Phi2 = 0, sinPhi1Phi2Phi3m = 0, sinPhi1Phi2Phi3p = 0; + Double_t phi = 0; Double_t multi = 0, multp = 0, mp = 0, mq = 0; - Double_t W2 = 0, W4 = 0, W2p = 0, W4p = 0; + Double_t w2 = 0, w4 = 0, w2p = 0, w4p = 0; // We loop over the data 1 time! - for (Int_t eta = 1; eta <= dNdetadphi->GetNbinsX(); eta++) { - // The values for each individual eta bins are reset + for (Int_t etaBin = 1; etaBin <= dNdetaBindphi->GetNbinsX(); etaBin++) { + // The values for each individual etaBin bins are reset mp = 0; pnRe = 0; p2nRe = 0; pnIm = 0; p2nIm = 0; - for (Int_t phii = 1; phii <= dNdphi->GetNbinsX()+1; phii++) { - Phii = dNdphi->GetXaxis()->GetBinCenter(phii); - multi = dNdphi->GetBinContent(phii); + for (Int_t phiN = 1; phiN <= dNdphi->GetNbinsX(); phiN++) { + phi = dNdphi->GetXaxis()->GetBinCenter(phiN); + multi = dNdphi->GetBinContent(phiN); - // In the phi loop on the first eta loop the integrated flow + // In the phi loop on the first etaBin loop the integrated flow // is calculated from the dNdphi histogram - if(eta == 1) { - QnRe += multi*TMath::Cos(n*Phii); - QnIm += multi*TMath::Sin(n*Phii); - Q2nRe += multi*TMath::Cos(2.*n*Phii); - Q2nIm += multi*TMath::Sin(2.*n*Phii); + if(etaBin == 1) { + 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); } - // For each eta bin the necesarry values for differential flow + // For each etaBin bin the necesarry values for differential flow // is calculated. Here is the loop over the phi's. - multp = dNdetadphi->GetBinContent(eta, phii); + multp = dNdetaBindphi->GetBinContent(etaBin, phiN); mp += multp; - pnRe += multp*TMath::Cos(n*Phii); - pnIm += multp*TMath::Sin(n*Phii); - p2nRe += multp*TMath::Cos(2.*n*Phii); - p2nIm += multp*TMath::Sin(2.*n*Phii); + 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); } - // The integrated flow is calculated - if (eta == 1) { - Double_t Eta = flowHist->GetXaxis()->GetBinCenter(0); + // The reference flow is calculated for the differential flow histograms + if (etaBin == 1) { + Double_t eta = flowHist->GetXaxis()->GetBinCenter(0); // 2-particle - W2 = Mult * (Mult - 1.); - avg2 = QnRe*QnRe + QnIm*QnIm - Mult; - avg2 /= W2; - w2avg2sq = W2 * avg2 * avg2; + w2 = mult * (mult - 1.); + avg2 = dQnRe*dQnRe + dQnIm*dQnIm - mult; + avg2 /= w2; + w2avg2sq = w2 * avg2 * avg2; - flowHist->Fill(Eta, 1, W2 * avg2); - flowHist->Fill(Eta, 3, W2); - flowHist->Fill(Eta, 4, w2avg2sq); + flowHist->Fill(eta, fCent, kRW2avg2, w2 * avg2); + flowHist->Fill(eta, fCent, kRW2, w2); + flowHist->Fill(eta, fCent, kRW2avg2sq, w2avg2sq); - flowHist->Fill(Eta, 26, QnRe); - flowHist->Fill(Eta, 27, QnIm); - flowHist->Fill(Eta, 28, Mult); + flowHist->Fill(eta, fCent, kQnRe, dQnRe); + flowHist->Fill(eta, fCent, kQnIm, dQnIm); + flowHist->Fill(eta, fCent, kRM, mult); // 4-particle - W4 = Mult * (Mult - 1.) * (Mult - 2.) * (Mult - 3.); - Double_t real = Q2nRe*QnRe*QnRe - Q2nRe*QnIm*QnIm + 2.*Q2nIm*QnRe*QnIm; + w4 = mult * (mult - 1.) * (mult - 2.) * (mult - 3.); + Double_t real = dQ2nRe*dQnRe*dQnRe - dQ2nRe*dQnIm*dQnIm + 2.*dQ2nIm*dQnRe*dQnIm; - avg4 = TMath::Power(QnRe*QnRe + QnIm*QnIm, 2); - avg4 += Q2nRe*Q2nRe + Q2nIm*Q2nIm - 2.*real; - avg4 -= 4.*(Mult - 2.)*(QnRe*QnRe + QnIm*QnIm) - 2.*Mult*(Mult - 3.); + avg4 = TMath::Power(dQnRe*dQnRe + dQnIm*dQnIm, 2); + avg4 += dQ2nRe*dQ2nRe + dQ2nIm*dQ2nIm - 2.*real; + avg4 -= 4.*(mult - 2.)*(dQnRe*dQnRe + dQnIm*dQnIm) - 2.*mult*(mult - 3.); - avg4 /= W4; - w4avg4sq = W4 * avg4 * avg4; - w2w4avg2avg4 = W2 * W4 * avg2 * avg4; - - flowHist->Fill(Eta, 9, W4 * avg4); - flowHist->Fill(Eta, 11, W4); - flowHist->Fill(Eta, 12, w4avg4sq); - flowHist->Fill(Eta, 14, W2 * W4); - flowHist->Fill(Eta, 15, w2w4avg2avg4); - - CosPhi1Phi2 = QnRe*QnRe - QnIm*QnIm - Q2nRe; - SinPhi1Phi2 = 2.*QnRe*QnIm - Q2nIm; + avg4 /= w4; + w4avg4sq = w4 * avg4 * avg4; + w2w4avg2avg4 = w2 * w4 * avg2 * avg4; + + flowHist->Fill(eta, fCent, kRW4avg4, w4 * avg4); + flowHist->Fill(eta, fCent, kRW4, w4); + flowHist->Fill(eta, fCent, kRW4avg4sq, w4avg4sq); + flowHist->Fill(eta, fCent, kRW2w4, w2 * w4); + flowHist->Fill(eta, fCent, kRW2w4avg2avg4, w2w4avg2avg4); + + cosPhi1Phi2 = dQnRe*dQnRe - dQnIm*dQnIm - dQ2nRe; + sinPhi1Phi2 = 2.*dQnRe*dQnIm - dQ2nIm; - CosPhi1Phi2Phi3m = TMath::Power(QnRe, 3) + QnRe*QnIm*QnIm; - CosPhi1Phi2Phi3m -= QnRe*Q2nRe + QnIm*Q2nIm + 2.*(Mult - 1.)*QnRe; - SinPhi1Phi2Phi3m = -TMath::Power(QnIm, 3) - QnRe*QnRe*QnIm; - SinPhi1Phi2Phi3m -= QnIm*Q2nRe - QnRe*Q2nIm + 2.*(Mult - 1.)*QnIm; + cosPhi1Phi2Phi3m = TMath::Power(dQnRe, 3) + dQnRe*dQnIm*dQnIm; + cosPhi1Phi2Phi3m -= dQnRe*dQ2nRe + dQnIm*dQ2nIm + 2.*(mult - 1.)*dQnRe; + + sinPhi1Phi2Phi3m = -TMath::Power(dQnIm, 3) - dQnRe*dQnRe*dQnIm; + sinPhi1Phi2Phi3m -= dQnIm*dQ2nRe - dQnRe*dQ2nIm - 2.*(mult - 1.)*dQnIm; - flowHist->Fill(Eta, 29, CosPhi1Phi2); - flowHist->Fill(Eta, 30, SinPhi1Phi2); - flowHist->Fill(Eta, 31, CosPhi1Phi2Phi3m); - flowHist->Fill(Eta, 32, SinPhi1Phi2Phi3m); - flowHist->Fill(Eta, 33, Mult*(Mult-1.)*(Mult-2.)); + flowHist->Fill(eta, fCent, kRCosphi1phi2, cosPhi1Phi2); + flowHist->Fill(eta, fCent, kRSinphi1phi2, sinPhi1Phi2); + flowHist->Fill(eta, fCent, kRCosphi1phi2phi3m, cosPhi1Phi2Phi3m); + flowHist->Fill(eta, fCent, kRSinphi1phi2phi3m, sinPhi1Phi2Phi3m); + flowHist->Fill(eta, fCent, kRMm1m2, mult*(mult-1.)*(mult-2.)); // Count number of events - flowHist->Fill(Eta, 0., 1.); - } // end of harmonics loop + flowHist->Fill(eta, fCent, 0., 1.); + } - // Differential flow calculations for each eta bin is done: + // Differential flow calculations for each etaBin bin is done: if (mp == 0) continue; - Double_t Eta = dNdetadphi->GetXaxis()->GetBinCenter(eta); + Double_t eta = dNdetaBindphi->GetXaxis()->GetBinCenter(etaBin); +// eta = TMath::Abs(eta); mq = mp; qnRe = pnRe; qnIm = pnIm; - qnnRe = p2nRe; - qnnIm = p2nIm; + q2nRe = p2nRe; + q2nIm = p2nIm; + if (type.Contains("SPD") || (type.Contains("") && eta >= 4.)) { + mq = 0; + qnRe = 0; + qnIm = 0; + q2nRe = 0; + q2nIm = 0; + } + + // Then the reference flow is calculated for each etaBin bin also, + // TODO: Find smart way to implement in above calculations... + + // 2-particle + w2p = mp * (mp - 1.); + avg2p = pnRe*pnRe + pnIm*pnIm - mp; + avg2p /= w2p; + w2pavg2psq = w2p * avg2p * avg2p; + + flowHist->Fill(eta, fCent, kRW2avg2, w2p * avg2p); + flowHist->Fill(eta, fCent, kRW2, w2p); + flowHist->Fill(eta, fCent, kRW2avg2sq, w2pavg2psq); + + flowHist->Fill(eta, fCent, kRM, mp); + + // 4-particle + w4p = mp * (mp - 1.) * (mp - 2.) * (mp - 3.); + Double_t real = p2nRe*pnRe*pnRe - p2nRe*pnIm*pnIm + 2.*p2nIm*pnRe*pnIm; + + avg4p = TMath::Power(pnRe*pnRe + pnIm*pnIm, 2); + avg4p += p2nRe*p2nRe + p2nIm*p2nIm - 2.*real; + avg4p -= 4.*(mp - 2.)*(pnRe*pnRe + pnIm*pnIm) - 2.*mp*(mp - 3.); + + avg4p /= w4p; + w4pavg4psq = w4p * avg4p * avg4p; + w2w4pavg2avg4p = w2p * w4p * avg2p * avg4p; + + flowHist->Fill(eta, fCent, kRW4avg4, w4p * avg4p); + flowHist->Fill(eta, fCent, kRW4, w4p); + flowHist->Fill(eta, fCent, kRW4avg4sq, w4pavg4psq); + flowHist->Fill(eta, fCent, kRW2w4, w2p * w4p); + flowHist->Fill(eta, fCent, kRW2w4avg2avg4, w2w4pavg2avg4p); + + cosPhi1Phi2 = pnRe*pnRe - pnIm*pnIm - p2nRe; + sinPhi1Phi2 = 2.*pnRe*pnIm - p2nIm; + + cosPhi1Phi2Phi3m = TMath::Power(pnRe, 3) + pnRe*pnIm*pnIm; + cosPhi1Phi2Phi3m -= pnRe*p2nRe + pnIm*p2nIm + 2.*(mp - 1.)*pnRe; + + sinPhi1Phi2Phi3m = -TMath::Power(pnIm, 3) - pnRe*pnRe*pnIm; + sinPhi1Phi2Phi3m -= pnIm*p2nRe - pnRe*p2nIm - 2.*(mp - 1.)*pnIm; + + flowHist->Fill(eta, fCent, kRCosphi1phi2, cosPhi1Phi2); + flowHist->Fill(eta, fCent, kRSinphi1phi2, sinPhi1Phi2); + flowHist->Fill(eta, fCent, kRCosphi1phi2phi3m, cosPhi1Phi2Phi3m); + flowHist->Fill(eta, fCent, kRSinphi1phi2phi3m, sinPhi1Phi2Phi3m); + flowHist->Fill(eta, fCent, kRMm1m2, mp*(mp-1.)*(mp-2.)); // 2-particle differential flow - W2p = mp * Mult - mq; - avg2p = pnRe*QnRe + pnIm*QnIm - mq; - avg2p /= W2p; - w2pavg2psq = W2p * avg2p * avg2p; - w2w2pavg2avg2p = W2 * W2p * avg2 * avg2p; + w2p = mp * mult - mq; + avg2p = pnRe*dQnRe + pnIm*dQnIm - mq; + avg2p /= w2p; + w2pavg2psq = w2p * avg2p * avg2p; + w2w2pavg2avg2p = w2 * w2p * avg2 * avg2p; - flowHist->Fill(Eta, 1, W2p * avg2p); - flowHist->Fill(Eta, 3, W2p); - flowHist->Fill(Eta, 4, w2pavg2psq); - flowHist->Fill(Eta, 6, w2w2pavg2avg2p); - flowHist->Fill(Eta, 8, W2 * W2p); + flowHist->Fill(eta, fCent, kW2avg2, w2p * avg2p); + flowHist->Fill(eta, fCent, kW2, w2p); + flowHist->Fill(eta, fCent, kW2avg2sq, w2pavg2psq); + flowHist->Fill(eta, fCent, kW2w2pavg2avg2p, w2w2pavg2avg2p); + flowHist->Fill(eta, fCent, kW2w2p, w2 * w2p); - flowHist->Fill(Eta, 26, pnRe); - flowHist->Fill(Eta, 27, pnIm); - flowHist->Fill(Eta, 28, mp); + flowHist->Fill(eta, fCent, kQnRe, pnRe); + flowHist->Fill(eta, fCent, kQnIm, pnIm); + flowHist->Fill(eta, fCent, kM, mp); // 4-particle differential flow - W4p = (mp * Mult - 3.*mq)*(Mult - 1.)*(Mult - 2.); - - avg4p = pnRe*QnRe*(QnRe*QnRe + QnIm*QnIm) + pnIm*QnIm*(QnRe*QnRe + QnIm*QnIm); - avg4p -= qnnRe*QnRe*QnRe - qnnRe*QnIm*QnIm + 2.*qnnIm*QnRe*QnIm; - avg4p -= pnRe*QnRe*Q2nRe - pnRe*QnIm*Q2nIm + pnIm*QnRe*Q2nIm + pnIm*QnIm*Q2nRe; - avg4p -= 2.*Mult*(pnRe*QnRe + pnIm*QnIm); - - avg4p += - 2.*mq*(QnRe*QnRe + QnIm*QnIm) + 7.*(qnRe*QnRe + qnIm*QnIm); - avg4p += - (QnRe*qnRe + QnIm*qnIm) + (qnnRe*Q2nRe + qnnIm*Q2nIm); - avg4p += 2.*(pnRe*QnRe + pnIm*QnIm) + 2.*mq*Mult - 6.*mq; - avg4p /= W4p; - - w4pavg4psq = W4p * avg4p * avg4p; - w2w4pavg2avg4p = W2 * W4p * avg2 * avg4p; - w4w2pavg4avg2p = W4 * W2p * avg4 * avg2p; - w4w4pavg4avg4p = W4 * W4p * avg4 * avg4p; - w2pw4pavg2pavg4p = W2p * W4p * avg2p * avg4p; - - flowHist->Fill(Eta, 9, W4p * avg4p); - flowHist->Fill(Eta, 11, W4p); - flowHist->Fill(Eta, 12, w4pavg4psq); - flowHist->Fill(Eta, 14, W2p * W4p); - flowHist->Fill(Eta, 15, w2pw4pavg2pavg4p); - flowHist->Fill(Eta, 17, W2 * W4p); - flowHist->Fill(Eta, 18, w2w4pavg2avg4p); - flowHist->Fill(Eta, 20, W4 * W2p); - flowHist->Fill(Eta, 21, w4w2pavg4avg2p); - flowHist->Fill(Eta, 23, W4 * W4p); - flowHist->Fill(Eta, 24, w4w4pavg4avg4p); - - CosPhi1Phi2 = pnRe*QnRe - pnIm*QnIm - qnnRe; - SinPhi1Phi2 = pnRe*QnIm + pnIm*QnRe - qnnIm; - - CosPhi1Phi2Phi3p = pnRe*(QnRe*QnRe + QnIm*QnIm - Mult); - CosPhi1Phi2Phi3p -= qnnRe*QnRe - qnnIm*QnIm + mq*QnRe - 2.*qnRe; - SinPhi1Phi2Phi3p = pnIm*(QnRe*QnRe + QnIm*QnIm - Mult); - SinPhi1Phi2Phi3p -= qnnIm*QnRe - qnnRe*QnIm + mq*QnIm - 2.*qnIm; - - CosPhi1Phi2Phi3m = pnRe*(QnRe*QnRe - QnIm*QnIm) + 2.*pnIm*QnRe*QnIm; - CosPhi1Phi2Phi3m -= pnRe*Q2nRe + pnIm*Q2nIm + 2.*mq*QnRe - 2.*qnRe; - SinPhi1Phi2Phi3m = pnIm*(QnRe*QnRe - QnIm*QnIm) - 2.*pnRe*QnRe*QnIm; - SinPhi1Phi2Phi3m += - pnIm*Q2nRe + pnRe*Q2nIm + 2.*mq*QnIm - 2.*qnIm; - - flowHist->Fill(Eta, 29, CosPhi1Phi2); - flowHist->Fill(Eta, 30, SinPhi1Phi2); - flowHist->Fill(Eta, 31, CosPhi1Phi2Phi3m); - flowHist->Fill(Eta, 32, SinPhi1Phi2Phi3m); - flowHist->Fill(Eta, 33, (mp*Mult-2.*mq)*(Mult-1.)); - flowHist->Fill(Eta, 34, CosPhi1Phi2Phi3p); - flowHist->Fill(Eta, 35, SinPhi1Phi2Phi3p); + w4p = (mp * mult - 3.*mq)*(mult - 1.)*(mult - 2.); + + avg4p = pnRe*dQnRe*(dQnRe*dQnRe + dQnIm*dQnIm) + pnIm*dQnIm*(dQnRe*dQnRe + dQnIm*dQnIm); + avg4p -= q2nRe*dQnRe*dQnRe - q2nRe*dQnIm*dQnIm + 2.*q2nIm*dQnRe*dQnIm; + avg4p -= pnRe*dQnRe*dQ2nRe - pnRe*dQnIm*dQ2nIm + pnIm*dQnRe*dQ2nIm + pnIm*dQnIm*dQ2nRe; + avg4p -= 2.*mult*(pnRe*dQnRe + pnIm*dQnIm); + + avg4p += - 2.*mq*(dQnRe*dQnRe + dQnIm*dQnIm) + 7.*(qnRe*dQnRe + qnIm*dQnIm); + avg4p += - (dQnRe*qnRe + dQnIm*qnIm) + (q2nRe*dQ2nRe + q2nIm*dQ2nIm); + avg4p += 2.*(pnRe*dQnRe + pnIm*dQnIm) + 2.*mq*mult - 6.*mq; + avg4p /= w4p; + + w4pavg4psq = w4p * avg4p * avg4p; + w2w4pavg2avg4p = w2 * w4p * avg2 * avg4p; + w4w2pavg4avg2p = w4 * w2p * avg4 * avg2p; + w4w4pavg4avg4p = w4 * w4p * avg4 * avg4p; + w2pw4pavg2pavg4p = w2p * w4p * avg2p * avg4p; + + flowHist->Fill(eta, fCent, kW4avg4, w4p * avg4p); + flowHist->Fill(eta, fCent, kW4, w4p); + flowHist->Fill(eta, fCent, kW4avg4sq, w4pavg4psq); + flowHist->Fill(eta, fCent, kW2w4, w2p * w4p); + flowHist->Fill(eta, fCent, kW2w4avg2avg4, w2pw4pavg2pavg4p); + flowHist->Fill(eta, fCent, kW2w4p, w2 * w4p); + flowHist->Fill(eta, fCent, kW2w4pavg2avg4p, w2w4pavg2avg4p); + flowHist->Fill(eta, fCent, kW4w2p, w4 * w2p); + flowHist->Fill(eta, fCent, kW4w2pavg4avg2p, w4w2pavg4avg2p); + flowHist->Fill(eta, fCent, kW4w4p, w4 * w4p); + flowHist->Fill(eta, fCent, kW4w4pavg4avg4p, w4w4pavg4avg4p); + + cosPhi1Phi2 = pnRe*dQnRe - pnIm*dQnIm - q2nRe; + sinPhi1Phi2 = pnRe*dQnIm + pnIm*dQnRe - q2nIm; + + cosPhi1Phi2Phi3p = pnRe*(dQnRe*dQnRe + dQnIm*dQnIm - mult); + cosPhi1Phi2Phi3p -= q2nRe*dQnRe - q2nIm*dQnIm + mq*dQnRe - 2.*qnRe; + sinPhi1Phi2Phi3p = pnIm*(dQnRe*dQnRe + dQnIm*dQnIm - mult); + sinPhi1Phi2Phi3p -= q2nIm*dQnRe - q2nRe*dQnIm + mq*dQnIm - 2.*qnIm; + + cosPhi1Phi2Phi3m = pnRe*(dQnRe*dQnRe - dQnIm*dQnIm) + 2.*pnIm*dQnRe*dQnIm; + cosPhi1Phi2Phi3m -= pnRe*dQ2nRe + pnIm*dQ2nIm + 2.*mq*dQnRe - 2.*qnRe; + sinPhi1Phi2Phi3m = pnIm*(dQnRe*dQnRe - dQnIm*dQnIm) - 2.*pnRe*dQnRe*dQnIm; + sinPhi1Phi2Phi3m += - pnIm*dQ2nRe + pnRe*dQ2nIm + 2.*mq*dQnIm - 2.*qnIm; + + flowHist->Fill(eta, fCent, kCosphi1phi2, cosPhi1Phi2); + flowHist->Fill(eta, fCent, kSinphi1phi2, sinPhi1Phi2); + flowHist->Fill(eta, fCent, kCosphi1phi2phi3m, cosPhi1Phi2Phi3m); + flowHist->Fill(eta, fCent, kSinphi1phi2phi3m, sinPhi1Phi2Phi3m); + flowHist->Fill(eta, fCent, kMm1m2, (mp*mult-2.*mq)*(mult-1.)); + flowHist->Fill(eta, fCent, kCosphi1phi2phi3p, cosPhi1Phi2Phi3p); + flowHist->Fill(eta, fCent, kSinphi1phi2phi3p, sinPhi1Phi2Phi3p); } @@ -458,243 +553,291 @@ void AliForwardFlowTaskQC::Terminate(Option_t */*option*/) // Parameters: // option Not used // +// Double_t x[6] = { 0., 1., 2., 3., 4.5, 6. }; + + TH3D* hcumulantsHist = 0; + TH2D* cumulantsHist = new TH2D("tmphist", "tmphist", fEtaBins, -4, 6, 40, 0.5, 40.5); + TH1D* cumulant2refHist = 0; + TH1D* cumulant4refHist = 0; + TH1D* cumulant2diffHist = 0; + TH1D* cumulant4diffHist = 0; + + // For flow calculations + Double_t two = 0, qc2 = 0, vTwo2 = 0, four = 0, qc4 = 0, vTwo4 = 0; + Double_t twoPrime = 0, qc2Prime = 0, vTwo2diff = 0, fourPrime = 0, qc4Prime = 0, vTwo4diff = 0; + Double_t w2 = 0, w4 = 0, w2p = 0, w4p = 0, sqrtW2sq = 0, sqrtW2psq = 0, w2W2p = 0; + Double_t w2W4 = 0, w2W4p = 0, w4W2p = 0, w4W4p = 0, w2pW4p = 0; + Double_t sqrtW4sq = 0, sqrtW4psq = 0; + Double_t w2avg2 = 0, w2pavg2p = 0, w4avg4 = 0, w4pavg4p = 0; + Double_t cosP1nPhi = 0, sinP1nPhi = 0, mult = 0, cosP1nPhi1P1nPhi2 = 0, sinP1nPhi1P1nPhi2 = 0; + Double_t cosP1nPhi1M1nPhi2M1nPhi3 = 0, sinP1nPhi1M1nPhi2M1nPhi3 = 0, multm1m2 = 0; + Double_t cosP1nPsi = 0, sinP1nPsi = 0, mp = 0, cosP1nPsi1P1nPhi2 = 0, sinP1nPsi1P1nPhi2 = 0; + Double_t cosP1nPsi1M1nPhi2M1nPhi3 = 0, sinP1nPsi1M1nPhi2M1nPhi3 = 0, mpqMult = 0; + Double_t cosP1nPsi1P1nPhi2M1nPhi3 = 0, sinP1nPsi1P1nPhi2M1nPhi3 = 0; + + // For error calculations + Double_t w2avg2sq = 0, w2W2pavg2avg2p = 0, w2pavg2psq = 0; + Double_t w4avg4sq = 0, w2W4avg2avg4 = 0, w4W4pavg4avg4p = 0, w4pavg4psq = 0; + Double_t w2W4pavg2avg4p = 0, w4W2pavg4avg2p = 0; + Double_t stwosq = 0, stwoPrimesq = 0, sfoursq = 0, sfourPrimesq = 0; + Double_t vTwo2err = 0, vTwo2diffErr = 0, vTwo4err = 0, vTwo4diffErr = 0; + Double_t cov22p = 0, cov24 = 0, cov24p = 0, cov42p = 0, cov44p = 0, cov2p2np = 0; + Double_t w2pW4pavg2pavg4p = 0; + - TH2D* cumulantsHist; - TH1D* cumulant2Hist; - TH1D* cumulant4Hist; - - Int_t nLoops = (fMC ? 3 : 1); + Int_t nLoops = (fMC ? 4 : 2); + TString type = ""; + Int_t y[11] = { 0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 100 }; + // Do a loop over the difference analysis types, calculating flow - // 1 loop for real data, 3 for MC data + // 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++) { - TString type; if (loop == 1) type = ""; - if (loop == 2) type = "MC"; - if (loop == 3) type = "TrRef"; + if (loop == 2) type = "SPD"; + if (loop == 3) type = "MC"; + if (loop == 4) type = "TrRef"; for (Int_t n = 1; n <= 4; n++) { if (!fv[n]) continue; - cumulantsHist = (TH2D*)fOutputList->FindObject(Form("hQ%dCumuHist%s", n, type.Data())); - cumulant2Hist = (TH1D*)fOutputList->FindObject(Form("hQ%dCumulant2Flow%s", n, type.Data())); - cumulant4Hist = (TH1D*)fOutputList->FindObject(Form("hQ%dCumulant4Flow%s", n, type.Data())); - - // For flow calculations - Double_t Avg2 = 0, c2 = 0, vTwo2 = 0, Avg4 = 0, c4 = 0, vTwo4 = 0; - Double_t Avg2p = 0, d2 = 0, vTwo2diff = 0, Avg4p = 0, d4 = 0, vTwo4diff = 0; - Double_t W2 = 0, W4 = 0, W2p = 0, W4p = 0, sqrtW2sq = 0, sqrtW2psq = 0, W2W2p = 0; - Double_t W2W4 = 0, W2W4p = 0, W4W2p = 0, W4W4p = 0, W2pW4p = 0; - Double_t sqrtW4sq = 0, sqrtW4psq = 0; - Double_t W2avg2 = 0, W2pavg2p = 0, W4avg4 = 0, W4pavg4p = 0; - Double_t AvgCos2Phi = 0, AvgSin2Phi = 0, Mult = 0, AvgCos2Phi1Phi2 = 0, AvgSin2Phi1Phi2 = 0; - Double_t AvgCos2Phi1Phi2Phi3m = 0, AvgSin2Phi1Phi2Phi3m = 0, Multm1m2 = 0; - Double_t AvgCos2Psi = 0, AvgSin2Psi = 0, mp = 0, AvgCos2Psi1Phi2 = 0, AvgSin2Psi1Phi2 = 0; - Double_t AvgCos2Psi1Phi2Phi3m = 0, AvgSin2Psi1Phi2Phi3m = 0, mpqMult = 0; - Double_t AvgCos2Psi1Phi2Phi3p = 0, AvgSin2Psi1Phi2Phi3p = 0; - - // For error calculations - Double_t W2avg2sq = 0, W2W2pavg2avg2p = 0, W2pavg2psq = 0; - Double_t W4avg4sq = 0, W2W4avg2avg4 = 0, W4W4pavg4avg4p = 0, W4pavg4psq = 0; - Double_t W2W4pavg2avg4p = 0, W4W2pavg4avg2p = 0; - Double_t sAvg2sq = 0, sAvg2psq = 0, sAvg4sq = 0, sAvg4psq = 0; - Double_t vTwo2err = 0, vTwo2diffErr = 0, vTwo4err = 0, vTwo4diffErr = 0; - Double_t Cov22p = 0, Cov24 = 0, Cov24p = 0, Cov42p = 0, Cov44p = 0, Cov2p2np = 0; - Double_t W2pW4pavg2pavg4p = 0; - - for (Int_t eta = 0; eta <= cumulantsHist->GetNbinsX(); eta++) { - if (eta == 0) { - // 2-particle reference flow - W2avg2 = cumulantsHist->GetBinContent(eta, 1); - if (!W2avg2) continue; - W2 = cumulantsHist->GetBinContent(eta, 3); - AvgCos2Phi = cumulantsHist->GetBinContent(eta, 26); - AvgSin2Phi = cumulantsHist->GetBinContent(eta, 27); - Mult = cumulantsHist->GetBinContent(eta, 28); - AvgCos2Phi /= Mult; - AvgSin2Phi /= Mult; - Avg2 = W2avg2 / W2; - c2 = Avg2 - TMath::Power(AvgCos2Phi, 2) - TMath::Power(AvgSin2Phi, 2); - vTwo2 = TMath::Sqrt(c2); - cumulant2Hist->SetBinContent(eta, vTwo2); - - // 4-particle reference flow - W4avg4 = cumulantsHist->GetBinContent(eta, 9); - W4 = cumulantsHist->GetBinContent(eta, 11); - AvgCos2Phi1Phi2 = cumulantsHist->GetBinContent(eta, 29); - AvgSin2Phi1Phi2 = cumulantsHist->GetBinContent(eta, 30); - AvgCos2Phi1Phi2 /= W2; - AvgSin2Phi1Phi2 /= W2; - AvgCos2Phi1Phi2Phi3m = cumulantsHist->GetBinContent(eta, 31); - AvgSin2Phi1Phi2Phi3m = cumulantsHist->GetBinContent(eta, 32); - Multm1m2 = cumulantsHist->GetBinContent(eta, 33); - AvgCos2Phi1Phi2Phi3m /= Multm1m2; - AvgSin2Phi1Phi2Phi3m /= Multm1m2; - Avg4 = W4avg4 / W4; - c4 = Avg4 - 2. * Avg2 * Avg2; - c4 -= 4.*AvgCos2Phi*AvgCos2Phi1Phi2Phi3m; - c4 += 4.*AvgSin2Phi*AvgSin2Phi1Phi2Phi3m; - c4 -= TMath::Power(AvgCos2Phi1Phi2, 2) + TMath::Power(AvgSin2Phi1Phi2 , 2); - c4 += 4.*AvgCos2Phi1Phi2*(TMath::Power(AvgCos2Phi, 2) - TMath::Power(AvgSin2Phi, 2)); - c4 += 8.*AvgSin2Phi1Phi2*AvgSin2Phi*AvgCos2Phi; - c4 += 8.*Avg2*(TMath::Power(AvgCos2Phi, 2) + TMath::Power(AvgSin2Phi, 2)); - c4 -= 6.*TMath::Power(TMath::Power(AvgCos2Phi, 2)+TMath::Power(AvgSin2Phi, 2), 2); - - vTwo4 = TMath::Power(-c4, 0.25); - cumulant4Hist->SetBinContent(eta, vTwo4); + hcumulantsHist = (TH3D*)fOutputList->FindObject(Form("hQ%dCumuHist%s", n, type.Data())); - // 2-particle reference flow error calculations - W2avg2sq = cumulantsHist->GetBinContent(eta, 4); - sqrtW2sq = cumulantsHist->GetBinError(eta, 3); - - sAvg2sq = VarSQ(W2avg2sq, Avg2, W2, W2avg2, sqrtW2sq); - vTwo2err = sqrtW2sq * TMath::Sqrt(sAvg2sq) / (2. * TMath::Sqrt(Avg2) * W2); - cumulant2Hist->SetBinError(eta, vTwo2err); - - // 4-particle reference flow error calculations - W4avg4sq = cumulantsHist->GetBinContent(eta, 12); - sqrtW4sq = cumulantsHist->GetBinError(eta, 11); - W2W4 = cumulantsHist->GetBinContent(eta, 14); - W2W4avg2avg4 = cumulantsHist->GetBinContent(eta, 15); - - sAvg4sq = VarSQ(W4avg4sq, Avg4, W4, W4avg4, sqrtW4sq); - Cov24 = CovXY(W2W4avg2avg4, W2W4, Avg2*Avg4, W2, W4); - - vTwo4err = Avg2*Avg2 * TMath::Power(sqrtW2sq, 2) * sAvg2sq / (W2*W2); - vTwo4err += TMath::Power(sqrtW4sq, 2) * sAvg4sq / (16. * W4*W4); - vTwo4err -= Avg2 * W2W4 * Cov24 / (2. * W2 * W4); - vTwo4err /= TMath::Power(2. * Avg2*Avg2 - Avg4, 1.5); - vTwo4err = TMath::Sqrt(vTwo4err); - cumulant4Hist->SetBinError(eta, vTwo4err); - - continue; + // Centrality loop + for (Int_t c = -2; c < hcumulantsHist->GetNbinsY(); c++) { + if (c == -2) { + hcumulantsHist->GetYaxis()->SetRange(0, 9); + cumulantsHist = (TH2D*)hcumulantsHist->Project3D("zx oe"); + + cumulant2refHist = (TH1D*)fOutputList->FindObject(Form("hQ%dCumulant2RefFlow%s_mb", n, type.Data())); + cumulant4refHist = (TH1D*)fOutputList->FindObject(Form("hQ%dCumulant4RefFlow%s_mb", n, type.Data())); + cumulant2diffHist = (TH1D*)fOutputList->FindObject(Form("hQ%dCumulant2DiffFlow%s_mb", n, type.Data())); + cumulant4diffHist = (TH1D*)fOutputList->FindObject(Form("hQ%dCumulant4DiffFlow%s_mb", n, type.Data())); + } + else if (c == -1) { + hcumulantsHist->GetYaxis()->SetRange(1, 5); + cumulantsHist = (TH2D*)hcumulantsHist->Project3D("zx oe"); + + cumulant2refHist = (TH1D*)fOutputList->FindObject(Form("hQ%dCumulant2RefFlow%s_0_40", n, type.Data())); + cumulant4refHist = (TH1D*)fOutputList->FindObject(Form("hQ%dCumulant4RefFlow%s_0_40", n, type.Data())); + cumulant2diffHist = (TH1D*)fOutputList->FindObject(Form("hQ%dCumulant2DiffFlow%s_0_40", n, type.Data())); + cumulant4diffHist = (TH1D*)fOutputList->FindObject(Form("hQ%dCumulant4DiffFlow%s_0_40", n, type.Data())); + } + else { + hcumulantsHist->GetYaxis()->SetRange(c+1, c+1); + cumulantsHist = (TH2D*)hcumulantsHist->Project3D("zx e"); + cumulant2refHist = (TH1D*)fOutputList->FindObject(Form("hQ%dCumulant2RefFlow%s_%d_%d", n, type.Data(), y[c], y[c+1])); + cumulant4refHist = (TH1D*)fOutputList->FindObject(Form("hQ%dCumulant4RefFlow%s_%d_%d", n, type.Data(), y[c], y[c+1])); + cumulant2diffHist = (TH1D*)fOutputList->FindObject(Form("hQ%dCumulant2DiffFlow%s_%d_%d", n, type.Data(), y[c], y[c+1])); + cumulant4diffHist = (TH1D*)fOutputList->FindObject(Form("hQ%dCumulant4DiffFlow%s_%d_%d", n, type.Data(), y[c], y[c+1])); } - // 2-particle differential flow - W2pavg2p = cumulantsHist->GetBinContent(eta, 1); - if (!W2pavg2p) continue; - W2p = cumulantsHist->GetBinContent(eta, 3); - AvgCos2Psi = cumulantsHist->GetBinContent(eta, 26); - AvgSin2Psi = cumulantsHist->GetBinContent(eta, 27); - mp = cumulantsHist->GetBinContent(eta, 28); - AvgCos2Psi /= mp; - AvgSin2Psi /= mp; - Avg2p = W2pavg2p / W2p; - d2 = Avg2p - AvgCos2Psi*AvgCos2Phi - AvgSin2Psi*AvgSin2Phi; - vTwo2diff = d2 / TMath::Sqrt(c2); - cumulant2Hist->SetBinContent(eta, vTwo2diff); - - // 4-particle differential flow - W4pavg4p = cumulantsHist->GetBinContent(eta, 9); - W4p = cumulantsHist->GetBinContent(eta, 11); - AvgCos2Psi1Phi2 = cumulantsHist->GetBinContent(eta, 29); - AvgSin2Psi1Phi2 = cumulantsHist->GetBinContent(eta, 30); - AvgCos2Psi1Phi2 /= W2p; - AvgSin2Psi1Phi2 /= W2p; - AvgCos2Psi1Phi2Phi3m = cumulantsHist->GetBinContent(eta, 31); - AvgSin2Psi1Phi2Phi3m = cumulantsHist->GetBinContent(eta, 32); - mpqMult = cumulantsHist->GetBinContent(eta, 33); - AvgCos2Psi1Phi2Phi3m /= mpqMult; - AvgSin2Psi1Phi2Phi3m /= mpqMult; - AvgCos2Psi1Phi2Phi3p = cumulantsHist->GetBinContent(eta, 34); - AvgSin2Psi1Phi2Phi3p = cumulantsHist->GetBinContent(eta, 35); - AvgCos2Psi1Phi2Phi3p /= mpqMult; - AvgSin2Psi1Phi2Phi3p /= mpqMult; - - Avg4p = W4pavg4p / W4p; - d4 = Avg4p - 2. * Avg2p * Avg2; - d4 -= AvgCos2Psi*AvgCos2Phi1Phi2Phi3m; - d4 += AvgSin2Psi*AvgSin2Phi1Phi2Phi3m; - d4 -= AvgCos2Phi*AvgCos2Psi1Phi2Phi3m; - d4 += AvgSin2Phi*AvgSin2Psi1Phi2Phi3m; - d4 -= 2.*AvgCos2Phi*AvgCos2Psi1Phi2Phi3p; - d4 -= 2.*AvgSin2Phi*AvgSin2Psi1Phi2Phi3p; - d4 -= AvgCos2Psi1Phi2*AvgCos2Phi1Phi2; - d4 -= AvgSin2Psi1Phi2*AvgSin2Phi1Phi2; - d4 += 2.*AvgCos2Phi1Phi2*(AvgCos2Psi*AvgCos2Phi - AvgSin2Psi*AvgSin2Phi); - d4 += 2.*AvgSin2Phi1Phi2*(AvgCos2Psi*AvgSin2Phi + AvgSin2Psi*AvgCos2Phi); - d4 += 4.*Avg2*(AvgCos2Psi*AvgCos2Phi + AvgSin2Psi*AvgSin2Phi); - d4 += 2.*AvgCos2Psi1Phi2*(TMath::Power(AvgCos2Phi, 2) - TMath::Power(AvgSin2Phi, 2)); - d4 += 4.*AvgSin2Psi1Phi2*AvgCos2Phi*AvgSin2Phi; - d4 += 4.*Avg2p*(TMath::Power(AvgCos2Phi, 2) + TMath::Power(AvgSin2Phi, 2)); - d4 -= 6.*(TMath::Power(AvgCos2Phi, 2) - TMath::Power(AvgSin2Phi, 2)) - *(AvgCos2Psi*AvgCos2Phi-AvgSin2Psi*AvgSin2Phi); - d4 -= 12.*AvgCos2Phi*AvgSin2Phi*(AvgSin2Psi*AvgCos2Phi+AvgCos2Psi*AvgSin2Phi); - - vTwo4diff = - d4 / TMath::Power(-c4, 0.75); - cumulant4Hist->SetBinContent(eta, vTwo4diff); - - // 2-particle differential flow error calculations - W2pavg2psq = cumulantsHist->GetBinContent(eta, 4); - sqrtW2psq = cumulantsHist->GetBinError(eta, 3); - W2W2pavg2avg2p = cumulantsHist->GetBinContent(eta, 6); - W2W2p = cumulantsHist->GetBinContent(eta, 8); + Bool_t refDone = kFALSE; + + for (Int_t etaBin = 1; etaBin <= cumulantsHist->GetNbinsX(); etaBin++) { + if (refDone == kFALSE || etaBin == 0) { + // 2-particle reference flow + w2avg2 = cumulantsHist->GetBinContent(etaBin, kRW2avg2); + if (etaBin == cumulantsHist->GetNbinsX()) { + refDone = kTRUE; + etaBin = -1; + continue; + } + if (!w2avg2) continue; + w2 = cumulantsHist->GetBinContent(etaBin, kRW2); + cosP1nPhi = cumulantsHist->GetBinContent(etaBin, kQnRe); + sinP1nPhi = cumulantsHist->GetBinContent(etaBin, kQnIm); + mult = cumulantsHist->GetBinContent(etaBin, kRM); + cosP1nPhi /= mult; + sinP1nPhi /= mult; + two = w2avg2 / w2; + qc2 = two - TMath::Power(cosP1nPhi, 2) - TMath::Power(sinP1nPhi, 2); + vTwo2 = TMath::Sqrt(qc2); + if (etaBin == 0) cumulant2diffHist->SetBinContent(etaBin, vTwo2); + if (etaBin > 0) cumulant2refHist->SetBinContent(etaBin, vTwo2); + + // 4-particle reference flow + w4avg4 = cumulantsHist->GetBinContent(etaBin, kRW4avg4); + w4 = cumulantsHist->GetBinContent(etaBin, kRW4); + cosP1nPhi1P1nPhi2 = cumulantsHist->GetBinContent(etaBin, kRCosphi1phi2); + sinP1nPhi1P1nPhi2 = cumulantsHist->GetBinContent(etaBin, kRSinphi1phi2); + cosP1nPhi1P1nPhi2 /= w2; + sinP1nPhi1P1nPhi2 /= w2; + cosP1nPhi1M1nPhi2M1nPhi3 = cumulantsHist->GetBinContent(etaBin, kRCosphi1phi2phi3m); + sinP1nPhi1M1nPhi2M1nPhi3 = cumulantsHist->GetBinContent(etaBin, kRSinphi1phi2phi3m); + multm1m2 = cumulantsHist->GetBinContent(etaBin, kRMm1m2); + cosP1nPhi1M1nPhi2M1nPhi3 /= multm1m2; + sinP1nPhi1M1nPhi2M1nPhi3 /= multm1m2; + four = w4avg4 / w4; + qc4 = four-2.*pow(two,2.) + - 4.*cosP1nPhi*cosP1nPhi1M1nPhi2M1nPhi3 + + 4.*sinP1nPhi*sinP1nPhi1M1nPhi2M1nPhi3-pow(cosP1nPhi1P1nPhi2,2.)-pow(sinP1nPhi1P1nPhi2,2.) + + 4.*cosP1nPhi1P1nPhi2*(pow(cosP1nPhi,2.)-pow(sinP1nPhi,2.))+8.*sinP1nPhi1P1nPhi2*sinP1nPhi*cosP1nPhi + + 8.*two*(pow(cosP1nPhi,2.)+pow(sinP1nPhi,2.))-6.*pow((pow(cosP1nPhi,2.)+pow(sinP1nPhi,2.)),2.); + + + vTwo4 = TMath::Power(-qc4, 0.25); + if (etaBin == 0) cumulant4diffHist->SetBinContent(etaBin, vTwo4); + if (etaBin > 0) cumulant4refHist->SetBinContent(etaBin, vTwo4); + + // 2-particle reference flow error calculations + w2avg2sq = cumulantsHist->GetBinContent(etaBin, kRW2avg2sq); + sqrtW2sq = cumulantsHist->GetBinError(etaBin, kRW2); + + stwosq = VarSQ(w2avg2sq, two, w2, w2avg2, sqrtW2sq); + vTwo2err = sqrtW2sq * TMath::Sqrt(stwosq) / (2. * TMath::Sqrt(two) * w2); + if (etaBin == 0) cumulant2diffHist->SetBinError(etaBin, vTwo2err); + if (etaBin > 0) cumulant2refHist->SetBinError(etaBin, vTwo2err); + + // 4-particle reference flow error calculations + w4avg4sq = cumulantsHist->GetBinContent(etaBin, kRW4avg4sq); + sqrtW4sq = cumulantsHist->GetBinError(etaBin, kRW4); + w2W4 = cumulantsHist->GetBinContent(etaBin, kRW2w4); + w2W4avg2avg4 = cumulantsHist->GetBinContent(etaBin, kRW2w4avg2avg4); + + sfoursq = VarSQ(w4avg4sq, four, w4, w4avg4, sqrtW4sq); + cov24 = CovXY(w2W4avg2avg4, w2W4, two*four, w2, w4); - Cov22p = CovXY(W2W2pavg2avg2p, W2W2p, Avg2*Avg2p, W2, W2p); - sAvg2psq = VarSQ(W2pavg2psq, Avg2p, W2p, W2pavg2p, sqrtW2psq); + vTwo4err = two*two * TMath::Power(sqrtW2sq, 2) * stwosq / (w2*w2); + vTwo4err += TMath::Power(sqrtW4sq, 2) * sfoursq / (16. * w4*w4); + vTwo4err -= two * w2W4 * cov24 / (2. * w2 * w4); + vTwo4err /= TMath::Power(2. * two*two - four, 1.5); + vTwo4err = TMath::Sqrt(vTwo4err); + if (etaBin == 0) cumulant4diffHist->SetBinError(etaBin, vTwo4err); + if (etaBin > 0) cumulant4refHist->SetBinError(etaBin, vTwo4err); + + continue; + } + + // 2-particle differential flow + w2pavg2p = cumulantsHist->GetBinContent(etaBin, kW2avg2); + if (!w2pavg2p) continue; + w2p = cumulantsHist->GetBinContent(etaBin, kW2); + cosP1nPsi = cumulantsHist->GetBinContent(etaBin, kQnRe); + sinP1nPsi = cumulantsHist->GetBinContent(etaBin, kQnIm); + mp = cumulantsHist->GetBinContent(etaBin, kM); + cosP1nPsi /= mp; + sinP1nPsi /= mp; + twoPrime = w2pavg2p / w2p; + qc2Prime = twoPrime - sinP1nPsi*sinP1nPhi - cosP1nPsi*cosP1nPhi; + vTwo2diff = qc2Prime / TMath::Sqrt(qc2); + cumulant2diffHist->SetBinContent(etaBin, vTwo2diff); - vTwo2diffErr = Avg2p*Avg2p*TMath::Power(sqrtW2psq, 2)*sAvg2sq/(W2*W2); - vTwo2diffErr += 4.*Avg2*Avg2*TMath::Power(sqrtW2psq, 2)*sAvg2psq/(W2p*W2p); - vTwo2diffErr -= 4.*Avg2*Avg2p*W2W2p*Cov22p/(W2*W2p); - vTwo2diffErr /= (4. * TMath::Power(Avg2, 3)); - vTwo2diffErr = TMath::Sqrt(vTwo2diffErr); - cumulant2Hist->SetBinError(eta, vTwo2diffErr); - - // 4-particle differential flow error calculations - sqrtW4psq = cumulantsHist->GetBinError(eta, 11); - W4pavg4psq = cumulantsHist->GetBinContent(eta, 12); - W2pW4p = cumulantsHist->GetBinContent(eta, 14); - W2pW4pavg2pavg4p = cumulantsHist->GetBinContent(eta, 15); - W2W4p = cumulantsHist->GetBinContent(eta, 17); - W2W4pavg2avg4p = cumulantsHist->GetBinContent(eta, 18); - W4W2p = cumulantsHist->GetBinContent(eta, 20); - W4W2pavg4avg2p = cumulantsHist->GetBinContent(eta, 21); - W4W4p = cumulantsHist->GetBinContent(eta, 23); - W4W4pavg4avg4p = cumulantsHist->GetBinContent(eta, 24); - - sAvg4psq = VarSQ(W4pavg4psq, Avg4p, W4p, W4pavg4p, sqrtW4psq); - Cov24p = CovXY(W2W4pavg2avg4p, W2W4p, Avg2*Avg4p, W2, W4p); - Cov42p = CovXY(W4W2pavg4avg2p, W4W2p, Avg4*Avg2p, W4, W2p); - Cov44p = CovXY(W4W4pavg4avg4p, W4W4p, Avg4*Avg4p, W4, W4p); - Cov2p2np = CovXY(W2pW4pavg2pavg4p, W2pW4p, Avg2p*Avg4p, W2p, W4p); + // 4-particle differential flow + w4pavg4p = cumulantsHist->GetBinContent(etaBin, kW4avg4); + w4p = cumulantsHist->GetBinContent(etaBin, kW4); + cosP1nPsi1P1nPhi2 = cumulantsHist->GetBinContent(etaBin, kCosphi1phi2); + sinP1nPsi1P1nPhi2 = cumulantsHist->GetBinContent(etaBin, kSinphi1phi2); + cosP1nPsi1P1nPhi2 /= w2p; + sinP1nPsi1P1nPhi2 /= w2p; + cosP1nPsi1M1nPhi2M1nPhi3 = cumulantsHist->GetBinContent(etaBin, kCosphi1phi2phi3m); + sinP1nPsi1M1nPhi2M1nPhi3 = cumulantsHist->GetBinContent(etaBin, kSinphi1phi2phi3m); + mpqMult = cumulantsHist->GetBinContent(etaBin, kMm1m2); + cosP1nPsi1M1nPhi2M1nPhi3 /= mpqMult; + sinP1nPsi1M1nPhi2M1nPhi3 /= mpqMult; + cosP1nPsi1P1nPhi2M1nPhi3 = cumulantsHist->GetBinContent(etaBin, kCosphi1phi2phi3p); + sinP1nPsi1P1nPhi2M1nPhi3 = cumulantsHist->GetBinContent(etaBin, kSinphi1phi2phi3p); + cosP1nPsi1P1nPhi2M1nPhi3 /= mpqMult; + sinP1nPsi1P1nPhi2M1nPhi3 /= mpqMult; + + fourPrime = w4pavg4p / 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*(pow(cosP1nPhi,2.)-pow(sinP1nPhi,2.)) + + 4.*sinP1nPsi1P1nPhi2*cosP1nPhi*sinP1nPhi + + 4.*twoPrime*(pow(cosP1nPhi,2.)+pow(sinP1nPhi,2.)) + - 6.*(pow(cosP1nPhi,2.)-pow(sinP1nPhi,2.)) + * (cosP1nPsi*cosP1nPhi-sinP1nPsi*sinP1nPhi) + - 12.*cosP1nPhi*sinP1nPhi + * (sinP1nPsi*cosP1nPhi+cosP1nPsi*sinP1nPhi); - // Numbers on the side reference term number in paper (cite needed) loosely - /*1*/ vTwo4diffErr = TMath::Power(2.*Avg2*Avg2*Avg2p - 3.*Avg2*Avg4p + 2.*Avg4*Avg2p, 2) - * TMath::Power(sqrtW2sq, 2) * sAvg2sq / (W2*W2); - /*2*/ vTwo4diffErr += 9. * TMath::Power(2.*Avg2*Avg2p - Avg4p, 2) * TMath::Power(sqrtW4sq, 2) - * sAvg4sq / (16. * W4*W4); - /*3*/ vTwo4diffErr += 4. * Avg2*Avg2 * TMath::Power(2.*Avg2*Avg2 - Avg4, 2) * TMath::Power(sqrtW2psq, 2) - * sAvg2psq / (W2p*W2p); - /*4*/ vTwo4diffErr += TMath::Power(2.*Avg2*Avg2 - Avg4, 2) * TMath::Power(sqrtW4psq, 2) * sAvg4psq - / (W4p*W4p); - /*5*/ vTwo4diffErr -= 1.5 * (2.*Avg2*Avg2p - Avg4p) * (2.*Avg2*Avg2*Avg2p - 3.*Avg2*Avg4p + 2.*Avg4*Avg2p) - * W2W4 * Cov24 / (W2*W4); - /*6*/ vTwo4diffErr -= 4. * Avg2 * (2.*Avg2*Avg2 - Avg4) - * (2.*Avg2*Avg2*Avg2p - 3.*Avg2*Avg4p + 2.*Avg4*Avg2p) - * W2W2p * Cov22p / (W2 * W2p); - /*7*/ vTwo4diffErr += 2. * (2.*Avg2*Avg2 - Avg4) - * (2.*Avg2*Avg2*Avg2p - 3.*Avg2*Avg4p + 2.*Avg4*Avg2p) - * W2W4p * Cov24p / (W2 * W4p); - /*8*/ vTwo4diffErr += 3.*Avg2*(2.*Avg2*Avg2 - Avg4)*(2.*Avg2*Avg2p - Avg4p) - * W4W2p * Cov42p / (W4*W2p); - /*9*/ vTwo4diffErr -= 1.5 * (2.*Avg2*Avg2 - Avg4)*(2.*Avg2*Avg2p - Avg4p) - * W4W4p * Cov44p / (W4 * W4p); - /*10*/vTwo4diffErr -= 4.*Avg2*TMath::Power(2.*Avg2*Avg2 - Avg4, 2) - * W2pW4p * Cov2p2np / (W2p * W4p); - /*11*/vTwo4diffErr /= TMath::Power(2.*Avg2*Avg2 - Avg4, 3.5); - vTwo4diffErr = TMath::Sqrt(vTwo4diffErr); - - cumulant4Hist->SetBinError(eta, vTwo4diffErr); - } // End of eta loop - - // Number of events: - Int_t nEv = (Int_t)cumulantsHist->GetBinContent(0,0); - cumulant2Hist->SetBinContent(cumulant2Hist->GetNbinsX() + 1, nEv); - cumulant4Hist->SetBinContent(cumulant4Hist->GetNbinsX() + 1, nEv); + vTwo4diff = - qc4Prime / TMath::Power(-qc4, 0.75); + cumulant4diffHist->SetBinContent(etaBin, vTwo4diff); + + // 2-particle differential flow error calculations + w2pavg2psq = cumulantsHist->GetBinContent(etaBin, kW2avg2sq); + sqrtW2psq = cumulantsHist->GetBinError(etaBin, kW2); + w2W2pavg2avg2p = cumulantsHist->GetBinContent(etaBin, kW2w2pavg2avg2p); + w2W2p = cumulantsHist->GetBinContent(etaBin, kW2w2p); + + cov22p = CovXY(w2W2pavg2avg2p, w2W2p, two*twoPrime, w2, w2p); + stwoPrimesq = VarSQ(w2pavg2psq, twoPrime, w2p, w2pavg2p, sqrtW2psq); + + vTwo2diffErr = twoPrime*twoPrime*TMath::Power(sqrtW2sq, 2)*stwosq/(w2*w2); + vTwo2diffErr += 4.*two*two*TMath::Power(sqrtW2psq, 2)*stwoPrimesq/(w2p*w2p); + vTwo2diffErr -= 4.*two*twoPrime*w2W2p*cov22p/(w2*w2p); + vTwo2diffErr /= (4. * TMath::Power(two, 3)); + vTwo2diffErr = TMath::Sqrt(vTwo2diffErr); + cumulant2diffHist->SetBinError(etaBin, vTwo2diffErr); + + // 4-particle differential flow error calculations + sqrtW4psq = cumulantsHist->GetBinError(etaBin, kW4); + w4pavg4psq = cumulantsHist->GetBinContent(etaBin, kW4avg4sq); + w2pW4p = cumulantsHist->GetBinContent(etaBin, kW2w4); + w2pW4pavg2pavg4p = cumulantsHist->GetBinContent(etaBin, kW2w4avg2avg4); + w2W4p = cumulantsHist->GetBinContent(etaBin, kW2w4p); + w2W4pavg2avg4p = cumulantsHist->GetBinContent(etaBin, kW2w4pavg2avg4p); + w4W2p = cumulantsHist->GetBinContent(etaBin, kW4w2p); + w4W2pavg4avg2p = cumulantsHist->GetBinContent(etaBin, kW4w2pavg4avg2p); + w4W4p = cumulantsHist->GetBinContent(etaBin, kW4w4p); + w4W4pavg4avg4p = cumulantsHist->GetBinContent(etaBin, kW4w4pavg4avg4p); + + sfourPrimesq = VarSQ(w4pavg4psq, fourPrime, w4p, w4pavg4p, sqrtW4psq); + cov24p = CovXY(w2W4pavg2avg4p, w2W4p, two*fourPrime, w2, w4p); + cov42p = CovXY(w4W2pavg4avg2p, w4W2p, four*twoPrime, w4, w2p); + cov44p = CovXY(w4W4pavg4avg4p, w4W4p, four*fourPrime, w4, w4p); + cov2p2np = CovXY(w2pW4pavg2pavg4p, w2pW4p, twoPrime*fourPrime, w2p, w4p); + + // Numbers on the side reference term number in paper (cite needed) loosely + /*1*/ vTwo4diffErr = TMath::Power(2.*two*two*twoPrime - 3.*two*fourPrime + 2.*four*twoPrime, 2) + * TMath::Power(sqrtW2sq, 2) * stwosq / (w2*w2); + /*2*/ vTwo4diffErr += 9. * TMath::Power(2.*two*twoPrime - fourPrime, 2) * TMath::Power(sqrtW4sq, 2) + * sfoursq / (16. * w4*w4); + /*3*/ vTwo4diffErr += 4. * two*two * TMath::Power(2.*two*two - four, 2) * TMath::Power(sqrtW2psq, 2) + * stwoPrimesq / (w2p*w2p); + /*4*/ vTwo4diffErr += TMath::Power(2.*two*two - four, 2) * TMath::Power(sqrtW4psq, 2) * sfourPrimesq + / (w4p*w4p); + /*5*/ vTwo4diffErr -= 1.5 * (2.*two*twoPrime - fourPrime) * (2.*two*two*twoPrime - 3.*two*fourPrime + 2.*four*twoPrime) + * w2W4 * cov24 / (w2*w4); + /*6*/ vTwo4diffErr -= 4. * two * (2.*two*two - four) + * (2.*two*two*twoPrime - 3.*two*fourPrime + 2.*four*twoPrime) + * w2W2p * cov22p / (w2 * w2p); + /*7*/ vTwo4diffErr += 2. * (2.*two*two - four) + * (2.*two*two*twoPrime - 3.*two*fourPrime + 2.*four*twoPrime) + * w2W4p * cov24p / (w2 * w4p); + /*8*/ vTwo4diffErr += 3.*two*(2.*two*two - four)*(2.*two*twoPrime - fourPrime) + * w4W2p * cov42p / (w4*w2p); + /*9*/ vTwo4diffErr -= 1.5 * (2.*two*two - four)*(2.*two*twoPrime - fourPrime) + * w4W4p * cov44p / (w4 * w4p); + /*10*/vTwo4diffErr -= 4.*two*TMath::Power(2.*two*two - four, 2) + * w2pW4p * cov2p2np / (w2p * w4p); + /*11*/vTwo4diffErr /= TMath::Power(2.*two*two - four, 3.5); + vTwo4diffErr = TMath::Sqrt(vTwo4diffErr); + + cumulant4diffHist->SetBinError(etaBin, vTwo4diffErr); + } // End of etaBin loop + // Number of events: + Int_t nEv = (Int_t)cumulantsHist->GetBinContent(0,0); + cumulant2diffHist->SetBinContent(cumulant2diffHist->GetNbinsX() + 1, nEv); + cumulant4diffHist->SetBinContent(cumulant4diffHist->GetNbinsX() + 1, nEv); + } // End of centrality loop } // End of harmonics loop } // End of type loop + + delete cumulantsHist; + } // _____________________________________________________________________ void AliForwardFlowTaskQC::ProcessPrimary() @@ -705,45 +848,29 @@ void AliForwardFlowTaskQC::ProcessPrimary() // can be run on it. // - // Histograms are loaded and reset - TH1D* dNdphi = (TH1D*)fOutputList->FindObject("hdNdphiSEMC"); - TH2D* dNdetadphi = (TH2D*)fOutputList->FindObject("hdNdetadphiSEMC"); - TH1D* dNdphiTrRef = (TH1D*)fOutputList->FindObject("hdNdphiSETrRef"); - TH2D* dNdetadphiTrRef = (TH2D*)fOutputList->FindObject("hdNdetadphiSETrRef"); - - dNdphi->Reset(); - dNdetadphi->Reset(); - dNdphiTrRef->Reset(); - dNdetadphiTrRef->Reset(); - - // Loads AliFMDFlowCommon and fills histograms and runs analysis. - // AOD events also get a TrackRef histogram - AliForwardFlowBase* common = new AliForwardFlowBase(fOutputList); - - if (fAOD) { - if (!common->LoopAODmc(fAOD)) return; - if (!common->LoopAODtrrefHits(fAOD)) return; - // if (!common->LoopMCaddptFlow(fAOD)) return; - // if (!common->LoopMCaddpdgFlow(fAOD)) return; - // if (!common->LoopMCaddetaFlow(fAOD)) return; + if (fFlowUtil->LoopAODtrrefHits(fAOD)) { + // Run analysis on TrackRefs + for (Int_t n = 1; n <= 4; n++) { + if (fv[n]) + CumulantsMethod("TrRef", n); + } } - // Run analysis on MC truth - for (Int_t n = 1; n <= 4; n++) { - if (fv[n]) - CumulantsMethod("MC", n); - } - - // Run analysis on TrackRefs - for (Int_t n = 1; n <= 4; n++) { - if (fv[n]) - CumulantsMethod("TrRef", n); + if (fFlowUtil->LoopAODmc(fAOD, fAddFlow, fAddType, fAddOrder)) { + // Run analysis on MC truth + for (Int_t n = 1; n <= 4; n++) { + if (fv[n]) + CumulantsMethod("MC", n); + } } -} + } //_____________________________________________________________________ -Double_t AliForwardFlowTaskQC::VarSQ(Double_t wxx2, Double_t x, Double_t wx, - Double_t wxx, Double_t sqrtwx2) +Double_t AliForwardFlowTaskQC::VarSQ(Double_t wxx2, + Double_t x, + Double_t wx, + Double_t wxx, + Double_t sqrtwx2) const { // // Small function to compute the variance squared - used by Terminte() @@ -757,20 +884,23 @@ Double_t AliForwardFlowTaskQC::VarSQ(Double_t wxx2, Double_t x, Double_t wx, return sx; } //_____________________________________________________________________ -Double_t AliForwardFlowTaskQC::CovXY(Double_t wxwyxy, Double_t wxwy, - Double_t xy, Double_t wx, Double_t wy) +Double_t AliForwardFlowTaskQC::CovXY(Double_t wxwyxy, + Double_t wxwy, + Double_t xy, + Double_t wx, + Double_t wy) const { // // Small function to compute the covariance between two numbers // - used by Terminate() // - Double_t Cov, denominator, numerator; + Double_t cov, denominator, numerator; denominator = (wxwyxy / wxwy) - xy; numerator = 1 - (wxwy / (wx * wy)); - Cov = denominator / numerator; - return Cov; + cov = denominator / numerator; + return cov; } //_____________________________________________________________________ // diff --git a/PWG2/FORWARD/analysis2/AliForwardFlowTaskQC.h b/PWG2/FORWARD/analysis2/AliForwardFlowTaskQC.h index 70532aca88a..dcee800a2d8 100644 --- a/PWG2/FORWARD/analysis2/AliForwardFlowTaskQC.h +++ b/PWG2/FORWARD/analysis2/AliForwardFlowTaskQC.h @@ -4,7 +4,8 @@ #ifndef ALIFORWARDFLOWTASKQC_H #define ALIFORWARDFLOWTASKQC_H #include "AliAnalysisTaskSE.h" -#include "AliAODEvent.h" +#include "AliForwardFlowUtil.h" +class AliAODEvent; /** * Calculate the flow in the forward regions using the Q cumulants method @@ -95,10 +96,25 @@ public: void SetDoHarmonics(Bool_t v1 = kTRUE, Bool_t v2 = kTRUE, Bool_t v3 = kTRUE, Bool_t v4 = kTRUE) { fv[1] = v1; fv[2] = v2; fv[3] = v3; fv[4] = v4; } + /* + * Set string to add flow to MC truth particles + */ + void AddFlow(TString type = "") { fAddFlow = type; } + /* + * Set which function fAddFlow should use + */ + void AddFlowType(Int_t number = 0) { fAddType = number; } + /* + * Set which order of flow to add + */ + void AddFlowOrder(Int_t order = 2) { fAddOrder = order; } + /** + * Set Z vertex range - Used by flow task + */ + void SetVertexRange(Int_t vertex = 2) { fZvertex = vertex; } /** * @} */ - virtual void SetDebugLevel(Int_t level) { fDebug = level; } protected: /* @@ -123,22 +139,28 @@ protected: * */ Double_t VarSQ(Double_t wxx2, Double_t x, Double_t wx, - Double_t wxx, Double_t sqrtwx2); + Double_t wxx, Double_t sqrtwx2) const ; /** * Caclulate the covariance between x and y - used to finalize * calculations in Terminate() * */ Double_t CovXY(Double_t wxwyxy, Double_t wxwy, Double_t XY, - Double_t wx, Double_t wy); + Double_t wx, Double_t wy) const; - Int_t fDebug; // Debug flag 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 Bool_t fv[5]; // Calculate v_{n} flag - + TString fAddFlow; // Add flow string + Int_t fAddType; // Add flow type # + Int_t fAddOrder; // Add flow order + Int_t fZvertex; // Z vertex range + Double_t fCent; // Centrality + + ClassDef(AliForwardFlowTaskQC, 1); // Analysis task for FMD analysis }; diff --git a/PWG2/FORWARD/analysis2/AliForwardFlowUtil.cxx b/PWG2/FORWARD/analysis2/AliForwardFlowUtil.cxx new file mode 100644 index 00000000000..1425dd07760 --- /dev/null +++ b/PWG2/FORWARD/analysis2/AliForwardFlowUtil.cxx @@ -0,0 +1,363 @@ +// +// 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" + +ClassImp(AliForwardFlowUtil) + +//_____________________________________________________________________ +AliForwardFlowUtil::AliForwardFlowUtil() : + fList(0), + fZvertex(0) {} + // + // Default Constructor + // +//_____________________________________________________________________ +AliForwardFlowUtil::AliForwardFlowUtil(TList* outputList) : + fList(0), + fZvertex(0) +{ + // + // Constructor + // + // Parameters: + // TList: list containing histograms for flow analysis + // + fList = outputList; +} +//_____________________________________________________________________ +Bool_t AliForwardFlowUtil::AODCheck(const AliAODForwardMult* aodfm) const +{ + // + // Function to check that and AOD event meets the cuts + // + // Parameters: + // AliAODForwardMult: forward mult object with trigger and vertex info + // + if (!aodfm->IsTriggerBits(AliAODForwardMult::kInel)) return kFALSE; + if (!aodfm->HasIpZ()) return kFALSE; + if (!aodfm->InRange(-fZvertex,fZvertex)) return kFALSE; + + return kTRUE; +} +//_____________________________________________________________________ +Bool_t AliForwardFlowUtil::LoopAODFMD(const AliAODEvent* aodevent) const +{ + // + // 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(); + TH1D* dNdphi = (TH1D*)fList->FindObject("hdNdphiSE"); + TH2D* dNdetadphi = (TH2D*)fList->FindObject("hdNdetadphiSE"); + 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); + for (Int_t phiBin = 1; phiBin<=hdNdetadphi.GetNbinsY(); phiBin++) { + phi = hdNdetadphi.GetYaxis()->GetBinCenter(phiBin); + weight = hdNdetadphi.GetBinContent(etaBin, phiBin); + dNdetadphi->Fill(eta, phi, weight); + if (eta < 4.) dNdphi->Fill(phi, weight); + mult += weight; + } + } + dNdphi->SetBinContent(0, mult); +// if (aodfmult->HasCentrality()) dNdphi->SetBinContent(dNdphi->GetNbinsX()+1, aodfmult->GetCentrality()); +// else dNdphi->SetBinContent(dNdphi->GetNbinsX()+1, -1); + + dNdphi->SetBinContent(dNdphi->GetNbinsX()+1, GetCentFromMC(aodevent)); + 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* dNdetadphi = (TH2D*)fList->FindObject("hdNdetadphiSESPD"); + dNdetadphi->Reset(); + + 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); + for (Int_t phiBin = 1; phiBin<=hdNdetadphi.GetNbinsY(); phiBin++) { + phi = hdNdetadphi.GetYaxis()->GetBinCenter(phiBin); + weight = hdNdetadphi.GetBinContent(etaBin, phiBin); + dNdetadphi->Fill(eta, phi, weight); + } + } + return kTRUE; +} +//_____________________________________________________________________ +Bool_t AliForwardFlowUtil::LoopAODtrrefHits(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; + // if (!AODCheck(aodfmult)) return kFALSE; + + TH2D hdNdetadphi = aodfmult->GetHistogram(); + TH1D* dNdphi = (TH1D*)fList->FindObject("hdNdphiSETrRef"); + TH2D* dNdetadphi = (TH2D*)fList->FindObject("hdNdetadphiSETrRef"); + 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); + 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(phi, weight); + mult += weight; + } + } + dNdphi->SetBinContent(0, mult); + + 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; + if (addFlow.Length() > 1) { + AliAODMCHeader* header = (AliAODMCHeader*)aodevent->FindListObject(AliAODMCHeader::StdBranchName()); + rp = header->GetReactionPlaneAngle(); + } + + TH2D* dNdetadphi = (TH2D*)fList->FindObject("hdNdetadphiSEMC"); + TH1D* dNdphi = (TH1D*)fList->FindObject("hdNdphiSEMC"); + 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 = 0; + 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; + if (particle->Eta() > -3.4 && particle->Eta() < 5) { + // Add flow if it is in the argument + if (addFlow.Length() > 1) { + if (addFlow.Contains("pt")) + weight += AddptFlow(particle->Pt(), type); + if (addFlow.Contains("pid")) + weight += AddpidFlow(particle->PdgCode(), type); + if (addFlow.Contains("eta")) + weight += AddetaFlow(particle->Eta(), type); + if (addFlow.Contains("flat")) + weight = 0.1*type; + weight *= 2.*TMath::Cos((Double_t)order*(particle->Phi()-rp)); + } + weight += 1; + + dNdphi->Fill(particle->Phi(), weight); + dNdetadphi->Fill(particle->Eta(), particle->Phi(), weight); + mult += weight; + } + } + + dNdphi->SetBinContent(0, mult); + + return kTRUE; +} +//_____________________________________________________________________ +Double_t AliForwardFlowUtil::AddptFlow(Double_t Pt = 0, Int_t type = 0) const +{ + // + // Add pt dependent flow factor + // + Double_t weight = 0; + + if (type == 1) weight = 0.125*Pt; + + if (type == 2) { + weight = 0.2*(Pt/2.0); + if (Pt > 2.0) weight = 0.2; + } + + if (type == 3) { + weight = 0.05*(Pt/2.0); + if (Pt < 2.0) weight = 0.05; + } + + if (type == 4) { + weight = 0.2*(Pt/1.0); + if (Pt > 1.0) weight = 0.2; + } + + if (type == 5) { + weight = 0.05*(Pt/1.0); + if (Pt < 1.0) weight = 0.05; + } + + if (type == 6) { + weight = 0.2*(Pt/3.0); + if (Pt > 3.0) weight = 0.2; + } + + return weight; +} +//_____________________________________________________________________ +Double_t AliForwardFlowUtil::AddpidFlow(Int_t ID = 0, Int_t type = 0) const +{ + // + // Add pid dependent flow factor + // + Double_t weight = 0; + + if (type == 1) { + weight = 0.07; + if (TMath::Abs(ID) == 211) // pion flow + weight = 0.1; + if (TMath::Abs(ID) == 2212) // proton flow + weight = 0.05; + } + + if (type == 2) { + weight = 0.06; + if (TMath::Abs(ID) == 211) // pion flow + weight = 0.1; + if (TMath::Abs(ID) == 2212) // proton flow + weight = 0.08; + } + + if (type == 3) { + weight = 0.05; + if (TMath::Abs(ID) == 211) // pion flow + weight = 0.1; + if (TMath::Abs(ID) == 2212) // proton flow + weight = 0.07; + } + + if (type == 4) { + weight = 0.07; + if (TMath::Abs(ID) == 211) // pion flow + weight = 0.1; + if (TMath::Abs(ID) == 2212) // proton flow + weight = 0.085; + } + + return weight; +} +//_____________________________________________________________________ +Double_t AliForwardFlowUtil::AddetaFlow(Double_t eta = 0, Int_t type = 0) const +{ + // + // Add eta dependent flow factor + // + Double_t weight = 0; + + if (type == 1) weight = 0.03 + 0.07 * (1 - TMath::Abs(eta) / 6.); + + if (type == 2) weight = 0.07 * (1 - TMath::Abs(eta) / 6.); + + if (type == 3) weight = 0.07 * (1 - TMath::Abs(eta) / 8.); + + if (type == 4) weight = 0.07 * (1 - TMath::Abs(eta) / 10.); + + if (type == 5) weight = 0.07 * (1 - TMath::Abs(eta) / 12.); + + if (type == 6) weight = 0.07 * (1 - TMath::Abs(eta) / 14.); + + if (type == 7) weight = 0.07 * (1 - TMath::Abs(eta) / 16.); + + 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(); + + if ( 0.00 <= b && b < 3.50) + cent = 2.5; + if ( 3.50 <= b && b < 4.95) + cent = 7.5; + if ( 4.95 <= b && b < 6.98) + cent = 15.; + if ( 6.98 <= b && b < 8.55) + cent = 25.; + if ( 8.55 <= b && b < 9.88) + cent = 35.; + if ( 9.88 <= b && b < 11.04) + cent = 45.; + if (11.04 <= b && b < 12.09) + cent = 55.; + if (12.09 <= b && b < 13.06) + cent = 65.; + if (13.06 <= b && b < 13.97) + cent = 75.; + if (13.97 <= b && b < 19.61) + cent = 90.; + + return cent; +} +//_____________________________________________________________________ +// +// +// EOF + diff --git a/PWG2/FORWARD/analysis2/AliForwardFlowUtil.h b/PWG2/FORWARD/analysis2/AliForwardFlowUtil.h new file mode 100644 index 00000000000..fecd7f5a898 --- /dev/null +++ b/PWG2/FORWARD/analysis2/AliForwardFlowUtil.h @@ -0,0 +1,102 @@ +// +// 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 +#include "TNamed.h" +class AliAODForwardMult; +class AliAODEvent; +class TList; + +/** + * + * Class used to handle the input from AODs and put it into histograms + * the Forward Flow tasks can run on. + * + * @ingroup pwg2_forward_tasks + */ +class AliForwardFlowUtil : public TNamed +{ +public: + /** + * Constructor + */ + AliForwardFlowUtil(); + /* + * Constructor + * + * @param fList list of histograms for flow analysis + */ + AliForwardFlowUtil(TList* fList); + /* + * Copy constructor + * + * @param o Object to copy from + */ + AliForwardFlowUtil(const AliForwardFlowUtil& o) : TNamed(), + fList(o.fList), + fZvertex(o.fZvertex) {} + /** + * Assignment operator + * + * @param o Object to assign from + * + * @return Reference to this object + */ + AliForwardFlowUtil& operator=(const AliForwardFlowUtil&) { return *this; } + /** + * Check that AOD event meet trigger requirements + */ + Bool_t AODCheck(const AliAODForwardMult* aodfm) const; + /** + * Loop over AliAODForwardMult object and fill flow histograms + */ + Bool_t LoopAODFMD(const AliAODEvent* AODevent) const; + /* + * Loop over AliAODForwardCentral object and fill flow histograms + */ + Bool_t LoopAODSPD(const AliAODEvent* AODevent) const; + /** + * Loop over AliAODForwardMult object and fill flow histograms from + * track refs + */ + Bool_t LoopAODtrrefHits(const AliAODEvent* AODevent) const; + /** + * Loop over AliAODMCParticle branch object and fill flow histograms + * add flow if arguments are set + */ + Bool_t LoopAODmc(const AliAODEvent* AODevent, TString addFlow, Int_t type, Int_t order) const; + /** + * Set Z vertex range - Used by flow task + */ + void SetVertexRange(Int_t vertex = 2) { fZvertex = vertex; } + +protected: + /** + * Add pt dependent flow factor + */ + Double_t AddptFlow(Double_t Pt, Int_t type) const; + /** + * Add pid dependent flow factor + */ + Double_t AddpidFlow(Int_t ID, Int_t type) const; + /** + * Add eta dependent flow factor + */ + Double_t AddetaFlow(Double_t Eta, Int_t type) const; + /** + * Get centrality form MC impact parameter + */ + Double_t GetCentFromMC(const AliAODEvent* AODevent) const; + + TList* fList; // List of flow histograms + Int_t fZvertex; // Z vertex range + + ClassDef(AliForwardFlowUtil, 1); +}; + +#endif +// Local Variables: +// mode: C++ +// End: diff --git a/PWG2/FORWARD/analysis2/MakeFlow.C b/PWG2/FORWARD/analysis2/MakeFlow.C index 1f401fbe6cd..a85d4132f1a 100644 --- a/PWG2/FORWARD/analysis2/MakeFlow.C +++ b/PWG2/FORWARD/analysis2/MakeFlow.C @@ -1,43 +1,62 @@ /** - * Script to analyse AOD input for flow + * Script to analyse AOD input for flow + * + * Takes either a single (AOD) .root file as input or a .txt + * The .txt file is expected to contain the path to the files + * from the current directory or the absolute path. * * @par Inputs: - * - + * * * @par Outputs: * - * */ -void MakeFlow(TString path = "", - Bool_t recursive = true, - Int_t nevents = 100, +void MakeFlow(TString data = "", + Int_t nevents = 0, TString type = "", - Int_t etabins) + Int_t etabins = 40, + Int_t zVertex = 2, + TString addFlow = "", + Int_t addFType = 0, + Int_t addFOrder = 0) { - Bool_t proof = false; + Bool_t proof = kFALSE; - // --- Load libs ------------------------------------------------- + // --- Load libs --------------------------------------------------- gROOT->Macro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C"); // --- Check for proof mode, and possibly upload pars -------------- if (proof> 0) { gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadPars.C"); if (!LoadPars(proof)) { - Error("MakeAOD", "Failed to load PARs"); + AliError("MakeFlow", "Failed to load PARs"); return; } } - // --- Add to chain either ESD or AOD ---------------------------- - gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/MakeChain.C"); - TChain* chain = MakeChain("AOD", path, recursive); - // If 0 or less events is select, choose all - if (nevents <= 0) nevents = chain->GetEntries(); + // --- Add to chain either AOD ------------------------------------ + if (data.Length() <= 1) { + AliError("You didn't add a data file"); + return; + } + TChain* chain = new TChain("aodTree"); + + if (data.Contains(".txt")) + MakeChain(data, chain); - // --- Initiate the event handlers ------------------------------- + if (data.Contains(".root")) { + if (!TFile::Open(data.Data())) { + AliError(Form("AOD file %s not found", data.Data())); + return; + } + chain->Add(data.Data()); + } + + // --- Initiate the event handlers -------------------------------- AliAnalysisManager *mgr = new AliAnalysisManager("Forward Flow", - "Flow in forward region"); - mgr->SetUseProgressBar(kTRUE); + "Flow in the forward region"); + mgr->SetUseProgressBar(kTRUE, 10); // --- AOD input handler ------------------------------------------- AliAODInputHandler *aodInputHandler = new AliAODInputHandler(); @@ -50,7 +69,7 @@ void MakeFlow(TString path = "", // --- Add the tasks --------------------------------------------- gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/AddTaskForwardFlow.C"); - AddTaskForwardFlow(type.Data(), etabins); + AddTaskForwardFlow(type, etabins, zVertex, addFlow, addFType, addFOrder); // --- Run the analysis -------------------------------------------- TStopwatch t; @@ -66,6 +85,32 @@ void MakeFlow(TString path = "", 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; + while(in.good()) + { + in >> line; + + if (line.Length() == 0) + continue; + + if (TFile::Open(line)) + chain->Add(line); + } + + in.close(); + + return; +} // // EOF // diff --git a/PWG2/PWG2forward2LinkDef.h b/PWG2/PWG2forward2LinkDef.h index 36f5d47cf3e..f9ed7e93304 100644 --- a/PWG2/PWG2forward2LinkDef.h +++ b/PWG2/PWG2forward2LinkDef.h @@ -71,7 +71,7 @@ #pragma link C++ class AliCentralCorrSecondaryMap+; #pragma link C++ class AliCentralCorrAcceptance+; #pragma link C++ class AliCentraldNdetaTask+; -#pragma link C++ class AliForwardFlowBase+; +#pragma link C++ class AliForwardFlowUtil+; #pragma link C++ class AliForwardFlowTaskQC+; #else -- 2.43.0