--- /dev/null
+#include "AliFMDMCCorrections.h"
+#include <AliESDFMD.h>
+#include <TAxis.h>
+#include <TList.h>
+#include <TMath.h>
+#include "AliForwardCorrectionManager.h"
+// #include "AliFMDAnaParameters.h"
+#include "AliLog.h"
+#include <TH2D.h>
+#include <TROOT.h>
+#include <TProfile2D.h>
+#include <iostream>
+#include <iomanip>
+
+ClassImp(AliFMDMCCorrections)
+#if 0
+; // For Emacs
+#endif
+
+
+//____________________________________________________________________
+AliFMDMCCorrections::~AliFMDMCCorrections()
+{
+ if (fComps) fComps->Clear();
+ if (fFMD1i) delete fFMD1i;
+ if (fFMD2i) delete fFMD2i;
+ if (fFMD2o) delete fFMD2o;
+ if (fFMD3i) delete fFMD3i;
+ if (fFMD3o) delete fFMD3o;
+}
+
+//____________________________________________________________________
+AliFMDMCCorrections&
+AliFMDMCCorrections::operator=(const AliFMDMCCorrections& o)
+{
+ AliFMDCorrections::operator=(o);
+
+ return *this;
+}
+
+//____________________________________________________________________
+Bool_t
+AliFMDMCCorrections::CorrectMC(AliForwardUtil::Histos& hists,
+ UShort_t vtxbin)
+{
+ AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
+
+ UShort_t uvb = vtxbin;
+ 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* bg = fcm.GetSecondaryMap()->GetCorrection(d,r,uvb);
+ TH2D* ef = fcm.GetVertexBias()->GetCorrection(r, uvb);
+ if (!bg) {
+ AliWarning(Form("No secondary correction for FMDM%d%c in vertex bin %d",
+ d, r, uvb));
+ continue;
+ }
+ if (!ef) {
+ AliWarning(Form("No event vertex bias correction in vertex bin %d",
+ uvb));
+ continue;
+ }
+
+ // Divide by primary/total ratio
+ h->Divide(bg);
+
+ // Divide by the event selection efficiency
+ h->Divide(ef);
+ }
+ }
+
+ return kTRUE;
+}
+
+//____________________________________________________________________
+void
+AliFMDMCCorrections::Init(const TAxis& eAxis)
+{
+ fFMD1i = Make(1,'I',eAxis);
+ fFMD2i = Make(2,'I',eAxis);
+ fFMD2o = Make(2,'O',eAxis);
+ fFMD3i = Make(3,'I',eAxis);
+ fFMD3o = Make(3,'O',eAxis);
+
+ fComps->Add(fFMD1i);
+ fComps->Add(fFMD2i);
+ fComps->Add(fFMD2o);
+ fComps->Add(fFMD3i);
+ fComps->Add(fFMD3o);
+}
+
+//____________________________________________________________________
+TProfile2D*
+AliFMDMCCorrections::Make(UShort_t d, Char_t r,
+ const TAxis& axis) const
+{
+ TProfile2D* ret = new TProfile2D(Form("FMD%d%c_esd_vs_mc", d, r),
+ Form("ESD/MC signal for FMD%d%c", d, r),
+ axis.GetNbins(),
+ axis.GetXmin(),
+ axis.GetXmax(),
+ (r == 'I' || r == 'i') ? 20 : 40,
+ 0, 2*TMath::Pi());
+ ret->GetXaxis()->SetTitle("#eta");
+ ret->GetYaxis()->SetTitle("#varphi [degrees]");
+ ret->GetZaxis()->SetTitle("#LT primary density ESD/MC#GT");
+ ret->SetDirectory(0);
+ return ret;
+}
+//____________________________________________________________________
+void
+AliFMDMCCorrections::Fill(UShort_t d, Char_t r, TH2* esd, TH2* mc)
+{
+ if (!esd || !mc) return;
+ TProfile2D* p = 0;
+ switch (d) {
+ case 1: p = fFMD1i; break;
+ case 2: p = (r == 'I' || r == 'i' ? fFMD2i : fFMD2o); break;
+ case 3: p = (r == 'I' || r == 'i' ? fFMD3i : fFMD3o); break;
+ }
+ if (!p) return;
+
+ for (Int_t iEta = 1; iEta <= esd->GetNbinsX(); iEta++) {
+ Double_t eta = esd->GetXaxis()->GetBinCenter(iEta);
+ for (Int_t iPhi = 1; iPhi <= esd->GetNbinsY(); iPhi++) {
+ Double_t phi = esd->GetYaxis()->GetBinCenter(iPhi);
+ Double_t mEsd = esd->GetBinContent(iEta,iPhi);
+ Double_t mMc = mc->GetBinContent(iEta,iPhi);
+
+ p->Fill(eta, phi, (mMc > 0 ? mEsd / mMc : 0));
+ }
+ }
+}
+
+//____________________________________________________________________
+Bool_t
+AliFMDMCCorrections::CompareResults(AliForwardUtil::Histos& esd,
+ AliForwardUtil::Histos& mc)
+{
+ Fill(1, 'I', esd.Get(1,'I'), mc.Get(1,'I'));
+ Fill(2, 'I', esd.Get(2,'I'), mc.Get(2,'I'));
+ Fill(2, 'O', esd.Get(2,'O'), mc.Get(2,'O'));
+ Fill(3, 'I', esd.Get(3,'I'), mc.Get(3,'I'));
+ Fill(3, 'O', esd.Get(3,'O'), mc.Get(3,'O'));
+
+ return kTRUE;
+}
+
+//____________________________________________________________________
+void
+AliFMDMCCorrections::DefineOutput(TList* dir)
+{
+ AliFMDCorrections::DefineOutput(dir);
+ TList* d = static_cast<TList*>(dir->FindObject(GetName()));
+
+ fComps = new TList;
+ fComps->SetName("esd_mc_comparison");
+ d->Add(fComps);
+}
+
+//____________________________________________________________________
+//
+// EOF
+//
+
+
+
--- /dev/null
+#ifndef ALIFMDMCCORRECTIONS_H
+#define ALIFMDMCCORRECTIONS_H
+#include "AliFMDCorrections.h"
+#include <TList.h>
+class TProfile2D;
+class TH2;
+
+/**
+ * @defgroup pwg2_forward_mc Monte-carlo code
+ *
+ * @ingroup pwg2_forward
+ */
+/**
+ * This class calculates the exclusive charged particle density
+ * in each for the 5 FMD rings.
+ *
+ * @par Input:
+ * - 5 RingHistos objects - each with a number of vertex dependent
+ * 2D histograms of the inclusive charge particle density
+ *
+ * @par Output:
+ * - 5 RingHistos objects - each with a number of vertex dependent
+ * 2D histograms of the exclusive charge particle density
+ *
+ * @par Corrections used:
+ * - AliFMDCorrSecondaryMap;
+ * - AliFMDCorrVertexBias
+ * - AliFMDCorrMergingEfficiency
+ *
+ * @ingroup pwg2_forward_algo
+ * @ingroup pwg2_forward_mc
+ */
+class AliFMDMCCorrections : public AliFMDCorrections
+{
+public:
+ /**
+ * Constructor
+ */
+ AliFMDMCCorrections()
+ : AliFMDCorrections(),
+ fFMD1i(0),
+ fFMD2i(0),
+ fFMD2o(0),
+ fFMD3i(0),
+ fFMD3o(0),
+ fComps(0)
+ {}
+ /**
+ * Constructor
+ *
+ * @param name Name of object
+ */
+ AliFMDMCCorrections(const char* name)
+ : AliFMDCorrections(name),
+ fFMD1i(0),
+ fFMD2i(0),
+ fFMD2o(0),
+ fFMD3i(0),
+ fFMD3o(0),
+ fComps(0)
+ {}
+ /**
+ * Copy constructor
+ *
+ * @param o Object to copy from
+ */
+ AliFMDMCCorrections(const AliFMDMCCorrections& o)
+ : AliFMDCorrections(o),
+ fFMD1i(o.fFMD1i),
+ fFMD2i(o.fFMD2i),
+ fFMD2o(o.fFMD2o),
+ fFMD3i(o.fFMD3i),
+ fFMD3o(o.fFMD3o),
+ fComps(0)
+ {}
+ /**
+ * Destructor
+ */
+ virtual ~AliFMDMCCorrections();
+ /**
+ * Assignement operator
+ *
+ * @param o Object to assign from
+ *
+ * @return Reference to this object
+ */
+ AliFMDMCCorrections& operator=(const AliFMDMCCorrections&);
+ /**
+ * Initialize this object
+ *
+ * @param etaAxis Eta axis to use
+ */
+ void Init(const TAxis& etaAxis);
+ /**
+ * Do the calculations
+ *
+ * @param hists Cache of histograms
+ * @param vtxBin Vertex bin
+ *
+ * @return true on successs
+ */
+ virtual Bool_t CorrectMC(AliForwardUtil::Histos& hists, UShort_t vtxBin);
+ /**
+ * Compare the result of analysing the ESD for
+ * the inclusive charged particle density to analysing
+ * MC truth
+ *
+ * @param esd
+ * @param mc
+ *
+ * @return
+ */
+ virtual Bool_t CompareResults(AliForwardUtil::Histos& esd,
+ AliForwardUtil::Histos& mc);
+ /**
+ * Output diagnostic histograms to directory
+ *
+ * @param dir List to write in
+ */
+ void DefineOutput(TList* dir);
+protected:
+ /**
+ * MAke comparison profiles
+ *
+ * @param d Detector
+ * @param r Ring
+ * @param axis Eta axis
+ *
+ * @return Newly allocated profile object
+ */
+ TProfile2D* Make(UShort_t d, Char_t r, const TAxis& axis) const;
+ /**
+ * Fill comparison profiles
+ *
+ * @param d Detector
+ * @param r Ring
+ * @param esd ESD histogram
+ * @param mc MC histogram
+ */
+ void Fill(UShort_t d, Char_t r, TH2* esd, TH2* mc);
+
+ TProfile2D* fFMD1i; // Comparison
+ TProfile2D* fFMD2i; // Comparison
+ TProfile2D* fFMD2o; // Comparison
+ TProfile2D* fFMD3i; // Comparison
+ TProfile2D* fFMD3o; // Comparison
+ TList* fComps; // List of comparisons
+
+ ClassDef(AliFMDMCCorrections,1); // Calculate Nch density
+};
+
+#endif
+// Local Variables:
+// mode: C++
+// End:
+