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

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