https://github.com/root-project/root
Tip revision: 1d3fcd975eb2957d69334911e994a25b7e1aa50c authored by Fons Rademakers on 05 November 2010, 14:42:46 UTC
tag patch release v5-27-06b.
tag patch release v5-27-06b.
Tip revision: 1d3fcd9
stressSpectrum.cxx
// @(#)root/test:$name: $:$id: stressSpectrum.cxx,v 1.15 2002/10/25 10:47:51 rdm exp $
// Author: Rene Brun 17/01/2006
/////////////////////////////////////////////////////////////////
//
// TSPectrum test suite
// ====================
//
// This stress program tests many elements of the TSpectrum, TSpectrum2 classes.
//
// To run in batch, do
// stressSpectrum : run 100 experiments with graphics (default)
// stressSpectrum 1000 : run 1000 experiments with graphics
// stressSpectrum -b 200 : run 200 experiments in batch mode
// stressSpectrum -b : run 100 experiments in batch mode
//
// To run interactively, do
// root -b
// Root > .x stressSpectrum.cxx : run 100 experiments with graphics (default)
// Root > .x stressSpectrum.cxx(20) : run 20 experiments
// Root > .x stressSpectrum.cxx+(30) : run 30 experiments via ACLIC
//
// Several tests are run sequentially. Each test will produce one line (Test OK or Test FAILED) .
// At the end of the test a table is printed showing the global results
// Real Time and Cpu Time.
// One single number (ROOTMARKS) is also calculated showing the relative
// performance of your machine compared to a reference machine
// a Pentium IV 3.0 Ghz) with 512 MBytes of memory
// and 120 GBytes IDE disk.
//
// An example of output when all the tests run OK is shown below:
//
//////////////////////////////////////////////////////////////////////////
// //
//****************************************************************************
//* Starting stress S P E C T R U M *
//****************************************************************************
//Peak1 : found = 70.21/ 73.75, good = 65.03/ 68.60, ghost = 8.54/ 8.39,--- OK
//Peak2 : found =163/300, good =163, ghost =8,---------------------------- OK
//****************************************************************************
//stressSpectrum: Real Time = 19.86 seconds Cpu Time = 19.04 seconds
//****************************************************************************
//* ROOTMARKS = 810.9 * Root5.09/01 20051216/1229
//****************************************************************************
#include <stdlib.h>
#include "TApplication.h"
#include "TBenchmark.h"
#include "TCanvas.h"
#include "TH2.h"
#include "TF2.h"
#include "TRandom.h"
#include "TSpectrum.h"
#include "TSpectrum2.h"
#include "TStyle.h"
#include "Riostream.h"
#include "TROOT.h"
#include "TMath.h"
Int_t npeaks;
Double_t fpeaks(Double_t *x, Double_t *par) {
Double_t result = par[0] + par[1]*x[0];
for (Int_t p=0;p<npeaks;p++) {
Double_t norm = par[3*p+2];
Double_t mean = par[3*p+3];
Double_t sigma = par[3*p+4];
result += norm*TMath::Gaus(x[0],mean,sigma);
}
return result;
}
Double_t fpeaks2(Double_t *x, Double_t *par) {
Double_t result = 0.1;
for (Int_t p=0;p<npeaks;p++) {
Double_t norm = par[5*p+0];
Double_t mean1 = par[5*p+1];
Double_t sigma1 = par[5*p+2];
Double_t mean2 = par[5*p+3];
Double_t sigma2 = par[5*p+4];
result += norm*TMath::Gaus(x[0],mean1,sigma1)*TMath::Gaus(x[1],mean2,sigma2);
}
return result;
}
void findPeaks(Int_t pmin, Int_t pmax, Int_t &nfound, Int_t &ngood, Int_t &nghost) {
npeaks = (Int_t)gRandom->Uniform(pmin,pmax);
Int_t nbins = 500;
Double_t dxbins = 2;
TH1F *h = new TH1F("h","test",nbins,0,nbins*dxbins);
//generate n peaks at random
Double_t par[3000];
par[0] = 0.8;
par[1] = -0.6/1000;
Int_t p,pf;
for (p=0;p<npeaks;p++) {
par[3*p+2] = 1;
par[3*p+3] = 10+gRandom->Rndm()*(nbins-20)*dxbins;
par[3*p+4] = 3+2*gRandom->Rndm();
}
TF1 *f = new TF1("f",fpeaks,0,nbins*dxbins,2+3*npeaks);
f->SetNpx(1000);
f->SetParameters(par);
h->FillRandom("f",200000);
TSpectrum *s = new TSpectrum(4*npeaks);
nfound = s->Search(h,2,"goff");
//Search found peaks
ngood = 0;
Float_t *xpeaks = s->GetPositionX();
for (p=0;p<npeaks;p++) {
for (Int_t pf=0;pf<nfound;pf++) {
Double_t dx = TMath::Abs(xpeaks[pf] - par[3*p+3]);
if (dx <dxbins) ngood++;
}
}
//Search ghost peaks
nghost = 0;
for (pf=0;pf<nfound;pf++) {
Int_t nf=0;
for (Int_t p=0;p<npeaks;p++) {
Double_t dx = TMath::Abs(xpeaks[pf] - par[3*p+3]);
if (dx <dxbins) nf++;
}
if (nf == 0) nghost++;
}
delete f;
delete h;
delete s;
}
void stress1(Int_t ntimes) {
Int_t pmin = 5;
Int_t pmax = 55;
TCanvas *c1 = new TCanvas("c1","Spectrum results",10,10,800,800);
c1->Divide(2,2);
gStyle->SetOptFit();
TH1F *hpeaks = new TH1F("hpeaks","Number of peaks",pmax-pmin,pmin,pmax);
TH1F *hfound = new TH1F("hfound","% peak founds",100,0,100);
TH1F *hgood = new TH1F("hgood", "% good peaks",100,0,100);
TH1F *hghost = new TH1F("hghost","% ghost peaks",100,0,100);
Int_t nfound,ngood,nghost;
for (Int_t i=0;i<ntimes;i++) {
findPeaks(pmin,pmax,nfound,ngood,nghost);
hpeaks->Fill(npeaks);
hfound->Fill(100*Double_t(nfound)/Double_t(npeaks));
hgood->Fill(100*Double_t(ngood)/Double_t(npeaks));
hghost->Fill(100*Double_t(nghost)/Double_t(npeaks));
//printf("npeaks = %d, nfound = %d, ngood = %d, nghost = %d\n",npeaks,nfound,ngood,nghost);
}
c1->cd(1);
hpeaks->Fit("pol1","lq");
c1->cd(2);
hfound->Fit("gaus","lq");
c1->cd(3);
hgood->Fit("gaus","lq");
c1->cd(4);
hghost->Fit("gaus","lq","",0,30);
c1->cd();
Double_t p1 = hfound->GetFunction("gaus")->GetParameter(1);
Double_t ep1 = hfound->GetFunction("gaus")->GetParError(1);
Double_t p2 = hgood->GetFunction("gaus")->GetParameter(1);
Double_t ep2 = hgood->GetFunction("gaus")->GetParError(1);
Double_t p3 = hghost->GetFunction("gaus")->GetParameter(1);
Double_t ep3 = hghost->GetFunction("gaus")->GetParError(1);
Double_t p1ref = 70.21; //ref numbers obtained with ntimes=1000
Double_t p2ref = 65.03;
Double_t p3ref = 8.54;
//printf("p1=%g+-%g, p2=%g+-%g, p3=%g+-%g\n",p1,ep1,p2,ep2,p3,ep3);
char sok[20];
if (TMath::Abs(p1ref-p1) < 2*ep1 && TMath::Abs(p2ref-p2) < 2*ep2 && TMath::Abs(p3ref-p3) < 2*ep3 ) {
sprintf(sok,"OK");
} else {
sprintf(sok,"failed");
}
printf("Peak1 : found =%6.2f/%6.2f, good =%6.2f/%6.2f, ghost =%5.2f/%5.2f,--- %s\n",
p1,p1ref,p2,p2ref,p3,p3ref,sok);
}
void stress2(Int_t np2) {
npeaks = np2;
TRandom r;
Int_t nbinsx = 200;
Int_t nbinsy = 200;
Double_t xmin = 0;
Double_t xmax = (Double_t)nbinsx;
Double_t ymin = 0;
Double_t ymax = (Double_t)nbinsy;
Double_t dx = (xmax-xmin)/nbinsx;
Double_t dy = (ymax-ymin)/nbinsy;
TH2F *h2 = new TH2F("h2","test",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
h2->SetStats(0);
//generate n peaks at random
Double_t par[3000];
Int_t p;
for (p=0;p<npeaks;p++) {
par[5*p+0] = r.Uniform(0.2,1);
par[5*p+1] = r.Uniform(xmin,xmax);
par[5*p+2] = r.Uniform(dx,5*dx);
par[5*p+3] = r.Uniform(ymin,ymax);
par[5*p+4] = r.Uniform(dy,5*dy);
}
TF2 *f2 = new TF2("f2",fpeaks2,xmin,xmax,ymin,ymax,5*npeaks);
f2->SetNpx(100);
f2->SetNpy(100);
f2->SetParameters(par);
h2->FillRandom("f2",500000);
//now the real stuff
TSpectrum2 *s = new TSpectrum2(2*npeaks);
Int_t nfound = s->Search(h2,2,"goff noMarkov");
//searching good and ghost peaks (approximation)
Int_t pf,ngood = 0;
Float_t *xpeaks = s->GetPositionX();
Float_t *ypeaks = s->GetPositionY();
for (p=0;p<npeaks;p++) {
for (Int_t pf=0;pf<nfound;pf++) {
Double_t diffx = TMath::Abs(xpeaks[pf] - par[5*p+1]);
Double_t diffy = TMath::Abs(ypeaks[pf] - par[5*p+3]);
if (diffx < 2*dx && diffy < 2*dy) ngood++;
}
}
if (ngood > nfound) ngood = nfound;
//Search ghost peaks (approximation)
Int_t nghost = 0;
for (pf=0;pf<nfound;pf++) {
Int_t nf=0;
for (Int_t p=0;p<npeaks;p++) {
Double_t diffx = TMath::Abs(xpeaks[pf] - par[5*p+1]);
Double_t diffy = TMath::Abs(ypeaks[pf] - par[5*p+3]);
if (diffx < 2*dx && diffy < 2*dy) nf++;
}
if (nf == 0) nghost++;
}
delete s;
delete f2;
delete h2;
Int_t nfoundRef = 163;
Int_t ngoodRef = 163;
Int_t nghostRef = 8;
char sok[20];
if ( TMath::Abs(nfound - nfoundRef) < 5
&& TMath::Abs(ngood - ngoodRef) < 5
&& TMath::Abs(nghost - nghostRef) < 5) {
sprintf(sok,"OK");
} else {
sprintf(sok,"failed");
}
printf("Peak2 : found =%d/%d, good =%d, ghost =%2d,---------------------------- %s\n",
nfound,npeaks,ngood,nghost,sok);
}
#ifndef __CINT__
void stressSpectrum(Int_t ntimes) {
#else
void stressSpectrum(Int_t ntimes=100) {
#endif
cout << "****************************************************************************" <<endl;
cout << "* Starting stress S P E C T R U M *" <<endl;
cout << "****************************************************************************" <<endl;
gBenchmark->Start("stressSpectrum");
stress1(ntimes);
stress2(300);
gBenchmark->Stop ("stressSpectrum");
Double_t reftime100 = 19.04; //pcbrun compiled
Double_t ct = gBenchmark->GetCpuTime("stressSpectrum");
const Double_t rootmarks = 800*reftime100*ntimes/(100*ct);
printf("****************************************************************************\n");
gBenchmark->Print("stressSpectrum");
printf("****************************************************************************\n");
printf("* ROOTMARKS =%6.1f * Root%-8s %d/%d\n",rootmarks,gROOT->GetVersion(),
gROOT->GetVersionDate(),gROOT->GetVersionTime());
printf("****************************************************************************\n");
}
#ifndef __CINT__
int main(int argc, char **argv)
{
TApplication theApp("App", &argc, argv);
gROOT->SetBatch();
gBenchmark = new TBenchmark();
Int_t ntimes = 100;
if (argc > 1) ntimes = atoi(argv[1]);
stressSpectrum(ntimes);
return 0;
}
#endif