#ifdef BUILD
#include "AliAnalysisTaskSE.h"
#include "AliAnalysisManager.h"
#include "AliAnalysisDataContainer.h"
#include "AliMCEvent.h"
#include "AliESDEvent.h"
#include "AliStack.h"
#include "AliMultiplicity.h"
#include "AliFMDMCEventInspector.h"
#include "AliAODForwardMult.h"
#include "AliLog.h"
#include <TH1I.h>
#include <TH2D.h>
#include <TAxis.h>
#include <TList.h>
#include <TObjArray.h>
#include <TParameter.h>
#include <TStopwatch.h>
#include <TROOT.h>
#include <THStack.h>
#include <TStyle.h>

//====================================================================
/**
 * Task to evaluate trigger bias in pp 
 * 
 */
class EvaluateTrigger : public AliAnalysisTaskSE
{
public:
  enum { 
    kNone = 0x0, 
    kESD  = 0x1, 
    kMC   = 0x2
  };
  enum { 
    kMaxN = 9
  };
  /** 
   * Constructor
   */
  EvaluateTrigger() 
    : AliAnalysisTaskSE(),
      fInel(),
      fInelGt0(),
      fNSD(),
      fNClusterGt0(),
      fInspector(), 
      fFirstEvent(true),
      fData(0), 
      fTriggers(0), 
      fTrackletRequirement(kESD),
      fVertexRequirement(kESD), 
      fVertexAxis(0, 0, 0), 
      fVertexESD(0),
      fVertexMC(0), 
      fM(0)
  {}
  /** 
   * Constructor 
   */
  EvaluateTrigger(const char* /*name*/) 
    : AliAnalysisTaskSE("evaluateTrigger"),
      fInel(AliAODForwardMult::kInel),
      fInelGt0(AliAODForwardMult::kInelGt0),
      fNSD(AliAODForwardMult::kNSD),
      fNClusterGt0(AliAODForwardMult::kNClusterGt0),
      fInspector("eventInspector"), 
      fFirstEvent(true), 
      fData(0), 
      fTriggers(0),
      fTrackletRequirement(kESD),
      fVertexRequirement(kESD), 
      fVertexAxis(10, -10, 10), 
      fVertexESD(0),
      fVertexMC(0), 
      fM(0)
  {
    DefineOutput(1, TList::Class());
    DefineOutput(2, TList::Class());
  }
  void SetVertexRequirement(UShort_t m) { fVertexRequirement = m; }
  void SetTrackletRequirement(UShort_t m) { fTrackletRequirement = m; }
  void SetVertexAxis(Int_t nBins, Double_t low, Double_t high) 
  {
    fVertexAxis.Set(nBins, low, high);
  }
  /** 
   * Intialize 
   */
  void Init() {}
  /** 
   * Create output objects 
   */
  void UserCreateOutputObjects()
  {
    fList = new TList;
    fList->SetOwner();
    fList->SetName("triggerSums");

    // Double_t mb[] = { 0, 1, 2, 3, 4, 5, 6, 8, 9, 10, 11 };
    // Int_t    nM   = 10;
    TAxis eAxis(200, -4, 6);
    TAxis pAxis(40, 0, 2*TMath::Pi());

    fData = new TH2D("data", "Cache", 
		     eAxis.GetNbins(), eAxis.GetXmin(), eAxis.GetXmax(), 
		     pAxis.GetNbins(), pAxis.GetXmin(), pAxis.GetXmax());
    fData->SetDirectory(0);
    fData->SetXTitle("#eta");
    fData->SetYTitle("#varphi [radians]");
    fData->SetZTitle("N_{ch}(#eta,#varphi)");
    fData->Sumw2();
    
    fM = new TH1D("m", "Distribution of N_{ch}|_{|#eta|<1}", kMaxN+1,0,kMaxN+1);
    fM->SetXTitle("N_{ch}|_{|#eta|<1}");
    fM->SetYTitle("Events");
    fM->SetFillColor(kRed+1);
    fM->SetFillStyle(3001);
    fM->SetDirectory(0);
    fList->Add(fM);

    for (Int_t i = 0; i <= kMaxN; i++) { 
      TString lbl;
      if (i == 0)          lbl = "all";
      else if (i == kMaxN) lbl = Form("%d+",i-1);
      else                 lbl = Form("<%d",i);
      fM->GetXaxis()->SetBinLabel(i+1, lbl);
    }

    fTriggers = new TH1I("triggers", "Triggers", 8, -.5, 7.5);
    fTriggers->SetDirectory(0);
    fTriggers->GetXaxis()->SetBinLabel(1, "INEL (MC)");
    fTriggers->GetXaxis()->SetBinLabel(2, "INEL (ESD)");
    fTriggers->GetXaxis()->SetBinLabel(3, "INEL & N_{cluster}>0 (MC)");
    fTriggers->GetXaxis()->SetBinLabel(4, "INEL & N_{cluster}>0 (ESD)");
    fTriggers->GetXaxis()->SetBinLabel(5, "INEL>0 (MC)");
    fTriggers->GetXaxis()->SetBinLabel(6, "INEL>0 (ESD)");
    fTriggers->GetXaxis()->SetBinLabel(7, "NSD (MC)");
    fTriggers->GetXaxis()->SetBinLabel(8, "NSD (ESD)");
    fTriggers->SetFillColor(kYellow+1);
    fTriggers->SetFillStyle(3001);
    fList->Add(fTriggers);

    fVertexESD = new TH1D("vertexESD", "ESD vertex distribution", 
			  fVertexAxis.GetNbins(), 
			  fVertexAxis.GetXmin(), 
			  fVertexAxis.GetXmax());
    fVertexESD->SetDirectory(0);
    fVertexESD->SetFillColor(kRed+1);
    fVertexESD->SetFillStyle(3001);
    fVertexESD->SetXTitle("v_{z} [cm]");
    fVertexESD->SetYTitle("P(v_{z})");
    fList->Add(fVertexESD);

    fVertexMC = new TH1D("vertexMC", "MC vertex distribution", 
			  fVertexAxis.GetNbins(), 
			  fVertexAxis.GetXmin(), 
			  fVertexAxis.GetXmax());
    fVertexMC->SetDirectory(0);
    fVertexMC->SetFillColor(kBlue+1);
    fVertexMC->SetFillStyle(3001);
    fVertexMC->SetXTitle("v_{z} [cm]");
    fVertexMC->SetYTitle("P(v_{z})");
    fList->Add(fVertexMC);

    fInel.CreateObjects(fList, fM, fData);
    fInelGt0.CreateObjects(fList, fM, fData);
    fNSD.CreateObjects(fList, fM, fData);
    fNClusterGt0.CreateObjects(fList, fM, fData);


    fInspector.DefineOutput(fList);
    fInspector.Init(fVertexAxis);

    PostData(1, fList);
  }
  /** 
   * Event processing 
   */
  void UserExec(Option_t*) 
  {
    // Get the input data - MC event
    AliMCEvent*  mcEvent = MCEvent();
    if (!mcEvent) { 
      AliWarning("No MC event found");
      return;
    }
    
    // Get the input data - ESD event
    AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
    if (!esd) { 
      AliWarning("No ESD event found for input event");
      return;
    }

    if (fFirstEvent && esd->GetESDRun()) {
      fInspector.ReadRunDetails(esd);

      AliInfo(Form("Initializing with parameters from the ESD:\n"
		   "         AliESDEvent::GetBeamEnergy()   ->%f\n"
		   "         AliESDEvent::GetBeamType()     ->%s\n"
		   "         AliESDEvent::GetCurrentL3()    ->%f\n"
		   "         AliESDEvent::GetMagneticField()->%f\n"
		   "         AliESDEvent::GetRunNumber()    ->%d\n",
		   esd->GetBeamEnergy(),
		   esd->GetBeamType(),
		   esd->GetCurrentL3(),
		   esd->GetMagneticField(),
		   esd->GetRunNumber()));
      
      fFirstEvent = false;
    }

    // Get the particle stack 
    AliStack* stack = mcEvent->Stack();

    // Some variables 
    UInt_t   triggers; // Trigger bits
    Bool_t   lowFlux;  // Low flux flag
    UShort_t iVz;      // Vertex bin from ESD
    Double_t vZ;       // Z coordinate from ESD
    Double_t cent;     // Centrality 
    UShort_t iVzMc;    // Vertex bin from MC
    Double_t vZMc;     // Z coordinate of IP vertex from MC
    Double_t b;        // Impact parameter
    Int_t    nPart;    // Number of participants 
    Int_t    nBin;     // Number of binary collisions 
    Double_t phiR;     // Reaction plane from MC
    UShort_t nClusters;// Number of clisters 
    // Process the data 
    Int_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, vZ, cent,
				      nClusters);
    Int_t retMC  = fInspector.ProcessMC(mcEvent, triggers, iVzMc, 
					vZMc, b, nPart, nBin, phiR);

