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 |
|
---|
8 | #include <stdint.h>
|
---|
9 | #include <cstdio>
|
---|
10 |
|
---|
11 | #define HAVE_ZLIB
|
---|
12 | #include "fits.h"
|
---|
13 | #include "FOpenDataFile.c"
|
---|
14 | #include "FOpenCalibFile.c"
|
---|
15 | #include "FPedestal.c"
|
---|
16 | #include "FOscilloscope.c"
|
---|
17 |
|
---|
18 | #define n_peaks 4
|
---|
19 |
|
---|
20 | int peakfinderscope(const char *name, const char *drsname, size_t eventnr, 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 | //Create the canvas, plot and title
|
---|
60 | //-------------------------------------------
|
---|
61 | char title[500];
|
---|
62 | std::sprintf(title,"Data: %s, DRS: %s, Px %i Ev %i",name,drsname,pixelnr,eventnr);
|
---|
63 | TCanvas *canv = new TCanvas( "canv", "Waveform", 200, 10, 700, 500 );
|
---|
64 | // TProfile *pulseshape = new TProfile("pulseshape", title, 70, -10.5, 59.5);
|
---|
65 | TH2F *pulseshape = new TH2F("pulseshape",title,70,-10.5,59.5,350,-10.5,59.5);
|
---|
66 | pulseshape->GetXaxis()->SetTitle("Time after trigger (bins @2 GHz)");
|
---|
67 | pulseshape->GetYaxis()->SetTitle("Amplitude (mV)");
|
---|
68 |
|
---|
69 | char title_base[500];
|
---|
70 | std::sprintf(title_base,"All samples of Px %i (data: %s, DRS: %s)",pixelnr,name,drsname);
|
---|
71 |
|
---|
72 | char title_spect[500];
|
---|
73 | std::sprintf(title_spect,"Spectrum of all pixels (data: %s, DRS: %s)",name,drsname);
|
---|
74 | TCanvas *canv_spect = new TCanvas( "canv_spect", "Spectrum", 950, 10, 700, 500 );
|
---|
75 | TH1F *spectrum = new TH1F("spectrum", title_spect,300,0,600);
|
---|
76 | spectrum->GetXaxis()->SetTitle("Pulse size (a.u.)");
|
---|
77 | spectrum->GetYaxis()->SetTitle("Counts");
|
---|
78 |
|
---|
79 | //-------------------------------------------
|
---|
80 | //Define the fit
|
---|
81 | //-------------------------------------------
|
---|
82 | int fit_start = 30;
|
---|
83 | int fit_end = 380;
|
---|
84 | 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-[11])/[12])",fit_start,fit_end);
|
---|
85 | Double_t par[5+2*n_peaks] = {70,90,500000,15,50000,15,5000,15,500,15,10000,30,30};
|
---|
86 | // 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);
|
---|
87 | // Double_t par[3+2*n_peaks] = {70,90,500000,15,50000,15,5000,15,500,15,800000,20,0};
|
---|
88 | f_peaks->SetParameters(par);
|
---|
89 | f_peaks->SetLineColor(2);
|
---|
90 |
|
---|
91 | TF1 *f_pedestal = new TF1("pedestal","[0]*TMath::Exp((-x-[1])/[2])",fit_start,fit_end);
|
---|
92 | f_pedestal->SetLineColor(2);
|
---|
93 | f_pedestal->SetLineWidth(1);
|
---|
94 |
|
---|
95 | TF1* f_peak[n_peaks];
|
---|
96 | for(int i=0; i<n_peaks; i++) {
|
---|
97 | f_peak[i] = new TF1("peak","gaus",fit_start,fit_end);
|
---|
98 | f_peak[i]->SetLineColor(2);
|
---|
99 | f_peak[i]->SetLineWidth(1);
|
---|
100 | }
|
---|
101 |
|
---|
102 | //-------------------------------------------
|
---|
103 | //Open the file
|
---|
104 | //-------------------------------------------
|
---|
105 | fits datafile(name);
|
---|
106 | if (!datafile)
|
---|
107 | {
|
---|
108 | cout << "Couldn't properly open the datafile." << endl;
|
---|
109 | return 1;
|
---|
110 | }
|
---|
111 | FOpenDataFile(datafile, data, data_offset, data_num, data_n, data_roi, data_px);
|
---|
112 |
|
---|
113 | //-------------------------------------------
|
---|
114 | //Find the baseline
|
---|
115 | //-------------------------------------------
|
---|
116 | float pedestal_mean, pedestal_rms;
|
---|
117 | // FPedestal(title_base, datafile, data, data_offset, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, pixelnr, pedestal_mean, pedestal_rms);
|
---|
118 | //Return values
|
---|
119 | //Value of maximal probability: -2.225 +- 4.3144
|
---|
120 | pedestal_mean = -2.225;
|
---|
121 | pedestal_rms = 4.3144;
|
---|
122 |
|
---|
123 | //-------------------------------------------
|
---|
124 | //Oscilloscope
|
---|
125 | //-------------------------------------------
|
---|
126 | // float threshold = pedestal_mean+pedestal_rms;
|
---|
127 | float threshold = 2.5;
|
---|
128 |
|
---|
129 | //***********************
|
---|
130 | FOscilloscope(datafile, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, threshold, pulseshape, spectrum);
|
---|
131 | canv->cd();
|
---|
132 | pulseshape->Draw();
|
---|
133 | canv->Modified();
|
---|
134 | canv->Update();
|
---|
135 |
|
---|
136 | canv_spect->cd();
|
---|
137 | spectrum->Draw();
|
---|
138 | canv_spect->Modified();
|
---|
139 | canv_spect->Update();
|
---|
140 | //***********************
|
---|
141 | fits datafile2("20110804_025.fits");
|
---|
142 | if (!datafile2)
|
---|
143 | {
|
---|
144 | cout << "Couldn't properly open the datafile2." << endl;
|
---|
145 | return 1;
|
---|
146 | }
|
---|
147 | FOpenDataFile(datafile2, data, data_offset, data_num, data_n, data_roi, data_px);
|
---|
148 | FOscilloscope(datafile2, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, threshold, pulseshape, spectrum);
|
---|
149 |
|
---|
150 | canv->cd();
|
---|
151 | pulseshape->Draw();
|
---|
152 | canv->Modified();
|
---|
153 | canv->Update();
|
---|
154 |
|
---|
155 | canv_spect->cd();
|
---|
156 | spectrum->Draw();
|
---|
157 | canv_spect->Modified();
|
---|
158 | canv_spect->Update();
|
---|
159 | //***********************
|
---|
160 | fits datafile3("20110804_026.fits");
|
---|
161 | if (!datafile3)
|
---|
162 | {
|
---|
163 | cout << "Couldn't properly open the datafile3." << endl;
|
---|
164 | return 1;
|
---|
165 | }
|
---|
166 | FOpenDataFile(datafile3, data, data_offset, data_num, data_n, data_roi, data_px);
|
---|
167 | FOscilloscope(datafile3, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, threshold, pulseshape, spectrum);
|
---|
168 |
|
---|
169 | canv->cd();
|
---|
170 | pulseshape->Draw();
|
---|
171 | canv->Modified();
|
---|
172 | canv->Update();
|
---|
173 |
|
---|
174 | canv_spect->cd();
|
---|
175 | spectrum->Draw();
|
---|
176 | canv_spect->Modified();
|
---|
177 | canv_spect->Update();
|
---|
178 | //***********************
|
---|
179 | fits datafile4("20110804_027.fits");
|
---|
180 | if (!datafile4)
|
---|
181 | {
|
---|
182 | cout << "Couldn't properly open the datafile4." << endl;
|
---|
183 | return 1;
|
---|
184 | }
|
---|
185 | FOpenDataFile(datafile4, data, data_offset, data_num, data_n, data_roi, data_px);
|
---|
186 | FOscilloscope(datafile4, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, threshold, pulseshape, spectrum);
|
---|
187 |
|
---|
188 | canv->cd();
|
---|
189 | pulseshape->Draw();
|
---|
190 | canv->Modified();
|
---|
191 | canv->Update();
|
---|
192 |
|
---|
193 | canv_spect->cd();
|
---|
194 | spectrum->Draw();
|
---|
195 | canv_spect->Modified();
|
---|
196 | canv_spect->Update();
|
---|
197 | //***********************
|
---|
198 | spectrum->Fit(f_peaks,"R+");
|
---|
199 | f_peaks->GetParameters(&par[0]);
|
---|
200 |
|
---|
201 | f_peak[0]->SetParameters(par[2],par[0],par[3]);
|
---|
202 | f_peak[1]->SetParameters(par[4],par[0]+par[1],par[5]);
|
---|
203 | f_peak[2]->SetParameters(par[6],par[0]+2*par[1],par[7]);
|
---|
204 | f_peak[3]->SetParameters(par[8],par[0]+3*par[1],par[9]);
|
---|
205 | // f_peak[4]->SetParameters(par[10],par[12],par[11]);
|
---|
206 | f_pedestal->SetParameters(par[10],par[11],par[12]);
|
---|
207 |
|
---|
208 | std::cout << "Crosstalk propability approx. " << ((par[4]*par[5])/(par[2]*par[3]+par[4]*par[5])) << std::endl;
|
---|
209 | //-------------------------------------------
|
---|
210 | //Draw the data
|
---|
211 | //-------------------------------------------
|
---|
212 | canv->cd();
|
---|
213 | pulseshape->Draw();
|
---|
214 | canv->Modified();
|
---|
215 | canv->Update();
|
---|
216 |
|
---|
217 | canv_spect->cd();
|
---|
218 | canv_spect->SetLogy();
|
---|
219 | spectrum->Draw();
|
---|
220 | f_peaks->Draw("same");
|
---|
221 |
|
---|
222 | for(int i=0; i<n_peaks; i++) {
|
---|
223 | f_peak[i]->Draw("same");
|
---|
224 | }
|
---|
225 |
|
---|
226 | f_pedestal->Draw("same");
|
---|
227 |
|
---|
228 | canv_spect->Modified();
|
---|
229 | canv_spect->Update();
|
---|
230 |
|
---|
231 | return 0;
|
---|
232 | }
|
---|