source: Sophya/trunk/SophyaLib/NTools/hisprof.cc@ 372

#include "machdefs.h"
#include <string.h>
#include <stdio.h>
#include <math.h>
#include "hisprof.h"
#include "perrors.h"
#include "cvector.h"

using namespace PlanckDPC;

//++
// Class        HProf
// Lib  Outils++ 
// include      hisprof.h
//
//      Classe de profil d'histogrammes.
//--

//++
// Titre        Constructeurs
//--

/********* Methode *********/
//++
HProf::HProf()
//
//      Constructeur par defaut.
//--
: Histo()
, SumY(NULL), SumY2(NULL), SumW(NULL), Ok(false), YMin(1.), YMax(-1.), Opt(0)
{
 END_CONSTRUCTOR
}

/********* Methode *********/
//++
HProf::HProf(float xMin, float xMax, int nBin, float yMin, float yMax)
//
//      Constructeur. Histogramme de profil de ``nBin'' bins entre ``xMin''
//      et ``xMax'' avec coupure d'acceptance sur y entre ``yMin'' et ``yMax''.
//      Si yMin>=yMax alors pas de coupure d'acceptance sur y.
//      Par defaut l'erreur du profil represente la dispersion dans le bin,
//      SetErrOpt(1) permet de demander de calculer l'erreur sur la moyenne.
//--
: Histo(xMin,xMax,nBin)
, SumY(new double[nBin]), SumY2(new double[nBin]), SumW(new double[nBin])
, Ok(false), YMin(yMin), YMax(yMax), Opt(0)
{
  Histo::Errors();
  Zero();
  END_CONSTRUCTOR
}

/********* Methode *********/
//++
HProf::HProf(const HProf& H)
//
//      Constructeur par copie.
//--
: Histo(H)
, SumY(new double[H.bins]), SumY2(new double[H.bins]), SumW(new double[H.bins])
, Ok(H.Ok), YMin(H.YMin), YMax(H.YMax), Opt(H.Opt)
{
  memcpy(SumY,  H.SumY,  bins*sizeof(double));
  memcpy(SumY2, H.SumY2, bins*sizeof(double));
  memcpy(SumW,  H.SumW,  bins*sizeof(double));
  END_CONSTRUCTOR
}

/********* Methode *********/
void HProf::Delete()
{
  if( SumY  != NULL ) { delete[] SumY;  SumY = NULL;}
  if( SumY2 != NULL ) { delete[] SumY2; SumY2 = NULL;}
  if( SumW  != NULL ) { delete[] SumW;  SumW = NULL;}
  Ok = false;
}

/********* Methode *********/
HProf::~HProf()
{
Delete();
}

//++
// Titre        Methodes
//--

/********* Methode *********/
//++
void HProf::Zero()
//
//      Remise a zero
//--
{
  memset(SumY,  0, bins*sizeof(double));
  memset(SumY2, 0, bins*sizeof(double));
  memset(SumW,  0, bins*sizeof(double));
  Ok = false;
  Histo::Zero();
}

/********* Methode *********/
//++
void HProf::SetErrOpt(bool spread)
//
//      Pour changer l'option de calcul de l'erreur du profile.
//      Si ``spread'' = true alors l'erreur represente la dispersion
//      des donnees dans le bin, si = false elle represente l'erreur
//      sur la moyenne du bin.
//| - Pour le bin (j):
//| H(j) = sum(y), E(j) = sum(y^2), L(j) = sum(w)
//| ->  s(j) = sqrt(E(j)/L(j) - (H(j)/L(j))^2)  dispersion
//| ->  e(j) = s(j)/sqrt(L(j)) erreur sur la moyenne
//| spread=true:  opt=0 : dispersion des donnees dans le bin = s(j)
//| spread=false: opt=1 : erreur sur la moyenne du bin = e(j)
//--
{
int opt = (spread) ? 0 : 1;
if(opt!=Opt) {Opt=opt; Ok=false;}
}

