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

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