Subversion Repositories f9daq

Rev

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

Rev 172 Rev 173
1
#include "../include/sipmscan.h"
1
#include "../include/sipmscan.h"
2
#include "../include/workstation.h"
2
#include "../include/workstation.h"
3
 
3
 
4
#include <stdio.h>
4
#include <stdio.h>
5
#include <stdlib.h>
5
#include <stdlib.h>
6
 
6
 
7
// Peak detection function
7
// Peak detection function
8
int npeaks = 20;
8
int npeaks = 20;
9
double FindPeaks(double *x, double *par)
9
double FindPeaks(double *x, double *par)
10
{
10
{
11
   double result = 0;
11
   double result = 0;
12
   for(int i = 0; i < npeaks; i++)
12
   for(int i = 0; i < npeaks; i++)
13
   {
13
   {
14
      double norm = par[3*i];
14
      double norm = par[3*i];
15
      double mean = par[3*i+1];
15
      double mean = par[3*i+1];
16
      double sigma = par[3*i+2];
16
      double sigma = par[3*i+2];
17
      result += norm*TMath::Gaus(x[0], mean, sigma);
17
      result += norm*TMath::Gaus(x[0], mean, sigma);
18
   }
18
   }
19
   return result;
19
   return result;
20
}
20
}
21
 
21
 
22
// File browser for selecting the dark run histogram
22
// File browser for selecting the dark run histogram
23
void TGAppMainFrame::SelectDarkHist()
23
void TGAppMainFrame::SelectDarkHist()
24
{
24
{
25
   TGFileInfo file_info;
25
   TGFileInfo file_info;
26
   const char *filetypes[] = {"Histograms",histextall,0,0};
26
   const char *filetypes[] = {"Histograms",histextall,0,0};
27
   char *cTemp;
27
   char *cTemp;
28
   file_info.fFileTypes = filetypes;
28
   file_info.fFileTypes = filetypes;
29
//   cTemp = new char[1024];
29
//   cTemp = new char[1024];
30
//   sprintf(cTemp, "%s/results", rootdir);
30
//   sprintf(cTemp, "%s/results", rootdir);
31
//   file_info.fIniDir = StrDup(cTemp);
31
//   file_info.fIniDir = StrDup(cTemp);
32
   file_info.fIniDir = StrDup(currentAnalDir);
32
   file_info.fIniDir = StrDup(currentAnalDir);
33
   file_info.fMultipleSelection = kFALSE;
33
   file_info.fMultipleSelection = kFALSE;
34
   new TGFileDialog(gClient->GetDefaultRoot(), fMain, kFDOpen, &file_info);
34
   new TGFileDialog(gClient->GetDefaultRoot(), fMain, kFDOpen, &file_info);
35
//   delete[] cTemp;
35
//   delete[] cTemp;
36
 
36
 
37
   if(file_info.fFilename != NULL)
37
   if(file_info.fFilename != NULL)
38
   {
38
   {
39
      darkRun->widgetTE->SetText(file_info.fFilename);
39
      darkRun->widgetTE->SetText(file_info.fFilename);
40
      fileList->AddEntry(file_info.fFilename, fileList->GetNumberOfEntries());
40
      fileList->AddEntry(file_info.fFilename, fileList->GetNumberOfEntries());
41
      fileList->Layout();
41
      fileList->Layout();
42
   }
42
   }
43
   else
43
   else
44
      darkRun->widgetTE->SetText("");
44
      darkRun->widgetTE->SetText("");
45
}
45
}
46
 
46
 
47
// Reset analysis defaults for the current analysis type (0 = Integrate spectrum, 1 = Relative PDE, 2 = Breakdown voltage, 3 = Surface scan, 4 = Timing analysis)
47
// Reset analysis defaults for the current analysis type (0 = Integrate spectrum, 1 = Relative PDE, 2 = Breakdown voltage, 3 = Surface scan, 4 = Timing analysis)
48
void TGAppMainFrame::AnalysisDefaults()
48
void TGAppMainFrame::AnalysisDefaults()
49
{
49
{
50
   printf("AnalysisDefaults(): Current analysis tab = %d\n", analTab->GetCurrent());
50
   printf("AnalysisDefaults(): Current analysis tab = %d\n", analTab->GetCurrent());
51
   if(analTab->GetCurrent() == 0)       // Integrate spectrum
51
   if(analTab->GetCurrent() == 0)       // Integrate spectrum
52
   {
52
   {
53
      intSpect->widgetChBox[0]->SetState(kButtonUp);
53
      intSpect->widgetChBox[0]->SetState(kButtonUp);
54
      intSpect->widgetChBox[1]->SetState(kButtonUp);
54
      intSpect->widgetChBox[1]->SetState(kButtonUp);
55
      intSpect->widgetChBox[2]->SetState(kButtonDown);
55
      intSpect->widgetChBox[2]->SetState(kButtonDown);
56
      resol2d->widgetNE[0]->SetNumber(40);
56
      resol2d->widgetNE[0]->SetNumber(40);
57
      resol2d->widgetNE[1]->SetNumber(40);
57
      resol2d->widgetNE[1]->SetNumber(40);
58
   }
58
   }
59
   else if(analTab->GetCurrent() == 1)  // Relative PDE
59
   else if(analTab->GetCurrent() == 1)  // Relative PDE
60
   {
60
   {
61
      relPde->widgetChBox[0]->SetState(kButtonDown);
61
      relPde->widgetChBox[0]->SetState(kButtonDown);
-
 
62
      relPde->widgetChBox[1]->SetState(kButtonUp);
62
      midPeak->widgetChBox[0]->SetState(kButtonUp);
63
      midPeak->widgetChBox[0]->SetState(kButtonUp);
63
      zeroAngle->widgetNE[0]->SetNumber(0.00);
64
      zeroAngle->widgetNE[0]->SetNumber(0.00);
64
   }
65
   }
65
   else if(analTab->GetCurrent() == 2)  // Breakdown voltage
66
   else if(analTab->GetCurrent() == 2)  // Breakdown voltage
66
   {
67
   {
67
      minPeak->widgetNE[0]->SetNumber(2);
68
      minPeak->widgetNE[0]->SetNumber(2);
68
      minPeak->widgetNE[0]->SetNumber(1);
69
      minPeak->widgetNE[0]->SetNumber(1);
69
   }
70
   }
70
   else if(analTab->GetCurrent() == 3)  // Surface scan
71
   else if(analTab->GetCurrent() == 3)  // Surface scan
71
   {
72
   {
72
      surfScanOpt->widgetChBox[0]->SetState(kButtonDown);
73
      surfScanOpt->widgetChBox[0]->SetState(kButtonDown);
73
      surfScanOpt->widgetChBox[1]->SetState(kButtonUp);
74
      surfScanOpt->widgetChBox[1]->SetState(kButtonUp);
74
      resolSurf->widgetNE[0]->SetNumber(40);
75
      resolSurf->widgetNE[0]->SetNumber(40);
75
      resolSurf->widgetNE[1]->SetNumber(40);
76
      resolSurf->widgetNE[1]->SetNumber(40);
76
   }
77
   }
77
}
78
}
78
 
79
 
79
// Analysis handle function
80
// Analysis handle function
80
void TGAppMainFrame::AnalysisHandle(int type)
81
void TGAppMainFrame::AnalysisHandle(int type)
81
{
82
{
82
   TList *files;
83
   TList *files;
83
   bool createTab = true;
84
   bool createTab = true;
84
   int tabid = -1;
85
   int tabid = -1;
85
   int analtab = analTab->GetCurrent();
86
   int analtab = analTab->GetCurrent();
86
 
87
 
87
   int analtype;
88
   int analtype;
88
   if( (analtab == 0) && (!intSpect->widgetChBox[0]->IsDown() && !intSpect->widgetChBox[1]->IsDown()) )
89
   if( (analtab == 0) && (!intSpect->widgetChBox[0]->IsDown() && !intSpect->widgetChBox[1]->IsDown()) )
89
      analtype = 0;
90
      analtype = 0;
90
   else if( (analtab == 0) && (intSpect->widgetChBox[0]->IsDown() || intSpect->widgetChBox[1]->IsDown()) )
91
   else if( (analtab == 0) && (intSpect->widgetChBox[0]->IsDown() || intSpect->widgetChBox[1]->IsDown()) )
91
      analtype = 1;
92
      analtype = 1;
92
   else if( analtab == 1 )
93
   else if( analtab == 1 )
93
      analtype = 2;
94
      analtype = 2;
94
   else if( analtab == 2 )
95
   else if( analtab == 2 )
95
      analtype = 3;
96
      analtype = 3;
96
   else if( analtab == 3 )
97
   else if( analtab == 3 )
97
      analtype = 4;
98
      analtype = 4;
-
 
99
 
-
 
100
   // Save analysis settings any time we run a new analysis
-
 
101
   SaveAnalSettings();
98
 
102
 
99
   // Only integrate spectrum or make relative PDE
103
   // Only integrate spectrum or make relative PDE
100
   if(type == 0)
104
   if(type == 0)
101
   {
105
   {
102
      files = new TList();
106
      files = new TList();
103
      fileList->GetSelectedEntries(files);
107
      fileList->GetSelectedEntries(files);
104
 
108
 
105
      if( analtype == 0 )
109
      if( analtype == 0 )
106
         IntegSpectrum(files, 0, 0);
110
         IntegSpectrum(files, 0, 0);
107
 
111
 
108
      if( intSpect->widgetChBox[0]->IsDown() )
112
      if( intSpect->widgetChBox[0]->IsDown() )
109
         IntegSpectrum(files, 1, 0);
113
         IntegSpectrum(files, 1, 0);
110
 
114
 
111
      if( intSpect->widgetChBox[1]->IsDown() )
115
      if( intSpect->widgetChBox[1]->IsDown() )
112
         IntegSpectrum(files, 2, 0);
116
         IntegSpectrum(files, 2, 0);
113
 
117
 
114
      if( analtype == 2 )
118
      if( analtype == 2 )
115
         PhotonMu(files, 0);
119
         PhotonMu(files, 0);
116
 
120
 
117
      if( analtype == 3 )
121
      if( analtype == 3 )
118
         BreakdownVolt(files, 0);
122
         BreakdownVolt(files, 0);
119
 
123
 
120
      if( analtype == 4 )
124
      if( analtype == 4 )
121
         SurfaceScan(files, 0);
125
         SurfaceScan(files, 0);
122
   }
126
   }
123
   // Integrate spectrum or make relative PDE and open edit window
127
   // Integrate spectrum or make relative PDE and open edit window
124
   else if(type == 1)
128
   else if(type == 1)
125
   {
129
   {
126
      files = new TList();
130
      files = new TList();
127
      fileList->GetSelectedEntries(files);
131
      fileList->GetSelectedEntries(files);
128
 
132
 
129
      // Prepare a new analysis edit tab, if we want to edit plots
133
      // Prepare a new analysis edit tab, if we want to edit plots
130
      for(int i = 0; i < fTab->GetNumberOfTabs(); i++)
134
      for(int i = 0; i < fTab->GetNumberOfTabs(); i++)
131
      {
135
      {
132
         if(strcmp("Analysis edit", fTab->GetTabTab(i)->GetString() ) == 0)
136
         if(strcmp("Analysis edit", fTab->GetTabTab(i)->GetString() ) == 0)
133
         {
137
         {
134
            createTab = false;
138
            createTab = false;
135
            tabid = i;
139
            tabid = i;
136
         }
140
         }
137
         
141
         
138
         if(DBGSIG) printf("AnalysisHandle(): Name of tab = %s\n", fTab->GetTabTab(i)->GetString() );
142
         if(DBGSIG) printf("AnalysisHandle(): Name of tab = %s\n", fTab->GetTabTab(i)->GetString() );
139
      }
143
      }
140
   
144
   
141
      if(files->GetSize() > 0)
145
      if(files->GetSize() > 0)
142
      {
146
      {
143
         TempAnalysisTab(fTab, createTab, &tabid, analtype);
147
         TempAnalysisTab(fTab, createTab, &tabid, analtype);
144
 
148
 
145
         // Integrate spectra
149
         // Integrate spectra
146
         if( analtype == 0 )
150
         if( analtype == 0 )
147
            IntegSpectrum(files, 0, 1);
151
            IntegSpectrum(files, 0, 1);
148
 
152
 
149
         if( intSpect->widgetChBox[0]->IsDown() )
153
         if( intSpect->widgetChBox[0]->IsDown() )
150
            IntegSpectrum(files, 1, 1);
154
            IntegSpectrum(files, 1, 1);
151
 
155
 
152
         if( intSpect->widgetChBox[1]->IsDown() )
156
         if( intSpect->widgetChBox[1]->IsDown() )
153
            IntegSpectrum(files, 2, 1);
157
            IntegSpectrum(files, 2, 1);
154
 
158
 
155
         if( analtype == 2 )
159
         if( analtype == 2 )
156
            PhotonMu(files, 1);
160
            PhotonMu(files, 1);
157
 
161
 
158
         if( analtype == 3 )
162
         if( analtype == 3 )
159
            BreakdownVolt(files, 1);
163
            BreakdownVolt(files, 1);
160
 
164
 
161
         if( analtype == 4 )
165
         if( analtype == 4 )
162
            SurfaceScan(files, 1);
166
            SurfaceScan(files, 1);
163
 
167
 
164
         fTab->SetTab(tabid);
168
         fTab->SetTab(tabid);
165
      }
169
      }
166
   }
170
   }
167
 
171
 
168
   delete files;
172
   delete files;
169
}
173
}
170
 
174
 
171
// Analysis functions ----------------------------------
175
// Analysis functions ----------------------------------
172
void TGAppMainFrame::IntegSpectrum(TList *files, int direction, int edit)
176
void TGAppMainFrame::IntegSpectrum(TList *files, int direction, int edit)
173
{
177
{
174
   unsigned int nrfiles = files->GetSize();
178
   unsigned int nrfiles = files->GetSize();
175
   char ctemp[1024];
179
   char ctemp[1024];
176
   char exportname[1024];
180
   char exportname[1024];
177
   int j, k = 0, m = 0;
181
   int j, k = 0, m = 0;
178
 
182
 
179
   TTree *header_data, *meas_data;
183
   TTree *header_data, *meas_data;
180
   double *integralCount, *integralAcc;
184
   double *integralCount, *integralAcc;
181
   integralCount = new double[nrfiles];
185
   integralCount = new double[nrfiles];
182
   integralAcc = new double[nrfiles];
186
   integralAcc = new double[nrfiles];
183
   double *surfxy, *surfz;
187
   double *surfxy, *surfz;
184
   surfxy = new double[nrfiles];
188
   surfxy = new double[nrfiles];
185
   surfz = new double[nrfiles];
189
   surfz = new double[nrfiles];
186
   double minInteg, maxInteg;
190
   double minInteg, maxInteg;
187
   bool norm = intSpect->widgetChBox[2]->IsDown();
191
   bool norm = intSpect->widgetChBox[2]->IsDown();
188
   double curzval;
192
   double curzval;
189
   bool edge2d = false;
193
   bool edge2d = false;
190
   TCanvas *gCanvas;
194
   TCanvas *gCanvas;
191
 
195
 
192
   float progVal = 0;
196
   float progVal = 0;
193
   analysisProgress->widgetPB->SetPosition(progVal);
197
   analysisProgress->widgetPB->SetPosition(progVal);
194
   gVirtualX->Update(1);
198
   gVirtualX->Update(1);
195
 
199
 
196
   // Zero the integral count and accumulated vaules
200
   // Zero the integral count and accumulated vaules
197
   for(int i = 0; i < (int)nrfiles; i++) {integralCount[i] = 0; integralAcc[i] = 0; }
201
   for(int i = 0; i < (int)nrfiles; i++) {integralCount[i] = 0; integralAcc[i] = 0; }
198
 
202
 
199
   // Start if we select at least one file
203
   // Start if we select at least one file
200
   if(nrfiles > 0)
204
   if(nrfiles > 0)
201
   {
205
   {
202
      for(int i = 0; i < (int)nrfiles; i++)
206
      for(int i = 0; i < (int)nrfiles; i++)
203
      {
207
      {
204
         if(files->At(i))
208
         if(files->At(i))
205
         {
209
         {
206
            // Read the X,Y and Z positions from header and ADC and TDC values from the measurements
210
            // Read the X,Y and Z positions from header and ADC and TDC values from the measurements
207
            sprintf(ctemp, "%s", files->At(i)->GetTitle());
211
            sprintf(ctemp, "%s", files->At(i)->GetTitle());
208
            inroot = new TFile(ctemp, "READ");
212
            inroot = new TFile(ctemp, "READ");
209
         
213
         
210
            inroot->GetObject("header_data", header_data);
214
            inroot->GetObject("header_data", header_data);
211
            inroot->GetObject("meas_data", meas_data);
215
            inroot->GetObject("meas_data", meas_data);
212
         
216
         
213
            header_data->SetBranchAddress("xpos", &evtheader.xpos);
217
            header_data->SetBranchAddress("xpos", &evtheader.xpos);
214
            header_data->GetEntry(0);
218
            header_data->GetEntry(0);
215
            header_data->SetBranchAddress("ypos", &evtheader.ypos);
219
            header_data->SetBranchAddress("ypos", &evtheader.ypos);
216
            header_data->GetEntry(0);
220
            header_data->GetEntry(0);
217
            header_data->SetBranchAddress("zpos", &evtheader.zpos);
221
            header_data->SetBranchAddress("zpos", &evtheader.zpos);
218
            header_data->GetEntry(0);
222
            header_data->GetEntry(0);
219
 
223
 
220
            char rdc[256];
224
            char rdc[256];
221
            j = selectCh->widgetNE[0]->GetNumber();
225
            j = selectCh->widgetNE[0]->GetNumber();
222
            double rangetdc[2];
226
            double rangetdc[2];
223
            rangetdc[0] = tdcRange->widgetNE[0]->GetNumber();
227
            rangetdc[0] = tdcRange->widgetNE[0]->GetNumber();
224
            rangetdc[1] = tdcRange->widgetNE[1]->GetNumber();
228
            rangetdc[1] = tdcRange->widgetNE[1]->GetNumber();
225
   
229
   
226
            k = 0;
230
            k = 0;
227
            m = 0;
231
            m = 0;
228
         
232
         
229
            // Reading the data
233
            // Reading the data
230
            for(int e = 0; e < meas_data->GetEntries(); e++)
234
            for(int e = 0; e < meas_data->GetEntries(); e++)
231
            {
235
            {
232
               sprintf(rdc, "ADC%d", j);
236
               sprintf(rdc, "ADC%d", j);
233
               meas_data->SetBranchAddress(rdc, &evtdata.adcdata[j]);
237
               meas_data->SetBranchAddress(rdc, &evtdata.adcdata[j]);
234
               meas_data->GetEntry(e);
238
               meas_data->GetEntry(e);
235
         
239
         
236
               sprintf(rdc, "TDC%d", j);
240
               sprintf(rdc, "TDC%d", j);
237
               meas_data->SetBranchAddress(rdc, &evtdata.tdcdata[j]);
241
               meas_data->SetBranchAddress(rdc, &evtdata.tdcdata[j]);
238
               meas_data->GetEntry(e);
242
               meas_data->GetEntry(e);
239
   
243
   
240
               // Use data point only if it is inside the TDC window
244
               // Use data point only if it is inside the TDC window
241
               if( ((double)evtdata.tdcdata[j]/tdctimeconversion >= rangetdc[0]) && ((double)evtdata.tdcdata[j]/tdctimeconversion <= rangetdc[1]) )
245
               if( ((double)evtdata.tdcdata[j]/tdctimeconversion >= rangetdc[0]) && ((double)evtdata.tdcdata[j]/tdctimeconversion <= rangetdc[1]) )
242
               {
246
               {
243
                  k++;
247
                  k++;
244
                  m += evtdata.adcdata[j];
248
                  m += evtdata.adcdata[j];
245
               }
249
               }
246
            }
250
            }
247
 
251
 
248
            // X, Y and Z values from each file (table units or microns)
252
            // X, Y and Z values from each file (table units or microns)
249
            if(posUnitsPlot->widgetCB->GetSelected() == 0)
253
            if(posUnitsPlot->widgetCB->GetSelected() == 0)
250
            {
254
            {
251
               if(direction == 1)
255
               if(direction == 1)
252
                  surfxy[i] = (double)(evtheader.xpos);
256
                  surfxy[i] = (double)(evtheader.xpos);
253
               else if(direction == 2)
257
               else if(direction == 2)
254
                  surfxy[i] = (double)(evtheader.ypos);
258
                  surfxy[i] = (double)(evtheader.ypos);
255
               surfz[i] = (double)(evtheader.zpos);
259
               surfz[i] = (double)(evtheader.zpos);
256
            }
260
            }
257
            else if(posUnitsPlot->widgetCB->GetSelected() == 1)
261
            else if(posUnitsPlot->widgetCB->GetSelected() == 1)
258
            {
262
            {
259
               if(direction == 1)
263
               if(direction == 1)
260
                  surfxy[i] = (double)(evtheader.xpos*lenconversion);
264
                  surfxy[i] = (double)(evtheader.xpos*lenconversion);
261
               else if(direction == 2)
265
               else if(direction == 2)
262
                  surfxy[i] = (double)(evtheader.ypos*lenconversion);
266
                  surfxy[i] = (double)(evtheader.ypos*lenconversion);
263
               surfz[i] = (double)(evtheader.zpos*lenconversion);
267
               surfz[i] = (double)(evtheader.zpos*lenconversion);
264
            }
268
            }
265
 
269
 
266
            // Check if we have different Z values: if no, just make the edge plots; if yes, save edge plots and make a 2d edge plot
270
            // Check if we have different Z values: if no, just make the edge plots; if yes, save edge plots and make a 2d edge plot
267
            if(i == 0) curzval = surfz[i];
271
            if(i == 0) curzval = surfz[i];
268
            else
272
            else
269
            {
273
            {
270
               if(surfz[i] != curzval)
274
               if(surfz[i] != curzval)
271
                  edge2d = true;
275
                  edge2d = true;
272
            }
276
            }
273
 
277
 
274
            // Print the calculated integral, if no X or Y direction edge plots are enabled; otherwise, just save for later plotting
278
            // Print the calculated integral, if no X or Y direction edge plots are enabled; otherwise, just save for later plotting
275
            if(direction == 0)
279
            if(direction == 0)
276
            {
280
            {
277
               if(norm)
281
               if(norm)
278
               {
282
               {
279
                  integralCount[i] += ((double)m)/((double)k);
283
                  integralCount[i] += ((double)m)/((double)k);
280
                  printf("IntegSpectrum(): %s: The integral is: %lf\n", ctemp, integralCount[i]);
284
                  printf("IntegSpectrum(): %s: The integral is: %lf\n", ctemp, integralCount[i]);
281
               }
285
               }
282
               else
286
               else
283
               {
287
               {
284
                  integralCount[i] += m;
288
                  integralCount[i] += m;
285
                  printf("IntegSpectrum(): %s: The integral is: %d\n", ctemp, (int)integralCount[i]);
289
                  printf("IntegSpectrum(): %s: The integral is: %d\n", ctemp, (int)integralCount[i]);
286
               }
290
               }
287
            }
291
            }
288
            else
292
            else
289
            {
293
            {
290
               if(norm)
294
               if(norm)
291
                  integralCount[i] += ((double)m)/((double)k);
295
                  integralCount[i] += ((double)m)/((double)k);
292
               else
296
               else
293
                  integralCount[i] += m;
297
                  integralCount[i] += m;
294
            }
298
            }
295
           
299
           
296
            inroot->Close();
300
            inroot->Close();
297
            delete inroot;
301
            delete inroot;
298
         }
302
         }
299
 
303
 
300
         // Update the progress bar
304
         // Update the progress bar
301
         progVal = (float)(75.00/nrfiles)*i;
305
         progVal = (float)(75.00/nrfiles)*i;
302
         analysisProgress->widgetPB->SetPosition(progVal);
306
         analysisProgress->widgetPB->SetPosition(progVal);
303
         gVirtualX->Update(1);
307
         gVirtualX->Update(1);
304
      }
308
      }
305
 
309
 
306
      printf("IntegSpectrum(): %d files were selected.\n", nrfiles);
310
      printf("IntegSpectrum(): %d files were selected.\n", nrfiles);
307
 
311
 
308
      // If only an integral is needed, do not plot and exit here
312
      // If only an integral is needed, do not plot and exit here
309
      if( direction == 0 )
313
      if( direction == 0 )
310
      {
314
      {
311
         delete[] integralCount;
315
         delete[] integralCount;
312
         delete[] surfxy;
316
         delete[] surfxy;
313
         delete[] surfz;
317
         delete[] surfz;
314
         return;
318
         return;
315
      }
319
      }
316
 
320
 
317
      // Current z value and the accumulated counter
321
      // Current z value and the accumulated counter
318
      curzval = surfz[0];
322
      curzval = surfz[0];
319
      j = 0;
323
      j = 0;
320
      int acc = 0;
324
      int acc = 0;
321
      int zb;
325
      int zb;
322
      for(int i = 0; i <= (int)nrfiles; i++)
326
      for(int i = 0; i <= (int)nrfiles; i++)
323
      {
327
      {
324
         // Collect the accumulated integral in order to produce a PDF from a CDF
328
         // Collect the accumulated integral in order to produce a PDF from a CDF
325
         // While we are at the same Z value, save under one set
329
         // While we are at the same Z value, save under one set
326
         if( (surfz[i] == curzval) && (acc != nrfiles) )
330
         if( (surfz[i] == curzval) && (acc != nrfiles) )
327
         {
331
         {
328
            integralAcc[j] = integralCount[i];
332
            integralAcc[j] = integralCount[i];
329
            if(DBGSIG) printf("IntegSpectrum(): Integral check 1 (i=%d,j=%d,z=%.2lf): %lf\t%lf\n", i, j, surfz[i], integralCount[i], integralAcc[j]);
333
            if(DBGSIG) printf("IntegSpectrum(): Integral check 1 (i=%d,j=%d,z=%.2lf): %lf\t%lf\n", i, j, surfz[i], integralCount[i], integralAcc[j]);
330
            j++;
334
            j++;
331
            acc++;
335
            acc++;
332
         }
336
         }
333
         // When we switch to a new set of Z values and at the end, we must save the previous ones to make 1D edge plots
337
         // When we switch to a new set of Z values and at the end, we must save the previous ones to make 1D edge plots
334
         else
338
         else
335
         {
339
         {
336
            // Find minimal and maximal integral values to subtract the offset and normate PDF to 1
340
            // Find minimal and maximal integral values to subtract the offset and normate PDF to 1
337
            NormateSet(i, j, &minInteg, &maxInteg, integralCount, integralAcc);
341
            NormateSet(i, j, &minInteg, &maxInteg, integralCount, integralAcc);
338
 
342
 
339
            if(acc != nrfiles)
343
            if(acc != nrfiles)
340
            {
344
            {
341
               curzval = surfz[i];
345
               curzval = surfz[i];
342
               // PDF and CDF plot
346
               // PDF and CDF plot
343
               PlotEdgeDistribution(files, i, j, &minInteg, &maxInteg, surfxy, integralAcc, direction, edge2d, edit);
347
               PlotEdgeDistribution(files, i, j, &minInteg, &maxInteg, surfxy, integralAcc, direction, edge2d, edit);
344
               i--;
348
               i--;
345
               j = 0;
349
               j = 0;
346
            }
350
            }
347
            else
351
            else
348
            {
352
            {
349
               // PDF and CDF plot
353
               // PDF and CDF plot
350
               PlotEdgeDistribution(files, i, j, &minInteg, &maxInteg, surfxy, integralAcc, direction, edge2d, edit);
354
               PlotEdgeDistribution(files, i, j, &minInteg, &maxInteg, surfxy, integralAcc, direction, edge2d, edit);
351
               i--;
355
               i--;
352
               break;
356
               break;
353
            }
357
            }
354
         }
358
         }
355
 
359
 
356
         // Update the progress bar
360
         // Update the progress bar
357
         progVal = (float)(15.00/nrfiles)*i+75.00;
361
         progVal = (float)(15.00/nrfiles)*i+75.00;
358
         analysisProgress->widgetPB->SetPosition(progVal);
362
         analysisProgress->widgetPB->SetPosition(progVal);
359
         gVirtualX->Update(1);
363
         gVirtualX->Update(1);
360
      }
364
      }
361
 
365
 
362
      // Make the 2D edge plot
366
      // Make the 2D edge plot
363
      if(edge2d)
367
      if(edge2d)
364
      {
368
      {
365
         if(edit == 0)
369
         if(edit == 0)
366
            gCanvas = new TCanvas("canv","canv",900,900);
370
            gCanvas = new TCanvas("canv","canv",900,900);
367
         else
371
         else
368
            gCanvas = tempAnalysisCanvas->GetCanvas();
372
            gCanvas = tempAnalysisCanvas->GetCanvas();
369
 
373
 
370
         double range[4];
374
         double range[4];
371
         TGraph2D *gScan2D;
375
         TGraph2D *gScan2D;
372
         gScan2D = new TGraph2D();
376
         gScan2D = new TGraph2D();
373
         gScan2D->SetName("edgescan");
-
 
374
         range[0] = TMath::MinElement(nrfiles, surfxy);
377
         range[0] = TMath::MinElement(nrfiles, surfxy);
375
         range[1] = TMath::MaxElement(nrfiles, surfxy);
378
         range[1] = TMath::MaxElement(nrfiles, surfxy);
376
         range[2] = TMath::MinElement(nrfiles, surfz);
379
         range[2] = TMath::MinElement(nrfiles, surfz);
377
         range[3] = TMath::MaxElement(nrfiles, surfz);
380
         range[3] = TMath::MaxElement(nrfiles, surfz);
378
 
381
 
379
         for(int i = 0; i < nrfiles; i++)
382
         for(int i = 0; i < nrfiles; i++)
380
         {
383
         {
381
            if(DBGSIG)
384
            if(DBGSIG)
382
               printf("IntegSpectrum(): %.2lf\t%.2lf\t%lf\n", surfxy[i], surfz[i], integralCount[i]);
385
               printf("IntegSpectrum(): %.2lf\t%.2lf\t%lf\n", surfxy[i], surfz[i], integralCount[i]);
383
            gScan2D->SetPoint(i, surfxy[i], surfz[i], integralCount[i]);
386
            gScan2D->SetPoint(i, surfxy[i], surfz[i], integralCount[i]);
384
 
387
 
385
            // Update the progress bar
388
            // Update the progress bar
386
            progVal = (float)(9.00/nrfiles)*i+90.00;
389
            progVal = (float)(9.00/nrfiles)*i+90.00;
387
            analysisProgress->widgetPB->SetPosition(progVal);
390
            analysisProgress->widgetPB->SetPosition(progVal);
388
            gVirtualX->Update(1);
391
            gVirtualX->Update(1);
389
         }
392
         }
390
 
393
 
391
         gCanvas->cd();
394
         gCanvas->cd();
392
         gStyle->SetPalette(1);
395
         gStyle->SetPalette(1);
393
         gCanvas->SetLeftMargin(0.15);
396
         gCanvas->SetLeftMargin(0.15);
394
         gCanvas->SetRightMargin(0.126);
397
         gCanvas->SetRightMargin(0.126);
395
         gCanvas->SetTopMargin(0.077);
398
         gCanvas->SetTopMargin(0.077);
396
         gScan2D->Draw("COLZ");
399
         gScan2D->Draw("COLZ");
397
         
400
         
398
         gCanvas->Modified();
401
         gCanvas->Modified();
399
         gCanvas->Update();
402
         gCanvas->Update();
400
 
403
 
401
         gScan2D->SetNpx((int)resol2d->widgetNE[0]->GetNumber());
404
         gScan2D->SetNpx((int)resol2d->widgetNE[0]->GetNumber());
402
         gScan2D->SetNpy((int)resol2d->widgetNE[1]->GetNumber());
405
         gScan2D->SetNpy((int)resol2d->widgetNE[1]->GetNumber());
403
         
406
         
404
         gCanvas->Modified();
407
         gCanvas->Modified();
405
         gCanvas->Update();
408
         gCanvas->Update();
406
         
409
         
407
         if(direction == 1)
410
         if(direction == 1)
408
            gScan2D->GetXaxis()->SetTitle("X [#mum]");
411
            gScan2D->GetXaxis()->SetTitle("X [#mum]");
409
         else if(direction == 2)
412
         else if(direction == 2)
410
            gScan2D->GetXaxis()->SetTitle("Y [#mum]");
413
            gScan2D->GetXaxis()->SetTitle("Y [#mum]");
411
 
414
 
412
 
415
 
413
         gScan2D->GetXaxis()->SetTitleOffset(1.3);
416
         gScan2D->GetXaxis()->SetTitleOffset(1.3);
414
         gScan2D->GetXaxis()->CenterTitle(kTRUE);
417
         gScan2D->GetXaxis()->CenterTitle(kTRUE);
415
         gScan2D->GetXaxis()->SetLabelSize(0.027);
418
         gScan2D->GetXaxis()->SetLabelSize(0.027);
416
         gScan2D->GetXaxis()->SetLabelOffset(0.02);
419
         gScan2D->GetXaxis()->SetLabelOffset(0.02);
417
         gScan2D->GetXaxis()->SetRangeUser(range[0], range[1]);
420
         gScan2D->GetXaxis()->SetRangeUser(range[0], range[1]);
418
         gScan2D->GetXaxis()->SetNoExponent(kTRUE);
421
         gScan2D->GetXaxis()->SetNoExponent(kTRUE);
419
         gScan2D->GetYaxis()->SetTitleOffset(1.9);
422
         gScan2D->GetYaxis()->SetTitleOffset(1.9);
420
         gScan2D->GetYaxis()->CenterTitle(kTRUE);
423
         gScan2D->GetYaxis()->CenterTitle(kTRUE);
421
         gScan2D->GetYaxis()->SetLabelSize(0.027);
424
         gScan2D->GetYaxis()->SetLabelSize(0.027);
422
         gScan2D->GetYaxis()->SetLabelOffset(0.02);
425
         gScan2D->GetYaxis()->SetLabelOffset(0.02);
423
         gScan2D->GetYaxis()->SetRangeUser(range[2], range[3]);
426
         gScan2D->GetYaxis()->SetRangeUser(range[2], range[3]);
424
         gScan2D->GetYaxis()->SetNoExponent(kTRUE);
427
         gScan2D->GetYaxis()->SetNoExponent(kTRUE);
425
 
428
 
426
/*         TGaxis *yax = (TGaxis*)gScan2D->GetYaxis();
429
/*         TGaxis *yax = (TGaxis*)gScan2D->GetYaxis();
427
         yax->SetMaxDigits(4);*/
430
         yax->SetMaxDigits(4);*/
428
 
431
 
429
         if(!cleanPlots)
432
         if(!cleanPlots)
430
         {
433
         {
431
            if(direction == 1)
434
            if(direction == 1)
432
            {
435
            {
433
               if(posUnitsPlot->widgetCB->GetSelected() == 0)
436
               if(posUnitsPlot->widgetCB->GetSelected() == 0)
434
                  gScan2D->SetTitle("Laser focal point;X [table units];Z [table units]");
437
                  gScan2D->SetTitle("Laser focal point;X [table units];Z [table units]");
435
               else if(posUnitsPlot->widgetCB->GetSelected() == 1)
438
               else if(posUnitsPlot->widgetCB->GetSelected() == 1)
436
                  gScan2D->SetTitle("Laser focal point;X [#mum];Z [#mum]");
439
                  gScan2D->SetTitle("Laser focal point;X [#mum];Z [#mum]");
437
            }
440
            }
438
            else if(direction == 2)
441
            else if(direction == 2)
439
            {
442
            {
440
               if(posUnitsPlot->widgetCB->GetSelected() == 0)
443
               if(posUnitsPlot->widgetCB->GetSelected() == 0)
441
                  gScan2D->SetTitle("Laser focal point;Y [table units];Z [table units]");
444
                  gScan2D->SetTitle("Laser focal point;Y [table units];Z [table units]");
442
               else if(posUnitsPlot->widgetCB->GetSelected() == 1)
445
               else if(posUnitsPlot->widgetCB->GetSelected() == 1)
443
                  gScan2D->SetTitle("Laser focal point;Y [#mum];Z [#mum]");
446
                  gScan2D->SetTitle("Laser focal point;Y [#mum];Z [#mum]");
444
            }
447
            }
445
         }
448
         }
446
         else
449
         else
447
         {
450
         {
448
            if(direction == 1)
451
            if(direction == 1)
449
            {
452
            {
450
               if(posUnitsPlot->widgetCB->GetSelected() == 0)
453
               if(posUnitsPlot->widgetCB->GetSelected() == 0)
451
                  gScan2D->SetTitle(";X [table units];Z [table units]");
454
                  gScan2D->SetTitle(";X [table units];Z [table units]");
452
               else if(posUnitsPlot->widgetCB->GetSelected() == 1)
455
               else if(posUnitsPlot->widgetCB->GetSelected() == 1)
453
                  gScan2D->SetTitle(";X [#mum];Z [#mum]");
456
                  gScan2D->SetTitle(";X [#mum];Z [#mum]");
454
            }
457
            }
455
            else if(direction == 2)
458
            else if(direction == 2)
456
            {
459
            {
457
               if(posUnitsPlot->widgetCB->GetSelected() == 0)
460
               if(posUnitsPlot->widgetCB->GetSelected() == 0)
458
                  gScan2D->SetTitle(";Y [table units];Z [table units]");
461
                  gScan2D->SetTitle(";Y [table units];Z [table units]");
459
               else if(posUnitsPlot->widgetCB->GetSelected() == 1)
462
               else if(posUnitsPlot->widgetCB->GetSelected() == 1)
460
                  gScan2D->SetTitle(";Y [#mum];Z [#mum]");
463
                  gScan2D->SetTitle(";Y [#mum];Z [#mum]");
461
            }
464
            }
462
         }
465
         }
463
 
466
 
464
         gCanvas->Modified();
467
         gCanvas->Modified();
465
         gCanvas->Update();
468
         gCanvas->Update();
466
 
469
 
467
         remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp);
470
         remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp);
468
         sprintf(exportname, "%s", ctemp);
471
         sprintf(exportname, "%s", ctemp);
469
         remove_from_last(exportname, '_', ctemp);
472
         remove_from_last(exportname, '_', ctemp);
470
         if(direction == 1)
473
         if(direction == 1)
471
            sprintf(exportname, "%s_xdir_focalpoint.pdf", ctemp);
474
            sprintf(exportname, "%s_xdir_focalpoint.pdf", ctemp);
472
         else if(direction == 2)
475
         else if(direction == 2)
473
            sprintf(exportname, "%s_ydir_focalpoint.pdf", ctemp);
476
            sprintf(exportname, "%s_ydir_focalpoint.pdf", ctemp);
474
         gCanvas->SaveAs(exportname);
477
         gCanvas->SaveAs(exportname);
475
 
478
 
476
         // Update the progress bar
479
         // Update the progress bar
477
         analysisProgress->widgetPB->SetPosition(100.0);
480
         analysisProgress->widgetPB->SetPosition(100.0);
478
         gVirtualX->Update(1);
481
         gVirtualX->Update(1);
479
 
482
 
480
         if(edit == 0)
483
         if(edit == 0)
481
         {
484
         {
482
            delete gScan2D;
485
            delete gScan2D;
483
            delete gCanvas;
486
            delete gCanvas;
484
         }
487
         }
485
      }
488
      }
486
      else
489
      else
487
      {
490
      {
488
         // Update the progress bar
491
         // Update the progress bar
489
         analysisProgress->widgetPB->SetPosition(100.0);
492
         analysisProgress->widgetPB->SetPosition(100.0);
490
         gVirtualX->Update(1);
493
         gVirtualX->Update(1);
491
      }
494
      }
492
   }
495
   }
493
}
496
}
494
 
497
 
495
void TGAppMainFrame::PlotEdgeDistribution(TList *files, int filenr, int points, double *min, double *max, double *xy, double *integAcc, int axis, bool edge2d, int edit)
498
void TGAppMainFrame::PlotEdgeDistribution(TList *files, int filenr, int points, double *min, double *max, double *xy, double *integAcc, int axis, bool edge2d, int edit)
496
{
499
{
497
   TGraph *gScan[2];
500
   TGraph *gScan[2];
498
   int pdfmax = -1;
501
   int pdfmax = -1;
499
   int count = 0;
502
   int count = 0;
500
   char ctemp[1024];
503
   char ctemp[1024];
501
   char exportname[1024];
504
   char exportname[1024];
502
   TCanvas *gCanvas;
505
   TCanvas *gCanvas;
503
 
506
 
504
   // Prepare the CDF plot
507
   // Prepare the CDF plot
505
   gScan[1] = new TGraph();
508
   gScan[1] = new TGraph();
506
   for(int i = 0; i < points; i++)
509
   for(int i = 0; i < points; i++)
507
   {
510
   {
508
      count = filenr - points + i;
511
      count = filenr - points + i;
509
      gScan[1]->SetPoint(i, (double)xy[count], (double)integAcc[i]/(*max));
512
      gScan[1]->SetPoint(i, (double)xy[count], (double)integAcc[i]/(*max));
510
      if(DBGSIG) printf("PlotEdgeDistribution(): CDF %d: %lf, %lf\n", i, (double)xy[count], (double)integAcc[i]/(*max));
513
      if(DBGSIG) printf("PlotEdgeDistribution(): CDF %d: %lf, %lf\n", i, (double)xy[count], (double)integAcc[i]/(*max));
511
 
514
 
512
      if( ((integAcc[i+1]-integAcc[i])/(*max) > pdfmax) && (i < points-1) )
515
      if( ((integAcc[i+1]-integAcc[i])/(*max) > pdfmax) && (i < points-1) )
513
         pdfmax = (integAcc[i+1]-integAcc[i])/(*max);
516
         pdfmax = (integAcc[i+1]-integAcc[i])/(*max);
514
   }
517
   }
515
 
518
 
516
   pdfmax = (TMath::Ceil(pdfmax*10))/10.;
519
   pdfmax = (TMath::Ceil(pdfmax*10))/10.;
517
 
520
 
518
   // Prepare the PDF plot
521
   // Prepare the PDF plot
519
   gScan[0] = new TGraph();
522
   gScan[0] = new TGraph();
520
   for(int i = points-1; i >= 0; i--)
523
   for(int i = points-1; i >= 0; i--)
521
   {
524
   {
522
      count = (filenr-1) - (points-1) + i;
525
      count = (filenr-1) - (points-1) + i;
523
      // Set any negative values of the PDF to 0
526
      // Set any negative values of the PDF to 0
524
      if( (integAcc[i]-integAcc[i-1])/(*max) < 0 )
527
      if( (integAcc[i]-integAcc[i-1])/(*max) < 0 )
525
         gScan[0]->SetPoint(i, (double)xy[count], 0);
528
         gScan[0]->SetPoint(i, (double)xy[count], 0);
526
      else
529
      else
527
         gScan[0]->SetPoint(i, (double)xy[count], (integAcc[i]-integAcc[i-1])/(*max));
530
         gScan[0]->SetPoint(i, (double)xy[count], (integAcc[i]-integAcc[i-1])/(*max));
528
      if(DBGSIG) printf("PlotEdgeDistribution(): PDF %d: %lf, %lf\n", i, (double)xy[count], (integAcc[i+1]-integAcc[i])/(*max));
531
      if(DBGSIG) printf("PlotEdgeDistribution(): PDF %d: %lf, %lf\n", i, (double)xy[count], (integAcc[i+1]-integAcc[i])/(*max));
529
   }
532
   }
530
 
533
 
531
   remove_from_last((char*)files->At(filenr-1)->GetTitle(), '_', ctemp);
534
   remove_from_last((char*)files->At(filenr-1)->GetTitle(), '_', ctemp);
532
   sprintf(exportname, "%s_edge.pdf", ctemp);
535
   sprintf(exportname, "%s_edge.pdf", ctemp);
533
 
536
 
534
   if(edit == 0)
537
   if(edit == 0)
535
      gCanvas = new TCanvas("canv1d","canv1d",1200,900);
538
      gCanvas = new TCanvas("canv1d","canv1d",1200,900);
536
   else
539
   else
537
      gCanvas = tempAnalysisCanvas->GetCanvas();
540
      gCanvas = tempAnalysisCanvas->GetCanvas();
538
 
541
 
539
   // Fit the PDF with a gaussian
542
   // Fit the PDF with a gaussian
540
   gScan[0]->Fit("gaus","Q");
543
   gScan[0]->Fit("gaus","Q");
541
   gScan[0]->GetFunction("gaus")->SetNpx(400);
544
   gScan[0]->GetFunction("gaus")->SetNpx(400);
542
 
545
 
543
   gStyle->SetOptFit(1);
546
   gStyle->SetOptFit(1);
544
 
547
 
545
   gScan[1]->Draw("AL");
548
   gScan[1]->Draw("AL");
546
   gPad->Update();
549
   gPad->Update();
547
   gScan[0]->Draw("LP");
550
   gScan[0]->Draw("LP");
548
 
551
 
549
   gCanvas->Modified();
552
   gCanvas->Modified();
550
   gCanvas->Update();
553
   gCanvas->Update();
551
 
554
 
552
   TPaveStats *stats = (TPaveStats*)gScan[0]->FindObject("stats");
555
   TPaveStats *stats = (TPaveStats*)gScan[0]->FindObject("stats");
553
   if(!cleanPlots)
556
   if(!cleanPlots)
554
   {
557
   {
555
      stats->SetX1NDC(0.86); stats->SetX2NDC(1.0);
558
      stats->SetX1NDC(0.86); stats->SetX2NDC(1.0);
556
      stats->SetY1NDC(0.87); stats->SetY2NDC(1.0);
559
      stats->SetY1NDC(0.87); stats->SetY2NDC(1.0);
557
   }
560
   }
558
   else
561
   else
559
   {
562
   {
560
      stats->SetX1NDC(1.1); stats->SetX2NDC(1.3);
563
      stats->SetX1NDC(1.1); stats->SetX2NDC(1.3);
561
      stats->SetY1NDC(1.1); stats->SetY2NDC(1.3);
564
      stats->SetY1NDC(1.1); stats->SetY2NDC(1.3);
562
   }
565
   }
563
 
566
 
564
   gCanvas->SetGridx(1);
567
   gCanvas->SetGridx(1);
565
   gCanvas->SetGridy(1);
568
   gCanvas->SetGridy(1);
566
   if(axis == 1)
569
   if(axis == 1)
567
      gScan[1]->GetXaxis()->SetTitle("X [#mum]");
570
      gScan[1]->GetXaxis()->SetTitle("X [#mum]");
568
   else if(axis == 2)
571
   else if(axis == 2)
569
      gScan[1]->GetXaxis()->SetTitle("Y [#mum]");
572
      gScan[1]->GetXaxis()->SetTitle("Y [#mum]");
570
   gScan[1]->GetXaxis()->SetTitleOffset(1.3);
573
   gScan[1]->GetXaxis()->SetTitleOffset(1.3);
571
   gScan[1]->GetXaxis()->CenterTitle(kTRUE);
574
   gScan[1]->GetXaxis()->CenterTitle(kTRUE);
572
   gScan[1]->GetXaxis()->SetLabelSize(0.027);
575
   gScan[1]->GetXaxis()->SetLabelSize(0.027);
573
   gScan[1]->GetXaxis()->SetLabelOffset(0.02);
576
   gScan[1]->GetXaxis()->SetLabelOffset(0.02);
574
   gScan[1]->GetXaxis()->SetNoExponent(kTRUE);
577
   gScan[1]->GetXaxis()->SetNoExponent(kTRUE);
575
   gScan[1]->GetYaxis()->SetTitle("Normalized ADC integral");
578
   gScan[1]->GetYaxis()->SetTitle("Normalized ADC integral");
576
 
579
 
577
   gScan[1]->GetYaxis()->CenterTitle(kTRUE);
580
   gScan[1]->GetYaxis()->CenterTitle(kTRUE);
578
   gScan[1]->GetYaxis()->SetLabelSize(0.027);
581
   gScan[1]->GetYaxis()->SetLabelSize(0.027);
579
   gScan[1]->GetYaxis()->SetLabelOffset(0.02);
582
   gScan[1]->GetYaxis()->SetLabelOffset(0.02);
580
   gScan[1]->GetYaxis()->SetRangeUser(0,1);
583
   gScan[1]->GetYaxis()->SetRangeUser(0,1);
581
   gScan[1]->GetYaxis()->SetTitleOffset(1.4);
584
   gScan[1]->GetYaxis()->SetTitleOffset(1.4);
582
   gScan[1]->GetYaxis()->SetTitleSize(0.030);
585
   gScan[1]->GetYaxis()->SetTitleSize(0.030);
583
 
586
 
584
   if(!cleanPlots)
587
   if(!cleanPlots)
585
   {
588
   {
586
      if(axis == 1)
589
      if(axis == 1)
587
      {
590
      {
588
         if(posUnitsPlot->widgetCB->GetSelected() == 0)
591
         if(posUnitsPlot->widgetCB->GetSelected() == 0)
589
            gScan[1]->SetTitle("SiPM edge detection;X [table units];Normalized ADC integral");
592
            gScan[1]->SetTitle("SiPM edge detection;X [table units];Normalized ADC integral");
590
         else if(posUnitsPlot->widgetCB->GetSelected() == 1)
593
         else if(posUnitsPlot->widgetCB->GetSelected() == 1)
591
            gScan[1]->SetTitle("SiPM edge detection;X [#mum];Normalized ADC integral");
594
            gScan[1]->SetTitle("SiPM edge detection;X [#mum];Normalized ADC integral");
592
      }
595
      }
593
      else if(axis == 2)
596
      else if(axis == 2)
594
      {
597
      {
595
         if(posUnitsPlot->widgetCB->GetSelected() == 0)
598
         if(posUnitsPlot->widgetCB->GetSelected() == 0)
596
            gScan[1]->SetTitle("SiPM edge detection;Y [table units];Normalized ADC integral");
599
            gScan[1]->SetTitle("SiPM edge detection;Y [table units];Normalized ADC integral");
597
         else if(posUnitsPlot->widgetCB->GetSelected() == 1)
600
         else if(posUnitsPlot->widgetCB->GetSelected() == 1)
598
            gScan[1]->SetTitle("SiPM edge detection;Y [#mum];Normalized ADC integral");
601
            gScan[1]->SetTitle("SiPM edge detection;Y [#mum];Normalized ADC integral");
599
      }
602
      }
600
   }
603
   }
601
   else
604
   else
602
   {
605
   {
603
      if(axis == 1)
606
      if(axis == 1)
604
      {
607
      {
605
         if(posUnitsPlot->widgetCB->GetSelected() == 0)
608
         if(posUnitsPlot->widgetCB->GetSelected() == 0)
606
            gScan[1]->SetTitle(";X [table units];Normalized ADC integral");
609
            gScan[1]->SetTitle(";X [table units];Normalized ADC integral");
607
         else if(posUnitsPlot->widgetCB->GetSelected() == 1)
610
         else if(posUnitsPlot->widgetCB->GetSelected() == 1)
608
            gScan[1]->SetTitle(";X [#mum];Normalized ADC integral");
611
            gScan[1]->SetTitle(";X [#mum];Normalized ADC integral");
609
      }
612
      }
610
      else if(axis == 2)
613
      else if(axis == 2)
611
      {
614
      {
612
         if(posUnitsPlot->widgetCB->GetSelected() == 0)
615
         if(posUnitsPlot->widgetCB->GetSelected() == 0)
613
            gScan[1]->SetTitle(";Y [table units];Normalized ADC integral");
616
            gScan[1]->SetTitle(";Y [table units];Normalized ADC integral");
614
         else if(posUnitsPlot->widgetCB->GetSelected() == 1)
617
         else if(posUnitsPlot->widgetCB->GetSelected() == 1)
615
            gScan[1]->SetTitle(";Y [#mum];Normalized ADC integral");
618
            gScan[1]->SetTitle(";Y [#mum];Normalized ADC integral");
616
      }
619
      }
617
   }
620
   }
618
   gScan[1]->SetLineColor(kBlue);
621
   gScan[1]->SetLineColor(kBlue);
619
   gScan[0]->SetLineWidth(2);
622
   gScan[0]->SetLineWidth(2);
620
   gScan[1]->SetLineWidth(2);
623
   gScan[1]->SetLineWidth(2);
621
 
624
 
622
   gCanvas->Modified();
625
   gCanvas->Modified();
623
   gCanvas->Update();
626
   gCanvas->Update();
624
 
627
 
625
   gCanvas->SaveAs(exportname);
628
   gCanvas->SaveAs(exportname);
626
 
629
 
627
   // If 2D edge plot, delete the 1D edge plots as we go
630
   // If 2D edge plot, delete the 1D edge plots as we go
628
   if(edge2d)
631
   if(edge2d)
629
   {
632
   {
630
      delete gScan[0];
633
      delete gScan[0];
631
      delete gScan[1];
634
      delete gScan[1];
632
      if(edit == 0)
635
      if(edit == 0)
633
         delete gCanvas;
636
         delete gCanvas;
634
   }
637
   }
635
   else
638
   else
636
   {
639
   {
637
      if(edit == 0)
640
      if(edit == 0)
638
      {
641
      {
639
         delete gScan[0];
642
         delete gScan[0];
640
         delete gScan[1];
643
         delete gScan[1];
641
         delete gCanvas;
644
         delete gCanvas;
642
      }
645
      }
643
   }
646
   }
644
}
647
}
645
 
648
 
646
void TGAppMainFrame::PhotonMu(TList *files, int edit)
649
void TGAppMainFrame::PhotonMu(TList *files, int edit)
647
{
650
{
648
   unsigned int nrfiles = files->GetSize();
651
   unsigned int nrfiles = files->GetSize();
649
   char ctemp[1024];
652
   char ctemp[1024];
650
   char exportname[1024];
653
   char exportname[1024];
651
   int j, k = 0, m = 0, n = 0, k2 = 0, m2 = 0;
654
   int j, k = 0, m = 0, n = 0, k2 = 0, m2 = 0;
652
 
655
 
653
   TCanvas *gCanvas;
656
   TCanvas *gCanvas;
654
   TTree *header_data, *meas_data;
657
   TTree *header_data, *meas_data;
655
   double *integralCount, *integralPedestal;
658
   double *integralCount, *integralPedestal;
656
   integralCount = new double[nrfiles];
659
   integralCount = new double[nrfiles];
657
   integralPedestal = new double[nrfiles];
660
   integralPedestal = new double[nrfiles];
658
   double *angle, *pdeval, *muval;
661
   double *angle, *pdeval, *muval;
659
   angle = new double[nrfiles];
662
   angle = new double[nrfiles];
660
   pdeval = new double[nrfiles];
663
   pdeval = new double[nrfiles];
661
   muval = new double[nrfiles];
664
   muval = new double[nrfiles];
662
 
665
 
663
   // Zero the integral count and accumulated vaules
666
   // Zero the integral count and accumulated vaules
664
   for(int i = 0; i < (int)nrfiles; i++) {integralCount[i] = 0; integralPedestal[i] = 0; }
667
   for(int i = 0; i < (int)nrfiles; i++) {integralCount[i] = 0; integralPedestal[i] = 0; }
665
 
668
 
666
   // Fitting variables
669
   // Fitting variables
667
   TSpectrum *spec;
670
   TSpectrum *spec;
668
   TH1F *histtemp;
671
   TH1F *histtemp;
669
   TH1 *histback;
672
   TH1 *histback;
670
   TH1F *h2;
673
   TH1F *h2;
671
   float *xpeaks;
674
   float *xpeaks;
672
   TF1 *fit;
675
   TF1 *fit;
673
   TF1 *fittingfunc;
676
   TF1 *fittingfunc;
674
   double *fparam, *fparamerr;
677
   double *fparam, *fparamerr;
675
   double meansel[20];
678
   double meansel[20];
676
   double sigmasel[20];
679
   double sigmasel[20];
677
   double meanparam, paramsigma;
680
   double meanparam, paramsigma;
678
   int sortindex[20];
681
   int sortindex[20];
679
   int adcpedestal[2];
682
   int adcpedestal[2];
680
   int zeromu = 0;
683
   int zeromu = 0;
681
   int darkhist = -1;
684
   int darkhist = -1;
-
 
685
   int nopeaks = -1;
682
 
686
 
683
   double pointest[12];
687
   double pointest[12];
684
   bool exclude = false;
688
   bool exclude = false;
685
 
689
 
686
   // Zero the parameter values
690
   // Zero the parameter values
687
   for(int i = 0; i < 20; i++) {meansel[i] = 0; sigmasel[i] = 0; }
691
   for(int i = 0; i < 20; i++) {meansel[i] = 0; sigmasel[i] = 0; }
688
 
692
 
689
   float progVal = 0;
693
   float progVal = 0;
690
   analysisProgress->widgetPB->SetPosition(progVal);
694
   analysisProgress->widgetPB->SetPosition(progVal);
691
   gVirtualX->Update(1);
695
   gVirtualX->Update(1);
-
 
696
 
-
 
697
   // Check if the checkbox for no peaks is selected - TODO: Still need the situation when we do not have a peaked ADC spectrum
-
 
698
   if(relPde->widgetChBox[1]->IsDown())
-
 
699
   {
-
 
700
      printf("PhotonMu(): ADC spectrum has no peak structure.\n");
-
 
701
      nopeaks = 1;
-
 
702
 
-
 
703
      // Error if there is no darkhist
-
 
704
      if(strcmp("", darkRun->widgetTE->GetText()) == 0)
-
 
705
      {
-
 
706
         printf("PhotonMu(): Error! The no peak structure option needs a dark histogram.\n");
-
 
707
         delete[] integralCount;
-
 
708
         delete[] integralPedestal;
-
 
709
         delete[] angle;
-
 
710
         delete[] pdeval;
-
 
711
         delete[] muval;
-
 
712
         return;
-
 
713
      }
-
 
714
 
-
 
715
   }
692
 
716
 
693
   // Start if we select at least one file
717
   // Start if we select at least one file
694
   if(nrfiles > 0)
718
   if(nrfiles > 0)
695
   {
719
   {
-
 
720
      // Find the pedestal peak for the dark histogram, and use it for all if there are no peaks
-
 
721
      if(nopeaks != -1)
-
 
722
      {
-
 
723
         // Replot the spectrum on analysisCanvas and do not close the input file
-
 
724
         DisplayHistogram( (char*)(darkRun->widgetTE->GetText()), 0, 1);
-
 
725
         analysisCanvas->GetCanvas()->Modified();
-
 
726
         analysisCanvas->GetCanvas()->Update();
-
 
727
       
-
 
728
         // Get the spectrum
-
 
729
         histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname);
-
 
730
         npeaks = 15;
-
 
731
         double par[300];
-
 
732
         spec = new TSpectrum(npeaks);
-
 
733
         // Find spectrum background
-
 
734
         histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same");
-
 
735
         // Clone histogram and subtract background from it if we select that option
-
 
736
         h2 = (TH1F*)histtemp->Clone("h2");
-
 
737
         if(fitChecks->widgetChBox[0]->IsDown())
-
 
738
            h2->Add(histback, -1);
-
 
739
         // Search for the peaks
-
 
740
         int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() );
-
 
741
         printf("PhotonMu(): Found %d candidates to fit.\n",found);
-
 
742
         npeaks = found;
-
 
743
   
-
 
744
         // Set initial peak parameters
-
 
745
         xpeaks = spec->GetPositionX();
-
 
746
         for(j = 0; j < found; j++)
-
 
747
         {
-
 
748
            float xp = xpeaks[j];
-
 
749
            int bin = h2->GetXaxis()->FindBin(xp);
-
 
750
            float yp = h2->GetBinContent(bin);
-
 
751
            par[3*j] = yp;
-
 
752
            par[3*j+1] = xp;
-
 
753
            par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber();
-
 
754
         }
-
 
755
         
-
 
756
         // Fit the histogram
-
 
757
         fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks);
-
 
758
         TVirtualFitter::Fitter(histtemp, 3*npeaks);
-
 
759
         fit->SetParameters(par);
-
 
760
         fit->SetNpx(300);
-
 
761
         h2->Fit("fit","Q");
-
 
762
         // Get the fitted parameters
-
 
763
         fittingfunc = h2->GetFunction("fit");
-
 
764
         fparam = fittingfunc->GetParameters();
-
 
765
         fparamerr = fittingfunc->GetParErrors();
-
 
766
   
-
 
767
         // Gather the parameters (mean peak value for now)
-
 
768
         int j = 1;
-
 
769
         int nrfit = 0;
-
 
770
         while(1)
-
 
771
         {
-
 
772
            if( (fparam[j] < 1.E-30) || (nrfit > 8) )
-
 
773
               break;
-
 
774
            else
-
 
775
            {
-
 
776
               // Check if pedestal is above the lower limit and sigma is smaller than the mean
-
 
777
               if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) )
-
 
778
               {
-
 
779
                  // With the additional ADC offset, we can shift the mean values slightly, so they are not close to the X.5, but to the X.0 values
-
 
780
                  meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber());
-
 
781
                  sigmasel[nrfit] = fparam[j+1];
-
 
782
                  nrfit++;
-
 
783
               }
-
 
784
            }
-
 
785
   
-
 
786
            j+=3;
-
 
787
         }
-
 
788
         TMath::Sort(nrfit, meansel, sortindex, kFALSE);
-
 
789
 
-
 
790
         meanparam = meansel[sortindex[0]];
-
 
791
         paramsigma = sigmasel[sortindex[0]];
-
 
792
 
-
 
793
         for(j = 0; j < nrfit; j++)
-
 
794
            if(DBGSIG)
-
 
795
               printf("PhotonMu(): %d: peak mean = %lf\n", j, meansel[sortindex[j]]);
-
 
796
       
-
 
797
         j = 0;
-
 
798
         adcpedestal[0] = 0;
-
 
799
         adcpedestal[1] = -1;
-
 
800
 
-
 
801
         while(1)
-
 
802
         {
-
 
803
            int bin = histtemp->GetXaxis()->FindBin((int)(j+meanparam+paramsigma));
-
 
804
            int yp = histtemp->GetBinContent(bin);
-
 
805
 
-
 
806
            // Check where we get to first minimum after pedestal peak or where we get to the half maximum of the pedestal peak (in case there is only a pedestal peak)
-
 
807
            if(adcpedestal[1] == -1)
-
 
808
            {
-
 
809
               adcpedestal[0] = j+meanparam+paramsigma;
-
 
810
               adcpedestal[1] = yp;
-
 
811
            }
-
 
812
            else
-
 
813
            {
-
 
814
               if( (npeaks > 1) && (adcpedestal[1] >= yp) )
-
 
815
               {
-
 
816
                  adcpedestal[0] = j+meanparam+paramsigma;
-
 
817
                  adcpedestal[1] = yp;
-
 
818
               }
-
 
819
               else if( (npeaks == 1) && (adcpedestal[0] < meanparam+5*paramsigma) ) // TODO -> Determining the pedestal when only one peak
-
 
820
               {
-
 
821
                  adcpedestal[0] = j+meanparam+paramsigma;
-
 
822
                  adcpedestal[1] = yp;
-
 
823
               }
-
 
824
               else
-
 
825
                  break;
-
 
826
            }
-
 
827
       
-
 
828
            j++;
-
 
829
            if(j > 50) break;
-
 
830
         }
-
 
831
 
-
 
832
         if(midPeak->widgetChBox[0]->IsDown())
-
 
833
         {
-
 
834
            if( (meanparam - (int)meanparam >= 0.) && (meanparam - (int)meanparam < 0.5) )
-
 
835
               m = TMath::Floor(meanparam);
-
 
836
            else if( (meanparam - (int)meanparam >= 0.5) && (meanparam - (int)meanparam < 1.) )
-
 
837
               m = TMath::Ceil(meanparam);
-
 
838
            int bin = histtemp->GetXaxis()->FindBin(m);
-
 
839
            adcpedestal[0] = m;
-
 
840
            printf("midpeak x = %d, ", adcpedestal[0]);
-
 
841
            adcpedestal[1] = histtemp->GetBinContent(bin);
-
 
842
         }
-
 
843
 
-
 
844
         // Option to show the fit
-
 
845
         fittingfunc->Draw("L SAME");
-
 
846
         analysisCanvas->GetCanvas()->Modified();
-
 
847
         analysisCanvas->GetCanvas()->Update();
-
 
848
       
-
 
849
         printf("Pedestal ends = %d and nr. of counts = %d\n", adcpedestal[0], adcpedestal[1]);
-
 
850
 
-
 
851
         // Delete the opened histogram and spectrum
-
 
852
         delete spec;
-
 
853
         delete inroot;
-
 
854
 
-
 
855
//       return;
-
 
856
      }
-
 
857
 
-
 
858
      printf("PhotonMu(): Continuing with the rest of the spectra.\n");
-
 
859
 
-
 
860
      // Check all histograms for pedestal peak values
696
      for(int i = 0; i < (int)nrfiles; i++)
861
      for(int i = 0; i < (int)nrfiles; i++)
697
      {
862
      {
698
         if( (nrfiles == 1) || (!multiSelect->widgetChBox[0]->IsDown()) )
863
         if( (nrfiles == 1) || (!multiSelect->widgetChBox[0]->IsDown()) )
699
         {
864
         {
700
            printf("PhotonMu(): Only one file selected. Not running analysis, just showing the fit.\n");
865
            printf("PhotonMu(): Only one file selected. Not running analysis, just showing the fit.\n");
701
 
866
 
702
            // Replot the spectrum on analysisCanvas and do not close the input file
867
            // Replot the spectrum on analysisCanvas and do not close the input file
703
            DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1);
868
            DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1);
704
            analysisCanvas->GetCanvas()->Modified();
869
            analysisCanvas->GetCanvas()->Modified();
705
            analysisCanvas->GetCanvas()->Update();
870
            analysisCanvas->GetCanvas()->Update();
706
       
871
       
707
            // Get the spectrum
872
            // Get the spectrum
708
            histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname);
873
            histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname);
709
            npeaks = 15;
874
            npeaks = 15;
710
            double par[300];
875
            double par[300];
711
            spec = new TSpectrum(npeaks);
876
            spec = new TSpectrum(npeaks);
712
            // Find spectrum background
877
            // Find spectrum background
713
            histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same");
878
            histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same");
714
            // Clone histogram and subtract background from it if we select that option
879
            // Clone histogram and subtract background from it if we select that option
715
            h2 = (TH1F*)histtemp->Clone("h2");
880
            h2 = (TH1F*)histtemp->Clone("h2");
716
            if(fitChecks->widgetChBox[0]->IsDown())
881
            if(fitChecks->widgetChBox[0]->IsDown())
717
               h2->Add(histback, -1);
882
               h2->Add(histback, -1);
718
            // Search for the peaks
883
            // Search for the peaks
719
            int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() );
884
            int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() );
720
            printf("PhotonMu(): Found %d candidates to fit.\n",found);
885
            printf("PhotonMu(): Found %d candidates to fit.\n",found);
721
            npeaks = found;
886
            npeaks = found;
722
   
887
   
723
            // Set initial peak parameters
888
            // Set initial peak parameters
724
            xpeaks = spec->GetPositionX();
889
            xpeaks = spec->GetPositionX();
725
            for(j = 0; j < found; j++)
890
            for(j = 0; j < found; j++)
726
            {
891
            {
727
               float xp = xpeaks[j];
892
               float xp = xpeaks[j];
728
               int bin = h2->GetXaxis()->FindBin(xp);
893
               int bin = h2->GetXaxis()->FindBin(xp);
729
               float yp = h2->GetBinContent(bin);
894
               float yp = h2->GetBinContent(bin);
730
               par[3*j] = yp;
895
               par[3*j] = yp;
731
               par[3*j+1] = xp;
896
               par[3*j+1] = xp;
732
               par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber();
897
               par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber();
733
            }
898
            }
734
         
899
         
735
            // Fit the histogram
900
            // Fit the histogram
736
            fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks);
901
            fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks);
737
            TVirtualFitter::Fitter(histtemp, 3*npeaks);
902
            TVirtualFitter::Fitter(histtemp, 3*npeaks);
738
            fit->SetParameters(par);
903
            fit->SetParameters(par);
739
            fit->SetNpx(300);
904
            fit->SetNpx(300);
740
            h2->Fit("fit","Q");
905
            h2->Fit("fit","Q");
741
            // Get the fitted parameters
906
            // Get the fitted parameters
742
            fittingfunc = h2->GetFunction("fit");
907
            fittingfunc = h2->GetFunction("fit");
743
            fparam = fittingfunc->GetParameters();
908
            fparam = fittingfunc->GetParameters();
744
            fparamerr = fittingfunc->GetParErrors();
909
            fparamerr = fittingfunc->GetParErrors();
745
   
910
   
746
            // Gather the parameters (mean peak value for now)
911
            // Gather the parameters (mean peak value for now)
747
            int j = 1;
912
            int j = 1;
748
            int nrfit = 0;
913
            int nrfit = 0;
749
            while(1)
914
            while(1)
750
            {
915
            {
751
               if( (fparam[j] < 1.E-30) || (nrfit > 8) )
916
               if( (fparam[j] < 1.E-30) || (nrfit > 8) )
752
                  break;
917
                  break;
753
               else
918
               else
754
               {
919
               {
755
                  // Check if pedestal is above the lower limit and sigma is smaller than the mean
920
                  // Check if pedestal is above the lower limit and sigma is smaller than the mean
756
                  if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) )
921
                  if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) )
757
                  {
922
                  {
758
                     // With the additional ADC offset, we can shift the mean values slightly, so they are not close to the X.5, but to the X.0 values
923
                     // With the additional ADC offset, we can shift the mean values slightly, so they are not close to the X.5, but to the X.0 values
759
                     meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber());
924
                     meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber());
760
                     sigmasel[nrfit] = fparam[j+1];
925
                     sigmasel[nrfit] = fparam[j+1];
761
                     nrfit++;
926
                     nrfit++;
762
                  }
927
                  }
763
               }
928
               }
764
   
929
   
765
               j+=3;
930
               j+=3;
766
            }
931
            }
767
            TMath::Sort(nrfit, meansel, sortindex, kFALSE);
932
            TMath::Sort(nrfit, meansel, sortindex, kFALSE);
768
 
933
 
769
            fittingfunc->Draw("SAME");
934
            fittingfunc->Draw("SAME");
770
            analysisCanvas->GetCanvas()->Modified();
-
 
771
            analysisCanvas->GetCanvas()->Update();
-
 
772
 
-
 
773
            meanparam = meansel[sortindex[0]];
-
 
774
            paramsigma = sigmasel[sortindex[0]];
-
 
775
 
-
 
776
            for(j = 0; j < nrfit; j++)
-
 
777
               printf("PhotonMu(): %d: peak mean = %lf\n", j, meansel[sortindex[j]]);
-
 
778
 
-
 
779
 
-
 
780
 
-
 
781
            return;
-
 
782
         }
-
 
783
         if(files->At(i))
-
 
784
         {
-
 
785
            if(strcmp(files->At(i)->GetTitle(),darkRun->widgetTE->GetText()) == 0)
-
 
786
            {
-
 
787
               printf("PhotonMu(): %s is the dark histogram file.\n", files->At(i)->GetTitle());
-
 
788
               darkhist = i;
-
 
789
            }
-
 
790
 
-
 
791
            // Replot the spectrum on analysisCanvas and do not close the input file
-
 
792
            DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1);
-
 
793
            analysisCanvas->GetCanvas()->Modified();
935
            analysisCanvas->GetCanvas()->Modified();
794
            analysisCanvas->GetCanvas()->Update();
936
            analysisCanvas->GetCanvas()->Update();
795
       
-
 
796
            // Get the spectrum
-
 
797
            histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname);
-
 
798
            npeaks = 15;
-
 
799
            double par[300];
-
 
800
            spec = new TSpectrum(npeaks);
-
 
801
            // Find spectrum background
-
 
802
            histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same");
-
 
803
            // Clone histogram and subtract background from it if we select that option
-
 
804
            h2 = (TH1F*)histtemp->Clone("h2");
-
 
805
            if(fitChecks->widgetChBox[0]->IsDown())
-
 
806
               h2->Add(histback, -1);
-
 
807
            // Search for the peaks
-
 
808
            int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() );
-
 
809
            printf("PhotonMu(): Found %d candidates to fit.\n",found);
-
 
810
            npeaks = found;
-
 
811
   
-
 
812
            // Set initial peak parameters
-
 
813
            xpeaks = spec->GetPositionX();
-
 
814
            for(j = 0; j < found; j++)
-
 
815
            {
-
 
816
               float xp = xpeaks[j];
-
 
817
               int bin = h2->GetXaxis()->FindBin(xp);
-
 
818
               float yp = h2->GetBinContent(bin);
-
 
819
               par[3*j] = yp;
-
 
820
               par[3*j+1] = xp;
-
 
821
               par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber();
-
 
822
            }
-
 
823
         
-
 
824
            // Fit the histogram
-
 
825
            fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks);
-
 
826
            TVirtualFitter::Fitter(histtemp, 3*npeaks);
-
 
827
            fit->SetParameters(par);
-
 
828
            fit->SetNpx(300);
-
 
829
            h2->Fit("fit","Q");
-
 
830
            // Get the fitted parameters
