source: fact/tools/rootmacros/FPulseTemplate.C@ 13428

Last change on this file since 13428 was 13422, checked in by Jens Buss, 14 years ago
now with median, mean and max of slice
File size: 66.4 KB
Line 
1/********************** FPulseTemplate ***********************
2* read FACT raw data
3* remove spikes
4* calculate baseline
5* find peaks (CFD and normal discriminator)
6* compute pulse height and pulse integral spektrum of the peaks
7+ search for Peaks in data
8+ put peaks into different histos depending und Amplitude
9+ draw histos for single, double, tripple, ....above 100mV Photonpulses
10+ draw Parameterdevelopment depending on photon quantity
11+ form a template (shape depending on number of photons)
12 so it can be used for detection of other peaks
13*************************************************************/
14
15//----------------------------------------------------------------------------
16// root libraries
17//----------------------------------------------------------------------------
18
19#include <TROOT.h>
20#include <TCanvas.h>
21#include <TProfile.h>
22#include <TTimer.h>
23#include <TH1F.h>
24#include <TH2F.h>
25#include <Getline.h>
26#include <TLine.h>
27#include <TBox.h>
28#include <TMath.h>
29#include <TFile.h>
30#include <TStyle.h>
31#include <TString.h>
32#include <TProfile.h>
33
34#include <stdio.h>
35#include <stdint.h>
36#include <cstdio>
37#include <iostream>
38#include <fstream>
39using namespace std;
40
41#define NPIX 1440
42#define NCELL 1024
43#define FAD_MAX_SAMPLES 1024
44#define MAX_PULS_ORDER 1
45#define HAVE_ZLIB
46
47//----------------------------------------------------------------------------
48// rootmacros
49//----------------------------------------------------------------------------
50
51#include "fits.h"
52
53#include "openFits.h"
54#include "openFits.c"
55
56#include "discriminator.h"
57#include "discriminator.C"
58#include "zerosearch.h"
59#include "zerosearch.C"
60#include "factfir.C"
61
62#include "DrsCalibration.C"
63#include "DrsCalibration.h"
64
65#include "SpikeRemoval.h"
66#include "SpikeRemoval.C"
67
68//----------------------------------------------------------------------------
69// Decleration of global variables
70//----------------------------------------------------------------------------
71
72bool breakout=false;
73
74int gNEvents;
75vector<int16_t> AllPixelDataVector;
76vector<int16_t> StartCellVector;
77unsigned int CurrentEventID;
78size_t PXLxROI;
79UInt_t gRegionOfInterest;
80UInt_t NumberOfPixels;
81TString histotitle;
82
83size_t TriggerOffsetROI, RC;
84size_t drs_n;
85vector<float> drs_basemean;
86vector<float> drs_gainmean;
87vector<float> drs_triggeroffsetmean;
88
89vector<float> Ameas(FAD_MAX_SAMPLES); // copy of the data (measured amplitude
90vector<float> Vcorr(FAD_MAX_SAMPLES); // corrected Values
91vector<float> Vslide(FAD_MAX_SAMPLES); // sliding average result
92vector<float> Vcfd(FAD_MAX_SAMPLES); // CDF result
93vector<float> Vcfd2(FAD_MAX_SAMPLES); // CDF result + 2nd sliding average
94
95float gGainMean = 9;
96float gBSLMean = 0;
97
98typedef struct{
99 double maxAmpl;
100 int countOfMax;
101 } OverlayMaximum;
102
103//typedef struct{
104// TString name;
105// TString title;
106// TString xTitle;
107// TString yTitle;
108// float yMax;
109// float yMin;
110//} histAttrib_t;
111
112//histAttrib_t* gHistoAttrib[MAX_PULS_ORDER];
113//histAttrib_t* gMaximumAttrib[MAX_PULS_ORDER];
114//histAttrib_t* gMaxGausAttrib[MAX_PULS_ORDER];
115//histAttrib_t* gProfileAttrib[MAX_PULS_ORDER];
116
117
118// histograms
119const int Number_Of_Debug_Histo_Types = 7;
120
121const unsigned int
122 Ameas_ = 0,
123 N1mean_ = 1,
124 Vcorr_ = 2,
125 Vtest_ = 3,
126 Vslide_ = 4,
127 Vcfd_ = 5,
128 Vcfd2_ = 6;
129
130//const char* gHistoTypes[8] = {
131// "hPixelOverlay",
132// "hPixelMax",
133// "hPixelEdgeOverlay",
134// "hPixelProfile",
135// "hAllPixelOverlay",
136// "hAllPixelMax",
137// "hAllPixelMaxGaus",
138// "hAllPixelProfile"
139//};
140
141//----------------------------------------------------------------------------
142// Initialisation of histograms
143//----------------------------------------------------------------------------
144
145// Temporary Objects
146TH1F* debugHistos = NULL;
147TH2F* hOverlayTemp = NULL;
148TH1F* hMaximumTemp = NULL;
149TH1D* hProjPeak = NULL;
150TH1F* hTesthisto = NULL;
151TH2F* hTesthisto2 = NULL;
152
153//Pixelwise Histograms
154TH2F* hPixelOverlay[MAX_PULS_ORDER]= {NULL};//histogrammm for overlay of detected Peaks
155TH2F* hPixelEdgeOverlay[MAX_PULS_ORDER] = {NULL};
156TProfile* hPixelProfile[MAX_PULS_ORDER] = {NULL};//histogrammm for Profile of detected Peaks
157TProfile* hPixelProfile2[MAX_PULS_ORDER] = {NULL};//histogrammm for Profile of detected Peaks
158TH1F* hPixelMax[MAX_PULS_ORDER] = {NULL};
159TH1F* hPixelMedian[MAX_PULS_ORDER] = {NULL};
160TH1F* hPixelMean[MAX_PULS_ORDER] = {NULL};
161
162//All Pixel Histograms
163TH2F* hAllPixelOverlay[MAX_PULS_ORDER] = {NULL};
164TProfile* hAllPixelProfile[MAX_PULS_ORDER] = {NULL};
165TProfile* hAllPixelProfile2[MAX_PULS_ORDER] = {NULL};
166TH1F* hAllPixelMax[MAX_PULS_ORDER] = {NULL};
167TH1F* hAllPixelMedian[MAX_PULS_ORDER] = {NULL};
168TH1F* hAllPixelMean[MAX_PULS_ORDER] = {NULL};
169TH2F* hAllPixelEdgeOverlay[MAX_PULS_ORDER] = {NULL};
170
171//Histogram Parameters
172Int_t gPixelOverlayXaxisLeft = 0;
173Int_t gPixelOverlayXaxisRight = 0;
174
175//Root-File Objects
176TObjArray* hList = NULL;
177TObjArray* hAllPixelList = NULL;
178TObjArray* hRootList = NULL;
179
180//Canvases
181TCanvas* cgpPixelPulses[MAX_PULS_ORDER] = {NULL};
182TCanvas* cgpAllPixelPulses[MAX_PULS_ORDER] = {NULL};
183TCanvas* cgpTestHistos = NULL;
184//TCanvas* gpcDevelopment = NULL;
185
186//----------------------------------------------------------------------------
187// Functions
188//----------------------------------------------------------------------------
189
190void BookHistos(int );
191void BookPixelHistos(int );
192void BookTestHistos( int );
193void DeletePixelHistos(int );
194void SaveHistograms( const char*, const char*, TObjArray*, int );
195void FillHistograms(vector<Region>*, int, int, int);
196void PlotPulseShapeOfMaxPropability(unsigned int, int, int);
197
198void DrawPulseHistograms(int, int);
199void FitMaxPropabilityPuls( TH1F* , int );
200void DrawMaximumHistograms( int, int );
201void DrawTestHistograms( int);
202Double_t MedianOfH1( TH1 *h1 );
203void PlotMedianEachSliceOfPulse(unsigned int, int, int);
204
205
206void UpdateCanvases( int, int);
207void DeletePixelCanvases( int );
208
209void CreateRootFile(const char*, int );
210TFile *OpenRootFile(const char*, int );
211void CloseRootFile(TFile *tf);
212TString CreateSubDirName(int); //creates a String containing the path to the subdirectory in root file
213TString CreateSubDirName(const char* );
214void WritePixelTemplateToCsv( TString, const char*, const char*, int, int);
215void WriteAllPixelTemplateToCsv( TString, const char*, const char*, int);
216//void SetHistogrammSettings( const char*, int, int , int);
217//void CreatePulseCanvas( TCanvas* , unsigned int, const char* , const char*, const char*, int);
218//void DeletePulseCanvas( TCanvas* , unsigned int);
219
220
221//----------------------------------------------------------------------------
222//----------------------------------------------------------------------------
223// MAIN - Funtion
224//----------------------------------------------------------------------------
225//----------------------------------------------------------------------------
226int FPulseTemplate(
227 char* datafilename = "../data/2011/11/09/20111109_006.fits.gz",
228 const char* drsfilename = "../data/2011/11/09/20111109_003.drs.fits.gz",
229 const char* OutRootFileName = "../analysis/FPulseTemplate/20111109_006/20111109_006.pulses.root",
230 const char* OutPutPath = "../analysis/FPulseTemplate/20111109_006/Templates/",
231 bool ProduceGraphic = false,
232 int refresh_rate = 500, //refresh rate for canvases
233 bool spikeDebug = false,
234 bool debugPixel = false,
235 bool fitdata = false,
236 int verbosityLevel = 0, // different verbosity levels can be implemented here
237 int firstevent = 0,
238 int nevents = -1,
239 int firstpixel = 0,
240 int npixel = -1,
241 int AmplWindowWidth = 14, //Width of Window for selection of pulses to histograms
242 float GainMean = 8,
243 float BSLMean = -1, //4 Histogramms will be drawn, decide how far pulsehights differ from eachother
244 int avg1 = 8,
245 int avg2 = 8,
246 int OverlayWindowLeft = 70,
247 int OverlayWindowRight = 230
248 )
249{
250 gGainMean = GainMean;
251 gBSLMean = BSLMean;
252//----------------------------------------------------------------------------
253// Save-Root-File Settings
254//----------------------------------------------------------------------------
255 CreateRootFile( OutRootFileName, verbosityLevel );
256
257
258//----------------------------------------------------------------------------
259// root global Settings
260//----------------------------------------------------------------------------
261
262 gPixelOverlayXaxisLeft = OverlayWindowLeft;
263 gPixelOverlayXaxisRight = OverlayWindowRight;
264
265 gStyle->SetPalette(1,0);
266 gROOT->SetStyle("Plain");
267// gPad->SetGrid();
268
269 for (
270 int pulse_order = MAX_PULS_ORDER - 1;
271 pulse_order >= 0 ;
272 pulse_order--
273 )
274 {
275 TString cName ="cgpPixelPulses";
276 cName += pulse_order;
277 TString cTitle ="Pulses of Order: ";
278 cTitle += pulse_order;
279 cgpPixelPulses[pulse_order] = new TCanvas(cName,cTitle, 0,pulse_order*20,800,800);
280 cgpPixelPulses[pulse_order]->Divide(2, 2);
281 cName = "cgpAllPixelPulses";
282 cName += pulse_order;
283 cTitle ="AllPixel average of puls shapes of Order: ";
284 cTitle += pulse_order;
285 cgpAllPixelPulses[pulse_order] = new TCanvas( cName, cTitle, 801, pulse_order*20, 800, 800 );
286 cgpAllPixelPulses[pulse_order]->Divide(2, 2);
287 }
288
289 // Create (pointer to) Canvases, which are used in every run,
290 // also in 'non-debug' runs
291 // Canvases only need if spike Debug, but I want to deklare
292 // the pointers anyway ...
293
294 TCanvas *cFiltered = NULL;
295 if (spikeDebug)
296 {
297 cFiltered = new TCanvas("cFiltered","filtered DRS Waveforms", 1,310,400,300);
298 cFiltered->Divide(1, 3);
299 }
300
301 cgpTestHistos = new TCanvas( "cgpTestHistos", "Test Histograms", 801, 0, 800, 800 );
302 cgpTestHistos->Divide(2,0);
303//-----------------------------------------------------------------------------
304// Filter-Settings
305//-----------------------------------------------------------------------------
306// CFD filter settings
307 int k_cfd = 10;
308 vector<double> a_cfd(k_cfd, 0);
309 double b_cfd = 1.;
310 a_cfd[0]=-0.75;
311 a_cfd[k_cfd-1]=1.;
312
313//-----------------------------------------------------------------------------
314// prepare datafile
315//-----------------------------------------------------------------------------
316
317// Open the data file
318
319 fits * datafile;
320 // Opens the raw data file and 'binds' the variables given as
321 // Parameters to the data file. So they are filled with
322 // raw data as soon as datafile->GetRow(int) is called.
323 gNEvents = openDataFits(
324 datafilename,
325 &datafile,
326 AllPixelDataVector,
327 StartCellVector,
328 CurrentEventID,
329 gRegionOfInterest,
330 NumberOfPixels,
331 PXLxROI,
332 verbosityLevel
333 );
334
335 if (gNEvents == 0)
336 {
337 cout << "return code of OpenDataFile:" << datafilename<< endl;
338 cout << "is zero -> aborting." << endl;
339 return 1;
340 }
341
342 if (verbosityLevel > 0)
343 {
344 cout << endl <<"number of events in file: "<< gNEvents << "\t";
345 }
346
347 if ( nevents == -1 || nevents > gNEvents ) nevents = gNEvents; // -1 means all!
348
349 if (verbosityLevel > 0)
350 {
351 cout <<"of, which "<< nevents << " will be processed"<< endl;
352 cout <<"Total # of Pixel: "<< NumberOfPixels << "\t";
353 }
354
355 if ( npixel == -1 || npixel > (int)NumberOfPixels ) npixel = (int)NumberOfPixels; // -1 means all!
356
357 if (verbosityLevel > 0)
358 {
359 cout <<"of, which "<< npixel << " will be processed"<< endl;
360 }
361
362 //Get the DRS calibration
363 RC = openCalibFits( drsfilename, drs_basemean, drs_gainmean, drs_triggeroffsetmean, TriggerOffsetROI);
364 if (RC == 0)
365 {
366 cout << "return code of openCalibFits:" << drsfilename << endl;
367 cout << "is zero -> aborting." << endl;
368 return 1;
369 }
370
371 // Book the histograms
372 BookHistos(verbosityLevel );
373 BookTestHistos( verbosityLevel );
374//-------------------------------------
375// Loop over every Pixel
376//-------------------------------------
377 for ( int pixel = firstpixel; pixel < firstpixel + npixel; pixel++ )
378 {
379 if (verbosityLevel >= 0)
380 {
381 cout << "------------------------------------------------" << endl
382 << "...processing Pixel: " << pixel << endl;
383 }
384
385 BookPixelHistos(verbosityLevel );
386
387//--------------------------------------------------------------------
388// Loops over Every Event of Pixel
389//--------------------------------------------------------------------
390 for ( int ev = firstevent; ev < firstevent + nevents; ev++)
391 {
392 // Get an Event --> consists of 1440 Pixel ...erm....data
393 datafile->GetRow( ev );
394
395 if (verbosityLevel > 0)
396 {
397 cout << "-------------------------------------" << endl
398 << "...processing Event: " << CurrentEventID
399 << "/" << nevents << " of Pixel: " << pixel << endl;
400 }
401
402//-------------------------------------
403// Apply Calibration
404//-------------------------------------
405 if (verbosityLevel > 2) cout << "...applying DrsCalibration";
406 applyDrsCalibration(
407 Ameas,
408 pixel,
409 12,
410 12,
411 drs_basemean,
412 drs_gainmean,
413 drs_triggeroffsetmean,
414 gRegionOfInterest,
415 AllPixelDataVector,
416 StartCellVector
417 );
418 if (verbosityLevel > 2) cout << "...done " << endl;
419
420//-------------------------------------
421// Apply Filters
422//-------------------------------------
423 // finds spikes in the raw data, and interpolates the value
424 // spikes are: 1 or 2 slice wide, positive non physical artifacts
425 if (verbosityLevel > 2) cout << "...removeing Spikes";
426 removeSpikes (Ameas, Vcorr);
427 if (verbosityLevel > 2) cout << "...done " << endl;
428
429 // filter Vcorr with sliding average using FIR filter function
430 if (verbosityLevel > 2) cout << "...applying sliding average filter";
431 sliding_avg(Vcorr, Vslide, avg1);
432 if (verbosityLevel > 2) cout << "...done " << endl;
433
434 // filter Vslide with CFD using FIR filter function
435 if (verbosityLevel > 2) cout << "...apllying factfir filter";
436 factfir(b_cfd , a_cfd, k_cfd, Vslide, Vcfd);
437 if (verbosityLevel > 2) cout << "...done " << endl;
438
439 // filter Vcfd with sliding average using FIR filter function
440 if (verbosityLevel > 2) cout << "...applying 2nd sliding average filter";
441 sliding_avg(Vcfd, Vcfd2, avg2);
442 if (verbosityLevel > 2) cout << "...done " << endl;
443
444//-------------------------------------
445// Search vor Zero crossings
446//-------------------------------------
447 if (verbosityLevel > 2) cout << endl << "...searching zero crossings" ;
448 // peaks in Ameas[] are found by searching for zero crossings
449 // in Vcfd2
450 // first Argument 1 means ... *rising* edge
451 // second Argument 1 means ... search with stepsize 1 ... 10 is okay as well
452 vector<Region>* pZXings = zerosearch( Vcfd2 , 1 , 8);
453 // pZXings means "zero cross ings"
454 EnlargeRegion(*pZXings, 10, 10);
455 findAbsMaxInRegions(*pZXings, Vslide);
456 removeMaximaBelow( *pZXings, 3.0);
457 removeRegionWithMaxOnEdge( *pZXings, 2);
458 removeRegionOnFallingEdge( *pZXings, 100);
459 findTimeOfHalfMaxLeft(*pZXings, Vcorr, gBSLMean, 5, 10, verbosityLevel );
460 if (verbosityLevel > 2) cout << "...done" << endl;
461
462//-----------------------------------------------------------------------------
463// Fill Overlay Histos
464//-----------------------------------------------------------------------------
465 FillHistograms(pZXings, AmplWindowWidth, ev, verbosityLevel );
466
467//-----------------------------------------------------------------------------
468// Spike Debug
469//-----------------------------------------------------------------------------
470 if ( spikeDebug )
471 {
472 // TODO do this correct. The vectors should be the rigt ones... this is just luck
473 debugHistos[Ameas_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
474 debugHistos[Vcorr_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
475 debugHistos[Vslide_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
476 debugHistos[Vcfd_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
477 debugHistos[Vcfd2_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
478
479 for ( unsigned int sl = 0; sl < gRegionOfInterest; sl++)
480 {
481 debugHistos[Ameas_].SetBinContent(sl, Ameas[sl]);
482 debugHistos[Vcorr_].SetBinContent(sl, Vcorr[sl]);
483 debugHistos[Vslide_].SetBinContent( sl, Vslide[sl] );
484 debugHistos[Vcfd_].SetBinContent( sl, Vcfd[sl] );
485 debugHistos[Vcfd2_].SetBinContent( sl, Vcfd2[sl] );
486 }
487
488
489 cFiltered->cd(1);
490 gPad->SetGrid();
491 debugHistos[Ameas_].Draw();
492
493 cFiltered->cd(2);
494 gPad->SetGrid();
495 debugHistos[Vcorr_].Draw();
496
497 cFiltered->cd(3);
498 gPad->SetGrid();
499 debugHistos[Vslide_].Draw();
500
501 TBox *OneBox;
502 vector<TBox*> MyBoxes;
503 for (unsigned int i=0; i<pZXings->size(); i++){
504 OneBox = new TBox(
505 pZXings->at(i).maxPos -10 ,
506 pZXings->at(i).maxVal -0.5,
507 pZXings->at(i).maxPos +10 ,
508 pZXings->at(i).maxVal +0.5);
509 OneBox->SetLineColor(kBlue);
510 OneBox->SetLineWidth(1);
511 OneBox->SetFillStyle(0);
512 OneBox->SetFillColor(kRed);
513 MyBoxes.push_back(OneBox);
514 OneBox->Draw();
515 }
516
517// cFiltered->cd(3);
518// gPad->SetGrid();
519// debugHistos[Vcfd2_].Draw();
520// TLine *zeroline = new TLine(0, 0, 1024, 0);
521// zeroline->SetLineColor(kBlue);
522// zeroline->Draw();
523
524 cFiltered->Update();
525
526
527 //Process gui events asynchronously during input
528 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
529 timer.TurnOn();
530 TString input = Getline("Type 'q' to exit, <return> to go on: ");
531 timer.TurnOff();
532 if (input=="q\n") {
533 breakout=true;
534 }
535
536 //TODO!!!!!!!!!
537 // do some Garbage collection ...
538 // all the Objects on the heap should be deleted here.
539
540 }// end of if(spikeDebug)
541
542 delete pZXings;
543 if (breakout) break;
544//-------------------------------------
545// Draw 1. Set of Pixel Histograms
546//-------------------------------------
547 if ((ev % refresh_rate) == 0)
548 {
549 DrawPulseHistograms(
550 verbosityLevel,
551 MAX_PULS_ORDER
552 );
553 DrawTestHistograms(
554 verbosityLevel
555 );
556 if (ProduceGraphic)
557 {
558 UpdateCanvases(
559 verbosityLevel,
560 MAX_PULS_ORDER
561 );
562 }
563 }
564
565 if (breakout) break;
566
567 if (verbosityLevel > 2)
568 {
569 cout << endl << "-------------------------------------"<< endl;
570 }
571 }//End of Loop over Events
572//-------------------------------------
573// end of Loops over Events
574//-------------------------------------
575
576
577
578
579//-------------------------------------
580// Histogramms of Maximas in Overlay Spectra
581//-------------------------------------
582
583 PlotPulseShapeOfMaxPropability(
584 MAX_PULS_ORDER,
585 fitdata,
586 verbosityLevel
587 );
588
589 PlotMedianEachSliceOfPulse(
590 MAX_PULS_ORDER,
591 fitdata,
592 verbosityLevel
593 );
594
595 WritePixelTemplateToCsv(
596 OutPutPath,
597 "PulseTemplate_PointSet",
598 "Maximum",
599 pixel,
600 verbosityLevel
601 );
602
603 if (verbosityLevel > 2) cout << "...done" << endl;
604 //here is what happends at the end of each loop over all pixels
605 //Save Histograms of Pixel into Output rootfile
606 DrawPulseHistograms(
607 verbosityLevel,
608 MAX_PULS_ORDER
609 );
610
611 DrawMaximumHistograms(
612 verbosityLevel,
613 MAX_PULS_ORDER
614 );
615
616 if (ProduceGraphic)
617 {
618 UpdateCanvases(
619 verbosityLevel,
620 MAX_PULS_ORDER
621 );
622 }
623
624// SaveHistograms(
625// OutRootFileName,
626// CreateSubDirName(pixel),
627// hList,
628// verbosityLevel
629// );
630
631 if (debugPixel)
632 {
633 //Process gui events asynchronously during input
634 cout << endl;
635 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
636 timer.TurnOn();
637 TString input = Getline("Type 'q' to exit, <return> to go on: ");
638 timer.TurnOff();
639 if (input=="q\n") {
640 break;
641 }
642 }
643
644 DeletePixelHistos(verbosityLevel);
645
646 if (verbosityLevel > 2)
647 {
648 cout << endl << "...End of Pixel"
649 << endl << "------------------------------------------------"
650 << endl;
651 }
652 }
653 // End of Loop over all Pixels
654
655//-------------------------------------
656// Draw Histograms
657//-------------------------------------
658
659 SaveHistograms( //save histograms of all pixel into output root file
660 OutRootFileName,
661 CreateSubDirName("All"),
662 hAllPixelList,
663 verbosityLevel
664 );
665
666// SaveHistograms( //save histograms of generell results into output root file
667// OutRootFileName,
668// "root",
669// hRootList,
670// verbosityLevel
671// );
672
673
674
675 if (ProduceGraphic)
676 {
677 UpdateCanvases(
678 verbosityLevel,
679 MAX_PULS_ORDER
680 );
681 }
682
683 WriteAllPixelTemplateToCsv(
684 OutPutPath,
685 "PulseTemplate_PointSet",
686 "Maximum",
687 verbosityLevel
688 );
689
690 DeletePixelCanvases( verbosityLevel );
691 return( 0 );
692}
693//----------------------------------------------------------------------------
694// end of main function
695//-----------------------------------------------------------------------------
696
697
698
699
700//-----------------------------------------------------------------------------
701// Funktions
702//-----------------------------------------------------------------------------
703
704
705void
706BookPixelHistos( int verbosityLevel)
707{
708 if (verbosityLevel > 2) cout << endl << "...book pixel histograms" << endl;
709 hList = new TObjArray;
710 TString histoname;
711 for (int order = 0; order < MAX_PULS_ORDER; order ++)
712 {
713// SetHistogramAttributes(
714// "PeakOverlay",
715// order,
716// gOverlayAttrib,
717// MaxAmplOfFirstPulse
718// );
719 histoname = "hPixelOverlay";
720 histoname += order;
721 if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
722 hPixelOverlay[order] = new TH2F(
723 histoname,
724 "Overlay of detected pulses of one pulse order for one Pixel",
725 gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
726 (-1*gPixelOverlayXaxisLeft)-0.5,
727 gPixelOverlayXaxisRight-0.5 ,
728 512,
729 -55.5,
730 200.5
731 );
732 hPixelOverlay[order]->SetAxisRange(
733 gBSLMean - 5,
734 (gGainMean*(order+1)) + 10,
735 "Y");
736 hPixelOverlay[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
737 hPixelOverlay[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
738 //hPixelProfile->SetBit(TH2F::kCanRebin);
739 hList->Add( hPixelOverlay[order] );
740
741 //------------------------------------------------------------------------
742// SetHistogramAttributes(
743// "PeakMaximum",
744// order,
745// gMaximumAttrib,
746// MaxAmplOfFirstPulse
747// );
748 histoname = "hPixelMax";
749 histoname += order;
750 if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
751 hPixelMax[order] = new TH1F(
752 histoname,
753 "Maximum value of each slice in overlay plot",
754 gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
755 (-1*gPixelOverlayXaxisLeft)-0.5,
756 gPixelOverlayXaxisRight-0.5
757 );
758 hPixelMax[order]->SetAxisRange(
759 gBSLMean - 5,
760 (gGainMean*(order+1)) + 10,
761 "Y");
762 hPixelMax[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
763 hPixelMax[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
764 hList->Add( hPixelMax[order] );
765
766//------------------------------------------------------------------------
767 // SetHistogramAttributes(
768 // "PeakMaximum",
769 // order,
770 // gMaximumAttrib,
771 // MaxAmplOfFirstPulse
772 // );
773 histoname = "hPixelMedian";
774 histoname += order;
775 if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
776 hPixelMedian[order] = new TH1F(
777 histoname,
778 "Median of each slice in overlay plot",
779 gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
780 (-1*gPixelOverlayXaxisLeft)-0.5,
781 gPixelOverlayXaxisRight-0.5
782 );
783 hPixelMedian[order]->SetAxisRange(
784 gBSLMean - 5,
785 (gGainMean*(order+1)) + 10,
786 "Y");
787 hPixelMedian[order]->SetLineColor(kRed);
788 hPixelMedian[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
789 hPixelMedian[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
790 hList->Add( hPixelMedian[order] );
791
792//------------------------------------------------------------------------
793 // SetHistogramAttributes(
794 // "PeakMaximum",
795 // order,
796 // gMaximumAttrib,
797 // MaxAmplOfFirstPulse
798 // );
799 histoname = "hPixelMean";
800 histoname += order;
801 if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
802 hPixelMean[order] = new TH1F(
803 histoname,
804 "Mean of each slice in overlay plot",
805 gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
806 (-1*gPixelOverlayXaxisLeft)-0.5,
807 gPixelOverlayXaxisRight-0.5
808 );
809 hPixelMean[order]->SetAxisRange(
810 gBSLMean - 5,
811 (gGainMean*(order+1)) + 10,
812 "Y");
813 hPixelMean[order]->SetLineColor(kBlue);
814 hPixelMean[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
815 hPixelMean[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
816 hList->Add( hPixelMean[order] );
817
818//------------------------------------------------------------------------
819// SetHistogramAttributes(
820// "PeakMaxGaus",
821// order,
822// gMaxGausAttrib,
823// MaxAmplOfFirstPulse
824// );
825 histoname = "hPixelEdgeOverlay";
826 histoname += order;
827 if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
828 hPixelEdgeOverlay[order] = new TH2F(
829 histoname,
830 "Overlay at rising edge of detected pulses of one pulse order for one Pixel ",
831 gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
832 (-1*gPixelOverlayXaxisLeft)-0.5,
833 gPixelOverlayXaxisRight-0.5 ,
834 512,
835 -55.5,
836 200.5
837 );
838
839 hPixelEdgeOverlay[order]->SetAxisRange(
840 gBSLMean - 5,
841 (gGainMean*(order+1)) + 10,
842 "Y");
843 hPixelEdgeOverlay[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
844 hPixelEdgeOverlay[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
845 hList->Add( hPixelEdgeOverlay[order] );
846
847 //------------------------------------------------------------------------
848// SetHistogramAttributes(
849// "PeakProfile",
850// order,
851// gProfileAttrib,
852// MaxAmplOfFirstPulse
853// );
854 histoname = "hPixelProfile";
855 histoname += order;
856 if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
857 hPixelProfile[order] = new TProfile(
858 histoname,
859 "Mean value of each slice in overlay plot (Tprofile)",
860 gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,//nbinsx
861 (-1*gPixelOverlayXaxisLeft)-0.5, //xlow
862 gPixelOverlayXaxisRight-0.5 , //xup
863 // 512, //nbinsy
864 -55.5, //ylow
865 300.5, //yup
866 "s"); //option
867 hPixelProfile[order]->SetAxisRange(
868 gBSLMean - 5,
869 (gGainMean*(order+1)) + 10,
870 "Y");
871 hPixelProfile[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
872 hPixelProfile[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
873 //hPixelProfile->SetBit(TH2F::kCanRebin);
874 hList->Add( hPixelProfile[order] );
875
876
877 //------------------------------------------------------------------------
878 // SetHistogramAttributes(
879 // "PeakProfile",
880 // order,
881 // gProfileAttrib,
882 // MaxAmplOfFirstPulse
883 // );
884 histoname = "hPixelProfile2";
885 histoname += order;
886 if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
887 hPixelProfile2[order] = new TProfile(
888 histoname,
889 "Mean value of each slice in overlay plot (Tprofile)",
890 gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,//nbinsx
891 (-1*gPixelOverlayXaxisLeft)-0.5, //xlow
892 gPixelOverlayXaxisRight-0.5 , //xup
893 // 512, //nbinsy
894 -55.5, //ylow
895 300.5, //yup
896 "s"); //option
897 hPixelProfile2[order]->SetLineColor(kRed);
898 hPixelProfile2[order]->SetAxisRange(
899 gBSLMean - 5,
900 (gGainMean*(order+1)) + 10,
901 "Y");
902 hPixelProfile2[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
903 hPixelProfile2[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
904 //hPixelProfile->SetBit(TH2F::kCanRebin);
905
906
907
908
909 }
910 if (verbosityLevel > 2) cout << "...done" << endl;
911}
912//end of BookPixelHistos
913//----------------------------------------------------------------------------
914
915
916void DeletePixelHistos( int verbosityLevel )
917{
918 if (verbosityLevel > 2) cout << endl << "...delete current pixel histograms" ;
919 for (int order = 0; order < MAX_PULS_ORDER; order ++)
920 {
921 if (verbosityLevel > 3) cout << endl << "...deleting hPixelOverlay" << order;
922 delete hPixelOverlay[order];
923 if (verbosityLevel > 3) cout << endl << "...deleting hPixelMax" << order;
924 delete hPixelMax[order];
925 if (verbosityLevel > 3) cout << endl << "...deleting hPixelMedian" << order;
926 delete hPixelMedian[order];
927 if (verbosityLevel > 3) cout << endl << "...deleting hPixelMean" << order;
928 delete hPixelMean[order];
929 if (verbosityLevel > 3) cout << endl << "...deleting hPixelEdgeOverlay" << order;
930 delete hPixelEdgeOverlay[order];
931 if (verbosityLevel > 3) cout << endl << "...deleting hPixelProfile" << order;
932 delete hPixelProfile[order];
933 if (verbosityLevel > 3) cout << endl << "...deleting hPixelProfile2" << order;
934 delete hPixelProfile2[order];
935 }
936 if (verbosityLevel > 3) cout << endl << "...deleting hList";
937 delete hList;
938 if (verbosityLevel > 2) cout << endl << "...done" << endl;
939}
940//end of DeletePixelHistos
941//----------------------------------------------------------------------------
942
943void BookTestHistos( int verbosityLevel )
944{
945 if (verbosityLevel > 2) cout << endl << "...book pixel histograms" << endl;
946
947 hTesthisto = new TH1F (
948 "hTesthisto",
949 "Deviation of rising edge and maximum",
950 600,
951 -10.1,
952 10.1
953 );
954 hTesthisto->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
955 hTesthisto->GetYaxis()->SetTitle( "counts" );
956
957 hTesthisto2 = new TH2F (
958 "hTesthisto2",
959 "Deviation of rising edge and maximum by event #",
960 gNEvents,
961 250,
962 800,
963 600,
964 -10.1,
965 10.1
966 );
967// hTesthisto2->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
968// hTesthisto2->GetYaxis()->SetTitle( "counts" );
969// hTesthisto2->SetMarkerStyle( 2 );
970
971 if (verbosityLevel > 2) cout << "...done" << endl;
972}
973//end of BookTestHistos
974//----------------------------------------------------------------------------
975
976void BookHistos( int verbosityLevel )
977{
978 if (verbosityLevel > 2) cout << endl << "...book histograms" << endl;
979 hAllPixelList = new TObjArray;
980 hRootList = new TObjArray;
981 debugHistos = new TH1F[ Number_Of_Debug_Histo_Types ];
982 for ( int type = 0; type < Number_Of_Debug_Histo_Types; type++){
983 debugHistos[ type ].SetBins(1024, 0, 1024);
984 debugHistos[ type ].SetLineColor(1);
985 debugHistos[ type ].SetLineWidth(2);
986 debugHistos[ type ].SetStats(false);
987
988 // set X axis paras
989 debugHistos[ type ].GetXaxis()->SetLabelSize(0.08);
990 debugHistos[ type ].GetXaxis()->SetTitleSize(0.08);
991 debugHistos[ type ].GetXaxis()->SetTitleOffset(1.0);
992 debugHistos[ type ].GetXaxis()->SetTitle(Form("Time slice [%.1f ns/slice]", 1./2.));
993
994 // set Y axis paras
995 debugHistos[ type ].GetYaxis()->SetLabelSize(0.08);
996 debugHistos[ type ].GetYaxis()->SetTitleSize(0.08);
997 debugHistos[ type ].GetYaxis()->SetTitleOffset(0.3);
998 debugHistos[ type ].GetYaxis()->SetTitle("Amplitude [mV]");
999 }
1000
1001 // All Pixel Histograms
1002 //------------------------------------------------------------------------
1003 TString histoname;
1004 for (int order = 0; order < MAX_PULS_ORDER; order ++)
1005 {
1006// SetHistogramAttributes(
1007// "AllPixelPeakOverlay",
1008// order,
1009// gOverlayAttrib,
1010// MaxAmplOfFirstPulse
1011// );
1012 histoname = "hAllPixelOverlay";
1013 histoname += order;
1014 if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
1015 hAllPixelOverlay[order] = new TH2F(
1016 histoname,
1017 "Overlay of detected Pulses of one Order for All Pixels",
1018 gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
1019 (-1*gPixelOverlayXaxisLeft)-0.5,
1020 gPixelOverlayXaxisRight-0.5 ,
1021 512,
1022 -55.5,
1023 300.5
1024 );
1025 hAllPixelOverlay[order]->SetAxisRange(
1026 gBSLMean - 5,
1027 (gGainMean*(order+1)) + 10,
1028 "Y");
1029 hAllPixelOverlay[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
1030 hAllPixelOverlay[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
1031 //hAllPixelOverlay->SetBit(TH2F::kCanRebin);
1032 hAllPixelList->Add( hAllPixelOverlay[order] );
1033
1034 //------------------------------------------------------------------------
1035
1036// SetHistogramAttributes(
1037// "AllPixelPeakMax",
1038// order,
1039// gMaximumAttrib,
1040// MaxAmplOfFirstPulse
1041// );
1042 histoname = "hAllPixelMax";
1043 histoname += order;
1044 if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
1045 hAllPixelMax[order] = new TH1F(
1046 histoname,
1047 "Maximum value of each slice in overlay plot for All Pixels",
1048 gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
1049 (-1*gPixelOverlayXaxisLeft)-0.5,
1050 gPixelOverlayXaxisRight-0.5
1051 );
1052 hAllPixelMax[order]->SetAxisRange(
1053 gBSLMean - 5,
1054 (gGainMean*(order+1)) + 10,
1055 "Y");
1056 hAllPixelMax[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
1057 hAllPixelMax[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
1058 hAllPixelList->Add( hAllPixelMax[order] );
1059
1060//------------------------------------------------------------------------
1061// SetHistogramAttributes(
1062// "PeakMaximum",
1063// order,
1064// gMaximumAttrib,
1065// MaxAmplOfFirstPulse
1066// );
1067 histoname = "hAllPixelMedian";
1068 histoname += order;
1069 if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
1070 hAllPixelMedian[order] = new TH1F(
1071 histoname,
1072 "Median of each slice in overlay plot for All Pixels",
1073 gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
1074 (-1*gPixelOverlayXaxisLeft)-0.5,
1075 gPixelOverlayXaxisRight-0.5
1076 );
1077 hAllPixelMedian[order]->SetAxisRange(
1078 gBSLMean - 5,
1079 (gGainMean*(order+1)) + 10,
1080 "Y");
1081 hAllPixelMedian[order]->SetLineColor(kRed);
1082 hAllPixelMedian[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
1083 hAllPixelMedian[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
1084 hAllPixelList->Add( hAllPixelMedian[order] );
1085 if (verbosityLevel > 3) cout << "...BREAKPOINT in " << histoname << endl;
1086//------------------------------------------------------------------------
1087// SetHistogramAttributes(
1088// "PeakMaximum",
1089// order,
1090// gMaximumAttrib,
1091// MaxAmplOfFirstPulse
1092// );
1093 histoname = "hAllPixelMean";
1094 histoname += order;
1095 if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
1096 hAllPixelMean[order] = new TH1F(
1097 histoname,
1098 "Mean of each slice in overlay plot for All Pixels",
1099 gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
1100 (-1*gPixelOverlayXaxisLeft)-0.5,
1101 gPixelOverlayXaxisRight-0.5
1102 );
1103 hAllPixelMean[order]->SetAxisRange(
1104 gBSLMean - 5,
1105 (gGainMean*(order+1)) + 10,
1106 "Y");
1107 hAllPixelMean[order]->SetLineColor(kBlue);
1108 hAllPixelMean[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
1109 hAllPixelMean[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
1110 hAllPixelList->Add( hAllPixelMean[order] );
1111
1112//------------------------------------------------------------------------
1113// SetHistogramAttributes(
1114// "PeakMaxGaus",
1115// order,
1116// gMaxGausAttrib,
1117// MaxAmplOfFirstPulse
1118// );
1119 histoname = "hAllPixelEdgeOverlay";
1120 histoname += order;
1121 if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
1122 hAllPixelEdgeOverlay[order] = new TH2F(
1123 histoname,
1124 "Overlay at rising edge of detected pulses of one pulse order for All Pixels ",
1125 gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
1126 (-1*gPixelOverlayXaxisLeft)-0.5,
1127 gPixelOverlayXaxisRight-0.5 ,
1128 512,
1129 -55.5,
1130 200.5
1131 );
1132
1133 hAllPixelEdgeOverlay[order]->SetAxisRange(
1134 gBSLMean - 5,
1135 (gGainMean*(order+1)) + 10,
1136 "Y");
1137 hAllPixelEdgeOverlay[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
1138 hAllPixelEdgeOverlay[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
1139 hAllPixelList->Add( hAllPixelEdgeOverlay[order] );
1140
1141//------------------------------------------------------------------------
1142// SetHistogramAttributes(
1143// "AllPixelPeakProfile",
1144// order,
1145// gProfileAttrib,
1146// MaxAmplOfFirstPulse
1147// );
1148 histoname = "hAllPixelProfile";
1149 histoname += order;
1150 if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
1151 hAllPixelProfile[order] = new TProfile(
1152 histoname,
1153 "Mean value of each slice in overlay plot for All Pixels (Tprofile)",
1154 gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,//nbinsx
1155 (-1*gPixelOverlayXaxisLeft)-0.5, //xlow
1156 gPixelOverlayXaxisRight-0.5 , //xup
1157 // 512, //nbinsy
1158 -55.5, //ylow
1159 300.5, //yup
1160 "s"); //option
1161 hAllPixelProfile[order]->SetAxisRange(
1162 gBSLMean - 5,
1163 (gGainMean*(order+1)) + 10,
1164 "Y");
1165 hAllPixelProfile[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
1166 hAllPixelProfile[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
1167 //hPixelProfile->SetBit(TH2F::kCanRebin);
1168 hAllPixelList->Add( hAllPixelProfile[order] );
1169
1170//------------------------------------------------------------------------
1171// SetHistogramAttributes(
1172// "AllPixelPeakProfile",
1173// order,
1174// gProfileAttrib,
1175// MaxAmplOfFirstPulse
1176// );
1177 histoname = "hAllPixelProfile2";
1178 histoname += order;
1179 if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
1180 hAllPixelProfile2[order] = new TProfile(
1181 histoname,
1182 "Mean value of each slice in overlay plot for All Pixels (Tprofile)",
1183 gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,//nbinsx
1184 (-1*gPixelOverlayXaxisLeft)-0.5, //xlow
1185 gPixelOverlayXaxisRight-0.5 , //xup
1186 // 512, //nbinsy
1187 -55.5, //ylow
1188 300.5, //yup
1189 "s"); //option
1190 hAllPixelProfile2[order]->SetAxisRange(
1191 gBSLMean - 5,
1192 (gGainMean*(order+1)) + 10,
1193 "Y");
1194 hAllPixelProfile2[order]->SetLineColor(kRed);
1195 hAllPixelProfile2[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
1196 hAllPixelProfile2[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
1197 //hPixelProfile->SetBit(TH2F::kCanRebin);
1198 hAllPixelList->Add( hAllPixelProfile2[order] );
1199
1200
1201
1202 }//end of for over order
1203 if (verbosityLevel > 2) cout << "...done" << endl;
1204}
1205// end of BookHistos
1206//----------------------------------------------------------------------------
1207
1208
1209void CreateRootFile(const char * loc_fname, int verbosityLevel)
1210{
1211 TFile* tf = new TFile(loc_fname,"UPDATE");
1212 if (tf->IsZombie())
1213 {
1214 cout << "Error opening file" << endl;
1215 exit(-1);
1216 }
1217 else {
1218 if (verbosityLevel > 1) cout << "creating root-file successfull" << endl;
1219 tf->Close();
1220 }
1221}
1222//end of CreateRootFile
1223//----------------------------------------------------------------------------
1224
1225
1226TFile* OpenRootFile(const char * loc_fname, int verbosityLevel)
1227{
1228 TFile* tf = new TFile(loc_fname,"UPDATE");
1229 if (tf->IsZombie())
1230 {
1231 cout << "Error opening file" << endl;
1232 exit(-1);
1233 }
1234 else {
1235 if (verbosityLevel > 1) cout << "...opening root-file:"
1236 << loc_fname
1237 << " successfull" << endl;
1238 }
1239 return tf;
1240}
1241//end of OpenRootFile
1242//----------------------------------------------------------------------------
1243
1244
1245void CloseRootFile(TFile* tf)
1246{
1247 if (tf->IsZombie())
1248 {
1249 cout << "Error closing file" << endl;
1250 exit(-1);
1251 } else {
1252 tf->Close();
1253 }
1254}
1255//end of CloseRootFile
1256//----------------------------------------------------------------------------
1257
1258
1259void SaveHistograms(
1260 const char* loc_fname,
1261 const char* subdirectory,
1262 TObjArray* histList,
1263 int verbosityLevel
1264 )
1265{
1266 if (verbosityLevel > 2) cout << endl
1267 << "...saving pixel histograms to subdirectory: "
1268 << subdirectory << endl ;
1269
1270 TFile* tf = OpenRootFile(loc_fname, verbosityLevel);
1271 if ( subdirectory != "root")
1272 {
1273 tf->cd();
1274 tf->mkdir(subdirectory,subdirectory);
1275 tf->cd(subdirectory);
1276 }
1277 else
1278 {
1279 tf->cd();
1280 }
1281 histList->Write(); // write the major histograms into the top level directory
1282 if (verbosityLevel > 3) tf->ls();
1283 CloseRootFile( tf ); // close the file
1284 if (verbosityLevel > 2) cout << "...done" << endl;
1285}
1286// end of SaveHistograms
1287//----------------------------------------------------------------------------
1288
1289
1290TString CreateSubDirName(int pixel)
1291{
1292 TString path = "Pixel_";
1293 path += pixel;
1294// cout << "path " << path << endl ;
1295 return path;
1296}
1297// end of CreateSubDirName
1298//----------------------------------------------------------------------------
1299
1300TString CreateSubDirName(const char* title)
1301{
1302 TString path = title;
1303 path += "_Pixel";
1304// cout << "path " << path << endl ;
1305 return path;
1306}
1307// end of CreateSubDirName
1308//----------------------------------------------------------------------------
1309
1310void FillHistograms(
1311 vector<Region>* pZXings,
1312 int AmplWindowWidth,
1313 int eventNumber,
1314 int verbosityLevel
1315 )
1316{
1317 if (verbosityLevel > 2) cout << endl << "...filling pulse histograms" ;
1318 bool use_this_peak=false;
1319 int order_of_pulse=0;
1320 vector<Region>::iterator reg;
1321
1322 //Loop over all found zerocrossings in Region
1323 for (reg = pZXings->begin() ; reg < pZXings->end() ; reg++)
1324 {
1325 //skip those who are at the Rim of the ROI
1326 if (reg->maxPos < 12 || (unsigned int) reg->maxPos > gRegionOfInterest-12)
1327 {
1328 if (verbosityLevel > 2) cout << endl << "\t...out of range" << endl;
1329 continue;
1330 }
1331 //domstest->Fill(reg->maxPos);
1332
1333 // Define axis range of Histogram
1334 int Left = reg->maxPos - gPixelOverlayXaxisLeft;
1335 int Right = reg->maxPos + gPixelOverlayXaxisRight;
1336 int EdgeLeft = reg->halfRisingEdgePos - gPixelOverlayXaxisLeft;
1337 int EdgeRight = reg->halfRisingEdgePos + gPixelOverlayXaxisRight;
1338
1339 if (Left < 0) Left =0;
1340 if (Right > (int)Ameas.size() ) Right=Ameas.size() ;
1341 if (EdgeLeft < 0) EdgeLeft =0;
1342 if (EdgeRight > (int)Ameas.size() ) EdgeRight=Ameas.size() ;
1343
1344
1345 hTesthisto->Fill( reg->slopeOfRisingEdge ) ;
1346 hTesthisto2->Fill( eventNumber, reg->slopeOfRisingEdge ) ;
1347
1348 if (verbosityLevel > 2) cout << endl << "\t...choosing Histogram" ;
1349
1350 //determine order of pulse and which histogram shall be filled
1351 if (verbosityLevel > 2) cout << endl << "\t...choosing Histogram" ;
1352 for(int order = 0; order < MAX_PULS_ORDER; order++)
1353 {
1354 if (Ameas[reg->maxPos] >= ((gGainMean*(order+1)) - (AmplWindowWidth/2))
1355 && Ameas[reg->maxPos] < ((gGainMean*(order+1)) + AmplWindowWidth/2))
1356 {
1357 //hOverlayTemp = hPixelOverlay[order];
1358 if (verbosityLevel > 2) cout << "...#" << order ;
1359 order_of_pulse = order;
1360 use_this_peak = true;
1361 break;
1362 }
1363 else if (order >= (MAX_PULS_ORDER - 1)){
1364 use_this_peak = false;
1365 if (verbosityLevel > 2) cout << "...NONE" << endl ;
1366 }
1367 }
1368
1369 //Fill overlay und profile histograms
1370 if (use_this_peak){
1371 if (verbosityLevel > 2) cout << "...filling" ;
1372 for ( int pos = Left; pos < Right; pos++){
1373// if ();
1374 hPixelOverlay[order_of_pulse]->Fill( pos - (reg->maxPos), Ameas[pos]) ;
1375 hPixelProfile[order_of_pulse]->Fill( pos - (reg->maxPos), Ameas[pos]) ;
1376 hAllPixelOverlay[order_of_pulse]->Fill( pos - (reg->maxPos), Ameas[pos]) ;
1377 hAllPixelProfile[order_of_pulse]->Fill( pos - (reg->maxPos), Ameas[pos]) ;
1378 }
1379 for ( int pos = EdgeLeft; pos < EdgeRight; pos++){
1380// cout << endl << "###filling edge histo###" << endl;
1381 hPixelEdgeOverlay[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), Ameas[pos]) ;
1382 hPixelProfile2[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), Ameas[pos]) ;
1383 hAllPixelEdgeOverlay[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), Ameas[pos]) ;
1384 hAllPixelProfile2[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), Ameas[pos]) ;
1385// cout << endl << "######" << endl;
1386 }
1387 if (verbosityLevel > 2) cout << "...done" << endl;
1388 }
1389 //Breakout option
1390 if (breakout)
1391 break;
1392 }
1393 if (verbosityLevel > 2) cout << "...done" << endl;
1394}
1395// end of FillHistograms
1396//----------------------------------------------------------------------------
1397
1398void DrawTestHistograms(
1399 int verbosityLevel
1400 )
1401{
1402 if (verbosityLevel > 2) cout << endl << "...drawing pulse histograms" ;
1403
1404 cgpTestHistos->cd(1);
1405 hTesthisto->Draw();
1406 cgpTestHistos->cd(2);
1407 hTesthisto2->Draw("BOX");
1408}
1409// end of DrawTestHistograms
1410//----------------------------------------------------------------------------
1411
1412void DrawPulseHistograms(
1413 int verbosityLevel,
1414 int max_pulse_order
1415 )
1416{
1417 if (verbosityLevel > 2) cout << endl << "...drawing pulse histograms" ;
1418 for (int pulse_order = 0; pulse_order < max_pulse_order; pulse_order++)
1419 {
1420 cgpPixelPulses[pulse_order]->cd(1);
1421 hPixelOverlay[pulse_order]->Draw("COLZ");
1422 cgpPixelPulses[pulse_order]->cd(3);
1423 hPixelProfile[pulse_order]->Draw("COLZ");
1424 hPixelProfile2[pulse_order]->Draw("SAME");
1425
1426 cgpPixelPulses[pulse_order]->cd(2);
1427 hPixelEdgeOverlay[pulse_order]->Draw("COLZ");
1428// cgpPixelPulses[pulse_order]->cd(4);
1429// hTesthisto->Draw();
1430
1431 cgpAllPixelPulses[pulse_order]->cd(1);
1432 hAllPixelOverlay[pulse_order]->Draw("COLZ");
1433 cgpAllPixelPulses[pulse_order]->cd(2);
1434 hAllPixelEdgeOverlay[pulse_order]->Draw("COLZ");
1435
1436 cgpAllPixelPulses[pulse_order]->cd(3);
1437 hAllPixelProfile[pulse_order]->Draw("COLZ");
1438 hAllPixelProfile2[pulse_order]->Draw("SAME");
1439
1440 }
1441}
1442// end of DrawPulseHistograms
1443//----------------------------------------------------------------------------
1444
1445
1446void DrawMaximumHistograms(
1447 int verbosityLevel,
1448 int max_pulse_order
1449 )
1450{
1451 if (verbosityLevel > 2) cout << endl << "...drawing maximum histograms" ;
1452 for (int pulse_order = 0; pulse_order < max_pulse_order; pulse_order++)
1453 {
1454 cgpPixelPulses[pulse_order]->cd(4);
1455 hPixelMax[pulse_order]->Draw();
1456 hPixelMedian[pulse_order]->Draw("SAME");
1457 hPixelMean[pulse_order]->Draw("SAME");
1458
1459 cgpAllPixelPulses[pulse_order]->cd(4);
1460 hAllPixelMax[pulse_order]->Draw();
1461 hAllPixelMedian[pulse_order]->Draw("SAME");
1462 hAllPixelMean[pulse_order]->Draw("SAME");
1463// cgpAllPixelPulses[pulse_order]->cd(3);
1464 // hAllPixelMaxGaus[pulse_order]->Draw("COLZ");
1465 }
1466}
1467// end of DrawMaximumHistograms
1468//----------------------------------------------------------------------------
1469
1470
1471void UpdateCanvases(
1472 int verbosityLevel,
1473 int max_pulse_order
1474 )
1475{
1476 if (verbosityLevel > 3) cout << endl << "...updating canvases" ;
1477 for (int pulse_order = 0; pulse_order < max_pulse_order; pulse_order++)
1478 {
1479 cgpPixelPulses[pulse_order]->Modified();
1480 cgpPixelPulses[pulse_order]->Update();
1481 cgpAllPixelPulses[pulse_order]->Modified();
1482 cgpAllPixelPulses[pulse_order]->Update();
1483 cgpTestHistos->Modified();
1484 cgpTestHistos->Update();
1485 }
1486}
1487// end of UpdateCanvases
1488//----------------------------------------------------------------------------
1489
1490
1491
1492void
1493PlotMedianEachSliceOfPulse(
1494 unsigned int max_pulse_order,
1495 int fitdata,
1496 int verbosityLevel
1497 )
1498{
1499 if (verbosityLevel > 2) cout << endl
1500 << "...calculating pulse shape of slice's Median"
1501 << endl;
1502 for (unsigned int pulse_order = 0 ; pulse_order < max_pulse_order ; pulse_order ++)
1503 {
1504 if (verbosityLevel > 2) cout << "\t...calculation of "
1505 << "pulse order " << pulse_order;
1506
1507 Int_t nbins = hPixelOverlay[pulse_order]->GetXaxis()->GetNbins();
1508
1509 for (Int_t TimeSlice=1;TimeSlice<=nbins;TimeSlice++) {
1510 TH1 *h1 = hPixelOverlay[pulse_order]->ProjectionY("",TimeSlice,TimeSlice);
1511 float median = MedianOfH1(h1);
1512 float mean = h1->GetMean();
1513 if (verbosityLevel > 2) printf("Median of Slice %d, Median=%g\n",TimeSlice,median);
1514 delete h1;
1515 hPixelMedian[pulse_order]->SetBinContent(TimeSlice, median );
1516 hPixelMean[pulse_order]->SetBinContent(TimeSlice, mean );
1517 hAllPixelMedian[pulse_order]->SetBinContent(TimeSlice, median );
1518 hAllPixelMean[pulse_order]->SetBinContent(TimeSlice, mean );
1519 }
1520
1521 if (verbosityLevel > 2) cout << "\t...done" << endl;
1522
1523 if (fitdata)
1524 {
1525 FitMaxPropabilityPuls(
1526 hPixelMedian[pulse_order],
1527 verbosityLevel
1528 );
1529 }
1530 }
1531}
1532// end of PlotMedianEachSliceOfPulse
1533//----------------------------------------------------------------------------
1534
1535Double_t MedianOfH1 (
1536 TH1* h1
1537 )
1538{
1539 //compute the median for 1-d histogram h1
1540 Int_t nbins = h1->GetXaxis()->GetNbins();
1541 Double_t *x = new Double_t[nbins];
1542 Double_t *y = new Double_t[nbins];
1543 for (Int_t i=0;i<nbins;i++) {
1544 x[i] = h1->GetXaxis()->GetBinCenter(i+1);
1545 y[i] = h1->GetBinContent(i+1);
1546 }
1547 Double_t median = TMath::Median(nbins,x,y);
1548 delete [] x;
1549 delete [] y;
1550 return median;
1551}
1552// end of PlotMedianEachSliceOfPulse
1553//----------------------------------------------------------------------------
1554
1555void
1556PlotPulseShapeOfMaxPropability(
1557 unsigned int max_pulse_order,
1558 int fitdata,
1559 int verbosityLevel
1560 )
1561{
1562 if (verbosityLevel > 2) cout << endl
1563 << "...calculating pulse shape of max. propability"
1564 << endl;
1565 for (unsigned int pulse_order = 0 ; pulse_order < max_pulse_order ; pulse_order ++)
1566 {
1567 if (verbosityLevel > 2) cout << "\t...calculation of "
1568 << "pulse order " << pulse_order;
1569 // vector max_value_of to number of timeslices in Overlay Spectra
1570 float max_value_of_slice;
1571
1572 Int_t nbins = hPixelOverlay[pulse_order]->GetXaxis()->GetNbins();
1573 if (verbosityLevel > 2) cout << "...generating projections" << endl;
1574 for (Int_t TimeSlice=1;TimeSlice<=nbins;TimeSlice++)
1575 {
1576 if (verbosityLevel > 2) cout << "...Timeslice: " << TimeSlice;
1577 //put maximumvalue of every Y-projection of every timeslice into vector
1578 TH1D *h1 = hPixelOverlay[pulse_order]->ProjectionY( "", TimeSlice,TimeSlice);
1579
1580 max_value_of_slice = h1->GetBinCenter( h1->GetMaximumBin() );
1581
1582 if (verbosityLevel > 2) cout << " with max value "
1583 << max_value_of_slice << endl;
1584 hPixelMax[pulse_order]->SetBinContent(TimeSlice, max_value_of_slice );
1585 hAllPixelMax[pulse_order]->SetBinContent(TimeSlice, max_value_of_slice );
1586 delete h1;
1587 }
1588
1589 if (verbosityLevel > 2) cout << "\t...done" << endl;
1590
1591 if (fitdata)
1592 {
1593 FitMaxPropabilityPuls(
1594 hPixelMax[pulse_order],
1595 verbosityLevel
1596 );
1597 }
1598 }
1599}
1600
1601
1602void
1603FitMaxPropabilityPuls(
1604 TH1F* hMaximumTemp,
1605 int verbosityLevel
1606 )
1607 {
1608 if (verbosityLevel > 2) cout << "...fit Landau in histograms" ;
1609 if (verbosityLevel > 2) cout << "\t...calculating Landaufit" ;
1610 hMaximumTemp->Fit("landau", "", "", -50, 250);
1611 if (verbosityLevel > 2) cout << "...done" << endl;
1612 }
1613
1614
1615//void SetHistogrammSettings(
1616// const char* histoType,
1617// int max_puls_order,
1618// int MaxAmplOfFirstPulse,
1619// int verbosityLevel
1620// )
1621//{
1622// if (verbosityLevel > 2) cout << endl << "...set histogram Settings" << endl;
1623// for (int pulse_order = 1; pulse_order <= max_puls_order; pulse_order++)
1624// {
1625// qHistSetting[pulse_order].yMax = (MaxAmplOfFirstPulse * pulse_order);
1626// cout << "\t Max @ " << gHistSetting[pulse_order].yMax << endl;
1627// gHistSetting[pulse_order].name = "h";
1628// gHistSetting[pulse_order].name += histoType;
1629// gHistSetting[pulse_order].name += "PeakOverlay";
1630// gHistSetting[pulse_order].name += pulse_order;
1631// gHistSetting[pulse_order].title = "PeakOverlay with Amplitude around ";
1632// gHistSetting[pulse_order].title += gHistSetting[pulse_order].yMax;
1633// gHistSetting[pulse_order].title += " mV";
1634// gHistSetting[pulse_order].Mname = "h";
1635// gHistSetting[pulse_order].Mname += histoType;
1636// gHistSetting[pulse_order].Mname += "PeakMaximum";
1637// gHistSetting[pulse_order].Mname += pulse_order;
1638// gHistSetting[pulse_order].Mtitle = "Peak (approx) derived with Maximum of above Spektrum with Max of ";
1639// gHistSetting[pulse_order].Mtitle += gHistSetting[pulse_order].yMax;
1640// gHistSetting[pulse_order].Mtitle += " mV";
1641// }
1642// if (verbosityLevel > 2) cout << "...done" << endl;
1643//}
1644
1645
1646//void CreatePulseCanvas(
1647// TCanvas* array_of_canvases,
1648// unsigned int max_pulse_order,
1649// const char* canvas_name,
1650// const char* canvas_title,
1651// const char* pixel,
1652// int verbosityLevel
1653// )
1654//{
1655// if (verbosityLevel > 2) cout << "...building Canvases" ;
1656
1657// for (
1658// unsigned int pulse_order = 1;
1659// pulse_order <= max_pulse_order;
1660// pulse_order++
1661// )
1662// {
1663// TString name = canvas_name;
1664// name += " ";
1665// name += pulse_order;
1666// name += " ";
1667
1668// TString title = "Order ";
1669// title += pulse_order;
1670// title += " ";
1671// title += canvas_title;
1672// title += " of ";
1673// title += pixel;
1674// title += ". Pixel";
1675
1676// array_of_canvases[pulse_order] = new TCanvas(name, title, 1, 1, 1400, 700);
1677// array_of_canvases[pulse_order].Divide(2,2);
1678// }
1679// if (verbosityLevel > 2) cout << "...done" << endl;
1680//}
1681
1682
1683//void
1684//DeletePulseCanvas(
1685// TCanvas* array_of_canvases,
1686// unsigned int max_pulse_order
1687// )
1688//{
1689// for (
1690// unsigned int pulse_order = 1;
1691// pulse_order <= max_pulse_order;
1692// pulse_order++
1693// )
1694// {
1695// delete array_of_canvases[pulse_order];
1696// }
1697//}
1698
1699
1700
1701
1702void
1703DeletePixelCanvases( int verbosityLevel )
1704{
1705 if (verbosityLevel > 2) cout << endl << "...delete pixel Canvas" ;
1706 for (int pulse_order = 0; pulse_order < MAX_PULS_ORDER; pulse_order ++)
1707 {
1708 delete cgpPixelPulses[pulse_order];
1709 }
1710}
1711
1712//----------------------------------------------------------------------------
1713//----------------------------------------------------------------------------
1714
1715//void SetHistogramAttributes(
1716// const char* histoType,
1717// int pulse_order,
1718// histAttrib_t* histoAttrib,
1719// int max_ampl_of_first_pulse
1720// )
1721//{
1722// histoAttrib[pulse_order].yMax = (max_ampl_of_first_pulse * pulse_order);
1723// histoAttrib[pulse_order].yMin = (histoAttrib[pulse_order].yMax/2);
1724// histoAttrib[pulse_order].name = "h";
1725// histoAttrib[pulse_order].name += histoType;
1726// histoAttrib[pulse_order].name = "_";
1727// histoAttrib[pulse_order].name += pulse_order;
1728
1729// histoAttrib[pulse_order].title += histoType;
1730// histoAttrib[pulse_order].title += " with Amplitude around ";
1731// histoAttrib[pulse_order].title += histoAttrib[pulse_order].yMax;
1732// histoAttrib[pulse_order].title += " mV";
1733
1734// histoAttrib[pulse_order].xTitle += "Timeslices [a.u.]";
1735// histoAttrib[pulse_order].yTitle += "Amplitude [mV]";
1736//}
1737
1738void WritePixelTemplateToCsv(
1739 TString path,
1740 const char* csv_file_name,
1741 const char* overlay_method,
1742 int pixel,
1743 int verbosityLevel
1744 )
1745{
1746 path += csv_file_name;
1747 path += "_";
1748 path += pixel;
1749 path += ".csv";
1750
1751 Int_t nbins = hPixelMax[0]->GetXaxis()->GetNbins();
1752 if (verbosityLevel > 0)
1753 {
1754 cout << "writing point-set to csv file: " ;
1755 cout << path << endl;
1756 cout << "...opening file" << endl;
1757 }
1758 if (verbosityLevel > 2) cout << "...number of bins " << nbins << endl;
1759 ofstream out;
1760 out.open( path );
1761 out << "### point-set of a single photon pulse template" << endl;
1762 out << "### template determined with pulse overlay at: "
1763 << overlay_method << endl;
1764 out << "### Slice's Amplitude determined by calculating the " << endl
1765 << "### value of maximum propability of slice -> Amplitude1 " << endl
1766 << "### mean of slice -> Amplitude2 " << endl
1767 << "### median of slice -> Amplitude3 " << endl
1768 << "### for each slice" << endl;
1769 out << "### Pixel number (CHid): " << pixel << endl
1770 << endl;
1771
1772 out << "time [slices],Amplitude1 [mV],Amplitude2 [mV],Amplitude3 [mV]" << endl;
1773
1774 for (int TimeSlice=1;TimeSlice<=nbins;TimeSlice++)
1775 {
1776 out << TimeSlice << "," ;
1777 out << hPixelMax[0]->GetBinContent(TimeSlice) << ",";
1778 out << hPixelMean[0]->GetBinContent(TimeSlice) << ",";
1779 out << hPixelMedian[0]->GetBinContent(TimeSlice) << endl;
1780 }
1781 out.close();
1782 if (verbosityLevel > 0) cout << "...file closed" << endl;
1783}
1784
1785void WriteAllPixelTemplateToCsv(
1786 TString path,
1787 const char* csv_file_name,
1788 const char* overlay_method,
1789 int verbosityLevel
1790 )
1791{
1792 path += csv_file_name;
1793 path += "_AllPixel.csv";
1794
1795 Int_t nbins = hAllPixelMax[0]->GetXaxis()->GetNbins();
1796 if (verbosityLevel > 0)
1797 {
1798 cout << "writing point-set to csv file: " ;
1799 cout << path << endl;
1800 cout << "...opening file" << endl;
1801 }
1802 if (verbosityLevel > 2) cout << "...number of bins " << nbins << endl;
1803 ofstream out;
1804 out.open( path );
1805 out << "### point-set of a single photon pulse template" << endl;
1806 out << "### template determined with pulse overlay at: "
1807 << overlay_method << "of all Pixels";
1808 out << "### Slice's Amplitude determined by calculating the " << endl
1809 << "### value of maximum propability of slice -> Amplitude1 " << endl
1810 << "### mean of slice -> Amplitude2 " << endl
1811 << "### median of slice -> Amplitude3 " << endl
1812 << "### for each slice" << endl
1813 << endl;
1814
1815 out << "time [slices],Amplitude1 [mV],Amplitude2 [mV],Amplitude3 [mV]" << endl;
1816
1817 for (int TimeSlice=1;TimeSlice<=nbins;TimeSlice++)
1818 {
1819 out << TimeSlice << "," ;
1820 out << hAllPixelMax[0]->GetBinContent(TimeSlice) << ",";
1821 out << hAllPixelMean[0]->GetBinContent(TimeSlice) << ",";
1822 out << hAllPixelMedian[0]->GetBinContent(TimeSlice) << endl;
1823 }
1824 out.close();
1825 if (verbosityLevel > 0) cout << "...file closed" << endl;
1826}
Note: See TracBrowser for help on using the repository browser.