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

Subversion Repositories f9daq

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