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

Subversion Repositories f9daq

(root)/lab/sipmscan/trunk/src/analysis.cpp - Rev 173