H2D.cpp

Go to the documentation of this file.

#ifdef __ICC
// disable icc remark #2259: non-pointer conversion from "X" to "Y" may lose significant bits
//   TODO: To be removed, since it comes from ROOT TMathBase.h
#pragma warning(disable:2259)
#endif
#ifdef WIN32
// Disable warning
//   warning C4996: 'sprintf': This function or variable may be unsafe.
// coming from TString.h
#pragma warning(disable:4996)
#endif
#include <array>
#include <GaudiKernel/ObjectFactory.h>
#include <GaudiCommonSvc/HistogramUtility.h>
#include <GaudiCommonSvc/H2D.h>
#include <GaudiCommonSvc/H1D.h>
#include <GaudiCommonSvc/P1D.h>
#include <TH2D.h>
#include <TProfile.h>
#include "GaudiPI.h"

namespace {
   using AIDA::IProfile1D;
   using AIDA::IHistogram1D;
   using AIDA::IHistogram2D;
}

std::pair<DataObject*,IHistogram2D*> Gaudi::createH2D(const std::string & title,int binsX,double iminX,double imaxX,int binsY,double iminY,double imaxY) {
  auto p = new Histogram2D(new TH2D(title.c_str(),title.c_str(),binsX, iminX, imaxX, binsY, iminY, imaxY));
  return { p, p };
}

std::pair<DataObject*,IHistogram2D*> Gaudi::createH2D(const std::string & title, const Edges& eX, const Edges& eY) {
  auto p = new Histogram2D(new TH2D(title.c_str(),title.c_str(),eX.size()-1,&eX.front(),eY.size()-1,&eY.front()));
  return { p, p };
}

std::pair<DataObject*,IHistogram2D*> Gaudi::createH2D(TH2D* rep) {
  auto p = new Histogram2D(rep);
  return { p, p };
}

std::pair<DataObject*,IHistogram2D*> Gaudi::createH2D(const IHistogram2D& hist)  {
  TH2D *h = getRepresentation<AIDA::IHistogram2D,TH2D>(hist);
  Histogram2D *n = h ? new Histogram2D(new TH2D(*h)) : nullptr;
  return { n, n };
}

std::pair<DataObject*,IHistogram1D*>
Gaudi::slice1DX(const std::string& nam,const IHistogram2D& hist,int first, int last)  {
  TH2 *r = getRepresentation<IHistogram2D,TH2>(hist);
  TH1D *t = r ? r->ProjectionX("_px",first,last,"e") : nullptr;
  if ( t ) t->SetName(nam.c_str());
  Histogram1D* p = ( t ? new Histogram1D(t) : nullptr );
  return { p, p };
}

std::pair<DataObject*,IHistogram1D*>
Gaudi::slice1DY(const std::string& nam,const IHistogram2D& hist,int first, int last)  {
  TH2  *r = getRepresentation<IHistogram2D,TH2>(hist);
  TH1D *t = r ? r->ProjectionY("_py",first,last,"e") : nullptr;
  if ( t ) t->SetName(nam.c_str());
  Histogram1D* p = ( t ? new Histogram1D(t) : nullptr );
  return { p, p };
}

std::pair<DataObject*,IHistogram1D*>
Gaudi::project1DY(const std::string& nam,const IHistogram2D& hist,int first,int last) {
  TH2 *r = getRepresentation<IHistogram2D,TH2>(hist);
  TH1D *t = r ? r->ProjectionY("_px",first,last,"e") : nullptr;
  if ( t ) t->SetName(nam.c_str());
  Histogram1D* p = ( t ? new Histogram1D(t) : nullptr );
  return { p, p };
}

std::pair<DataObject*,IProfile1D*>
Gaudi::profile1DX(const std::string& nam,const IHistogram2D& hist,int first,int last)  {
  TH2 *r = Gaudi::getRepresentation<IHistogram2D,TH2>(hist);
  TProfile *t = r ? r->ProfileX("_pfx",first,last,"e") : nullptr;
  if ( t ) t->SetName(nam.c_str());
  Profile1D* p = ( t ? new Profile1D(t) : nullptr );
  return { p, p };
}

std::pair<DataObject*,IProfile1D*>
Gaudi::profile1DY(const std::string& nam, const IHistogram2D& hist,int first, int last) {
  TH2 *r = getRepresentation<IHistogram2D,TH2>(hist);
  TProfile *t = r ? r->ProfileY("_pfx",first,last,"e") : nullptr;
  if ( t ) t->SetName(nam.c_str());
  Profile1D* p = ( t ? new Profile1D(t) : nullptr );
  return { p, p };
}


namespace Gaudi {
  template <>
  void * Generic2D<IHistogram2D,TH2D>::cast(const std::string & className) const {
    if (className == "AIDA::IHistogram2D")
      return (IHistogram2D*)this;
    else if (className == "AIDA::IHistogram")
      return (IHistogram*)this;
    return nullptr;
  }

  template <>
  int Generic2D<IHistogram2D,TH2D>::binEntries(int indexX,int indexY) const {
    if (binHeight(indexX, indexY)<=0) return 0;
    double xx =  binHeight(indexX, indexY)/binError(indexX, indexY);
    return int(xx*xx+0.5);
  }

  template <>
  void Generic2D<AIDA::IHistogram2D,TH2D>::adoptRepresentation(TObject* rep)  {
    TH2D* imp = dynamic_cast<TH2D*>(rep);
    if ( !imp )  throw std::runtime_error("Cannot adopt native histogram representation.");
    m_rep.reset( imp );
    m_xAxis.initialize(m_rep->GetXaxis(),true);
    m_yAxis.initialize(m_rep->GetYaxis(),true);
    const TArrayD* a = m_rep->GetSumw2();
    if ( !a || (a && a->GetSize()==0) ) m_rep->Sumw2();
    setTitle(m_rep->GetTitle());
  }
}

