source: fact/tools/rootmacros/PulseTemplates/FCalcPulseTemplate.C@ 13636

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