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

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