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

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