| # Gaudi Framework, version v22r4 |

Home |  | Generated: Fri Sep 2 2011

# H2D.cpp

Go to the documentation of this file.

    00001 // $Id: H2D.cpp,v 1.16 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 "H2D.h"
    00014 #include "H1D.h"
    00015 #include "P1D.h"
    00016 #include "TH1D.h"
    00017 #include "TProfile.h"
    00018 #include "GaudiKernel/ObjectFactory.h"
    00019
    00020 std::pair<DataObject*,IHistogram2D*> Gaudi::createH2D(const std::string & title,int binsX,double iminX,double imaxX,int binsY,double iminY,double imaxY) {
    00021   Histogram2D* p = new Histogram2D(new TH2D(title.c_str(),title.c_str(),binsX, iminX, imaxX, binsY, iminY, imaxY));
    00022   return std::pair<DataObject*,IHistogram2D*>(p,p);
    00023 }
    00024
    00025 std::pair<DataObject*,IHistogram2D*> Gaudi::createH2D(const std::string & title, const Edges& eX, const Edges& eY) {
    00026   Histogram2D* p = new Histogram2D(new TH2D(title.c_str(),title.c_str(),eX.size()-1,&eX.front(),eY.size()-1,&eY.front()));
    00027   return std::pair<DataObject*,IHistogram2D*>(p,p);
    00028 }
    00029
    00030 std::pair<DataObject*,IHistogram2D*> Gaudi::createH2D(TH2D* rep) {
    00031   Histogram2D* p = new Histogram2D(rep);
    00032   return std::pair<DataObject*,IHistogram2D*>(p,p);
    00033 }
    00034
    00035 std::pair<DataObject*,IHistogram2D*> Gaudi::createH2D(const IHistogram2D& hist)  {
    00036   TH2D *h = getRepresentation<AIDA::IHistogram2D,TH2D>(hist);
    00037   Histogram2D *n = h ? new Histogram2D(new TH2D(*h)) : 0;
    00038   return std::pair<DataObject*,IHistogram2D*>(n,n);
    00039 }
    00040
    00041 std::pair<DataObject*,IHistogram1D*>
    00042 Gaudi::slice1DX(const std::string& nam,const IHistogram2D& hist,int first, int last)  {
    00043   TH2 *r = getRepresentation<IHistogram2D,TH2>(hist);
    00044   TH1D *t = r ? r->ProjectionX("_px",first,last,"e") : 0;
    00045   if ( t ) t->SetName(nam.c_str());
    00046   Histogram1D* p = t ? new Histogram1D(t) : 0;
    00047   return std::pair<DataObject*,IHistogram1D*>(p,p);
    00048 }
    00049
    00050 std::pair<DataObject*,IHistogram1D*>
    00051 Gaudi::slice1DY(const std::string& nam,const IHistogram2D& hist,int first, int last)  {
    00052   TH2  *r = getRepresentation<IHistogram2D,TH2>(hist);
    00053   TH1D *t = r ? r->ProjectionY("_py",first,last,"e") : 0;
    00054   if ( t ) t->SetName(nam.c_str());
    00055   Histogram1D* p = t ? new Histogram1D(t) : 0;
    00056   return std::pair<DataObject*,IHistogram1D*>(p,p);
    00057 }
    00058
    00059 std::pair<DataObject*,IHistogram1D*>
    00060 Gaudi::project1DY(const std::string& nam,const IHistogram2D& hist,int first,int last) {
    00061   TH2 *r = getRepresentation<IHistogram2D,TH2>(hist);
    00062   TH1D *t = r ? r->ProjectionY("_px",first,last,"e") : 0;
    00063   if ( t ) t->SetName(nam.c_str());
    00064   Histogram1D* p = t ? new Histogram1D(t) : 0;
    00065   return std::pair<DataObject*,IHistogram1D*>(p,p);
    00066 }
    00067
    00068 std::pair<DataObject*,IProfile1D*>
    00069 Gaudi::profile1DX(const std::string& nam,const IHistogram2D& hist,int first,int last)  {
    00070   TH2 *r = Gaudi::getRepresentation<IHistogram2D,TH2>(hist);
    00071   TProfile *t = r ? r->ProfileX("_pfx",first,last,"e") : 0;
    00072   if ( t ) t->SetName(nam.c_str());
    00073   Profile1D* p = t ? new Profile1D(t) : 0;
    00074   return std::pair<DataObject*,IProfile1D*>(p,p);
    00075 }
    00076
    00077 std::pair<DataObject*,IProfile1D*>
    00078 Gaudi::profile1DY(const std::string& nam, const IHistogram2D& hist,int first, int last) {
    00079   TH2 *r = getRepresentation<IHistogram2D,TH2>(hist);
    00080   TProfile *t = r ? r->ProfileY("_pfx",first,last,"e") : 0;
    00081   if ( t ) t->SetName(nam.c_str());
    00082   Profile1D* p = t ? new Profile1D(t) : 0;
    00083   return std::pair<DataObject*,IProfile1D*>(p,p);
    00084 }
    00085
    00086
    00087 namespace Gaudi {
    00088   template <>
    00089   void * Generic2D<IHistogram2D,TH2D>::cast(const std::string & className) const {
    00090     if (className == "AIDA::IHistogram2D")
    00091       return (IHistogram2D*)this;
    00092     else if (className == "AIDA::IHistogram")
    00093       return (IHistogram*)this;
    00094     return 0;
    00095   }
    00096
    00097   template <>
    00098   int Generic2D<IHistogram2D,TH2D>::binEntries(int indexX,int indexY) const {
    00099     if (binHeight(indexX, indexY)<=0) return 0;
    00100     double xx =  binHeight(indexX, indexY)/binError(indexX, indexY);
    00101     return int(xx*xx+0.5);
    00102   }
    00103
    00104   template <>
    00105   void Generic2D<AIDA::IHistogram2D,TH2D>::adoptRepresentation(TObject* rep)  {
    00106     TH2D* imp = dynamic_cast<TH2D*>(rep);
    00107     if ( imp )  {
    00108       if ( m_rep ) delete m_rep;
    00109       m_rep = imp;
    00110       m_xAxis.initialize(m_rep->GetXaxis(),true);
    00111       m_yAxis.initialize(m_rep->GetYaxis(),true);
    00112       const TArrayD* a = m_rep->GetSumw2();
    00113       if ( 0 == a || (a && a->GetSize()==0) ) m_rep->Sumw2();
    00114       setTitle(m_rep->GetTitle());
    00115       return;
    00116     }
    00117     throw std::runtime_error("Cannot adopt native histogram representation.");
    00118   }
    00119 }
    00120
    00121 Gaudi::Histogram2D::Histogram2D() {
    00122   m_rep = new TH2D();
    00123   m_rep->Sumw2();
    00124   m_sumEntries = 0;
    00125   m_sumwx = m_sumwy = 0;
    00126   setTitle("");
    00127   m_rep->SetDirectory(0);
    00128 }
    00129
    00130 Gaudi::Histogram2D::Histogram2D(TH2D* rep)  {
    00131   m_rep = 0;
    00132   adoptRepresentation(rep);
    00133   m_sumEntries = 0;
    00134   m_sumwx = m_sumwy = 0;
    00135   m_rep->SetDirectory(0);
    00136 }
    00137
    00138 bool Gaudi::Histogram2D::setBinContents( int i,int j,int entries,double height,double error,double centreX,double centreY) {
    00139   m_rep->SetBinContent(rIndexX(i), rIndexY(j), height);
    00140   m_rep->SetBinError(rIndexX(i), rIndexY(j), error);
    00141   // accumulate sumwx for in range bins
    00142   if (i >=0 && j >= 0) {
    00143     m_sumwx += centreX*height;
    00144     m_sumwy += centreY*height;
    00145   }
    00146   m_sumEntries += entries;
    00147   return true;
    00148 }
    00149
    00150 bool  Gaudi::Histogram2D::reset() {
    00151   m_sumwx = 0;
    00152   m_sumwy = 0;
    00153   return Base::reset();
    00154 }
    00155
    00156 #ifdef __ICC
    00157 // disable icc remark #1572: floating-point equality and inequality comparisons are unreliable
    00158 //   The comparison are meant
    00159 #pragma warning(push)
    00160 #pragma warning(disable:1572)
    00161 #endif
    00162 bool Gaudi::Histogram2D::fill ( double x,double y,double weight) {
    00163   (weight==1.) ? m_rep->Fill(x,y) : m_rep->Fill(x,y,weight );
    00164   return true;
    00165 }
    00166
    00167 bool Gaudi::Histogram2D::setRms(double rmsX,double rmsY) {
    00168   m_rep->SetEntries(m_sumEntries);
    00169   std::vector<double> stat(11);
    00170   stat[0] =  sumBinHeights();
    00171   stat[1] = 0;
    00172   if(equivalentBinEntries() != 0)
    00173     stat[1] = (sumBinHeights() * sumBinHeights()) / equivalentBinEntries();
    00174   stat[2] = m_sumwx;
    00175   double meanX = 0;
    00176   if(sumBinHeights() != 0) meanX =  m_sumwx/ sumBinHeights();
    00177   stat[3] = (meanX*meanX + rmsX*rmsX) * sumBinHeights();
    00178   stat[4] = m_sumwy;
    00179   double meanY = 0;
    00180   if(sumBinHeights() != 0) meanY =  m_sumwy/ sumBinHeights();
    00181   stat[5] = (meanY*meanY + rmsY*rmsY) * sumBinHeights();
    00182   stat[6] = 0;
    00183   m_rep->PutStats(&stat.front());
    00184   return true;
    00185 }
    00186
    00187 void Gaudi::Histogram2D::copyFromAida(const IHistogram2D& h) {
    00188   // implement here the copy
    00189   delete m_rep;
    00190   const char* tit = h.title().c_str();
    00191   if (h.xAxis().isFixedBinning() && h.yAxis().isFixedBinning()  )
    00192     m_rep = new TH2D(tit,tit,
    00193                      h.xAxis().bins(),h.xAxis().lowerEdge(),h.xAxis().upperEdge(),
    00194                      h.yAxis().bins(),h.yAxis().lowerEdge(),h.yAxis().upperEdge() );
    00195   else {
    00196     Edges eX, eY;
    00197     for (int i =0; i < h.xAxis().bins(); ++i)
    00198       eX.push_back(h.xAxis().binLowerEdge(i));
    00199     // add also upperedges at the end
    00200     eX.push_back(h.xAxis().upperEdge() );
    00201     for (int i =0; i < h.yAxis().bins(); ++i)
    00202       eY.push_back(h.yAxis().binLowerEdge(i));
    00203     // add also upperedges at the end
    00204     eY.push_back(h.yAxis().upperEdge() );
    00205     m_rep = new TH2D(tit,tit,eX.size()-1,&eX.front(),eY.size()-1,&eY.front());
    00206   }
    00207   m_xAxis.initialize(m_rep->GetXaxis(),true);
    00208   m_yAxis.initialize(m_rep->GetYaxis(),true);
    00209   m_rep->Sumw2();
    00210   m_sumEntries = 0;
    00211   m_sumwx = 0;
    00212   m_sumwy = 0;
    00213   // statistics
    00214   double sumw = h.sumBinHeights();
    00215   double sumw2 = 0;
    00216   if (h.equivalentBinEntries() != 0)
    00217     sumw2 = ( sumw * sumw ) /h.equivalentBinEntries();
    00218   double sumwx = h.meanX()*h.sumBinHeights();
    00219   double sumwx2 = (h.meanX()*h.meanX() + h.rmsX()*h.rmsX() )*h.sumBinHeights();
    00220   double sumwy = h.meanY()*h.sumBinHeights();
    00221   double sumwy2 = (h.meanY()*h.meanY() + h.rmsY()*h.rmsY() )*h.sumBinHeights();
    00222   double sumwxy = 0;
    00223
    00224   // copy the contents in  (AIDA underflow/overflow are -2,-1)
    00225   for (int i=-2; i < xAxis().bins(); ++i) {
    00226     for (int j=-2; j < yAxis().bins(); ++j) {
    00227       // root binning starts from one !
    00228       m_rep->SetBinContent(rIndexX(i), rIndexY(j), h.binHeight(i,j) );
    00229       m_rep->SetBinError(rIndexX(i), rIndexY(j), h.binError(i,j) );
    00230       // calculate statistics
    00231       if ( i >= 0 && j >= 0) {
    00232         sumwxy += h.binHeight(i,j)*h.binMeanX(i,j)*h.binMeanY(i,j);
    00233       }
    00234     }
    00235   }
    00236   // need to do set entries after setting contents otherwise root will recalulate them
    00237   // taking into account how many time  SetBinContents() has been called
    00238   m_rep->SetEntries(h.allEntries());
    00239   // fill stat vector
    00240   std::vector<double> stat(11);
    00241   stat[0] = sumw;
    00242   stat[1] = sumw2;
    00243   stat[2] = sumwx;
    00244   stat[3] = sumwx2;
    00245   stat[4] = sumwy;
    00246   stat[5] = sumwy2;
    00247   stat[6] = sumwxy;
    00248   m_rep->PutStats(&stat.front());
    00249 }
    00250
    00251 typedef Gaudi::Histogram2D H2D;
    00252 DECLARE_DATAOBJECT_FACTORY(H2D)
    00253
    00254 #ifdef __ICC
    00255 // re-enable icc remark #1572
    00256 #pragma warning(pop)
    00257 #endif

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