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

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