    Bool_t hasESDVtx = retESD == AliFMDEventInspector::kOk;
    Bool_t hasMCVtx  = retMC  == AliFMDEventInspector::kOk;
    if (hasESDVtx) fVertexESD->Fill(vZ);
    if (hasMCVtx)  fVertexMC->Fill(vZMc);

    Bool_t isMcInel = true; // (triggers & AliAODForwardMult::kB);
    Bool_t isMcNSD  = (triggers & AliAODForwardMult::kMCNSD);

    Int_t mESD = 0;
    const AliMultiplicity* spdmult = esd->GetMultiplicity();
    if (!spdmult) {
      AliWarning("No SPD multiplicity");
    }
    else { 
      // Check if we have one or more tracklets 
      // in the range -1 < eta < 1 to set the INEL>0 
      // trigger flag. 
      Int_t n = spdmult->GetNumberOfTracklets();
      for (Int_t j = 0; j < n; j++) 
	if(TMath::Abs(spdmult->GetEta(j)) < 1) mESD++;
    }

    // Reset cache 
    fData->Reset();
    Int_t mMC = 0; // Number of particles in |eta|<1

    // Loop over all tracks 
    Int_t nTracks = mcEvent->GetNumberOfTracks();
    for (Int_t iTr = 0; iTr < nTracks; iTr++) { 
      AliMCParticle* particle = 
	static_cast<AliMCParticle*>(mcEvent->GetTrack(iTr));
    
      // Check the returned particle 
      if (!particle) continue;

      // Check if this charged and a primary 
      Bool_t isCharged = particle->Charge() != 0;
      Bool_t isPrimary = stack->IsPhysicalPrimary(iTr);

      if (!isCharged || !isPrimary) continue;

      
      // Fill (eta,phi) of the particle into histograsm for b
      Double_t eta = particle->Eta();
      Double_t phi = particle->Phi();
      
      fData->Fill(eta, phi);
      if (TMath::Abs(eta) <= 1) mMC++;
    }
    Int_t m = mESD;
    if (fTrackletRequirement == kMC) m = mMC;
    fM->Fill(m);

