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

Last change on this file since 12204 was 12166, checked in by neise, 13 years ago
initial commit
File size: 9.9 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 datafile and plot the DRS-calibrated data
24//ATTENTION: only works for ROI=1024
25// (array indices of the calibration wrong otherwise)
26//******************************************************************************
27
28 gROOT->SetStyle("Plain");
29
30//-------------------------------------------
31//Define the data variables
32//-------------------------------------------
33 vector<int16_t> data;
34 vector<int16_t> data_offset;
35 unsigned int data_num;
36 size_t data_n;
37 UInt_t data_px;
38 UInt_t data_roi;
39
40//-------------------------------------------
41//Get the DRS calibration
42//-------------------------------------------
43 size_t drs_n;
44 vector<float> drs_basemean;
45 vector<float> drs_gainmean;
46 vector<float> drs_triggeroffsetmean;
47 FOpenCalibFile(drsname, drs_basemean, drs_gainmean, drs_triggeroffsetmean, drs_n);
48
49//-------------------------------------------
50//Check the sizes of the data columns
51//-------------------------------------------
52// if(drs_n!=data_n)
53// {
54// cout << "Data and DRS file incompatible (Px*ROI disagree)" << endl;
55// return 1;
56// }
57
58//-------------------------------------------
59//Open the file
60//-------------------------------------------
61 fits datafile(name);
62 if (!datafile)
63 {
64 cout << "Couldn't properly open the datafile." << endl;
65 return 1;
66 }
67 FOpenDataFile(datafile, data, data_offset, data_num, data_n, data_roi, data_px);
68
69//-------------------------------------------
70//Create the canvas, plot and title
71//-------------------------------------------
72 char title_spect[500];
73 char name_spect[50];
74 std::sprintf(title_spect,"Spectrum of Px %i (data: %s, DRS: %s)",pixelnr, name,drsname);
75 TCanvas *canv_spect = new TCanvas( "canv_spect", "Spectrum", 950, 10, 700, 500 );
76 TH1F* spectrum[data_px];
77 for(int i=0; i<data_px; i++) {
78 std::sprintf(title_spect,"Spectrum of Px %i (data: %s, DRS: %s)",i,name,drsname);
79 std::sprintf(name_spect,"spect%i",i);
80 spectrum[i] = new TH1F(name_spect,title_spect,300,0,600);
81 spectrum[i]->GetXaxis()->SetTitle("Pulse size (a.u.)");
82 spectrum[i]->GetYaxis()->SetTitle("Counts");
83 }
84
85 TCanvas *canv_gain= new TCanvas( "canv_gain", "Gain in arbitrary units", 50, 500, 1800, 800 );
86//-------------------------------------------
87//Define the fit
88//-------------------------------------------
89 int fit_start = 20;
90 int fit_end = 200;
91 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);
92 Double_t par_init[5+2*n_peaks] = {70,90,100,15,10,15,10,30,30};
93 Double_t parread_tmp[5+2*n_peaks];
94 Double_t parread[data_px][5+2*n_peaks];
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
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 Double_t parread_dist[data_px];
177 Double_t parread_x[data_px];
178 for(int i=0; i<data_px; i++) {
179 f_peaks->SetParameters(par_init);
180 if(!(spectrum[i]->Fit(f_peaks,"NR"))) {
181// std::cout << "Fit successful." << std::endl;
182 f_peaks->GetParameters(&parread_tmp[0]);
183 for(int j=0; j<5+2*n_peaks; j++) {
184 parread[i][j] = parread_tmp[j];
185 }
186 parread_dist[i] = parread_tmp[1];
187 parread_x[i] = i;
188 }
189 else {
190 for(int j=0; j<5+2*n_peaks; j++) {
191 parread[i][j] = 0;
192 }
193 parread_dist[i] = 0;
194 parread_x[i] = i;
195 }
196 }
197
198 canv_gain->cd();
199 TGraph* plot_gain = new TGraph(data_px,parread_x,parread_dist);
200 plot_gain->GetYaxis()->SetRangeUser(0,150);
201 plot_gain->Draw("A*");
202 canv_gain->Modified();
203 canv_gain->Update();
204
205 canv_spect->cd();
206 canv_spect->SetLogy();
207
208//-------------------------------------------
209//Draw the data where the gain is out of limits
210//-------------------------------------------
211 char temp = 'x';
212 std::cout << "Plot the spectra out of limits: Enter for next, 'a' to abort." << std::endl;
213 for(int i=0; i<data_px; i++) {
214 if((parread_dist[i]<80)||(parread_dist[i]>105)) {
215 std::cout << "Gain out of limits: Px " << i << " at " << parread_dist[i] << std::endl;
216 spectrum[i]->Draw();
217
218 f_peaks->SetParameters(&parread[i][0]);
219 f_peak[0]->SetParameters(parread[i][2],parread[i][0],parread[i][3]);
220 f_peak[1]->SetParameters(parread[i][4],parread[i][0]+parread[i][1],parread[i][5]);
221// f_peak[2]->SetParameters(parread[i][6],parread[i][0]+2*parread[i][1],parread[i][7]);
222// f_peak[3]->SetParameters(parread[i][8],parread[i][0]+3*parread[i][1],parread[i][9]);
223// f_peak[4]->SetParameters(parread[i][10],parread[i][12],parread[i][11]);
224 f_pedestal->SetParameters(parread[i][6],parread[i][7],parread[i][8]);
225
226 f_peaks->Draw("same");
227 for(int i=0; i<n_peaks; i++) {
228 f_peak[i]->Draw("same");
229 }
230 f_pedestal->Draw("same");
231
232 canv_spect->Modified();
233 canv_spect->Update();
234 temp='x';
235 while ((temp!='\n') && (temp!='a'))
236 cin.get(temp);
237 if(temp=='a') break;
238 }
239 }
240//-------------------------------------------
241//Draw the data
242//-------------------------------------------
243 canv_spect->cd();
244 Int_t plotspect_nr;
245 std::cout << "Enter the number of the spectrum to plot (-1 to quit): ";
246 std::cin >> plotspect_nr;
247 while(plotspect_nr>=0) {
248 plotspect_nr = plotspect_nr % data_px;
249 spectrum[plotspect_nr]->Draw();
250
251 f_peaks->SetParameters(&parread[plotspect_nr][0]);
252 f_peak[0]->SetParameters(parread[plotspect_nr][2],parread[plotspect_nr][0],parread[plotspect_nr][3]);
253 f_peak[1]->SetParameters(parread[plotspect_nr][4],parread[plotspect_nr][0]+parread[plotspect_nr][1],parread[plotspect_nr][5]);
254// f_peak[2]->SetParameters(parread[plotspect_nr][6],parread[plotspect_nr][0]+2*parread[plotspect_nr][1],parread[plotspect_nr][7]);
255// f_peak[3]->SetParameters(parread[plotspect_nr][8],parread[plotspect_nr][0]+3*parread[plotspect_nr][1],parread[plotspect_nr][9]);
256// f_peak[4]->SetParameters(parread[plotspect_nr][10],parread[plotspect_nr][12],parread[plotspect_nr][11]);
257 f_pedestal->SetParameters(parread[plotspect_nr][6],parread[plotspect_nr][7],parread[plotspect_nr][8]);
258
259 f_peaks->Draw("same");
260 for(int i=0; i<n_peaks; i++) {
261 f_peak[i]->Draw("same");
262 }
263 f_pedestal->Draw("same");
264
265 canv_spect->Modified();
266 canv_spect->Update();
267
268 std::cout << "Enter the number of the spectrum to plot (-1 to quit): ";
269 std::cin >> plotspect_nr;
270 }
271
272 return 0;
273}
Note: See TracBrowser for help on using the repository browser.