/********* Methode *********/
//++
void HProf::UpdateHisto() const
//
//      Pour mettre a jour l'histogramme de profil.
//      Il est important d'appeler cette methode si on veut
//      utiliser les methodes de la classe Histo qui ne sont
//      pas redefinies dans la classe HProf.
//      En effet, pour des raisons de precision la classe
//      HProf travaille avec des tableaux en double precision
//      et seulement au moment ou l'on a besoin de l'histogramme
//      ce dernier est rempli avec les valeurs demandees (moyenne
//      et dispersion/erreur sur la moyenne).
//--
{

float m,e2;
for(int i=0;i<bins;i++) {
  if(SumW[i]<=0.) {
    m = e2 = 0.;
  } else {
    m  = SumY[i]/SumW[i];
    e2 = SumY2[i]/SumW[i] - m*m;
    if(Opt) e2 /= SumW[i];
  }
  data[i] = m;
  err2[i] = e2;
}
// Attention, a cause de "WriteSelf const" UpdateHisto doit etre "const".
// Comme on veut modifier Ok, on est oblige de faire cette entourloupe:
HProf *buff = (HProf *) this;
buff->Ok = true;
}

/********* Methode *********/
//++
void HProf::Add(float x, float y, float w)
//
//      Addition au contenu de l'histo pour abscisse x de la valeur y et poids w
//--
{
  if(YMax>YMin && (y<YMin || YMax<y)) return;
  int numBin = FindBin(x);
  if (numBin<0) under += w;
  else if (numBin>=bins) over += w;
  else {
    Ok = false;
    SumY[numBin]  += y;
    SumY2[numBin] += y*y;
    SumW[numBin]  += w;
    nHist += w;
    nEntries++;
  }
}

/********* Methode *********/
//++
void HProf::AddBin(int numBin, float y, float w)
//
//      Addition au contenu de l'histo bin numBin de la valeur y poids w
//--
{
  if(YMax>YMin && (y<YMin || YMax<y)) return;
  if (numBin<0) under += w;
  else if (numBin>=bins) over += w;
  else {
    Ok = false;
    SumY[numBin]  += y;
    SumY2[numBin] += y*y;
    SumW[numBin]  += w;
    nHist += w;
    nEntries++;
  }
}

/********* Methode *********/
//++
HProf& HProf::operator = (const HProf& h)
//
//--
{
if(this == &h) return *this;
if( h.bins > bins ) Delete();
Histo *hthis = (Histo *) this;
*hthis = (Histo) h;
if(!SumY)  SumY  = new double[bins];
if(!SumY2) SumY2 = new double[bins];
if(!SumW)  SumW  = new double[bins];
memcpy(SumY,  h.SumY,  bins*sizeof(double));
memcpy(SumY2, h.SumY2, bins*sizeof(double));
memcpy(SumW,  h.SumW,  bins*sizeof(double));
Ok = h.Ok;
YMin = h.YMin;
YMax = h.YMax;
Opt = h.Opt;

return *this;
}

/********* Methode *********/
//++
HProf& HProf::operator += (const HProf& a)
//
//      Attention dans cette addition il n'y a pas de gestion
//      des YMin et YMax. L'addition est faite meme si les limites
//      en Y de ``a'' sont differentes de celles de ``this''.
//--
{
if(bins!=a.bins) THROW(sizeMismatchErr);
Histo *hthis = (Histo *) this;
*hthis += (Histo) a;
for(int i=0;i<bins;i++) {
  SumY[i]  += a.SumY[i];
  SumY2[i] += a.SumY2[i];
  SumW[i]  += a.SumW[i];
}
Ok = false;

return *this;
}

/********* Methode *********/
//++
HProf operator + (const HProf& a, const HProf& b)
//
//      Meme commentaire que pour l'operateur +=
//--
{
if (b.NBins()!=a.NBins()) THROW(sizeMismatchErr);
HProf c(a);
return (c += b);
}

