CountHist.cxx

Go to the documentation of this file.

// $Id: CountHist.cxx,v 1.6 2008/05/27 18:22:48 rustem Exp $

//C++
#include <cmath>
#include <cassert>
#include <iostream>

//ROOT
#include "TH1.h"

//Local
#include "PhysicsNtuple/Hist1d.h"
#include "CountHist.h"

using namespace std;

//---------------------------------------------------------------------------------------
Anp::CountHist::CountHist(int iactn_, int reson_, Particle::Particle_t particle_)
   :hall(0),
    hsel(0),
    hrat(0),
    nmin(1000.0),
    particle(particle_),
    iactn(iactn_),
    reson(reson_),
    verbose(false)
{
}

//------------------------------------------------------------------------------------------
Anp::CountHist::~CountHist()
{
}

//------------------------------------------------------------------------------------------
void Anp::CountHist::Fill(const double value, const double weight_all, const double weight_sel)
{
   if(!CountHist::Valid())
   {
      return;
   }

   if(weight_all < 0.0 || weight_sel < 0.0)
   {
      cerr << "CountHist::Fill - event with negative weight" << endl;
   }

   hall -> Fill(value, weight_all);
   hsel -> Fill(value, weight_sel);
}

//------------------------------------------------------------------------------------------
TH1* Anp::CountHist::ComputeRatio(const string &option)
{
   if(!Valid())
   {
      return 0;
   }

   CountHist::SetRatio(hsel, hall, hrat);

   if(option.find("merge") != string::npos && nmin > 0.0)
   {
      Anp::Hist1d<double> hist(*hsel);
      hist.Merge(nmin);

      TH1 *HAll = CountHist::Refill(hist, hall);
      TH1 *HSel = CountHist::Refill(hist, hsel);
      TH1 *HRat = CountHist::Refill(hist, hrat);
      
      const string name = std::string(hrat -> GetName()) + "_merge";

      HRat -> Reset();
      Anp::SetDir(HRat, hrat -> GetDirectory(), name);

      SetRatio(HSel, HAll, HRat);

      delete HAll;
      delete HSel;
   }  

   return hrat;
}

//------------------------------------------------------------------------------------------
bool Anp::CountHist::SetRatio(TH1 *hN, TH1 *hD, TH1 *hR)
{
   if(!hN || !hD || !hR)
   {
      return false;
   }

   const int nbin = hN -> GetNbinsX();
   
   if(hD -> GetNbinsX() != nbin)
   {
      cerr << "CountHist::SetRatio - mismatched number of bins" << endl;
      return false;
   }

   for(int ibin = 1; ibin <= nbin; ++ibin)
   {      
      const double valueN = hN -> GetBinContent(ibin);
      const double valueD = hD -> GetBinContent(ibin);

      if(!(valueN > 0.0) || !(valueD > 0.0))
      {
         continue;
      }

      if(valueD < valueN)
      {
         if(verbose)
         {
         cerr << "CountHist::SetRatio - denominator is smaller than nominator " 
              << valueN << "/" << valueD << " = " << valueN/valueD << endl;
         }

         continue;
      }

      const double ratio = valueN/valueD;           
      const double error = std::sqrt(valueN*(1.0 - valueN/valueD))/valueD;

      hR -> SetBinContent(ibin, ratio);
      hR -> SetBinError(ibin, error);
   }   

   return true;
}

//------------------------------------------------------------------------------------------
TH1* Anp::CountHist::Refill(Hist1d<double> &hist, TH1 *hO)
{
   if(!hO)
   {
      return 0;
   }

   hist.Reset();
   hist.Fill(*hO);
   
   TH1 *hN = Anp::CreateTH1<double>(hist, hO -> GetName());
   Anp::SetDir(hN, 0, hO -> GetName());

   hN -> SetTitle(hO -> GetTitle());
   hN -> GetXaxis() -> SetTitle(hO -> GetXaxis() -> GetTitle());
   hN -> GetYaxis() -> SetTitle(hO -> GetYaxis() -> GetTitle());
   hN -> GetXaxis() -> CenterTitle();
   hN -> GetYaxis() -> CenterTitle();

   return hN;
}

//---------------------------------------------------------------------------------------
bool Anp::CountHist::Valid() const
{
   if(hall && hsel && hrat)
   {
      return true;
   }

   return false;
}

//---------------------------------------------------------------------------------------
Particle::Particle_t Anp::CountHist::GetParticle() const
{
   return particle;
}

//---------------------------------------------------------------------------------------
int Anp::CountHist::GetInteraction() const
{
   return iactn;
}

//---------------------------------------------------------------------------------------
int Anp::CountHist::GetReson() const
{
   return reson;
}

//---------------------------------------------------------------------------------------
bool Anp::operator==(const Handle<CountHist> &lhs, const Handle<CountHist> &rhs)
{
   if(!lhs.valid() || !rhs.valid())
   {
      return false;
   }

   if(lhs -> GetInteraction() != rhs -> GetInteraction())
   {
      return false;
   }

   if(lhs -> GetReson() != rhs -> GetReson())
   {
      return false;
   }

   if(lhs -> GetParticle() != rhs -> GetParticle())
   {
      return false;
   }   
   
   return true;
}

---