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

#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), mCorrel(0)
{
}

/********* 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(), mNData(NULL)
{
  CreateOrResize(xMin,xMax,nBin);
}

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

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

/********* Methode *********/
/*! Gestion de l'allocation */
void HistoErr::CreateOrResize(r_8 xMin, r_8 xMax, int_4 nBin)
{
  //cout<<"HistoErr::CreateOrResize()"<<endl;
  bool samelen = (nBin==mBins)? true: false;
  Histo::CreateOrResize(xMin,xMax,nBin);
  Histo::Errors();
  if(mNData!=NULL && !samelen) {delete [] mNData; mNData=NULL;}
  if(nBin>0 && mNData==NULL) mNData = new r_8[nBin];
  if(mNData) memset(mNData,0,nBin*sizeof(r_8));
  mCorrel = 0;
}

/********* Methode *********/
/*! Delete des tableaux */
void HistoErr::Delete(void)
{
  //cout<<"HistoErr::Delete()"<<endl;
  Histo::Delete();
  if(mNData) {delete [] mNData; mNData=NULL;}
  mCorrel = 0;
}

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

/********* 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 *********/
/*!
  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;
 mCorrel++;
 for(int_4 i=0;i<mBins;i++) {
   if(mNData[i]<1.) continue;
   mData[i] /= mNData[i];
   mErr2[i] /= mNData[i];
 }
 return;
}

/********* Methode *********/
/*!
 Recompute back the original HistoErr before ToCorrel action
*/
void HistoErr::FromCorrel(void)
{
 if(mBins<1) return;
 mCorrel--;
 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++;
   }
 }
 mCorrel = hfrom.mCorrel;

}

/********* Methode *********/
/*!
  Operateur egal HistoErr = HistoErr
*/
HistoErr& HistoErr::operator = (const HistoErr& h)
{
  if(this==&h) return *this;
  CreateOrResize(h.mMin,h.mMax,h.mBins);
  if(mData) memcpy(mData,h.mData,mBins*sizeof(r_8));
  if(mErr2) memcpy(mErr2,h.mErr2,mBins*sizeof(r_8));
  mUnder = h.mUnder; mOver = h.mOver;
  nHist = h.nHist; nEntries = h.nEntries;
  if(mNData) memcpy(mNData,h.mNData,mBins*sizeof(r_8));
  mCorrel = h.mCorrel;

  return *this;
}

/********* Methode *********/
/*!
  Print info
*/
void HistoErr::Show(ostream & os) const
{
  Histo::Show(os);
  os <<"       ncorrel="<<mCorrel<< endl;
}

///////////////////////////////////////////////////////////
// --------------------------------------------------------
//   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);

// Nombre d'appels a ToCorrel/FromCorrel
is.Get(dobj->mCorrel);

// Lecture des donnees HistoErr
is.Get(dobj->mBins);
if(dobj->mBins>0) dobj->mNData = new r_8[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;

// Ecriture entete
strg = "HistErr";
os.PutStr(strg);

// Nombre d'appels a ToCorrel/FromCorrel
os.Put(dobj->mCorrel);

// 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 SOPHYA::ObjFileIO<HistoErr>;
#endif