// Rappel des inlines functions pour commentaires
//++
//   inline Histo GetHisto()
//      Retourne l'histogramme de profil.
//--
//++
// inline void GetMean(Vector& v)
//      Retourne le contenu de la moyenne dans le vecteur v
//--
//++
// inline void GetError2(Vector& v)
//      Retourne le contenu au carre de la dispersion/erreur dans le vecteur v
//--
//++
// inline float operator()(int i) const
//      Retourne le contenu du bin i
//--
//++
// inline float Error2(int i) const
//      Retourne le carre de la dispersion/erreur du bin i
//--
//++
// inline float Error(int i) const
//      Retourne la dispersion/erreur du bin i
//--
//++
// inline int Fit(GeneralFit& gfit)
//      Fit du profile par ``gfit''.
//--
//++
// inline Histo* FitResidus(GeneralFit& gfit)
//      Retourne l'Histogramme des residus par ``gfit''.
//--
//++
// inline Histo* FitFunction(GeneralFit& gfit)
//      Retourne l'Histogramme de la fonction fittee par ``gfit''.
//--
//++
// inline void Print(int dyn,float hmin,float hmax,int pflag,int il,int ih)
//      Print, voir detail dans Histo::Print
//--


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

FIO_HProf::FIO_HProf()
{
dobj=new HProf;
ownobj=true;
}

FIO_HProf::FIO_HProf(string const & filename) 
{
dobj=new HProf;
ownobj=true; 
Read(filename);
}

FIO_HProf::FIO_HProf(const HProf & obj) 
{ 
dobj = new HProf(obj);
ownobj=true; 
}

FIO_HProf::FIO_HProf(HProf * obj) 
{ 
dobj = obj;
ownobj=false; 
}

FIO_HProf::~FIO_HProf()
{
if (ownobj && dobj) delete dobj;
}

AnyDataObj* FIO_HProf::DataObj()
{
return(dobj);
}

void FIO_HProf::ReadSelf(PInPersist& is)
{
char strg[256];

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

// Lecture entete
is.GetLine(strg,255);

// Ecriture des valeurs
is.Get(dobj->bins);
is.Get(dobj->YMin);
is.Get(dobj->YMax);
is.Get(dobj->Opt);
dobj->Ok = true;

// Ecriture des donnees propres a l'histogramme de profil.
is.GetLine(strg,255);
dobj->SumY  = new double[dobj->bins];
dobj->SumY2 = new double[dobj->bins];
dobj->SumW  = new double[dobj->bins];
is.Get(dobj->SumY,  dobj->bins);
is.Get(dobj->SumY2, dobj->bins);
is.Get(dobj->SumW,  dobj->bins);

// Ecriture de l'histogramme
FIO_Histo fio_h((Histo&)*dobj);
fio_h.Read(is);

return;
}

void FIO_HProf::WriteSelf(POutPersist& os) const
{
if (dobj == NULL)   return;
char strg[256];

dobj->UpdateHisto();

// Ecriture entete pour identifier facilement
sprintf(strg,"HProf: YMin=%f  YMax=%f  Opt=%1d",dobj->YMin,dobj->YMax,dobj->Opt);
os.PutLine(strg);

// Ecriture des valeurs
os.Put(dobj->bins);
os.Put(dobj->YMin);
os.Put(dobj->YMax);
os.Put(dobj->Opt);

// Ecriture des donnees propres a l'histogramme de profil.
sprintf(strg,"HProf: SumY SumY2 SumW");
os.PutLine(strg);
os.Put(dobj->SumY,  dobj->bins);
os.Put(dobj->SumY2, dobj->bins);
os.Put(dobj->SumW,  dobj->bins);

// Ecriture de l'histogramme
FIO_Histo fio_h((Histo&)*dobj);
fio_h.Write(os);

return;
}