Gaudi::Histogram2D::Histogram2D()
    : Base( new TH2D() )
{
  m_rep->Sumw2();
  m_sumwx = m_sumwy = 0;
  setTitle("");
  m_rep->SetDirectory(nullptr);
}

Gaudi::Histogram2D::Histogram2D(TH2D* rep)
{
  adoptRepresentation(rep);
  m_sumwx = m_sumwy = 0;
  m_rep->SetDirectory(nullptr);
}

bool Gaudi::Histogram2D::setBinContents( int i,int j,int entries,double height,double error,double centreX,double centreY) {
  m_rep->SetBinContent(rIndexX(i), rIndexY(j), height);
  m_rep->SetBinError(rIndexX(i), rIndexY(j), error);
  // accumulate sumwx for in range bins
  if (i >=0 && j >= 0) {
    m_sumwx += centreX*height;
    m_sumwy += centreY*height;
  }
  m_sumEntries += entries;
  return true;
}

bool  Gaudi::Histogram2D::reset() {
  m_sumwx = 0;
  m_sumwy = 0;
  return Base::reset();
}

#ifdef __ICC
// disable icc remark #1572: floating-point equality and inequality comparisons are unreliable
//   The comparison are meant
#pragma warning(push)
#pragma warning(disable:1572)
#endif
bool Gaudi::Histogram2D::fill ( double x,double y,double weight) {
  (weight==1.) ? m_rep->Fill(x,y) : m_rep->Fill(x,y,weight );
  return true;
}

bool Gaudi::Histogram2D::setRms(double rmsX,double rmsY) {
  m_rep->SetEntries(m_sumEntries);
  std::vector<double> stat(11);
  stat[0] =  sumBinHeights();
  stat[1] = 0;
  if(equivalentBinEntries() != 0)
    stat[1] = (sumBinHeights() * sumBinHeights()) / equivalentBinEntries();
  stat[2] = m_sumwx;
  double meanX = 0;
  if(sumBinHeights() != 0) meanX =  m_sumwx/ sumBinHeights();
  stat[3] = (meanX*meanX + rmsX*rmsX) * sumBinHeights();
  stat[4] = m_sumwy;
  double meanY = 0;
  if(sumBinHeights() != 0) meanY =  m_sumwy/ sumBinHeights();
  stat[5] = (meanY*meanY + rmsY*rmsY) * sumBinHeights();
  stat[6] = 0;
  m_rep->PutStats(&stat.front());
  return true;
}

void Gaudi::Histogram2D::copyFromAida(const IHistogram2D& h) {
  // implement here the copy
  const char* tit = h.title().c_str();
  if (h.xAxis().isFixedBinning() && h.yAxis().isFixedBinning()  )
    m_rep.reset( new TH2D(tit,tit,
                     h.xAxis().bins(),h.xAxis().lowerEdge(),h.xAxis().upperEdge(),
                     h.yAxis().bins(),h.yAxis().lowerEdge(),h.yAxis().upperEdge() ) );
  else {
    Edges eX, eY;
    for (int i =0; i < h.xAxis().bins(); ++i)
      eX.push_back(h.xAxis().binLowerEdge(i));
    // add also upperedges at the end
    eX.push_back(h.xAxis().upperEdge() );
    for (int i =0; i < h.yAxis().bins(); ++i)
      eY.push_back(h.yAxis().binLowerEdge(i));
    // add also upperedges at the end
    eY.push_back(h.yAxis().upperEdge() );
    m_rep.reset( new TH2D(tit,tit,eX.size()-1,&eX.front(),eY.size()-1,&eY.front()) );
  }
  m_xAxis.initialize(m_rep->GetXaxis(),true);
  m_yAxis.initialize(m_rep->GetYaxis(),true);
  m_rep->Sumw2();
  m_sumEntries = 0;
  m_sumwx = 0;
  m_sumwy = 0;
  // statistics
  double sumw = h.sumBinHeights();
  double sumw2 = 0;
  if (h.equivalentBinEntries() != 0)
    sumw2 = ( sumw * sumw ) /h.equivalentBinEntries();
  double sumwx = h.meanX()*h.sumBinHeights();
  double sumwx2 = (h.meanX()*h.meanX() + h.rmsX()*h.rmsX() )*h.sumBinHeights();
  double sumwy = h.meanY()*h.sumBinHeights();
  double sumwy2 = (h.meanY()*h.meanY() + h.rmsY()*h.rmsY() )*h.sumBinHeights();
  double sumwxy = 0;

  // copy the contents in  (AIDA underflow/overflow are -2,-1)
  for (int i=-2; i < xAxis().bins(); ++i) {
    for (int j=-2; j < yAxis().bins(); ++j) {
      // root binning starts from one !
      m_rep->SetBinContent(rIndexX(i), rIndexY(j), h.binHeight(i,j) );
      m_rep->SetBinError(rIndexX(i), rIndexY(j), h.binError(i,j) );
      // calculate statistics
      if ( i >= 0 && j >= 0) {
        sumwxy += h.binHeight(i,j)*h.binMeanX(i,j)*h.binMeanY(i,j);
      }
    }
  }
  // need to do set entries after setting contents otherwise root will recalculate them
  // taking into account how many time  SetBinContents() has been called
  m_rep->SetEntries(h.allEntries());
  // fill stat vector
  std::array<double,11> stat = {{ sumw,  sumw2,
                                  sumwx, sumwx2,
                                  sumwy, sumwy2,
                                  sumwxy }};
  m_rep->PutStats(stat.data());
}

#ifdef __ICC
// re-enable icc remark #1572
#pragma warning(pop)
#endif