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

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