ROOT logo
// @(#)root/hist:$Id: TH1K.cxx 20882 2007-11-19 11:31:26Z rdm $
// 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 "TH1K.h"
#include "TMath.h"

//______________________________________________________________________________
// The TH1K 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()
{
   // Constructor.

   fDimension = 1;
   fNIn   = 0;
   fReady = 0;
   fKOrd  = 3;
}


//______________________________________________________________________________
TH1K::TH1K(const char *name,const char *title,Int_t nbins,Double_t xlow,Double_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()
{
   // Destructor.
}


//______________________________________________________________________________
Int_t TH1K::Fill(Double_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()) {
      if (!fgStatOverflows) return -1;
   }
   ++fTsumw;
   ++fTsumw2;
   fTsumwx  += x;
   fTsumwx2 += x*x;
   fReady = 0;
   if (fNIn == fN) Set(fN*2);
   AddAt(x,fNIn++);
   return bin;
}


//______________________________________________________________________________
Double_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;
   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;
}


//______________________________________________________________________________
Double_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)
{
   // Reset.

   fNIn   =0;
   fReady = 0;
   TH1::Reset(option);
}


//______________________________________________________________________________
void TH1K::SavePrimitive(ostream &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)
{
   // Compare.

   if (*((float*)f1) < *((float*)f2)) return -1;
   if (*((float*)f1) > *((float*)f2)) return  1;
   return 0;
}


