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

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