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
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
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
/**
* @file AddTaskForwardFlow.C
- * @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk>
- * @date Wed Mar 23 12:14:17 2011
+ * @author Alexander Hansen alexander.hansen@cern.ch
+ * @date Wed Sep 07 12:14:17 2011
*
* @brief
*
* @ingroup pwglf_forward_flow
*/
void AddTaskForwardFlow(TString type = "",
- Int_t etabins = 48,
Bool_t mc = kFALSE,
TString addFlow = "",
Int_t addFType = 0,
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_cast<AliForwardMCFlowTaskQC*>task;
+ 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;
}
--- /dev/null
+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;
+}
--- /dev/null
+//
+// Class that contains results from FMD eventplane calculations
+//
+#include "AliAODForwardEP.h"
+#include <TBrowser.h>
+#include <iostream>
+#include <TMath.h>
+#include <TObjString.h>
+#include <TObjArray.h>
+#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 <<std::endl;
+ std::cout << "FMD1+2 EP: \t" << fEpA <<std::endl;
+ std::cout << "FMD3 EP: \t" << fEpC <<std::endl;
+ std::cout << "Random Subep 1: \t" << fEp1 <<std::endl;
+ std::cout << "Random Subep 2: \t" << fEp2 <<std::endl;
+}
+
+//____________________________________________________________________
+//
+// EOF
+//
--- /dev/null
+
+//
+// Class that contains results from FMD eventplane calculations
+//
+#ifndef ALIAODFORWARDEP_H
+#define ALIAODFORWARDEP_H
+/**
+ * @file AliAODForwardMult.h
+ * @author Alexander Hansen
+ * @date Tue Feb 14 2012
+ *
+ * @brief
+ *
+ *
+ * @ingroup pwglf_forward
+ */
+#include <TObject.h>
+#include <TH1D.h>
+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:
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");
--- /dev/null
+// This task finds the eventplane
+// using the FMD
+//
+#include "AliFMDEventPlaneFinder.h"
+#include <TList.h>
+#include <TMath.h>
+#include "AliLog.h"
+#include <TH2D.h>
+#include <TH3D.h>
+#include "TROOT.h"
+#include <iostream>
+#include <iomanip>
+#include <TString.h>
+#include <TFile.h>
+#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; q<nr; q++) {
+ Char_t r = (q == 0 ? 'I' : 'O');
+ CalcQVectors(hists->Get(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
+//
+
+
+
--- /dev/null
+// 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 <TNamed.h>
+#include <TVector2.h>
+#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:
+
#include <TList.h>
#include <TROOT.h>
#include <iostream>
+#include "AliGenHijingEventHeader.h"
+#include <TF1.h>
+#include <TGraph.h>
//____________________________________________________________________
AliFMDMCTrackDensity::AliFMDMCTrackDensity()
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;
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();
return 0;
}
+// #define USE_FLOW_WEIGHTS 1
//____________________________________________________________________
Bool_t
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<AliGenHijingEventHeader*>
+ (event.GenEventHeader());
+ Double_t rp = (hd ? hd->ReactionPlaneAngle() : 0.);
+ Double_t b = (hd ? hd->ImpactParameter() : -1 );
+#endif
+
AliStack* stack = const_cast<AliMCEvent&>(event).Stack();
Int_t nTracks = stack->GetNtrack();//event.GetNumberOfTracks();
Int_t nPrim = stack->GetNprimary();//event.GetNumberOfPrimary();
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
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
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
*
* @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
#include <TList.h>
#include <iostream>
#include <TMath.h>
-#include "THnSparse.h"
#include "TH3D.h"
#include "TProfile2D.h"
#include "AliLog.h"
#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
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
// 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
// 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<VertexBin*>(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<VertexBin*>(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<AliAODEvent*>(InputEvent());
- if (!fAOD) return;
+ if (!fAOD) return kFALSE;
- // Fill histograms
- if (!fFlowUtil->LoopAODFMD(fAOD)) return;
- if (!fFlowUtil->LoopAODSPD(fAOD)) return;
+ AliAODForwardMult* aodfmult = static_cast<AliAODForwardMult*>(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<AliAODCentralMult*>(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<VertexBin*>(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<VertexBin*>(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<TList*> (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<TProfile2D*>(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<VertexBin*>(nextFMD()))) {
+ bin->CumulantsTerminate(fSumList, fOutputList);
+ }
+ TIter nextSPD(&fBinsSPD);
+ while ((bin = static_cast<VertexBin*>(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;
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;
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.);
-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;
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;
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.);
+ 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;
- 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;
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);
- }
- }
-
-}
//_____________________________________________________________________
//
//
#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
* @ingroup pwglf_forward_tasks_flow
* @ingroup pwglf_forward_flow
*
- * @todo Add centrality stuff
*
*/
class AliForwardFlowTaskQC : public AliAnalysisTaskSE
*/
/**
* Create output objects
- *
*/
virtual void UserCreateOutputObjects();
/**
* Initialize the task
- *
*/
- virtual void Init() {}
+ virtual void Init() {}
/**
* Process each event
*
*/
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
*
* @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
+++ /dev/null
-//
-// 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 <iostream>
-#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<AliAODForwardMult*>(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<AliAODCentralMult*>(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<AliAODForwardMult*>(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<AliAODCentralMult*>(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<AliAODMCHeader*>(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
-
+++ /dev/null
-//
-// 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:
--- /dev/null
+//
+// 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<VertexBin*>(nextFMDTR()))) {
+ bin->AddOutput(fSumList);
+ }
+ TIter nextSPDTR(&fBinsSPDTR);
+ while ((bin = static_cast<VertexBin*>(nextSPDTR()))) {
+ bin->AddOutput(fSumList);
+ }
+ TIter nextMC(&fBinsMC);
+ while ((bin = static_cast<VertexBin*>(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<AliAODForwardMult*>(fAOD->FindListObject("ForwardMC"));
+ if (!aodfmult) return kFALSE;
+ TH2D fmdTRdNdetadphi = aodfmult->GetHistogram();
+
+ AliAODCentralMult* aodcmult = static_cast<AliAODCentralMult*>(fAOD->FindListObject("CentralClustersMC"));
+ if (!aodcmult) return kFALSE;
+ TH2D spdTRdNdetadphi = aodcmult->GetHistogram();
+
+ TIter nextFMDTR(&fBinsFMDTR);
+ VertexBin* bin = 0;
+ while ((bin = static_cast<VertexBin*>(nextFMDTR()))) {
+ if (bin->CheckVertex(fZvertex)) {
+ if (!bin->FillHists(fmdTRdNdetadphi)) return kFALSE;
+ bin->CumulantsAccumulate(fCent);
+ }
+ }
+
+ TIter nextSPDTR(&fBinsSPDTR);
+ while ((bin = static_cast<VertexBin*>(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<VertexBin*>(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<VertexBin*>(nextFMDTR()))) {
+ bin->CumulantsTerminate(fSumList, fOutputList);
+ }
+ TIter nextSPDTR(&fBinsSPDTR);
+ while ((bin = static_cast<VertexBin*>(nextSPDTR()))) {
+ bin->CumulantsTerminate(fSumList, fOutputList);
+ }
+ TIter nextMC(&fBinsMC);
+ while ((bin = static_cast<VertexBin*>(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<AliAODMCHeader*>(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
+
--- /dev/null
+//
+// 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 <TH2D.h>
+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:
fESDFMD(),
fHistos(),
fAODFMD(),
+ fAODEP(),
fMCESDFMD(),
fMCHistos(),
fMCAODFMD(),
fDensityCalculator(),
fCorrections(),
fHistCollector(),
+ fEventPlaneFinder(),
fList(0)
{
//
fESDFMD(),
fHistos(),
fAODFMD(kFALSE),
+ fAODEP(kFALSE),
fMCESDFMD(),
fMCHistos(),
fMCAODFMD(kTRUE),
fDensityCalculator("density"),
fCorrections("corrections"),
fHistCollector("collector"),
+ fEventPlaneFinder("eventplane"),
fList(0)
{
//
fESDFMD(o.fESDFMD),
fHistos(o.fHistos),
fAODFMD(o.fAODFMD),
+ fAODEP(o.fAODEP),
fMCESDFMD(o.fMCESDFMD),
fMCHistos(o.fMCHistos),
fMCAODFMD(o.fMCAODFMD),
fDensityCalculator(o.fDensityCalculator),
fCorrections(o.fCorrections),
fHistCollector(o.fHistCollector),
+ fEventPlaneFinder(o.fEventPlaneFinder),
fList(o.fList)
{
//
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;
fDensityCalculator.SetDebug(dbg);
fCorrections.SetDebug(dbg);
fHistCollector.SetDebug(dbg);
+ fEventPlaneFinder.SetDebug(dbg);
}
//____________________________________________________________________
void
fHistos.Init(*pe);
fAODFMD.Init(*pe);
+ fAODEP.Init(*pe);
fMCHistos.Init(*pe);
fMCAODFMD.Init(*pe);
fRingSums.Init(*pe);
fDensityCalculator.Init(*pe);
fCorrections.Init(*pe);
fHistCollector.Init(*pv,*pe);
+ fEventPlaneFinder.Init(*pe);
this->Print();
}
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");
fDensityCalculator.DefineOutput(fList);
fCorrections.DefineOutput(fList);
fHistCollector.DefineOutput(fList);
+ fEventPlaneFinder.DefineOutput(fList);
PostData(1, fList);
}
fHistos.Clear();
fESDFMD.Clear();
fAODFMD.Clear();
+ fAODEP.Clear();
fMCHistos.Clear();
fMCESDFMD.Clear();
fMCAODFMD.Clear();
}
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");
#include "AliFMDMCCorrector.h"
#include "AliFMDHistCollector.h"
#include "AliAODForwardMult.h"
+#include "AliAODForwardEP.h"
#include "AliFMDEnergyFitter.h"
+#include "AliFMDEventPlaneFinder.h"
#include <AliESDFMD.h>
class AliESDEvent;
class TH2D;
* @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; }
/**
* @}
*/
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
AliFMDMCDensityCalculator fDensityCalculator; // Algorithm
AliFMDMCCorrector fCorrections; // Algorithm
AliFMDHistCollector fHistCollector; // Algorithm
+ AliFMDEventPlaneFinder fEventPlaneFinder; // Algorithm
TList* fList; // Output list
#include "AliFMDDensityCalculator.h"
#include "AliFMDCorrector.h"
#include "AliFMDHistCollector.h"
+#include "AliFMDEventPlaneFinder.h"
#include "AliESDEvent.h"
#include <TROOT.h>
#include <TAxis.h>
fFirstEvent = false;
+ GetEventPlaneFinder().SetRunNumber(esd->GetRunNumber());
InitializeSubs();
}
return esd;
GetDensityCalculator().Print(option);
GetCorrections() .Print(option);
GetHistCollector() .Print(option);
+ GetEventPlaneFinder() .Print(option);
gROOT->DecreaseDirLevel();
}
class AliFMDCorrector;
class AliFMDHistCollector;
class AliForwardCorrectionManager;
+class AliFMDEventPlaneFinder;
class AliESDEvent;
class TH2D;
class TList;
* @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;
+
/**
* @}
*/
fESDFMD(),
fHistos(),
fAODFMD(),
+ fAODEP(),
fRingSums(),
fEventInspector(),
fSharingFilter(),
fDensityCalculator(),
fCorrections(),
fHistCollector(),
+ fEventPlaneFinder(),
fList(0)
{
//
fESDFMD(),
fHistos(),
fAODFMD(false),
+ fAODEP(false),
fRingSums(),
fEventInspector("event"),
fSharingFilter("sharing"),
fDensityCalculator("density"),
fCorrections("corrections"),
fHistCollector("collector"),
+ fEventPlaneFinder("eventplane"),
fList(0)
{
//
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)
{
//
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;
fDensityCalculator.SetDebug(dbg);
fCorrections.SetDebug(dbg);
fHistCollector.SetDebug(dbg);
+ fEventPlaneFinder.SetDebug(dbg);
}
//____________________________________________________________________
fHistos.Init(*pe);
fAODFMD.Init(*pe);
+ fAODEP.Init(*pe);
fRingSums.Init(*pe);
fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
fDensityCalculator.Init(*pe);
fCorrections.Init(*pe);
fHistCollector.Init(*pv,*pe);
+ fEventPlaneFinder.Init(*pe);
this->Print();
}
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);
}
fHistos.Clear();
fESDFMD.Clear();
fAODFMD.Clear();
+ fAODEP.Clear();
Bool_t lowFlux = kFALSE;
UInt_t triggers = 0;
fAODFMD.SetCentrality(cent);
fAODFMD.SetNClusters(nClusters);
MarkEventForStore();
-
+
if (found & AliFMDEventInspector::kNoSPD) return;
if (found & AliFMDEventInspector::kNoFMD) return;
if (found & AliFMDEventInspector::kNoVertex) return;
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)) {
#include "AliFMDCorrector.h"
#include "AliFMDHistCollector.h"
#include "AliAODForwardMult.h"
+#include "AliAODForwardEP.h"
#include "AliFMDEnergyFitter.h"
+#include "AliFMDEventPlaneFinder.h"
#include <AliESDFMD.h>
class AliESDEvent;
class TH2D;
* @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
*
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
AliFMDDensityCalculator fDensityCalculator; // Algorithm
AliFMDCorrector fCorrections; // Algorithm
AliFMDHistCollector fHistCollector; // Algorithm
+ AliFMDEventPlaneFinder fEventPlaneFinder; // Algorithm
TList* fList; // Output list
#include <TList.h>
#include <TROOT.h>
#include <iostream>
+#include "AliGenHijingEventHeader.h"
+#include <TF1.h>
+#include <TGraph.h>
//____________________________________________________________________
AliSPDMCTrackDensity::AliSPDMCTrackDensity()
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;
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);
}
return 0;
}
+// #define USE_FLOW_WEIGHTS
//____________________________________________________________________
Bool_t
AliSPDMCTrackDensity::Calculate(const AliMCEvent& event,
// True on succes, false otherwise
//
+#ifdef USE_FLOW_WEIGHTS
+ AliGenHijingEventHeader* hd = dynamic_cast<AliGenHijingEventHeader*>
+ (event.GenEventHeader());
+ Double_t rp = (hd ? hd->ReactionPlaneAngle() : 0.);
+ Double_t b = (hd ? hd->ImpactParameter() : -1 );
+#endif
+
AliStack* stack = const_cast<AliMCEvent&>(event).Stack();
Int_t nTracks = stack->GetNtrack();
Int_t nPrim = stack->GetNtrack();
// 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
const AliMCParticle* mother,
Int_t refNo,
Double_t vz,
- TH2D& output) const;
+ TH2D& output,
+ Double_t w) const;
/**
* Get ultimate mother of a track
*
*/
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
// 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;
*
* @brief
*
- * @ingroup pwglf_forward_scripts_makers
+ * @ingroup pwg2_forward_scripts_makers
*
*/
/**
* @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 ------------------------------------------
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",
// --- 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;
}
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
//
+#include "TrainSetup.C"
+
/**
* @file MakeAODTrain.C
* @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk>
* @ingroup pwglf_forward_trains
*/
//====================================================================
-#include "TrainSetup.C"
/**
* Analysis train to make Forward and Central multiplicity
// --- 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");
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");
+
}
//__________________________________________________________________
/**
// --- 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);
+ }*/
}
//__________________________________________________________________
/**
--- /dev/null
+#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
+//
#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+;
#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+;
#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+;
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff