source: Sophya/trunk/SophyaLib/HiStats/hisprof.cc@ 1089

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

/*!
  \class SOPHYA::HProf
  \ingroup HiStats
  Classe de profil d'histogrammes.
*/

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

/********* Methode *********/
/*!
  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.
*/
HProf::HProf(float xMin, float xMax, int nBin, float yMin, float yMax)
: 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 *********/
/*!
  Constructeur par copie.
*/
HProf::HProf(const HProf& H)
: Histo(H)
, SumY((H.bins>0)  ? new double[H.bins] : NULL)
, SumY2((H.bins>0) ? new double[H.bins] : NULL)
, SumW((H.bins>0)  ? new double[H.bins] : NULL)
, Ok(H.Ok), YMin(H.YMin), YMax(H.YMax), Opt(H.Opt)
{
  if(bins>0) {
    memcpy(SumY,  H.SumY,  bins*sizeof(double));
    memcpy(SumY2, H.SumY2, bins*sizeof(double));
    memcpy(SumW,  H.SumW,  bins*sizeof(double));
  }
  UpdateHisto();
  END_CONSTRUCTOR
}

/********* Methode *********/
/*!
  Des-allocation
*/
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 *********/
/*!
  Destructeur
*/
HProf::~HProf()
{
Delete();
}

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

/********* Methode *********/
/*!
  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.
  \verbatim
    - 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)
  \endverbatim
*/
void HProf::SetErrOpt(bool spread)
{
int opt = (spread) ? 0 : 1;
if(opt!=Opt) {Opt=opt; Ok=false;}
}

/********* Methode *********/
/*!
  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).
*/
void HProf::updatehisto() const
{
float m,e2;
if(bins>0) 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;
}
Ok = true;
// 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 *********/
/*!
  Addition au contenu de l'histo pour abscisse x de la valeur y et poids w
*/
void HProf::Add(float x, float y, float 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 *********/
/*!
  Addition au contenu de l'histo bin numBin de la valeur y poids w
*/
void HProf::AddBin(int numBin, float y, float 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 *********/
/*!
  Operateur HProf H = H1
*/
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 *********/
/*!
  Operateur HProf H += H1

  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''.
*/
HProf& HProf::operator += (const HProf& a)
{
if(bins!=a.bins) THROW(sizeMismatchErr);
Histo *hthis = (Histo *) this;
*hthis += (Histo) a;
if(bins>0) for(int i=0;i<bins;i++) {
  SumY[i]  += a.SumY[i];
  SumY2[i] += a.SumY2[i];
  SumW[i]  += a.SumW[i];
}
updatehisto();

return *this;
}

/********* Methode *********/
/*!
  Pour rebinner l'histogramme de profile sur nbinew bins
*/
void HProf::HRebin(int nbinew)
{
  if( bins <= 0 ) return; // createur par default
  if( nbinew <= 0 ) return;

  // memorisation de l'HProf original
  HProf H(*this);

  // Rebin de la partie Histo
  int binold = bins;
  Histo::HRebin(nbinew);

  // Le nombre de bins est il plus grand ??
  if( bins > binold ) {
    delete [] SumY;  SumY  = new double[nbinew]; memset(SumY, 0,bins*sizeof(double));
    delete [] SumY2; SumY2 = new double[nbinew]; memset(SumY2,0,bins*sizeof(double));
    delete [] SumW;  SumW  = new double[nbinew]; memset(SumW, 0,bins*sizeof(double));
  }

  // Remplissage des parties propres au HPprof
  for(int i=0;i<bins;i++) {
    float xmi = BinLowEdge(i);
    float xma = BinHighEdge(i);
    int imi =  H.FindBin(xmi);
    if( imi < 0 ) imi = 0;
    int ima =  H.FindBin(xma);
    if( ima >= H.bins ) ima = H.bins-1;
    double wY = 0., wY2 = 0., wW = 0.;
    if( ima == imi ) {
      wY  = H.SumY[imi]  * binWidth/H.binWidth;
      wY2 = H.SumY2[imi] * binWidth/H.binWidth;
      wW  = H.SumW[imi]  * binWidth/H.binWidth;
    } else {
      wY  += H.SumY[imi]  * (H.BinHighEdge(imi)-xmi)/H.binWidth;
      wY  += H.SumY[ima]  * (xma -H.BinLowEdge(ima))/H.binWidth;
      wY2 += H.SumY2[imi] * (H.BinHighEdge(imi)-xmi)/H.binWidth;
      wY2 += H.SumY2[ima] * (xma -H.BinLowEdge(ima))/H.binWidth;
      wW  += H.SumW[imi]  * (H.BinHighEdge(imi)-xmi)/H.binWidth;
      wW  += H.SumW[ima]  * (xma -H.BinLowEdge(ima))/H.binWidth;
      if(ima>imi+1) for(int ii=imi+1;ii<ima;ii++)
        {wY += H.SumY[ii]; wY2 += H.SumY2[ii]; wW += H.SumW[ii];}
    }
    SumY[i] += wY; SumY2[i] += wY2; SumW[i] += wW;
  }
  // On synchronise les tableaux Sum?? et l'Histogramme
  updatehisto();
}

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


void ObjFileIO<HProf>::ReadSelf(PInPersist& is)
{
char strg[256];

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

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

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

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

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

void ObjFileIO<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);
os << (Histo&)(*dobj);

return;
}

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

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