--- /dev/null
+/**
+ * This is the macro to include the Forward multiplicity in a train.
+ *
+ * @ingroup pwg2_forward_analysis_scripts
+ */
+AliForwardMultiplicity*
+AddTaskFMD(Int_t nCutBins=1, Float_t correctionCut=0.1)
+{
+ AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+ if (!mgr) {
+ Error("AddTaskFMD", "No analysis manager to connect to.");
+ return NULL;
+ }
+
+ AliForwardMultiplicity* task = new AliForwardMultiplicity("FMD");
+ task->GetHistCollector().SetNCutBins(nCutBins);
+ task->GetHistCollector().SetCorrectionCut(correctionCut);
+ mgr->AddTask(task);
+
+ AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+ AliMCEventHandler* mcHandler =
+ dynamic_cast<AliMCEventHandler*>(mgr->GetMCtruthEventHandler());
+ Info("AddTaskFMD", "MC handler %p", mcHandler);
+ if(mcHandler) {
+ pars->SetRealData(kFALSE);
+ pars->SetProcessPrimary(kTRUE);
+ pars->SetProcessHits(kFALSE);
+ }
+ else {
+ pars->SetRealData(kTRUE);
+ pars->SetProcessPrimary(kFALSE);
+ pars->SetProcessHits(kFALSE);
+ }
+ pars->Init();
+
+ TString outputfile = AliAnalysisManager::GetCommonFileName();
+ outputfile += Form(":%s",pars->GetDndetaAnalysisName());
+ AliAnalysisDataContainer* histOut =
+ mgr->CreateContainer("Forward", TList::Class(),
+ AliAnalysisManager::kOutputContainer,outputfile);
+
+
+ mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
+ mgr->ConnectOutput(task, 1, histOut);
+
+ return task;
+}
--- /dev/null
+#include "AliAODForwardMult.h"
+#include <TBrowser.h>
+#include <iostream>
+#include <TMath.h>
+#include <TObjString.h>
+
+ClassImp(AliAODForwardMult)
+#if 0
+; // For Emacs
+#endif
+
+//____________________________________________________________________
+const Float_t AliAODForwardMult::fgkInvalidIpZ = 1e6;
+
+//____________________________________________________________________
+AliAODForwardMult::AliAODForwardMult()
+ : fHist(),
+ fTriggers(0),
+ fIpZ(fgkInvalidIpZ)
+{}
+
+//____________________________________________________________________
+AliAODForwardMult::AliAODForwardMult(Bool_t)
+ : fHist("forwardMult", "d^{2}N_{ch}/d#etad#varphi in the forward regions",
+ 200, -4, 6, 20, 0, 2*TMath::Pi()),
+ fTriggers(0),
+ fIpZ(fgkInvalidIpZ)
+{
+ fHist.SetXTitle("#eta");
+ fHist.SetYTitle("#varphi [radians]");
+ fHist.SetZTitle("#frac{d^{2}N_{ch}}{d#etad#varphi}");
+ fHist.Sumw2();
+}
+
+//____________________________________________________________________
+void
+AliAODForwardMult::Init(const TAxis& etaAxis)
+{
+ fHist.SetBins(etaAxis.GetNbins(), etaAxis.GetXmin(), etaAxis.GetXmax(),
+ 20, 0, 2*TMath::Pi());
+}
+
+//____________________________________________________________________
+void
+AliAODForwardMult::Clear(Option_t* option)
+{
+ fHist.Reset(option);
+ fTriggers = 0;
+ fIpZ = fgkInvalidIpZ;
+}
+//____________________________________________________________________
+Bool_t
+AliAODForwardMult::HasIpZ() const
+{
+ return TMath::Abs(fIpZ - fgkInvalidIpZ) > 1;
+}
+
+//____________________________________________________________________
+void
+AliAODForwardMult::Browse(TBrowser* b)
+{
+ static TObjString ipz;
+ static TObjString trg;
+ ipz = Form("ip_z=%fcm", fIpZ);
+ trg = GetTriggerString(fTriggers);
+ b->Add(&fHist);
+ b->Add(&ipz);
+ b->Add(&trg);
+}
+
+//____________________________________________________________________
+const Char_t*
+AliAODForwardMult::GetTriggerString(UInt_t mask)
+{
+ static TString trg;
+ trg = "";
+ if ((mask & kInel) != 0x0) trg.Append("INEL ");
+ if ((mask & kInelGt0) != 0x0) trg.Append("INEL>0 ");
+ if ((mask & kNSD) != 0x0) trg.Append("NSD ");
+ if ((mask & kA) != 0x0) trg.Append("A ");
+ if ((mask & kB) != 0x0) trg.Append("B ");
+ if ((mask & kC) != 0x0) trg.Append("C ");
+ if ((mask & kE) != 0x0) trg.Append("E ");
+ return trg.Data();
+}
+
+//____________________________________________________________________
+void
+AliAODForwardMult::Print(Option_t* option) const
+{
+ fHist.Print(option);
+ std::cout << "Ipz: " << fIpZ << "cm " << (HasIpZ() ? "" : "in")
+ << "valid\n"
+ << "Triggers: " << GetTriggerString(fTriggers) << std::endl;
+}
+
+//____________________________________________________________________
+//
+// EOF
+//
--- /dev/null
+#ifndef ALIROOT_PWG2_FORWARD_ANALYSIS_ALIAODFORWARDMULT_H
+#define ALIROOT_PWG2_FORWARD_ANALYSIS_ALIAODFORWARDMULT_H
+#include <TObject.h>
+#include <TH2D.h>
+class TBrowser;
+/**
+ * Class that contains the forward multiplicity data per event
+ *
+ * This class contains a histogram of
+ * @f[
+ * \frac{d^2N_{ch}}{d\eta d\phi}\quad,
+ * @f]
+ * as well as a trigger mask for each analysed event.
+ *
+ * The eta acceptance of the event is stored in the underflow bins of
+ * the histogram. So to build the final histogram, one needs to
+ * correct for this acceptance (properly weighted by the events), and
+ * the vertex efficiency. This simply boils down to defining a 2D
+ * histogram and summing the event histograms in that histogram. One
+ * should of course also do proper book-keeping of the accepted event.
+ *
+ * @code
+ * TH2D* sum = 0; // Summed hist
+ * TTree* tree = GetAODTree(); // AOD tree
+ * AliAODForwardMult* mult = 0; // AOD object
+ * Int_t nTriggered = 0; // # of triggered ev.
+ * Int_t nAccepted = 0; // # of ev. w/vertex
+ * Int_t nAvailable = tree->GetEntries(); // How many entries
+ * Float_t vzLow = -10; // Lower ip cut
+ * Float_t vzHigh = 10; // Upper ip cut
+ * Int_t mask = AliAODForward::kINEL;// Trigger mask
+ * tree->SetBranchAddress("forward", &forward); // Set the address
+ *
+ * for (int i = 0; i < nAvailable; i++) {
+ * // Create sum histogram on first event - to match binning to input
+ * if (!sum) sum = static_cast<TH2D*>(mult->Clone("d2ndetadphi"));
+ *
+ * tree->GetEntry(i);
+ *
+ * // Other trigger/event requirements could be defined
+ * if (!mult->IsTriggerBits(mask)) continue;
+ * nTriggered++;
+ *
+ * // Select vertex range (in centimeters)
+ * if (!mult->InRange(vzLow, vzHigh) continue;
+ * nAccepted++;
+ *
+ * // Add contribution from this event
+ * sum->Add(&(mult->GetHistogram()));
+ * }
+ *
+ * // Get acceptance normalisation from underflow bins
+ * TH1D* norm = sum->Projection("norm", 0, 1, "");
+ * // Project onto eta axis - _ignoring_underflow_bins_!
+ * TH1D* dndeta = sum->Projection("dndeta", 1, -1, "e");
+ * // Normalize to the acceptance
+ * dndeta->Divide(norm);
+ * // Scale by the vertex efficiency
+ * dndeta->Scale(Double_t(nAccepted)/nTriggered, "width");
+ * // And draw the result
+ * dndeta->Draw();
+ * @endcode
+ *
+ * The above code will draw the final @f$ dN_{ch}/d\eta@f$ for the
+ * selected event class and vertex range
+ *
+ * The histogram can be used as input for other kinds of analysis too,
+ * like flow, event-plane, centrality, and so on.
+ *
+ * @ingroup pwg2_forward_analysis
+ */
+class AliAODForwardMult : public TObject
+{
+public:
+ /**
+ * Bits of the trigger pattern
+ */
+ enum {
+ /** In-elastic collision */
+ kInel = 0x001,
+ /** In-elastic collision with at least one SPD tracklet */
+ kInelGt0 = 0x002,
+ /** Non-single diffractive collision */
+ kNSD = 0x004,
+ /** Empty bunch crossing */
+ kEmpty = 0x008,
+ /** A-side trigger */
+ kA = 0x010,
+ /** B(arrel) trigger */
+ kB = 0x020,
+ /** C-side trigger */
+ kC = 0x080,
+ /** Empty trigger */
+ kE = 0x100
+ };
+ /**
+ * Default constructor
+ *
+ * Used by ROOT I/O sub-system - do not use
+ */
+ AliAODForwardMult();
+ /**
+ * Constructor
+ *
+ */
+ AliAODForwardMult(Bool_t);
+ /**
+ * Destructor
+ */
+ ~AliAODForwardMult() {}
+ /**
+ * Initialize
+ *
+ * @param etaAxis Pseudo-rapidity axis
+ */
+ void Init(const TAxis& etaAxis);
+ /**
+ * Get the @f$ d^2N_{ch}/d\eta d\phi@f$ histogram,
+ *
+ * @return @f$ d^2N_{ch}/d\eta d\phi@f$ histogram,
+ */
+ const TH2D& GetHistogram() const { return fHist; }
+ /**
+ * Get the @f$ d^2N_{ch}/d\eta d\phi@f$ histogram,
+ *
+ * @return @f$ d^2N_{ch}/d\eta d\phi@f$ histogram,
+ */
+ TH2D& GetHistogram() { return fHist; }
+ /**
+ * Get the trigger mask
+ *
+ * @return Trigger mask
+ */
+ UInt_t GetTriggerMask() const { return fTriggers; }
+ /**
+ * Set the trigger mask
+ *
+ * @param trg Trigger mask
+ */
+ void SetTriggerMask(UInt_t trg) { fTriggers = trg; }
+ /**
+ * Set bit(s) in trigger mask
+ *
+ * @param bits bit(s) to set
+ */
+ void SetTriggerBits(UInt_t bits) { fTriggers |= bits; }
+ /**
+ * Check if bit(s) are set in the trigger mask
+ *
+ * @param bits Bits to test for
+ *
+ * @return
+ */
+ Bool_t IsTriggerBits(UInt_t bits) const;
+ /**
+ * Whether we have any trigger bits
+ */
+ Bool_t HasTrigger() const { return fTriggers != 0; }
+ /**
+ * 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; }
+ /**
+ * Print content
+ *
+ * @param option Passed verbatim to TH2::Print
+ */
+ void Print(Option_t* option="") const;
+ /**
+ * Set the z coordinate of the interaction point
+ *
+ * @param ipZ Interaction point z coordinate
+ */
+ void SetIpZ(Float_t ipZ) { fIpZ = ipZ; }
+ /**
+ * Set the z coordinate of the interaction point
+ *
+ * @return Interaction point z coordinate
+ */
+ Float_t GetIpZ() const { return fIpZ; }
+ /**
+ * Check if we have a valid z coordinate of the interaction point
+ *
+ * @return True if we have a valid interaction point z coordinate
+ */
+ Bool_t HasIpZ() const;
+ /**
+ * Check if the z coordinate of the interaction point is within the
+ * given limits. Note that the convention used corresponds to the
+ * convention used in ROOTs TAxis.
+ *
+ * @param low Lower cut (inclusive)
+ * @param high Upper cut (exclusive)
+ *
+ * @return true if @f$ low \ge ipz < high@f$
+ */
+ Bool_t InRange(Float_t low, Float_t high) const;
+ /**
+ * Get the name of the object
+ *
+ * @return Name of object
+ */
+ const Char_t* GetName() const { return "Forward"; }
+ /**
+ * Get a string correspondig to the trigger mask
+ *
+ * @param mask Trigger mask
+ *
+ * @return Static string (copy before use)
+ */
+ static const Char_t* GetTriggerString(UInt_t mask);
+protected:
+ TH2D fHist; // Histogram of d^2N_{ch}/(deta dphi) for this event
+ UInt_t fTriggers; // Trigger bit mask
+ Float_t fIpZ; // Z coordinate of the interaction point
+
+ static const Float_t fgkInvalidIpZ; // Invalid IpZ value
+ ClassDef(AliAODForwardMult,1); // AOD forward multiplicity
+};
+
+//____________________________________________________________________
+inline Bool_t
+AliAODForwardMult::InRange(Float_t low, Float_t high) const
+{
+ return HasIpZ() && fIpZ >= low && fIpZ < high;
+}
+
+//____________________________________________________________________
+inline Bool_t
+AliAODForwardMult::IsTriggerBits(UInt_t bits) const
+{
+ return HasTrigger() && ((fTriggers & bits) != 0);
+}
+
+
+#endif
+// Local Variables:
+// mode: C++
+// End:
+
--- /dev/null
+#include "AliFMDCorrections.h"
+#include <AliESDFMD.h>
+#include <TAxis.h>
+#include <TList.h>
+#include <TMath.h>
+#include "AliFMDAnaParameters.h"
+#include "AliLog.h"
+#include <TH2D.h>
+
+ClassImp(AliFMDCorrections)
+#if 0
+; // For Emacs
+#endif
+
+//____________________________________________________________________
+AliFMDCorrections::AliFMDCorrections()
+ : TNamed(),
+ fRingHistos(),
+ fMultCut(0.3)
+{}
+
+//____________________________________________________________________
+AliFMDCorrections::AliFMDCorrections(const char* title)
+ : TNamed("fmdCorrections", title),
+ fRingHistos(),
+ fMultCut(0.3)
+{
+ fRingHistos.SetName(GetName());
+ fRingHistos.Add(new RingHistos(1, 'I'));
+ fRingHistos.Add(new RingHistos(2, 'I'));
+ fRingHistos.Add(new RingHistos(2, 'O'));
+ fRingHistos.Add(new RingHistos(3, 'I'));
+ fRingHistos.Add(new RingHistos(3, 'O'));
+}
+
+//____________________________________________________________________
+AliFMDCorrections::AliFMDCorrections(const AliFMDCorrections& o)
+ : TNamed(o),
+ fRingHistos(),
+ fMultCut(o.fMultCut)
+{
+ TIter next(&o.fRingHistos);
+ TObject* obj = 0;
+ while ((obj = next())) fRingHistos.Add(obj);
+}
+
+//____________________________________________________________________
+AliFMDCorrections::~AliFMDCorrections()
+{
+ fRingHistos.Delete();
+}
+
+//____________________________________________________________________
+AliFMDCorrections&
+AliFMDCorrections::operator=(const AliFMDCorrections& o)
+{
+ SetName(o.GetName());
+ SetTitle(o.GetTitle());
+
+ fMultCut = o.fMultCut;
+
+ fRingHistos.Delete();
+ TIter next(&o.fRingHistos);
+ TObject* obj = 0;
+ while ((obj = next())) fRingHistos.Add(obj);
+
+ return *this;
+}
+
+//____________________________________________________________________
+AliFMDCorrections::RingHistos*
+AliFMDCorrections::GetRingHistos(UShort_t d, Char_t r) const
+{
+ Int_t idx = -1;
+ switch (d) {
+ case 1: idx = 0; break;
+ case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
+ case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
+ }
+ if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0;
+
+ return static_cast<RingHistos*>(fRingHistos.At(idx));
+}
+
+//____________________________________________________________________
+Bool_t
+AliFMDCorrections::Correct(AliForwardUtil::Histos& hists,
+ Int_t vtxbin)
+{
+ AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+
+ 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');
+ TH2D* h = hists.Get(d,r);
+ RingHistos* rh= GetRingHistos(d,r);
+ TH2F* bg= pars->GetBackgroundCorrection(d, r, vtxbin);
+ TH2F* ef= pars->GetEventSelectionEfficiency("INEL",vtxbin,r);
+ if (!bg) {
+ AliWarning(Form("No secondary correction for FMDM%d%c in vertex bin %d",
+ d, r, vtxbin));
+ continue;
+ }
+ if (!ef) {
+ AliWarning(Form("No event %s selection efficiency in vertex bin %d",
+ "INEL", vtxbin));
+ continue;
+ }
+
+ // Divide by primary/total ratio
+ h->Divide(bg);
+
+ // Divide by the event selection efficiency
+ h->Divide(ef);
+
+ if(!pars->SharingEffPresent()) {
+ AliWarning("No sharing efficiencies");
+ continue;
+ }
+
+ TH1F* sf = pars->GetSharingEfficiencyTrVtx(d,r,vtxbin);
+
+ for (Int_t ieta = 1; ieta <= h->GetNbinsX(); ieta++) {
+ Float_t c = sf->GetBinContent(ieta);
+ Float_t ec = sf->GetBinError(ieta);
+
+ if (c == 0) continue;
+
+ for (Int_t iphi = 1; iphi <= h->GetNbinsY(); iphi++) {
+ Double_t m = h->GetBinContent(ieta, iphi) / c;
+ Double_t em = h->GetBinError(ieta, iphi);
+
+ Double_t e = TMath::Sqrt(em * em + (m * ec) * (m * ec));
+
+ h->SetBinContent(ieta,iphi,m);
+ h->SetBinError(ieta,iphi,e);
+ }
+ }
+
+ rh->fDensity->Add(h);
+ }
+ }
+
+ return kTRUE;
+}
+
+//____________________________________________________________________
+void
+AliFMDCorrections::ScaleHistograms(Int_t nEvents)
+{
+ if (nEvents <= 0) return;
+
+ TIter next(&fRingHistos);
+ RingHistos* o = 0;
+ while ((o = static_cast<RingHistos*>(next()))) {
+ o->fDensity->Scale(1. / nEvents);
+ }
+}
+
+//____________________________________________________________________
+void
+AliFMDCorrections::Output(TList* dir)
+{
+ TList* d = new TList;
+ d->SetName(GetName());
+ dir->Add(d);
+ TIter next(&fRingHistos);
+ RingHistos* o = 0;
+ while ((o = static_cast<RingHistos*>(next()))) {
+ o->Output(d);
+ }
+}
+
+//====================================================================
+AliFMDCorrections::RingHistos::RingHistos()
+ : fDet(0),
+ fRing('\0'),
+ fDensity(0)
+{}
+
+//____________________________________________________________________
+AliFMDCorrections::RingHistos::RingHistos(UShort_t d, Char_t r)
+ : fDet(d),
+ fRing(r),
+ fDensity(0)
+{
+ fDensity = new TH2D(Form("FMD%d%c_Primary_Density", d, r),
+ Form("in FMD%d%c", d, r),
+ 200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40),
+ 0, 2*TMath::Pi());
+ fDensity->SetDirectory(0);
+ fDensity->SetXTitle("#eta");
+ fDensity->SetYTitle("#phi [radians]");
+ fDensity->SetZTitle("Primary N_{ch} density");
+}
+//____________________________________________________________________
+AliFMDCorrections::RingHistos::RingHistos(const RingHistos& o)
+ : TObject(o),
+ fDet(o.fDet),
+ fRing(o.fRing),
+ fDensity(o.fDensity)
+{}
+
+//____________________________________________________________________
+AliFMDCorrections::RingHistos&
+AliFMDCorrections::RingHistos::operator=(const RingHistos& o)
+{
+ fDet = o.fDet;
+ fRing = o.fRing;
+
+ if (fDensity) delete fDensity;
+
+ fDensity = static_cast<TH2D*>(o.fDensity->Clone());
+
+ return *this;
+}
+//____________________________________________________________________
+AliFMDCorrections::RingHistos::~RingHistos()
+{
+ if (fDensity) delete fDensity;
+}
+
+//____________________________________________________________________
+void
+AliFMDCorrections::RingHistos::Output(TList* dir)
+{
+ TList* d = new TList;
+ d->SetName(Form("FMD%d%c", fDet, fRing));
+ d->Add(fDensity);
+ dir->Add(d);
+}
+
+//____________________________________________________________________
+//
+// EOF
+//
+
+
+
--- /dev/null
+#ifndef ALIROOT_PWG2_FORWARD_ANALYSIS2_ALIFMDCORRECTIONS_H
+#define ALIROOT_PWG2_FORWARD_ANALYSIS2_ALIFMDCORRECTIONS_H
+#include <TNamed.h>
+#include <TList.h>
+#include "AliForwardUtil.h"
+class AliESDFMD;
+class TH2D;
+
+/**
+ * This class calculates the inclusive charged particle density
+ * in each for the 5 FMD rings.
+ *
+ * @par Input:
+ * - AliESDFMD object possibly corrected for sharing
+ *
+ * @par Output:
+ * - 5 RingHistos objects - each with a number of vertex dependent
+ * 2D histograms of the inclusive charge particle density
+ *
+ * @par Corrections used:
+ * - AliFMDAnaCalibBackgroundCorrection
+ * - AliFMDAnaCalibEventSelectionEfficiency
+ * - AliFMDAnaCalibSharingEfficiency
+ *
+ * @ingroup pwg2_forward_analysis
+ */
+class AliFMDCorrections : public TNamed
+{
+public:
+ /**
+ * Constructor
+ */
+ AliFMDCorrections();
+ /**
+ * Constructor
+ *
+ * @param name Name of object
+ */
+ AliFMDCorrections(const char* name);
+ /**
+ * Copy constructor
+ *
+ * @param o Object to copy from
+ */
+ AliFMDCorrections(const AliFMDCorrections& o);
+ /**
+ * Destructor
+ */
+ virtual ~AliFMDCorrections();
+ /**
+ * Assignement operator
+ *
+ * @param o Object to assign from
+ *
+ * @return Reference to this object
+ */
+ AliFMDCorrections& operator=(const AliFMDCorrections&);
+ /**
+ * Do the calculations
+ *
+ * @param hists Cache of histograms
+ * @param vtxBin Vertex bin
+ *
+ * @return true on successs
+ */
+ virtual Bool_t Correct(AliForwardUtil::Histos& hists, Int_t vtxBin);
+ /**
+ * Scale the histograms to the total number of events
+ *
+ * @param nEvents Number of events
+ */
+ void ScaleHistograms(Int_t nEvents);
+ /**
+ * Output diagnostic histograms to directory
+ *
+ * @param dir List to write in
+ */
+ void Output(TList* dir);
+protected:
+ /**
+ * Internal data structure to keep track of the histograms
+ */
+ struct RingHistos : public TObject
+ {
+ /**
+ * Default CTOR
+ */
+ RingHistos();
+ /**
+ * Constructor
+ *
+ * @param d detector
+ * @param r ring
+ */
+ RingHistos(UShort_t d, Char_t r);
+ /**
+ * Copy constructor
+ *
+ * @param o Object to copy from
+ */
+ RingHistos(const RingHistos& o);
+ /**
+ * Assignment operator
+ *
+ * @param o Object to assign from
+ *
+ * @return Reference to this
+ */
+ RingHistos& operator=(const RingHistos& o);
+ /**
+ * Destructor
+ */
+ ~RingHistos();
+ void Output(TList* dir);
+ UShort_t fDet; // Detector
+ Char_t fRing; // Ring
+ TH2D* fDensity; // Distribution primary Nch
+ };
+ /**
+ * Get the ring histogram container
+ *
+ * @param d Detector
+ * @param r Ring
+ *
+ * @return Ring histogram container
+ */
+ RingHistos* GetRingHistos(UShort_t d, Char_t r) const;
+
+ TList fRingHistos; // List of histogram containers
+ Double_t fMultCut; // Low cut on scaled energy loss
+
+ ClassDef(AliFMDCorrections,1); // Calculate Nch density
+};
+
+#endif
+// Local Variables:
+// mode: C++
+// End:
+
--- /dev/null
+#include "AliFMDDensityCalculator.h"
+#include <AliESDFMD.h>
+#include <TAxis.h>
+#include <TList.h>
+#include <TMath.h>
+#include "AliFMDAnaParameters.h"
+#include "AliLog.h"
+#include <TH2D.h>
+
+ClassImp(AliFMDDensityCalculator)
+#if 0
+; // For Emacs
+#endif
+
+//____________________________________________________________________
+AliFMDDensityCalculator::AliFMDDensityCalculator()
+ : TNamed(),
+ fRingHistos(),
+ fMultCut(0.3)
+{}
+
+//____________________________________________________________________
+AliFMDDensityCalculator::AliFMDDensityCalculator(const char* title)
+ : TNamed("fmdDensityCalculator", title),
+ fRingHistos(),
+ fMultCut(0.3)
+{
+ fRingHistos.SetName(GetName());
+ fRingHistos.Add(new RingHistos(1, 'I'));
+ fRingHistos.Add(new RingHistos(2, 'I'));
+ fRingHistos.Add(new RingHistos(2, 'O'));
+ fRingHistos.Add(new RingHistos(3, 'I'));
+ fRingHistos.Add(new RingHistos(3, 'O'));
+}
+
+//____________________________________________________________________
+AliFMDDensityCalculator::AliFMDDensityCalculator(const
+ AliFMDDensityCalculator& o)
+ : TNamed(o),
+ fRingHistos(),
+ fMultCut(o.fMultCut)
+{
+ TIter next(&o.fRingHistos);
+ TObject* obj = 0;
+ while ((obj = next())) fRingHistos.Add(obj);
+}
+
+//____________________________________________________________________
+AliFMDDensityCalculator::~AliFMDDensityCalculator()
+{
+ fRingHistos.Delete();
+}
+
+//____________________________________________________________________
+AliFMDDensityCalculator&
+AliFMDDensityCalculator::operator=(const AliFMDDensityCalculator& o)
+{
+ SetName(o.GetName());
+ SetTitle(o.GetTitle());
+
+ fMultCut = o.fMultCut;
+
+ fRingHistos.Delete();
+ TIter next(&o.fRingHistos);
+ TObject* obj = 0;
+ while ((obj = next())) fRingHistos.Add(obj);
+
+ return *this;
+}
+
+//____________________________________________________________________
+AliFMDDensityCalculator::RingHistos*
+AliFMDDensityCalculator::GetRingHistos(UShort_t d, Char_t r) const
+{
+ Int_t idx = -1;
+ switch (d) {
+ case 1: idx = 0; break;
+ case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
+ case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
+ }
+ if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0;
+
+ return static_cast<RingHistos*>(fRingHistos.At(idx));
+}
+
+//____________________________________________________________________
+Bool_t
+AliFMDDensityCalculator::Calculate(const AliESDFMD& fmd,
+ AliForwardUtil::Histos& hists,
+ Int_t vtxbin,
+ Bool_t lowFlux)
+{
+ 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');
+ UShort_t ns= (q == 0 ? 20 : 40);
+ UShort_t nt= (q == 0 ? 512 : 256);
+ TH2D* h = hists.Get(d,r);
+ RingHistos* rh= GetRingHistos(d,r);
+
+ for (UShort_t s=0; s<ns; s++) {
+ for (UShort_t t=0; t<nt; t++) {
+ Float_t mult = fmd.Multiplicity(d,r,s,t);
+
+ if (mult == 0 || mult > 20) continue;
+
+ Float_t phi = fmd.Phi(d,r,s,t) / 180 * TMath::Pi();
+ Float_t eta = fmd.Eta(d,r,s,t);
+
+ Float_t n = NParticles(mult,d,r,s,t,vtxbin,eta,lowFlux);
+ rh->fEvsN->Fill(mult,n);
+
+ Float_t c = Correction(d,r,s,t,vtxbin,eta,lowFlux);
+ if (c > 0) n /= c;
+ rh->fEvsM->Fill(mult,n);
+
+ h->Fill(eta,phi,n);
+ rh->fDensity->Fill(eta,phi,n);
+ } // for t
+ } // for s
+ } // for q
+ } // for d
+
+ return kTRUE;
+}
+
+//_____________________________________________________________________
+Float_t
+AliFMDDensityCalculator::NParticles(Float_t mult,
+ UShort_t d,
+ Char_t r,
+ UShort_t /*s*/,
+ UShort_t /*t*/,
+ Int_t /*v*/,
+ Float_t eta,
+ Bool_t lowFlux) const
+{
+ if (mult <= fMultCut) return 0;
+ if (lowFlux) return 1;
+
+ AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+
+ Float_t mpv = pars->GetMPV(d,r,eta);
+ Float_t w = pars->GetSigma(d,r,eta);
+ Float_t w2 = pars->Get2MIPWeight(d,r,eta);
+ Float_t w3 = pars->Get3MIPWeight(d,r,eta);
+ Float_t mpv2 = 2*mpv+2*w*TMath::Log(2);
+ Float_t mpv3 = 3*mpv+3*w*TMath::Log(3);
+
+ Float_t sum = (TMath::Landau(mult,mpv,w,kTRUE) +
+ w2 * TMath::Landau(mult,mpv2,2*w,kTRUE) +
+ w3 * TMath::Landau(mult,mpv3,3*w,kTRUE));
+ Float_t wsum = (TMath::Landau(mult,mpv,w,kTRUE) +
+ 2*w2 * TMath::Landau(mult,mpv2,2*w,kTRUE) +
+ 3*w3 * TMath::Landau(mult,mpv3,3*w,kTRUE));
+
+ return (sum > 0) ? wsum / sum : 1;
+}
+
+//_____________________________________________________________________
+Float_t
+AliFMDDensityCalculator::Correction(UShort_t d, Char_t r, UShort_t /*s*/,
+ UShort_t t, Int_t v, Float_t eta,
+ Bool_t lowFlux) const
+{
+ AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+
+ Float_t correction = AcceptanceCorrection(r,t);
+ if (lowFlux) {
+ TH1F* dblHitCor = pars->GetDoubleHitCorrection(d,r);
+ if (dblHitCor) {
+ Float_t dblC = dblHitCor->GetBinContent(dblHitCor->FindBin(eta));
+ if (dblC > 0)
+ correction *= dblC;
+ }
+ else
+ AliWarning(Form("Missing double hit correction for FMD%d%c",d,r));
+ }
+
+ TH1F* deadCor = pars->GetFMDDeadCorrection(v);
+ if (deadCor) {
+ Float_t deadC = deadCor->GetBinContent(deadCor->FindBin(eta));
+ if (deadC > 0)
+ correction *= deadC;
+ }
+
+ return correction;
+}
+
+
+//_____________________________________________________________________
+Float_t
+AliFMDDensityCalculator::AcceptanceCorrection(Char_t r, UShort_t t) const
+{
+ AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+
+ //AliFMDRing fmdring(ring);
+ //fmdring.Init();
+ Float_t rad = pars->GetMaxR(r)-pars->GetMinR(r);
+ Float_t slen = pars->GetStripLength(r,t);
+ Float_t sblen = pars->GetBaseStripLength(r,t);
+ Float_t nstrips = (r == 'I' ? 512 : 256);
+ Float_t segment = rad / nstrips;
+ Float_t radius = pars->GetMinR(r) + segment*t;
+
+ Float_t basearea1 = 0.5*sblen*TMath::Power(radius,2);
+ Float_t basearea2 = 0.5*sblen*TMath::Power((radius-segment),2);
+ Float_t basearea = basearea1 - basearea2;
+
+ Float_t area1 = 0.5*slen*TMath::Power(radius,2);
+ Float_t area2 = 0.5*slen*TMath::Power((radius-segment),2);
+ Float_t area = area1 - area2;
+
+ return area/basearea;
+}
+
+//____________________________________________________________________
+void
+AliFMDDensityCalculator::ScaleHistograms(Int_t nEvents)
+{
+ if (nEvents <= 0) return;
+
+ TIter next(&fRingHistos);
+ RingHistos* o = 0;
+ while ((o = static_cast<RingHistos*>(next()))) {
+ o->fDensity->Scale(1. / nEvents);
+ }
+}
+
+//____________________________________________________________________
+void
+AliFMDDensityCalculator::Output(TList* dir)
+{
+ TList* d = new TList;
+ d->SetName(GetName());
+ dir->Add(d);
+ TIter next(&fRingHistos);
+ RingHistos* o = 0;
+ while ((o = static_cast<RingHistos*>(next()))) {
+ o->Output(d);
+ }
+}
+
+//====================================================================
+AliFMDDensityCalculator::RingHistos::RingHistos()
+ : fDet(0),
+ fRing('\0'),
+ fEvsN(0),
+ fEvsM(0),
+ fDensity(0)
+{}
+
+//____________________________________________________________________
+AliFMDDensityCalculator::RingHistos::RingHistos(UShort_t d, Char_t r)
+ : fDet(d),
+ fRing(r),
+ fEvsN(0),
+ fEvsM(0),
+ fDensity(0)
+{
+ fEvsN = new TH2D(Form("FMD%d%c_Eloss_N_nocorr", d, r),
+ Form("Energy loss vs uncorrected inclusive "
+ "N_{ch} in FMD%d%c", d, r),
+ 100, -.5, 24.5, 100, -.5, 24.5);
+ fEvsM = new TH2D(Form("FMD%d%c_Eloss_N_corr", d, r),
+ Form("Energy loss vs corrected inclusive N_{ch} in FMD%d%c",
+ d, r), 100, -.5, 24.5, 100, -.5, 24.5);
+ fEvsN->SetXTitle("#Delta E/#Delta E_{mip}");
+ fEvsN->SetYTitle("Inclusive N_{ch} (uncorrected)");
+ fEvsN->Sumw2();
+ fEvsN->SetDirectory(0);
+ fEvsN->SetXTitle("#Delta E/#Delta E_{mip}");
+ fEvsN->SetYTitle("Inclusive N_{ch} (corrected)");
+ fEvsM->Sumw2();
+ fEvsM->SetDirectory(0);
+
+ fDensity = new TH2D(Form("FMD%d%c_Incl_Density", d, r),
+ Form("in FMD%d%c", d, r),
+ 200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40),
+ 0, 2*TMath::Pi());
+ fDensity->SetDirectory(0);
+ fDensity->SetXTitle("#eta");
+ fDensity->SetYTitle("#phi [radians]");
+ fDensity->SetZTitle("Inclusive N_{ch} density");
+}
+//____________________________________________________________________
+AliFMDDensityCalculator::RingHistos::RingHistos(const RingHistos& o)
+ : TObject(o),
+ fDet(o.fDet),
+ fRing(o.fRing),
+ fEvsN(o.fEvsN),
+ fEvsM(o.fEvsM),
+ fDensity(o.fDensity)
+{}
+
+//____________________________________________________________________
+AliFMDDensityCalculator::RingHistos&
+AliFMDDensityCalculator::RingHistos::operator=(const RingHistos& o)
+{
+ fDet = o.fDet;
+ fRing = o.fRing;
+
+ if (fEvsN) delete fEvsN;
+ if (fEvsM) delete fEvsM;
+ if (fDensity) delete fDensity;
+
+ fEvsN = static_cast<TH2D*>(o.fEvsN->Clone());
+ fEvsM = static_cast<TH2D*>(o.fEvsM->Clone());
+ fDensity = static_cast<TH2D*>(o.fDensity->Clone());
+
+ return *this;
+}
+//____________________________________________________________________
+AliFMDDensityCalculator::RingHistos::~RingHistos()
+{
+ if (fEvsN) delete fEvsN;
+ if (fEvsM) delete fEvsM;
+ if (fDensity) delete fDensity;
+}
+
+//____________________________________________________________________
+void
+AliFMDDensityCalculator::RingHistos::Output(TList* dir)
+{
+ TList* d = new TList;
+ d->SetName(Form("FMD%d%c", fDet, fRing));
+ d->Add(fEvsN);
+ d->Add(fEvsM);
+ d->Add(fDensity);
+ dir->Add(d);
+}
+
+//____________________________________________________________________
+//
+// EOF
+//
+
+
+
--- /dev/null
+#ifndef ALIROOT_PWG2_FORWARD_ANALYSIS2_ALIFMDDENSITYCALCULATOR_H
+#define ALIROOT_PWG2_FORWARD_ANALYSIS2_ALIFMDDENSITYCALCULATOR_H
+#include <TNamed.h>
+#include <TList.h>
+#include "AliForwardUtil.h"
+class AliESDFMD;
+class TH2D;
+
+/**
+ * This class calculates the inclusive charged particle density
+ * in each for the 5 FMD rings.
+ *
+ * @par Input:
+ * - AliESDFMD object possibly corrected for sharing
+ *
+ * @par Output:
+ * - 5 RingHistos objects - each with a number of vertex dependent
+ * 2D histograms of the inclusive charge particle density
+ *
+ * @par Corrections used:
+ * - AliFMDAnaCalibEnergyDistribution
+ * - AliFMDDoubleHitCorrection
+ * - AliFMDDeadCorrection
+ *
+ * @ingroup pwg2_forward_analysis
+ */
+class AliFMDDensityCalculator : public TNamed
+{
+public:
+ /**
+ * Constructor
+ */
+ AliFMDDensityCalculator();
+ /**
+ * Constructor
+ *
+ * @param name Name of object
+ */
+ AliFMDDensityCalculator(const char* name);
+ /**
+ * Copy constructor
+ *
+ * @param o Object to copy from
+ */
+ AliFMDDensityCalculator(const AliFMDDensityCalculator& o);
+ /**
+ * Destructor
+ */
+ virtual ~AliFMDDensityCalculator();
+ /**
+ * Assignement operator
+ *
+ * @param o Object to assign from
+ *
+ * @return Reference to this object
+ */
+ AliFMDDensityCalculator& operator=(const AliFMDDensityCalculator&);
+ /**
+ * Do the calculations
+ *
+ * @param fmd AliESDFMD object (possibly) corrected for sharing
+ * @param hists Histogram cache
+ * @param vtxBin Vertex bin
+ * @param lowFlux Low flux flag.
+ *
+ * @return true on successs
+ */
+ virtual Bool_t Calculate(const AliESDFMD& fmd,
+ AliForwardUtil::Histos& hists,
+ Int_t vtxBin, Bool_t lowFlux);
+ /**
+ * Scale the histograms to the total number of events
+ *
+ * @param nEvents Number of events
+ */
+ void ScaleHistograms(Int_t nEvents);
+ /**
+ * Output diagnostic histograms to directory
+ *
+ * @param dir List to write in
+ */
+ void Output(TList* dir);
+protected:
+ /**
+ * Get the number of particles corresponding to the signal mult
+ *
+ * @param mult Signal
+ * @param d Detector
+ * @param r Ring
+ * @param s Sector
+ * @param t Strip (not used)
+ * @param v Vertex bin
+ * @param eta Pseudo-rapidity
+ * @param lowFlux Low-flux flag
+ *
+ * @return The number of particles
+ */
+ virtual Float_t NParticles(Float_t mult,
+ UShort_t d, Char_t r, UShort_t s, UShort_t t,
+ Int_t v, Float_t eta, Bool_t lowFlux) const;
+ /**
+ * Get the inverse correction factor. This consist of
+ *
+ * - acceptance correction (corners of sensors)
+ * - double hit correction (for low-flux events)
+ * - dead strip correction
+ *
+ * @param d Detector
+ * @param r Ring
+ * @param s Sector
+ * @param t Strip (not used)
+ * @param v Vertex bin
+ * @param eta Pseudo-rapidity
+ * @param lowFlux Low-flux flag
+ *
+ * @return
+ */
+ virtual Float_t Correction(UShort_t d, Char_t r, UShort_t s, UShort_t t,
+ Int_t v, Float_t eta, Bool_t lowFlux) const;
+ /**
+ * Get the acceptance correction for strip @a t in an ring of type @a r
+ *
+ * @param r Ring type ('I' or 'O')
+ * @param t Strip number
+ *
+ * @return Inverse acceptance correction
+ */
+ virtual Float_t AcceptanceCorrection(Char_t r, UShort_t t) const;
+
+ /**
+ * Internal data structure to keep track of the histograms
+ */
+ struct RingHistos : public TObject
+ {
+ /**
+ * Default CTOR
+ */
+ RingHistos();
+ /**
+ * Constructor
+ *
+ * @param d detector
+ * @param r ring
+ */
+ RingHistos(UShort_t d, Char_t r);
+ /**
+ * Copy constructor
+ *
+ * @param o Object to copy from
+ */
+ RingHistos(const RingHistos& o);
+ /**
+ * Assignment operator
+ *
+ * @param o Object to assign from
+ *
+ * @return Reference to this
+ */
+ RingHistos& operator=(const RingHistos& o);
+ /**
+ * Destructor
+ */
+ ~RingHistos();
+ void Output(TList* dir);
+ UShort_t fDet; // Detector
+ Char_t fRing; // Ring
+ TH2D* fEvsN; // Correlation of Eloss vs uncorrected Nch
+ TH2D* fEvsM; // Correlation of Eloss vs corrected Nch
+ TH2D* fDensity; // Distribution inclusive Nch
+ };
+ /**
+ * Get the ring histogram container
+ *
+ * @param d Detector
+ * @param r Ring
+ *
+ * @return Ring histogram container
+ */
+ RingHistos* GetRingHistos(UShort_t d, Char_t r) const;
+
+ TList fRingHistos; // List of histogram containers
+ Double_t fMultCut; // Low cut on scaled energy loss
+
+ ClassDef(AliFMDDensityCalculator,1); // Calculate Nch density
+};
+
+#endif
+// Local Variables:
+// mode: C++
+// End:
+
--- /dev/null
+
+#include "AliFMDHistCollector.h"
+#include <AliESDFMD.h>
+#include <TAxis.h>
+#include <TList.h>
+#include <TMath.h>
+#include "AliFMDAnaParameters.h"
+#include "AliLog.h"
+#include <TH2D.h>
+#include <TH1I.h>
+#include <TProfile.h>
+#include <TProfile2D.h>
+#include <TArrayI.h>
+
+ClassImp(AliFMDHistCollector)
+#if 0
+; // For Emacs
+#endif
+
+
+//____________________________________________________________________
+AliFMDHistCollector&
+AliFMDHistCollector::operator=(const AliFMDHistCollector& o)
+{
+ SetName(o.GetName());
+ SetTitle(o.GetTitle());
+
+ fNCutBins = o.fNCutBins;
+ fCorrectionCut = o.fCorrectionCut;
+ fFirstBins = o.fFirstBins;
+ fLastBins = o.fLastBins;
+
+ return *this;
+}
+
+//____________________________________________________________________
+void
+AliFMDHistCollector::Init(const TAxis& vtxAxis)
+{
+ AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+
+ Int_t nVz = vtxAxis.GetNbins();
+ fFirstBins.Set(5*nVz);
+ fLastBins.Set(5*nVz);
+
+ // Find the eta bin ranges
+ for (Int_t iVz = 0; iVz < nVz; iVz++) {
+ // Find the first and last eta bin to use for each ring for
+ // each vertex bin. This is instead of using the methods
+ // provided by AliFMDAnaParameters
+ for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
+ UShort_t d = 0;
+ Char_t r = 0;
+ GetDetRing(iIdx, d, r);
+
+ // Get the background object
+ TH2F* bg = pars->GetBackgroundCorrection(d,r,iVz);
+ Int_t nEta = bg->GetNbinsX();
+ Int_t first = nEta+1;
+ Int_t last = 0;
+
+ // Loop over the eta bins
+ for (Int_t ie = 1; ie <= nEta; ie++) {
+ if (bg->GetBinContent(ie,1) < fCorrectionCut) continue;
+
+ // Loop over the phi bins to make sure that we
+ // do not have holes in the coverage
+ bool ok = true;
+ for (Int_t ip = 1; ip <= bg->GetNbinsY(); ip++) {
+ if (bg->GetBinContent(ie,ip) < 0.001) {
+ ok = false;
+ continue;
+ }
+ }
+ if (!ok) continue;
+
+ first = TMath::Min(ie, first);
+ last = TMath::Max(ie, last);
+ }
+
+ // Store the result for later use
+ fFirstBins[iVz*5+iIdx] = first;
+ fLastBins[iVz*5+iIdx] = last;
+ } // for j
+ }
+}
+
+//____________________________________________________________________
+Int_t
+AliFMDHistCollector::GetIdx(UShort_t d, Char_t r) const
+{
+ Int_t idx = -1;
+ switch (d) {
+ case 1: idx = 0; break;
+ case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
+ case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
+ }
+ return idx;
+}
+//____________________________________________________________________
+void
+AliFMDHistCollector::GetDetRing(Int_t idx, UShort_t& d, Char_t& r) const
+{
+ d = 0;
+ r = '\0';
+ switch (idx) {
+ case 0: d = 1; r = 'I'; break;
+ case 1: d = 2; r = 'I'; break;
+ case 2: d = 2; r = 'O'; break;
+ case 3: d = 3; r = 'I'; break;
+ case 4: d = 3; r = 'O'; break;
+ }
+}
+
+//____________________________________________________________________
+void
+AliFMDHistCollector::GetFirstAndLast(Int_t idx, Int_t vtxbin,
+ Int_t& first, Int_t& last) const
+{
+ first = 0;
+ last = 0;
+
+ if (idx < 0) return;
+ idx += vtxbin * 5;
+
+ if (idx < 0 || idx >= fFirstBins.GetSize()) return;
+
+ first = fFirstBins.At(idx)+fNCutBins;
+ last = fLastBins.At(idx)-fNCutBins;
+}
+
+//____________________________________________________________________
+Int_t
+AliFMDHistCollector::GetFirst(Int_t idx, Int_t v) const
+{
+ Int_t f, l;
+ GetFirstAndLast(idx,v,f,l);
+ return f;
+}
+
+
+//____________________________________________________________________
+Int_t
+AliFMDHistCollector::GetLast(Int_t idx, Int_t v) const
+{
+ Int_t f, l;
+ GetFirstAndLast(idx,v,f,l);
+ return l;
+}
+
+//____________________________________________________________________
+Int_t
+AliFMDHistCollector::GetOverlap(UShort_t d, Char_t r,
+ Int_t bin, Int_t vtxbin) const
+{
+ // Get the possibly overlapping histogram
+ Int_t other = -1;
+ if (d == 1) {
+ if (bin <= GetLast(2,'I',vtxbin)) other = GetIdx(2,'I');
+ }
+ else if (d == 2 && r == 'I') {
+ if (bin <= GetLast(2, 'O', vtxbin)) other = GetIdx(2, 'O');
+ else if (bin >= GetFirst(1, 'I', vtxbin)) other = GetIdx(1, 'I');
+ }
+ else if (d == 2 && r == 'O') {
+ if (bin >= GetFirst(2, 'I', vtxbin)) other = GetIdx(2,'I');
+ }
+ else if (d == 3 && r == 'O') {
+ if (bin <= GetLast(3, 'I', vtxbin)) other = GetIdx(3, 'I');
+ }
+ else if (d == 3 && r == 'I') {
+ if (bin >= GetFirst(3, 'O', vtxbin)) other = GetIdx(3, 'O');
+ }
+ // AliInfo(Form("FMD%d%c (%d) -> %d", d, r, GetIdx(d,r), other));
+ return other;
+}
+//____________________________________________________________________
+Int_t
+AliFMDHistCollector::GetOverlap(Int_t idx, Int_t bin, Int_t vtxbin) const
+{
+ UShort_t d = 0;
+ Char_t r = '\0';
+ GetDetRing(idx, d, r);
+ if (d == 0 || r == '\0') return 0;
+
+ return GetOverlap(d, r, bin, vtxbin);
+}
+
+
+//____________________________________________________________________
+Bool_t
+AliFMDHistCollector::Collect(AliForwardUtil::Histos& hists,
+ Int_t vtxbin,
+ TH2D& out)
+{
+ 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');
+ TH2D* h = hists.Get(d,r);
+ TH2D* t = static_cast<TH2D*>(h->Clone(Form("FMD%d%c_tmp",d,r)));
+
+ // Outer rings have better phi segmentation - rebin to same as inner.
+ if (q == 1) t->RebinY(2);
+
+ Int_t first = 0;
+ Int_t last = 0;
+ GetFirstAndLast(d, r, vtxbin, first, last);
+
+ // Now update profile output
+ for (Int_t iEta = first; iEta <= last; iEta++) {
+
+ // Get the possibly overlapping histogram
+ Int_t overlap = GetOverlap(d,r,iEta,vtxbin);
+
+ // Fill eta acceptance for this event into the phi underlow bin
+ Float_t ooc = out.GetBinContent(iEta,0);
+ Float_t noc = overlap >= 0? 0.5 : 1;
+ out.SetBinContent(iEta, 0, ooc + noc);
+
+ for (Int_t iPhi = 1; iPhi <= h->GetNbinsY(); iPhi++) {
+ Double_t c = t->GetBinContent(iEta,iPhi);
+ Double_t e = t->GetBinError(iEta,iPhi);
+
+ // If there's no signal, continue
+ // if (e <= 0) continue;
+ if (c <= 0 || e <= 0) continue;
+
+ // If there's no overlapping histogram (ring), then
+ // fill in data and continue to the next phi bin
+ if (overlap < 0) {
+ out.SetBinContent(iEta,iPhi,c);
+ out.SetBinError(iEta,iPhi,e);
+ continue;
+ }
+
+ // Get the current bin content and error
+ Double_t oc = out.GetBinContent(iEta,iPhi);
+ Double_t oe = out.GetBinError(iEta,iPhi);
+
+#define USE_STRAIGHT_MEAN
+// #define USE_STRAIGHT_MEAN_NONZERO
+// #define USE_WEIGHTED_MEAN
+// #define USE_MOST_CERTAIN
+#if defined(USE_STRAIGHT_MEAN)
+ // calculate the average of old value (half the original),
+ // and this value, as well as the summed squared errors
+ // of the existing content (sqrt((e_1/2)^2=sqrt(e_1^2/4)=e_1/2)
+ // and half the error of this.
+ //
+ // So, on the first overlapping histogram we get
+ //
+ // c = c_1 / 2
+ // e = sqrt((e_1 / 2)^2) = e_1/2
+ //
+ // On the second we get
+ //
+ // c' = c_2 / 2 + c = c_2 / 2 + c_1 / 2 = (c_1+c_2)/2
+ // e' = sqrt(e^2 + (e_2/2)^2)
+ // = sqrt(e_1^2/4 + e_2^2/4)
+ // = sqrt(1/4 * (e_1^2+e_2^2))
+ // = 1/2 * sqrt(e_1^2 + e_2^2)
+ out.SetBinContent(iEta,iPhi,oc + c/2);
+ out.SetBinError(iEta,iPhi,TMath::Sqrt(oe*oe+(e*e)/4));
+#elif defined(USE_STRAIGHT_MEAN_NONZERO)
+# define ZERO_OTHER
+ // If there's data in the overlapping histogram,
+ // calculate the average and add the errors in
+ // quadrature.
+ //
+ // If there's no data in the overlapping histogram,
+ // then just fill in the data
+ if (oe > 0) {
+ out.SetBinContent(iEta,iPhi,(oc + c)/2);
+ out.SetBinError(iEta,iPhi,TMath::Sqrt(oe*oe + e*e)/2);
+ }
+ else {
+ out.SetBinContent(iEta,iPhi,c);
+ out.SetBinError(iEta,iPhi,e);
+ }
+#elif defined(USE_WEIGHTED_MEAN)
+ // Calculate the weighted mean
+ Double_t w = 1/(e*e);
+ Double_t sc = w * c;
+ Double_t sw = w;
+ if (oe > 0) {
+ Double_t ow = 1/(oe*oe);
+ sc += ow * oc;
+ sw += ow;
+ }
+ Double_t nc = sc / sw;
+ Double_t ne = TMath::Sqrt(1 / sw);
+ out.SetBinContent(iEta,iPhi,nc,ne);
+#elif defined(USE_MOST_CERTAIN)
+# define ZERO_OTHER
+ if (e < oe) {
+ out.SetBinContent(iEta,iPhi,c);
+ out.SetBinError(iEta,iPhi,e);
+ }
+ else {
+ out.SetBinContent(iEta,iPhi,oc);
+ out.SetBinError(iEta,iPhi,oe);
+ }
+#else
+# error No method for defining content of overlapping bins defined
+#endif
+#if defined(ZERO_OTHER)
+ // Get the content of the overlapping histogram,
+ // and zero the content so that we won't use it
+ // again
+ UShort_t od; Char_t oR;
+ GetDetRing(overlap, od, oR);
+ TH2D* other = hists.Get(od,oR);
+ other->SetBinContent(iEta,iPhi,0);
+ other->SetBinError(iEta,iPhi,0);
+#endif
+ }
+ }
+ // Remove temporary histogram
+ delete t;
+ } // for r
+ } // for d
+ return kTRUE;
+}
+
+
+//____________________________________________________________________
+//
+// EOF
+//
+
+
+
--- /dev/null
+#ifndef ALIROOT_PWG2_FORWARD_ANALYSIS2_ALIFMDHISTCOLLECTOR_H
+#define ALIROOT_PWG2_FORWARD_ANALYSIS2_ALIFMDHISTCOLLECTOR_H
+#include <TNamed.h>
+#include <TList.h>
+#include <TArrayI.h>
+#include "AliForwardUtil.h"
+class AliESDFMD;
+class TH2D;
+
+/**
+ * This class collects the event histograms into single histograms,
+ * one for each ring in each vertex bin.
+ *
+ * @par Input:
+ * - AliESDFMD object possibly corrected for sharing
+ *
+ * @par Output:
+ * - 5 RingHistos objects - each with a number of vertex dependent
+ * 2D histograms of the inclusive charge particle density
+ *
+ * @par HistCollector used:
+ * - AliFMDAnaCalibBackgroundCorrection
+ *
+ * @ingroup pwg2_forward_analysis
+ */
+class AliFMDHistCollector : public TNamed
+{
+public:
+ /**
+ * Constructor
+ */
+ AliFMDHistCollector()
+ : fNCutBins(0), fCorrectionCut(0), fFirstBins(), fLastBins()
+ {}
+ /**
+ * Constructor
+ *
+ * @param title Name of object
+ */
+ AliFMDHistCollector(const char* title)
+ : TNamed("fmdHistCollector", title),
+ fNCutBins(1), fCorrectionCut(0.5), fFirstBins(1), fLastBins(1)
+ {}
+ /**
+ * Copy constructor
+ *
+ * @param o Object to copy from
+ */
+ AliFMDHistCollector(const AliFMDHistCollector& o)
+ : TNamed(o),
+ fNCutBins(o.fNCutBins), fCorrectionCut(o.fCorrectionCut),
+ fFirstBins(o.fFirstBins), fLastBins(o.fLastBins)
+ {}
+
+ /**
+ * Destructor
+ */
+ virtual ~AliFMDHistCollector() {}
+ /**
+ * Assignement operator
+ *
+ * @return Reference to this object
+ */
+ AliFMDHistCollector& operator=(const AliFMDHistCollector&);
+ /**
+ * Intialise
+ *
+ * @param vtxAxis Vertex axis
+ */
+ virtual void Init(const TAxis& vtxAxis);
+ /**
+ * Do the calculations
+ *
+ * @param hists Cache of histograms
+ * @param vtxBin Vertex bin
+ * @param out Output histogram
+ *
+ * @return true on successs
+ */
+ virtual Bool_t Collect(AliForwardUtil::Histos& hists, Int_t vtxBin,
+ TH2D& out);
+ /**
+ * Set the number of extra bins (beyond the secondary map border)
+ * to cut away.
+ *
+ * @param n Number of bins
+ */
+ void SetNCutBins(UInt_t n=2) { fNCutBins = n; }
+ /**
+ * Set the correction cut, that is, when bins in the secondary
+ * correction maps have a value lower than this cut, they are
+ * considered uncertain and not used
+ *
+ * @param cut Cut-off
+ */
+ void SetCorrectionCut(Float_t cut=0.5) { fCorrectionCut = cut; }
+protected:
+ /**
+ * Get the first and last eta bin to use for a given ring and vertex
+ *
+ * @param d Detector
+ * @param r Ring
+ * @param vtxBin Vertex bin
+ * @param first On return, the first eta bin to use
+ * @param last On return, the last eta bin to use
+ */
+ virtual void GetFirstAndLast(UShort_t d, Char_t r, Int_t vtxBin,
+ Int_t& first, Int_t& last) const;
+ /**
+ * Get the first and last eta bin to use for a given ring and vertex
+ *
+ * @param idx Ring index as given by GetIdx
+ * @param vtxBin Vertex bin
+ * @param first On return, the first eta bin to use
+ * @param last On return, the last eta bin to use
+ */
+ virtual void GetFirstAndLast(Int_t idx, Int_t vtxBin,
+ Int_t& first, Int_t& last) const;
+ /**
+ * Get the first eta bin to use for a given ring and vertex
+ *
+ * @param d Detector
+ * @param r Ring
+ * @param v vertex bin
+ *
+ * @return First eta bin to use, or -1 in case of problems
+ */
+ Int_t GetFirst(UShort_t d, Char_t r, Int_t v) const;
+ /**
+ * Get the first eta bin to use for a given ring and vertex
+ *
+ * @param idx Ring index as given by GetIdx
+ * @param v vertex bin
+ *
+ * @return First eta bin to use, or -1 in case of problems
+ */
+ Int_t GetFirst(Int_t idx, Int_t v) const;
+ /**
+ * Get the last eta bin to use for a given ring and vertex
+ *
+ * @param d Detector
+ * @param r Ring
+ * @param v vertex bin
+ *
+ * @return Last eta bin to use, or -1 in case of problems
+ */
+ Int_t GetLast(UShort_t d, Char_t r, Int_t v) const;
+ /**
+ * Get the last eta bin to use for a given ring and vertex
+ *
+ * @param idx Ring index as given by GetIdx
+ * @param v vertex bin
+ *
+ * @return Last eta bin to use, or -1 in case of problems
+ */
+ Int_t GetLast(Int_t idx, Int_t v) const;
+ /**
+ * Get the detector and ring from the ring index
+ *
+ * @param idx Ring index
+ * @param d On return, the detector or 0 in case of errors
+ * @param r On return, the ring id or '0' in case of errors
+ */
+ void GetDetRing(Int_t idx, UShort_t& d, Char_t& r) const;
+ /**
+ * Get the ring index from detector number and ring identifier
+ *
+ * @param d Detector
+ * @param r Ring identifier
+ *
+ * @return ring index or -1 in case of problems
+ */
+ Int_t GetIdx(UShort_t d, Char_t r) const;
+ /**
+ * Get the possibly overlapping histogram of eta bin @a e in
+ * detector and ring
+ *
+ * @param d Detector
+ * @param r Ring
+ * @param e Eta bin
+ * @param v Vertex bin
+ *
+ * @return Overlapping histogram index or -1
+ */
+ Int_t GetOverlap(UShort_t d, Char_t r, Int_t e, Int_t v) const;
+ /**
+ * Get the possibly overlapping histogram of eta bin @a e in
+ * detector and ring
+ *
+ * @param i Ring index
+ * @param e Eta bin
+ * @param v Vertex bin
+ *
+ * @return Overlapping histogram index or -1
+ */
+ Int_t GetOverlap(Int_t i, Int_t e, Int_t v) const;
+ /**
+ * Check if there's an overlapping histogram with this eta bin of
+ * the detector and ring
+ *
+ * @param d Detector
+ * @param r Ring
+ * @param e eta bin
+ * @param v Vertex bin
+ *
+ * @return True if there's an overlapping histogram
+ */
+ Bool_t HasOverlap(UShort_t d, Char_t r, Int_t e, Int_t v) const;
+ /**
+ * Check if there's an overlapping histogram with this eta bin of
+ * ring
+ *
+ * @param i Ring index
+ * @param e eta bin
+ * @param v Vertex bin
+ *
+ * @return True if there's an overlapping histogram
+ */
+ Bool_t HasOverlap(Int_t i, Int_t e, Int_t v) const;
+
+
+ Int_t fNCutBins; // Number of additional bins to cut away
+ Float_t fCorrectionCut; // Cut-off on secondary corrections
+ TArrayI fFirstBins; // Array of first eta bins
+ TArrayI fLastBins; // Array of last eta bins
+
+ ClassDef(AliFMDHistCollector,1); // Calculate Nch density
+};
+
+//____________________________________________________________________
+inline void
+AliFMDHistCollector::GetFirstAndLast(UShort_t d, Char_t r, Int_t vtxbin,
+ Int_t& first, Int_t& last) const
+{
+ GetFirstAndLast(GetIdx(d,r), vtxbin, first, last);
+}
+//____________________________________________________________________
+inline Int_t
+AliFMDHistCollector::GetFirst(UShort_t d, Char_t r, Int_t v) const
+{
+ return GetFirst(GetIdx(d,r), v);
+}
+//____________________________________________________________________
+inline Int_t
+AliFMDHistCollector::GetLast(UShort_t d, Char_t r, Int_t v) const
+{
+ return GetLast(GetIdx(d, r), v);
+}
+//____________________________________________________________________
+inline Bool_t
+AliFMDHistCollector::HasOverlap(UShort_t d, Char_t r, Int_t e, Int_t v) const
+{
+ return GetOverlap(d,r,e,v) >= 0;
+}
+//____________________________________________________________________
+inline Bool_t
+AliFMDHistCollector::HasOverlap(Int_t i, Int_t e, Int_t v) const
+{
+ return GetOverlap(i,e,v) >= 0;
+}
+
+#endif
+// Local Variables:
+// mode: C++
+// End:
+
--- /dev/null
+//
+// Class to do the sharing correction of FMD ESD data
+//
+#include "AliFMDSharingFilter.h"
+#include <AliESDFMD.h>
+#include <TAxis.h>
+#include <TList.h>
+#include <TH1.h>
+#include <TMath.h>
+#include "AliFMDAnaParameters.h"
+#include <AliLog.h>
+
+ClassImp(AliFMDSharingFilter)
+#if 0
+; // This is for Emacs
+#endif
+
+
+//____________________________________________________________________
+AliFMDSharingFilter::AliFMDSharingFilter()
+ : TNamed(),
+ fRingHistos(),
+ fLowCut(0.3),
+ fCorrectAngles(kFALSE),
+ fEtaCorr(0)
+{}
+
+//____________________________________________________________________
+AliFMDSharingFilter::AliFMDSharingFilter(const char* title)
+ : TNamed("fmdSharingFilter", title),
+ fRingHistos(),
+ fLowCut(0.3),
+ fCorrectAngles(kFALSE),
+ fEtaCorr(0)
+{
+ fRingHistos.Add(new RingHistos(1, 'I'));
+ fRingHistos.Add(new RingHistos(2, 'I'));
+ fRingHistos.Add(new RingHistos(2, 'O'));
+ fRingHistos.Add(new RingHistos(3, 'I'));
+ fRingHistos.Add(new RingHistos(3, 'O'));
+ fEtaCorr = new TH2F("etaCorr", "Correction of #eta",
+ 200, -4, 6, 200, -4, 6);
+ fEtaCorr->SetXTitle("#eta from ESD");
+ fEtaCorr->SetYTitle("#eta from AnaParameters");
+}
+
+//____________________________________________________________________
+AliFMDSharingFilter::AliFMDSharingFilter(const AliFMDSharingFilter& o)
+ : TNamed(o),
+ fRingHistos(),
+ fLowCut(o.fLowCut),
+ fCorrectAngles(o.fCorrectAngles),
+ fEtaCorr(o.fEtaCorr)
+{
+ TIter next(&o.fRingHistos);
+ TObject* obj = 0;
+ while ((obj = next())) fRingHistos.Add(obj);
+}
+
+//____________________________________________________________________
+AliFMDSharingFilter::~AliFMDSharingFilter()
+{
+ fRingHistos.Delete();
+}
+
+//____________________________________________________________________
+AliFMDSharingFilter&
+AliFMDSharingFilter::operator=(const AliFMDSharingFilter& o)
+{
+ SetName(o.GetName());
+ SetTitle(o.GetTitle());
+
+ fLowCut = o.fLowCut;
+ fCorrectAngles = o.fCorrectAngles;
+
+ fRingHistos.Delete();
+ TIter next(&o.fRingHistos);
+ TObject* obj = 0;
+ while ((obj = next())) fRingHistos.Add(obj);
+
+ return *this;
+}
+
+//____________________________________________________________________
+AliFMDSharingFilter::RingHistos*
+AliFMDSharingFilter::GetRingHistos(UShort_t d, Char_t r) const
+{
+ Int_t idx = -1;
+ switch (d) {
+ case 1: idx = 0; break;
+ case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
+ case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
+ }
+ if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0;
+
+ return static_cast<RingHistos*>(fRingHistos.At(idx));
+}
+
+//____________________________________________________________________
+Bool_t
+AliFMDSharingFilter::Filter(const AliESDFMD& input,
+ Bool_t lowFlux,
+ AliESDFMD& output,
+ Double_t vz)
+{
+ output.Clear();
+ TIter next(&fRingHistos);
+ RingHistos* o = 0;
+ while ((o = static_cast<RingHistos*>(next())))
+ o->Clear();
+
+ AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+
+ for(UShort_t d = 1; d <= 3; d++) {
+ Int_t nRings = (d == 1 ? 1 : 2);
+ for (UShort_t q = 0; q < nRings; q++) {
+ Char_t r = (q == 0 ? 'I' : 'O');
+ UShort_t nsec = (q == 0 ? 20 : 40);
+ UShort_t nstr = (q == 0 ? 512 : 256);
+
+ RingHistos* histos = GetRingHistos(d, r);
+
+ for(UShort_t s = 0; s < nsec; s++) {
+ Bool_t usedThis = kFALSE;
+ Bool_t usedPrev = kFALSE;
+
+ for(UShort_t t = 0; t < nstr; t++) {
+ output.SetMultiplicity(d,r,s,t,0.);
+ Float_t mult = SignalInStrip(input,d,r,s,t);
+
+ // Keep dead-channel information.
+ if(mult == AliESDFMD::kInvalidMult)
+ output.SetMultiplicity(d,r,s,t,AliESDFMD::kInvalidMult);
+
+ // If no signal or dead strip, go on.
+ if (mult == AliESDFMD::kInvalidMult || mult == 0) continue;
+
+ // Get the pseudo-rapidity
+ Double_t eta1 = input.Eta(d,r,s,t);
+ Double_t eta2 = pars->GetEtaFromStrip(d,r,s,t,vz);
+ fEtaCorr->Fill(eta1, eta2);
+ Double_t eta = eta2;
+
+ // Fill the diagnostics histogram
+ histos->fBefore->Fill(mult);
+
+ // Get next and previous signal - if any
+ Double_t prevE = 0;
+ Double_t nextE = 0;
+ if (t != 0) {
+ prevE = SignalInStrip(input,d,r,s,t-1);
+ if (prevE == AliESDFMD::kInvalidMult) prevE = 0;
+ }
+ if (t != nstr - 1) {
+ nextE = SignalInStrip(input,d,r,s,t+1);
+ if (nextE == AliESDFMD::kInvalidMult) nextE = 0;
+ }
+
+ Double_t mergedEnergy = MultiplicityOfStrip(mult,eta,prevE,nextE,
+ lowFlux,d,r,s,t,
+ usedPrev,usedThis);
+ if (mergedEnergy > fLowCut) histos->Incr();
+ histos->fAfter->Fill(mergedEnergy);
+
+ output.SetMultiplicity(d,r,s,t,mergedEnergy);
+ output.SetEta(d,r,s,t,eta);
+ } // for strip
+ } // for sector
+ } // for ring
+ } // for detector
+
+ next.Reset();
+ while ((o = static_cast<RingHistos*>(next())))
+ o->Finish();
+
+ return kTRUE;
+}
+
+//_____________________________________________________________________
+Double_t
+AliFMDSharingFilter::SignalInStrip(const AliESDFMD& input,
+ UShort_t d,
+ Char_t r,
+ UShort_t s,
+ UShort_t t) const
+{
+ Double_t mult = input.Multiplicity(d,r,s,t);
+ if (mult == AliESDFMD::kInvalidMult || mult == 0) return mult;
+
+ if (fCorrectAngles && !input.IsAngleCorrected())
+ mult = AngleCorrect(mult, input.Eta(d,r,s,t));
+ else if (!fCorrectAngles && input.IsAngleCorrected())
+ mult = DeAngleCorrect(mult, input.Eta(d,r,s,t));
+ return mult;
+}
+
+//_____________________________________________________________________
+Double_t
+AliFMDSharingFilter::GetHighCut(UShort_t d, Char_t r, Double_t eta) const
+{
+ AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+
+ // Get the high cut. The high cut is defined as the
+ // most-probably-value peak found from the energy distributions, minus
+ // 2 times the width of the corresponding Landau.
+ Double_t mpv = pars->GetMPV(d,r,eta);
+ Double_t w = pars->GetSigma(d,r,eta);
+ if (mpv > 100)
+ AliWarning(Form("FMD%d%c, eta=%f, MPV=%f, W=%f", d, r, eta, mpv, w));
+ Double_t highCut = (mpv - 2 * w);
+ return highCut;
+}
+
+//_____________________________________________________________________
+Double_t
+AliFMDSharingFilter::MultiplicityOfStrip(Double_t mult,
+ Double_t eta,
+ Double_t prevE,
+ Double_t nextE,
+ Bool_t lowFlux,
+ UShort_t d,
+ Char_t r,
+ UShort_t /*s*/,
+ UShort_t /*t*/,
+ Bool_t& usedPrev,
+ Bool_t& usedThis)
+{
+ // IF the multiplicity is very large, or below the cut, reject it, and
+ // flag it as candidate for sharing
+ if(mult > 12 || mult < fLowCut) {
+ usedThis = kFALSE;
+ usedPrev = kFALSE;
+ return 0;
+ }
+
+ // If this was shared with the previous signal, return 0 and mark it as
+ // not a candiate for sharing
+ if(usedThis) {
+ usedThis = kFALSE;
+ usedPrev = kTRUE;
+ return 0.;
+ }
+
+ //analyse and perform sharing on one strip
+
+ // Calculate the total energy loss
+ Double_t highCut = GetHighCut(d, r, eta);
+
+ // If this signal is smaller than the next, and the next signal is smaller
+ // than then the high cut, and it's a low-flux event, then mark this and
+ // the next signal as candidates for sharing
+ //
+ // This is the test if the signal is the smaller of two possibly
+ // shared signals.
+ if (mult < nextE &&
+ nextE > highCut &&
+ lowFlux) {
+ usedThis = kFALSE;
+ usedPrev = kFALSE;
+ return 0;
+ }
+
+ Double_t totalE = mult;
+
+ // If the previous signal was larger than the low cut, and
+ // the previous was smaller than high cut, and the previous signal
+ // isn't marked as used, then add it's energy to this signal
+ if (prevE > fLowCut &&
+ prevE < highCut &&
+ !usedPrev)
+ totalE += prevE;
+
+ // If the next signal is larger than the low cut, and
+ // the next signal is smaller than the low cut, then add the next signal
+ // to this, and marked the next signal as used
+ if(nextE > fLowCut && nextE < highCut ) {
+ totalE += nextE;
+ usedThis = kTRUE;
+ }
+
+ // If we're processing on non-angle corrected data, we should do the
+ // angle correction here
+ if (!fCorrectAngles)
+ totalE = AngleCorrect(totalE, eta);
+
+ // Fill post processing histogram
+ // if(totalE > fLowCut)
+ // GetRingHistos(d,r)->fAfter->Fill(totalE);
+
+ // Return value
+ Double_t mergedEnergy = 0;
+
+ if(totalE > 0) {
+ // If we have a signal, then this is used (sharedPrev=true) and
+ // the signal is set to the result
+ mergedEnergy = totalE;
+ usedPrev = kTRUE;
+ }
+ else {
+ // If we do not have a signal (too low), then this is not shared,
+ // and the next strip is not shared either
+ usedThis = kFALSE;
+ usedPrev = kFALSE;
+ }
+
+ return mergedEnergy;
+}
+//____________________________________________________________________
+Double_t
+AliFMDSharingFilter::AngleCorrect(Double_t mult, Double_t eta) const
+{
+ Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
+ if (eta < 0) theta -= TMath::Pi();
+ return mult * TMath::Cos(theta);
+}
+//____________________________________________________________________
+Double_t
+AliFMDSharingFilter::DeAngleCorrect(Double_t mult, Double_t eta) const
+{
+ Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
+ if (eta < 0) theta -= TMath::Pi();
+ return mult / TMath::Cos(theta);
+}
+
+//____________________________________________________________________
+void
+AliFMDSharingFilter::ScaleHistograms(Int_t nEvents)
+{
+ if (nEvents <= 0) return;
+
+ TIter next(&fRingHistos);
+ RingHistos* o = 0;
+ while ((o = static_cast<RingHistos*>(next()))) {
+ o->fBefore->Scale(1. / nEvents);
+ o->fAfter->Scale(1. / nEvents);
+ }
+}
+
+//____________________________________________________________________
+void
+AliFMDSharingFilter::Output(TList* dir)
+{
+ TList* d = new TList;
+ d->SetName(GetName());
+ dir->Add(d);
+ d->Add(fEtaCorr);
+ TIter next(&fRingHistos);
+ RingHistos* o = 0;
+ while ((o = static_cast<RingHistos*>(next()))) {
+ o->Output(d);
+ }
+}
+
+//====================================================================
+AliFMDSharingFilter::RingHistos::RingHistos()
+ : fDet(0),
+ fRing('\0'),
+ fBefore(0),
+ fAfter(0),
+ fHits(0),
+ fNHits(0)
+{}
+
+//____________________________________________________________________
+AliFMDSharingFilter::RingHistos::RingHistos(UShort_t d, Char_t r)
+ : fDet(d),
+ fRing(r),
+ fBefore(0),
+ fAfter(0),
+ fHits(0),
+ fNHits(0)
+{
+ fBefore = new TH1D(Form("FMD%d%c_ESD_Eloss", d, r),
+ Form("Energy loss in FMD%d%c (reconstruction)", d, r),
+ 1000, 0, 25);
+ fAfter = new TH1D(Form("FMD%d%c_Ana_Eloss", d, r),
+ Form("Energy loss in FMD%d%c (sharing corrected)", d, r),
+ 1000, 0, 25);
+ fBefore->SetXTitle("#Delta E/#Delta E_{mip}");
+ fBefore->SetYTitle("P(#Delta E/#Delta E_{mip})");
+ fBefore->SetFillColor(kRed+1);
+ fBefore->SetFillStyle(3001);
+ // fBefore->Sumw2();
+ fBefore->SetDirectory(0);
+ fAfter->SetXTitle("#Delta E/#Delta E_{mip}");
+ fAfter->SetYTitle("P(#Delta E/#Delta E_{mip})");
+ fAfter->SetFillColor(kBlue+1);
+ fAfter->SetFillStyle(3001);
+ // fAfter->Sumw2();
+ fAfter->SetDirectory(0);
+
+ fHits = new TH1D(Form("FMD%d%c_Hits", d, r),
+ Form("Number of hits in FMD%d%c", d, r),
+ 200, 0, 200000);
+ fHits->SetDirectory(0);
+ fHits->SetXTitle("# of hits");
+ fHits->SetFillColor(kGreen+1);
+}
+//____________________________________________________________________
+AliFMDSharingFilter::RingHistos::RingHistos(const RingHistos& o)
+ : TObject(o),
+ fDet(o.fDet),
+ fRing(o.fRing),
+ fBefore(o.fBefore),
+ fAfter(o.fAfter),
+ fHits(o.fHits),
+ fNHits(o.fNHits)
+{}
+
+//____________________________________________________________________
+AliFMDSharingFilter::RingHistos&
+AliFMDSharingFilter::RingHistos::operator=(const RingHistos& o)
+{
+ fDet = o.fDet;
+ fRing = o.fRing;
+
+ if (fBefore) delete fBefore;
+ if (fAfter) delete fAfter;
+ if (fHits) delete fHits;
+
+ fBefore = static_cast<TH1D*>(o.fBefore->Clone());
+ fAfter = static_cast<TH1D*>(o.fAfter->Clone());
+ fHits = static_cast<TH1D*>(o.fHits->Clone());
+
+ return *this;
+}
+//____________________________________________________________________
+AliFMDSharingFilter::RingHistos::~RingHistos()
+{
+ if (fBefore) delete fBefore;
+ if (fAfter) delete fAfter;
+ if (fHits) delete fHits;
+}
+//____________________________________________________________________
+void
+AliFMDSharingFilter::RingHistos::Finish()
+{
+ fHits->Fill(fNHits);
+}
+
+//____________________________________________________________________
+void
+AliFMDSharingFilter::RingHistos::Output(TList* dir)
+{
+ TList* d = new TList;
+ d->SetName(Form("FMD%d%c", fDet, fRing));
+ d->Add(fBefore);
+ d->Add(fAfter);
+ d->Add(fHits);
+ dir->Add(d);
+}
+
+//____________________________________________________________________
+//
+// EOF
+//
--- /dev/null
+#ifndef ALIROOT_PWG2_FORWARD_ALIFMDSHARINGFILTER_H
+#define ALIROOT_PWG2_FORWARD_ALIFMDSHARINGFILTER_H
+#include <TNamed.h>
+#include <TH2.h>
+#include <TList.h>
+class AliESDFMD;
+class TAxis;
+class TList;
+class TH2;
+
+/**
+ * Class to do the sharing correction. That is, a filter that merges
+ * adjacent strip signals presumably originating from a single particle
+ * that impinges on the detector in such a way that it deposite energy
+ * into two or more strips.
+ *
+ * @par Input:
+ * - AliESDFMD object - from reconstruction
+ *
+ * @par Output:
+ * - AliESDFMD object - copy of input, but with signals merged
+ *
+ * @par Corrections used:
+ * - AliFMDAnaCalibEnergyDistribution objects
+ *
+ * @par Histograms:
+ * - For each ring (FMD1i, FMD2i, FMD2o, FMD3i, FMD3o) the distribution of
+ * signals before and after the filter.
+ * - For each ring (see above), an array of distributions of number of
+ * hit strips for each vertex bin (if enabled - see Init method)
+ *
+ *
+ * @ingroup pwg2_forward_analysis
+ */
+class AliFMDSharingFilter : public TNamed
+{
+public:
+ /**
+ * Destructor
+ */
+ virtual ~AliFMDSharingFilter();
+ /**
+ * Default Constructor - do not use
+ */
+ AliFMDSharingFilter();
+ /**
+ * Constructor
+ *
+ * @param title Title of object - not significant
+ */
+ AliFMDSharingFilter(const char* title);
+ /**
+ * Copy constructor
+ *
+ * @param o Object to copy from
+ */
+ AliFMDSharingFilter(const AliFMDSharingFilter& o);
+ /**
+ * Assignment operator
+ *
+ * @param o Object to assign from
+ *
+ * @return Reference to this
+ */
+ AliFMDSharingFilter& operator=(const AliFMDSharingFilter& o);
+
+ /**
+ * Initialise the filter
+ *
+ */
+ void Init() {}
+ /**
+ * Set the low cut used for sharing
+ *
+ * @param lowCut Low cut
+ */
+ void SetLowCut(Double_t lowCut=0.3) { fLowCut = lowCut; }
+
+ /**
+ * Enable use of angle corrected signals in the algorithm
+ *
+ * @param use If true, use angle corrected signals,
+ * otherwise use de-corrected signals. In the final output, the
+ * signals are always angle corrected.
+ */
+ void UseAngleCorrectedSignals(Bool_t use) { fCorrectAngles = use; }
+ /**
+ * Filter the input AliESDFMD object
+ *
+ * @param input Input
+ * @param lowFlux If this is a low-flux event
+ * @param output Output AliESDFMD object
+ * @param vz Current vertex position
+ *
+ * @return True on success, false otherwise
+ */
+ Bool_t Filter(const AliESDFMD& input,
+ Bool_t lowFlux,
+ AliESDFMD& output,
+ Double_t vz);
+ /**
+ * Scale the histograms to the total number of events
+ *
+ * @param nEvents Number of events
+ */
+ void ScaleHistograms(Int_t nEvents);
+
+ void Output(TList* dir);
+protected:
+ /**
+ * Internal data structure to keep track of the histograms
+ */
+ struct RingHistos : public TObject
+ {
+ /**
+ * Default CTOR
+ */
+ RingHistos();
+ /**
+ * Constructor
+ *
+ * @param d detector
+ * @param r ring
+ */
+ RingHistos(UShort_t d, Char_t r);
+ /**
+ * Copy constructor
+ *
+ * @param o Object to copy from
+ */
+ RingHistos(const RingHistos& o);
+ /**
+ * Assignment operator
+ *
+ * @param o Object to assign from
+ *
+ * @return Reference to this
+ */
+ RingHistos& operator=(const RingHistos& o);
+ /**
+ * Destructor
+ */
+ ~RingHistos();
+ /**
+ * Initialise this object
+ */
+ void Init() {}
+ void Clear(const Option_t* ="") { fNHits = 0; }
+ void Incr() { fNHits++; }
+ void Finish();
+ void Output(TList* dir);
+ UShort_t fDet; // Detector
+ Char_t fRing; // Ring
+ TH1D* fBefore; // Distribution of signals before filter
+ TH1D* fAfter; // Distribution of signals after filter
+ TH1D* fHits; // Distribution of hit strips.
+ Int_t fNHits; // Number of hit strips per event
+ };
+ /**
+ * Get the ring histogram container
+ *
+ * @param d Detector
+ * @param r Ring
+ *
+ * @return Ring histogram container
+ */
+ RingHistos* GetRingHistos(UShort_t d, Char_t r) const;
+ /**
+ * Get the signal in a strip
+ *
+ * @param fmd ESD object
+ * @param d Detector
+ * @param r Ring
+ * @param s Sector
+ * @param t Strip
+ *
+ * @return The energy signal
+ */
+ Double_t SignalInStrip(const AliESDFMD& fmd,
+ UShort_t d,
+ Char_t r,
+ UShort_t s,
+ UShort_t t) const;
+ /**
+ * The actual algorithm
+ *
+ * @param mult The unfiltered signal in the strip
+ * @param eta Psuedo rapidity
+ * @param prevE Previous strip signal (or 0)
+ * @param nextE Next strip signal (or 0)
+ * @param lowFlux Whether this is a low flux event
+ * @param d Detector
+ * @param r Ring
+ * @param s Sector
+ * @param t Strip
+ * @param usedPrev Whether the previous strip was used in sharing or not
+ * @param usedThis Wether this strip was used in sharing or not.
+ *
+ * @return The filtered signal in the strip
+ */
+ Double_t MultiplicityOfStrip(Double_t mult,
+ Double_t eta,
+ Double_t prevE,
+ Double_t nextE,
+ Bool_t lowFlux,
+ UShort_t d,
+ Char_t r,
+ UShort_t s,
+ UShort_t t,
+ Bool_t& usedPrev,
+ Bool_t& usedThis);
+ /**
+ * Angle correct the signal
+ *
+ * @param mult Angle Un-corrected Signal
+ * @param eta Pseudo-rapidity
+ *
+ * @return Angle corrected signal
+ */
+ Double_t AngleCorrect(Double_t mult, Double_t eta) const;
+ /**
+ * Angle de-correct the signal
+ *
+ * @param mult Angle corrected Signal
+ * @param eta Pseudo-rapidity
+ *
+ * @return Angle un-corrected signal
+ */
+ Double_t DeAngleCorrect(Double_t mult, Double_t eta) const;
+ /**
+ * Get the high cut. The high cut is defined as the
+ * most-probably-value peak found from the energy distributions, minus
+ * 2 times the width of the corresponding Landau.
+ */
+ virtual Double_t GetHighCut(UShort_t d, Char_t r, Double_t eta) const;
+
+ TList fRingHistos; // List of histogram containers
+ Double_t fLowCut; // Low cut on sharing
+ Bool_t fCorrectAngles; // Whether to work on angle corrected signals
+ TH2* fEtaCorr;
+
+ ClassDef(AliFMDSharingFilter,1); //
+};
+
+#endif
+// Local Variables:
+// mode: C++
+// End:
--- /dev/null
+#include "AliForwardMultiplicity.h"
+#include "AliTriggerAnalysis.h"
+#include "AliPhysicsSelection.h"
+#include "AliLog.h"
+#include "AliFMDAnaParameters.h"
+#include "AliESDEvent.h"
+#include "AliAODHandler.h"
+#include "AliMultiplicity.h"
+#include "AliInputEventHandler.h"
+#include <TH1.h>
+#include <TDirectory.h>
+#include <TTree.h>
+
+//====================================================================
+AliForwardMultiplicity::AliForwardMultiplicity()
+ : AliAnalysisTaskSE(),
+ fHEventsTr(0),
+ fHEventsTrVtx(0),
+ fHTriggers(0),
+ fHData(0),
+ fFirstEvent(true),
+ fLowFluxCut(1000),
+ fESDFMD(),
+ fHistos(),
+ fAODFMD(),
+ fSharingFilter(),
+ fDensityCalculator(),
+ fCorrections(),
+ fHistCollector(),
+ fList(0),
+ fTree(0)
+{
+}
+
+//____________________________________________________________________
+AliForwardMultiplicity::AliForwardMultiplicity(const char* name)
+ : AliAnalysisTaskSE(name),
+ fHEventsTr(0),
+ fHEventsTrVtx(0),
+ fHTriggers(0),
+ fHData(0),
+ fFirstEvent(true),
+ fLowFluxCut(1000),
+ fESDFMD(),
+ fHistos(),
+ fAODFMD(kTRUE),
+ fSharingFilter("sharing"),
+ fDensityCalculator("density"),
+ fCorrections("corrections"),
+ fHistCollector("collector"),
+ fList(0),
+ fTree(0)
+{
+ DefineOutput(1, TList::Class());
+ // DefineOutput(2, TTree::Class());
+}
+
+//____________________________________________________________________
+AliForwardMultiplicity::AliForwardMultiplicity(const AliForwardMultiplicity& o)
+ : AliAnalysisTaskSE(o),
+ fHEventsTr(o.fHEventsTr),
+ fHEventsTrVtx(o.fHEventsTrVtx),
+ fHTriggers(o.fHTriggers),
+ fHData(o.fHData),
+ fFirstEvent(true),
+ fLowFluxCut(1000),
+ fESDFMD(o.fESDFMD),
+ fHistos(o.fHistos),
+ fAODFMD(o.fAODFMD),
+ fSharingFilter(o.fSharingFilter),
+ fDensityCalculator(o.fDensityCalculator),
+ fCorrections(o.fCorrections),
+ fHistCollector(o.fHistCollector),
+ fList(o.fList),
+ fTree(o.fTree)
+{
+}
+
+//____________________________________________________________________
+AliForwardMultiplicity&
+AliForwardMultiplicity::operator=(const AliForwardMultiplicity& o)
+{
+ fHEventsTr = o.fHEventsTr;
+ fHEventsTrVtx = o.fHEventsTrVtx;
+ fHTriggers = o.fHTriggers;
+ fHData = o.fHData;
+ fFirstEvent = o.fFirstEvent;
+ fSharingFilter = o.fSharingFilter;
+ fDensityCalculator = o.fDensityCalculator;
+ fCorrections = o.fCorrections;
+ fHistCollector = o.fHistCollector;
+ fHistos = o.fHistos;
+ fAODFMD = o.fAODFMD;
+ fList = o.fList;
+ fTree = o.fTree;
+
+ return *this;
+}
+
+//____________________________________________________________________
+void
+AliForwardMultiplicity::Init()
+{
+ fFirstEvent = true;
+}
+
+//____________________________________________________________________
+void
+AliForwardMultiplicity::InitializeSubs()
+{
+ AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+ pars->Init(kTRUE);
+
+ fHEventsTr = new TH1I("nEvents", "Number of events w/trigger",
+ pars->GetNvtxBins(),
+ -pars->GetVtxCutZ(),
+ pars->GetVtxCutZ());
+ fHEventsTr->SetXTitle("v_{z} [cm]");
+ fHEventsTr->SetYTitle("# of events");
+ fHEventsTr->SetFillColor(kRed+1);
+ fHEventsTr->SetFillStyle(3001);
+ fHEventsTr->SetDirectory(0);
+ // fHEventsTr->Sumw2();
+
+ fHEventsTrVtx = new TH1I("nEventsTrVtx",
+ "Number of events w/trigger and vertex",
+ pars->GetNvtxBins(),
+ -pars->GetVtxCutZ(),
+ pars->GetVtxCutZ());
+ fHEventsTrVtx->SetXTitle("v_{z} [cm]");
+ fHEventsTrVtx->SetYTitle("# of events");
+ fHEventsTrVtx->SetFillColor(kBlue+1);
+ fHEventsTrVtx->SetFillStyle(3001);
+ fHEventsTrVtx->SetDirectory(0);
+ // fHEventsTrVtx->Sumw2();
+
+
+ fHTriggers = new TH1I("triggers", "Triggers", 10, 0, 10);
+ fHTriggers->SetFillColor(kRed+1);
+ fHTriggers->SetFillStyle(3001);
+ fHTriggers->SetStats(0);
+ fHTriggers->SetDirectory(0);
+ fHTriggers->GetXaxis()->SetBinLabel(1,"INEL");
+ fHTriggers->GetXaxis()->SetBinLabel(2,"INEL>0");
+ fHTriggers->GetXaxis()->SetBinLabel(3,"NSD");
+ fHTriggers->GetXaxis()->SetBinLabel(4,"Empty");
+ fHTriggers->GetXaxis()->SetBinLabel(5,"A");
+ fHTriggers->GetXaxis()->SetBinLabel(6,"B");
+ fHTriggers->GetXaxis()->SetBinLabel(7,"C");
+ fHTriggers->GetXaxis()->SetBinLabel(8,"E");
+
+ TAxis e(pars->GetNetaBins(), pars->GetEtaMin(), pars->GetEtaMax());
+ fHistos.Init(e);
+ fAODFMD.Init(e);
+
+ fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
+ fHData->SetStats(0);
+ fHData->SetDirectory(0);
+ fSharingFilter.Init();
+ fHistCollector.Init(*(fHEventsTr->GetXaxis()));
+}
+
+//____________________________________________________________________
+void
+AliForwardMultiplicity::UserCreateOutputObjects()
+{
+ fList = new TList;
+
+ AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
+ AliAODHandler* ah =
+ dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
+ if (!ah)
+ AliFatal("No AOD output handler set in analysis manager");
+
+
+ TObject* obj = &fAODFMD;
+ ah->AddBranch("AliAODForwardMult", &obj);
+
+ // fTree = new TTree("T", "T");
+ // fTree->Branch("forward", &fAODFMD);
+
+ PostData(1, fList);
+ // PostData(2, fTree);
+}
+//____________________________________________________________________
+void
+AliForwardMultiplicity::UserExec(Option_t*)
+{
+ // Get the input data
+ AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
+ if (!esd) {
+ AliWarning("No ESD event found for input event");
+ return;
+ }
+
+#if 0
+ static Int_t nEvents = 0;
+ nEvents++;
+ if (nEvents % 100 == 0) AliInfo(Form("Event # %6d", nEvents));
+#endif
+
+ // On the first event, initialize the parameters
+ if (fFirstEvent) {
+ AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+ pars->SetParametersFromESD(esd);
+ pars->PrintStatus();
+ fFirstEvent = false;
+
+ InitializeSubs();
+ }
+ // Clear stuff
+ fHistos.Clear();
+ fESDFMD.Clear();
+ fAODFMD.Clear();
+
+ // Read trigger information from the ESD and store in AOD object
+ if (!ReadTriggers(esd)) {
+#ifdef VERBOSE
+ AliWarning("Failed to read triggers from ESD");
+#endif
+ return;
+ }
+
+ // Mark this event for storage
+ MarkEventForStore();
+
+ // Check if this is a high-flux event
+ const AliMultiplicity* testmult = esd->GetMultiplicity();
+ if (!testmult) {
+#ifdef VERBOSE
+ AliWarning("No central multiplicity object found");
+#endif
+ return;
+ }
+ Bool_t lowFlux = testmult->GetNumberOfTracklets() < fLowFluxCut;
+
+ // Get the FMD ESD data
+ AliESDFMD* esdFMD = esd->GetFMDData();
+ if (!esdFMD) {
+#ifdef VERBOSE
+ AliWarning("No FMD data found in ESD");
+#endif
+ return;
+ }
+
+ // Get the vertex information
+ Double_t vz = 0;
+ Bool_t vzOk = ReadVertex(esd, vz);
+
+ fHEventsTr->Fill(vz);
+ if (!vzOk) {
+#ifdef VERBOSE
+ AliWarning("Failed to read vertex from ESD");
+#endif
+ return;
+ }
+ fHEventsTrVtx->Fill(vz);
+
+ // Get the vertex bin
+ Int_t ivz = fHEventsTr->GetXaxis()->FindBin(vz)-1;
+ fAODFMD.SetIpZ(vz);
+ if (ivz < 0 || ivz >= fHEventsTr->GetXaxis()->GetNbins()) {
+#if 0
+ AliWarning(Form("Vertex @ %f outside of range [%f,%f]",
+ vz, fHEventsTr->GetXaxis()->GetXmin(),
+ fHEventsTr->GetXaxis()->GetXmax()));
+#endif
+ return;
+ }
+
+ // Apply the sharing filter (or hit merging or clustering if you like)
+ if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, vz)) {
+#ifdef VERBOSE
+ AliWarning("Sharing filter failed!");
+#endif
+ return;
+ }
+
+ // Calculate the inclusive charged particle density
+ if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux)) {
+ AliWarning("Density calculator failed!");
+ return;
+ }
+
+ // Do the secondary and other corrections.
+ if (!fCorrections.Correct(fHistos, ivz)) {
+ AliWarning("Corrections failed");
+ return;
+ }
+
+ if (!fHistCollector.Collect(fHistos, ivz, fAODFMD.GetHistogram())) {
+ AliWarning("Histogram collector failed");
+ return;
+ }
+
+ if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
+ fHData->Add(&(fAODFMD.GetHistogram()));
+}
+
+//____________________________________________________________________
+void
+AliForwardMultiplicity::Terminate(Option_t*)
+{
+ TList* list = dynamic_cast<TList*>(GetOutputData(1));
+ if (!list) {
+ AliError("No output list defined");
+ return;
+ }
+ // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
+ TH1D* dNdeta = fHData->ProjectionX("dNdeta", 1, -1, "e");
+ TH1D* norm = fHData->ProjectionX("dNdeta", 0, 1, "");
+ dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
+ dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
+ dNdeta->Divide(norm);
+ dNdeta->SetStats(0);
+ dNdeta->Scale(Double_t(fHEventsTrVtx->GetEntries())/fHEventsTr->GetEntries(),
+ "width");
+
+ list->Add(fHEventsTr);
+ list->Add(fHEventsTrVtx);
+ list->Add(fHTriggers);
+ list->Add(fHData);
+ list->Add(dNdeta);
+
+ TList* last = new TList;
+ last->SetName("LastEvent");
+ list->Add(last);
+ last->Add(&fAODFMD.GetHistogram());
+ last->Add(fHistos.fFMD1i);
+ last->Add(fHistos.fFMD2i);
+ last->Add(fHistos.fFMD2o);
+ last->Add(fHistos.fFMD3i);
+ last->Add(fHistos.fFMD3o);
+
+
+ fSharingFilter.ScaleHistograms(fHEventsTr->Integral());
+ fSharingFilter.Output(list);
+
+ fDensityCalculator.ScaleHistograms(fHEventsTrVtx->Integral());
+ fDensityCalculator.Output(list);
+
+ fCorrections.ScaleHistograms(fHEventsTrVtx->Integral());
+ fCorrections.Output(list);
+}
+
+//____________________________________________________________________
+void
+AliForwardMultiplicity::MarkEventForStore() const
+{
+ // Make sure the AOD tree is filled
+ AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
+ AliAODHandler* ah =
+ dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
+ if (!ah)
+ AliFatal("No AOD output handler set in analysis manager");
+
+ ah->SetFillAOD(kTRUE);
+}
+//____________________________________________________________________
+Bool_t
+AliForwardMultiplicity::ReadTriggers(AliESDEvent* esd)
+{
+ // Get the analysis manager - should always be there
+ AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
+ if (!am) {
+ AliWarning("No analysis manager defined!");
+ return kFALSE;
+ }
+
+ // Get the input handler - should always be there
+ AliInputEventHandler* ih =
+ static_cast<AliInputEventHandler*>(am->GetInputEventHandler());
+ if (!ih) {
+ AliWarning("No input handler");
+ return kFALSE;
+ }
+
+ // Get the physics selection - add that by using the macro
+ // AddTaskPhysicsSelection.C
+ AliPhysicsSelection* ps =
+ static_cast<AliPhysicsSelection*>(ih->GetEventSelection());
+ if (!ps) {
+ AliWarning("No physics selection");
+ return kFALSE;
+ }
+
+ // Check if this is a collision candidate (INEL)
+ Bool_t inel = ps->IsCollisionCandidate(esd);
+ if (inel) {
+ fAODFMD.SetTriggerBits(AliAODForwardMult::kInel);
+ fHTriggers->Fill(.5);
+ }
+
+
+ // IF this is inel, see if we have a tracklet
+ if (inel) {
+ const AliMultiplicity* spdmult = esd->GetMultiplicity();
+ if (!spdmult) {
+ AliWarning("No SPD multiplicity");
+ }
+ else {
+ Int_t n = spdmult->GetNumberOfTracklets();
+ for (Int_t j = 0; j < n; j++) {
+ if(TMath::Abs(spdmult->GetEta(j)) < 1) {
+ fAODFMD.SetTriggerBits(AliAODForwardMult::kInelGt0);
+ fHTriggers->Fill(1.5);
+ break;
+ }
+ }
+ }
+ }
+
+ // Analyse some trigger stuff
+ AliTriggerAnalysis ta;
+ if (ta.IsOfflineTriggerFired(esd, AliTriggerAnalysis::kNSD1)) {
+ fAODFMD.SetTriggerBits(AliAODForwardMult::kNSD);
+ fHTriggers->Fill(2.5);
+ }
+
+ // Get trigger stuff
+ TString triggers = esd->GetFiredTriggerClasses();
+ if (triggers.Contains("CBEAMB-ABCE-NOPF-ALL")) {
+ fAODFMD.SetTriggerBits(AliAODForwardMult::kEmpty);
+ fHTriggers->Fill(3.5);
+ }
+
+ if (triggers.Contains("CINT1A-ABCE-NOPF-ALL")) {
+ fAODFMD.SetTriggerBits(AliAODForwardMult::kA);
+ fHTriggers->Fill(4.5);
+ }
+
+ if (triggers.Contains("CINT1B-ABCE-NOPF-ALL")) {
+ fAODFMD.SetTriggerBits(AliAODForwardMult::kB);
+ fHTriggers->Fill(5.5);
+ }
+
+
+ if (triggers.Contains("CINT1C-ABCE-NOPF-ALL")) {
+ fAODFMD.SetTriggerBits(AliAODForwardMult::kC);
+ fHTriggers->Fill(6.5);
+ }
+
+ if (triggers.Contains("CINT1-E-NOPF-ALL")) {
+ fAODFMD.SetTriggerBits(AliAODForwardMult::kE);
+ fHTriggers->Fill(7.5);
+ }
+
+ return kTRUE;
+}
+//____________________________________________________________________
+Bool_t
+AliForwardMultiplicity::ReadVertex(AliESDEvent* esd, Double_t& vz)
+{
+ // Get the vertex
+ const AliESDVertex* vertex = esd->GetPrimaryVertexSPD();
+ if (!vertex) {
+#ifdef VERBOSE
+ AliWarning("No SPD vertex found in ESD");
+#endif
+ return kFALSE;
+ }
+
+ // Check that enough tracklets contributed
+ if(vertex->GetNContributors() <= 0) {
+#ifdef VERBOSE
+ AliWarning(Form("Number of contributors to vertex is %d<=0",
+ vertex->GetNContributors()));
+#endif
+ return kFALSE;
+ }
+
+ // Check that the uncertainty isn't too large
+ if (vertex->GetZRes() > 0.1) {
+#ifdef VERBOSE
+ AliWarning(Form("Uncertaintity in Z of vertex is too large %f > 0.1",
+ vertex->GetZRes()));
+#endif
+ return kFALSE;
+ }
+
+ // Get the z coordiante
+ vz = vertex->GetZ();
+
+ return kTRUE;
+}
+
+//____________________________________________________________________
+void
+AliForwardMultiplicity::Print(Option_t*) const
+{
+}
+
+//
+// EOF
+//
--- /dev/null
+#ifndef ALIROOT_PWG2_FORWARD_ALIFORWARDMULTIPLICTY_H
+#define ALIROOT_PWG2_FORWARD_ALIFORWARDMULTIPLICTY_H
+#include <AliAnalysisTaskSE.h>
+#include "AliForwardUtil.h"
+#include "AliFMDSharingFilter.h"
+#include "AliFMDDensityCalculator.h"
+#include "AliFMDCorrections.h"
+#include "AliFMDHistCollector.h"
+#include "AliAODForwardMult.h"
+#include <AliESDFMD.h>
+#include <TH1I.h>
+class AliFMDAnaParameters;
+class AliESDEvent;
+class TH2D;
+class TList;
+class TTree;
+
+
+/**
+ * @mainpage ALICE PWG2 Forward Multiplcity Analysis
+ */
+/**
+ * @defgroup pwg2_forward_analysis PWG2 Forward analysis
+ *
+ * Code to do the multiplicity analysis in the forward psuedo-rapidity
+ * regions
+ *
+ */
+/**
+ * Calculate the multiplicity in the forward regions event-by-event
+ *
+ * @par Inputs:
+ * - AliESDEvent
+ *
+ * @par Outputs:
+ * - AliAODForwardMult
+ *
+ * @par Histograms
+ *
+ * @par Corrections used
+ *
+ * @ingroup pwg2_forward_analysis
+ *
+ */
+class AliForwardMultiplicity : public AliAnalysisTaskSE
+{
+public:
+ /**
+ * Constructor
+ *
+ * @param name Name of task
+ */
+ AliForwardMultiplicity(const char* name);
+ /**
+ * Constructor
+ */
+ AliForwardMultiplicity();
+ /**
+ * Copy constructor
+ *
+ * @param o Object to copy from
+ */
+ AliForwardMultiplicity(const AliForwardMultiplicity& o);
+ /**
+ * Assignment operator
+ *
+ * @param o Object to assign from
+ *
+ * @return Reference to this object
+ */
+ AliForwardMultiplicity& operator=(const AliForwardMultiplicity& o);
+ /**
+ * @{
+ * @name Interface methods
+ */
+ /**
+ * Initialize the task
+ *
+ */
+ virtual void Init();
+ /**
+ * Create output objects
+ *
+ */
+ virtual void UserCreateOutputObjects();
+ /**
+ * Process each event
+ *
+ * @param option Not used
+ */
+ virtual void UserExec(Option_t* option);
+ /**
+ * End of job
+ *
+ * @param option Not used
+ */
+ virtual void Terminate(Option_t* option);
+ /**
+ * @}
+ */
+ void Print(Option_t* option="") const;
+
+ /**
+ * Get reference to the SharingFilter algorithm
+ *
+ * @return Reference to AliFMDSharingFilter object
+ */
+ AliFMDSharingFilter& GetSharingFilter() { return fSharingFilter; }
+ /**
+ * Get reference to the DensityCalculator algorithm
+ *
+ * @return Reference to AliFMDDensityCalculator object
+ */
+ AliFMDDensityCalculator& GetDensityCalculator() { return fDensityCalculator; }
+ /**
+ * Get reference to the Corrections algorithm
+ *
+ * @return Reference to AliFMDCorrections object
+ */
+ AliFMDCorrections& GetCorrections() { return fCorrections; }
+ /**
+ * Get reference to the HistCollector algorithm
+ *
+ * @return Reference to AliFMDHistCollector object
+ */
+ AliFMDHistCollector& GetHistCollector() { return fHistCollector; }
+protected:
+ /**
+ * Initialise the sub objects and stuff. Called on first event
+ *
+ */
+ virtual void InitializeSubs();
+ /**
+ * Mark this event as one to store in the AOD
+ *
+ */
+ virtual void MarkEventForStore() const;
+ /**
+ * Read the trigger information and store in output
+ *
+ * @param esd ESD event object
+ *
+ * @return true in case of success
+ */
+ virtual Bool_t ReadTriggers(AliESDEvent* esd);
+ /**
+ * Read the vertex information, and return the z coordinate
+ *
+ * @param esd ESD event object
+ * @param vz On return, the z coordinate of the primary vertex
+ *
+ * @return true on success
+ */
+ virtual Bool_t ReadVertex(AliESDEvent* esd, Double_t& vz);
+
+ TH1I* fHEventsTr; // Histogram of events w/trigger
+ TH1I* fHEventsTrVtx; // Events w/trigger and vertex
+ TH1I* fHTriggers; // Triggers
+ TH2D* fHData; // Summed 1/Nd^2N_{ch}/dphideta
+ Bool_t fFirstEvent; // Whether the event is the first seen
+ Int_t fLowFluxCut; // Low flux cut
+ AliESDFMD fESDFMD; // Sharing corrected ESD object
+ AliForwardUtil::Histos fHistos; // Cache histograms
+ AliAODForwardMult fAODFMD; // Output object
+
+ AliFMDSharingFilter fSharingFilter; // Algorithm
+ AliFMDDensityCalculator fDensityCalculator; // Algorithm
+ AliFMDCorrections fCorrections; // Algorithm
+ AliFMDHistCollector fHistCollector; // Algorithm
+
+ TList* fList; // Output list
+ TTree* fTree; // Output tree
+
+ ClassDef(AliForwardMultiplicity,1) // Forward multiplicity class
+};
+
+#endif
+// Local Variables:
+// mode: C++
+// End:
+
--- /dev/null
+#include "AliForwardUtil.h"
+#include <TH2D.h>
+#include <TMath.h>
+
+//====================================================================
+AliForwardUtil::Histos::~Histos()
+{
+ if (fFMD1i) delete fFMD1i;
+ if (fFMD2i) delete fFMD2i;
+ if (fFMD2o) delete fFMD2o;
+ if (fFMD3i) delete fFMD3i;
+ if (fFMD3o) delete fFMD3o;
+}
+
+//____________________________________________________________________
+TH2D*
+AliForwardUtil::Histos::Make(UShort_t d, Char_t r,
+ const TAxis& etaAxis) const
+{
+ Int_t ns = (r == 'I' || r == 'i') ? 20 : 40;
+ TH2D* hist = new TH2D(Form("FMD%d%c_cache", d, r),
+ Form("FMD%d%c cache", d, r),
+ etaAxis.GetNbins(), etaAxis.GetXmin(),
+ etaAxis.GetXmax(), ns, 0, 2*TMath::Pi());
+ hist->SetXTitle("#eta");
+ hist->SetYTitle("#phi [radians]");
+ hist->SetZTitle("d^{2}N_{ch}/d#etad#phi");
+ hist->Sumw2();
+ hist->SetDirectory(0);
+
+ return hist;
+}
+//____________________________________________________________________
+void
+AliForwardUtil::Histos::Init(const TAxis& etaAxis)
+{
+ fFMD1i = Make(1, 'I', etaAxis);
+ fFMD2i = Make(2, 'I', etaAxis);
+ fFMD2o = Make(2, 'O', etaAxis);
+ fFMD3i = Make(3, 'I', etaAxis);
+ fFMD3o = Make(3, 'O', etaAxis);
+}
+//____________________________________________________________________
+void
+AliForwardUtil::Histos::Clear(Option_t* option)
+{
+ fFMD1i->Reset(option);
+ fFMD2i->Reset(option);
+ fFMD2o->Reset(option);
+ fFMD3i->Reset(option);
+ fFMD3o->Reset(option);
+}
+
+//____________________________________________________________________
+TH2D*
+AliForwardUtil::Histos::Get(UShort_t d, Char_t r) const
+{
+ switch (d) {
+ case 1: return fFMD1i;
+ case 2: return (r == 'I' || r == 'i' ? fFMD2i : fFMD2o);
+ case 3: return (r == 'I' || r == 'i' ? fFMD3i : fFMD3o);
+ }
+ return 0;
+}
+#if 0
+//____________________________________________________________________
+TH2D*
+AliForwardUtil::Histos::Get(UShort_t d, Char_t r)
+{
+ switch (d) {
+ case 1: return fFMD1i;
+ case 2: return (r == 'I' || r == 'i' ? fFMD2i : fFMD2o);
+ case 3: return (r == 'I' || r == 'i' ? fFMD3i : fFMD3o);
+ }
+ default: return 0;
+}
+#endif
+
+//
+// EOF
+//
--- /dev/null
+#ifndef ALIROOT_PWG2_FORWARD_ALIFORWARDUTIL_H
+#define ALIROOT_PWG2_FORWARD_ALIFORWARDUTIL_H
+#include <TObject.h>
+class TH2D;
+class TAxis;
+
+/**
+ * Utilities used in the forward multiplcity analysis
+ *
+ * @ingroup pwg2_forward_analysis
+ */
+class AliForwardUtil : public TObject
+{
+ public:
+ /**
+ * Structure to hold histograms
+ *
+ * @ingroup pwg2_forward_analysis
+ */
+ struct Histos : public TObject
+ {
+ /**
+ * Constructor
+ *
+ *
+ */
+ Histos() : fFMD1i(0), fFMD2i(0), fFMD2o(0), fFMD3i(0), fFMD3o(0) {}
+ /**
+ * Copy constructor
+ *
+ * @param o Object to copy from
+ */
+ Histos(const Histos& o)
+ : TObject(o),
+ fFMD1i(o.fFMD1i),
+ fFMD2i(o.fFMD2i),
+ fFMD2o(o.fFMD2o),
+ fFMD3i(o.fFMD3i),
+ fFMD3o(o.fFMD3o)
+ {}
+ /**
+ * Assignement operator
+ *
+ * @return Reference to this
+ */
+ Histos& operator=(const Histos&) { return *this;}
+ /**
+ * Destructor
+ */
+ ~Histos();
+ /**
+ * Initialize the object
+ *
+ * @param etaAxis Eta axis to use
+ */
+ void Init(const TAxis& etaAxis);
+ /**
+ * Make a histogram
+ *
+ * @param d Detector
+ * @param r Ring
+ * @param etaAxis Eta axis to use
+ *
+ * @return Newly allocated histogram
+ */
+ TH2D* Make(UShort_t d, Char_t r, const TAxis& etaAxis) const;
+ /**
+ * Clear data
+ *
+ * @param option Not used
+ */
+ void Clear(Option_t* option="");
+ // const TH2D* Get(UShort_t d, Char_t r) const;
+ /**
+ * Get the histogram for a particular detector,ring
+ *
+ * @param d Detector
+ * @param r Ring
+ *
+ * @return Histogram for detector,ring or nul
+ */
+ TH2D* Get(UShort_t d, Char_t r) const;
+ TH2D* fFMD1i; // Histogram for FMD1i
+ TH2D* fFMD2i; // Histogram for FMD2i
+ TH2D* fFMD2o; // Histogram for FMD2o
+ TH2D* fFMD3i; // Histogram for FMD3i
+ TH2D* fFMD3o; // Histogram for FMD3o
+ };
+
+};
+
+#endif
+// Local Variables:
+// mode: C++
+// End:
+
--- /dev/null
+//____________________________________________________________________
+//
+// $Id: Compile.C 30305 2008-12-09 05:45:53Z cholm $
+//
+// Script to compile (using ACLic) and load a script. It sets the
+// include path to contain the relevant directories.
+//
+/**
+ * Compile an FMD script using ACLic
+ *
+ * @param script Script to compile
+ * @param option Compile option
+ *
+ * @ingroup pwg2_forward_analysis_scripts
+ */
+Bool_t
+Compile(const char* script, Option_t* option="g")
+{
+ if (!script || script[0] == '\0') {
+ std::cerr << "No script to compile!" << std::endl;
+ return kFALSE;
+ }
+ gSystem->Load("libANALYSIS.so");
+ gSystem->Load("libANALYSISalice.so");
+ gSystem->Load("libPWG2forward2.so");
+ TString macroPath(gROOT->GetMacroPath());
+ macroPath.Append(":${ALICE_ROOT}/FMD/scripts");
+ gROOT->SetMacroPath(macroPath.Data());
+ gSystem->SetIncludePath("-I`root-config --incdir` "
+ "-I${ALICE_ROOT} "
+ "-I${ALICE_ROOT}/include "
+ "-I${ALICE_ROOT}/PWG2/FORWARD/analysis2 ");
+ Long_t ret = gROOT->ProcessLine(Form(".L %s+%s", script, option));
+ return ret == 0;
+}
+
+//____________________________________________________________________
+//
+// EOF
+//
--- /dev/null
+#include <TH1D.h>
+#include <TH2D.h>
+#include <TCanvas.h>
+#include <TPad.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <TError.h>
+#include <TStyle.h>
+#include <THStack.h>
+#include <TLegend.h>
+#include <TMath.h>
+#include <TMultiGraph.h>
+#include <TGraph.h>
+#include <TGraphErrors.h>
+#include <TLatex.h>
+#include "AliAODForwardMult.h"
+#include "OtherData.C"
+
+/**
+ * Draw the data stored in the AOD
+ *
+ * To use this, do
+ *
+ * @code
+ * Root> .L $ALICE_ROOT/PWG2/FORWARD/analysis2/Compile.C
+ * Root> Compile("$ALICE_ROOT/PWG2/FORWARD/analysis2/DrawRes.C")
+ * Root> DrawRes dr
+ * Root> dr.Run("AliAODs.root",-10,10,5,AliAODForwardMult::kInel,900)
+ * @endcode
+ *
+ * The output is stored in a ROOT file
+ *
+ * See also the script Pass2.C
+ *
+ * @ingroup pwg2_forward_analysis_scripts
+ */
+struct DrawRes
+{
+public:
+ /** AOD tree */
+ TTree* fTree;
+ /** AOD object */
+ AliAODForwardMult* fAOD;
+ /** Output file */
+ TFile* fOut;
+ /** Summed histogram */
+ TH2D* fSum;
+ /** Vertex efficiency */
+ Double_t fVtxEff;
+ /** Title to put on the plot */
+ TString fTitle;
+
+ //__________________________________________________________________
+ /**
+ * Constructor
+ *
+ */
+ DrawRes()
+ : fTree(0),
+ fAOD(0),
+ fOut(0),
+ fSum(0),
+ fVtxEff(0),
+ fTitle("")
+ {}
+ //__________________________________________________________________
+ /**
+ * Reset internal structures
+ *
+ */
+ void Clear(Option_t* )
+ {
+ if (fTree && fTree->GetCurrentFile()) {
+ fTree->GetCurrentFile()->Close();
+ delete fTree;
+ }
+ if (fOut) {
+ fOut->Close();
+ delete fOut;
+ }
+ if (fSum) delete fSum;
+ fTree = 0;
+ fOut = 0;
+ fSum = 0;
+ fVtxEff = 0;
+
+ }
+ //__________________________________________________________________
+ /**
+ * Run
+ *
+ * @param file Input file with AOD tree
+ * @param vzMin Minimum interaction point z coordinate
+ * @param vzMax Maximum interaction point z coordinate
+ * @param rebin How many times to re-bin the @f$ dN_{ch}/d\eta@f$
+ * @param mask Trigger mask
+ * @param energy Collision energy @f$ \sqrt{s}@f$ or @f$\sqrt{s_{NN}}@f$
+ * @param title Title to put on the plot
+ *
+ * @return True on success, false otherwise
+ */
+ Bool_t Run(const char* file="AliAODs.root",
+ Double_t vzMin=-10, Double_t vzMax=10, Int_t rebin=1,
+ Int_t mask=AliAODForwardMult::kInel, Int_t energy=900,
+ const char* title="")
+ {
+ TString trgName;
+ if (mask & AliAODForwardMult::kInel) trgName.Append("inel-");
+ if (mask & AliAODForwardMult::kInelGt0) trgName.Append("inelgt0-");
+ if (mask & AliAODForwardMult::kNSD) trgName.Append("nsd-");
+ if (trgName.EndsWith("-")) trgName.Remove(trgName.Length()-1);
+ if (trgName.IsNull()) trgName = "unknown";
+ TString outName =
+ TString::Format("dndeta_%04dGeV_%c%02d-%c%02dcm_rb%02d_%s.root",
+ energy,
+ vzMin < 0 ? 'm' : 'p', Int_t(TMath::Abs(vzMin)+.5),
+ vzMax < 0 ? 'm' : 'p', Int_t(TMath::Abs(vzMax)+.5),
+ rebin, trgName.Data());
+ fTitle = title;
+ if (!Open(file, outName)) return kFALSE;
+ if (!Process(vzMin,vzMax,mask)) return kFALSE;
+ if (!Finish(rebin, mask, energy)) return kFALSE;
+
+ return kTRUE;
+ }
+ //__________________________________________________________________
+ /**
+ * Open the files @a fname containg the tree with AliAODEvent objects,
+ * and also open the output file @a outname for writting.
+ *
+ * @param fname Input file name
+ * @param outname Output file name
+ *
+ * @return true on success, false otherwise
+ */
+ Bool_t Open(const char* fname, const char* outname)
+ {
+ Clear("");
+
+ TFile* file = TFile::Open(fname, "READ");
+ if (!file) {
+ Error("Open", "Cannot open file %s", fname);
+ return kFALSE;
+ }
+
+ // Get the AOD tree
+ fTree = static_cast<TTree*>(file->Get("aodTree"));
+ if (!fTree) {
+ Error("Init", "Couldn't get the tree");
+ return kFALSE;
+ }
+
+ // Set the branch pointer
+ fTree->SetBranchAddress("Forward", &fAOD);
+
+ fOut = TFile::Open(outname, "RECREATE");
+ if (!fOut) {
+ Error("Open", "Couldn't open %s", outname);
+ return kFALSE;
+ }
+ return kTRUE;
+ }
+ //__________________________________________________________________
+ /**
+ * Process the events
+ *
+ * @param vzMin Minimum interaction point z coordinate
+ * @param vzMax Maximum interaction point z coordinate
+ * @param mask Trigger mask
+ *
+ * @return true on success, false otherwise
+ */
+ Bool_t Process(Double_t vzMin, Double_t vzMax, Int_t mask)
+ {
+ Int_t nTriggered = 0; // # of triggered ev.
+ Int_t nAccepted = 0; // # of ev. w/vertex
+ Int_t nAvailable = fTree->GetEntries(); // How many entries
+
+ for (int i = 0; i < nAvailable; i++) {
+ fTree->GetEntry(i);
+
+ // Create sum histogram on first event - to match binning to input
+ if (!fSum) {
+ fSum = static_cast<TH2D*>(fAOD->GetHistogram().Clone("d2ndetadphi"));
+ Info("Process", "Created sum histogram (%d,%f,%f)x(%d,%f,%f)",
+ fSum->GetNbinsX(),
+ fSum->GetXaxis()->GetXmin(),
+ fSum->GetXaxis()->GetXmax(),
+ fSum->GetNbinsY(),
+ fSum->GetYaxis()->GetXmin(),
+ fSum->GetYaxis()->GetXmax());
+ }
+
+ // fAOD->Print();
+ // Other trigger/event requirements could be defined
+ if (!fAOD->IsTriggerBits(mask)) continue;
+ nTriggered++;
+
+ // Select vertex range (in centimeters)
+ if (!fAOD->InRange(vzMin, vzMax)) continue;
+ nAccepted++;
+
+ // Add contribution from this event
+ fSum->Add(&(fAOD->GetHistogram()));
+ }
+ fVtxEff = Double_t(nAccepted)/nTriggered;
+
+ Info("Process", "Got %6d triggers, accepted %6d out of %6d events",
+ nTriggered, nAccepted, nAvailable);
+
+ return kTRUE;
+ }
+ //__________________________________________________________________
+ /**
+ * Finish the stuff and draw
+ *
+ * @param rebin How many times to rebin
+ * @param energy Collision energy
+ * @param mask Trigger mask
+ *
+ * @return true on success, false otherwise
+ */
+ Bool_t Finish(Int_t rebin, Int_t mask, Int_t energy)
+ {
+ fOut->cd();
+
+ // Get acceptance normalisation from underflow bins
+ TH1D* norm = fSum->ProjectionX("norm", 0, 1, "");
+ // Project onto eta axis - _ignoring_underflow_bins_!
+ TH1D* dndeta = fSum->ProjectionX("dndeta", 1, -1, "e");
+ dndeta->SetTitle("1/N dN_{ch}/d#eta|_{forward}");
+ // Normalize to the acceptance
+ dndeta->Divide(norm);
+ // Scale by the vertex efficiency
+ dndeta->Scale(fVtxEff, "width");
+ dndeta->SetMarkerColor(kRed+1);
+ dndeta->SetMarkerStyle(20);
+ dndeta->SetMarkerSize(1);
+ dndeta->SetFillStyle(0);
+
+ TH1* sym = Symmetrice(dndeta);
+
+ Rebin(dndeta, rebin);
+ Rebin(sym, rebin);
+
+ THStack* stack = new THStack("results", "Results");
+ stack->Add(dndeta);
+ stack->Add(sym);
+
+ TString hhdName(fOut->GetName());
+ hhdName.ReplaceAll("dndeta", "hhd");
+ TH1* hhd = GetHHD(hhdName.Data(), mask & AliAODForwardMult::kNSD);
+ TH1* hhdsym = 0;
+ if (hhd) {
+ hhdsym = Symmetrice(hhd);
+ stack->Add(hhd);
+ stack->Add(hhdsym);
+ }
+
+ TMultiGraph* mg = GetData(energy, mask);
+
+ gStyle->SetOptTitle(0);
+ Double_t yd = (mg->GetListOfGraphs()->GetEntries() ? 0.3 : 0);
+ TCanvas* c = new TCanvas("c", "C", 900, 800);
+ c->SetFillColor(0);
+ c->SetBorderMode(0);
+ c->SetBorderSize(0);
+ c->cd();
+
+ TPad* p1 = new TPad("p1", "p1", 0, yd, 1.0, 1.0, 0, 0);
+ p1->SetTopMargin(0.05);
+ p1->SetBottomMargin(0.001);
+ p1->SetRightMargin(0.05);
+ p1->SetGridx();
+ p1->SetTicks(1,1);
+ p1->Draw();
+ p1->cd();
+ stack->SetMinimum(-0.1);
+ FixAxis(stack, 1/(1-yd)/1.5, "#frac{1}{N} #frac{dN_{ch}}{#eta}");
+ p1->Clear();
+ stack->DrawClone("nostack e1");
+
+
+ THStack* ratios = new THStack("ratios", "Ratios");
+ TGraphAsymmErrors* o = 0;
+ TIter nextG(mg->GetListOfGraphs());
+ while ((o = static_cast<TGraphAsymmErrors*>(nextG()))) {
+ o->Draw("same p");
+ ratios->Add(Ratio(dndeta, o));
+ ratios->Add(Ratio(sym, o));
+ ratios->Add(Ratio(hhd, o));
+ ratios->Add(Ratio(hhdsym, o));
+ }
+ // Make a legend
+ TString trigs(AliAODForwardMult::GetTriggerString(mask));
+ TLegend* l = p1->BuildLegend(.3,p1->GetBottomMargin()+.01,.7,.4,
+ Form("#sqrt{s}=%dGeV, %s",
+ energy, trigs.Data()));
+ l->SetFillColor(0);
+ l->SetFillStyle(0);
+ l->SetBorderSize(0);
+ p1->cd();
+ c->cd();
+
+ if (!ratios->GetHists() || ratios->GetHists()->GetEntries() <= 0) {
+ p1->SetPad(0, 0, 1, 1);
+ p1->SetBottomMargin(0.1);
+ l->SetY1(0.11);
+ FixAxis(stack, (1-yd)/1, "#frac{1}{N} #frac{dN_{ch}}{#eta}",false);
+ p1->cd();
+ c->cd();
+ }
+ else {
+ TPad* p2 = new TPad("p2", "p2", 0, 0.0, 1.0, yd, 0, 0);
+ p2->SetTopMargin(0.001);
+ p2->SetRightMargin(0.05);
+ p2->SetBottomMargin(1/yd * 0.07);
+ p2->SetGridx();
+ p2->SetTicks(1,1);
+ p2->Draw();
+ p2->cd();
+ FixAxis(ratios, 1/yd/1.5, "Ratios");
+ p2->Clear();
+ ratios->DrawClone("nostack e1");
+
+ TLegend* l2 = p2->BuildLegend(.3,p2->GetBottomMargin()+.01,.7,.6);
+ l2->SetFillColor(0);
+ l2->SetFillStyle(0);
+ l2->SetBorderSize(0);
+
+ p2->cd();
+ TGraphErrors* band = new TGraphErrors(2);
+ band->SetPoint(0, sym->GetXaxis()->GetXmin(), 1);
+ band->SetPoint(1, dndeta->GetXaxis()->GetXmax(), 1);
+ band->SetPointError(0, 0, .1);
+ band->SetPointError(1, 0, .1);
+ band->SetFillColor(kYellow+2);
+ band->SetFillStyle(3002);
+ band->Draw("3 same");
+ ratios->DrawClone("nostack e1 same");
+
+ c->cd();
+ }
+ p1->cd();
+ TLatex* tit = new TLatex(0.10, 0.95, fTitle.Data());
+ tit->SetNDC();
+ tit->SetTextSize(0.05);
+ tit->Draw();
+
+ TString imgName(fOut->GetName());
+ imgName.ReplaceAll(".root", ".png");
+ c->SaveAs(imgName.Data());
+
+ stack->Write();
+ mg->Write();
+ ratios->Write();
+
+ fOut->Close();
+
+ return kTRUE;
+ }
+ //__________________________________________________________________
+ /**
+ * Get the result from previous analysis code
+ *
+ * @param fn File to open
+ *
+ * @return null or result of previous analysis code
+ */
+ TH1* GetHHD(const char* fn="fmd_dNdeta_mult.root", bool nsd=false)
+ {
+ TDirectory* savdir = gDirectory;
+ TFile* file = TFile::Open(fn, "READ");
+ if (!file) {
+ Warning("GetHHD", "couldn't open HHD file %s", fn);
+ savdir->cd();
+ 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->ls();
+ file->Close();
+ savdir->cd();
+ return 0;
+ }
+ TH1* r = static_cast<TH1*>(h->Clone("dndeta_hhd"));
+ r->SetTitle("1/N dN_{ch}/d#eta (HHD)");
+ r->SetFillStyle(0);
+ r->SetFillColor(0);
+ r->SetMarkerStyle(21);
+ r->SetMarkerColor(kPink+1);
+ r->SetDirectory(savdir);
+
+ file->Close();
+ savdir->cd();
+ return r;
+ }
+ //__________________________________________________________________
+ /**
+ * 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
+ */
+ void FixAxis(THStack* stack, Double_t s, const char* ytitle,
+ 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->SetTicks("+-");
+ ya->SetNdivisions(10);
+ ya->SetTitleSize(s*ya->GetTitleSize());
+ ya->SetLabelSize(s*ya->GetLabelSize());
+ }
+ }
+ //__________________________________________________________________
+ /**
+ * 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) 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(h->GetMarkerStyle());
+ ret->SetMarkerColor(g->GetMarkerColor());
+ ret->SetMarkerSize(0.7*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; }
+ return ret;
+ }
+
+ //__________________________________________________________________
+ /**
+ * 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 = new TH1D(Form("%s_mirror", h->GetName()),
+ Form("%s (mirrored)", h->GetTitle()),
+ 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());
+
+ // 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));
+ }
+ return s;
+ }
+ //__________________________________________________________________
+ /**
+ * Rebin a histogram
+ *
+ * @param h Histogram to rebin
+ * @param rebin Rebinning factor
+ *
+ * @return
+ */
+ virtual void Rebin(TH1* h, Int_t rebin) const
+ {
+ if (rebin <= 1) return;
+
+ Int_t nBins = h->GetNbinsX();
+ if(nBins % rebin != 0) {
+ Warning("Rebin", "Rebin factor %d is not a devisor of current number "
+ "of bins %d in the histogram %s", rebin, nBins, h->GetName());
+ return;
+ }
+
+ // Make a copy
+ TH1* tmp = static_cast<TH1*>(h->Clone("tmp"));
+ tmp->Rebin(rebin);
+ tmp->SetDirectory(0);
+
+ // The new number of bins
+ Int_t nBinsNew = nBins / rebin;
+ 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<=rebin;j++) {
+ Int_t bin = (i-1)*rebin + j;
+ if(h->GetBinContent(bin) <= 0) continue;
+ Double_t c = h->GetBinContent(bin);
+ Double_t w = 1 / TMath::Power(c,2);
+ content += c;
+ sumw += w;
+ wsum += w * c;
+ nbins++;
+ }
+
+ if(content > 0 ) {
+ tmp->SetBinContent(i, wsum / sumw);
+ tmp->SetBinError(i,TMath::Sqrt(sumw));
+ }
+ }
+
+ // Finally, rebin the histogram, and set new content
+ h->Rebin(rebin);
+ for(Int_t i =1;i<=nBinsNew; i++) {
+ h->SetBinContent(i,tmp->GetBinContent(i));
+ // h->SetBinError(i,tmp->GetBinError(i));
+ }
+
+ delete tmp;
+ }
+};
+
+//____________________________________________________________________
+//
+// EOF
+//
+
--- /dev/null
+#include <TGraphAsymmErrors.h>
+#include <TMultiGraph.h>
+#include <TStyle.h>
+#include <TMath.h>
+#include <TCanvas.h>
+#include <TLegend.h>
+
+//____________________________________________________________________
+/**
+ * @defgroup pwg2_forward_analysis_otherdata External data
+ *
+ * @ingroup pwg2_forward_analysis_scripts
+ */
+//____________________________________________________________________
+/**
+ * Values used
+ *
+ * @ingroup pwg2_forward_analysis_otherdata
+ */
+enum {
+ /** Color used for UA5 data */
+ UA5Color = kBlue+1,
+ /** Color used for CMS data */
+ CMSColor = kGreen+1,
+ /** Color used for ALICE data */
+ ALICEColor = kMagenta+1,
+ /** Marker style INEL data */
+ INELStyle = 20,
+ /** Marker style INEL>0 data */
+ INELGtStyle= 22,
+ /** Marker style NSD data */
+ NSDStyle = 21,
+ /** Colour offset for mirror data */
+ MirrorOff = 4
+};
+
+//____________________________________________________________________
+/**
+ * Set graph attributes
+ *
+ * @param g Graph
+ * @param marker Marker style
+ * @param color Marker and line color
+ * @param name Name of graph
+ * @param title Title of graph
+ *
+ * @ingroup pwg2_forward_analysis_otherdata
+ */
+void
+SetGraphAttributes(TGraph* g, Int_t marker, Int_t color,
+ const Char_t* name, const Char_t* title)
+{
+ g->SetName(name);
+ g->SetTitle(title);
+ g->SetMarkerStyle(marker);
+ g->SetMarkerColor(color);
+ g->SetLineColor(color);
+ g->SetFillColor(0);
+ g->SetFillStyle(0);
+ g->GetHistogram()->SetStats(kFALSE);
+ g->GetHistogram()->SetXTitle("#eta");
+ g->GetHistogram()->SetYTitle("#frac{1}{N} #frac{dN_{ch}}{#eta}");
+}
+
+//____________________________________________________________________
+/**
+ * Get the UA5 NSD data for pp at @f$ \sqrt{s} = 900GeV@f$
+ * p7886_d1x1y4 - Z.Phys.C33:1-6,1986.
+ *
+ * @param mirrored Wether to produce the mirrored plot
+ *
+ * @return graph of data
+ *
+ * @ingroup pwg2_forward_analysis_otherdata
+ */
+TGraphAsymmErrors* UA5Nsd(Bool_t mirrored=false)
+{
+ //UA5 data NSD - p7886_d1x1y4 - Z.Phys.C33:1-6,1986.
+ double x[] = { 0.125, 0.375, 0.625, 0.875, 1.125, 1.375, 1.625, 1.875, 2.125,
+ 2.375, 2.625, 2.875, 3.125, 3.375, 3.625, 3.875, 4.125, 4.375,
+ 4.625 };
+ double exm[] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
+ double exp[] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
+ double y[] = { 3.48, 3.38, 3.52, 3.68, 3.71, 3.86, 3.76, 3.66, 3.72,
+ 3.69, 3.56, 3.41, 3.15, 3.09, 2.74, 2.73, 2.32, 1.99, 1.69 };
+ double eym[] = { 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07,
+ 0.07, 0.07, 0.07, 0.07, 0.08, 0.08, 0.09, 0.09, 0.1, 0.13 };
+ double eyp[] = { 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07,
+ 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];
+
+ TGraphAsymmErrors* g = new TGraphAsymmErrors(19,x, y,exm,exp,eym,eyp);
+ TGraphAsymmErrors* gm = new TGraphAsymmErrors(19,xm,y,exm,exp,eym,eyp);
+ SetGraphAttributes(g, NSDStyle, UA5Color,"ua5_nsd", "UA5 NSD");
+ SetGraphAttributes(gm, NSDStyle+MirrorOff, UA5Color,"ua5_nsd_mirrored",
+ "UA5 NSD (mirrored)");
+
+ return (mirrored ? gm : g);
+}
+
+//____________________________________________________________________
+/**
+ * Get the UA5 INEL data for pp at @f$ \sqrt{s} = 900GeV@f$
+ * p7886_d2x1y2 - Z.Phys.C33:1-6,1986.
+ *
+ * @param mirrored Wether to produce the mirrored plot
+ *
+ * @return graph of data
+ *
+ * @ingroup pwg2_forward_analysis_otherdata
+ */
+TGraphAsymmErrors* UA5Inel(Bool_t mirrored=false)
+{
+ //UA5 data INEL - p7886_d2x1y2 - Z.Phys.C33:1-6,1986.
+ double x[] = { 0.125, 0.375, 0.625, 0.875, 1.125, 1.375, 1.625, 1.875, 2.125,
+ 2.375, 2.625, 2.875, 3.125, 3.375, 3.625, 3.875, 4.125, 4.375,
+ 4.625 };
+ double exm[] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
+ double exp[] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
+ double y[] = { 3.14, 3.04, 3.17, 3.33, 3.33, 3.53, 3.46, 3.41, 3.45,
+ 3.39, 3.07, 3.07, 2.93, 2.93, 2.55, 2.48, 2.18, 1.91, 1.52 };
+ double eym[] = { 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07,
+ 0.07, 0.07, 0.07, 0.07, 0.08, 0.08, 0.09, 0.09, 0.1, 0.13 };
+ double eyp[] = { 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07,
+ 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);
+
+ SetGraphAttributes(g, INELStyle, UA5Color, "ua5_inel", "UA5 INEL");
+ SetGraphAttributes(gm, INELStyle+MirrorOff, UA5Color, "ua5_inel_mirrored",
+ "UA5 INEL (mirrored)");
+
+ return (mirrored ? gm : g);
+}
+
+//____________________________________________________________________
+/**
+ * Get the ALICE INEL data in @f$ |\eta|<1.3@f$ for pp at @f$ \sqrt{s}
+ * = 900GeV@f$
+ * p7742_d1x1y1 - Eur.Phys.J.C68:89-108,2010.
+ *
+ * @return graph of data
+ *
+ * @ingroup pwg2_forward_analysis_otherdata
+ */
+TGraphAsymmErrors* AliceCentralInel900()
+{
+ // INEL - p7742_d1x1y1 - Eur.Phys.J.C68:89-108,2010.
+ TGraphAsymmErrors* g = 0;
+ double x[] = { -1.3, -1.1, -0.9, -0.7, -0.5, -0.3,
+ -0.1, 0.1, 0.3, 0.5, 0.7, 0.9, 1.1, 1.3 };
+ double exm[] = { 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
+ 0.1, 0.1, 0.1, 0.1, 0.1 };
+ double exp[] = { 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
+ 0.1, 0.1, 0.1, 0.1, 0.1 };
+ double y[] = { 3.28, 3.28, 3.22, 3.12, 3.06, 3.02, 2.98, 3.02, 3.02,
+ 3.05, 3.15, 3.21, 3.26, 3.33 };
+ double eym[] = { 0.06324555320336758, 0.06324555320336758,
+ 0.06324555320336758, 0.06324555320336758,
+ 0.06324555320336758, 0.05385164807134505,
+ 0.05385164807134505, 0.05385164807134505,
+ 0.05385164807134505, 0.06324555320336758,
+ 0.06324555320336758, 0.06324555320336758,
+ 0.06324555320336758, 0.06324555320336758 };
+ double eyp[] = { 0.08246211251235322, 0.08246211251235322,
+ 0.08246211251235322, 0.08246211251235322,
+ 0.08246211251235322, 0.08246211251235322,
+ 0.07280109889280519, 0.08246211251235322,
+ 0.08246211251235322, 0.08246211251235322,
+ 0.08246211251235322, 0.08246211251235322,
+ 0.08246211251235322, 0.08246211251235322 };
+ const int np = 14;
+ for (int i = 0; i < np; i++) {
+ eym[i] += 0.02;
+ eyp[i] += 0.02;
+ }
+ g = new TGraphAsymmErrors(np, x, y, exm, exp, eym, eyp);
+ SetGraphAttributes(g, INELStyle, ALICEColor, "alice_inel",
+ "ALICE INEL (publ.)");
+
+ return g;
+}
+
+//____________________________________________________________________
+/**
+ * Get the ALICE INEL>0 data in @f$ |\eta|<1.3@f$ for pp at @f$
+ * \sqrt{s} = 900GeV@f$
+ *
+ * p7741_d4x1y1 - Eur.Phys.J.C68:345-354,2010.
+ *
+ * @return graph of data
+ *
+ * @ingroup pwg2_forward_analysis_otherdata
+ */
+TGraphAsymmErrors* AliceCentralInelGt900()
+{
+ // INEL>0 - p7741_d4x1y1 - Eur.Phys.J.C68:345-354,2010.
+ double x[] = { -0.9, -0.7, -0.5, -0.3, -0.1, 0.1, 0.3, 0.5, 0.7,
+ 0.9 };
+ double exm[] = { 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
+ 0.1 };
+ double exp[] = { 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
+ 0.1 };
+ double y[] = { 4.0, 3.87, 3.8, 3.7, 3.67, 3.73, 3.72, 3.77, 3.92,
+ 4.01 };
+ double eym[] = { 0.07615773105863909, 0.07615773105863909,
+ 0.07615773105863909, 0.06324555320336758,
+ 0.06324555320336758, 0.06324555320336758,
+ 0.0670820393249937, 0.07615773105863909,
+ 0.07615773105863909, 0.07615773105863909 };
+ double eyp[] = { 0.08544003745317531, 0.07615773105863909,
+ 0.07615773105863909, 0.07280109889280519,
+ 0.07280109889280519, 0.07280109889280519,
+ 0.07615773105863909, 0.07615773105863909,
+ 0.08544003745317531, 0.08544003745317531 };
+ const int np = 10;
+ for (int i = 0; i < np; i++) {
+ double stat = (i >= 3 && i<=5) ? 0.02 : 0.03;
+ eym[i] += stat;
+ eyp[i] += stat;
+ }
+
+ TGraphAsymmErrors* g = new TGraphAsymmErrors(np, x, y, exm, exp, eym, eyp);
+ SetGraphAttributes(g, INELGtStyle, ALICEColor, "alice_inelgt900",
+ "ALICE INEL>0 (publ.)");
+ return g;
+
+}
+
+//____________________________________________________________________
+/**
+ * Get the ALICE INEL>0 data in @f$ |\eta|<0.9@f$ for pp at @f$
+ * \sqrt{s} = 2.36TeV@f$
+ *
+ * p7741_d5x1y1 - Eur.Phys.J.C68:345-354,2010.
+ *
+ * @return graph of data
+ *
+ * @ingroup pwg2_forward_analysis_otherdata
+ */
+TGraphAsymmErrors* AliceCentralInelGt2360()
+{
+ // INEL>0 - p7741_d5x1y1 - Eur.Phys.J.C68:345-354,2010.
+ double x[] = { -0.9, -0.7, -0.5, -0.3, -0.1, 0.1, 0.3, 0.5, 0.7, 0.9 };
+ double exm[] = { 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1 };
+ double exp[] = { 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1 };
+ double y[] = { 4.91, 4.76, 4.63, 4.64, 4.55, 4.55, 4.64, 4.66, 4.82, 4.88 };
+ double eym[] = { 0.08544003745317531, 0.08544003745317531,
+ 0.08544003745317531, 0.08544003745317531,
+ 0.08544003745317531, 0.08544003745317531,
+ 0.08544003745317531, 0.08544003745317531,
+ 0.08544003745317531, 0.08544003745317531 };
+ double eyp[] = { 0.11401754250991379, 0.11401754250991379,
+ 0.1044030650891055, 0.1044030650891055,
+ 0.1044030650891055, 0.1044030650891055,
+ 0.1044030650891055, 0.1044030650891055,
+ 0.11401754250991379, 0.11401754250991379 };
+ const int np = 10;
+ for (int i = 0; i < np; i++) {
+ double stat = 0.3;
+ eym[i] += stat;
+ eyp[i] += stat;
+ }
+
+ TGraphAsymmErrors* g = new TGraphAsymmErrors(np, x, y, exm, exp, eym, eyp);
+ SetGraphAttributes(g, INELGtStyle, ALICEColor, "alice_inelgt2360",
+ "ALICE INEL>0 (publ.)");
+ return g;
+}
+
+//____________________________________________________________________
+/**
+ * Get the ALICE INEL>0 data in @f$ |\eta|<0.9@f$ for pp at @f$
+ * \sqrt{s} = 7TeV@f$
+ *
+ * p7741_d6x1y1 - Eur.Phys.J.C68:345-354,2010.
+ *
+ * @return graph of data
+ *
+ * @ingroup pwg2_forward_analysis_otherdata
+ */
+TGraphAsymmErrors* AliceCentralInelGt7000()
+{
+ // INEL>0 - p7741_d6x1y1 - Eur.Phys.J.C68:345-354,2010.
+// Plot: p7741_d6x1y1
+ double x[] = { -0.9, -0.7, -0.5, -0.3, -0.1, 0.1, 0.3, 0.5, 0.7, 0.9 };
+ double exm[] = { 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1 };
+ double exp[] = { 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1 };
+ double y[] = { 6.22, 6.07, 6.01, 5.84, 5.85, 5.85, 5.91, 6.01, 6.17, 6.26 };
+ double eym[] = { 0.1216552506059644, 0.1216552506059644,
+ 0.1216552506059644, 0.11180339887498948,
+ 0.11180339887498948, 0.11180339887498948,
+ 0.11180339887498948, 0.1216552506059644,
+ 0.1216552506059644, 0.1216552506059644 };
+ double eyp[] = { 0.21095023109728983, 0.21095023109728983,
+ 0.2009975124224178, 0.2009975124224178,
+ 0.2009975124224178, 0.2009975124224178,
+ 0.2009975124224178, 0.2009975124224178,
+ 0.21095023109728983, 0.21095023109728983 };
+ const int np = 10;
+ for (int i = 0; i < np; i++) {
+ double stat = 0.2;
+ eym[i] += stat;
+ eyp[i] += stat;
+ }
+
+ TGraphAsymmErrors* g = new TGraphAsymmErrors(np, x, y, exm, exp, eym, eyp);
+ SetGraphAttributes(g, INELGtStyle, ALICEColor, "alice_inelgt7000",
+ "ALICE INEL>0 (publ.)");
+ return g;
+}
+
+//____________________________________________________________________
+/**
+ * Get the ALICE NSD data in @f$ |\eta|<1.3@f$ for pp at @f$
+ * \sqrt{s} = 900GeV@f$
+ *
+ * p7742_d2x1y1 - Eur.Phys.J.C68:89-108,2010.
+ *
+ * @return graph of data
+ *
+ * @ingroup pwg2_forward_analysis_otherdata
+ */
+TGraphAsymmErrors* AliceCentralNsd900()
+{
+ //NSD - p7742_d2x1y1 - Eur.Phys.J.C68:89-108,2010.
+ double x[] = { -1.3, -1.1, -0.9, -0.7, -0.5, -0.3, -0.1, 0.1, 0.3,
+ 0.5, 0.7, 0.9, 1.1, 1.3 };
+ double exm[] = { 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
+ 0.1, 0.1, 0.1, 0.1, 0.1 };
+ double exp[] = { 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
+ 0.1, 0.1, 0.1, 0.1, 0.1 };
+ double y[] = { 3.9, 3.89, 3.81, 3.7, 3.64, 3.59, 3.53, 3.58, 3.59,
+ 3.61, 3.74, 3.8, 3.87, 3.95 };
+ double eym[] = { 0.13341664064126335, 0.13152946437965907,
+ 0.13152946437965907, 0.1216552506059644,
+ 0.1216552506059644, 0.1216552506059644,
+ 0.1216552506059644, 0.1216552506059644,
+ 0.1216552506059644, 0.1216552506059644,
+ 0.1216552506059644, 0.13152946437965907,
+ 0.13152946437965907, 0.13341664064126335 };
+ double eyp[] = { 0.13341664064126335, 0.13152946437965907,
+ 0.13152946437965907, 0.1216552506059644,
+ 0.1216552506059644, 0.1216552506059644,
+ 0.1216552506059644, 0.1216552506059644,
+ 0.1216552506059644, 0.1216552506059644,
+ 0.1216552506059644, 0.13152946437965907,
+ 0.13152946437965907, 0.13341664064126335 };
+ const int np = 14;
+ for (int i = 0; i < np; i++) {
+ double stat = (i == 0 || i == np-1) ? 0.03 : 0.02;
+ eym[i] += stat;
+ eyp[i] += stat;
+ }
+
+ TGraphAsymmErrors* g = new TGraphAsymmErrors(np, x, y, exm, exp, eym, eyp);
+ SetGraphAttributes(g, NSDStyle, ALICEColor, "alice_nsd", "ALICE NSD (publ.)");
+
+ return g;
+}
+
+//____________________________________________________________________
+/**
+ * Get the ALICE INEL data in @f$ |\eta|<1.3@f$ for pp at @f$
+ * \sqrt{s} = 2.36TeV@f$
+ *
+ * p7742_d3x1y1 - Eur.Phys.J.C68:89-108,2010.
+ *
+ * @return graph of data
+ *
+ * @ingroup pwg2_forward_analysis_otherdata
+ */
+TGraphAsymmErrors* AliceCentralInel2360()
+{
+ // INEL - p7742_d3x1y1 - Eur.Phys.J.C68:89-108,2010.
+ double x[] = { -1.3, -1.1, -0.9, -0.7, -0.5, -0.3, -0.1, 0.1, 0.3,
+ 0.5, 0.7, 0.9, 1.1, 1.3 };
+ double exm[] = { 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
+ 0.1, 0.1, 0.1, 0.1, 0.1 };
+ double exp[] = { 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
+ 0.1, 0.1, 0.1, 0.1, 0.1 };
+ double y[] = { 4.08, 4.01, 4.0, 3.88, 3.77, 3.8, 3.73, 3.71, 3.79,
+ 3.82, 3.94, 4.02, 4.05, 4.16 };
+ double eym[] = { 0.13341664064126335, 0.12369316876852982,
+ 0.12369316876852982, 0.1216552506059644,
+ 0.1216552506059644, 0.1216552506059644,
+ 0.1216552506059644, 0.11180339887498948,
+ 0.1216552506059644, 0.1216552506059644,
+ 0.12369316876852982, 0.12369316876852982,
+ 0.13341664064126335, 0.13341664064126335 };
+ double eyp[] = { 0.2716615541441225, 0.2716615541441225,
+ 0.2716615541441225, 0.260768096208106,
+ 0.25079872407968906, 0.25079872407968906,
+ 0.25079872407968906, 0.25079872407968906,
+ 0.25079872407968906, 0.260768096208106,
+ 0.261725046566048, 0.2716615541441225,
+ 0.2716615541441225, 0.2816025568065745 };
+ const int np = 14;
+ for (int i = 0; i < np; i++) {
+ double stat = (i < 3 || i > np-1-4) ? 0.03 : 0.02;
+ eym[i] += stat;
+ eyp[i] += stat;
+ }
+
+ TGraphAsymmErrors* g = new TGraphAsymmErrors(np, x, y, exm, exp, eym, eyp);
+ SetGraphAttributes(g, NSDStyle, ALICEColor, "alice_inel2360",
+ "ALICE INEL (publ.)");
+ return g;
+}
+
+//____________________________________________________________________
+/**
+ * Get the ALICE NSD data in @f$ |\eta|<1.3@f$ for pp at @f$
+ * \sqrt{s} = 2.36TeV@f$
+ *
+ * p7742_d4x1y1 - Eur.Phys.J.C68:89-108,2010.
+ *
+ * @return graph of data
+ *
+ * @ingroup pwg2_forward_analysis_otherdata
+ */
+TGraphAsymmErrors* AliceCentralNsd2360()
+{
+ // NSD - p7742_d4x1y1 - Eur.Phys.J.C68:89-108,2010.
+ double x[] = { -1.3, -1.1, -0.9, -0.7, -0.5, -0.3, -0.1, 0.1, 0.3,
+ 0.5, 0.7, 0.9, 1.1, 1.3 };
+ double exm[] = { 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
+ 0.1, 0.1, 0.1, 0.1, 0.1 };
+ double exp[] = { 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
+ 0.1, 0.1, 0.1, 0.1, 0.1 };
+ double y[] = { 4.79, 4.72, 4.7, 4.56, 4.44, 4.47, 4.39, 4.37, 4.47,
+ 4.5, 4.64, 4.73, 4.76, 4.9 };
+ double eym[] = { 0.13601470508735444, 0.13341664064126335,
+ 0.13341664064126335, 0.12369316876852982,
+ 0.12369316876852982, 0.12369316876852982,
+ 0.12369316876852982, 0.12369316876852982,
+ 0.12369316876852982, 0.12369316876852982,
+ 0.12369316876852982, 0.13341664064126335,
+ 0.13341664064126335, 0.13341664064126335 };
+ double eyp[] = { 0.18439088914585774, 0.18248287590894657,
+ 0.18248287590894657, 0.1726267650163207,
+ 0.1726267650163207, 0.1726267650163207,
+ 0.16278820596099708, 0.16278820596099708,
+ 0.1726267650163207, 0.1726267650163207,
+ 0.1726267650163207, 0.18248287590894657,
+ 0.18248287590894657, 0.18248287590894657 };
+ const int np = 14;
+
+ for (int i = 0; i < np; i++) {
+ double stat = (i < 1) ? 0.03 : 0.02;
+ eym[i] += stat;
+ eyp[i] += stat;
+ }
+
+ TGraphAsymmErrors* g = new TGraphAsymmErrors(np, x, y, exm, exp, eym, eyp);
+ SetGraphAttributes(g, NSDStyle, ALICEColor, "alice_nsd2360",
+ "ALICE NSD (publ.)");
+ return g;
+}
+
+
+//____________________________________________________________________
+/**
+ * Get the CMS NSD data in @f$ |\eta|<2.25@f$ for pp at @f$
+ * \sqrt{s} = 900GeV@f$
+ *
+ * p7743_d8x1y1
+ *
+ * @return graph of data
+ *
+ * @ingroup pwg2_forward_analysis_otherdata
+ */
+TGraphAsymmErrors* CMSNsd900()
+{
+ // CMS published NSD data - p7743_d8x1y1
+ double x[] = { -2.25, -1.75, -1.25, -0.75, -0.25, 0.25, 0.75, 1.25, 1.75,
+ 2.25 };
+ double exm[] = { 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25 };
+ double exp[] = { 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25 };
+ double y[] = { 3.6, 3.73, 3.62, 3.54, 3.48, 3.48, 3.54, 3.62, 3.73, 3.6 };
+ double eym[] = { 0.13, 0.14, 0.13, 0.13, 0.13, 0.13, 0.13, 0.13, 0.14,0.13 };
+ double eyp[] = { 0.13, 0.14, 0.13, 0.13, 0.13, 0.13, 0.13, 0.13, 0.14, 0.13 };
+ const int np = 10;
+ TGraphAsymmErrors* g = new TGraphAsymmErrors(np, x, y, exm, exp, eym, eyp);
+ SetGraphAttributes(g, NSDStyle, CMSColor, "cms_nsd900", "CMS NSD");
+
+ return g;
+}
+
+
+//____________________________________________________________________
+/**
+ * Get the CMS NSD data in @f$ |\eta|<2.25@f$ for pp at @f$
+ * \sqrt{s} = 2.36GeV@f$
+ *
+ * p7743_d8x1y2
+ *
+ * @return graph of data
+ *
+ * @ingroup pwg2_forward_analysis_otherdata
+ */
+TGraphAsymmErrors* CMSNsd2360()
+{
+ // CMS NSD 2360 - p7743_d8x1y2
+ double x[] = { -2.25, -1.75, -1.25, -0.75, -0.25, 0.25, 0.75,1.25,1.75,2.25 };
+ double exm[] = { 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25 };
+ double exp[] = { 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25 };
+ double y[] = { 4.78, 4.81, 4.66, 4.61, 4.47, 4.47, 4.61, 4.66, 4.81, 4.78 };
+ double eym[] = { 0.17, 0.18, 0.17, 0.17, 0.16, 0.16, 0.17, 0.17, 0.18, 0.17 };
+ double eyp[] = { 0.17, 0.18, 0.17, 0.17, 0.16, 0.16, 0.17, 0.17, 0.18, 0.17 };
+ const int np = 10;
+ TGraphAsymmErrors* g = new TGraphAsymmErrors(np, x, y, exm, exp, eym, eyp);
+ SetGraphAttributes(g, NSDStyle, CMSColor, "cms_nsd2360", "CMS NSD");
+ return g;
+}
+
+
+//____________________________________________________________________
+/**
+ * Get the CMS NSD data in @f$ |\eta|<2.25@f$ for pp at @f$
+ * \sqrt{s} = 7TeV@f$
+ *
+ * p7838_d5x1y1
+ *
+ * @return graph of data
+ *
+ * @ingroup pwg2_forward_analysis_otherdata
+ */
+TGraphAsymmErrors* CMSNsd7000()
+{
+ // CMS NSD 7000 - Plot: p7838_d5x1y1
+ double x[] = { -2.25, -1.75, -1.25, -0.75, -0.25, 0.25, 0.75,1.25,1.75,2.25 };
+ double exm[] = { 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25 };
+ double exp[] = { 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25 };
+ double y[] = { 6.18, 6.26, 6.14, 6.01, 5.78, 5.78, 6.01, 6.14, 6.26, 6.18 };
+ double eym[] = { 0.25, 0.25, 0.24, 0.24, 0.23, 0.23, 0.24, 0.24, 0.25, 0.25 };
+ double eyp[] = { 0.25, 0.25, 0.24, 0.24, 0.23, 0.23, 0.24, 0.24, 0.25, 0.25 };
+ const int np = 10;
+ TGraphAsymmErrors* g = new TGraphAsymmErrors(np, x, y, exm, exp, eym, eyp);
+ SetGraphAttributes(g, NSDStyle, CMSColor, "cms_nsd7000", "CMS NSD");
+ return g;
+}
+
+//____________________________________________________________________
+/**
+ * Get a multi graph of data for a given energy and trigger type
+ *
+ * @param energy Energy in GeV (900, 2360, 7000)
+ * @param type Bit pattern of trigger type
+ * - 0x1 INEL
+ * - 0x2 INEL>0
+ * - 0x4 NSD
+ *
+ * @return A multi graph with the selected data.
+ *
+ * @ingroup pwg2_forward_analysis_otherdata
+ */
+TMultiGraph*
+GetData(Int_t energy, Int_t type)
+{
+ TMultiGraph* mp = new TMultiGraph(Form("dndeta_%dGeV_%d", energy, type),"");
+ TString tn;
+ TString en;
+ if (TMath::Abs(energy-900) < 10) {
+ en.Append(" #sqrt{s}=900GeV");
+ if (type & 0x1) {
+ tn.Append(" INEL");
+ mp->Add(UA5Inel(false));
+ mp->Add(UA5Inel(true));
+ mp->Add(AliceCentralInel900());
+ }
+ if (type & 0x4) {
+ tn.Append(" NSD");
+ mp->Add(UA5Nsd(false));
+ mp->Add(UA5Nsd(true));
+ mp->Add(AliceCentralNsd900());
+ mp->Add(CMSNsd900());
+ }
+ if (type & 0x2) {
+ tn.Append(" INEL>0");
+ mp->Add(AliceCentralInelGt900());
+ }
+ }
+ if (TMath::Abs(energy-2360) < 10) {
+ en.Append(" #sqrt{s}=2.36TeV");
+ if (type & 0x1) {
+ tn.Append(" INEL");
+ mp->Add(AliceCentralInel2360());
+ }
+ if (type & 0x4) {
+ tn.Append(" NSD");
+ mp->Add(AliceCentralNsd2360());
+ mp->Add(CMSNsd2360());
+ }
+ if (type & 0x1) {
+ tn.Append(" INEL>0");
+ mp->Add(AliceCentralInelGt2360());
+ }
+ }
+ if (TMath::Abs(energy-7000) < 10) {
+ en.Append(" #sqrt{s}=7TeV");
+ if (type & 0x1) {
+ tn.Append(" INEL");
+ }
+ if (type & 0x4) {
+ tn.Append(" NSD");
+ mp->Add(CMSNsd7000());
+ }
+ if (type & 0x1) {
+ tn.Append(" INEL>0");
+ mp->Add(AliceCentralInelGt7000());
+ }
+ }
+ mp->SetTitle(Form("1/N dN_{ch}/d#eta, pp(p#bar{p}), %s, %s",
+ en.Data(), tn.Data()));
+
+ return mp;
+}
+
+//____________________________________________________________________
+/**
+ * Plot external data for a given selection of energy and trigger type
+ * (see GetData)
+ *
+ * @param energy Energy in GeV
+ * @param type Trigger type bit mask
+ *
+ * @ingroup pwg2_forward_analysis_otherdata
+ */
+void
+OtherData(Int_t energy=900, Int_t type=0x1)
+{
+ gStyle->SetTitleX(0.1);
+ gStyle->SetTitleY(1.0);
+ gStyle->SetTitleW(0.85);
+ gStyle->SetTitleH(0.05);
+ gStyle->SetTitleBorderSize(0);
+ gStyle->SetTitleTextColor(kWhite);
+ gStyle->SetTitleFillColor(kBlack);
+ gStyle->SetTitleFontSize(0.02);
+
+ gStyle->SetOptTitle(1);
+ gStyle->SetOptStat(0);
+
+ TCanvas* c = new TCanvas("c", "dN/deta", 800, 600);
+ c->SetFillColor(0);
+ c->SetBorderSize(0);
+ c->SetBorderMode(0);
+ c->SetRightMargin(0.05);
+ c->SetTopMargin(0.05);
+
+ TMultiGraph* mp = GetData(energy, type);
+ mp->SetMinimum(0);
+ mp->Draw("ap");
+ if (mp->GetXaxis())
+ mp->GetXaxis()->SetTitle("#eta");
+ if (mp->GetYaxis())
+ mp->GetYaxis()->SetTitle("#frac{1}{N} #frac{dN_{ch}}{#eta}");
+
+ TLegend* l = c->BuildLegend(0.3, 0.15, 0.7, 0.5);
+ l->SetFillColor(0);
+ l->SetBorderSize(0);
+
+ c->cd();
+}
+
+//____________________________________________________________________
+//
+// EOF
+//
--- /dev/null
+/**
+ * Run first pass of the analysis - that is read in ESD and produce AOD
+ *
+ * @param file ESD input file
+ * @param nEvents Number of events to process
+ * @param nCutBins Number of additional bins to cut off
+ * @param correctionCut Threshold for when to use secondary map
+ *
+ * @ingroup pwg2_forward_analysis_scripts
+ */
+void
+Pass1(const char* file="AliESDs.root",
+ Int_t nEvents=1000,
+ Int_t nCutBins=1,
+ Int_t correctionCut=0.1)
+{
+ gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/RunManager.C");
+
+ RunManager(file, kFALSE, nEvents, nCutBins, correctionCut);
+}
+//
+// EOF
+//
--- /dev/null
+/**
+ * Read in AOD and generate @f$ dN/d\eta@f$ for the selected
+ * trigger classes and vertex ranges
+ *
+ * @param file Input file (AOD)
+ * @param triggers Triggers to investigate
+ * @param energy Energy (only used for comparisons)
+ * @param vzMin Minimum interaction point z coordinate
+ * @param vzMax Maximum interaction point z coordinate
+ * @param rebin How many bins to group
+ * @param title Title to put on the plot
+ *
+ * @ingroup pwg2_forward_analysis_scripts
+ */
+void
+Pass2(const char* file="AliAODs.root",
+ const char* triggers="INEL",
+ Int_t energy=900,
+ Double_t vzMin=-10,
+ Double_t vzMax=10,
+ Int_t rebin=5,
+ const char* title="")
+{
+ gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/Compile.C");
+ Compile("$ALICE_ROOT/PWG2/FORWARD/analysis2/DrawRes.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");
+ }
+
+ 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"
+ "--------------------------------------\n",
+ file, vzMin, vzMax, rebin, trgMask, trgs.Data(), energy, title);
+
+ DrawRes dr;
+ dr.Run(file, vzMin, vzMax, rebin, trgMask, energy, title);
+}
+//
+// EOF
+//
+
+
--- /dev/null
+/**
+ * Script to set-up a train
+ *
+ * @param esd ESD file
+ * @param mc Whether to do MC or not
+ * @param nEvents Number of events
+ * @param nCutBins Bins to cut away
+ * @param correctionCut
+ *
+ * @ingroup pwg2_forward_analysis_scripts
+ */
+void RunManager(const char* esd, Bool_t mc=kFALSE, Int_t nEvents=1000,
+ Int_t nCutBins=1, Float_t correctionCut=0.1)
+{
+ gSystem->Load("libVMC");
+ gSystem->Load("libTree");
+
+ gSystem->Load("libSTEERBase");
+
+ gSystem->Load("libESD") ;
+ gSystem->Load("libAOD") ;
+ gSystem->Load("libANALYSIS");
+ gSystem->Load("libANALYSISalice");
+
+ gSystem->Load("libPhysics");
+ gSystem->Load("libPWG0base");
+ gSystem->Load("libPWG0dep");
+ gSystem->Load("libPWG2forward");
+ gSystem->Load("libPWG2forward2");
+
+ //You can expand this chain if you have more data :-)
+ TChain* chain = new TChain("esdTree");
+ chain->Add(esd);
+
+ //Creating the manager and handlers
+ AliAnalysisManager *mgr = new AliAnalysisManager("Analysis Train",
+ "FMD analysis train");
+ AliESDInputHandler *esdHandler = new AliESDInputHandler();
+ esdHandler->SetInactiveBranches("AliESDACORDE "
+ "AliRawDataErrorLogs "
+ "CaloClusters "
+ "Cascades "
+ "EMCALCells "
+ "EMCALTrigger "
+ "Kinks "
+ "Cascades "
+ "MuonTracks "
+ "TrdTracks "
+ "CaloClusters");
+ mgr->SetInputEventHandler(esdHandler);
+
+
+ // Monte Carlo handler
+ // AliMCEventHandler* mcHandler = new AliMCEventHandler();
+ // mgr->SetMCtruthEventHandler(mcHandler);
+ // mcHandler->SetReadTR(readTR);
+
+ // AOD output handler
+ AliAODHandler* aodHandler = new AliAODHandler();
+ mgr->SetOutputEventHandler(aodHandler);
+ aodHandler->SetFillAOD(kTRUE);
+ aodHandler->SetFillAODforRun(kTRUE);
+ aodHandler->SetOutputFileName("AliAODs.root");
+
+ gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/AddTaskFMD.C");
+ gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
+ AddTaskFMD(nCutBins, correctionCut);
+ AddTaskPhysicsSelection(mc, kTRUE, kTRUE);
+
+ // Run the analysis
+
+ TStopwatch t;
+ if (mgr->InitAnalysis()) {
+ mgr->PrintStatus();
+ mgr->SetUseProgressBar(kTRUE);
+ t.Start();
+ mgr->StartAnalysis("local", chain, nEvents);
+ t.Stop();
+ t.Print();
+ }
+}
+//
+// EOF
+//
--- /dev/null
+
+#include "AliFMDHistCollector.h"
+#include <AliESDFMD.h>
+#include <TAxis.h>
+#include <TList.h>
+#include <TMath.h>
+#include "AliFMDAnaParameters.h"
+#include "AliLog.h"
+#include <TH2D.h>
+#include <TH1I.h>
+#include <TProfile.h>
+#include <TProfile2D.h>
+#include <TArrayI.h>
+
+ClassImp(AliFMDHistCollector)
+#if 0
+; // For Emacs
+#endif
+
+//____________________________________________________________________
+AliFMDHistCollector::AliFMDHistCollector()
+ : TNamed(),
+ fList(),
+ fNEvents(0),
+ fNCutBins(0),
+ fCorrectionCut(0),
+ fFirstBins(),
+ fLastBins(),
+ fUseEtaFromData(kFALSE),
+ fEtaNorm(0),
+ fOutput()
+{}
+
+//____________________________________________________________________
+AliFMDHistCollector::AliFMDHistCollector(const char* title)
+ : TNamed("fmdHistCollector", title),
+ fList(),
+ fNEvents(0),
+ fNCutBins(1),
+ fCorrectionCut(0.5),
+ fFirstBins(1),
+ fLastBins(1),
+ fUseEtaFromData(kFALSE),
+ fEtaNorm(0),
+ fOutput()
+{
+ fList.SetName(GetName());
+ fOutput.SetName(GetName());
+}
+
+//____________________________________________________________________
+AliFMDHistCollector::AliFMDHistCollector(const AliFMDHistCollector& o)
+ : TNamed(o),
+ fList(),
+ fNEvents(o.fNEvents),
+ fNCutBins(o.fNCutBins),
+ fCorrectionCut(o.fCorrectionCut),
+ fFirstBins(o.fFirstBins),
+ fLastBins(o.fLastBins),
+ fUseEtaFromData(o.fUseEtaFromData),
+ fEtaNorm(o.fEtaNorm),
+ fOutput()
+{
+
+ TObject* obj = 0;
+ TIter nextL(&o.fList);
+ while ((obj = nextL())) fList.Add(obj->Clone());
+ TIter nextO(&o.fOutput);
+ while ((obj = nextO())) fOutput.Add(obj->Clone());
+}
+
+//____________________________________________________________________
+AliFMDHistCollector::~AliFMDHistCollector()
+{
+ fList.Delete();
+ fOutput.Delete();
+ if (fEtaNorm) delete fEtaNorm;
+}
+
+//____________________________________________________________________
+AliFMDHistCollector&
+AliFMDHistCollector::operator=(const AliFMDHistCollector& o)
+{
+ SetName(o.GetName());
+ SetTitle(o.GetTitle());
+ if (fEtaNorm) delete fEtaNorm;
+ fList.Delete();
+ fOutput.Delete();
+
+ fNEvents = o.fNEvents;
+ fNCutBins = o.fNCutBins;
+ fCorrectionCut = o.fCorrectionCut;
+ fFirstBins = o.fFirstBins;
+ fLastBins = o.fLastBins;
+ fUseEtaFromData = o.fUseEtaFromData;
+
+ TObject* obj = 0;
+ TIter nextL(&o.fList);
+ while ((obj = nextL())) fList.Add(obj->Clone());
+ TIter nextO(&o.fOutput);
+ while ((obj = nextO())) fOutput.Add(obj->Clone());
+
+ return *this;
+}
+
+//____________________________________________________________________
+void
+AliFMDHistCollector::Init(const TAxis& vtxAxis, const TAxis& etaAxis)
+{
+ AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+
+ Int_t nVz = vtxAxis.GetNbins();
+ fFirstBins.Set(5*nVz);
+ fLastBins.Set(5*nVz);
+
+ // Find the eta bin ranges
+ for (Int_t iVz = 0; iVz < nVz; iVz++) {
+ AliForwardUtil::Histos* h = new AliForwardUtil::Histos;
+ h->Init(etaAxis);
+ fList.AddAt(h, iVz);
+
+ // Find the first and last eta bin to use for each ring for
+ // each vertex bin. This is instead of using the methods
+ // provided by AliFMDAnaParameters
+ for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
+ UShort_t d = 0;
+ Char_t r = 0;
+ GetDetRing(iIdx, d, r);
+
+ // Get the background object
+ TH2F* bg = pars->GetBackgroundCorrection(d,r,iVz);
+ Int_t nEta = bg->GetNbinsX();
+ Int_t first = nEta+1;
+ Int_t last = 0;
+
+ // Loop over the eta bins
+ for (Int_t ie = 1; ie <= nEta; ie++) {
+ if (bg->GetBinContent(ie,1) < fCorrectionCut) continue;
+
+ // Loop over the phi bins to make sure that we
+ // do not have holes in the coverage
+ bool ok = true;
+ for (Int_t ip = 1; ip <= bg->GetNbinsY(); ip++) {
+ if (bg->GetBinContent(ie,ip) < 0.001) {
+ ok = false;
+ continue;
+ }
+ }
+ if (!ok) continue;
+
+ first = TMath::Min(ie, first);
+ last = TMath::Max(ie, last);
+ }
+
+ // Store the result for later use
+ fFirstBins[iVz*5+iIdx] = first;
+ fLastBins[iVz*5+iIdx] = last;
+ } // for j
+ }
+
+ // Next, using the ranges found above, do the eta acceptance
+ // normalisation.
+ fEtaNorm = new TH1F("etaAcceptance", "Acceptance in eta",
+ etaAxis.GetNbins(),
+ etaAxis.GetXmin(),
+ etaAxis.GetXmax());
+ fEtaNorm->SetDirectory(0);
+ fEtaNorm->SetFillColor(kRed+1);
+ fEtaNorm->SetFillStyle(3001);
+
+ TH1F* tmp = new TH1F("tmp", "tmp",
+ etaAxis.GetNbins(),
+ etaAxis.GetXmin(),
+ etaAxis.GetXmax());
+ tmp->SetDirectory(0);
+ TH2D* bg = new TH2D("bg", "Background map",
+ etaAxis.GetNbins(), etaAxis.GetXmin(),etaAxis.GetXmax(),
+ 20, 0, 2*TMath::Pi());
+ bg->SetDirectory(0);
+
+ Int_t colors[] = { kRed, kPink, kMagenta, kViolet, kBlue,
+ kAzure, kCyan, kTeal, kGreen, kSpring };
+
+
+ // Loop over the vertex bins
+ for (Int_t iVz = 0; iVz < nVz; iVz++) {
+ // Clear cache
+ tmp->Reset();
+ tmp->SetFillColor(colors[iVz % 10] + iVz/10);
+ tmp->SetFillStyle(3001);
+
+ bg->Reset();
+ TList* vzList = new TList;
+ vzList->SetName(Form("vtxbin%02d", iVz));
+ fOutput.AddAt(vzList, iVz);
+
+ // Loop over rings
+ for (Int_t d = 1; d <= 3; d++) {
+ Int_t nq = (d == 1 ? 1 : 2);
+ for (Int_t q = 0; q < nq; q++) {
+ Char_t r = (q == 0 ? 'I' : 'O');
+
+ // Get the first and last bin to use
+ Int_t first, last;
+ GetFirstAndLast(d, r, iVz, first, last);
+
+ TH2F* ibg = static_cast<TH2F*>(pars->GetBackgroundCorrection(d,r,iVz)
+ ->Clone("bg"));
+
+ // Zero out-side of selected range
+ if (q == 1) { ibg->RebinY(2); ibg->Scale(0.5); }
+ for (Int_t iEta = 1; iEta < first; iEta++)
+ for (Int_t iPhi = 1; iPhi <= 20; iPhi++)
+ ibg->SetBinContent(iEta,iPhi,0);
+ for (Int_t iEta = last+1; iEta <= etaAxis.GetNbins(); iEta++)
+ for (Int_t iPhi = 1; iPhi <= 20; iPhi++)
+ ibg->SetBinContent(iEta,iPhi,0);
+ // Add to output
+ bg->Add(ibg);
+ delete ibg;
+
+ // Loop over the eta bins
+ for (Int_t iEta = first; iEta <= last; iEta++) {
+ Int_t overlap = GetOverlap(d,r,iEta,iVz);
+ Float_t oc = tmp->GetBinContent(iEta);
+ Float_t nc = overlap >= 0? 0.5 : 1;
+ tmp->SetBinContent(iEta, oc + nc);
+ }
+ }
+ }
+ vzList->Add(tmp->Clone("etaAcceptance"));
+ vzList->Add(bg->Clone("secondaryMap"));
+ fEtaNorm->Add(tmp);
+ }
+ fEtaNorm->Scale(1. / nVz);
+ fOutput.Add(fEtaNorm);
+ delete tmp;
+ delete bg;
+
+}
+
+//____________________________________________________________________
+Int_t
+AliFMDHistCollector::GetIdx(UShort_t d, Char_t r) const
+{
+ Int_t idx = -1;
+ switch (d) {
+ case 1: idx = 0; break;
+ case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
+ case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
+ }
+ return idx;
+}
+//____________________________________________________________________
+void
+AliFMDHistCollector::GetDetRing(Int_t idx, UShort_t& d, Char_t& r) const
+{
+ d = 0;
+ r = '\0';
+ switch (idx) {
+ case 0: d = 1; r = 'I'; break;
+ case 1: d = 2; r = 'I'; break;
+ case 2: d = 2; r = 'O'; break;
+ case 3: d = 3; r = 'I'; break;
+ case 4: d = 3; r = 'O'; break;
+ }
+}
+
+//____________________________________________________________________
+void
+AliFMDHistCollector::GetFirstAndLast(Int_t idx, Int_t vtxbin,
+ Int_t& first, Int_t& last) const
+{
+ first = 0;
+ last = 0;
+
+ if (idx < 0) return;
+ idx += vtxbin * 5;
+
+ if (idx < 0 || idx >= fFirstBins.GetSize()) return;
+
+ first = fFirstBins.At(idx)+fNCutBins;
+ last = fLastBins.At(idx)-fNCutBins;
+}
+
+//____________________________________________________________________
+Int_t
+AliFMDHistCollector::GetFirst(Int_t idx, Int_t v) const
+{
+ Int_t f, l;
+ GetFirstAndLast(idx,v,f,l);
+ return f;
+}
+
+
+//____________________________________________________________________
+Int_t
+AliFMDHistCollector::GetLast(Int_t idx, Int_t v) const
+{
+ Int_t f, l;
+ GetFirstAndLast(idx,v,f,l);
+ return l;
+}
+
+//____________________________________________________________________
+Int_t
+AliFMDHistCollector::GetOverlap(UShort_t d, Char_t r,
+ Int_t bin, Int_t vtxbin) const
+{
+ // Get the possibly overlapping histogram
+ Int_t other = -1;
+ if (d == 1) {
+ if (bin <= GetLast(2,'I',vtxbin)) other = GetIdx(2,'I');
+ }
+ else if (d == 2 && r == 'I') {
+ if (bin <= GetLast(2, 'O', vtxbin)) other = GetIdx(2, 'O');
+ else if (bin >= GetFirst(1, 'I', vtxbin)) other = GetIdx(1, 'I');
+ }
+ else if (d == 2 && r == 'O') {
+ if (bin >= GetFirst(2, 'I', vtxbin)) other = GetIdx(2,'I');
+ }
+ else if (d == 3 && r == 'O') {
+ if (bin <= GetLast(3, 'I', vtxbin)) other = GetIdx(3, 'I');
+ }
+ else if (d == 3 && r == 'I') {
+ if (bin >= GetFirst(3, 'O', vtxbin)) other = GetIdx(3, 'O');
+ }
+ // AliInfo(Form("FMD%d%c (%d) -> %d", d, r, GetIdx(d,r), other));
+ return other;
+}
+//____________________________________________________________________
+Int_t
+AliFMDHistCollector::GetOverlap(Int_t idx, Int_t bin, Int_t vtxbin) const
+{
+ UShort_t d = 0;
+ Char_t r = '\0';
+ GetDetRing(idx, d, r);
+ if (d == 0 || r == '\0') return 0;
+
+ return GetOverlap(d, r, bin, vtxbin);
+}
+
+
+//____________________________________________________________________
+Bool_t
+AliFMDHistCollector::Collect(AliForwardUtil::Histos& hists,
+ Int_t vtxbin,
+ TH2D& out)
+{
+ if (fList.GetEntries() <= vtxbin) {
+ AliWarning(Form("Vertex bin %d out of range [0,%d]",
+ vtxbin, fList.GetEntries()));
+ return kFALSE;
+ }
+
+ Merge(hists, vtxbin, out);
+ Store(hists, vtxbin);
+
+ return kTRUE;
+}
+
+//____________________________________________________________________
+void
+AliFMDHistCollector::Store(AliForwardUtil::Histos& hists,
+ Int_t vtxbin)
+{
+ AliForwardUtil::Histos* histos =
+ static_cast<AliForwardUtil::Histos*>(fList.At(vtxbin));
+ if (!histos) {
+ AliWarning(Form("No histogram for vertex bin %d", vtxbin));
+ return;
+ }
+ if (!histos->fFMD1i ||
+ !histos->fFMD2i ||
+ !histos->fFMD2o ||
+ !histos->fFMD3i ||
+ !histos->fFMD3o) {
+ AliWarning(Form("Histograms for vertex bin %d not initialised "
+ "(%p,%p,%p,%p,%p)", vtxbin,
+ histos->fFMD1i, histos->fFMD2i, histos->fFMD2o,
+ histos->fFMD3i, histos->fFMD3o));
+ return;
+ }
+ if (!hists.fFMD1i ||
+ !hists.fFMD2i ||
+ !hists.fFMD2o ||
+ !hists.fFMD3i ||
+ !hists.fFMD3o) {
+ AliWarning(Form("Cache histograms not initialised (%p,%p,%p,%p,%p)",
+ hists.fFMD1i, hists.fFMD2i, hists.fFMD2o,
+ hists.fFMD3i, hists.fFMD3o));
+ return;
+ }
+
+ histos->fFMD1i->Add(hists.fFMD1i);
+ histos->fFMD2i->Add(hists.fFMD2i);
+ histos->fFMD2o->Add(hists.fFMD2o);
+ histos->fFMD3i->Add(hists.fFMD3i);
+ histos->fFMD3o->Add(hists.fFMD3o);
+}
+
+
+//____________________________________________________________________
+void
+AliFMDHistCollector::Merge(AliForwardUtil::Histos& hists,
+ Int_t vtxbin,
+ TH2D& out)
+{
+ // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+ Bool_t isProfile = (out.InheritsFrom(TProfile2D::Class()));
+
+ 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');
+ TH2D* h = hists.Get(d,r);
+ TH2D* t = static_cast<TH2D*>(h->Clone(Form("FMD%d%c_tmp",d,r)));
+
+ // Outer rings have better phi segmentation - rebin to same as inner.
+ if (q == 1) t->RebinY(2);
+
+ Int_t first = 0;
+ Int_t last = 0;
+ GetFirstAndLast(d, r, vtxbin, first, last);
+
+ // Now update profile output
+ for (Int_t iEta = first; iEta <= last; iEta++) {
+
+ // Get the possibly overlapping histogram
+ Int_t overlap = GetOverlap(d,r,iEta,vtxbin);
+
+ // Fill eta acceptance for this event into the phi underlow bin
+ Float_t ooc = out.GetBinContent(iEta,0);
+ Float_t noc = overlap >= 0? 0.5 : 1;
+ out.SetBinContent(iEta, 0, ooc + noc);
+
+ for (Int_t iPhi = 1; iPhi <= h->GetNbinsY(); iPhi++) {
+ Double_t c = t->GetBinContent(iEta,iPhi);
+ Double_t e = t->GetBinError(iEta,iPhi);
+
+ // If there's no signal, continue
+ // if (e <= 0) continue;
+ if (c <= 0 || e <= 0) continue;
+
+ if (isProfile) {
+ static_cast<TProfile2D&>(out).Fill(h->GetXaxis()->GetBinCenter(iEta),
+ h->GetYaxis()->GetBinCenter(iPhi),
+ c, e);
+ continue;
+ }
+
+ // If there's no overlapping histogram (ring), then
+ // fill in data and continue to the next phi bin
+ if (overlap < 0) {
+ out.SetBinContent(iEta,iPhi,c);
+ out.SetBinError(iEta,iPhi,e);
+ continue;
+ }
+
+ // Get the current bin content and error
+ Double_t oc = out.GetBinContent(iEta,iPhi);
+ Double_t oe = out.GetBinError(iEta,iPhi);
+
+#define USE_STRAIGHT_MEAN
+// #define USE_STRAIGHT_MEAN_NONZERO
+// #define USE_WEIGHTED_MEAN
+// #define USE_MOST_CERTAIN
+#if defined(USE_STRAIGHT_MEAN)
+ // calculate the average of old value (half the original),
+ // and this value, as well as the summed squared errors
+ // of the existing content (sqrt((e_1/2)^2=sqrt(e_1^2/4)=e_1/2)
+ // and half the error of this.
+ //
+ // So, on the first overlapping histogram we get
+ //
+ // c = c_1 / 2
+ // e = sqrt((e_1 / 2)^2) = e_1/2
+ //
+ // On the second we get
+ //
+ // c' = c_2 / 2 + c = c_2 / 2 + c_1 / 2 = (c_1+c_2)/2
+ // e' = sqrt(e^2 + (e_2/2)^2)
+ // = sqrt(e_1^2/4 + e_2^2/4)
+ // = sqrt(1/4 * (e_1^2+e_2^2))
+ // = 1/2 * sqrt(e_1^2 + e_2^2)
+ out.SetBinContent(iEta,iPhi,oc + c/2);
+ out.SetBinError(iEta,iPhi,TMath::Sqrt(oe*oe+(e*e)/4));
+#elif defined(USE_STRAIGHT_MEAN_NONZERO)
+# define ZERO_OTHER
+ // If there's data in the overlapping histogram,
+ // calculate the average and add the errors in
+ // quadrature.
+ //
+ // If there's no data in the overlapping histogram,
+ // then just fill in the data
+ if (oe > 0) {
+ out.SetBinContent(iEta,iPhi,(oc + c)/2);
+ out.SetBinError(iEta,iPhi,TMath::Sqrt(oe*oe + e*e)/2);
+ }
+ else {
+ out.SetBinContent(iEta,iPhi,c);
+ out.SetBinError(iEta,iPhi,e);
+ }
+#elif defined(USE_WEIGHTED_MEAN)
+ // Calculate the weighted mean
+ Double_t w = 1/(e*e);
+ Double_t sc = w * c;
+ Double_t sw = w;
+ if (oe > 0) {
+ Double_t ow = 1/(oe*oe);
+ sc += ow * oc;
+ sw += ow;
+ }
+ Double_t nc = sc / sw;
+ Double_t ne = TMath::Sqrt(1 / sw);
+ out.SetBinContent(iEta,iPhi,nc,ne);
+#elif defined(USE_MOST_CERTAIN)
+# define ZERO_OTHER
+ if (e < oe) {
+ out.SetBinContent(iEta,iPhi,c);
+ out.SetBinError(iEta,iPhi,e);
+ }
+ else {
+ out.SetBinContent(iEta,iPhi,oc);
+ out.SetBinError(iEta,iPhi,oe);
+ }
+#else
+# error No method for defining content of overlapping bins defined
+#endif
+#if defined(ZERO_OTHER)
+ // Get the content of the overlapping histogram,
+ // and zero the content so that we won't use it
+ // again
+ UShort_t od; Char_t oR;
+ GetDetRing(overlap, od, oR);
+ TH2D* other = hists.Get(od,oR);
+ other->SetBinContent(iEta,iPhi,0);
+ other->SetBinError(iEta,iPhi,0);
+#endif
+ }
+ }
+ // Remove temporary histogram
+ delete t;
+ } // for r
+ } // for d
+
+ // Scale the histogram to the bin size
+ // Int_t nEta = out.GetNbinsX();
+ // Int_t nPhi = out.GetNbinsY();
+ // Double_t rEta = out.GetXaxis()->GetXmax()-out.GetXaxis()->GetXmin();
+ // Double_t rPhi = out.GetYaxis()->GetXmax()-out.GetYaxis()->GetXmin();
+ // out.Scale(1. / (rEta/nEta*rPhi/nPhi));
+ // out.Scale(1., "width"); //
+}
+
+//____________________________________________________________________
+void
+AliFMDHistCollector::ScaleHistograms(const TH1I& nEvents)
+{
+ fNEvents = &nEvents;
+
+#if 0
+ TIter next(&fList);
+ AliForwardUtil::Histos* histos = 0;
+ Int_t i = 1;
+ while ((histos = static_cast<AliForwardUtil::Histos*>(next()))) {
+ Int_t nev = nEvents.GetBinContent(i++);
+ if (nev <= 0) continue;
+ histos->fFMD1i->Scale(1. / nev);
+ histos->fFMD2i->Scale(1. / nev);
+ histos->fFMD2o->Scale(1. / nev);
+ histos->fFMD3i->Scale(1. / nev);
+ histos->fFMD3o->Scale(1. / nev);
+ }
+#endif
+}
+
+//____________________________________________________________________
+void
+AliFMDHistCollector::Output(TList* dir)
+{
+ dir->Add(&fOutput);
+
+ // The axis
+ TAxis* eAxis = static_cast<AliForwardUtil::Histos*>(fList.At(0))
+ ->fFMD1i->GetXaxis();
+ TAxis* vAxis = fNEvents->GetXaxis();
+ Int_t nVtx = vAxis->GetNbins();
+
+ // Create profiles of each vertex bin
+ TList mults;
+ for (Int_t i = 0; i < nVtx; i++) {
+ TProfile* mult = new TProfile(Form("dndeta_vtx%02d", i),
+ Form("dN_{ch}/d#eta %f<v_{z}<%f",
+ fNEvents->GetXaxis()->GetBinLowEdge(i+1),
+ fNEvents->GetXaxis()->GetBinUpEdge(i+1)),
+ eAxis->GetNbins(),eAxis->GetXmin(),
+ eAxis->GetXmax());
+ mult->Sumw2();
+ mults.AddAt(mult, i);
+ fOutput.Add(mult);
+ }
+ // Loop over vertex bins
+ for (Int_t iVtx = 0; iVtx < nVtx; iVtx++) {
+ AliForwardUtil::Histos* histos =
+ static_cast<AliForwardUtil::Histos*>(fList.At(iVtx));
+ // Output list for this vertex
+ TList* l =
+ static_cast<TList*>(fOutput.FindObject(Form("vtxbin%02d", iVtx)));
+
+ // Number of events in this vertex bin
+ Int_t nev = fNEvents->GetBinContent(iVtx+1);
+ if (nev <= 0) {
+ continue;
+ }
+
+ // One of the outputs
+ TProfile* mult = static_cast<TProfile*>(mults.At(iVtx));
+ if (!mult) {
+ AliWarning(Form("No multiplicity for vertex %d", iVtx));
+ continue;
+ }
+
+ // Array of histograms and iterator
+ TH2D* hh[] = { histos->fFMD1i, histos->fFMD2i, histos->fFMD2o,
+ histos->fFMD3i, histos->fFMD3o, 0 };
+ TH2D** pp = hh;
+ Int_t idx = 0;
+ while (*pp) {
+ // (eta,phi) histogram of the ring in this vertex
+ (*pp)->SetName(Form("d2ndetadphi_%s", (*pp)->GetName()));
+ // Projection (sum) over phi of each vertex histogram
+ TH1D* h = (*pp)->ProjectionX(Form("%s_proj", (*pp)->GetName()));
+ // Scale to number of events and bin width
+ h->Scale(1./nev, "width");
+ // add to output
+ l->Add(h);
+
+ Int_t firstNonZero = h->GetNbinsX();
+ Int_t lastNonZero = 0;
+ if (!fUseEtaFromData) {
+ GetFirstAndLast(idx, iVtx, firstNonZero, lastNonZero);
+ // firstNonZero += 2;
+ // lastNonZero -= 2;
+ }
+ else {
+ // Find first and las non-zerio bins
+ for (Int_t j = 1; j <= h->GetNbinsX(); j++) {
+ if (h->GetBinContent(j) != 0) {
+ firstNonZero = TMath::Min(j, firstNonZero);
+ lastNonZero = TMath::Max(j, lastNonZero);
+ }
+ }
+ firstNonZero += fNCutBins;
+ lastNonZero -= fNCutBins;
+ }
+
+ // Fill profile histogram
+ UShort_t d; Char_t r;
+ GetDetRing(idx, d, r);
+ for (Int_t iEta = firstNonZero; iEta <= lastNonZero; iEta++) {
+ Double_t c = h->GetBinContent(iEta);
+ Double_t e = h->GetBinError(iEta);
+ if (e <= 0) {
+ AliWarning(Form("error=%f in bin %d of %s/vtx=%d",
+ e, iEta, h->GetName(), iVtx));
+ continue;
+ }
+ mult->Fill(h->GetBinCenter(iEta), c, e);
+ }
+ idx++;
+ pp++;
+ } // while (*pp)
+
+ // Loop again to output 2D hists
+ pp = hh;
+ while (*pp) {
+ // (*pp)->Scale(1./nev, "width");
+ l->Add(*pp);
+ pp++;
+ }
+ } // for histos
+
+ // Result
+ TProfile* total = new TProfile("dndeta", "#frac{1}{N}#frac{dN_{ch}}{d#eta}",
+ eAxis->GetNbins(),eAxis->GetXmin(),
+ eAxis->GetXmax());
+
+
+ for (Int_t iVtx = 0; iVtx < vAxis->GetNbins(); iVtx++) {
+ TProfile* mult = static_cast<TProfile*>(mults.At(iVtx));
+ total->Add(mult);
+ }
+ fOutput.Add(total);
+}
+
+//____________________________________________________________________
+//
+// EOF
+//
+
+
+
--- /dev/null
+#ifndef ALIROOT_PWG2_FORWARD_ANALYSIS2_ALIFMDHISTCOLLECTOR_H
+#define ALIROOT_PWG2_FORWARD_ANALYSIS2_ALIFMDHISTCOLLECTOR_H
+#include <TNamed.h>
+#include <TList.h>
+#include <TArrayI.h>
+#include "AliForwardUtil.h"
+class AliESDFMD;
+class TH2D;
+class TH1I;
+class TH1F;
+
+/**
+ * This class collects the event histograms into single histograms,
+ * one for each ring in each vertex bin.
+ *
+ * @par Input:
+ * - AliESDFMD object possibly corrected for sharing
+ *
+ * @par Output:
+ * - 5 RingHistos objects - each with a number of vertex dependent
+ * 2D histograms of the inclusive charge particle density
+ *
+ * @par HistCollector used:
+ * - AliFMDAnaCalibBackgroundCorrection
+ *
+ * @ingroup pwg2_forward_analysis
+ */
+class AliFMDHistCollector : public TNamed
+{
+public:
+ /**
+ * Constructor
+ */
+ AliFMDHistCollector();
+ /**
+ * Constructor
+ *
+ * @param name Name of object
+ */
+ AliFMDHistCollector(const char* name);
+ /**
+ * Copy constructor
+ *
+ * @param o Object to copy from
+ */
+ AliFMDHistCollector(const AliFMDHistCollector& o);
+ /**
+ * Destructor
+ */
+ virtual ~AliFMDHistCollector();
+ /**
+ * Assignement operator
+ *
+ * @param o Object to assign from
+ *
+ * @return Reference to this object
+ */
+ AliFMDHistCollector& operator=(const AliFMDHistCollector&);
+ /**
+ * Intialise
+ *
+ * @param vtxAxis Vertex axis
+ * @param etaAxis Eta axis
+ */
+ virtual void Init(const TAxis& vtxAxis, const TAxis& etaAxis);
+ /**
+ * Do the calculations
+ *
+ * @param hists Cache of histograms
+ * @param vtxBin Vertex bin
+ * @param out Output histogram
+ *
+ * @return true on successs
+ */
+ virtual Bool_t Collect(AliForwardUtil::Histos& hists, Int_t vtxBin,
+ TH2D& out);
+ /**
+ * Scale the histograms to the total number of events
+ *
+ * @param nEvents Number of events
+ */
+ void ScaleHistograms(const TH1I& nEvents);
+ /**
+ * Output diagnostic histograms to directory
+ *
+ * @param dir List to write in
+ */
+ void Output(TList* dir);
+ /**
+ * Set the number of extra bins (beyond the secondary map border)
+ * to cut away.
+ *
+ * @param n Number of bins
+ */
+ void SetNCutBins(UInt_t n=2) { fNCutBins = n; }
+ /**
+ * Set the correction cut, that is, when bins in the secondary
+ * correction maps have a value lower than this cut, they are
+ * considered uncertain and not used
+ *
+ * @param cut Cut-off
+ */
+ void SetCorrectionCut(Float_t cut=0.5) { fCorrectionCut = cut; }
+ /**
+ * Whether to use the eta range from the data
+ *
+ * @param use
+ */
+ void UseEtaFromData(Bool_t use=kTRUE) { fUseEtaFromData = use; }
+protected:
+ /**
+ * Add the 5 input histograms to our internal sum of vertex
+ * dependent histograms
+ *
+ * @param hists Result
+ * @param vtxBin Vertex bin
+ */
+ virtual void Store(AliForwardUtil::Histos& hists, Int_t vtxBin);
+ /**
+ * Merge the 5 input histograms into a single histogram
+ *
+ * @param hists Result
+ * @param vtxBin Vertex bin
+ * @param out Output histogram
+ */
+ virtual void Merge(AliForwardUtil::Histos& hists, Int_t vtxBin, TH2D& out);
+ /**
+ * Get the first and last eta bin to use for a given ring and vertex
+ *
+ * @param d Detector
+ * @param r Ring
+ * @param vtxBin Vertex bin
+ * @param first On return, the first eta bin to use
+ * @param last On return, the last eta bin to use
+ */
+ virtual void GetFirstAndLast(UShort_t d, Char_t r, Int_t vtxBin,
+ Int_t& first, Int_t& last) const;
+ /**
+ * Get the first and last eta bin to use for a given ring and vertex
+ *
+ * @param idx Ring index as given by GetIdx
+ * @param vtxBin Vertex bin
+ * @param first On return, the first eta bin to use
+ * @param last On return, the last eta bin to use
+ */
+ virtual void GetFirstAndLast(Int_t idx, Int_t vtxBin,
+ Int_t& first, Int_t& last) const;
+ /**
+ * Get the first eta bin to use for a given ring and vertex
+ *
+ * @param d Detector
+ * @param r Ring
+ * @param v vertex bin
+ *
+ * @return First eta bin to use, or -1 in case of problems
+ */
+ Int_t GetFirst(UShort_t d, Char_t r, Int_t v) const;
+ /**
+ * Get the first eta bin to use for a given ring and vertex
+ *
+ * @param i Ring index as given by GetIdx
+ * @param v vertex bin
+ *
+ * @return First eta bin to use, or -1 in case of problems
+ */
+ Int_t GetFirst(Int_t idx, Int_t v) const;
+ /**
+ * Get the last eta bin to use for a given ring and vertex
+ *
+ * @param d Detector
+ * @param r Ring
+ * @param v vertex bin
+ *
+ * @return Last eta bin to use, or -1 in case of problems
+ */
+ Int_t GetLast(UShort_t d, Char_t r, Int_t v) const;
+ /**
+ * Get the last eta bin to use for a given ring and vertex
+ *
+ * @param i Ring index as given by GetIdx
+ * @param v vertex bin
+ *
+ * @return Last eta bin to use, or -1 in case of problems
+ */
+ Int_t GetLast(Int_t idx, Int_t v) const;
+ /**
+ * Get the detector and ring from the ring index
+ *
+ * @param idx Ring index
+ * @param d On return, the detector or 0 in case of errors
+ * @param r On return, the ring id or '\0' in case of errors
+ */
+ void GetDetRing(Int_t idx, UShort_t& d, Char_t& r) const;
+ /**
+ * Get the ring index from detector number and ring identifier
+ *
+ * @param d Detector
+ * @param r Ring identifier
+ *
+ * @return ring index or -1 in case of problems
+ */
+ Int_t GetIdx(UShort_t d, Char_t r) const;
+ /**
+ * Get the possibly overlapping histogram of eta bin @a e in
+ * detector and ring
+ *
+ * @param d Detector
+ * @param r Ring
+ * @param e Eta bin
+ * @param v Vertex bin
+ *
+ * @return Overlapping histogram index or -1
+ */
+ Int_t GetOverlap(UShort_t d, Char_t r, Int_t e, Int_t v) const;
+ /**
+ * Get the possibly overlapping histogram of eta bin @a e in
+ * detector and ring
+ *
+ * @param i Ring index
+ * @param e Eta bin
+ * @param v Vertex bin
+ *
+ * @return Overlapping histogram index or -1
+ */
+ Int_t GetOverlap(Int_t i, Int_t e, Int_t v) const;
+ /**
+ * Check if there's an overlapping histogram with this eta bin of
+ * the detector and ring
+ *
+ * @param d Detector
+ * @param r Ring
+ * @param e eta bin
+ * @param v Vertex bin
+ *
+ * @return True if there's an overlapping histogram
+ */
+ Bool_t HasOverlap(UShort_t d, Char_t r, Int_t e, Int_t v) const;
+ /**
+ * Check if there's an overlapping histogram with this eta bin of
+ * ring
+ *
+ * @param i Ring index
+ * @param e eta bin
+ * @param v Vertex bin
+ *
+ * @return True if there's an overlapping histogram
+ */
+ Bool_t HasOverlap(Int_t i, Int_t e, Int_t v) const;
+
+
+ TList fList; // List of histogram containers
+ const TH1I* fNEvents; // Reference event histogram
+ Int_t fNCutBins; // Number of additional bins to cut away
+ Float_t fCorrectionCut; // Cut-off on secondary corrections
+ TArrayI fFirstBins; // Array of first eta bins
+ TArrayI fLastBins; // Array of last eta bins
+ Bool_t fUseEtaFromData; // Wether to use the data for the limits
+ TH1F* fEtaNorm; // Normalisation in eta
+ TList fOutput;
+
+ ClassDef(AliFMDHistCollector,1); // Calculate Nch density
+};
+
+//____________________________________________________________________
+inline void
+AliFMDHistCollector::GetFirstAndLast(UShort_t d, Char_t r, Int_t vtxbin,
+ Int_t& first, Int_t& last) const
+{
+ GetFirstAndLast(GetIdx(d,r), vtxbin, first, last);
+}
+//____________________________________________________________________
+inline Int_t
+AliFMDHistCollector::GetFirst(UShort_t d, Char_t r, Int_t v) const
+{
+ return GetFirst(GetIdx(d,r), v);
+}
+//____________________________________________________________________
+inline Int_t
+AliFMDHistCollector::GetLast(UShort_t d, Char_t r, Int_t v) const
+{
+ return GetLast(GetIdx(d, r), v);
+}
+//____________________________________________________________________
+inline Bool_t
+AliFMDHistCollector::HasOverlap(UShort_t d, Char_t r, Int_t e, Int_t v) const
+{
+ return GetOverlap(d,r,e,v) >= 0;
+}
+//____________________________________________________________________
+inline Bool_t
+AliFMDHistCollector::HasOverlap(Int_t i, Int_t e, Int_t v) const
+{
+ return GetOverlap(i,e,v) >= 0;
+}
+
+#endif
+// Local Variables:
+// mode: C++
+// End:
+
--- /dev/null
+#include "DrawResBase.h"
+
+/**
+ * @defgroup pwg2_forward_analysis_scripts PWG2 Forward analysis - scripts
+ *
+ * @ingroup pwg2_forward_analysis
+ */
+/**
+ * Example macro to loop over the event-by-event 2D histogram of
+ * @f[
+ * \frac{d^{2}N_{ch}}{d\eta\,d\phi}
+ * @f]
+ * stored in an AOD.
+ *
+ * The class needs the files <<i>base</i>><tt>_hists.root</tt>
+ * containing the histograms generated by AliForwardMultiplicity and
+ * the file <<i>base</i>><tt>_aods.root</tt> containing the tree
+ * with AliAODEvent objects where AliAODForwardMult objects have been
+ * added to in the branch <tt>Forward</tt>
+ *
+ * @ingroup pwg2_forward_analysis_scripts
+ */
+class DrawRes1D : public DrawResBase
+{
+public:
+ /**
+ * Constructor
+ *
+ * @param special If true, add to the list of 'specials'
+ */
+ DrawRes1D()
+ : fTotal1D(0)
+ {}
+ //__________________________________________________________________
+ virtual void Clear(Option_t* option)
+ {
+ DrawResBase::Clear(option);
+ if (fTotal1D) delete fTotal1D;
+ fTotal1D = 0;
+ }
+ //__________________________________________________________________
+ virtual Bool_t IsInit() const { return fTotal1D; }
+ //__________________________________________________________________
+ /**
+ * Utility function to set-up histograms based on the input
+ * @f$ dd^{2}N_{ch}/d\eta\,d\phi@f$ histogram. This member function
+ * is called on the first event so that we have the proper binning
+ *
+ * @param templ Input histogram
+ *
+ * @return true on succcess
+ */
+ virtual Bool_t FirstEvent(const TH2D& templ)
+ {
+ if (!DrawResBase::FirstEvent(templ)) return kFALSE;
+ const TAxis* etaAxis = templ.GetXaxis();
+
+ // Generate sum histograms.
+ // - fTotal1D will be the sum of projections on the X axis of the input
+ // histograms
+ fTotal1D = new TH1D("dndeta",
+ "#frac{1}{N}#frac{dN_{ch}}{d#eta}",
+ etaAxis->GetNbins(),
+ etaAxis->GetXmin(),
+ etaAxis->GetXmax());
+ fTotal1D->SetMarkerStyle(29);
+ fTotal1D->SetMarkerSize(1);
+ fTotal1D->SetStats(0);
+ fTotal1D->SetDirectory(0);
+ fTotal1D->Sumw2();
+
+ return kTRUE;
+ }
+
+
+ //__________________________________________________________________
+ /**
+ * Process the events
+ *
+ *
+ * @return true on success, false otherwise
+ */
+ virtual void AddContrib(const TH2D& hist, Int_t)
+ {
+ // Make a projection on the X axis of the input histogram
+ TH1D* proj = hist.ProjectionX("_px", 0, -1, "e");
+
+ // Add to 1D summed histogram
+ fTotal1D->Add(proj); // , scale);
+
+ // Remove the projection
+ delete proj;
+ }
+ //__________________________________________________________________
+ /**
+ * Process the events
+ *
+ *
+ * @return true on success, false otherwise
+ */
+ virtual TH1D* GetResult()
+ {
+ // Scale our direct sum of the projects of the input histograms to
+ // the number of vertex bins and the bin width. If we do rebinning,
+ // we must scale it one more time.
+ // fTotal1D->Scale(1. / fNAccepted, "width");
+ fTotal1D->Divide(fNorm);
+ fTotal1D->Scale(1., "width");
+ return fTotal1D;
+ }
+ /**
+ * Browse this object
+ *
+ * @param b Browser to use
+ */
+ virtual void Browse(TBrowser* b)
+ {
+ if (fTotal1D) b->Add(fTotal1D);
+ DrawResBase::Browse(b);
+ }
+ /**
+ * Create a new object of this class, and add it to the list of
+ * specials, and create a browse and browse to this object
+ *
+ * @return Newly created object
+ */
+ static DrawRes1D* Create()
+ {
+ DrawRes1D* dr = new DrawRes1D;
+ gROOT->GetListOfSpecials()->Add(dr);
+ TBrowser* b = new TBrowser("b");
+ b->BrowseObject(gROOT->GetListOfSpecials());
+ return dr;
+ }
+protected:
+ TH1D* fTotal1D; // Direct sum of input histograms
+
+ ClassDef(DrawRes1D,0)
+};
+
+//
+// EOF
+//
--- /dev/null
+#include "DrawResBase.h"
+
+/**
+ * @defgroup pwg2_forward_analysis_scripts PWG2 Forward analysis - scripts
+ *
+ * @ingroup pwg2_forward_analysis
+ */
+/**
+ * Example macro to loop over the event-by-event 2D histogram of
+ * @f[
+ * \frac{d^{2}N_{ch}}{d\eta\,d\phi}
+ * @f]
+ * stored in an AOD.
+ *
+ * The class needs the files <<i>base</i>><tt>_hists.root</tt>
+ * containing the histograms generated by AliForwardMultiplicity and
+ * the file <<i>base</i>><tt>_aods.root</tt> containing the tree
+ * with AliAODEvent objects where AliAODForwardMult objects have been
+ * added to in the branch <tt>Forward</tt>
+ *
+ * @ingroup pwg2_forward_analysis_scripts
+ */
+class DrawRes2D : public DrawResBase
+{
+public:
+ /**
+ * Constructor
+ *
+ * @param special If true, add to the list of 'specials'
+ */
+ DrawRes2D()
+ : fTotal2D(0)
+ {
+ }
+ //__________________________________________________________________
+ virtual void Clear(Option_t* option)
+ {
+ DrawResBase::Clear(option);
+ if (fTotal2D) delete fTotal2D;
+ fTotal2D = 0;
+ }
+ //__________________________________________________________________
+ Bool_t IsInit() const { return fTotal2D; }
+ //__________________________________________________________________
+ /**
+ * Utility function to set-up histograms based on the input
+ * @f$ dd^{2}N_{ch}/d\eta\,d\phi@f$ histogram. This member function
+ * is called on the first event so that we have the proper binning
+ *
+ * @param templ Input histogram
+ *
+ * @return true on succcess
+ */
+ Bool_t FirstEvent(const TH2D& templ)
+ {
+ if (!DrawResBase::FirstEvent(templ)) return kFALSE;
+
+ const TAxis* etaAxis = templ.GetXaxis();
+ const TAxis* phiAxis = templ.GetYaxis();
+
+ // Generate sum histograms.
+ // - fTotal2D will be the direct sum of the input histograms.
+ fTotal2D = new TH2D("d2ndetadphi", "1/N dN^{2}_{ch}/d#etad#phi",
+ etaAxis->GetNbins(),
+ etaAxis->GetXmin(),
+ etaAxis->GetXmax(),
+ phiAxis->GetNbins(),
+ phiAxis->GetXmin(),
+ phiAxis->GetXmax());
+ fTotal2D->SetXTitle("#eta");
+ fTotal2D->SetYTitle("#varphi [radians]");
+ fTotal2D->SetZTitle(fTotal2D->GetTitle());
+ fTotal2D->Sumw2();
+ fTotal2D->SetStats(0);
+ fTotal2D->SetDirectory(0);
+
+ return kTRUE;
+ }
+
+
+ //__________________________________________________________________
+ /**
+ * Process the events
+ *
+ *
+ * @return true on success, false otherwise
+ */
+ void AddContrib(const TH2D& hist, Int_t)
+ {
+ fTotal2D->Add(&hist);
+ }
+ //__________________________________________________________________
+ /**
+ * Process the events
+ *
+ *
+ * @return true on success, false otherwise
+ */
+ TH1D* GetResult()
+ {
+ // Do not sum underflow bins!
+ TH1D* proj = fTotal2D->ProjectionX("dndeta_proj", 1, -1, "e");
+ TH1D* norm = fTotal2D->ProjectionX("norm", 0, 1, "");
+ proj->SetTitle("1/N dN_{ch}/d#eta");
+ proj->Divide(norm);
+ proj->Scale(1., "width");
+ return proj;
+ }
+ /**
+ * Browse this object
+ *
+ * @param b Browser to use
+ */
+ void Browse(TBrowser* b)
+ {
+ if (fTotal2D) b->Add(fTotal2D);
+ DrawResBase::Browse(b);
+ }
+ /**
+ * Create a new object of this class, and add it to the list of
+ * specials, and create a browse and browse to this object
+ *
+ * @return Newly created object
+ */
+ static DrawResBase* Create()
+ {
+ DrawRes2D* dr = new DrawRes2D;
+ gROOT->GetListOfSpecials()->Add(dr);
+ TBrowser* b = new TBrowser("b");
+ b->BrowseObject(gROOT->GetListOfSpecials());
+ return dr;
+ }
+protected:
+ TH2D* fTotal2D; // Direct sum of input histograms
+
+ ClassDef(DrawRes2D,0)
+};
+
+//
+// EOF
+//
--- /dev/null
+#include <TH1D.h>
+#include <TH2D.h>
+#include <TBrowser.h>
+#include <TROOT.h>
+#include <TH1I.h>
+#include <TProfile.h>
+#include <TList.h>
+#include <TAxis.h>
+#include <TCanvas.h>
+#include <TPad.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <TError.h>
+#include <TStyle.h>
+#include <THStack.h>
+#include <TLegend.h>
+#include <TMath.h>
+#include <TParameter.h>
+#include "AliAODForwardMult.h"
+
+/**
+ * @defgroup pwg2_forward_analysis_scripts PWG2 Forward analysis - scripts
+ *
+ * @ingroup pwg2_forward_analysis
+ */
+/**
+ * Example macro to loop over the event-by-event 2D histogram of
+ * @f[
+ * \frac{d^{2}N_{ch}}{d\eta\,d\phi}
+ * @f]
+ * stored in an AOD.
+ *
+ * The class needs the files <<i>base</i>><tt>_hists.root</tt>
+ * containing the histograms generated by AliForwardMultiplicity and
+ * the file <<i>base</i>><tt>_aods.root</tt> containing the tree
+ * with AliAODEvent objects where AliAODForwardMult objects have been
+ * added to in the branch <tt>Forward</tt>
+ *
+ * @ingroup pwg2_forward_analysis_scripts
+ */
+class DrawResBase : public TObject
+{
+public:
+ /**
+ * Constructor
+ *
+ * @param special If true, add to the list of 'specials'
+ */
+ DrawResBase()
+ : fBinVzMin(0),
+ fBinVzMax(0),
+ fTree(0),
+ fAOD(0),
+ fNorm(0),
+ fNorm2(0),
+ fVtx(0),
+ fNAccepted(0),
+ fNTriggered(0),
+ fOut(0),
+ fEtaNorms(),
+ fDirect(0)
+ {
+ }
+ virtual ~DrawResBase()
+ {
+ Clear();
+ }
+ /**
+ * Clear internal structures
+ *
+ */
+ virtual void Clear(Option_t* ="")
+ {
+ // Clear previously created data objects
+ if (fTree && fTree->GetCurrentFile()) {
+ fTree->GetCurrentFile()->Close();
+ delete fTree;
+ }
+ if (fOut) {
+ fOut->Close();
+ delete fOut;
+ }
+
+ if (fAOD) delete fAOD;
+ if (fNorm) delete fNorm;
+ if (fNorm2) delete fNorm2;
+ if (fVtx) delete fVtx;
+ fTree = 0;
+ fAOD = 0;
+ fNorm = 0;
+ fNorm2 = 0;
+ fVtx = 0;
+ fOut = 0;
+
+ fEtaNorms.Clear();
+ fNAccepted = 0;
+ fNTriggered = 0;
+ }
+ //__________________________________________________________________
+ /**
+ * Open the files <<i>base</i>><tt>_hists.root</tt>
+ * containing the histograms generated by AliForwardMultiplicity and
+ * the file <<i>base</i>><tt>_aods.root</tt> containing the tree
+ * with AliAODEvent objects.
+ *
+ * @param base Base name of files
+ * @param vzMin Minimum collision vertex z position to use
+ * @param vzMax Maximum collision vertex z position to use
+ * @param rebin Rebinning factor
+ *
+ * @return true on success, false otherwise
+ */
+ virtual Bool_t Open(const char* base,
+ Double_t vzMin=-10,
+ Double_t vzMax=10,
+ Bool_t getOld=false)
+ {
+ // Set our cuts etc.
+ Double_t tVzMin = vzMin;
+ Double_t tVzMax = vzMax;
+ if (tVzMax < tVzMin && tVzMin < 0) tVzMax = -tVzMin;
+
+ // Clear previously created data objects
+ Clear();
+
+ // Open the AOD file
+ TString fn = TString::Format("%s_aods.root", base);
+ TFile* file = TFile::Open(fn.Data(), "READ");
+ if (!file) {
+ Error("Init", "Couldn't open %s", fn.Data());
+ return kFALSE;
+ }
+
+ // Get the AOD tree
+ fTree = static_cast<TTree*>(file->Get("aodTree"));
+ if (!fTree) {
+ Error("Init", "Couldn't get the tree");
+ return kFALSE;
+ }
+
+ // Set the branch pointer
+ fTree->SetBranchAddress("Forward", &fAOD);
+
+ // Open the histogram file
+ fn = TString::Format("%s_hists.root", base);
+ if (!ReadAuxObjects(fn.Data(), tVzMin, tVzMax, getOld)) {
+ Error("Init", "Failed to read auxillary objects from %s", fn.Data());
+ return kFALSE;
+ }
+
+ // Open the output file
+ fn = TString::Format("%s_out.root", base);
+ fOut = TFile::Open(fn.Data(), "RECREATE");
+ if (!fOut) {
+ Error("Init", "Couldn't open %s", fn.Data());
+ return kFALSE;
+ }
+
+ Info("Open", "Selected vertex bins are [%d,%d]", fBinVzMin, fBinVzMax);
+
+ return kTRUE;
+ }
+ /**
+ * Read in auxillary objects from file, and initialise the acceptance
+ * histograms, and find the vertex bins to use.
+ *
+ * This member function expects the following structure in the file
+ *
+ * PWG2forwardDnDeta TDirectory
+ * +- Forward TList
+ * +- nEventsTrVtx TH1I
+ * +- fmdHistCollector TList
+ * +- vtxbin00 TList
+ * | +- etaAcceptance TH1F
+ * +- vtxbin01 TList
+ * | +- etaAcceptance TH1F
+ * |...
+ *
+ * where nEventsTrVtx has nVtx bins, and the subscript on the vtxbin
+ * lists run from 0 to nVtx-1.
+ *
+ * - nEventsTrVtx is a 1D histogram of events per vertex bin. It is
+ * only used to find the minimum and the maximum vertex bins to
+ * use.
+ * - The etaAcceptance histograms are 1D histograms of the covered
+ * eta bins for each vertex bin.
+ *
+ * A file with the correct output is normally generated when
+ * executing the AliForwardMultiplicity task, but can also be
+ * generated outside of that analysis task. In that case, one has
+ * to be certain, that the secondary map cuts used are the same, and
+ * that the vertex bins are the same as used in the analysis.
+ *
+ * @param filename File to read from
+ * @param vzMin Least vertex Z coordinate to accept
+ * @param vzMax Most vertex Z coordinate to accept
+ *
+ * @return true on succes, error otherwise
+ */
+ virtual Bool_t ReadAuxObjects(const char* filename,
+ Double_t vzMin,
+ Double_t vzMax,
+ Bool_t getOld)
+ {
+ TFile* file = TFile::Open(filename, "READ");
+ if (!file) {
+ Error("Init", "Couldn't open %s", filename);
+ return kFALSE;
+ }
+
+ // Get the list stored in the file
+ TList* forward =
+ static_cast<TList*>(file->Get("PWG2forwardDnDeta/Forward"));
+ if (!forward) {
+ Error("ReadAuxObjects", "Couldn't get forward list");
+ return kFALSE;
+ }
+
+ // Get the list from the collector
+ TList* collect =
+ static_cast<TList*>(forward->FindObject("fmdHistCollector"));
+ if (!collect) {
+ Error("Init", "Couldn't get collector list");
+ return kFALSE;
+ }
+
+ // Get the event (by vertex bin) histogram
+ TH1* events = static_cast<TH1I*>(forward->FindObject("nEventsTrVtx"));
+ if (!events) {
+ Error("ReadAuxObjects", "Couldn't get the event histogram");
+ return kFALSE;
+ }
+ // Find the min/max bins to use based on the cuts given
+ fBinVzMin = events->FindBin(vzMin);
+ fBinVzMax = events->FindBin(vzMax-.0000001);
+
+ // Make the normalisation
+ for (Int_t iVz = fBinVzMin; iVz <= fBinVzMax; iVz++) {
+
+ // Get the acceptance from the input file
+ TList* l = static_cast<TList*>(collect
+ ->FindObject(Form("vtxbin%02d",iVz-1)));
+ if (!l) {
+ Error("ReadAuxObjets", "List vtxbin%02d not found in %s",
+ iVz-1, collect->GetName());
+ continue;
+ }
+ TH1F* ve = static_cast<TH1F*>(l->FindObject("etaAcceptance"));
+ if (!ve){
+ Error("ReadAuxObjects", "No eta acceptance histogram found in "
+ "vtxbin%02d/etaAcceptance", iVz-1);
+ continue;
+ }
+ TH1F* te = static_cast<TH1F*>(ve->Clone(Form("etaAcceptance%0d",iVz)));
+ te->SetDirectory(0);
+ fEtaNorms.AddAtAndExpand(te, iVz);
+ }
+ fEtaNorms.SetOwner(kTRUE);
+
+ Bool_t ret = ReadMoreAuxObjects(file,forward,collect,getOld);
+ file->Close();
+ return ret;
+ }
+ //__________________________________________________________________
+ virtual Bool_t ReadMoreAuxObjects(TFile* file,
+ TList* forward,
+ TList* collect,
+ Bool_t getOld)
+ {
+ if (getOld) {
+ fDirect = static_cast<TH1D*>(collect->FindObject("dndeta"));
+ if (fDirect) {
+ fDirect->SetTitle(Form("%s directly", fDirect->GetTitle()));
+ fDirect->SetLineColor(kBlue+1);
+ fDirect->SetMarkerColor(kBlue+1);
+ fDirect->SetMarkerStyle(21);
+ fDirect->SetMarkerSize(0.8);
+ }
+ else
+ Warning("Finish", "Couldn't get the old dN/deta histogram");
+ }
+
+ return kTRUE;
+ }
+ //__________________________________________________________________
+ /**
+ * Check if the passed vertex bin number [1,nVtxBins] is within our
+ * cut.
+ *
+ * @param bin Vertex bin [1,nVtxBins]
+ *
+ * @return true if within cut, false otherwise
+ */
+ Bool_t IsInsideVtxCut(Int_t bin) const
+ {
+ return bin >= fBinVzMin && bin <= fBinVzMax;
+ }
+ virtual Bool_t IsInit() const = 0;
+ //__________________________________________________________________
+ /**
+ * Utility function to set-up histograms based on the input
+ * @f$ dd^{2}N_{ch}/d\eta\,d\phi@f$ histogram. This member function
+ * is called on the first event so that we have the proper binning
+ *
+ * @param templ Input histogram
+ *
+ * @return true on succcess
+ */
+ virtual Bool_t FirstEvent(const TH2D& templ)
+ {
+ Info("FirstEvent", "Generating histograms");
+
+ const TAxis* etaAxis = templ.GetXaxis();
+ // - fNorm is the normalisation
+ fNorm = new TH1D("norm", "Normalisation",
+ etaAxis->GetNbins(),
+ etaAxis->GetXmin(),
+ etaAxis->GetXmax());
+ fNorm->SetFillColor(kRed);
+ fNorm->SetFillStyle(3001);
+ fNorm->SetXTitle("#eta");
+ fNorm->SetYTitle("Normalisation");
+ fNorm->SetStats(0);
+ fNorm->SetDirectory(0);
+
+ fNorm2 = new TH1D("norm2", "Normalisation",
+ etaAxis->GetNbins(),
+ etaAxis->GetXmin(),
+ etaAxis->GetXmax());
+ fNorm2->SetFillColor(kBlue);
+ fNorm2->SetFillStyle(3003);
+ fNorm2->SetXTitle("#eta");
+ fNorm2->SetYTitle("Normalisation");
+ fNorm2->SetStats(0);
+ fNorm2->SetDirectory(0);
+
+ // - fVtx is the vertex distribution
+ fVtx = new TH1D("vtx", "Events per vertex bin",
+ fBinVzMax-fBinVzMin+1,
+ fBinVzMin, fBinVzMax);
+ fVtx->SetXTitle("v_{z} bin");
+ fVtx->SetYTitle("Events");
+ fVtx->SetDirectory(0);
+ fVtx->SetFillColor(kRed+1);
+ fVtx->SetFillStyle(3001);
+ fVtx->SetStats(0);
+
+ return kTRUE;
+ }
+
+
+ //__________________________________________________________________
+ /**
+ * Process the events
+ *
+ *
+ * @return true on success, false otherwise
+ */
+ virtual Bool_t Process()
+ {
+ // Get the number of events in the tree
+ Int_t nEntries = fTree->GetEntries();
+
+ // Loop over the events in the tree
+ for (Int_t event = 0; event < nEntries; event++) {
+ fTree->GetEntry(event);
+ if ((event+1) % 1000 == 0)
+ Info("Process", "Event # %6d of %d", event+1, nEntries);
+
+ // Get our input histogram
+ const TH2D& hist = fAOD->GetHistogram();
+
+
+ // If fTotal2D is not made yet, do so (first event)
+ if (!IsInit()) {
+ if (!FirstEvent(hist)) {
+ Error("Process", "Failed to initialize on first event");
+ return kFALSE;
+ }
+ }
+
+ // Check the trigger
+ if (!fAOD->IsTriggerBits(AliAODForwardMult::kInel)) {
+ // Info("Process", "Not an INEL event");
+ continue;
+ }
+ fNTriggered++;
+
+ // Get the vertex bin - add 1 as we are using histogram bin
+ // numbers in this class
+ Int_t vtxBin = fAOD->GetVtxBin()+1;
+
+ // Check if we're within vertex cut
+ if (!IsInsideVtxCut(vtxBin)) continue;
+
+ // Increment our vertex event and total accepted counters
+ fVtx->AddBinContent(vtxBin);
+ fNAccepted++;
+
+ for (Int_t iEta = 1; iEta <= hist.GetNbinsX(); iEta++)
+ if (hist.GetBinContent(iEta,0))
+ fNorm2->AddBinContent(iEta);
+
+ AddContrib(hist, vtxBin);
+ }
+ FindNormalization();
+
+ Info("Process", "Got %6d trigger, and accepted %6d out of %6d events",
+ fNTriggered, fNAccepted, nEntries);
+
+ return kTRUE;
+ }
+ //__________________________________________________________________
+ /**
+ * Normalize the result
+ *
+ */
+ void FindNormalization()
+ {
+ fEtaNorms.ls();
+ // Make the normalisation
+ for (Int_t iVz = fBinVzMin; iVz <= fBinVzMax; iVz++) {
+
+ // Get the acceptance from the input file
+ TH1F* ve = static_cast<TH1F*>(fEtaNorms.At(iVz));
+ if (!ve){
+ Error("FindNormalization",
+ "No eta acceptance histogram found for vertex bin %d", iVz);
+ continue;
+ }
+ // Add the vertex-dependent acceptance weighted by the number
+ // of events in this vertex bin
+ fNorm->Add(ve, fVtx->GetBinContent(iVz));
+ }
+ // Set all errors to zero (not used)
+ for (Int_t i = 1; i <= fNorm->GetNbinsX(); i++) {
+ fNorm->SetBinError(i, 0);
+ fNorm2->SetBinError(i, 0);
+ }
+ }
+ //__________________________________________________________________
+ /**
+ * Add a contribtion from @a hist to result
+ *
+ * @param hist Contribution to add
+ * @param vtxBin Vertex bin.
+ */
+ virtual void AddContrib(const TH2D& hist, Int_t vtxBin) = 0;
+
+ //__________________________________________________________________
+ /**
+ * Called at the end of the job.
+ *
+ * It plots the result of the analysis in tree canvases.
+ * - One shows the per-vertex accumalted histograms and compares them
+ * to the per-vertex histograms made in the analysis.
+ * - Another shows the final @f$ dN_{ch}/d\eta@f$ calculated in
+ * different ways and compare those to the @f$ dN_{ch}/d\eta@f$ made
+ * during the analysis.
+ * - The last canvas shows the @f$ dd^{2}N_{ch}/d\eta\,d\phi@f$ histogram.
+ *
+ * @return true on succes, false otherwise
+ */
+ Bool_t Finish(Int_t rebin)
+ {
+ // Set the style
+ gStyle->SetOptTitle(0);
+ gStyle->SetPadColor(0);
+ gStyle->SetPadBorderSize(0);
+ gStyle->SetPadBorderMode(0);
+ gStyle->SetPadRightMargin(0.05);
+ gStyle->SetPadTopMargin(0.05);
+ gStyle->SetPalette(1);
+
+ // Make a stack to 'auto-scale' when plotting more than 1 histogram.
+ THStack* stack = new THStack("results","Results for #frac{dN_{ch}}{d#eta}");
+
+ // Generate the projection from the direct sum of the input histograms,
+ // and scale it to the number of vertex bins and the bin width
+ TH1D* tmp = fNorm;
+ fNorm = fNorm2;
+ fNorm2 = tmp;
+ TH1D* h = GetResult();
+ h->SetMarkerColor(kRed+1);
+ h->SetMarkerStyle(20);
+ h->SetMarkerSize(1);
+ h->Scale(Float_t(fNAccepted)/fNTriggered);
+ TH1D* s = Symmetrice(h);
+ Rebin(h, rebin);
+ Rebin(s, rebin);
+ stack->Add(h);
+ stack->Add(s);
+
+ // Get the result from the analysis and plit that too (after modifying
+ // the style and possible rebinning)
+ Float_t y1 = 0;
+ if (fDirect) {
+ fDirect->Scale(Float_t(fNAccepted)/fNTriggered);
+ Rebin(fDirect, rebin);
+ stack->Add(fDirect, "");
+ y1 = 0.3;
+ }
+
+ // The second canvas
+ TCanvas* cTotal1D = new TCanvas("total", "Result", 800, 800);
+ cTotal1D->SetFillColor(0);
+ cTotal1D->SetBorderMode(0);
+ cTotal1D->SetBorderSize(0);
+ cTotal1D->cd();
+
+ // Draw stack in pad
+ TPad* p1 = new TPad("p1", "p1", 0, y1, 1.0, 1.0, 0, 0);
+ p1->Draw();
+ p1->cd();
+ stack->DrawClone("nostack e1 p");
+ stack->Write();
+
+ // Make a legend
+ TLegend* l = p1->BuildLegend(0.31, 0.15, 0.5, 0.6);
+ l->SetFillColor(0);
+ l->SetBorderSize(0);
+ l->SetTextSize(0.04); // l->GetTextSize()/3);
+ cTotal1D->cd();
+
+ if (fDirect) {
+ // Create a new stack
+ stack = new THStack("ratios", "Ratios to old method");
+
+ // Calculate the ratio of our direct sum of input histograms
+ // result from the analysis and draw it
+ TH1D* ratio = static_cast<TH1D*>(h->Clone("ratio"));
+ ratio->SetDirectory(0);
+ ratio->Divide(fDirect);
+ TH1* z = Symmetrice(ratio);
+ stack->Add(ratio, "");
+ stack->Add(z, "");
+
+ // Generate our second pad and draw
+ TPad* p2 = new TPad("p2", "p2", 0, 0.0, 1.0, y1, 0, 0);
+ p2->Draw();
+ p2->cd();
+ stack->Draw("nostack e1");
+ stack->DrawClone("nostack e1");
+ stack->Write();
+ }
+
+ // update
+ cTotal1D->cd();
+
+ // Generate the last canvas and show the summed input histogram
+ TCanvas* other = new TCanvas("other", "Other", 800, 600);
+ other->SetFillColor(0);
+ other->SetBorderMode(0);
+ other->SetBorderSize(0);
+ other->Divide(2,1);
+
+ other->cd(1);
+ fNorm->DrawCopy("HIST");
+ fNorm->Write();
+
+ other->cd(2);
+ fNorm2->DrawCopy("HIST");
+ fNorm2->Write();
+ // fVtx->DrawCopy();
+ // fVtx->Write();
+
+ // out->Write();
+ fOut->Close();
+ return kTRUE;
+ }
+ //__________________________________________________________________
+ /**
+ * Get the result
+ *
+ * @return The result
+ */
+ virtual TH1D* GetResult() { return 0; }
+
+ //__________________________________________________________________
+ /**
+ * Rebin a histogram
+ *
+ * @param h Histogram to rebin
+ * @param rebin Rebinning factor
+ *
+ * @return
+ */
+ virtual void Rebin(TH1D* h, Int_t rebin) const
+ {
+ if (rebin <= 1) return;
+
+ Int_t nBins = h->GetNbinsX();
+ if(nBins % rebin != 0) {
+ Warning("Rebin", "Rebin factor %d is not a devisor of current number "
+ "of bins %d in the histogram %s", rebin, nBins, h->GetName());
+ return;
+ }
+ if (h->IsA()->InheritsFrom(TProfile::Class())){
+ h->Rebin(rebin);
+ return;
+ }
+
+
+ // Make a copy
+ TH1D* tmp = static_cast<TH1D*>(h->Clone("tmp"));
+ tmp->Rebin(rebin);
+ tmp->SetDirectory(0);
+
+ // The new number of bins
+ Int_t nBinsNew = nBins / rebin;
+ 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<=rebin;j++) {
+ Int_t bin = (i-1)*rebin + j;
+ if(h->GetBinContent(bin) <= 0) continue;
+ Double_t c = h->GetBinContent(bin);
+ Double_t w = 1 / TMath::Power(c,2);
+ content += c;
+ sumw += w;
+ wsum += w * c;
+ nbins++;
+ }
+
+ if(content > 0 ) {
+ tmp->SetBinContent(i, wsum / sumw);
+ tmp->SetBinError(i,TMath::Sqrt(sumw));
+ }
+ }
+
+ // Finally, rebin the histogram, and set new content
+ h->Rebin(rebin);
+ for(Int_t i =1;i<=nBinsNew; i++) {
+ h->SetBinContent(i,tmp->GetBinContent(i));
+ // h->SetBinError(i,tmp->GetBinError(i));
+ }
+
+ delete tmp;
+ }
+ //__________________________________________________________________
+ TH1D* Symmetrice(const TH1D* h) const
+ {
+ Int_t nBins = h->GetNbinsX();
+ TH1D* s = new TH1D(Form("%s_mirror", h->GetName()),
+ Form("%s (mirrored)", h->GetTitle()),
+ nBins,
+ -h->GetXaxis()->GetXmax(),
+ -h->GetXaxis()->GetXmin());
+ s->SetDirectory(0);
+ s->SetMarkerColor(h->GetMarkerColor());
+ s->SetMarkerSize(h->GetMarkerSize());
+ s->SetMarkerStyle(h->GetMarkerStyle()+4);
+
+ // 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);
+ Double_t xlast = h->GetBinCenter(last);
+ Int_t f1 = h->GetXaxis()->FindBin(-xfirst);
+ Int_t f2 = s->GetXaxis()->FindBin(-xlast);
+ Int_t l2 = s->GetXaxis()->FindBin(xfirst);
+ Info("Symmetrice", "Data in [%d,%d], copying data from [%d,%d] to [%d,%d]",
+ first, last, f1, last, f2, l2);
+ for (Int_t i = f1, j=l2; i <= last; i++,j--) {
+ s->SetBinContent(j, h->GetBinContent(i));
+ s->SetBinError(j, h->GetBinError(i));
+ }
+ return s;
+ }
+
+
+ //__________________________________________________________________
+ /**
+ * Run the post-processing.
+ *
+ * This will open the files <<i>base</i>><tt>_hists.root</tt>
+ * containing the histograms generated by AliForwardMultiplicity and
+ * the file <<i>base</i>><tt>_aods.root</tt> containing the
+ * tree with AliAODEvent objects.
+ *
+ * Then it will loop over the events, accepting only INEL events
+ * that have a primary collision point along z within @a vzMin and
+ * @a vzMax.
+ *
+ * After the processing, the result will be shown in 3 canvases,
+ * possibly rebinning the result by the factor @a rebin.
+ *
+ * @param base Base name of files
+ * @param vzMin Minimum collision vertex z position to use
+ * @param vzMax Maximum collision vertex z position to use
+ * @param rebin Rebinning factor
+ *
+ * @return true on success, false otherwise
+ */
+ Bool_t Run(const char* base,
+ Double_t vzMin=-10,
+ Double_t vzMax=10,
+ Int_t rebin=1,
+ Bool_t getOld=false)
+ {
+ if (!Open(base,vzMin,vzMax,getOld)) return kFALSE;
+ if (!Process()) return kFALSE;
+ if (!Finish(rebin)) return kFALSE;
+
+ return kTRUE;
+ } // *MENU*
+ /**
+ * Get name of object
+ *
+ * @return Object name
+ */
+ const char* GetName() const { return "drawRes"; }
+ /**
+ * This is a folder
+ *
+ * @return always true
+ */
+ Bool_t IsFolder() const { return kTRUE; }
+ /**
+ * Browse this object
+ *
+ * @param b Browser to use
+ */
+ virtual void Browse(TBrowser* b)
+ {
+ b->Add(&fEtaNorms, "#eta Acceptances");
+ if (fTree) b->Add(fTree);
+ if (fNorm) b->Add(fNorm);
+ if (fNorm2) b->Add(fNorm2);
+ if (fVtx) b->Add(fVtx);
+ if (fDirect) b->Add(fDirect);
+ b->Add(new TParameter<Int_t> ("binVzMin", fBinVzMin));
+ b->Add(new TParameter<Int_t> ("binVzMax", fBinVzMax));
+ b->Add(new TParameter<Int_t> ("fNAccepted", fNAccepted));
+ b->Add(new TParameter<Int_t> ("fNTriggered", fNTriggered));
+ }
+protected:
+ Int_t fBinVzMin; // Corresponding bin to min vertex
+ Int_t fBinVzMax; // Corresponding bin to max vertex
+ TTree* fTree; // Event tree
+ AliAODForwardMult* fAOD; // Our event-by-event data
+ TH1D* fNorm; // Histogram of # events with data per bin
+ TH1D* fNorm2; // Histogram of # events with data per bin
+ TH1D* fVtx; // Histogram of # events with data per bin
+ Int_t fNAccepted; // # of accepted events with vertex
+ Int_t fNTriggered;// # of events with trigger
+ TFile* fOut; // Output file
+ TObjArray fEtaNorms; // Array of eta acceptances per vertex
+ TH1D* fDirect;
+
+ ClassDef(DrawResBase,0)
+};
+
+// Local Variables:
+// mode: C++
+// End:
--- /dev/null
+#include <TH1D.h>
+#include <TH2D.h>
+#include <TH1I.h>
+#include <TProfile.h>
+#include <TList.h>
+#include <TAxis.h>
+#include <TCanvas.h>
+#include <TPad.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <TError.h>
+#include <TStyle.h>
+#include <THStack.h>
+#include <TLegend.h>
+#include <TMath.h>
+#include "AliAODForwardMult.h"
+
+/**
+ * @defgroup pwg2_forward_analysis_scripts PWG2 Forward analysis - scripts
+ *
+ * @ingroup pwg2_forward_analysis
+ */
+/**
+ * Example macro to loop over the event-by-event 2D histogram of
+ * @f[
+ * \frac{d^{2}N_{ch}}{d\eta\,d\phi}
+ * @f]
+ * stored in an AOD.
+ *
+ * The class needs the files <<i>base</i>><tt>_hists.root</tt>
+ * containing the histograms generated by AliForwardMultiplicity and
+ * the file <<i>base</i>><tt>_aods.root</tt> containing the tree
+ * with AliAODEvent objects where AliAODForwardMult objects have been
+ * added to in the branch <tt>Forward</tt>
+ *
+ * @ingroup pwg2_forward_analysis_scripts
+ */
+class DrawRes
+{
+public:
+ /**
+ * Constructor
+ *
+ * @param special If true, add to the list of 'specials'
+ */
+ DrawRes(Bool_t special=true)
+ : fVzMin(-10),
+ fVzMax(10),
+ fBinVzMin(0),
+ fBinVzMax(0),
+ fRebin(1),
+ fTree(0),
+ fAOD(0),
+ fEvents(0),
+ fCollect(0),
+ fByVtx1D(0),
+ fTotal1D(0),
+ fTotal2D(0),
+ fFromVtx(0),
+ fNorm(0)
+ {}
+ //__________________________________________________________________
+ /**
+ * Open the files <<i>base</i>><tt>_hists.root</tt>
+ * containing the histograms generated by AliForwardMultiplicity and
+ * the file <<i>base</i>><tt>_aods.root</tt> containing the tree
+ * with AliAODEvent objects.
+ *
+ * @param base Base name of files
+ * @param vzMin Minimum collision vertex z position to use
+ * @param vzMax Maximum collision vertex z position to use
+ * @param rebin Rebinning factor
+ *
+ * @return true on success, false otherwise
+ */
+ Bool_t Open(const char* base,
+ Double_t vzMin=-10,
+ Double_t vzMax=10,
+ Int_t rebin=1)
+ {
+ // Set our cuts etc.
+ fVzMin = vzMin;
+ fVzMax = vzMax;
+ if (fVzMax < fVzMin && fVzMin < 0) fVzMax = -fVzMin;
+ fRebin = rebin;
+
+ // Clear previously created data objects
+ fByVtx1D.Delete();
+ if (fTotal1D) { delete fTotal1D; fTotal1D = 0; }
+ if (fTotal2D) { delete fTotal2D; fTotal2D = 0; }
+ if (fTree && fTree->GetCurrentFile()) {
+ fTree->GetCurrentFile()->Close();
+ }
+ fCollect = 0;
+ fTree = 0;
+
+ // Open the AOD file
+ TString fn = TString::Format("%s_aods.root", base);
+ TFile* file = TFile::Open(fn.Data(), "READ");
+ if (!file) {
+ Error("Init", "Couldn't open %s", fn.Data());
+ return kFALSE;
+ }
+
+ // Get the AOD tree
+ fTree = static_cast<TTree*>(file->Get("aodTree"));
+ if (!fTree) {
+ Error("Init", "Couldn't get the tree");
+ return kFALSE;
+ }
+
+ // Set the branch pointer
+ fTree->SetBranchAddress("Forward", &fAOD);
+
+ // Open the histogram file
+ fn = TString::Format("%s_hists.root", base);
+ file = TFile::Open(fn.Data(), "READ");
+ if (!file) {
+ Error("Init", "Couldn't open %s", fn.Data());
+ return kFALSE;
+ }
+
+ // Get the list stored in the file
+ TList* forward =
+ static_cast<TList*>(file->Get("PWG2forwardDnDeta/Forward"));
+ if (!forward) {
+ Error("Init", "Couldn't get forward list");
+ return kFALSE;
+ }
+
+ // Get the list from the collector
+ fCollect =
+ static_cast<TList*>(forward->FindObject("fmdHistCollector"));
+ if (!fCollect) {
+ Error("Init", "Couldn't get collector list");
+ return kFALSE;
+ }
+
+ // Get the event (by vertex bin) histogram
+ fEvents = static_cast<TH1I*>(forward->FindObject("nEventsTrVtx"));
+ if (!fEvents) {
+ Error("Init", "Couldn't get the event histogram");
+ return kFALSE;
+ }
+
+ // Find the min/max bins to use based on the cuts given
+ fBinVzMin = fEvents->FindBin(fVzMin);
+ fBinVzMax = fEvents->FindBin(fVzMax-.0000001);
+ Info("Open", "Selected vertex bins are [%d,%d]", fBinVzMin, fBinVzMax);
+
+ return kTRUE;
+ }
+ //__________________________________________________________________
+ /**
+ * Check if the passed vertex bin number [1,nVtxBins] is within our
+ * cut.
+ *
+ * @param bin Vertex bin [1,nVtxBins]
+ *
+ * @return true if within cut, false otherwise
+ */
+ Bool_t IsInsideVtxCut(Int_t bin) const
+ {
+ return bin >= fBinVzMin && bin <= fBinVzMax;
+ }
+
+ //__________________________________________________________________
+ /**
+ * Make a 1D histogram of @f$ dN_{ch}/d\eta@f$ with the given name and
+ * title
+ *
+ * @param name Name of histogram
+ * @param title Title of histogram
+ * @param a1 Eta axis
+ *
+ * @return Newly allocated 1D histogram object
+ */
+ TH1D* Make1D(const char* name, const char* title, const TAxis* a1)
+ {
+ TH1D* ret = new TH1D(name, title,
+ a1->GetNbins(), a1->GetXmin(), a1->GetXmax());
+ ret->SetXTitle("#eta");
+ ret->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
+ ret->Sumw2();
+ ret->SetMarkerColor(kRed+1);
+ ret->SetLineColor(kRed+1);
+ ret->SetMarkerStyle(24);
+ ret->SetMarkerSize(1);
+ ret->SetStats(0);
+ ret->SetDirectory(0);
+
+ return ret;
+ }
+ //__________________________________________________________________
+ /**
+ * Make a 2D histogram of @f$ d^{2}N_{ch}/d\eta\,d\phi@f$
+ *
+ * @param name Name of histogram
+ * @param title Title of histogram
+ * @param a1 Eta axis
+ * @param a2 Phi axis
+ *
+ * @return
+ */
+ TH2D* Make2D(const char* name, const char* title,
+ const TAxis* a1, const TAxis* a2)
+ {
+ TH2D* ret = new TH2D(name, title,
+ a1->GetNbins(), a1->GetXmin(), a1->GetXmax(),
+ a2->GetNbins(), a2->GetXmin(), a2->GetXmax());
+ ret->SetXTitle("#eta");
+ ret->SetYTitle("#varphi");
+ ret->SetZTitle("#frac{1}{N}#frac{dN^{2}_{ch}}{d#etad#varphi}");
+ ret->Sumw2();
+ ret->SetStats(0);
+ ret->SetDirectory(0);
+
+ return ret;
+ }
+ //__________________________________________________________________
+ /**
+ * Utility function to set-up histograms based on the input
+ * @f$ dd^{2}N_{ch}/d\eta\,d\phi@f$ histogram. This member function
+ * is called on the first event so that we have the proper binning
+ *
+ * @param templ Input histogram
+ *
+ * @return true on succcess
+ */
+ Bool_t FirstEvent(const TH2D& templ)
+ {
+ const TAxis* etaAxis = templ.GetXaxis();
+ const TAxis* phiAxis = templ.GetYaxis();
+
+ // Generate sum histograms.
+ // - fTotal1D will be the sum of projections on the X axis of the input
+ // histograms
+ // - fTotal2D will be the direct sum of the input histograms.
+ // - fFromVtx will be the sum of the per vertex histograms after
+ // processing all events
+ fTotal1D = Make1D("dndeta",
+ "#frac{1}{N}#frac{dN_{ch}}{d#eta} direct sum", etaAxis);
+ fTotal2D = Make2D("d2ndetadphi",
+ "#frac{1}{N}#frac{dN^{2}_{ch}}{d#etad#phi} direct sum",
+ etaAxis, phiAxis);
+ fFromVtx = Make1D("dndeta_test", "#frac{1}{N}#frac{dN_{ch}}{d#eta} "
+ "from VTX", etaAxis);
+ fTotal1D->SetMarkerStyle(25);
+ fTotal1D->SetMarkerSize(1.1);
+ fNorm = new TH1D("norm", "Normalisation",
+ etaAxis->GetNbins(),
+ etaAxis->GetXmin(),
+ etaAxis->GetXmax());
+ fNorm->SetFillColor(kRed);
+ fNorm->SetFillStyle(3001);
+ fNorm->SetXTitle("#eta");
+ fNorm->SetYTitle("Normalisation");
+ fNorm->Sumw2();
+ fNorm->SetDirectory(0);
+ fVtx = new TH1D("vtx", "Events per vertex bin",
+ fEvents->GetXaxis()->GetNbins(),
+ fEvents->GetXaxis()->GetXmin(),
+ fEvents->GetXaxis()->GetXmax());
+ fVtx->SetXTitle("v_{z} [cm]");
+ fVtx->SetYTitle("Events");
+ fVtx->SetDirectory(0);
+ fVtx->SetFillColor(kRed+1);
+ fVtx->SetFillStyle(3001);
+
+
+ // Generate the per-vertex bin histograms. These will be the sum
+ // of each of the per-event input histograms for a given vertex range.
+ for (Int_t i = fBinVzMin; i <= fBinVzMax; i++) {
+ TH1* h1 = Make1D(Form("dndeta_vz%02d", i),
+ Form("#frac{1}{N}#frac{dN_{ch}}{d#eta}, vtxbin=%d", i),
+ etaAxis);
+ fByVtx1D.AddAtAndExpand(h1, i);
+ }
+ return kTRUE;
+ }
+
+
+ //__________________________________________________________________
+ /**
+ * Process the events
+ *
+ *
+ * @return true on success, false otherwise
+ */
+ Bool_t Process()
+ {
+ // Get the number of events in the tree
+ Int_t nEntries = fTree->GetEntries();
+ fNAccepted = 0;
+
+ // Loop over the events in the tree
+ for (Int_t event = 0; event < nEntries; event++) {
+ fTree->GetEntry(event);
+
+ // Get our input histogram
+ const TH2D& hist = fAOD->GetHistogram();
+
+
+ // If fTotal1D or fTotal2D are not made yet, do so (first event)
+ if (!fTotal2D || !fTotal1D) {
+ if (!FirstEvent(hist)) {
+ Error("Process", "Failed to initialize on first event");
+ return kFALSE;
+ }
+ }
+
+ // Check the trigger
+ if (!fAOD->IsTriggerBits(AliAODForwardMult::kInel)) {
+ // Info("Process", "Not an INEL event");
+ continue;
+ }
+
+ // Get the vertex bin - add 1 as we are using histogram bin
+ // numbers in this class
+ Int_t vtxBin = fAOD->GetVtxBin()+1;
+
+ // Check if we're within vertex cut
+ if (!IsInsideVtxCut(vtxBin)) {
+ Info("Process", "In event # %d, %d not within vertex cut "
+ "[%d,%d] (%f,%f)",
+ event, vtxBin, fBinVzMin, fBinVzMax, fVzMin, fVzMax);
+ continue;
+ }
+ fVtx->AddBinContent(vtxBin);
+
+ // Increment our accepted counter
+ fNAccepted++;
+
+ // Get the scale of this vertex range
+ Double_t vtxScale = fEvents->GetBinContent(vtxBin);
+
+ // Get 'by-vertex' histograms
+ TH1* tmp1 = static_cast<TH1D*>(fByVtx1D.At(vtxBin));
+ if (!tmp1) {
+ Warning("Process", "No histogram for vertex bin %d)",vtxBin);
+ continue;
+ }
+
+ Double_t scale = 1. / vtxScale;
+ if (tmp1->InheritsFrom(TProfile::Class()))
+ scale = 1;
+
+ if (!AddContrib(hist, *tmp1, scale)) return kFALSE;
+
+ }
+ for (Int_t iVz = fBinVzMin; iVz <= fBinVzMax; iVz++) {
+ TList* l = static_cast<TList*>(fCollect
+ ->FindObject(Form("vtxbin%02d",iVz-1)));
+ if (!l) {
+ Error("Process", "List vtxbin%02d not found in %s",
+ iVz-1, fCollect->GetName());
+ continue;
+ }
+ TH1F* ve = static_cast<TH1F*>(l->FindObject("etaAcceptance"));
+ if (!ve){
+ Error("Process", "No eta acceptance histogram found in "
+ "vtxbin%02d/etaAcceptance", iVz-1);
+ continue;
+ }
+ ve->Sumw2();
+ fNorm->Add(ve, fVtx->GetBinContent(iVz));
+ }
+ for (Int_t i = 1; i <= fNorm->GetNbinsX(); i++)
+ fNorm->SetBinError(i, 0);
+ // fNorm->Scale(1.);
+
+ Info("Process", "Accepted %6d out of %6d events", fNAccepted, nEntries);
+ return kTRUE;
+ }
+ //__________________________________________________________________
+ Bool_t AddContrib(const TH2D& hist, TH1& vtx1D, Double_t scale)
+ {
+ // Make a projection on the X axis of the input histogram
+ TH1D* proj = hist.ProjectionX("_px", 0, -1, "e");
+
+ // Add to per-vertex histogram
+ vtx1D.Add(proj); // , scale);
+
+#if defined(USE_NORM_HIST)
+ scale = 1;
+#endif
+
+#if defined(USE_WEIGHTED_MEAN)
+ Int_t nBinPhi = hist.GetXaxis()->GetNbins();
+ for (Int_t binEta = 1; binEta <= nBinEta; binEta++) {
+ AddBinContrib(binEta, *proj, *fTotal1D);
+ for (Int_t binPhi = 1; binPhi < nBinPhi; binPhi++)
+ AddBinContrib(hist.GetBin(binEta, binPhi), hist, *fTotal2D);
+ }
+#else
+ // Add to 1D summed histogram
+ fTotal1D->Add(proj); // , scale);
+
+ // Add to 2D summed histogram
+ fTotal2D->Add(&hist); // , scale);
+#endif
+
+ // Remove the projection
+ delete proj;
+
+ return kTRUE;
+ }
+ //__________________________________________________________________
+ Bool_t AddBinContrib(Int_t bin, const TH1& in, TH1& out)
+ {
+ Double_t ic = in.GetBinContent(bin);
+ Double_t ie = in.GetBinError(bin);
+ if (ie <= 0.00001) return kTRUE;
+
+ Double_t iw = 1 / (ie*ie);
+ Double_t oc = out.GetBinContent(bin);
+ Double_t oe = out.GetBinError(bin);
+
+ // Store weighted sum in histogram
+ out.SetBinContent(bin, oc + iw * ic);
+ // Store sum of weights in histogram
+ out.SetBinError(bin, oe + iw);
+
+ return kTRUE;
+ }
+ //__________________________________________________________________
+ Bool_t SetResult(TH1*)
+ {
+#if defined(USE_WEIGHTED_MEAN)
+ Double_t scale = 1.;
+ Int_t nBinEta = fTotal2D->GetXaxis()->GetNbins();
+ Int_t nBinPhi = fTotal1D->GetXaxis()->GetNbins();
+ for (Int_t binEta = 1; binEta <= nBinEta; binEta++) {
+ SetBinResult(binEta, *fTotal1D);
+ for (Int_t binPhi = 1; binPhi < nBinPhi; binPhi++)
+ SetBinResult(fTotal2D->GetBin(binEta, binPhi), *fTotal2D);
+ }
+#elif defined (USE_NORM_HIST)
+ fTotal1D->Divide(fNorm);
+ Int_t nBinEta = fTotal2D->GetXaxis()->GetNbins();
+ Int_t nBinPhi = fTotal1D->GetXaxis()->GetNbins();
+ for (Int_t binEta = 1; binEta <= nBinEta; binEta++) {
+ Double_t n = fNorm->GetBinContent(binEta);
+ for (Int_t binPhi = 1; binPhi < nBinPhi; binPhi++) {
+ if (n <= 0) {
+ fTotal2D->SetBinContent(binEta,binPhi,0);
+ continue;
+ }
+ Double_t c = fTotal2D->GetBinContent(binEta,binPhi);
+ fTotal2D->SetBinContent(binEta,binPhi,c/n);
+ }
+ }
+#else
+#endif
+ // fTotal2D->Scale(1, "width");
+ return kTRUE;
+ }
+
+ //__________________________________________________________________
+ Bool_t SetBinResult(Int_t bin, TH1& in)
+ {
+ Double_t ic = in.GetBinContent(bin);
+ Double_t ie = in.GetBinError(bin);
+ if (ie <= 0.00001) return kTRUE;
+
+ Double_t av = ic / ie;
+ Double_t er = TMath::Sqrt(1/ie);
+
+ // Set bin to be
+ // @f[
+ // \frac{\sum_i w_i c_i}{\sum_i w_i}
+ // @f]
+ // where @f$ w_i = 1/e_i^2@f$, and set the error to be
+ // @f[
+ // \sqrt{\frac{1}{\sum_i w_i}}
+ // @f]
+ in.SetBinContent(bin, av);
+ in.SetBinError(bin, er);
+
+ return kTRUE;
+ }
+ //__________________________________________________________________
+ /**
+ * Called at the end of the job.
+ *
+ * It plots the result of the analysis in tree canvases.
+ * - One shows the per-vertex accumalted histograms and compares them
+ * to the per-vertex histograms made in the analysis.
+ * - Another shows the final @f$ dN_{ch}/d\eta@f$ calculated in
+ * different ways and compare those to the @f$ dN_{ch}/d\eta@f$ made
+ * during the analysis.
+ * - The last canvas shows the @f$ dd^{2}N_{ch}/d\eta\,d\phi@f$ histogram.
+ *
+ * @return true on succes, false otherwise
+ */
+ Bool_t Finish()
+ {
+ TFile* out = TFile::Open("out.root", "RECREATE");
+
+ // Set the style
+ gStyle->SetOptTitle(0);
+ gStyle->SetPadColor(0);
+ gStyle->SetPadBorderSize(0);
+ gStyle->SetPadBorderMode(0);
+ gStyle->SetPadRightMargin(0.05);
+ gStyle->SetPadTopMargin(0.05);
+ gStyle->SetPalette(1);
+
+ // Get number of bins
+ Int_t nBin = fBinVzMax-fBinVzMin+1;
+
+ // Make first canvas
+ TCanvas* cVtx = new TCanvas("cVtx", "By Vertex", 1000, 700);
+ cVtx->SetBorderSize(0);
+ cVtx->SetBorderMode(0);
+ cVtx->Divide((nBin+.5)/2, 2, 0.0005, 0.0005);
+
+ // Loop over vertex histograms
+ Int_t nHists = 0;
+ for (Int_t i = fBinVzMin; i <= fBinVzMax; i++) {
+ TDirectory* vtxDir = out->mkdir(Form("vtx%02d", i));
+ vtxDir->cd();
+
+ TH1D* vh1 = static_cast<TH1D*>(fByVtx1D.At(i));
+ if (!vh1) {
+ Error("Finish", "No histogram at %d", i);
+ continue;
+ }
+
+ fFromVtx->Add(vh1);
+
+ // Scale and add to output
+ if (fVtx->GetBinContent(i) > 0)
+ vh1->Scale(1. / fVtx->GetBinContent(i), "width");
+
+ // Write to output file
+ vh1->Write();
+
+ // Divide the pad
+ TVirtualPad* pad = cVtx->cd(i-fBinVzMin+1);
+ pad->Divide(1,2,0.005,0);
+
+ // Draw the per-vertex histogram
+ pad->cd(1);
+ vh1->Draw();
+
+ // Get the same histogram from the analysis and draw it
+ TProfile* p =
+ static_cast<TProfile*>(fCollect
+ ->FindObject(Form("dndeta_vtx%02d",i-1)));
+ p->SetMarkerColor(kBlue+1);
+ p->SetLineColor(kBlue+1);
+ p->SetMarkerSize(0.8);
+ p->SetMarkerStyle(20);
+ p->Draw("same");
+ p->Write();
+
+ // Make the ratio of the two histograms and draw it
+ pad->cd(2);
+ TH1* ratio = static_cast<TH1*>(vh1->Clone(Form("ratio_vtx%02d", i-1)));
+ // ratio->Scale(1./fEvents->GetBinContent(i));
+ ratio->Divide(p);
+ ratio->SetMarkerSize(.8);
+ ratio->SetLineColor(kGray);
+ ratio->Draw("P");
+ }
+ cVtx->cd();
+ out->cd();
+
+ // Get the number of events selected
+ // Int_t nEvSelected = fEvents->Integral(fBinVzMin, fBinVzMax);
+
+ // The second canvas
+ TCanvas* cTotal1D = new TCanvas("total", "Result", 800, 800);
+ cTotal1D->SetFillColor(0);
+ cTotal1D->SetBorderMode(0);
+ cTotal1D->SetBorderSize(0);
+ cTotal1D->cd();
+
+ // Make our main pad
+ TPad* p1 = new TPad("p1", "p1", 0, 0.3, 1.0, 1.0, 0, 0);
+ p1->Draw();
+ p1->cd();
+
+ // Make a stack to 'auto-scale' when plotting more than 1 histogram.
+ THStack* stack = new THStack("results", "Results for #frac{dN_{ch}{d#eta}");
+
+ // Scale the histogram summed over the vertex bins. This must be
+ // divided by the number of bins we have summed over. If we do
+ // rebinning, we should scale it again.
+ fFromVtx->Divide(fNorm);
+ // fFromVtx->Scale(1.);
+ fFromVtx->Scale(1, "width");
+
+
+ if (fRebin > 1) { fFromVtx->Rebin(fRebin); fFromVtx->Scale(1. / fRebin); }
+ stack->Add(fFromVtx, "");
+
+ // Generate the projection from the direct sum of the input histograms,
+ // and scale it to the number of vertex bins and the bin width
+ TH1* proj = fTotal2D->ProjectionX("dndeta_proj", 0, -1, "e");
+ proj->SetLineColor(kGreen+1);
+ proj->SetMarkerColor(kGreen+1);
+ proj->SetMarkerStyle(24);
+ proj->SetMarkerSize(1.5);
+ proj->SetTitle(Form("%s directly", proj->GetTitle()));
+ proj->Divide(fNorm);
+ proj->Scale(1., "width");
+ stack->Add(proj, "");
+
+ // Scale our direct sum of the projects of the input histograms to
+ // the number of vertex bins and the bin width. If we do rebinning,
+ // we must scale it one more time.
+ // fTotal1D->Scale(1. / fNAccepted, "width");
+ fTotal1D->Divide(fNorm);
+ fTotal1D->Scale(1., "width");
+ if (fRebin > 1) { fTotal1D->Rebin(fRebin); fTotal1D->Scale(1. / fRebin); }
+ stack->Add(fTotal1D);
+
+ // Get the result from the analysis and plit that too (after modifying
+ // the style and possible rebinning)
+ TProfile* dtotal = static_cast<TProfile*>(fCollect->FindObject("dndeta"));
+ if (!dtotal) {
+ Error("Finish", "Couldn't get the event histogram");
+ return kFALSE;
+ }
+ dtotal->SetTitle(Form("%s directly", dtotal->GetTitle()));
+ dtotal->SetLineColor(kBlue+1);
+ dtotal->SetMarkerColor(kBlue+1);
+ dtotal->SetMarkerStyle(20);
+ dtotal->SetMarkerSize(0.8);
+ if (fRebin > 1) { dtotal->Rebin(fRebin); }
+ stack->Add(dtotal, "");
+
+ // Draw stack
+ stack->Draw("nostack e1");
+ stack->Write();
+
+ // Make a legend
+ TLegend* l = p1->BuildLegend(0.31, 0.15, 0.5, 0.6);
+ l->SetFillColor(0);
+ l->SetBorderSize(0);
+ l->SetTextSize(0.04); // l->GetTextSize()/3);
+ cTotal1D->cd();
+
+ // Generate our second pad
+ TPad* p2 = new TPad("p2", "p2", 0, 0.0, 1.0, 0.3, 0, 0);
+ p2->Draw();
+ p2->cd();
+
+ // Create a new stack
+ stack = new THStack("ratios", "Ratios to old method");
+ // Calculate the ratio of our summed over vertex bins to the
+ // result from the analysis and draw it
+ TH1* ratioVtx = static_cast<TH1*>(fFromVtx->Clone("ratioVtx"));
+ ratioVtx->SetDirectory(0);
+ ratioVtx->Divide(dtotal);
+ stack->Add(ratioVtx, "");
+
+ // Calculate the ratio of our direct sum of input histograms
+ // result from the analysis and draw it
+ TH1* ratioSum = static_cast<TH1*>(fTotal1D->Clone("ratioSum"));
+ ratioSum->SetDirectory(0);
+ ratioSum->Divide(dtotal);
+ stack->Add(ratioSum, "");
+
+ // Calculate the ratio of our direct sum of input histograms
+ // result from the analysis and draw it
+ TH1* ratioProj = static_cast<TH1*>(proj->Clone("ratioProj"));
+ ratioProj->SetDirectory(0);
+ ratioProj->Divide(dtotal);
+ stack->Add(ratioProj, "");
+
+ stack->Draw("nostack e1");
+ stack->Write();
+
+ // update
+ cTotal1D->cd();
+
+ // Generate the last canvas and show the summed input histogram
+ TCanvas* other = new TCanvas("other", "Other", 800, 600);
+ other->SetFillColor(0);
+ other->SetBorderMode(0);
+ other->SetBorderSize(0);
+ other->Divide(2,1);
+
+ other->cd(1);
+ fNorm->Draw("HIST");
+ fNorm->Write();
+
+ other->cd(2);
+ fEvents->Draw();
+ fVtx->Draw("same");
+ fEvents->Write();
+ // fTotal2D->Draw("lego2 e");
+
+ return kTRUE;
+ }
+ //__________________________________________________________________
+ /**
+ * Run the post-processing.
+ *
+ * This will open the files <<i>base</i>><tt>_hists.root</tt>
+ * containing the histograms generated by AliForwardMultiplicity and
+ * the file <<i>base</i>><tt>_aods.root</tt> containing the
+ * tree with AliAODEvent objects.
+ *
+ * Then it will loop over the events, accepting only INEL events
+ * that have a primary collision point along z within @a vzMin and
+ * @a vzMax.
+ *
+ * After the processing, the result will be shown in 3 canvases,
+ * possibly rebinning the result by the factor @a rebin.
+ *
+ * @param base Base name of files
+ * @param vzMin Minimum collision vertex z position to use
+ * @param vzMax Maximum collision vertex z position to use
+ * @param rebin Rebinning factor
+ *
+ * @return true on success, false otherwise
+ */
+ Bool_t Run(const char* base,
+ Double_t vzMin=-10,
+ Double_t vzMax=10,
+ Int_t rebin=1)
+ {
+ if (!Open(base,vzMin,vzMax,rebin)) return kFALSE;
+ if (!Process()) return kFALSE;
+ if (!Finish()) return kFALSE;
+
+ return kTRUE;
+ }
+protected:
+ Double_t fVzMin; // Minimum vertex
+ Double_t fVzMax; // Maximum vertex
+ Int_t fBinVzMin; // Corresponding bin to min vertex
+ Int_t fBinVzMax; // Corresponding bin to max vertex
+ Int_t fRebin; // Rebin factor
+ TTree* fTree; // Event tree
+ AliAODForwardMult* fAOD; // Our event-by-event data
+ TH1I* fEvents; // histogram of event counts per vertex
+ TList* fCollect; // List of analysis histograms
+ TObjArray fByVtx1D; // List of per-vertex 1D histograms
+ TH1D* fTotal1D; // Direct sum of input histograms
+ TH2D* fTotal2D; // Direct sum of input histograms
+ TH1D* fFromVtx; // Sum of per-vertex histograms
+ TH1D* fNorm; // Histogram of # events with data per bin
+ TH1D* fVtx;
+ Int_t fNAccepted; // # of accepted events
+};
+
+//
+// EOF
+//
--- /dev/null
+#include "DrawRes1D.C"
+
+/**
+ * @defgroup pwg2_forward_analysis_scripts PWG2 Forward analysis - scripts
+ *
+ * @ingroup pwg2_forward_analysis
+ */
+/**
+ * Example macro to loop over the event-by-event 2D histogram of
+ * @f[
+ * \frac{d^{2}N_{ch}}{d\eta\,d\phi}
+ * @f]
+ * stored in an AOD.
+ *
+ * The class needs the files <<i>base</i>><tt>_hists.root</tt>
+ * containing the histograms generated by AliForwardMultiplicity and
+ * the file <<i>base</i>><tt>_aods.root</tt> containing the tree
+ * with AliAODEvent objects where AliAODForwardMult objects have been
+ * added to in the branch <tt>Forward</tt>
+ *
+ * @ingroup pwg2_forward_analysis_scripts
+ */
+class DrawResVtx : public DrawRes1D
+{
+public:
+ /**
+ * Constructor
+ *
+ * @param special If true, add to the list of 'specials'
+ */
+ DrawResVtx()
+ : fByVtx1D(0)
+ {}
+ //__________________________________________________________________
+ /**
+ * Open the files <<i>base</i>><tt>_hists.root</tt>
+ * containing the histograms generated by AliForwardMultiplicity and
+ * the file <<i>base</i>><tt>_aods.root</tt> containing the tree
+ * with AliAODEvent objects.
+ *
+ * @param base Base name of files
+ * @param vzMin Minimum collision vertex z position to use
+ * @param vzMax Maximum collision vertex z position to use
+ * @param rebin Rebinning factor
+ *
+ * @return true on success, false otherwise
+ */
+ Bool_t Open(const char* base,
+ Double_t vzMin=-10,
+ Double_t vzMax=10)
+ {
+ // Clear previously created data objects
+ fByVtx1D.Delete();
+ return DrawRes1D::Open(base, vzMin, vzMax);
+ }
+
+ //__________________________________________________________________
+ /**
+ * Make a 1D histogram of @f$ dN_{ch}/d\eta@f$ with the given name and
+ * title
+ *
+ * @param name Name of histogram
+ * @param title Title of histogram
+ * @param a1 Eta axis
+ *
+ * @return Newly allocated 1D histogram object
+ */
+ TH1D* Make1D(const char* name, const char* title, const TAxis* a1)
+ {
+ TH1D* ret = new TH1D(name, title,
+ a1->GetNbins(), a1->GetXmin(), a1->GetXmax());
+ ret->SetXTitle("#eta");
+ ret->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
+ ret->Sumw2();
+ ret->SetMarkerColor(kRed+1);
+ ret->SetLineColor(kRed+1);
+ ret->SetMarkerStyle(20);
+ ret->SetMarkerSize(1);
+ ret->SetStats(0);
+ ret->SetDirectory(0);
+
+ return ret;
+ }
+ //__________________________________________________________________
+ /**
+ * Utility function to set-up histograms based on the input
+ * @f$ dd^{2}N_{ch}/d\eta\,d\phi@f$ histogram. This member function
+ * is called on the first event so that we have the proper binning
+ *
+ * @param templ Input histogram
+ *
+ * @return true on succcess
+ */
+ Bool_t FirstEvent(const TH2D& templ)
+ {
+ if (!DrawRes1D::FirstEvent(templ)) return kFALSE;
+
+ const TAxis* etaAxis = templ.GetXaxis();
+
+ // Generate the per-vertex bin histograms. These will be the sum
+ // of each of the per-event input histograms for a given vertex range.
+ for (Int_t i = fBinVzMin; i <= fBinVzMax; i++) {
+ TH1* h1 = Make1D(Form("dndeta_vz%02d", i),
+ Form("#frac{1}{N}#frac{dN_{ch}}{d#eta}, vtxbin=%d", i),
+ etaAxis);
+ fByVtx1D.AddAtAndExpand(h1, i);
+ }
+ return kTRUE;
+ }
+
+
+ //__________________________________________________________________
+ /**
+ * Process the events
+ *
+ *
+ * @return true on success, false otherwise
+ */
+ void AddContrib(const TH2D& hist, Int_t vtxBin)
+ {
+ // Get 'by-vertex' histograms
+ TH1* tmp1 = static_cast<TH1D*>(fByVtx1D.At(vtxBin));
+ if (!tmp1) {
+ Warning("Process", "No histogram for vertex bin %d)",vtxBin);
+ return;
+ }
+ // Make a projection on the X axis of the input histogram
+ TH1D* proj = hist.ProjectionX("_px", 0, -1, "e");
+
+ // Add to per-vertex histogram
+ tmp1->Add(proj);
+
+ delete proj;
+ }
+ //__________________________________________________________________
+ /**
+ * Called at the end of the job.
+ *
+ * It plots the result of the analysis in tree canvases.
+ * - One shows the per-vertex accumalted histograms and compares them
+ * to the per-vertex histograms made in the analysis.
+ * - Another shows the final @f$ dN_{ch}/d\eta@f$ calculated in
+ * different ways and compare those to the @f$ dN_{ch}/d\eta@f$ made
+ * during the analysis.
+ * - The last canvas shows the @f$ dd^{2}N_{ch}/d\eta\,d\phi@f$ histogram.
+ *
+ * @return true on succes, false otherwise
+ */
+ TH1D* GetResult()
+ {
+ // Get number of bins
+ Int_t nBin = fBinVzMax-fBinVzMin+1;
+
+ // Make first canvas
+ TCanvas* cVtx = new TCanvas("cVtx", "By Vertex", 1000, 700);
+ cVtx->SetBorderSize(0);
+ cVtx->SetBorderMode(0);
+ cVtx->Divide((nBin+.5)/2, 2, 0.0005, 0.0005);
+
+ // Loop over vertex histograms
+ TDirectory* savdir = gDirectory;
+ for (Int_t i = fBinVzMin; i <= fBinVzMax; i++) {
+ TDirectory* vtxDir = savdir->mkdir(Form("vtx%02d", i));
+ vtxDir->cd();
+
+ TH1D* vh1 = static_cast<TH1D*>(fByVtx1D.At(i));
+ if (!vh1) {
+ Error("Finish", "No histogram at %d", i);
+ continue;
+ }
+
+ fTotal1D->Add(vh1);
+
+ // Scale and add to output
+ if (fVtx->GetBinContent(i) > 0)
+ vh1->Scale(1. / fVtx->GetBinContent(i), "width");
+
+ // Write to output file
+ vh1->Write();
+
+ // Divide the pad
+ TVirtualPad* pad = cVtx->cd(i-fBinVzMin+1);
+ pad->Divide(1,2,0.005,0);
+
+ // Draw the per-vertex histogram
+ pad->cd(1);
+ vh1->DrawClone();
+
+ // Get the same histogram from the analysis and draw it
+ TProfile* p =
+ static_cast<TProfile*>(fCollect
+ ->FindObject(Form("dndeta_vtx%02d",i-1)));
+ p->SetMarkerColor(kBlue+1);
+ p->SetLineColor(kBlue+1);
+ p->SetMarkerSize(0.8);
+ p->SetMarkerStyle(20);
+ p->DrawClone("same");
+ p->Write();
+
+ // Make the ratio of the two histograms and draw it
+ pad->cd(2);
+ TH1* ratio = static_cast<TH1*>(vh1->Clone(Form("ratio_vtx%02d", i-1)));
+ // ratio->Scale(1./fEvents->GetBinContent(i));
+ ratio->Divide(p);
+ ratio->SetMarkerSize(.8);
+ ratio->SetLineColor(kGray);
+ ratio->DrawClone("P");
+ }
+ cVtx->cd();
+ savdir->cd();
+
+ // Scale the histogram summed over the vertex bins. This must be
+ // divided by the number of bins we have summed over. If we do
+ // rebinning, we should scale it again.
+ fTotal1D->Divide(fNorm);
+ fTotal1D->Scale(1, "width");
+ return fTotal1D;
+ }
+ /**
+ * Browse this object
+ *
+ * @param b Browser to use
+ */
+ void Browse(TBrowser* b)
+ {
+ b->Add(&fByVtx1D);
+ DrawRes1D::Browse(b);
+ }
+ /**
+ * Create a new object of this class, and add it to the list of
+ * specials, and create a browse and browse to this object
+ *
+ * @return Newly created object
+ */
+ static DrawResVtx* Create()
+ {
+ DrawResVtx* dr = new DrawResVtx;
+ gROOT->GetListOfSpecials()->Add(dr);
+ TBrowser* b = new TBrowser("b");
+ b->BrowseObject(gROOT->GetListOfSpecials());
+ return dr;
+ }
+
+protected:
+ TObjArray fByVtx1D; // List of per-vertex 1D histograms
+
+ ClassDef(DrawResVtx,0)
+};
+
+//
+// EOF
+//