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

Last change on this file since 13602 was 13602, checked in by Jens Buss, 13 years ago
line that could be uncommented to also save debug histos
File size: 38.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);
180void 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 const char* datafilename = "/fact/raw/2011/11/09/20111109_006.fits.gz",
192 const char* drsfilename = "/fact/raw/2011/11/09/20111109_003.drs.fits.gz",
193 const char* OutRootFileName = "/home_nfs/isdc/jbbuss/analysis/FPulseTemplate/test.root",
194 bool ProduceGraphic = true,
195 bool spikeDebug = true,
196 bool debugPixel = true,
197 bool testmode = false,
198 bool saveResults = true,
199 int verbosityLevel = 1, // different verbosity levels can be implemented here
200 int firstevent = 1, // Hast to be between 1 and infinity NOT 0!!!!
201 int nevents = 100,
202 int firstpixel = 15,
203 int npixel = 3,
204 int pixelSetSize = 40,
205 int maxPulseOrder = 3,
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 if (spikeDebug)
389 BookDebugHistos(verbosityLevel );
390
391 if (testmode)
392 BookTestHistos( verbosityLevel );
393
394
395//-----------------------------------------------------------------------------
396// Main Cycle
397//-----------------------------------------------------------------------------
398
399//----------------------------------------------------------------------------
400// Initialize Pixel
401//----------------------------------------------------------------------------
402 Pixel** pixel = new Pixel*[NPIX];
403
404 for (int i = 0 ; i < NPIX; i++)
405 {
406 pixel[i] = NULL;
407 }
408//-------------------------------------
409// Loop over Pixel Sets
410//-------------------------------------
411 for ( int firstPixelOfSet = firstpixel;
412 firstPixelOfSet < firstpixel + npixel;
413 firstPixelOfSet += pixelSetSize )
414 {
415
416 if (verbosityLevel == 0)
417 {
418 cout << "------------------------------------------------" << endl
419 << "...processing Set from Pixel "
420 << firstPixelOfSet
421 << " to Pixel "
422 << firstPixelOfSet+pixelSetSize-1 << endl;
423 }
424
425//--------------------------------------------------------------------
426// Loop over every Event of Pixel Set
427//--------------------------------------------------------------------
428 for ( int ev = firstevent ; ev < firstevent + nevents; ev++)
429 {
430 // Get an Event --> consists of 1440 Pixel ...erm....data
431 datafile->GetRow( ev - 1 );
432
433 if (verbosityLevel == 1)
434 {
435 cout << "-------------------------------------" << endl
436 << "...processing Set from Pixel "
437 << firstPixelOfSet
438 << " to Pixel "
439 << firstPixelOfSet+pixelSetSize-1
440 << "... Event: " << CurrentEventID
441 << "/" << nevents << endl;
442 }
443
444//--------------------------------------------------------------------
445// Loops over every Pixel of a Set of Pixels
446//--------------------------------------------------------------------
447 for ( int pixelID = firstPixelOfSet;
448 pixelID < firstPixelOfSet + pixelSetSize
449 && pixelID < firstpixel + npixel;
450 pixelID++ )
451 {
452 if (verbosityLevel > 1)
453 {
454 cout << "-------------------------------------" << endl
455 << "...processing Set from Pixel "
456 << firstPixelOfSet
457 << " to Pixel "
458 << firstPixelOfSet+pixelSetSize-1
459 << "... Event: " << CurrentEventID
460 << "/" << nevents << endl
461 << " Pixel: " << pixelID
462 << "/" << firstpixel + npixel -1 << endl;
463 }
464
465//-------------------------------------
466// Create individual Pixel
467//-------------------------------------
468 if (ev == firstevent)
469 {
470 pixel[pixelID] = new Pixel(
471 pixelID,
472 maxPulseOrder,
473 verbosityLevel,
474 false,
475 NULL,
476 gPixelOverlayXaxisLeft,
477 gPixelOverlayXaxisRight,
478 gBSLMean,
479 gGainMean,
480 histoOptions
481 );
482 }
483
484//-------------------------------------
485// Apply Calibration
486//-------------------------------------
487 if (verbosityLevel > 2) cout << "...applying DrsCalibration";
488 applyDrsCalibration(
489 Ameas,
490 pixelID,
491 12,
492 12,
493 drs_basemean,
494 drs_gainmean,
495 drs_triggeroffsetmean,
496 gRegionOfInterest,
497 AllPixelDataVector,
498 StartCellVector
499 );
500 if (verbosityLevel > 2) cout << "...done " << endl;
501
502//-------------------------------------
503// Apply Filters
504//-------------------------------------
505 // finds spikes in the raw data, and interpolates the value
506 // spikes are: 1 or 2 slice wide, positive non physical artifacts
507 if (verbosityLevel > 2) cout << "...removeing Spikes";
508 removeSpikes (Ameas, Vcorr);
509 if (verbosityLevel > 2) cout << "...done " << endl;
510
511 // filter Vcorr with sliding average using FIR filter function
512 if (verbosityLevel > 2) cout << "...applying sliding average filter";
513 sliding_avg(Vcorr, Vslide, avg1);
514 if (verbosityLevel > 2) cout << "...done " << endl;
515
516 // filter Vslide with CFD using FIR filter function
517 if (verbosityLevel > 2) cout << "...apllying factfir filter";
518 factfir(b_cfd , a_cfd, k_cfd, Vslide, Vcfd);
519 if (verbosityLevel > 2) cout << "...done " << endl;
520
521 // filter Vcfd with sliding average using FIR filter function
522 if (verbosityLevel > 2) cout << "...applying 2nd sliding average filter";
523 sliding_avg(Vcfd, Vcfd2, avg2);
524 if (verbosityLevel > 2) cout << "...done " << endl;
525
526//-------------------------------------
527// Search vor Zero crossings
528//-------------------------------------
529 if (verbosityLevel > 2) cout << endl << "...searching zero crossings" ;
530 // peaks in Ameas[] are found by searching for zero crossings
531 // in Vcfd2
532 // first Argument 1 means ... *rising* edge
533 // second Argument 1 means ... search with stepsize 1 ... 10 is okay as well
534 vector<Region>* pZXings = zerosearch( Vcfd2 , 1 , 8);
535 // pZXings means "zero cross ings"
536 EnlargeRegion(*pZXings, 10, 10);
537 findAbsMaxInRegions(*pZXings, Vslide);
538 removeMaximaBelow( *pZXings, 3.0);
539 removeRegionWithMaxOnEdge( *pZXings, 2);
540 removeRegionOnFallingEdge( *pZXings, 100);
541 findTimeOfHalfMaxLeft(*pZXings, Vcorr, gBSLMean, 5, 10, verbosityLevel );
542 removeRegionWithToFlatSlope(*pZXings, 0.5);
543 if (verbosityLevel > 2) cout << "...done" << endl;
544
545//-----------------------------------------------------------------------------
546// Fill Overlay Histos
547//-----------------------------------------------------------------------------
548 if (verbosityLevel > 2)
549 {
550 cout << endl << "Filling Histograms of Pixel# [Set] " << pixelID
551 << " Pixel# [Got] " << pixel[pixelID]->mChid;
552 }
553 FillHistograms(
554 pixel[pixelID],
555 pZXings,
556 AmplWindowWidth,
557 ev,
558// histoOptions,
559 maxPulseOrder,
560 verbosityLevel
561 );
562//-----------------------------------------------------------------------------
563// Spike Debug
564//-----------------------------------------------------------------------------
565 if ( spikeDebug )
566 {
567 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
1027void
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 bool has_negative_slope = false;
1048 cFiltered->cd(1);
1049 gPad->SetGrid();
1050 debugHistos[Ameas_].Draw();
1051
1052 TF1 *OneEdge;
1053// vector<TF1*> MyEdges;
1054 int left = debugHistos[Ameas_].GetXaxis()->GetFirst();
1055 int right = debugHistos[Ameas_].GetXaxis()->GetLast();
1056 vector<Region>::iterator reg;
1057 for (reg = pZXings->begin() ; reg < pZXings->end() ; reg++)
1058 {
1059 if (reg->slopeOfRisingEdge < 0)
1060 {
1061 has_negative_slope = true;
1062 cout << "Slope is: " << reg->slopeOfRisingEdge << endl;
1063
1064// xBegin = (reg->halfRisingEdgePos - reg->distanceEdgeToMax);
1065// yBegin = reg->interceptRisingEdge + reg->slopeOfRisingEdge * (xBegin)
1066
1067 OneEdge = new TF1(
1068 "OneEdge",
1069 "[0]*x+[1]",
1070 left,
1071 right
1072 );
1073
1074 OneEdge->SetParameters(
1075 reg->slopeOfRisingEdge,
1076 reg->interceptRisingEdge
1077 );
1078
1079 OneEdge->SetLineColor(kRed);
1080 OneEdge->SetLineWidth(1);
1081// MyEdges.push_back(OneEdge);
1082
1083 OneEdge->Draw("SAME");
1084 }
1085 }
1086
1087 TLine *OneLine;
1088 TBox *OneBox;
1089// vector<TBox*> MyBoxes;
1090 for (unsigned int i=0; i<pZXings->size(); i++){
1091 OneBox = new TBox(
1092 pZXings->at(i).maxPos -10 ,
1093 pZXings->at(i).maxVal -0.5,
1094 pZXings->at(i).maxPos +10 ,
1095 pZXings->at(i).maxVal +0.5);
1096 OneBox->SetLineColor(kBlue);
1097 OneBox->SetLineWidth(1);
1098 OneBox->SetFillStyle(0);
1099 OneBox->SetFillColor(kRed);
1100// MyBoxes.push_back(OneBox);
1101 OneBox->Draw();
1102
1103 OneLine = new TLine(
1104 pZXings->at(i).halfRisingEdgePos,
1105 Ameas[pZXings->at(i).halfRisingEdgePos],
1106 pZXings->at(i).halfRisingEdgePos,
1107 pZXings->at(i).halfRisingEdgePos * pZXings->at(i).slopeOfRisingEdge + pZXings->at(i).interceptRisingEdge
1108 );
1109 OneLine->SetLineColor(kBlue);
1110 OneLine->SetLineWidth(3);
1111// MyBoxes.push_back(OneBox);
1112 OneLine->Draw();
1113 }
1114
1115 cFiltered->cd(2);
1116 gPad->SetGrid();
1117 debugHistos[Vcorr_].Draw();
1118
1119 cFiltered->cd(3);
1120 gPad->SetGrid();
1121 debugHistos[Vslide_].Draw();
1122
1123
1124// cFiltered->cd(3);
1125// gPad->SetGrid();
1126// debugHistos[Vcfd2_].Draw();
1127// TLine *zeroline = new TLine(0, 0, 1024, 0);
1128// zeroline->SetLineColor(kBlue);
1129// zeroline->Draw();
1130 if (has_negative_slope)
1131 {
1132 cFiltered->Update();
1133
1134
1135 //Process gui events asynchronously during input
1136 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
1137 timer.TurnOn();
1138 TString input = Getline("Type 'q' to exit, <return> to go on: ");
1139 timer.TurnOff();
1140 if (input=="q\n") {
1141 breakout=true;
1142 }
1143
1144
1145// for (unsigned int i=0; i<pZXings->size(); i++)
1146// {
1147// delete MyEdges.at(i)
1148// delete MyBoxes.at(i);
1149// }
1150 delete OneBox;
1151 delete OneEdge;
1152 delete OneLine;
1153
1154// MyEdges.erase();
1155// MyBoxes.erase();
1156 }
1157
1158 //TODO!!!!!!!!!
1159 // do some Garbage collection ...
1160 // all the Objects on the heap should be deleted here.
1161
1162}// end of if(spikeDebug)
1163
1164
1165int main()
1166{
1167
1168FPulseOverlay();
1169return 0;
1170
1171}
Note: See TracBrowser for help on using the repository browser.