    bool isMcInelGt0 = isMcInel && (mMC > 0);
    
    bool hasVertex   = true;
    if (fVertexRequirement & kMC)  hasVertex = hasVertex && hasMCVtx;
    if (fVertexRequirement & kESD) hasVertex = hasVertex && hasESDVtx;

    if (isMcInel) {
      fTriggers->Fill(0);
      bool triggered = (triggers & AliAODForwardMult::kInel);
      if (triggered) fTriggers->Fill(1);
      fInel.AddEvent(triggered, hasVertex, m, fData);
    }
    if (isMcInel) { // && nClusters > 0) {
      fTriggers->Fill(2);
      bool triggered = (triggers & AliAODForwardMult::kNClusterGt0);
      if (triggered) fTriggers->Fill(3);
      fNClusterGt0.AddEvent(triggered, hasVertex, m, fData);
    }
    if (isMcInelGt0) {
      fTriggers->Fill(4);
      bool triggered = (triggers & AliAODForwardMult::kInelGt0);
      if (triggered) fTriggers->Fill(5);
      fInelGt0.AddEvent(triggered, hasVertex, m, fData);
    }
    if (isMcNSD) {
      fTriggers->Fill(6);
      bool triggered = (triggers & AliAODForwardMult::kNSD);
      if (triggered) fTriggers->Fill(7);
      fNSD.AddEvent(triggered, hasVertex, m, fData);
    }
    PostData(1, fList);
  }
  /** 
   * End of job processing 
   */
  void Terminate(Option_t*)
  {
    fList = dynamic_cast<TList*>(GetOutputData(1));
    if (!fList) {
      AliError(Form("No output list defined (%p)", GetOutputData(1)));
      if (GetOutputData(1)) GetOutputData(1)->Print();
      return;
    }


    TList* output = new TList;
    output->SetName("triggerResults");
    output->SetOwner();

    fVertexMC = static_cast<TH1D*>(fList->FindObject("vertexMC"));
    fVertexESD = static_cast<TH1D*>(fList->FindObject("vertexESD"));
    fM         = static_cast<TH1D*>(fList->FindObject("m"));
    if (fVertexMC) { 
      TH1D* vtxMC = static_cast<TH1D*>(fVertexMC->Clone("vertexMC"));
      vtxMC->SetDirectory(0);
      if (vtxMC->GetEntries() > 0)
	vtxMC->Scale(1. / vtxMC->GetEntries());
      else 
	vtxMC->Scale(0);
      output->Add(vtxMC);
    }
    if (fVertexESD) { 
      TH1D* vtxESD = static_cast<TH1D*>(fVertexESD->Clone("vertexESD"));
      vtxESD->SetDirectory(0);
      if (vtxESD->GetEntries() > 0)
	vtxESD->Scale(1. / vtxESD->GetEntries());
      else 
	vtxESD->Scale(0);
      output->Add(vtxESD);
    }
    if (fM) { 
      TH1D* m = static_cast<TH1D*>(fM->Clone("m"));
      m->SetDirectory(0);
      m->SetYTitle("P(N_{ch}|_{|#eta|<1} < X)");
      if (m->GetBinContent(1) > 0)
	m->Scale(1. / m->GetBinContent(1));
      else 
	m->Scale(0);
      output->Add(m);
    }      

    TString vtxReq;
    if (fVertexRequirement & kMC)  vtxReq.Append("MC ");
    if (fVertexRequirement & kESD) vtxReq.Append("ESD ");
    output->Add(new TNamed("vtxReq", vtxReq.Data()));
    output->Add(new TNamed("trkReq",
			   fTrackletRequirement == kMC ? "MC" : "ESD"));

    fInel.Finish(fList, output);
    fInelGt0.Finish(fList, output);
    fNSD.Finish(fList, output);
    fNClusterGt0.Finish(fList, output);

    PostData(2, output);
  }
    
