source: Sophya/trunk/SophyaLib/HiStats/histerr.cc@ 2628

#include "sopnamsp.h"
#include "machdefs.h"
#include <string.h>
#include <stdio.h>
#include <math.h>
#include "histinit.h"
#include "histerr.h"
#include "perrors.h"

/*!
  \class SOPHYA::HistoErr
  \ingroup HiStats
  Classe d'histogrammes 1D avec erreurs donnees par l'utilisateur
*/

/********* Methode *********/
/*! Constructeur par defaut */
HistoErr::HistoErr(void)
: Histo(), mNData(NULL)
{
}

/********* Methode *********/
/*! Constructeur d'un histo de nBin bins allant de xMin a xMax */
HistoErr::HistoErr(r_8 xMin, r_8 xMax, int_4 nBin)
: Histo(xMin,xMax,nBin), mNData(NULL)
{
 this->Errors();
 allocate_mNData(nBin);
 Zero();
}

/********* Methode *********/
/*! Constructeur par copie */
HistoErr::HistoErr(const HistoErr& H)
: Histo(H), mNData(NULL)
{
 allocate_mNData(H.mBins);
 if(mBins>0) memcpy(mNData,H.mNData,mBins*sizeof(r_8));
}

/********* Methode *********/
/*! Destructeur */
HistoErr::~HistoErr(void)
{
 Delete();
}

/********* Methode *********/
/*! Allocation du tableau mNData */
void HistoErr::allocate_mNData(int nbin)
{
 if(nbin<=0) return;
 if(mNData) {delete [] mNData; mNData=NULL;}
 mNData = new r_8[nbin];
}

/********* Methode *********/
/*! Delete des tableaux */
void HistoErr::Delete(void)
{
 Histo::Delete();
 if(mNData) {delete [] mNData; mNData=NULL;}
}

/********* Methode *********/
/*!
  Remise a zero
*/
void HistoErr::Zero(void)
{
 Histo::Zero();
 if(mNData) memset(mNData,0,mBins*sizeof(r_8));
}

/********* Methode *********/
/*!
 Addition du contenu de l'histo pour abscisse x poids w et l'erreur e
*/
void HistoErr::Add(r_8 x, r_8 w, r_8 e)
{
 int_4 numBin = (int_4)floor((x-mMin)/binWidth);
 if(numBin<0) mUnder += w;
 else if(numBin>=mBins) mOver += w;
 else {
   mData[numBin] += w; mNData[numBin] += 1.; mErr2[numBin] += e*e;
   nHist += w; nEntries++;
 }
}

/********* Methode *********/
/*!
 Addition du contenu de l'histo pour le bin numBin poids w et l'erreur e
*/
void HistoErr::AddBin(int_4 numBin, r_8 w, r_8 e)
{
 if(numBin<0) mUnder += w;
 else if(numBin>=mBins) mOver += w;
 else {
   mData[numBin] += w; mNData[numBin] += 1.; mErr2[numBin] += e*e;
   nHist += w; nEntries++;
 }
}

/********* Methode *********/
/*!
 Remplissage du contenu de l'histo pour le bin numBin poids w et l'erreur e
*/
void HistoErr::SetBin(int_4 numBin, r_8 w, r_8 e, r_8 nb)
{
Histo::SetBin(numBin,w);
Histo::SetErr2(numBin,e*e);
SetNentB(numBin,nb);
return;
}

/********* Methode *********/
/*!
 Remplissage nombre d'entrees pour le bin numBin
*/
void HistoErr::SetNentB(int_4 numBin, r_8 nb)
{
 if(numBin>=0 && numBin<mBins) mNData[numBin] = nb;
return;
}

/*!
  Remplissage d'un tableau avec les nombres d'entrees dans le bin de l'histo
*/
void HistoErr::GetNBin(TVector<r_8>& v) const
{
v.Realloc(mBins);
for(int_4 i=0;i<mBins;i++) v(i) = mNData[i];
return;
}

