| # Gaudi Framework, version v21r6 |

Home |  | Generated: 11 Nov 2009

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