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

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