// @(#)root/hist:$Name:  $:$Id: TH1K.cxx,v 1.5 2002/01/24 11:39:29 rdm Exp $
// Author: Victor Perevoztchikov <perev@bnl.gov>  21/02/2001

/*************************************************************************
 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

#include <stdlib.h>

#include "Riostream.h"
#include "TROOT.h"
#include "TH1K.h"
#include "TMath.h"

//*-**-*-*-*-*-*-*-*-*-*-*The T H 1 K   Class*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*--*-*-*-*
//*-*                  ===============================
//*-*
//*-*  TH1K class supports the nearest K Neighbours method,
//*-*       widely used in cluster analysis.
//*-*       This method is especially useful for small statistics.
//*-*
//*-*       In this method :
//*-*         DensityOfProbability ~ 1/DistanceToNearestKthNeighbour
//*-*
//*-*      Ctr TH1K::TH1K(name,title,nbins,xlow,xup,K=0)
//*-*      differs from TH1F only by "K"
//*-*      K - is the order of K Neighbours method, usually >=3
//*-*      K = 0, means default, where K is selected by TH1K in such a way
//*-*             that DistanceToNearestKthNeighbour > BinWidth and K >=3
//*-*
//*-*  This class has been implemented by Victor Perevoztchikov <perev@bnl.gov>
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

ClassImp(TH1K)

//______________________________________________________________________________
 TH1K::TH1K(): TH1(), TArrayF()
{
   fDimension = 1;
   fNIn   = 0;
   fReady = 0;
   fKOrd  = 3;
}

//______________________________________________________________________________
 TH1K::TH1K(const char *name,const char *title,Int_t nbins,Axis_t xlow,Axis_t xup,Int_t k)
     : TH1(name,title,nbins,xlow,xup), TArrayF(100)
{
//
//    Create a 1-Dim histogram with fix bins of type float
//    ====================================================
//           (see TH1K::TH1 for explanation of parameters)
//
   fDimension = 1;
   fNIn = 0;
   fReady = 0;
   fKOrd = k;
}

//______________________________________________________________________________
 TH1K::~TH1K()
{
}

//______________________________________________________________________________
 Int_t TH1K::Fill(Axis_t x)
{
//*-*-*-*-*-*-*-*-*-*Increment bin with abscissa X by 1*-*-*-*-*-*-*-*-*-*-*
//*-*                ==================================
//*-*
//*-* if x is less than the low-edge of the first bin, the Underflow bin is incremented
//*-* if x is greater than the upper edge of last bin, the Overflow bin is incremented
//*-*
//*-* If the storage of the sum of squares of weights has been triggered,
//*-* via the function Sumw2, then the sum of the squares of weights is incremented
//*-* by 1 in the bin corresponding to x.
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

   fReady = 0;
   Int_t bin;
   fEntries++;
   bin =fXaxis.FindBin(x);
   if (bin == 0 || bin > fXaxis.GetNbins()) return -1;
   ++fTsumw;
   ++fTsumw2;
   fTsumwx  += x;
   fTsumwx2 += x*x;
   fReady = 0;
   if (fNIn == fN) Set(fN*2);
   AddAt(x,fNIn++);
   return bin;
}

//______________________________________________________________________________
 Stat_t TH1K::GetBinContent(Int_t bin) const
{
//*-*-*-*-*-*-*Return content of global bin number bin*-*-*-*-*-*-*-*-*-*-*-*
//*-*          =======================================
   if (!fReady) {
      ((TH1K*)this)->Sort();
      ((TH1K*)this)->fReady=1;
   }
   if (!fNIn)   return 0.;
   float x = GetBinCenter(bin);
   int left = TMath::BinarySearch(fNIn,fArray,x);
   int jl=left,jr=left+1,nk,nkmax =fKOrd;
//   assert(jl==-1   || fArray[jl] < x);
//   assert(jr==fNIn || fArray[jr] > x);
   float fl,fr,ff=0.,ffmin=1.e-6;
   if (!nkmax) {nkmax = 3; ffmin = GetBinWidth(bin);}
   if (nkmax >= fNIn) nkmax = fNIn-1;
   for (nk = 1; nk <= nkmax || ff <= ffmin; nk++) {
     fl = (jl>=0  ) ? TMath::Abs(fArray[jl]-x) : 1.e+20;
     fr = (jr<fNIn) ? TMath::Abs(fArray[jr]-x) : 1.e+20;
     if (jl<0 && jr>=fNIn) break;
     if (fl < fr) { ff = fl; jl--;}
     else         { ff = fr; jr++;}
   }
   ((TH1K*)this)->fKCur = nk - 1;
   return fNIn * 0.5*fKCur/((float)(fNIn+1))*GetBinWidth(bin)/ff;
}

//______________________________________________________________________________
 Stat_t TH1K::GetBinError(Int_t bin) const
{
//*-*-*-*-*-*-*Return content of global bin error *-*-*-*-*-*-*-*-*-*-*-*
//*-*          ===================================
   return TMath::Sqrt(((double)(fNIn-fKCur+1))/((fNIn+1)*(fKCur-1)))*GetBinContent(bin);
}

//______________________________________________________________________________
 void   TH1K::Reset(Option_t *option)
{
  fNIn   =0;
  fReady = 0;
  TH1::Reset(option);
}

//______________________________________________________________________________
 void TH1K::SavePrimitive(ofstream &out, Option_t *option)
{
    // Save primitive as a C++ statement(s) on output stream out

   //Note the following restrictions in the code generated:
   // - variable bin size not implemented
   // - Objects in list of functions not saved (fits)
   // - Contours not saved

   char quote = '"';
   out<<"   "<<endl;
   out<<"   "<<"TH1 *";

   out<<GetName()<<" = new "<<ClassName()<<"("<<quote<<GetName()<<quote<<","<<quote<<GetTitle()<<quote
                 <<","<<GetXaxis()->GetNbins()
                 <<","<<GetXaxis()->GetXmin()
                 <<","<<GetXaxis()->GetXmax()
                 <<","<<fKOrd;

              out<<");"<<endl;
   if (fDirectory == 0) {
      out<<"   "<<GetName()<<"->SetDirectory(0);"<<endl;
   }
   if (TestBit(kNoStats)) {
      out<<"   "<<GetName()<<"->SetStats(0);"<<endl;
   }
   if (fOption.Length() != 0) {
      out<<"   "<<GetName()<<"->SetOption("<<quote<<fOption.Data()<<quote<<");"<<endl;
   }

   if (fNIn) {
      out<<"   Float_t Arr[]={"<<endl;
      for (int i=0; i<fNIn; i++) {
         out<< fArray[i];
         if (i != fNIn-1) {out<< ",";} else { out<< "};";}
         if (i%10 == 9)   {out<< endl;}
      }
      out<< endl;
      out<<"   for(int i=0;i<" << fNIn << ";i++)"<<GetName()<<"->Fill(Arr[i]);";
      out<< endl;
   }
   SaveFillAttributes(out,GetName(),0,1001);
   SaveLineAttributes(out,GetName(),1,1,1);
   SaveMarkerAttributes(out,GetName(),1,1,1);
   fXaxis.SaveAttributes(out,GetName(),"->GetXaxis()");
   fYaxis.SaveAttributes(out,GetName(),"->GetYaxis()");
   fZaxis.SaveAttributes(out,GetName(),"->GetZaxis()");
   TString opt = option;
   opt.ToLower();
   if (!opt.Contains("nodraw")) {
      out<<"   "<<GetName()<<"->Draw("
      <<quote<<option<<quote<<");"<<endl;
   }
}


//______________________________________________________________________________
static int TH1K_fcompare(const void *f1,const void *f2)
{
  if (*((float*)f1) < *((float*)f2)) return -1;
  if (*((float*)f1) > *((float*)f2)) return  1;
  return 0;
}

//______________________________________________________________________________
 void TH1K::Sort()
{
  if (fNIn<2) return;
  qsort(GetArray(),fNIn,sizeof(Float_t),&TH1K_fcompare);
}