source: fact/tools/rootmacros/gainanalysis.C@ 18481

Last change on this file since 18481 was 12303, checked in by kraehenb, 13 years ago
Save spectra of gainanalysis to the file gainspectra.root.
File size: 10.7 KB
Line 
1#include <TROOT.h>
2#include <TCanvas.h>
3#include <TProfile.h>
4#include <TH1.h>
5#include <TH2.h>
6#include <TF1.h>
7#include <TGraph.h>
8#include <TFile.h>
9
10#include <stdint.h>
11#include <cstdio>
12
13#define HAVE_ZLIB
14#include "fits.h"
15#include "FOpenDataFile.c"
16#include "FOpenCalibFile.c"
17#include "FPedestal.c"
18#include "FOscilloscopeAllPx.c"
19
20#define n_peaks 2
21
22int gainanalysis(const char *name, const char *drsname, size_t pixelnr)
23{
24//******************************************************************************
25//Read a file with G-APD singles, analyse them and plot the spectra
26//This source code is best read with a wide editor (>100 characters per line) or a horizontal scroll bar.
27//ATTENTION: this script only works for ROI=1024
28// (array indices of the calibration wrong otherwise)
29//******************************************************************************
30
31 gROOT->SetStyle("Plain");
32 TFile outfile("gainspectra.root","RECREATE");
33//-------------------------------------------
34//Define the data variables
35//-------------------------------------------
36 vector<int16_t> data;
37 vector<int16_t> data_offset;
38 unsigned int data_num;
39 size_t data_n;
40 UInt_t data_px;
41 UInt_t data_roi;
42
43//-------------------------------------------
44//Get the DRS calibration
45//-------------------------------------------
46 size_t drs_n;
47 vector<float> drs_basemean;
48 vector<float> drs_gainmean;
49 vector<float> drs_triggeroffsetmean;
50 FOpenCalibFile(drsname, drs_basemean, drs_gainmean, drs_triggeroffsetmean, drs_n);
51
52//-------------------------------------------
53//Check the sizes of the data columns
54//-------------------------------------------
55// if(drs_n!=data_n)
56// {
57// cout << "Data and DRS file incompatible (Px*ROI disagree)" << endl;
58// return 1;
59// }
60
61//-------------------------------------------
62//Open the file
63//-------------------------------------------
64 fits datafile(name);
65 if (!datafile)
66 {
67 cout << "Couldn't properly open the datafile." << endl;
68 return 1;
69 }
70 FOpenDataFile(datafile, data, data_offset, data_num, data_n, data_roi, data_px);
71
72//-------------------------------------------
73//Create the canvas, plot and title
74//-------------------------------------------
75 char title_spect[500];
76 char name_spect[50];
77 std::sprintf(title_spect,"Spectrum of Px %i (data: %s, DRS: %s)",pixelnr, name,drsname);
78 TCanvas *canv_spect = new TCanvas( "canv_spect", "Spectrum", 950, 10, 700, 500 );
79 TH1F* spectrum[data_px];
80 for(int i=0; i<data_px; i++) {
81 std::sprintf(title_spect,"Spectrum of Px %i (data: %s, DRS: %s)",i,name,drsname);
82 std::sprintf(name_spect,"spect%i",i);
83 spectrum[i] = new TH1F(name_spect,title_spect,350,-100,600);
84 spectrum[i]->GetXaxis()->SetTitle("Pulse size (a.u.)");
85 spectrum[i]->GetYaxis()->SetTitle("Counts");
86 }
87
88 TCanvas *canv_gain= new TCanvas( "canv_gain", "Gain in arbitrary units", 50, 500, 1800, 800 );
89//-------------------------------------------
90//Define the fit
91//-------------------------------------------
92 int fit_start = 20;
93 int fit_end = 200;
94 TF1 *f_peaks = new TF1("total","[2]*TMath::Exp(-(x-[0])*(x-[0])/(2*[3]*[3])) + [4]*TMath::Exp(-(x-[0]-[1])*(x-[0]-[1])/(2*[5]*[5])) + [6]*TMath::Exp((-x-[7])/[8])",fit_start,fit_end);
95 Double_t par_init[5+2*n_peaks] = {70,90,100,15,10,15,10,30,30};
96// TF1 *f_peaks = new TF1("total","[2]*TMath::Exp(-(x-[0])*(x-[0])/(2*[3]*[3])) + [4]*TMath::Exp(-(x-[0]-[1])*(x-[0]-[1])/(2*[5]*[5])) + [6]*TMath::Exp(-(x-[0]-2*[1])*(x-[0]-2*[1])/(2*[7]*[7])) + [8]*TMath::Exp(-(x-[0]-3*[1])*(x-[0]-3*[1])/(2*[9]*[9])) + [10]*TMath::Exp(-(x-[12])*(x-[12])/(2*[11]*[11]))",fit_start,fit_end);
97// Double_t par_init[3+2*n_peaks] = {70,90,500000,15,50000,15,5000,15,500,15,800000,20,0};
98 f_peaks->SetParameters(par_init);
99 f_peaks->SetLineColor(2);
100
101 TF1 *f_pedestal = new TF1("pedestal","[0]*TMath::Exp((-x-[1])/[2])",fit_start,fit_end);
102 f_pedestal->SetLineColor(2);
103 f_pedestal->SetLineWidth(1);
104
105 TF1* f_peak[n_peaks];
106 for(int i=0; i<n_peaks; i++) {
107 f_peak[i] = new TF1("peak","gaus",fit_start,fit_end);
108 f_peak[i]->SetLineColor(2);
109 f_peak[i]->SetLineWidth(1);
110 }
111
112//-------------------------------------------
113//Find the baseline
114//-------------------------------------------
115 float pedestal_mean, pedestal_rms;
116// FPedestal(title_base, datafile, data, data_offset, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, pixelnr, pedestal_mean, pedestal_rms);
117 //Return values
118 //Value of maximal probability: -2.225 +- 4.3144
119 pedestal_mean = -2.225;
120 pedestal_rms = 4.3144;
121
122//-------------------------------------------
123//Oscilloscope (search for peaks and fill integral into the histograms)
124//-------------------------------------------
125// float threshold = pedestal_mean+pedestal_rms;
126 float threshold = 2.5;
127
128//***********************
129 FOscilloscopeAllPx(datafile, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, threshold, spectrum);
130
131 canv_spect->cd();
132 spectrum[0]->Draw();
133 canv_spect->Modified();
134 canv_spect->Update();
135//***********************
136// fits datafile2("20110804_025.fits");
137// if (!datafile2)
138// {
139// cout << "Couldn't properly open the datafile2." << endl;
140// return 1;
141// }
142// FOpenDataFile(datafile2, data, data_offset, data_num, data_n, data_roi, data_px);
143// FOscilloscopeAllPx(datafile2, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, threshold, spectrum);//
144// canv_spect->cd();
145// spectrum[0]->Draw();
146// canv_spect->Modified();
147// canv_spect->Update();
148////***********************
149// fits datafile3("20110804_026.fits");
150// if (!datafile3)
151// {
152// cout << "Couldn't properly open the datafile3." << endl;
153// return 1;
154// }
155// FOpenDataFile(datafile3, data, data_offset, data_num, data_n, data_roi, data_px);
156// FOscilloscopeAllPx(datafile3, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, threshold, spectrum);
157//
158// canv_spect->cd();
159// spectrum[0]->Draw();
160// canv_spect->Modified();
161// canv_spect->Update();
162////***********************
163// fits datafile4("20110804_027.fits");
164// if (!datafile4)
165// {
166// cout << "Couldn't properly open the datafile4." << endl;
167// return 1;
168// }
169// FOpenDataFile(datafile4, data, data_offset, data_num, data_n, data_roi, data_px);
170// FOscilloscopeAllPx(datafile4, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, threshold, spectrum);
171//
172// canv_spect->cd();
173// spectrum[0]->Draw();
174// canv_spect->Modified();
175// canv_spect->Update();
176//***********************
177
178//-------------------------------------------
179//Fit the histograms and store the parameters into parread
180//parread_dist: peak difference (gain) per pixel, parread_x: x-Values for the TGraph
181//-------------------------------------------
182 Double_t parread[data_px][5+2*n_peaks];
183 Double_t parread_tmp[5+2*n_peaks];
184 Double_t parread_dist[data_px];
185 Double_t parread_x[data_px];
186 for(int i=0; i<data_px; i++) {
187 f_peaks->SetParameters(par_init);
188 if(!(spectrum[i]->Fit(f_peaks,"NR"))) {
189// std::cout << "Fit successful." << std::endl;
190 f_peaks->GetParameters(&parread_tmp[0]);
191 for(int j=0; j<5+2*n_peaks; j++) {
192 parread[i][j] = parread_tmp[j];
193 }
194 parread_dist[i] = parread_tmp[1];
195 parread_x[i] = i;
196 }
197 else {
198 for(int j=0; j<5+2*n_peaks; j++) {
199 parread[i][j] = 0;
200 }
201 parread_dist[i] = 0;
202 parread_x[i] = i;
203 }
204 }
205
206//-------------------------------------------
207//Plot the gain of all pixel
208//-------------------------------------------
209 canv_gain->cd();
210 TGraph* plot_gain = new TGraph(data_px,parread_x,parread_dist);
211 plot_gain->GetYaxis()->SetRangeUser(0,150);
212 plot_gain->Draw("A*");
213 canv_gain->Modified();
214 canv_gain->Update();
215
216 canv_spect->cd();
217 canv_spect->SetLogy();
218
219//-------------------------------------------
220//Draw the spectrum of pixel where the gain is out of limits
221//Drawn are the fitted function (f_peaks) and its components (f_peak[i], f_pedestal)
222//-------------------------------------------
223 char temp = 'x';
224 std::cout << "Plot the spectra out of limits: Enter for next, 'a' to abort." << std::endl;
225 for(int i=0; i<data_px; i++) {
226 if((parread_dist[i]<80)||(parread_dist[i]>105)) {
227 std::cout << "Gain out of limits: Px " << i << " at " << parread_dist[i] << std::endl;
228 spectrum[i]->Draw();
229
230 f_peaks->SetParameters(&parread[i][0]);
231 f_peak[0]->SetParameters(parread[i][2],parread[i][0],parread[i][3]);
232 f_peak[1]->SetParameters(parread[i][4],parread[i][0]+parread[i][1],parread[i][5]);
233// f_peak[2]->SetParameters(parread[i][6],parread[i][0]+2*parread[i][1],parread[i][7]);
234// f_peak[3]->SetParameters(parread[i][8],parread[i][0]+3*parread[i][1],parread[i][9]);
235// f_peak[4]->SetParameters(parread[i][10],parread[i][12],parread[i][11]);
236 f_pedestal->SetParameters(parread[i][6],parread[i][7],parread[i][8]);
237
238 f_peaks->Draw("same");
239 for(int i=0; i<n_peaks; i++) {
240 f_peak[i]->Draw("same");
241 }
242 f_pedestal->Draw("same");
243
244 canv_spect->Modified();
245 canv_spect->Update();
246 temp='x';
247 while ((temp!='\n') && (temp!='a'))
248 cin.get(temp);
249 if(temp=='a') break;
250 }
251 }
252
253//-------------------------------------------
254//Give the possibility to draw the spectrum of any pixel
255//-------------------------------------------
256 canv_spect->cd();
257 Int_t plotspect_nr;
258 std::cout << "Enter the number of the spectrum to plot (-1 to quit): ";
259 std::cin >> plotspect_nr;
260 while(plotspect_nr>=0) {
261 plotspect_nr = plotspect_nr % data_px;
262 spectrum[plotspect_nr]->Draw();
263
264 f_peaks->SetParameters(&parread[plotspect_nr][0]);
265 f_peak[0]->SetParameters(parread[plotspect_nr][2],parread[plotspect_nr][0],parread[plotspect_nr][3]);
266 f_peak[1]->SetParameters(parread[plotspect_nr][4],parread[plotspect_nr][0]+parread[plotspect_nr][1],parread[plotspect_nr][5]);
267// f_peak[2]->SetParameters(parread[plotspect_nr][6],parread[plotspect_nr][0]+2*parread[plotspect_nr][1],parread[plotspect_nr][7]);
268// f_peak[3]->SetParameters(parread[plotspect_nr][8],parread[plotspect_nr][0]+3*parread[plotspect_nr][1],parread[plotspect_nr][9]);
269// f_peak[4]->SetParameters(parread[plotspect_nr][10],parread[plotspect_nr][12],parread[plotspect_nr][11]);
270 f_pedestal->SetParameters(parread[plotspect_nr][6],parread[plotspect_nr][7],parread[plotspect_nr][8]);
271
272 f_peaks->Draw("same");
273 for(int i=0; i<n_peaks; i++) {
274 f_peak[i]->Draw("same");
275 }
276 f_pedestal->Draw("same");
277
278 canv_spect->Modified();
279 canv_spect->Update();
280
281 std::cout << "Enter the number of the spectrum to plot (-1 to quit): ";
282 std::cin >> plotspect_nr;
283 }
284 outfile.Write();
285 outfile.Close();
286 return 0;
287}
Note: See TracBrowser for help on using the repository browser.