source: fact/tools/rootmacros/PulseTemplates/FPulseOverlay_filelist.C@ 16318

Last change on this file since 16318 was 15372, checked in by Jens Buss, 13 years ago
add some new histos
File size: 61.6 KB
Line 
1/********************** FPulseOverlay ***********************
2* read FACT raw data
3* remove spikes
4* calculate baseline
5* find peaks (CFD and normal discriminator)
6+ seagRCh 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 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.h"
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#include "configfile.h"
73
74
75//----------------------------------------------------------------------------
76// Decleration of global variables
77//----------------------------------------------------------------------------
78
79bool gBreakout=false;
80
81vector<int16_t> gAllPixelDataVector;
82vector<int16_t> gStartCellVector;
83unsigned int gCurrentEventID;
84size_t gPXLxROI;
85UInt_t gRegionOfInterest;
86UInt_t gNumberOfPixels;
87TString ghistotitle;
88
89size_t gTriggerOffsetROI;
90size_t gRC;
91size_t gDrs_n;
92vector<float> gDrsBaseMean;
93vector<float> gDrsGainMean;
94vector<float> gDrsTriggeroffsetMean;
95
96vector<float> gAmeas(FAD_MAX_SAMPLES); // copy of the data (measured amplitude
97vector<float> gVcorr(FAD_MAX_SAMPLES); // corrected Values
98vector<float> gVslide(FAD_MAX_SAMPLES); // sliding average result
99vector<float> gVcfd(FAD_MAX_SAMPLES); // CDF result
100vector<float> gVcfd2(FAD_MAX_SAMPLES); // CDF result + 2nd sliding average
101
102int gNEvents=0;
103float gGainMean = 9;
104float gBSLMean = 0;
105
106//typedef struct{
107// double maxAmpl;
108// int countOfMax;
109// } OverlayMaximum;
110
111//typedef struct{
112// TString name;
113// TString title;
114// TString xTitle;
115// TString yTitle;
116// float yMax;
117// float yMin;
118//} histAttrib_t;
119
120//histAttrib_t* gHistoAttrib[maxPulseOrder];
121//histAttrib_t* gProfileAttrib[maxPulseOrder];
122
123
124// histograms
125const int Number_Of_Debug_Histo_Types = 7;
126
127const unsigned int
128 gAmeas_ = 0,
129 N1mean_ = 1,
130 gVcorr_ = 2,
131 Vtest_ = 3,
132 gVslide_ = 4,
133 gVcfd_ = 5,
134 gVcfd2_ = 6;
135
136//const char* gHistoTypes[8] = {
137// "hMaxOverlay",
138// "hPixelMax",
139// "hEdgeOverlay",
140// "hMaxProfile",
141// "hAllPixelOverlay",
142// "hAllPixelMax",
143// "hAllPixelMaxGaus",
144// "hAllPixelProfile"
145//};
146
147//----------------------------------------------------------------------------
148// Initialisation of histograms
149//----------------------------------------------------------------------------
150
151 // Temporary Objects
152 TH1F* debugHistos = NULL;
153 //TH2F* hOverlayTemp = NULL;
154 //TH1D* hProjPeak = NULL;
155 TH1F* hTesthisto = NULL;
156 TH2F* hTesthisto2 = NULL;
157
158 //Histogram Parameters
159 Int_t gPixelOverlayXaxisLeft = 0;
160 Int_t gPixelOverlayXaxisRight = 0;
161
162 //Root-File Objects
163// TObjArray* hList[pixelSetSize] = {NULL};
164 TList* hRootList = NULL;
165
166 TCanvas** cgpPixelPulses = NULL;
167 TCanvas** cgpDistributions = NULL;
168 TCanvas* cFiltered = NULL;
169 TCanvas* cgpTestHistos = NULL;
170
171//----------------------------------------------------------------------------
172// Functions
173//----------------------------------------------------------------------------
174
175void BookDebugHistos(int );
176void BookTestHistos( int );
177
178//void DeletePixelHistos(int );
179//void SaveHistograms( const char*, const char*, TObjArray*, int );
180
181void FillHistograms(Pixel*, vector<Region>*, int, int, int , bool , int verbosityLevel );
182void DrawTestHistograms( int);
183bool ProduceDebugHistos( vector<Region> *pZXings);
184bool UseThisPulse( int, int, int, float, float bslMean, int maxPulseOrder, int verbosityLevel);
185void UpdateCanvases( int, int, bool);
186void DeletePixelCanvases( int, int );
187
188int HandleFitsFile(
189 fits * datafile,
190 TString datafilename,
191 TString drsfilename,
192 int nevents,
193 int npixel,
194 int verbosityLevel
195 );
196
197void ReadSequenzFile(
198 TString fileName,
199 vector<TString> *fitsFiles,
200 vector<TString> *drsFiles
201 );
202
203TString gDataDirectory;
204TString gHomeDirectory;
205TString ghostname;
206
207//----------------------------------------------------------------------------
208//----------------------------------------------------------------------------
209// MAIN - Funtion
210//----------------------------------------------------------------------------
211//----------------------------------------------------------------------------
212int FPulseOverlay_filelist(
213 //---------------{Process Files}---------------
214 TString inputPath = "raw/2012/08/02/",
215 TString sequenzFileName = "/home_nfs/isdc/jbbuss/analysis/pulseShape/lists/20120802templ_vs_100mV.csv",
216 TString OutputPath = "analysis/PulsTemplateAnalysis/20120802/",
217 //TString OutRootFileName = "20120802_100mV_Stacking_G12B-2AWW8.root",
218 TString OutRootFileName = "new_dist_test.root",
219 //---------------{Process parameters}---------------
220 int nRuns = 5,
221 int firstevent = 1, // Hast to be between 1 and infinity NOT 0!!!!
222 int nevents = -1,
223 int firstpixel = 0,
224 int npixel = 1,
225 int pixelSetSize = 150,
226 int maxPulseOrder = 6,
227 //---------------{stacking parameters}---------------
228 TString histoOptions = "SRMTILDAP",
229 int AmplWindowWidth = 8, //Width of Window for selection of pulses to histograms
230 float GainMean = 12,
231 float BSLMean = -2, //4 Histogramms will be drawn, decide how far pulsehights differ from eachother
232 int avg1 = 8,
233 int avg2 = 8,
234 int OverlayWindowLeft = 100,
235 int OverlayWindowRight = 1024,
236
237 // pulses within this distance behind pulse max will be discarted
238 int beginRisingEdge = 5,
239 int endRisingEdge = 10,
240
241 // pulses within this distance behind pulse max will be discarted
242 float risingEdgeMinSlope = 0.5,
243 //---------------{graphical options}---------------
244 bool select_pulses = false,
245 bool ProduceGraphic = false,
246 bool spikeDebug = false,
247 bool debugPixel = false,
248 bool stats = true,
249 bool saveResults = true,
250 bool testmode = false,
251 //---------------{verbosity}---------------
252 int refresh_rate = 500, //refresh rate for canvases
253 int verbosityLevel = 0 // different verbosity levels can be implemented here
254 )
255{
256//-----------------------------------------------------------------------------
257// configuration
258//-----------------------------------------------------------------------------
259
260 // Amplitude below which pulses will be discarted
261 float minPulseAmplitude = 0.3*(GainMean)+BSLMean;
262 // Amplitude above which pulses will be discarted
263 float maxPulseAmplitude = (maxPulseOrder+1)*(GainMean)+BSLMean;
264 // pulses within this distance to timeline edge be discarted
265 int timeLineEdgeWidth = 2;
266 // pulses within this distance behind pulse max will be discarted
267 int FallingEdgeWidth = 200;
268
269
270//-----------------------------------------------------------------------------
271// prepare datafiles
272//-----------------------------------------------------------------------------
273
274 inputPath = SetHostsPaths(false, inputPath );
275
276 TString fitsFileSuffix = ".fits.gz";
277 TString drsFileSuffix = ".drs.fits.gz";
278
279 // declare the data file
280 fits * datafile = NULL;
281
282 // Opens the raw data file and 'binds' the variables given as
283 // Parameters to the data file. So they are filled with
284 // raw data as soon as datafile->GetRow(int) is called.
285
286 //vector of fits files to read
287 vector<TString> fitsFileNames;
288 vector<TString> drsFileNames;
289
290 cout << "...reading sequnezfile" << endl;
291 ReadSequenzFile(
292 sequenzFileName,
293 &fitsFileNames,
294 &drsFileNames
295 );
296
297 cout << "filelist\n" ;
298
299 //fill vector of fits files to read
300 cout << "generating file list:" << endl;
301 if(nRuns < (int)fitsFileNames.size() && nRuns != -1)
302 {
303 fitsFileNames.resize( nRuns);
304 drsFileNames.resize( nRuns);
305 }
306
307 for (unsigned int fileNr = 0; fileNr < fitsFileNames.size(); fileNr++)
308 {
309 if ( fileNr >= drsFileNames.size() )
310 {
311 cout << "too less drs-files in vector, skipping rest" << endl;
312 break;
313
314 drsFileNames.resize(fileNr);
315 fitsFileNames.resize(fileNr);
316 }
317
318 cout << "...checking raw files";
319// cout << "\t fits:\t" << fitsFileNames.at(fileNr);
320// cout << "\n\t\t drs:\t" << drsFileNames.at(fileNr) << endl;
321
322 //filename/path construction
323 TString fits_filename = inputPath;
324 fits_filename.Append(fitsFileNames.at(fileNr));
325 fits_filename.Append(fitsFileSuffix);
326 TString drs_filename = inputPath;
327 drs_filename.Append(drsFileNames.at(fileNr));
328 drs_filename.Append(drsFileSuffix);
329 cout << "\n\t fits:\t" << fits_filename;
330 cout << "\n\t drs:\t" << drs_filename << endl;
331
332 //check if files are readable
333 ifstream fitsFile(fits_filename);
334 ifstream drsFile(drs_filename);
335 if (fitsFile.good() && drsFile.good())
336 {
337 //fill vector
338 cout << "...added\n" << endl;
339 fitsFileNames.at(fileNr) = fits_filename;
340 drsFileNames.at(fileNr) = drs_filename;
341 }
342 else {
343 cout << "...FAILED:";
344 cout << "\t deleting line"
345// << fitsFileNames.at(fileNr)
346// << " and "
347// << drsFileNames.at(fileNr)
348 << " from list\n" << endl;
349
350 fitsFileNames.erase(fitsFileNames.begin()+fileNr);
351 drsFileNames.erase(drsFileNames.begin()+fileNr);
352 fileNr--;
353// continue;
354 }
355 }
356 cout << endl << "size of file vectors: fits: "
357 << fitsFileNames.size()
358 << " drs: "
359 << drsFileNames.size() << endl;
360
361 if (fitsFileNames.size() < 1) return 1;
362
363//----------------------------------------------------------------------------
364// Save-Root-File Settings
365//----------------------------------------------------------------------------
366
367 OutputPath = SetHostsPaths(true, OutputPath );
368 OutRootFileName.Prepend(OutputPath);
369// drsfilename = SetHostsPaths(false, drsfilename );
370// OutRootFileName = SetHostsPaths(true, OutRootFileName );
371 cout << endl;
372// cout << "drs filename:\t" << drsfilename << endl;
373// cout << "data filename:\t" << datafilename << endl;
374 cout << "out filename:\t" << OutRootFileName << endl;
375
376
377 if (saveResults)
378 {
379 CreateRootFile( OutRootFileName, true, verbosityLevel );
380 }
381
382
383
384
385
386//----------------------------------------------------------------------------
387// root global Settings
388//----------------------------------------------------------------------------
389
390 gGainMean = GainMean;
391 gBSLMean = BSLMean;
392 gPixelOverlayXaxisLeft = OverlayWindowLeft;
393 gPixelOverlayXaxisRight = OverlayWindowRight;
394
395 gStyle->SetPalette(1,0);
396 gROOT->SetStyle("Plain");
397// gPad->SetGrid();
398
399//----------------------------------------------------------------------------
400// root Canvas Settings
401//----------------------------------------------------------------------------
402 //Canvas Pad numbering
403 int pulsesCanvasFrameNrs[4] = {
404 1, // Top left
405 2, // Top right
406 3, // bottom left
407 4 // bootom right
408 };
409
410 //Canvas Pad numbering
411 int distributionCanvasFrameNrs[4] = {
412 1, // Top left
413 2, // Top right
414 3, // bottom left
415 4 // bootom right
416 };
417
418 if (ProduceGraphic)
419 {
420
421 //Canvases
422 cgpPixelPulses = new TCanvas*[maxPulseOrder];
423 cgpDistributions = new TCanvas*[maxPulseOrder];
424
425 //TCanvas* gpcDevelopment = NULL;
426 TString cName = "";
427 TString cTitle = "";
428
429 //naming of pulse canvases
430 for (
431 int pulse_order = maxPulseOrder - 1;
432 pulse_order >= 0 ;
433 pulse_order--
434 )
435 {
436 cName ="cgpDistributions";
437 cName += pulse_order;
438 cTitle ="Distributions of Pulses with Order of: ";
439 cTitle += pulse_order;
440 cgpDistributions[pulse_order] = new TCanvas(cName,cTitle, 720,pulse_order*20,720,720);
441 cgpDistributions[pulse_order]->Divide(2, 2);
442 cName ="cgpPixelPulses";
443 cName += pulse_order;
444 cTitle ="Overlays of Pulses with Order of: ";
445 cTitle += pulse_order;
446 cgpPixelPulses[pulse_order] = new TCanvas(cName,cTitle, 0,pulse_order*20,720,720);
447 cgpPixelPulses[pulse_order]->Divide(2, 2);
448 }
449
450
451 // Create (pointer to) Canvases, which are used in every run,
452 // also in 'non-debug' runs
453 // Canvases only need if spike Debug, but I want to deklare
454 // the pointers anyway ...
455 if (spikeDebug)
456 {
457 cFiltered = new TCanvas("cFiltered","filtered DRS Waveforms", 1080,420,360,360);
458 cFiltered->Divide(1, 3);
459 }
460
461 if (testmode)
462 {
463 //additional Test histograms
464 cgpTestHistos = new TCanvas( "cgpTestHistos", "Test Histograms", 360, 420, 360, 360 );
465 cgpTestHistos->Divide(2,0);
466 }
467 }
468
469//-----------------------------------------------------------------------------
470// Filter-Settings
471//-----------------------------------------------------------------------------
472// CFD filter settings
473 int k_cfd = 10;
474 vector<double> a_cfd(k_cfd, 0);
475 double b_cfd = 1.;
476 a_cfd[0]=-0.75;
477 a_cfd[k_cfd-1]=1.;
478
479//-----------------------------------------------------------------------------
480// Book the histograms
481//-----------------------------------------------------------------------------
482
483
484
485 if (testmode)
486 {
487 BookTestHistos( verbosityLevel );
488 }
489//----------------------------------------------------------------------------
490// Initialize Pixel
491//----------------------------------------------------------------------------
492 Pixel** pixel = new Pixel*[NPIX];
493
494 for (int i = 0 ; i < NPIX; i++)
495 {
496 pixel[i] = NULL;
497 }
498
499
500//-----------------------------------------------------------------------------
501// Main Cycle
502//-----------------------------------------------------------------------------
503 cout << "------------------Main cylcle------------------------------" << endl;
504
505 if(npixel == -1) npixel = 1439;
506
507 if ( pixelSetSize == -1 )
508 {
509 pixelSetSize = firstpixel +npixel;
510 }
511 int lastPixelOfSet = 0;
512 cout << "check " << npixel << endl;
513//-------------------------------------
514// Loop over Pixel Sets
515//-------------------------------------
516 for ( int firstPixelOfSet = firstpixel;
517 firstPixelOfSet < firstpixel + npixel;
518 firstPixelOfSet += pixelSetSize )
519 {
520 lastPixelOfSet = firstPixelOfSet + pixelSetSize-1;
521 if (lastPixelOfSet > firstpixel + npixel)
522 {
523 lastPixelOfSet = firstpixel + npixel;
524 }
525 if (lastPixelOfSet > 1439)
526 {
527 lastPixelOfSet = 1439;
528 }
529 if (verbosityLevel == 0)
530 {
531 cout << "------------------------------------------------" << endl
532 << "...processing Set from Pixel "
533 << firstPixelOfSet
534 << " to Pixel "
535 << lastPixelOfSet << endl;
536 }
537
538 //-------------------------------------
539 // Loop over Fits Files
540 //-------------------------------------
541
542 cout << "loop over files" << endl;
543 for( unsigned int fileId = 0; fileId < fitsFileNames.size(); fileId++ )
544 {
545 datafile = NULL;
546// int return_val = 0;
547// return_val = HandleFitsFile(
548// datafile,
549// datafiles.at(fileId),
550// drsfilename,
551// nevents,
552// npixel,
553// verbosityLevel
554// );
555// if (return_val)
556// {
557// cerr << "ERROR" << endl;
558// return return_val;
559// }
560//==================================================
561
562 // Opens the raw data file and 'binds' the variables given as
563 // Parameters to the data file. So they are filled with
564 // raw data as soon as datafile->GetRow(int) is called.
565 cout << "\nopening " << fitsFileNames.at(fileId) << endl;
566
567 gNEvents = openDataFits(
568 fitsFileNames.at(fileId),
569 &datafile,
570 gAllPixelDataVector,
571 gStartCellVector,
572 gCurrentEventID,
573 gRegionOfInterest,
574 gNumberOfPixels,
575 gPXLxROI,
576 verbosityLevel
577 );
578
579 if (gNEvents == 0)
580 {
581 cout << "return code of OpenDataFile:" << fitsFileNames.at(fileId)<< endl;
582 cout << "is zero -> aborting." << endl;
583 return 1;
584 }
585
586 if (verbosityLevel > 0)
587 {
588 cout << endl <<"number of events in file: "<< gNEvents << "\t";
589 }
590
591 if ( nevents == -1 || nevents > gNEvents ) nevents = gNEvents; // -1 means all!
592
593 if (verbosityLevel > 0)
594 {
595 cout <<"of, which "<< nevents << " will be processed"<< endl;
596 cout <<"Total # of Pixel: "<< gNumberOfPixels << "\t";
597 }
598
599 if ( npixel == -1 || npixel > (int)gNumberOfPixels ) npixel = (int)gNumberOfPixels; // -1 means all!
600
601 if (verbosityLevel > 0)
602 {
603 cout <<"of, which "<< npixel << " will be processed"<< endl;
604 }
605
606 //Get the DRS calibration
607 gRC = openCalibFits(
608 drsFileNames.at(fileId),
609 gDrsBaseMean,
610 gDrsGainMean,
611 gDrsTriggeroffsetMean,
612 gTriggerOffsetROI
613 );
614
615 if (gRC == 0)
616 {
617 cout << "return code of openCalibFits:" << drsFileNames.at(fileId) << endl;
618 cout << "is zero -> aborting." << endl;
619 return 1;
620 }
621
622 //==================================================================
623
624//--------------------------------------------------------------------
625// Loop over every Event of Pixel Set
626//--------------------------------------------------------------------
627 for ( int ev = firstevent ; ev < firstevent + nevents; ev++)
628 {
629 // Get an Event --> consists of 1440 Pixel ...erm....data
630 datafile->GetRow( ev - 1 );
631
632 if (verbosityLevel == 1)
633 {
634 cout << "-------------------------------------" << endl
635 << "...processing Set from Pixel "
636 << firstPixelOfSet
637 << " to Pixel "
638 << lastPixelOfSet
639 << "... Event: " << gCurrentEventID
640 << "/" << nevents << endl;
641 }
642
643 //--------------------------------------------------------------------
644 // Loops over every Pixel of a Set of Pixels
645 //--------------------------------------------------------------------
646 for ( int pixelID = firstPixelOfSet;
647 pixelID < lastPixelOfSet + 1
648 && pixelID < firstpixel + npixel;
649 pixelID++ )
650 {
651 if (verbosityLevel > 1)
652 {
653 cout << "-------------------------------------" << endl
654 << "...processing Set from Pixel "
655 << firstPixelOfSet
656 << " to Pixel "
657 << lastPixelOfSet
658 << "... Event: " << gCurrentEventID
659 << "/" << nevents << endl
660 << " Pixel: " << pixelID
661 << "/" << firstpixel + npixel -1 << endl;
662 }
663
664 //-------------------------------------
665 // Create individual Pixel
666 //-------------------------------------
667 if (ev == firstevent && fileId == 0)
668 {
669 pixel[pixelID] = new Pixel(
670 pixelID,
671 maxPulseOrder,
672 verbosityLevel,
673 stats,
674 histoOptions,
675 gPixelOverlayXaxisLeft,
676 gPixelOverlayXaxisRight,
677 gBSLMean,
678 gGainMean
679 );
680 }
681
682 //-------------------------------------
683 // Apply Calibration
684 //-------------------------------------
685 if (verbosityLevel > 2) cout << "...applying DrsCalibration";
686 applyDrsCalibration(
687 gAmeas,
688 pixelID,
689 12,
690 12,
691 gDrsBaseMean,
692 gDrsGainMean,
693 gDrsTriggeroffsetMean,
694 gRegionOfInterest,
695 gAllPixelDataVector,
696 gStartCellVector
697 );
698 if (verbosityLevel > 2) cout << "...done " << endl;
699
700 //-------------------------------------
701 // Apply Filters
702 //-------------------------------------
703 // finds spikes in the raw data, and interpolates the value
704 // spikes are: 1 or 2 slice wide, positive non physical artifacts
705 if (verbosityLevel > 2) cout << "...removeing Spikes";
706 removeSpikes (gAmeas, gVcorr);
707 if (verbosityLevel > 2) cout << "...done " << endl;
708
709 // filter gVcorr with sliding average using FIR filter function
710 if (verbosityLevel > 2) cout << "...applying sliding average filter";
711 sliding_avg(gVcorr, gVslide, avg1);
712 if (verbosityLevel > 2) cout << "...done " << endl;
713
714 // filter gVslide with CFD using FIR filter function
715 if (verbosityLevel > 2) cout << "...apllying factfir filter";
716 factfir(b_cfd , a_cfd, k_cfd, gVslide, gVcfd);
717 if (verbosityLevel > 2) cout << "...done " << endl;
718
719 // filter gVcfd with sliding average using FIR filter function
720 if (verbosityLevel > 2) cout << "...applying 2nd sliding average filter";
721 sliding_avg(gVcfd, gVcfd2, avg2);
722 if (verbosityLevel > 2) cout << "...done " << endl;
723
724 //-------------------------------------
725 // Search vor Zero crossings
726 //-------------------------------------
727 if (verbosityLevel > 2) cout << endl << "...seagRChing zero crossings" ;
728 // peaks in gAmeas[] are found by seagRChing for zero crossings
729 // in gVcfd2
730 // first Argument 1 means ... *rising* edge
731 // second Argument 1 means ... seagRCh with stepsize 1 ... 10 is okay as well
732 vector<Region>* pZXings = zerosearch( gVcfd2 , 1 , 8, verbosityLevel);
733 // pZXings means "zero cross ings"
734
735 float init_reg_size = pZXings->size();
736 float curr_reg_size = init_reg_size;
737 int step = 0;
738
739 pixel[pixelID]->hDiscartedPulses->Fill(step, 1 - curr_reg_size/init_reg_size);
740
741 //enlare the search window
742 curr_reg_size = EnlargeRegion(*pZXings, 10, 10);
743 pixel[pixelID]->hDiscartedPulses->Fill(step++, 1 - curr_reg_size/init_reg_size);
744
745 // find the the max position of a peak
746 curr_reg_size = findAbsMaxInRegions(*pZXings, gVslide);
747 pixel[pixelID]->hDiscartedPulses->Fill(step++, 1 - curr_reg_size/init_reg_size);
748
749// curr_reg_size = calcCFDPositions(*pZXings, gVcfd2);
750// pixel[pixelID]->hDiscartedPulses->Fill(step++, 1 - curr_reg_size/init_reg_size);
751
752
753 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++
754 //do a small after pulse analysis, fill in delayed pulses
755
756 vector<Region>::iterator reg_it = pZXings->begin();
757 int first_pos = pZXings->begin()->maxPos;
758 while( reg_it != pZXings->end() )
759 {
760 //float minimum = 0.0;
761 for (int sl = reg_it->maxPos; sl < 1; sl--)
762 {
763 if (gVcfd2[sl] > 0 && gVcfd2[sl-2] <= 0 && gVcfd2[sl-4] < 0)
764 {
765 //minimum = gVcfd2[sl-2];
766 break;
767 }
768 }
769 pixel[pixelID]->hAfterPulses->Fill(
770 reg_it->maxPos - first_pos,
771 reg_it->maxVal// - minimum
772 );
773 reg_it++;
774 }
775 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++
776
777 // remove peaks below 1 p.e.
778 curr_reg_size = removeMaximaBelow( *pZXings, minPulseAmplitude);
779 pixel[pixelID]->hDiscartedPulses->Fill(step++, 1 - curr_reg_size/init_reg_size);
780
781 // Calculate the position (sample) of a pulses half maximum
782 curr_reg_size = findTimeOfHalfMaxLeft(*pZXings, gVslide, gBSLMean, beginRisingEdge, endRisingEdge, verbosityLevel );
783 pixel[pixelID]->hDiscartedPulses->Fill(step++, 1 - curr_reg_size/init_reg_size);
784
785
786 // Fill arrival times into histograms
787 vector<Region>::iterator it = pZXings->begin();
788 while( it != pZXings->end() )
789 {
790 pixel[pixelID]->hMaxPos->Fill( it->maxPos);
791 pixel[pixelID]->hEdgePos->Fill( it->halfRisingEdgePos);
792 pixel[pixelID]->hEdgeSlope->Fill( it->slopeOfRisingEdge);
793 pixel[pixelID]->hMaxEdgeSlope->Fill( it->maxSlope);
794 pixel[pixelID]->hIntercept->Fill( it->interceptRisingEdge);
795 pixel[pixelID]->hEdgeLength->Fill( it->distanceEdgeToMax);
796 pixel[pixelID]->hMaxAmpl->Fill( it->maxVal);
797 ++it;
798 }
799
800 // remove pulses at the pipeline fringes
801 curr_reg_size = removeRegionWithMaxOnEdge( *pZXings, timeLineEdgeWidth);
802 pixel[pixelID]->hDiscartedPulses->Fill(step++, 1 - curr_reg_size/init_reg_size);
803
804 //sometimes pulses where found with a very flat slope, DISCARD them
805 curr_reg_size = removeRegionWithToFlatSlope(*pZXings, risingEdgeMinSlope);
806 pixel[pixelID]->hDiscartedPulses->Fill(step++, 1 - curr_reg_size/init_reg_size);
807 if (verbosityLevel > 2) cout << "...done" << endl;
808
809 // remove peaks above that are larger tha the largest amplitude to be studied
810 curr_reg_size = removeMaximaAbove( *pZXings, maxPulseAmplitude);
811 pixel[pixelID]->hDiscartedPulses->Fill(step++, 1 - curr_reg_size/init_reg_size);
812
813
814 // remove "Afterpulses" from pulse list
815 curr_reg_size = removeRegionOnFallingEdge( *pZXings, FallingEdgeWidth);
816 pixel[pixelID]->hDiscartedPulses->Fill(step++, 1 - curr_reg_size/init_reg_size);
817
818
819
820 //-----------------------------------------------------------------------------
821 // Fill Overlay Histos
822 //-----------------------------------------------------------------------------
823 if (verbosityLevel > 2)
824 {
825 cout << endl << "Filling Histograms of Pixel# [Set] " << pixelID
826 << " Pixel# [Got] " << pixel[pixelID]->mChid;
827 }
828 FillHistograms(
829 pixel[pixelID],
830 pZXings,
831 AmplWindowWidth,
832 ev,
833 // histoOptions,
834 maxPulseOrder,
835 select_pulses,
836 verbosityLevel
837 );
838 //-----------------------------------------------------------------------------
839 // Spike Debug
840 //-----------------------------------------------------------------------------
841 if ( spikeDebug )
842 {
843 BookDebugHistos(verbosityLevel );
844 gBreakout = ProduceDebugHistos( pZXings );
845 }// end of if(spikeDebug)
846
847 delete pZXings;
848 if (gBreakout) break;
849 //-------------------------------------
850 // Draw 1. Set of Pixel Histograms
851 //-------------------------------------
852 if ((ev % refresh_rate) == 0)
853 {
854 if (ProduceGraphic)
855 {
856 if (debugPixel)
857 {
858 // for ( int order = 0;
859 // order < maxPulseOrder;
860 // order++
861 // )
862 // {
863 pixel[firstPixelOfSet]->DrawOverlayHistograms(
864 cgpPixelPulses,
865 pulsesCanvasFrameNrs
866 );
867
868 pixel[firstPixelOfSet]->DrawDistributionHistograms(
869 cgpDistributions,
870 distributionCanvasFrameNrs
871 );
872 // }
873 if (testmode)
874 {
875 DrawTestHistograms(
876 verbosityLevel
877 );
878 }
879
880 UpdateCanvases(
881 verbosityLevel,
882 maxPulseOrder,
883 testmode
884 );
885 }
886 }
887 }
888
889 if (gBreakout) break;
890
891 if (verbosityLevel > 1)
892 {
893 cout << endl << "...End of Pixel"
894 << endl << "-------------------------------------"<< endl;
895 }
896 }//End of Loop over Set
897 //-------------------------------------
898 // end of Loops over Set
899 //-------------------------------------
900
901 }//End of Loop over Events
902 //-------------------------------------
903 // end of Loops over Events
904 //-------------------------------------
905
906 }//End of Loop over Files
907 //-------------------------------------
908 // end of Loops over Files
909 //-------------------------------------
910
911
912//-------------------------------------
913// Draw Pixel Histogramms of Overlay Spectra
914//-------------------------------------
915
916 for ( int pixelID = firstPixelOfSet;
917 pixelID < firstPixelOfSet + pixelSetSize
918 && pixelID < firstpixel + npixel;
919 pixelID++ )
920 {
921 //here is what happends at the end of each loop over all Events
922
923 if (ProduceGraphic)
924 {
925 pixel[pixelID]->DrawOverlayHistograms(
926 cgpPixelPulses,
927 pulsesCanvasFrameNrs
928 );
929
930 pixel[pixelID]->DrawDistributionHistograms(
931 cgpDistributions,
932 distributionCanvasFrameNrs
933 );
934
935 UpdateCanvases(
936 verbosityLevel,
937 maxPulseOrder,
938 testmode
939 );
940 }
941
942 //Save Histograms of Pixels into Output rootfile
943 if (verbosityLevel > 2)
944 {
945 cout << endl << "Saving Pixel# [Set] " << pixelID
946 << " Pixel# [Got] " << pixel[pixelID]->mChid;
947 }
948 pixel[pixelID]->SavePixelHistograms( OutRootFileName, saveResults );
949
950
951 if (debugPixel)
952 {
953 //Process gui events asynchronously during input
954 cout << endl;
955 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
956 timer.TurnOn();
957 TString input = Getline("Type 'q' to exit, <return> to go on: ");
958 timer.TurnOff();
959 if (input=="q\n") {
960 break;
961 }
962 }
963
964 if (verbosityLevel > 2)
965 {
966 cout << endl << "Deleting Pixel# [Set] " << pixelID
967 << " Pixel# [Got] " << pixel[pixelID]->mChid;
968 }
969 delete pixel[pixelID];
970
971 if (verbosityLevel > 2)
972 {
973 cout << endl << "...End of Pixel-Set"
974 << endl << "------------------------------------------------"
975 << endl;
976 }
977 }
978 }
979 // End of Loop over all Pixels
980
981//-------------------------------------
982// Draw Histograms
983//-------------------------------------
984 TString test = "root";
985 SaveList( //save histograms of generell results into output root file
986 OutRootFileName,
987 test,
988 hRootList,
989 saveResults,
990 verbosityLevel
991 );
992
993 if (ProduceGraphic)
994 {
995 UpdateCanvases(
996 verbosityLevel,
997 maxPulseOrder,
998 testmode
999 );
1000
1001 //Process gui events asynchronously during input
1002 cout << endl;
1003 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
1004 timer.TurnOn();
1005 Getline("Press <return> to go on");
1006 timer.TurnOff();
1007
1008 DeletePixelCanvases(
1009 maxPulseOrder,
1010 verbosityLevel
1011 );
1012
1013 if (testmode)
1014 delete cgpTestHistos;
1015
1016 if (spikeDebug)
1017 delete cFiltered;
1018 }
1019
1020 delete datafile;
1021 delete[] pixel;
1022
1023 return( 0 );
1024}
1025//----------------------------------------------------------------------------
1026// end of main function
1027//-----------------------------------------------------------------------------
1028
1029
1030
1031
1032//-----------------------------------------------------------------------------
1033// Funktions
1034//-----------------------------------------------------------------------------
1035
1036/*
1037ich erzeuge sowas wie 36 pixel mit ca 20 Histogrammen pro pixel
1038wie viel speicher braucht das?
1039Ich brauche eine möglichkeit jedes dieser histogramme zu identifizieren
1040am besten ein array oder eine liste oder einen abstracten datentyp mit den histogrammen als attribute.
1041ne klasse wäre super
1042
1043*/
1044
1045void BookTestHistos( int verbosityLevel )
1046{
1047 if (verbosityLevel > 2) cout << endl << "...book pixel histograms" << endl;
1048
1049 hTesthisto = new TH1F (
1050 "hTesthisto",
1051 "Deviation of rising edge and maximum",
1052 600,
1053 -10.1,
1054 10.1
1055 );
1056 hTesthisto->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
1057 hTesthisto->GetYaxis()->SetTitle( "counts" );
1058
1059 hTesthisto2 = new TH2F (
1060 "hTesthisto2",
1061 "Deviation of rising edge and maximum by event #",
1062 gNEvents,
1063 250,
1064 800,
1065 600,
1066 -10.1,
1067 10.1
1068 );
1069// hTesthisto2->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
1070// hTesthisto2->GetYaxis()->SetTitle( "counts" );
1071// hTesthisto2->SetMarkerStyle( 2 );
1072
1073 if (verbosityLevel > 2) cout << "...done" << endl;
1074}
1075//end of BookTestHistos
1076//----------------------------------------------------------------------------
1077
1078
1079void
1080BookDebugHistos( int verbosityLevel )
1081{
1082 if (verbosityLevel > 2) cout << endl << "...book histograms" << endl;
1083 hRootList = new TList();
1084 debugHistos = new TH1F[ Number_Of_Debug_Histo_Types ];
1085 for ( int type = 0; type < Number_Of_Debug_Histo_Types; type++){
1086 debugHistos[ type ].SetBins(1024, 0, 1024);
1087 debugHistos[ type ].SetLineColor(1);
1088 debugHistos[ type ].SetLineWidth(2);
1089 debugHistos[ type ].SetStats(false);
1090
1091 // set X axis paras
1092 debugHistos[ type ].GetXaxis()->SetLabelSize(0.08);
1093 debugHistos[ type ].GetXaxis()->SetTitleSize(0.08);
1094 debugHistos[ type ].GetXaxis()->SetTitleOffset(1.0);
1095 debugHistos[ type ].GetXaxis()->SetTitle(Form("Time slice [%.1f ns/slice]", 1./2.));
1096
1097 // set Y axis paras
1098 debugHistos[ type ].GetYaxis()->SetLabelSize(0.08);
1099 debugHistos[ type ].GetYaxis()->SetTitleSize(0.08);
1100 debugHistos[ type ].GetYaxis()->SetTitleOffset(0.3);
1101 debugHistos[ type ].GetYaxis()->SetTitle("Amplitude [mV]");
1102// hRootList->Add( &debugHistos[ type ] );
1103 }
1104 if (verbosityLevel > 2) cout << "...done" << endl;
1105}
1106// end of BookDebugHistos
1107//----------------------------------------------------------------------------
1108
1109
1110void
1111FillHistograms(
1112 Pixel* CurrentPixel,
1113 vector<Region>* pZXings,
1114 int AmplWindowWidth,
1115 int eventNumber,
1116 int maxPulseOrder,
1117 bool select_pulses,
1118 int verbosityLevel
1119 )
1120{
1121 //entry counters
1122 int maxOverlayEntries = 0;
1123 int edgeOverlayEntries = 0;
1124 int maxProfileEntries = 0;
1125 int edgeProfileEntries = 0;
1126
1127
1128
1129
1130 if (verbosityLevel > 2) cout << endl << "...filling pulse histograms" ;
1131 if (verbosityLevel > 3) cout << endl << "...EventNR " << eventNumber ;
1132 bool use_this_peak=false;
1133 int order_of_pulse=0;
1134 vector<Region>::iterator reg;
1135 //Loop over all found zerocrossings in Region
1136 for (reg = pZXings->begin() ; reg < pZXings->end() ; reg++)
1137 {
1138 //skip those who are at the Rim of the ROI
1139 if (reg->maxPos < 12 || (unsigned int) reg->maxPos > gRegionOfInterest-12)
1140 {
1141 if (verbosityLevel > 2) cout << endl << "\t...out of range" << endl;
1142 continue;
1143 }
1144
1145 // Define axis range of Histogram
1146 int Left = reg->maxPos - gPixelOverlayXaxisLeft;
1147 int Right = reg->maxPos + gPixelOverlayXaxisRight;
1148 int EdgeLeft = reg->halfRisingEdgePos - gPixelOverlayXaxisLeft;
1149 int EdgeRight = reg->halfRisingEdgePos + gPixelOverlayXaxisRight;
1150
1151 if (Left < 0) Left = 0;
1152 if (EdgeLeft < 0) EdgeLeft = 0;
1153 if (Right > (int)gVcorr.size() ) Right = gVcorr.size();
1154 if (EdgeRight > (int)gVcorr.size() ) EdgeRight = gVcorr.size();
1155
1156
1157 if (verbosityLevel > 2) cout << endl << "\t...choosing Histogram" ;
1158 //determine order of pulse and which histogram shall be filled
1159 for ( int order = 0; order < maxPulseOrder; order++)
1160 {
1161 use_this_peak = UseThisPulse(
1162 reg->maxPos,
1163 order,
1164 AmplWindowWidth,
1165 gGainMean,
1166 gBSLMean,
1167 maxPulseOrder,
1168 verbosityLevel
1169 );
1170 if ( use_this_peak)
1171 {
1172 order_of_pulse = order;
1173 break;
1174 }
1175 }
1176 if (select_pulses){
1177 order_of_pulse = 0;
1178 use_this_peak = true;
1179 }
1180 //Fill histograms
1181 if (use_this_peak)
1182 {
1183 //get number of entries of used histograms
1184 maxOverlayEntries = CurrentPixel->hMaxOverlay[order_of_pulse]->GetEntries();
1185 edgeOverlayEntries = CurrentPixel->hEdgeOverlay[order_of_pulse]->GetEntries();
1186 maxProfileEntries = CurrentPixel->hMaxProfile[order_of_pulse]->GetEntries();
1187 edgeProfileEntries = CurrentPixel->hEdgeProfile[order_of_pulse]->GetEntries();
1188
1189
1190 //Histograms for Distributions
1191 if (CurrentPixel->mOptions.Contains("S") )
1192 {
1193 if (verbosityLevel > 2)
1194 cout << endl << "\t...filling Edge Histogram" ;
1195 CurrentPixel->hSlopeRisingEdge[order_of_pulse]->Fill( reg->slopeOfRisingEdge ) ;
1196 }
1197
1198 if (CurrentPixel->mOptions.Contains("R") )
1199 {
1200 if (verbosityLevel > 2)
1201 cout << endl << "\t...filling RisingEdgeToMax Histogram" ;
1202 CurrentPixel->hRisingEdgeToMax[order_of_pulse]->Fill( reg->distanceEdgeToMax ) ;
1203 }
1204
1205 if (CurrentPixel->mOptions.Contains("M") )
1206 {
1207 if (verbosityLevel > 2)
1208 cout << endl << "\t...filling PosOfMax Histogram" ;
1209 CurrentPixel->hPosOfMax[order_of_pulse]->Fill( reg->maxPos ) ;
1210 }
1211
1212 float average_of_max = gVcorr[reg->maxPos];
1213 average_of_max += gVcorr[reg->maxPos + 1];
1214 average_of_max += gVcorr[reg->maxPos - 1];
1215 average_of_max /= 3;
1216
1217 float slope_of_follow = 0.0;
1218
1219 //Histograms for Maximum Overlay
1220 for ( int pos = Left; pos < Right; pos++)
1221 {
1222 if (pos > reg->maxPos + 1)
1223 {
1224 if ( pos+5 >= Right )
1225 {
1226 break;
1227 }
1228 if (gVcorr[pos+1] > average_of_max)
1229 {
1230 break;
1231 }
1232 slope_of_follow = (gVcorr[pos+5] - gVcorr[pos+1])/5;
1233 if (slope_of_follow > 1)
1234 {
1235 break;
1236 }
1237 }
1238 CurrentPixel->hMaxOverlay[order_of_pulse]->Fill( pos - (reg->maxPos), gVcorr[pos]) ;
1239 CurrentPixel->hMaxProfile[order_of_pulse]->Fill( pos - (reg->maxPos), gVcorr[pos]) ;
1240 }
1241
1242 //Histograms for Edge Overlay
1243 int last_pos = 0;
1244 bool ap_cutoff = false;
1245 slope_of_follow = 0.0;
1246
1247 for ( int pos = EdgeLeft; pos < EdgeRight; pos++)
1248 {
1249 if (pos > reg->maxPos + 1)
1250 {
1251 last_pos = pos;
1252 // check if the current slice is close to the end of the overlay window
1253 if ( pos+5 >= Right )
1254 {
1255 ap_cutoff = false;
1256 break;
1257 }
1258 // check if the next slice is larger than the pulse maximum
1259 // i guess it makes more sense to compare the average amplitude of the slices before
1260 //to the following pulse
1261 if (gVcorr[pos+1] > average_of_max)
1262 {
1263 ap_cutoff = true;
1264 break;
1265 }
1266
1267 // check if the next 5 slices have a steep slope
1268 slope_of_follow = (gVcorr[pos+5] - gVcorr[pos+1])/5;
1269 if (slope_of_follow > 1)
1270 {
1271 ap_cutoff = true;
1272 break;
1273 }
1274 }
1275// if (pos > reg->maxPos
1276// && gVcorr[pos] < average_of_max
1277// && gVcorr[pos+3] > average_of_max
1278// && gVcorr[pos+10] > average_of_max
1279// )
1280// {
1281// break;
1282// }
1283 CurrentPixel->hEdgeOverlay[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), gVcorr[pos]) ;
1284 CurrentPixel->hEdgeProfile[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), gVcorr[pos]) ;
1285 }
1286
1287 CurrentPixel->hPulseLength->Fill(last_pos - (reg->halfRisingEdgePos));
1288 if (ap_cutoff)
1289 {
1290 CurrentPixel->hPulseLengthAPcutoff->Fill(last_pos - (reg->halfRisingEdgePos));
1291 }
1292 else
1293 {
1294 CurrentPixel->hPulseLengthTLcutoff->Fill(last_pos - (reg->halfRisingEdgePos));
1295 }
1296
1297 //set histogram entries
1298 CurrentPixel->hMaxOverlay[order_of_pulse]->SetEntries( maxOverlayEntries +1 );
1299 CurrentPixel->hEdgeOverlay[order_of_pulse]->SetEntries( edgeOverlayEntries +1 );
1300 CurrentPixel->hMaxProfile[order_of_pulse]->SetEntries( maxProfileEntries +1 );
1301 CurrentPixel->hEdgeProfile[order_of_pulse]->SetEntries( edgeProfileEntries +1 );
1302
1303 if (verbosityLevel > 2) cout << "...done" << endl;
1304 }
1305
1306 //gBreakout option
1307 if (gBreakout)
1308 {
1309 break;
1310 }
1311 }
1312 if (verbosityLevel > 2) cout << "...done" << endl;
1313}
1314// end of FillHistograms
1315//----------------------------------------------------------------------------
1316
1317bool UseThisPulse(
1318 int maxPos,
1319 int order,
1320 int AmplWindowWidth,
1321 float gainMean,
1322 float bslMean,
1323 int maxPulseOrder,
1324 int verbosityLevel
1325 )
1326{
1327 //determine if pulse is of given order
1328bool use_this_peak = false;
1329 if (gVcorr[maxPos] >= ((gainMean*(order+1)) - bslMean - (AmplWindowWidth/2))
1330 && gVcorr[maxPos] < ((gainMean*(order+1)) -bslMean + AmplWindowWidth/2))
1331 {
1332 if (verbosityLevel > 3) cout << "...#" << order ;
1333 use_this_peak = true;
1334 }
1335 else
1336 {
1337 if (verbosityLevel > 3)
1338 {
1339 if (order >= (maxPulseOrder - 1))
1340 {
1341 cout << "...NONE" << endl ;
1342 }
1343 }
1344 use_this_peak = false;
1345 }
1346
1347return use_this_peak;
1348}
1349
1350
1351void
1352DrawTestHistograms(
1353 int verbosityLevel
1354 )
1355{
1356 if (verbosityLevel > 2) cout << endl << "...drawing pulse histograms TODO CANVAS" ;
1357
1358// cgpTestHistos->cd(1);
1359// hTesthisto->Draw();
1360// cgpTestHistos->cd(2);
1361// hTesthisto2->Draw("BOX");
1362}
1363// end of DrawTestHistograms
1364//----------------------------------------------------------------------------
1365
1366void
1367UpdateCanvases(
1368 int verbosityLevel,
1369 int max_pulse_order,
1370 bool testmode
1371 )
1372{
1373 if (verbosityLevel > 3) cout << endl << "...updating canvases" ;
1374 for (int pulse_order = 0; pulse_order < max_pulse_order; pulse_order++)
1375 {
1376 cgpPixelPulses[pulse_order]->Modified();
1377 cgpPixelPulses[pulse_order]->Update();
1378 cgpDistributions[pulse_order]->Modified();
1379 cgpDistributions[pulse_order]->Update();
1380 if ( testmode )
1381 {
1382 cgpTestHistos->Modified();
1383 cgpTestHistos->Update();
1384 }
1385 }
1386}
1387// end of UpdateCanvases
1388//----------------------------------------------------------------------------
1389
1390void
1391DeletePixelCanvases(
1392 int maxPulseOrder,
1393 int verbosityLevel
1394 )
1395{
1396 if (verbosityLevel > 2) cout << endl << "...delete pixel Canvases" ;
1397 for (int pulse_order = 0; pulse_order < maxPulseOrder; pulse_order++ )
1398 {
1399 delete cgpPixelPulses[pulse_order];
1400 delete cgpDistributions[pulse_order];
1401 }
1402 delete[] cgpPixelPulses;
1403 delete[] cgpDistributions;
1404}
1405
1406//----------------------------------------------------------------------------
1407//----------------------------------------------------------------------------
1408
1409bool
1410ProduceDebugHistos(
1411 vector<Region>* pZXings
1412 )
1413{
1414 // TODO do this correct. The vectors should be the rigt ones... this is just luck
1415 debugHistos[gAmeas_].GetXaxis()->Set(gAmeas.size() , -0.5 , gAmeas.size()-0.5);
1416 debugHistos[gVcorr_].GetXaxis()->Set(gAmeas.size() , -0.5 , gAmeas.size()-0.5);
1417 debugHistos[gVslide_].GetXaxis()->Set(gAmeas.size() , -0.5 , gAmeas.size()-0.5);
1418 debugHistos[gVcfd_].GetXaxis()->Set(gAmeas.size() , -0.5 , gAmeas.size()-0.5);
1419 debugHistos[gVcfd2_].GetXaxis()->Set(gAmeas.size() , -0.5 , gAmeas.size()-0.5);
1420
1421 for ( unsigned int sl = 0; sl < gRegionOfInterest; sl++)
1422 {
1423 debugHistos[gAmeas_].SetBinContent(sl, gAmeas[sl]);
1424 debugHistos[gVcorr_].SetBinContent(sl, gVcorr[sl]);
1425 debugHistos[gVslide_].SetBinContent( sl, gVslide[sl] );
1426 debugHistos[gVcfd_].SetBinContent( sl, gVcfd[sl] );
1427 debugHistos[gVcfd2_].SetBinContent( sl, gVcfd2[sl] );
1428 }
1429
1430 bool has_negative_slope = false;
1431 cFiltered->cd(1);
1432 gPad->SetGrid();
1433 debugHistos[gAmeas_].Draw();
1434
1435 TF1 *OneEdge;
1436 TLine *OneLine;
1437 TBox *OneBox;
1438 vector<TF1*> MyEdges;
1439 vector<TBox*> MyBoxes;
1440 vector<TBox*> MyLines;
1441
1442 int left = debugHistos[gAmeas_].GetXaxis()->GetFirst();
1443 int right = debugHistos[gAmeas_].GetXaxis()->GetLast();
1444 vector<Region>::iterator region;
1445 for (region = pZXings->begin() ; region < pZXings->end() ; region++)
1446 {
1447 if (region->slopeOfRisingEdge < 0)
1448 {
1449 has_negative_slope = true;
1450 }
1451 }
1452
1453 if (has_negative_slope)
1454 {
1455 for (unsigned int i=0; i<pZXings->size(); i++){
1456 Region* reg = &pZXings->at(i);
1457 cout << "Slope is: " << reg->slopeOfRisingEdge << endl;
1458
1459// xBegin = (reg->halfRisingEdgePos - reg->distanceEdgeToMax);
1460// yBegin = reg->integRCeptRisingEdge + reg->slopeOfRisingEdge * (xBegin)
1461
1462 OneEdge = new TF1(
1463 "OneEdge",
1464 "[0]*x+[1]",
1465 left,
1466 right
1467 );
1468
1469 OneEdge->SetParameters(
1470 reg->slopeOfRisingEdge,
1471 reg->interceptRisingEdge
1472 );
1473
1474 OneEdge->SetLineColor(kRed);
1475 OneEdge->SetLineWidth(1);
1476 MyEdges.push_back(OneEdge);
1477
1478 OneEdge->Draw("SAME");
1479// delete OneEdge;
1480
1481 OneBox = new TBox(
1482 pZXings->at(i).maxPos -10 ,
1483 pZXings->at(i).maxVal -0.5,
1484 pZXings->at(i).maxPos +10 ,
1485 pZXings->at(i).maxVal +0.5);
1486 OneBox->SetLineColor(kBlue);
1487 OneBox->SetLineWidth(1);
1488 OneBox->SetFillStyle(0);
1489 OneBox->SetFillColor(kRed);
1490 MyBoxes.push_back(OneBox);
1491 OneBox->Draw();
1492
1493 OneLine = new TLine(
1494 pZXings->at(i).halfRisingEdgePos,
1495 gAmeas[pZXings->at(i).halfRisingEdgePos],
1496 pZXings->at(i).halfRisingEdgePos,
1497 pZXings->at(i).halfRisingEdgePos * pZXings->at(i).slopeOfRisingEdge + pZXings->at(i).interceptRisingEdge
1498 );
1499 OneLine->SetLineColor(kBlue);
1500 OneLine->SetLineWidth(3);
1501 MyLines.push_back(OneBox);
1502 OneLine->Draw();
1503 }
1504
1505 cFiltered->cd(2);
1506 gPad->SetGrid();
1507 debugHistos[gVslide_].Draw();
1508 debugHistos[gVcorr_].Draw("SAME");
1509
1510 for (unsigned int i=0; i<MyEdges.size(); i++){
1511 MyEdges.at(i)->Draw("SAME");
1512 }
1513 for (unsigned int i=0; i<MyBoxes.size(); i++){
1514 MyBoxes.at(i)->Draw();
1515 }
1516 for (unsigned int i=0; i<MyLines.size(); i++){
1517 MyLines.at(i)->Draw();
1518 }
1519
1520 cFiltered->cd(3);
1521 gPad->SetGrid();
1522 debugHistos[gVcorr_].SetLineColor(kBlue);
1523 debugHistos[gVcorr_].Draw();
1524 debugHistos[gAmeas_].Draw("SAME");
1525 debugHistos[gVslide_].SetLineColor(kRed);
1526 debugHistos[gVslide_].Draw("SAME");
1527 debugHistos[gVcfd_].SetLineColor(kGreen);
1528 debugHistos[gVcfd_].Draw("SAME");
1529 debugHistos[gVcfd2_].SetLineColor(kViolet);
1530 debugHistos[gVcfd2_].Draw("SAME");
1531 for (unsigned int i=0; i<MyEdges.size(); i++){
1532 MyEdges.at(i)->Draw("SAME");
1533 }
1534 for (unsigned int i=0; i<MyBoxes.size(); i++){
1535 MyBoxes.at(i)->Draw();
1536 }
1537 for (unsigned int i=0; i<MyLines.size(); i++){
1538 MyLines.at(i)->Draw();
1539 }
1540
1541
1542// cFiltered->cd(3);
1543// gPad->SetGrid();
1544// debugHistos[gVcfd2_].Draw();
1545// TLine *zeroline = new TLine(0, 0, 1024, 0);
1546// zeroline->SetLineColor(kBlue);
1547// zeroline->Draw();
1548
1549 cFiltered->Update();
1550
1551
1552 //Process gui events asynchronously during input
1553 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
1554 timer.TurnOn();
1555 TString input = Getline("Type 'q' to exit, <return> to go on: ");
1556 timer.TurnOff();
1557 if (input=="q\n") {
1558 return true;
1559 }
1560
1561 MyLines.clear();
1562 MyBoxes.clear();
1563 MyEdges.clear();
1564
1565// for (unsigned int i=0; i<pZXings->size(); i++)
1566// {
1567// delete MyEdges.at(i)
1568// delete MyBoxes.at(i);
1569// }
1570// delete OneBox;
1571// delete OneEdge;
1572// delete OneLine;
1573
1574// MyEdges.erase();
1575// MyBoxes.erase();
1576 }
1577
1578 //TODO!!!!!!!!!
1579 // do some Garbage collection ...
1580 // all the Objects on the heap should be deleted here.
1581return false;
1582}// end of if(spikeDebug)
1583
1584int HandleFitsFile(
1585 fits * datafile,
1586 TString datafilename,
1587 TString drsfilename,
1588 int nevents,
1589 int npixel,
1590 int verbosityLevel
1591 )
1592{
1593 // Opens the raw data file and 'binds' the variables given as
1594 // Parameters to the data file. So they are filled with
1595 // raw data as soon as datafile->GetRow(int) is called.
1596 cout << "\nopening " << datafilename << endl;
1597
1598 gNEvents = openDataFits(
1599 datafilename,
1600 &datafile,
1601 gAllPixelDataVector,
1602 gStartCellVector,
1603 gCurrentEventID,
1604 gRegionOfInterest,
1605 gNumberOfPixels,
1606 gPXLxROI,
1607 verbosityLevel
1608 );
1609
1610 if (gNEvents == 0)
1611 {
1612 cout << "return code of OpenDataFile:" << datafilename<< endl;
1613 cout << "is zero -> aborting." << endl;
1614 return 1;
1615 }
1616
1617 if (verbosityLevel > 0)
1618 {
1619 cout << endl <<"number of events in file: "<< gNEvents << "\t";
1620 }
1621
1622 if ( nevents == -1 || nevents > gNEvents ) nevents = gNEvents; // -1 means all!
1623
1624 if (verbosityLevel > 0)
1625 {
1626 cout <<"of, which "<< nevents << " will be processed"<< endl;
1627 cout <<"Total # of Pixel: "<< gNumberOfPixels << "\t";
1628 }
1629
1630 if ( npixel == -1 || npixel > (int)gNumberOfPixels ) npixel = (int)gNumberOfPixels; // -1 means all!
1631
1632 if (verbosityLevel > 0)
1633 {
1634 cout <<"of, which "<< npixel << " will be processed"<< endl;
1635 }
1636
1637 //Get the DRS calibration
1638 gRC = openCalibFits(
1639 drsfilename,
1640 gDrsBaseMean,
1641 gDrsGainMean,
1642 gDrsTriggeroffsetMean,
1643 gTriggerOffsetROI
1644 );
1645
1646 if (gRC == 0)
1647 {
1648 cout << "return code of openCalibFits:" << drsfilename << endl;
1649 cout << "is zero -> aborting." << endl;
1650 return 1;
1651 }
1652 return 0;
1653}
1654
1655void ReadSequenzFile(TString fileName, vector<TString> *fitsFiles, vector<TString> *drsFiles)
1656{
1657 ifstream seqFile;
1658 seqFile.open(fileName,ios::in);
1659
1660 TString fitsName, drsName;
1661
1662 if(!seqFile)
1663 {
1664 cerr << "***** no config-file in the user defined directory.\n"
1665 << "Check the path .\n" << endl;
1666 return;
1667 }
1668
1669 if(!seqFile.good()) {
1670 cerr << "***** not able to read config-file." << endl;
1671 return;
1672 }
1673
1674 while(!seqFile.eof())
1675 {
1676 seqFile >> fitsName >> drsName;
1677
1678 if (!fitsName.IsWhitespace() && !fitsName.IsNull())
1679 {
1680 fitsFiles->push_back(fitsName);
1681
1682 if (!drsName.IsWhitespace() && !drsName.IsNull())
1683 {
1684 drsFiles->push_back(drsName);
1685 }
1686 else drsFiles->push_back("0");
1687 }
1688 }
1689
1690 seqFile.close();
1691
1692 return;
1693}
1694
1695
1696//int main(int argc,char *argv[])
1697//{
1698
1699// TString test;
1700// TString gRCFileName;
1701// TString processType = "overlay";
1702// bool gRCFileNameCmdSet = false;
1703// int verbLevel = 0; // different verbosity levels can be implemented here
1704// bool verbLevelCmdSet = false;
1705// bool save = false;
1706// bool produceGraphic = false;
1707
1708// // decode arguments
1709// if(argc < 2)
1710// {
1711// printf("no arguements given, using standard arguments\n");
1712// }
1713
1714// // set conditions for functions arguments
1715// for (int i=1;i<argc;i++)
1716// {
1717// test = argv[i];
1718
1719// if (test.Contains("--config_file") || test.Contains("-c"))
1720// {
1721// cout << "gRC-File: \"" << argv[i + 1] << "\"" << endl;
1722// gRCFileName = argv[i + 1];
1723// gRCFileNameCmdSet = true;
1724// continue;
1725// }
1726
1727// if (test.Contains("--verbosity") || test.Contains("-v"))
1728// {
1729// cout << "Verbosity Level: \"" << argv[i + 1] << "\"" << endl;
1730// verbLevel = atoi(argv[i + 1]);
1731// continue;
1732// }
1733
1734// if (test.Contains("--save") || test.Contains("-s"))
1735// {
1736// cout << "will save results" << endl;
1737// save = true;
1738// continue;
1739// }
1740
1741// if (test.Contains("--graphics") || test.Contains("-g"))
1742// {
1743// cout << "will produce graphics" << endl;
1744// produceGraphic = true;
1745// continue;
1746// }
1747// }
1748
1749// // reading gRC-File:
1750// configfile gRCfile( gRCFileName, processType );
1751
1752// if (save)
1753// {
1754// gRCfile.mSave = save;
1755// }
1756
1757// if (verbLevelCmdSet)
1758// {
1759// gRCfile.mVerbLevel = verbLevel;
1760// }
1761
1762// if (produceGraphic)
1763// {
1764// gRCfile.mProduceGraphic = produceGraphic;
1765// }
1766
1767// FPulseOverlay(
1768// gRCfile.mInputPath,
1769// gRCfile.mDataFileName,
1770// gRCfile.mDrsFileName,
1771// gRCfile.mOutputPath,
1772// gRCfile.mOutputFile,
1773// gRCfile.mFirstEvent,
1774// gRCfile.mNumEvents,
1775// gRCfile.mFirstPixel,
1776// gRCfile.mNumPixel,
1777// gRCfile.mPixelSetSize,
1778// gRCfile.mMaxOrder,
1779// gRCfile.mVerbLevel,
1780// gRCfile.mHistoOptions,
1781// gRCfile.mAmplWindowWidth,
1782// gRCfile.mGainMean,
1783// gRCfile.mBSLMean,
1784// gRCfile.mAvg1,
1785// gRCfile.mAvg2,
1786// gRCfile.mOverlayWindowLeft,
1787// gRCfile.mOverlayWindowRight,
1788// gRCfile.mProduceGraphic,
1789// gRCfile.mSpikeDebug,
1790// gRCfile.mDbgPixel,
1791// gRCfile.mPrintStats,
1792// gRCfile.mSave,
1793// gRCfile.mTestMode,
1794// gRCfile.mRefreshRate
1795// );
1796// return 0;
1797
1798//}
Note: See TracBrowser for help on using the repository browser.