StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
FindPeaks.C
1 // This script generates a random number of gaussian peaks
2 // on top of a linear background.
3 // The position of the peaks is found via TSpectrum and injected
4 // as initial values of parameters to make a global fit.
5 // The background is computed and drawn on top of the original histogram.
6 //
7 // To execute this example, do
8 // root > .x peaks.C (generate 10 peaks by default)
9 // root > .x peaks.C++ (use the compiler)
10 // root > .x peaks.C++(30) (generates 30 peaks)
11 //
12 // To execute only the first part of teh script (without fitting)
13 // specify a negative value for the number of peaks, eg
14 // root > .x peaks.C(-20)
15 //
16 //Author: Rene Brun
17 #include "Riostream.h"
18 #include <stdio.h>
19 #include "TSystem.h"
20 #include "TCanvas.h"
21 #include "TH1.h"
22 #include "TF1.h"
23 #include "TRandom.h"
24 #include "TSpectrum.h"
25 #include "TVirtualFitter.h"
26 #include "TMath.h"
27 #include "TDirectory.h"
28 #include "TKey.h"
29 ofstream out;
30 Int_t npeaks = 30;
31 TCanvas *c1 = 0;
32 TCanvas *c2 = 0;
33 //________________________________________________________________________________
34 Double_t fpeaks(Double_t *x, Double_t *par) {
35  Double_t result = par[0] + par[1]*x[0];
36  for (Int_t p=0;p<npeaks;p++) {
37  Double_t norm = par[3*p+2];
38  Double_t mean = par[3*p+3];
39  Double_t sigma = par[3*p+4];
40  result += norm*TMath::Gaus(x[0],mean,sigma);
41  }
42  return result;
43 }
44 //________________________________________________________________________________
45 void peaks(TH1* h, Int_t np=10, Int_t ij=0) {
46  if (! h) return;
47  npeaks = TMath::Abs(np);
48  if (! c1) c1 = new TCanvas();
49  else c1->Clear();
50  if (c2 && ij > 0) {c2->cd(ij); h->Draw(); c2->Update();}
51  c1->Divide(1,2);
52  c1->cd(1);
53  h->Draw();
54  TH1F *h2 = (TH1F*)h->Clone("h2");
55  //Use TSpectrum to find the peak candidates
56  TSpectrum *s = new TSpectrum(2*npeaks);
57  Int_t nfound = s->Search(h,5,"",0.05);
58  printf("Found %d candidate peaks to fit\n",nfound);
59  if (! nfound) return;
60  //Estimate background using TSpectrum::Background
61  TH1 *hb = s->Background(h,20,"same");
62  hb->Draw("same");
63  c1->Update();
64  if (c2 && ij > 0) {c2->cd(ij); h->Draw(); c2->Update();}
65  if (np <0) return;
66 
67  //estimate linear background using a fitting method
68  c1->cd(2);
69  TF1 *fline = new TF1("fline","pol1",0,1000);
70  fline->FixParameter(1,0.);
71  h->Fit("fline","qnlw");
72  if (c2 && ij > 0) {c2->cd(ij+1); h->Draw(); c2->Update(); c1->cd(2);}
73  //Loop on all found peaks. Eliminate peaks at the background level
74  Double_t par[3000];
75  par[0] = fline->GetParameter(0);
76  par[1] = fline->GetParameter(1);
77  npeaks = 0;
78  Float_t *xpeaks = s->GetPositionX();
79  Float_t ymax = 0;
80  for (Int_t p=0;p<nfound;p++) {
81  Float_t xp = xpeaks[p];
82  Int_t bin = h->GetXaxis()->FindBin(xp);
83  Float_t yp = h->GetBinContent(bin);
84  if (yp-3*TMath::Sqrt(yp) < fline->Eval(xp)) continue;
85  par[3*npeaks+2] = yp;
86  if (yp > ymax) ymax = yp;
87  par[3*npeaks+3] = xp;
88  par[3*npeaks+4] = 3;
89  npeaks++;
90  }
91  cout << "Found " << npeaks << " useful peaks to fit" << endl;
92  Int_t NP = 0;
93  Int_t nbad = 0;
94  TString LineH(" {\""); LineH += h->GetName(); LineH += "\"";
95  TString Line("");
96  struct ParErr_t {Double_t par, err;};
97  ParErr_t parErr[10];
98  TF1 *fit = 0;
99  if (ymax > 2*par[0]) {
100  cout << "Now fitting: Be patient" << endl;
101  fit = new TF1("fit",fpeaks,0,1000,2+3*npeaks);
102  TVirtualFitter::Fitter(h2,10+3*npeaks); //we may have more than the default 25 parameters
103  fit->SetParameter(0,par[0]);
104  fit->FixParameter(1,0.);
105  for (Int_t p = 0; p < npeaks; p++) {
106  fit->SetParName(3*p+2,Form("A%i",p));
107  fit->SetParLimits(3*p+2,0,1e6);
108  fit->SetParameter(3*p+2,par[3*p+2]);
109  fit->SetParName(3*p+3,Form("#mu%i",p));
110  fit->SetParLimits(3*p+3,TMath::Max(0.,par[3*p+3]-2), TMath::Min(240.,par[3*p+3]+2));
111  fit->SetParameter(3*p+3,par[3*p+3]);
112  fit->SetParName(3*p+4,Form("#sigma%i",p));
113  fit->SetParLimits(3*p+4,0.01,20);
114  fit->SetParameter(3*p+4,par[3*p+4]);
115  }
116  fit->SetNpx(1000);
117  h2->SetStats(0);
118  h2->Fit("fit","l");
119  if (c2 && ij > 0) {c2->cd(ij); h2->Draw("same"); c2->Update(); c1->cd(2);}
120  fit->GetParameters(par);
121  for (Int_t p = 0; p<np;p++) {
122  if (p < npeaks && par[3*p+2] > 5*fit->GetParError(3*p+2) &&
123  par[3*p+2] > par[0]) {
124  if (TMath::Abs(par[3*p+4]) > 5.0) nbad++;
125  // Line += Form(",%f,%f,%7.2f,%5.2f",par[3*p+2],fit->GetParError(3*p+2),par[3*p+3],TMath::Abs(par[3*p+4]));
126  parErr[NP].par = par[3*p+3];
127  parErr[NP].err = TMath::Abs(par[3*p+4]);
128  for (Int_t i = 0; i < NP; i++) {
129  if (parErr[i].par > parErr[NP].par) {
130  ParErr_t temp = parErr[i];
131  parErr[i] = parErr[NP];
132  parErr[NP] = temp;
133  }
134  }
135  NP++;
136  }
137  }
138  }
139  for (Int_t p = 0; p < np; p++) {
140  if (p < NP) Line += Form(",%7.2f,%5.2f",parErr[p].par,parErr[p].err);
141  else Line += ",0,0";
142  }
143  Line += "},";
144  if (nbad > 1) NP = -NP; // reject whole hybrid
145  LineH += Form(",%3i",NP);
146  cout << LineH << Line << endl;
147  out << LineH << Line << endl;
148  c1->Update();
149  if (c2) c2->Update();
150  delete fit;
151  delete s;
152 }
153 //________________________________________________________________________________
154 void FindPeaks(Int_t np=10) {
155  TString Out("Results.");
156  Out += gSystem->BaseName(gDirectory->GetName());
157  Out.ReplaceAll(".root","");
158  Out.ReplaceAll(" ","");
159  out.open(Out, ios::out);
160  TIter nextkey(gDirectory->GetListOfKeys() );
161  TKey *key;
162  Int_t i = 0;
163  Int_t nhists = 0;
164  while ((key = (TKey*) nextkey())) {
165  TObject *obj = key->ReadObj();
166  if (! obj->IsA()->InheritsFrom( "TH1" ) ) continue;
167  TH1 *h1 = (TH1*)obj;
168  if (h1->GetEntries() < 1000) continue;
169  nhists++;
170  }
171  Int_t nxy = nhists;
172  Int_t nx = (Int_t) TMath::Sqrt(nxy) + 1;
173  Int_t ny = nxy/nx + 1;
174  c2 = new TCanvas("Anodes","Anodes");
175  c2->Divide(nx,ny);
176  nextkey.Reset();
177  Int_t histNo = 1;
178  while ((key = (TKey*) nextkey())) {
179  TObject *obj = key->ReadObj();
180  if (! obj->IsA()->InheritsFrom( "TH1" ) ) continue;
181  if ( obj->IsA()->InheritsFrom( "TH1C" ) ||
182  obj->IsA()->InheritsFrom( "TH1S" ) ||
183  obj->IsA()->InheritsFrom( "TH1I" ) ||
184  obj->IsA()->InheritsFrom( "TH1F" ) ) {
185  cout << "Found histogram " << obj->GetName() << endl;
186  TH1 *h1 = (TH1*)obj;
187  if (h1->GetEntries() < 1000) {
188  TString LineH(" {\""); LineH += h1->GetName(); LineH += "\"";
189  TString Line("");
190  LineH += ",-99";
191  for (Int_t p = 0; p < np; p++) Line += ",0,0";
192  Line += "},// dead";
193  cout << LineH << Line << endl;
194  out << LineH << Line << endl;
195  continue;
196  }
197  peaks(h1,np,histNo);
198  histNo++;
199  // break;
200  if (i < 99) {
201  cout << "Type something" << endl;
202  cin >> i;
203  if (i <= 0) break;
204  }
205  }
206  }
207  out.close();
208 }
Line
STAR includes.
Definition: Line.hh:23
Generated by   doxygen 1.8.5