Subversion Repositories f9daq

Rev

Rev 173 | Go to most recent revision | 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
}