Hi Rene,
as I promised you during ACAT in Chicago, I send you and everybody
new TH1K histogram class.
//*-*
//*-* TH1K class supports the nearest K Neighbors 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 Neighbors method, usually >=3
//*-* K = 0, means default, where K is selected by TH1K in such a way
//*-* that DistanceToNearestKthNeighbour > BinWidth and K >=3
//*-*
//*-*
I hope it will be useful, especially in low statistics case.
It was tested in STAR.
All comments and proposals from everybody are welcome
Victor
--
Victor M. Perevoztchikov perev@bnl.gov perev@vxcern.cern.ch
Brookhaven National Laboratory MS 510A PO Box 5000 Upton NY 11973-5000
tel office : 631-344-7894; fax 631-344-4206; home 631-345-2690
// @(#)root/hist:$Name: $:$Id: TH1K.cxx,v 1.1 2001/02/10 02:01:09 perev Exp $
// Author: Rene Brun 26/12/94
/*************************************************************************
* 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 <string.h>
#include <stdio.h>
#include <ctype.h>
#include <fstream.h>
#include <iostream.h>
#include <float.h>
#include <assert.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 claster analisys.
//*-* 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
//*-*
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
ClassImp(TH1K)
//______________________________________________________________________________
// TH1K methods
//______________________________________________________________________________
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()
{
}
//______________________________________________________________________________
void TH1K::Reset(Option_t *option)
{
fNIn =0; fReady = 0;
TH1::Reset(option);
}
//______________________________________________________________________________
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 ) ? fabs(fArray[jl]-x) : 1.e+20;
fr = (jr<fNIn) ? fabs(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);
}
//______________________________________________________________________________
static int 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),&fcompare);
}
// @(#)root/hist:$Name: $:$Id: TH1K.h,v 1.1 2001/02/10 02:00:42 perev Exp $
// Author: Rene Brun 26/12/94
/*************************************************************************
* 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. *
*************************************************************************/
#ifndef ROOT_TH1K
#define ROOT_TH1K
//////////////////////////////////////////////////////////////////////////
// //
// TH1K //
// //
// 1-Dim histogram nearest K Neighbour class. //
// //
//////////////////////////////////////////////////////////////////////////
#include "TH1.h"
//________________________________________________________________________
class TH1K : public TH1, public TArrayF {
public:
TH1K();
TH1K(const char *name,const char *title,Int_t nbinsx,Axis_t xlow,Axis_t xup,Int_t k=0);
virtual ~TH1K();
virtual Stat_t GetBinContent(Int_t bin) const;
virtual Stat_t GetBinContent(Int_t bin,Int_t) const
{return GetBinContent(bin);}
virtual Stat_t GetBinContent(Int_t bin,Int_t,Int_t) const
{return GetBinContent(bin);}
virtual Stat_t GetBinError(Int_t bin) const;
virtual Stat_t GetBinError(Int_t bin,Int_t) const
{return GetBinError(bin);}
virtual Stat_t GetBinError(Int_t bin,Int_t,Int_t) const
{return GetBinError(bin);}
virtual void Reset(Option_t *option="");
virtual Int_t Fill(Axis_t x);
virtual Int_t Fill(Axis_t x,Stat_t w){return TH1::Fill(x,w);}
void SetKOrd(Int_t k){fKOrd=k;}
protected:
Int_t fReady; //!
Int_t fNIn;
Int_t fKOrd; //!
Int_t fKCur; //!
private:
void Sort();
ClassDef(TH1K,1) //1-Dim Nearest Kth neighbour method
};
#endif
void padRefresh(TPad *pad,int flag=0);
void hksimple()
{
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
//*-*
//*-* This program creates :
//*-* - a one dimensional histogram
//*-* - a two dimensional histogram
//*-* - a profile histogram
//*-* - a memory-resident ntuple
//*-*
//*-* These objects are filled with some random numbers and saved on a file.
//*-*
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
gROOT->Reset();
TH1 *hpx[3];
int j;
// Create a new canvas.
c1 = new TCanvas("c1","Dynamic Filling Example",200,10,700,500);
c1->SetFillColor(42);
c1->GetFrame()->SetFillColor(21);
c1->GetFrame()->SetBorderSize(6);
c1->GetFrame()->SetBorderMode(-1);
// Create some histograms, a profile histogram and an ntuple
hpx[0] = new TH1F("hp0","This is the px distribution",1000,-4,4);
hpx[1] = new TH1K("hk1","This is the px distribution",1000,-4,4);
hpx[2] = new TH1K("hk2","This is the px distribution",1000,-4,4,16);
c1->Divide(1,3);
for (j=0;j<3;j++) {c1->cd(j+1); hpx[j]->Draw();hpx[j]->SetFillColor(48);}
// Set canvas/frame attributes (save old attributes)
gBenchmark->Start("hksimple");
// Fill histograms randomly
gRandom->SetSeed();
Float_t px, py, pz;
const Int_t kUPDATE = 10;
for (Int_t i = 0; i <= 300; i++) {
gRandom->Rannor(px,py);
for (j=0;j<3;j++) {hpx[j]->Fill(px);}
if (i && (i%kUPDATE) == 0) {
padRefresh(c1);
}
}
gBenchmark->Show("hksimple");
}
void padRefresh(TPad *pad,int flag)
{
if (!pad) return;
pad->Modified();
pad->Update();
TList *tl = pad->GetListOfPrimitives();
if (!tl) return;
TListIter next(tl);
TObject *to;
while ((to=next())) {
if (to->InheritsFrom(TPad::Class())) padRefresh((TPad*)to,1);}
if (flag) return;
gSystem->ProcessEvents();
}
This archive was generated by hypermail 2b29 : Tue Jan 01 2002 - 17:50:36 MET