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

Subversion Repositories f9daq

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