-
 
831
            fittingfunc = h2->GetFunction("fit");
-
 
832
            fparam = fittingfunc->GetParameters();
-
 
833
            fparamerr = fittingfunc->GetParErrors();
-
 
834
   
-
 
835
            // Gather the parameters (mean peak value for now)
-
 
836
            int j = 1;
-
 
837
            int nrfit = 0;
-
 
838
            while(1)
-
 
839
            {
-
 
840
               if( (fparam[j] < 1.E-30) || (nrfit > 8) )
-
 
841
                  break;
-
 
842
               else
-
 
843
               {
-
 
844
                  // Check if pedestal is above the lower limit and sigma is smaller than the mean
-
 
845
                  if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) )
-
 
846
                  {
-
 
847
                     // With the additional ADC offset, we can shift the mean values slightly, so they are not close to the X.5, but to the X.0 values
-
 
848
                     meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber());
-
 
849
                     sigmasel[nrfit] = fparam[j+1];
-
 
850
                     nrfit++;
-
 
851
                  }
-
 
852
               }
-
 
853
   
-
 
854
               j+=3;
-
 
855
            }
-
 
856
            TMath::Sort(nrfit, meansel, sortindex, kFALSE);
-
 
857
 
937
 
858
            meanparam = meansel[sortindex[0]];
938
            meanparam = meansel[sortindex[0]];
859
            paramsigma = sigmasel[sortindex[0]];
939
            paramsigma = sigmasel[sortindex[0]];
860
 
940
 
861
            for(j = 0; j < nrfit; j++)
941
            for(j = 0; j < nrfit; j++)
862
               if(DBGSIG)
-
 
863
                  printf("PhotonMu(): %d: peak mean = %lf\n", j, meansel[sortindex[j]]);
942
               printf("PhotonMu(): %d: peak mean = %lf\n", j, meansel[sortindex[j]]);
-
 
943
 
-
 
944
 
-
 
945
 
864
       
946
            return;
-
 
947
         }
865
            j = 0;
948
         if(files->At(i))
-
 
949
         {
866
            adcpedestal[0] = 0;
950
            if(strcmp(files->At(i)->GetTitle(),darkRun->widgetTE->GetText()) == 0)
-
 
951
            {
-
 
952
               printf("PhotonMu(): %s is the dark histogram file.\n", files->At(i)->GetTitle());
867
            adcpedestal[1] = -1;
953
               darkhist = i;
-
 
954
            }
868
 
955
 
869
            while(1)
956
            if(nopeaks == -1)
870
            {
957
            {
-
 
958
               // Replot the spectrum on analysisCanvas and do not close the input file
-
 
959
               DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1);
871
               int bin = histtemp->GetXaxis()->FindBin((int)(j+meanparam+paramsigma));
960
               analysisCanvas->GetCanvas()->Modified();
872
               int yp = histtemp->GetBinContent(bin);
961
               analysisCanvas->GetCanvas()->Update();
873
 
962
       
-
 
963
               // Get the spectrum
874
               // 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)
964
               histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname);
875
               if(adcpedestal[1] == -1)
965
               npeaks = 15;
876
               {
966
               double par[300];
877
                  adcpedestal[0] = j+meanparam+paramsigma;
967
               spec = new TSpectrum(npeaks);
-
 
968
               // Find spectrum background
-
 
969
               histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same");
-
 
970
               // Clone histogram and subtract background from it if we select that option
-
 
971
               h2 = (TH1F*)histtemp->Clone("h2");
-
 
972
               if(fitChecks->widgetChBox[0]->IsDown())
878
                  adcpedestal[1] = yp;
973
                  h2->Add(histback, -1);
879
               }
974
               // Search for the peaks
-
 
975
               int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() );
-
 
976
               printf("PhotonMu(): Found %d candidates to fit.\n",found);
880
               else
977
               npeaks = found;
-
 
978
   
-
 
979
               // Set initial peak parameters
-
 
980
               xpeaks = spec->GetPositionX();
-
 
981
               for(j = 0; j < found; j++)
881
               {
982
               {
-
 
983
                  float xp = xpeaks[j];
-
 
984
                  int bin = h2->GetXaxis()->FindBin(xp);
-
 
985
                  float yp = h2->GetBinContent(bin);
-
 
986
                  par[3*j] = yp;
-
 
987
                  par[3*j+1] = xp;
-
 
988
                  par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber();
-
 
989
               }
-
 
990
         
-
 
991
               // Fit the histogram
-
 
992
               fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks);
-
 
993
               TVirtualFitter::Fitter(histtemp, 3*npeaks);
-
 
994
               fit->SetParameters(par);
-
 
995
               fit->SetNpx(300);
-
 
996
               h2->Fit("fit","Q");
-
 
997
               // Get the fitted parameters
-
 
998
               fittingfunc = h2->GetFunction("fit");
-
 
999
               fparam = fittingfunc->GetParameters();
-
 
1000
               fparamerr = fittingfunc->GetParErrors();
-
 
1001
   
-
 
1002
               // Gather the parameters (mean peak value for now)
-
 
1003
               int j = 1;
-
 
1004
               int nrfit = 0;
-
 
1005
               while(1)
-
 
1006
               {
882
                  if( (npeaks > 1) && (adcpedestal[1] >= yp) )
1007
                  if( (fparam[j] < 1.E-30) || (nrfit > 8) )
-
 
1008
                     break;
-
 
1009
                  else
883
                  {
1010
                  {
-
 
1011
                     // Check if pedestal is above the lower limit and sigma is smaller than the mean
-
 
1012
                     if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) )
-
 
1013
                     {
-
 
1014
                        // With the additional ADC offset, we can shift the mean values slightly, so they are not close to the X.5, but to the X.0 values
-
 
1015
                        meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber());
884
                     adcpedestal[0] = j+meanparam+paramsigma;
1016
                        sigmasel[nrfit] = fparam[j+1];
885
                     adcpedestal[1] = yp;
1017
                        nrfit++;
-
 
1018
                     }
886
                  }
1019
                  }
-
 
1020
   
-
 
1021
                  j+=3;
-
 
1022
               }
-
 
1023
               TMath::Sort(nrfit, meansel, sortindex, kFALSE);
-
 
1024
 
-
 
1025
               meanparam = meansel[sortindex[0]];
-
 
1026
               paramsigma = sigmasel[sortindex[0]];
-
 
1027
 
-
 
1028
               for(j = 0; j < nrfit; j++)
-
 
1029
                  if(DBGSIG)
-
 
1030
                     printf("PhotonMu(): %d: peak mean = %lf\n", j, meansel[sortindex[j]]);
-
 
1031
       
-
 
1032
               j = 0;
-
 
1033
               adcpedestal[0] = 0;
-
 
1034
               adcpedestal[1] = -1;
-
 
1035
 
-
 
1036
               while(1)
-
 
1037
               {
-
 
1038
                  int bin = histtemp->GetXaxis()->FindBin((int)(j+meanparam+paramsigma));
-
 
1039
                  int yp = histtemp->GetBinContent(bin);
-
 
1040
 
887
                  else if( (npeaks == 1) && (adcpedestal[0] < meanparam+5*paramsigma) ) // TODO -> Determining the pedestal when only one peak
1041
                  // Check where we get to first minimum after pedestal peak or where we get to the half maximum of the pedestal peak (in case there is only a pedestal peak)
-
 
1042
                  if(adcpedestal[1] == -1)
888
                  {
1043
                  {
889
                     adcpedestal[0] = j+meanparam+paramsigma;
1044
                     adcpedestal[0] = j+meanparam+paramsigma;
890
                     adcpedestal[1] = yp;
1045
                     adcpedestal[1] = yp;
891
                  }
1046
                  }
892
                  else
1047
                  else
-
 
1048
                  {
-
 
1049
                     if( (npeaks > 1) && (adcpedestal[1] >= yp) )
-
 
1050
                     {
-
 
1051
                        adcpedestal[0] = j+meanparam+paramsigma;
-
 
1052
                        adcpedestal[1] = yp;
-
 
1053
                     }
-
 
1054
                     else if( (npeaks == 1) && (adcpedestal[0] < meanparam+5*paramsigma) ) // TODO -> Determining the pedestal when only one peak
-
 
1055
                     {
-
 
1056
                        adcpedestal[0] = j+meanparam+paramsigma;
-
 
1057
                        adcpedestal[1] = yp;
-
 
1058
                     }
-
 
1059
                     else
893
                     break;
1060
                        break;
-
 
1061
                  }
-
 
1062
       
-
 
1063
                  j++;
-
 
1064
                  if(j > 50) break;
894
               }
1065
               }
895
       
-
 
896
               j++;
-
 
897
               if(j > 50) break;
-
 
898
            }
-
 
899
 
-
 
900
            if( (npeaks > 1) && (nrfit > 1) )
-
 
901
            {
-
 
902
               int bin = histtemp->GetXaxis()->FindBin((int)(meanparam+meansel[sortindex[1]])/2);
-
 
903
               adcpedestal[0] = (meanparam+meansel[sortindex[1]])/2;
-
 
904
               printf("PhotonMu(): multipeak x = %d, ", adcpedestal[0]);
-
 
905
               adcpedestal[1] = histtemp->GetBinContent(bin);
-
 
906
            }
-
 
907
 
-
 
908
            if(midPeak->widgetChBox[0]->IsDown())
-
 
909
            {
-
 
910
               if( (meanparam - (int)meanparam >= 0.) && (meanparam - (int)meanparam < 0.5) )
-
 
911
                  m = TMath::Floor(meanparam);
-
 
912
               else if( (meanparam - (int)meanparam >= 0.5) && (meanparam - (int)meanparam < 1.) )
-
 
913
                  m = TMath::Ceil(meanparam);
-
 
914
               int bin = histtemp->GetXaxis()->FindBin(m);
-
 
915
               adcpedestal[0] = m;
-
 
916
               printf("midpeak x = %d, ", adcpedestal[0]);
-
 
917
               adcpedestal[1] = histtemp->GetBinContent(bin);
-
 
918
            }
-
 
919
 
1066
 
-
 
1067
               if(npeaks > 1)
-
 
1068
               {
-
 
1069
                  int bin = histtemp->GetXaxis()->FindBin((int)(meanparam+meansel[sortindex[1]])/2);
-
 
1070
                  adcpedestal[0] = (meanparam+meansel[sortindex[1]])/2;
-
 
1071
                  printf("PhotonMu(): multipeak x = %d, ", adcpedestal[0]);
-
 
1072
                  adcpedestal[1] = histtemp->GetBinContent(bin);
-
 
1073
               }
-
 
1074
 
-
 
1075
               if(midPeak->widgetChBox[0]->IsDown())
-
 
1076
               {
-
 
1077
                  if( (meanparam - (int)meanparam >= 0.) && (meanparam - (int)meanparam < 0.5) )
-
 
1078
                     m = TMath::Floor(meanparam);
-
 
1079
                  else if( (meanparam - (int)meanparam >= 0.5) && (meanparam - (int)meanparam < 1.) )
-
 
1080
                     m = TMath::Ceil(meanparam);
-
 
1081
                  int bin = histtemp->GetXaxis()->FindBin(m);
-
 
1082
                  adcpedestal[0] = m;
-
 
1083
                  printf("midpeak x = %d, ", adcpedestal[0]);
-
 
1084
                  adcpedestal[1] = histtemp->GetBinContent(bin);
-
 
1085
               }
-
 
1086
 
920
/*          // Option to show the fit
1087
/*             // Option to show the fit
921
            fittingfunc->Draw("L SAME");
1088
               fittingfunc->Draw("L SAME");
922
            analysisCanvas->GetCanvas()->Modified();
1089
               analysisCanvas->GetCanvas()->Modified();
923
            analysisCanvas->GetCanvas()->Update();*/
1090
               analysisCanvas->GetCanvas()->Update();*/
924
       
1091
       
925
            printf("Pedestal ends = %d and nr. of counts = %d\n", adcpedestal[0], adcpedestal[1]);
1092
               printf("Pedestal ends = %d and nr. of counts = %d\n", adcpedestal[0], adcpedestal[1]);
926
 
1093
 
927
            // Delete the opened histogram and spectrum
1094
               // Delete the opened histogram and spectrum
928
            delete spec;
1095
               delete spec;
929
            delete inroot;
1096
               delete inroot;
-
 
1097
            }
930
 
1098
 
931
            // Open the input file and read header, ADC and TDC values
1099
            // Open the input file and read header, ADC and TDC values
932
            sprintf(ctemp, "%s", files->At(i)->GetTitle());
1100
            sprintf(ctemp, "%s", files->At(i)->GetTitle());
933
            inroot = new TFile(ctemp, "READ");
1101
            inroot = new TFile(ctemp, "READ");
934
         
1102
         
935
            inroot->GetObject("header_data", header_data);
1103
            inroot->GetObject("header_data", header_data);
936
            inroot->GetObject("meas_data", meas_data);
1104
            inroot->GetObject("meas_data", meas_data);
937
         
1105
         
938
            // Reading the header
1106
            // Reading the header
939
            if( header_data->FindBranch("angle") )
1107
            if( header_data->FindBranch("angle") )
940
            {
1108
            {
941
               header_data->SetBranchAddress("angle", &evtheader.angle);
1109
               header_data->SetBranchAddress("angle", &evtheader.angle);
942
               header_data->GetEntry(0);
1110
               header_data->GetEntry(0);
943
            }
1111
            }
944
            else
1112
            else
945
            {
1113
            {
946
               printf("PhotonMu(): Error! Selected file has no angle header value. Please edit header to add the angle header value.\n");
1114
               printf("PhotonMu(): Error! Selected file has no angle header value. Please edit header to add the angle header value.\n");
947
               break;
1115
               break;
948
            }
1116
            }
949
   
1117
   
950
            char rdc[256];
1118
            char rdc[256];
951
            j = selectCh->widgetNE[0]->GetNumber();
1119
            j = selectCh->widgetNE[0]->GetNumber();
952
            double rangetdc[2];
1120
            double rangetdc[2];
953
            rangetdc[0] = tdcRange->widgetNE[0]->GetNumber();
1121
            rangetdc[0] = tdcRange->widgetNE[0]->GetNumber();
954
            rangetdc[1] = tdcRange->widgetNE[1]->GetNumber();
1122
            rangetdc[1] = tdcRange->widgetNE[1]->GetNumber();
955
   
1123
   
956
            k = 0;
1124
            k = 0;
957
            k2 = 0;
1125
            k2 = 0;
958
            m = 0;
1126
            m = 0;
959
            m2 = 0;
1127
            m2 = 0;
960
 
1128
 
961
            // Reading the data
1129
            // Reading the data
962
            for(int e = 0; e < meas_data->GetEntries(); e++)
1130
            for(int e = 0; e < meas_data->GetEntries(); e++)
