| # Gaudi Framework, version v21r9 |

Home |  | Generated: 3 May 2010

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