Rev 173 |
Blame |
Compare with Previous |
Last modification |
View Log
| RSS feed
#include "../include/sipmscan.h"
#include "../include/workstation.h"
#include <stdio.h>
#include <stdlib.h>
// Peak detection function
int npeaks = 20;
double FindPeaks(double *x, double *par)
{
double result = 0;
for(int i = 0; i < npeaks; i++)
{
double norm = par[3*i];
double mean = par[3*i+1];
double sigma = par[3*i+2];
result += norm*TMath::Gaus(x[0], mean, sigma);
}
return result;
}
// File browser for selecting the dark run histogram
void TGAppMainFrame::SelectDarkHist()
{
TGFileInfo file_info;
const char *filetypes[] = {"Histograms",histextall,0,0};
char *cTemp;
file_info.fFileTypes = filetypes;
// cTemp = new char[1024];
// sprintf(cTemp, "%s/results", rootdir);
// file_info.fIniDir = StrDup(cTemp);
file_info.fIniDir = StrDup(currentAnalDir);
file_info.fMultipleSelection = kFALSE;
new TGFileDialog(gClient->GetDefaultRoot(), fMain, kFDOpen, &file_info);
// delete[] cTemp;
if(file_info.fFilename != NULL)
{
darkRun->widgetTE->SetText(file_info.fFilename);
fileList->AddEntry(file_info.fFilename, fileList->GetNumberOfEntries());
fileList->Layout();
}
else
darkRun->widgetTE->SetText("");
}
// Reset analysis defaults for the current analysis type (0 = Integrate spectrum, 1 = Relative PDE, 2 = Breakdown voltage, 3 = Surface scan, 4 = Timing analysis)
void TGAppMainFrame::AnalysisDefaults()
{
printf("AnalysisDefaults(): Current analysis tab = %d\n", analTab->GetCurrent());
if(analTab->GetCurrent() == 0) // Integrate spectrum
{
intSpect->widgetChBox[0]->SetState(kButtonUp);
intSpect->widgetChBox[1]->SetState(kButtonUp);
intSpect->widgetChBox[2]->SetState(kButtonDown);
resol2d->widgetNE[0]->SetNumber(40);
resol2d->widgetNE[1]->SetNumber(40);
}
else if(analTab->GetCurrent() == 1) // Relative PDE
{
relPde->widgetChBox[0]->SetState(kButtonDown);
midPeak->widgetChBox[0]->SetState(kButtonUp);
zeroAngle->widgetNE[0]->SetNumber(0.00);
}
else if(analTab->GetCurrent() == 2) // Breakdown voltage
{
minPeak->widgetNE[0]->SetNumber(2);
minPeak->widgetNE[0]->SetNumber(1);
}
else if(analTab->GetCurrent() == 3) // Surface scan
{
surfScanOpt->widgetChBox[0]->SetState(kButtonDown);
surfScanOpt->widgetChBox[1]->SetState(kButtonUp);
resolSurf->widgetNE[0]->SetNumber(40);
resolSurf->widgetNE[1]->SetNumber(40);
}
}
// Analysis handle function
void TGAppMainFrame::AnalysisHandle(int type)
{
TList *files;
bool createTab = true;
int tabid = -1;
int analtab = analTab->GetCurrent();
int analtype;
if( (analtab == 0) && (!intSpect->widgetChBox[0]->IsDown() && !intSpect->widgetChBox[1]->IsDown()) )
analtype = 0;
else if( (analtab == 0) && (intSpect->widgetChBox[0]->IsDown() || intSpect->widgetChBox[1]->IsDown()) )
analtype = 1;
else if( analtab == 1 )
analtype = 2;
else if( analtab == 2 )
analtype = 3;
else if( analtab == 3 )
analtype = 4;
// Save analysis settings any time we run a new analysis
SaveAnalSettings();
// Only integrate spectrum or make relative PDE
if(type == 0)
{
files = new TList();
fileList->GetSelectedEntries(files);
if( analtype == 0 )
IntegSpectrum(files, 0, 0);
if( intSpect->widgetChBox[0]->IsDown() )
IntegSpectrum(files, 1, 0);
if( intSpect->widgetChBox[1]->IsDown() )
IntegSpectrum(files, 2, 0);
if( analtype == 2 )
PhotonMu(files, 0);
if( analtype == 3 )
BreakdownVolt(files, 0);
if( analtype == 4 )
SurfaceScan(files, 0);
}
// Integrate spectrum or make relative PDE and open edit window
else if(type == 1)
{
files = new TList();
fileList->GetSelectedEntries(files);
// Prepare a new analysis edit tab, if we want to edit plots
for(int i = 0; i < fTab->GetNumberOfTabs(); i++)
{
if(strcmp("Analysis edit", fTab->GetTabTab(i)->GetString() ) == 0)
{
createTab = false;
tabid = i;
}
if(DBGSIG) printf("AnalysisHandle(): Name of tab = %s\n", fTab->GetTabTab(i)->GetString() );
}
if(files->GetSize() > 0)
{
TempAnalysisTab(fTab, createTab, &tabid, analtype);
// Integrate spectra
if( analtype == 0 )
IntegSpectrum(files, 0, 1);
if( intSpect->widgetChBox[0]->IsDown() )
IntegSpectrum(files, 1, 1);
if( intSpect->widgetChBox[1]->IsDown() )
IntegSpectrum(files, 2, 1);
if( analtype == 2 )
PhotonMu(files, 1);
if( analtype == 3 )
BreakdownVolt(files, 1);
if( analtype == 4 )
SurfaceScan(files, 1);
fTab->SetTab(tabid);
}
}
delete files;
}
// Analysis functions ----------------------------------
void TGAppMainFrame::IntegSpectrum(TList *files, int direction, int edit)
{
unsigned int nrfiles = files->GetSize();
char ctemp[1024];
char exportname[1024];
int j, k = 0, m = 0;
TTree *header_data, *meas_data;
double *integralCount, *integralAcc;
integralCount = new double[nrfiles];
integralAcc = new double[nrfiles];
double *surfxy, *surfz;
surfxy = new double[nrfiles];
surfz = new double[nrfiles];
double minInteg, maxInteg;
bool norm = intSpect->widgetChBox[2]->IsDown();
double curzval;
bool edge2d = false;
TCanvas *gCanvas;
float progVal = 0;
analysisProgress->widgetPB->SetPosition(progVal);
gVirtualX->Update(1);
// Zero the integral count and accumulated vaules
for(int i = 0; i < (int)nrfiles; i++) {integralCount[i] = 0; integralAcc[i] = 0; }
// Start if we select at least one file
if(nrfiles > 0)
{
for(int i = 0; i < (int)nrfiles; i++)
{
if(files->At(i))
{
// Read the X,Y and Z positions from header and ADC and TDC values from the measurements
sprintf(ctemp, "%s", files->At(i)->GetTitle());
inroot = new TFile(ctemp, "READ");
inroot->GetObject("header_data", header_data);
inroot->GetObject("meas_data", meas_data);
header_data->SetBranchAddress("xpos", &evtheader.xpos);
header_data->GetEntry(0);
header_data->SetBranchAddress("ypos", &evtheader.ypos);
header_data->GetEntry(0);
header_data->SetBranchAddress("zpos", &evtheader.zpos);
header_data->GetEntry(0);
char rdc[256];
j = selectCh->widgetNE[0]->GetNumber();
double rangetdc[2];
rangetdc[0] = tdcRange->widgetNE[0]->GetNumber();
rangetdc[1] = tdcRange->widgetNE[1]->GetNumber();
k = 0;
m = 0;
// Reading the data
for(int e = 0; e < meas_data->GetEntries(); e++)
{
sprintf(rdc, "ADC%d", j);
meas_data->SetBranchAddress(rdc, &evtdata.adcdata[j]);
meas_data->GetEntry(e);
sprintf(rdc, "TDC%d", j);
meas_data->SetBranchAddress(rdc, &evtdata.tdcdata[j]);
meas_data->GetEntry(e);
// Use data point only if it is inside the TDC window
if( ((double)evtdata.tdcdata[j]/tdctimeconversion >= rangetdc[0]) && ((double)evtdata.tdcdata[j]/tdctimeconversion <= rangetdc[1]) )
{
k++;
m += evtdata.adcdata[j];
}
}
// X, Y and Z values from each file (table units or microns)
if(posUnitsPlot->widgetCB->GetSelected() == 0)
{
if(direction == 1)
surfxy[i] = (double)(evtheader.xpos);
else if(direction == 2)
surfxy[i] = (double)(evtheader.ypos);
surfz[i] = (double)(evtheader.zpos);
}
else if(posUnitsPlot->widgetCB->GetSelected() == 1)
{
if(direction == 1)
surfxy[i] = (double)(evtheader.xpos*lenconversion);
else if(direction == 2)
surfxy[i] = (double)(evtheader.ypos*lenconversion);
surfz[i] = (double)(evtheader.zpos*lenconversion);
}
// 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
if(i == 0) curzval = surfz[i];
else
{
if(surfz[i] != curzval)
edge2d = true;
}
// Print the calculated integral, if no X or Y direction edge plots are enabled; otherwise, just save for later plotting
if(direction == 0)
{
if(norm)
{
integralCount[i] += ((double)m)/((double)k);
printf("IntegSpectrum(): %s: The integral is: %lf\n", ctemp, integralCount[i]);
}
else
{
integralCount[i] += m;
printf("IntegSpectrum(): %s: The integral is: %d\n", ctemp, (int)integralCount[i]);
}
}
else
{
if(norm)
integralCount[i] += ((double)m)/((double)k);
else
integralCount[i] += m;
}
inroot->Close();
delete inroot;
}
// Update the progress bar
progVal = (float)(75.00/nrfiles)*i;
analysisProgress->widgetPB->SetPosition(progVal);
gVirtualX->Update(1);
}
printf("IntegSpectrum(): %d files were selected.\n", nrfiles);
// If only an integral is needed, do not plot and exit here
if( direction == 0 )
{
delete[] integralCount;
delete[] surfxy;
delete[] surfz;
return;
}
// Current z value and the accumulated counter
curzval = surfz[0];
j = 0;
int acc = 0;
int zb;
for(int i = 0; i <= (int)nrfiles; i++)
{
// Collect the accumulated integral in order to produce a PDF from a CDF
// While we are at the same Z value, save under one set
if( (surfz[i] == curzval) && (acc != nrfiles) )
{
integralAcc[j] = integralCount[i];
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]);
j++;
acc++;
}
// 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
else
{
// Find minimal and maximal integral values to subtract the offset and normate PDF to 1
NormateSet(i, j, &minInteg, &maxInteg, integralCount, integralAcc);
if(acc != nrfiles)
{
curzval = surfz[i];
// PDF and CDF plot
PlotEdgeDistribution(files, i, j, &minInteg, &maxInteg, surfxy, integralAcc, direction, edge2d, edit);
i--;
j = 0;
}
else
{
// PDF and CDF plot
PlotEdgeDistribution(files, i, j, &minInteg, &maxInteg, surfxy, integralAcc, direction, edge2d, edit);
i--;
break;
}
}
// Update the progress bar
progVal = (float)(15.00/nrfiles)*i+75.00;
analysisProgress->widgetPB->SetPosition(progVal);
gVirtualX->Update(1);
}
// Make the 2D edge plot
if(edge2d)
{
if(edit == 0)
gCanvas = new TCanvas("canv","canv",900,900);
else
gCanvas = tempAnalysisCanvas->GetCanvas();
double range[4];
TGraph2D *gScan2D;
gScan2D = new TGraph2D();
range[0] = TMath::MinElement(nrfiles, surfxy);
range[1] = TMath::MaxElement(nrfiles, surfxy);
range[2] = TMath::MinElement(nrfiles, surfz);
range[3] = TMath::MaxElement(nrfiles, surfz);
for(int i = 0; i < nrfiles; i++)
{
if(DBGSIG)
printf("IntegSpectrum(): %.2lf\t%.2lf\t%lf\n", surfxy[i], surfz[i], integralCount[i]);
gScan2D->SetPoint(i, surfxy[i], surfz[i], integralCount[i]);
// Update the progress bar
progVal = (float)(9.00/nrfiles)*i+90.00;
analysisProgress->widgetPB->SetPosition(progVal);
gVirtualX->Update(1);
}
gCanvas->cd();
gStyle->SetPalette(1);
gCanvas->SetLeftMargin(0.15);
gCanvas->SetRightMargin(0.126);
gCanvas->SetTopMargin(0.077);
gScan2D->Draw("COLZ");
gCanvas->Modified();
gCanvas->Update();
gScan2D->SetNpx((int)resol2d->widgetNE[0]->GetNumber());
gScan2D->SetNpy((int)resol2d->widgetNE[1]->GetNumber());
gCanvas->Modified();
gCanvas->Update();
if(direction == 1)
gScan2D->GetXaxis()->SetTitle("X [#mum]");
else if(direction == 2)
gScan2D->GetXaxis()->SetTitle("Y [#mum]");
gScan2D->GetXaxis()->SetTitleOffset(1.3);
gScan2D->GetXaxis()->CenterTitle(kTRUE);
gScan2D->GetXaxis()->SetLabelSize(0.027);
gScan2D->GetXaxis()->SetLabelOffset(0.02);
gScan2D->GetXaxis()->SetRangeUser(range[0], range[1]);
gScan2D->GetXaxis()->SetNoExponent(kTRUE);
gScan2D->GetYaxis()->SetTitleOffset(1.9);
gScan2D->GetYaxis()->CenterTitle(kTRUE);
gScan2D->GetYaxis()->SetLabelSize(0.027);
gScan2D->GetYaxis()->SetLabelOffset(0.02);
gScan2D->GetYaxis()->SetRangeUser(range[2], range[3]);
gScan2D->GetYaxis()->SetNoExponent(kTRUE);
/* TGaxis *yax = (TGaxis*)gScan2D->GetYaxis();
yax->SetMaxDigits(4);*/
if(!cleanPlots)
{
if(direction == 1)
{
if(posUnitsPlot->widgetCB->GetSelected() == 0)
gScan2D->SetTitle("Laser focal point;X [table units];Z [table units]");
else if(posUnitsPlot->widgetCB->GetSelected() == 1)
gScan2D->SetTitle("Laser focal point;X [#mum];Z [#mum]");
}
else if(direction == 2)
{
if(posUnitsPlot->widgetCB->GetSelected() == 0)
gScan2D->SetTitle("Laser focal point;Y [table units];Z [table units]");
else if(posUnitsPlot->widgetCB->GetSelected() == 1)
gScan2D->SetTitle("Laser focal point;Y [#mum];Z [#mum]");
}
}
else
{
if(direction == 1)
{
if(posUnitsPlot->widgetCB->GetSelected() == 0)
gScan2D->SetTitle(";X [table units];Z [table units]");
else if(posUnitsPlot->widgetCB->GetSelected() == 1)
gScan2D->SetTitle(";X [#mum];Z [#mum]");
}
else if(direction == 2)
{
if(posUnitsPlot->widgetCB->GetSelected() == 0)
gScan2D->SetTitle(";Y [table units];Z [table units]");
else if(posUnitsPlot->widgetCB->GetSelected() == 1)
gScan2D->SetTitle(";Y [#mum];Z [#mum]");
}
}
gCanvas->Modified();
gCanvas->Update();
remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp);
sprintf(exportname, "%s", ctemp);
remove_from_last(exportname, '_', ctemp);
if(direction == 1)
sprintf(exportname, "%s_xdir_focalpoint.pdf", ctemp);
else if(direction == 2)
sprintf(exportname, "%s_ydir_focalpoint.pdf", ctemp);
gCanvas->SaveAs(exportname);
// Update the progress bar
analysisProgress->widgetPB->SetPosition(100.0);
gVirtualX->Update(1);
if(edit == 0)
{
delete gScan2D;
delete gCanvas;
}
}
else
{
// Update the progress bar
analysisProgress->widgetPB->SetPosition(100.0);
gVirtualX->Update(1);
}
}
}
void TGAppMainFrame::PlotEdgeDistribution(TList *files, int filenr, int points, double *min, double *max, double *xy, double *integAcc, int axis, bool edge2d, int edit)
{
TGraph *gScan[2];
int pdfmax = -1;
int count = 0;
char ctemp[1024];
char exportname[1024];
TCanvas *gCanvas;
// Prepare the CDF plot
gScan[1] = new TGraph();
for(int i = 0; i < points; i++)
{
count = filenr - points + i;
gScan[1]->SetPoint(i, (double)xy[count], (double)integAcc[i]/(*max));
if(DBGSIG) printf("PlotEdgeDistribution(): CDF %d: %lf, %lf\n", i, (double)xy[count], (double)integAcc[i]/(*max));
if( ((integAcc[i+1]-integAcc[i])/(*max) > pdfmax) && (i < points-1) )
pdfmax = (integAcc[i+1]-integAcc[i])/(*max);
}
pdfmax = (TMath::Ceil(pdfmax*10))/10.;
// Prepare the PDF plot
gScan[0] = new TGraph();
for(int i = points-1; i >= 0; i--)
{
count = (filenr-1) - (points-1) + i;
// Set any negative values of the PDF to 0
if( (integAcc[i]-integAcc[i-1])/(*max) < 0 )
gScan[0]->SetPoint(i, (double)xy[count], 0);
else
gScan[0]->SetPoint(i, (double)xy[count], (integAcc[i]-integAcc[i-1])/(*max));
if(DBGSIG) printf("PlotEdgeDistribution(): PDF %d: %lf, %lf\n", i, (double)xy[count], (integAcc[i+1]-integAcc[i])/(*max));
}
remove_from_last((char*)files->At(filenr-1)->GetTitle(), '_', ctemp);
sprintf(exportname, "%s_edge.pdf", ctemp);
if(edit == 0)
gCanvas = new TCanvas("canv1d","canv1d",1200,900);
else
gCanvas = tempAnalysisCanvas->GetCanvas();
// Fit the PDF with a gaussian
gScan[0]->Fit("gaus","Q");
gScan[0]->GetFunction("gaus")->SetNpx(400);
gStyle->SetOptFit(1);
gScan[1]->Draw("AL");
gPad->Update();
gScan[0]->Draw("LP");
gCanvas->Modified();
gCanvas->Update();
TPaveStats *stats = (TPaveStats*)gScan[0]->FindObject("stats");
if(!cleanPlots)
{
stats->SetX1NDC(0.86); stats->SetX2NDC(1.0);
stats->SetY1NDC(0.87); stats->SetY2NDC(1.0);
}
else
{
stats->SetX1NDC(1.1); stats->SetX2NDC(1.3);
stats->SetY1NDC(1.1); stats->SetY2NDC(1.3);
}
gCanvas->SetGridx(1);
gCanvas->SetGridy(1);
if(axis == 1)
gScan[1]->GetXaxis()->SetTitle("X [#mum]");
else if(axis == 2)
gScan[1]->GetXaxis()->SetTitle("Y [#mum]");
gScan[1]->GetXaxis()->SetTitleOffset(1.3);
gScan[1]->GetXaxis()->CenterTitle(kTRUE);
gScan[1]->GetXaxis()->SetLabelSize(0.027);
gScan[1]->GetXaxis()->SetLabelOffset(0.02);
gScan[1]->GetXaxis()->SetNoExponent(kTRUE);
gScan[1]->GetYaxis()->SetTitle("Normalized ADC integral");
gScan[1]->GetYaxis()->CenterTitle(kTRUE);
gScan[1]->GetYaxis()->SetLabelSize(0.027);
gScan[1]->GetYaxis()->SetLabelOffset(0.02);
gScan[1]->GetYaxis()->SetRangeUser(0,1);
gScan[1]->GetYaxis()->SetTitleOffset(1.4);
gScan[1]->GetYaxis()->SetTitleSize(0.030);
if(!cleanPlots)
{
if(axis == 1)
{
if(posUnitsPlot->widgetCB->GetSelected() == 0)
gScan[1]->SetTitle("SiPM edge detection;X [table units];Normalized ADC integral");
else if(posUnitsPlot->widgetCB->GetSelected() == 1)
gScan[1]->SetTitle("SiPM edge detection;X [#mum];Normalized ADC integral");
}
else if(axis == 2)
{
if(posUnitsPlot->widgetCB->GetSelected() == 0)
gScan[1]->SetTitle("SiPM edge detection;Y [table units];Normalized ADC integral");
else if(posUnitsPlot->widgetCB->GetSelected() == 1)
gScan[1]->SetTitle("SiPM edge detection;Y [#mum];Normalized ADC integral");
}
}
else
{
if(axis == 1)
{
if(posUnitsPlot->widgetCB->GetSelected() == 0)
gScan[1]->SetTitle(";X [table units];Normalized ADC integral");
else if(posUnitsPlot->widgetCB->GetSelected() == 1)
gScan[1]->SetTitle(";X [#mum];Normalized ADC integral");
}
else if(axis == 2)
{
if(posUnitsPlot->widgetCB->GetSelected() == 0)
gScan[1]->SetTitle(";Y [table units];Normalized ADC integral");
else if(posUnitsPlot->widgetCB->GetSelected() == 1)
gScan[1]->SetTitle(";Y [#mum];Normalized ADC integral");
}
}
gScan[1]->SetLineColor(kBlue);
gScan[0]->SetLineWidth(2);
gScan[1]->SetLineWidth(2);
gCanvas->Modified();
gCanvas->Update();
gCanvas->SaveAs(exportname);
// If 2D edge plot, delete the 1D edge plots as we go
if(edge2d)
{
delete gScan[0];
delete gScan[1];
if(edit == 0)
delete gCanvas;
}
else
{
if(edit == 0)
{
delete gScan[0];
delete gScan[1];
delete gCanvas;
}
}
}
void TGAppMainFrame::PhotonMu(TList *files, int edit)
{
unsigned int nrfiles = files->GetSize();
char ctemp[1024];
char exportname[1024];
int j, k = 0, m = 0, n = 0, k2 = 0, m2 = 0;
TCanvas *gCanvas;
TTree *header_data, *meas_data;
double *integralCount, *integralPedestal;
integralCount = new double[nrfiles];
integralPedestal = new double[nrfiles];
double *angle, *pdeval, *muval;
angle = new double[nrfiles];
pdeval = new double[nrfiles];
muval = new double[nrfiles];
// Zero the integral count and accumulated vaules
for(int i = 0; i < (int)nrfiles; i++) {integralCount[i] = 0; integralPedestal[i] = 0; }
// Fitting variables
TSpectrum *spec;
TH1F *histtemp;
TH1 *histback;
TH1F *h2;
float *xpeaks;
TF1 *fit;
TF1 *fittingfunc;
double *fparam, *fparamerr;
double meansel[20];
double sigmasel[20];
double meanparam, paramsigma;
int sortindex[20];
int adcpedestal[2];
int zeromu = 0;
int darkhist = -1;
int nopeaks = -1;
double pointest[12];
bool exclude = false;
// Zero the parameter values
for(int i = 0; i < 20; i++) {meansel[i] = 0; sigmasel[i] = 0; }
float progVal = 0;
analysisProgress->widgetPB->SetPosition(progVal);
gVirtualX->Update(1);
// Check if the checkbox for no peaks is selected - TODO: Still need the situation when we do not have a peaked ADC spectrum
/* if(relPde->widgetChBox[1]->IsDown())
{
printf("PhotonMu(): ADC spectrum has no peak structure.\n");
nopeaks = 1;
// Error if there is no darkhist
if(strcmp("", darkRun->widgetTE->GetText()) == 0)
{
printf("PhotonMu(): Error! The no peak structure option needs a dark histogram.\n");
delete[] integralCount;
delete[] integralPedestal;
delete[] angle;
delete[] pdeval;
delete[] muval;
return;
}
}*/
// Start if we select at least one file
if(nrfiles > 0)
{
// Find the pedestal peak for the dark histogram, and use it for all if there are no peaks
if(nopeaks != -1)
{
// Replot the spectrum on analysisCanvas and do not close the input file
DisplayHistogram( (char*)(darkRun->widgetTE->GetText()), 0, 1);
analysisCanvas->GetCanvas()->Modified();
analysisCanvas->GetCanvas()->Update();
// Get the spectrum
histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname);
npeaks = 15;
double par[300];
spec = new TSpectrum(npeaks);
// Find spectrum background
histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same");
// Clone histogram and subtract background from it if we select that option
h2 = (TH1F*)histtemp->Clone("h2");
if(fitChecks->widgetChBox[0]->IsDown())
h2->Add(histback, -1);
// Search for the peaks
int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() );
printf("PhotonMu(): Found %d candidates to fit.\n",found);
npeaks = found;
// Set initial peak parameters
xpeaks = spec->GetPositionX();
for(j = 0; j < found; j++)
{
float xp = xpeaks[j];
int bin = h2->GetXaxis()->FindBin(xp);
float yp = h2->GetBinContent(bin);
par[3*j] = yp;
par[3*j+1] = xp;
par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber();
}
// Fit the histogram
fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks);
TVirtualFitter::Fitter(histtemp, 3*npeaks);
fit->SetParameters(par);
fit->SetNpx(300);
h2->Fit("fit","Q");
// Get the fitted parameters
fittingfunc = h2->GetFunction("fit");
fparam = fittingfunc->GetParameters();
fparamerr = fittingfunc->GetParErrors();
// Gather the parameters (mean peak value for now)
int j = 1;
int nrfit = 0;
while(1)
{
if( (fparam[j] < 1.E-30) || (nrfit > 8) )
break;
else
{
// Check if pedestal is above the lower limit and sigma is smaller than the mean
if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) )
{
// 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
meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber());
sigmasel[nrfit] = fparam[j+1];
nrfit++;
}
}
j+=3;
}
TMath::Sort(nrfit, meansel, sortindex, kFALSE);
meanparam = meansel[sortindex[0]];
paramsigma = sigmasel[sortindex[0]];
for(j = 0; j < nrfit; j++)
if(DBGSIG)
printf("PhotonMu(): %d: peak mean = %lf\n", j, meansel[sortindex[j]]);
j = 0;
adcpedestal[0] = 0;
adcpedestal[1] = -1;
while(1)
{
int bin = histtemp->GetXaxis()->FindBin((int)(j+meanparam+paramsigma));
int yp = histtemp->GetBinContent(bin);
// 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)
if(adcpedestal[1] == -1)
{
adcpedestal[0] = j+meanparam+paramsigma;
adcpedestal[1] = yp;
}
else
{
if( (npeaks > 1) && (adcpedestal[1] >= yp) )
{
adcpedestal[0] = j+meanparam+paramsigma;
adcpedestal[1] = yp;
}
else if( (npeaks == 1) && (adcpedestal[0] < meanparam+5*paramsigma) ) // TODO -> Determining the pedestal when only one peak
{
adcpedestal[0] = j+meanparam+paramsigma;
adcpedestal[1] = yp;
}
else
break;
}
j++;
if(j > 50) break;
}
if(midPeak->widgetChBox[0]->IsDown())
{
if( (meanparam - (int)meanparam >= 0.) && (meanparam - (int)meanparam < 0.5) )
m = TMath::Floor(meanparam);
else if( (meanparam - (int)meanparam >= 0.5) && (meanparam - (int)meanparam < 1.) )
m = TMath::Ceil(meanparam);
int bin = histtemp->GetXaxis()->FindBin(m);
adcpedestal[0] = m;
printf("midpeak x = %d, ", adcpedestal[0]);
adcpedestal[1] = histtemp->GetBinContent(bin);
}
// Option to show the fit
fittingfunc->Draw("L SAME");
analysisCanvas->GetCanvas()->Modified();
analysisCanvas->GetCanvas()->Update();
printf("Pedestal ends = %d and nr. of counts = %d\n", adcpedestal[0], adcpedestal[1]);
// Delete the opened histogram and spectrum
delete spec;
delete inroot;
// return;
}
printf("PhotonMu(): Continuing with the rest of the spectra.\n");
// Check all histograms for pedestal peak values
for(int i = 0; i < (int)nrfiles; i++)
{
if( (nrfiles == 1) || (!multiSelect->widgetChBox[0]->IsDown()) )
{
printf("PhotonMu(): Only one file selected. Not running analysis, just showing the fit.\n");
// Replot the spectrum on analysisCanvas and do not close the input file
DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1);
analysisCanvas->GetCanvas()->Modified();
analysisCanvas->GetCanvas()->Update();
// Get the spectrum
histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname);
npeaks = 15;
double par[300];
spec = new TSpectrum(npeaks);
// Find spectrum background
histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same");
// Clone histogram and subtract background from it if we select that option
h2 = (TH1F*)histtemp->Clone("h2");
if(fitChecks->widgetChBox[0]->IsDown())
h2->Add(histback, -1);
// Search for the peaks
int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() );
printf("PhotonMu(): Found %d candidates to fit.\n",found);
npeaks = found;
// Set initial peak parameters
xpeaks = spec->GetPositionX();
for(j = 0; j < found; j++)
{
float xp = xpeaks[j];
int bin = h2->GetXaxis()->FindBin(xp);
float yp = h2->GetBinContent(bin);
par[3*j] = yp;
par[3*j+1] = xp;
par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber();
}
// Fit the histogram
fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks);
TVirtualFitter::Fitter(histtemp, 3*npeaks);
fit->SetParameters(par);
fit->SetNpx(300);
h2->Fit("fit","Q");
// Get the fitted parameters
fittingfunc = h2->GetFunction("fit");
fparam = fittingfunc->GetParameters();
fparamerr = fittingfunc->GetParErrors();
// Gather the parameters (mean peak value for now)
int j = 1;
int nrfit = 0;
while(1)
{
if( (fparam[j] < 1.E-30) || (nrfit > 8) )
break;
else
{
// Check if pedestal is above the lower limit and sigma is smaller than the mean
if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) )
{
// 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
meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber());
sigmasel[nrfit] = fparam[j+1];
nrfit++;
}
}
j+=3;
}
TMath::Sort(nrfit, meansel, sortindex, kFALSE);
fittingfunc->Draw("SAME");
analysisCanvas->GetCanvas()->Modified();
analysisCanvas->GetCanvas()->Update();
meanparam = meansel[sortindex[0]];
paramsigma = sigmasel[sortindex[0]];
for(j = 0; j < nrfit; j++)
printf("PhotonMu(): %d: peak mean = %lf\n", j, meansel[sortindex[j]]);
return;
}
if(files->At(i))
{
if(strcmp(files->At(i)->GetTitle(),darkRun->widgetTE->GetText()) == 0)
{
printf("PhotonMu(): %s is the dark histogram file.\n", files->At(i)->GetTitle());
darkhist = i;
}
if(nopeaks == -1)
{
// Replot the spectrum on analysisCanvas and do not close the input file
DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1);
analysisCanvas->GetCanvas()->Modified();
analysisCanvas->GetCanvas()->Update();
// Get the spectrum
histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname);
npeaks = 15;
double par[300];
spec = new TSpectrum(npeaks);
// Find spectrum background
histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same");
// Clone histogram and subtract background from it if we select that option
h2 = (TH1F*)histtemp->Clone("h2");
if(fitChecks->widgetChBox[0]->IsDown())
h2->Add(histback, -1);
// Search for the peaks
int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() );
printf("PhotonMu(): Found %d candidates to fit.\n",found);
npeaks = found;
// Set initial peak parameters
xpeaks = spec->GetPositionX();
for(j = 0; j < found; j++)
{
float xp = xpeaks[j];
int bin = h2->GetXaxis()->FindBin(xp);
float yp = h2->GetBinContent(bin);
par[3*j] = yp;
par[3*j+1] = xp;
par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber();
}
// Fit the histogram
fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks);
TVirtualFitter::Fitter(histtemp, 3*npeaks);
fit->SetParameters(par);
fit->SetNpx(300);
h2->Fit("fit","Q");
// Get the fitted parameters
fittingfunc = h2->GetFunction("fit");
fparam = fittingfunc->GetParameters();
fparamerr = fittingfunc->GetParErrors();
// Gather the parameters (mean peak value for now)
int j = 1;
int nrfit = 0;
while(1)
{
if( (fparam[j] < 1.E-30) || (nrfit > 8) )
break;
else
{
// Check if pedestal is above the lower limit and sigma is smaller than the mean
if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) )
{
// 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
meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber());
sigmasel[nrfit] = fparam[j+1];
nrfit++;
}
}
j+=3;
}
TMath::Sort(nrfit, meansel, sortindex, kFALSE);
meanparam = meansel[sortindex[0]];
paramsigma = sigmasel[sortindex[0]];
for(j = 0; j < nrfit; j++)
if(DBGSIG)
printf("PhotonMu(): %d: peak mean = %lf\n", j, meansel[sortindex[j]]);
j = 0;
adcpedestal[0] = 0;
adcpedestal[1] = -1;
while(1)
{
int bin = histtemp->GetXaxis()->FindBin((int)(j+meanparam+paramsigma));
int yp = histtemp->GetBinContent(bin);
// 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)
if(adcpedestal[1] == -1)
{
adcpedestal[0] = j+meanparam+paramsigma;
adcpedestal[1] = yp;
}
else
{
if( (npeaks > 1) && (adcpedestal[1] >= yp) )
{
adcpedestal[0] = j+meanparam+paramsigma;
adcpedestal[1] = yp;
}
else if( (npeaks == 1) && (adcpedestal[0] < meanparam+5*paramsigma) ) // TODO -> Determining the pedestal when only one peak
{
adcpedestal[0] = j+meanparam+paramsigma;
adcpedestal[1] = yp;
}
else
break;
}
j++;
if(j > 50) break;
}
if(npeaks > 1)
{
int bin = histtemp->GetXaxis()->FindBin((int)(meanparam+meansel[sortindex[1]])/2);
adcpedestal[0] = (meanparam+meansel[sortindex[1]])/2;
printf("PhotonMu(): multipeak x = %d, ", adcpedestal[0]);
adcpedestal[1] = histtemp->GetBinContent(bin);
}
if(midPeak->widgetChBox[0]->IsDown())
{
if( (meanparam - (int)meanparam >= 0.) && (meanparam - (int)meanparam < 0.5) )
m = TMath::Floor(meanparam);
else if( (meanparam - (int)meanparam >= 0.5) && (meanparam - (int)meanparam < 1.) )
m = TMath::Ceil(meanparam);
int bin = histtemp->GetXaxis()->FindBin(m);
adcpedestal[0] = m;
printf("midpeak x = %d, ", adcpedestal[0]);
adcpedestal[1] = histtemp->GetBinContent(bin);
}
/* // Option to show the fit
fittingfunc->Draw("L SAME");
analysisCanvas->GetCanvas()->Modified();
analysisCanvas->GetCanvas()->Update();*/
printf("Pedestal ends = %d and nr. of counts = %d\n", adcpedestal[0], adcpedestal[1]);
// Delete the opened histogram and spectrum
delete spec;
delete inroot;
}
// Open the input file and read header, ADC and TDC values
sprintf(ctemp, "%s", files->At(i)->GetTitle());
inroot = new TFile(ctemp, "READ");
inroot->GetObject("header_data", header_data);
inroot->GetObject("meas_data", meas_data);
// Reading the header
if( header_data->FindBranch("angle") )
{
header_data->SetBranchAddress("angle", &evtheader.angle);
header_data->GetEntry(0);
}
else
{
printf("PhotonMu(): Error! Selected file has no angle header value. Please edit header to add the angle header value.\n");
break;
}
char rdc[256];
j = selectCh->widgetNE[0]->GetNumber();
double rangetdc[2];
rangetdc[0] = tdcRange->widgetNE[0]->GetNumber();
rangetdc[1] = tdcRange->widgetNE[1]->GetNumber();
k = 0;
k2 = 0;
m = 0;
m2 = 0;
// Reading the data
for(int e = 0; e < meas_data->GetEntries(); e++)
{
sprintf(rdc, "ADC%d", j);
meas_data->SetBranchAddress(rdc, &evtdata.adcdata[j]);
meas_data->GetEntry(e);
sprintf(rdc, "TDC%d", j);
meas_data->SetBranchAddress(rdc, &evtdata.tdcdata[j]);
meas_data->GetEntry(e);
// If our data point is inside the TDC window
if( ((double)evtdata.tdcdata[j]/tdctimeconversion >= rangetdc[0]) && ((double)evtdata.tdcdata[j]/tdctimeconversion <= rangetdc[1]) )
{
// Gather only the integral of the pedestal
if((double)evtdata.adcdata[j] < (double)adcpedestal[0]+0.5 )
{
k2++;
m2 += evtdata.adcdata[j];
}
// Gather the complete integral
k++;
m += evtdata.adcdata[j];
}
}
// Determine the angle, mu and relative PDE
angle[i] = (double)(evtheader.angle);
integralCount[i] += (double)m;
printf("PhotonMu(): %lf: Complete integral (%d evts) = %lf\n", angle[i], k, integralCount[i]);
integralPedestal[i] += (double)m2;
printf("PhotonMu(): %lf: Pedestal integral (%d evts) = %lf\n", angle[i], k2, integralPedestal[i]);
if( (angle[i] == zeroAngle->widgetNE[0]->GetNumber()) && (darkhist != i) )
zeromu = i;
muval[i] = -TMath::Log((double)k2/(double)k);
printf("PhotonMu(): %lf: muval = %lf\n", angle[i], muval[i]);
inroot->Close();
delete inroot;
}
// Update the progress bar
progVal = (float)(90.00/nrfiles)*i;
analysisProgress->widgetPB->SetPosition(progVal);
gVirtualX->Update(1);
}
printf("PhotonMu(): %d files were selected.\n", nrfiles);
printf("PhotonMu(): angle\tmu\trelative PDE\n");
m = 0;
// Set the 0 degree muval, reuse meansel[1]
meansel[1] = muval[zeromu];
printf("Zero value (id=%d, angle=%lf) = %lf\n", zeromu, angle[zeromu], meansel[1]);
// TODO - point estimation still not working correctly!
for(int i = 0; i < (int)nrfiles; i++)
{
// Estimate next point and check error (5 point least square fit estimation)
if( ((i > 4) && (m == i)) || ((i > 5) && (m < i)) )
{
// Set exclude signal to false
exclude = false;
// Get next point values (if zero value -> need to add the dark hist value again)
pointest[10] = angle[i];
pointest[11] = muval[i];
// printf("Last point = [%lf,%lf]\n", pointest[10], pointest[11]);
// Check if next point has larger error than acceptable (if yes, set exclude signal to true), reuse meansel[0]
meansel[0] = PointEstimate(5, pointest); // PointEstimate only works with very small step size
if(meansel[0] > accError->widgetNE[0]->GetNumber())
{
printf("PhotonMu(): Point (%lf, %lf) excluded due to error: %lf\n", pointest[10], pointest[11], meansel[0]);
exclude = true;
}
if(isinf(meansel[0]))
{
printf("PhotonMu(): Point (%lf, %lf) excluded due to being infinite: %lf\n", pointest[10], pointest[11], meansel[0]);
exclude = true;
}
// Value with 0 angle and dark histogram are always needed, so should not be excluded
if(i == darkhist)
exclude = false;
// If nothing excluded, pass the points in pointest variable like in a FIFO
if(!exclude)
{
for(int j = 0; j < 10; j++)
{
if(DBGSIG) printf("PhotonMu(): j = %d: Old X value = %lf\n", j, pointest[j]);
pointest[j] = pointest[j+2];
}
}
else
{
for(int j = 0; j < 10; j++)
if(DBGSIG) printf("PhotonMu(): No j = %d: Old X value = %lf\n", j, pointest[j]);
}
}
else
{
// First 5 points act as estimator points for next one
pointest[2*m] = angle[i];
pointest[2*m+1] = muval[i];
}
printf("Exclude signal = %d\n", (int)exclude);
// Run only if we have a dark run histogram and middle pedestal peak estimation
if( (darkhist != -1) && midPeak->widgetChBox[0]->IsDown() )
{
if(DBGSIG) printf("PhotonMu(): m = %d, i = %d: muval = %lf, ", m, i, muval[i]);
// Subtract the dark value from all values
angle[m] = angle[i];
muval[m] = muval[i] - muval[darkhist];
if(DBGSIG) printf("angle = %lf, newmuval = %lf, darkmuval = %lf, ", angle[m], muval[m], muval[darkhist]);
// Calculate relative PDE
// pdeval[m] = muval[m]/(muval[zeromu]*TMath::Cos(angle[m]*TMath::ACos(-1.)/180.));
pdeval[m] = muval[m]/((meansel[1]-muval[darkhist])*TMath::Cos(angle[m]*TMath::ACos(-1.)/180.));
if(DBGSIG) printf("pdeval = %lf\n", pdeval[m]);
// Only increase counter if error is not too big
if( (darkhist != i) && (!exclude) )
m++;
}
else
{
// Relative PDE calculation
angle[m] = angle[i];
muval[m] = muval[i];
// pdeval[m] = muval[m]/(muval[zeromu]*TMath::Cos(angle[m]*TMath::ACos(-1.)/180.));
pdeval[m] = muval[m]/(meansel[1]*TMath::Cos(angle[m]*TMath::ACos(-1.)/180.));
// Only increase counter if error is not too big
if(!exclude)
m++;
}
printf("PhotonMu(): %lf\t%lf\t%lf\n", angle[i], muval[i], pdeval[i]);
}
// Check for range of values to plot
double plotMax = 0.;
for(int i = 0; i < (int)nrfiles; i++)
{
plotMax = TMath::Max(plotMax, muval[i]);
plotMax = TMath::Max(plotMax, pdeval[i]);
}
if(plotMax <= 0.)
plotMax = 1.1;
printf("PhotonMu(): Maximum value: %lf\n", plotMax);
if(DBGSIG) printf("\n");
if(darkhist != -1)
printf("PhotonMu(): Number of excluded points: %d\n", (nrfiles-1-m));
else
printf("PhotonMu(): Number of excluded points: %d\n", (nrfiles-m));
// Plot mu and PDE angle dependance plots
if(edit == 0)
gCanvas = new TCanvas("canv","canv",1200,900);
else if(edit == 1)
gCanvas = tempAnalysisCanvas->GetCanvas();
gCanvas->cd();
gCanvas->SetGrid();
TGraph *pde = new TGraph(m, angle, pdeval);
pde->SetMarkerStyle(21);
pde->SetMarkerSize(0.7);
pde->SetMarkerColor(2);
pde->SetLineWidth(1);
pde->SetLineColor(2);
pde->GetXaxis()->SetLabelSize(0.030);
pde->GetXaxis()->CenterTitle();
// pde->GetXaxis()->SetRange(angle[0],angle[nrfiles-1]);
// pde->GetXaxis()->SetRangeUser(angle[0],angle[nrfiles-1]);
pde->GetXaxis()->SetRange(-90.0,90.0);
pde->GetXaxis()->SetRangeUser(-90.0,90.0);
pde->GetXaxis()->SetLimits(-90.0,90.0);
pde->GetYaxis()->SetTitleOffset(1.2);
pde->GetYaxis()->SetLabelSize(0.030);
pde->GetYaxis()->CenterTitle();
pde->GetYaxis()->SetRange(0., 1.1*plotMax);
pde->GetYaxis()->SetRangeUser(0., 1.1*plotMax);
pde->GetYaxis()->SetLimits(0., 1.1*plotMax);
pde->SetName("pde");
pde->Draw("ALP");
pde->SetTitle(";Incidence angle (#circ);Relative PDE(#theta) / #mu(#theta)");
TGraph *mugr = new TGraph(m, angle, muval);
mugr->SetMarkerStyle(20);
mugr->SetMarkerSize(0.7);
mugr->SetMarkerColor(4);
mugr->SetLineWidth(1);
mugr->SetLineColor(4);
mugr->SetName("muval");
mugr->Draw("SAME;LP");
gCanvas->Modified();
gCanvas->Update();
if(edit == 0)
{
remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp);
sprintf(exportname, "%s_pde.pdf", ctemp);
gCanvas->SaveAs(exportname);
delete mugr;
delete pde;
delete gCanvas;
}
else if(edit == 1)
{
gCanvas->Modified();
gCanvas->Update();
}
// Update the progress bar
analysisProgress->widgetPB->SetPosition(100.);
gVirtualX->Update(1);
}
}
void TGAppMainFrame::BreakdownVolt(TList *files, int edit)
{
unsigned int nrfiles = files->GetSize();
char ctemp[1024];
char exportname[1024];
char paramname[1024];
int j, k = 0;
TCanvas *gCanvas;
TTree *header_data, *meas_data;
// Fitting variables
TSpectrum *spec;
TH1F *histtemp;
TH1 *histback;
TH1F *h2;
float *xpeaks;
TF1 *fit;
TF1 *fittingfunc;
TLatex *latex;
double *fparam, *fparamerr;
double meansel[20];
double meanselerr[20];
double sigmasel[20];
double meanparam, meanparamerr, paramsigma;
int sortindex[20];
bool exclude = false;
int p = 0;
double volt[nrfiles], volterr[nrfiles], sep[3][nrfiles], seperr[3][nrfiles];
int first = 1;
FILE *fp;
remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp);
sprintf(paramname, "%s_fitresult.txt", ctemp);
fp = fopen(paramname, "w");
fclose(fp);
int peaklimit = minPeak->widgetNE[0]->GetNumber()+1; // +1 to account for the pedestal peak
// Zero the parameter values
for(int i = 0; i < nrfiles; i++)
{
volt[i] = 0; volterr[i] = 0;
sep[0][i] = 0; sep[1][i] = 0; sep[2][i] = 0;
seperr[0][i] = 0; seperr[1][i] = 0; seperr[2][i] = 0;
if(i < 20) { meansel[i] = 0; meanselerr[i] = 0; sigmasel[i] = 0; }
}
float progVal = 0;
analysisProgress->widgetPB->SetPosition(progVal);
gVirtualX->Update(1);
// Start if we select at least one file
if(nrfiles > 0)
{
for(int i = 0; i < (int)nrfiles; i++)
{
if( (nrfiles == 1) || (!multiSelect->widgetChBox[0]->IsDown()) )
{
printf("BreakdownVolt(): Only one file selected. Not running analysis, just showing the fit.\n");
// Replot the spectrum on analysisCanvas and do not close the input file
DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1);
analysisCanvas->GetCanvas()->Modified();
analysisCanvas->GetCanvas()->Update();
// Get the spectrum
histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname);
npeaks = 15;
double par[300];
spec = new TSpectrum(npeaks);
// Find spectrum background
histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same");
// Clone histogram and subtract background from it if we select that option
h2 = (TH1F*)histtemp->Clone("h2");
if(fitChecks->widgetChBox[0]->IsDown())
h2->Add(histback, -1);
// Search for the peaks
int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() );
printf("PhotonMu(): Found %d candidates to fit.\n",found);
npeaks = found;
// Set initial peak parameters
xpeaks = spec->GetPositionX();
for(j = 0; j < found; j++)
{
float xp = xpeaks[j];
int bin = h2->GetXaxis()->FindBin(xp);
float yp = h2->GetBinContent(bin);
par[3*j] = yp;
par[3*j+1] = xp;
par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber();
}
// Fit the histogram
fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks);
TVirtualFitter::Fitter(histtemp, 3*npeaks);
fit->SetParameters(par);
fit->SetNpx(300);
h2->Fit("fit","Q");
// Get the fitted parameters
fittingfunc = h2->GetFunction("fit");
fparam = fittingfunc->GetParameters();
fparamerr = fittingfunc->GetParErrors();
// Gather the parameters (mean peak value for now)
int j = 1;
int nrfit = 0;
while(1)
{
if( (fparam[j] < 1.E-30) || (nrfit > 8) )
break;
else
{
// Check if pedestal is above the lower limit and sigma is smaller than the mean
if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) )
{
// 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
meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber());
meanselerr[nrfit] = fparamerr[j];
sigmasel[nrfit] = fparam[j+1];
nrfit++;
}
}
j+=3;
}
TMath::Sort(nrfit, meansel, sortindex, kFALSE);
fittingfunc->Draw("SAME");
analysisCanvas->GetCanvas()->Modified();
analysisCanvas->GetCanvas()->Update();
meanparam = meansel[sortindex[0]];
meanparamerr = meanselerr[sortindex[0]];
paramsigma = sigmasel[sortindex[0]];
for(j = 0; j < nrfit; j++)
printf("BreakdownVolt(): %d: peak mean = %lf, peak err = %lf\n", j, meansel[sortindex[j]], meanselerr[sortindex[j]]);
return;
}
if(files->At(i))
{
// Replot the spectrum on analysisCanvas and do not close the input file
DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1);
analysisCanvas->GetCanvas()->Modified();
analysisCanvas->GetCanvas()->Update();
// Get the spectrum
histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname);
npeaks = 15;
double par[300];
spec = new TSpectrum(npeaks);
// Find spectrum background
histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same");
// Clone histogram and subtract background from it if we select that option
h2 = (TH1F*)histtemp->Clone("h2");
if(fitChecks->widgetChBox[0]->IsDown())
h2->Add(histback, -1);
// Search for the peaks
int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() );
printf("BreakdownVolt(): Found %d candidates to fit.\n",found);
npeaks = found;
// Set initial peak parameters
xpeaks = spec->GetPositionX();
for(j = 0; j < found; j++)
{
float xp = xpeaks[j];
int bin = h2->GetXaxis()->FindBin(xp);
float yp = h2->GetBinContent(bin);
par[3*j] = yp;
par[3*j+1] = xp;
par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber();
}
// Fit the histogram
fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks);
TVirtualFitter::Fitter(histtemp, 3*npeaks);
fit->SetParameters(par);
fit->SetNpx(300);
h2->Fit("fit","Q");
// Get the fitted parameters
fittingfunc = h2->GetFunction("fit");
if(nrfiles == 1) // TODO: Show the fit if only one file is selected
fittingfunc->Draw("SAME");
fparam = fittingfunc->GetParameters();
fparamerr = fittingfunc->GetParErrors();
// Gather the parameters (mean peak value for now)
int j = 1;
int nrfit = 0;
while(1)
{
if( (fparam[j] < 1.E-30) || (fparamerr[j] < 1.E-10) )// TODO: Maybe not correct for the error
break;
else
{
// Check if pedestal is above the lower limit and sigma is smaller than the mean
if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) )
{
meansel[nrfit] = fparam[j];
meanselerr[nrfit] = fparamerr[j];
sigmasel[nrfit] = fparam[j+1];
nrfit++;
}
}
j+=3;
}
TMath::Sort(nrfit, meansel, sortindex, kFALSE);
meanparam = meansel[sortindex[0]];
meanparamerr = meanselerr[sortindex[0]];
paramsigma = sigmasel[sortindex[0]];
for(j = 0; j < nrfit; j++)
{
if(DBGSIG)
printf("BreakdownVolt(): %d: peak mean = %lf, peak err = %lf\n", j, meansel[sortindex[j]], meanselerr[sortindex[j]]);
}
// Delete the opened histogram and spectrum
delete spec;
delete inroot;
// Open the input file and read header, ADC and TDC values
sprintf(ctemp, "%s", files->At(i)->GetTitle());
inroot = new TFile(ctemp, "READ");
inroot->GetObject("header_data", header_data);
inroot->GetObject("meas_data", meas_data);
// Reading the header
if( header_data->FindBranch("biasvolt") )
{
header_data->SetBranchAddress("biasvolt", &evtheader.biasvolt);
header_data->GetEntry(0);
}
else
{
printf("BreakdownVolt(): Error! Selected file has no bias voltage header value. Please edit header to add the bias voltage header value.\n");
break;
}
// h2->SetStats(0);
analysisCanvas->GetCanvas()->Modified();
analysisCanvas->GetCanvas()->Update();
// Save each fitting plot
if(fitChecks->widgetChBox[1]->IsDown()) // TODO: Check if this works
{
remove_ext((char*)files->At(i)->GetTitle(), ctemp);
sprintf(exportname, "%s_fit.pdf", ctemp);
analysisCanvas->GetCanvas()->SaveAs(exportname);
}
// Calculate the separation between peaks
volt[p] = evtheader.biasvolt;
volterr[p] = 1.e-5;
if(nrfit == 3)
{
sep[0][p] = meansel[sortindex[2]] - meansel[sortindex[1]];
seperr[0][p] = TMath::Abs(meanselerr[sortindex[2]]) + TMath::Abs(meanselerr[sortindex[1]]);
exclude = (seperr[0][p]/sep[0][p] < accError->widgetNE[0]->GetNumber());
}
else if(nrfit == 4)
{
sep[0][p] = meansel[sortindex[2]] - meansel[sortindex[1]];
sep[1][p] = meansel[sortindex[3]] - meansel[sortindex[2]];
seperr[0][p] = TMath::Abs(meanselerr[sortindex[2]]) + TMath::Abs(meanselerr[sortindex[1]]);
seperr[1][p] = TMath::Abs(meanselerr[sortindex[3]]) + TMath::Abs(meanselerr[sortindex[2]]);
exclude = ((seperr[0][p]/sep[0][p] < accError->widgetNE[0]->GetNumber()) && (seperr[1][p]/sep[1][p] < accError->widgetNE[0]->GetNumber()));
}
else if(nrfit > 4)
{
sep[0][p] = meansel[sortindex[2]] - meansel[sortindex[1]];
sep[1][p] = meansel[sortindex[3]] - meansel[sortindex[2]];
sep[2][p] = meansel[sortindex[4]] - meansel[sortindex[3]];
seperr[0][p] = TMath::Abs(meanselerr[sortindex[2]]) + TMath::Abs(meanselerr[sortindex[1]]);
seperr[1][p] = TMath::Abs(meanselerr[sortindex[3]]) + TMath::Abs(meanselerr[sortindex[2]]);
seperr[2][p] = TMath::Abs(meanselerr[sortindex[4]]) + TMath::Abs(meanselerr[sortindex[3]]);
exclude = ((seperr[0][p]/sep[0][p] < accError->widgetNE[0]->GetNumber()) && (seperr[1][p]/sep[1][p] < accError->widgetNE[0]->GetNumber()) && (seperr[2][p]/sep[2][p] < accError->widgetNE[0]->GetNumber()));
}
else
{
printf("BreakdownVolt(): The current separation measurements does not fall within acceptable errors.\n");
exclude = false;
}
printf("BreakdownVolt(): Calculated separation between peaks %d and %d: sep = %lf, seperr = %lf\n", (int)peakSepCalc->widgetNE[0]->GetNumber()-1, (int)peakSepCalc->widgetNE[0]->GetNumber(), sep[(int)peakSepCalc->widgetNE[0]->GetNumber()-1][p], seperr[(int)peakSepCalc->widgetNE[0]->GetNumber()-1][p]);
// Write out parameters to a file
fp = fopen(paramname, "a");
if(exclude)
{
if(first == 1)
{
fprintf(fp, "%le\t%d\t", evtheader.biasvolt, nrfit);
for(j = 0; j < nrfit; j++)
fprintf(fp, "%le\t%le\t", meansel[sortindex[j]], meanselerr[sortindex[j]]);
fprintf(fp,"\n");
first = 0;
}
p++;
}
else
{
if(nrfit == 3)
printf("Point (at %.2lfV) rejected due to too large errors: %lf\n", volt[p], seperr[0][p]/sep[0][p]);
else if(nrfit == 4)
printf("Point (at %.2lfV) rejected due to too large errors: %lf, %lf\n", volt[p], seperr[0][p]/sep[0][p], seperr[1][p]/sep[1][p]);
else if(nrfit > 4)
printf("Point (at %.2lfV) rejected due to too large errors: %lf, %lf, %lf\n", volt[p], seperr[0][p]/sep[0][p], seperr[1][p]/sep[1][p], seperr[2][p]/sep[2][p]);
}
fclose(fp);
}
if(nrfiles == 1) break;
first = 1;
// Update the progress bar
progVal = (float)(90.00/nrfiles)*i;
analysisProgress->widgetPB->SetPosition(progVal);
gVirtualX->Update(1);
}
printf("BreakdownVolt(): %d files were selected.\n", nrfiles);
printf("BreakdownVolt(): Number of points to plot is %d\n", p);
// Plot and fit breakdown voltage plot
if(edit == 0)
gCanvas = new TCanvas("canv","canv",1200,900);
else if(edit == 1)
gCanvas = tempAnalysisCanvas->GetCanvas();
gCanvas->cd();
gCanvas->SetGrid();
TGraphErrors* bdplot;
k = peakSepCalc->widgetNE[0]->GetNumber();
if(k < 4)
bdplot = new TGraphErrors(p, volt, sep[k-1], volterr, seperr[k-1]);
else
{
printf("BreakdownVolt(): Unsupported peak separation selected (%d).\n", k);
return;
}
bdplot->SetMarkerStyle(20);
bdplot->SetMarkerSize(0.4);
bdplot->SetMarkerColor(kBlue);
bdplot->SetLineWidth(1);
bdplot->SetLineColor(kBlue);
bdplot->GetXaxis()->SetLabelSize(0.030);
bdplot->GetXaxis()->CenterTitle();
// bdplot->GetXaxis()->SetRange(-90,90);
// bdplot->GetXaxis()->SetRangeUser(-90,90);
// bdplot->GetXaxis()->SetLimits(-90,90);
bdplot->GetYaxis()->SetTitleOffset(1.2);
bdplot->GetYaxis()->SetLabelSize(0.030);
bdplot->GetYaxis()->CenterTitle();
// bdplot->GetYaxis()->SetRangeUser(0., 1.2);
bdplot->SetName("bdplot");
bdplot->Draw("AP");
bdplot->SetTitle(";Bias voltage (V);Peak separation");
// Fit the breakdown voltage plot with a line
bdplot->Fit("pol1","Q");
TF1 *fit1 = bdplot->GetFunction("pol1");
meansel[0] = fit1->GetParameter(0);
meanselerr[0] = fit1->GetParError(0);
meansel[1] = fit1->GetParameter(1);
meanselerr[1] = fit1->GetParError(1);
meansel[2] = -meansel[0]/meansel[1];
if(!cleanPlots)
{
sprintf(ctemp, "#splitline{#Delta_{p}(U) = (%.2lf #pm %.2lf)#timesU + (%.2lf #pm %.3lf)}{U_{0} = %.2lf #pm %.3lf}", meansel[0], meanselerr[0], meansel[1], meanselerr[1], meansel[2], meansel[2]*(TMath::Abs(meanselerr[0]/meansel[0]) + TMath::Abs(meanselerr[1]/meansel[1])) );
latex = new TLatex();
latex->SetTextSize(0.039);
latex->DrawLatex(volt[0], 0.97*sep[0][sortindex[1]], ctemp);
printf("#Delta_{p}(U) = (%.6lf #pm %.8lf)#timesU + (%.6lf #pm %.8lf)}{U_{0} = %.6lf #pm %.8lf\n", meansel[0], meanselerr[0], meansel[1], meanselerr[1], meansel[2], meansel[2]*(TMath::Abs(meanselerr[0]/meansel[0]) + TMath::Abs(meanselerr[1]/meansel[1])) );
}
else
printf("#Delta_{p}(U) = (%.2lf #pm %.2lf)#timesU + (%.2lf #pm %.3lf)}{U_{0} = %.2lf #pm %.3lf\n", meansel[0], meanselerr[0], meansel[1], meanselerr[1], meansel[2], meansel[2]*(TMath::Abs(meanselerr[0]/meansel[0]) + TMath::Abs(meanselerr[1]/meansel[1])) );
if(edit == 0)
{
remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp);
sprintf(exportname, "%s_breakdown.pdf", ctemp);
gCanvas->SaveAs(exportname);
delete bdplot;
delete gCanvas;
}
else if(edit == 1)
{
gCanvas->Modified();
gCanvas->Update();
}
// Update the progress bar
analysisProgress->widgetPB->SetPosition(100.);
gVirtualX->Update(1);
}
}
void TGAppMainFrame::SurfaceScan(TList *files, int edit)
{
unsigned int nrfiles = files->GetSize();
char ctemp[1024];
char exportname[1024];
int j, k = 0, m = 0, n = 0;
TTree *header_data, *meas_data;
double *integralCount;
integralCount = new double[nrfiles];
double *surfx, *surfy;
surfx = new double[nrfiles];
surfy = new double[nrfiles];
double xsurfmin = 0, ysurfmin = 0;
int nrentries;
double minInteg, maxInteg;
bool norm = surfScanOpt->widgetChBox[0]->IsDown();
double curyval;
// bool edge2d = false;
TCanvas *gCanvas;
float progVal = 0;
analysisProgress->widgetPB->SetPosition(progVal);
gVirtualX->Update(1);
// Zero the integral count and accumulated vaules
for(int i = 0; i < (int)nrfiles; i++) integralCount[i] = 0;
// Start if we select at more than one file
if( (nrfiles > 1) && (multiSelect->widgetChBox[0]->IsOn()) )
{
printf("SurfaceScan(): Creating a surface plot. Please wait...\n");
for(int i = 0; i < (int)nrfiles; i++)
{
if(files->At(i))
{
n++;
// Read the X and Y positions from header and ADC and TDC values from the measurements
sprintf(ctemp, "%s", files->At(i)->GetTitle());
inroot = new TFile(ctemp, "READ");
inroot->GetObject("header_data", header_data);
inroot->GetObject("meas_data", meas_data);
header_data->SetBranchAddress("xpos", &evtheader.xpos);
header_data->GetEntry(0);
header_data->SetBranchAddress("ypos", &evtheader.ypos);
header_data->GetEntry(0);
char rdc[256];
j = selectCh->widgetNE[0]->GetNumber();
double rangetdc[2];
rangetdc[0] = tdcRange->widgetNE[0]->GetNumber();
rangetdc[1] = tdcRange->widgetNE[1]->GetNumber();
k = 0;
m = 0;
// Reading the data
for(int e = 0; e < meas_data->GetEntries(); e++)
{
sprintf(rdc, "ADC%d", j);
meas_data->SetBranchAddress(rdc, &evtdata.adcdata[j]);
meas_data->GetEntry(e);
sprintf(rdc, "TDC%d", j);
meas_data->SetBranchAddress(rdc, &evtdata.tdcdata[j]);
meas_data->GetEntry(e);
// Use data point only if it is inside the TDC window
if( ((double)evtdata.tdcdata[j]/tdctimeconversion >= rangetdc[0]) && ((double)evtdata.tdcdata[j]/tdctimeconversion <= rangetdc[1]) )
{
k++;
m += evtdata.adcdata[j];
}
}
// X, Y and Z values from each file (table units or microns)
if(posUnitsPlot->widgetCB->GetSelected() == 0)
{
if(n == 1)
{
xsurfmin = (double)(evtheader.xpos);
ysurfmin = (double)(evtheader.ypos);
}
surfx[i] = (double)(evtheader.xpos);
surfy[i] = (double)(evtheader.ypos);
}
else if(posUnitsPlot->widgetCB->GetSelected() == 1)
{
if(n == 1)
{
xsurfmin = (double)(evtheader.xpos*lenconversion);
ysurfmin = (double)(evtheader.ypos*lenconversion);
}
surfx[i] = (double)(evtheader.xpos*lenconversion);
surfy[i] = (double)(evtheader.ypos*lenconversion);
}
// Save result for later plotting
if(norm)
integralCount[i] += ((double)m)/((double)k);
else
integralCount[i] += m;
inroot->Close();
delete inroot;
}
// Update the progress bar
progVal = (float)(75.00/nrfiles)*i;
analysisProgress->widgetPB->SetPosition(progVal);
gVirtualX->Update(1);
}
nrentries = n;
printf("SurfaceScan(): %d files were selected.\n", nrfiles);
// Check for Minimum and Maximum values and normalize to 1, if normalization is selected
if(norm)
{
minInteg = TMath::MinElement(nrfiles, integralCount);
for(int i = 0; i < nrfiles; i++)
{
integralCount[i] -= minInteg;
if(DBGSIG) printf("Subtraction: %lf\n", integralCount[i]);
}
maxInteg = TMath::MaxElement(nrfiles, integralCount);
for(int i = 0; i < nrfiles; i++)
{
integralCount[i] = integralCount[i]/maxInteg;
if(DBGSIG) printf("Normalization: %lf\n", integralCount[i]);
}
}
// Make the 2D surface plot
if(edit == 0)
gCanvas = new TCanvas("canv","canv",1100,900);
else
gCanvas = tempAnalysisCanvas->GetCanvas();
double range[4];
TGraph2D *gScan2D;
gScan2D = new TGraph2D();
range[0] = TMath::MinElement(nrfiles, surfx);
range[1] = TMath::MaxElement(nrfiles, surfx);
range[2] = TMath::MinElement(nrfiles, surfy);
range[3] = TMath::MaxElement(nrfiles, surfy);
for(int i = 0; i < nrfiles; i++)
{
if(DBGSIG)
{
if(surfScanOpt->widgetChBox[1]->IsDown())
printf("SurfaceScan(): %.2lf\t%.2lf\t%lf\n", surfx[i]-range[0], surfy[i]-range[2], integralCount[i]);
else
printf("SurfaceScan(): %.2lf\t%.2lf\t%lf\n", surfx[i], surfy[i], integralCount[i]);
}
if(surfScanOpt->widgetChBox[1]->IsDown())
gScan2D->SetPoint(i, surfx[i]-range[0], surfy[i]-range[2], integralCount[i]);
else
gScan2D->SetPoint(i, surfx[i], surfy[i], integralCount[i]);
// Update the progress bar
progVal = (float)(9.00/nrfiles)*i+90.00;
analysisProgress->widgetPB->SetPosition(progVal);
gVirtualX->Update(1);
}
if(surfScanOpt->widgetChBox[1]->IsDown())
{
range[1] -= range[0];
range[3] -= range[2];
range[0] -= range[0];
range[2] -= range[2];
}
gCanvas->cd();
gStyle->SetPalette(1);
gCanvas->SetLeftMargin(0.15);
gCanvas->SetRightMargin(0.126);
gCanvas->SetTopMargin(0.077);
gScan2D->Draw("COLZ");
gCanvas->Modified();
gCanvas->Update();
if(!cleanPlots)
{
if(posUnitsPlot->widgetCB->GetSelected() == 0)
gScan2D->SetTitle("Surface scan;X [table units];Y [table units]");
else if(posUnitsPlot->widgetCB->GetSelected() == 1)
gScan2D->SetTitle("Surface scan;X [#mum];Y [#mum]");
}
else
{
if(posUnitsPlot->widgetCB->GetSelected() == 0)
gScan2D->SetTitle(";X [table units];Y [table units]");
else if(posUnitsPlot->widgetCB->GetSelected() == 1)
gScan2D->SetTitle(";X [#mum];Y [#mum]");
}
/* TGaxis *xax = (TGaxis*)gScan2D->GetXaxis();
xax->SetMaxDigits(4);
TGaxis *yax = (TGaxis*)gScan2D->GetYaxis();
yax->SetMaxDigits(4);*/
gScan2D->SetNpx((int)resolSurf->widgetNE[0]->GetNumber());
gScan2D->SetNpy((int)resolSurf->widgetNE[1]->GetNumber());
gCanvas->Modified();
gCanvas->Update();
gScan2D->GetXaxis()->SetTitleOffset(1.3);
gScan2D->GetXaxis()->CenterTitle(kTRUE);
gScan2D->GetXaxis()->SetLabelSize(0.027);
gScan2D->GetXaxis()->SetLabelOffset(0.02);
gScan2D->GetXaxis()->SetRangeUser(range[0], range[1]);
gScan2D->GetXaxis()->SetNoExponent(kTRUE);
gScan2D->GetYaxis()->SetTitleOffset(1.9);
gScan2D->GetYaxis()->CenterTitle(kTRUE);
gScan2D->GetYaxis()->SetLabelSize(0.027);
gScan2D->GetYaxis()->SetLabelOffset(0.02);
gScan2D->GetYaxis()->SetRangeUser(range[2], range[3]);
gScan2D->GetYaxis()->SetNoExponent(kTRUE);
gCanvas->Modified();
gCanvas->Update();
remove_from_last((char*)files->At(0)->GetTitle(), '_', ctemp);
sprintf(exportname, "%s_surfscan.pdf", ctemp);
gCanvas->SaveAs(exportname);
// Update the progress bar
analysisProgress->widgetPB->SetPosition(100.0);
gVirtualX->Update(1);
if(edit == 0)
{
delete gScan2D;
delete gCanvas;
}
else if(edit == 1)
{
gCanvas->Modified();
gCanvas->Update();
}
}
}