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