analysis from AODs. Added script to add to manager.
Updated Run.sh and Pass2.C to use this new task.
Added script DrawdNdeta.C to draw the result of the
task. Renamed old MakeESDChain.C MakeChain.C
FORWARD/analysis2/AliForwardMultiplicityTask.cxx
FORWARD/analysis2/AliForwardMCMultiplicityTask.cxx
FORWARD/analysis2/AliForwardMCCorrectionsTask.cxx
+ FORWARD/analysis2/AliForwarddNdetaTask.cxx
FORWARD/analysis2/AliForwardUtil.cxx
)
set ( EINCLUDE ANALYSIS )
-set ( EXPORT FORWARD/analysis2/AliAODForwardMult.h )
+set ( EXPORT FORWARD/analysis2/AliAODForwardMult.h
+ FORWARD/analysis2/AliForwardUtil.h )
set ( DHDR PWG2forward2LinkDef.h)
+# --------------------------------------------------------------------
+# Extra targets
+#
+add_custom_command( OUTPUT FORWARD/doc/doc.pdf
+ COMMAND pdflatex doc
+ COMMAND pdflatex doc
+ COMMAND pdflatex doc
+ COMMAND mkdir -p ${CMAKE_CURRENT_BINARY_DIR}/FORWARD/doc
+ COMMAND mv doc.pdf ${CMAKE_CURRENT_BINARY_DIR}/FORWARD/doc/
+ WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/FORWARD/doc
+ MAIN_DEPENDENCY FORWARD/doc/doc.tex )
+
+add_custom_command( OUTPUT FORWARD/doc/html/index.html
+ COMMAND mkdir -p ${CMAKE_CURRENT_BINARY_DIR}/FORWARD/doc
+ COMMAND sed -e 's,\\\(INPUT\\|EXAMPLE_PATH\\\) *= ../analysis2,\\1 = ${CMAKE_CURRENT_SOURCE_DIR}/FORWARD/analysis2,' -e 's,OUTPUT_DIRECTORY *=.*,OUTPUT_DIRECTORY = FORWARD/doc/,' < ${CMAKE_CURRENT_SOURCE_DIR}/FORWARD/doc/Doxyfile > ${CMAKE_CURRENT_BINARY_DIR}/FORWARD/doc/Doxyfile
+ COMMAND doxygen ${CMAKE_CURRENT_BINARY_DIR}/FORWARD/doc/Doxyfile
+ DEPENDS ${SRCS} ${HDRS} )
+
+add_custom_target( PWG2forward-doc DEPENDS FORWARD/doc/doc.pdf )
+add_custom_target( PWG2forward-doxy DEPENDS FORWARD/doc/html/index.html )
+
# --------------------------------------------------------------------
# Extra installation targets
#
FORWARD/analysis2/AddTaskCopyHeader.C
FORWARD/analysis2/AddTaskFMD.C
FORWARD/analysis2/AddTaskFMDELoss.C
- FORWARD/analysis2/AnalyseAOD.C
+ FORWARD/analysis2/AddTaskForwarddNdeta.C
+ FORWARD/analysis2/DrawdNdeta.C
FORWARD/analysis2/MakeAOD.C
+ FORWARD/analysis2/MakedNdeta.C
FORWARD/analysis2/MakeELossFits.C
FORWARD/analysis2/OtherData.C
FORWARD/analysis2/Pass1.C
--- /dev/null
+/**
+ * @file AddTaskForwarddNdeta.C
+ * @author Christian Holm Christensen <cholm@nbi.dk>
+ * @date Fri Jan 28 10:22:26 2011
+ *
+ * @brief Script to add a multiplicity task for the central
+ * @f$\eta@f$ region
+ *
+ *
+ */
+AliAnalysisTask*
+AddTaskForwarddNdeta(const char* trig="INEL", Double_t vzMin=-10, Double_t vzMax=10)
+{
+ // analysis manager
+ AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
+
+ // Make our object. 2nd argumenent is absolute max Eta
+ // 3rd argument is absolute max Vz
+ AliForwarddNdetaTask* task = new AliForwarddNdetaTask("Forward");
+ task->SetVertexRange(vzMin, vzMax);
+ task->SetTriggerMask(trig);
+ mgr->AddTask(task);
+
+ // create containers for input/output
+ AliAnalysisDataContainer *sums =
+ mgr->CreateContainer("ForwardSums", TList::Class(),
+ AliAnalysisManager::kOutputContainer,
+ AliAnalysisManager::GetCommonFileName());
+ AliAnalysisDataContainer *output =
+ mgr->CreateContainer("ForwardResults", TList::Class(),
+ AliAnalysisManager::kParamContainer,
+ AliAnalysisManager::GetCommonFileName());
+
+ // connect input/output
+ mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
+ mgr->ConnectOutput(task, 1, sums);
+ mgr->ConnectOutput(task, 2, output);
+
+ return task;
+}
+
+
+//________________________________________________________________________
+//
+// EOF
+//
fTriggers = 0;
fIpZ = fgkInvalidIpZ;
}
+//____________________________________________________________________
+void
+AliAODForwardMult::SetSNN(UShort_t snn)
+{
+ // set the center of mass energy per nucleon pair (GeV).
+ // This is stored in bin (0,0) of the histogram
+ //
+ // Parameters:
+ // sNN Center of mass energy per nuclean
+ fHist.SetBinContent(0,0,snn);
+}
+//____________________________________________________________________
+void
+AliAODForwardMult::SetSystem(UShort_t sys)
+{
+ // set the center of mass energy per nucleon pair (GeV).
+ // This is stored in bin (N+1,0) of the histogram
+ //
+ // Parameters:
+ // sys Collision system number
+ fHist.SetBinContent(fHist.GetNbinsX()+1,0,sys);
+}
+
//____________________________________________________________________
Bool_t
AliAODForwardMult::HasIpZ() const
return TMath::Abs(fIpZ - fgkInvalidIpZ) > 1;
}
+//____________________________________________________________________
+UShort_t
+AliAODForwardMult::GetSNN() const
+{
+ // set the center of mass energy per nucleon pair (GeV).
+ // This is stored in bin (0,0) of the histogram
+ //
+ // Parameters:
+ // sNN Center of mass energy per nuclean
+ return fHist.GetBinContent(0,0);
+}
+
+//____________________________________________________________________
+UShort_t
+AliAODForwardMult::GetSystem() const
+{
+ // set the center of mass energy per nucleon pair (GeV).
+ // This is stored in bin (N+1,0) of the histogram
+ //
+ // Parameters:
+ // sNN Center of mass energy per nuclean
+ return fHist.GetBinContent(fHist.GetNbinsX()+1,0);
+}
+
//____________________________________________________________________
void
AliAODForwardMult::Browse(TBrowser* b)
* @param ipZ Interaction point z coordinate
*/
void SetIpZ(Float_t ipZ) { fIpZ = ipZ; } // Set Ip's Z coordinate
+ /**
+ * Set the center of mass energy per nucleon-pair. This is stored
+ * in the (0,0) of the histogram
+ *
+ * @param sNN Center of mass energy per nucleon pair (GeV)
+ */
+ void SetSNN(UShort_t sNN);
+ /**
+ * Get the collision system number
+ * - 0: Unknown
+ * - 1: pp
+ * - 2: PbPb
+ *
+ * @param sys Collision system number
+ */
+ void SetSystem(UShort_t sys);
/**
* Set the z coordinate of the interaction point
*
* @return True if we have a valid interaction point z coordinate
*/
Bool_t HasIpZ() const;
+ /**
+ * Get the center of mass energy per nucleon pair (GeV)
+ *
+ * @return Center of mass energy per nucleon pair (GeV)
+ */
+ UShort_t GetSNN() const;
+ /**
+ * Get the collision system number
+ * - 0: Unknown
+ * - 1: pp
+ * - 2: PbPb
+ *
+ * @return Collision system number
+ */
+ UShort_t GetSystem() const;
/**
* Check if the z coordinate of the interaction point is within the
* given limits. Note that the convention used corresponds to the
// Set trigger bits, and mark this event for storage
fAODFMD.SetTriggerBits(triggers);
+ fAODFMD.SetSNN(fEventInspector.GetEnergy());
+ fAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
MarkEventForStore();
if (found & AliFMDEventInspector::kNoSPD) return;
--- /dev/null
+//====================================================================
+#include "AliForwarddNdetaTask.h"
+#include <TMath.h>
+#include <TH2D.h>
+#include <TH1D.h>
+#include <THStack.h>
+#include <TList.h>
+#include <AliAnalysisManager.h>
+#include <AliAODEvent.h>
+#include <AliAODHandler.h>
+#include <AliAODInputHandler.h>
+#include "AliForwardUtil.h"
+
+//____________________________________________________________________
+AliForwarddNdetaTask::AliForwarddNdetaTask()
+ : AliAnalysisTaskSE(),
+ fSumForward(0), // Sum of histograms
+ fSumForwardMC(0), // Sum of MC histograms (if any)
+ fSumPrimary(0), // Sum of primary histograms
+ fSumCentral(0), // Sum of central histograms
+ fCentral(0), // Cache of central histogram
+ fPrimary(0), // Cache of primary histogram
+ fSums(0), // Container of sums
+ fOutput(0), // Container of outputs
+ fTriggers(0), // Histogram of triggers
+ fVtxMin(0), // Minimum v_z
+ fVtxMax(0), // Maximum v_z
+ fTriggerMask(0), // Trigger mask
+ fRebin(0), // Rebinning factor
+ fCutEdges(false),
+ fSNNString(0),
+ fSysString(0)
+{}
+
+//____________________________________________________________________
+AliForwarddNdetaTask::AliForwarddNdetaTask(const char* /* name */)
+ : AliAnalysisTaskSE("Forward"),
+ fSumForward(0), // Sum of histograms
+ fSumForwardMC(0), // Sum of MC histograms (if any)
+ fSumPrimary(0), // Sum of primary histograms
+ fSumCentral(0), // Sum of central histograms
+ fCentral(0), // Cache of central histogram
+ fPrimary(0), // Cache of primary histogram
+ fSums(0), // Container of sums
+ fOutput(0), // Container of outputs
+ fTriggers(0), // Histogram of triggers
+ fVtxMin(-10), // Minimum v_z
+ fVtxMax(10), // Maximum v_z
+ fTriggerMask(AliAODForwardMult::kInel),
+ fRebin(5), // Rebinning factor
+ fCutEdges(false),
+ fSNNString(0),
+ fSysString(0)
+{
+ // Output slot #1 writes into a TH1 container
+ DefineOutput(1, TList::Class());
+ DefineOutput(2, TList::Class());
+}
+
+//____________________________________________________________________
+AliForwarddNdetaTask::AliForwarddNdetaTask(const AliForwarddNdetaTask& o)
+ : AliAnalysisTaskSE(o),
+ fSumForward(o.fSumForward), // TH2D* - Sum of histograms
+ fSumForwardMC(o.fSumForwardMC), // TH2D* - Sum of MC histograms (if any)
+ fSumPrimary(o.fSumPrimary), // TH2D* - Sum of primary histograms
+ fSumCentral(o.fSumCentral), // TH2D* - Sum of central histograms
+ fCentral(o.fCentral), //! Cache of central histogram
+ fPrimary(o.fPrimary), //! Cache of primary histogram
+ fSums(o.fSums), // TList* - Container of sums
+ fOutput(o.fOutput), // TList* - Container of outputs
+ fTriggers(o.fTriggers), // TH1D* - Histogram of triggers
+ fVtxMin(o.fVtxMin), // Double_t - Minimum v_z
+ fVtxMax(o.fVtxMax), // Double_t - Maximum v_z
+ fTriggerMask(o.fTriggerMask), // Int_t - Trigger mask
+ fRebin(o.fRebin), // Int_t - Rebinning factor
+ fCutEdges(o.fCutEdges), // Bool_t - Whether to cut edges when rebinning
+ fSNNString(o.fSNNString), // TNamed* -
+ fSysString(o.fSysString) // TNamed* -
+{}
+
+//____________________________________________________________________
+AliForwarddNdetaTask::~AliForwarddNdetaTask()
+{
+ if (fSums) {
+ fSums->Delete();
+ delete fSums;
+ fSums = 0;
+ }
+ if (fOutputs) {
+ fOutputs->Delete();
+ delete fOutputs;
+ fOutputs = 0;
+ }
+ if (fCentral) delete fCentral;
+ if (fPrimary) delete fPrimary;
+}
+
+//________________________________________________________________________
+void
+AliForwarddNdetaTask::SetTriggerMask(const char* mask)
+{
+ UShort_t trgMask = 0;
+ TString trgs(mask);
+ trgs.ToUpper();
+ TObjString* trg;
+ TIter next(trgs.Tokenize(" ,|"));
+ while ((trg = static_cast<TObjString*>(next()))) {
+ TString s(trg->GetString());
+ if (s.IsNull()) continue;
+ if (s.CompareTo("INEL") == 0) trgMask = AliAODForwardMult::kInel;
+ else if (s.CompareTo("INEL>0")== 0) trgMask = AliAODForwardMult::kInelGt0;
+ else if (s.CompareTo("NSD") == 0) trgMask = AliAODForwardMult::kNSD;
+ else
+ Warning("SetTriggerMask", "Unknown trigger %s", s.Data());
+ }
+ if (trgMask == 0) trgMask = 1;
+ SetTriggerMask(trgMask);
+}
+
+//________________________________________________________________________
+void
+AliForwarddNdetaTask::UserCreateOutputObjects()
+{
+ // Create histograms
+ // Called once (on the worker node)
+
+ fOutput = new TList;
+ fOutput->SetName(Form("%s_result", GetName()));
+ fOutput->SetOwner();
+
+ fSums = new TList;
+ fSums->SetName(Form("%s_sums", GetName()));
+ fSums->SetOwner();
+
+
+ fTriggers = new TH1D("triggers", "Number of triggers",
+ kAccepted, 1, kAccepted);
+ fTriggers->SetYTitle("# of events");
+ fTriggers->GetXaxis()->SetBinLabel(kAll, "All events");
+ fTriggers->GetXaxis()->SetBinLabel(kB, "w/B trigger");
+ fTriggers->GetXaxis()->SetBinLabel(kA, "w/A trigger");
+ fTriggers->GetXaxis()->SetBinLabel(kC, "w/C trigger");
+ fTriggers->GetXaxis()->SetBinLabel(kE, "w/E trigger");
+ fTriggers->GetXaxis()->SetBinLabel(kMB, "w/Collision trigger");
+ fTriggers->GetXaxis()->SetBinLabel(kWithVertex, "w/Vertex");
+ fTriggers->GetXaxis()->SetBinLabel(kWithTrigger, "w/Selected trigger");
+ fTriggers->GetXaxis()->SetBinLabel(kAccepted, "Accepted by cut");
+ fTriggers->GetXaxis()->SetNdivisions(kAccepted, false);
+ fTriggers->SetFillColor(kRed+1);
+ fTriggers->SetFillStyle(3001);
+ fTriggers->SetStats(0);
+ fSums->Add(fTriggers);
+
+ fSNNString = new TNamed("sNN", "");
+ fSysString = new TNamed("sys", "");
+ fSums->Add(fSNNString);
+ fSums->Add(fSysString);
+
+ // Check that we have an AOD input handler
+ AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
+ AliAODInputHandler* ah =
+ dynamic_cast<AliAODInputHandler*>(am->GetInputEventHandler());
+ if (!ah) AliFatal("No AOD input handler set in analysis manager");
+
+ // Post data for ALL output slots >0 here, to get at least an empty histogram
+ PostData(1, fSums);
+}
+
+//____________________________________________________________________
+TH2D*
+AliForwarddNdetaTask::CloneHist(TH2D* in, const char* name)
+{
+ if (!in) return 0;
+ TH2D* ret = static_cast<TH2D*>(in->Clone(name));
+ ret->SetDirectory(0);
+ ret->Sumw2();
+ fSums->Add(ret);
+
+ return ret;
+}
+
+//____________________________________________________________________
+void
+AliForwarddNdetaTask::UserExec(Option_t *)
+{
+ // Main loop
+ AliAODEvent* aod = dynamic_cast<AliAODEvent*>(InputEvent());
+ if (!aod) {
+ AliError("Cannot get the AOD event");
+ return;
+ }
+
+ // Get objects from the event structure
+ TObject* oForward = aod->FindListObject("Forward");
+ TObject* oForwardMC = aod->FindListObject("ForwardMC");
+ TObject* oPrimary = aod->FindListObject("primary");
+ TObject* oCentral = aod->FindListObject("Central");
+
+ // We should have a forward object at least
+ if (!oForward) {
+ AliWarning("No Forward object found AOD");
+ return;
+ }
+
+ // Cast to good types
+ AliAODForwardMult* forward = static_cast<AliAODForwardMult*>(oForward);
+ AliAODForwardMult* forwardMC = static_cast<AliAODForwardMult*>(oForwardMC);
+ TH2D* primary = static_cast<TH2D*>(oPrimary);
+ TH2D* central = static_cast<TH2D*>(oCentral);
+
+ static bool first = true;
+ if (first) {
+ UShort_t sNN = forward->GetSNN();
+ UShort_t sys = forward->GetSystem();
+ fSNNString->SetTitle(AliForwardUtil::CenterOfMassEnergyString(sNN));
+ fSNNString->SetUniqueID(sNN);
+ fSysString->SetTitle(AliForwardUtil::CollisionSystemString(sys));
+ fSysString->SetUniqueID(sys);
+ first = false;
+ }
+
+ // Create our sum histograms
+ if (!fSumForward) fSumForward = CloneHist(&forward->GetHistogram(),"forward");
+ if (!fSumForwardMC && forwardMC)
+ fSumForwardMC = CloneHist(&forwardMC->GetHistogram(),"forwardMC");
+ if (!fSumPrimary && primary) fSumPrimary = CloneHist(primary, "primary");
+ if (!fSumCentral && central) fSumCentral = CloneHist(central, "central");
+
+ // Add contribtion from MC
+ if (primary) fSumPrimary->Add(primary);
+
+ // Count event
+ fTriggers->AddBinContent(kAll);
+ if (forward->IsTriggerBits(AliAODForwardMult::kB))
+ fTriggers->AddBinContent(kB);
+ if (forward->IsTriggerBits(AliAODForwardMult::kA))
+ fTriggers->AddBinContent(kA);
+ if (forward->IsTriggerBits(AliAODForwardMult::kC))
+ fTriggers->AddBinContent(kC);
+ if (forward->IsTriggerBits(AliAODForwardMult::kE))
+ fTriggers->AddBinContent(kE);
+ if (forward->IsTriggerBits(AliAODForwardMult::kInel))
+ fTriggers->AddBinContent(kMB);
+
+ // Check if we have an event of interest.
+ if (!forward->IsTriggerBits(fTriggerMask)) return;
+ fTriggers->AddBinContent(kWithTrigger);
+
+ // Check that we have a valid vertex
+ if (!forward->HasIpZ()) return;
+ fTriggers->AddBinContent(kWithVertex);
+
+ // Check that vertex is within cuts
+ if (!forward->InRange(fVtxMin, fVtxMax)) return;
+ fTriggers->AddBinContent(kAccepted);
+
+ // Add contribution
+ fSumForward->Add(&(forward->GetHistogram()));
+ if (fSumForwardMC) fSumForwardMC->Add(&(forwardMC->GetHistogram()));
+ if (fSumPrimary) fSumPrimary->Add(primary);
+ if (fSumCentral) fSumCentral->Add(central);
+
+ PostData(1, fSums);
+}
+
+//________________________________________________________________________
+void
+AliForwarddNdetaTask::SetHistogramAttributes(TH1D* h, Int_t colour, Int_t marker,
+ const char* title, const char* ytitle)
+{
+ h->SetTitle(title);
+ h->SetMarkerColor(colour);
+ h->SetMarkerStyle(marker);
+ h->SetMarkerSize(1);
+ h->SetFillStyle(0);
+ h->SetYTitle(ytitle);
+ h->GetXaxis()->SetTitleFont(132);
+ h->GetXaxis()->SetLabelFont(132);
+ h->GetXaxis()->SetNdivisions(10);
+ h->GetYaxis()->SetTitleFont(132);
+ h->GetYaxis()->SetLabelFont(132);
+ h->GetYaxis()->SetNdivisions(10);
+ h->GetYaxis()->SetDecimals();
+ h->SetStats(0);
+}
+
+//________________________________________________________________________
+void
+AliForwarddNdetaTask::Terminate(Option_t *)
+{
+ // Draw result to screen, or perform fitting, normalizations Called
+ // once at the end of the query
+
+ fSums = dynamic_cast<TList*> (GetOutputData(1));
+ if(!fOutput) {
+ AliError("Could not retrieve TList fSums");
+ return;
+ }
+
+ if (!fOutput) {
+ fOutput = new TList;
+ fOutput->SetName(Form("%s_result", GetName()));
+ fOutput->SetOwner();
+ }
+
+ fSumForward = static_cast<TH2D*>(fSums->FindObject("forward"));
+ fSumForwardMC = static_cast<TH2D*>(fSums->FindObject("forwardMC"));
+ fSumPrimary = static_cast<TH2D*>(fSums->FindObject("primary"));
+ fSumCentral = static_cast<TH2D*>(fSums->FindObject("central"));
+ fTriggers = static_cast<TH1D*>(fSums->FindObject("triggers"));
+ fSNNString = static_cast<TNamed*>(fSums->FindObject("sNN"));
+ fSysString = static_cast<TNamed*>(fSums->FindObject("sys"));
+
+ if (!fTriggers) {
+ AliError("Couldn't find histogram 'triggers' in list");
+ return;
+ }
+ if (!fSumForward) {
+ AliError("Couldn't find histogram 'forward' in list");
+ return;
+ }
+
+ Int_t nAll = Int_t(fTriggers->GetBinContent(kAll));
+ Int_t nB = Int_t(fTriggers->GetBinContent(kB));
+ Int_t nA = Int_t(fTriggers->GetBinContent(kA));
+ Int_t nC = Int_t(fTriggers->GetBinContent(kC));
+ Int_t nE = Int_t(fTriggers->GetBinContent(kE));
+ Int_t nMB = Int_t(fTriggers->GetBinContent(kMB));
+ Int_t nTriggered = Int_t(fTriggers->GetBinContent(kWithTrigger));
+ Int_t nWithVertex = Int_t(fTriggers->GetBinContent(kWithVertex));
+ Int_t nAccepted = Int_t(fTriggers->GetBinContent(kAccepted));
+ Int_t nGood = nB - nA - nC + 2 * nE;
+ Double_t vtxEff = Double_t(nMB) / nTriggered * Double_t(nAccepted) / nGood;
+ AliInfo(Form("Total of %9d events\n"
+ " of these %9d are minimum bias\n"
+ " of these %9d has a %s trigger\n"
+ " of these %9d has a vertex\n"
+ " of these %9d were in [%+4.1f,%+4.1f]cm\n"
+ " Triggers by type:\n"
+ " B = %9d\n"
+ " A|C = %9d (%9d+%-9d)\n"
+ " E = %9d\n"
+ " Implies %9d good triggers\n"
+ " Vertex efficiency: %f",
+ nAll, nMB, nTriggered,
+ AliAODForwardMult::GetTriggerString(fTriggerMask),
+ nWithVertex, nAccepted,
+ fVtxMin, fVtxMax,
+ nB, nA+nC, nA, nC, nE, nGood, vtxEff));
+
+ // Get acceptance normalisation from underflow bins
+ TH1D* normForward = fSumForward->ProjectionX("normForward", 0, 1, "");
+ // Project onto eta axis - _ignoring_underflow_bins_!
+ TH1D* dndetaForward = fSumForward->ProjectionX("dndetaForward", 1, -1, "e");
+ // Normalize to the acceptance
+ dndetaForward->Divide(normForward);
+ // Scale by the vertex efficiency
+ dndetaForward->Scale(vtxEff, "width");
+
+ SetHistogramAttributes(dndetaForward, kRed+1, 20, "ALICE Forward");
+ SetHistogramAttributes(normForward, kRed+1, 20, "ALICE Forward", "Normalisation");
+
+ fOutput->Add(fTriggers->Clone());
+ fOutput->Add(fSNNString->Clone());
+ fOutput->Add(fSysString->Clone());
+ fOutput->Add(dndetaForward);
+ fOutput->Add(normForward);
+ fOutput->Add(Rebin(dndetaForward));
+
+ if (fSumForwardMC) {
+ // Get acceptance normalisation from underflow bins
+ TH1D* normForwardMC = fSumForwardMC->ProjectionX("normForwardMC", 0, 1, "");
+ // Project onto eta axis - _ignoring_underflow_bins_!
+ TH1D* dndetaForwardMC = fSumForwardMC->ProjectionX("dndetaForwardMC", 1, -1, "e");
+ // Normalize to the acceptance
+ dndetaForwardMC->Divide(normForwardMC);
+ // Scale by the vertex efficiency
+ dndetaForwardMC->Scale(vtxEff, "width");
+
+ SetHistogramAttributes(dndetaForwardMC, kRed+3, 21, "ALICE Forward (MC)");
+ SetHistogramAttributes(normForwardMC, kRed+3, 21, "ALICE Forward (MC)",
+ "Normalisation");
+
+ fOutput->Add(dndetaForwardMC);
+ fOutput->Add(normForwardMC);
+ fOutput->Add(Rebin(dndetaForwardMC));
+ }
+
+ if (fSumPrimary) {
+ TH1D* dndetaTruth = fSumPrimary->ProjectionX("dndetaTruth",-1,-1,"e");
+ dndetaTruth->Scale(1./nAll, "width");
+
+ SetHistogramAttributes(dndetaTruth, kGray+3, 22, "Monte-Carlo truth");
+
+ fOutput->Add(dndetaTruth);
+ fOutput->Add(Rebin(dndetaTruth));
+ }
+ if (fSumCentral) {
+ TH1D* dndetaCentral = fSumCentral->ProjectionX("dndetaCentral",-1,-1,"e");
+ dndetaCentral->Scale(1./nGood, "width");
+
+ SetHistogramAttributes(dndetaCentral, kGray+3, 22, "ALICE Central - track(let)s");
+
+ dndetaCentral->GetXaxis()->SetRangeUser(-1,1);
+ fOutput->Add(dndetaCentral);
+ fOutput->Add(Rebin(dndetaCentral));
+ }
+
+ TNamed* trigString =
+ new TNamed("trigString", AliAODForwardMult::GetTriggerString(fTriggerMask));
+ trigString->SetUniqueID(fTriggerMask);
+ fOutput->Add(trigString);
+
+ TAxis* vtxAxis = new TAxis(1,fVtxMin,fVtxMax);
+ vtxAxis->SetName("vtxAxis");
+ vtxAxis->SetTitle(Form("v_{z}#in[%+5.1f,%+5.1f]cm", fVtxMin,fVtxMax));
+ fOutput->Add(vtxAxis);
+
+ // If only there was a away to get sqrt{s_NN} and beam type
+
+
+ PostData(2, fOutput);
+}
+
+//________________________________________________________________________
+TH1D*
+AliForwarddNdetaTask::Rebin(const TH1D* h) const
+{
+ if (fRebin <= 1) return 0;
+
+ Int_t nBins = h->GetNbinsX();
+ if(nBins % fRebin != 0) {
+ Warning("Rebin", "Rebin factor %d is not a devisor of current number "
+ "of bins %d in the histogram %s", fRebin, nBins, h->GetName());
+ return 0;
+ }
+
+ // Make a copy
+ TH1D* tmp = static_cast<TH1D*>(h->Clone(Form("%s_rebin%02d",
+ h->GetName(), fRebin)));
+ tmp->Rebin(fRebin);
+ tmp->SetDirectory(0);
+
+ // The new number of bins
+ Int_t nBinsNew = nBins / fRebin;
+ for(Int_t i = 1;i<= nBinsNew; i++) {
+ Double_t content = 0;
+ Double_t sumw = 0;
+ Double_t wsum = 0;
+ Int_t nbins = 0;
+ for(Int_t j = 1; j<=fRebin;j++) {
+ Int_t bin = (i-1)*fRebin + j;
+ Double_t c = h->GetBinContent(bin);
+
+ if (c <= 0) continue;
+
+ if (fCutEdges) {
+ if (h->GetBinContent(bin+1)<=0 ||
+ h->GetBinContent(bin-1)) {
+ Warning("Rebin", "removing bin %d=%f of %s (%d=%f,%d=%f)",
+ bin, c, h->GetName(),
+ bin+1, h->GetBinContent(bin+1),
+ bin-1, h->GetBinContent(bin-1));
+ continue;
+ }
+ }
+ Double_t e = h->GetBinError(bin);
+ Double_t w = 1 / (e*e); // 1/c/c
+ content += c;
+ sumw += w;
+ wsum += w * c;
+ nbins++;
+ }
+
+ if(content > 0 && nbins > 0) {
+ tmp->SetBinContent(i, wsum / sumw);
+ tmp->SetBinError(i,1./TMath::Sqrt(sumw));
+ }
+ }
+
+ return tmp;
+}
--- /dev/null
+//
+// Task to analyse the AOD for for dN/deta in the forward regions
+//
+#ifndef ALIFORWARDDNDETATASK_H
+#define ALIFORWARDDNDETATASK_H
+#include <AliAnalysisTaskSE.h>
+#include <AliAODForwardMult.h>
+class TList;
+class TH2D;
+class TH1D;
+
+/**
+ * Task to determine the
+ */
+class AliForwarddNdetaTask : public AliAnalysisTaskSE
+{
+public:
+ /**
+ * Constructor
+ *
+ */
+ AliForwarddNdetaTask();
+ /**
+ * Constructor
+ *
+ * @param name Name of task
+ * @param maxVtx Set @f$v_z@f$ range
+ */
+ AliForwarddNdetaTask(const char* name);
+
+ void SetVertexRange(Double_t min, Double_t max) { fVtxMin=min; fVtxMax=max; }
+ void SetRebinning(Int_t rebin) { fRebin = rebin; }
+ void SetTriggerMask(UShort_t mask) { fTriggerMask = mask; }
+ void SetTriggerMask(const char* mask);
+ /**
+ * Destructor
+ *
+ */
+ virtual ~AliForwarddNdetaTask();
+ /**
+ * Initialise on master - does nothing
+ *
+ */
+ virtual void Init() {}
+ /**
+ * Create output objects.
+ *
+ * This is called once per slave process
+ */
+ virtual void UserCreateOutputObjects();
+ /**
+ * Process a single event
+ *
+ * @param option Not used
+ */
+ virtual void UserExec(Option_t* option);
+ /**
+ * Called at end of event processing..
+ *
+ * This is called once in the master
+ *
+ * @param option Not used
+ */
+ virtual void Terminate(Option_t* option);
+protected:
+ AliForwarddNdetaTask(const AliForwarddNdetaTask&);
+ AliForwarddNdetaTask& operator=(const AliForwarddNdetaTask&) { return *this; }
+ /**
+ * Clone a 2D histogram
+ *
+ * @param in Histogram to clone.
+ * @param name New name of clone.
+ *
+ * @return The clone
+ */
+ TH2D* CloneHist(TH2D* in, const char* name);
+ /**
+ * Make a copy of the input histogram and rebin that histogram
+ *
+ * @param h Histogram to rebin
+ *
+ * @return New (rebinned) histogram
+ */
+ TH1D* Rebin(const TH1D* h) const;
+ /**
+ * Set histogram graphical options, etc.
+ *
+ * @param h Histogram to modify
+ * @param colour Marker color
+ * @param marker Marker style
+ * @param title Title of histogram
+ * @param ytitle Title on y-axis.
+ */
+ void SetHistogramAttributes(TH1D* h, Int_t colour, Int_t marker, const char* title,
+ const char* ytitle="#frac{1}{N} #frac{dN_{ch}}{d#eta}");
+ /**
+ * Trigger histogram bins
+ */
+ enum {
+ kAll = 1,
+ kB = 2,
+ kA = 3,
+ kC = 4,
+ kE = 5,
+ kMB = 6,
+ kWithTrigger= 7,
+ kWithVertex = 8,
+ kAccepted = 9
+ };
+
+ TH2D* fSumForward; // Sum of histograms
+ TH2D* fSumForwardMC; // Sum of MC histograms (if any)
+ TH2D* fSumPrimary; // Sum of primary histograms
+ TH2D* fSumCentral; // Sum of central histograms
+ TH2D* fCentral; //! Cache of central histogram
+ TH2D* fPrimary; //! Cache of primary histogram
+
+ TList* fSums; // Container of sums
+ TList* fOutput; // Container of outputs
+
+ TH1D* fTriggers; // Histogram of triggers
+
+ Double_t fVtxMin; // Minimum v_z
+ Double_t fVtxMax; // Maximum v_z
+ Int_t fTriggerMask; // Trigger mask
+ Int_t fRebin; // Rebinning factor
+ Bool_t fCutEdges; // Whether to cut edges when rebinning
+ TNamed* fSNNString; //
+ TNamed* fSysString; //
+
+ ClassDef(AliForwarddNdetaTask,1); // Determine multiplicity in central area
+};
+
+#endif
--- /dev/null
+#include <TH1.h>
+#include <THStack.h>
+#include <TGraphErrors.h>
+#include <TGraphAsymmErrors.h>
+#include <TMultiGraph.h>
+#include <TFile.h>
+#include <TList.h>
+#include <TString.h>
+#include <TError.h>
+#include <TSystem.h>
+#include <TROOT.h>
+#include <TMath.h>
+#include <TCanvas.h>
+#include <TPad.h>
+#include <TStyle.h>
+#include <TLegend.h>
+#include <TLegendEntry.h>
+#include <TLatex.h>
+#include <TImage.h>
+
+struct dNdetaDrawer
+{
+ struct RangeParam
+ {
+ TAxis* fMasterAxis; // Master axis
+ TAxis* fSlave1Axis; // First slave axis
+ TVirtualPad* fSlave1Pad; // First slave pad
+ TAxis* fSlave2Axis; // Second slave axis
+ TVirtualPad* fSlave2Pad; // Second slave pad
+ };
+ //__________________________________________________________________
+ dNdetaDrawer()
+ : fShowOthers(false), // Bool_t
+ fShowAlice(false), // Bool_t
+ fShowRatios(false), // Bool_t
+ fShowLeftRight(false), // Bool_t
+ fRebin(5), // UShort_t
+ fCutEdges(false), // Bool_t
+ fTitle(""), // TString
+ fHHDFile(""), // TString
+ fTrigString(0), // TNamed*
+ fSNNString(0), // TNamed*
+ fSysString(0), // TNamed*
+ fVtxAxis(0), // TAxis*
+ fForward(0), // TH1*
+ fForwardMC(0), // TH1*
+ fForwardHHD(0), // TH1*
+ fTruth(0), // TH1*
+ fCentral(0), // TH1*
+ fForwardSym(0), // TH1*
+ fForwardMCSym(0), // TH1*
+ fForwardHHDSym(0), // TH1*
+ fTriggers(0), // TH1*
+ fRangeParam(0)
+
+ {
+ fRangeParam = new RangeParam;
+ fRangeParam->fMasterAxis = 0;
+ fRangeParam->fSlave1Axis = 0;
+ fRangeParam->fSlave1Pad = 0;
+ fRangeParam->fSlave2Axis = 0;
+ fRangeParam->fSlave2Pad = 0;
+ }
+ void SetShowOthers(Bool_t x) { fShowOthers = x; }
+ void SetShowAlice(Bool_t x) { fShowAlice = x; }
+ void SetShowRatios(Bool_t x) { fShowRatios = x; }
+ void SetShowLeftRight(Bool_t x) { fShowLeftRight = x; }
+ void SetRebin(UShort_t x) { fRebin = x; }
+ void SetCutEdges(Bool_t x) { fCutEdges = x; }
+ void SetTitle(TString x) { fTitle = x; }
+ void SetHHDFile(const char* fn) { fHHDFile = fn; }
+
+ //__________________________________________________________________
+ void Run(const char* filename)
+ {
+ if (!Open(filename)) return;
+
+ Double_t max = 0;
+
+ // Create our stack of results
+ THStack* results = StackResults(max);
+
+ // Create our stack of other results
+ TMultiGraph* other = 0;
+ if (fShowOthers || fShowAlice) other = StackOther(max);
+
+ Double_t smax = 0;
+ THStack* ratios = 0;
+ if (fShowRatios) ratios = StackRatios(other, smax);
+
+ Double_t amax = 0;
+ THStack* leftright = 0;
+ if (fShowLeftRight) leftright = StackLeftRight(amax);
+
+ Plot(results, other, max, ratios, smax, leftright, amax);
+ }
+
+ //__________________________________________________________________
+ Bool_t Open(const char* filename)
+ {
+ TFile* file = TFile::Open(filename, "READ");
+ if (!file) {
+ Error("Open", "Cannot open %s", filename);
+ return false;
+ }
+
+ TList* results = static_cast<TList*>(file->Get("ForwardResults"));
+ if (!results) {
+ Error("Open", "Couldn't find list ForwardResults");
+ return false;
+ }
+
+ fForward = GetResult(results, "dndetaForward");
+ fForwardMC = GetResult(results, "dndetaForwardMC");
+ fTruth = GetResult(results, "dndetaTruth");
+ fCentral = GetResult(results, "dndetaCentral");
+
+ fTrigString = static_cast<TNamed*>(results->FindObject("trigString"));
+ fSNNString = static_cast<TNamed*>(results->FindObject("sNN"));
+ fSysString = static_cast<TNamed*>(results->FindObject("sys"));
+ fVtxAxis = static_cast<TAxis*>(results->FindObject("vtxAxis"));
+
+ TList* sums = static_cast<TList*>(file->Get("ForwardSums"));
+ if (sums)
+ fTriggers = GetResult(sums, "triggers");
+
+ if (!fForward) {
+ Error("Open", "Couldn't find the result of the forward analysis");
+ return false;
+ }
+ file->Close();
+
+
+ fForwardHHD = GetHHD();
+
+ return true;
+ }
+ //__________________________________________________________________
+ TH1* GetResult(TList* list, const char* name) const
+ {
+ if (!list) return 0;
+
+ TH1* ret = static_cast<TH1*>(list->FindObject(name));
+ if (!ret)
+ Warning("GetResult", "Histogram %s not found", name);
+
+ return ret;
+ }
+ //__________________________________________________________________
+ /**
+ * Get the result from previous analysis code
+ *
+ * @param fn File to open
+ * @param nsd Whether this is NSD
+ *
+ * @return null or result of previous analysis code
+ */
+ TH1* GetHHD()
+ {
+ if (fHHDFile.IsNull()) return 0;
+ const char* fn = fHHDFile.Data();
+ Bool_t nsd = (fTrigString ? fTrigString->GetUniqueID() & 0x4 : false);
+ TDirectory* savdir = gDirectory;
+ if (gSystem->AccessPathName(fn)) {
+ Warning("GetHHD", "Output of HHD analysis (%s) not available", fn);
+ return 0;
+ }
+ TFile* file = TFile::Open(fn, "READ");
+ if (!file) {
+ Warning("GetHHD", "couldn't open HHD file %s", fn);
+ return 0;
+ }
+ TString hist(Form("dNdeta_dNdeta%s", (nsd ? "NSD" : "")));
+ TH1* h = static_cast<TH1*>(file->Get(hist.Data()));
+ if (!h) {
+ Warning("GetHHD", "Couldn't find HHD histogram %s in %s",
+ hist.Data(), fn);
+ file->Close();
+ savdir->cd();
+ return 0;
+ }
+ TH1* r = static_cast<TH1*>(h->Clone("dndeta_hhd"));
+ r->SetTitle("ALICE Forward (HHD)");
+ r->SetFillStyle(0);
+ r->SetFillColor(0);
+ r->SetMarkerStyle(21);
+ r->SetMarkerColor(kPink+1);
+ r->SetDirectory(0);
+
+ file->Close();
+ savdir->cd();
+ return r;
+ }
+ //__________________________________________________________________
+ THStack* StackResults(Double_t& max)
+ {
+ THStack* stack = new THStack("results", "Stack of Results");
+ max = TMath::Max(max, AddHistogram(stack, fTruth, "e5 p"));
+ max = TMath::Max(max, AddHistogram(stack, fForwardHHD, "", fForwardHHDSym));
+ max = TMath::Max(max, AddHistogram(stack, fForwardMC, "", fForwardMCSym));
+ max = TMath::Max(max, AddHistogram(stack, fCentral, ""));
+ max = TMath::Max(max, AddHistogram(stack, fForward, "", fForwardSym));
+ return stack;
+ }
+ //__________________________________________________________________
+ TMultiGraph* StackOther(Double_t& max) const
+ {
+ gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/OtherData.C");
+ Int_t error = 0;
+ Bool_t onlya = (fShowOthers ? false : true);
+ Int_t trg = (fTrigString ? fTrigString->GetUniqueID() : 0);
+ UShort_t snn = (fSNNString ? fSNNString->GetUniqueID() : 0);
+ Long_t ret = gROOT->ProcessLine(Form("GetData(%d,%d,%d);",
+ snn,trg,onlya));
+ if (error) {
+ Error("StackOther", "Failed to execute GetData(%d,%d,%d)",
+ snn, trg, onlya);
+ return 0;
+ }
+ if (!ret) {
+ Warning("StackOther", "No other data found for sNN=%d, trigger=%d",
+ snn, trg);
+ return 0;
+ }
+ TMultiGraph* other = reinterpret_cast<TMultiGraph*>(ret);
+
+ TGraphAsymmErrors* o = 0;
+ TIter next(other->GetListOfGraphs());
+ while ((o = static_cast<TGraphAsymmErrors*>(next())))
+ max = TMath::Max(max, TMath::MaxElement(o->GetN(), o->GetY()));
+
+ return other;
+ }
+ //__________________________________________________________________
+ THStack* StackRatios(TMultiGraph* others, Double_t& max)
+ {
+ THStack* ratios = new THStack("ratios", "Ratios");
+
+ if (others) {
+ TGraphAsymmErrors* ua5_1 = 0;
+ TGraphAsymmErrors* ua5_2 = 0;
+ TGraphAsymmErrors* alice = 0;
+ TGraphAsymmErrors* cms = 0;
+ TGraphAsymmErrors* o = 0;
+ TIter nextG(others->GetListOfGraphs());
+ while ((o = static_cast<TGraphAsymmErrors*>(nextG()))) {
+ ratios->Add(Ratio(fForward, o, max));
+ ratios->Add(Ratio(fForwardSym, o, max));
+ ratios->Add(Ratio(fForwardHHD, o, max));
+ ratios->Add(Ratio(fForwardHHDSym, o, max));
+ ratios->Add(Ratio(fCentral, o, max));
+ TString oName(o->GetName());
+ oName.ToLower();
+ if (oName.Contains("ua5")) { if (ua5_1) ua5_2 = o; else ua5_1 = o; }
+ if (oName.Contains("alice")) alice = o;
+ if (oName.Contains("cms")) cms = o;
+ }
+ if (ua5_1 && alice) ratios->Add(Ratio(alice, ua5_1, max));
+ if (ua5_2 && alice) ratios->Add(Ratio(alice, ua5_2, max));
+ if (cms && alice) ratios->Add(Ratio(alice, cms, max));
+ }
+
+ // If we have data from HHD's analysis, then do the ratio of
+ // our result to that data.
+ if (fForwardHHD) {
+ ratios->Add(Ratio(fForward, fForwardHHD, max));
+ ratios->Add(Ratio(fForwardSym, fForwardHHDSym, max));
+ }
+
+ // Do comparison to MC
+ if (fForwardMC) {
+ ratios->Add(Ratio(fForward, fForwardMC, max));
+ ratios->Add(Ratio(fForwardSym, fForwardMCSym, max));
+ }
+
+ // Check if we have ratios
+ if (!ratios->GetHists() ||
+ (ratios->GetHists()->GetEntries() <= 0)) {
+ delete ratios;
+ ratios = 0;
+ }
+ return ratios;
+ }
+ //__________________________________________________________________
+ THStack* StackLeftRight(Double_t& max)
+ {
+ THStack* ret = new THStack("leftright", "Left-right asymmetry");
+ ret->Add(Asymmetry(fForward, max));
+ ret->Add(Asymmetry(fForwardHHD, max));
+ ret->Add(Asymmetry(fForwardMC, max));
+
+ if (!ret->GetHists() ||
+ (ret->GetHists()->GetEntries() <= 0)) {
+ delete ret;
+ ret = 0;
+ }
+ return ret;
+ }
+ //__________________________________________________________________
+ TAxis* FindXAxis(TVirtualPad* p, const char* name)
+ {
+ TObject* o = p->GetListOfPrimitives()->FindObject(name);
+ if (!o) {
+ Warning("FindXAxis", "%s not found in pad", name);
+ return 0;
+ }
+ THStack* stack = dynamic_cast<THStack*>(o);
+ if (!stack) {
+ Warning("FindXAxis", "%s is not a THStack", name);
+ return 0;
+ }
+ if (!stack->GetHistogram()) {
+ Warning("FindXAxis", "%s has no histogram", name);
+ return 0;
+ }
+ TAxis* ret = stack->GetHistogram()->GetXaxis();
+ return ret;
+ }
+ //__________________________________________________________________
+ void Plot(THStack* results,
+ TMultiGraph* others,
+ Double_t max,
+ THStack* ratios,
+ Double_t rmax,
+ THStack* leftright,
+ Double_t amax)
+ {
+ gStyle->SetOptTitle(0);
+ gStyle->SetTitleFont(132, "xyz");
+ gStyle->SetLabelFont(132, "xyz");
+
+ Int_t h = 800;
+ Int_t w = 800; // h / TMath::Sqrt(2);
+ if (!ratios) w *= 1.4;
+ if (!leftright) w *= 1.4;
+ TCanvas* c = new TCanvas("Results", "Results", w, h);
+ c->SetFillColor(0);
+ c->SetBorderSize(0);
+ c->SetBorderMode(0);
+
+ Double_t y1 = 0;
+ Double_t y2 = 0;
+ Double_t y3 = 0;
+ if (ratios) y1 = 0.3;
+ if (leftright) {
+ if (y1 > 0.0001) {
+ y2 = 0.2;
+ y1 = 0.4;
+ }
+ else {
+ y1 = 0.2;
+ y2 = y1;
+ }
+ }
+ PlotResults(results, others, max, y1);
+ c->cd();
+
+ PlotRatios(ratios, rmax, y2, y1);
+ c->cd( );
+
+ PlotLeftRight(leftright, amax, y3, y2);
+ c->cd();
+ }
+ //__________________________________________________________________
+ void PlotResults(THStack* results, TMultiGraph* others,
+ Double_t max, Double_t yd)
+ {
+ // Make a sub-pad for the result itself
+ TPad* p1 = new TPad("p1", "p1", 0, yd, 1.0, 1.0, 0, 0, 0);
+ p1->SetTopMargin(0.05);
+ p1->SetBorderSize(0);
+ p1->SetBorderMode(0);
+ p1->SetBottomMargin(yd > 0.001 ? 0.001 : 0.1);
+ p1->SetRightMargin(0.05);
+ p1->SetGridx();
+ p1->SetTicks(1,1);
+ p1->SetNumber(1);
+ p1->Draw();
+ p1->cd();
+
+ results->SetMaximum(1.15*max);
+ results->SetMinimum(yd > 0.00001 ? -0.1 : 0);
+
+ FixAxis(results, 1/(1-yd)/1.7, "#frac{1}{N} #frac{dN_{ch}}{d#eta}");
+
+ p1->Clear();
+ results->DrawClone("nostack e1");
+
+ fRangeParam->fSlave1Axis = results->GetXaxis();
+ fRangeParam->fSlave1Pad = p1;
+
+ // Draw other data
+ if (others) {
+ TGraphAsymmErrors* o = 0;
+ TIter nextG(others->GetListOfGraphs());
+ while ((o = static_cast<TGraphAsymmErrors*>(nextG())))
+ o->DrawClone("same p");
+ }
+
+ // Make a legend in the main result pad
+ TLegend* l = p1->BuildLegend(.15,p1->GetBottomMargin()+.01,.90,.35);
+ l->SetNColumns(2);
+ l->SetFillColor(0);
+ l->SetFillStyle(0);
+ l->SetBorderSize(0);
+ l->SetTextFont(132);
+
+ // Put a title on top
+ TLatex* tit = new TLatex(0.10, 0.95, fTitle.Data());
+ tit->SetNDC();
+ tit->SetTextFont(132);
+ tit->SetTextSize(0.05);
+ tit->Draw();
+
+ // Put a nice label in the plot
+ TString eS;
+ UShort_t snn = (fSNNString ? fSNNString->GetUniqueID() : 0);
+ const char* sys = (fSysString ? fSysString->GetTitle() : "?");
+ if (snn > 1000) eS = Form("%4.2fTeV", float(snn)/1000);
+ else eS = Form("%3dGeV", snn);
+ TLatex* tt = new TLatex(.93, .93, Form("%s #sqrt{s}=%s, %s",
+ sys,
+ eS.Data(),
+ fTrigString ?
+ fTrigString->GetTitle() : "?"));
+ tt->SetNDC();
+ tt->SetTextFont(132);
+ tt->SetTextAlign(33);
+ tt->Draw();
+
+ // Put number of accepted events on the plot
+ Int_t nev = fTriggers->GetBinContent(fTriggers->GetNbinsX());
+ TLatex* et = new TLatex(.93, .83, Form("%d events", nev));
+ et->SetNDC();
+ et->SetTextFont(132);
+ et->SetTextAlign(33);
+ et->Draw();
+
+ // Put number of accepted events on the plot
+ if (fVtxAxis) {
+ TLatex* vt = new TLatex(.93, .88, fVtxAxis->GetTitle());
+ vt->SetNDC();
+ vt->SetTextFont(132);
+ vt->SetTextAlign(33);
+ vt->Draw();
+ }
+ // results->Draw("nostack e1 same");
+
+ fRangeParam->fSlave1Axis = FindXAxis(p1, results->GetName());
+ fRangeParam->fSlave1Pad = p1;
+
+
+ // Mark the plot as preliminary
+ TLatex* pt = new TLatex(.12, .93, "Preliminary");
+ pt->SetNDC();
+ pt->SetTextFont(22);
+ pt->SetTextSize(0.07);
+ pt->SetTextColor(kRed+1);
+ pt->SetTextAlign(13);
+ pt->Draw();
+
+ if (!gSystem->AccessPathName("ALICE.png")) {
+ TPad* logo = new TPad("logo", "logo", .12, .7, .32, .90, 0, 0, 0);
+ logo->SetFillStyle(0);
+ logo->Draw();
+ logo->cd();
+ TImage* i = TImage::Create();
+ i->ReadImage("ALICE.png");
+ i->Draw();
+ }
+ p1->cd();
+ }
+ //__________________________________________________________________
+ void PlotRatios(THStack* ratios, Double_t max, Double_t y1, Double_t y2)
+ {
+ if (!ratios) return;
+ bool isBottom = (y1 < 0.0001);
+ Double_t yd = y2 - y1;
+ // Make a sub-pad for the result itself
+ TPad* p2 = new TPad("p2", "p2", 0, y1, 1.0, y2, 0, 0, 0);
+ p2->SetTopMargin(0.001);
+ p2->SetRightMargin(0.05);
+ p2->SetBottomMargin(isBottom ? 1/yd * 0.07 : 0.0001);
+ p2->SetGridx();
+ p2->SetTicks(1,1);
+ p2->SetNumber(2);
+ p2->Draw();
+ p2->cd();
+
+ // Fix up axis
+ FixAxis(ratios, 1/yd/1.7, "Ratios", 7);
+
+ ratios->SetMaximum(1+TMath::Max(.22,1.05*max));
+ ratios->SetMinimum(1-TMath::Max(.32,1.05*max));
+ p2->Clear();
+ ratios->DrawClone("nostack e1");
+
+
+ // Make a legend
+ TLegend* l2 = p2->BuildLegend(.15,p2->GetBottomMargin()+.01,.9,
+ isBottom ? .6 : .4);
+ l2->SetNColumns(2);
+ l2->SetFillColor(0);
+ l2->SetFillStyle(0);
+ l2->SetBorderSize(0);
+ l2->SetTextFont(132);
+
+ // Make a nice band from 0.9 to 1.1
+ TGraphErrors* band = new TGraphErrors(2);
+ band->SetPoint(0, fForwardSym->GetXaxis()->GetXmin(), 1);
+ band->SetPoint(1, fForward->GetXaxis()->GetXmax(), 1);
+ band->SetPointError(0, 0, .1);
+ band->SetPointError(1, 0, .1);
+ band->SetFillColor(kYellow+2);
+ band->SetFillStyle(3002);
+ band->SetLineStyle(2);
+ band->SetLineWidth(1);
+ band->Draw("3 same");
+ band->DrawClone("X L same");
+
+ // Replot the ratios on top
+ ratios->DrawClone("nostack e1 same");
+
+ if (isBottom) {
+ fRangeParam->fMasterAxis = FindXAxis(p2, ratios->GetName());
+ p2->AddExec("range", Form("RangeExec((dNdetaDrawer::RangeParam*)%p)",
+ fRangeParam));
+ }
+ else {
+ fRangeParam->fSlave2Axis = FindXAxis(p2, ratios->GetName());
+ fRangeParam->fSlave2Pad = p2;
+ }
+ }
+ //__________________________________________________________________
+ void PlotLeftRight(THStack* leftright, Double_t max,
+ Double_t y1, Double_t y2)
+ {
+ if (!leftright) return;
+ bool isBottom = (y1 < 0.0001);
+ Double_t yd = y2 - y1;
+ // Make a sub-pad for the result itself
+ TPad* p3 = new TPad("p3", "p3", 0, y1, 1.0, y2, 0, 0, 0);
+ p3->SetTopMargin(0.001);
+ p3->SetRightMargin(0.05);
+ p3->SetBottomMargin(isBottom ? 1/yd * 0.07 : 0.0001);
+ p3->SetGridx();
+ p3->SetTicks(1,1);
+ p3->SetNumber(2);
+ p3->Draw();
+ p3->cd();
+
+ TH1* dummy = 0;
+ if (leftright->GetHists()->GetEntries() == 1) {
+ // Add dummy histogram
+ dummy = new TH1F("dummy","", 10, -6, 6);
+ dummy->SetLineColor(0);
+ dummy->SetFillColor(0);
+ dummy->SetMarkerColor(0);
+ leftright->Add(dummy);
+ }
+
+ // Fix up axis
+ FixAxis(leftright, 1/yd/1.7, "Right/Left", 4);
+
+ leftright->SetMaximum(1+TMath::Max(.12,1.05*max));
+ leftright->SetMinimum(1-TMath::Max(.15,1.05*max));
+ p3->Clear();
+ leftright->DrawClone("nostack e1");
+
+
+ // Make a legend
+ TLegend* l2 = p3->BuildLegend(.15,p3->GetBottomMargin()+.01,.9,.5);
+ l2->SetNColumns(2);
+ l2->SetFillColor(0);
+ l2->SetFillStyle(0);
+ l2->SetBorderSize(0);
+ l2->SetTextFont(132);
+#ifndef __CINT__
+ if (dummy) {
+ TList* prims = l2->GetListOfPrimitives();
+ TIter next(prims);
+ TLegendEntry* o = 0;
+ while ((o = static_cast<TLegendEntry*>(next()))) {
+ TString lbl(o->GetLabel());
+ if (lbl != "dummy") continue;
+ prims->Remove(o);
+ break;
+ }
+ }
+#endif
+ // Make a nice band from 0.9 to 1.1
+ TGraphErrors* band = new TGraphErrors(2);
+ band->SetPoint(0, fForwardSym->GetXaxis()->GetXmin(), 1);
+ band->SetPoint(1, fForward->GetXaxis()->GetXmax(), 1);
+ band->SetPointError(0, 0, .05);
+ band->SetPointError(1, 0, .05);
+ band->SetFillColor(kYellow+2);
+ band->SetFillStyle(3002);
+ band->SetLineStyle(2);
+ band->SetLineWidth(1);
+ band->Draw("3 same");
+ band->DrawClone("X L same");
+
+ leftright->DrawClone("nostack e1 same");
+ if (isBottom) {
+ fRangeParam->fMasterAxis = FindXAxis(p3, leftright->GetName());
+ p3->AddExec("range", Form("RangeExec((dNdetaDrawer::RangeParam*)%p)",
+ fRangeParam));
+ }
+ else {
+ fRangeParam->fSlave2Axis = FindXAxis(p3, leftright->GetName());
+ fRangeParam->fSlave2Pad = p3;
+ }
+ }
+
+ //__________________________________________________________________
+ /**
+ * Fix the apperance of the axis in a stack
+ *
+ * @param stack stack of histogram
+ * @param s Scaling factor
+ * @param ytitle Y axis title
+ * @param force Whether to draw the stack first or not
+ * @param ynDiv Divisions on Y axis
+ */
+ void FixAxis(THStack* stack, Double_t s, const char* ytitle,
+ Int_t ynDiv=210, Bool_t force=true)
+ {
+ if (force) stack->Draw("nostack e1");
+
+ TH1* h = stack->GetHistogram();
+ if (!h) return;
+
+ h->SetXTitle("#eta");
+ h->SetYTitle(ytitle);
+ TAxis* xa = h->GetXaxis();
+ TAxis* ya = h->GetYaxis();
+ if (xa) {
+ xa->SetTitle("#eta");
+ // xa->SetTicks("+-");
+ xa->SetTitleSize(s*xa->GetTitleSize());
+ xa->SetLabelSize(s*xa->GetLabelSize());
+ xa->SetTickLength(s*xa->GetTickLength());
+ }
+ if (ya) {
+ ya->SetTitle(ytitle);
+ ya->SetDecimals();
+ // ya->SetTicks("+-");
+ ya->SetNdivisions(ynDiv);
+ ya->SetTitleSize(s*ya->GetTitleSize());
+ ya->SetTitleOffset(ya->GetTitleOffset()/s);
+ ya->SetLabelSize(s*ya->GetLabelSize());
+ }
+ }
+
+ //__________________________________________________________________
+ Double_t AddHistogram(THStack* stack, TH1* hist, Option_t* option) const
+ {
+ // Check if we have input
+ if (!hist) return 0;
+
+ // Rebin if needed
+ Rebin(hist);
+
+ stack->Add(hist, option);
+ return hist->GetMaximum();
+ }
+ //__________________________________________________________________
+ Double_t AddHistogram(THStack* stack, TH1* hist, Option_t* option,
+ TH1*& sym) const
+ {
+ // Check if we have input
+ if (!hist) return 0;
+
+ // Rebin if needed
+ Rebin(hist);
+ stack->Add(hist, option);
+
+ // Now symmetrice the histogram
+ sym = Symmetrice(hist);
+ stack->Add(sym, option);
+
+ return hist->GetMaximum();
+ }
+
+ //__________________________________________________________________
+ /**
+ * Rebin a histogram
+ *
+ * @param h Histogram to rebin
+ * @param rebin Rebinning factor
+ *
+ * @return
+ */
+ virtual void Rebin(TH1* h) const
+ {
+ if (fRebin <= 1) return;
+
+ Int_t nBins = h->GetNbinsX();
+ if(nBins % fRebin != 0) {
+ Warning("Rebin", "Rebin factor %d is not a devisor of current number "
+ "of bins %d in the histogram %s", fRebin, nBins, h->GetName());
+ return;
+ }
+
+ // Make a copy
+ TH1* tmp = static_cast<TH1*>(h->Clone("tmp"));
+ tmp->Rebin(fRebin);
+ tmp->SetDirectory(0);
+
+ // The new number of bins
+ Int_t nBinsNew = nBins / fRebin;
+ for(Int_t i = 1;i<= nBinsNew; i++) {
+ Double_t content = 0;
+ Double_t sumw = 0;
+ Double_t wsum = 0;
+ Int_t nbins = 0;
+ for(Int_t j = 1; j<=fRebin;j++) {
+ Int_t bin = (i-1)*fRebin + j;
+ Double_t c = h->GetBinContent(bin);
+
+ if (c <= 0) continue;
+
+ if (fCutEdges) {
+ if (h->GetBinContent(bin+1)<=0 ||
+ h->GetBinContent(bin-1)) {
+ Warning("Rebin", "removing bin %d=%f of %s (%d=%f,%d=%f)",
+ bin, c, h->GetName(),
+ bin+1, h->GetBinContent(bin+1),
+ bin-1, h->GetBinContent(bin-1));
+ continue;
+ }
+ }
+ Double_t e = h->GetBinError(bin);
+ Double_t w = 1 / (e*e); // 1/c/c
+ content += c;
+ sumw += w;
+ wsum += w * c;
+ nbins++;
+ }
+
+ if(content > 0 && nbins > 0) {
+ tmp->SetBinContent(i, wsum / sumw);
+ tmp->SetBinError(i,1./TMath::Sqrt(sumw));
+ }
+ }
+
+ // Finally, rebin the histogram, and set new content
+ h->Rebin(fRebin);
+ h->Reset();
+ for(Int_t i = 1; i<= nBinsNew; i++) {
+ h->SetBinContent(i,tmp->GetBinContent(i));
+ h->SetBinError(i, tmp->GetBinError(i));
+ }
+
+ delete tmp;
+ }
+ //__________________________________________________________________
+ /**
+ * Make an extension of @a h to make it symmetric about 0
+ *
+ * @param h Histogram to symmertrice
+ *
+ * @return Symmetric extension of @a h
+ */
+ TH1* Symmetrice(const TH1* h) const
+ {
+ Int_t nBins = h->GetNbinsX();
+ TH1* s = static_cast<TH1*>(h->Clone(Form("%s_mirror", h->GetName())));
+ s->SetTitle(Form("%s (mirrored)", h->GetTitle()));
+ s->Reset();
+ s->SetBins(nBins, -h->GetXaxis()->GetXmax(), -h->GetXaxis()->GetXmin());
+ s->SetMarkerColor(h->GetMarkerColor());
+ s->SetMarkerSize(h->GetMarkerSize());
+ s->SetMarkerStyle(h->GetMarkerStyle()+4);
+ s->SetFillColor(h->GetFillColor());
+ s->SetFillStyle(h->GetFillStyle());
+ s->SetDirectory(0);
+
+ // Find the first and last bin with data
+ Int_t first = nBins+1;
+ Int_t last = 0;
+ for (Int_t i = 1; i <= nBins; i++) {
+ if (h->GetBinContent(i) <= 0) continue;
+ first = TMath::Min(first, i);
+ last = TMath::Max(last, i);
+ }
+
+ Double_t xfirst = h->GetBinCenter(first-1);
+ Int_t f1 = h->GetXaxis()->FindBin(-xfirst);
+ Int_t l2 = s->GetXaxis()->FindBin(xfirst);
+ for (Int_t i = f1, j=l2; i <= last; i++,j--) {
+ s->SetBinContent(j, h->GetBinContent(i));
+ s->SetBinError(j, h->GetBinError(i));
+ }
+ // Fill in overlap bin
+ s->SetBinContent(l2+1, h->GetBinContent(first));
+ s->SetBinError(l2+1, h->GetBinError(first));
+ return s;
+ }
+ //__________________________________________________________________
+ Double_t RatioMax(TH1* h) const
+ {
+ Int_t nBins = h->GetNbinsX();
+ Double_t ret = 0;
+ for (Int_t i = 1; i <= nBins; i++) {
+ Double_t c = h->GetBinContent(i);
+ if (c == 0) continue;
+ Double_t e = h->GetBinError(i);
+ Double_t d = TMath::Abs(1-c-e);
+ ret = TMath::Max(d, ret);
+ }
+ return ret;
+ }
+ //__________________________________________________________________
+ /**
+ * Compute the ratio of @a h to @a g. @a g is evaluated at the bin
+ * centers of @a h
+ *
+ * @param h Numerator
+ * @param g Divisor
+ *
+ * @return h/g
+ */
+ TH1* Ratio(const TH1* h, const TGraph* g, Double_t& max) const
+ {
+ if (!h || !g) return 0;
+
+ TH1* ret = static_cast<TH1*>(h->Clone("tmp"));
+ ret->SetName(Form("%s_over_%s", h->GetName(), g->GetName()));
+ ret->SetTitle(Form("%s / %s", h->GetTitle(), g->GetTitle()));
+ ret->Reset();
+ ret->SetMarkerStyle(g->GetMarkerStyle());
+ ret->SetMarkerColor(h->GetMarkerColor());
+ ret->SetMarkerSize(0.9*g->GetMarkerSize());
+ Double_t xlow = g->GetX()[0];
+ Double_t xhigh = g->GetX()[g->GetN()-1];
+ if (xlow > xhigh) { Double_t t = xhigh; xhigh = xlow; xlow = t; }
+
+ for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
+ Double_t c = h->GetBinContent(i);
+ if (c <= 0) continue;
+
+ Double_t x = h->GetBinCenter(i);
+ if (x < xlow || x > xhigh) continue;
+
+ Double_t f = g->Eval(x);
+ if (f <= 0) continue;
+
+ ret->SetBinContent(i, c / f);
+ ret->SetBinError(i, h->GetBinError(i) / f);
+ }
+ if (ret->GetEntries() <= 0) {
+ delete ret;
+ ret = 0;
+ }
+ else
+ max = TMath::Max(RatioMax(ret), max);
+ return ret;
+ }
+ //__________________________________________________________________
+ /**
+ * Make the ratio of h1 to h2
+ *
+ * @param h1 First histogram (numerator)
+ * @param h2 Second histogram (denominator)
+ *
+ * @return h1 / h2
+ */
+ TH1* Ratio(const TH1* h1, const TH1* h2, Double_t& max) const
+ {
+ if (!h1 || !h2) return 0;
+ TH1* t1 = static_cast<TH1*>(h1->Clone(Form("%s_%s",
+ h1->GetName(),
+ h2->GetName())));
+ t1->SetTitle(Form("%s / %s", h1->GetTitle(), h2->GetTitle()));
+ t1->Divide(h2);
+ t1->SetMarkerColor(h1->GetMarkerColor());
+ t1->SetMarkerStyle(h2->GetMarkerStyle());
+ max = TMath::Max(RatioMax(t1), max);
+ return t1;
+ }
+ //__________________________________________________________________
+ /**
+ * Calculate the ratio of two graphs - g1 / g2
+ *
+ * @param g1 Numerator
+ * @param g2 Denominator
+ *
+ * @return g1 / g2 in a histogram
+ */
+ TH1* Ratio(const TGraphAsymmErrors* g1,
+ const TGraphAsymmErrors* g2, Double_t& max) const
+ {
+ Int_t nBins = g1->GetN();
+ TArrayF bins(nBins+1);
+ Double_t dx = 0;
+ for (Int_t i = 0; i < nBins; i++) {
+ Double_t x = g1->GetX()[i];
+ Double_t exl = g1->GetEXlow()[i];
+ Double_t exh = g1->GetEXhigh()[i];
+ bins.fArray[i] = x-exl;
+ bins.fArray[i+1] = x+exh;
+ Double_t dxi = exh+exl;
+ if (i == 0) dx = dxi;
+ else if (dxi != dx) dx = 0;
+ }
+ TString name(Form("%s_%s", g1->GetName(), g2->GetName()));
+ TString title(Form("%s / %s", g1->GetTitle(), g2->GetTitle()));
+ TH1* h = 0;
+ if (dx != 0) {
+ h = new TH1F(name.Data(), title.Data(), nBins, bins[0], bins[nBins]);
+ }
+ else {
+ h = new TH1F(name.Data(), title.Data(), nBins, bins.fArray);
+ }
+ h->SetMarkerStyle(g2->GetMarkerStyle());
+ h->SetMarkerColor(g1->GetMarkerColor());
+ h->SetMarkerSize(0.9*g2->GetMarkerSize());
+
+ Double_t low = g2->GetX()[0];
+ Double_t high = g2->GetX()[g2->GetN()-1];
+ if (low > high) { Double_t t = low; low = high; high = t; }
+ for (Int_t i = 0; i < nBins; i++) {
+ Double_t x = g1->GetX()[i];
+ if (x < low-dx || x > high+dx) continue;
+ Double_t c1 = g1->GetY()[i];
+ Double_t e1 = g1->GetErrorY(i);
+ Double_t c2 = g2->Eval(x);
+
+ h->SetBinContent(i+1, c1 / c2);
+ h->SetBinError(i+1, e1 / c2);
+ }
+ max = TMath::Max(RatioMax(h), max);
+ return h;
+ }
+
+ //__________________________________________________________________
+ TH1* Graph2Hist(const TGraphAsymmErrors* g) const
+ {
+ Int_t nBins = g->GetN();
+ TArrayF bins(nBins+1);
+ Double_t dx = 0;
+ for (Int_t i = 0; i < nBins; i++) {
+ Double_t x = g->GetX()[i];
+ Double_t exl = g->GetEXlow()[i];
+ Double_t exh = g->GetEXhigh()[i];
+ bins.fArray[i] = x-exl;
+ bins.fArray[i+1] = x+exh;
+ Double_t dxi = exh+exl;
+ if (i == 0) dx = dxi;
+ else if (dxi != dx) dx = 0;
+ }
+ TString name(g->GetName());
+ TString title(g->GetTitle());
+ TH1D* h = 0;
+ if (dx != 0) {
+ h = new TH1D(name.Data(), title.Data(), nBins, bins[0], bins[nBins]);
+ }
+ else {
+ h = new TH1D(name.Data(), title.Data(), nBins, bins.fArray);
+ }
+ h->SetMarkerStyle(g->GetMarkerStyle());
+ h->SetMarkerColor(g->GetMarkerColor());
+ h->SetMarkerSize(g->GetMarkerSize());
+
+ return h;
+ }
+ //__________________________________________________________________
+ TH1* Asymmetry(TH1* h, Double_t& max)
+ {
+ if (!h) return 0;
+
+ TH1* ret = static_cast<TH1*>(h->Clone(Form("%s_leftright", h->GetName())));
+ // Int_t oBins = h->GetNbinsX();
+ // Double_t high = h->GetXaxis()->GetXmax();
+ // Double_t low = h->GetXaxis()->GetXmin();
+ // Double_t dBin = (high - low) / oBins;
+ // Int_t tBins = Int_t(2*high/dBin+.5);
+ // ret->SetBins(tBins, -high, high);
+ ret->Reset();
+ ret->SetTitle(Form("%s (+/-)", h->GetTitle()));
+ ret->SetYTitle("Right/Left");
+ Int_t nBins = h->GetNbinsX();
+ for (Int_t i = 1; i <= nBins; i++) {
+ Double_t x = h->GetBinCenter(i);
+ if (x > 0) break;
+
+ Double_t c1 = h->GetBinContent(i);
+ Double_t e1 = h->GetBinError(i);
+ if (c1 <= 0) continue;
+
+ Int_t j = h->FindBin(-x);
+ if (j <= 0 || j > nBins) continue;
+
+ Double_t c2 = h->GetBinContent(j);
+ Double_t e2 = h->GetBinError(j);
+
+ Double_t c12 = c1*c1;
+ Double_t e = TMath::Sqrt((e2*e2*c1*c1+e1*e1*c2*c2)/(c12*c12));
+
+ Int_t k = ret->FindBin(x);
+ ret->SetBinContent(k, c2/c1);
+ ret->SetBinError(k, e);
+ }
+ max = TMath::Max(max, RatioMax(ret));
+
+ return ret;
+ }
+
+
+ //__________________________________________________________________
+ Bool_t fShowOthers;
+ Bool_t fShowAlice;
+ Bool_t fShowRatios;
+ Bool_t fShowLeftRight;
+ UShort_t fRebin;
+ Bool_t fCutEdges;
+ TString fTitle;
+ TString fHHDFile;
+ TNamed* fTrigString;
+ TNamed* fSNNString;
+ TNamed* fSysString;
+ TAxis* fVtxAxis;
+ TH1* fForward;
+ TH1* fForwardMC;
+ TH1* fForwardHHD;
+ TH1* fTruth;
+ TH1* fCentral;
+ TH1* fForwardSym;
+ TH1* fForwardMCSym;
+ TH1* fForwardHHDSym;
+ TH1* fTriggers;
+ RangeParam* fRangeParam;
+
+};
+
+//=== Stuff for auto zooming =========================================
+void UpdateRange(dNdetaDrawer::RangeParam* p)
+{
+ if (!p) {
+ Warning("UpdateRange", "No parameters %p", p);
+ return;
+ }
+ if (!p->fMasterAxis) {
+ Warning("UpdateRange", "No master axis %p", p->fMasterAxis);
+ return;
+ }
+ Int_t first = p->fMasterAxis->GetFirst();
+ Int_t last = p->fMasterAxis->GetLast();
+ Double_t x1 = p->fMasterAxis->GetBinCenter(first);
+ Double_t x2 = p->fMasterAxis->GetBinCenter(last);
+ //Info("UpdateRange", "Range set to [%3d,%3d]->[%f,%f]", first, last, x1,x2);
+
+ if (p->fSlave1Axis) {
+ Int_t i1 = p->fSlave1Axis->FindBin(x1);
+ Int_t i2 = p->fSlave1Axis->FindBin(x2);
+ p->fSlave1Axis->SetRange(i1, i2);
+ p->fSlave1Pad->Modified();
+ p->fSlave1Pad->Update();
+ }
+ if (p->fSlave2Axis) {
+ Int_t i1 = p->fSlave2Axis->FindBin(x1);
+ Int_t i2 = p->fSlave2Axis->FindBin(x2);
+ p->fSlave2Axis->SetRange(i1, i2);
+ p->fSlave2Pad->Modified();
+ p->fSlave2Pad->Update();
+ }
+ TCanvas* c = gPad->GetCanvas();
+ c->cd();
+}
+
+//____________________________________________________________________
+void RangeExec(dNdetaDrawer::RangeParam* p)
+{
+ // Event types:
+ // 51: Mouse move
+ // 53:
+ // 1: Button down
+ // 21: Mouse drag
+ // 11: Mouse release
+ // dNdetaDrawer::RangeParam* p =
+ // reinterpret_cast<dNdetaDrawer::RangeParam*>(addr);
+ Int_t event = gPad->GetEvent();
+ TObject *select = gPad->GetSelected();
+ if (event == 53) {
+ UpdateRange(p);
+ return;
+ }
+ if (event != 11 || !select || select != p->fMasterAxis) return;
+ UpdateRange(p);
+}
+
+//=== Steering function ==============================================
+void
+DrawdNdeta(const char* filename="forward_dndeta.root",
+ Int_t flags=0xf,
+ const char* title=""
+ UShort_t rebin=5,
+ Bool_t cutEdges=false)
+{
+ dNdetaDrawer d;
+ d.SetRebin(rebin);
+ d.SetCutEdges(cutEdges);
+ d.SetTitle(title);
+ d.SetHHDFile("");
+ d.SetShowOthers(flags & 0x1);
+ d.SetShowAlice(flags & 0x2);
+ d.SetShowRatios(flags & 0x4);
+ d.SetShowLeftRight(flags & 0x8);
+ d.Run(filename);
+}
}
// --- Our data chain ----------------------------------------------
- gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/MakeESDChain.C");
- TChain* chain = MakeESDChain(esddir,true);
+ gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/MakeChain.C");
+ TChain* chain = MakeChain("ESD", esddir,true);
// If 0 or less events is select, choose all
if (nEvents <= 0) nEvents = chain->GetEntries();
// Central multiplicity
- Compile("$ALICE_ROOT/PWG2/FORWARD/analysis2/AddTaskCentralMult.C","");
- AddTaskCentralMult();
+ // Compile("$ALICE_ROOT/PWG2/FORWARD/analysis2/AddTaskCentralMult.C","");
+ // AddTaskCentralMult();
#endif
// FMD
--- /dev/null
+/**
+ * @file
+ *
+ * @ingroup pwg2_forward_scripts
+ */
+
+/**
+ * Run first pass of the analysis - that is read in ESD and produce AOD
+ *
+ * @param esddir ESD input directory. Any file matching the pattern
+ * *AliESDs*.root are added to the chain
+ * @param nEvents Number of events to process. If 0 or less, then
+ * all events are analysed
+ * @param proof Run in proof mode
+ * @param mc Run over MC data
+ *
+ * If PROOF mode is selected, then Terminate will be run on the master node
+ * in any case.
+ *
+ *
+ * @ingroup pwg2_forward_scripts
+ */
+void MakedNdeta(const char* aoddir=".",
+ Int_t nEvents=-1,
+ const char* trig="INEL",
+ Double_t vzMin=-10,
+ Double_t vzMax=10,
+ Int_t proof=0)
+{
+ // --- Libraries to load -------------------------------------------
+ gROOT->Macro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C");
+
+ // --- Check for proof mode, and possibly upload pars --------------
+ if (proof> 0) {
+ gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadPars.C");
+ LoadPars(proof);
+ }
+
+ // --- Our data chain ----------------------------------------------
+ gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/MakeChain.C");
+ TChain* chain = MakeChain("AOD", aoddir,true);
+ // If 0 or less events is select, choose all
+ if (nEvents <= 0) nEvents = chain->GetEntries();
+
+ // --- Creating the manager and handlers ---------------------------
+ AliAnalysisManager *mgr = new AliAnalysisManager("Forward Train",
+ "Forward dN/deta");
+ AliAnalysisManager::SetCommonFileName("forward_dndeta.root");
+
+ // --- ESD input handler -------------------------------------------
+ AliAODInputHandler *aodInputHandler = new AliAODInputHandler();
+ mgr->SetInputEventHandler(aodInputHandler);
+
+ // --- Add tasks ---------------------------------------------------
+ // Centrality
+ gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/AddTaskForwarddNdeta.C");
+ AddTaskForwarddNdeta(trig, vzMin, vzMax);
+
+ // --- Run the analysis --------------------------------------------
+ TStopwatch t;
+ if (!mgr->InitAnalysis()) {
+ Error("MakedNdeta", "Failed to initialize analysis train!");
+ return;
+ }
+ // Skip terminate if we're so requested and not in Proof or full mode
+ mgr->SetSkipTerminate(false);
+ // Some informative output
+ mgr->PrintStatus();
+ // mgr->SetDebugLevel(3);
+ if (mgr->GetDebugLevel() < 1 && !proof)
+ mgr->SetUseProgressBar(kTRUE);
+
+ // Run the train
+ t.Start();
+ Printf("=== RUNNING ANALYSIS ==================================");
+ mgr->StartAnalysis(proof ? "proof" : "local", chain, nEvents);
+ t.Stop();
+ t.Print();
+}
+//
+// EOF
+//
0.07, 0.07, 0.07, 0.07, 0.08, 0.08, 0.09, 0.09, 0.1, 0.13 };
const int np = 20;
double xm[np];
- for (int i = 0; i < np; i++) xm[i] = -x[i];
+ double exmm[np];
+ double expm[np];
+ double ym[np];
+ double eymm[np];
+ double eypm[np];
+ for (int i = 0; i < np; i++) {
+ int j = np-1-i;
+ xm[i] = -x[j];
+ exmm[i] = exm[j];
+ expm[i] = exp[j];
+ ym[i] = y[j];
+ eymm[i] = eym[j];
+ eypm[i] = eyp[j];
+ }
- TGraphAsymmErrors* g = new TGraphAsymmErrors(19,x, y,exm,exp,eym,eyp);
- TGraphAsymmErrors* gm = new TGraphAsymmErrors(19,xm,y,exm,exp,eym,eyp);
+ TGraphAsymmErrors* g = new TGraphAsymmErrors(19,x, y, exm, exp, eym, eyp);
+ TGraphAsymmErrors* gm = new TGraphAsymmErrors(19,xm,ym,exmm,expm,eymm,eypm);
SetGraphAttributes(g, NSDStyle, UA5Color,"ua5_nsd", "UA5 NSD");
SetGraphAttributes(gm, NSDStyle+MirrorOff, UA5Color,"ua5_nsd_mirrored",
"UA5 NSD (mirrored)");
0.07, 0.07, 0.07, 0.07, 0.08, 0.08, 0.09, 0.09, 0.1, 0.13 };
const int np = 19;
double xm[np];
- for (int i = 0; i < np; i++) xm[i] = -x[i];
- TGraphAsymmErrors* g = new TGraphAsymmErrors(np,x, y,exm,exp,eym,eyp);
- TGraphAsymmErrors* gm = new TGraphAsymmErrors(np,xm,y,exm,exp,eym,eyp);
+ double exmm[np];
+ double expm[np];
+ double ym[np];
+ double eymm[np];
+ double eypm[np];
+ for (int i = 0; i < np; i++) {
+ int j = np-1-i;
+ xm[i] = -x[j];
+ exmm[i] = exm[j];
+ expm[i] = exp[j];
+ ym[i] = y[j];
+ eymm[i] = eym[j];
+ eypm[i] = eyp[j];
+ }
+ TGraphAsymmErrors* g = new TGraphAsymmErrors(np,x, y, exm, exp, eym, eyp);
+ TGraphAsymmErrors* gm = new TGraphAsymmErrors(np,xm,ym,exmm,expm,eymm,eypm);
SetGraphAttributes(g, INELStyle, UA5Color, "ua5_inel", "UA5 INEL");
SetGraphAttributes(gm, INELStyle+MirrorOff, UA5Color, "ua5_inel_mirrored",
* @ingroup pwg2_forward_otherdata
*/
TMultiGraph*
-GetData(Int_t energy, Int_t type)
+GetData(Int_t energy, Int_t type, bool aliceOnly=false)
{
TMultiGraph* mp = new TMultiGraph(Form("dndeta_%dGeV_%d", energy, type),"");
TString tn;
en.Append(" #sqrt{s}=900GeV");
if (type & 0x1) {
tn.Append(" INEL");
- mp->Add(UA5Inel(false));
- mp->Add(UA5Inel(true));
+ if (!aliceOnly) mp->Add(UA5Inel(false));
+ if (!aliceOnly) mp->Add(UA5Inel(true));
mp->Add(AliceCentralInel900());
}
if (type & 0x4) {
tn.Append(" NSD");
- mp->Add(UA5Nsd(false));
- mp->Add(UA5Nsd(true));
+ if (!aliceOnly) mp->Add(UA5Nsd(false));
+ if (!aliceOnly) mp->Add(UA5Nsd(true));
mp->Add(AliceCentralNsd900());
- mp->Add(CMSNsd900());
+ if (!aliceOnly) mp->Add(CMSNsd900());
}
if (type & 0x2) {
tn.Append(" INEL>0");
if (type & 0x4) {
tn.Append(" NSD");
mp->Add(AliceCentralNsd2360());
- mp->Add(CMSNsd2360());
+ if (!aliceOnly) mp->Add(CMSNsd2360());
}
if (type & 0x1) {
tn.Append(" INEL>0");
}
if (type & 0x4) {
tn.Append(" NSD");
- mp->Add(CMSNsd7000());
+ if (!aliceOnly) mp->Add(CMSNsd7000());
}
if (type & 0x1) {
tn.Append(" INEL>0");
* @ingroup pwg2_forward_otherdata
*/
void
-OtherData(Int_t energy=900, Int_t type=0x1)
+OtherData(Int_t energy=900, Int_t type=0x1, bool aliceOnly=false)
{
- TMultiGraph* mp = GetData(energy, type);
+ TMultiGraph* mp = GetData(energy, type, aliceOnly);
if (!mp) return;
gStyle->SetTitleX(0.1);
* @ingroup pwg2_forward_scripts
*/
void
-Pass2(const char* file=".",
+Pass2(const char* aoddir=".",
+ Int_t nEvents=-1,
const char* triggers="INEL",
- Int_t energy=900,
Double_t vzMin=-10,
- Double_t vzMax=10,
- Int_t rebin=5,
- const char* title="",
- bool hhd=true,
- bool comp=true)
+ Double_t vzMax=10,
+ Int_t proof=0)
{
- gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/Compile.C");
- Compile("$ALICE_ROOT/PWG2/FORWARD/analysis2/AnalyseAOD.C","+g");
-
- Int_t trgMask;
- TString trgs(triggers);
- trgs.ToUpper();
- TObjString* trg;
- TIter next(trgs.Tokenize(" ,|"));
- while ((trg = static_cast<TObjString*>(next()))) {
- TString s(trg->GetString());
- if (s.IsNull()) continue;
- if (s.CompareTo("INEL") == 0) trgMask = AliAODForwardMult::kInel;
- else if (s.CompareTo("INEL>0")== 0) trgMask = AliAODForwardMult::kInelGt0;
- else if (s.CompareTo("NSD") == 0) trgMask = AliAODForwardMult::kNSD;
- else
- Warning("Pass2", "Unknown trigger %s", s.Data());
- }
- if (trgMask == 0) {
- trgMask = 1;
- trgs.Append("INEL");
- }
- TString tit(title);
- tit.ReplaceAll("@", " ");
+ gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/MakedNdeta.C");
- printf("--------------------------------------\n"
- "Settings for this:\n"
- " Input AOD: %s\n"
- " Vertex range: %+4.1f -> %+4.1f cm\n"
- " Rebinning: %d\n"
- " Trigger mask: 0x%02x (%s)\n"
- " Energy: %dGeV\n"
- " Title: %s\n"
- " HHD comp.: %s\n"
- " Other comp.: %s\n"
- "--------------------------------------\n",
- file, vzMin, vzMax, rebin, trgMask, trgs.Data(), energy, tit.Data(),
- hhd ? "yes" : "no", comp ? "yes" : "no");
-
- AnalyseAOD dr;
- TStopwatch t;
- t.Start();
- dr.Run(file, vzMin, vzMax, rebin, trgMask, energy, tit.Data(), hhd, comp);
- t.Stop();
- t.Print();
+ MakedNdeta(aoddir, nEvents, triggers, vzMin, vzMax, proof);
}
//
// EOF
-n,--events N Number of events ($nev)
-1,--pass1 Run only pass 1, no draw ($nodraw)
-2,--pass2 Run only pass 2, just draw ($noanal)
- -r,--rebin N Rebin by N ($rebin)
-v,--vz-min CM Minimum value of vz ($vzmin)
-V,--vz-max CM Maximum value of vz ($vzmax)
-t,--trigger TYPE Select trigger TYPE ($type)
- -e,--energy CMS Center of mass energy ($cms)
-b,--batch Do batch processing ($batch)
-P,--proof NWORKERS Run in PROOF(Lite) mode ($proof)
-M,--mc Run over MC data ($mc)
- -S,--title STRING Set the title string ($tit)
-g,--gdb Run in GDB mode ($gdb)
- -H,--hhd Do comparison to HHD ($hhd)
- -O,--other Do comparison to other ($comp)
-E,--eloss Run energy loss script
TYPE is a comma or space separated list of
-P|--proof) proof=$2 ; shift ;;
-M|--mc) mc=`toggle $mc` ;;
-g|--gdb) gdb=`toggle $gdb` ;;
- -H|--hhd) hhd=`toggle $hhd` ;;
- -O|--other) comp=`toggle $comp` ;;
- -r|--rebin) rebin=$2 ; shift ;;
-v|--vz-min) vzmin=$2 ; shift ;;
-V|--vz-max) vzmax=$2 ; shift ;;
- -e|--energy) cms=$2 ; shift ;;
- -S|--title) tit="$2" ; shift ;;
-E|--eloss) pass1=MakeELossFits.C ;
output=forward_eloss.root ;
nodraw=1 ;;
aliroot $opts $opts1 ${ana}/${pass1}\(\".\",$nev,$proof,$mc\)
fi
fail=$?
- rm -f event_stat.root \
- EventStat_temp.root \
- outputs_valid \
+ rm -f outputs_valid \
`printf %09d.stat $nev`
if test $fail -gt 0 ; then
echo "Return value $fail not 0" ; exit $fail
if test $nodraw -lt 1 ; then
rm -f result.root
tit=`echo $tit | tr ' ' '@'`
- echo "Running aliroot ${opts} ${opts1} ${ana}/Pass2.C\(\".\",\"$type\",$cms,$vzmin,$vzmax,$rebin,$hhd,$comp\)"
- aliroot ${opts} ${ana}/Pass2.C\(\".\",\"$type\",$cms,$vzmin,$vzmax,$rebin,\"$tit\",$hhd,$comp\)
+ echo "Running aliroot ${opts} ${opts1} ${ana}/Pass2.C\(\".\",$nev,\"$type\",$vzmin,$vzmax,$proof\)"
+ aliroot ${opts} ${ana}/Pass2.C\(\".\",$nev,\"$type\",$vzmin,$vzmax,$proof\)
fi
--- /dev/null
+/**
+ * Scan a directory (optionally recursive) for data files to add to
+ * the chain. Only ROOT files, and files which name contain the
+ * passed pattern are considered.
+ *
+ * @param dir Directory to scan
+ * @param chain Chain to add data to
+ * @param pattern Pattern that the file name must contain
+ * @param recursive Whether to scan recursively
+ */
+void
+ScanDirectory(TSystemDirectory* dir, TChain* chain,
+ const char* pattern, bool recursive)
+{
+ // Get list of files, and go back to old working directory
+ TString oldDir(gSystem->WorkingDirectory());
+ TList* files = dir->GetListOfFiles();
+ gSystem->ChangeDirectory(oldDir);
+
+ // Sort list of files and check if we should add it
+ files->Sort();
+ TIter next(files);
+ TSystemFile* file = 0;
+ while ((file = static_cast<TSystemFile*>(next()))) {
+ TString name(file->GetName());
+
+ // Ignore special links
+ if (name == "." || name == "..") continue;
+
+ // Check if this is a directory
+ if (file->IsDirectory()) {
+ if (recursive)
+ ScanDirectory(static_cast<TSystemDirectory*>(file),chain,recursive);
+ continue;
+ }
+
+ // If this is not a root file, ignore
+ if (!name.EndsWith(".root")) continue;
+
+ // If this file does not contain the pattern, ignore
+ if (!name.Contains(pattern)) continue;
+
+ // Get the path
+ TString data(Form("%s/%s", file->GetTitle(), name.Data()));
+
+ chain->Add(data);
+
+ }
+}
+
+/**
+ * Make a chain of specified data
+ *
+ * @param what What data to chain. Possible values are
+ * - ESD Event summary data (AliESD)
+ * - AOD Analysis object data (AliAOD)
+ * - MC Simulation data (galice)
+ * @param datadir Data directory to scan
+ * @param recursive Whether to recurse into sub-directories
+ *
+ * @return Pointer to newly create chain, or null
+ */
+TChain*
+MakeChain(const char* what, const char* datadir, bool recursive=false)
+{
+ TString w(what);
+ w.ToUpper();
+ const char* treeName = 0;
+ const char* pattern = 0;
+ if (w.Contains("ESD")) { treeName = "esdTree"; pattern = "AliESD"; }
+ else if (w.Contains("AOD")) { treeName = "aodTree"; pattern = "AliAOD"; }
+ else if (w.Contains("MC")) { treeName = "TE"; pattern = "galice"; }
+ else {
+ Error("MakeChain", "Unknown mode '%s' (not one of ESD,AOD, or MC)", what);
+ return 0;
+ }
+
+ // --- Our data chain ----------------------------------------------
+ TChain* chain = new TChain(treeName);
+
+ // --- Get list of ESDs --------------------------------------------
+ // Open source directory, and make sure we go back to were we were
+ TString oldDir(gSystem->WorkingDirectory());
+ TSystemDirectory d(datadir, datadir);
+ ScanDirectory(&d, chain, pattern, recursive);
+
+ return chain;
+}
+++ /dev/null
-void
-ScanDirectory(TSystemDirectory* dir, TChain* chain, bool recursive)
-{
- // gROOT->IndentLevel();
- // Printf("Scanning %s ...", dir->GetName());
- // gROOT->IncreaseDirLevel();
-
-
- // Get list of files, and go back to old working directory
- TString oldDir(gSystem->WorkingDirectory());
- TList* files = dir->GetListOfFiles();
- gSystem->ChangeDirectory(oldDir);
-
- // Sort list of files and check if we should add it
- files->Sort();
- TIter next(files);
- TSystemFile* file = 0;
- while ((file = static_cast<TSystemFile*>(next()))) {
- TString name(file->GetName());
-
- // Ignore special links
- if (name == "." || name == "..") continue;
-
- // Check if this is a directory
- if (file->IsDirectory()) {
- if (recursive)
- ScanDirectory(static_cast<TSystemDirectory*>(file),chain,recursive);
- continue;
- }
-
- // If this is not a root file, ignore
- if (!name.EndsWith(".root")) continue;
-
- // If this file does not contain AliESDs, ignore
- if (!name.Contains("AliESDs")) continue;
-
- // Get the path
- TString esd(Form("%s/%s", file->GetTitle(), name.Data()));
-
- // Print and add
- // gROOT->IndentLevel();
- // Printf("adding %s", esd.Data());
- chain->Add(esd);
-
- }
- // gROOT->DecreaseDirLevel();
-}
-
-TChain*
-MakeESDChain(const char* esddir, bool recursive=false)
-{
- // --- Our data chain ----------------------------------------------
- TChain* chain = new TChain("esdTree");
-
- // --- Get list of ESDs --------------------------------------------
- // Open source directory, and make sure we go back to were we were
- TString oldDir(gSystem->WorkingDirectory());
- TSystemDirectory d(esddir, esddir);
- ScanDirectory(&d, chain, recursive);
-
- // chain->GetListOfFiles()->ls();
-
- return chain;
-}
// Note: custom streamer to ensure singleton consistency!
#pragma link C++ class AliForwardCorrectionManager-;
#pragma link C++ class AliForwardMCCorrectionsTask+;
+#pragma link C++ class AliForwarddNdetaTask+;
#else
# error Not for compilation
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff