Subversion Repositories f9daq

Rev

Rev 146 | Rev 172 | Go to most recent revision | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 146 Rev 167
1
#include "../include/sipmscan.h"
1
#include "../include/sipmscan.h"
2
#include "../include/workstation.h"
2
#include "../include/workstation.h"
3
 
3
 
4
#include <stdio.h>
4
#include <stdio.h>
5
#include <stdlib.h>
5
#include <stdlib.h>
6
 
6
 
7
// Peak detection function
7
// Peak detection function
8
int npeaks = 20;
8
int npeaks = 20;
9
double FindPeaks(double *x, double *par)
9
double FindPeaks(double *x, double *par)
10
{
10
{
11
   double result = 0;
11
   double result = 0;
12
   for(int i = 0; i < npeaks; i++)
12
   for(int i = 0; i < npeaks; i++)
13
   {
13
   {
14
      double norm = par[3*i];
14
      double norm = par[3*i];
15
      double mean = par[3*i+1];
15
      double mean = par[3*i+1];
16
      double sigma = par[3*i+2];
16
      double sigma = par[3*i+2];
17
      result += norm*TMath::Gaus(x[0], mean, sigma);
17
      result += norm*TMath::Gaus(x[0], mean, sigma);
18
   }
18
   }
19
   return result;
19
   return result;
20
}
20
}
21
 
21
 
22
// File browser for selecting the dark run histogram
22
// File browser for selecting the dark run histogram
23
void TGAppMainFrame::SelectDarkHist()
23
void TGAppMainFrame::SelectDarkHist()
24
{
24
{
25
   TGFileInfo file_info;
25
   TGFileInfo file_info;
26
   const char *filetypes[] = {"Histograms",histextall,0,0};
26
   const char *filetypes[] = {"Histograms",histextall,0,0};
27
   char *cTemp;
27
   char *cTemp;
28
   file_info.fFileTypes = filetypes;
28
   file_info.fFileTypes = filetypes;
29
   cTemp = new char[1024];
29
//   cTemp = new char[1024];
30
   sprintf(cTemp, "%s/results", rootdir);
30
//   sprintf(cTemp, "%s/results", rootdir);
31
   file_info.fIniDir = StrDup(cTemp);
31
//   file_info.fIniDir = StrDup(cTemp);
-
 
32
   file_info.fIniDir = StrDup(currentOpenDir);
32
   file_info.fMultipleSelection = kFALSE;
33
   file_info.fMultipleSelection = kFALSE;
33
   new TGFileDialog(gClient->GetDefaultRoot(), fMain, kFDOpen, &file_info);
34
   new TGFileDialog(gClient->GetDefaultRoot(), fMain, kFDOpen, &file_info);
34
   delete[] cTemp;
35
//   delete[] cTemp;
35
 
36
 
36
   if(file_info.fFilename != NULL)
37
   if(file_info.fFilename != NULL)
37
   {
38
   {
38
      darkRun->widgetTE->SetText(file_info.fFilename);
39
      darkRun->widgetTE->SetText(file_info.fFilename);
39
      fileList->AddEntry(file_info.fFilename, fileList->GetNumberOfEntries());
40
      fileList->AddEntry(file_info.fFilename, fileList->GetNumberOfEntries());
40
      fileList->Layout();
41
      fileList->Layout();
41
   }
42
   }
42
   else
43
   else
43
      darkRun->widgetTE->SetText("");
44
      darkRun->widgetTE->SetText("");
44
}
45
}
45
 
46
 
46
// Reset analysis defaults for the current analysis type (0 = Integrate spectrum, 1 = Relative PDE, 2 = Breakdown voltage, 3 = Surface scan, 4 = Timing analysis)
47
// Reset analysis defaults for the current analysis type (0 = Integrate spectrum, 1 = Relative PDE, 2 = Breakdown voltage, 3 = Surface scan, 4 = Timing analysis)
47
void TGAppMainFrame::AnalysisDefaults()
48
void TGAppMainFrame::AnalysisDefaults()
48
{
49
{
49
   printf("AnalysisDefaults(): Current analysis tab = %d\n", analTab->GetCurrent());
50
   printf("AnalysisDefaults(): Current analysis tab = %d\n", analTab->GetCurrent());
50
   if(analTab->GetCurrent() == 0)       // Integrate spectrum
51
   if(analTab->GetCurrent() == 0)       // Integrate spectrum
51
   {
52
   {
52
      intSpect->widgetChBox[0]->SetState(kButtonUp);
53
      intSpect->widgetChBox[0]->SetState(kButtonUp);
53
      intSpect->widgetChBox[1]->SetState(kButtonUp);
54
      intSpect->widgetChBox[1]->SetState(kButtonUp);
54
      intSpect->widgetChBox[2]->SetState(kButtonDown);
55
      intSpect->widgetChBox[2]->SetState(kButtonDown);
55
      resol2d->widgetNE[0]->SetNumber(40);
56
      resol2d->widgetNE[0]->SetNumber(40);
56
      resol2d->widgetNE[1]->SetNumber(40);
57
      resol2d->widgetNE[1]->SetNumber(40);
57
   }
58
   }
58
   else if(analTab->GetCurrent() == 1)  // Relative PDE
59
   else if(analTab->GetCurrent() == 1)  // Relative PDE
59
   {
60
   {
60
      relPde->widgetChBox[1]->SetState(kButtonDown);
61
      relPde->widgetChBox[1]->SetState(kButtonDown);
61
      midPeak->widgetChBox[0]->SetState(kButtonUp);
62
      midPeak->widgetChBox[0]->SetState(kButtonUp);
62
      zeroAngle->widgetNE[0]->SetNumber(0.00);
63
      zeroAngle->widgetNE[0]->SetNumber(0.00);
63
   }
64
   }
64
   else if(analTab->GetCurrent() == 2)  // Breakdown voltage
65
   else if(analTab->GetCurrent() == 2)  // Breakdown voltage
65
   {
66
   {
66
      minPeak->widgetNE[0]->SetNumber(2);
67
      minPeak->widgetNE[0]->SetNumber(2);
67
      minPeak->widgetNE[0]->SetNumber(1);
68
      minPeak->widgetNE[0]->SetNumber(1);
68
   }
69
   }
69
   else if(analTab->GetCurrent() == 3)  // Surface scan
70
   else if(analTab->GetCurrent() == 3)  // Surface scan
70
   {
71
   {
71
      surfScanOpt->widgetChBox[0]->SetState(kButtonDown);
72
      surfScanOpt->widgetChBox[0]->SetState(kButtonDown);
72
      surfScanOpt->widgetChBox[1]->SetState(kButtonUp);
73
      surfScanOpt->widgetChBox[1]->SetState(kButtonUp);
73
      resolSurf->widgetNE[0]->SetNumber(40);
74
      resolSurf->widgetNE[0]->SetNumber(40);
74
      resolSurf->widgetNE[1]->SetNumber(40);
75
      resolSurf->widgetNE[1]->SetNumber(40);
75
   }
76
   }
76
}
77
}
77
 
78
 
78
// Analysis handle function
79
// Analysis handle function
79
void TGAppMainFrame::AnalysisHandle(int type)
80
void TGAppMainFrame::AnalysisHandle(int type)
80
{
81
{
81
   TList *files;
82
   TList *files;
82
   bool createTab = true;
83
   bool createTab = true;
83
   int tabid = -1;
84
   int tabid = -1;
84
   int analtab = analTab->GetCurrent();
85
   int analtab = analTab->GetCurrent();
85
 
86
 
86
   int analtype;
87
   int analtype;
87
   if( (analtab == 0) && (!intSpect->widgetChBox[0]->IsDown() && !intSpect->widgetChBox[1]->IsDown()) )
88
   if( (analtab == 0) && (!intSpect->widgetChBox[0]->IsDown() && !intSpect->widgetChBox[1]->IsDown()) )
88
      analtype = 0;
89
      analtype = 0;
89
   else if( (analtab == 0) && (intSpect->widgetChBox[0]->IsDown() || intSpect->widgetChBox[1]->IsDown()) )
90
   else if( (analtab == 0) && (intSpect->widgetChBox[0]->IsDown() || intSpect->widgetChBox[1]->IsDown()) )
90
      analtype = 1;
91
      analtype = 1;
91
   else if( analtab == 1 )
92
   else if( analtab == 1 )
92
      analtype = 2;
93
      analtype = 2;
93
   else if( analtab == 2 )
94
   else if( analtab == 2 )
94
      analtype = 3;
95
      analtype = 3;
95
   else if( analtab == 3 )
96
   else if( analtab == 3 )
96
      analtype = 4;
97
      analtype = 4;
97
 
98
 
98
   // Only integrate spectrum or make relative PDE
99
   // Only integrate spectrum or make relative PDE
99
   if(type == 0)
100
   if(type == 0)
100
   {
101
   {
101
      files = new TList();
102
      files = new TList();
102
      fileList->GetSelectedEntries(files);
103
      fileList->GetSelectedEntries(files);
103
 
104
 
104
      if( analtype == 0 )
105
      if( analtype == 0 )
105
         IntegSpectrum(files, 0, 0);
106
         IntegSpectrum(files, 0, 0);
106
 
107
 
107
      if( intSpect->widgetChBox[0]->IsDown() )
108
      if( intSpect->widgetChBox[0]->IsDown() )
108
         IntegSpectrum(files, 1, 0);
109
         IntegSpectrum(files, 1, 0);
109
 
110
 
110
      if( intSpect->widgetChBox[1]->IsDown() )
111
      if( intSpect->widgetChBox[1]->IsDown() )
111
         IntegSpectrum(files, 2, 0);
112
         IntegSpectrum(files, 2, 0);
112
 
113
 
113
      if( analtype == 2 )
114
      if( analtype == 2 )
114
         PhotonMu(files, 0);
115
         PhotonMu(files, 0);
115
 
116
 
116
      if( analtype == 3 )
117
      if( analtype == 3 )
117
         BreakdownVolt(files, 0);
118
         BreakdownVolt(files, 0);
118
 
119
 
119
      if( analtype == 4 )
120
      if( analtype == 4 )
120
         SurfaceScan(files, 0);
121
         SurfaceScan(files, 0);
121
   }
122
   }
122
   // Integrate spectrum or make relative PDE and open edit window
123
   // Integrate spectrum or make relative PDE and open edit window
123
   else if(type == 1)
124
   else if(type == 1)
124
   {
125
   {
125
      files = new TList();
126
      files = new TList();
126
      fileList->GetSelectedEntries(files);
127
      fileList->GetSelectedEntries(files);
127
 
128
 
128
      // Prepare a new analysis edit tab, if we want to edit plots
129
      // Prepare a new analysis edit tab, if we want to edit plots
129
      for(int i = 0; i < fTab->GetNumberOfTabs(); i++)
130
      for(int i = 0; i < fTab->GetNumberOfTabs(); i++)
130
      {
131
      {
131
         if(strcmp("Analysis edit", fTab->GetTabTab(i)->GetString() ) == 0)
132
         if(strcmp("Analysis edit", fTab->GetTabTab(i)->GetString() ) == 0)
132
         {
133
         {
133
            createTab = false;
134
            createTab = false;
134
            tabid = i;
135
            tabid = i;
135
         }
136
         }
136
         
137
         
137
         if(DBGSIG) printf("AnalysisHandle(): Name of tab = %s\n", fTab->GetTabTab(i)->GetString() );
138
         if(DBGSIG) printf("AnalysisHandle(): Name of tab = %s\n", fTab->GetTabTab(i)->GetString() );
138
      }
139
      }
139
   
140
   
140
      if(files->GetSize() > 0)
141
      if(files->GetSize() > 0)
141
      {
142
      {
142
         TempAnalysisTab(fTab, createTab, &tabid, analtype);
143
         TempAnalysisTab(fTab, createTab, &tabid, analtype);
143
 
144
 
144
         // Integrate spectra
145
         // Integrate spectra
145
         if( analtype == 0 )
146
         if( analtype == 0 )
146
            IntegSpectrum(files, 0, 1);
147
            IntegSpectrum(files, 0, 1);
147
 
148
 
148
         if( intSpect->widgetChBox[0]->IsDown() )
149
         if( intSpect->widgetChBox[0]->IsDown() )
149
            IntegSpectrum(files, 1, 1);
150
            IntegSpectrum(files, 1, 1);
150
 
151
 
151
         if( intSpect->widgetChBox[1]->IsDown() )
152
         if( intSpect->widgetChBox[1]->IsDown() )
152
            IntegSpectrum(files, 2, 1);
153
            IntegSpectrum(files, 2, 1);
153
 
154
 
154
         if( analtype == 2 )
155
         if( analtype == 2 )
155
            PhotonMu(files, 1);
156
            PhotonMu(files, 1);
156
 
157
 
157
         if( analtype == 3 )
158
         if( analtype == 3 )
158
            BreakdownVolt(files, 1);
159
            BreakdownVolt(files, 1);
159
 
160
 
160
         if( analtype == 4 )
161
         if( analtype == 4 )
161
            SurfaceScan(files, 1);
162
            SurfaceScan(files, 1);
162
 
163
 
163
         fTab->SetTab(tabid);
164
         fTab->SetTab(tabid);
164
      }
165
      }
165
   }
166
   }
166
 
167
 
167
   delete files;
168
   delete files;
168
}
169
}
169
 
170
 
170
// Analysis functions ----------------------------------
171
// Analysis functions ----------------------------------
171
void TGAppMainFrame::IntegSpectrum(TList *files, int direction, int edit)
172
void TGAppMainFrame::IntegSpectrum(TList *files, int direction, int edit)
172
{
173
{
173
   unsigned int nrfiles = files->GetSize();
174
   unsigned int nrfiles = files->GetSize();
174
   char ctemp[1024];
175
   char ctemp[1024];
175
   char exportname[1024];
176
   char exportname[1024];
176
   int j, k = 0, m = 0;
177
   int j, k = 0, m = 0;
177
 
178
 
178
   TTree *header_data, *meas_data;
179
   TTree *header_data, *meas_data;
179
   double *integralCount, *integralAcc;
180
   double *integralCount, *integralAcc;
180
   integralCount = new double[nrfiles];
181
   integralCount = new double[nrfiles];
181
   integralAcc = new double[nrfiles];
182
   integralAcc = new double[nrfiles];
182
   double *surfxy, *surfz;
183
   double *surfxy, *surfz;
183
   surfxy = new double[nrfiles];
184
   surfxy = new double[nrfiles];
184
   surfz = new double[nrfiles];
185
   surfz = new double[nrfiles];
185
   double minInteg, maxInteg;
186
   double minInteg, maxInteg;
186
   bool norm = intSpect->widgetChBox[2]->IsDown();
187
   bool norm = intSpect->widgetChBox[2]->IsDown();
187
   double curzval;
188
   double curzval;
188
   bool edge2d = false;
189
   bool edge2d = false;
189
   TCanvas *gCanvas;
190
   TCanvas *gCanvas;
190
 
191
 
191
   float progVal = 0;
192
   float progVal = 0;
192
   analysisProgress->widgetPB->SetPosition(progVal);
193
   analysisProgress->widgetPB->SetPosition(progVal);
193
   gVirtualX->Update(1);
194
   gVirtualX->Update(1);
194
 
195
 
195
   // Zero the integral count and accumulated vaules
196
   // Zero the integral count and accumulated vaules
196
   for(int i = 0; i < (int)nrfiles; i++) {integralCount[i] = 0; integralAcc[i] = 0; }
197
   for(int i = 0; i < (int)nrfiles; i++) {integralCount[i] = 0; integralAcc[i] = 0; }
197
 
198
 
198
   // Start if we select at least one file
199
   // Start if we select at least one file
199
   if(nrfiles > 0)
200
   if(nrfiles > 0)
200
   {
201
   {
201
      for(int i = 0; i < (int)nrfiles; i++)
202
      for(int i = 0; i < (int)nrfiles; i++)
202
      {
203
      {
203
         if(files->At(i))
204
         if(files->At(i))
204
         {
205
         {
205
            // Read the X,Y and Z positions from header and ADC and TDC values from the measurements
206
            // Read the X,Y and Z positions from header and ADC and TDC values from the measurements
206
            sprintf(ctemp, "%s", files->At(i)->GetTitle());
207
            sprintf(ctemp, "%s", files->At(i)->GetTitle());
207
            inroot = new TFile(ctemp, "READ");
208
            inroot = new TFile(ctemp, "READ");
208
         
209
         
209
            inroot->GetObject("header_data", header_data);
210
            inroot->GetObject("header_data", header_data);
210
            inroot->GetObject("meas_data", meas_data);
211
            inroot->GetObject("meas_data", meas_data);
211
         
212
         
212
            header_data->SetBranchAddress("xpos", &evtheader.xpos);
213
            header_data->SetBranchAddress("xpos", &evtheader.xpos);
213
            header_data->GetEntry(0);
214
            header_data->GetEntry(0);
214
            header_data->SetBranchAddress("ypos", &evtheader.ypos);
215
            header_data->SetBranchAddress("ypos", &evtheader.ypos);
215
            header_data->GetEntry(0);
216
            header_data->GetEntry(0);
216
            header_data->SetBranchAddress("zpos", &evtheader.zpos);
217
            header_data->SetBranchAddress("zpos", &evtheader.zpos);
217
            header_data->GetEntry(0);
218
            header_data->GetEntry(0);
218
 
219
 
219
            char rdc[256];
220
            char rdc[256];
220
            j = selectCh->widgetNE[0]->GetNumber();
221
            j = selectCh->widgetNE[0]->GetNumber();
221
            double rangetdc[2];
222
            double rangetdc[2];
222
            rangetdc[0] = tdcRange->widgetNE[0]->GetNumber();
223
            rangetdc[0] = tdcRange->widgetNE[0]->GetNumber();
223
            rangetdc[1] = tdcRange->widgetNE[1]->GetNumber();
224
            rangetdc[1] = tdcRange->widgetNE[1]->GetNumber();
224
   
225
   
225
            k = 0;
226
            k = 0;
226
            m = 0;
227
            m = 0;
227
         
228
         
228
            // Reading the data
229
            // Reading the data
229
            for(int e = 0; e < meas_data->GetEntries(); e++)
230
            for(int e = 0; e < meas_data->GetEntries(); e++)
230
            {
231
            {
231
               sprintf(rdc, "ADC%d", j);
232
               sprintf(rdc, "ADC%d", j);
232
               meas_data->SetBranchAddress(rdc, &evtdata.adcdata[j]);
233
               meas_data->SetBranchAddress(rdc, &evtdata.adcdata[j]);
233
               meas_data->GetEntry(e);
234
               meas_data->GetEntry(e);
234
         
235
         
235
               sprintf(rdc, "TDC%d", j);
236
               sprintf(rdc, "TDC%d", j);
236
               meas_data->SetBranchAddress(rdc, &evtdata.tdcdata[j]);
237
               meas_data->SetBranchAddress(rdc, &evtdata.tdcdata[j]);
237
               meas_data->GetEntry(e);
238
               meas_data->GetEntry(e);
238
   
239
   
239
               // Use data point only if it is inside the TDC window
240
               // Use data point only if it is inside the TDC window
240
               if( ((double)evtdata.tdcdata[j]/tdctimeconversion >= rangetdc[0]) && ((double)evtdata.tdcdata[j]/tdctimeconversion <= rangetdc[1]) )
241
               if( ((double)evtdata.tdcdata[j]/tdctimeconversion >= rangetdc[0]) && ((double)evtdata.tdcdata[j]/tdctimeconversion <= rangetdc[1]) )
241
               {
242
               {
242
                  k++;
243
                  k++;
243
                  m += evtdata.adcdata[j];
244
                  m += evtdata.adcdata[j];
244
               }
245
               }
245
            }
246
            }
246
 
247
 
247
            // X, Y and Z values from each file (table units or microns)
248
            // X, Y and Z values from each file (table units or microns)
248
            if(posUnits->widgetCB->GetSelected() == 0)
249
            if(posUnitsPlot->widgetCB->GetSelected() == 0)
249
            {
250
            {
250
               if(direction == 1)
251
               if(direction == 1)
251
                  surfxy[i] = (double)(evtheader.xpos);
252
                  surfxy[i] = (double)(evtheader.xpos);
252
               else if(direction == 2)
253
               else if(direction == 2)
253
                  surfxy[i] = (double)(evtheader.ypos);
254
                  surfxy[i] = (double)(evtheader.ypos);
254
               surfz[i] = (double)(evtheader.zpos);
255
               surfz[i] = (double)(evtheader.zpos);
255
            }
256
            }
256
            else if(posUnits->widgetCB->GetSelected() == 1)
257
            else if(posUnitsPlot->widgetCB->GetSelected() == 1)
257
            {
258
            {
258
               if(direction == 1)
259
               if(direction == 1)
259
                  surfxy[i] = (double)(evtheader.xpos*lenconversion);
260
                  surfxy[i] = (double)(evtheader.xpos*lenconversion);
260
               else if(direction == 2)
261
               else if(direction == 2)
261
                  surfxy[i] = (double)(evtheader.ypos*lenconversion);
262
                  surfxy[i] = (double)(evtheader.ypos*lenconversion);
262
               surfz[i] = (double)(evtheader.zpos*lenconversion);
263
               surfz[i] = (double)(evtheader.zpos*lenconversion);
263
            }
264
            }
264
 
265
 
265
            // 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
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
266
            if(i == 0) curzval = surfz[i];
267
            if(i == 0) curzval = surfz[i];
267
            else
268
            else
268
            {
269
            {
269
               if(surfz[i] != curzval)
270
               if(surfz[i] != curzval)
270
                  edge2d = true;
271
                  edge2d = true;
271
            }
272
            }
272
 
273
 
273
            // Print the calculated integral, if no X or Y direction edge plots are enabled; otherwise, just save for later plotting
274
            // Print the calculated integral, if no X or Y direction edge plots are enabled; otherwise, just save for later plotting
274
            if(direction == 0)
275
            if(direction == 0)
275
            {
276
            {
276
               if(norm)
277
               if(norm)
277
               {
278
               {
278
                  integralCount[i] += ((double)m)/((double)k);
279
                  integralCount[i] += ((double)m)/((double)k);
279
                  printf("IntegSpectrum(): %s: The integral is: %lf\n", ctemp, integralCount[i]);
280
                  printf("IntegSpectrum(): %s: The integral is: %lf\n", ctemp, integralCount[i]);
280
               }
281
               }
281
               else
282
               else
282
               {
283
               {
283
                  integralCount[i] += m;
284
                  integralCount[i] += m;
284
                  printf("IntegSpectrum(): %s: The integral is: %d\n", ctemp, (int)integralCount[i]);
285
                  printf("IntegSpectrum(): %s: The integral is: %d\n", ctemp, (int)integralCount[i]);
285
               }
286
               }
286
            }
287
            }
287
            else
288
            else
288
            {
289
            {
289
               if(norm)
290
               if(norm)
290
                  integralCount[i] += ((double)m)/((double)k);
291
                  integralCount[i] += ((double)m)/((double)k);
291
               else
292
               else
292
                  integralCount[i] += m;
293
                  integralCount[i] += m;
293
            }
294
            }
294
           
295
           
295
            inroot->Close();
296
            inroot->Close();
296
            delete inroot;
297
            delete inroot;
297
         }
298
         }
298
 
299
 
299
         // Update the progress bar
300
         // Update the progress bar
300
         progVal = (float)(75.00/nrfiles)*i;
301
         progVal = (float)(75.00/nrfiles)*i;
301
         analysisProgress->widgetPB->SetPosition(progVal);
302
         analysisProgress->widgetPB->SetPosition(progVal);
302
         gVirtualX->Update(1);
303
         gVirtualX->Update(1);
303
      }
304
      }
304
 
305
 
305
      printf("IntegSpectrum(): %d files were selected.\n", nrfiles);
306
      printf("IntegSpectrum(): %d files were selected.\n", nrfiles);
306
 
307
 
307
      // If only an integral is needed, do not plot and exit here
308
      // If only an integral is needed, do not plot and exit here
308
      if( direction == 0 )
309
      if( direction == 0 )
309
      {
310
      {
310
         delete[] integralCount;
311
         delete[] integralCount;
311
         delete[] surfxy;
312
         delete[] surfxy;
312
         delete[] surfz;
313
         delete[] surfz;
313
         return;
314
         return;
314
      }
315
      }
315
 
316
 
316
      // Current z value and the accumulated counter
317
      // Current z value and the accumulated counter
317
      curzval = surfz[0];
318
      curzval = surfz[0];
318
      j = 0;
319
      j = 0;
319
      int acc = 0;
320
      int acc = 0;
320
      int zb;
321
      int zb;
321
      for(int i = 0; i <= (int)nrfiles; i++)
322
      for(int i = 0; i <= (int)nrfiles; i++)
322
      {
323
      {
323
         // Collect the accumulated integral in order to produce a PDF from a CDF
324
         // Collect the accumulated integral in order to produce a PDF from a CDF
324
         // While we are at the same Z value, save under one set
325
         // While we are at the same Z value, save under one set
325
         if( (surfz[i] == curzval) && (acc != nrfiles) )
326
         if( (surfz[i] == curzval) && (acc != nrfiles) )
326
         {
327
         {
327
            integralAcc[j] = integralCount[i];
328
            integralAcc[j] = integralCount[i];
328
            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]);
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]);
329
            j++;
330
            j++;
330
            acc++;
331
            acc++;
331
         }
332
         }
332
         // 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
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
333
         else
334
         else
334
         {
335
         {
335
            // Find minimal and maximal integral values to subtract the offset and normate PDF to 1
336
            // Find minimal and maximal integral values to subtract the offset and normate PDF to 1
336
            NormateSet(i, j, &minInteg, &maxInteg, integralCount, integralAcc);
337
            NormateSet(i, j, &minInteg, &maxInteg, integralCount, integralAcc);
337
 
338
 
338
            if(acc != nrfiles)
339
            if(acc != nrfiles)
339
            {
340
            {
340
               curzval = surfz[i];
341
               curzval = surfz[i];
341
               // PDF and CDF plot
342
               // PDF and CDF plot
342
               PlotEdgeDistribution(files, i, j, &minInteg, &maxInteg, surfxy, integralAcc, direction, edge2d, edit);
343
               PlotEdgeDistribution(files, i, j, &minInteg, &maxInteg, surfxy, integralAcc, direction, edge2d, edit);
343
               i--;
344
               i--;
344
               j = 0;
345
               j = 0;
345
            }
346
            }
346
            else
347
            else
347
            {
348
            {
348
               // PDF and CDF plot
349
               // PDF and CDF plot
349
               PlotEdgeDistribution(files, i, j, &minInteg, &maxInteg, surfxy, integralAcc, direction, edge2d, edit);
350
               PlotEdgeDistribution(files, i, j, &minInteg, &maxInteg, surfxy, integralAcc, direction, edge2d, edit);
350
               i--;
351
               i--;
351
               break;
352
               break;
352
            }
353
            }
353
         }
354
         }
354
 
355
 
355
         // Update the progress bar
356
         // Update the progress bar
356
         progVal = (float)(15.00/nrfiles)*i+75.00;
357
         progVal = (float)(15.00/nrfiles)*i+75.00;
357
         analysisProgress->widgetPB->SetPosition(progVal);
358
         analysisProgress->widgetPB->SetPosition(progVal);
358
         gVirtualX->Update(1);
359
         gVirtualX->Update(1);
359
      }
360
      }
360
 
361
 
361
      // Make the 2D edge plot
362
      // Make the 2D edge plot
362
      if(edge2d)
363
      if(edge2d)
363
      {
364
      {
364
         if(edit == 0)
365
         if(edit == 0)
365
            gCanvas = new TCanvas("canv","canv",900,900);
366
            gCanvas = new TCanvas("canv","canv",900,900);
366
         else
367
         else
367
            gCanvas = tempAnalysisCanvas->GetCanvas();
368
            gCanvas = tempAnalysisCanvas->GetCanvas();
368
 
369
 
369
         double range[4];
370
         double range[4];
370
         TGraph2D *gScan2D;
371
         TGraph2D *gScan2D;
371
         gScan2D = new TGraph2D();
372
         gScan2D = new TGraph2D();
372
         range[0] = TMath::MinElement(nrfiles, surfxy);
373
         range[0] = TMath::MinElement(nrfiles, surfxy);
373
         range[1] = TMath::MaxElement(nrfiles, surfxy);
374
         range[1] = TMath::MaxElement(nrfiles, surfxy);
374
         range[2] = TMath::MinElement(nrfiles, surfz);
375
         range[2] = TMath::MinElement(nrfiles, surfz);
375
         range[3] = TMath::MaxElement(nrfiles, surfz);
376
         range[3] = TMath::MaxElement(nrfiles, surfz);
376
 
377
 
377
         for(int i = 0; i < nrfiles; i++)
378
         for(int i = 0; i < nrfiles; i++)
378
         {
379
         {
379
            if(DBGSIG)
380
            if(DBGSIG)
380
               printf("IntegSpectrum(): %.2lf\t%.2lf\t%lf\n", surfxy[i], surfz[i], integralCount[i]);
381
               printf("IntegSpectrum(): %.2lf\t%.2lf\t%lf\n", surfxy[i], surfz[i], integralCount[i]);
381
            gScan2D->SetPoint(i, surfxy[i], surfz[i], integralCount[i]);
382
            gScan2D->SetPoint(i, surfxy[i], surfz[i], integralCount[i]);
382
 
383
 
383
            // Update the progress bar
384
            // Update the progress bar
384
            progVal = (float)(9.00/nrfiles)*i+90.00;
385
            progVal = (float)(9.00/nrfiles)*i+90.00;
385
            analysisProgress->widgetPB->SetPosition(progVal);
386
            analysisProgress->widgetPB->SetPosition(progVal);
386
            gVirtualX->Update(1);
387
            gVirtualX->Update(1);
387
         }
388
         }
388
 
389
 
389
         gCanvas->cd();
390
         gCanvas->cd();
390
         gStyle->SetPalette(1);
391
         gStyle->SetPalette(1);
391
         gCanvas->SetLeftMargin(0.15);
392
         gCanvas->SetLeftMargin(0.15);
392
         gCanvas->SetRightMargin(0.126);
393
         gCanvas->SetRightMargin(0.126);
393
         gCanvas->SetTopMargin(0.077);
394
         gCanvas->SetTopMargin(0.077);
394
         gScan2D->Draw("COLZ");
395
         gScan2D->Draw("COLZ");
395
         
396
         
396
         gCanvas->Modified();
397
         gCanvas->Modified();
397
         gCanvas->Update();
398
         gCanvas->Update();
398
 
399
 
399
         gScan2D->SetNpx((int)resol2d->widgetNE[0]->GetNumber());
400
         gScan2D->SetNpx((int)resol2d->widgetNE[0]->GetNumber());
400
         gScan2D->SetNpy((int)resol2d->widgetNE[1]->GetNumber());
401
         gScan2D->SetNpy((int)resol2d->widgetNE[1]->GetNumber());
401
         
402
         
402
         gCanvas->Modified();
403
         gCanvas->Modified();
403
         gCanvas->Update();
404
         gCanvas->Update();
404
         
405
         
405
         if(direction == 1)
406
         if(direction == 1)
406
            gScan2D->GetXaxis()->SetTitle("X [#mum]");
407
            gScan2D->GetXaxis()->SetTitle("X [#mum]");
407
         else if(direction == 2)
408
         else if(direction == 2)
408
            gScan2D->GetXaxis()->SetTitle("Y [#mum]");
409
            gScan2D->GetXaxis()->SetTitle("Y [#mum]");
409
 
410
 
410
 
411
 
411
         gScan2D->GetXaxis()->SetTitleOffset(1.3);
412
         gScan2D->GetXaxis()->SetTitleOffset(1.3);
412
         gScan2D->GetXaxis()->CenterTitle(kTRUE);
413
         gScan2D->GetXaxis()->CenterTitle(kTRUE);
413
         gScan2D->GetXaxis()->SetLabelSize(0.027);
414
         gScan2D->GetXaxis()->SetLabelSize(0.027);
414
         gScan2D->GetXaxis()->SetLabelOffset(0.02);
415
         gScan2D->GetXaxis()->SetLabelOffset(0.02);
415
         gScan2D->GetXaxis()->SetRangeUser(range[0], range[1]);
416
         gScan2D->GetXaxis()->SetRangeUser(range[0], range[1]);
416
         gScan2D->GetXaxis()->SetNoExponent(kTRUE);
417
         gScan2D->GetXaxis()->SetNoExponent(kTRUE);
417
         gScan2D->GetYaxis()->SetTitleOffset(1.9);
418
         gScan2D->GetYaxis()->SetTitleOffset(1.9);
418
         gScan2D->GetYaxis()->CenterTitle(kTRUE);
419
         gScan2D->GetYaxis()->CenterTitle(kTRUE);
419
         gScan2D->GetYaxis()->SetLabelSize(0.027);
420
         gScan2D->GetYaxis()->SetLabelSize(0.027);
420
         gScan2D->GetXaxis()->SetLabelOffset(0.02);
421
         gScan2D->GetYaxis()->SetLabelOffset(0.02);
421
         gScan2D->GetYaxis()->SetRangeUser(range[2], range[3]);
422
         gScan2D->GetYaxis()->SetRangeUser(range[2], range[3]);
422
         gScan2D->GetYaxis()->SetNoExponent(kTRUE);
423
         gScan2D->GetYaxis()->SetNoExponent(kTRUE);
423
 
424
 
424
/*         TGaxis *yax = (TGaxis*)gScan2D->GetYaxis();
425
/*         TGaxis *yax = (TGaxis*)gScan2D->GetYaxis();
425
         yax->SetMaxDigits(4);*/
426
         yax->SetMaxDigits(4);*/
426
 
427
 
427
         if(!cleanPlots)
428
         if(!cleanPlots)
428
         {
429
         {
429
            if(direction == 1)
430
            if(direction == 1)
430
            {
431
            {
431
               if(posUnits->widgetCB->GetSelected() == 0)
432
               if(posUnitsPlot->widgetCB->GetSelected() == 0)
432
                  gScan2D->SetTitle("Laser focal point;X [table units];Z [table units]");
433
                  gScan2D->SetTitle("Laser focal point;X [table units];Z [table units]");
433
               else if(posUnits->widgetCB->GetSelected() == 1)
434
               else if(posUnitsPlot->widgetCB->GetSelected() == 1)
434
                  gScan2D->SetTitle("Laser focal point;X [#mum];Z [#mum]");
435
                  gScan2D->SetTitle("Laser focal point;X [#mum];Z [#mum]");
435
            }
436
            }
436
            else if(direction == 2)
437
            else if(direction == 2)
437
            {
438
            {
438
               if(posUnits->widgetCB->GetSelected() == 0)
439
               if(posUnitsPlot->widgetCB->GetSelected() == 0)
439
                  gScan2D->SetTitle("Laser focal point;Y [table units];Z [table units]");
440
                  gScan2D->SetTitle("Laser focal point;Y [table units];Z [table units]");
440
               else if(posUnits->widgetCB->GetSelected() == 1)
441
               else if(posUnitsPlot->widgetCB->GetSelected() == 1)
441
                  gScan2D->SetTitle("Laser focal point;Y [#mum];Z [#mum]");
442
                  gScan2D->SetTitle("Laser focal point;Y [#mum];Z [#mum]");
442
            }
443
            }
443
         }
444
         }
444
         else
445
         else
445
         {
446
         {
446
            if(direction == 1)
447
            if(direction == 1)
447
            {
448
            {
448
               if(posUnits->widgetCB->GetSelected() == 0)
449
               if(posUnitsPlot->widgetCB->GetSelected() == 0)
449
                  gScan2D->SetTitle(";X [table units];Z [table units]");
450
                  gScan2D->SetTitle(";X [table units];Z [table units]");
450
               else if(posUnits->widgetCB->GetSelected() == 1)
451
               else if(posUnitsPlot->widgetCB->GetSelected() == 1)
451
                  gScan2D->SetTitle(";X [#mum];Z [#mum]");
452
                  gScan2D->SetTitle(";X [#mum];Z [#mum]");
452
            }
453
            }
453
            else if(direction == 2)
454
            else if(direction == 2)
454
            {
455
            {
455
               if(posUnits->widgetCB->GetSelected() == 0)
456
               if(posUnitsPlot->widgetCB->GetSelected() == 0)
456
                  gScan2D->SetTitle(";Y [table units];Z [table units]");
457
                  gScan2D->SetTitle(";Y [table units];Z [table units]");
457
               else if(posUnits->widgetCB->GetSelected() == 1)
458
               else if(posUnitsPlot->widgetCB->GetSelected() == 1)
458
                  gScan2D->SetTitle(";Y [#mum];Z [#mum]");
459
                  gScan2D->SetTitle(";Y [#mum];Z [#mum]");
459
            }
460
            }
460
         }
461
         }
461
 
462
 
462
         gCanvas->Modified();
463
         gCanvas->Modified();
463
         gCanvas->Update();
464
         gCanvas->Update();
464
 
465
 
465
         remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp);
466
         remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp);
466
         sprintf(exportname, "%s", ctemp);
467
         sprintf(exportname, "%s", ctemp);
467
         remove_from_last(exportname, '_', ctemp);
468
         remove_from_last(exportname, '_', ctemp);
468
         if(direction == 1)
469
         if(direction == 1)
469
            sprintf(exportname, "%s_xdir_focalpoint.pdf", ctemp);
470
            sprintf(exportname, "%s_xdir_focalpoint.pdf", ctemp);
470
         else if(direction == 2)
471
         else if(direction == 2)
471
            sprintf(exportname, "%s_ydir_focalpoint.pdf", ctemp);
472
            sprintf(exportname, "%s_ydir_focalpoint.pdf", ctemp);
472
         gCanvas->SaveAs(exportname);
473
         gCanvas->SaveAs(exportname);
473
 
474
 
474
         // Update the progress bar
475
         // Update the progress bar
475
         analysisProgress->widgetPB->SetPosition(100.0);
476
         analysisProgress->widgetPB->SetPosition(100.0);
476
         gVirtualX->Update(1);
477
         gVirtualX->Update(1);
477
 
478
 
478
         if(edit == 0)
479
         if(edit == 0)
479
         {
480
         {
480
            delete gScan2D;
481
            delete gScan2D;
481
            delete gCanvas;
482
            delete gCanvas;
482
         }
483
         }
483
      }
484
      }
484
      else
485
      else
485
      {
486
      {
486
         // Update the progress bar
487
         // Update the progress bar
487
         analysisProgress->widgetPB->SetPosition(100.0);
488
         analysisProgress->widgetPB->SetPosition(100.0);
488
         gVirtualX->Update(1);
489
         gVirtualX->Update(1);
489
      }
490
      }
490
   }
491
   }
491
}
492
}
492
 
493
 
493
void TGAppMainFrame::PlotEdgeDistribution(TList *files, int filenr, int points, double *min, double *max, double *xy, double *integAcc, int axis, bool edge2d, int edit)
494
void TGAppMainFrame::PlotEdgeDistribution(TList *files, int filenr, int points, double *min, double *max, double *xy, double *integAcc, int axis, bool edge2d, int edit)
494
{
495
{
495
   TGraph *gScan[2];
496
   TGraph *gScan[2];
496
   int pdfmax = -1;
497
   int pdfmax = -1;
497
   int count = 0;
498
   int count = 0;
498
   char ctemp[1024];
499
   char ctemp[1024];
499
   char exportname[1024];
500
   char exportname[1024];
500
   TCanvas *gCanvas;
501
   TCanvas *gCanvas;
501
 
502
 
502
   // Prepare the CDF plot
503
   // Prepare the CDF plot
503
   gScan[1] = new TGraph();
504
   gScan[1] = new TGraph();
504
   for(int i = 0; i < points; i++)
505
   for(int i = 0; i < points; i++)
505
   {
506
   {
506
      count = filenr - points + i;
507
      count = filenr - points + i;
507
      gScan[1]->SetPoint(i, (double)xy[count], (double)integAcc[i]/(*max));
508
      gScan[1]->SetPoint(i, (double)xy[count], (double)integAcc[i]/(*max));
508
      if(DBGSIG) printf("PlotEdgeDistribution(): CDF %d: %lf, %lf\n", 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));
509
 
510
 
510
      if( ((integAcc[i+1]-integAcc[i])/(*max) > pdfmax) && (i < points-1) )
511
      if( ((integAcc[i+1]-integAcc[i])/(*max) > pdfmax) && (i < points-1) )
511
         pdfmax = (integAcc[i+1]-integAcc[i])/(*max);
512
         pdfmax = (integAcc[i+1]-integAcc[i])/(*max);
512
   }
513
   }
513
 
514
 
514
   pdfmax = (TMath::Ceil(pdfmax*10))/10.;
515
   pdfmax = (TMath::Ceil(pdfmax*10))/10.;
515
 
516
 
516
   // Prepare the PDF plot
517
   // Prepare the PDF plot
517
   gScan[0] = new TGraph();
518
   gScan[0] = new TGraph();
518
   for(int i = points-1; i >= 0; i--)
519
   for(int i = points-1; i >= 0; i--)
519
   {
520
   {
520
      count = (filenr-1) - (points-1) + i;
521
      count = (filenr-1) - (points-1) + i;
521
      // Set any negative values of the PDF to 0
522
      // Set any negative values of the PDF to 0
522
      if( (integAcc[i]-integAcc[i-1])/(*max) < 0 )
523
      if( (integAcc[i]-integAcc[i-1])/(*max) < 0 )
523
         gScan[0]->SetPoint(i, (double)xy[count], 0);
524
         gScan[0]->SetPoint(i, (double)xy[count], 0);
524
      else
525
      else
525
         gScan[0]->SetPoint(i, (double)xy[count], (integAcc[i]-integAcc[i-1])/(*max));
526
         gScan[0]->SetPoint(i, (double)xy[count], (integAcc[i]-integAcc[i-1])/(*max));
526
      if(DBGSIG) printf("PlotEdgeDistribution(): PDF %d: %lf, %lf\n", i, (double)xy[count], (integAcc[i+1]-integAcc[i])/(*max));
527
      if(DBGSIG) printf("PlotEdgeDistribution(): PDF %d: %lf, %lf\n", i, (double)xy[count], (integAcc[i+1]-integAcc[i])/(*max));
527
   }
528
   }
528
 
529
 
529
   remove_from_last((char*)files->At(filenr-1)->GetTitle(), '_', ctemp);
530
   remove_from_last((char*)files->At(filenr-1)->GetTitle(), '_', ctemp);
530
   sprintf(exportname, "%s_edge.pdf", ctemp);
531
   sprintf(exportname, "%s_edge.pdf", ctemp);
531
 
532
 
532
   if(edit == 0)
533
   if(edit == 0)
533
      gCanvas = new TCanvas("canv1d","canv1d",1200,900);
534
      gCanvas = new TCanvas("canv1d","canv1d",1200,900);
534
   else
535
   else
535
      gCanvas = tempAnalysisCanvas->GetCanvas();
536
      gCanvas = tempAnalysisCanvas->GetCanvas();
536
 
537
 
537
   // Fit the PDF with a gaussian
538
   // Fit the PDF with a gaussian
538
   gScan[0]->Fit("gaus","Q");
539
   gScan[0]->Fit("gaus","Q");
539
   gScan[0]->GetFunction("gaus")->SetNpx(400);
540
   gScan[0]->GetFunction("gaus")->SetNpx(400);
540
 
541
 
541
   gStyle->SetOptFit(1);
542
   gStyle->SetOptFit(1);
542
 
543
 
543
   gScan[1]->Draw("AL");
544
   gScan[1]->Draw("AL");
544
   gPad->Update();
545
   gPad->Update();
545
   gScan[0]->Draw("LP");
546
   gScan[0]->Draw("LP");
546
 
547
 
547
   gCanvas->Modified();
548
   gCanvas->Modified();
548
   gCanvas->Update();
549
   gCanvas->Update();
549
 
550
 
550
   TPaveStats *stats = (TPaveStats*)gScan[0]->FindObject("stats");
551
   TPaveStats *stats = (TPaveStats*)gScan[0]->FindObject("stats");
551
   if(!cleanPlots)
552
   if(!cleanPlots)
552
   {
553
   {
553
      stats->SetX1NDC(0.86); stats->SetX2NDC(1.0);
554
      stats->SetX1NDC(0.86); stats->SetX2NDC(1.0);
554
      stats->SetY1NDC(0.87); stats->SetY2NDC(1.0);
555
      stats->SetY1NDC(0.87); stats->SetY2NDC(1.0);
555
   }
556
   }
556
   else
557
   else
557
   {
558
   {
558
      stats->SetX1NDC(1.1); stats->SetX2NDC(1.3);
559
      stats->SetX1NDC(1.1); stats->SetX2NDC(1.3);
559
      stats->SetY1NDC(1.1); stats->SetY2NDC(1.3);
560
      stats->SetY1NDC(1.1); stats->SetY2NDC(1.3);
560
   }
561
   }
561
 
562
 
562
   gCanvas->SetGridx(1);
563
   gCanvas->SetGridx(1);
563
   gCanvas->SetGridy(1);
564
   gCanvas->SetGridy(1);
564
   if(axis == 1)
565
   if(axis == 1)
565
      gScan[1]->GetXaxis()->SetTitle("X [#mum]");
566
      gScan[1]->GetXaxis()->SetTitle("X [#mum]");
566
   else if(axis == 2)
567
   else if(axis == 2)
567
      gScan[1]->GetXaxis()->SetTitle("Y [#mum]");
568
      gScan[1]->GetXaxis()->SetTitle("Y [#mum]");
568
   gScan[1]->GetXaxis()->SetTitleOffset(1.3);
569
   gScan[1]->GetXaxis()->SetTitleOffset(1.3);
569
   gScan[1]->GetXaxis()->CenterTitle(kTRUE);
570
   gScan[1]->GetXaxis()->CenterTitle(kTRUE);
570
   gScan[1]->GetXaxis()->SetLabelSize(0.027);
571
   gScan[1]->GetXaxis()->SetLabelSize(0.027);
571
   gScan[1]->GetXaxis()->SetLabelOffset(0.02);
572
   gScan[1]->GetXaxis()->SetLabelOffset(0.02);
572
   gScan[1]->GetXaxis()->SetNoExponent(kTRUE);
573
   gScan[1]->GetXaxis()->SetNoExponent(kTRUE);
573
   gScan[1]->GetYaxis()->SetTitle("Normalized ADC integral");
574
   gScan[1]->GetYaxis()->SetTitle("Normalized ADC integral");
574
 
575
 
575
   gScan[1]->GetYaxis()->CenterTitle(kTRUE);
576
   gScan[1]->GetYaxis()->CenterTitle(kTRUE);
576
   gScan[1]->GetYaxis()->SetLabelSize(0.027);
577
   gScan[1]->GetYaxis()->SetLabelSize(0.027);
577
   gScan[1]->GetYaxis()->SetLabelOffset(0.02);
578
   gScan[1]->GetYaxis()->SetLabelOffset(0.02);
578
   gScan[1]->GetYaxis()->SetRangeUser(0,1);
579
   gScan[1]->GetYaxis()->SetRangeUser(0,1);
579
   gScan[1]->GetYaxis()->SetTitleOffset(1.4);
580
   gScan[1]->GetYaxis()->SetTitleOffset(1.4);
580
   gScan[1]->GetYaxis()->SetTitleSize(0.030);
581
   gScan[1]->GetYaxis()->SetTitleSize(0.030);
581
 
582
 
582
   if(!cleanPlots)
583
   if(!cleanPlots)
583
   {
584
   {
584
      if(axis == 1)
585
      if(axis == 1)
585
      {
586
      {
586
         if(posUnits->widgetCB->GetSelected() == 0)
587
         if(posUnitsPlot->widgetCB->GetSelected() == 0)
587
            gScan[1]->SetTitle("SiPM edge detection;X [table units];Normalized ADC integral");
588
            gScan[1]->SetTitle("SiPM edge detection;X [table units];Normalized ADC integral");
588
         else if(posUnits->widgetCB->GetSelected() == 1)
589
         else if(posUnitsPlot->widgetCB->GetSelected() == 1)
589
            gScan[1]->SetTitle("SiPM edge detection;X [#mum];Normalized ADC integral");
590
            gScan[1]->SetTitle("SiPM edge detection;X [#mum];Normalized ADC integral");
590
      }
591
      }
591
      else if(axis == 2)
592
      else if(axis == 2)
592
      {
593
      {
593
         if(posUnits->widgetCB->GetSelected() == 0)
594
         if(posUnitsPlot->widgetCB->GetSelected() == 0)
594
            gScan[1]->SetTitle("SiPM edge detection;Y [table units];Normalized ADC integral");
595
            gScan[1]->SetTitle("SiPM edge detection;Y [table units];Normalized ADC integral");
595
         else if(posUnits->widgetCB->GetSelected() == 1)
596
         else if(posUnitsPlot->widgetCB->GetSelected() == 1)
596
            gScan[1]->SetTitle("SiPM edge detection;Y [#mum];Normalized ADC integral");
597
            gScan[1]->SetTitle("SiPM edge detection;Y [#mum];Normalized ADC integral");
597
      }
598
      }
598
   }
599
   }
599
   else
600
   else
600
   {
601
   {
601
      if(axis == 1)
602
      if(axis == 1)
602
      {
603
      {
603
         if(posUnits->widgetCB->GetSelected() == 0)
604
         if(posUnitsPlot->widgetCB->GetSelected() == 0)
604
            gScan[1]->SetTitle(";X [table units];Normalized ADC integral");
605
            gScan[1]->SetTitle(";X [table units];Normalized ADC integral");
605
         else if(posUnits->widgetCB->GetSelected() == 1)
606
         else if(posUnitsPlot->widgetCB->GetSelected() == 1)
606
            gScan[1]->SetTitle(";X [#mum];Normalized ADC integral");
607
            gScan[1]->SetTitle(";X [#mum];Normalized ADC integral");
607
      }
608
      }
608
      else if(axis == 2)
609
      else if(axis == 2)
609
      {
610
      {
610
         if(posUnits->widgetCB->GetSelected() == 0)
611
         if(posUnitsPlot->widgetCB->GetSelected() == 0)
611
            gScan[1]->SetTitle(";Y [table units];Normalized ADC integral");
612
            gScan[1]->SetTitle(";Y [table units];Normalized ADC integral");
612
         else if(posUnits->widgetCB->GetSelected() == 1)
613
         else if(posUnitsPlot->widgetCB->GetSelected() == 1)
613
            gScan[1]->SetTitle(";Y [#mum];Normalized ADC integral");
614
            gScan[1]->SetTitle(";Y [#mum];Normalized ADC integral");
614
      }
615
      }
615
   }
616
   }
616
   gScan[1]->SetLineColor(kBlue);
617
   gScan[1]->SetLineColor(kBlue);
617
   gScan[0]->SetLineWidth(2);
618
   gScan[0]->SetLineWidth(2);
618
   gScan[1]->SetLineWidth(2);
619
   gScan[1]->SetLineWidth(2);
619
 
620
 
620
   gCanvas->Modified();
621
   gCanvas->Modified();
621
   gCanvas->Update();
622
   gCanvas->Update();
622
 
623
 
623
   gCanvas->SaveAs(exportname);
624
   gCanvas->SaveAs(exportname);
624
 
625
 
625
   // If 2D edge plot, delete the 1D edge plots as we go
626
   // If 2D edge plot, delete the 1D edge plots as we go
626
   if(edge2d)
627
   if(edge2d)
627
   {
628
   {
628
      delete gScan[0];
629
      delete gScan[0];
629
      delete gScan[1];
630
      delete gScan[1];
630
      if(edit == 0)
631
      if(edit == 0)
631
         delete gCanvas;
632
         delete gCanvas;
632
   }
633
   }
633
   else
634
   else
634
   {
635
   {
635
      if(edit == 0)
636
      if(edit == 0)
636
      {
637
      {
637
         delete gScan[0];
638
         delete gScan[0];
638
         delete gScan[1];
639
         delete gScan[1];
639
         delete gCanvas;
640
         delete gCanvas;
640
      }
641
      }
641
   }
642
   }
642
}
643
}
643
 
644
 
644
void TGAppMainFrame::PhotonMu(TList *files, int edit)
645
void TGAppMainFrame::PhotonMu(TList *files, int edit)
645
{
646
{
646
   unsigned int nrfiles = files->GetSize();
647
   unsigned int nrfiles = files->GetSize();
647
   char ctemp[1024];
648
   char ctemp[1024];
648
   char exportname[1024];
649
   char exportname[1024];
649
   int j, k = 0, m = 0, n = 0, k2 = 0, m2 = 0;
650
   int j, k = 0, m = 0, n = 0, k2 = 0, m2 = 0;
650
 
651
 
651
   TCanvas *gCanvas;
652
   TCanvas *gCanvas;
652
   TTree *header_data, *meas_data;
653
   TTree *header_data, *meas_data;
653
   double *integralCount, *integralPedestal;
654
   double *integralCount, *integralPedestal;
654
   integralCount = new double[nrfiles];
655
   integralCount = new double[nrfiles];
655
   integralPedestal = new double[nrfiles];
656
   integralPedestal = new double[nrfiles];
656
   double *angle, *pdeval, *muval;
657
   double *angle, *pdeval, *muval;
657
   angle = new double[nrfiles];
658
   angle = new double[nrfiles];
658
   pdeval = new double[nrfiles];
659
   pdeval = new double[nrfiles];
659
   muval = new double[nrfiles];
660
   muval = new double[nrfiles];
660
 
661
 
661
   // Zero the integral count and accumulated vaules
662
   // Zero the integral count and accumulated vaules
662
   for(int i = 0; i < (int)nrfiles; i++) {integralCount[i] = 0; integralPedestal[i] = 0; }
663
   for(int i = 0; i < (int)nrfiles; i++) {integralCount[i] = 0; integralPedestal[i] = 0; }
663
 
664
 
664
   // Fitting variables
665
   // Fitting variables
665
   TSpectrum *spec;
666
   TSpectrum *spec;
666
   TH1F *histtemp;
667
   TH1F *histtemp;
667
   TH1 *histback;
668
   TH1 *histback;
668
   TH1F *h2;
669
   TH1F *h2;
669
   float *xpeaks;
670
   float *xpeaks;
670
   TF1 *fit;
671
   TF1 *fit;
671
   TF1 *fittingfunc;
672
   TF1 *fittingfunc;
672
   double *fparam, *fparamerr;
673
   double *fparam, *fparamerr;
673
   double meansel[20];
674
   double meansel[20];
674
   double sigmasel[20];
675
   double sigmasel[20];
675
   double meanparam, paramsigma;
676
   double meanparam, paramsigma;
676
   int sortindex[20];
677
   int sortindex[20];
677
   int adcpedestal[2];
678
   int adcpedestal[2];
678
   int zeromu = 0;
679
   int zeromu = 0;
679
   int darkhist = -1;
680
   int darkhist = -1;
680
 
681
 
681
   double pointest[12];
682
   double pointest[12];
682
   bool exclude = false;
683
   bool exclude = false;
683
 
684
 
684
   // Zero the parameter values
685
   // Zero the parameter values
685
   for(int i = 0; i < 20; i++) {meansel[i] = 0; sigmasel[i] = 0; }
686
   for(int i = 0; i < 20; i++) {meansel[i] = 0; sigmasel[i] = 0; }
686
 
687
 
687
   float progVal = 0;
688
   float progVal = 0;
688
   analysisProgress->widgetPB->SetPosition(progVal);
689
   analysisProgress->widgetPB->SetPosition(progVal);
689
   gVirtualX->Update(1);
690
   gVirtualX->Update(1);
690
 
691
 
691
   // Start if we select at least one file
692
   // Start if we select at least one file
692
   if(nrfiles > 0)
693
   if(nrfiles > 0)
693
   {
694
   {
694
      for(int i = 0; i < (int)nrfiles; i++)
695
      for(int i = 0; i < (int)nrfiles; i++)
695
      {
696
      {
696
         if(files->At(i))
697
         if( (nrfiles == 1) || (!multiSelect->widgetChBox[0]->IsDown()) )
697
         {
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
         }
-
 
782
         if(files->At(i))
-
 
783
         {
698
            if(strcmp(files->At(i)->GetTitle(),darkRun->widgetTE->GetText()) == 0)
784
            if(strcmp(files->At(i)->GetTitle(),darkRun->widgetTE->GetText()) == 0)
699
            {
785
            {
700
               printf("PhotonMu(): %s is the dark histogram file.\n", files->At(i)->GetTitle());
786
               printf("PhotonMu(): %s is the dark histogram file.\n", files->At(i)->GetTitle());
701
               darkhist = i;
787
               darkhist = i;
702
            }
788
            }
703
 
789
 
704
            // Replot the spectrum on analysisCanvas and do not close the input file
790
            // Replot the spectrum on analysisCanvas and do not close the input file
705
            DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1);
791
            DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1);
706
            analysisCanvas->GetCanvas()->Modified();
792
            analysisCanvas->GetCanvas()->Modified();
707
            analysisCanvas->GetCanvas()->Update();
793
            analysisCanvas->GetCanvas()->Update();
708
       
794
       
709
            // Get the spectrum
795
            // Get the spectrum
710
            histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname);
796
            histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname);
711
            npeaks = 15;
797
            npeaks = 15;
712
            double par[300];
798
            double par[300];
713
            spec = new TSpectrum(npeaks);
799
            spec = new TSpectrum(npeaks);
714
            // Find spectrum background
800
            // Find spectrum background
715
            histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same");
801
            histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same");
716
            // Clone histogram and subtract background from it if we select that option
802
            // Clone histogram and subtract background from it if we select that option
717
            h2 = (TH1F*)histtemp->Clone("h2");
803
            h2 = (TH1F*)histtemp->Clone("h2");
718
            if(fitChecks->widgetChBox[0]->IsDown())
804
            if(fitChecks->widgetChBox[0]->IsDown())
719
               h2->Add(histback, -1);
805
               h2->Add(histback, -1);
720
            // Search for the peaks
806
            // Search for the peaks
721
            int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() );
807
            int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() );
722
            printf("PhotonMu(): Found %d candidates to fit.\n",found);
808
            printf("PhotonMu(): Found %d candidates to fit.\n",found);
723
            npeaks = found;
809
            npeaks = found;
724
   
810
   
725
            // Set initial peak parameters
811
            // Set initial peak parameters
726
            xpeaks = spec->GetPositionX();
812
            xpeaks = spec->GetPositionX();
727
            for(j = 0; j < found; j++)
813
            for(j = 0; j < found; j++)
728
            {
814
            {
729
               float xp = xpeaks[j];
815
               float xp = xpeaks[j];
730
               int bin = h2->GetXaxis()->FindBin(xp);
816
               int bin = h2->GetXaxis()->FindBin(xp);
731
               float yp = h2->GetBinContent(bin);
817
               float yp = h2->GetBinContent(bin);
732
               par[3*j] = yp;
818
               par[3*j] = yp;
733
               par[3*j+1] = xp;
819
               par[3*j+1] = xp;
734
               par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber();
820
               par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber();
735
            }
821
            }
736
         
822
         
737
            // Fit the histogram
823
            // Fit the histogram
738
            fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks);
824
            fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks);
739
            TVirtualFitter::Fitter(histtemp, 3*npeaks);
825
            TVirtualFitter::Fitter(histtemp, 3*npeaks);
740
            fit->SetParameters(par);
826
            fit->SetParameters(par);
741
            fit->SetNpx(300);
827
            fit->SetNpx(300);
742
            h2->Fit("fit","Q");
828
            h2->Fit("fit","Q");
743
            // Get the fitted parameters
829
            // Get the fitted parameters
744
            fittingfunc = h2->GetFunction("fit");
830
            fittingfunc = h2->GetFunction("fit");
745
            fparam = fittingfunc->GetParameters();
831
            fparam = fittingfunc->GetParameters();
746
            fparamerr = fittingfunc->GetParErrors();
832
            fparamerr = fittingfunc->GetParErrors();
747
   
833
   
748
            // Gather the parameters (mean peak value for now)
834
            // Gather the parameters (mean peak value for now)
749
            int j = 1;
835
            int j = 1;
750
            int nrfit = 0;
836
            int nrfit = 0;
751
            while(1)
837
            while(1)
752
            {
838
            {
753
               if( (fparam[j] < 1.E-30) || (nrfit > 8) )
839
               if( (fparam[j] < 1.E-30) || (nrfit > 8) )
754
                  break;
840
                  break;
755
               else
841
               else
756
               {
842
               {
757
                  // Check if pedestal is above the lower limit and sigma is smaller than the mean
843
                  // Check if pedestal is above the lower limit and sigma is smaller than the mean
758
                  if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) )
844
                  if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) )
759
                  {
845
                  {
760
                     // 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
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
761
                     meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber());
847
                     meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber());
762
                     sigmasel[nrfit] = fparam[j+1];
848
                     sigmasel[nrfit] = fparam[j+1];
763
                     nrfit++;
849
                     nrfit++;
764
                  }
850
                  }
765
               }
851
               }
766
   
852
   
767
               j+=3;
853
               j+=3;
768
            }
854
            }
769
            TMath::Sort(nrfit, meansel, sortindex, kFALSE);
855
            TMath::Sort(nrfit, meansel, sortindex, kFALSE);
770
 
856
 
771
            meanparam = meansel[sortindex[0]];
857
            meanparam = meansel[sortindex[0]];
772
            paramsigma = sigmasel[sortindex[0]];
858
            paramsigma = sigmasel[sortindex[0]];
773
 
859
 
774
            for(j = 0; j < nrfit; j++)
860
            for(j = 0; j < nrfit; j++)
775
               if(DBGSIG)
861
               if(DBGSIG)
776
                  printf("PhotonMu(): %d: peak mean = %lf\n", j, meansel[sortindex[j]]);
862
                  printf("PhotonMu(): %d: peak mean = %lf\n", j, meansel[sortindex[j]]);
777
       
863
       
778
            j = 0;
864
            j = 0;
779
            adcpedestal[0] = 0;
865
            adcpedestal[0] = 0;
780
            adcpedestal[1] = -1;
866
            adcpedestal[1] = -1;
781
 
867
 
782
            while(1)
868
            while(1)
783
            {
869
            {
784
               int bin = histtemp->GetXaxis()->FindBin((int)(j+meanparam+paramsigma));
870
               int bin = histtemp->GetXaxis()->FindBin((int)(j+meanparam+paramsigma));
785
               int yp = histtemp->GetBinContent(bin);
871
               int yp = histtemp->GetBinContent(bin);
786
 
872
 
787
               // 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)
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)
788
               if(adcpedestal[1] == -1)
874
               if(adcpedestal[1] == -1)
789
               {
875
               {
790
                  adcpedestal[0] = j+meanparam+paramsigma;
876
                  adcpedestal[0] = j+meanparam+paramsigma;
791
                  adcpedestal[1] = yp;
877
                  adcpedestal[1] = yp;
792
               }
878
               }
793
               else
879
               else
794
               {
880
               {
795
                  if( (npeaks > 1) && (adcpedestal[1] >= yp) )
881
                  if( (npeaks > 1) && (adcpedestal[1] >= yp) )
796
                  {
882
                  {
797
                     adcpedestal[0] = j+meanparam+paramsigma;
883
                     adcpedestal[0] = j+meanparam+paramsigma;
798
                     adcpedestal[1] = yp;
884
                     adcpedestal[1] = yp;
799
                  }
885
                  }
800
                  else if( (npeaks == 1) && (adcpedestal[0] < meanparam+5*paramsigma) ) // TODO -> Determining the pedestal when only one peak
886
                  else if( (npeaks == 1) && (adcpedestal[0] < meanparam+5*paramsigma) ) // TODO -> Determining the pedestal when only one peak
801
                  {
887
                  {
802
                     adcpedestal[0] = j+meanparam+paramsigma;
888
                     adcpedestal[0] = j+meanparam+paramsigma;
803
                     adcpedestal[1] = yp;
889
                     adcpedestal[1] = yp;
804
                  }
890
                  }
805
                  else
891
                  else
806
                     break;
892
                     break;
807
               }
893
               }
808
       
894
       
809
               j++;
895
               j++;
810
               if(j > 50) break;
896
               if(j > 50) break;
811
            }
897
            }
812
 
898
 
813
            if(npeaks > 1)
899
            if(npeaks > 1)
814
            {
900
            {
815
               int bin = histtemp->GetXaxis()->FindBin((int)(meanparam+meansel[sortindex[1]])/2);
901
               int bin = histtemp->GetXaxis()->FindBin((int)(meanparam+meansel[sortindex[1]])/2);
816
               adcpedestal[0] = (meanparam+meansel[sortindex[1]])/2;
902
               adcpedestal[0] = (meanparam+meansel[sortindex[1]])/2;
817
               printf("PhotonMu(): multipeak x = %d, ", adcpedestal[0]);
903
               printf("PhotonMu(): multipeak x = %d, ", adcpedestal[0]);
818
               adcpedestal[1] = histtemp->GetBinContent(bin);
904
               adcpedestal[1] = histtemp->GetBinContent(bin);
819
            }
905
            }
820
 
906
 
821
            if(midPeak->widgetChBox[0]->IsDown())
907
            if(midPeak->widgetChBox[0]->IsDown())
822
            {
908
            {
823
               if( (meanparam - (int)meanparam >= 0.) && (meanparam - (int)meanparam < 0.5) )
909
               if( (meanparam - (int)meanparam >= 0.) && (meanparam - (int)meanparam < 0.5) )
824
                  m = TMath::Floor(meanparam);