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

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

Subversion Repositories f9daq

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