[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / older / DrawResKeep.C

#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 &lt;<i>base</i>&gt;<tt>_hists.root</tt> 
 * containing the histograms generated by AliForwardMultiplicity and 
 * the file &lt;<i>base</i>&gt;<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 &lt;<i>base</i>&gt;<tt>_hists.root</tt> 
   * containing the histograms generated by AliForwardMultiplicity and 
   * the file &lt;<i>base</i>&gt;<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 &lt;<i>base</i>&gt;<tt>_hists.root</tt>
   * containing the histograms generated by AliForwardMultiplicity and
   * the file &lt;<i>base</i>&gt;<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
//