963
            {
1131
            {
964
               sprintf(rdc, "ADC%d", j);
1132
               sprintf(rdc, "ADC%d", j);
965
               meas_data->SetBranchAddress(rdc, &evtdata.adcdata[j]);
1133
               meas_data->SetBranchAddress(rdc, &evtdata.adcdata[j]);
966
               meas_data->GetEntry(e);
1134
               meas_data->GetEntry(e);
967
         
1135
         
968
               sprintf(rdc, "TDC%d", j);
1136
               sprintf(rdc, "TDC%d", j);
969
               meas_data->SetBranchAddress(rdc, &evtdata.tdcdata[j]);
1137
               meas_data->SetBranchAddress(rdc, &evtdata.tdcdata[j]);
970
               meas_data->GetEntry(e);
1138
               meas_data->GetEntry(e);
971
   
1139
   
972
               // If our data point is inside the TDC window
1140
               // If our data point is inside the TDC window
973
               if( ((double)evtdata.tdcdata[j]/tdctimeconversion >= rangetdc[0]) && ((double)evtdata.tdcdata[j]/tdctimeconversion <= rangetdc[1]) )
1141
               if( ((double)evtdata.tdcdata[j]/tdctimeconversion >= rangetdc[0]) && ((double)evtdata.tdcdata[j]/tdctimeconversion <= rangetdc[1]) )
974
               {
1142
               {
975
                  // Gather only the integral of the pedestal
1143
                  // Gather only the integral of the pedestal
976
                  if((double)evtdata.adcdata[j] < (double)adcpedestal[0]+0.5 )
1144
                  if((double)evtdata.adcdata[j] < (double)adcpedestal[0]+0.5 )
977
                  {
1145
                  {
978
                     k2++;
1146
                     k2++;
979
                     m2 += evtdata.adcdata[j];
1147
                     m2 += evtdata.adcdata[j];
980
                  }
1148
                  }
981
 
1149
 
982
                  // Gather the complete integral
1150
                  // Gather the complete integral
983
                  k++;
1151
                  k++;
984
                  m += evtdata.adcdata[j];
1152
                  m += evtdata.adcdata[j];
985
               }
1153
               }
986
            }
1154
            }
987
 
1155
 
988
            // Determine the angle, mu and relative PDE
1156
            // Determine the angle, mu and relative PDE
989
            angle[i] = (double)(evtheader.angle);
1157
            angle[i] = (double)(evtheader.angle);
990
            integralCount[i] += (double)m;
1158
            integralCount[i] += (double)m;
991
            printf("PhotonMu(): %lf: Complete integral (%d evts) = %lf\n", angle[i], k, integralCount[i]);
1159
            printf("PhotonMu(): %lf: Complete integral (%d evts) = %lf\n", angle[i], k, integralCount[i]);
992
            integralPedestal[i] += (double)m2;
1160
            integralPedestal[i] += (double)m2;
993
            printf("PhotonMu(): %lf: Pedestal integral (%d evts) = %lf\n", angle[i], k2, integralPedestal[i]);
1161
            printf("PhotonMu(): %lf: Pedestal integral (%d evts) = %lf\n", angle[i], k2, integralPedestal[i]);
994
            if( (angle[i] == zeroAngle->widgetNE[0]->GetNumber()) && (darkhist != i) )
1162
            if( (angle[i] == zeroAngle->widgetNE[0]->GetNumber()) && (darkhist != i) )
995
               zeromu = i;
1163
               zeromu = i;
996
 
1164
 
997
            // Checking for errors when fitting a histogram
-
 
998
            if(k2 == 0)
-
 
999
            {
-
 
1000
               printf("PhotonMu(): No pedestal entries found. Check the fitting results.\n");
-
 
1001
               muval[i] = -1;
-
 
1002
            }
-
 
1003
            else
-
 
1004
               muval[i] = -TMath::Log((double)k2/(double)k);
1165
            muval[i] = -TMath::Log((double)k2/(double)k);
1005
            printf("PhotonMu(): %lf: muval = %lf\n", angle[i], muval[i]);
1166
            printf("PhotonMu(): %lf: muval = %lf\n", angle[i], muval[i]);
1006
 
1167
 
1007
            inroot->Close();
1168
            inroot->Close();
1008
            delete inroot;
1169
            delete inroot;
1009
         }
1170
         }
1010
 
1171
 
1011
         // Update the progress bar
1172
         // Update the progress bar
1012
         progVal = (float)(90.00/nrfiles)*i;
1173
         progVal = (float)(90.00/nrfiles)*i;
1013
         analysisProgress->widgetPB->SetPosition(progVal);
1174
         analysisProgress->widgetPB->SetPosition(progVal);
1014
         gVirtualX->Update(1);
1175
         gVirtualX->Update(1);
1015
      }
1176
      }
1016
 
1177
 
1017
      printf("PhotonMu(): %d files were selected.\n", nrfiles);
1178
      printf("PhotonMu(): %d files were selected.\n", nrfiles);
1018
 
1179
 
1019
      printf("PhotonMu(): angle\tmu\trelative PDE\n");
1180
      printf("PhotonMu(): angle\tmu\trelative PDE\n");
1020
      m = 0;
1181
      m = 0;
1021
     
1182
 
1022
      // Set the 0 degree muval, reuse meansel[1]
1183
      // Set the 0 degree muval, reuse meansel[1]
1023
      meansel[1] = muval[zeromu];
1184
      meansel[1] = muval[zeromu];
1024
      printf("Zero value (id=%d, angle=%lf) = %lf\n", zeromu, angle[zeromu], meansel[1]);
1185
      printf("Zero value (id=%d, angle=%lf) = %lf\n", zeromu, angle[zeromu], meansel[1]);
1025
 
1186
 
1026
      // TODO - point estimation still not working correctly!
1187
      // TODO - point estimation still not working correctly!
1027
      for(int i = 0; i < (int)nrfiles; i++)
1188
      for(int i = 0; i < (int)nrfiles; i++)
1028
      {
1189
      {
1029
         // Estimate next point and check error (5 point least square fit estimation)
1190
         // Estimate next point and check error (5 point least square fit estimation)
1030
         if( ((i > 4) && (m == i)) || ((i > 5) && (m < i)) )
1191
         if( ((i > 4) && (m == i)) || ((i > 5) && (m < i)) )
1031
         {
1192
         {
1032
            // Set exclude signal to false
1193
            // Set exclude signal to false
1033
            exclude = false;
1194
            exclude = false;
1034
 
1195
 
1035
            // Get next point values (if zero value -> need to add the dark hist value again)
1196
            // Get next point values (if zero value -> need to add the dark hist value again)
1036
            pointest[10] = angle[i];
1197
            pointest[10] = angle[i];
1037
            pointest[11] = muval[i];
1198
            pointest[11] = muval[i];
1038
 
1199
 
1039
            // Check if next point has larger error than acceptable (if yes, set exclude signal to true), reuse meansel[0]
1200
            // Check if next point has larger error than acceptable (if yes, set exclude signal to true), reuse meansel[0]
1040
            meansel[0] = PointEstimate(5, pointest);    // PointEstimate only works with very small step size
1201
            meansel[0] = PointEstimate(5, pointest);    // PointEstimate only works with very small step size
1041
            if(meansel[0] > accError->widgetNE[0]->GetNumber())
1202
            if(meansel[0] > accError->widgetNE[0]->GetNumber())
1042
            {
1203
            {
1043
               printf("PhotonMu(): Point (%lf, %lf) excluded due to error: %lf\n", pointest[10], pointest[11], meansel[0]);
1204
               printf("PhotonMu(): Point (%lf, %lf) excluded due to error: %lf\n", pointest[10], pointest[11], meansel[0]);
1044
               exclude = true;
1205
               exclude = true;
1045
            }
1206
            }
1046
 
1207
 
1047
            // Value with 0 angle and dark histogram are always needed, so should not be excluded
1208
            // Value with 0 angle and dark histogram are always needed, so should not be excluded
1048
            if(i == darkhist)
1209
            if(i == darkhist)
1049
               exclude = false;
1210
               exclude = false;
1050
 
-
 
1051
            // Wrong fit
-
 
1052
            if(muval[i] == -1)
-
 
1053
               exclude = true;
-
 
1054
 
1211
 
1055
            // If nothing excluded, pass the points in pointest variable like in a FIFO
1212
            // If nothing excluded, pass the points in pointest variable like in a FIFO
1056
            if(!exclude)
1213
            if(!exclude)
1057
            {
1214
            {
1058
               for(int j = 0; j < 10; j++)
1215
               for(int j = 0; j < 10; j++)
1059
               {
1216
               {
1060
                  if(DBGSIG) printf("PhotonMu(): j = %d: Old X value = %lf\n", j, pointest[j]);
1217
                  if(DBGSIG) printf("PhotonMu(): j = %d: Old X value = %lf\n", j, pointest[j]);
1061
                  pointest[j] = pointest[j+2];
1218
                  pointest[j] = pointest[j+2];
1062
               }
1219
               }
1063
            }
1220
            }
1064
            else
1221
            else
1065
            {
1222
            {
1066
               for(int j = 0; j < 10; j++)
1223
               for(int j = 0; j < 10; j++)
1067
                  if(DBGSIG) printf("PhotonMu(): No j = %d: Old X value = %lf\n", j, pointest[j]);
1224
                  if(DBGSIG) printf("PhotonMu(): No j = %d: Old X value = %lf\n", j, pointest[j]);
1068
            }
1225
            }
1069
         }
1226
         }
1070
         else
1227
         else
1071
         {
1228
         {
1072
            // First 5 points act as estimator points for next one
1229
            // First 5 points act as estimator points for next one
1073
            pointest[2*m] = angle[i];
1230
            pointest[2*m] = angle[i];
1074
            pointest[2*m+1] = muval[i];
1231
            pointest[2*m+1] = muval[i];
1075
         }
1232
         }
1076
 
1233
 
1077
         // Run only if we have a dark run histogram and middle pedestal peak estimation
1234
         // Run only if we have a dark run histogram and middle pedestal peak estimation
1078
         if( (darkhist != -1) && midPeak->widgetChBox[0]->IsDown() )
1235
         if( (darkhist != -1) && midPeak->widgetChBox[0]->IsDown() )
1079
         {
1236
         {
1080
            if(DBGSIG) printf("PhotonMu(): m = %d, i = %d: muval = %lf, ", m, i, muval[i]);
1237
            if(DBGSIG) printf("PhotonMu(): m = %d, i = %d: muval = %lf, ", m, i, muval[i]);
1081
 
1238
 
1082
            // Subtract the dark value from all values
1239
            // Subtract the dark value from all values
1083
            angle[m] = angle[i];
1240
            angle[m] = angle[i];
1084
            muval[m] = muval[i] - muval[darkhist];
1241
            muval[m] = muval[i] - muval[darkhist];
1085
 
1242
 
1086
            if(DBGSIG) printf("angle = %lf, newmuval = %lf, darkmuval = %lf, ", angle[m], muval[m], muval[darkhist]);
1243
            if(DBGSIG) printf("angle = %lf, newmuval = %lf, darkmuval = %lf, ", angle[m], muval[m], muval[darkhist]);
1087
 
1244
 
1088
            // Calculate relative PDE
1245
            // Calculate relative PDE
1089
//          pdeval[m] = muval[m]/(muval[zeromu]*TMath::Cos(angle[m]*TMath::ACos(-1.)/180.));
1246
//          pdeval[m] = muval[m]/(muval[zeromu]*TMath::Cos(angle[m]*TMath::ACos(-1.)/180.));
1090
            pdeval[m] = muval[m]/((meansel[1]-muval[darkhist])*TMath::Cos(angle[m]*TMath::ACos(-1.)/180.));
1247
            pdeval[m] = muval[m]/((meansel[1]-muval[darkhist])*TMath::Cos(angle[m]*TMath::ACos(-1.)/180.));
1091
            if(DBGSIG) printf("pdeval = %lf\n", pdeval[m]);
1248
            if(DBGSIG) printf("pdeval = %lf\n", pdeval[m]);
1092
 
1249
 
1093
            // Only increase counter if error is not too big
1250
            // Only increase counter if error is not too big
1094
            if( (darkhist != i) && (!exclude) )
1251
            if( (darkhist != i) && (!exclude) )
1095
               m++;
1252
               m++;
1096
         }
1253
         }
1097
         else
1254
         else
1098
         {
1255
         {
1099
            // Relative PDE calculation
1256
            // Relative PDE calculation
1100
            angle[m] = angle[i];
1257
            angle[m] = angle[i];
1101
            muval[m] = muval[i];
1258
            muval[m] = muval[i];
1102
//            pdeval[m] = muval[m]/(muval[zeromu]*TMath::Cos(angle[m]*TMath::ACos(-1.)/180.));
1259
//            pdeval[m] = muval[m]/(muval[zeromu]*TMath::Cos(angle[m]*TMath::ACos(-1.)/180.));
1103
            pdeval[m] = muval[m]/(meansel[1]*TMath::Cos(angle[m]*TMath::ACos(-1.)/180.));
1260
            pdeval[m] = muval[m]/(meansel[1]*TMath::Cos(angle[m]*TMath::ACos(-1.)/180.));
1104
 
1261
 
1105
            // Only increase counter if error is not too big
1262
            // Only increase counter if error is not too big
1106
            if(!exclude)
1263
            if(!exclude)
1107
               m++;
1264
               m++;
1108
         }
1265
         }
1109
         printf("PhotonMu(): %lf\t%lf\t%lf\n", angle[i], muval[i], pdeval[i]);
1266
         printf("PhotonMu(): %lf\t%lf\t%lf\n", angle[i], muval[i], pdeval[i]);
1110
      }
1267
      }
1111
 
1268
 
1112
      // Check for range of values to plot
1269
      // Check for range of values to plot
1113
      double plotMax = 0.;
1270
      double plotMax = 0.;
1114
      for(int i = 0; i < (int)nrfiles; i++)
1271
      for(int i = 0; i < (int)nrfiles; i++)
1115
      {
1272
      {
1116
         plotMax = TMath::Max(plotMax, muval[i]);
1273
         plotMax = TMath::Max(plotMax, muval[i]);
1117
         plotMax = TMath::Max(plotMax, pdeval[i]);
1274
         plotMax = TMath::Max(plotMax, pdeval[i]);
1118
      }
1275
      }
1119
      if(plotMax <= 0.)
1276
      if(plotMax <= 0.)
1120
         plotMax = 1.1;
1277
         plotMax = 1.1;
1121
      printf("PhotonMu(): Maximum value: %lf\n", plotMax);
1278
      printf("PhotonMu(): Maximum value: %lf\n", plotMax);
1122
 
1279
 
1123
      if(DBGSIG) printf("\n");
1280
      if(DBGSIG) printf("\n");
1124
      if(darkhist != -1)
1281
      if(darkhist != -1)
1125
         printf("PhotonMu(): Number of excluded points: %d\n", (nrfiles-1-m));
