WebSVN – f9daq – Blame – //lab/sipmscan/trunk/src/analysis.cpp

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

Subversion Repositories f9daq

(root)//lab/sipmscan/trunk/src/analysis.cpp – Rev 181