Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
peeks.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file peeks.C
1 // Illustrates how to find peaks in histograms.
2 // This script generates a random number of gaussian peaks
3 // on top of a linear background.
4 // The position of the peaks is found via TSpectrum and injected
5 // as initial values of parameters to make a global fit.
6 // The background is computed and drawn on top of the original histogram.
7 //
8 // To execute this example, do
9 // root > .x peaks.C (generate 10 peaks by default)
10 // root > .x peaks.C++ (use the compiler)
11 // root > .x peaks.C++(30) (generates 30 peaks)
12 //
13 // To execute only the first part of the script (without fitting)
14 // specify a negative value for the number of peaks, eg
15 // root > .x peaks.C(-20)
16 //
17 //Author: Rene Brun
18 
19 #include "TCanvas.h"
20 #include "TMath.h"
21 #include "TH1.h"
22 #include "TF1.h"
23 #include "TRandom.h"
24 #include "TSpectrum.h"
25 #include "TVirtualFitter.h"
26 
27 Int_t npeaks = 30;
28 Double_t fpeaks(Double_t *x, Double_t *par) {
29  Double_t result = par[0] + par[1]*x[0];
30  for (Int_t p=0;p<npeaks;p++) {
31  Double_t norm = par[3*p+2];
32  Double_t mean = par[3*p+3];
33  Double_t sigma = par[3*p+4];
34  result += norm*TMath::Gaus(x[0],mean,sigma);
35  }
36  return result;
37 }
38 void peaks(Int_t np=10) {
39  npeaks = TMath::Abs(np);
40  TH1F *h = new TH1F("h","test",500,0,1000);
41  //generate n peaks at random
42  Double_t par[3000];
43  par[0] = 0.8;
44  par[1] = -0.6/1000;
45  Int_t p;
46  for (p=0;p<npeaks;p++) {
47  par[3*p+2] = 1;
48  par[3*p+3] = 10+gRandom->Rndm()*980;
49  par[3*p+4] = 3+2*gRandom->Rndm();
50  }
51  TF1 *f = new TF1("f",fpeaks,0,1000,2+3*npeaks);
52  f->SetNpx(1000);
53  f->SetParameters(par);
54  TCanvas *c1 = new TCanvas("c1","c1",10,10,1000,900);
55  c1->Divide(1,2);
56  c1->cd(1);
57  h->FillRandom("f",200000);
58  h->Draw();
59  TH1F *h2 = (TH1F*)h->Clone("h2");
60  //Use TSpectrum to find the peak candidates
61  TSpectrum *s = new TSpectrum(2*npeaks);
62  Int_t nfound = s->Search(h,2,"",0.10);
63  printf("Found %d candidate peaks to fit\n",nfound);
64  //Estimate background using TSpectrum::Background
65  TH1 *hb = s->Background(h,20,"same");
66  if (hb) c1->Update();
67  if (np <0) return;
68 
69  //estimate linear background using a fitting method
70  c1->cd(2);
71  TF1 *fline = new TF1("fline","pol1",0,1000);
72  h->Fit("fline","qn");
73  //Loop on all found peaks. Eliminate peaks at the background level
74  par[0] = fline->GetParameter(0);
75  par[1] = fline->GetParameter(1);
76  npeaks = 0;
77  Double_t *xpeaks = s->GetPositionX();
78  for (p=0;p<nfound;p++) {
79  Double_t xp = xpeaks[p];
80  Int_t bin = h->GetXaxis()->FindBin(xp);
81  Double_t yp = h->GetBinContent(bin);
82  if (yp-TMath::Sqrt(yp) < fline->Eval(xp)) continue;
83  par[3*npeaks+2] = yp;
84  par[3*npeaks+3] = xp;
85  par[3*npeaks+4] = 3;
86  npeaks++;
87  }
88  printf("Found %d useful peaks to fit\n",npeaks);
89  printf("Now fitting: Be patient\n");
90  TF1 *fit = new TF1("fit",fpeaks,0,1000,2+3*npeaks);
91  //we may have more than the default 25 parameters
92  TVirtualFitter::Fitter(h2,10+3*npeaks);
93  fit->SetParameters(par);
94  fit->SetNpx(1000);
95  h2->Fit("fit");
96 }
97 
98 
99 
100 
101 
102 // ;; This buffer is for notes you don't want to save, and for Lisp evaluation.
103 // ;; If you want to create a file, visit that file with C-x C-f,
104 // ;; then enter the text in that file's own buffer.
105 
106 // **************************************************************************
107 // assume background is linear what is most likely not true ....
108 // par[0] background constant
109 // [1] background slope
110 // assume the signal is the sequence of equally spaced gaussian peaks
111 // par[2] peak to peak separation
112 // [3] normalization constant for the first peak (pedestal)
113 // [4] location of the first peak
114 // [5] rms of the first peak - at this stage we will allow it to vary peak-to-peak
115 // then 2 parameters per peak (normalization and RMS)
116 Double_t sipmpeaks(Double_t *x, Double_t *par) {
117  Double_t result = par[0] + par[1]*x[0];
118  // add pedestal peak first
119  result += par[3]*TMath::Gaus(x[0],par[4],par[5]);
120  for (Int_t p=0; p<npeaks; p++) {
121  Double_t norm = par[6 + p*2];
122  Double_t mean = par[4] + par[2]*(p+1);
123  Double_t sigma = par[8 + p*2];
124  result += norm*TMath::Gaus(x[0],mean,sigma);
125  }
126  return result;
127 }
129 // do single cell calibration using low end data stored in rv_# and fm_# channel histograms
130 // the assumption - we have beautiful picture with at least 20 peaks !!!!!!
131 Int_t spmcalib(TString hName)
132 {
133  TH1F * data = (TH1F*)(gSystem->FindObject(hName));
134  if(!data){
135  cout<<"spmCalib Data Histogramm "<<hName<<" not found"<<endl;
136  exit();
137  }
138  Int_t maxPeaks(20); // maximum number of peaks we have chance to resolve
139  Double_t range(150.); // range of amplitudes for peak search
140  Double_t Double_t norm0(1.), slope0(-1./range);
141  Double_t peak_to_peak(10.);
142  Double_t pLoc(0.), pRMS(1.5);
143  Double_t par[6+maxPeaks*2]; // backgrownd, pedestal and nPeaks of gaussian peaks with constant peak-to-peak separation
144 
145  par[0] = norm0;
146  par[1] = slope0;
147  par[2] = peak_to_peak; // peak-to-peak separation
148  par[3] = norm0; // normalization for pedestal peak
149  par[4] = pLoc; // pedestal mean
150  par[5] = pRMS; // pedestal rms
151  for (int p=0; p< maxPeaks; p++) {
152  par[6 + p*3] = norm0; // norm of the peak p
153  par[7 + p*3] = par[4] + par[2]*(p+1); // mean of the peak p
154  par[8 + p*3] = 2.; // rms of the peak p
155  }
156  TString scId = "sc_"; scId += hName;
157  TF1 * scF = (TF1*) (gROOT->FindObject(scId)); if(scF) delete scF;
158  scF = new TF1(scId,&sipmpeaks,0,range,6+maxPeaks*2);
159  TH1F * dataCl = (TH1F*)(data->Clone("dataCl"));
160  dataCl -> SetAxisRange(0., range);
161  Int_t binsToUse = dataCl->GetXaxis()->FindBin(range);
162  TCanvas * peaks = new TCanvas("peaks","peaks",10,10,1000,900);
163  peaks -> Divide(1,2);
164  peaks -> cd(1);
165  dataCl -> Draw();
166  TH1F * dataCl2 = (TH1F*)data->Clone("dataCl2");
167 
168  // Use TSpectrum to find the peak candidates
169  TSpectrum *s = new TSpectrum(maxPeaks+1);
170  Int_t nfound = s->Search(dataCl, 2, "", 0.10);
171  printf("TSpectrum: Found %d candidate peaks to fit\n",nfound);
172 
173  //Estimate background using TSpectrum::Background
174  TH1 *hb = s->Background(h,20,"same");
175  if (hb) peaks->Update();
176  // if (np <0) return;
177  // function to estimate linear background
178  peaks -> cd(2);
179  TF1 *fline = new TF1("fline","pol1",0,range);
180  dataCl -> Fit("fline","qn");
181  //Loop on all found peaks. Eliminate peaks at the background level
182  par[0] = fline->GetParameter(0);
183  par[1] = fline->GetParameter(1);
184  Int_t npeaks = 0;
185  Double_t *xpeaks = s->GetPositionX();
186  Double_t spacing = 0.;
187  for (p=0; p<nfound; p++) {
188  Double_t xp = xpeaks[p];
189  Int_t bin = dataCl->GetXaxis()->FindBin(xp);
190  Double_t yp = dataCl->GetBinContent(bin);
191  cout<<"Testing Peak "<<p<<" xp "<<xp<<" bin "<<bin<<" yP "<<yp<<" BCKG "<<fline->Eval(xp)<<endl;
192  if (yp-TMath::Sqrt(yp) < fline->Eval(xp)) continue;
193  if(npeaks) spacing += xp - par[4 + (npeaks-1)+3];
194  par[3 + npeaks*3] = yp;
195  par[4 + npeaks*3] = xp;
196  par[5 + npeaks*3] = 3.;
197  npeaks++;
198  }
199  printf("Found %d useful peaks to fit\n",npeaks);
200  printf("Now fitting: Be patient\n");
201  TF1 *fit = new TF1("fit",&sipmpeaks,0,1000,3+3*npeaks); // first peak is the pedestal
202  //we may have more than the default 25 parameters
203  TVirtualFitter::Fitter(h2,10+3*npeaks);
204  fit->SetParameters(par);
205  fit->SetNpx(1000);
206  h2->Fit("fit");
207 }
208 
209 
210 
211 
212 
214 // do single cell calibration using low end data stored in rv_# and fm_# channel histograms
215 // the assumption - we have beautiful picture with at least 20 peaks !!!!!!
216 void sccalib(nPeaks = 20)
217 {
218  Double_t rcpar[2+3+2*nPeaks]; // backgrownd and nPeaks of uncorrelated gaussian peaks
219  Double_t fcpar[2+3+2*nPeaks]; // backgrownd and nPeaks of uncorrelated gaussian peaks
220 
221  int nx_c = 2;
222  int ny_c = NCH;
223  TString cSC = "Single Cell Calibration (black - raw, red - fit) Run "; cSC += runId.Data();
224  scCalibDisplay = new TCanvas("scCalibDisplay",cSC,200*nx_c,100*ny_c);
225  scCalibDisplay->Divide(nx_c, ny_c);
226  gStyle->SetOptStat(1110);
227  gStyle->SetOptFit(111);
228  Int_t binsToUse(200); Double_t slope0(-1./(Double_t))binsToUse); Double_t norm0(1.);
229  Double_t peak_to_peak(10.);
230  for (int ich = 0; ich<NCH; ich++){
231  TH1F *hrv = (TH1F*)rpeak[ich];
232  hrv -> SetLineColor(1);
233  hrv -> SetAxisRange(0., binsToUse);
234  TH1F *fmv = (TH1F*)fm[ich];
235  fmv -> SetLineColor(2);
236  fmv -> SetAxisRange(0., binsToUse);
237 
238  // estimate exponential backgrownd we go to 500 ADC counts
239  rcpar[0] = norm0;
240  rcpar[1] = slope0;
241  rcpar[2] = peak_to_peak; // peak-to-peak separation
242  rcpar[3] = norm0; // normalization for pedestal peak
243  rcpar[4] = 2.; // pedestal mean
244  rcpar[5] = 1.5; // pedestal rms
245  // we fit peak-to-peak not the individual independent peak positions
246  for (int p=1; p<= nPeaks; p++) {
247  rcpar[2*p+4] = norm0; // norm of the peak p
248  rcpar[2*p+5] = 2.; // rms of the peak p
249  }
250  fcpar[0] = 1.;
251  fcpar[1] = slope0;
252  fcpar[2] = peak_to_peak; // peak-to-peak separation
253  fcpar[3] = norm0; // normalization for pedestal peak
254  fcpar[4] = 2.; // pedestal mean
255  fcpar[5] = 1.5; // pedestal rms
256  // we fit peak-to-peak not the individual independent peak positions
257  for (int p=1; p<= nPeaks; p++) {
258  fcpar[2*p+4] = norm0; // norm of the peak p
259  fcpar[2*p+5] = 2.; // rms of the peak p
260  }
261 
262  TString rvscId = "rvsc_"; rvscId += ich;
263  TString fmscId = "fmsc_"; fmscId += ich;
264  rvsc[ich] = (TF1*) (gROOT->FindObject(rvscId)); if(rvsc[ich]) delete rvsc[ich];
265  fmsc[ich] = (TF1*) (gROOT->FindObject(fmscId)); if(fmsc[ich]) delete fmsc[ich];
266  rvsc[ich] = new TF1(rvscId,&scpeaks,0,500,5+2*nPeaks); fmsc[ich] = new TF1(fmscId,&scpeaks,0,500,5+2*nPeaks);
267  rvsc[ich] -> SetNpx(binsToUse); fmsc[ich] -> SetNpx(binsToUse);
268  rvsc[ich] -> SetParameters(rcpar); fmsc[ich] -> SetParameters(fcpar);
269  rvsc[ich] -> SetLineColor(1); fmsc[ich] -> SetLineColor(2);
270  TH1F * rh = (TH1F*)(hrv->Clone("rh")); TH1F * fh = (TH1F*)(fmv->Clone("fh"));
271 
272  // function to estimate linear background
273  TF1 *fline = new TF1("fline","pol1",0,binsToUse);
274  // Use TSpectrum to find the peak candidates
275  TSpectrum *s = new TSpectrum(2*nPeaks);
276  scCalibDisplay -> cd(ich*2+1); hrv -> Draw();
277  Int_t rfound = s->Search(hrv,1,"new");
278  printf("Found %d candidate peaks in Raw hist ",rfound);
279  scCalibDisplay -> Update();
280  hrv -> Fit("fline","qn");
281  rcpar[0] = fline->GetParameter(0);
282  rcpar[1] = fline->GetParameter(1);
283  //Loop on all found peaks. Eliminate peaks at the background level
284  nrpeaks = 0;
285  Float_t *xrpeaks = s->GetPositionX();
286  for (p=0;p<nfound;p++) {
287  Float_t xp = xpeaks[p];
288  Int_t bin = hrv->GetXaxis()->FindBin(xp);
289  Float_t yp = h->GetBinContent(bin);
290  if (yp-TMath::Sqrt(yp) < fline->Eval(xp)) continue;
291  par[2*npeaks+2] = yp;
292  par[3*npeaks+3] = xp;
293  par[3*npeaks+4] = 3;
294  nrpeaks++;
295  }
296 
297 
298 
299 
300  scCalibDisplay -> cd(ich*2+2); fmv -> Draw();
301  int_t ffound = s->Search(fmv,1,"new");
302  printf("Found %d candidate peaks in Fit hist ",ffound);
303  scCalibDisplay -> Update();
304  fmv -> Fit("fline","qn");
305  fcpar[0] = fline->GetParameter(0);
306  fcpar[1] = fline->GetParameter(1);
307 
308  //Loop on all found peaks. Eliminate peaks at the background level
309  npeaks = 0;
310  Float_t *xpeaks = s->GetPositionX();
311  for (p=0;p<nfound;p++) {
312  Float_t xp = xpeaks[p];
313  Int_t bin = h->GetXaxis()->FindBin(xp);
314  Float_t yp = h->GetBinContent(bin);
315  if (yp-TMath::Sqrt(yp) < fline->Eval(xp)) continue;
316  par[3*npeaks+2] = yp;
317  par[3*npeaks+3] = xp;
318  par[3*npeaks+4] = 3;
319  npeaks++;
320  }
321 
322 
323 
324 
325 
326 
327 
328 
329  // rvsc[ich] -> SetParNames ("#mu","g","#sigma","A"); fmsc[ich] -> SetParNames ("#mu","g","#sigma","A");
330  scCalibDisplay -> cd(ich+1); hrv->SetAxisRange(0., 200); hrv ->Draw(); fmv->Draw("same");
331  hrv -> Fit(rvsc[ich],"rnQ"); fmv -> Fit(fmsc[ich],"rnq");
332  rvsc[ich] -> Draw("same"); fmsc[ich] -> Draw("same");
333  }
334  scCalibDisplay -> Update();
335 }