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

Last change on this file since 14943 was 14943, checked in by Jens Buss, 12 years ago
disttribution histogramms added to hist list and saved to root file
File size: 43.2 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_model2.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 = 6,
116 // ------------{ analysis parameters }------------
117 float meanGain = 11.,
118 float meanBsl = -2.5,
119 int bootstrapIt = 10,
120 int pulseModell = 0,
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 // Calculate Max Prop. Value of each slice
550 //-------------------------------------
551
552// pixel[pixelID]->Normalize2Dhistos(pulse_order);
553
554 //from Maximum Overlay
555 if (verbosityLevel > 2)
556 {
557 cout << "...extracting templates from Maximum Overlay "
558 << endl;
559 }
560
561 ExtractPulseTemplate(
562 pixel[pixelID],
563 "Maximum",
564 pulse_order,
565 10,
566 verbosityLevel
567 );
568
569 //from Edge Overlay
570 if (verbosityLevel > 2)
571 {
572 cout << "...extracting templates from Edge Overlay "
573 << endl;
574 }
575 ExtractPulseTemplate(
576 pixel[pixelID],
577 "Edge",
578 pulse_order,
579 10,
580 verbosityLevel
581 );
582
583 pixel[pixelID]->SetRangeUser(pixelOverlayXaxisLeft, pixelOverlayXaxisRight, pulse_order);
584
585
586 PixelCsv.WritePointSet(
587 pixel[pixelID],
588 "Maximum",
589 pulse_order
590 );
591
592 PixelCsv.WritePointSet(
593 pixel[pixelID],
594 "Edge",
595 pulse_order
596 );
597
598 if ( saveResults )
599 {
600 WritePixelTemplateToCsv(
601 pixel[pixelID],
602 OutPutPath,
603 "Maximum",
604 pulse_order,
605 verbosityLevel
606 );
607
608 WritePixelTemplateToCsv(
609 pixel[pixelID],
610 OutPutPath,
611 "Edge",
612 pulse_order,
613 verbosityLevel
614 );
615 }
616
617 if (ProduceGraphic && debugPixel)
618 {
619
620 pixel[pixelID]->DrawTemplateHistograms(
621 cgpPixelPulses,
622 PixelCanvasFrameNrs
623 );
624
625 pixel[pixelID]->DrawEdgeTemplateHistograms(
626 cgpPixelPulses,
627 PixelCanvasFrameNrs
628 );
629
630// cgpPixelPulses[pulse_order]->cd(6);
631
632
633 //-------------------------------------
634 //-------------------------------------
635 // Test Area
636 //-------------------------------------
637 //-------------------------------------
638
639 cgpPixelPulses[pulse_order]->cd(8);
640
641 //-------------------------------------
642 //-------------------------------------
643 // EOF Test Area
644 //-------------------------------------
645 //-------------------------------------
646 }
647
648 // =======================================================
649 // Fit of pulse function to template
650
651 Pulse* pulseFits[6];
652
653// pulses.cd(1);
654 Pulse maxMaxPulse("maxMaxPulse", pixel[pixelID]->hPixelMax[pulse_order], fitOption, pulseModell,pulse_order);
655 pulseFits[0] = &maxMaxPulse;
656
657// pulses.cd(2);
658 Pulse maxMedianPulse("maxMedianPulse", pixel[pixelID]->hPixelMedian[pulse_order], fitOption,pulseModell,pulse_order);
659 pulseFits[1] = &maxMedianPulse;
660
661// pulses.cd(3);
662 Pulse maxMeanPulse("maxMeanPulse", pixel[pixelID]->hPixelMean[pulse_order], fitOption,pulseModell,pulse_order);
663 pulseFits[2] = &maxMeanPulse;
664
665 Pulse edgeMaxPulse("edgeMaxPulse", pixel[pixelID]->hPixelEdgeMax[pulse_order], fitOption,pulseModell,pulse_order);
666 pulseFits[3] = &edgeMaxPulse;
667
668 Pulse edgeMedianPulse("edgeMedianPulse", pixel[pixelID]->hPixelEdgeMedian[pulse_order], fitOption,pulseModell,pulse_order);
669 pulseFits[4] = &edgeMedianPulse;
670
671 Pulse edgeMeanPulse("edgeMeanPulse", pixel[pixelID]->hPixelEdgeMean[pulse_order], fitOption,pulseModell,pulse_order);
672 pulseFits[5] = &edgeMeanPulse;
673
674
675 //--------------------------------------------------------
676 //Write parameters from pulse fit to CSV
677
678 for (int i = 0; i < 6; i++)
679 {
680 PixelModelCsv.WritePulseAttributes(
681 pixel[pixelID],
682 pulseFits[i],
683 "Edge_pix",
684 pulse_order
685 );
686 }
687 //-------------------------------------
688 // Fill Error Histogramms
689 //-------------------------------------
690
691 for (int i = pixel[pixelID]->hPixelEdgeMedian[pulse_order]->GetXaxis()->GetFirst();
692 i <= pixel[pixelID]->hPixelEdgeMedian[pulse_order]->GetXaxis()->GetLast();
693 i++)
694 {
695
696 hMedianErrors[pulse_order].SetBinContent(
697 i,
698 pixel[pixelID]->hPixelEdgeMedian[pulse_order]->GetBinError(i)
699 );
700
701 hMeanErrors[pulse_order].SetBinContent(
702 i,
703 pixel[pixelID]->hPixelEdgeMean[pulse_order]->GetBinError(i)
704 );
705
706 hMaxErrors[pulse_order].SetBinContent(
707 i,
708 pixel[pixelID]->hPixelEdgeMax[pulse_order]->GetBinError(i)
709 );
710 }
711
712 if (ProduceGraphic && debugPixel)
713 {
714 int current_order = 0;
715 pulses.cd(1);
716 hMaxErrors[current_order].Draw();
717
718 pulses.cd(2);
719 hMedianErrors[current_order].Draw();
720
721 pulses.cd(3);
722 hMeanErrors[current_order].Draw();
723 pulses.Update();
724 }
725
726 //-------------------------------------
727 // Fill Distribution Histogramms
728 //-------------------------------------
729
730 hBsl[pulse_order].Fill( edgeMedianPulse.GetBsl() );
731 hAsymptote[pulse_order].Fill( edgeMedianPulse.GetHeight() );
732 hT0[pulse_order].Fill( edgeMedianPulse.GetT0() );
733 hTau1[pulse_order].Fill( edgeMedianPulse.GetTau1() );
734 hTau2[pulse_order].Fill( edgeMedianPulse.GetTau2() );
735 hIntegral[pulse_order].Fill( edgeMedianPulse.GetIntegral());
736 hAmplitude[pulse_order].Fill( edgeMedianPulse.GetAmplitude());
737 hRiseTime[pulse_order].Fill( edgeMedianPulse.GetRiseTime());
738
739 //-------------------------------------
740 // Fill Distribution Graphs
741 //-------------------------------------
742
743 grPheVsAmpl.Set(grPheVsAmpl.GetN()+1);
744 grPheVsIntegral.Set(grPheVsIntegral.GetN()+1);
745
746
747 grPheVsAmpl.SetPoint(
748 grPheVsAmpl.GetN()-1,
749 edgeMedianPulse.GetPhE(),
750 edgeMedianPulse.GetAmplitude()
751 );
752
753 grPheVsAmpl.SetPointError(
754 grPheVsAmpl.GetN()-1,
755 edgeMedianPulse.GetPhEErr(),
756 edgeMedianPulse.GetAmplitudeErr()
757 );
758
759 grPheVsIntegral.SetPoint(
760 grPheVsIntegral.GetN()-1,
761 edgeMedianPulse.GetPhE(),
762 edgeMedianPulse.GetIntegral()
763 );
764
765 grPheVsAmpl.SetPointError(
766 grPheVsIntegral.GetN()-1,
767 edgeMedianPulse.GetPhEErr(),
768 edgeMedianPulse.GetIntegralErr()
769 );
770
771 grPheVsTau2.SetPoint(
772 grPheVsIntegral.GetN()-1,
773 edgeMedianPulse.GetPhE(),
774 edgeMedianPulse.GetTau2()
775 );
776
777 grPheVsTau2.SetPointError(
778 grPheVsIntegral.GetN()-1,
779 edgeMedianPulse.GetPhE(),
780 edgeMedianPulse.GetTau2Err()
781 );
782
783// cgpDistributions[pulse_order].cd(1);
784
785
786// FitMaxPropabilityPuls(
787// pixel[pixelID]->hPixelEdgeMean[pulse_order],
788// verbosityLevel
789// );
790
791 //Preparing Camera
792 if (first_pass)
793 {
794 if (verbosityLevel > 0)
795 {
796 cout << endl << "...preparing camera" << endl;
797 }
798
799 wholeCamera = new PixelSum(
800 "AllPixel",
801 1440,
802 maxPulseOrder,
803 verbosityLevel,
804 stats,
805 "C",
806 pixel[pixelID]->mPixelOverlayXaxisLeft,
807 pixel[pixelID]->mPixelOverlayXaxisRight ,
808 pixel[pixelID]->mBSLMean ,
809 pixel[pixelID]->mGainMean ,
810 outputRootFile
811 );
812
813 first_pass = false;
814 }
815
816 //-------------------------------------
817 // Fill Histogramms of Camera
818 //-------------------------------------
819
820 wholeCamera->hMaxOverlay[pulse_order]->Add(
821 pixel[pixelID]->hMaxOverlay[pulse_order]
822 );
823
824 wholeCamera->hMaxProfile[pulse_order]->Add(
825 pixel[pixelID]->hMaxProfile[pulse_order]
826 );
827
828 wholeCamera->hEdgeOverlay[pulse_order]->Add(
829 pixel[pixelID]->hEdgeOverlay[pulse_order]
830 );
831
832 wholeCamera->hEdgeProfile[pulse_order]->Add(
833 pixel[pixelID]->hEdgeProfile[pulse_order]
834 );
835
836 //-------------------------------------
837 // Comparisons
838 //-------------------------------------
839
840 //chi2 test
841// float chi2 =
842// pixel[pixelID]->hPixelEdgeMean[pulse_order]->Chi2Test(
843// pixel[pixelID]->hPixelMean[pulse_order],
844// "UUPCHI2"
845// );
846// cout << "p-Value :" << chi2 << endl;
847// wholeCamera->hChi2EdgetoMax[pulse_order]->Fill(chi2);
848// cgpDistributions[pulse_order]->cd();
849// wholeCamera->hChi2EdgetoMax[pulse_order]->Draw();
850
851 if (verbosityLevel > 2)
852 {
853 cout << endl
854 << "...End of pulseorder "
855 << pulse_order
856 << endl;
857 }
858
859 }
860 // End of Loop over pulsorder of current Pixel
861
862 if (ProduceGraphic && debugPixel)
863 {
864 UpdateCanvases(
865 verbosityLevel,
866 MAX_PULS_ORDER,
867 false
868 );
869 }
870
871 if (debugPixel)
872 {
873 //Process gui events asynchronously during input
874 cout << endl;
875 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
876 timer.TurnOn();
877 TString input = Getline("Type 'q' to exit, Type 's' to exit set, <return> to go on: ");
878 timer.TurnOff();
879 if (input=="q\n")
880 {
881 return(0);
882 }
883 if (input=="s\n")
884 {
885 break;
886 }
887 }
888
889 if ( saveResults )
890 {
891 pixel[pixelID]->SavePixelHistograms(
892 outFile,
893 saveResults
894 );
895 }
896
897 //deleteCurrent Pixel from Heap
898 delete pixel[pixelID];
899 pixel[pixelID] = NULL;
900
901 if (verbosityLevel > 2)
902 {
903 cout << endl
904 << "...End of Pixel"
905 << endl
906 << "------------------------------------------------"
907 << endl;
908 }
909 }
910 // End of Loop over all Pixels of set
911
912 if (verbosityLevel > 1)
913 {
914 cout << endl
915 << "...End of Loop over all Pixels of set"
916 << endl
917 << "------------------------------------------------"
918 << endl;
919 }
920 }
921 // End of Loop over all Pixelsets
922
923 if (verbosityLevel > 0)
924 {
925 cout << endl
926 << "...End of Loop over all Pixelsets"
927 << endl
928 << "------------------------------------------------"
929 << endl;
930 }
931
932 delete[] pixel;
933 pixel = NULL;
934
935//-------------------------------------
936// Draw All Pixel Histograms
937//-------------------------------------
938
939 for ( int pulse_order = 0 ;
940 pulse_order < wholeCamera->mMaxPulseOrder ;
941 pulse_order ++)
942 {
943 if (verbosityLevel > 2)
944 {
945 cout << "-------------------------------------"
946 << endl
947 << "...processing Pulse-Order: "
948 << pulse_order;
949 }
950
951// Calculate Max Prop. Value of each slice
952//-------------------------------------
953
954 //from Maximum Overlay
955 ExtractPulseTemplate(
956 wholeCamera,
957 "Maximum",
958 pulse_order,
959 bootstrapIt,
960 verbosityLevel
961 );
962 if ( saveResults )
963 {
964 WritePixelTemplateToCsv(
965 wholeCamera,
966 OutPutPath,
967 "Maximum",
968 pulse_order,
969 verbosityLevel
970 );
971 }
972
973 //from Edge Overlay
974 ExtractPulseTemplate(
975 wholeCamera,
976 "Edge",
977 pulse_order,
978 bootstrapIt,
979 verbosityLevel
980 );
981 if ( saveResults )
982 {
983 WritePixelTemplateToCsv(
984 wholeCamera,
985 OutPutPath,
986 "Edge",
987 pulse_order,
988 verbosityLevel
989 );
990 }
991
992 if (ProduceGraphic)
993 {
994 wholeCamera->DrawTemplateHistograms(
995 cgpPixelPulses,
996 PixelCanvasFrameNrs
997 );
998
999 wholeCamera->DrawEdgeTemplateHistograms(
1000 cgpPixelPulses,
1001 PixelCanvasFrameNrs
1002 );
1003 }
1004
1005 ShiftStartOfHistoInXTo(
1006 wholeCamera->hEdgeOverlay[0],
1007 0
1008 );
1009
1010 ShiftStartOfHistoInXTo(
1011 wholeCamera->hMaxOverlay[0],
1012 0
1013 );
1014
1015 // =======================================================
1016 // Fit of pulse function to template
1017
1018 Pulse* allPulseFits[6];
1019
1020 pulses.cd(1);
1021 Pulse maxCamMaxPulse("maxMaxPulse", wholeCamera->hPixelMax[pulse_order], "Q", pulseModell);
1022 allPulseFits[0] = &maxCamMaxPulse;
1023
1024 pulses.cd(2);
1025 Pulse maxCamMedianPulse("maxMedianPulse", wholeCamera->hPixelMedian[pulse_order], "Q", pulseModell);
1026 allPulseFits[1] = &maxCamMedianPulse;
1027
1028 pulses.cd(3);
1029 Pulse maxCamMeanPulse("maxMeanPulse", wholeCamera->hPixelMean[pulse_order],"Q", pulseModell);
1030 allPulseFits[2] = &maxCamMeanPulse;
1031
1032 Pulse edgeCamMaxPulse("edgeMaxPulse", wholeCamera->hPixelEdgeMax[pulse_order], "Q", pulseModell);
1033 allPulseFits[3] = &edgeCamMaxPulse;
1034
1035 Pulse edgeCamMedianPulse("edgeMedianPulse", wholeCamera->hPixelEdgeMedian[pulse_order],"Q", pulseModell);
1036 allPulseFits[4] = &edgeCamMedianPulse;
1037
1038 Pulse edgeCamMeanPulse("edgeMeanPulse", wholeCamera->hPixelEdgeMean[pulse_order], "Q", pulseModell);
1039 allPulseFits[5] = &edgeCamMeanPulse;
1040
1041
1042 //--------------------------------------------------------
1043 //Write parameters from pulse fit to CSV
1044
1045 for (int i = 0; i < 6; i++)
1046 {
1047 PixelModelCsv.WritePulseAttributes(
1048 wholeCamera,
1049 allPulseFits[i],
1050 "Edge",
1051 pulse_order
1052 );
1053 }
1054 } //EOF: Draw All Pixel Histograms
1055
1056
1057 //-------------------------------------
1058 // Draw Parameter Distributions
1059 //-------------------------------------
1060
1061 if (ProduceGraphic)
1062 {
1063 cgpGraphs->cd(1);
1064
1065 if (verbosityLevel > 0) cout << "...drawing grPheVsAmpl" << endl;
1066 grPheVsAmpl.Draw("AP");
1067
1068 cgpGraphs->cd(2);
1069 if (verbosityLevel > 0) cout << "...drawing grPheVsIntegral" << endl;
1070 grPheVsIntegral.Draw("AP");
1071
1072 cgpGraphs->cd(3);
1073 if (verbosityLevel > 0) cout << "...drawing grPheVsTau2" << endl;
1074 grPheVsTau2.Draw("AP");
1075
1076 cgpGraphs->Update();
1077
1078 for (int i = 0; i < maxPulseOrder; i++){
1079 cgpDistributions[i]->Divide(4,2);
1080
1081 cgpDistributions[i]->cd(1);
1082 hBsl[i].Draw();
1083 cgpDistributions[i]->cd(2);
1084 hAsymptote[i].Draw();
1085 cgpDistributions[i]->cd(3);
1086 hT0[i].Draw();
1087 cgpDistributions[i]->cd(4);
1088 hTau1[i].Draw();
1089 cgpDistributions[i]->cd(5);
1090 hTau2[i].Draw();
1091 cgpDistributions[i]->cd(6);
1092 hIntegral[i].Draw();
1093 cgpDistributions[i]->cd(7);
1094 hAmplitude[i].Draw();
1095 cgpDistributions[i]->cd(8);
1096 hRiseTime[i].Draw();
1097
1098
1099
1100 }
1101
1102
1103
1104
1105
1106 }
1107
1108
1109
1110
1111
1112
1113
1114//-------------------------------------
1115// Save All Pixel Histograms
1116//-------------------------------------
1117 if ( saveResults )
1118 {
1119 wholeCamera->SavePixelHistograms(
1120 outFile,
1121 saveResults
1122 );
1123
1124 SaveList(
1125 outFile,
1126 "",
1127 histoList,
1128 saveResults,
1129 verbosityLevel
1130 );
1131
1132
1133 }
1134
1135 if (ProduceGraphic)
1136 {
1137 UpdateCanvases(
1138 verbosityLevel,
1139 MAX_PULS_ORDER,
1140 false
1141 );
1142
1143 //Process gui events asynchronously during input
1144 cout << endl;
1145 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
1146 timer.TurnOn();
1147 TString input = Getline("press <return> to exit: ");
1148 timer.TurnOff();
1149 cout << input ;
1150 }
1151
1152//-------------------------------------
1153// Delete Objects on Heap
1154//-------------------------------------
1155 delete wholeCamera;
1156 if (ProduceGraphic)
1157 {
1158 DeletePixelCanvases( maxPulseOrder ,verbosityLevel );
1159 }
1160 delete inputRootFile;
1161 delete outputRootFile;
1162 delete histoList;
1163
1164 return( 0 );
1165}
1166//----------------------------------------------------------------------------
1167// end of main function
1168//-----------------------------------------------------------------------------
1169
1170
1171
1172
1173//-----------------------------------------------------------------------------
1174// Funktions
1175//-----------------------------------------------------------------------------
1176void
1177DeletePixelCanvases(
1178 int maxPulseOrder,
1179 int verbosityLevel
1180 )
1181{
1182 if (verbosityLevel > 2)
1183 {
1184 cout << endl
1185 << "...delete pixel Canvases"
1186 << endl;
1187 }
1188
1189 for ( int pulse_order = 0; pulse_order < maxPulseOrder; pulse_order++ )
1190 {
1191 delete cgpPixelPulses[pulse_order];
1192 cgpPixelPulses[pulse_order] = NULL;
1193
1194 delete cgpDistributions[pulse_order];
1195 cgpDistributions[pulse_order] = NULL;
1196 }
1197
1198 delete[] cgpPixelPulses;
1199 cgpPixelPulses = NULL;
1200
1201 delete[] cgpDistributions;
1202 cgpDistributions = NULL;
1203}
1204
1205void
1206UpdateCanvases(
1207 int verbosityLevel,
1208 int max_pulse_order,
1209 bool testmode
1210 )
1211{
1212 if (verbosityLevel > 3)
1213 {
1214 cout << endl << "...updating canvases" << endl;
1215 }
1216 for (int pulse_order = 0; pulse_order < max_pulse_order; pulse_order++)
1217 {
1218 cgpPixelPulses[pulse_order]->Modified();
1219 cgpPixelPulses[pulse_order]->Update();
1220
1221 cgpDistributions[pulse_order]->Modified();
1222 cgpDistributions[pulse_order]->Update();
1223// if ( testmode )
1224// {
1225// cgpTestHistos->Modified();
1226// cgpTestHistos->Update();
1227// }
1228 }
1229}
1230//-----------------------------------------------------------------------------
1231// MAIN Funktion for C compilation
1232//-----------------------------------------------------------------------------
1233
1234int main(int argc,char *argv[])
1235{
1236// Int_t index;
1237 TString test;
1238 TString rcFileName;
1239 TString processType = "template";
1240 bool rcFileNameCmdSet = false;
1241 int verbLevel = 0; // different verbosity levels can be implemented here
1242 bool verbLevelCmdSet = false;
1243 bool save = false;
1244 bool produceGraphic = false;
1245
1246// TString inputRootFile = "test.root";
1247// TString inputPath = "";
1248// TString outputRootFile = "test.root";
1249// TString outPutPath = "";
1250// int firstPixel = 0;
1251// int nPixel = -1;
1252// int pixelSetSize = 40;
1253// int maxOrder = 3;
1254// bool dbgPixel = false;
1255// bool fitdata = false;
1256// bool printStats = false;
1257
1258 // decode arguments
1259 if(argc < 2)
1260 {
1261 printf("no arguements given, using standard arguments\n");
1262 }
1263
1264 // set conditions for functions arguments
1265 for (int i=1;i<argc;i++)
1266 {
1267 test = argv[i];
1268
1269 if (test.Contains("--config-file") || test.Contains("-c"))
1270 {
1271 cout << "RC-File: \"" << argv[i + 1] << "\"" << endl;
1272 rcFileName = argv[i + 1];
1273 rcFileNameCmdSet = true;
1274 continue;
1275 }
1276
1277 if (test.Contains("--verbosity") || test.Contains("-v"))
1278 {
1279 cout << "Verbosity Level: \"" << argv[i + 1] << "\"" << endl;
1280 verbLevel = atoi(argv[i + 1]);
1281 continue;
1282 }
1283
1284 if (test.Contains("--save") || test.Contains("-s"))
1285 {
1286 cout << "will save results" << endl;
1287 save = true;
1288 continue;
1289 }
1290
1291 if (test.Contains("--graphics") || test.Contains("-g"))
1292 {
1293 cout << "will produce graphics" << endl;
1294 produceGraphic = true;
1295 continue;
1296 }
1297 }
1298
1299 // reading rc-File:
1300 if (rcFileNameCmdSet)
1301 {
1302 configfile rcfile( rcFileName, processType );
1303
1304 if (save)
1305 {
1306 rcfile.mSave = save;
1307 }
1308
1309 if (verbLevelCmdSet)
1310 {
1311 rcfile.mVerbLevel = verbLevel;
1312 }
1313
1314 if (produceGraphic)
1315 {
1316 rcfile.mProduceGraphic = produceGraphic;
1317 }
1318
1319 //
1320 FCalcPulseTemplate(
1321 rcfile.mInputFile,
1322 rcfile.mInputPath,
1323 rcfile.mOutputFile,
1324 rcfile.mOutputPath,
1325 rcfile.mFirstPixel,
1326 rcfile.mNumPixel,
1327 rcfile.mPixelSetSize,
1328 rcfile.mMaxOrder,
1329// rcfile.mHistoOptions,
1330 rcfile.mProduceGraphic,
1331 rcfile.mPrintStats,
1332 rcfile.mSave,
1333 rcfile.mDbgPixel,
1334 rcfile.mVerbLevel
1335 );
1336 }
1337 else
1338 {
1339 cout << "user: check if configfile is set correctly in cmd" << endl;
1340 }
1341
1342 return 0;
1343}
1344
1345
Note: See TracBrowser for help on using the repository browser.