protected:
  //__________________________________________________________________
  /** 
   * Structure to hold per trigger type information 
   */
  struct TriggerType : public TNamed
  {
    //________________________________________________________________
    /** 
     * Structure to hold per multiplicity bin information 
     */
    struct MBin : public TNamed
    {
      TH2D*     fTruth;
      TH2D*     fTriggered; 
      TH2D*     fAccepted;
      TH1I*     fCounts;
      UShort_t  fLow;
      UShort_t  fHigh;
      Bool_t IsAll() const { return fLow > fHigh; }
      Bool_t IsLast() const { return fHigh >= kMaxN; }
      /** 
       * Constructor 
       */
      MBin() 
	: fTruth(0), fTriggered(0), fAccepted(0), 
	  fCounts(0), fLow(0), fHigh(1000) {}
      /** 
       * Constructor 
       * 
       * @param p      Parent list 
       * @param low    Low cut 
       * @param high   High cut 
       * @param eAxis  Eta axis 
       * @param pAxis  Phi axis 
       */
      MBin(TList* p, UShort_t low, UShort_t high, const TH2D* dHist) 
	: fTruth(0), 
	  fTriggered(0), 
	  fAccepted(0),
	  fCounts(0), 
	  fLow(low), 
	  fHigh(high)
      {
	if (IsAll()) {
	  SetName("all");
	  SetTitle("All");
	}
	else if (IsLast()) { 
	  SetName(Form("m%03dplus", fLow));
	  SetTitle(Form("%d #leq N_{tracklets}|_{|#eta|<1}", fLow));
	}
	else {
	  SetName(Form("m%03d_%03d", fLow, fHigh));
	  SetTitle(Form("%d #leq N_{tracklets}|_{|#eta|<1} < %d", fLow,fHigh));
	}

	TList* l = new TList;
	l->SetOwner();
	l->SetName(GetName());
	p->Add(l);

	fTruth = static_cast<TH2D*>(dHist->Clone(("truth")));
	fTruth->SetTitle("MC truth");
	fTruth->SetDirectory(0);
	fTruth->SetZTitle("#sum_i^{N_X} N_{ch}(#eta,#varphi)");
	// fTruth->Sumw2();
	fTruth->Reset();
	l->Add(fTruth);

	fTriggered = static_cast<TH2D*>(fTruth->Clone(("triggered")));
	fTriggered->SetTitle("Triggered");
	fTriggered->SetDirectory(0);
	fTriggered->SetZTitle("#sum_i^{N_T} N_{ch}(#eta,#varphi)");
	// fTriggered->Sumw2();
	fTriggered->Reset();
	l->Add(fTriggered);

	fAccepted = static_cast<TH2D*>(fTruth->Clone(("accepted")));
	fAccepted->SetTitle("Accepted");
	fAccepted->SetDirectory(0);
	fAccepted->SetZTitle("#sum_i^{N_T} N_{ch}(#eta,#varphi)");
	// fAccepted->Sumw2();
	fAccepted->Reset();
	l->Add(fAccepted);
	
	fCounts = new TH1I("counts", "Event counts", 3, -.5, 2.5);
	fCounts->SetDirectory(0);
	fCounts->GetXaxis()->SetBinLabel(1, "Truth");
	fCounts->GetXaxis()->SetBinLabel(2, "Triggered");
	fCounts->GetXaxis()->SetBinLabel(3, "Accepted");
	fCounts->SetYTitle("# events");
	l->Add(fCounts);
      }
      /** 
       * Add event observation
       * 
       * @param triggered Whether the event was triggered
       * @param event     Data for this event 
       */
      void AddEvent(Bool_t triggered, Bool_t hasVtx, const TH2D* event) 
      {
	fCounts->Fill(0);
	fTruth->Add(event);
	if (triggered) { 
	  fCounts->Fill(1);
	  fTriggered->Add(event);
	  if (hasVtx) {
	    fCounts->Fill(2);
	    fAccepted->Add(event);
	  }
	}
      }
      /** 
       * End of job processing 
       * 
       * @param p      Parent list
       * @param o      Output parent list
       * @param stack  Stack of histograms
       * 
       * @return Trigger efficiency
       */ 
     Double_t Finish(const TList* p, TList* o, THStack* stack) 
      {
	TList* l = dynamic_cast<TList*>(p->FindObject(GetName()));
	if (!l) { 
	  Warning("Finish", "Cannot find %s in %s", GetName(), p->GetName());
	  return 0;
	}
	fTruth     = static_cast<TH2D*>(l->FindObject("truth"));
	fTriggered = static_cast<TH2D*>(l->FindObject("triggered"));
	fAccepted  = static_cast<TH2D*>(l->FindObject("accepted"));
	fCounts    = static_cast<TH1I*>(l->FindObject("counts"));
	
	Int_t    nTruth     = fCounts->GetBinContent(1);
	Int_t    nTriggered = fCounts->GetBinContent(2);
	Int_t    nAccepted  = fCounts->GetBinContent(3);
	Double_t eff        = 0;
	if (nTruth > 0) eff = double(nTriggered) / nTruth;
	else if (nTriggered == nTruth) eff = 1;

	if (nTruth > 0)     fTruth->Scale(1. / nTruth);
	if (nTriggered > 0) fTriggered->Scale(1. / nTriggered);
	if (nAccepted > 0)  fAccepted->Scale(1. / nAccepted);

	if (IsAll()) 
	  Info("Finish", "%-12s  [all]  E_X=N_T/N_X=%9d/%-9d=%f "
	       "E_V=N_A/N_T=%9d/%-9d=%f", 
	       p->GetName(), nTriggered, nTruth, eff, nAccepted, nTriggered, 
	       (nTriggered > 0 ? double(nAccepted) / nTriggered: 0));
	else if (IsLast()) 
	  Info("Finish", "%-12s  [%2d+]  E_X=N_T/N_X=%9d/%-9d=%f "
	       "E_V=N_A/N_T=%9d/%-9d=%f", 
	       p->GetName(), fLow, nTriggered, nTruth, eff, 
	       nAccepted, nTriggered, 
	       (nTriggered > 0 ? double(nAccepted) / nTriggered: 0));
	else
	  Info("Finish", "%-12s [%2d-%2d] E_X=N_T/N_X=%9d/%-9d=%f "
	       "E_V=N_A/N_T=%9d/%-9d=%f", 
	       p->GetName(), fLow, fHigh, nTriggered, nTruth, eff, 
	       nAccepted, nTriggered, 
	       (nTriggered > 0 ? double(nAccepted) / nTriggered: 0));
	
	TList* out = new TList;
	out->SetName(GetName());
	out->SetOwner();
	o->Add(out);
	
	out->Add(fTruth);
	out->Add(fTriggered);
	out->Add(new TParameter<double>("eff", eff));
	
	TH2D* bias = static_cast<TH2D*>(fAccepted->Clone("bias"));
	bias->Divide(fTruth);
	bias->SetDirectory(0);
	bias->SetZTitle("Trigger bias (accepted/truth)");
	out->Add(bias);

	TString title = GetTitle();
	TH1D* truth_px = static_cast<TH1D*>(fTruth->ProjectionX("truth_eta"));
	truth_px->SetTitle(title);
	truth_px->Scale(1. / fTruth->GetNbinsY());
	truth_px->SetDirectory(0);
	truth_px->SetLineColor(kRed+1);
	truth_px->SetMarkerColor(kRed+1);
	truth_px->SetFillColor(kRed+1);
	truth_px->SetFillStyle(0);
	out->Add(truth_px);

	TH1D* triggered_px = 
	  static_cast<TH1D*>(fTriggered->ProjectionX("triggered_eta"));
	triggered_px->SetTitle(title);
	triggered_px->Scale(1. / fTriggered->GetNbinsY());
	triggered_px->SetDirectory(0);
 	triggered_px->SetLineColor(kGreen+1);
	triggered_px->SetMarkerColor(kGreen+1);
	triggered_px->SetFillColor(kGreen+1);
	triggered_px->SetFillStyle(0);
	out->Add(triggered_px);

	TH1D* accepted_px = 
	  static_cast<TH1D*>(fAccepted->ProjectionX("accepted_eta"));
	accepted_px->SetTitle(title);
	accepted_px->Scale(1. / fAccepted->GetNbinsY());
 	accepted_px->SetLineColor(kBlue+1);
	accepted_px->SetMarkerColor(kBlue+1);
	accepted_px->SetFillColor(kBlue+1);
	accepted_px->SetDirectory(0);
	out->Add(accepted_px);

	THStack* data = new THStack("data", "Data distributions");
	data->Add(truth_px);
	data->Add(triggered_px);
	data->Add(accepted_px);
	out->Add(data);

#if 0
	TH1D* bias_px = static_cast<TH1D*>(bias->ProjectionX("bias_eta"));
	bias_px->SetTitle(title);
	bias_px->Scale(1. / bias->GetNbinsY());
#else
	TH1D* bias_px = static_cast<TH1D*>(accepted_px->Clone("bias_px"));
	bias_px->Divide(truth_px);
	bias_px->SetYTitle("Trigger bias (triggered/truth)");
#endif
	bias_px->SetDirectory(0);
	bias_px->SetMarkerStyle(20);
	bias_px->SetFillStyle(0);
	bias_px->SetMinimum(0);
	out->Add(bias_px);

	stack->Add(bias_px);

	return eff;
      }	
      ClassDef(MBin,1);
    };

    /** 
     * Constructor 
     */
    TriggerType() 
      : TNamed(), 
	fMask(0),
	fM(0),
	fBins(0)
    {}
    //--- Constructor ------------------------------------------------
    /** 
     * Constructor 
     * 
     * @param mask  Trigger mask
     */
    TriggerType(UShort_t mask) 
      : TNamed(AliAODForwardMult::GetTriggerString(mask), ""),
	fMask(mask), 
	fM(0), 
	fBins(0)
    {
      TString n(GetName());
      n = n.Strip(TString::kBoth);
      SetName(n);
    }
    /** 
     * Create objects 
     * 
     * @param list   PArent list
     * @param mAxis  Multiplicity axis 
     * @param eAxis  Eta axis
     * @param pAxis  Phi axis
     */
    void CreateObjects(TList* list, 
		       const TH1D* mHist, 
		       const TH2D* dHist)
    {
      TList* ours  = new TList;
      ours->SetName(GetName());
      ours->SetOwner();
      list->Add(ours);

      fM = static_cast<TH1D*>(mHist->Clone("m"));
      fM->SetDirectory(0);
      ours->Add(fM);
      
      fBins = new TObjArray;
      fBins->AddAt(new MBin(ours, 1000, 0, dHist), 0);
      for (Int_t i = 1; i < fM->GetNbinsX(); i++) { 
	Double_t low  = fM->GetXaxis()->GetBinLowEdge(i);
	Double_t high = fM->GetXaxis()->GetBinUpEdge(i);
	Info("CreatePbjects", "Adding bin [%3d,%3d] @ %d", 
	     int(low), int(high), i);
	fBins->AddAt(new MBin(ours, UShort_t(low), UShort_t(high), dHist), i);
      }
    }
    /** 
     * Find bin corresponding to m
     * 
     * @param m Multiplicity 
     * 
     * @return Bin. 
     */
    MBin* FindBin(UShort_t m) 
    { 
#if 0
      for (Int_t i = 1; i < fBins->GetEntries(); i++) {
	MBin* b = static_cast<MBin*>(fBins->At(i));
	if (m >= b->fLow && m < b->fHigh) return b;
      }
      Warning("FindBin", "Found no bin for m=%d", m);
      return 0;
#else
      Int_t b = fM->GetXaxis()->FindBin(m);
      if (b <= 0) return 0;
      if (b >= fM->GetNbinsX()+1) b = fM->GetNbinsX();
      MBin* bin = static_cast<MBin*>(fBins->At(b));
      return bin;
#endif
    }
    /** 
     * Add event observation 
     * 
     * @param triggered  IF this is triggered
     * @param m          Multiplicity 
     * @param data       Observation 
     */
    void AddEvent(Bool_t triggered, Bool_t hasVtx, UShort_t m, const TH2D* data)
    {
      fM->AddBinContent(1);
      fM->AddBinContent(TMath::Min(fM->GetNbinsX(), m+2));

      MBin* all = static_cast<MBin*>(fBins->At(0));
      all->AddEvent(triggered, hasVtx, data);
      
      MBin* bin = FindBin(m);
      if (!bin) return;
      bin->AddEvent(triggered, hasVtx, data);      
    }      
    /** 
     * End of job processing 
     * 
     * @param p Parent list 
     * @param o Parent output list 
     */
    void Finish(const TList* p, TList* o)
    {
      TList* l = dynamic_cast<TList*>(p->FindObject(GetName()));
      if (!l) { 
	Warning("Finish", "Cannot find %s in %s", GetName(), p->GetName());
	return;
      }
      
      TList* ours  = new TList;
      ours->SetName(GetName());
      ours->SetOwner();
      o->Add(ours);

      fM = static_cast<TH1D*>(l->FindObject("m"));
      if (!fM) { 
	Warning("Finish", "Didn't find object 'm' in %s", l->GetName());
	return;
      }
      TH1D* m = static_cast<TH1D*>(fM->Clone("m"));
      m->SetDirectory(0);
      m->SetYTitle("P(N_{ch}|_{|#eta|<1} < X)");
      if (m->GetBinContent(1) > 0) 
	m->Scale(1. / m->GetBinContent(1));
      ours->Add(m);

      Int_t nBin = fM->GetNbinsX();
      TH1D* effs = static_cast<TH1D*>(fM->Clone("effs"));
      effs->SetYTitle("#epsilon_{X}");
      effs->SetFillColor(kRed+1);
      effs->SetDirectory(0);
      effs->SetMinimum(0);

      gStyle->SetPalette(1);
      Int_t   nCol = gStyle->GetNumberOfColors();
      THStack* stack = new THStack("biases", "Trigger biases");
      for (Int_t i = 0; i < nBin; i++) { 
	MBin* bin = static_cast<MBin*>(fBins->At(i));
	effs->SetBinContent(i+1, bin->Finish(l, ours, stack));
	TH1* h = static_cast<TH1*>(stack->GetHists()->At(i));
	Int_t col = kBlack;
	if (i != 0) { 
	  Int_t icol = TMath::Min(nCol-1,int(double(i)/nBin * nCol + .5));
	  col        = gStyle->GetColorPalette(icol);
	}
	h->SetMarkerColor(col);
	h->SetFillColor(col);
	h->SetLineColor(col);
      }

      ours->Add(stack);
      ours->Add(effs);
    } 
    UShort_t   fMask;
    TH1D*      fM;
    TObjArray* fBins;
    ClassDef(TriggerType,1);
  };
  TriggerType            fInel;
  TriggerType            fInelGt0;
  TriggerType            fNSD;
  TriggerType            fNClusterGt0;
  AliFMDMCEventInspector fInspector;
  TList*                 fList;
  Bool_t                 fFirstEvent;
  TH2D*                  fData;
  TH1I*                  fTriggers;
  UInt_t                 fTrackletRequirement;
  UInt_t                 fVertexRequirement;
  TAxis                  fVertexAxis;
  TH1D*                  fVertexESD;
  TH1D*                  fVertexMC;
  TH1D*                  fM;
  ClassDef(EvaluateTrigger,1);
};
#else 
//====================================================================
void MakeEvaluateTriggers(const char* esddir, 
			  Int_t       nEvents    = -1, 
			  UInt_t      vtx        = 0x2,
			  UInt_t      trk        = 0x1,
			  Int_t       vz         = 10,
			  Int_t       proof      = 0)
{
  // --- Libraries to load -------------------------------------------
  gROOT->Macro("$ALICE_ROOT/PWGLF/FORWARD/analysis2/scripts/LoadLibs.C");

  // --- Check for proof mode, and possibly upload pars --------------
  if (proof> 0) { 
    gROOT->LoadMacro("$ALICE_ROOT/PWGLF/FORWARD/analysis2/scripts/LoadPars.C");
    if (!LoadPars(proof)) { 
      Error("MakeAOD", "Failed to load PARs");
      return;
    }
  }
  
  // --- Our data chain ----------------------------------------------
  gROOT->LoadMacro("$ALICE_ROOT/PWGLF/FORWARD/analysis2/scripts/MakeChain.C");
  TChain* chain = MakeChain("ESD", esddir,true);
  // If 0 or less events is select, choose all 
  if (nEvents <= 0) nEvents = chain->GetEntries();
  
  // --- Set the macro path ------------------------------------------
  gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWGLF/FORWARD/analysis2:"
			   "$ALICE_ROOT/ANALYSIS/macros",
			   gROOT->GetMacroPath()));

  // --- Creating the manager and handlers ---------------------------
  AliAnalysisManager *mgr  = new AliAnalysisManager("Triggers", 
						    "Forward multiplicity");
  AliAnalysisManager::SetCommonFileName("triggers.root");

  // --- ESD input handler -------------------------------------------
  AliESDInputHandler *esdHandler = new AliESDInputHandler();
  mgr->SetInputEventHandler(esdHandler);      
       
  // --- Monte Carlo handler -----------------------------------------
  AliMCEventHandler* mcHandler = new AliMCEventHandler();
  mgr->SetMCtruthEventHandler(mcHandler);
  mcHandler->SetReadTR(true);    

  // --- Add tasks ---------------------------------------------------
  // Physics selection 
  gROOT->Macro(Form("AddTaskPhysicsSelection.C(1,1,0)"));

