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

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