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

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