source: Sophya/trunk/Cosmo/RadioBeam/qhist.cc@ 3843

//  Simple specific histo-like classes 
//  R. Ansari - Avril-Mai 2010 

#include "qhist.h"

// -----------------------------------
// -- Classe ressemblant a un histo 2D 
// -----------------------------------

/* --Methode-- */
QHis1D::QHis1D()
  : nx_(0),xmin_(0.),xmax_(0.),
    nhis_(0), nunder_(0), nover_(0),
    swunder_(0.), swover_(0.)
{
}

/* --Methode-- */
QHis1D::QHis1D(QHis1D const& a)
  : nx_(a.nx_),xmin_(a.xmin_),xmax_(a.xmax_), aw_(a.aw_), 
    nhis_(a.nhis_), nunder_(a.nunder_), nover_(a.nover_),
    swunder_(a.swunder_), swover_(a.swover_)
{
}

/* --Methode-- */
QHis1D::QHis1D(r_8 xMin,r_8 xMax,int_4 nxb)
  : nx_(0),xmin_(0.),xmax_(0.),
    nhis_(0), nunder_(0), nover_(0),
    swunder_(0.), swover_(0.)
{
  Define(xMin, xMax, nxb);
}

/* --Methode-- */
void QHis1D::Define(r_8 xMin,r_8 xMax,int_4 nxb)
{
  nx_=nxb; 
  xmin_=xMin;  xmax_=xMax;  
  dxb_=(xmax_-xmin_)/(double)nx_;
  sa_size_t sz[5];  sz[0]=nx_; 
  aw_.ReSize(1,sz);
  nhis_=nunder_=nover_=0;
  swunder_=swover_=0.;
  return;
}

/* --Methode-- */
sa_size_t QHis1D::Add(r_8 x, r_8 w)
{
  sa_size_t ix = (sa_size_t)((x-xmin_)/dxb_);
  if (ix<0) {
    nunder_++;  swunder_+=w;
    return 0.;
  }
  if (ix>=nx_) {
    nover_++;  swover_+=w;
    return 0.;
  }
  nhis_++;
  aw_(ix) += w;
  return nhis_;
}

Histo QHis1D::Convert()
{
  Histo h(xmin_,xmax_,nx_);
  for(sa_size_t i=0; i<nx_; i++)  h(i) = aw_(i);
  cout << "QHis1D::Convert()/Info: Nx" << nx_ << " NUnder=" << nunder_ << " NHis=" << nhis_ 
       << " NOver=" << nover_ << "\n" 
       << " SumWUnder=" << swunder_ << " SumW=" << aw_.Sum() << " SumWOver=" << swover_ << endl;
  return h;
}

// -----------------------------------
// -- Classe ressemblant a un histo 2D 
// -----------------------------------
QHis2D::QHis2D()
  : nx(0),ny(0),xmin(0),xmax(0),ymin(0),ymax(0),sumw0(0.)
{
  ixb0 = jyb0 = 0;
}
QHis2D::QHis2D(r_8 xMin,r_8 xMax,int_4 nxb,r_8 yMin,r_8 yMax,int_4 nyb)
  : nx(0),ny(0),xmin(0),xmax(0),ymin(0),ymax(0),sumw0(0.)
{
  Define(xMin, xMax, nxb, yMin, yMax, nyb);
}
void QHis2D::Define(r_8 xMin,r_8 xMax,int_4 nxb,r_8 yMin,r_8 yMax,int_4 nyb)
{
  nx=nxb; ny=nyb;
  xmin=xMin;  xmax=xMax;  
  ymin=yMin;  ymax=yMax;  
  dxb=(xmax-xmin)/(double)nx;
  dyb=(ymax-ymin)/(double)ny;
  sa_size_t sz[5];  sz[0]=nx;  sz[1]=ny;
  aw.ReSize(2,sz);
  SetZeroBin();
  sumw0=0.;
  return;
}
double QHis2D::Add(r_8 x, r_8 y, r_8 w, bool fgfh)
{
  sa_size_t ix = (sa_size_t)((x-xmin)/dxb);
  sa_size_t jy = (sa_size_t)((y-ymin)/dyb);
  if ((ix<0)||(ix>=nx)||(jy<0)||(jy>=ny))  return 0.;
  double rw = ((ix==ixb0)&&(jy==jyb0)) ? w : 0.;
  sumw0 += rw; 
  if (fgfh) aw(ix,jy) += w;
  return rw;
}
void QHis2D::SetZeroBin(r_8 x, r_8 y)
{
  ixb0 = (sa_size_t)((x-xmin)/dxb);
  jyb0 = (sa_size_t)((y-ymin)/dyb);
}
Histo2D QHis2D::Convert()
{
  int_4 imn,jmn,imx,jmx;
  r_8 min = aw(0,0);
  r_8 max = aw(0,0);
  imn=jmn=imx=jmx=0;
  Histo2D h2(xmin,xmax,nx,ymin,ymax,ny);
  for(int_4 j=0; j<ny; j++) 
    for(int_4 i=0; i<nx; i++) {
      h2(i,j) = aw(i,j);
      if (aw(i,j)>max) {
        imx=i;  jmx=j;  max=aw(i,j);
      }
      if (aw(i,j)<min) {
        imn=i;  jmn=j;  min=aw(i,j);
      }
    }
  cout << "QHis2D::Convert()/Info: Nx,Ny=" << nx << "," << ny << " SumW=" << sumw0 
       << "\n ... Max:" << imx << "," << jmx << " ->" << max 
       << " @" << imx*dxb+xmin << "," << jmx*dyb+ymin
       << "\n ...Min:" << imn << "," << jmn << " ->" << min 
       << " @" << imn*dxb+xmin << "," << jmn*dyb+ymin << endl;
  return h2;
}