#if 0
  // --- Fix up physics selection to give proper A,C, and E triggers -
  AliInputEventHandler* ih =
    static_cast<AliInputEventHandler*>(mgr->GetInputEventHandler());
  AliPhysicsSelection* ps = 
    static_cast<AliPhysicsSelection*>(ih->GetEventSelection());
  // Ignore trigger class when selecting events.  This mean that we
  // get offline+(A,C,E) events too
  ps->SetSkipTriggerClassSelection(true);
#endif

  // --- compile our code --------------------------------------------
  gSystem->AddIncludePath("-I${ALICE_ROOT}/PWGLF/FORWARD/analysis2 "
			  "-I${ALICE_ROOT}/ANALYSIS "
			  "-I${ALICE_ROOT}/include -DBUILD=1");
  gROOT->LoadMacro("${ALICE_ROOT}/PWGLF/FORWARD/analysis2/MakeEvaluateTriggers.C++g");
  
  // --- Make our object ---------------------------------------------
  EvaluateTrigger* task = new EvaluateTrigger("triggers");
  mgr->AddTask(task);
  task->SetVertexRequirement(vtx);
  task->SetTrackletRequirement(trk);
  task->SetVertexAxis(10, -vz, vz);

  // --- create containers for input/output --------------------------
  AliAnalysisDataContainer *sums = 
    mgr->CreateContainer("triggerSums", TList::Class(), 
			 AliAnalysisManager::kOutputContainer, 
			 AliAnalysisManager::GetCommonFileName());
  AliAnalysisDataContainer *output = 
    mgr->CreateContainer("triggerResults", TList::Class(), 
			 AliAnalysisManager::kParamContainer, 
			 AliAnalysisManager::GetCommonFileName());
  
  // --- connect input/output ----------------------------------------
  mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
  mgr->ConnectOutput(task, 1, sums);
  mgr->ConnectOutput(task, 2, output);
  
  // --- Run the analysis --------------------------------------------
  TStopwatch t;
  if (!mgr->InitAnalysis()) {
    Error("MakeAOD", "Failed to initialize analysis train!");
    return;
  }
  // Skip terminate if we're so requested and not in Proof or full mode
  mgr->SetSkipTerminate(false);
  // Some informative output 
  mgr->PrintStatus();
  if (proof) mgr->SetDebugLevel(3);
  if (mgr->GetDebugLevel() < 1 && !proof) 
    mgr->SetUseProgressBar(kTRUE,100);

  // Run the train 
  t.Start();
  Printf("=== RUNNING ANALYSIS on %9d events ==========================",
	 nEvents);
  mgr->StartAnalysis(proof ? "proof" : "local", chain, nEvents);
  t.Stop();
  t.Print();
}
//====================================================================
void DrawEvaluateTriggers(const char* filename="triggers.root")
{
  TFile* file = TFile::Open(filename, "READ");
  if (!file) { 
    Error("DrawEvaluteTriggers", "Failed to open %s", filename);
    return;
  }

  TList* list = static_cast<TList*>(file->Get("triggerResults"));
  if (!list) { 
    Error("DrawEvaluateTriggers", "Faield to get triggerResults from %s", 
	  list->GetName());
    return;
  }

  TCanvas* c = new TCanvas("c", "c", 1200, 900);
  c->SetFillColor(0);
  c->SetBorderSize(0);
  c->SetBorderMode(0);

  TPad* p = new TPad("top", "Top", 0, .9, 1, 1, 0, 0, 0); 
  p->SetFillColor(0);
  p->SetFillStyle(0);
  p->SetBorderSize(0);
  p->SetBorderMode(0);
  c->cd();
  p->Draw();
  p->cd();
  TObject* vtxReq = list->FindObject("vtxReq");
  TObject* trkReq = list->FindObject("trkReq");
  if (!vtxReq) vtxReq = new TNamed("vtxReq", "Unknown");
  if (!trkReq) trkReq = new TNamed("trkReq", "Unknown");
  TLatex* ltx = new TLatex(0.5, 0.5, 
			   Form("Trigger bias and efficiency. " 
				"Vertex requirement: %s, "
				"Tracklet requirement %s",
				vtxReq->GetTitle(), trkReq->GetTitle()));
  ltx->Draw();
  ltx->SetTextAlign(22);
  ltx->SetTextSize(0.4);

  const char*  trigs[] = { "INEL", "NCluster>0", "INEL>0", "NSD", 0 };
  const char** trig    = trigs;
  Int_t n = 4;
  Int_t i = 0;
  while (*trig) { 
    TList* lt = static_cast<TList*>(list->FindObject(*trig));
    if (!lt) { 
      Warning("DrawEvaluteTriggers", "No list for '%s' in '%s'", 
	      *trig, list->GetName());
      // list->ls();
      TIter next(triggerResults);
      TObject* o ;
      while ((o = next())) Printf("'%s'", o->GetName());
      trig++;
      i++;
      continue;
    }

    Double_t y1 = 1-double(i+1)/n;
    Double_t y2 = 1-double(i)/n;
    Double_t x1 = double(i)/n;
    Double_t x2 = double(i+1)/n;
    p = new TPad(Form("p%d", 2*i), Form("%s biases", *trig), 
		 x1, .3, x2, .9, 0, 0);
    p->SetFillColor(0);
    p->SetFillStyle(0);
    p->SetBorderSize(0);
    p->SetBorderMode(0);
    p->SetTopMargin(0.01);
    p->SetBottomMargin(0.05);
    p->SetRightMargin(0.01);
    c->cd();
    p->Draw();
    p->cd();
    
    THStack* biases = static_cast<THStack*>(lt->FindObject("biases"));
    biases->SetTitle(*trig);
    TLegend* l = new TLegend(.1, .15, .95, .35, 
			     Form("1/N_{T}#sum N_{ch}}/"
				  "1/N_{X}#sum N_{ch} - %s", *trig));
    l->SetFillColor(0);
    l->SetFillStyle(0);
    l->SetBorderSize(0);
    l->SetNColumns(2);
    
    gStyle->SetOptTitle(0);
    TIter next(biases->GetHists());
    TH1*  frame = 0;
    TH1*  hist = 0;
    Int_t j = 1;
    while ((hist = static_cast<TH1*>(next()))) { 
      if (!frame) { 
	
	frame = new TH2D("frame", "", 
			 hist->GetXaxis()->GetNbins(),
			 hist->GetXaxis()->GetXmin(),
			 hist->GetXaxis()->GetXmax(),
			 100, 0, 4);
	frame->SetDirectory(0);
	frame->SetXTitle(hist->GetXaxis()->GetTitle());
	frame->SetYTitle(hist->GetYaxis()->GetTitle());
	frame->SetStats(0);
	frame->Draw();
      }
      // hist->Scale(j);
      // hist->Draw("same p hist");
      // if (j == 1) hist->SetTitle("all");
      l->AddEntry(hist, hist->GetTitle(), "p");
      j++;
    }
    // biases->SetMaximum(3.5);
    biases->Draw("p hist nostack");
    l->Draw();
    
    p = new TPad(Form("p%d", 2*i+1), Form("%s efficiencies", *trig), 
		 x1, 0, x2, .3, 0, 0, 0);
    p->SetFillColor(0);
    p->SetFillStyle(0);
    p->SetBorderSize(0);
    p->SetBorderMode(0);
    p->SetTopMargin(0.01);
    p->SetRightMargin(0.01);
    c->cd();
    p->Draw();
    p->cd();
    
    TH1* effs = static_cast<TH1*>(lt->FindObject("effs"));
    effs->SetTitle(*trig);
    effs->SetStats(0);
    effs->SetMinimum(0);
    effs->SetMaximum(1.4);
    effs->SetMarkerSize(2.5);
    effs->Draw("hist text");


    i++;
    trig++;
  }
  c->cd();

  c->SaveAs("triggers.png");
}

#endif