| # Gaudi Framework, version v22r2 |

Home |  | Generated: Tue May 10 2011

# H3D.cpp

Go to the documentation of this file.

    00001 #ifdef __ICC
    00002 // disable icc remark #2259: non-pointer conversion from "X" to "Y" may lose significant bits
    00003 //   TODO: To be removed, since it comes from ROOT TMathBase.h
    00004 #pragma warning(disable:2259)
    00005 #endif
    00006 #ifdef WIN32
    00007 // Disable warning
    00008 //   warning C4996: 'sprintf': This function or variable may be unsafe.
    00009 // coming from TString.h
    00010 #pragma warning(disable:4996)
    00011 #endif
    00012
    00013 #include "GaudiKernel/DataObject.h"
    00014 #include "GaudiKernel/ObjectFactory.h"
    00015 #include "GaudiPI.h"
    00016 #include "Generic3D.h"
    00017 #include "TH3D.h"
    00018
    00019 namespace Gaudi {
    00020
    00027   class GAUDI_API Histogram3D : public DataObject, public Generic3D<AIDA::IHistogram3D,TH3D>   {
    00028   public:
    00030     Histogram3D();
    00032     Histogram3D(TH3D* rep);
    00034     virtual ~Histogram3D() {}
    00036     bool fill ( double x, double y, double z, double weight);
    00038     virtual bool setBinContents( int i, int j, int k, int entries,double height,double error,double centreX, double centreY, double centreZ );
    00040     virtual bool setRms(double rmsX, double rmsY, double rmsZ);
    00041     // overwrite reset
    00042     virtual bool reset();
    00044     void* cast(const std::string & className) const;
    00046     void  copyFromAida(const AIDA::IHistogram3D & h);
    00048     virtual const CLID& clID() const { return Gaudi::Histogram3D::classID(); }
    00049     static const CLID& classID()     { return CLID_H3D; }
    00050
    00051   protected:
    00052     // cache sumwx and sumwy  when setting contents since I don't have bin mean
    00053     double m_sumwx;
    00054     double m_sumwy;
    00055     double m_sumwz;
    00056   };
    00057 }
    00058
    00059 namespace Gaudi {
    00060   template <>
    00061   void Generic3D<AIDA::IHistogram3D,TH3D>::adoptRepresentation(TObject* rep)  {
    00062     TH3D* imp = dynamic_cast<TH3D*>(rep);
    00063     if ( imp )  {
    00064       if ( m_rep ) delete m_rep;
    00065       m_rep = imp;
    00066       m_xAxis.initialize(m_rep->GetXaxis(),true);
    00067       m_yAxis.initialize(m_rep->GetYaxis(),true);
    00068       m_zAxis.initialize(m_rep->GetZaxis(),true);
    00069       const TArrayD* a = m_rep->GetSumw2();
    00070       if ( 0 == a || (a && a->GetSize()==0) ) m_rep->Sumw2();
    00071       setTitle(m_rep->GetTitle());
    00072       return;
    00073     }
    00074     throw std::runtime_error("Cannot adopt native histogram representation.");
    00075   }
    00076 }
    00077
    00079 std::pair<DataObject*,AIDA::IHistogram3D*>
    00080 Gaudi::createH3D(const std::string& title,int nBinsX,double xlow,double xup,int nBinsY,double ylow,double yup,int nBinsZ,double zlow,double zup) {
    00081   Histogram3D* p = new Histogram3D(new TH3D(title.c_str(),title.c_str(),nBinsX,xlow,xup,nBinsY,ylow,yup,nBinsZ,zlow,zup));
    00082   return std::pair<DataObject*,AIDA::IHistogram3D*>(p,p);
    00083 }
    00084
    00086 std::pair<DataObject*,AIDA::IHistogram3D*>
    00087 Gaudi::createH3D(const std::string& title, const Edges& eX,const Edges& eY,const Edges& eZ) {
    00088   Histogram3D* p = new Histogram3D(new TH3D(title.c_str(),title.c_str(),eX.size()-1,&eX.front(), eY.size()-1,&eY.front(), eZ.size()-1,&eZ.front()));
    00089   return std::pair<DataObject*,AIDA::IHistogram3D*>(p,p);
    00090 }
    00091
    00092 std::pair<DataObject*,AIDA::IHistogram3D*> Gaudi::createH3D(const AIDA::IHistogram3D& hist)  {
    00093   TH3D *h = getRepresentation<AIDA::IHistogram3D,TH3D>(hist);
    00094   Histogram3D *n = h ? new Histogram3D(new TH3D(*h)) : 0;
    00095   return std::pair<DataObject*,AIDA::IHistogram3D*>(n,n);
    00096 }
    00097
    00098 Gaudi::Histogram3D::Histogram3D() {
    00099   m_rep = new TH3D();
    00100   setTitle("");
    00101   m_rep->Sumw2();
    00102   m_sumEntries = 0;
    00103   m_sumwx = 0;
    00104   m_sumwy = 0;
    00105   m_sumwz = 0;
    00106   m_rep->SetDirectory(0);
    00107 }
    00108
    00109 Gaudi::Histogram3D::Histogram3D(TH3D* rep) {
    00110   m_rep = 0;
    00111   adoptRepresentation(rep);
    00112   m_sumEntries = 0;
    00113   m_sumwx = 0;
    00114   m_sumwy = 0;
    00115   m_sumwz = 0;
    00116   m_rep->SetDirectory(0);
    00117 }
    00118
    00119 // set bin content (entries and centre are not used )
    00120 bool Gaudi::Histogram3D::setBinContents(int i,int j,int k,int entries,double height,double error,double centreX, double centreY, double centreZ ) {
    00121   m_rep->SetBinContent(rIndexX(i), rIndexY(j), rIndexZ(k), height);
    00122   m_rep->SetBinError(rIndexX(i), rIndexY(j), rIndexZ(k), error);
    00123   // accumulate sum bin centers
    00124   if (i >=0 && j >= 0 && k >= 0) {
    00125     m_sumwx += centreX*height;
    00126     m_sumwy += centreY*height;
    00127     m_sumwz += centreZ*height;
    00128   }
    00129   m_sumEntries += entries;
    00130   return true;
    00131 }
    00132
    00133 bool Gaudi::Histogram3D::reset (  )   {
    00134   m_sumwx = 0;
    00135   m_sumwy = 0;
    00136   m_sumwz = 0;
    00137   m_sumEntries = 0;
    00138   m_rep->Reset ( );
    00139   return true;
    00140 }
    00141
    00142 bool Gaudi::Histogram3D::fill ( double x, double y, double z, double weight)  {
    00143   m_rep->Fill ( x , y, z, weight );
    00144   return true;
    00145 }
    00146
    00147 void* Gaudi::Histogram3D::cast(const std::string & className) const   {
    00148   if (className == "AIDA::IHistogram3D")   {
    00149     return (AIDA::IHistogram3D*)this;
    00150   }
    00151   else if (className == "AIDA::IHistogram")   {
    00152     return (AIDA::IHistogram*)this;
    00153   }
    00154   return 0;
    00155 }
    00156
    00157 #ifdef __ICC
    00158 // disable icc remark #1572: floating-point equality and inequality comparisons are unreliable
    00159 //   The comparison are meant
    00160 #pragma warning(push)
    00161 #pragma warning(disable:1572)
    00162 #endif
    00163 bool Gaudi::Histogram3D::setRms(double rmsX, double rmsY, double rmsZ   ) {
    00164   m_rep->SetEntries(m_sumEntries);
    00165   std::vector<double> stat(11);
    00166   // sum weights
    00167   stat[0] =  sumBinHeights();
    00168   stat[1] = 0;
    00169   if (equivalentBinEntries() != 0)
    00170     stat[1] = (  sumBinHeights() * sumBinHeights() ) / equivalentBinEntries();
    00171   stat[2] = m_sumwx;
    00172   double meanX = 0;
    00173   if ( sumBinHeights() != 0 ) meanX =  m_sumwx/ sumBinHeights();
    00174   stat[3] = ( meanX*meanX + rmsX*rmsX )* sumBinHeights();
    00175   stat[4] = m_sumwy;
    00176   double meanY = 0;
    00177   if ( sumBinHeights() != 0 ) meanY =  m_sumwy/ sumBinHeights();
    00178   stat[5] = ( meanY*meanY + rmsY*rmsY )* sumBinHeights();
    00179   stat[6] = 0;
    00180   stat[7] = m_sumwz;
    00181   double meanZ = 0;
    00182   if ( sumBinHeights() != 0 ) meanZ =  m_sumwz/ sumBinHeights();
    00183   stat[8] = ( meanZ*meanZ + rmsZ*rmsZ )* sumBinHeights();
    00184   // do not need to use sumwxy sumwxz and sumwyz
    00185
    00186   m_rep->PutStats(&stat.front());
    00187   return true;
    00188 }
    00189
    00190 void Gaudi::Histogram3D::copyFromAida(const AIDA::IHistogram3D & h) {
    00191   delete m_rep;
    00192   // implement here the copy
    00193   const char* tit = h.title().c_str();
    00194   if (h.xAxis().isFixedBinning() && h.yAxis().isFixedBinning() &&  h.zAxis().isFixedBinning() )  {
    00195     m_rep = new TH3D(tit,tit,
    00196       h.xAxis().bins(), h.xAxis().lowerEdge(), h.xAxis().upperEdge(),
    00197       h.yAxis().bins(), h.yAxis().lowerEdge(), h.yAxis().upperEdge(),
    00198       h.zAxis().bins(), h.zAxis().lowerEdge(), h.zAxis().upperEdge() );
    00199   }
    00200   else {
    00201     Edges eX, eY, eZ;
    00202     for (int i =0; i < h.xAxis().bins(); ++i)
    00203       eX.push_back(h.xAxis().binLowerEdge(i));
    00204     // add also upperedges at the end
    00205     eX.push_back(h.xAxis().upperEdge() );
    00206     for (int i =0; i < h.yAxis().bins(); ++i)
    00207       eY.push_back(h.yAxis().binLowerEdge(i));
    00208     // add also upperedges at the end
    00209     eY.push_back(h.yAxis().upperEdge() );
    00210     for (int i =0; i < h.zAxis().bins(); ++i)
    00211       eZ.push_back(h.zAxis().binLowerEdge(i));
    00212     // add also upperedges at the end
    00213     eZ.push_back(h.zAxis().upperEdge() );
    00214     m_rep = new TH3D(tit,tit,eX.size()-1,&eX.front(),eY.size()-1,&eY.front(),eZ.size()-1,&eZ.front());
    00215   }
    00216   m_xAxis.initialize(m_rep->GetXaxis(),true);
    00217   m_yAxis.initialize(m_rep->GetYaxis(),true);
    00218   m_zAxis.initialize(m_rep->GetZaxis(),true);
    00219   const TArrayD* a = m_rep->GetSumw2();
    00220   if ( 0 == a || (a && a->GetSize()==0) ) m_rep->Sumw2();
    00221   m_sumEntries = 0;
    00222   m_sumwx = 0;
    00223   m_sumwy = 0;
    00224   m_sumwz = 0;
    00225
    00226   // statistics
    00227   double sumw   = h.sumBinHeights();
    00228   double sumw2  = 0;
    00229   if (h.equivalentBinEntries() != 0)
    00230     sumw2 = ( sumw * sumw ) /h.equivalentBinEntries();
    00231   double sumwx  = h.meanX()*h.sumBinHeights();
    00232   double sumwx2 = (h.meanX()*h.meanX() + h.rmsX()*h.rmsX() )*h.sumBinHeights();
    00233   double sumwy  = h.meanY()*h.sumBinHeights();
    00234   double sumwy2 = (h.meanY()*h.meanY() + h.rmsY()*h.rmsY() )*h.sumBinHeights();
    00235   double sumwz  = h.meanZ()*h.sumBinHeights();
    00236   double sumwz2 = (h.meanZ()*h.meanZ() + h.rmsZ()*h.rmsZ() )*h.sumBinHeights();
    00237   double sumwxy = 0;
    00238   double sumwxz = 0;
    00239   double sumwyz = 0;
    00240
    00241   // copy the contents in  (AIDA underflow/overflow are -2,-1)
    00242   for (int i=-2; i < xAxis().bins(); ++i) {
    00243     for (int j=-2; j < yAxis().bins(); ++j) {
    00244       for (int k=-2; k < zAxis().bins(); ++k) {
    00245         m_rep->SetBinContent(rIndexX(i), rIndexY(j), rIndexZ(k), h.binHeight(i,j,k) );
    00246         m_rep->SetBinError(rIndexX(i), rIndexY(j), rIndexZ(k), h.binError(i,j,k) );
    00247         // calculate statistics
    00248         if ( i >= 0 && j >= 0 && k >= 0) {
    00249           sumwxy += h.binHeight(i,j,k)*h.binMeanX(i,j,k)*h.binMeanY(i,j,k);
    00250           sumwxz += h.binHeight(i,j,k)*h.binMeanX(i,j,k)*h.binMeanZ(i,j,k);
    00251           sumwyz += h.binHeight(i,j,k)*h.binMeanY(i,j,k)*h.binMeanZ(i,j,k);
    00252         }
    00253       }
    00254     }
    00255   }
    00256   // need to do set entries after setting contents otherwise root will recalulate them
    00257   // taking into account how many time  SetBinContents() has been called
    00258   m_rep->SetEntries(h.allEntries());
    00259
    00260   // fill stat vector
    00261   std::vector<double> stat(11);
    00262   stat[0] = sumw;
    00263   stat[1] = sumw2;
    00264   stat[2] = sumwx;
    00265   stat[3] = sumwx2;
    00266   stat[4] = sumwy;
    00267   stat[5] = sumwy2;
    00268   stat[6] = sumwxy;
    00269   stat[7] = sumwz;
    00270   stat[8] = sumwz2;
    00271   stat[9] = sumwxz;
    00272   stat[10] = sumwyz;
    00273   m_rep->PutStats(&stat.front());
    00274 }
    00275
    00276 #ifdef __ICC
    00277 // re-enable icc remark #1572
    00278 #pragma warning(pop)
    00279 #endif
    00280
    00281 typedef Gaudi::Histogram3D H3D;
    00282 DECLARE_DATAOBJECT_FACTORY(H3D)

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