//______________________________________________________________________________
void TH1K::Sort()
{
   // Sort.

   if (fNIn<2) return;
   qsort(GetArray(),fNIn,sizeof(Float_t),&TH1K_fcompare);
}
 TH1K.cxx:1
 TH1K.cxx:2
 TH1K.cxx:3
 TH1K.cxx:4
 TH1K.cxx:5
 TH1K.cxx:6
 TH1K.cxx:7
 TH1K.cxx:8
 TH1K.cxx:9
 TH1K.cxx:10
 TH1K.cxx:11
 TH1K.cxx:12
 TH1K.cxx:13
 TH1K.cxx:14
 TH1K.cxx:15
 TH1K.cxx:16
 TH1K.cxx:17
 TH1K.cxx:18
 TH1K.cxx:19
 TH1K.cxx:20
 TH1K.cxx:21
 TH1K.cxx:22
 TH1K.cxx:23
 TH1K.cxx:24
 TH1K.cxx:25
 TH1K.cxx:26
 TH1K.cxx:27
 TH1K.cxx:28
 TH1K.cxx:29
 TH1K.cxx:30
 TH1K.cxx:31
 TH1K.cxx:32
 TH1K.cxx:33
 TH1K.cxx:34
 TH1K.cxx:35
 TH1K.cxx:36
 TH1K.cxx:37
 TH1K.cxx:38
 TH1K.cxx:39
 TH1K.cxx:40
 TH1K.cxx:41
 TH1K.cxx:42
 TH1K.cxx:43
 TH1K.cxx:44
 TH1K.cxx:45
 TH1K.cxx:46
 TH1K.cxx:47
 TH1K.cxx:48
 TH1K.cxx:49
 TH1K.cxx:50
 TH1K.cxx:51
 TH1K.cxx:52
 TH1K.cxx:53
 TH1K.cxx:54
 TH1K.cxx:55
 TH1K.cxx:56
 TH1K.cxx:57
 TH1K.cxx:58
 TH1K.cxx:59
 TH1K.cxx:60
 TH1K.cxx:61
 TH1K.cxx:62
 TH1K.cxx:63
 TH1K.cxx:64
 TH1K.cxx:65
 TH1K.cxx:66
 TH1K.cxx:67
 TH1K.cxx:68
 TH1K.cxx:69
 TH1K.cxx:70
 TH1K.cxx:71
 TH1K.cxx:72
 TH1K.cxx:73
 TH1K.cxx:74
 TH1K.cxx:75
 TH1K.cxx:76
 TH1K.cxx:77
 TH1K.cxx:78
 TH1K.cxx:79
 TH1K.cxx:80
 TH1K.cxx:81
 TH1K.cxx:82
 TH1K.cxx:83
 TH1K.cxx:84
 TH1K.cxx:85
 TH1K.cxx:86
 TH1K.cxx:87
 TH1K.cxx:88
 TH1K.cxx:89
 TH1K.cxx:90
 TH1K.cxx:91
 TH1K.cxx:92
 TH1K.cxx:93
 TH1K.cxx:94
 TH1K.cxx:95
 TH1K.cxx:96
 TH1K.cxx:97
 TH1K.cxx:98
 TH1K.cxx:99
 TH1K.cxx:100
 TH1K.cxx:101
 TH1K.cxx:102
 TH1K.cxx:103
 TH1K.cxx:104
 TH1K.cxx:105
 TH1K.cxx:106
 TH1K.cxx:107
 TH1K.cxx:108
 TH1K.cxx:109
 TH1K.cxx:110
 TH1K.cxx:111
 TH1K.cxx:112
 TH1K.cxx:113
 TH1K.cxx:114
 TH1K.cxx:115
 TH1K.cxx:116
 TH1K.cxx:117
 TH1K.cxx:118
 TH1K.cxx:119
 TH1K.cxx:120
 TH1K.cxx:121
 TH1K.cxx:122
 TH1K.cxx:123
 TH1K.cxx:124
 TH1K.cxx:125
 TH1K.cxx:126
 TH1K.cxx:127
 TH1K.cxx:128
 TH1K.cxx:129
 TH1K.cxx:130
 TH1K.cxx:131
 TH1K.cxx:132
 TH1K.cxx:133
 TH1K.cxx:134
 TH1K.cxx:135
 TH1K.cxx:136
 TH1K.cxx:137
 TH1K.cxx:138
 TH1K.cxx:139
 TH1K.cxx:140
 TH1K.cxx:141
 TH1K.cxx:142
 TH1K.cxx:143
 TH1K.cxx:144
 TH1K.cxx:145
 TH1K.cxx:146
 TH1K.cxx:147
 TH1K.cxx:148
 TH1K.cxx:149
 TH1K.cxx:150
 TH1K.cxx:151
 TH1K.cxx:152
 TH1K.cxx:153
 TH1K.cxx:154
 TH1K.cxx:155
 TH1K.cxx:156
 TH1K.cxx:157
 TH1K.cxx:158
 TH1K.cxx:159
 TH1K.cxx:160
 TH1K.cxx:161
 TH1K.cxx:162
 TH1K.cxx:163
 TH1K.cxx:164
 TH1K.cxx:165
 TH1K.cxx:166
 TH1K.cxx:167
 TH1K.cxx:168
 TH1K.cxx:169
 TH1K.cxx:170
 TH1K.cxx:171
 TH1K.cxx:172
 TH1K.cxx:173
 TH1K.cxx:174
 TH1K.cxx:175
 TH1K.cxx:176
 TH1K.cxx:177
 TH1K.cxx:178
 TH1K.cxx:179
 TH1K.cxx:180
 TH1K.cxx:181
 TH1K.cxx:182
 TH1K.cxx:183
 TH1K.cxx:184
 TH1K.cxx:185
 TH1K.cxx:186
 TH1K.cxx:187
 TH1K.cxx:188
 TH1K.cxx:189
 TH1K.cxx:190
 TH1K.cxx:191
 TH1K.cxx:192
 TH1K.cxx:193
 TH1K.cxx:194
 TH1K.cxx:195
 TH1K.cxx:196
 TH1K.cxx:197
 TH1K.cxx:198
 TH1K.cxx:199
 TH1K.cxx:200
 TH1K.cxx:201
 TH1K.cxx:202
 TH1K.cxx:203
 TH1K.cxx:204
 TH1K.cxx:205
 TH1K.cxx:206
 TH1K.cxx:207
 TH1K.cxx:208
 TH1K.cxx:209
 TH1K.cxx:210
 TH1K.cxx:211
 TH1K.cxx:212
 TH1K.cxx:213
 TH1K.cxx:214
 TH1K.cxx:215
 TH1K.cxx:216
 TH1K.cxx:217
 TH1K.cxx:218
 TH1K.cxx:219
 TH1K.cxx:220
 TH1K.cxx:221
 TH1K.cxx:222
 TH1K.cxx:223
 TH1K.cxx:224