/********* Methode *********/
/*!
  Remplissage du nombre d'entrees dans les bins de l'histo avec les valeurs d'un vecteur
*/
void HistoErr::PutNBin(TVector<r_8> &v)
{
int_4 n = (v.NElts()<mBins) ? v.NElts(): mBins;
if(n>0) for(int_4 i=0;i<n;i++) mNData[i] = v(i);
return;
}

/********* Methode *********/
/*!
 Recompute XMin and XMax so that
 the CENTER of the first bin is exactly XMin and
 the CENTER of the last bin is exactly XMax.
 Remember that otherwise
 XMin is the beginning of the first bin
 and XMax is the end of the last bin
*/
void HistoErr::ReCenterBin(void)
{
 if(mBins<=1) return;
 double dx = (mMax-mMin)/(mBins-1);
 mMin -= dx/2.;
 mMax += dx/2.;
 binWidth = (mMax-mMin)/mBins;
}

/********* Methode *********/
/*!
  Compute the correlation histogram.
  Each bin content is divided by the number of entries in that bin.
  Each squared error is divided by the number of entries in that bin.
  The number of entries by bin is NOT set to 1 (calling ToCorrel
    many time will change the histogram !)
*/
void HistoErr::ToCorrel(void)
{
 if(mBins<1) return;
 for(int_4 i=0;i<mBins;i++) {
   if(mNData[i]<1.) continue;
   mData[i] /= mNData[i];
   mErr2[i] /= mNData[i];
 }
 return;
}

/********* Methode *********/
/*!
 Fill the histogram with an other histogram
*/
void HistoErr::FillFrHErr(HistoErr& hfrom)
{
 if(mBins<=0) return;
 if(hfrom.mBins<=0) return;

 Zero();

 for(int_4 i=0;i<hfrom.mBins;i++) {
   int numBin = FindBin(hfrom.BinCenter(i));
   if(numBin<0) mUnder += hfrom.mData[i];
   else if(numBin>=mBins) mOver += hfrom.mData[i];
   else {
     mData[numBin]  += hfrom.mData[i];
     mNData[numBin] += hfrom.mNData[i];
     mErr2[numBin]  += hfrom.mErr2[i];
     nHist          += hfrom.mData[i];
     nEntries++;
   }
 }

}

/********* Methode *********/
/*!
  Operateur egal HistoErr = HistoErr
*/
HistoErr& HistoErr::operator = (const HistoErr& h)
{
  if(this==&h) return *this;
  Delete();
  if(h.mBins<=0) return *this;

  // Copy the "Histo" part
  (Histo)(*this) = Histo::operator=(h);
  // Copy the "entries by bin" table
  allocate_mNData(h.mBins);
  memcpy(mNData,h.mNData,mBins*sizeof(r_8));

  return *this;
}

///////////////////////////////////////////////////////////
// --------------------------------------------------------
//   Les objets delegues pour la gestion de persistance 
// --------------------------------------------------------
///////////////////////////////////////////////////////////

DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
void ObjFileIO<HistoErr>::ReadSelf(PInPersist& is)
{
string strg;

if(dobj==NULL) dobj=new HistoErr;
  else         dobj->Delete();

// Lecture entete
is.GetStr(strg);

// Lecture des donnees HistoErr
is.Get(dobj->mBins);
  dobj->allocate_mNData(dobj->mBins);
is.Get(dobj->mNData, dobj->mBins);

// Lecture de l'histogramme
is >> (Histo&)(*dobj);

return;
}

DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
void ObjFileIO<HistoErr>::WriteSelf(POutPersist& os) const
{
if(dobj == NULL)   return;
string strg;

strg = "HistErr";
os.PutStr(strg);

// Ecriture des valeurs
os.Put(dobj->mBins);

// Ecriture des donnees HistoErr nombre d entree par bin
os.Put(dobj->mNData, dobj->mBins);

// Ecriture de l'histogramme
os << (Histo&)(*dobj);

return;
}

#ifdef __CXX_PRAGMA_TEMPLATES__
#pragma define_template ObjFileIO<HistoErr>
#endif

#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
template class ObjFileIO<HistoErr>;
#endif