ROOT   Reference Guide
Loading...
Searching...
No Matches
peaks2.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_spectrum
3/// \notebook
4/// Example to illustrate the 2-d peak finder (class TSpectrum2).
5///
6/// This script generates a random number of 2-d gaussian peaks
7/// The position of the peaks is found via TSpectrum2
8/// To execute this example, do:
9///
10/// ~~~{.cpp}
11/// root > .x peaks2.C (generate up to 50 peaks by default)
12/// root > .x peaks2.C(10) (generate up to 10 peaks)
13/// root > .x peaks2.C+(200) (generate up to 200 peaks via ACLIC)
14/// ~~~
15///
16/// The script will iterate generating a new histogram having
17/// between 5 and the maximun number of peaks specified.
18/// Double Click on the bottom right corner of the pad to go to a new spectrum
19/// To Quit, select the "quit" item in the canvas "File" menu
20///
21/// \macro_image
22/// \macro_code
23///
24/// \author Rene Brun
25
26#include "TSpectrum2.h"
27#include "TCanvas.h"
28#include "TRandom.h"
29#include "TH2.h"
30#include "TF2.h"
31#include "TMath.h"
32#include "TROOT.h"
33
34TSpectrum2 *s;
35TH2F *h2 = nullptr;
36Int_t npeaks = 30;
37Double_t fpeaks2(Double_t *x, Double_t *par) {
38 Double_t result = 0.1;
39 for (Int_t p=0;p<npeaks;p++) {
40 Double_t norm = par[5*p+0];
41 Double_t mean1 = par[5*p+1];
42 Double_t sigma1 = par[5*p+2];
43 Double_t mean2 = par[5*p+3];
44 Double_t sigma2 = par[5*p+4];
45 result += norm*TMath::Gaus(x[0],mean1,sigma1)*TMath::Gaus(x[1],mean2,sigma2);
46 }
47 return result;
48}
49void findPeak2() {
50 printf("Generating histogram with %d peaks\n",npeaks);
51 Int_t nbinsx = 200;
52 Int_t nbinsy = 200;
53 Double_t xmin = 0;
54 Double_t xmax = (Double_t)nbinsx;
55 Double_t ymin = 0;
56 Double_t ymax = (Double_t)nbinsy;
57 Double_t dx = (xmax-xmin)/nbinsx;
58 Double_t dy = (ymax-ymin)/nbinsy;
59 delete h2;
60 h2 = new TH2F("h2","test",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
61 h2->SetStats(false);
62 //generate n peaks at random
63 Double_t par[3000];
64 Int_t p;
65 for (p=0;p<npeaks;p++) {
66 par[5*p+0] = gRandom->Uniform(0.2,1);
67 par[5*p+1] = gRandom->Uniform(xmin,xmax);
68 par[5*p+2] = gRandom->Uniform(dx,5*dx);
69 par[5*p+3] = gRandom->Uniform(ymin,ymax);
70 par[5*p+4] = gRandom->Uniform(dy,5*dy);
71 }
72 TF2 *f2 = new TF2("f2",fpeaks2,xmin,xmax,ymin,ymax,5*npeaks);
73 f2->SetNpx(100);
74 f2->SetNpy(100);
75 f2->SetParameters(par);
76 TCanvas *c1 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("c1");
77 if (!c1) c1 = new TCanvas("c1","c1",10,10,1000,700);
78 h2->FillRandom("f2",500000);
79
80 //now the real stuff: Finding the peaks
81 Int_t nfound = s->Search(h2,2,"col");
82
83 //searching good and ghost peaks (approximation)
84 Int_t pf,ngood = 0;
85 Double_t *xpeaks = s->GetPositionX();
86 Double_t *ypeaks = s->GetPositionY();
87 for (p=0;p<npeaks;p++) {
88 for (pf=0;pf<nfound;pf++) {
89 Double_t diffx = TMath::Abs(xpeaks[pf] - par[5*p+1]);
90 Double_t diffy = TMath::Abs(ypeaks[pf] - par[5*p+3]);
91 if (diffx < 2*dx && diffy < 2*dy) ngood++;
92 }
93 }
94 if (ngood > nfound) ngood = nfound;
95 //Search ghost peaks (approximation)
96 Int_t nghost = 0;
97 for (pf=0;pf<nfound;pf++) {
98 Int_t nf=0;
99 for (p=0;p<npeaks;p++) {
100 Double_t diffx = TMath::Abs(xpeaks[pf] - par[5*p+1]);
101 Double_t diffy = TMath::Abs(ypeaks[pf] - par[5*p+3]);
102 if (diffx < 2*dx && diffy < 2*dy) nf++;
103 }
104 if (nf == 0) nghost++;
105 }
106 c1->Update();
107
108 s->Print();
109 printf("Gener=%d, Found=%d, Good=%d, Ghost=%d\n",npeaks,nfound,ngood,nghost);
110 if (!gROOT->IsBatch()) {
111 printf("\nDouble click in the bottom right corner of the pad to continue\n");
112 c1->WaitPrimitive();
113 }
114}
115void peaks2(Int_t maxpeaks=50) {
116 s = new TSpectrum2(2*maxpeaks);
117 for (int i=0; i<10; ++i) {
118 npeaks = (Int_t)gRandom->Uniform(5,maxpeaks);
119 findPeak2();
120 }
121}
122
123
int Int_t
Definition RtypesCore.h:45
double Double_t
Definition RtypesCore.h:59
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
float xmin
float ymin
float xmax
float ymax
#define gROOT
Definition TROOT.h:407
R__EXTERN TRandom * gRandom
Definition TRandom.h:62
The Canvas class.
Definition TCanvas.h:23
virtual void SetNpx(Int_t npx=100)
Set the number of points used to draw the function.
Definition TF1.cxx:3433
virtual void SetParameters(const Double_t *params)
Definition TF1.h:670
A 2-Dim function with parameters.
Definition TF2.h:29
virtual void SetNpy(Int_t npy=100)
Set the number of points used to draw the function.
Definition TF2.cxx:927
2-D histogram with a float per channel (see TH1 documentation)
Definition TH2.h:307
TObject * FindObject(const char *name) const override
Search if object named name is inside this pad or in pads inside this pad.
Definition TPad.cxx:2647
virtual Double_t Uniform(Double_t x1=1)
Returns a uniform deviate on the interval (0, x1).
Definition TRandom.cxx:682
Advanced 2-dimensional spectra processing.
Definition TSpectrum2.h:18
Double_t * GetPositionY() const
Definition TSpectrum2.h:45
Double_t * GetPositionX() const
Definition TSpectrum2.h:44
void Print(Option_t *option="") const override
Print the array of positions.
virtual Int_t Search(const TH1 *hist, Double_t sigma=2, Option_t *option="", Double_t threshold=0.05)
This function searches for peaks in source spectrum in hin The number of found peaks and their positi...
return c1
Definition legend1.C:41
Double_t x[n]
Definition legend1.C:17
Double_t Gaus(Double_t x, Double_t mean=0, Double_t sigma=1, Bool_t norm=kFALSE)
Calculates a gaussian function with mean and sigma.
Definition TMath.cxx:471
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:123