Subversion Repositories f9daq

Rev

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