Rev 172 | Go to most recent revision | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed
Rev 172 | Rev 173 | ||
---|---|---|---|
1 | #include "../include/sipmscan.h" |
1 | #include "../include/sipmscan.h" |
2 | #include "../include/workstation.h" |
2 | #include "../include/workstation.h" |
3 | 3 | ||
4 | #include <stdio.h> |
4 | #include <stdio.h> |
5 | #include <stdlib.h> |
5 | #include <stdlib.h> |
6 | 6 | ||
7 | // Peak detection function |
7 | // Peak detection function |
8 | int npeaks = 20; |
8 | int npeaks = 20; |
9 | double FindPeaks(double *x, double *par) |
9 | double FindPeaks(double *x, double *par) |
10 | { |
10 | { |
11 | double result = 0; |
11 | double result = 0; |
12 | for(int i = 0; i < npeaks; i++) |
12 | for(int i = 0; i < npeaks; i++) |
13 | { |
13 | { |
14 | double norm = par[3*i]; |
14 | double norm = par[3*i]; |
15 | double mean = par[3*i+1]; |
15 | double mean = par[3*i+1]; |
16 | double sigma = par[3*i+2]; |
16 | double sigma = par[3*i+2]; |
17 | result += norm*TMath::Gaus(x[0], mean, sigma); |
17 | result += norm*TMath::Gaus(x[0], mean, sigma); |
18 | } |
18 | } |
19 | return result; |
19 | return result; |
20 | } |
20 | } |
21 | 21 | ||
22 | // File browser for selecting the dark run histogram |
22 | // File browser for selecting the dark run histogram |
23 | void TGAppMainFrame::SelectDarkHist() |
23 | void TGAppMainFrame::SelectDarkHist() |
24 | { |
24 | { |
25 | TGFileInfo file_info; |
25 | TGFileInfo file_info; |
26 | const char *filetypes[] = {"Histograms",histextall,0,0}; |
26 | const char *filetypes[] = {"Histograms",histextall,0,0}; |
27 | char *cTemp; |
27 | char *cTemp; |
28 | file_info.fFileTypes = filetypes; |
28 | file_info.fFileTypes = filetypes; |
29 | // cTemp = new char[1024]; |
29 | // cTemp = new char[1024]; |
30 | // sprintf(cTemp, "%s/results", rootdir); |
30 | // sprintf(cTemp, "%s/results", rootdir); |
31 | // file_info.fIniDir = StrDup(cTemp); |
31 | // file_info.fIniDir = StrDup(cTemp); |
32 | file_info.fIniDir = StrDup(currentAnalDir); |
32 | file_info.fIniDir = StrDup(currentAnalDir); |
33 | file_info.fMultipleSelection = kFALSE; |
33 | file_info.fMultipleSelection = kFALSE; |
34 | new TGFileDialog(gClient->GetDefaultRoot(), fMain, kFDOpen, &file_info); |
34 | new TGFileDialog(gClient->GetDefaultRoot(), fMain, kFDOpen, &file_info); |
35 | // delete[] cTemp; |
35 | // delete[] cTemp; |
36 | 36 | ||
37 | if(file_info.fFilename != NULL) |
37 | if(file_info.fFilename != NULL) |
38 | { |
38 | { |
39 | darkRun->widgetTE->SetText(file_info.fFilename); |
39 | darkRun->widgetTE->SetText(file_info.fFilename); |
40 | fileList->AddEntry(file_info.fFilename, fileList->GetNumberOfEntries()); |
40 | fileList->AddEntry(file_info.fFilename, fileList->GetNumberOfEntries()); |
41 | fileList->Layout(); |
41 | fileList->Layout(); |
42 | } |
42 | } |
43 | else |
43 | else |
44 | darkRun->widgetTE->SetText(""); |
44 | darkRun->widgetTE->SetText(""); |
45 | } |
45 | } |
46 | 46 | ||
47 | // Reset analysis defaults for the current analysis type (0 = Integrate spectrum, 1 = Relative PDE, 2 = Breakdown voltage, 3 = Surface scan, 4 = Timing analysis) |
47 | // Reset analysis defaults for the current analysis type (0 = Integrate spectrum, 1 = Relative PDE, 2 = Breakdown voltage, 3 = Surface scan, 4 = Timing analysis) |
48 | void TGAppMainFrame::AnalysisDefaults() |
48 | void TGAppMainFrame::AnalysisDefaults() |
49 | { |
49 | { |
50 | printf("AnalysisDefaults(): Current analysis tab = %d\n", analTab->GetCurrent()); |
50 | printf("AnalysisDefaults(): Current analysis tab = %d\n", analTab->GetCurrent()); |
51 | if(analTab->GetCurrent() == 0) // Integrate spectrum |
51 | if(analTab->GetCurrent() == 0) // Integrate spectrum |
52 | { |
52 | { |
53 | intSpect->widgetChBox[0]->SetState(kButtonUp); |
53 | intSpect->widgetChBox[0]->SetState(kButtonUp); |
54 | intSpect->widgetChBox[1]->SetState(kButtonUp); |
54 | intSpect->widgetChBox[1]->SetState(kButtonUp); |
55 | intSpect->widgetChBox[2]->SetState(kButtonDown); |
55 | intSpect->widgetChBox[2]->SetState(kButtonDown); |
56 | resol2d->widgetNE[0]->SetNumber(40); |
56 | resol2d->widgetNE[0]->SetNumber(40); |
57 | resol2d->widgetNE[1]->SetNumber(40); |
57 | resol2d->widgetNE[1]->SetNumber(40); |
58 | } |
58 | } |
59 | else if(analTab->GetCurrent() == 1) // Relative PDE |
59 | else if(analTab->GetCurrent() == 1) // Relative PDE |
60 | { |
60 | { |
61 | relPde->widgetChBox[0]->SetState(kButtonDown); |
61 | relPde->widgetChBox[0]->SetState(kButtonDown); |
- | 62 | relPde->widgetChBox[1]->SetState(kButtonUp); |
|
62 | midPeak->widgetChBox[0]->SetState(kButtonUp); |
63 | midPeak->widgetChBox[0]->SetState(kButtonUp); |
63 | zeroAngle->widgetNE[0]->SetNumber(0.00); |
64 | zeroAngle->widgetNE[0]->SetNumber(0.00); |
64 | } |
65 | } |
65 | else if(analTab->GetCurrent() == 2) // Breakdown voltage |
66 | else if(analTab->GetCurrent() == 2) // Breakdown voltage |
66 | { |
67 | { |
67 | minPeak->widgetNE[0]->SetNumber(2); |
68 | minPeak->widgetNE[0]->SetNumber(2); |
68 | minPeak->widgetNE[0]->SetNumber(1); |
69 | minPeak->widgetNE[0]->SetNumber(1); |
69 | } |
70 | } |
70 | else if(analTab->GetCurrent() == 3) // Surface scan |
71 | else if(analTab->GetCurrent() == 3) // Surface scan |
71 | { |
72 | { |
72 | surfScanOpt->widgetChBox[0]->SetState(kButtonDown); |
73 | surfScanOpt->widgetChBox[0]->SetState(kButtonDown); |
73 | surfScanOpt->widgetChBox[1]->SetState(kButtonUp); |
74 | surfScanOpt->widgetChBox[1]->SetState(kButtonUp); |
74 | resolSurf->widgetNE[0]->SetNumber(40); |
75 | resolSurf->widgetNE[0]->SetNumber(40); |
75 | resolSurf->widgetNE[1]->SetNumber(40); |
76 | resolSurf->widgetNE[1]->SetNumber(40); |
76 | } |
77 | } |
77 | } |
78 | } |
78 | 79 | ||
79 | // Analysis handle function |
80 | // Analysis handle function |
80 | void TGAppMainFrame::AnalysisHandle(int type) |
81 | void TGAppMainFrame::AnalysisHandle(int type) |
81 | { |
82 | { |
82 | TList *files; |
83 | TList *files; |
83 | bool createTab = true; |
84 | bool createTab = true; |
84 | int tabid = -1; |
85 | int tabid = -1; |
85 | int analtab = analTab->GetCurrent(); |
86 | int analtab = analTab->GetCurrent(); |
86 | 87 | ||
87 | int analtype; |
88 | int analtype; |
88 | if( (analtab == 0) && (!intSpect->widgetChBox[0]->IsDown() && !intSpect->widgetChBox[1]->IsDown()) ) |
89 | if( (analtab == 0) && (!intSpect->widgetChBox[0]->IsDown() && !intSpect->widgetChBox[1]->IsDown()) ) |
89 | analtype = 0; |
90 | analtype = 0; |
90 | else if( (analtab == 0) && (intSpect->widgetChBox[0]->IsDown() || intSpect->widgetChBox[1]->IsDown()) ) |
91 | else if( (analtab == 0) && (intSpect->widgetChBox[0]->IsDown() || intSpect->widgetChBox[1]->IsDown()) ) |
91 | analtype = 1; |
92 | analtype = 1; |
92 | else if( analtab == 1 ) |
93 | else if( analtab == 1 ) |
93 | analtype = 2; |
94 | analtype = 2; |
94 | else if( analtab == 2 ) |
95 | else if( analtab == 2 ) |
95 | analtype = 3; |
96 | analtype = 3; |
96 | else if( analtab == 3 ) |
97 | else if( analtab == 3 ) |
97 | analtype = 4; |
98 | analtype = 4; |
- | 99 | ||
- | 100 | // Save analysis settings any time we run a new analysis |
|
- | 101 | SaveAnalSettings(); |
|
98 | 102 | ||
99 | // Only integrate spectrum or make relative PDE |
103 | // Only integrate spectrum or make relative PDE |
100 | if(type == 0) |
104 | if(type == 0) |
101 | { |
105 | { |
102 | files = new TList(); |
106 | files = new TList(); |
103 | fileList->GetSelectedEntries(files); |
107 | fileList->GetSelectedEntries(files); |
104 | 108 | ||
105 | if( analtype == 0 ) |
109 | if( analtype == 0 ) |
106 | IntegSpectrum(files, 0, 0); |
110 | IntegSpectrum(files, 0, 0); |
107 | 111 | ||
108 | if( intSpect->widgetChBox[0]->IsDown() ) |
112 | if( intSpect->widgetChBox[0]->IsDown() ) |
109 | IntegSpectrum(files, 1, 0); |
113 | IntegSpectrum(files, 1, 0); |
110 | 114 | ||
111 | if( intSpect->widgetChBox[1]->IsDown() ) |
115 | if( intSpect->widgetChBox[1]->IsDown() ) |
112 | IntegSpectrum(files, 2, 0); |
116 | IntegSpectrum(files, 2, 0); |
113 | 117 | ||
114 | if( analtype == 2 ) |
118 | if( analtype == 2 ) |
115 | PhotonMu(files, 0); |
119 | PhotonMu(files, 0); |
116 | 120 | ||
117 | if( analtype == 3 ) |
121 | if( analtype == 3 ) |
118 | BreakdownVolt(files, 0); |
122 | BreakdownVolt(files, 0); |
119 | 123 | ||
120 | if( analtype == 4 ) |
124 | if( analtype == 4 ) |
121 | SurfaceScan(files, 0); |
125 | SurfaceScan(files, 0); |
122 | } |
126 | } |
123 | // Integrate spectrum or make relative PDE and open edit window |
127 | // Integrate spectrum or make relative PDE and open edit window |
124 | else if(type == 1) |
128 | else if(type == 1) |
125 | { |
129 | { |
126 | files = new TList(); |
130 | files = new TList(); |
127 | fileList->GetSelectedEntries(files); |
131 | fileList->GetSelectedEntries(files); |
128 | 132 | ||
129 | // Prepare a new analysis edit tab, if we want to edit plots |
133 | // Prepare a new analysis edit tab, if we want to edit plots |
130 | for(int i = 0; i < fTab->GetNumberOfTabs(); i++) |
134 | for(int i = 0; i < fTab->GetNumberOfTabs(); i++) |
131 | { |
135 | { |
132 | if(strcmp("Analysis edit", fTab->GetTabTab(i)->GetString() ) == 0) |
136 | if(strcmp("Analysis edit", fTab->GetTabTab(i)->GetString() ) == 0) |
133 | { |
137 | { |
134 | createTab = false; |
138 | createTab = false; |
135 | tabid = i; |
139 | tabid = i; |
136 | } |
140 | } |
137 | 141 | ||
138 | if(DBGSIG) printf("AnalysisHandle(): Name of tab = %s\n", fTab->GetTabTab(i)->GetString() ); |
142 | if(DBGSIG) printf("AnalysisHandle(): Name of tab = %s\n", fTab->GetTabTab(i)->GetString() ); |
139 | } |
143 | } |
140 | 144 | ||
141 | if(files->GetSize() > 0) |
145 | if(files->GetSize() > 0) |
142 | { |
146 | { |
143 | TempAnalysisTab(fTab, createTab, &tabid, analtype); |
147 | TempAnalysisTab(fTab, createTab, &tabid, analtype); |
144 | 148 | ||
145 | // Integrate spectra |
149 | // Integrate spectra |
146 | if( analtype == 0 ) |
150 | if( analtype == 0 ) |
147 | IntegSpectrum(files, 0, 1); |
151 | IntegSpectrum(files, 0, 1); |
148 | 152 | ||
149 | if( intSpect->widgetChBox[0]->IsDown() ) |
153 | if( intSpect->widgetChBox[0]->IsDown() ) |
150 | IntegSpectrum(files, 1, 1); |
154 | IntegSpectrum(files, 1, 1); |
151 | 155 | ||
152 | if( intSpect->widgetChBox[1]->IsDown() ) |
156 | if( intSpect->widgetChBox[1]->IsDown() ) |
153 | IntegSpectrum(files, 2, 1); |
157 | IntegSpectrum(files, 2, 1); |
154 | 158 | ||
155 | if( analtype == 2 ) |
159 | if( analtype == 2 ) |
156 | PhotonMu(files, 1); |
160 | PhotonMu(files, 1); |
157 | 161 | ||
158 | if( analtype == 3 ) |
162 | if( analtype == 3 ) |
159 | BreakdownVolt(files, 1); |
163 | BreakdownVolt(files, 1); |
160 | 164 | ||
161 | if( analtype == 4 ) |
165 | if( analtype == 4 ) |
162 | SurfaceScan(files, 1); |
166 | SurfaceScan(files, 1); |
163 | 167 | ||
164 | fTab->SetTab(tabid); |
168 | fTab->SetTab(tabid); |
165 | } |
169 | } |
166 | } |
170 | } |
167 | 171 | ||
168 | delete files; |
172 | delete files; |
169 | } |
173 | } |
170 | 174 | ||
171 | // Analysis functions ---------------------------------- |
175 | // Analysis functions ---------------------------------- |
172 | void TGAppMainFrame::IntegSpectrum(TList *files, int direction, int edit) |
176 | void TGAppMainFrame::IntegSpectrum(TList *files, int direction, int edit) |
173 | { |
177 | { |
174 | unsigned int nrfiles = files->GetSize(); |
178 | unsigned int nrfiles = files->GetSize(); |
175 | char ctemp[1024]; |
179 | char ctemp[1024]; |
176 | char exportname[1024]; |
180 | char exportname[1024]; |
177 | int j, k = 0, m = 0; |
181 | int j, k = 0, m = 0; |
178 | 182 | ||
179 | TTree *header_data, *meas_data; |
183 | TTree *header_data, *meas_data; |
180 | double *integralCount, *integralAcc; |
184 | double *integralCount, *integralAcc; |
181 | integralCount = new double[nrfiles]; |
185 | integralCount = new double[nrfiles]; |
182 | integralAcc = new double[nrfiles]; |
186 | integralAcc = new double[nrfiles]; |
183 | double *surfxy, *surfz; |
187 | double *surfxy, *surfz; |
184 | surfxy = new double[nrfiles]; |
188 | surfxy = new double[nrfiles]; |
185 | surfz = new double[nrfiles]; |
189 | surfz = new double[nrfiles]; |
186 | double minInteg, maxInteg; |
190 | double minInteg, maxInteg; |
187 | bool norm = intSpect->widgetChBox[2]->IsDown(); |
191 | bool norm = intSpect->widgetChBox[2]->IsDown(); |
188 | double curzval; |
192 | double curzval; |
189 | bool edge2d = false; |
193 | bool edge2d = false; |
190 | TCanvas *gCanvas; |
194 | TCanvas *gCanvas; |
191 | 195 | ||
192 | float progVal = 0; |
196 | float progVal = 0; |
193 | analysisProgress->widgetPB->SetPosition(progVal); |
197 | analysisProgress->widgetPB->SetPosition(progVal); |
194 | gVirtualX->Update(1); |
198 | gVirtualX->Update(1); |
195 | 199 | ||
196 | // Zero the integral count and accumulated vaules |
200 | // Zero the integral count and accumulated vaules |
197 | for(int i = 0; i < (int)nrfiles; i++) {integralCount[i] = 0; integralAcc[i] = 0; } |
201 | for(int i = 0; i < (int)nrfiles; i++) {integralCount[i] = 0; integralAcc[i] = 0; } |
198 | 202 | ||
199 | // Start if we select at least one file |
203 | // Start if we select at least one file |
200 | if(nrfiles > 0) |
204 | if(nrfiles > 0) |
201 | { |
205 | { |
202 | for(int i = 0; i < (int)nrfiles; i++) |
206 | for(int i = 0; i < (int)nrfiles; i++) |
203 | { |
207 | { |
204 | if(files->At(i)) |
208 | if(files->At(i)) |
205 | { |
209 | { |
206 | // Read the X,Y and Z positions from header and ADC and TDC values from the measurements |
210 | // Read the X,Y and Z positions from header and ADC and TDC values from the measurements |
207 | sprintf(ctemp, "%s", files->At(i)->GetTitle()); |
211 | sprintf(ctemp, "%s", files->At(i)->GetTitle()); |
208 | inroot = new TFile(ctemp, "READ"); |
212 | inroot = new TFile(ctemp, "READ"); |
209 | 213 | ||
210 | inroot->GetObject("header_data", header_data); |
214 | inroot->GetObject("header_data", header_data); |
211 | inroot->GetObject("meas_data", meas_data); |
215 | inroot->GetObject("meas_data", meas_data); |
212 | 216 | ||
213 | header_data->SetBranchAddress("xpos", &evtheader.xpos); |
217 | header_data->SetBranchAddress("xpos", &evtheader.xpos); |
214 | header_data->GetEntry(0); |
218 | header_data->GetEntry(0); |
215 | header_data->SetBranchAddress("ypos", &evtheader.ypos); |
219 | header_data->SetBranchAddress("ypos", &evtheader.ypos); |
216 | header_data->GetEntry(0); |
220 | header_data->GetEntry(0); |
217 | header_data->SetBranchAddress("zpos", &evtheader.zpos); |
221 | header_data->SetBranchAddress("zpos", &evtheader.zpos); |
218 | header_data->GetEntry(0); |
222 | header_data->GetEntry(0); |
219 | 223 | ||
220 | char rdc[256]; |
224 | char rdc[256]; |
221 | j = selectCh->widgetNE[0]->GetNumber(); |
225 | j = selectCh->widgetNE[0]->GetNumber(); |
222 | double rangetdc[2]; |
226 | double rangetdc[2]; |
223 | rangetdc[0] = tdcRange->widgetNE[0]->GetNumber(); |
227 | rangetdc[0] = tdcRange->widgetNE[0]->GetNumber(); |
224 | rangetdc[1] = tdcRange->widgetNE[1]->GetNumber(); |
228 | rangetdc[1] = tdcRange->widgetNE[1]->GetNumber(); |
225 | 229 | ||
226 | k = 0; |
230 | k = 0; |
227 | m = 0; |
231 | m = 0; |
228 | 232 | ||
229 | // Reading the data |
233 | // Reading the data |
230 | for(int e = 0; e < meas_data->GetEntries(); e++) |
234 | for(int e = 0; e < meas_data->GetEntries(); e++) |
231 | { |
235 | { |
232 | sprintf(rdc, "ADC%d", j); |
236 | sprintf(rdc, "ADC%d", j); |
233 | meas_data->SetBranchAddress(rdc, &evtdata.adcdata[j]); |
237 | meas_data->SetBranchAddress(rdc, &evtdata.adcdata[j]); |
234 | meas_data->GetEntry(e); |
238 | meas_data->GetEntry(e); |
235 | 239 | ||
236 | sprintf(rdc, "TDC%d", j); |
240 | sprintf(rdc, "TDC%d", j); |
237 | meas_data->SetBranchAddress(rdc, &evtdata.tdcdata[j]); |
241 | meas_data->SetBranchAddress(rdc, &evtdata.tdcdata[j]); |
238 | meas_data->GetEntry(e); |
242 | meas_data->GetEntry(e); |
239 | 243 | ||
240 | // Use data point only if it is inside the TDC window |
244 | // Use data point only if it is inside the TDC window |
241 | if( ((double)evtdata.tdcdata[j]/tdctimeconversion >= rangetdc[0]) && ((double)evtdata.tdcdata[j]/tdctimeconversion <= rangetdc[1]) ) |
245 | if( ((double)evtdata.tdcdata[j]/tdctimeconversion >= rangetdc[0]) && ((double)evtdata.tdcdata[j]/tdctimeconversion <= rangetdc[1]) ) |
242 | { |
246 | { |
243 | k++; |
247 | k++; |
244 | m += evtdata.adcdata[j]; |
248 | m += evtdata.adcdata[j]; |
245 | } |
249 | } |
246 | } |
250 | } |
247 | 251 | ||
248 | // X, Y and Z values from each file (table units or microns) |
252 | // X, Y and Z values from each file (table units or microns) |
249 | if(posUnitsPlot->widgetCB->GetSelected() == 0) |
253 | if(posUnitsPlot->widgetCB->GetSelected() == 0) |
250 | { |
254 | { |
251 | if(direction == 1) |
255 | if(direction == 1) |
252 | surfxy[i] = (double)(evtheader.xpos); |
256 | surfxy[i] = (double)(evtheader.xpos); |
253 | else if(direction == 2) |
257 | else if(direction == 2) |
254 | surfxy[i] = (double)(evtheader.ypos); |
258 | surfxy[i] = (double)(evtheader.ypos); |
255 | surfz[i] = (double)(evtheader.zpos); |
259 | surfz[i] = (double)(evtheader.zpos); |
256 | } |
260 | } |
257 | else if(posUnitsPlot->widgetCB->GetSelected() == 1) |
261 | else if(posUnitsPlot->widgetCB->GetSelected() == 1) |
258 | { |
262 | { |
259 | if(direction == 1) |
263 | if(direction == 1) |
260 | surfxy[i] = (double)(evtheader.xpos*lenconversion); |
264 | surfxy[i] = (double)(evtheader.xpos*lenconversion); |
261 | else if(direction == 2) |
265 | else if(direction == 2) |
262 | surfxy[i] = (double)(evtheader.ypos*lenconversion); |
266 | surfxy[i] = (double)(evtheader.ypos*lenconversion); |
263 | surfz[i] = (double)(evtheader.zpos*lenconversion); |
267 | surfz[i] = (double)(evtheader.zpos*lenconversion); |
264 | } |
268 | } |
265 | 269 | ||
266 | // Check if we have different Z values: if no, just make the edge plots; if yes, save edge plots and make a 2d edge plot |
270 | // Check if we have different Z values: if no, just make the edge plots; if yes, save edge plots and make a 2d edge plot |
267 | if(i == 0) curzval = surfz[i]; |
271 | if(i == 0) curzval = surfz[i]; |
268 | else |
272 | else |
269 | { |
273 | { |
270 | if(surfz[i] != curzval) |
274 | if(surfz[i] != curzval) |
271 | edge2d = true; |
275 | edge2d = true; |
272 | } |
276 | } |
273 | 277 | ||
274 | // Print the calculated integral, if no X or Y direction edge plots are enabled; otherwise, just save for later plotting |
278 | // Print the calculated integral, if no X or Y direction edge plots are enabled; otherwise, just save for later plotting |
275 | if(direction == 0) |
279 | if(direction == 0) |
276 | { |
280 | { |
277 | if(norm) |
281 | if(norm) |
278 | { |
282 | { |
279 | integralCount[i] += ((double)m)/((double)k); |
283 | integralCount[i] += ((double)m)/((double)k); |
280 | printf("IntegSpectrum(): %s: The integral is: %lf\n", ctemp, integralCount[i]); |
284 | printf("IntegSpectrum(): %s: The integral is: %lf\n", ctemp, integralCount[i]); |
281 | } |
285 | } |
282 | else |
286 | else |
283 | { |
287 | { |
284 | integralCount[i] += m; |
288 | integralCount[i] += m; |
285 | printf("IntegSpectrum(): %s: The integral is: %d\n", ctemp, (int)integralCount[i]); |
289 | printf("IntegSpectrum(): %s: The integral is: %d\n", ctemp, (int)integralCount[i]); |
286 | } |
290 | } |
287 | } |
291 | } |
288 | else |
292 | else |
289 | { |
293 | { |
290 | if(norm) |
294 | if(norm) |
291 | integralCount[i] += ((double)m)/((double)k); |
295 | integralCount[i] += ((double)m)/((double)k); |
292 | else |
296 | else |
293 | integralCount[i] += m; |
297 | integralCount[i] += m; |
294 | } |
298 | } |
295 | 299 | ||
296 | inroot->Close(); |
300 | inroot->Close(); |
297 | delete inroot; |
301 | delete inroot; |
298 | } |
302 | } |
299 | 303 | ||
300 | // Update the progress bar |
304 | // Update the progress bar |
301 | progVal = (float)(75.00/nrfiles)*i; |
305 | progVal = (float)(75.00/nrfiles)*i; |
302 | analysisProgress->widgetPB->SetPosition(progVal); |
306 | analysisProgress->widgetPB->SetPosition(progVal); |
303 | gVirtualX->Update(1); |
307 | gVirtualX->Update(1); |
304 | } |
308 | } |
305 | 309 | ||
306 | printf("IntegSpectrum(): %d files were selected.\n", nrfiles); |
310 | printf("IntegSpectrum(): %d files were selected.\n", nrfiles); |
307 | 311 | ||
308 | // If only an integral is needed, do not plot and exit here |
312 | // If only an integral is needed, do not plot and exit here |
309 | if( direction == 0 ) |
313 | if( direction == 0 ) |
310 | { |
314 | { |
311 | delete[] integralCount; |
315 | delete[] integralCount; |
312 | delete[] surfxy; |
316 | delete[] surfxy; |
313 | delete[] surfz; |
317 | delete[] surfz; |
314 | return; |
318 | return; |
315 | } |
319 | } |
316 | 320 | ||
317 | // Current z value and the accumulated counter |
321 | // Current z value and the accumulated counter |
318 | curzval = surfz[0]; |
322 | curzval = surfz[0]; |
319 | j = 0; |
323 | j = 0; |
320 | int acc = 0; |
324 | int acc = 0; |
321 | int zb; |
325 | int zb; |
322 | for(int i = 0; i <= (int)nrfiles; i++) |
326 | for(int i = 0; i <= (int)nrfiles; i++) |
323 | { |
327 | { |
324 | // Collect the accumulated integral in order to produce a PDF from a CDF |
328 | // Collect the accumulated integral in order to produce a PDF from a CDF |
325 | // While we are at the same Z value, save under one set |
329 | // While we are at the same Z value, save under one set |
326 | if( (surfz[i] == curzval) && (acc != nrfiles) ) |
330 | if( (surfz[i] == curzval) && (acc != nrfiles) ) |
327 | { |
331 | { |
328 | integralAcc[j] = integralCount[i]; |
332 | integralAcc[j] = integralCount[i]; |
329 | if(DBGSIG) printf("IntegSpectrum(): Integral check 1 (i=%d,j=%d,z=%.2lf): %lf\t%lf\n", i, j, surfz[i], integralCount[i], integralAcc[j]); |
333 | if(DBGSIG) printf("IntegSpectrum(): Integral check 1 (i=%d,j=%d,z=%.2lf): %lf\t%lf\n", i, j, surfz[i], integralCount[i], integralAcc[j]); |
330 | j++; |
334 | j++; |
331 | acc++; |
335 | acc++; |
332 | } |
336 | } |
333 | // When we switch to a new set of Z values and at the end, we must save the previous ones to make 1D edge plots |
337 | // When we switch to a new set of Z values and at the end, we must save the previous ones to make 1D edge plots |
334 | else |
338 | else |
335 | { |
339 | { |
336 | // Find minimal and maximal integral values to subtract the offset and normate PDF to 1 |
340 | // Find minimal and maximal integral values to subtract the offset and normate PDF to 1 |
337 | NormateSet(i, j, &minInteg, &maxInteg, integralCount, integralAcc); |
341 | NormateSet(i, j, &minInteg, &maxInteg, integralCount, integralAcc); |
338 | 342 | ||
339 | if(acc != nrfiles) |
343 | if(acc != nrfiles) |
340 | { |
344 | { |
341 | curzval = surfz[i]; |
345 | curzval = surfz[i]; |
342 | // PDF and CDF plot |
346 | // PDF and CDF plot |
343 | PlotEdgeDistribution(files, i, j, &minInteg, &maxInteg, surfxy, integralAcc, direction, edge2d, edit); |
347 | PlotEdgeDistribution(files, i, j, &minInteg, &maxInteg, surfxy, integralAcc, direction, edge2d, edit); |
344 | i--; |
348 | i--; |
345 | j = 0; |
349 | j = 0; |
346 | } |
350 | } |
347 | else |
351 | else |
348 | { |
352 | { |
349 | // PDF and CDF plot |
353 | // PDF and CDF plot |
350 | PlotEdgeDistribution(files, i, j, &minInteg, &maxInteg, surfxy, integralAcc, direction, edge2d, edit); |
354 | PlotEdgeDistribution(files, i, j, &minInteg, &maxInteg, surfxy, integralAcc, direction, edge2d, edit); |
351 | i--; |
355 | i--; |
352 | break; |
356 | break; |
353 | } |
357 | } |
354 | } |
358 | } |
355 | 359 | ||
356 | // Update the progress bar |
360 | // Update the progress bar |
357 | progVal = (float)(15.00/nrfiles)*i+75.00; |
361 | progVal = (float)(15.00/nrfiles)*i+75.00; |
358 | analysisProgress->widgetPB->SetPosition(progVal); |
362 | analysisProgress->widgetPB->SetPosition(progVal); |
359 | gVirtualX->Update(1); |
363 | gVirtualX->Update(1); |
360 | } |
364 | } |
361 | 365 | ||
362 | // Make the 2D edge plot |
366 | // Make the 2D edge plot |
363 | if(edge2d) |
367 | if(edge2d) |
364 | { |
368 | { |
365 | if(edit == 0) |
369 | if(edit == 0) |
366 | gCanvas = new TCanvas("canv","canv",900,900); |
370 | gCanvas = new TCanvas("canv","canv",900,900); |
367 | else |
371 | else |
368 | gCanvas = tempAnalysisCanvas->GetCanvas(); |
372 | gCanvas = tempAnalysisCanvas->GetCanvas(); |
369 | 373 | ||
370 | double range[4]; |
374 | double range[4]; |
371 | TGraph2D *gScan2D; |
375 | TGraph2D *gScan2D; |
372 | gScan2D = new TGraph2D(); |
376 | gScan2D = new TGraph2D(); |
373 | gScan2D->SetName("edgescan"); |
- | |
374 | range[0] = TMath::MinElement(nrfiles, surfxy); |
377 | range[0] = TMath::MinElement(nrfiles, surfxy); |
375 | range[1] = TMath::MaxElement(nrfiles, surfxy); |
378 | range[1] = TMath::MaxElement(nrfiles, surfxy); |
376 | range[2] = TMath::MinElement(nrfiles, surfz); |
379 | range[2] = TMath::MinElement(nrfiles, surfz); |
377 | range[3] = TMath::MaxElement(nrfiles, surfz); |
380 | range[3] = TMath::MaxElement(nrfiles, surfz); |
378 | 381 | ||
379 | for(int i = 0; i < nrfiles; i++) |
382 | for(int i = 0; i < nrfiles; i++) |
380 | { |
383 | { |
381 | if(DBGSIG) |
384 | if(DBGSIG) |
382 | printf("IntegSpectrum(): %.2lf\t%.2lf\t%lf\n", surfxy[i], surfz[i], integralCount[i]); |
385 | printf("IntegSpectrum(): %.2lf\t%.2lf\t%lf\n", surfxy[i], surfz[i], integralCount[i]); |
383 | gScan2D->SetPoint(i, surfxy[i], surfz[i], integralCount[i]); |
386 | gScan2D->SetPoint(i, surfxy[i], surfz[i], integralCount[i]); |
384 | 387 | ||
385 | // Update the progress bar |
388 | // Update the progress bar |
386 | progVal = (float)(9.00/nrfiles)*i+90.00; |
389 | progVal = (float)(9.00/nrfiles)*i+90.00; |
387 | analysisProgress->widgetPB->SetPosition(progVal); |
390 | analysisProgress->widgetPB->SetPosition(progVal); |
388 | gVirtualX->Update(1); |
391 | gVirtualX->Update(1); |
389 | } |
392 | } |
390 | 393 | ||
391 | gCanvas->cd(); |
394 | gCanvas->cd(); |
392 | gStyle->SetPalette(1); |
395 | gStyle->SetPalette(1); |
393 | gCanvas->SetLeftMargin(0.15); |
396 | gCanvas->SetLeftMargin(0.15); |
394 | gCanvas->SetRightMargin(0.126); |
397 | gCanvas->SetRightMargin(0.126); |
395 | gCanvas->SetTopMargin(0.077); |
398 | gCanvas->SetTopMargin(0.077); |
396 | gScan2D->Draw("COLZ"); |
399 | gScan2D->Draw("COLZ"); |
397 | 400 | ||
398 | gCanvas->Modified(); |
401 | gCanvas->Modified(); |
399 | gCanvas->Update(); |
402 | gCanvas->Update(); |
400 | 403 | ||
401 | gScan2D->SetNpx((int)resol2d->widgetNE[0]->GetNumber()); |
404 | gScan2D->SetNpx((int)resol2d->widgetNE[0]->GetNumber()); |
402 | gScan2D->SetNpy((int)resol2d->widgetNE[1]->GetNumber()); |
405 | gScan2D->SetNpy((int)resol2d->widgetNE[1]->GetNumber()); |
403 | 406 | ||
404 | gCanvas->Modified(); |
407 | gCanvas->Modified(); |
405 | gCanvas->Update(); |
408 | gCanvas->Update(); |
406 | 409 | ||
407 | if(direction == 1) |
410 | if(direction == 1) |
408 | gScan2D->GetXaxis()->SetTitle("X [#mum]"); |
411 | gScan2D->GetXaxis()->SetTitle("X [#mum]"); |
409 | else if(direction == 2) |
412 | else if(direction == 2) |
410 | gScan2D->GetXaxis()->SetTitle("Y [#mum]"); |
413 | gScan2D->GetXaxis()->SetTitle("Y [#mum]"); |
411 | 414 | ||
412 | 415 | ||
413 | gScan2D->GetXaxis()->SetTitleOffset(1.3); |
416 | gScan2D->GetXaxis()->SetTitleOffset(1.3); |
414 | gScan2D->GetXaxis()->CenterTitle(kTRUE); |
417 | gScan2D->GetXaxis()->CenterTitle(kTRUE); |
415 | gScan2D->GetXaxis()->SetLabelSize(0.027); |
418 | gScan2D->GetXaxis()->SetLabelSize(0.027); |
416 | gScan2D->GetXaxis()->SetLabelOffset(0.02); |
419 | gScan2D->GetXaxis()->SetLabelOffset(0.02); |
417 | gScan2D->GetXaxis()->SetRangeUser(range[0], range[1]); |
420 | gScan2D->GetXaxis()->SetRangeUser(range[0], range[1]); |
418 | gScan2D->GetXaxis()->SetNoExponent(kTRUE); |
421 | gScan2D->GetXaxis()->SetNoExponent(kTRUE); |
419 | gScan2D->GetYaxis()->SetTitleOffset(1.9); |
422 | gScan2D->GetYaxis()->SetTitleOffset(1.9); |
420 | gScan2D->GetYaxis()->CenterTitle(kTRUE); |
423 | gScan2D->GetYaxis()->CenterTitle(kTRUE); |
421 | gScan2D->GetYaxis()->SetLabelSize(0.027); |
424 | gScan2D->GetYaxis()->SetLabelSize(0.027); |
422 | gScan2D->GetYaxis()->SetLabelOffset(0.02); |
425 | gScan2D->GetYaxis()->SetLabelOffset(0.02); |
423 | gScan2D->GetYaxis()->SetRangeUser(range[2], range[3]); |
426 | gScan2D->GetYaxis()->SetRangeUser(range[2], range[3]); |
424 | gScan2D->GetYaxis()->SetNoExponent(kTRUE); |
427 | gScan2D->GetYaxis()->SetNoExponent(kTRUE); |
425 | 428 | ||
426 | /* TGaxis *yax = (TGaxis*)gScan2D->GetYaxis(); |
429 | /* TGaxis *yax = (TGaxis*)gScan2D->GetYaxis(); |
427 | yax->SetMaxDigits(4);*/ |
430 | yax->SetMaxDigits(4);*/ |
428 | 431 | ||
429 | if(!cleanPlots) |
432 | if(!cleanPlots) |
430 | { |
433 | { |
431 | if(direction == 1) |
434 | if(direction == 1) |
432 | { |
435 | { |
433 | if(posUnitsPlot->widgetCB->GetSelected() == 0) |
436 | if(posUnitsPlot->widgetCB->GetSelected() == 0) |
434 | gScan2D->SetTitle("Laser focal point;X [table units];Z [table units]"); |
437 | gScan2D->SetTitle("Laser focal point;X [table units];Z [table units]"); |
435 | else if(posUnitsPlot->widgetCB->GetSelected() == 1) |
438 | else if(posUnitsPlot->widgetCB->GetSelected() == 1) |
436 | gScan2D->SetTitle("Laser focal point;X [#mum];Z [#mum]"); |
439 | gScan2D->SetTitle("Laser focal point;X [#mum];Z [#mum]"); |
437 | } |
440 | } |
438 | else if(direction == 2) |
441 | else if(direction == 2) |
439 | { |
442 | { |
440 | if(posUnitsPlot->widgetCB->GetSelected() == 0) |
443 | if(posUnitsPlot->widgetCB->GetSelected() == 0) |
441 | gScan2D->SetTitle("Laser focal point;Y [table units];Z [table units]"); |
444 | gScan2D->SetTitle("Laser focal point;Y [table units];Z [table units]"); |
442 | else if(posUnitsPlot->widgetCB->GetSelected() == 1) |
445 | else if(posUnitsPlot->widgetCB->GetSelected() == 1) |
443 | gScan2D->SetTitle("Laser focal point;Y [#mum];Z [#mum]"); |
446 | gScan2D->SetTitle("Laser focal point;Y [#mum];Z [#mum]"); |
444 | } |
447 | } |
445 | } |
448 | } |
446 | else |
449 | else |
447 | { |
450 | { |
448 | if(direction == 1) |
451 | if(direction == 1) |
449 | { |
452 | { |
450 | if(posUnitsPlot->widgetCB->GetSelected() == 0) |
453 | if(posUnitsPlot->widgetCB->GetSelected() == 0) |
451 | gScan2D->SetTitle(";X [table units];Z [table units]"); |
454 | gScan2D->SetTitle(";X [table units];Z [table units]"); |
452 | else if(posUnitsPlot->widgetCB->GetSelected() == 1) |
455 | else if(posUnitsPlot->widgetCB->GetSelected() == 1) |
453 | gScan2D->SetTitle(";X [#mum];Z [#mum]"); |
456 | gScan2D->SetTitle(";X [#mum];Z [#mum]"); |
454 | } |
457 | } |
455 | else if(direction == 2) |
458 | else if(direction == 2) |
456 | { |
459 | { |
457 | if(posUnitsPlot->widgetCB->GetSelected() == 0) |
460 | if(posUnitsPlot->widgetCB->GetSelected() == 0) |
458 | gScan2D->SetTitle(";Y [table units];Z [table units]"); |
461 | gScan2D->SetTitle(";Y [table units];Z [table units]"); |
459 | else if(posUnitsPlot->widgetCB->GetSelected() == 1) |
462 | else if(posUnitsPlot->widgetCB->GetSelected() == 1) |
460 | gScan2D->SetTitle(";Y [#mum];Z [#mum]"); |
463 | gScan2D->SetTitle(";Y [#mum];Z [#mum]"); |
461 | } |
464 | } |
462 | } |
465 | } |
463 | 466 | ||
464 | gCanvas->Modified(); |
467 | gCanvas->Modified(); |
465 | gCanvas->Update(); |
468 | gCanvas->Update(); |
466 | 469 | ||
467 | remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp); |
470 | remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp); |
468 | sprintf(exportname, "%s", ctemp); |
471 | sprintf(exportname, "%s", ctemp); |
469 | remove_from_last(exportname, '_', ctemp); |
472 | remove_from_last(exportname, '_', ctemp); |
470 | if(direction == 1) |
473 | if(direction == 1) |
471 | sprintf(exportname, "%s_xdir_focalpoint.pdf", ctemp); |
474 | sprintf(exportname, "%s_xdir_focalpoint.pdf", ctemp); |
472 | else if(direction == 2) |
475 | else if(direction == 2) |
473 | sprintf(exportname, "%s_ydir_focalpoint.pdf", ctemp); |
476 | sprintf(exportname, "%s_ydir_focalpoint.pdf", ctemp); |
474 | gCanvas->SaveAs(exportname); |
477 | gCanvas->SaveAs(exportname); |
475 | 478 | ||
476 | // Update the progress bar |
479 | // Update the progress bar |
477 | analysisProgress->widgetPB->SetPosition(100.0); |
480 | analysisProgress->widgetPB->SetPosition(100.0); |
478 | gVirtualX->Update(1); |
481 | gVirtualX->Update(1); |
479 | 482 | ||
480 | if(edit == 0) |
483 | if(edit == 0) |
481 | { |
484 | { |
482 | delete gScan2D; |
485 | delete gScan2D; |
483 | delete gCanvas; |
486 | delete gCanvas; |
484 | } |
487 | } |
485 | } |
488 | } |
486 | else |
489 | else |
487 | { |
490 | { |
488 | // Update the progress bar |
491 | // Update the progress bar |
489 | analysisProgress->widgetPB->SetPosition(100.0); |
492 | analysisProgress->widgetPB->SetPosition(100.0); |
490 | gVirtualX->Update(1); |
493 | gVirtualX->Update(1); |
491 | } |
494 | } |
492 | } |
495 | } |
493 | } |
496 | } |
494 | 497 | ||
495 | void TGAppMainFrame::PlotEdgeDistribution(TList *files, int filenr, int points, double *min, double *max, double *xy, double *integAcc, int axis, bool edge2d, int edit) |
498 | void TGAppMainFrame::PlotEdgeDistribution(TList *files, int filenr, int points, double *min, double *max, double *xy, double *integAcc, int axis, bool edge2d, int edit) |
496 | { |
499 | { |
497 | TGraph *gScan[2]; |
500 | TGraph *gScan[2]; |
498 | int pdfmax = -1; |
501 | int pdfmax = -1; |
499 | int count = 0; |
502 | int count = 0; |
500 | char ctemp[1024]; |
503 | char ctemp[1024]; |
501 | char exportname[1024]; |
504 | char exportname[1024]; |
502 | TCanvas *gCanvas; |
505 | TCanvas *gCanvas; |
503 | 506 | ||
504 | // Prepare the CDF plot |
507 | // Prepare the CDF plot |
505 | gScan[1] = new TGraph(); |
508 | gScan[1] = new TGraph(); |
506 | for(int i = 0; i < points; i++) |
509 | for(int i = 0; i < points; i++) |
507 | { |
510 | { |
508 | count = filenr - points + i; |
511 | count = filenr - points + i; |
509 | gScan[1]->SetPoint(i, (double)xy[count], (double)integAcc[i]/(*max)); |
512 | gScan[1]->SetPoint(i, (double)xy[count], (double)integAcc[i]/(*max)); |
510 | if(DBGSIG) printf("PlotEdgeDistribution(): CDF %d: %lf, %lf\n", i, (double)xy[count], (double)integAcc[i]/(*max)); |
513 | if(DBGSIG) printf("PlotEdgeDistribution(): CDF %d: %lf, %lf\n", i, (double)xy[count], (double)integAcc[i]/(*max)); |
511 | 514 | ||
512 | if( ((integAcc[i+1]-integAcc[i])/(*max) > pdfmax) && (i < points-1) ) |
515 | if( ((integAcc[i+1]-integAcc[i])/(*max) > pdfmax) && (i < points-1) ) |
513 | pdfmax = (integAcc[i+1]-integAcc[i])/(*max); |
516 | pdfmax = (integAcc[i+1]-integAcc[i])/(*max); |
514 | } |
517 | } |
515 | 518 | ||
516 | pdfmax = (TMath::Ceil(pdfmax*10))/10.; |
519 | pdfmax = (TMath::Ceil(pdfmax*10))/10.; |
517 | 520 | ||
518 | // Prepare the PDF plot |
521 | // Prepare the PDF plot |
519 | gScan[0] = new TGraph(); |
522 | gScan[0] = new TGraph(); |
520 | for(int i = points-1; i >= 0; i--) |
523 | for(int i = points-1; i >= 0; i--) |
521 | { |
524 | { |
522 | count = (filenr-1) - (points-1) + i; |
525 | count = (filenr-1) - (points-1) + i; |
523 | // Set any negative values of the PDF to 0 |
526 | // Set any negative values of the PDF to 0 |
524 | if( (integAcc[i]-integAcc[i-1])/(*max) < 0 ) |
527 | if( (integAcc[i]-integAcc[i-1])/(*max) < 0 ) |
525 | gScan[0]->SetPoint(i, (double)xy[count], 0); |
528 | gScan[0]->SetPoint(i, (double)xy[count], 0); |
526 | else |
529 | else |
527 | gScan[0]->SetPoint(i, (double)xy[count], (integAcc[i]-integAcc[i-1])/(*max)); |
530 | gScan[0]->SetPoint(i, (double)xy[count], (integAcc[i]-integAcc[i-1])/(*max)); |
528 | if(DBGSIG) printf("PlotEdgeDistribution(): PDF %d: %lf, %lf\n", i, (double)xy[count], (integAcc[i+1]-integAcc[i])/(*max)); |
531 | if(DBGSIG) printf("PlotEdgeDistribution(): PDF %d: %lf, %lf\n", i, (double)xy[count], (integAcc[i+1]-integAcc[i])/(*max)); |
529 | } |
532 | } |
530 | 533 | ||
531 | remove_from_last((char*)files->At(filenr-1)->GetTitle(), '_', ctemp); |
534 | remove_from_last((char*)files->At(filenr-1)->GetTitle(), '_', ctemp); |
532 | sprintf(exportname, "%s_edge.pdf", ctemp); |
535 | sprintf(exportname, "%s_edge.pdf", ctemp); |
533 | 536 | ||
534 | if(edit == 0) |
537 | if(edit == 0) |
535 | gCanvas = new TCanvas("canv1d","canv1d",1200,900); |
538 | gCanvas = new TCanvas("canv1d","canv1d",1200,900); |
536 | else |
539 | else |
537 | gCanvas = tempAnalysisCanvas->GetCanvas(); |
540 | gCanvas = tempAnalysisCanvas->GetCanvas(); |
538 | 541 | ||
539 | // Fit the PDF with a gaussian |
542 | // Fit the PDF with a gaussian |
540 | gScan[0]->Fit("gaus","Q"); |
543 | gScan[0]->Fit("gaus","Q"); |
541 | gScan[0]->GetFunction("gaus")->SetNpx(400); |
544 | gScan[0]->GetFunction("gaus")->SetNpx(400); |
542 | 545 | ||
543 | gStyle->SetOptFit(1); |
546 | gStyle->SetOptFit(1); |
544 | 547 | ||
545 | gScan[1]->Draw("AL"); |
548 | gScan[1]->Draw("AL"); |
546 | gPad->Update(); |
549 | gPad->Update(); |
547 | gScan[0]->Draw("LP"); |
550 | gScan[0]->Draw("LP"); |
548 | 551 | ||
549 | gCanvas->Modified(); |
552 | gCanvas->Modified(); |
550 | gCanvas->Update(); |
553 | gCanvas->Update(); |
551 | 554 | ||
552 | TPaveStats *stats = (TPaveStats*)gScan[0]->FindObject("stats"); |
555 | TPaveStats *stats = (TPaveStats*)gScan[0]->FindObject("stats"); |
553 | if(!cleanPlots) |
556 | if(!cleanPlots) |
554 | { |
557 | { |
555 | stats->SetX1NDC(0.86); stats->SetX2NDC(1.0); |
558 | stats->SetX1NDC(0.86); stats->SetX2NDC(1.0); |
556 | stats->SetY1NDC(0.87); stats->SetY2NDC(1.0); |
559 | stats->SetY1NDC(0.87); stats->SetY2NDC(1.0); |
557 | } |
560 | } |
558 | else |
561 | else |
559 | { |
562 | { |
560 | stats->SetX1NDC(1.1); stats->SetX2NDC(1.3); |
563 | stats->SetX1NDC(1.1); stats->SetX2NDC(1.3); |
561 | stats->SetY1NDC(1.1); stats->SetY2NDC(1.3); |
564 | stats->SetY1NDC(1.1); stats->SetY2NDC(1.3); |
562 | } |
565 | } |
563 | 566 | ||
564 | gCanvas->SetGridx(1); |
567 | gCanvas->SetGridx(1); |
565 | gCanvas->SetGridy(1); |
568 | gCanvas->SetGridy(1); |
566 | if(axis == 1) |
569 | if(axis == 1) |
567 | gScan[1]->GetXaxis()->SetTitle("X [#mum]"); |
570 | gScan[1]->GetXaxis()->SetTitle("X [#mum]"); |
568 | else if(axis == 2) |
571 | else if(axis == 2) |
569 | gScan[1]->GetXaxis()->SetTitle("Y [#mum]"); |
572 | gScan[1]->GetXaxis()->SetTitle("Y [#mum]"); |
570 | gScan[1]->GetXaxis()->SetTitleOffset(1.3); |
573 | gScan[1]->GetXaxis()->SetTitleOffset(1.3); |
571 | gScan[1]->GetXaxis()->CenterTitle(kTRUE); |
574 | gScan[1]->GetXaxis()->CenterTitle(kTRUE); |
572 | gScan[1]->GetXaxis()->SetLabelSize(0.027); |
575 | gScan[1]->GetXaxis()->SetLabelSize(0.027); |
573 | gScan[1]->GetXaxis()->SetLabelOffset(0.02); |
576 | gScan[1]->GetXaxis()->SetLabelOffset(0.02); |
574 | gScan[1]->GetXaxis()->SetNoExponent(kTRUE); |
577 | gScan[1]->GetXaxis()->SetNoExponent(kTRUE); |
575 | gScan[1]->GetYaxis()->SetTitle("Normalized ADC integral"); |
578 | gScan[1]->GetYaxis()->SetTitle("Normalized ADC integral"); |
576 | 579 | ||
577 | gScan[1]->GetYaxis()->CenterTitle(kTRUE); |
580 | gScan[1]->GetYaxis()->CenterTitle(kTRUE); |
578 | gScan[1]->GetYaxis()->SetLabelSize(0.027); |
581 | gScan[1]->GetYaxis()->SetLabelSize(0.027); |
579 | gScan[1]->GetYaxis()->SetLabelOffset(0.02); |
582 | gScan[1]->GetYaxis()->SetLabelOffset(0.02); |
580 | gScan[1]->GetYaxis()->SetRangeUser(0,1); |
583 | gScan[1]->GetYaxis()->SetRangeUser(0,1); |
581 | gScan[1]->GetYaxis()->SetTitleOffset(1.4); |
584 | gScan[1]->GetYaxis()->SetTitleOffset(1.4); |
582 | gScan[1]->GetYaxis()->SetTitleSize(0.030); |
585 | gScan[1]->GetYaxis()->SetTitleSize(0.030); |
583 | 586 | ||
584 | if(!cleanPlots) |
587 | if(!cleanPlots) |
585 | { |
588 | { |
586 | if(axis == 1) |
589 | if(axis == 1) |
587 | { |
590 | { |
588 | if(posUnitsPlot->widgetCB->GetSelected() == 0) |
591 | if(posUnitsPlot->widgetCB->GetSelected() == 0) |
589 | gScan[1]->SetTitle("SiPM edge detection;X [table units];Normalized ADC integral"); |
592 | gScan[1]->SetTitle("SiPM edge detection;X [table units];Normalized ADC integral"); |
590 | else if(posUnitsPlot->widgetCB->GetSelected() == 1) |
593 | else if(posUnitsPlot->widgetCB->GetSelected() == 1) |
591 | gScan[1]->SetTitle("SiPM edge detection;X [#mum];Normalized ADC integral"); |
594 | gScan[1]->SetTitle("SiPM edge detection;X [#mum];Normalized ADC integral"); |
592 | } |
595 | } |
593 | else if(axis == 2) |
596 | else if(axis == 2) |
594 | { |
597 | { |
595 | if(posUnitsPlot->widgetCB->GetSelected() == 0) |
598 | if(posUnitsPlot->widgetCB->GetSelected() == 0) |
596 | gScan[1]->SetTitle("SiPM edge detection;Y [table units];Normalized ADC integral"); |
599 | gScan[1]->SetTitle("SiPM edge detection;Y [table units];Normalized ADC integral"); |
597 | else if(posUnitsPlot->widgetCB->GetSelected() == 1) |
600 | else if(posUnitsPlot->widgetCB->GetSelected() == 1) |
598 | gScan[1]->SetTitle("SiPM edge detection;Y [#mum];Normalized ADC integral"); |
601 | gScan[1]->SetTitle("SiPM edge detection;Y [#mum];Normalized ADC integral"); |
599 | } |
602 | } |
600 | } |
603 | } |
601 | else |
604 | else |
602 | { |
605 | { |
603 | if(axis == 1) |
606 | if(axis == 1) |
604 | { |
607 | { |
605 | if(posUnitsPlot->widgetCB->GetSelected() == 0) |
608 | if(posUnitsPlot->widgetCB->GetSelected() == 0) |
606 | gScan[1]->SetTitle(";X [table units];Normalized ADC integral"); |
609 | gScan[1]->SetTitle(";X [table units];Normalized ADC integral"); |
607 | else if(posUnitsPlot->widgetCB->GetSelected() == 1) |
610 | else if(posUnitsPlot->widgetCB->GetSelected() == 1) |
608 | gScan[1]->SetTitle(";X [#mum];Normalized ADC integral"); |
611 | gScan[1]->SetTitle(";X [#mum];Normalized ADC integral"); |
609 | } |
612 | } |
610 | else if(axis == 2) |
613 | else if(axis == 2) |
611 | { |
614 | { |
612 | if(posUnitsPlot->widgetCB->GetSelected() == 0) |
615 | if(posUnitsPlot->widgetCB->GetSelected() == 0) |
613 | gScan[1]->SetTitle(";Y [table units];Normalized ADC integral"); |
616 | gScan[1]->SetTitle(";Y [table units];Normalized ADC integral"); |
614 | else if(posUnitsPlot->widgetCB->GetSelected() == 1) |
617 | else if(posUnitsPlot->widgetCB->GetSelected() == 1) |
615 | gScan[1]->SetTitle(";Y [#mum];Normalized ADC integral"); |
618 | gScan[1]->SetTitle(";Y [#mum];Normalized ADC integral"); |
616 | } |
619 | } |
617 | } |
620 | } |
618 | gScan[1]->SetLineColor(kBlue); |
621 | gScan[1]->SetLineColor(kBlue); |
619 | gScan[0]->SetLineWidth(2); |
622 | gScan[0]->SetLineWidth(2); |
620 | gScan[1]->SetLineWidth(2); |
623 | gScan[1]->SetLineWidth(2); |
621 | 624 | ||
622 | gCanvas->Modified(); |
625 | gCanvas->Modified(); |
623 | gCanvas->Update(); |
626 | gCanvas->Update(); |
624 | 627 | ||
625 | gCanvas->SaveAs(exportname); |
628 | gCanvas->SaveAs(exportname); |
626 | 629 | ||
627 | // If 2D edge plot, delete the 1D edge plots as we go |
630 | // If 2D edge plot, delete the 1D edge plots as we go |
628 | if(edge2d) |
631 | if(edge2d) |
629 | { |
632 | { |
630 | delete gScan[0]; |
633 | delete gScan[0]; |
631 | delete gScan[1]; |
634 | delete gScan[1]; |
632 | if(edit == 0) |
635 | if(edit == 0) |
633 | delete gCanvas; |
636 | delete gCanvas; |
634 | } |
637 | } |
635 | else |
638 | else |
636 | { |
639 | { |
637 | if(edit == 0) |
640 | if(edit == 0) |
638 | { |
641 | { |
639 | delete gScan[0]; |
642 | delete gScan[0]; |
640 | delete gScan[1]; |
643 | delete gScan[1]; |
641 | delete gCanvas; |
644 | delete gCanvas; |
642 | } |
645 | } |
643 | } |
646 | } |
644 | } |
647 | } |
645 | 648 | ||
646 | void TGAppMainFrame::PhotonMu(TList *files, int edit) |
649 | void TGAppMainFrame::PhotonMu(TList *files, int edit) |
647 | { |
650 | { |
648 | unsigned int nrfiles = files->GetSize(); |
651 | unsigned int nrfiles = files->GetSize(); |
649 | char ctemp[1024]; |
652 | char ctemp[1024]; |
650 | char exportname[1024]; |
653 | char exportname[1024]; |
651 | int j, k = 0, m = 0, n = 0, k2 = 0, m2 = 0; |
654 | int j, k = 0, m = 0, n = 0, k2 = 0, m2 = 0; |
652 | 655 | ||
653 | TCanvas *gCanvas; |
656 | TCanvas *gCanvas; |
654 | TTree *header_data, *meas_data; |
657 | TTree *header_data, *meas_data; |
655 | double *integralCount, *integralPedestal; |
658 | double *integralCount, *integralPedestal; |
656 | integralCount = new double[nrfiles]; |
659 | integralCount = new double[nrfiles]; |
657 | integralPedestal = new double[nrfiles]; |
660 | integralPedestal = new double[nrfiles]; |
658 | double *angle, *pdeval, *muval; |
661 | double *angle, *pdeval, *muval; |
659 | angle = new double[nrfiles]; |
662 | angle = new double[nrfiles]; |
660 | pdeval = new double[nrfiles]; |
663 | pdeval = new double[nrfiles]; |
661 | muval = new double[nrfiles]; |
664 | muval = new double[nrfiles]; |
662 | 665 | ||
663 | // Zero the integral count and accumulated vaules |
666 | // Zero the integral count and accumulated vaules |
664 | for(int i = 0; i < (int)nrfiles; i++) {integralCount[i] = 0; integralPedestal[i] = 0; } |
667 | for(int i = 0; i < (int)nrfiles; i++) {integralCount[i] = 0; integralPedestal[i] = 0; } |
665 | 668 | ||
666 | // Fitting variables |
669 | // Fitting variables |
667 | TSpectrum *spec; |
670 | TSpectrum *spec; |
668 | TH1F *histtemp; |
671 | TH1F *histtemp; |
669 | TH1 *histback; |
672 | TH1 *histback; |
670 | TH1F *h2; |
673 | TH1F *h2; |
671 | float *xpeaks; |
674 | float *xpeaks; |
672 | TF1 *fit; |
675 | TF1 *fit; |
673 | TF1 *fittingfunc; |
676 | TF1 *fittingfunc; |
674 | double *fparam, *fparamerr; |
677 | double *fparam, *fparamerr; |
675 | double meansel[20]; |
678 | double meansel[20]; |
676 | double sigmasel[20]; |
679 | double sigmasel[20]; |
677 | double meanparam, paramsigma; |
680 | double meanparam, paramsigma; |
678 | int sortindex[20]; |
681 | int sortindex[20]; |
679 | int adcpedestal[2]; |
682 | int adcpedestal[2]; |
680 | int zeromu = 0; |
683 | int zeromu = 0; |
681 | int darkhist = -1; |
684 | int darkhist = -1; |
- | 685 | int nopeaks = -1; |
|
682 | 686 | ||
683 | double pointest[12]; |
687 | double pointest[12]; |
684 | bool exclude = false; |
688 | bool exclude = false; |
685 | 689 | ||
686 | // Zero the parameter values |
690 | // Zero the parameter values |
687 | for(int i = 0; i < 20; i++) {meansel[i] = 0; sigmasel[i] = 0; } |
691 | for(int i = 0; i < 20; i++) {meansel[i] = 0; sigmasel[i] = 0; } |
688 | 692 | ||
689 | float progVal = 0; |
693 | float progVal = 0; |
690 | analysisProgress->widgetPB->SetPosition(progVal); |
694 | analysisProgress->widgetPB->SetPosition(progVal); |
691 | gVirtualX->Update(1); |
695 | gVirtualX->Update(1); |
- | 696 | ||
- | 697 | // Check if the checkbox for no peaks is selected - TODO: Still need the situation when we do not have a peaked ADC spectrum |
|
- | 698 | if(relPde->widgetChBox[1]->IsDown()) |
|
- | 699 | { |
|
- | 700 | printf("PhotonMu(): ADC spectrum has no peak structure.\n"); |
|
- | 701 | nopeaks = 1; |
|
- | 702 | ||
- | 703 | // Error if there is no darkhist |
|
- | 704 | if(strcmp("", darkRun->widgetTE->GetText()) == 0) |
|
- | 705 | { |
|
- | 706 | printf("PhotonMu(): Error! The no peak structure option needs a dark histogram.\n"); |
|
- | 707 | delete[] integralCount; |
|
- | 708 | delete[] integralPedestal; |
|
- | 709 | delete[] angle; |
|
- | 710 | delete[] pdeval; |
|
- | 711 | delete[] muval; |
|
- | 712 | return; |
|
- | 713 | } |
|
- | 714 | ||
- | 715 | } |
|
692 | 716 | ||
693 | // Start if we select at least one file |
717 | // Start if we select at least one file |
694 | if(nrfiles > 0) |
718 | if(nrfiles > 0) |
695 | { |
719 | { |
- | 720 | // Find the pedestal peak for the dark histogram, and use it for all if there are no peaks |
|
- | 721 | if(nopeaks != -1) |
|
- | 722 | { |
|
- | 723 | // Replot the spectrum on analysisCanvas and do not close the input file |
|
- | 724 | DisplayHistogram( (char*)(darkRun->widgetTE->GetText()), 0, 1); |
|
- | 725 | analysisCanvas->GetCanvas()->Modified(); |
|
- | 726 | analysisCanvas->GetCanvas()->Update(); |
|
- | 727 | ||
- | 728 | // Get the spectrum |
|
- | 729 | histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname); |
|
- | 730 | npeaks = 15; |
|
- | 731 | double par[300]; |
|
- | 732 | spec = new TSpectrum(npeaks); |
|
- | 733 | // Find spectrum background |
|
- | 734 | histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same"); |
|
- | 735 | // Clone histogram and subtract background from it if we select that option |
|
- | 736 | h2 = (TH1F*)histtemp->Clone("h2"); |
|
- | 737 | if(fitChecks->widgetChBox[0]->IsDown()) |
|
- | 738 | h2->Add(histback, -1); |
|
- | 739 | // Search for the peaks |
|
- | 740 | int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() ); |
|
- | 741 | printf("PhotonMu(): Found %d candidates to fit.\n",found); |
|
- | 742 | npeaks = found; |
|
- | 743 | ||
- | 744 | // Set initial peak parameters |
|
- | 745 | xpeaks = spec->GetPositionX(); |
|
- | 746 | for(j = 0; j < found; j++) |
|
- | 747 | { |
|
- | 748 | float xp = xpeaks[j]; |
|
- | 749 | int bin = h2->GetXaxis()->FindBin(xp); |
|
- | 750 | float yp = h2->GetBinContent(bin); |
|
- | 751 | par[3*j] = yp; |
|
- | 752 | par[3*j+1] = xp; |
|
- | 753 | par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber(); |
|
- | 754 | } |
|
- | 755 | ||
- | 756 | // Fit the histogram |
|
- | 757 | fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks); |
|
- | 758 | TVirtualFitter::Fitter(histtemp, 3*npeaks); |
|
- | 759 | fit->SetParameters(par); |
|
- | 760 | fit->SetNpx(300); |
|
- | 761 | h2->Fit("fit","Q"); |
|
- | 762 | // Get the fitted parameters |
|
- | 763 | fittingfunc = h2->GetFunction("fit"); |
|
- | 764 | fparam = fittingfunc->GetParameters(); |
|
- | 765 | fparamerr = fittingfunc->GetParErrors(); |
|
- | 766 | ||
- | 767 | // Gather the parameters (mean peak value for now) |
|
- | 768 | int j = 1; |
|
- | 769 | int nrfit = 0; |
|
- | 770 | while(1) |
|
- | 771 | { |
|
- | 772 | if( (fparam[j] < 1.E-30) || (nrfit > 8) ) |
|
- | 773 | break; |
|
- | 774 | else |
|
- | 775 | { |
|
- | 776 | // Check if pedestal is above the lower limit and sigma is smaller than the mean |
|
- | 777 | if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) ) |
|
- | 778 | { |
|
- | 779 | // With the additional ADC offset, we can shift the mean values slightly, so they are not close to the X.5, but to the X.0 values |
|
- | 780 | meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber()); |
|
- | 781 | sigmasel[nrfit] = fparam[j+1]; |
|
- | 782 | nrfit++; |
|
- | 783 | } |
|
- | 784 | } |
|
- | 785 | ||
- | 786 | j+=3; |
|
- | 787 | } |
|
- | 788 | TMath::Sort(nrfit, meansel, sortindex, kFALSE); |
|
- | 789 | ||
- | 790 | meanparam = meansel[sortindex[0]]; |
|
- | 791 | paramsigma = sigmasel[sortindex[0]]; |
|
- | 792 | ||
- | 793 | for(j = 0; j < nrfit; j++) |
|
- | 794 | if(DBGSIG) |
|
- | 795 | printf("PhotonMu(): %d: peak mean = %lf\n", j, meansel[sortindex[j]]); |
|
- | 796 | ||
- | 797 | j = 0; |
|
- | 798 | adcpedestal[0] = 0; |
|
- | 799 | adcpedestal[1] = -1; |
|
- | 800 | ||
- | 801 | while(1) |
|
- | 802 | { |
|
- | 803 | int bin = histtemp->GetXaxis()->FindBin((int)(j+meanparam+paramsigma)); |
|
- | 804 | int yp = histtemp->GetBinContent(bin); |
|
- | 805 | ||
- | 806 | // Check where we get to first minimum after pedestal peak or where we get to the half maximum of the pedestal peak (in case there is only a pedestal peak) |
|
- | 807 | if(adcpedestal[1] == -1) |
|
- | 808 | { |
|
- | 809 | adcpedestal[0] = j+meanparam+paramsigma; |
|
- | 810 | adcpedestal[1] = yp; |
|
- | 811 | } |
|
- | 812 | else |
|
- | 813 | { |
|
- | 814 | if( (npeaks > 1) && (adcpedestal[1] >= yp) ) |
|
- | 815 | { |
|
- | 816 | adcpedestal[0] = j+meanparam+paramsigma; |
|
- | 817 | adcpedestal[1] = yp; |
|
- | 818 | } |
|
- | 819 | else if( (npeaks == 1) && (adcpedestal[0] < meanparam+5*paramsigma) ) // TODO -> Determining the pedestal when only one peak |
|
- | 820 | { |
|
- | 821 | adcpedestal[0] = j+meanparam+paramsigma; |
|
- | 822 | adcpedestal[1] = yp; |
|
- | 823 | } |
|
- | 824 | else |
|
- | 825 | break; |
|
- | 826 | } |
|
- | 827 | ||
- | 828 | j++; |
|
- | 829 | if(j > 50) break; |
|
- | 830 | } |
|
- | 831 | ||
- | 832 | if(midPeak->widgetChBox[0]->IsDown()) |
|
- | 833 | { |
|
- | 834 | if( (meanparam - (int)meanparam >= 0.) && (meanparam - (int)meanparam < 0.5) ) |
|
- | 835 | m = TMath::Floor(meanparam); |
|
- | 836 | else if( (meanparam - (int)meanparam >= 0.5) && (meanparam - (int)meanparam < 1.) ) |
|
- | 837 | m = TMath::Ceil(meanparam); |
|
- | 838 | int bin = histtemp->GetXaxis()->FindBin(m); |
|
- | 839 | adcpedestal[0] = m; |
|
- | 840 | printf("midpeak x = %d, ", adcpedestal[0]); |
|
- | 841 | adcpedestal[1] = histtemp->GetBinContent(bin); |
|
- | 842 | } |
|
- | 843 | ||
- | 844 | // Option to show the fit |
|
- | 845 | fittingfunc->Draw("L SAME"); |
|
- | 846 | analysisCanvas->GetCanvas()->Modified(); |
|
- | 847 | analysisCanvas->GetCanvas()->Update(); |
|
- | 848 | ||
- | 849 | printf("Pedestal ends = %d and nr. of counts = %d\n", adcpedestal[0], adcpedestal[1]); |
|
- | 850 | ||
- | 851 | // Delete the opened histogram and spectrum |
|
- | 852 | delete spec; |
|
- | 853 | delete inroot; |
|
- | 854 | ||
- | 855 | // return; |
|
- | 856 | } |
|
- | 857 | ||
- | 858 | printf("PhotonMu(): Continuing with the rest of the spectra.\n"); |
|
- | 859 | ||
- | 860 | // Check all histograms for pedestal peak values |
|
696 | for(int i = 0; i < (int)nrfiles; i++) |
861 | for(int i = 0; i < (int)nrfiles; i++) |
697 | { |
862 | { |
698 | if( (nrfiles == 1) || (!multiSelect->widgetChBox[0]->IsDown()) ) |
863 | if( (nrfiles == 1) || (!multiSelect->widgetChBox[0]->IsDown()) ) |
699 | { |
864 | { |
700 | printf("PhotonMu(): Only one file selected. Not running analysis, just showing the fit.\n"); |
865 | printf("PhotonMu(): Only one file selected. Not running analysis, just showing the fit.\n"); |
701 | 866 | ||
702 | // Replot the spectrum on analysisCanvas and do not close the input file |
867 | // Replot the spectrum on analysisCanvas and do not close the input file |
703 | DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1); |
868 | DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1); |
704 | analysisCanvas->GetCanvas()->Modified(); |
869 | analysisCanvas->GetCanvas()->Modified(); |
705 | analysisCanvas->GetCanvas()->Update(); |
870 | analysisCanvas->GetCanvas()->Update(); |
706 | 871 | ||
707 | // Get the spectrum |
872 | // Get the spectrum |
708 | histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname); |
873 | histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname); |
709 | npeaks = 15; |
874 | npeaks = 15; |
710 | double par[300]; |
875 | double par[300]; |
711 | spec = new TSpectrum(npeaks); |
876 | spec = new TSpectrum(npeaks); |
712 | // Find spectrum background |
877 | // Find spectrum background |
713 | histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same"); |
878 | histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same"); |
714 | // Clone histogram and subtract background from it if we select that option |
879 | // Clone histogram and subtract background from it if we select that option |
715 | h2 = (TH1F*)histtemp->Clone("h2"); |
880 | h2 = (TH1F*)histtemp->Clone("h2"); |
716 | if(fitChecks->widgetChBox[0]->IsDown()) |
881 | if(fitChecks->widgetChBox[0]->IsDown()) |
717 | h2->Add(histback, -1); |
882 | h2->Add(histback, -1); |
718 | // Search for the peaks |
883 | // Search for the peaks |
719 | int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() ); |
884 | int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() ); |
720 | printf("PhotonMu(): Found %d candidates to fit.\n",found); |
885 | printf("PhotonMu(): Found %d candidates to fit.\n",found); |
721 | npeaks = found; |
886 | npeaks = found; |
722 | 887 | ||
723 | // Set initial peak parameters |
888 | // Set initial peak parameters |
724 | xpeaks = spec->GetPositionX(); |
889 | xpeaks = spec->GetPositionX(); |
725 | for(j = 0; j < found; j++) |
890 | for(j = 0; j < found; j++) |
726 | { |
891 | { |
727 | float xp = xpeaks[j]; |
892 | float xp = xpeaks[j]; |
728 | int bin = h2->GetXaxis()->FindBin(xp); |
893 | int bin = h2->GetXaxis()->FindBin(xp); |
729 | float yp = h2->GetBinContent(bin); |
894 | float yp = h2->GetBinContent(bin); |
730 | par[3*j] = yp; |
895 | par[3*j] = yp; |
731 | par[3*j+1] = xp; |
896 | par[3*j+1] = xp; |
732 | par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber(); |
897 | par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber(); |
733 | } |
898 | } |
734 | 899 | ||
735 | // Fit the histogram |
900 | // Fit the histogram |
736 | fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks); |
901 | fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks); |
737 | TVirtualFitter::Fitter(histtemp, 3*npeaks); |
902 | TVirtualFitter::Fitter(histtemp, 3*npeaks); |
738 | fit->SetParameters(par); |
903 | fit->SetParameters(par); |
739 | fit->SetNpx(300); |
904 | fit->SetNpx(300); |
740 | h2->Fit("fit","Q"); |
905 | h2->Fit("fit","Q"); |
741 | // Get the fitted parameters |
906 | // Get the fitted parameters |
742 | fittingfunc = h2->GetFunction("fit"); |
907 | fittingfunc = h2->GetFunction("fit"); |
743 | fparam = fittingfunc->GetParameters(); |
908 | fparam = fittingfunc->GetParameters(); |
744 | fparamerr = fittingfunc->GetParErrors(); |
909 | fparamerr = fittingfunc->GetParErrors(); |
745 | 910 | ||
746 | // Gather the parameters (mean peak value for now) |
911 | // Gather the parameters (mean peak value for now) |
747 | int j = 1; |
912 | int j = 1; |
748 | int nrfit = 0; |
913 | int nrfit = 0; |
749 | while(1) |
914 | while(1) |
750 | { |
915 | { |
751 | if( (fparam[j] < 1.E-30) || (nrfit > 8) ) |
916 | if( (fparam[j] < 1.E-30) || (nrfit > 8) ) |
752 | break; |
917 | break; |
753 | else |
918 | else |
754 | { |
919 | { |
755 | // Check if pedestal is above the lower limit and sigma is smaller than the mean |
920 | // Check if pedestal is above the lower limit and sigma is smaller than the mean |
756 | if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) ) |
921 | if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) ) |
757 | { |
922 | { |
758 | // With the additional ADC offset, we can shift the mean values slightly, so they are not close to the X.5, but to the X.0 values |
923 | // With the additional ADC offset, we can shift the mean values slightly, so they are not close to the X.5, but to the X.0 values |
759 | meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber()); |
924 | meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber()); |
760 | sigmasel[nrfit] = fparam[j+1]; |
925 | sigmasel[nrfit] = fparam[j+1]; |
761 | nrfit++; |
926 | nrfit++; |
762 | } |
927 | } |
763 | } |
928 | } |
764 | 929 | ||
765 | j+=3; |
930 | j+=3; |
766 | } |
931 | } |
767 | TMath::Sort(nrfit, meansel, sortindex, kFALSE); |
932 | TMath::Sort(nrfit, meansel, sortindex, kFALSE); |
768 | 933 | ||
769 | fittingfunc->Draw("SAME"); |
934 | fittingfunc->Draw("SAME"); |
770 | analysisCanvas->GetCanvas()->Modified(); |
- | |
771 | analysisCanvas->GetCanvas()->Update(); |
- | |
772 | - | ||
773 | meanparam = meansel[sortindex[0]]; |
- | |
774 | paramsigma = sigmasel[sortindex[0]]; |
- | |
775 | - | ||
776 | for(j = 0; j < nrfit; j++) |
- | |
777 | printf("PhotonMu(): %d: peak mean = %lf\n", j, meansel[sortindex[j]]); |
- | |
778 | - | ||
779 | - | ||
780 | - | ||
781 | return; |
- | |
782 | } |
- | |
783 | if(files->At(i)) |
- | |
784 | { |
- | |
785 | if(strcmp(files->At(i)->GetTitle(),darkRun->widgetTE->GetText()) == 0) |
- | |
786 | { |
- | |
787 | printf("PhotonMu(): %s is the dark histogram file.\n", files->At(i)->GetTitle()); |
- | |
788 | darkhist = i; |
- | |
789 | } |
- | |
790 | - | ||
791 | // Replot the spectrum on analysisCanvas and do not close the input file |
- | |
792 | DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1); |
- | |
793 | analysisCanvas->GetCanvas()->Modified(); |
935 | analysisCanvas->GetCanvas()->Modified(); |
794 | analysisCanvas->GetCanvas()->Update(); |
936 | analysisCanvas->GetCanvas()->Update(); |
795 | - | ||
796 | // Get the spectrum |
- | |
797 | histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname); |
- | |
798 | npeaks = 15; |
- | |
799 | double par[300]; |
- | |
800 | spec = new TSpectrum(npeaks); |
- | |
801 | // Find spectrum background |
- | |
802 | histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same"); |
- | |
803 | // Clone histogram and subtract background from it if we select that option |
- | |
804 | h2 = (TH1F*)histtemp->Clone("h2"); |
- | |
805 | if(fitChecks->widgetChBox[0]->IsDown()) |
- | |
806 | h2->Add(histback, -1); |
- | |
807 | // Search for the peaks |
- | |
808 | int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() ); |
- | |
809 | printf("PhotonMu(): Found %d candidates to fit.\n",found); |
- | |
810 | npeaks = found; |
- | |
811 | - | ||
812 | // Set initial peak parameters |
- | |
813 | xpeaks = spec->GetPositionX(); |
- | |
814 | for(j = 0; j < found; j++) |
- | |
815 | { |
- | |
816 | float xp = xpeaks[j]; |
- | |
817 | int bin = h2->GetXaxis()->FindBin(xp); |
- | |
818 | float yp = h2->GetBinContent(bin); |
- | |
819 | par[3*j] = yp; |
- | |
820 | par[3*j+1] = xp; |
- | |
821 | par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber(); |
- | |
822 | } |
- | |
823 | - | ||
824 | // Fit the histogram |
- | |
825 | fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks); |
- | |
826 | TVirtualFitter::Fitter(histtemp, 3*npeaks); |
- | |
827 | fit->SetParameters(par); |
- | |
828 | fit->SetNpx(300); |
- | |
829 | h2->Fit("fit","Q"); |
- | |
830 | // Get the fitted parameters |
- | |
831 | fittingfunc = h2->GetFunction("fit"); |
- | |
832 | fparam = fittingfunc->GetParameters(); |
- | |
833 | fparamerr = fittingfunc->GetParErrors(); |
- | |
834 | - | ||
835 | // Gather the parameters (mean peak value for now) |
- | |
836 | int j = 1; |
- | |
837 | int nrfit = 0; |
- | |
838 | while(1) |
- | |
839 | { |
- | |
840 | if( (fparam[j] < 1.E-30) || (nrfit > 8) ) |
- | |
841 | break; |
- | |
842 | else |
- | |
843 | { |
- | |
844 | // Check if pedestal is above the lower limit and sigma is smaller than the mean |
- | |
845 | if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) ) |
- | |
846 | { |
- | |
847 | // With the additional ADC offset, we can shift the mean values slightly, so they are not close to the X.5, but to the X.0 values |
- | |
848 | meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber()); |
- | |
849 | sigmasel[nrfit] = fparam[j+1]; |
- | |
850 | nrfit++; |
- | |
851 | } |
- | |
852 | } |
- | |
853 | - | ||
854 | j+=3; |
- | |
855 | } |
- | |
856 | TMath::Sort(nrfit, meansel, sortindex, kFALSE); |
- | |
857 | 937 | ||
858 | meanparam = meansel[sortindex[0]]; |
938 | meanparam = meansel[sortindex[0]]; |
859 | paramsigma = sigmasel[sortindex[0]]; |
939 | paramsigma = sigmasel[sortindex[0]]; |
860 | 940 | ||
861 | for(j = 0; j < nrfit; j++) |
941 | for(j = 0; j < nrfit; j++) |
862 | if(DBGSIG) |
- | |
863 |
|
942 | printf("PhotonMu(): %d: peak mean = %lf\n", j, meansel[sortindex[j]]); |
- | 943 | ||
- | 944 | ||
- | 945 | ||
864 | 946 | return; |
|
- | 947 | } |
|
865 |
|
948 | if(files->At(i)) |
- | 949 | { |
|
866 |
|
950 | if(strcmp(files->At(i)->GetTitle(),darkRun->widgetTE->GetText()) == 0) |
- | 951 | { |
|
- | 952 | printf("PhotonMu(): %s is the dark histogram file.\n", files->At(i)->GetTitle()); |
|
867 |
|
953 | darkhist = i; |
- | 954 | } |
|
868 | 955 | ||
869 |
|
956 | if(nopeaks == -1) |
870 | { |
957 | { |
- | 958 | // Replot the spectrum on analysisCanvas and do not close the input file |
|
- | 959 | DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1); |
|
871 |
|
960 | analysisCanvas->GetCanvas()->Modified(); |
872 |
|
961 | analysisCanvas->GetCanvas()->Update(); |
873 | 962 | ||
- | 963 | // Get the spectrum |
|
874 |
|
964 | histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname); |
875 |
|
965 | npeaks = 15; |
876 | |
966 | double par[300]; |
877 |
|
967 | spec = new TSpectrum(npeaks); |
- | 968 | // Find spectrum background |
|
- | 969 | histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same"); |
|
- | 970 | // Clone histogram and subtract background from it if we select that option |
|
- | 971 | h2 = (TH1F*)histtemp->Clone("h2"); |
|
- | 972 | if(fitChecks->widgetChBox[0]->IsDown()) |
|
878 |
|
973 | h2->Add(histback, -1); |
879 |
|
974 | // Search for the peaks |
- | 975 | int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() ); |
|
- | 976 | printf("PhotonMu(): Found %d candidates to fit.\n",found); |
|
880 |
|
977 | npeaks = found; |
- | 978 | ||
- | 979 | // Set initial peak parameters |
|
- | 980 | xpeaks = spec->GetPositionX(); |
|
- | 981 | for(j = 0; j < found; j++) |
|
881 | { |
982 | { |
- | 983 | float xp = xpeaks[j]; |
|
- | 984 | int bin = h2->GetXaxis()->FindBin(xp); |
|
- | 985 | float yp = h2->GetBinContent(bin); |
|
- | 986 | par[3*j] = yp; |
|
- | 987 | par[3*j+1] = xp; |
|
- | 988 | par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber(); |
|
- | 989 | } |
|
- | 990 | ||
- | 991 | // Fit the histogram |
|
- | 992 | fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks); |
|
- | 993 | TVirtualFitter::Fitter(histtemp, 3*npeaks); |
|
- | 994 | fit->SetParameters(par); |
|
- | 995 | fit->SetNpx(300); |
|
- | 996 | h2->Fit("fit","Q"); |
|
- | 997 | // Get the fitted parameters |
|
- | 998 | fittingfunc = h2->GetFunction("fit"); |
|
- | 999 | fparam = fittingfunc->GetParameters(); |
|
- | 1000 | fparamerr = fittingfunc->GetParErrors(); |
|
- | 1001 | ||
- | 1002 | // Gather the parameters (mean peak value for now) |
|
- | 1003 | int j = 1; |
|
- | 1004 | int nrfit = 0; |
|
- | 1005 | while(1) |
|
- | 1006 | { |
|
882 | if( ( |
1007 | if( (fparam[j] < 1.E-30) || (nrfit > 8) ) |
- | 1008 | break; |
|
- | 1009 | else |
|
883 | { |
1010 | { |
- | 1011 | // Check if pedestal is above the lower limit and sigma is smaller than the mean |
|
- | 1012 | if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) ) |
|
- | 1013 | { |
|
- | 1014 | // With the additional ADC offset, we can shift the mean values slightly, so they are not close to the X.5, but to the X.0 values |
|
- | 1015 | meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber()); |
|
884 |
|
1016 | sigmasel[nrfit] = fparam[j+1]; |
885 |
|
1017 | nrfit++; |
- | 1018 | } |
|
886 | } |
1019 | } |
- | 1020 | ||
- | 1021 | j+=3; |
|
- | 1022 | } |
|
- | 1023 | TMath::Sort(nrfit, meansel, sortindex, kFALSE); |
|
- | 1024 | ||
- | 1025 | meanparam = meansel[sortindex[0]]; |
|
- | 1026 | paramsigma = sigmasel[sortindex[0]]; |
|
- | 1027 | ||
- | 1028 | for(j = 0; j < nrfit; j++) |
|
- | 1029 | if(DBGSIG) |
|
- | 1030 | printf("PhotonMu(): %d: peak mean = %lf\n", j, meansel[sortindex[j]]); |
|
- | 1031 | ||
- | 1032 | j = 0; |
|
- | 1033 | adcpedestal[0] = 0; |
|
- | 1034 | adcpedestal[1] = -1; |
|
- | 1035 | ||
- | 1036 | while(1) |
|
- | 1037 | { |
|
- | 1038 | int bin = histtemp->GetXaxis()->FindBin((int)(j+meanparam+paramsigma)); |
|
- | 1039 | int yp = histtemp->GetBinContent(bin); |
|
- | 1040 | ||
887 |
|
1041 | // Check where we get to first minimum after pedestal peak or where we get to the half maximum of the pedestal peak (in case there is only a pedestal peak) |
- | 1042 | if(adcpedestal[1] == -1) |
|
888 | { |
1043 | { |
889 | adcpedestal[0] = j+meanparam+paramsigma; |
1044 | adcpedestal[0] = j+meanparam+paramsigma; |
890 | adcpedestal[1] = yp; |
1045 | adcpedestal[1] = yp; |
891 | } |
1046 | } |
892 | else |
1047 | else |
- | 1048 | { |
|
- | 1049 | if( (npeaks > 1) && (adcpedestal[1] >= yp) ) |
|
- | 1050 | { |
|
- | 1051 | adcpedestal[0] = j+meanparam+paramsigma; |
|
- | 1052 | adcpedestal[1] = yp; |
|
- | 1053 | } |
|
- | 1054 | else if( (npeaks == 1) && (adcpedestal[0] < meanparam+5*paramsigma) ) // TODO -> Determining the pedestal when only one peak |
|
- | 1055 | { |
|
- | 1056 | adcpedestal[0] = j+meanparam+paramsigma; |
|
- | 1057 | adcpedestal[1] = yp; |
|
- | 1058 | } |
|
- | 1059 | else |
|
893 | break; |
1060 | break; |
- | 1061 | } |
|
- | 1062 | ||
- | 1063 | j++; |
|
- | 1064 | if(j > 50) break; |
|
894 | } |
1065 | } |
895 | - | ||
896 | j++; |
- | |
897 | if(j > 50) break; |
- | |
898 | } |
- | |
899 | - | ||
900 | if( (npeaks > 1) && (nrfit > 1) ) |
- | |
901 | { |
- | |
902 | int bin = histtemp->GetXaxis()->FindBin((int)(meanparam+meansel[sortindex[1]])/2); |
- | |
903 | adcpedestal[0] = (meanparam+meansel[sortindex[1]])/2; |
- | |
904 | printf("PhotonMu(): multipeak x = %d, ", adcpedestal[0]); |
- | |
905 | adcpedestal[1] = histtemp->GetBinContent(bin); |
- | |
906 | } |
- | |
907 | - | ||
908 | if(midPeak->widgetChBox[0]->IsDown()) |
- | |
909 | { |
- | |
910 | if( (meanparam - (int)meanparam >= 0.) && (meanparam - (int)meanparam < 0.5) ) |
- | |
911 | m = TMath::Floor(meanparam); |
- | |
912 | else if( (meanparam - (int)meanparam >= 0.5) && (meanparam - (int)meanparam < 1.) ) |
- | |
913 | m = TMath::Ceil(meanparam); |
- | |
914 | int bin = histtemp->GetXaxis()->FindBin(m); |
- | |
915 | adcpedestal[0] = m; |
- | |
916 | printf("midpeak x = %d, ", adcpedestal[0]); |
- | |
917 | adcpedestal[1] = histtemp->GetBinContent(bin); |
- | |
918 | } |
- | |
919 | 1066 | ||
- | 1067 | if(npeaks > 1) |
|
- | 1068 | { |
|
- | 1069 | int bin = histtemp->GetXaxis()->FindBin((int)(meanparam+meansel[sortindex[1]])/2); |
|
- | 1070 | adcpedestal[0] = (meanparam+meansel[sortindex[1]])/2; |
|
- | 1071 | printf("PhotonMu(): multipeak x = %d, ", adcpedestal[0]); |
|
- | 1072 | adcpedestal[1] = histtemp->GetBinContent(bin); |
|
- | 1073 | } |
|
- | 1074 | ||
- | 1075 | if(midPeak->widgetChBox[0]->IsDown()) |
|
- | 1076 | { |
|
- | 1077 | if( (meanparam - (int)meanparam >= 0.) && (meanparam - (int)meanparam < 0.5) ) |
|
- | 1078 | m = TMath::Floor(meanparam); |
|
- | 1079 | else if( (meanparam - (int)meanparam >= 0.5) && (meanparam - (int)meanparam < 1.) ) |
|
- | 1080 | m = TMath::Ceil(meanparam); |
|
- | 1081 | int bin = histtemp->GetXaxis()->FindBin(m); |
|
- | 1082 | adcpedestal[0] = m; |
|
- | 1083 | printf("midpeak x = %d, ", adcpedestal[0]); |
|
- | 1084 | adcpedestal[1] = histtemp->GetBinContent(bin); |
|
- | 1085 | } |
|
- | 1086 | ||
920 | /* // Option to show the fit |
1087 | /* // Option to show the fit |
921 | fittingfunc->Draw("L SAME"); |
1088 | fittingfunc->Draw("L SAME"); |
922 | analysisCanvas->GetCanvas()->Modified(); |
1089 | analysisCanvas->GetCanvas()->Modified(); |
923 | analysisCanvas->GetCanvas()->Update();*/ |
1090 | analysisCanvas->GetCanvas()->Update();*/ |
924 | 1091 | ||
925 | printf("Pedestal ends = %d and nr. of counts = %d\n", adcpedestal[0], adcpedestal[1]); |
1092 | printf("Pedestal ends = %d and nr. of counts = %d\n", adcpedestal[0], adcpedestal[1]); |
926 | 1093 | ||
927 | // Delete the opened histogram and spectrum |
1094 | // Delete the opened histogram and spectrum |
928 | delete spec; |
1095 | delete spec; |
929 | delete inroot; |
1096 | delete inroot; |
- | 1097 | } |
|
930 | 1098 | ||
931 | // Open the input file and read header, ADC and TDC values |
1099 | // Open the input file and read header, ADC and TDC values |
932 | sprintf(ctemp, "%s", files->At(i)->GetTitle()); |
1100 | sprintf(ctemp, "%s", files->At(i)->GetTitle()); |
933 | inroot = new TFile(ctemp, "READ"); |
1101 | inroot = new TFile(ctemp, "READ"); |
934 | 1102 | ||
935 | inroot->GetObject("header_data", header_data); |
1103 | inroot->GetObject("header_data", header_data); |
936 | inroot->GetObject("meas_data", meas_data); |
1104 | inroot->GetObject("meas_data", meas_data); |
937 | 1105 | ||
938 | // Reading the header |
1106 | // Reading the header |
939 | if( header_data->FindBranch("angle") ) |
1107 | if( header_data->FindBranch("angle") ) |
940 | { |
1108 | { |
941 | header_data->SetBranchAddress("angle", &evtheader.angle); |
1109 | header_data->SetBranchAddress("angle", &evtheader.angle); |
942 | header_data->GetEntry(0); |
1110 | header_data->GetEntry(0); |
943 | } |
1111 | } |
944 | else |
1112 | else |
945 | { |
1113 | { |
946 | printf("PhotonMu(): Error! Selected file has no angle header value. Please edit header to add the angle header value.\n"); |
1114 | printf("PhotonMu(): Error! Selected file has no angle header value. Please edit header to add the angle header value.\n"); |
947 | break; |
1115 | break; |
948 | } |
1116 | } |
949 | 1117 | ||
950 | char rdc[256]; |
1118 | char rdc[256]; |
951 | j = selectCh->widgetNE[0]->GetNumber(); |
1119 | j = selectCh->widgetNE[0]->GetNumber(); |
952 | double rangetdc[2]; |
1120 | double rangetdc[2]; |
953 | rangetdc[0] = tdcRange->widgetNE[0]->GetNumber(); |
1121 | rangetdc[0] = tdcRange->widgetNE[0]->GetNumber(); |
954 | rangetdc[1] = tdcRange->widgetNE[1]->GetNumber(); |
1122 | rangetdc[1] = tdcRange->widgetNE[1]->GetNumber(); |
955 | 1123 | ||
956 | k = 0; |
1124 | k = 0; |
957 | k2 = 0; |
1125 | k2 = 0; |
958 | m = 0; |
1126 | m = 0; |
959 | m2 = 0; |
1127 | m2 = 0; |
960 | 1128 | ||
961 | // Reading the data |
1129 | // Reading the data |
962 | for(int e = 0; e < meas_data->GetEntries(); e++) |
1130 | for(int e = 0; e < meas_data->GetEntries(); e++) |
963 | { |
1131 | { |
964 | sprintf(rdc, "ADC%d", j); |
1132 | sprintf(rdc, "ADC%d", j); |
965 | meas_data->SetBranchAddress(rdc, &evtdata.adcdata[j]); |
1133 | meas_data->SetBranchAddress(rdc, &evtdata.adcdata[j]); |
966 | meas_data->GetEntry(e); |
1134 | meas_data->GetEntry(e); |
967 | 1135 | ||
968 | sprintf(rdc, "TDC%d", j); |
1136 | sprintf(rdc, "TDC%d", j); |
969 | meas_data->SetBranchAddress(rdc, &evtdata.tdcdata[j]); |
1137 | meas_data->SetBranchAddress(rdc, &evtdata.tdcdata[j]); |
970 | meas_data->GetEntry(e); |
1138 | meas_data->GetEntry(e); |
971 | 1139 | ||
972 | // If our data point is inside the TDC window |
1140 | // If our data point is inside the TDC window |
973 | if( ((double)evtdata.tdcdata[j]/tdctimeconversion >= rangetdc[0]) && ((double)evtdata.tdcdata[j]/tdctimeconversion <= rangetdc[1]) ) |
1141 | if( ((double)evtdata.tdcdata[j]/tdctimeconversion >= rangetdc[0]) && ((double)evtdata.tdcdata[j]/tdctimeconversion <= rangetdc[1]) ) |
974 | { |
1142 | { |
975 | // Gather only the integral of the pedestal |
1143 | // Gather only the integral of the pedestal |
976 | if((double)evtdata.adcdata[j] < (double)adcpedestal[0]+0.5 ) |
1144 | if((double)evtdata.adcdata[j] < (double)adcpedestal[0]+0.5 ) |
977 | { |
1145 | { |
978 | k2++; |
1146 | k2++; |
979 | m2 += evtdata.adcdata[j]; |
1147 | m2 += evtdata.adcdata[j]; |
980 | } |
1148 | } |
981 | 1149 | ||
982 | // Gather the complete integral |
1150 | // Gather the complete integral |
983 | k++; |
1151 | k++; |
984 | m += evtdata.adcdata[j]; |
1152 | m += evtdata.adcdata[j]; |
985 | } |
1153 | } |
986 | } |
1154 | } |
987 | 1155 | ||
988 | // Determine the angle, mu and relative PDE |
1156 | // Determine the angle, mu and relative PDE |
989 | angle[i] = (double)(evtheader.angle); |
1157 | angle[i] = (double)(evtheader.angle); |
990 | integralCount[i] += (double)m; |
1158 | integralCount[i] += (double)m; |
991 | printf("PhotonMu(): %lf: Complete integral (%d evts) = %lf\n", angle[i], k, integralCount[i]); |
1159 | printf("PhotonMu(): %lf: Complete integral (%d evts) = %lf\n", angle[i], k, integralCount[i]); |
992 | integralPedestal[i] += (double)m2; |
1160 | integralPedestal[i] += (double)m2; |
993 | printf("PhotonMu(): %lf: Pedestal integral (%d evts) = %lf\n", angle[i], k2, integralPedestal[i]); |
1161 | printf("PhotonMu(): %lf: Pedestal integral (%d evts) = %lf\n", angle[i], k2, integralPedestal[i]); |
994 | if( (angle[i] == zeroAngle->widgetNE[0]->GetNumber()) && (darkhist != i) ) |
1162 | if( (angle[i] == zeroAngle->widgetNE[0]->GetNumber()) && (darkhist != i) ) |
995 | zeromu = i; |
1163 | zeromu = i; |
996 | 1164 | ||
997 | // Checking for errors when fitting a histogram |
- | |
998 | if(k2 == 0) |
- | |
999 | { |
- | |
1000 | printf("PhotonMu(): No pedestal entries found. Check the fitting results.\n"); |
- | |
1001 | muval[i] = -1; |
- | |
1002 | } |
- | |
1003 | else |
- | |
1004 |
|
1165 | muval[i] = -TMath::Log((double)k2/(double)k); |
1005 | printf("PhotonMu(): %lf: muval = %lf\n", angle[i], muval[i]); |
1166 | printf("PhotonMu(): %lf: muval = %lf\n", angle[i], muval[i]); |
1006 | 1167 | ||
1007 | inroot->Close(); |
1168 | inroot->Close(); |
1008 | delete inroot; |
1169 | delete inroot; |
1009 | } |
1170 | } |
1010 | 1171 | ||
1011 | // Update the progress bar |
1172 | // Update the progress bar |
1012 | progVal = (float)(90.00/nrfiles)*i; |
1173 | progVal = (float)(90.00/nrfiles)*i; |
1013 | analysisProgress->widgetPB->SetPosition(progVal); |
1174 | analysisProgress->widgetPB->SetPosition(progVal); |
1014 | gVirtualX->Update(1); |
1175 | gVirtualX->Update(1); |
1015 | } |
1176 | } |
1016 | 1177 | ||
1017 | printf("PhotonMu(): %d files were selected.\n", nrfiles); |
1178 | printf("PhotonMu(): %d files were selected.\n", nrfiles); |
1018 | 1179 | ||
1019 | printf("PhotonMu(): angle\tmu\trelative PDE\n"); |
1180 | printf("PhotonMu(): angle\tmu\trelative PDE\n"); |
1020 | m = 0; |
1181 | m = 0; |
1021 | 1182 | ||
1022 | // Set the 0 degree muval, reuse meansel[1] |
1183 | // Set the 0 degree muval, reuse meansel[1] |
1023 | meansel[1] = muval[zeromu]; |
1184 | meansel[1] = muval[zeromu]; |
1024 | printf("Zero value (id=%d, angle=%lf) = %lf\n", zeromu, angle[zeromu], meansel[1]); |
1185 | printf("Zero value (id=%d, angle=%lf) = %lf\n", zeromu, angle[zeromu], meansel[1]); |
1025 | 1186 | ||
1026 | // TODO - point estimation still not working correctly! |
1187 | // TODO - point estimation still not working correctly! |
1027 | for(int i = 0; i < (int)nrfiles; i++) |
1188 | for(int i = 0; i < (int)nrfiles; i++) |
1028 | { |
1189 | { |
1029 | // Estimate next point and check error (5 point least square fit estimation) |
1190 | // Estimate next point and check error (5 point least square fit estimation) |
1030 | if( ((i > 4) && (m == i)) || ((i > 5) && (m < i)) ) |
1191 | if( ((i > 4) && (m == i)) || ((i > 5) && (m < i)) ) |
1031 | { |
1192 | { |
1032 | // Set exclude signal to false |
1193 | // Set exclude signal to false |
1033 | exclude = false; |
1194 | exclude = false; |
1034 | 1195 | ||
1035 | // Get next point values (if zero value -> need to add the dark hist value again) |
1196 | // Get next point values (if zero value -> need to add the dark hist value again) |
1036 | pointest[10] = angle[i]; |
1197 | pointest[10] = angle[i]; |
1037 | pointest[11] = muval[i]; |
1198 | pointest[11] = muval[i]; |
1038 | 1199 | ||
1039 | // Check if next point has larger error than acceptable (if yes, set exclude signal to true), reuse meansel[0] |
1200 | // Check if next point has larger error than acceptable (if yes, set exclude signal to true), reuse meansel[0] |
1040 | meansel[0] = PointEstimate(5, pointest); // PointEstimate only works with very small step size |
1201 | meansel[0] = PointEstimate(5, pointest); // PointEstimate only works with very small step size |
1041 | if(meansel[0] > accError->widgetNE[0]->GetNumber()) |
1202 | if(meansel[0] > accError->widgetNE[0]->GetNumber()) |
1042 | { |
1203 | { |
1043 | printf("PhotonMu(): Point (%lf, %lf) excluded due to error: %lf\n", pointest[10], pointest[11], meansel[0]); |
1204 | printf("PhotonMu(): Point (%lf, %lf) excluded due to error: %lf\n", pointest[10], pointest[11], meansel[0]); |
1044 | exclude = true; |
1205 | exclude = true; |
1045 | } |
1206 | } |
1046 | 1207 | ||
1047 | // Value with 0 angle and dark histogram are always needed, so should not be excluded |
1208 | // Value with 0 angle and dark histogram are always needed, so should not be excluded |
1048 | if(i == darkhist) |
1209 | if(i == darkhist) |
1049 | exclude = false; |
1210 | exclude = false; |
1050 | - | ||
1051 | // Wrong fit |
- | |
1052 | if(muval[i] == -1) |
- | |
1053 | exclude = true; |
- | |
1054 | 1211 | ||
1055 | // If nothing excluded, pass the points in pointest variable like in a FIFO |
1212 | // If nothing excluded, pass the points in pointest variable like in a FIFO |
1056 | if(!exclude) |
1213 | if(!exclude) |
1057 | { |
1214 | { |
1058 | for(int j = 0; j < 10; j++) |
1215 | for(int j = 0; j < 10; j++) |
1059 | { |
1216 | { |
1060 | if(DBGSIG) printf("PhotonMu(): j = %d: Old X value = %lf\n", j, pointest[j]); |
1217 | if(DBGSIG) printf("PhotonMu(): j = %d: Old X value = %lf\n", j, pointest[j]); |
1061 | pointest[j] = pointest[j+2]; |
1218 | pointest[j] = pointest[j+2]; |
1062 | } |
1219 | } |
1063 | } |
1220 | } |
1064 | else |
1221 | else |
1065 | { |
1222 | { |
1066 | for(int j = 0; j < 10; j++) |
1223 | for(int j = 0; j < 10; j++) |
1067 | if(DBGSIG) printf("PhotonMu(): No j = %d: Old X value = %lf\n", j, pointest[j]); |
1224 | if(DBGSIG) printf("PhotonMu(): No j = %d: Old X value = %lf\n", j, pointest[j]); |
1068 | } |
1225 | } |
1069 | } |
1226 | } |
1070 | else |
1227 | else |
1071 | { |
1228 | { |
1072 | // First 5 points act as estimator points for next one |
1229 | // First 5 points act as estimator points for next one |
1073 | pointest[2*m] = angle[i]; |
1230 | pointest[2*m] = angle[i]; |
1074 | pointest[2*m+1] = muval[i]; |
1231 | pointest[2*m+1] = muval[i]; |
1075 | } |
1232 | } |
1076 | 1233 | ||
1077 | // Run only if we have a dark run histogram and middle pedestal peak estimation |
1234 | // Run only if we have a dark run histogram and middle pedestal peak estimation |
1078 | if( (darkhist != -1) && midPeak->widgetChBox[0]->IsDown() ) |
1235 | if( (darkhist != -1) && midPeak->widgetChBox[0]->IsDown() ) |
1079 | { |
1236 | { |
1080 | if(DBGSIG) printf("PhotonMu(): m = %d, i = %d: muval = %lf, ", m, i, muval[i]); |
1237 | if(DBGSIG) printf("PhotonMu(): m = %d, i = %d: muval = %lf, ", m, i, muval[i]); |
1081 | 1238 | ||
1082 | // Subtract the dark value from all values |
1239 | // Subtract the dark value from all values |
1083 | angle[m] = angle[i]; |
1240 | angle[m] = angle[i]; |
1084 | muval[m] = muval[i] - muval[darkhist]; |
1241 | muval[m] = muval[i] - muval[darkhist]; |
1085 | 1242 | ||
1086 | if(DBGSIG) printf("angle = %lf, newmuval = %lf, darkmuval = %lf, ", angle[m], muval[m], muval[darkhist]); |
1243 | if(DBGSIG) printf("angle = %lf, newmuval = %lf, darkmuval = %lf, ", angle[m], muval[m], muval[darkhist]); |
1087 | 1244 | ||
1088 | // Calculate relative PDE |
1245 | // Calculate relative PDE |
1089 | // pdeval[m] = muval[m]/(muval[zeromu]*TMath::Cos(angle[m]*TMath::ACos(-1.)/180.)); |
1246 | // pdeval[m] = muval[m]/(muval[zeromu]*TMath::Cos(angle[m]*TMath::ACos(-1.)/180.)); |
1090 | pdeval[m] = muval[m]/((meansel[1]-muval[darkhist])*TMath::Cos(angle[m]*TMath::ACos(-1.)/180.)); |
1247 | pdeval[m] = muval[m]/((meansel[1]-muval[darkhist])*TMath::Cos(angle[m]*TMath::ACos(-1.)/180.)); |
1091 | if(DBGSIG) printf("pdeval = %lf\n", pdeval[m]); |
1248 | if(DBGSIG) printf("pdeval = %lf\n", pdeval[m]); |
1092 | 1249 | ||
1093 | // Only increase counter if error is not too big |
1250 | // Only increase counter if error is not too big |
1094 | if( (darkhist != i) && (!exclude) ) |
1251 | if( (darkhist != i) && (!exclude) ) |
1095 | m++; |
1252 | m++; |
1096 | } |
1253 | } |
1097 | else |
1254 | else |
1098 | { |
1255 | { |
1099 | // Relative PDE calculation |
1256 | // Relative PDE calculation |
1100 | angle[m] = angle[i]; |
1257 | angle[m] = angle[i]; |
1101 | muval[m] = muval[i]; |
1258 | muval[m] = muval[i]; |
1102 | // pdeval[m] = muval[m]/(muval[zeromu]*TMath::Cos(angle[m]*TMath::ACos(-1.)/180.)); |
1259 | // pdeval[m] = muval[m]/(muval[zeromu]*TMath::Cos(angle[m]*TMath::ACos(-1.)/180.)); |
1103 | pdeval[m] = muval[m]/(meansel[1]*TMath::Cos(angle[m]*TMath::ACos(-1.)/180.)); |
1260 | pdeval[m] = muval[m]/(meansel[1]*TMath::Cos(angle[m]*TMath::ACos(-1.)/180.)); |
1104 | 1261 | ||
1105 | // Only increase counter if error is not too big |
1262 | // Only increase counter if error is not too big |
1106 | if(!exclude) |
1263 | if(!exclude) |
1107 | m++; |
1264 | m++; |
1108 | } |
1265 | } |
1109 | printf("PhotonMu(): %lf\t%lf\t%lf\n", angle[i], muval[i], pdeval[i]); |
1266 | printf("PhotonMu(): %lf\t%lf\t%lf\n", angle[i], muval[i], pdeval[i]); |
1110 | } |
1267 | } |
1111 | 1268 | ||
1112 | // Check for range of values to plot |
1269 | // Check for range of values to plot |
1113 | double plotMax = 0.; |
1270 | double plotMax = 0.; |
1114 | for(int i = 0; i < (int)nrfiles; i++) |
1271 | for(int i = 0; i < (int)nrfiles; i++) |
1115 | { |
1272 | { |
1116 | plotMax = TMath::Max(plotMax, muval[i]); |
1273 | plotMax = TMath::Max(plotMax, muval[i]); |
1117 | plotMax = TMath::Max(plotMax, pdeval[i]); |
1274 | plotMax = TMath::Max(plotMax, pdeval[i]); |
1118 | } |
1275 | } |
1119 | if(plotMax <= 0.) |
1276 | if(plotMax <= 0.) |
1120 | plotMax = 1.1; |
1277 | plotMax = 1.1; |
1121 | printf("PhotonMu(): Maximum value: %lf\n", plotMax); |
1278 | printf("PhotonMu(): Maximum value: %lf\n", plotMax); |
1122 | 1279 | ||
1123 | if(DBGSIG) printf("\n"); |
1280 | if(DBGSIG) printf("\n"); |
1124 | if(darkhist != -1) |
1281 | if(darkhist != -1) |
1125 | printf("PhotonMu(): Number of excluded points: %d\n", (nrfiles-1-m)); |
1282 | printf("PhotonMu(): Number of excluded points: %d\n", (nrfiles-1-m)); |
1126 | else |
1283 | else |
1127 | printf("PhotonMu(): Number of excluded points: %d\n", (nrfiles-m)); |
1284 | printf("PhotonMu(): Number of excluded points: %d\n", (nrfiles-m)); |
1128 | 1285 | ||
1129 | // Plot mu and PDE angle dependance plots |
1286 | // Plot mu and PDE angle dependance plots |
1130 | if(edit == 0) |
1287 | if(edit == 0) |
1131 | gCanvas = new TCanvas("canv","canv",1200,900); |
1288 | gCanvas = new TCanvas("canv","canv",1200,900); |
1132 | else if(edit == 1) |
1289 | else if(edit == 1) |
1133 | gCanvas = tempAnalysisCanvas->GetCanvas(); |
1290 | gCanvas = tempAnalysisCanvas->GetCanvas(); |
1134 | gCanvas->cd(); |
1291 | gCanvas->cd(); |
1135 | gCanvas->SetGrid(); |
1292 | gCanvas->SetGrid(); |
1136 | 1293 | ||
1137 | TGraph *pde = new TGraph(m, angle, pdeval); |
1294 | TGraph *pde = new TGraph(m, angle, pdeval); |
1138 | pde->SetMarkerStyle(21); |
1295 | pde->SetMarkerStyle(21); |
1139 | pde->SetMarkerSize(0.7); |
1296 | pde->SetMarkerSize(0.7); |
1140 | pde->SetMarkerColor(2); |
1297 | pde->SetMarkerColor(2); |
1141 | pde->SetLineWidth(1); |
1298 | pde->SetLineWidth(1); |
1142 | pde->SetLineColor(2); |
1299 | pde->SetLineColor(2); |
1143 | pde->GetXaxis()->SetLabelSize(0.030); |
1300 | pde->GetXaxis()->SetLabelSize(0.030); |
1144 | pde->GetXaxis()->CenterTitle(); |
1301 | pde->GetXaxis()->CenterTitle(); |
1145 | // pde->GetXaxis()->SetRange(angle[0],angle[nrfiles-1]); |
1302 | // pde->GetXaxis()->SetRange(angle[0],angle[nrfiles-1]); |
1146 | // pde->GetXaxis()->SetRangeUser(angle[0],angle[nrfiles-1]); |
1303 | // pde->GetXaxis()->SetRangeUser(angle[0],angle[nrfiles-1]); |
1147 | pde->GetXaxis()->SetRange(-90.0,90.0); |
1304 | pde->GetXaxis()->SetRange(-90.0,90.0); |
1148 | pde->GetXaxis()->SetRangeUser(-90.0,90.0); |
1305 | pde->GetXaxis()->SetRangeUser(-90.0,90.0); |
1149 | pde->GetXaxis()->SetLimits(-90.0,90.0); |
1306 | pde->GetXaxis()->SetLimits(-90.0,90.0); |
1150 | pde->GetYaxis()->SetTitleOffset(1.2); |
1307 | pde->GetYaxis()->SetTitleOffset(1.2); |
1151 | pde->GetYaxis()->SetLabelSize(0.030); |
1308 | pde->GetYaxis()->SetLabelSize(0.030); |
1152 | pde->GetYaxis()->CenterTitle(); |
1309 | pde->GetYaxis()->CenterTitle(); |
1153 | pde->GetYaxis()->SetRange(0., 1.1*plotMax); |
1310 | pde->GetYaxis()->SetRange(0., 1.1*plotMax); |
1154 | pde->GetYaxis()->SetRangeUser(0., 1.1*plotMax); |
1311 | pde->GetYaxis()->SetRangeUser(0., 1.1*plotMax); |
1155 | pde->GetYaxis()->SetLimits(0., 1.1*plotMax); |
1312 | pde->GetYaxis()->SetLimits(0., 1.1*plotMax); |
1156 | pde->SetName("pde"); |
1313 | pde->SetName("pde"); |
1157 | pde->Draw("ALP"); |
1314 | pde->Draw("ALP"); |
1158 | 1315 | ||
1159 | pde->SetTitle(";Incidence angle (#circ);Relative PDE(#theta) / #mu(#theta)"); |
1316 | pde->SetTitle(";Incidence angle (#circ);Relative PDE(#theta) / #mu(#theta)"); |
1160 | 1317 | ||
1161 | TGraph *mugr = new TGraph(m, angle, muval); |
1318 | TGraph *mugr = new TGraph(m, angle, muval); |
1162 | mugr->SetMarkerStyle(20); |
1319 | mugr->SetMarkerStyle(20); |
1163 | mugr->SetMarkerSize(0.7); |
1320 | mugr->SetMarkerSize(0.7); |
1164 | mugr->SetMarkerColor(4); |
1321 | mugr->SetMarkerColor(4); |
1165 | mugr->SetLineWidth(1); |
1322 | mugr->SetLineWidth(1); |
1166 | mugr->SetLineColor(4); |
1323 | mugr->SetLineColor(4); |
1167 | mugr->SetName("muval"); |
1324 | mugr->SetName("muval"); |
1168 | mugr->Draw("SAME;LP"); |
1325 | mugr->Draw("SAME;LP"); |
1169 | 1326 | ||
1170 | gCanvas->Modified(); |
1327 | gCanvas->Modified(); |
1171 | gCanvas->Update(); |
1328 | gCanvas->Update(); |
1172 | 1329 | ||
1173 | if(edit == 0) |
1330 | if(edit == 0) |
1174 | { |
1331 | { |
1175 | remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp); |
1332 | remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp); |
1176 | sprintf(exportname, "%s_pde.pdf", ctemp); |
1333 | sprintf(exportname, "%s_pde.pdf", ctemp); |
1177 | gCanvas->SaveAs(exportname); |
1334 | gCanvas->SaveAs(exportname); |
1178 | 1335 | ||
1179 | delete mugr; |
1336 | delete mugr; |
1180 | delete pde; |
1337 | delete pde; |
1181 | delete gCanvas; |
1338 | delete gCanvas; |
1182 | } |
1339 | } |
1183 | else if(edit == 1) |
1340 | else if(edit == 1) |
1184 | { |
1341 | { |
1185 | gCanvas->Modified(); |
1342 | gCanvas->Modified(); |
1186 | gCanvas->Update(); |
1343 | gCanvas->Update(); |
1187 | } |
1344 | } |
1188 | 1345 | ||
1189 | // Update the progress bar |
1346 | // Update the progress bar |
1190 | analysisProgress->widgetPB->SetPosition(100.); |
1347 | analysisProgress->widgetPB->SetPosition(100.); |
1191 | gVirtualX->Update(1); |
1348 | gVirtualX->Update(1); |
1192 | } |
1349 | } |
1193 | } |
1350 | } |
1194 | 1351 | ||
1195 | void TGAppMainFrame::BreakdownVolt(TList *files, int edit) |
1352 | void TGAppMainFrame::BreakdownVolt(TList *files, int edit) |
1196 | { |
1353 | { |
1197 | unsigned int nrfiles = files->GetSize(); |
1354 | unsigned int nrfiles = files->GetSize(); |
1198 | char ctemp[1024]; |
1355 | char ctemp[1024]; |
1199 | char exportname[1024]; |
1356 | char exportname[1024]; |
1200 | char paramname[1024]; |
1357 | char paramname[1024]; |
1201 | int j, k = 0; |
1358 | int j, k = 0; |
1202 | 1359 | ||
1203 | TCanvas *gCanvas; |
1360 | TCanvas *gCanvas; |
1204 | TTree *header_data, *meas_data; |
1361 | TTree *header_data, *meas_data; |
1205 | 1362 | ||
1206 | // Fitting variables |
1363 | // Fitting variables |
1207 | TSpectrum *spec; |
1364 | TSpectrum *spec; |
1208 | TH1F *histtemp; |
1365 | TH1F *histtemp; |
1209 | TH1 *histback; |
1366 | TH1 *histback; |
1210 | TH1F *h2; |
1367 | TH1F *h2; |
1211 | float *xpeaks; |
1368 | float *xpeaks; |
1212 | TF1 *fit; |
1369 | TF1 *fit; |
1213 | TF1 *fittingfunc; |
1370 | TF1 *fittingfunc; |
1214 | TLatex *latex; |
1371 | TLatex *latex; |
1215 | double *fparam, *fparamerr; |
1372 | double *fparam, *fparamerr; |
1216 | double meansel[20]; |
1373 | double meansel[20]; |
1217 | double meanselerr[20]; |
1374 | double meanselerr[20]; |
1218 | double sigmasel[20]; |
1375 | double sigmasel[20]; |
1219 | double meanparam, meanparamerr, paramsigma; |
1376 | double meanparam, meanparamerr, paramsigma; |
1220 | int sortindex[20]; |
1377 | int sortindex[20]; |
1221 | bool exclude = false; |
1378 | bool exclude = false; |
1222 | 1379 | ||
1223 | int p = 0; |
1380 | int p = 0; |
1224 | double volt[nrfiles], volterr[nrfiles], sep[3][nrfiles], seperr[3][nrfiles]; |
1381 | double volt[nrfiles], volterr[nrfiles], sep[3][nrfiles], seperr[3][nrfiles]; |
1225 | int first = 1; |
1382 | int first = 1; |
1226 | 1383 | ||
1227 | FILE *fp; |
1384 | FILE *fp; |
1228 | remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp); |
1385 | remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp); |
1229 | sprintf(paramname, "%s_fitresult.txt", ctemp); |
1386 | sprintf(paramname, "%s_fitresult.txt", ctemp); |
1230 | fp = fopen(paramname, "w"); |
1387 | fp = fopen(paramname, "w"); |
1231 | fclose(fp); |
1388 | fclose(fp); |
1232 | 1389 | ||
1233 | int peaklimit = minPeak->widgetNE[0]->GetNumber()+1; // +1 to account for the pedestal peak |
1390 | int peaklimit = minPeak->widgetNE[0]->GetNumber()+1; // +1 to account for the pedestal peak |
1234 | 1391 | ||
1235 | // Zero the parameter values |
1392 | // Zero the parameter values |
1236 | for(int i = 0; i < nrfiles; i++) |
1393 | for(int i = 0; i < nrfiles; i++) |
1237 | { |
1394 | { |
1238 | volt[i] = 0; volterr[i] = 0; |
1395 | volt[i] = 0; volterr[i] = 0; |
1239 | sep[0][i] = 0; sep[1][i] = 0; sep[2][i] = 0; |
1396 | sep[0][i] = 0; sep[1][i] = 0; sep[2][i] = 0; |
1240 | seperr[0][i] = 0; seperr[1][i] = 0; seperr[2][i] = 0; |
1397 | seperr[0][i] = 0; seperr[1][i] = 0; seperr[2][i] = 0; |
1241 | if(i < 20) { meansel[i] = 0; meanselerr[i] = 0; sigmasel[i] = 0; } |
1398 | if(i < 20) { meansel[i] = 0; meanselerr[i] = 0; sigmasel[i] = 0; } |
1242 | } |
1399 | } |
1243 | 1400 | ||
1244 | float progVal = 0; |
1401 | float progVal = 0; |
1245 | analysisProgress->widgetPB->SetPosition(progVal); |
1402 | analysisProgress->widgetPB->SetPosition(progVal); |
1246 | gVirtualX->Update(1); |
1403 | gVirtualX->Update(1); |
1247 | 1404 | ||
1248 | // Start if we select at least one file |
1405 | // Start if we select at least one file |
1249 | if(nrfiles > 0) |
1406 | if(nrfiles > 0) |
1250 | { |
1407 | { |
1251 | for(int i = 0; i < (int)nrfiles; i++) |
1408 | for(int i = 0; i < (int)nrfiles; i++) |
1252 | { |
1409 | { |
1253 | if( (nrfiles == 1) || (!multiSelect->widgetChBox[0]->IsDown()) ) |
1410 | if( (nrfiles == 1) || (!multiSelect->widgetChBox[0]->IsDown()) ) |
1254 | { |
1411 | { |
1255 | printf("BreakdownVolt(): Only one file selected. Not running analysis, just showing the fit.\n"); |
1412 | printf("BreakdownVolt(): Only one file selected. Not running analysis, just showing the fit.\n"); |
1256 | 1413 | ||
1257 | // Replot the spectrum on analysisCanvas and do not close the input file |
1414 | // Replot the spectrum on analysisCanvas and do not close the input file |
1258 | DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1); |
1415 | DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1); |
1259 | analysisCanvas->GetCanvas()->Modified(); |
1416 | analysisCanvas->GetCanvas()->Modified(); |
1260 | analysisCanvas->GetCanvas()->Update(); |
1417 | analysisCanvas->GetCanvas()->Update(); |
1261 | 1418 | ||
1262 | // Get the spectrum |
1419 | // Get the spectrum |
1263 | histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname); |
1420 | histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname); |
1264 | npeaks = 15; |
1421 | npeaks = 15; |
1265 | double par[300]; |
1422 | double par[300]; |
1266 | spec = new TSpectrum(npeaks); |
1423 | spec = new TSpectrum(npeaks); |
1267 | // Find spectrum background |
1424 | // Find spectrum background |
1268 | histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same"); |
1425 | histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same"); |
1269 | // Clone histogram and subtract background from it if we select that option |
1426 | // Clone histogram and subtract background from it if we select that option |
1270 | h2 = (TH1F*)histtemp->Clone("h2"); |
1427 | h2 = (TH1F*)histtemp->Clone("h2"); |
1271 | if(fitChecks->widgetChBox[0]->IsDown()) |
1428 | if(fitChecks->widgetChBox[0]->IsDown()) |
1272 | h2->Add(histback, -1); |
1429 | h2->Add(histback, -1); |
1273 | // Search for the peaks |
1430 | // Search for the peaks |
1274 | int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() ); |
1431 | int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() ); |
1275 | printf("PhotonMu(): Found %d candidates to fit.\n",found); |
1432 | printf("PhotonMu(): Found %d candidates to fit.\n",found); |
1276 | npeaks = found; |
1433 | npeaks = found; |
1277 | 1434 | ||
1278 | // Set initial peak parameters |
1435 | // Set initial peak parameters |
1279 | xpeaks = spec->GetPositionX(); |
1436 | xpeaks = spec->GetPositionX(); |
1280 | for(j = 0; j < found; j++) |
1437 | for(j = 0; j < found; j++) |
1281 | { |
1438 | { |
1282 | float xp = xpeaks[j]; |
1439 | float xp = xpeaks[j]; |
1283 | int bin = h2->GetXaxis()->FindBin(xp); |
1440 | int bin = h2->GetXaxis()->FindBin(xp); |
1284 | float yp = h2->GetBinContent(bin); |
1441 | float yp = h2->GetBinContent(bin); |
1285 | par[3*j] = yp; |
1442 | par[3*j] = yp; |
1286 | par[3*j+1] = xp; |
1443 | par[3*j+1] = xp; |
1287 | par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber(); |
1444 | par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber(); |
1288 | } |
1445 | } |
1289 | 1446 | ||
1290 | // Fit the histogram |
1447 | // Fit the histogram |
1291 | fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks); |
1448 | fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks); |
1292 | TVirtualFitter::Fitter(histtemp, 3*npeaks); |
1449 | TVirtualFitter::Fitter(histtemp, 3*npeaks); |
1293 | fit->SetParameters(par); |
1450 | fit->SetParameters(par); |
1294 | fit->SetNpx(300); |
1451 | fit->SetNpx(300); |
1295 | h2->Fit("fit","Q"); |
1452 | h2->Fit("fit","Q"); |
1296 | // Get the fitted parameters |
1453 | // Get the fitted parameters |
1297 | fittingfunc = h2->GetFunction("fit"); |
1454 | fittingfunc = h2->GetFunction("fit"); |
1298 | fparam = fittingfunc->GetParameters(); |
1455 | fparam = fittingfunc->GetParameters(); |
1299 | fparamerr = fittingfunc->GetParErrors(); |
1456 | fparamerr = fittingfunc->GetParErrors(); |
1300 | 1457 | ||
1301 | // Gather the parameters (mean peak value for now) |
1458 | // Gather the parameters (mean peak value for now) |
1302 | int j = 1; |
1459 | int j = 1; |
1303 | int nrfit = 0; |
1460 | int nrfit = 0; |
1304 | while(1) |
1461 | while(1) |
1305 | { |
1462 | { |
1306 | if( (fparam[j] < 1.E-30) || (nrfit > 8) ) |
1463 | if( (fparam[j] < 1.E-30) || (nrfit > 8) ) |
1307 | break; |
1464 | break; |
1308 | else |
1465 | else |
1309 | { |
1466 | { |
1310 | // Check if pedestal is above the lower limit and sigma is smaller than the mean |
1467 | // Check if pedestal is above the lower limit and sigma is smaller than the mean |
1311 | if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) ) |
1468 | if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) ) |
1312 | { |
1469 | { |
1313 | // With the additional ADC offset, we can shift the mean values slightly, so they are not close to the X.5, but to the X.0 values |
1470 | // With the additional ADC offset, we can shift the mean values slightly, so they are not close to the X.5, but to the X.0 values |
1314 | meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber()); |
1471 | meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber()); |
- | 1472 | meanselerr[nrfit] = fparamerr[j]; |
|
1315 | sigmasel[nrfit] = fparam[j+1]; |
1473 | sigmasel[nrfit] = fparam[j+1]; |
1316 | nrfit++; |
1474 | nrfit++; |
1317 | } |
1475 | } |
1318 | } |
1476 | } |
1319 | 1477 | ||
1320 | j+=3; |
1478 | j+=3; |
1321 | } |
1479 | } |
1322 | TMath::Sort(nrfit, meansel, sortindex, kFALSE); |
1480 | TMath::Sort(nrfit, meansel, sortindex, kFALSE); |
1323 | 1481 | ||
1324 | fittingfunc->Draw("SAME"); |
1482 | fittingfunc->Draw("SAME"); |
1325 | analysisCanvas->GetCanvas()->Modified(); |
1483 | analysisCanvas->GetCanvas()->Modified(); |
1326 | analysisCanvas->GetCanvas()->Update(); |
1484 | analysisCanvas->GetCanvas()->Update(); |
1327 | 1485 | ||
1328 | meanparam = meansel[sortindex[0]]; |
1486 | meanparam = meansel[sortindex[0]]; |
1329 | meanparamerr = meanselerr[sortindex[0]]; |
1487 | meanparamerr = meanselerr[sortindex[0]]; |
1330 | paramsigma = sigmasel[sortindex[0]]; |
1488 | paramsigma = sigmasel[sortindex[0]]; |
1331 | 1489 | ||
1332 | for(j = 0; j < nrfit; j++) |
1490 | for(j = 0; j < nrfit; j++) |
1333 | printf("BreakdownVolt(): %d: peak mean = %lf, peak err = %lf\n", j, meansel[sortindex[j]], meanselerr[sortindex[j]]); |
1491 | printf("BreakdownVolt(): %d: peak mean = %lf, peak err = %lf\n", j, meansel[sortindex[j]], meanselerr[sortindex[j]]); |
1334 | 1492 | ||
1335 | return; |
1493 | return; |
1336 | } |
1494 | } |
1337 | if(files->At(i)) |
1495 | if(files->At(i)) |
1338 | { |
1496 | { |
1339 | // Replot the spectrum on analysisCanvas and do not close the input file |
1497 | // Replot the spectrum on analysisCanvas and do not close the input file |
1340 | DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1); |
1498 | DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1); |
1341 | analysisCanvas->GetCanvas()->Modified(); |
1499 | analysisCanvas->GetCanvas()->Modified(); |
1342 | analysisCanvas->GetCanvas()->Update(); |
1500 | analysisCanvas->GetCanvas()->Update(); |
1343 | 1501 | ||
1344 | // Get the spectrum |
1502 | // Get the spectrum |
1345 | histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname); |
1503 | histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname); |
1346 | npeaks = 15; |
1504 | npeaks = 15; |
1347 | double par[300]; |
1505 | double par[300]; |
1348 | spec = new TSpectrum(npeaks); |
1506 | spec = new TSpectrum(npeaks); |
1349 | // Find spectrum background |
1507 | // Find spectrum background |
1350 | histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same"); |
1508 | histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same"); |
1351 | // Clone histogram and subtract background from it if we select that option |
1509 | // Clone histogram and subtract background from it if we select that option |
1352 | h2 = (TH1F*)histtemp->Clone("h2"); |
1510 | h2 = (TH1F*)histtemp->Clone("h2"); |
1353 | if(fitChecks->widgetChBox[0]->IsDown()) |
1511 | if(fitChecks->widgetChBox[0]->IsDown()) |
1354 | h2->Add(histback, -1); |
1512 | h2->Add(histback, -1); |
1355 | // Search for the peaks |
1513 | // Search for the peaks |
1356 | int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() ); |
1514 | int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() ); |
1357 | printf("BreakdownVolt(): Found %d candidates to fit.\n",found); |
1515 | printf("BreakdownVolt(): Found %d candidates to fit.\n",found); |
1358 | npeaks = found; |
1516 | npeaks = found; |
1359 | 1517 | ||
1360 | // Set initial peak parameters |
1518 | // Set initial peak parameters |
1361 | xpeaks = spec->GetPositionX(); |
1519 | xpeaks = spec->GetPositionX(); |
1362 | for(j = 0; j < found; j++) |
1520 | for(j = 0; j < found; j++) |
1363 | { |
1521 | { |
1364 | float xp = xpeaks[j]; |
1522 | float xp = xpeaks[j]; |
1365 | int bin = h2->GetXaxis()->FindBin(xp); |
1523 | int bin = h2->GetXaxis()->FindBin(xp); |
1366 | float yp = h2->GetBinContent(bin); |
1524 | float yp = h2->GetBinContent(bin); |
1367 | par[3*j] = yp; |
1525 | par[3*j] = yp; |
1368 | par[3*j+1] = xp; |
1526 | par[3*j+1] = xp; |
1369 | par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber(); |
1527 | par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber(); |
1370 | } |
1528 | } |
1371 | 1529 | ||
1372 | // Fit the histogram |
1530 | // Fit the histogram |
1373 | fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks); |
1531 | fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks); |
1374 | TVirtualFitter::Fitter(histtemp, 3*npeaks); |
1532 | TVirtualFitter::Fitter(histtemp, 3*npeaks); |
1375 | fit->SetParameters(par); |
1533 | fit->SetParameters(par); |
1376 | fit->SetNpx(300); |
1534 | fit->SetNpx(300); |
1377 | h2->Fit("fit","Q"); |
1535 | h2->Fit("fit","Q"); |
1378 | // Get the fitted parameters |
1536 | // Get the fitted parameters |
1379 | fittingfunc = h2->GetFunction("fit"); |
1537 | fittingfunc = h2->GetFunction("fit"); |
1380 | if(nrfiles == 1) // TODO: Show the fit if only one file is selected |
1538 | if(nrfiles == 1) // TODO: Show the fit if only one file is selected |
1381 | fittingfunc->Draw("SAME"); |
1539 | fittingfunc->Draw("SAME"); |
1382 | fparam = fittingfunc->GetParameters(); |
1540 | fparam = fittingfunc->GetParameters(); |
1383 | fparamerr = fittingfunc->GetParErrors(); |
1541 | fparamerr = fittingfunc->GetParErrors(); |
1384 | 1542 | ||
1385 | // Gather the parameters (mean peak value for now) |
1543 | // Gather the parameters (mean peak value for now) |
1386 | int j = 1; |
1544 | int j = 1; |
1387 | int nrfit = 0; |
1545 | int nrfit = 0; |
1388 | while(1) |
1546 | while(1) |
1389 | { |
1547 | { |
1390 | if( (fparam[j] < 1.E-30) || (fparamerr[j] < 1.E-10) )// TODO: Maybe not correct for the error |
1548 | if( (fparam[j] < 1.E-30) || (fparamerr[j] < 1.E-10) )// TODO: Maybe not correct for the error |
1391 | break; |
1549 | break; |
1392 | else |
1550 | else |
1393 | { |
1551 | { |
1394 | // Check if pedestal is above the lower limit and sigma is smaller than the mean |
1552 | // Check if pedestal is above the lower limit and sigma is smaller than the mean |
1395 | if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) ) |
1553 | if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) ) |
1396 | { |
1554 | { |
1397 | meansel[nrfit] = fparam[j]; |
1555 | meansel[nrfit] = fparam[j]; |
1398 | meanselerr[nrfit] = fparamerr[j]; |
1556 | meanselerr[nrfit] = fparamerr[j]; |
1399 | sigmasel[nrfit] = fparam[j+1]; |
1557 | sigmasel[nrfit] = fparam[j+1]; |
1400 | nrfit++; |
1558 | nrfit++; |
1401 | } |
1559 | } |
1402 | } |
1560 | } |
1403 | 1561 | ||
1404 | j+=3; |
1562 | j+=3; |
1405 | } |
1563 | } |
1406 | TMath::Sort(nrfit, meansel, sortindex, kFALSE); |
1564 | TMath::Sort(nrfit, meansel, sortindex, kFALSE); |
1407 | 1565 | ||
1408 | meanparam = meansel[sortindex[0]]; |
1566 | meanparam = meansel[sortindex[0]]; |
1409 | meanparamerr = meanselerr[sortindex[0]]; |
1567 | meanparamerr = meanselerr[sortindex[0]]; |
1410 | paramsigma = sigmasel[sortindex[0]]; |
1568 | paramsigma = sigmasel[sortindex[0]]; |
1411 | 1569 | ||
1412 | for(j = 0; j < nrfit; j++) |
1570 | for(j = 0; j < nrfit; j++) |
1413 | { |
1571 | { |
1414 | if(DBGSIG) |
1572 | if(DBGSIG) |
1415 | printf("BreakdownVolt(): %d: peak mean = %lf, peak err = %lf\n", j, meansel[sortindex[j]], meanselerr[sortindex[j]]); |
1573 | printf("BreakdownVolt(): %d: peak mean = %lf, peak err = %lf\n", j, meansel[sortindex[j]], meanselerr[sortindex[j]]); |
1416 | } |
1574 | } |
1417 | 1575 | ||
1418 | // Delete the opened histogram and spectrum |
1576 | // Delete the opened histogram and spectrum |
1419 | delete spec; |
1577 | delete spec; |
1420 | delete inroot; |
1578 | delete inroot; |
1421 | 1579 | ||
1422 | // Open the input file and read header, ADC and TDC values |
1580 | // Open the input file and read header, ADC and TDC values |
1423 | sprintf(ctemp, "%s", files->At(i)->GetTitle()); |
1581 | sprintf(ctemp, "%s", files->At(i)->GetTitle()); |
1424 | inroot = new TFile(ctemp, "READ"); |
1582 | inroot = new TFile(ctemp, "READ"); |
1425 | 1583 | ||
1426 | inroot->GetObject("header_data", header_data); |
1584 | inroot->GetObject("header_data", header_data); |
1427 | inroot->GetObject("meas_data", meas_data); |
1585 | inroot->GetObject("meas_data", meas_data); |
1428 | 1586 | ||
1429 | // Reading the header |
1587 | // Reading the header |
1430 | if( header_data->FindBranch("biasvolt") ) |
1588 | if( header_data->FindBranch("biasvolt") ) |
1431 | { |
1589 | { |
1432 | header_data->SetBranchAddress("biasvolt", &evtheader.biasvolt); |
1590 | header_data->SetBranchAddress("biasvolt", &evtheader.biasvolt); |
1433 | header_data->GetEntry(0); |
1591 | header_data->GetEntry(0); |
1434 | } |
1592 | } |
1435 | else |
1593 | else |
1436 | { |
1594 | { |
1437 | printf("BreakdownVolt(): Error! Selected file has no bias voltage header value. Please edit header to add the bias voltage header value.\n"); |
1595 | printf("BreakdownVolt(): Error! Selected file has no bias voltage header value. Please edit header to add the bias voltage header value.\n"); |
1438 | break; |
1596 | break; |
1439 | } |
1597 | } |
1440 | 1598 | ||
1441 | // h2->SetStats(0); |
1599 | // h2->SetStats(0); |
1442 | 1600 | ||
1443 | analysisCanvas->GetCanvas()->Modified(); |
1601 | analysisCanvas->GetCanvas()->Modified(); |
1444 | analysisCanvas->GetCanvas()->Update(); |
1602 | analysisCanvas->GetCanvas()->Update(); |
1445 | 1603 | ||
1446 | // Save each fitting plot |
1604 | // Save each fitting plot |
1447 | if(fitChecks->widgetChBox[1]->IsDown()) // TODO: Check if this works |
1605 | if(fitChecks->widgetChBox[1]->IsDown()) // TODO: Check if this works |
1448 | { |
1606 | { |
1449 | remove_ext((char*)files->At(i)->GetTitle(), ctemp); |
1607 | remove_ext((char*)files->At(i)->GetTitle(), ctemp); |
1450 | sprintf(exportname, "%s_fit.pdf", ctemp); |
1608 | sprintf(exportname, "%s_fit.pdf", ctemp); |
1451 | analysisCanvas->GetCanvas()->SaveAs(exportname); |
1609 | analysisCanvas->GetCanvas()->SaveAs(exportname); |
1452 | } |
1610 | } |
1453 | 1611 | ||
1454 | // Calculate the separation between peaks |
1612 | // Calculate the separation between peaks |
1455 | volt[p] = evtheader.biasvolt; |
1613 | volt[p] = evtheader.biasvolt; |
1456 | volterr[p] = 1.e-5; |
1614 | volterr[p] = 1.e-5; |
1457 | 1615 | ||
1458 | if(nrfit == 3) |
1616 | if(nrfit == 3) |
1459 | { |
1617 | { |
1460 | sep[0][p] = meansel[sortindex[2]] - meansel[sortindex[1]]; |
1618 | sep[0][p] = meansel[sortindex[2]] - meansel[sortindex[1]]; |
1461 | seperr[0][p] = TMath::Abs(meanselerr[sortindex[2]]) + TMath::Abs(meanselerr[sortindex[1]]); |
1619 | seperr[0][p] = TMath::Abs(meanselerr[sortindex[2]]) + TMath::Abs(meanselerr[sortindex[1]]); |
1462 | 1620 | ||
1463 | exclude = (seperr[0][p]/sep[0][p] < accError->widgetNE[0]->GetNumber()); |
1621 | exclude = (seperr[0][p]/sep[0][p] < accError->widgetNE[0]->GetNumber()); |
1464 | } |
1622 | } |
1465 | else if(nrfit == 4) |
1623 | else if(nrfit == 4) |
1466 | { |
1624 | { |
1467 | sep[0][p] = meansel[sortindex[2]] - meansel[sortindex[1]]; |
1625 | sep[0][p] = meansel[sortindex[2]] - meansel[sortindex[1]]; |
1468 | sep[1][p] = meansel[sortindex[3]] - meansel[sortindex[2]]; |
1626 | sep[1][p] = meansel[sortindex[3]] - meansel[sortindex[2]]; |
1469 | seperr[0][p] = TMath::Abs(meanselerr[sortindex[2]]) + TMath::Abs(meanselerr[sortindex[1]]); |
1627 | seperr[0][p] = TMath::Abs(meanselerr[sortindex[2]]) + TMath::Abs(meanselerr[sortindex[1]]); |
1470 | seperr[1][p] = TMath::Abs(meanselerr[sortindex[3]]) + TMath::Abs(meanselerr[sortindex[2]]); |
1628 | seperr[1][p] = TMath::Abs(meanselerr[sortindex[3]]) + TMath::Abs(meanselerr[sortindex[2]]); |
1471 | 1629 | ||
1472 | exclude = ((seperr[0][p]/sep[0][p] < accError->widgetNE[0]->GetNumber()) && (seperr[1][p]/sep[1][p] < accError->widgetNE[0]->GetNumber())); |
1630 | exclude = ((seperr[0][p]/sep[0][p] < accError->widgetNE[0]->GetNumber()) && (seperr[1][p]/sep[1][p] < accError->widgetNE[0]->GetNumber())); |
1473 | } |
1631 | } |
1474 | else if(nrfit > 4) |
1632 | else if(nrfit > 4) |
1475 | { |
1633 | { |
1476 | sep[0][p] = meansel[sortindex[2]] - meansel[sortindex[1]]; |
1634 | sep[0][p] = meansel[sortindex[2]] - meansel[sortindex[1]]; |
1477 | sep[1][p] = meansel[sortindex[3]] - meansel[sortindex[2]]; |
1635 | sep[1][p] = meansel[sortindex[3]] - meansel[sortindex[2]]; |
1478 | sep[2][p] = meansel[sortindex[4]] - meansel[sortindex[3]]; |
1636 | sep[2][p] = meansel[sortindex[4]] - meansel[sortindex[3]]; |
1479 | seperr[0][p] = TMath::Abs(meanselerr[sortindex[2]]) + TMath::Abs(meanselerr[sortindex[1]]); |
1637 | seperr[0][p] = TMath::Abs(meanselerr[sortindex[2]]) + TMath::Abs(meanselerr[sortindex[1]]); |
1480 | seperr[1][p] = TMath::Abs(meanselerr[sortindex[3]]) + TMath::Abs(meanselerr[sortindex[2]]); |
1638 | seperr[1][p] = TMath::Abs(meanselerr[sortindex[3]]) + TMath::Abs(meanselerr[sortindex[2]]); |
1481 | seperr[2][p] = TMath::Abs(meanselerr[sortindex[4]]) + TMath::Abs(meanselerr[sortindex[3]]); |
1639 | seperr[2][p] = TMath::Abs(meanselerr[sortindex[4]]) + TMath::Abs(meanselerr[sortindex[3]]); |
1482 | 1640 | ||
1483 | exclude = ((seperr[0][p]/sep[0][p] < accError->widgetNE[0]->GetNumber()) && (seperr[1][p]/sep[1][p] < accError->widgetNE[0]->GetNumber()) && (seperr[2][p]/sep[2][p] < accError->widgetNE[0]->GetNumber())); |
1641 | exclude = ((seperr[0][p]/sep[0][p] < accError->widgetNE[0]->GetNumber()) && (seperr[1][p]/sep[1][p] < accError->widgetNE[0]->GetNumber()) && (seperr[2][p]/sep[2][p] < accError->widgetNE[0]->GetNumber())); |
1484 | } |
1642 | } |
1485 | else |
1643 | else |
1486 | { |
1644 | { |
1487 | printf("BreakdownVolt(): The current separation measurements does not fall within acceptable errors.\n"); |
1645 | printf("BreakdownVolt(): The current separation measurements does not fall within acceptable errors.\n"); |
1488 | exclude = false; |
1646 | exclude = false; |
1489 | } |
1647 | } |
- | 1648 | ||
- | 1649 | printf("BreakdownVolt(): Calculated separation between peaks %d and %d: sep = %lf, seperr = %lf\n", (int)peakSepCalc->widgetNE[0]->GetNumber()-1, (int)peakSepCalc->widgetNE[0]->GetNumber(), sep[(int)peakSepCalc->widgetNE[0]->GetNumber()-1][p], seperr[(int)peakSepCalc->widgetNE[0]->GetNumber()-1][p]); |
|
1490 | 1650 | ||
1491 | // Write out parameters to a file |
1651 | // Write out parameters to a file |
1492 | fp = fopen(paramname, "a"); |
1652 | fp = fopen(paramname, "a"); |
1493 | if(exclude) |
1653 | if(exclude) |
1494 | { |
1654 | { |
1495 | if(first == 1) |
1655 | if(first == 1) |
1496 | { |
1656 | { |
1497 | fprintf(fp, "%le\t%d\t", evtheader.biasvolt, nrfit); |
1657 | fprintf(fp, "%le\t%d\t", evtheader.biasvolt, nrfit); |
1498 | for(j = 0; j < nrfit; j++) |
1658 | for(j = 0; j < nrfit; j++) |
1499 | fprintf(fp, "%le\t%le\t", meansel[sortindex[j]], meanselerr[sortindex[j]]); |
1659 | fprintf(fp, "%le\t%le\t", meansel[sortindex[j]], meanselerr[sortindex[j]]); |
1500 | fprintf(fp,"\n"); |
1660 | fprintf(fp,"\n"); |
1501 | first = 0; |
1661 | first = 0; |
1502 | } |
1662 | } |
1503 | p++; |
1663 | p++; |
1504 | } |
1664 | } |
1505 | else |
1665 | else |
1506 | { |
1666 | { |
1507 | if(nrfit == 3) |
1667 | if(nrfit == 3) |
1508 | printf("Point (at %.2lfV) rejected due to too large errors: %lf\n", volt[p], seperr[0][p]/sep[0][p]); |
1668 | printf("Point (at %.2lfV) rejected due to too large errors: %lf\n", volt[p], seperr[0][p]/sep[0][p]); |
1509 | else if(nrfit == 4) |
1669 | else if(nrfit == 4) |
1510 | printf("Point (at %.2lfV) rejected due to too large errors: %lf, %lf\n", volt[p], seperr[0][p]/sep[0][p], seperr[1][p]/sep[1][p]); |
1670 | printf("Point (at %.2lfV) rejected due to too large errors: %lf, %lf\n", volt[p], seperr[0][p]/sep[0][p], seperr[1][p]/sep[1][p]); |
1511 | else if(nrfit > 4) |
1671 | else if(nrfit > 4) |
1512 | printf("Point (at %.2lfV) rejected due to too large errors: %lf, %lf, %lf\n", volt[p], seperr[0][p]/sep[0][p], seperr[1][p]/sep[1][p], seperr[2][p]/sep[2][p]); |
1672 | printf("Point (at %.2lfV) rejected due to too large errors: %lf, %lf, %lf\n", volt[p], seperr[0][p]/sep[0][p], seperr[1][p]/sep[1][p], seperr[2][p]/sep[2][p]); |
1513 | } |
1673 | } |
1514 | fclose(fp); |
1674 | fclose(fp); |
1515 | } |
1675 | } |
1516 | 1676 | ||
1517 | if(nrfiles == 1) break; |
1677 | if(nrfiles == 1) break; |
1518 | first = 1; |
1678 | first = 1; |
1519 | 1679 | ||
1520 | // Update the progress bar |
1680 | // Update the progress bar |
1521 | progVal = (float)(90.00/nrfiles)*i; |
1681 | progVal = (float)(90.00/nrfiles)*i; |
1522 | analysisProgress->widgetPB->SetPosition(progVal); |
1682 | analysisProgress->widgetPB->SetPosition(progVal); |
1523 | gVirtualX->Update(1); |
1683 | gVirtualX->Update(1); |
1524 | } |
1684 | } |
1525 | 1685 | ||
1526 | printf("BreakdownVolt(): %d files were selected.\n", nrfiles); |
1686 | printf("BreakdownVolt(): %d files were selected.\n", nrfiles); |
1527 | printf("BreakdownVolt(): Number of points to plot is %d\n", p); |
1687 | printf("BreakdownVolt(): Number of points to plot is %d\n", p); |
1528 | 1688 | ||
1529 | // Plot and fit breakdown voltage plot |
1689 | // Plot and fit breakdown voltage plot |
1530 | if(edit == 0) |
1690 | if(edit == 0) |
1531 | gCanvas = new TCanvas("canv","canv",1200,900); |
1691 | gCanvas = new TCanvas("canv","canv",1200,900); |
1532 | else if(edit == 1) |
1692 | else if(edit == 1) |
1533 | gCanvas = tempAnalysisCanvas->GetCanvas(); |
1693 | gCanvas = tempAnalysisCanvas->GetCanvas(); |
1534 | gCanvas->cd(); |
1694 | gCanvas->cd(); |
1535 | gCanvas->SetGrid(); |
1695 | gCanvas->SetGrid(); |
1536 | 1696 | ||
1537 | TGraphErrors* bdplot; |
1697 | TGraphErrors* bdplot; |
1538 | k = peakSepCalc->widgetNE[0]->GetNumber(); |
1698 | k = peakSepCalc->widgetNE[0]->GetNumber(); |
1539 | if(k < 4) |
1699 | if(k < 4) |
1540 | bdplot = new TGraphErrors(p, volt, sep[k-1], volterr, seperr[k-1]); |
1700 | bdplot = new TGraphErrors(p, volt, sep[k-1], volterr, seperr[k-1]); |
1541 | else |
1701 | else |
1542 | { |
1702 | { |
1543 | printf(" |
1703 | printf("BreakdownVolt(): Unsupported peak separation selected (%d).\n", k); |
1544 | return; |
1704 | return; |
1545 | } |
1705 | } |
1546 | 1706 | ||
1547 | bdplot->SetMarkerStyle(20); |
1707 | bdplot->SetMarkerStyle(20); |
1548 | bdplot->SetMarkerSize(0.4); |
1708 | bdplot->SetMarkerSize(0.4); |
1549 | bdplot->SetMarkerColor(kBlue); |
1709 | bdplot->SetMarkerColor(kBlue); |
1550 | bdplot->SetLineWidth(1); |
1710 | bdplot->SetLineWidth(1); |
1551 | bdplot->SetLineColor(kBlue); |
1711 | bdplot->SetLineColor(kBlue); |
1552 | bdplot->GetXaxis()->SetLabelSize(0.030); |
1712 | bdplot->GetXaxis()->SetLabelSize(0.030); |
1553 | bdplot->GetXaxis()->CenterTitle(); |
1713 | bdplot->GetXaxis()->CenterTitle(); |
1554 | // bdplot->GetXaxis()->SetRange(-90,90); |
1714 | // bdplot->GetXaxis()->SetRange(-90,90); |
1555 | // bdplot->GetXaxis()->SetRangeUser(-90,90); |
1715 | // bdplot->GetXaxis()->SetRangeUser(-90,90); |
1556 | // bdplot->GetXaxis()->SetLimits(-90,90); |
1716 | // bdplot->GetXaxis()->SetLimits(-90,90); |
1557 | bdplot->GetYaxis()->SetTitleOffset(1.2); |
1717 | bdplot->GetYaxis()->SetTitleOffset(1.2); |
1558 | bdplot->GetYaxis()->SetLabelSize(0.030); |
1718 | bdplot->GetYaxis()->SetLabelSize(0.030); |
1559 | bdplot->GetYaxis()->CenterTitle(); |
1719 | bdplot->GetYaxis()->CenterTitle(); |
1560 | // bdplot->GetYaxis()->SetRangeUser(0., 1.2); |
1720 | // bdplot->GetYaxis()->SetRangeUser(0., 1.2); |
1561 | bdplot->SetName("bdplot"); |
1721 | bdplot->SetName("bdplot"); |
1562 | bdplot->Draw("AP"); |
1722 | bdplot->Draw("AP"); |
1563 | bdplot->SetTitle(";Bias voltage (V);Peak separation"); |
1723 | bdplot->SetTitle(";Bias voltage (V);Peak separation"); |
1564 | 1724 | ||
1565 | // Fit the breakdown voltage plot with a line |
1725 | // Fit the breakdown voltage plot with a line |
1566 | bdplot->Fit("pol1","Q"); |
1726 | bdplot->Fit("pol1","Q"); |
1567 | 1727 | ||
1568 | TF1 *fit1 = bdplot->GetFunction("pol1"); |
1728 | TF1 *fit1 = bdplot->GetFunction("pol1"); |
1569 | meansel[0] = fit1->GetParameter(0); |
1729 | meansel[0] = fit1->GetParameter(0); |
1570 | meanselerr[0] = fit1->GetParError(0); |
1730 | meanselerr[0] = fit1->GetParError(0); |
1571 | meansel[1] = fit1->GetParameter(1); |
1731 | meansel[1] = fit1->GetParameter(1); |
1572 | meanselerr[1] = fit1->GetParError(1); |
1732 | meanselerr[1] = fit1->GetParError(1); |
1573 | 1733 | ||
1574 | meansel[2] = -meansel[0]/meansel[1]; |
1734 | meansel[2] = -meansel[0]/meansel[1]; |
1575 | if(!cleanPlots) |
1735 | if(!cleanPlots) |
1576 | { |
1736 | { |
1577 | sprintf(ctemp, "#splitline{#Delta_{p}(U) = (%.2lf #pm %.2lf)#timesU + (%.2lf #pm %.3lf)}{U_{0} = %.2lf #pm %.3lf}", meansel[0], meanselerr[0], meansel[1], meanselerr[1], meansel[2], meansel[2]*(TMath::Abs(meanselerr[0]/meansel[0]) + TMath::Abs(meanselerr[1]/meansel[1])) ); |
1737 | sprintf(ctemp, "#splitline{#Delta_{p}(U) = (%.2lf #pm %.2lf)#timesU + (%.2lf #pm %.3lf)}{U_{0} = %.2lf #pm %.3lf}", meansel[0], meanselerr[0], meansel[1], meanselerr[1], meansel[2], meansel[2]*(TMath::Abs(meanselerr[0]/meansel[0]) + TMath::Abs(meanselerr[1]/meansel[1])) ); |
1578 | latex = new TLatex(); |
1738 | latex = new TLatex(); |
1579 | latex->SetTextSize(0.039); |
1739 | latex->SetTextSize(0.039); |
1580 | latex->DrawLatex(volt[0], 0.97*sep[0][sortindex[ |
1740 | latex->DrawLatex(volt[0], 0.97*sep[0][sortindex[1]], ctemp); |
1581 | printf("#Delta_{p}(U) = (%.6lf #pm %.8lf)#timesU + (%.6lf #pm %.8lf)}{U_{0} = %.6lf #pm %.8lf\n", meansel[0], meanselerr[0], meansel[1], meanselerr[1], meansel[2], meansel[2]*(TMath::Abs(meanselerr[0]/meansel[0]) + TMath::Abs(meanselerr[1]/meansel[1])) ); |
1741 | printf("#Delta_{p}(U) = (%.6lf #pm %.8lf)#timesU + (%.6lf #pm %.8lf)}{U_{0} = %.6lf #pm %.8lf\n", meansel[0], meanselerr[0], meansel[1], meanselerr[1], meansel[2], meansel[2]*(TMath::Abs(meanselerr[0]/meansel[0]) + TMath::Abs(meanselerr[1]/meansel[1])) ); |
1582 | } |
1742 | } |
1583 | else |
1743 | else |
1584 | printf("#Delta_{p}(U) = (%.2lf #pm %.2lf)#timesU + (%.2lf #pm %.3lf)}{U_{0} = %.2lf #pm %.3lf\n", meansel[0], meanselerr[0], meansel[1], meanselerr[1], meansel[2], meansel[2]*(TMath::Abs(meanselerr[0]/meansel[0]) + TMath::Abs(meanselerr[1]/meansel[1])) ); |
1744 | printf("#Delta_{p}(U) = (%.2lf #pm %.2lf)#timesU + (%.2lf #pm %.3lf)}{U_{0} = %.2lf #pm %.3lf\n", meansel[0], meanselerr[0], meansel[1], meanselerr[1], meansel[2], meansel[2]*(TMath::Abs(meanselerr[0]/meansel[0]) + TMath::Abs(meanselerr[1]/meansel[1])) ); |
1585 | 1745 | ||
1586 | if(edit == 0) |
1746 | if(edit == 0) |
1587 | { |
1747 | { |
1588 | remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp); |
1748 | remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp); |
1589 | sprintf(exportname, "%s_breakdown.pdf", ctemp); |
1749 | sprintf(exportname, "%s_breakdown.pdf", ctemp); |
1590 | gCanvas->SaveAs(exportname); |
1750 | gCanvas->SaveAs(exportname); |
1591 | 1751 | ||
1592 | delete bdplot; |
1752 | delete bdplot; |
1593 | delete gCanvas; |
1753 | delete gCanvas; |
1594 | } |
1754 | } |
1595 | else if(edit == 1) |
1755 | else if(edit == 1) |
1596 | { |
1756 | { |
1597 | gCanvas->Modified(); |
1757 | gCanvas->Modified(); |
1598 | gCanvas->Update(); |
1758 | gCanvas->Update(); |
1599 | } |
1759 | } |
1600 | 1760 | ||
1601 | // Update the progress bar |
1761 | // Update the progress bar |
1602 | analysisProgress->widgetPB->SetPosition(100.); |
1762 | analysisProgress->widgetPB->SetPosition(100.); |
1603 | gVirtualX->Update(1); |
1763 | gVirtualX->Update(1); |
1604 | } |
1764 | } |
1605 | } |
1765 | } |
1606 | 1766 | ||
1607 | void TGAppMainFrame::SurfaceScan(TList *files, int edit) |
1767 | void TGAppMainFrame::SurfaceScan(TList *files, int edit) |
1608 | { |
1768 | { |
1609 | unsigned int nrfiles = files->GetSize(); |
1769 | unsigned int nrfiles = files->GetSize(); |
1610 | char ctemp[1024]; |
1770 | char ctemp[1024]; |
1611 | char exportname[1024]; |
1771 | char exportname[1024]; |
1612 | int j, k = 0, m = 0, n = 0; |
1772 | int j, k = 0, m = 0, n = 0; |
1613 | 1773 | ||
1614 | TTree *header_data, *meas_data; |
1774 | TTree *header_data, *meas_data; |
1615 | double *integralCount; |
1775 | double *integralCount; |
1616 | integralCount = new double[nrfiles]; |
1776 | integralCount = new double[nrfiles]; |
1617 | double *surfx, *surfy; |
1777 | double *surfx, *surfy; |
1618 | surfx = new double[nrfiles]; |
1778 | surfx = new double[nrfiles]; |
1619 | surfy = new double[nrfiles]; |
1779 | surfy = new double[nrfiles]; |
1620 | double xsurfmin = 0, ysurfmin = 0; |
1780 | double xsurfmin = 0, ysurfmin = 0; |
1621 | int nrentries; |
1781 | int nrentries; |
1622 | double minInteg, maxInteg; |
1782 | double minInteg, maxInteg; |
1623 | bool norm = surfScanOpt->widgetChBox[0]->IsDown(); |
1783 | bool norm = surfScanOpt->widgetChBox[0]->IsDown(); |
1624 | double curyval; |
1784 | double curyval; |
1625 | // bool edge2d = false; |
1785 | // bool edge2d = false; |
1626 | 1786 | ||
1627 | TCanvas *gCanvas; |
1787 | TCanvas *gCanvas; |
1628 | 1788 | ||
1629 | float progVal = 0; |
1789 | float progVal = 0; |
1630 | analysisProgress->widgetPB->SetPosition(progVal); |
1790 | analysisProgress->widgetPB->SetPosition(progVal); |
1631 | gVirtualX->Update(1); |
1791 | gVirtualX->Update(1); |
1632 | 1792 | ||
1633 | // Zero the integral count and accumulated vaules |
1793 | // Zero the integral count and accumulated vaules |
1634 | for(int i = 0; i < (int)nrfiles; i++) integralCount[i] = 0; |
1794 | for(int i = 0; i < (int)nrfiles; i++) integralCount[i] = 0; |
1635 | 1795 | ||
1636 | // Start if we select at more than one file |
1796 | // Start if we select at more than one file |
1637 | if( (nrfiles > 1) && (multiSelect->widgetChBox[0]->IsOn()) ) |
1797 | if( (nrfiles > 1) && (multiSelect->widgetChBox[0]->IsOn()) ) |
1638 | { |
1798 | { |
1639 | printf("SurfaceScan(): Creating a surface plot. Please wait...\n"); |
1799 | printf("SurfaceScan(): Creating a surface plot. Please wait...\n"); |
1640 | for(int i = 0; i < (int)nrfiles; i++) |
1800 | for(int i = 0; i < (int)nrfiles; i++) |
1641 | { |
1801 | { |
1642 | if(files->At(i)) |
1802 | if(files->At(i)) |
1643 | { |
1803 | { |
1644 | n++; |
1804 | n++; |
1645 | // Read the X and Y positions from header and ADC and TDC values from the measurements |
1805 | // Read the X and Y positions from header and ADC and TDC values from the measurements |
1646 | sprintf(ctemp, "%s", files->At(i)->GetTitle()); |
1806 | sprintf(ctemp, "%s", files->At(i)->GetTitle()); |
1647 | inroot = new TFile(ctemp, "READ"); |
1807 | inroot = new TFile(ctemp, "READ"); |
1648 | 1808 | ||
1649 | inroot->GetObject("header_data", header_data); |
1809 | inroot->GetObject("header_data", header_data); |
1650 | inroot->GetObject("meas_data", meas_data); |
1810 | inroot->GetObject("meas_data", meas_data); |
1651 | 1811 | ||
1652 | header_data->SetBranchAddress("xpos", &evtheader.xpos); |
1812 | header_data->SetBranchAddress("xpos", &evtheader.xpos); |
1653 | header_data->GetEntry(0); |
1813 | header_data->GetEntry(0); |
1654 | header_data->SetBranchAddress("ypos", &evtheader.ypos); |
1814 | header_data->SetBranchAddress("ypos", &evtheader.ypos); |
1655 | header_data->GetEntry(0); |
1815 | header_data->GetEntry(0); |
1656 | 1816 | ||
1657 | char rdc[256]; |
1817 | char rdc[256]; |
1658 | j = selectCh->widgetNE[0]->GetNumber(); |
1818 | j = selectCh->widgetNE[0]->GetNumber(); |
1659 | double rangetdc[2]; |
1819 | double rangetdc[2]; |
1660 | rangetdc[0] = tdcRange->widgetNE[0]->GetNumber(); |
1820 | rangetdc[0] = tdcRange->widgetNE[0]->GetNumber(); |
1661 | rangetdc[1] = tdcRange->widgetNE[1]->GetNumber(); |
1821 | rangetdc[1] = tdcRange->widgetNE[1]->GetNumber(); |
1662 | 1822 | ||
1663 | k = 0; |
1823 | k = 0; |
1664 | m = 0; |
1824 | m = 0; |
1665 | 1825 | ||
1666 | // Reading the data |
1826 | // Reading the data |
1667 | for(int e = 0; e < meas_data->GetEntries(); e++) |
1827 | for(int e = 0; e < meas_data->GetEntries(); e++) |
1668 | { |
1828 | { |
1669 | sprintf(rdc, "ADC%d", j); |
1829 | sprintf(rdc, "ADC%d", j); |
1670 | meas_data->SetBranchAddress(rdc, &evtdata.adcdata[j]); |
1830 | meas_data->SetBranchAddress(rdc, &evtdata.adcdata[j]); |
1671 | meas_data->GetEntry(e); |
1831 | meas_data->GetEntry(e); |
1672 | 1832 | ||
1673 | sprintf(rdc, "TDC%d", j); |
1833 | sprintf(rdc, "TDC%d", j); |
1674 | meas_data->SetBranchAddress(rdc, &evtdata.tdcdata[j]); |
1834 | meas_data->SetBranchAddress(rdc, &evtdata.tdcdata[j]); |
1675 | meas_data->GetEntry(e); |
1835 | meas_data->GetEntry(e); |
1676 | 1836 | ||
1677 | // Use data point only if it is inside the TDC window |
1837 | // Use data point only if it is inside the TDC window |
1678 | if( ((double)evtdata.tdcdata[j]/tdctimeconversion >= rangetdc[0]) && ((double)evtdata.tdcdata[j]/tdctimeconversion <= rangetdc[1]) ) |
1838 | if( ((double)evtdata.tdcdata[j]/tdctimeconversion >= rangetdc[0]) && ((double)evtdata.tdcdata[j]/tdctimeconversion <= rangetdc[1]) ) |
1679 | { |
1839 | { |
1680 | k++; |
1840 | k++; |
1681 | m += evtdata.adcdata[j]; |
1841 | m += evtdata.adcdata[j]; |
1682 | } |
1842 | } |
1683 | } |
1843 | } |
1684 | 1844 | ||
1685 | // X, Y and Z values from each file (table units or microns) |
1845 | // X, Y and Z values from each file (table units or microns) |
1686 | if(posUnitsPlot->widgetCB->GetSelected() == 0) |
1846 | if(posUnitsPlot->widgetCB->GetSelected() == 0) |
1687 | { |
1847 | { |
1688 | if(n == 1) |
1848 | if(n == 1) |
1689 | { |
1849 | { |
1690 | xsurfmin = (double)(evtheader.xpos); |
1850 | xsurfmin = (double)(evtheader.xpos); |
1691 | ysurfmin = (double)(evtheader.ypos); |
1851 | ysurfmin = (double)(evtheader.ypos); |
1692 | } |
1852 | } |
1693 | 1853 | ||
1694 | surfx[i] = (double)(evtheader.xpos); |
1854 | surfx[i] = (double)(evtheader.xpos); |
1695 | surfy[i] = (double)(evtheader.ypos); |
1855 | surfy[i] = (double)(evtheader.ypos); |
1696 | } |
1856 | } |
1697 | else if(posUnitsPlot->widgetCB->GetSelected() == 1) |
1857 | else if(posUnitsPlot->widgetCB->GetSelected() == 1) |
1698 | { |
1858 | { |
1699 | if(n == 1) |
1859 | if(n == 1) |
1700 | { |
1860 | { |
1701 | xsurfmin = (double)(evtheader.xpos*lenconversion); |
1861 | xsurfmin = (double)(evtheader.xpos*lenconversion); |
1702 | ysurfmin = (double)(evtheader.ypos*lenconversion); |
1862 | ysurfmin = (double)(evtheader.ypos*lenconversion); |
1703 | } |
1863 | } |
1704 | 1864 | ||
1705 | surfx[i] = (double)(evtheader.xpos*lenconversion); |
1865 | surfx[i] = (double)(evtheader.xpos*lenconversion); |
1706 | surfy[i] = (double)(evtheader.ypos*lenconversion); |
1866 | surfy[i] = (double)(evtheader.ypos*lenconversion); |
1707 | } |
1867 | } |
1708 | 1868 | ||
1709 | // Save result for later plotting |
1869 | // Save result for later plotting |
1710 | if(norm) |
1870 | if(norm) |
1711 | integralCount[i] += ((double)m)/((double)k); |
1871 | integralCount[i] += ((double)m)/((double)k); |
1712 | else |
1872 | else |
1713 | integralCount[i] += m; |
1873 | integralCount[i] += m; |
1714 | 1874 | ||
1715 | inroot->Close(); |
1875 | inroot->Close(); |
1716 | delete inroot; |
1876 | delete inroot; |
1717 | } |
1877 | } |
1718 | 1878 | ||
1719 | // Update the progress bar |
1879 | // Update the progress bar |
1720 | progVal = (float)(75.00/nrfiles)*i; |
1880 | progVal = (float)(75.00/nrfiles)*i; |
1721 | analysisProgress->widgetPB->SetPosition(progVal); |
1881 | analysisProgress->widgetPB->SetPosition(progVal); |
1722 | gVirtualX->Update(1); |
1882 | gVirtualX->Update(1); |
1723 | } |
1883 | } |
1724 | 1884 | ||
1725 | nrentries = n; |
1885 | nrentries = n; |
1726 | printf("SurfaceScan(): %d files were selected.\n", nrfiles); |
1886 | printf("SurfaceScan(): %d files were selected.\n", nrfiles); |
1727 | 1887 | ||
1728 | // Check for Minimum and Maximum values and normalize to 1, if normalization is selected |
1888 | // Check for Minimum and Maximum values and normalize to 1, if normalization is selected |
1729 | if(norm) |
1889 | if(norm) |
1730 | { |
1890 | { |
1731 | minInteg = TMath::MinElement(nrfiles, integralCount); |
1891 | minInteg = TMath::MinElement(nrfiles, integralCount); |
1732 | for(int i = 0; i < nrfiles; i++) |
1892 | for(int i = 0; i < nrfiles; i++) |
1733 | { |
1893 | { |
1734 | integralCount[i] -= minInteg; |
1894 | integralCount[i] -= minInteg; |
1735 | if(DBGSIG) printf("Subtraction: %lf\n", integralCount[i]); |
1895 | if(DBGSIG) printf("Subtraction: %lf\n", integralCount[i]); |
1736 | } |
1896 | } |
1737 | 1897 | ||
1738 | maxInteg = TMath::MaxElement(nrfiles, integralCount); |
1898 | maxInteg = TMath::MaxElement(nrfiles, integralCount); |
1739 | for(int i = 0; i < nrfiles; i++) |
1899 | for(int i = 0; i < nrfiles; i++) |
1740 | { |
1900 | { |
1741 | integralCount[i] = integralCount[i]/maxInteg; |
1901 | integralCount[i] = integralCount[i]/maxInteg; |
1742 | if(DBGSIG) printf("Normalization: %lf\n", integralCount[i]); |
1902 | if(DBGSIG) printf("Normalization: %lf\n", integralCount[i]); |
1743 | } |
1903 | } |
1744 | } |
1904 | } |
1745 | 1905 | ||
1746 | // Make the 2D surface plot |
1906 | // Make the 2D surface plot |
1747 | if(edit == 0) |
1907 | if(edit == 0) |
1748 | gCanvas = new TCanvas("canv","canv",1100,900); |
1908 | gCanvas = new TCanvas("canv","canv",1100,900); |
1749 | else |
1909 | else |
1750 | gCanvas = tempAnalysisCanvas->GetCanvas(); |
1910 | gCanvas = tempAnalysisCanvas->GetCanvas(); |
1751 | 1911 | ||
1752 | double range[4]; |
1912 | double range[4]; |
1753 | TGraph2D *gScan2D; |
1913 | TGraph2D *gScan2D; |
1754 | gScan2D = new TGraph2D(); |
1914 | gScan2D = new TGraph2D(); |
1755 | gScan2D->SetName("surfscan"); |
- | |
1756 | range[0] = TMath::MinElement(nrfiles, surfx); |
1915 | range[0] = TMath::MinElement(nrfiles, surfx); |
1757 | range[1] = TMath::MaxElement(nrfiles, surfx); |
1916 | range[1] = TMath::MaxElement(nrfiles, surfx); |
1758 | range[2] = TMath::MinElement(nrfiles, surfy); |
1917 | range[2] = TMath::MinElement(nrfiles, surfy); |
1759 | range[3] = TMath::MaxElement(nrfiles, surfy); |
1918 | range[3] = TMath::MaxElement(nrfiles, surfy); |
1760 | 1919 | ||
1761 | for(int i = 0; i < nrfiles; i++) |
1920 | for(int i = 0; i < nrfiles; i++) |
1762 | { |
1921 | { |
1763 | if(DBGSIG) |
1922 | if(DBGSIG) |
1764 | { |
1923 | { |
1765 | if(surfScanOpt->widgetChBox[1]->IsDown()) |
1924 | if(surfScanOpt->widgetChBox[1]->IsDown()) |
1766 | printf("SurfaceScan(): %.2lf\t%.2lf\t%lf\n", surfx[i]-range[0], surfy[i]-range[2], integralCount[i]); |
1925 | printf("SurfaceScan(): %.2lf\t%.2lf\t%lf\n", surfx[i]-range[0], surfy[i]-range[2], integralCount[i]); |
1767 | else |
1926 | else |
1768 | printf("SurfaceScan(): %.2lf\t%.2lf\t%lf\n", surfx[i], surfy[i], integralCount[i]); |
1927 | printf("SurfaceScan(): %.2lf\t%.2lf\t%lf\n", surfx[i], surfy[i], integralCount[i]); |
1769 | } |
1928 | } |
1770 | 1929 | ||
1771 | if(surfScanOpt->widgetChBox[1]->IsDown()) |
1930 | if(surfScanOpt->widgetChBox[1]->IsDown()) |
1772 | gScan2D->SetPoint(i, surfx[i]-range[0], surfy[i]-range[2], integralCount[i]); |
1931 | gScan2D->SetPoint(i, surfx[i]-range[0], surfy[i]-range[2], integralCount[i]); |
1773 | else |
1932 | else |
1774 | gScan2D->SetPoint(i, surfx[i], surfy[i], integralCount[i]); |
1933 | gScan2D->SetPoint(i, surfx[i], surfy[i], integralCount[i]); |
1775 | 1934 | ||
1776 | // Update the progress bar |
1935 | // Update the progress bar |
1777 | progVal = (float)(9.00/nrfiles)*i+90.00; |
1936 | progVal = (float)(9.00/nrfiles)*i+90.00; |
1778 | analysisProgress->widgetPB->SetPosition(progVal); |
1937 | analysisProgress->widgetPB->SetPosition(progVal); |
1779 | gVirtualX->Update(1); |
1938 | gVirtualX->Update(1); |
1780 | } |
1939 | } |
1781 | 1940 | ||
1782 | if(surfScanOpt->widgetChBox[1]->IsDown()) |
1941 | if(surfScanOpt->widgetChBox[1]->IsDown()) |
1783 | { |
1942 | { |
1784 | range[1] -= range[0]; |
1943 | range[1] -= range[0]; |
1785 | range[3] -= range[2]; |
1944 | range[3] -= range[2]; |
1786 | range[0] -= range[0]; |
1945 | range[0] -= range[0]; |
1787 | range[2] -= range[2]; |
1946 | range[2] -= range[2]; |
1788 | } |
1947 | } |
1789 | 1948 | ||
1790 | gCanvas->cd(); |
1949 | gCanvas->cd(); |
1791 | gStyle->SetPalette(1); |
1950 | gStyle->SetPalette(1); |
1792 | gCanvas->SetLeftMargin(0.15); |
1951 | gCanvas->SetLeftMargin(0.15); |
1793 | gCanvas->SetRightMargin(0.126); |
1952 | gCanvas->SetRightMargin(0.126); |
1794 | gCanvas->SetTopMargin(0.077); |
1953 | gCanvas->SetTopMargin(0.077); |
1795 | gScan2D->Draw("COLZ"); |
1954 | gScan2D->Draw("COLZ"); |
1796 | 1955 | ||
1797 | gCanvas->Modified(); |
1956 | gCanvas->Modified(); |
1798 | gCanvas->Update(); |
1957 | gCanvas->Update(); |
1799 | 1958 | ||
1800 | if(!cleanPlots) |
1959 | if(!cleanPlots) |
1801 | { |
1960 | { |
1802 | if(posUnitsPlot->widgetCB->GetSelected() == 0) |
1961 | if(posUnitsPlot->widgetCB->GetSelected() == 0) |
1803 | gScan2D->SetTitle("Surface scan;X [table units];Y [table units]"); |
1962 | gScan2D->SetTitle("Surface scan;X [table units];Y [table units]"); |
1804 | else if(posUnitsPlot->widgetCB->GetSelected() == 1) |
1963 | else if(posUnitsPlot->widgetCB->GetSelected() == 1) |
1805 | gScan2D->SetTitle("Surface scan;X [#mum];Y [#mum]"); |
1964 | gScan2D->SetTitle("Surface scan;X [#mum];Y [#mum]"); |
1806 | } |
1965 | } |
1807 | else |
1966 | else |
1808 | { |
1967 | { |
1809 | if(posUnitsPlot->widgetCB->GetSelected() == 0) |
1968 | if(posUnitsPlot->widgetCB->GetSelected() == 0) |
1810 | gScan2D->SetTitle(";X [table units];Y [table units]"); |
1969 | gScan2D->SetTitle(";X [table units];Y [table units]"); |
1811 | else if(posUnitsPlot->widgetCB->GetSelected() == 1) |
1970 | else if(posUnitsPlot->widgetCB->GetSelected() == 1) |
1812 | gScan2D->SetTitle(";X [#mum];Y [#mum]"); |
1971 | gScan2D->SetTitle(";X [#mum];Y [#mum]"); |
1813 | } |
1972 | } |
1814 | /* TGaxis *xax = (TGaxis*)gScan2D->GetXaxis(); |
1973 | /* TGaxis *xax = (TGaxis*)gScan2D->GetXaxis(); |
1815 | xax->SetMaxDigits(4); |
1974 | xax->SetMaxDigits(4); |
1816 | TGaxis *yax = (TGaxis*)gScan2D->GetYaxis(); |
1975 | TGaxis *yax = (TGaxis*)gScan2D->GetYaxis(); |
1817 | yax->SetMaxDigits(4);*/ |
1976 | yax->SetMaxDigits(4);*/ |
1818 | 1977 | ||
1819 | gScan2D->SetNpx((int)resolSurf->widgetNE[0]->GetNumber()); |
1978 | gScan2D->SetNpx((int)resolSurf->widgetNE[0]->GetNumber()); |
1820 | gScan2D->SetNpy((int)resolSurf->widgetNE[1]->GetNumber()); |
1979 | gScan2D->SetNpy((int)resolSurf->widgetNE[1]->GetNumber()); |
1821 | 1980 | ||
1822 | gCanvas->Modified(); |
1981 | gCanvas->Modified(); |
1823 | gCanvas->Update(); |
1982 | gCanvas->Update(); |
1824 | 1983 | ||
1825 | gScan2D->GetXaxis()->SetTitleOffset(1.3); |
1984 | gScan2D->GetXaxis()->SetTitleOffset(1.3); |
1826 | gScan2D->GetXaxis()->CenterTitle(kTRUE); |
1985 | gScan2D->GetXaxis()->CenterTitle(kTRUE); |
1827 | gScan2D->GetXaxis()->SetLabelSize(0.027); |
1986 | gScan2D->GetXaxis()->SetLabelSize(0.027); |
1828 | gScan2D->GetXaxis()->SetLabelOffset(0.02); |
1987 | gScan2D->GetXaxis()->SetLabelOffset(0.02); |
1829 | gScan2D->GetXaxis()->SetRangeUser(range[0], range[1]); |
1988 | gScan2D->GetXaxis()->SetRangeUser(range[0], range[1]); |
1830 | gScan2D->GetXaxis()->SetNoExponent(kTRUE); |
1989 | gScan2D->GetXaxis()->SetNoExponent(kTRUE); |
1831 | gScan2D->GetYaxis()->SetTitleOffset(1.9); |
1990 | gScan2D->GetYaxis()->SetTitleOffset(1.9); |
1832 | gScan2D->GetYaxis()->CenterTitle(kTRUE); |
1991 | gScan2D->GetYaxis()->CenterTitle(kTRUE); |
1833 | gScan2D->GetYaxis()->SetLabelSize(0.027); |
1992 | gScan2D->GetYaxis()->SetLabelSize(0.027); |
1834 | gScan2D->GetYaxis()->SetLabelOffset(0.02); |
1993 | gScan2D->GetYaxis()->SetLabelOffset(0.02); |
1835 | gScan2D->GetYaxis()->SetRangeUser(range[2], range[3]); |
1994 | gScan2D->GetYaxis()->SetRangeUser(range[2], range[3]); |
1836 | gScan2D->GetYaxis()->SetNoExponent(kTRUE); |
1995 | gScan2D->GetYaxis()->SetNoExponent(kTRUE); |
1837 | 1996 | ||
1838 | gCanvas->Modified(); |
1997 | gCanvas->Modified(); |
1839 | gCanvas->Update(); |
1998 | gCanvas->Update(); |
1840 | 1999 | ||
1841 | remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp); |
2000 | remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp); |
1842 | sprintf(exportname, "%s_surfscan.pdf", ctemp); |
2001 | sprintf(exportname, "%s_surfscan.pdf", ctemp); |
1843 | gCanvas->SaveAs(exportname); |
2002 | gCanvas->SaveAs(exportname); |
1844 | 2003 | ||
1845 | // Update the progress bar |
2004 | // Update the progress bar |
1846 | analysisProgress->widgetPB->SetPosition(100.0); |
2005 | analysisProgress->widgetPB->SetPosition(100.0); |
1847 | gVirtualX->Update(1); |
2006 | gVirtualX->Update(1); |
1848 | 2007 | ||
1849 | if(edit == 0) |
2008 | if(edit == 0) |
1850 | { |
2009 | { |
1851 | delete gScan2D; |
2010 | delete gScan2D; |
1852 | delete gCanvas; |
2011 | delete gCanvas; |
1853 | } |
2012 | } |
1854 | else if(edit == 1) |
2013 | else if(edit == 1) |
1855 | { |
2014 | { |
1856 | gCanvas->Modified(); |
2015 | gCanvas->Modified(); |
1857 | gCanvas->Update(); |
2016 | gCanvas->Update(); |
1858 | - | ||
1859 | UpdateIntegrateSurface(-1); |
- | |
1860 | } |
2017 | } |
1861 | } |
2018 | } |
1862 | } |
2019 | } |
1863 | 2020 |