source: fact/tools/rootmacros/PulseTemplates/FPulseOverlay.C@ 16858

Last change on this file since 16858 was 15033, checked in by Jens Buss, 12 years ago
uncommented save histograms
File size: 43.7 KB
Line 
1/********************** FPulseOverlay ***********************
2* read FACT raw data
3* remove spikes
4* calculate baseline
5* find peaks (CFD and normal discriminator)
6+ search for Peaks in data
7+ put peaks into different histos depending on Amplitude
8+ draw histos for single, double, tripple, ....Photonpulses
9+ draw Parameterdevelopment
10 so it can be used for detection of other peaks
11*************************************************************/
12//----------------------------------------------------------------------------
13// root libraries
14//----------------------------------------------------------------------------
15
16#include <TROOT.h>
17#include <TCanvas.h>
18#include <TProfile.h>
19#include <TTimer.h>
20#include <TH1F.h>
21#include <TH2F.h>
22#include <Getline.h>
23#include <TLine.h>
24#include <TBox.h>
25#include <TMath.h>
26#include <TFile.h>
27#include <TStyle.h>
28#include <TString.h>
29#include <TObjArray.h>
30#include <TF1.h>
31#include<TSystem.h>
32
33#include <stdio.h>
34#include <stdint.h>
35#include <cstdio>
36#include <iostream>
37#include <fstream>
38using namespace std;
39
40#define NPIX 1440
41#define FAD_MAX_SAMPLES 1024
42#define HAVE_ZLIB
43
44//----------------------------------------------------------------------------
45// rootmacros
46//----------------------------------------------------------------------------
47
48//#include "../fits.h"
49
50#include "../openFits.h"
51//#include "../openFits.c"
52
53#include "../discriminator.h"
54//#include "../discriminator.C"
55#include "../zerosearch.h"
56//#include "../zerosearch.C"
57# include "../factfir.h"
58//#include "../factfir.C"
59
60//#include "../DrsCalibration.C"
61#include "../DrsCalibration.h"
62
63#include "../SpikeRemoval.h"
64//#include "../SpikeRemoval.C"
65
66#include "rootfilehandler.h"
67//#include "rootfilehandler.C"
68
69#include "pixel.h"
70//#include "pixel.C"
71
72#include "configfile.h"
73
74
75//----------------------------------------------------------------------------
76// Decleration of global variables
77//----------------------------------------------------------------------------
78
79bool breakout=false;
80
81vector<int16_t> AllPixelDataVector;
82vector<int16_t> StartCellVector;
83unsigned int CurrentEventID;
84size_t PXLxROI;
85UInt_t gRegionOfInterest;
86UInt_t NumberOfPixels;
87TString histotitle;
88
89size_t TriggerOffsetROI, RC;
90size_t drs_n;
91vector<float> drs_basemean;
92vector<float> drs_gainmean;
93vector<float> drs_triggeroffsetmean;
94
95vector<float> Ameas(FAD_MAX_SAMPLES); // copy of the data (measured amplitude
96vector<float> Vcorr(FAD_MAX_SAMPLES); // corrected Values
97vector<float> Vslide(FAD_MAX_SAMPLES); // sliding average result
98vector<float> Vcfd(FAD_MAX_SAMPLES); // CDF result
99vector<float> Vcfd2(FAD_MAX_SAMPLES); // CDF result + 2nd sliding average
100
101int gNEvents=0;
102float gGainMean = 9;
103float gBSLMean = 0;
104
105//typedef struct{
106// double maxAmpl;
107// int countOfMax;
108// } OverlayMaximum;
109
110//typedef struct{
111// TString name;
112// TString title;
113// TString xTitle;
114// TString yTitle;
115// float yMax;
116// float yMin;
117//} histAttrib_t;
118
119//histAttrib_t* gHistoAttrib[maxPulseOrder];
120//histAttrib_t* gProfileAttrib[maxPulseOrder];
121
122
123// histograms
124const int Number_Of_Debug_Histo_Types = 7;
125
126const unsigned int
127 Ameas_ = 0,
128 N1mean_ = 1,
129 Vcorr_ = 2,
130 Vtest_ = 3,
131 Vslide_ = 4,
132 Vcfd_ = 5,
133 Vcfd2_ = 6;
134
135//const char* gHistoTypes[8] = {
136// "hMaxOverlay",
137// "hPixelMax",
138// "hEdgeOverlay",
139// "hMaxProfile",
140// "hAllPixelOverlay",
141// "hAllPixelMax",
142// "hAllPixelMaxGaus",
143// "hAllPixelProfile"
144//};
145
146//----------------------------------------------------------------------------
147// Initialisation of histograms
148//----------------------------------------------------------------------------
149
150 // Temporary Objects
151 TH1F* debugHistos = NULL;
152 //TH2F* hOverlayTemp = NULL;
153 //TH1D* hProjPeak = NULL;
154 TH1F* hTesthisto = NULL;
155 TH2F* hTesthisto2 = NULL;
156
157 //Histogram Parameters
158 Int_t gPixelOverlayXaxisLeft = 0;
159 Int_t gPixelOverlayXaxisRight = 0;
160
161 //Root-File Objects
162// TObjArray* hList[pixelSetSize] = {NULL};
163 TList* hRootList = NULL;
164
165 TCanvas** cgpPixelPulses = NULL;
166 TCanvas** cgpDistributions = NULL;
167 TCanvas* cFiltered = NULL;
168 TCanvas* cgpTestHistos = NULL;
169
170//----------------------------------------------------------------------------
171// Functions
172//----------------------------------------------------------------------------
173
174void BookDebugHistos(int );
175void BookTestHistos( int );
176
177//void DeletePixelHistos(int );
178//void SaveHistograms( const char*, const char*, TObjArray*, int );
179
180void FillHistograms(Pixel*, vector<Region>*, int, int, int , int );
181void DrawTestHistograms( int);
182bool ProduceDebugHistos( vector<Region> *pZXings);
183bool UseThisPulse( int, int, int, float, int, int);
184void UpdateCanvases( int, int, bool);
185void DeletePixelCanvases( int, int );
186
187TString gDataDirectory;
188TString gHomeDirectory;
189TString ghostname;
190
191//----------------------------------------------------------------------------
192//----------------------------------------------------------------------------
193// MAIN - Funtion
194//----------------------------------------------------------------------------
195//----------------------------------------------------------------------------
196int FPulseOverlay(
197 TString inputPath = "raw/2012/03/09/",
198 TString datafilename = "20120309_018.fits.gz",
199 TString drsfilename = "20120309_012.drs.fits.gz",
200 TString OutputPath = "analysis/FPulseTemplate/20120309_018/Overlay/",
201 TString OutRootFileName = "20120309_018.root",
202 int firstevent = 1, // Hast to be between 1 and infinity NOT 0!!!!
203 int nevents = -1,
204 int firstpixel = 0,
205 int npixel = 2,
206 int pixelSetSize = 150,
207 int maxPulseOrder = 5,
208 int verbosityLevel = 1, // different verbosity levels can be implemented here
209 TString histoOptions = "SRM",
210 int AmplWindowWidth = 16, //Width of Window for selection of pulses to histograms
211 float GainMean = 9,
212 float BSLMean = -1, //4 Histogramms will be drawn, decide how far pulsehights differ from eachother
213 int avg1 = 8,
214 int avg2 = 8,
215 int OverlayWindowLeft = 70,
216 int OverlayWindowRight = 230,
217 bool ProduceGraphic = false,
218 bool spikeDebug = false,
219 bool debugPixel = false,
220 bool stats = false,
221 bool saveResults = true,
222 bool testmode = false,
223 int refresh_rate = 500 //refresh rate for canvases
224 )
225{
226
227//----------------------------------------------------------------------------
228// Save-Root-File Settings
229//----------------------------------------------------------------------------
230
231 inputPath = SetHostsPaths(false, inputPath );
232 OutputPath = SetHostsPaths(true, OutputPath );
233 datafilename.Prepend(inputPath);
234 drsfilename.Prepend(inputPath);
235 OutRootFileName.Prepend(OutputPath);
236// drsfilename = SetHostsPaths(false, drsfilename );
237// OutRootFileName = SetHostsPaths(true, OutRootFileName );
238 cout << endl;
239 cout << "drs filename:\t" << drsfilename << endl;
240 cout << "data filename:\t" << datafilename << endl;
241 cout << "out filename:\t" << OutRootFileName << endl;
242
243
244 if (saveResults)
245 {
246 CreateRootFile( OutRootFileName, true, verbosityLevel );
247 }
248
249//----------------------------------------------------------------------------
250// root global Settings
251//----------------------------------------------------------------------------
252
253 gGainMean = GainMean;
254 gBSLMean = BSLMean;
255 gPixelOverlayXaxisLeft = OverlayWindowLeft;
256 gPixelOverlayXaxisRight = OverlayWindowRight;
257
258 gStyle->SetPalette(1,0);
259 gROOT->SetStyle("Plain");
260// gPad->SetGrid();
261
262//----------------------------------------------------------------------------
263// root Canvas Settings
264//----------------------------------------------------------------------------
265 //Canvas Pad numbering
266 int pulsesCanvasFrameNrs[4] = {
267 1, // Top left
268 2, // Top right
269 3, // bottom left
270 4 // bootom right
271 };
272
273 //Canvas Pad numbering
274 int distributionCanvasFrameNrs[4] = {
275 1, // Top left
276 2, // Top right
277 3, // bottom left
278 4 // bootom right
279 };
280
281 if (ProduceGraphic)
282 {
283
284 //Canvases
285 cgpPixelPulses = new TCanvas*[maxPulseOrder];
286 cgpDistributions = new TCanvas*[maxPulseOrder];
287
288 //TCanvas* gpcDevelopment = NULL;
289 TString cName = "";
290 TString cTitle = "";
291
292 //naming of pulse canvases
293 for (
294 int pulse_order = maxPulseOrder - 1;
295 pulse_order >= 0 ;
296 pulse_order--
297 )
298 {
299 cName ="cgpDistributions";
300 cName += pulse_order;
301 cTitle ="Distributions of Pulses with Order of: ";
302 cTitle += pulse_order;
303 cgpDistributions[pulse_order] = new TCanvas(cName,cTitle, 720,pulse_order*20,720,720);
304 cgpDistributions[pulse_order]->Divide(2, 2);
305 cName ="cgpPixelPulses";
306 cName += pulse_order;
307 cTitle ="Overlays of Pulses with Order of: ";
308 cTitle += pulse_order;
309 cgpPixelPulses[pulse_order] = new TCanvas(cName,cTitle, 0,pulse_order*20,720,720);
310 cgpPixelPulses[pulse_order]->Divide(2, 2);
311 }
312
313
314 // Create (pointer to) Canvases, which are used in every run,
315 // also in 'non-debug' runs
316 // Canvases only need if spike Debug, but I want to deklare
317 // the pointers anyway ...
318 if (spikeDebug)
319 {
320 cFiltered = new TCanvas("cFiltered","filtered DRS Waveforms", 1080,420,360,360);
321 cFiltered->Divide(1, 3);
322 }
323
324 if (testmode)
325 {
326 //additional Test histograms
327 cgpTestHistos = new TCanvas( "cgpTestHistos", "Test Histograms", 360, 420, 360, 360 );
328 cgpTestHistos->Divide(2,0);
329 }
330 }
331
332//-----------------------------------------------------------------------------
333// Filter-Settings
334//-----------------------------------------------------------------------------
335// CFD filter settings
336 int k_cfd = 10;
337 vector<double> a_cfd(k_cfd, 0);
338 double b_cfd = 1.;
339 a_cfd[0]=-0.75;
340 a_cfd[k_cfd-1]=1.;
341
342//-----------------------------------------------------------------------------
343// prepare and read datafile
344//-----------------------------------------------------------------------------
345
346// Open the data file
347
348 fits * datafile = NULL;
349 // Opens the raw data file and 'binds' the variables given as
350 // Parameters to the data file. So they are filled with
351 // raw data as soon as datafile->GetRow(int) is called.
352 gNEvents = openDataFits(
353 datafilename,
354 &datafile,
355 AllPixelDataVector,
356 StartCellVector,
357 CurrentEventID,
358 gRegionOfInterest,
359 NumberOfPixels,
360 PXLxROI,
361 verbosityLevel
362 );
363
364 if (gNEvents == 0)
365 {
366 cout << "return code of OpenDataFile:" << datafilename<< endl;
367 cout << "is zero -> aborting." << endl;
368 return 1;
369 }
370
371 if (verbosityLevel > 0)
372 {
373 cout << endl <<"number of events in file: "<< gNEvents << "\t";
374 }
375
376 if ( nevents == -1 || nevents > gNEvents ) nevents = gNEvents; // -1 means all!
377
378 if (verbosityLevel > 0)
379 {
380 cout <<"of, which "<< nevents << " will be processed"<< endl;
381 cout <<"Total # of Pixel: "<< NumberOfPixels << "\t";
382 }
383
384 if ( npixel == -1 || npixel > (int)NumberOfPixels ) npixel = (int)NumberOfPixels; // -1 means all!
385
386 if (verbosityLevel > 0)
387 {
388 cout <<"of, which "<< npixel << " will be processed"<< endl;
389 }
390
391 //Get the DRS calibration
392 RC = openCalibFits(
393 drsfilename,
394 drs_basemean,
395 drs_gainmean,
396 drs_triggeroffsetmean,
397 TriggerOffsetROI
398 );
399
400 if (RC == 0)
401 {
402 cout << "return code of openCalibFits:" << drsfilename << endl;
403 cout << "is zero -> aborting." << endl;
404 return 1;
405 }
406
407//-----------------------------------------------------------------------------
408// Book the histograms
409//-----------------------------------------------------------------------------
410
411
412
413 if (testmode)
414 BookTestHistos( verbosityLevel );
415
416
417//-----------------------------------------------------------------------------
418// Main Cycle
419//-----------------------------------------------------------------------------
420 if ( pixelSetSize == -1 )
421 {
422 pixelSetSize = firstpixel +npixel;
423 }
424//----------------------------------------------------------------------------
425// Initialize Pixel
426//----------------------------------------------------------------------------
427 Pixel** pixel = new Pixel*[NPIX];
428
429 for (int i = 0 ; i < NPIX; i++)
430 {
431 pixel[i] = NULL;
432 }
433 int lastPixelOfSet = 0;
434//-------------------------------------
435// Loop over Pixel Sets
436//-------------------------------------
437 for ( int firstPixelOfSet = firstpixel;
438 firstPixelOfSet < firstpixel + npixel;
439 firstPixelOfSet += pixelSetSize )
440 {
441 lastPixelOfSet = firstPixelOfSet + pixelSetSize-1;
442 if (lastPixelOfSet > firstpixel + npixel)
443 {
444 lastPixelOfSet = firstpixel + npixel;
445 }
446 if (lastPixelOfSet > 1439)
447 {
448 lastPixelOfSet = 1439;
449 }
450 if (verbosityLevel == 0)
451 {
452 cout << "------------------------------------------------" << endl
453 << "...processing Set from Pixel "
454 << firstPixelOfSet
455 << " to Pixel "
456 << lastPixelOfSet << endl;
457 }
458
459//--------------------------------------------------------------------
460// Loop over every Event of Pixel Set
461//--------------------------------------------------------------------
462 for ( int ev = firstevent ; ev < firstevent + nevents; ev++)
463 {
464 // Get an Event --> consists of 1440 Pixel ...erm....data
465 datafile->GetRow( ev - 1 );
466
467 if (verbosityLevel == 1)
468 {
469 cout << "-------------------------------------" << endl
470 << "...processing Set from Pixel "
471 << firstPixelOfSet
472 << " to Pixel "
473 << lastPixelOfSet
474 << "... Event: " << CurrentEventID
475 << "/" << nevents << endl;
476 }
477
478//--------------------------------------------------------------------
479// Loops over every Pixel of a Set of Pixels
480//--------------------------------------------------------------------
481 for ( int pixelID = firstPixelOfSet;
482 pixelID < lastPixelOfSet + 1
483 && pixelID < firstpixel + npixel;
484 pixelID++ )
485 {
486 if (verbosityLevel > 1)
487 {
488 cout << "-------------------------------------" << endl
489 << "...processing Set from Pixel "
490 << firstPixelOfSet
491 << " to Pixel "
492 << lastPixelOfSet
493 << "... Event: " << CurrentEventID
494 << "/" << nevents << endl
495 << " Pixel: " << pixelID
496 << "/" << firstpixel + npixel -1 << endl;
497 }
498
499//-------------------------------------
500// Create individual Pixel
501//-------------------------------------
502 if (ev == firstevent)
503 {
504 pixel[pixelID] = new Pixel(
505 pixelID,
506 maxPulseOrder,
507 verbosityLevel,
508 stats,
509 histoOptions,
510 gPixelOverlayXaxisLeft,
511 gPixelOverlayXaxisRight,
512 gBSLMean,
513 gGainMean
514 );
515 }
516
517//-------------------------------------
518// Apply Calibration
519//-------------------------------------
520 if (verbosityLevel > 2) cout << "...applying DrsCalibration";
521 applyDrsCalibration(
522 Ameas,
523 pixelID,
524 12,
525 12,
526 drs_basemean,
527 drs_gainmean,
528 drs_triggeroffsetmean,
529 gRegionOfInterest,
530 AllPixelDataVector,
531 StartCellVector
532 );
533 if (verbosityLevel > 2) cout << "...done " << endl;
534
535//-------------------------------------
536// Apply Filters
537//-------------------------------------
538 // finds spikes in the raw data, and interpolates the value
539 // spikes are: 1 or 2 slice wide, positive non physical artifacts
540 if (verbosityLevel > 2) cout << "...removeing Spikes";
541 removeSpikes (Ameas, Vcorr);
542 if (verbosityLevel > 2) cout << "...done " << endl;
543
544 // filter Vcorr with sliding average using FIR filter function
545 if (verbosityLevel > 2) cout << "...applying sliding average filter";
546 sliding_avg(Vcorr, Vslide, avg1);
547 if (verbosityLevel > 2) cout << "...done " << endl;
548
549 // filter Vslide with CFD using FIR filter function
550 if (verbosityLevel > 2) cout << "...apllying factfir filter";
551 factfir(b_cfd , a_cfd, k_cfd, Vslide, Vcfd);
552 if (verbosityLevel > 2) cout << "...done " << endl;
553
554 // filter Vcfd with sliding average using FIR filter function
555 if (verbosityLevel > 2) cout << "...applying 2nd sliding average filter";
556 sliding_avg(Vcfd, Vcfd2, avg2);
557 if (verbosityLevel > 2) cout << "...done " << endl;
558
559//-------------------------------------
560// Search vor Zero crossings
561//-------------------------------------
562 if (verbosityLevel > 2) cout << endl << "...searching zero crossings" ;
563 // peaks in Ameas[] are found by searching for zero crossings
564 // in Vcfd2
565 // first Argument 1 means ... *rising* edge
566 // second Argument 1 means ... search with stepsize 1 ... 10 is okay as well
567 vector<Region>* pZXings = zerosearch( Vcfd2 , 1 , 8);
568 // pZXings means "zero cross ings"
569 EnlargeRegion(*pZXings, 10, 10);
570 findAbsMaxInRegions(*pZXings, Vslide);
571 removeMaximaBelow( *pZXings, 3.0);
572 removeRegionWithMaxOnEdge( *pZXings, 2);
573 removeRegionOnFallingEdge( *pZXings, 100);
574 findTimeOfHalfMaxLeft(*pZXings, Vslide, gBSLMean, 5, 10, verbosityLevel );
575 removeRegionWithToFlatSlope(*pZXings, 0.5);
576 if (verbosityLevel > 2) cout << "...done" << endl;
577
578//-----------------------------------------------------------------------------
579// Fill Overlay Histos
580//-----------------------------------------------------------------------------
581 if (verbosityLevel > 2)
582 {
583 cout << endl << "Filling Histograms of Pixel# [Set] " << pixelID
584 << " Pixel# [Got] " << pixel[pixelID]->mChid;
585 }
586 FillHistograms(
587 pixel[pixelID],
588 pZXings,
589 AmplWindowWidth,
590 ev,
591// histoOptions,
592 maxPulseOrder,
593 verbosityLevel
594 );
595//-----------------------------------------------------------------------------
596// Spike Debug
597//-----------------------------------------------------------------------------
598 if ( spikeDebug )
599 {
600 BookDebugHistos(verbosityLevel );
601 breakout = ProduceDebugHistos( pZXings );
602 }// end of if(spikeDebug)
603
604 delete pZXings;
605 if (breakout) break;
606//-------------------------------------
607// Draw 1. Set of Pixel Histograms
608//-------------------------------------
609 if ((ev % refresh_rate) == 0)
610 {
611 if (ProduceGraphic)
612 {
613 if (debugPixel)
614 {
615// for ( int order = 0;
616// order < maxPulseOrder;
617// order++
618// )
619// {
620 pixel[firstPixelOfSet]->DrawOverlayHistograms(
621 cgpPixelPulses,
622 pulsesCanvasFrameNrs
623 );
624
625 pixel[firstPixelOfSet]->DrawDistributionHistograms(
626 cgpDistributions,
627 distributionCanvasFrameNrs
628 );
629// }
630 if (testmode)
631 {
632 DrawTestHistograms(
633 verbosityLevel
634 );
635 }
636
637 UpdateCanvases(
638 verbosityLevel,
639 maxPulseOrder,
640 testmode
641 );
642 }
643 }
644 }
645
646 if (breakout) break;
647
648 if (verbosityLevel > 1)
649 {
650 cout << endl << "...End of Pixel"
651 << endl << "-------------------------------------"<< endl;
652 }
653 }//End of Loop over Set
654//-------------------------------------
655// end of Loops over Set
656//-------------------------------------
657
658 }//End of Loop over Events
659//-------------------------------------
660// end of Loops over Events
661//-------------------------------------
662
663
664
665//-------------------------------------
666// Draw Pixel Histogramms of Overlay Spectra
667//-------------------------------------
668
669 for ( int pixelID = firstPixelOfSet;
670 pixelID < firstPixelOfSet + pixelSetSize
671 && pixelID < firstpixel + npixel;
672 pixelID++ )
673 {
674 //here is what happends at the end of each loop over all Events
675
676 if (ProduceGraphic)
677 {
678 pixel[pixelID]->DrawOverlayHistograms(
679 cgpPixelPulses,
680 pulsesCanvasFrameNrs
681 );
682
683 pixel[pixelID]->DrawDistributionHistograms(
684 cgpDistributions,
685 distributionCanvasFrameNrs
686 );
687
688 UpdateCanvases(
689 verbosityLevel,
690 maxPulseOrder,
691 testmode
692 );
693 }
694
695 //Save Histograms of Pixels into Output rootfile
696 if (verbosityLevel > 2)
697 {
698 cout << endl << "Saving Pixel# [Set] " << pixelID
699 << " Pixel# [Got] " << pixel[pixelID]->mChid;
700 }
701 pixel[pixelID]->SavePixelHistograms( OutRootFileName, saveResults );
702
703
704 if (debugPixel)
705 {
706 //Process gui events asynchronously during input
707 cout << endl;
708 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
709 timer.TurnOn();
710 TString input = Getline("Type 'q' to exit, <return> to go on: ");
711 timer.TurnOff();
712 if (input=="q\n") {
713 break;
714 }
715 }
716
717 if (verbosityLevel > 2)
718 {
719 cout << endl << "Deleting Pixel# [Set] " << pixelID
720 << " Pixel# [Got] " << pixel[pixelID]->mChid;
721 }
722 delete pixel[pixelID];
723
724 if (verbosityLevel > 2)
725 {
726 cout << endl << "...End of Pixel-Set"
727 << endl << "------------------------------------------------"
728 << endl;
729 }
730 }
731 }
732 // End of Loop over all Pixels
733
734//-------------------------------------
735// Draw Histograms
736//-------------------------------------
737 TString test = "root";
738// SaveHistograms( //save histograms of generell results into output root file
739// OutRootFileName,
740// test,
741// hRootList,
742// saveResults,
743// verbosityLevel
744// );
745
746 if (ProduceGraphic)
747 {
748 UpdateCanvases(
749 verbosityLevel,
750 maxPulseOrder,
751 testmode
752 );
753
754 //Process gui events asynchronously during input
755 cout << endl;
756 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
757 timer.TurnOn();
758 Getline("Press <return> to go on");
759 timer.TurnOff();
760
761 DeletePixelCanvases(
762 maxPulseOrder,
763 verbosityLevel
764 );
765
766 if (testmode)
767 delete cgpTestHistos;
768
769 if (spikeDebug)
770 delete cFiltered;
771 }
772
773 delete datafile;
774 delete[] pixel;
775
776 return( 0 );
777}
778//----------------------------------------------------------------------------
779// end of main function
780//-----------------------------------------------------------------------------
781
782
783
784
785//-----------------------------------------------------------------------------
786// Funktions
787//-----------------------------------------------------------------------------
788
789/*
790ich erzeuge sowas wie 36 pixel mit ca 20 Histogrammen pro pixel
791wie viel speicher braucht das?
792Ich brauche eine möglichkeit jedes dieser histogramme zu identifizieren
793am besten ein array oder eine liste oder einen abstracten datentyp mit den histogrammen als attribute.
794ne klasse wäre super
795
796*/
797
798void BookTestHistos( int verbosityLevel )
799{
800 if (verbosityLevel > 2) cout << endl << "...book pixel histograms" << endl;
801
802 hTesthisto = new TH1F (
803 "hTesthisto",
804 "Deviation of rising edge and maximum",
805 600,
806 -10.1,
807 10.1
808 );
809 hTesthisto->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
810 hTesthisto->GetYaxis()->SetTitle( "counts" );
811
812 hTesthisto2 = new TH2F (
813 "hTesthisto2",
814 "Deviation of rising edge and maximum by event #",
815 gNEvents,
816 250,
817 800,
818 600,
819 -10.1,
820 10.1
821 );
822// hTesthisto2->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
823// hTesthisto2->GetYaxis()->SetTitle( "counts" );
824// hTesthisto2->SetMarkerStyle( 2 );
825
826 if (verbosityLevel > 2) cout << "...done" << endl;
827}
828//end of BookTestHistos
829//----------------------------------------------------------------------------
830
831
832void
833BookDebugHistos( int verbosityLevel )
834{
835 if (verbosityLevel > 2) cout << endl << "...book histograms" << endl;
836 hRootList = new TList();
837 debugHistos = new TH1F[ Number_Of_Debug_Histo_Types ];
838 for ( int type = 0; type < Number_Of_Debug_Histo_Types; type++){
839 debugHistos[ type ].SetBins(1024, 0, 1024);
840 debugHistos[ type ].SetLineColor(1);
841 debugHistos[ type ].SetLineWidth(2);
842 debugHistos[ type ].SetStats(false);
843
844 // set X axis paras
845 debugHistos[ type ].GetXaxis()->SetLabelSize(0.08);
846 debugHistos[ type ].GetXaxis()->SetTitleSize(0.08);
847 debugHistos[ type ].GetXaxis()->SetTitleOffset(1.0);
848 debugHistos[ type ].GetXaxis()->SetTitle(Form("Time slice [%.1f ns/slice]", 1./2.));
849
850 // set Y axis paras
851 debugHistos[ type ].GetYaxis()->SetLabelSize(0.08);
852 debugHistos[ type ].GetYaxis()->SetTitleSize(0.08);
853 debugHistos[ type ].GetYaxis()->SetTitleOffset(0.3);
854 debugHistos[ type ].GetYaxis()->SetTitle("Amplitude [mV]");
855// hRootList->Add( &debugHistos[ type ] );
856 }
857 if (verbosityLevel > 2) cout << "...done" << endl;
858}
859// end of BookDebugHistos
860//----------------------------------------------------------------------------
861
862
863void
864FillHistograms(
865 Pixel* CurrentPixel,
866 vector<Region>* pZXings,
867 int AmplWindowWidth,
868 int eventNumber,
869 int maxPulseOrder,
870 int verbosityLevel
871 )
872{
873 if (verbosityLevel > 2) cout << endl << "...filling pulse histograms" ;
874 if (verbosityLevel > 3) cout << endl << "...EventNR " << eventNumber ;
875 bool use_this_peak=false;
876 int order_of_pulse=0;
877 vector<Region>::iterator reg;
878 //Loop over all found zerocrossings in Region
879 for (reg = pZXings->begin() ; reg < pZXings->end() ; reg++)
880 {
881 //skip those who are at the Rim of the ROI
882 if (reg->maxPos < 12 || (unsigned int) reg->maxPos > gRegionOfInterest-12)
883 {
884 if (verbosityLevel > 2) cout << endl << "\t...out of range" << endl;
885 continue;
886 }
887
888 // Define axis range of Histogram
889 int Left = reg->maxPos - gPixelOverlayXaxisLeft;
890 int Right = reg->maxPos + gPixelOverlayXaxisRight;
891 int EdgeLeft = reg->halfRisingEdgePos - gPixelOverlayXaxisLeft;
892 int EdgeRight = reg->halfRisingEdgePos + gPixelOverlayXaxisRight;
893
894 if (Left < 0) Left = 0;
895 if (EdgeLeft < 0) EdgeLeft = 0;
896 if (Right > (int)Ameas.size() ) Right = Ameas.size();
897 if (EdgeRight > (int)Ameas.size() ) EdgeRight = Ameas.size();
898
899
900 if (verbosityLevel > 2) cout << endl << "\t...choosing Histogram" ;
901 //determine order of pulse and which histogram shall be filled
902 for ( int order = 0; order < maxPulseOrder; order++)
903 {
904 use_this_peak = UseThisPulse(
905 reg->maxPos,
906 order,
907 AmplWindowWidth,
908 gGainMean,
909 maxPulseOrder,
910 verbosityLevel
911 );
912 if ( use_this_peak)
913 {
914 order_of_pulse = order;
915 break;
916 }
917 }
918 //Fill histograms
919 if (use_this_peak)
920 {
921 //Histograms for Distributions
922 if (CurrentPixel->mOptions.Contains("S") )
923 {
924 if (verbosityLevel > 2)
925 cout << endl << "\t...filling Edge Histogram" ;
926 CurrentPixel->hSlopeRisingEdge[order_of_pulse]->Fill( reg->slopeOfRisingEdge ) ;
927 }
928
929 if (CurrentPixel->mOptions.Contains("R") )
930 {
931 if (verbosityLevel > 2)
932 cout << endl << "\t...filling RisingEdgeToMax Histogram" ;
933 CurrentPixel->hRisingEdgeToMax[order_of_pulse]->Fill( reg->distanceEdgeToMax ) ;
934 }
935
936 if (CurrentPixel->mOptions.Contains("M") )
937 {
938 if (verbosityLevel > 2)
939 cout << endl << "\t...filling PosOfMax Histogram" ;
940 CurrentPixel->hPosOfMax[order_of_pulse]->Fill( reg->maxPos ) ;
941 }
942
943 //Histograms for Maximum Overlay
944 for ( int pos = Left; pos < Right; pos++)
945 {
946 CurrentPixel->hMaxOverlay[order_of_pulse]->Fill( pos - (reg->maxPos), Ameas[pos]) ;
947 CurrentPixel->hMaxProfile[order_of_pulse]->Fill( pos - (reg->maxPos), Ameas[pos]) ;
948 }
949 //Histograms for Edge Overlay
950 for ( int pos = EdgeLeft; pos < EdgeRight; pos++)
951 {
952 CurrentPixel->hEdgeOverlay[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), Ameas[pos]) ;
953 CurrentPixel->hEdgeProfile[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), Ameas[pos]) ;
954 }
955
956 if (verbosityLevel > 2) cout << "...done" << endl;
957 }
958
959 //Breakout option
960 if (breakout)
961 {
962 break;
963 }
964 }
965 if (verbosityLevel > 2) cout << "...done" << endl;
966}
967// end of FillHistograms
968//----------------------------------------------------------------------------
969
970bool UseThisPulse(
971 int maxPos,
972 int order,
973 int AmplWindowWidth,
974 float gGainMean,
975 int maxPulseOrder,
976 int verbosityLevel
977 )
978{
979 //determine if pulse is of given order
980bool use_this_peak = false;
981 if (Ameas[maxPos] >= ((gGainMean*(order+1)) - (AmplWindowWidth/2))
982 && Ameas[maxPos] < ((gGainMean*(order+1)) + AmplWindowWidth/2))
983 {
984 if (verbosityLevel > 3) cout << "...#" << order ;
985 use_this_peak = true;
986 }
987 else
988 {
989 if (verbosityLevel > 3)
990 {
991 if (order >= (maxPulseOrder - 1))
992 {
993 cout << "...NONE" << endl ;
994 }
995 }
996 use_this_peak = false;
997 }
998
999return use_this_peak;
1000}
1001
1002
1003void
1004DrawTestHistograms(
1005 int verbosityLevel
1006 )
1007{
1008 if (verbosityLevel > 2) cout << endl << "...drawing pulse histograms TODO CANVAS" ;
1009
1010// cgpTestHistos->cd(1);
1011// hTesthisto->Draw();
1012// cgpTestHistos->cd(2);
1013// hTesthisto2->Draw("BOX");
1014}
1015// end of DrawTestHistograms
1016//----------------------------------------------------------------------------
1017
1018void
1019UpdateCanvases(
1020 int verbosityLevel,
1021 int max_pulse_order,
1022 bool testmode
1023 )
1024{
1025 if (verbosityLevel > 3) cout << endl << "...updating canvases" ;
1026 for (int pulse_order = 0; pulse_order < max_pulse_order; pulse_order++)
1027 {
1028 cgpPixelPulses[pulse_order]->Modified();
1029 cgpPixelPulses[pulse_order]->Update();
1030 cgpDistributions[pulse_order]->Modified();
1031 cgpDistributions[pulse_order]->Update();
1032 if ( testmode )
1033 {
1034 cgpTestHistos->Modified();
1035 cgpTestHistos->Update();
1036 }
1037 }
1038}
1039// end of UpdateCanvases
1040//----------------------------------------------------------------------------
1041
1042void
1043DeletePixelCanvases(
1044 int maxPulseOrder,
1045 int verbosityLevel
1046 )
1047{
1048 if (verbosityLevel > 2) cout << endl << "...delete pixel Canvases" ;
1049 for (int pulse_order = 0; pulse_order < maxPulseOrder; pulse_order++ )
1050 {
1051 delete cgpPixelPulses[pulse_order];
1052 delete cgpDistributions[pulse_order];
1053 }
1054 delete[] cgpPixelPulses;
1055 delete[] cgpDistributions;
1056}
1057
1058//----------------------------------------------------------------------------
1059//----------------------------------------------------------------------------
1060
1061bool
1062ProduceDebugHistos(
1063 vector<Region>* pZXings
1064 )
1065{
1066 // TODO do this correct. The vectors should be the rigt ones... this is just luck
1067 debugHistos[Ameas_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
1068 debugHistos[Vcorr_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
1069 debugHistos[Vslide_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
1070 debugHistos[Vcfd_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
1071 debugHistos[Vcfd2_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
1072
1073 for ( unsigned int sl = 0; sl < gRegionOfInterest; sl++)
1074 {
1075 debugHistos[Ameas_].SetBinContent(sl, Ameas[sl]);
1076 debugHistos[Vcorr_].SetBinContent(sl, Vcorr[sl]);
1077 debugHistos[Vslide_].SetBinContent( sl, Vslide[sl] );
1078 debugHistos[Vcfd_].SetBinContent( sl, Vcfd[sl] );
1079 debugHistos[Vcfd2_].SetBinContent( sl, Vcfd2[sl] );
1080 }
1081
1082 bool has_negative_slope = false;
1083 cFiltered->cd(1);
1084 gPad->SetGrid();
1085 debugHistos[Ameas_].Draw();
1086
1087 TF1 *OneEdge;
1088 TLine *OneLine;
1089 TBox *OneBox;
1090 vector<TF1*> MyEdges;
1091 vector<TBox*> MyBoxes;
1092 vector<TBox*> MyLines;
1093
1094 int left = debugHistos[Ameas_].GetXaxis()->GetFirst();
1095 int right = debugHistos[Ameas_].GetXaxis()->GetLast();
1096 vector<Region>::iterator reg;
1097 for (reg = pZXings->begin() ; reg < pZXings->end() ; reg++)
1098 {
1099 if (reg->slopeOfRisingEdge < 0)
1100 {
1101 has_negative_slope = true;
1102 }
1103 }
1104
1105 if (has_negative_slope)
1106 {
1107 for (unsigned int i=0; i<pZXings->size(); i++){
1108 Region* reg = &pZXings->at(i);
1109 cout << "Slope is: " << reg->slopeOfRisingEdge << endl;
1110
1111// xBegin = (reg->halfRisingEdgePos - reg->distanceEdgeToMax);
1112// yBegin = reg->interceptRisingEdge + reg->slopeOfRisingEdge * (xBegin)
1113
1114 OneEdge = new TF1(
1115 "OneEdge",
1116 "[0]*x+[1]",
1117 left,
1118 right
1119 );
1120
1121 OneEdge->SetParameters(
1122 reg->slopeOfRisingEdge,
1123 reg->interceptRisingEdge
1124 );
1125
1126 OneEdge->SetLineColor(kRed);
1127 OneEdge->SetLineWidth(1);
1128 MyEdges.push_back(OneEdge);
1129
1130 OneEdge->Draw("SAME");
1131// delete OneEdge;
1132
1133 OneBox = new TBox(
1134 pZXings->at(i).maxPos -10 ,
1135 pZXings->at(i).maxVal -0.5,
1136 pZXings->at(i).maxPos +10 ,
1137 pZXings->at(i).maxVal +0.5);
1138 OneBox->SetLineColor(kBlue);
1139 OneBox->SetLineWidth(1);
1140 OneBox->SetFillStyle(0);
1141 OneBox->SetFillColor(kRed);
1142 MyBoxes.push_back(OneBox);
1143 OneBox->Draw();
1144
1145 OneLine = new TLine(
1146 pZXings->at(i).halfRisingEdgePos,
1147 Ameas[pZXings->at(i).halfRisingEdgePos],
1148 pZXings->at(i).halfRisingEdgePos,
1149 pZXings->at(i).halfRisingEdgePos * pZXings->at(i).slopeOfRisingEdge + pZXings->at(i).interceptRisingEdge
1150 );
1151 OneLine->SetLineColor(kBlue);
1152 OneLine->SetLineWidth(3);
1153 MyLines.push_back(OneBox);
1154 OneLine->Draw();
1155 }
1156
1157 cFiltered->cd(2);
1158 gPad->SetGrid();
1159 debugHistos[Vslide_].Draw();
1160 debugHistos[Vcorr_].Draw("SAME");
1161
1162 for (unsigned int i=0; i<MyEdges.size(); i++){
1163 MyEdges.at(i)->Draw("SAME");
1164 }
1165 for (unsigned int i=0; i<MyBoxes.size(); i++){
1166 MyBoxes.at(i)->Draw();
1167 }
1168 for (unsigned int i=0; i<MyLines.size(); i++){
1169 MyLines.at(i)->Draw();
1170 }
1171
1172 cFiltered->cd(3);
1173 gPad->SetGrid();
1174 debugHistos[Vcorr_].SetLineColor(kBlue);
1175 debugHistos[Vcorr_].Draw();
1176 debugHistos[Ameas_].Draw("SAME");
1177 debugHistos[Vslide_].SetLineColor(kRed);
1178 debugHistos[Vslide_].Draw("SAME");
1179 debugHistos[Vcfd_].SetLineColor(kGreen);
1180 debugHistos[Vcfd_].Draw("SAME");
1181 debugHistos[Vcfd2_].SetLineColor(kViolet);
1182 debugHistos[Vcfd2_].Draw("SAME");
1183 for (unsigned int i=0; i<MyEdges.size(); i++){
1184 MyEdges.at(i)->Draw("SAME");
1185 }
1186 for (unsigned int i=0; i<MyBoxes.size(); i++){
1187 MyBoxes.at(i)->Draw();
1188 }
1189 for (unsigned int i=0; i<MyLines.size(); i++){
1190 MyLines.at(i)->Draw();
1191 }
1192
1193
1194// cFiltered->cd(3);
1195// gPad->SetGrid();
1196// debugHistos[Vcfd2_].Draw();
1197// TLine *zeroline = new TLine(0, 0, 1024, 0);
1198// zeroline->SetLineColor(kBlue);
1199// zeroline->Draw();
1200
1201 cFiltered->Update();
1202
1203
1204 //Process gui events asynchronously during input
1205 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
1206 timer.TurnOn();
1207 TString input = Getline("Type 'q' to exit, <return> to go on: ");
1208 timer.TurnOff();
1209 if (input=="q\n") {
1210 return true;
1211 }
1212
1213 MyLines.clear();
1214 MyBoxes.clear();
1215 MyEdges.clear();
1216
1217// for (unsigned int i=0; i<pZXings->size(); i++)
1218// {
1219// delete MyEdges.at(i)
1220// delete MyBoxes.at(i);
1221// }
1222// delete OneBox;
1223// delete OneEdge;
1224// delete OneLine;
1225
1226// MyEdges.erase();
1227// MyBoxes.erase();
1228 }
1229
1230 //TODO!!!!!!!!!
1231 // do some Garbage collection ...
1232 // all the Objects on the heap should be deleted here.
1233return false;
1234}// end of if(spikeDebug)
1235
1236
1237
1238
1239
1240
1241int main(int argc,char *argv[])
1242{
1243
1244 TString test;
1245 TString rcFileName;
1246 TString processType = "overlay";
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 // decode arguments
1254 if(argc < 2)
1255 {
1256 printf("no arguements given, using standard arguments\n");
1257 }
1258
1259 // set conditions for functions arguments
1260 for (int i=1;i<argc;i++)
1261 {
1262 test = argv[i];
1263
1264 if (test.Contains("--config_file") || test.Contains("-c"))
1265 {
1266 cout << "RC-File: \"" << argv[i + 1] << "\"" << endl;
1267 rcFileName = argv[i + 1];
1268 rcFileNameCmdSet = true;
1269 continue;
1270 }
1271
1272 if (test.Contains("--verbosity") || test.Contains("-v"))
1273 {
1274 cout << "Verbosity Level: \"" << argv[i + 1] << "\"" << endl;
1275 verbLevel = atoi(argv[i + 1]);
1276 continue;
1277 }
1278
1279 if (test.Contains("--save") || test.Contains("-s"))
1280 {
1281 cout << "will save results" << endl;
1282 save = true;
1283 continue;
1284 }
1285
1286 if (test.Contains("--graphics") || test.Contains("-g"))
1287 {
1288 cout << "will produce graphics" << endl;
1289 produceGraphic = true;
1290 continue;
1291 }
1292 }
1293
1294 // reading rc-File:
1295 configfile rcfile( rcFileName, processType );
1296
1297 if (save)
1298 {
1299 rcfile.mSave = save;
1300 }
1301
1302 if (verbLevelCmdSet)
1303 {
1304 rcfile.mVerbLevel = verbLevel;
1305 }
1306
1307 if (produceGraphic)
1308 {
1309 rcfile.mProduceGraphic = produceGraphic;
1310 }
1311
1312 FPulseOverlay(
1313 rcfile.mInputPath,
1314 rcfile.mDataFileName,
1315 rcfile.mDrsFileName,
1316 rcfile.mOutputPath,
1317 rcfile.mOutputFile,
1318 rcfile.mFirstEvent,
1319 rcfile.mNumEvents,
1320 rcfile.mFirstPixel,
1321 rcfile.mNumPixel,
1322 rcfile.mPixelSetSize,
1323 rcfile.mMaxOrder,
1324 rcfile.mVerbLevel,
1325 rcfile.mHistoOptions,
1326 rcfile.mAmplWindowWidth,
1327 rcfile.mGainMean,
1328 rcfile.mBSLMean,
1329 rcfile.mAvg1,
1330 rcfile.mAvg2,
1331 rcfile.mOverlayWindowLeft,
1332 rcfile.mOverlayWindowRight,
1333 rcfile.mProduceGraphic,
1334 rcfile.mSpikeDebug,
1335 rcfile.mDbgPixel,
1336 rcfile.mPrintStats,
1337 rcfile.mSave,
1338 rcfile.mTestMode,
1339 rcfile.mRefreshRate
1340 );
1341 return 0;
1342
1343}
Note: See TracBrowser for help on using the repository browser.