Context Navigation

source: fact/tools/rootmacros/PulseTemplates/FCalcPulseTemplate.C@ 14782

Visit:

Last change on this file since 14782 was 14781, checked in by Jens Buss, 12 years ago
distribution histos filled
File size: 47.7 KB

Line
1	/******************** FPulseTemplate *********************
2	* open a root-file with pulse overlay histograms
3	* calculate most propable Pulse-Shape
4	* compute pulse height and pulse integral spektrum of the peaks
5	+ draw histos for single, double, tripple, ....above 100mV Photonpulses
6	+ draw Parameterdevelopment depending on photon quantity
7	+ form a template (shape depending on number of photons)
8	so it can be used for detection of other peaks
9	*************************************************************/
10
11	//----------------------------------------------------------------------------
12	// root libraries
13	//----------------------------------------------------------------------------
14
15	#include <TROOT.h>
16	#include <TCanvas.h>
17	#include <TTimer.h>
18	#include <Getline.h>
19	#include <TMath.h>
20	#include <TFile.h>
21	#include <TGraph.h>
22	#include <TGraphErrors.h>
23	#include <TStyle.h>
24	#include <TString.h>
25	#include <TF1.h>
26
27	#include <stdio.h>
28	#include <stdint.h>
29	#include <cstdio>
30	#include <iostream>
31	#include <fstream>
32
33	#include <vector>
34	#include <utility>
35
36	using namespace std;
37
38	#define NPIX 1440
39	#define NCELL 1024
40	#define FAD_MAX_SAMPLES 1024
41	#define MAX_PULS_ORDER 1
42	#define HAVE_ZLIB
43
44	//----------------------------------------------------------------------------
45	// Classes
46	//----------------------------------------------------------------------------
47	#include "rootfilehandler.h"
48	#include "pixel.h"
49	#include "pixelsum.h"
50	#include "pulse.h"
51	#include "templateextractors.h"
52	#include "configfile.h"
53	#include "csv.h"
54
55	//----------------------------------------------------------------------------
56	// Decleration of global variables
57	//----------------------------------------------------------------------------
58	bool breakout=false;
59
60	//----------------------------------------------------------------------------
61	// Initialisation of Root Objects
62	//----------------------------------------------------------------------------
63
64	// Temporary Objects
65
66	//Pixelwise Histograms
67
68	//All Pixel Histograms
69
70	//Root-File Objects
71	TObjArray* hAllPixelList = NULL;
72	TObjArray* hRootList = NULL;
73
74	//Canvases
75	TCanvas** cgpPixelPulses = NULL;
76	TCanvas** cgpAllPixelPulses = NULL;
77	TCanvas** cgpDistributions = NULL;
78	TCanvas* cgpGraphs = NULL;
79	//TCanvas* cgpTestHistos = NULL;
80
81	//----------------------------------------------------------------------------
82	// Functions
83	//----------------------------------------------------------------------------
84	void DeletePixelCanvases( int, int );
85	void
86	UpdateCanvases(
87	int verbosityLevel,
88	int max_pulse_order,
89	bool testmode
90	);
91
92	//----------------------------------------------------------------------------
93	//----------------------------------------------------------------------------
94	// MAIN - Funtion
95	//----------------------------------------------------------------------------
96	//----------------------------------------------------------------------------
97	int FCalcPulseTemplate(
98	TString InRootFileName = "20120309_017.root",
99	// TString InRootFileName = "20120309_032.root",
100	TString InputPath = "analysis/analysis/FPulseTemplate/20120309_017/Overlay/",
101	// TString InputPath = "analysis/analysis/FPulseTemplate/20120309_032/Overlay/",
102	TString OutputRootFileName = "test.root",
103	TString OutPutPath = "analysis/FPulseTemplate/20120309_017/",
104	int firstpixel = 0,
105	int npixel = 100,
106	int pixelSetSize = 200,
107	int maxPulseOrder = 1,
108	// TString histoOptions = "SRM",
109	bool ProduceGraphic = true,
110	bool stats = true,
111	bool saveResults = false,
112	// bool fitdata = false,
113	bool debugPixel = false,
114	// int refresh_rate = 500, //refresh rate for canvases
115	int verbosityLevel = 1 // different verbosity levels can be implemented here
116	)
117	{
118
119	InputPath = SetHostsPaths(true, InputPath );
120	OutPutPath = SetHostsPaths(true, OutPutPath );
121
122	TString outFile = OutPutPath;
123	outFile += OutputRootFileName;
124
125	//----------------------------------------------------------------------------
126	// Open-Root-File Settings
127	//----------------------------------------------------------------------------
128	if (verbosityLevel > 0)
129	{
130	cout << endl << "...load root file" << endl;
131	}
132
133	TFile * inputRootFile
134	= LoadRootFile( InputPath, InRootFileName, verbosityLevel );
135
136	if (inputRootFile == NULL)
137	{
138	cerr << "input file not readable, check file path!!!" << endl;
139	return(0);
140	}
141	if (saveResults)
142	{
143	CreateRootFile( outFile, true, verbosityLevel );
144	}
145
146	TFile * outputRootFile
147	= OpenRootFile( OutPutPath, OutputRootFileName, verbosityLevel );
148	//----------------------------------------------------------------------------
149	// Define operation range
150	//----------------------------------------------------------------------------
151	if ( npixel == -1 )
152	{
153	npixel = NPIX - firstpixel;
154	}
155
156
157	if ( pixelSetSize == -1 )
158	{
159	pixelSetSize = firstpixel +npixel;
160	}
161
162	// float GainMean = GainMean; // this has to be extracted from root files
163	// float BSLMean = BSLMean; // this has to be extracted from root files
164
165	//----------------------------------------------------------------------------
166	// root global Settings
167	//----------------------------------------------------------------------------
168	if (verbosityLevel > 0)
169	{
170	cout << endl << "...setting up root environment" ;
171	}
172
173	gStyle->SetPalette(1,0);
174	gROOT->SetStyle("Plain");
175	// gPad->SetGrid();
176
177
178
179	//----------------------------------------------------------------------------
180	// root Canvas Settings
181	//----------------------------------------------------------------------------
182	if (verbosityLevel > 0)
183	{
184	cout << endl << "...preparing canvases" << endl;
185	}
186
187	//Canvas Pad numbering
188	int PixelCanvasFrameNrs[8] =
189	{
190	1, // Top left
191	2, // Top mid left
192	3, // Top mid right
193	4, // Top right
194	5, // bootom left
195	6, // bootom mid left
196	7, // bottom mid right
197	8 // bootom right
198	};
199
200	// //Canvas Pad numbering
201	// int DistributionCanvasFrameNrs[4] = {
202	// 1, // Top left
203	// 2, // Top right
204	// 3, // bottom left
205	// 4 // bootom right
206	// };
207
208	if (ProduceGraphic)
209	{
210
211	//Canvases
212	cgpPixelPulses = new TCanvas*[maxPulseOrder];
213	cgpDistributions = new TCanvas*[maxPulseOrder];
214	cgpGraphs = new TCanvas("graphs", "Dependences to pe" );
215
216	//TCanvas* gpcDevelopment = NULL;
217	TString cName = "";
218	TString cTitle = "";
219
220	//naming of pulse canvases
221	for (
222	int pulse_order = maxPulseOrder - 1;
223	pulse_order >= 0 ;
224	pulse_order--
225	)
226	{
227	cName ="cgpDistributions";
228	cName += pulse_order;
229
230	cTitle ="Distributions of Pulses with Order of: ";
231	cTitle += pulse_order;
232
233	cgpDistributions[pulse_order]
234	= new TCanvas(cName,cTitle, 720,pulse_order*20,720,720);
235	// cgpDistributions[pulse_order]->Divide(2, 2);
236
237	cName ="cgpPixelPulses";
238	cName += pulse_order;
239
240	cTitle ="Overlays of Pulses with Order of: ";
241	cTitle += pulse_order;
242
243	cgpPixelPulses[pulse_order]
244	= new TCanvas(cName,cTitle, 0,pulse_order*20,800,800);
245	cgpPixelPulses[pulse_order]->Divide(4, 2);
246	}
247
248	cgpGraphs->Divide(1,3);
249	}
250
251	//-----------------------------------------------------------------------------
252	// Filter-Settings
253	//-----------------------------------------------------------------------------
254
255	//-----------------------------------------------------------------------------
256	// Distrribution Histograms
257	//-----------------------------------------------------------------------------
258	TH1D hBsl[maxPulseOrder];
259	TH1D hAsymptote[maxPulseOrder];
260	TH1D hT0[maxPulseOrder];
261	TH1D hTau1[maxPulseOrder];
262	TH1D hTau2[maxPulseOrder];
263	TH1D hIntegral[maxPulseOrder];
264	TH1D hAmplitude[maxPulseOrder];
265	TH1D hRiseTime[maxPulseOrder];
266
267	//Lopp over pulse order configure distribution histograms
268	for (int i = 0; i<maxPulseOrder; i++)
269	{
270	//BSL
271	TString title;
272	title = "hBSL";
273	title.Append(i);
274	hBsl[i].SetNameTitle(title, title);
275	hBsl[i].SetBins(400, -2-0.005, 2-0.005);
276
277	//Asymptote
278	title = "hAsymptote";
279	title.Append(i);
280	hAsymptote[i].SetNameTitle(title, title);
281	hAsymptote[i].SetBins(60, -10.5, 49.5);
282
283	//T0
284	title = "hT0";
285	title.Append(i);
286	hT0[i].SetNameTitle(title, title);
287	hT0[i].SetBins(80, 39.5, 119.5);
288
289	//Tau1
290	title = "hTau1";
291	title.Append(i);
292	hTau1[i].SetNameTitle(title, title);
293	hTau1[i].SetBins(10, -10.5, 9.5);
294
295	//Tau2
296	title = "hTau2";
297	title.Append(i);
298	hTau2[i].SetNameTitle(title, title);
299	hTau2[i].SetBins(40, 19.5, 59.5);
300
301	//Integral
302	title = "hIntegral";
303	title.Append(i);
304	hIntegral[i].SetNameTitle(title, title);
305	hIntegral[i].SetBins(500, (i+1)-0.5, (i+1)499.5);
306
307	//Amplitude
308	title = "hAmplitude";
309	title.Append(i);
310	hAmplitude[i].SetNameTitle(title, title);
311	hAmplitude[i].SetBins(400, (i+1)4.95, (i+1)24.95);
312
313	//hRiseTime
314	title = "hRiseTime";
315	title.Append(i);
316	hRiseTime[i].SetNameTitle(title, title);
317	hRiseTime[i].SetBins(20, -0, 10);
318	}
319
320	// Error Histograms
321
322	TH1D hMedianErrors[maxPulseOrder];
323	TH1D hMeanErrors[maxPulseOrder];
324	TH1D hMaxErrors[maxPulseOrder];
325
326	//Lopp over pulse order configure distribution histograms
327	for (int i = 0; i<maxPulseOrder; i++)
328	{
329	//hMedianErrors
330	TString title;
331	title = "hMedianErrors";
332	title.Append(i);
333	hMedianErrors[i].SetNameTitle(title, title);
334	hMedianErrors[i].SetBins(300, -0.5, 300.5);
335
336	//hMeanErrors
337	title = "hMeanErrors";
338	title.Append(i);
339	hMeanErrors[i].SetNameTitle(title, title);
340	hMeanErrors[i].SetBins(300, -0.5, 300.5);
341
342	//hMaxErrors
343	title = "hMaxErrors";
344	title.Append(i);
345	hMaxErrors[i].SetNameTitle(title, title);
346	hMaxErrors[i].SetBins(300, -0.5, 300.5);
347	}
348
349
350	//-----------------------------------------------------------------------------
351	// Graphs
352	//-----------------------------------------------------------------------------
353	TGraphErrors grPheVsAmpl;
354	grPheVsAmpl.SetNameTitle("grPheVsAmpl", "Distribution of Pulse Amplitude vs. Photon equivalent");
355	grPheVsAmpl.SetMarkerStyle(21);
356
357	TGraphErrors grPheVsIntegral;
358	grPheVsIntegral.SetNameTitle("grPheVsIntegral", "Distribution of Pulse Integral vs. Photon equivalent");
359	grPheVsIntegral.SetMarkerStyle(21);
360
361	TGraphErrors grPheVsTau2;
362	grPheVsTau2.SetNameTitle("grPheVsTau2", "Distribution of Tau2 vs. Photon equivalent");
363	grPheVsTau2.SetMarkerStyle(21);
364	//-----------------------------------------------------------------------------
365	// Variables
366	//-----------------------------------------------------------------------------
367
368
369	//----------------------------------------------------------------------------
370	// Initialize Pixels
371	//----------------------------------------------------------------------------
372	if (verbosityLevel > 0)
373	{
374	cout << endl << "...preparing pixels" << endl;
375	}
376
377	Pixel** pixel = new Pixel*[NPIX];
378
379	for (int i = 0 ; i < NPIX; i++)
380	{
381	pixel[i] = NULL;
382	}
383
384	PixelSum* wholeCamera = NULL;
385
386	bool first_pass = true;
387
388	//----------------------------------------------------------------------------
389	// Initialize Outputfiles
390	//----------------------------------------------------------------------------
391
392	Csv PixelCsv(OutPutPath, InRootFileName.Remove(12, 5), "points", verbosityLevel);
393	PixelCsv.WritePointSetHeader();
394	Csv PixelModelCsv(OutPutPath, InRootFileName.Remove(12, 5), "fitResults", verbosityLevel);
395	PixelModelCsv.WritePulseAttributesHeader();
396	// Csv AllPixelCsv(OutPutPath, InRootFileName.Remove(12, 5), -1, verbosityLevel);
397	// Csv AllPixelModelCsv(OutPutPath, InRootFileName.Remove(12, 5), "AllPixelModel", verbosityLevel);
398
399	TCanvas pulses("Templates","Templates", 0,0,1200,600);
400	pulses.Divide(3,1);
401	//-------------------------------------
402	// Loop over Pixel Sets
403	//-------------------------------------
404	for ( int firstPixelOfSet = firstpixel;
405	firstPixelOfSet < firstpixel + npixel;
406	firstPixelOfSet += pixelSetSize )
407	{
408	if (verbosityLevel >= 0)
409	{
410	cout << "------------------------------------------------"
411	<< endl
412	<< "...processing Pixel: "
413	<< firstPixelOfSet
414	<< " to Pixel: "
415	<< firstPixelOfSet+pixelSetSize-1
416	<< endl;
417	}
418
419	//--------------------------------------------------------------------
420	// Loop over every Pixel of a Set of Pixels
421	//--------------------------------------------------------------------
422	for ( int pixelID = firstPixelOfSet;
423	pixelID < firstPixelOfSet + pixelSetSize
424	&& pixelID < firstpixel + npixel;
425	pixelID++ )
426	{
427	if (verbosityLevel >= 0)
428	{
429	cout << "------------------------------------------------"
430	<< endl
431	<< "...processing Pixel: "
432	<< pixelID << "/"
433	<< firstPixelOfSet + pixelSetSize - 1
434	<< endl;
435	}
436	if (verbosityLevel > 1)
437	{
438	cout << "-------------------------------------"
439	<< endl
440	<< "...processing Set from Pixel "
441	<< firstPixelOfSet
442	<< " to Pixel "
443	<< firstPixelOfSet+pixelSetSize-1
444	<< " Pixel: " << pixelID
445	<< "/" << firstpixel + npixel -1
446	<< endl;
447	}
448
449	//-------------------------------------
450	// Create individual Pixel
451	//-------------------------------------
452	if (verbosityLevel > 3)
453	{
454	cout << "...creating pixel: " << pixelID << endl;
455	}
456
457	pixel[pixelID] = new Pixel(
458	pixelID,
459	maxPulseOrder, ///One should get this from the file
460	verbosityLevel,
461	stats,
462	"L",
463	70, //pixelOverlayXaxisLeft /TODO: get it from the root file
464	230, //pixelOverlayXaxisRight /TODO: get it from the root file
465	-1, //bslMean /TODO: get it from the root file
466	9, //gain mean /TODO: get it from the root file
467	inputRootFile,
468	outputRootFile
469	);
470
471	if (breakout) break;
472
473	//Preparing Camera
474	if (first_pass)
475	{
476	if (verbosityLevel > 0)
477	{
478	cout << endl << "...preparing camera" << endl;
479	}
480
481	wholeCamera = new PixelSum(
482	"AllPixel",
483	1440,
484	maxPulseOrder,
485	verbosityLevel,
486	stats,
487	"C",
488	pixel[pixelID]->mPixelOverlayXaxisLeft,
489	pixel[pixelID]->mPixelOverlayXaxisRight ,
490	pixel[pixelID]->mBSLMean ,
491	pixel[pixelID]->mGainMean ,
492	outputRootFile
493	);
494
495	first_pass = false;
496	}
497	//-------------------------------------
498	// Histogramms of Maximas in Overlay Spectra
499	//-------------------------------------
500
501	for ( int pulse_order = 0 ;
502	pulse_order < pixel[pixelID]->mMaxPulseOrder ;
503	pulse_order ++)
504	{
505	if (verbosityLevel > 2)
506	{
507	cout << "-------------------------------------"
508	<< endl
509	<< "...processing Set from Pixel "
510	<< firstPixelOfSet
511	<< " to Pixel "
512	<< firstPixelOfSet+pixelSetSize-1
513	<< " Pixel: " << pixelID
514	<< "/" << firstpixel + npixel -1
515	<< " Pulse-Order: " << pulse_order
516	<< endl;
517	}
518
519	// Calculate Max Prop. Value of each slice
520	//-------------------------------------
521
522	//from Maximum Overlay
523	if (verbosityLevel > 2)
524	{
525	cout << "...extracting templates from Maximum Overlay "
526	<< endl;
527	}
528
529	ExtractPulseTemplate(
530	pixel[pixelID],
531	"Maximum",
532	pulse_order,
533	10,
534	verbosityLevel
535	);
536
537	//from Edge Overlay
538	if (verbosityLevel > 2)
539	{
540	cout << "...extracting templates from Edge Overlay "
541	<< endl;
542	}
543	ExtractPulseTemplate(
544	pixel[pixelID],
545	"Edge",
546	pulse_order,
547	10,
548	verbosityLevel
549	);
550
551
552	// ShiftStartOfHistoInXTo(
553	// pixel[pixelID]->hEdgeOverlay[pulse_order],
554	// 0
555	// );
556
557	// ShiftStartOfHistoInXTo(
558	// pixel[pixelID]->hMaxOverlay[pulse_order],
559	// 0
560	// );
561
562	PixelCsv.WritePointSet(
563	pixel[pixelID],
564	"Maximum",
565	pulse_order
566	);
567
568	PixelCsv.WritePointSet(
569	pixel[pixelID],
570	"Edge",
571	pulse_order
572	);
573
574	if ( saveResults )
575	{
576	WritePixelTemplateToCsv(
577	pixel[pixelID],
578	OutPutPath,
579	"Maximum",
580	pulse_order,
581	verbosityLevel
582	);
583
584	WritePixelTemplateToCsv(
585	pixel[pixelID],
586	OutPutPath,
587	"Edge",
588	pulse_order,
589	verbosityLevel
590	);
591	}
592
593	if (ProduceGraphic && debugPixel)
594	{
595	// pixel[pixelID]->hPixelEdgeMax[pulse_order]->GetXaxis()->SetLimits(
596	// 0,
597	// 300
598	// );
599	// pixel[pixelID]->hPixelEdgeMean[pulse_order]->GetXaxis()->SetLimits(
600	// 0,
601	// 300
602	// );
603	// pixel[pixelID]->hPixelEdgeMedian[pulse_order]->GetXaxis()->SetLimits(
604	// 0,
605	// 300
606	// );
607	// pixel[pixelID]->hPixelMax[pulse_order]->GetXaxis()->SetLimits(
608	// 0,
609	// 300
610	// );
611	// pixel[pixelID]->hPixelMean[pulse_order]->GetXaxis()->SetLimits(
612	// 0,
613	// 300
614	// );
615	// pixel[pixelID]->hPixelMedian[pulse_order]->GetXaxis()->SetLimits(
616	// 0,
617	// 300
618	// );
619
620	// pixel[pixelID]->hEdgeProfile[pulse_order]->GetXaxis()->SetLimits(
621	// 0,
622	// 300
623	// );
624
625
626
627
628	// pixel[pixelID]->ShiftHistoInY(
629	// pixel[pixelID]->hPixelEdgeMax[pulse_order],
630	// 0.25
631	// );
632
633
634	pixel[pixelID]->DrawTemplateHistograms(
635	cgpPixelPulses,
636	PixelCanvasFrameNrs
637	);
638
639	pixel[pixelID]->DrawEdgeTemplateHistograms(
640	cgpPixelPulses,
641	PixelCanvasFrameNrs
642	);
643
644	// cgpPixelPulses[pulse_order]->cd(6);
645
646
647	//-------------------------------------
648	//-------------------------------------
649	// Test Area
650	//-------------------------------------
651	//-------------------------------------
652
653	cgpPixelPulses[pulse_order]->cd(8);
654
655	// pixel[pixelID]->hEdgeProfile[pulse_order]->SetLineColor(kBlack);
656	// Pulse pulseMean("PulseMean", pixel[pixelID]->hEdgeProfile[pulse_order], 0);
657
658
659
660
661	// double fit_parameters5[3];
662	// fit_parameters5[0] = fit_parameters[0];
663	// fit_parameters5[0] = fit_parameters[0];
664	// FitPulse(
665	// "hugo3",
666	// pixel[pixelID]->hPixelEdgeMax[pulse_order],
667	// pixel[pixelID]->hPixelEdgeMax[pulse_order]->GetMaximumBin()-12,
668	// pixel[pixelID]->hPixelEdgeMax[pulse_order]->GetMaximumBin(),
669	// fit_parameters5
670	// );
671	// cout << "Parameters\t" << fit_parameters[0] << "\n"
672	// << fit_parameters[1] << "\n"
673	// << fit_parameters[2] << "\n" << endl;
674
675	// TF1 *func = new TF1("func", template_function, 0, 300, 10);
676
677	// func->SetParameters(
678	// -0.5, // bsl
679	// pixel[pixelID]->hPixelEdgeMax[pulse_order]->GetXaxis()->GetFirst() + 50,
680	// pixel[pixelID]->hPixelEdgeMax[pulse_order]->GetMaximumBin()
681	// ); // range
682
683
684
685	// func->SetParNames("baseline",
686	// "begin of pol3", "begin of exp",
687	// "exp-factor", "exp-tau", "ext-t0",
688	// "pol3_0", "pol3_1", "pol3_2", "pol3_3"
689	// );
690	// pixel[pixelID]->hPixelEdgeMax[pulse_order]->Fit(func);
691
692	// TH1F* hTest = new TH1F(
693	// "hTest",
694	// "Test",
695	// pixel[pixelID]->mPixelOverlayXaxisLeft
696	// + pixel[pixelID]->mPixelOverlayXaxisRight ,
697	// (-1*pixel[pixelID]->mPixelOverlayXaxisLeft)-0.5,
698	// pixel[pixelID]->mPixelOverlayXaxisRight-0.5
699	// );
700
701	// hTest->GetXaxis()->SetLimits(
702	// 0,
703	// 300
704	// );
705
706	// for (int bin = 0;
707	// bin < pixel[pixelID]->hPixelEdgeMax[pulse_order]->GetXaxis()->GetLast();
708	// bin++
709	// )
710	// {
711	// hTest->SetBinContent( bin, (-1)(fit_parameters[0]+TMath::Exp(fit_parameters[1]+fit_parameters[2](bin))) );
712	// }
713
714	// hTest->Add(pixel[pixelID]->hPixelEdgeMax[pulse_order], 1);
715
716	// cgpPixelPulses[pulse_order]->cd(7);
717	// hTest->Draw();
718
719	// for (int bin = 0;
720	// bin < pixel[pixelID]->hPixelEdgeMax[pulse_order]->GetXaxis()->GetLast();
721	// bin++
722	// )
723	// {
724	// hTest->SetBinContent( bin, (-1)*(hTest->GetBinContent(bin)) );
725	// }
726
727
728	// double fit_parameters2[3];
729	// FitFallingEdge(
730	// "hugo2",
731	// hTest,
732	// hTest->GetXaxis()->GetFirst()+68,
733	// hTest->GetXaxis()->GetFirst()+80,
734	// fit_parameters2
735	// );
736
737	// cgpPixelPulses[pulse_order]->cd(8);
738
739	// TH1F* hTest2 = new TH1F(
740	// "hTest2",
741	// "Test",
742	// pixel[pixelID]->mPixelOverlayXaxisLeft
743	// + pixel[pixelID]->mPixelOverlayXaxisRight ,
744	// (-1*pixel[pixelID]->mPixelOverlayXaxisLeft)-0.5,
745	// pixel[pixelID]->mPixelOverlayXaxisRight-0.5
746	// );
747
748	// hTest2->GetXaxis()->SetLimits(
749	// 0,
750	// 300
751	// );
752
753	// for (int bin = 0;
754	// bin < pixel[pixelID]->hPixelEdgeMax[pulse_order]->GetXaxis()->GetLast();
755	// bin++
756	// )
757	// {
758	// hTest2->SetBinContent( bin, (-1)*(hTest->GetBinContent(bin)) );
759	// }
760
761	// double fit_parameters3[3];
762	// FitRisingEdge(
763	// "hugo3",
764	// hTest2,
765	// hTest2->GetXaxis()->GetFirst()+68,
766	// hTest2->GetXaxis()->GetFirst()+80,
767	// fit_parameters3
768	// );
769
770
771	// hTest2->Draw();
772
773
774	//-------------------------------------
775	//-------------------------------------
776	// EOF Test Area
777	//-------------------------------------
778	//-------------------------------------
779	}
780
781	// =======================================================
782	// Fit of pulse function to template
783
784	Pulse* pulseFits[6];
785
786	// pulses.cd(1);
787	Pulse maxMaxPulse("maxMaxPulse", pixel[pixelID]->hPixelMax[pulse_order], "Q", 1);
788	pulseFits[0] = &maxMaxPulse;
789
790	// pulses.cd(2);
791	Pulse maxMedianPulse("maxMedianPulse", pixel[pixelID]->hPixelMedian[pulse_order], "Q",1);
792	pulseFits[1] = &maxMedianPulse;
793
794	// pulses.cd(3);
795	Pulse maxMeanPulse("maxMeanPulse", pixel[pixelID]->hPixelMean[pulse_order], "Q",1);
796	pulseFits[2] = &maxMeanPulse;
797
798	Pulse edgeMaxPulse("edgeMaxPulse", pixel[pixelID]->hPixelEdgeMax[pulse_order], "Q",1);
799	pulseFits[3] = &edgeMaxPulse;
800
801	Pulse edgeMedianPulse("edgeMedianPulse", pixel[pixelID]->hPixelEdgeMedian[pulse_order], "Q",1);
802	pulseFits[4] = &edgeMedianPulse;
803
804	Pulse edgeMeanPulse("edgeMeanPulse", pixel[pixelID]->hPixelEdgeMean[pulse_order], "Q",1);
805	pulseFits[5] = &edgeMeanPulse;
806
807
808	//--------------------------------------------------------
809	//Write parameters from pulse fit to CSV
810
811	for (int i = 0; i < 6; i++)
812	{
813	PixelModelCsv.WritePulseAttributes(
814	pixel[pixelID],
815	pulseFits[i],
816	"Edge_pix",
817	pulse_order
818	);
819	}
820	//-------------------------------------
821	// Fill Error Histogramms
822	//-------------------------------------
823
824	for (int i = pixel[pixelID]->hPixelEdgeMedian[pulse_order]->GetXaxis()->GetFirst();
825	i <= pixel[pixelID]->hPixelEdgeMedian[pulse_order]->GetXaxis()->GetLast();
826	i++)
827	{
828
829	hMedianErrors[pulse_order].SetBinContent(
830	i,
831	pixel[pixelID]->hPixelEdgeMedian[pulse_order]->GetBinError(i)
832	);
833
834	hMeanErrors[pulse_order].SetBinContent(
835	i,
836	pixel[pixelID]->hPixelEdgeMean[pulse_order]->GetBinError(i)
837	);
838
839	hMaxErrors[pulse_order].SetBinContent(
840	i,
841	pixel[pixelID]->hPixelEdgeMax[pulse_order]->GetBinError(i)
842	);
843	}
844
845	if (ProduceGraphic && debugPixel)
846	{
847	int current_order = 0;
848	pulses.cd(1);
849	hMaxErrors[current_order].Draw();
850
851	pulses.cd(2);
852	hMedianErrors[current_order].Draw();
853
854	pulses.cd(3);
855	hMeanErrors[current_order].Draw();
856	pulses.Update();
857	}
858
859	//-------------------------------------
860	// Fill Distribution Histogramms
861	//-------------------------------------
862
863	hBsl[pulse_order].Fill( edgeMedianPulse.GetBsl() );
864	hAsymptote[pulse_order].Fill( edgeMedianPulse.GetHeight() );
865	hT0[pulse_order].Fill( edgeMedianPulse.GetT0() );
866	hTau1[pulse_order].Fill( edgeMedianPulse.GetTau1() );
867	hTau2[pulse_order].Fill( edgeMedianPulse.GetTau2() );
868	hIntegral[pulse_order].Fill( edgeMedianPulse.GetIntegral());
869	hAmplitude[pulse_order].Fill( edgeMedianPulse.GetAmplitude());
870	hRiseTime[pulse_order].Fill( edgeMedianPulse.GetRiseTime());
871
872	//-------------------------------------
873	// Fill Distribution Graphs
874	//-------------------------------------
875
876	grPheVsAmpl.Set(grPheVsAmpl.GetN()+1);
877	grPheVsIntegral.Set(grPheVsIntegral.GetN()+1);
878
879
880	grPheVsAmpl.SetPoint(
881	grPheVsAmpl.GetN()-1,
882	edgeMedianPulse.GetPhE(),
883	edgeMedianPulse.GetAmplitude()
884	);
885
886	grPheVsAmpl.SetPointError(
887	grPheVsAmpl.GetN()-1,
888	edgeMedianPulse.GetPhEErr(),
889	edgeMedianPulse.GetAmplitudeErr()
890	);
891
892	grPheVsIntegral.SetPoint(
893	grPheVsIntegral.GetN()-1,
894	edgeMedianPulse.GetPhE(),
895	edgeMedianPulse.GetIntegral()
896	);
897
898	grPheVsAmpl.SetPointError(
899	grPheVsIntegral.GetN()-1,
900	edgeMedianPulse.GetPhEErr(),
901	edgeMedianPulse.GetIntegralErr()
902	);
903
904	grPheVsTau2.SetPoint(
905	grPheVsIntegral.GetN()-1,
906	edgeMedianPulse.GetPhE(),
907	edgeMedianPulse.GetTau2()
908	);
909
910	grPheVsTau2.SetPointError(
911	grPheVsIntegral.GetN()-1,
912	edgeMedianPulse.GetPhE(),
913	edgeMedianPulse.GetTau2Err()
914	);
915
916	// cgpDistributions[pulse_order].cd(1);
917
918
919	// FitMaxPropabilityPuls(
920	// pixel[pixelID]->hPixelEdgeMean[pulse_order],
921	// verbosityLevel
922	// );
923	//-------------------------------------
924	// Fill Histogramms of Camera
925	//-------------------------------------
926
927	wholeCamera->hMaxOverlay[pulse_order]->Add(
928	pixel[pixelID]->hMaxOverlay[pulse_order]
929	);
930
931	wholeCamera->hMaxProfile[pulse_order]->Add(
932	pixel[pixelID]->hMaxProfile[pulse_order]
933	);
934
935	wholeCamera->hEdgeOverlay[pulse_order]->Add(
936	pixel[pixelID]->hEdgeOverlay[pulse_order]
937	);
938
939	wholeCamera->hEdgeProfile[pulse_order]->Add(
940	pixel[pixelID]->hEdgeProfile[pulse_order]
941	);
942
943	//-------------------------------------
944	// Comparisons
945	//-------------------------------------
946
947	//chi2 test
948	// float chi2 =
949	// pixel[pixelID]->hPixelEdgeMean[pulse_order]->Chi2Test(
950	// pixel[pixelID]->hPixelMean[pulse_order],
951	// "UUPCHI2"
952	// );
953	// cout << "p-Value :" << chi2 << endl;
954	// wholeCamera->hChi2EdgetoMax[pulse_order]->Fill(chi2);
955	// cgpDistributions[pulse_order]->cd();
956	// wholeCamera->hChi2EdgetoMax[pulse_order]->Draw();
957
958	if (verbosityLevel > 2)
959	{
960	cout << endl
961	<< "...End of pulseorder "
962	<< pulse_order
963	<< endl;
964	}
965
966	}
967	// End of Loop over pulsorder of current Pixel
968
969	if (ProduceGraphic && debugPixel)
970	{
971	UpdateCanvases(
972	verbosityLevel,
973	MAX_PULS_ORDER,
974	false
975	);
976	}
977
978	if (debugPixel)
979	{
980	//Process gui events asynchronously during input
981	cout << endl;
982	TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
983	timer.TurnOn();
984	TString input = Getline("Type 'q' to exit, Type 's' to exit set, <return> to go on: ");
985	timer.TurnOff();
986	if (input=="q\n")
987	{
988	return(0);
989	}
990	if (input=="s\n")
991	{
992	break;
993	}
994	}
995
996	if ( saveResults )
997	{
998	pixel[pixelID]->SavePixelHistograms(
999	outFile,
1000	saveResults
1001	);
1002	}
1003
1004	//deleteCurrent Pixel from Heap
1005	delete pixel[pixelID];
1006	pixel[pixelID] = NULL;
1007
1008	if (verbosityLevel > 2)
1009	{
1010	cout << endl
1011	<< "...End of Pixel"
1012	<< endl
1013	<< "------------------------------------------------"
1014	<< endl;
1015	}
1016	}
1017	// End of Loop over all Pixels of set
1018
1019	if (verbosityLevel > 1)
1020	{
1021	cout << endl
1022	<< "...End of Loop over all Pixels of set"
1023	<< endl
1024	<< "------------------------------------------------"
1025	<< endl;
1026	}
1027	}
1028	// End of Loop over all Pixelsets
1029
1030	if (verbosityLevel > 0)
1031	{
1032	cout << endl
1033	<< "...End of Loop over all Pixelsets"
1034	<< endl
1035	<< "------------------------------------------------"
1036	<< endl;
1037	}
1038
1039	delete[] pixel;
1040	pixel = NULL;
1041
1042	//-------------------------------------
1043	// Draw All Pixel Histograms
1044	//-------------------------------------
1045
1046	for ( int pulse_order = 0 ;
1047	pulse_order < wholeCamera->mMaxPulseOrder ;
1048	pulse_order ++)
1049	{
1050	if (verbosityLevel > 2)
1051	{
1052	cout << "-------------------------------------"
1053	<< endl
1054	<< "...processing Pulse-Order: "
1055	<< pulse_order;
1056	}
1057
1058	// Calculate Max Prop. Value of each slice
1059	//-------------------------------------
1060
1061	//from Maximum Overlay
1062	ExtractPulseTemplate(
1063	wholeCamera,
1064	"Maximum",
1065	pulse_order,
1066	10,
1067	verbosityLevel
1068	);
1069	if ( saveResults )
1070	{
1071	WritePixelTemplateToCsv(
1072	wholeCamera,
1073	OutPutPath,
1074	"Maximum",
1075	pulse_order,
1076	verbosityLevel
1077	);
1078	}
1079
1080	//from Edge Overlay
1081	ExtractPulseTemplate(
1082	wholeCamera,
1083	"Edge",
1084	pulse_order,
1085	10,
1086	verbosityLevel
1087	);
1088	if ( saveResults )
1089	{
1090	WritePixelTemplateToCsv(
1091	wholeCamera,
1092	OutPutPath,
1093	"Edge",
1094	pulse_order,
1095	verbosityLevel
1096	);
1097	}
1098
1099	if (ProduceGraphic)
1100	{
1101	wholeCamera->DrawTemplateHistograms(
1102	cgpPixelPulses,
1103	PixelCanvasFrameNrs
1104	);
1105
1106	wholeCamera->DrawEdgeTemplateHistograms(
1107	cgpPixelPulses,
1108	PixelCanvasFrameNrs
1109	);
1110	}
1111
1112	ShiftStartOfHistoInXTo(
1113	wholeCamera->hEdgeOverlay[0],
1114	0
1115	);
1116
1117	ShiftStartOfHistoInXTo(
1118	wholeCamera->hMaxOverlay[0],
1119	0
1120	);
1121
1122	// =======================================================
1123	// Fit of pulse function to template
1124
1125	Pulse* allPulseFits[6];
1126
1127	pulses.cd(1);
1128	Pulse maxCamMaxPulse("maxMaxPulse", wholeCamera->hPixelMax[pulse_order], "Q", 1);
1129	allPulseFits[0] = &maxCamMaxPulse;
1130
1131	pulses.cd(2);
1132	Pulse maxCamMedianPulse("maxMedianPulse", wholeCamera->hPixelMedian[pulse_order], "Q", 1);
1133	allPulseFits[1] = &maxCamMedianPulse;
1134
1135	pulses.cd(3);
1136	Pulse maxCamMeanPulse("maxMeanPulse", wholeCamera->hPixelMean[pulse_order],"Q", 1);
1137	allPulseFits[2] = &maxCamMeanPulse;
1138
1139	Pulse edgeCamMaxPulse("edgeMaxPulse", wholeCamera->hPixelEdgeMax[pulse_order], "Q", 1);
1140	allPulseFits[3] = &edgeCamMaxPulse;
1141
1142	Pulse edgeCamMedianPulse("edgeMedianPulse", wholeCamera->hPixelEdgeMedian[pulse_order],"Q", 1);
1143	allPulseFits[4] = &edgeCamMedianPulse;
1144
1145	Pulse edgeCamMeanPulse("edgeMeanPulse", wholeCamera->hPixelEdgeMean[pulse_order], "Q", 1);
1146	allPulseFits[5] = &edgeCamMeanPulse;
1147
1148
1149	//--------------------------------------------------------
1150	//Write parameters from pulse fit to CSV
1151
1152	for (int i = 0; i < 6; i++)
1153	{
1154	PixelModelCsv.WritePulseAttributes(
1155	wholeCamera,
1156	allPulseFits[i],
1157	"Edge",
1158	pulse_order
1159	);
1160	}
1161	} //EOF: Draw All Pixel Histograms
1162
1163
1164	//-------------------------------------
1165	// Draw Parameter Distributions
1166	//-------------------------------------
1167
1168	if (ProduceGraphic)
1169	{
1170	cgpGraphs->cd(1);
1171
1172	if (verbosityLevel > 0) cout << "...drawing grPheVsAmpl" << endl;
1173	grPheVsAmpl.Draw("AP");
1174
1175	cgpGraphs->cd(2);
1176	if (verbosityLevel > 0) cout << "...drawing grPheVsIntegral" << endl;
1177	grPheVsIntegral.Draw("AP");
1178
1179	cgpGraphs->cd(3);
1180	if (verbosityLevel > 0) cout << "...drawing grPheVsTau2" << endl;
1181	grPheVsTau2.Draw("AP");
1182
1183	cgpGraphs->Update();
1184
1185	for (int i = 0; i < maxPulseOrder; i++){
1186	cgpDistributions[i]->Divide(4,2);
1187
1188	cgpDistributions[i]->cd(1);
1189	hBsl[i].Draw();
1190	cgpDistributions[i]->cd(2);
1191	hAsymptote[i].Draw();
1192	cgpDistributions[i]->cd(3);
1193	hT0[i].Draw();
1194	cgpDistributions[i]->cd(4);
1195	hTau1[i].Draw();
1196	cgpDistributions[i]->cd(5);
1197	hTau2[i].Draw();
1198	cgpDistributions[i]->cd(6);
1199	hIntegral[i].Draw();
1200	cgpDistributions[i]->cd(7);
1201	hAmplitude[i].Draw();
1202	cgpDistributions[i]->cd(8);
1203	hRiseTime[i].Draw();
1204
1205
1206
1207	}
1208
1209
1210
1211
1212
1213	}
1214
1215
1216
1217
1218
1219
1220
1221	//-------------------------------------
1222	// Save All Pixel Histograms
1223	//-------------------------------------
1224	if ( saveResults )
1225	{
1226	wholeCamera->SavePixelHistograms(
1227	outFile,
1228	saveResults
1229	);
1230	}
1231
1232	if (ProduceGraphic)
1233	{
1234	UpdateCanvases(
1235	verbosityLevel,
1236	MAX_PULS_ORDER,
1237	false
1238	);
1239
1240	//Process gui events asynchronously during input
1241	cout << endl;
1242	TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
1243	timer.TurnOn();
1244	TString input = Getline("press <return> to exit: ");
1245	timer.TurnOff();
1246	cout << input ;
1247	}
1248
1249	//-------------------------------------
1250	// Delete Objects on Heap
1251	//-------------------------------------
1252	delete wholeCamera;
1253	if (ProduceGraphic)
1254	{
1255	DeletePixelCanvases( maxPulseOrder ,verbosityLevel );
1256	}
1257	delete inputRootFile;
1258	delete outputRootFile;
1259
1260	return( 0 );
1261	}
1262	//----------------------------------------------------------------------------
1263	// end of main function
1264	//-----------------------------------------------------------------------------
1265
1266
1267
1268
1269	//-----------------------------------------------------------------------------
1270	// Funktions
1271	//-----------------------------------------------------------------------------
1272	void
1273	DeletePixelCanvases(
1274	int maxPulseOrder,
1275	int verbosityLevel
1276	)
1277	{
1278	if (verbosityLevel > 2)
1279	{
1280	cout << endl
1281	<< "...delete pixel Canvases"
1282	<< endl;
1283	}
1284
1285	for ( int pulse_order = 0; pulse_order < maxPulseOrder; pulse_order++ )
1286	{
1287	delete cgpPixelPulses[pulse_order];
1288	cgpPixelPulses[pulse_order] = NULL;
1289
1290	delete cgpDistributions[pulse_order];
1291	cgpDistributions[pulse_order] = NULL;
1292	}
1293
1294	delete[] cgpPixelPulses;
1295	cgpPixelPulses = NULL;
1296
1297	delete[] cgpDistributions;
1298	cgpDistributions = NULL;
1299	}
1300
1301	void
1302	UpdateCanvases(
1303	int verbosityLevel,
1304	int max_pulse_order,
1305	bool testmode
1306	)
1307	{
1308	if (verbosityLevel > 3)
1309	{
1310	cout << endl << "...updating canvases" << endl;
1311	}
1312	for (int pulse_order = 0; pulse_order < max_pulse_order; pulse_order++)
1313	{
1314	cgpPixelPulses[pulse_order]->Modified();
1315	cgpPixelPulses[pulse_order]->Update();
1316
1317	cgpDistributions[pulse_order]->Modified();
1318	cgpDistributions[pulse_order]->Update();
1319	// if ( testmode )
1320	// {
1321	// cgpTestHistos->Modified();
1322	// cgpTestHistos->Update();
1323	// }
1324	}
1325	}
1326	//-----------------------------------------------------------------------------
1327	// MAIN Funktion for C compilation
1328	//-----------------------------------------------------------------------------
1329
1330	int main(int argc,char *argv[])
1331	{
1332	// Int_t index;
1333	TString test;
1334	TString rcFileName;
1335	TString processType = "template";
1336	bool rcFileNameCmdSet = false;
1337	int verbLevel = 0; // different verbosity levels can be implemented here
1338	bool verbLevelCmdSet = false;
1339	bool save = false;
1340	bool produceGraphic = false;
1341
1342	// TString inputRootFile = "test.root";
1343	// TString inputPath = "";
1344	// TString outputRootFile = "test.root";
1345	// TString outPutPath = "";
1346	// int firstPixel = 0;
1347	// int nPixel = -1;
1348	// int pixelSetSize = 40;
1349	// int maxOrder = 3;
1350	// bool dbgPixel = false;
1351	// bool fitdata = false;
1352	// bool printStats = false;
1353
1354	// decode arguments
1355	if(argc < 2)
1356	{
1357	printf("no arguements given, using standard arguments\n");
1358	}
1359
1360	// set conditions for functions arguments
1361	for (int i=1;i<argc;i++)
1362	{
1363	test = argv[i];
1364
1365	if (test.Contains("--config-file") \|\| test.Contains("-c"))
1366	{
1367	cout << "RC-File: \"" << argv[i + 1] << "\"" << endl;
1368	rcFileName = argv[i + 1];
1369	rcFileNameCmdSet = true;
1370	continue;
1371	}
1372
1373	if (test.Contains("--verbosity") \|\| test.Contains("-v"))
1374	{
1375	cout << "Verbosity Level: \"" << argv[i + 1] << "\"" << endl;
1376	verbLevel = atoi(argv[i + 1]);
1377	continue;
1378	}
1379
1380	if (test.Contains("--save") \|\| test.Contains("-s"))
1381	{
1382	cout << "will save results" << endl;
1383	save = true;
1384	continue;
1385	}
1386
1387	if (test.Contains("--graphics") \|\| test.Contains("-g"))
1388	{
1389	cout << "will produce graphics" << endl;
1390	produceGraphic = true;
1391	continue;
1392	}
1393	}
1394
1395	// reading rc-File:
1396	if (rcFileNameCmdSet)
1397	{
1398	configfile rcfile( rcFileName, processType );
1399
1400	if (save)
1401	{
1402	rcfile.mSave = save;
1403	}
1404
1405	if (verbLevelCmdSet)
1406	{
1407	rcfile.mVerbLevel = verbLevel;
1408	}
1409
1410	if (produceGraphic)
1411	{
1412	rcfile.mProduceGraphic = produceGraphic;
1413	}
1414
1415	//
1416	FCalcPulseTemplate(
1417	rcfile.mInputFile,
1418	rcfile.mInputPath,
1419	rcfile.mOutputFile,
1420	rcfile.mOutputPath,
1421	rcfile.mFirstPixel,
1422	rcfile.mNumPixel,
1423	rcfile.mPixelSetSize,
1424	rcfile.mMaxOrder,
1425	// rcfile.mHistoOptions,
1426	rcfile.mProduceGraphic,
1427	rcfile.mPrintStats,
1428	rcfile.mSave,
1429	rcfile.mDbgPixel,
1430	rcfile.mVerbLevel
1431	);
1432	}
1433	else
1434	{
1435	cout << "user: check if configfile is set correctly in cmd" << endl;
1436	}
1437
1438	return 0;
1439	}
1440
1441

Note: See TracBrowser for help on using the repository browser.

Download in other formats: