| # Gaudi Framework, version v21r11 |

Home |  | Generated: 30 Sep 2010

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