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

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