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

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