| # Gaudi Framework, version v21r6 |

Home |  | Generated: 11 Nov 2009

# H1D.cpp

Go to the documentation of this file.

    00001 // $Id: H1D.cpp,v 1.14 2007/07/16 13:36:17 hmd Exp $
    00002 #include "H1D.h"
    00003 #include "GaudiPI.h"
    00004 #include "GaudiKernel/StreamBuffer.h"
    00005 #include "GaudiKernel/ObjectFactory.h"
    00006 typedef Gaudi::Histogram1D H1D;
    00007 DECLARE_DATAOBJECT_FACTORY(H1D)
    00008
    00009 std::pair<DataObject*,AIDA::IHistogram1D*> Gaudi::createH1D(const std::string& title,int nBins,double xlow, double xup)  {
    00010   Histogram1D* p = new Histogram1D(new TH1D(title.c_str(),title.c_str(),nBins,xlow,xup));
    00011   return std::pair<DataObject*,AIDA::IHistogram1D*>(p,p);
    00012 }
    00013
    00014 std::pair<DataObject*,AIDA::IHistogram1D*> Gaudi::createH1D(const std::string& title, const Edges& e)  {
    00015   Histogram1D* p = new Histogram1D(new TH1D(title.c_str(),title.c_str(),e.size()-1,&e.front()));
    00016   return std::pair<DataObject*,AIDA::IHistogram1D*>(p,p);
    00017 }
    00018
    00019 std::pair<DataObject*,AIDA::IHistogram1D*> Gaudi::createH1D(const AIDA::IHistogram1D& hist)  {
    00020   TH1D *h = getRepresentation<AIDA::IHistogram1D,TH1D>(hist);
    00021   Histogram1D *n = h ? new Histogram1D(new TH1D(*h)) : 0;
    00022   return std::pair<DataObject*,AIDA::IHistogram1D*>(n,n);
    00023 }
    00024 namespace Gaudi {
    00025   template<> void *Generic1D<AIDA::IHistogram1D,TH1D>::cast(const std::string & className) const  {
    00026     if (className == "AIDA::IHistogram1D")
    00027       return const_cast<AIDA::IHistogram1D*>((AIDA::IHistogram1D*)this);
    00028     if (className == "AIDA::IHistogram")
    00029       return const_cast<AIDA::IHistogram*>((AIDA::IHistogram*)this);
    00030     return 0;
    00031   }
    00032
    00033   template<> int Generic1D<AIDA::IHistogram1D,TH1D>::binEntries (int index) const  {
    00034     if (binHeight(index)<=0) return 0;
    00035     double xx =  binHeight(index)/binError(index);
    00036     return int(xx*xx+0.5);
    00037   }
    00038
    00039   template <>
    00040   void Generic1D<AIDA::IHistogram1D,TH1D>::adoptRepresentation(TObject* rep)  {
    00041     TH1D* imp = dynamic_cast<TH1D*>(rep);
    00042     if ( imp )  {
    00043       if ( m_rep ) delete m_rep;
    00044       m_rep = imp;
    00045       return;
    00046     }
    00047     throw std::runtime_error("Cannot adopt native histogram representation.");
    00048   }
    00049 }
    00050
    00051 Gaudi::Histogram1D::Histogram1D()  {
    00052   m_rep = new TH1D();
    00053   init("",false);
    00054 }
    00055
    00056 Gaudi::Histogram1D::Histogram1D(TH1D* rep)  {
    00057   m_rep = rep;
    00058   init(m_rep->GetTitle());
    00059   initSums();
    00060 }
    00061
    00062 void Gaudi::Histogram1D::init(const std::string& title, bool initialize_axis)  {
    00063   m_classType = "IHistogram1D";
    00064   if ( initialize_axis )  {
    00065     m_axis.initialize(m_rep->GetXaxis(),false);
    00066   }
    00067   const TArrayD* a = m_rep->GetSumw2();
    00068   if ( 0 == a || (a && a->GetSize()==0) ) m_rep->Sumw2();
    00069   setTitle(title);
    00070   m_rep->SetDirectory(0);
    00071   m_sumEntries = 0;
    00072   m_sumwx = 0;
    00073 }
    00074
    00075 void Gaudi::Histogram1D::initSums()  {
    00076   m_sumwx = 0;
    00077   m_sumEntries = 0;
    00078   for(int i=1, n=m_rep->GetNbinsX(); i<=n; ++i)    {
    00079     m_sumwx += m_rep->GetBinContent(i)*m_rep->GetBinCenter(i);
    00080     m_sumEntries += (int)m_rep->GetBinContent(i);
    00081   }
    00082 }
    00083
    00084 bool Gaudi::Histogram1D::reset()   {
    00085   m_sumwx = 0;
    00086   m_sumEntries = 0;
    00087   return Base::reset();
    00088 }
    00089
    00091 void Gaudi::Histogram1D::adoptRepresentation(TObject* rep)  {
    00092   Gaudi::Generic1D<AIDA::IHistogram1D,TH1D>::adoptRepresentation(rep);
    00093   if ( m_rep )  {
    00094     init(m_rep->GetTitle());
    00095     initSums();
    00096   }
    00097 }
    00098
    00099 bool Gaudi::Histogram1D::setBinContents(int i,int entries ,double height,double error,double centre) {
    00100   m_rep->SetBinContent(rIndex(i),height);
    00101   m_rep->SetBinError(rIndex(i),error);
    00102   // accumulate sumwx for in range bins
    00103   if (i != AIDA::IAxis::UNDERFLOW_BIN && i !=  AIDA::IAxis::OVERFLOW_BIN )
    00104     m_sumwx  += centre*height;
    00105   m_sumEntries += entries;
    00106   return true;
    00107 }
    00108
    00109 bool Gaudi::Histogram1D::setRms(double rms) {
    00110   m_rep->SetEntries(m_sumEntries);
    00111   std::vector<double> stat(11);
    00112   // sum weights
    00113   stat[0] =  sumBinHeights();
    00114   stat[1] = 0;
    00115   if (equivalentBinEntries() != 0)
    00116     stat[1] = (  sumBinHeights() * sumBinHeights() ) / equivalentBinEntries();
    00117   stat[2] = m_sumwx;
    00118   double mean = 0;
    00119   if ( sumBinHeights() != 0 ) mean =  m_sumwx/ sumBinHeights();
    00120   stat[3] = ( mean*mean + rms*rms )* sumBinHeights();
    00121   m_rep->PutStats(&stat.front());
    00122   return true;
    00123 }
    00124
    00125 // set histogram statistics
    00126 bool Gaudi::Histogram1D::setStatistics(int allEntries,double eqBinEntries,double mean,double rms) {
    00127   m_rep->SetEntries(allEntries);
    00128   // fill statistcal vector for Root
    00129   std::vector<double> stat(11);
    00130   // sum weights
    00131   stat[0] =  sumBinHeights();
    00132   // sum weights **2
    00133   stat[1] = 0;
    00134   if (eqBinEntries != 0)
    00135     stat[1] = (  sumBinHeights() * sumBinHeights() ) / eqBinEntries;
    00136   // sum weights * x
    00137   stat[2] = mean*sumBinHeights();
    00138   // sum weight * x **2
    00139   stat[3] = ( mean*mean + rms*rms )* sumBinHeights();
    00140   m_rep->PutStats(&stat.front());
    00141   return true;
    00142 }
    00143
    00144 bool Gaudi::Histogram1D::fill ( double x,double weight )  {
    00145   (weight == 1.) ? m_rep->Fill(x) : m_rep->Fill(x,weight);
    00146   return true;
    00147 }
    00148
    00149 void Gaudi::Histogram1D::copyFromAida(const AIDA::IHistogram1D & h) {
    00150  // implement here the copy
    00151   std::string tit = h.title()+"Copy";
    00152   delete m_rep;
    00153   if (h.axis().isFixedBinning() )  {
    00154     m_rep = new TH1D(tit.c_str(),tit.c_str(),h.axis().bins(),h.axis().lowerEdge(),h.axis().upperEdge());
    00155   }
    00156   else {
    00157     Edges e;
    00158     for (int i =0; i < h.axis().bins(); ++i)  {
    00159       e.push_back(h.axis().binLowerEdge(i));
    00160     }
    00161     // add also upperedges at the end
    00162     e.push_back(h.axis().upperEdge() );
    00163     m_rep = new TH1D(tit.c_str(),tit.c_str(),e.size()-1,&e.front());
    00164   }
    00165   m_axis.initialize(m_rep->GetXaxis(),false);
    00166   m_rep->Sumw2();
    00167   m_sumEntries = 0;
    00168   m_sumwx = 0;
    00169   // sumw
    00170   double sumw = h.sumBinHeights();
    00171   // sumw2
    00172   double sumw2 = 0;
    00173   if (h.equivalentBinEntries() != 0)
    00174     sumw2 = ( sumw * sumw ) /h.equivalentBinEntries();
    00175
    00176   double sumwx = h.mean()*h.sumBinHeights();
    00177   double sumwx2 = (h.mean()*h.mean() + h.rms()*h.rms() )*h.sumBinHeights();
    00178
    00179   // copy the contents in
    00180   for (int i=-2; i < axis().bins(); ++i) {
    00181     // root binning starts from one !
    00182     m_rep->SetBinContent(rIndex(i),h.binHeight(i) );
    00183     m_rep->SetBinError(rIndex(i),h.binError(i) );
    00184   }
    00185   // need to do set entries after setting contents otherwise root will recalulate them
    00186   // taking into account how many time  SetBinContents() has been called
    00187   m_rep->SetEntries(h.allEntries());
    00188     // stat vector
    00189   std::vector<double> stat(11);
    00190   stat[0] = sumw;
    00191   stat[1] = sumw2;
    00192   stat[2] = sumwx;
    00193   stat[3] = sumwx2;
    00194   m_rep->PutStats(&stat.front());
    00195 }
    00196
    00197 StreamBuffer& Gaudi::Histogram1D::serialize(StreamBuffer& s) {
    00198   DataObject::serialize(s);
    00199   std::string title;
    00200   int size;
    00201   s >> size;
    00202   for (int j = 0; j < size; j++) {
    00203     std::string key, value;
    00204     s >> key >> value;
    00205     annotation().addItem (key, value);
    00206     if ("Title" == key) {
    00207       title = value;
    00208     }
    00209   }
    00210   double lowerEdge, upperEdge, binHeight, binError;
    00211   int    isFixedBinning, bins;
    00212   s >> isFixedBinning >> bins;
    00213
    00214   if ( m_rep ) delete m_rep;
    00215   if ( isFixedBinning ) {
    00216     s >> lowerEdge >> upperEdge;
    00217     m_rep = new TH1D(title.c_str(),title.c_str(),bins,lowerEdge,upperEdge);
    00218   } else {
    00219     Edges edges;
    00220     edges.resize(bins);
    00221     for ( int i = 0; i <= bins; ++i )
    00222       s >> *(double*)&edges[i];
    00223     m_rep = new TH1D(title.c_str(),title.c_str(),edges.size()-1,&edges.front());
    00224   }
    00225   m_axis.initialize(m_rep->GetXaxis(),true);
    00226   m_rep->Sumw2();
    00227   m_sumEntries = 0;
    00228   m_sumwx = 0;
    00229
    00230   for ( int i = 0; i <= bins + 1; ++i ) {
    00231     s >> binHeight >> binError;
    00232     m_rep->SetBinContent( i, binHeight );
    00233     m_rep->SetBinError( i, binError );
    00234   }
    00235   Stat_t allEntries;
    00236   s >> allEntries;
    00237   m_rep->SetEntries( allEntries );
    00238   Stat_t stats[4];                                    // stats array
    00239   s >> stats[0] >> stats[1] >> stats[2] >> stats[3];
    00240   m_rep->PutStats( stats );
    00241   return s;
    00242 }
    00243
    00244 StreamBuffer& Gaudi::Histogram1D::serialize(StreamBuffer& s) const {
    00245   DataObject::serialize(s);
    00246   s << static_cast<int>( annotation().size() );
    00247   for (int i = 0; i < annotation().size(); i++) {
    00248     s << annotation().key(i);
    00249     s << annotation().value(i);
    00250   }
    00251   const AIDA::IAxis & axis( this->axis() );
    00252   const int isFixedBinning = axis.isFixedBinning();
    00253   const int bins = axis.bins();
    00254   s << isFixedBinning << bins;
    00255   if ( isFixedBinning ) {
    00256     s << axis.lowerEdge();
    00257   }
    00258   else {
    00259     for ( int i = 0; i < bins; ++i )
    00260       s << axis.binLowerEdge(i);
    00261   }
    00262   s << axis.upperEdge();
    00263   for ( int i = 0; i <= bins + 1; ++i )
    00264     s << m_rep->GetBinContent(i) << m_rep->GetBinError( i );
    00265
    00266   s << m_rep->GetEntries();
    00267   Stat_t stats[4];  // stats array
    00268   m_rep->GetStats( stats );
    00269   s << stats[0] << stats[1] << stats[2] << stats[3];
    00270   return s;
    00271 }
    00272