| # Gaudi Framework, version v23r1 |

Home |  | Generated: Wed Feb 29 2012

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

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Defines