1282
         printf("PhotonMu(): Number of excluded points: %d\n", (nrfiles-1-m));
1126
      else
1283
      else
1127
         printf("PhotonMu(): Number of excluded points: %d\n", (nrfiles-m));
1284
         printf("PhotonMu(): Number of excluded points: %d\n", (nrfiles-m));
1128
 
1285
 
1129
      // Plot mu and PDE angle dependance plots
1286
      // Plot mu and PDE angle dependance plots
1130
      if(edit == 0)
1287
      if(edit == 0)
1131
         gCanvas = new TCanvas("canv","canv",1200,900);
1288
         gCanvas = new TCanvas("canv","canv",1200,900);
1132
      else if(edit == 1)
1289
      else if(edit == 1)
1133
         gCanvas = tempAnalysisCanvas->GetCanvas();
1290
         gCanvas = tempAnalysisCanvas->GetCanvas();
1134
      gCanvas->cd();
1291
      gCanvas->cd();
1135
      gCanvas->SetGrid();
1292
      gCanvas->SetGrid();
1136
 
1293
 
1137
      TGraph *pde = new TGraph(m, angle, pdeval);
1294
      TGraph *pde = new TGraph(m, angle, pdeval);
1138
      pde->SetMarkerStyle(21);
1295
      pde->SetMarkerStyle(21);
1139
      pde->SetMarkerSize(0.7);
1296
      pde->SetMarkerSize(0.7);
1140
      pde->SetMarkerColor(2);
1297
      pde->SetMarkerColor(2);
1141
      pde->SetLineWidth(1);
1298
      pde->SetLineWidth(1);
1142
      pde->SetLineColor(2);
1299
      pde->SetLineColor(2);
1143
      pde->GetXaxis()->SetLabelSize(0.030);
1300
      pde->GetXaxis()->SetLabelSize(0.030);
1144
      pde->GetXaxis()->CenterTitle();
1301
      pde->GetXaxis()->CenterTitle();
1145
//      pde->GetXaxis()->SetRange(angle[0],angle[nrfiles-1]);
1302
//      pde->GetXaxis()->SetRange(angle[0],angle[nrfiles-1]);
1146
//      pde->GetXaxis()->SetRangeUser(angle[0],angle[nrfiles-1]);
1303
//      pde->GetXaxis()->SetRangeUser(angle[0],angle[nrfiles-1]);
1147
      pde->GetXaxis()->SetRange(-90.0,90.0);
1304
      pde->GetXaxis()->SetRange(-90.0,90.0);
1148
      pde->GetXaxis()->SetRangeUser(-90.0,90.0);
1305
      pde->GetXaxis()->SetRangeUser(-90.0,90.0);
1149
      pde->GetXaxis()->SetLimits(-90.0,90.0);
1306
      pde->GetXaxis()->SetLimits(-90.0,90.0);
1150
      pde->GetYaxis()->SetTitleOffset(1.2);
1307
      pde->GetYaxis()->SetTitleOffset(1.2);
1151
      pde->GetYaxis()->SetLabelSize(0.030);
1308
      pde->GetYaxis()->SetLabelSize(0.030);
1152
      pde->GetYaxis()->CenterTitle();
1309
      pde->GetYaxis()->CenterTitle();
1153
      pde->GetYaxis()->SetRange(0., 1.1*plotMax);
1310
      pde->GetYaxis()->SetRange(0., 1.1*plotMax);
1154
      pde->GetYaxis()->SetRangeUser(0., 1.1*plotMax);
1311
      pde->GetYaxis()->SetRangeUser(0., 1.1*plotMax);
1155
      pde->GetYaxis()->SetLimits(0., 1.1*plotMax);
1312
      pde->GetYaxis()->SetLimits(0., 1.1*plotMax);
1156
      pde->SetName("pde");
1313
      pde->SetName("pde");
1157
      pde->Draw("ALP");
1314
      pde->Draw("ALP");
1158
 
1315
 
1159
      pde->SetTitle(";Incidence angle (#circ);Relative PDE(#theta) / #mu(#theta)");
1316
      pde->SetTitle(";Incidence angle (#circ);Relative PDE(#theta) / #mu(#theta)");
1160
 
1317
 
1161
      TGraph *mugr = new TGraph(m, angle, muval);
1318
      TGraph *mugr = new TGraph(m, angle, muval);
1162
      mugr->SetMarkerStyle(20);
1319
      mugr->SetMarkerStyle(20);
1163
      mugr->SetMarkerSize(0.7);
1320
      mugr->SetMarkerSize(0.7);
1164
      mugr->SetMarkerColor(4);
1321
      mugr->SetMarkerColor(4);
1165
      mugr->SetLineWidth(1);
1322
      mugr->SetLineWidth(1);
1166
      mugr->SetLineColor(4);
1323
      mugr->SetLineColor(4);
1167
      mugr->SetName("muval");
1324
      mugr->SetName("muval");
1168
      mugr->Draw("SAME;LP");
1325
      mugr->Draw("SAME;LP");
1169
 
1326
 
1170
      gCanvas->Modified();
1327
      gCanvas->Modified();
1171
      gCanvas->Update();
1328
      gCanvas->Update();
1172
 
1329
 
1173
      if(edit == 0)
1330
      if(edit == 0)
1174
      {
1331
      {
1175
         remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp);
1332
         remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp);
1176
         sprintf(exportname, "%s_pde.pdf", ctemp);
1333
         sprintf(exportname, "%s_pde.pdf", ctemp);
1177
         gCanvas->SaveAs(exportname);
1334
         gCanvas->SaveAs(exportname);
1178
 
1335
 
1179
         delete mugr;
1336
         delete mugr;
1180
         delete pde;
1337
         delete pde;
1181
         delete gCanvas;
1338
         delete gCanvas;
1182
      }
1339
      }
1183
      else if(edit == 1)
1340
      else if(edit == 1)
1184
      {
1341
      {
1185
         gCanvas->Modified();
1342
         gCanvas->Modified();
1186
         gCanvas->Update();
1343
         gCanvas->Update();
1187
      }
1344
      }
1188
 
1345
 
1189
      // Update the progress bar
1346
      // Update the progress bar
1190
      analysisProgress->widgetPB->SetPosition(100.);
1347
      analysisProgress->widgetPB->SetPosition(100.);
1191
      gVirtualX->Update(1);
1348
      gVirtualX->Update(1);
1192
   }
1349
   }
1193
}
1350
}
1194
 
1351
 
1195
void TGAppMainFrame::BreakdownVolt(TList *files, int edit)
1352
void TGAppMainFrame::BreakdownVolt(TList *files, int edit)
1196
{
1353
{
1197
   unsigned int nrfiles = files->GetSize();
1354
   unsigned int nrfiles = files->GetSize();
1198
   char ctemp[1024];
1355
   char ctemp[1024];
1199
   char exportname[1024];
1356
   char exportname[1024];
1200
   char paramname[1024];
1357
   char paramname[1024];
1201
   int j, k = 0;
1358
   int j, k = 0;
1202
 
1359
 
1203
   TCanvas *gCanvas;
1360
   TCanvas *gCanvas;
1204
   TTree *header_data, *meas_data;
1361
   TTree *header_data, *meas_data;
1205
 
1362
 
1206
   // Fitting variables
1363
   // Fitting variables
1207
   TSpectrum *spec;
1364
   TSpectrum *spec;
1208
   TH1F *histtemp;
1365
   TH1F *histtemp;
1209
   TH1 *histback;
1366
   TH1 *histback;
1210
   TH1F *h2;
1367
   TH1F *h2;
1211
   float *xpeaks;
1368
   float *xpeaks;
1212
   TF1 *fit;
1369
   TF1 *fit;
1213
   TF1 *fittingfunc;
1370
   TF1 *fittingfunc;
1214
   TLatex *latex;
1371
   TLatex *latex;
1215
   double *fparam, *fparamerr;
1372
   double *fparam, *fparamerr;
1216
   double meansel[20];
1373
   double meansel[20];
1217
   double meanselerr[20];
1374
   double meanselerr[20];
1218
   double sigmasel[20];
1375
   double sigmasel[20];
1219
   double meanparam, meanparamerr, paramsigma;
1376
   double meanparam, meanparamerr, paramsigma;
1220
   int sortindex[20];
1377
   int sortindex[20];
1221
   bool exclude = false;
1378
   bool exclude = false;
1222
 
1379
 
1223
   int p = 0;
1380
   int p = 0;
1224
   double volt[nrfiles], volterr[nrfiles], sep[3][nrfiles], seperr[3][nrfiles];
1381
   double volt[nrfiles], volterr[nrfiles], sep[3][nrfiles], seperr[3][nrfiles];
1225
   int first = 1;
1382
   int first = 1;
1226
 
1383
 
1227
   FILE *fp;
1384
   FILE *fp;
1228
   remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp);
1385
   remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp);
1229
   sprintf(paramname, "%s_fitresult.txt", ctemp);
1386
   sprintf(paramname, "%s_fitresult.txt", ctemp);
1230
   fp = fopen(paramname, "w");
1387
   fp = fopen(paramname, "w");
1231
   fclose(fp);
1388
   fclose(fp);
1232
 
1389
 
1233
   int peaklimit = minPeak->widgetNE[0]->GetNumber()+1; // +1 to account for the pedestal peak
1390
   int peaklimit = minPeak->widgetNE[0]->GetNumber()+1; // +1 to account for the pedestal peak
1234
 
1391
 
1235
   // Zero the parameter values
1392
   // Zero the parameter values
1236
   for(int i = 0; i < nrfiles; i++)
1393
   for(int i = 0; i < nrfiles; i++)
1237
   {
1394
   {
1238
      volt[i] = 0; volterr[i] = 0;
1395
      volt[i] = 0; volterr[i] = 0;
1239
      sep[0][i] = 0; sep[1][i] = 0; sep[2][i] = 0;
1396
      sep[0][i] = 0; sep[1][i] = 0; sep[2][i] = 0;
1240
      seperr[0][i] = 0; seperr[1][i] = 0; seperr[2][i] = 0;
1397
      seperr[0][i] = 0; seperr[1][i] = 0; seperr[2][i] = 0;
1241
      if(i < 20) { meansel[i] = 0; meanselerr[i] = 0; sigmasel[i] = 0; }
1398
      if(i < 20) { meansel[i] = 0; meanselerr[i] = 0; sigmasel[i] = 0; }
1242
   }
1399
   }
1243
 
1400
 
1244
   float progVal = 0;
1401
   float progVal = 0;
1245
   analysisProgress->widgetPB->SetPosition(progVal);
1402
   analysisProgress->widgetPB->SetPosition(progVal);
1246
   gVirtualX->Update(1);
1403
   gVirtualX->Update(1);
1247
 
1404
 
1248
   // Start if we select at least one file
1405
   // Start if we select at least one file
1249
   if(nrfiles > 0)
1406
   if(nrfiles > 0)
1250
   {
1407
   {
1251
      for(int i = 0; i < (int)nrfiles; i++)
1408
      for(int i = 0; i < (int)nrfiles; i++)
1252
      {
1409
      {
1253
         if( (nrfiles == 1) || (!multiSelect->widgetChBox[0]->IsDown()) )
1410
         if( (nrfiles == 1) || (!multiSelect->widgetChBox[0]->IsDown()) )
1254
         {
1411
         {
1255
            printf("BreakdownVolt(): Only one file selected. Not running analysis, just showing the fit.\n");
1412
            printf("BreakdownVolt(): Only one file selected. Not running analysis, just showing the fit.\n");
1256
 
1413
 
1257
            // Replot the spectrum on analysisCanvas and do not close the input file
1414
            // Replot the spectrum on analysisCanvas and do not close the input file
1258
            DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1);
1415
            DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1);
1259
            analysisCanvas->GetCanvas()->Modified();
1416
            analysisCanvas->GetCanvas()->Modified();
1260
            analysisCanvas->GetCanvas()->Update();
1417
            analysisCanvas->GetCanvas()->Update();
1261
       
1418
       
1262
            // Get the spectrum
1419
            // Get the spectrum
1263
            histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname);
1420
            histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname);
1264
            npeaks = 15;
1421
            npeaks = 15;
1265
            double par[300];
1422
            double par[300];
1266
            spec = new TSpectrum(npeaks);
1423
            spec = new TSpectrum(npeaks);
1267
            // Find spectrum background
1424
            // Find spectrum background
1268
            histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same");
1425
            histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same");
1269
            // Clone histogram and subtract background from it if we select that option
1426
            // Clone histogram and subtract background from it if we select that option
1270
            h2 = (TH1F*)histtemp->Clone("h2");
1427
            h2 = (TH1F*)histtemp->Clone("h2");
1271
            if(fitChecks->widgetChBox[0]->IsDown())
1428
            if(fitChecks->widgetChBox[0]->IsDown())
1272
               h2->Add(histback, -1);
1429
               h2->Add(histback, -1);
1273
            // Search for the peaks
1430
            // Search for the peaks
1274
            int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() );
1431
            int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() );
1275
            printf("PhotonMu(): Found %d candidates to fit.\n",found);
1432
            printf("PhotonMu(): Found %d candidates to fit.\n",found);
1276
            npeaks = found;
1433
            npeaks = found;
1277
   
1434
   
1278
            // Set initial peak parameters
1435
            // Set initial peak parameters
1279
            xpeaks = spec->GetPositionX();
1436
            xpeaks = spec->GetPositionX();
1280
            for(j = 0; j < found; j++)
1437
            for(j = 0; j < found; j++)
1281
            {
1438
            {
1282
               float xp = xpeaks[j];
1439
               float xp = xpeaks[j];
1283
               int bin = h2->GetXaxis()->FindBin(xp);
1440
               int bin = h2->GetXaxis()->FindBin(xp);
1284
               float yp = h2->GetBinContent(bin);
1441
               float yp = h2->GetBinContent(bin);
1285
               par[3*j] = yp;
1442
               par[3*j] = yp;
1286
               par[3*j+1] = xp;
1443
               par[3*j+1] = xp;
1287
               par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber();
1444
               par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber();
1288
            }
1445
            }
1289
         
1446
         
1290
            // Fit the histogram
1447
            // Fit the histogram
1291
            fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks);
1448
            fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks);
1292
            TVirtualFitter::Fitter(histtemp, 3*npeaks);
1449
            TVirtualFitter::Fitter(histtemp, 3*npeaks);
1293
            fit->SetParameters(par);
1450
            fit->SetParameters(par);
1294
            fit->SetNpx(300);
1451
            fit->SetNpx(300);
1295
            h2->Fit("fit","Q");
1452
            h2->Fit("fit","Q");
1296
            // Get the fitted parameters
1453
            // Get the fitted parameters
1297
            fittingfunc = h2->GetFunction("fit");
1454
            fittingfunc = h2->GetFunction("fit");
1298
            fparam = fittingfunc->GetParameters