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

Last change on this file since 12919 was 12919, checked in by Jens Buss, 14 years ago
Bug Fixes
File size: 32.9 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
33#include <stdio.h>
34#include <stdint.h>
35#include <cstdio>
36
37#define NPIX 1440
38#define NCELL 1024
39#define FAD_MAX_SAMPLES 1024
40
41#define HAVE_ZLIB
42
43//----------------------------------------------------------------------------
44// rootmacros
45//----------------------------------------------------------------------------
46
47#include "fits.h"
48
49#include "openFits.h"
50#include "openFits.c"
51
52#include "discriminator.h"
53#include "discriminator.C"
54#include "zerosearch.h"
55#include "zerosearch.C"
56#include "factfir.C"
57
58#include "DrsCalibration.C"
59#include "DrsCalibration.h"
60
61#include "SpikeRemoval.h"
62#include "SpikeRemoval.C"
63
64//----------------------------------------------------------------------------
65// Decleration of global variables
66//----------------------------------------------------------------------------
67
68bool breakout=false;
69bool UseThisPeak=false;
70int NEvents;
71vector<int16_t> AllPixelDataVector;
72vector<int16_t> StartCellVector;
73unsigned int CurrentEventID;
74size_t PXLxROI;
75UInt_t RegionOfInterest;
76UInt_t NumberOfPixels;
77TString histotitle;
78
79size_t TriggerOffsetROI, RC;
80size_t drs_n;
81vector<float> drs_basemean;
82vector<float> drs_gainmean;
83vector<float> drs_triggeroffsetmean;
84
85vector<float> Ameas(FAD_MAX_SAMPLES); // copy of the data (measured amplitude
86vector<float> Vcorr(FAD_MAX_SAMPLES); // corrected Values
87vector<float> Vslide(FAD_MAX_SAMPLES); // sliding average result
88vector<float> Vcfd(FAD_MAX_SAMPLES); // CDF result
89vector<float> Vcfd2(FAD_MAX_SAMPLES); // CDF result + 2nd sliding average
90
91typedef struct{
92 float maxAmpl;
93 float countOfMax;
94 } OverlayMaximum;
95
96typedef struct{
97 TString name;
98 TString title;
99 TString Mname;
100 TString Mtitle;
101 float yMax;
102} hist_t;
103
104hist_t histSetting[4];
105
106vector<OverlayMaximum> SingleMaximum;
107// histograms
108const int NumberOfDebugHistoTypes = 7;
109
110const unsigned int
111 Ameas_ = 0,
112 N1mean_ = 1,
113 Vcorr_ = 2,
114 Vtest_ = 3,
115 Vslide_ = 4,
116 Vcfd_ = 5,
117 Vcfd2_ = 6;
118
119//----------------------------------------------------------------------------
120// Initialisation of histograms
121//----------------------------------------------------------------------------
122
123TH1F* h= NULL;
124TH1F *debugHistos= NULL;
125TH2F *hOverlayTemp = NULL;
126TH2F *hPeakOverlay[4];//histogrammm for overlay of detected Peaks
127TH2F *hSinglePeakOverlay;//histogrammm for overlay of detected Peaks
128TH2F *hDoublePeakOverlay;
129TH2F *hTripplePeakOverlay;
130TH2F *hLargePeakOverlay;
131
132TH1F *hMaximumTemp = NULL;
133TH1F *hPeakMaximum[4];
134TH1F *hSinglePeakMaximum;
135TH1F *hDoublePeakMaximum;
136TH1F *hTripplePeakMaximum;
137TH1F *hLargePeakMaximum;
138TH1D *hProjPeak = NULL;
139int gPeakOverlayXaxisLeft;
140int hPeakOverlayXaxisRight;
141
142TObjArray hList;
143TObjArray hListBaseline;
144
145//----------------------------------------------------------------------------
146// Functions
147//----------------------------------------------------------------------------
148
149void BookHistos( );
150void SaveHistograms( const char * );
151
152//----------------------------------------------------------------------------
153//----------------------------------------------------------------------------
154// MAIN - Funtion
155//----------------------------------------------------------------------------
156//----------------------------------------------------------------------------
157int FPulseTemplate(
158 char* datafilename = "../data/2011/11/10/20111110_005.fits.gz",
159 const char* drsfilename = "../data/2011/11/10/20111110_003.drs.fits.gz",
160 const char* OutRootFileName = "../analysis/20111110_003.test.root",
161 bool ProduceGraphic = true,
162 bool spikeDebug = false,
163 bool fitdata = false,
164 int verbosityLevel = 1, // different verbosity levels can be implemented here
165 int firstevent = 0,
166 int nevents = -1,
167 int firstpixel = 400,
168 int npixel = 1,
169 int PintervalWidth = 5,
170 int PintervalMin = 8,
171 int PintervalMax = 32, //4 Histogramms will be drawn, decide how far pulsehights differ from eachother
172 int avg1 = 14,
173 int avg2 = 8,
174 int OverlayWindowLeft = 70,
175 int OverlayWindowRight = 230)
176{
177//-----------------------------------------------------------------------------
178// histogramm Settings
179//-----------------------------------------------------------------------------
180
181 gPeakOverlayXaxisLeft = OverlayWindowLeft;
182 hPeakOverlayXaxisRight = OverlayWindowRight;
183
184 gStyle->SetPalette(1,0);
185 gROOT->SetStyle("Plain");
186
187// Create (pointer to) Canvases, which are used in every run,
188// also in 'non-debug' runs
189 // Canvases only need if spike Debug, but I want to deklare
190 // the pointers anyway ...
191 TCanvas *cFiltered = NULL;
192 if (spikeDebug){
193 cFiltered = new TCanvas("cFiltered","filtered DRS Waveforms", 1,310,400,300);
194 cFiltered->Divide(1, 3);
195 }
196 // Canvases to show the peak template
197 TCanvas *cPulsetypes = NULL;
198 cPulsetypes = new TCanvas("cPulsetypes", "Overlay of detected Peaks", 1, 1, 1400, 700);
199 cPulsetypes->Divide(4,2);
200
201
202 for (int stepper = 0; stepper < 4; stepper++){
203 if (stepper == 0) histSetting[stepper].yMax = PintervalMin;
204 else if (stepper == 3) histSetting[stepper].yMax = PintervalMax;
205 else histSetting[stepper].yMax = ((PintervalMax-PintervalMin)*(stepper)/3)+PintervalMin;
206 cout << "Max @ " << histSetting[stepper].yMax << endl;
207 histSetting[stepper].name = "hPeakOverlay";
208 histSetting[stepper].name += stepper;
209 histSetting[stepper].title = "PeakOverlay with Amplitude around ";
210 histSetting[stepper].title += histSetting[stepper].yMax;
211 histSetting[stepper].title += " mV";
212 histSetting[stepper].Mname = "hPeakMaximum";
213 histSetting[stepper].Mname += stepper;
214 histSetting[stepper].Mtitle = "Peak (approx) derived with Maximum of above Spektrum with Max of ";
215 histSetting[stepper].Mtitle += histSetting[stepper].yMax;
216 histSetting[stepper].Mtitle += " mV";
217 }
218
219//-----------------------------------------------------------------------------
220// Filter-Settings
221//-----------------------------------------------------------------------------
222// CFD filter settings
223 int k_cfd = 10;
224 vector<double> a_cfd(k_cfd, 0);
225 double b_cfd = 1.;
226 a_cfd[0]=-0.75;
227 a_cfd[k_cfd-1]=1.;
228
229//-----------------------------------------------------------------------------
230// prepare datafile
231//-----------------------------------------------------------------------------
232
233// Open the data file
234
235 fits * datafile;
236 // Opens the raw data file and 'binds' the variables given as
237 // Parameters to the data file. So they are filled with
238 // raw data as soon as datafile->GetRow(int) is called.
239 NEvents = openDataFits( datafilename, &datafile,
240 AllPixelDataVector, StartCellVector, CurrentEventID,
241 RegionOfInterest, NumberOfPixels, PXLxROI, verbosityLevel);
242 if (NEvents == 0){
243 cout << "return code of OpenDataFile:" << datafilename<< endl;
244 cout << "is zero -> aborting." << endl;
245 return 1;
246 }
247
248 if (verbosityLevel > 0)
249 cout <<"number of events in file: "<< NEvents << "\t";
250 if ( nevents == -1 || nevents > NEvents ) nevents = NEvents; // -1 means all!
251 if (verbosityLevel > 0)
252 cout <<"of, which "<< nevents << "will be processed"<< endl;
253
254 if (verbosityLevel > 0)
255 cout <<"Total # of Pixel: "<< NumberOfPixels << "\t";
256 if ( npixel == -1 || npixel > (int)NumberOfPixels ) npixel = (int)NumberOfPixels; // -1 means all!
257 if (verbosityLevel > 0)
258 cout <<"of, which "<< npixel << "will be processed"<< endl;
259//Get the DRS calibration
260 RC = openCalibFits( drsfilename, drs_basemean, drs_gainmean, drs_triggeroffsetmean, TriggerOffsetROI);
261 if (RC == 0){
262 cout << "return code of openCalibFits:" << drsfilename << endl;
263 cout << "is zero -> aborting." << endl;
264 return 1;
265 }
266// Book the histograms
267 BookHistos( );
268
269//-----------------------------------------------------------------------------
270// Loops over Every Event and Pixel
271//-----------------------------------------------------------------------------
272 for ( int ev = firstevent; ev < firstevent + nevents; ev++) {
273 // Get an Event --> consists of 1440 Pixel ...erm....data
274 datafile->GetRow( ev );
275
276//-------------------------------------
277// Loop over every Pixel of Event
278//-------------------------------------
279
280 for ( int pix = firstpixel; pix < firstpixel + npixel; pix++ ){
281 if (verbosityLevel > 0){
282 cout << endl;
283 if (pix == firstpixel){
284 cout << "...processing Event: " << CurrentEventID << "/" << nevents << endl;
285 }
286 }
287
288//-------------------------------------
289// Apply Filters
290//-------------------------------------
291 cout << "...applying DrsCalibration";
292 applyDrsCalibration( Ameas,pix,12,12,
293 drs_basemean, drs_gainmean, drs_triggeroffsetmean,
294 RegionOfInterest, AllPixelDataVector, StartCellVector);
295 cout << "...done " << endl;
296
297 // finds spikes in the raw data, and interpolates the value
298 // spikes are: 1 or 2 slice wide, positive non physical artifacts
299 cout << "...removeing Spikes";
300 removeSpikes (Ameas, Vcorr);
301 cout << "...done " << endl;
302 // filter Vcorr with sliding average using FIR filter function
303 cout << "...applying sliding average filter";
304 sliding_avg(Vcorr, Vslide, avg1);
305 cout << "...done " << endl;
306 // filter Vslide with CFD using FIR filter function
307 cout << "...apllying factfir filter";
308 factfir(b_cfd , a_cfd, k_cfd, Vslide, Vcfd);
309 cout << "...done " << endl;
310 // filter Vcfd with sliding average using FIR filter function
311 cout << "...applying 2nd sliding average filter";
312 sliding_avg(Vcfd, Vcfd2, avg2);
313 cout << "...done " << endl;
314
315//-------------------------------------
316// Search vor Zero crossings
317//-------------------------------------
318 // peaks in Ameas[] are found by searching for zero crossings
319 // in Vcfd2
320 // first Argument 1 means ... *rising* edge
321 // second Argument 1 means ... search with stepsize 1 ... 10 is okay as well
322 vector<Region> * zXings = zerosearch( Vcfd2 , 1 , 8);
323 // zXings means "zero cross ings"
324 EnlargeRegion(*zXings, 10, 10);
325 findAbsMaxInRegions(*zXings, Vslide);
326 removeMaximaBelow( *zXings, 3.0);
327 removeRegionWithMaxOnEdge( *zXings, 2);
328 removeRegionOnFallingEdge( *zXings, 100);
329
330//-----------------------------------------------------------------------------
331// Fill Overlay Histos
332//-----------------------------------------------------------------------------
333 vector<Region>::iterator it;
334 for (it = zXings->begin() ; it < zXings->end() ; it++){
335 if (it->maxPos < 12 || it->maxPos > RegionOfInterest-12)
336 continue;
337 //domstest->Fill(it->maxPos);
338 int Left = it->maxPos - OverlayWindowLeft;
339 int Right = it->maxPos + OverlayWindowRight;
340 if (Left < 0)
341 Left =0;
342 if (Right > (int)Vcorr.size() )
343 Right=Vcorr.size() ;
344
345 cout << "...choosing Histogram" ;
346 for(int stepper = 0; stepper < 4; stepper++){
347 if (Vslide[it->maxPos] >= (histSetting[stepper].yMax - (PintervalWidth/2)) && Vslide[it->maxPos] < (histSetting[stepper].yMax + PintervalWidth/2)){
348 hOverlayTemp = hPeakOverlay[stepper];
349 cout << "...#" << stepper ;
350 UseThisPeak = true;
351 break;
352 }
353 else if (stepper > 2){
354 UseThisPeak = false;
355 cout << "...NONE" << endl ;
356 }
357 }
358
359// if (Vslide[it->maxPos] >= 5 && Vslide[it->maxPos] < 13) hOverlayTemp = &*hSinglePeakOverlay;
360// else if (Vslide[it->maxPos] >= 13 && Vslide[it->maxPos] < 23) hOverlayTemp = hDoublePeakOverlay;
361// else if (Vslide[it->maxPos] >= 23 && Vslide[it->maxPos] < 33) hOverlayTemp = hTripplePeakOverlay;
362// else if (Vslide[it->maxPos] >= 33) hOverlayTemp = hLargePeakOverlay;
363 if (UseThisPeak){
364 cout << "...filling" ;
365 for ( int pos = Left; pos < Right; pos++){
366 // if (Vslide[it->maxPos] >= 5 && Vslide[it->maxPos] < 15) hSinglePeakOverlay->Fill( pos - (it->maxPos), Vslide[pos]);
367 // else if (Vslide[it->maxPos] >= 15 && Vslide[it->maxPos] < 25) hDoublePeakOverlay->Fill( pos - (it->maxPos), Vslide[pos]);
368 // else if (Vslide[it->maxPos] >= 25 && Vslide[it->maxPos] < 35) hTripplePeakOverlay->Fill( pos - (it->maxPos), Vslide[pos]);
369 // else if (Vslide[it->maxPos] >= 35) hOverlayTemp = hLargePeakOverlay;
370
371 hOverlayTemp->Fill( pos - (it->maxPos), Vslide[pos]) ;
372 }
373 cout << "...done" << endl;
374 }
375
376 }
377
378//-----------------------------------------------------------------------------
379// Spike Debug
380//-----------------------------------------------------------------------------
381 if ( spikeDebug ){
382 // TODO do this correct. The vectors should be the rigt ones... this is just luck
383 debugHistos[Ameas_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
384 debugHistos[Vcorr_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
385 debugHistos[Vslide_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
386 debugHistos[Vcfd_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
387 debugHistos[Vcfd2_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
388
389 for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
390 debugHistos[Ameas_].SetBinContent(sl, Ameas[sl]);
391 debugHistos[Vcorr_].SetBinContent(sl, Vcorr[sl]);
392 debugHistos[Vslide_].SetBinContent( sl, Vslide[sl] );
393 debugHistos[Vcfd_].SetBinContent( sl, Vcfd[sl] );
394 debugHistos[Vcfd2_].SetBinContent( sl, Vcfd2[sl] );
395 }
396
397
398 cFiltered->cd(1);
399 gPad->SetGrid();
400 debugHistos[Vcorr_].Draw();
401
402 cFiltered->cd(2);
403 gPad->SetGrid();
404 debugHistos[Vslide_].Draw();
405
406 TBox *OneBox;
407 vector<TBox*> MyBoxes;
408 for (unsigned int i=0; i<zXings->size(); i++){
409 OneBox = new TBox(
410 zXings->at(i).maxPos -10 ,
411 zXings->at(i).maxVal -0.5,
412 zXings->at(i).maxPos +10 ,
413 zXings->at(i).maxVal +0.5);
414 OneBox->SetLineColor(kBlue);
415 OneBox->SetLineWidth(1);
416 OneBox->SetFillStyle(0);
417 OneBox->SetFillColor(kRed);
418 MyBoxes.push_back(OneBox);
419 OneBox->Draw();
420 }
421
422 cFiltered->cd(3);
423 gPad->SetGrid();
424 debugHistos[Vcfd2_].Draw();
425 TLine *zeroline = new TLine(0, 0, 1024, 0);
426 zeroline->SetLineColor(kBlue);
427 zeroline->Draw();
428
429 cFiltered->Update();
430
431 //OVERLAY PEAKS
432 for(int stepper = 0; stepper < 4; stepper++){
433 cPulsetypes->cd(stepper+1);
434 gPad->SetGrid();
435 hPeakOverlay[stepper]->Draw("COLZ");
436 }
437
438
439// cPulsetypes->cd(1);
440// gPad->SetGrid();
441// hSinglePeakOverlay->Draw("COLZ");
442//
443// cPulsetypes->cd(2);
444// gPad->SetGrid();
445// hDoublePeakOverlay->Draw("COLZ");
446//
447// cPulsetypes->cd(3);
448// gPad->SetGrid();
449// hTripplePeakOverlay->Draw("COLZ");
450//
451// cPulsetypes->cd(4);
452// gPad->SetGrid();
453// hLargePeakOverlay->Draw("COLZ");
454
455 cPulsetypes->Modified();
456 cPulsetypes->Update();
457
458 //Process gui events asynchronously during input
459 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
460 timer.TurnOn();
461 TString input = Getline("Type 'q' to exit, <return> to go on: ");
462 timer.TurnOff();
463 if (input=="q\n") {
464 breakout=true;
465 }
466
467 //TODO!!!!!!!!!
468 // do some Garbage collection ...
469 // all the Objects on the heap should be deleted here.
470
471 }// end of if(spikeDebug)
472
473 delete zXings;
474 if (breakout)
475 break;
476
477 }
478 //-------------------------------------
479 // end of loop over pixels
480 //-------------------------------------
481 if (ProduceGraphic){
482 if (ev % 50 == 0){
483
484 //OVERLAY PEAKS
485 cout << "...drawing histograms after Loop over Pixels" ;
486 for(int stepper = 0; stepper < 4; stepper++){
487 cPulsetypes->cd(stepper+1);
488 gPad->SetGrid();
489 hPeakOverlay[stepper]->Draw("COLZ");
490 cout << "..." << stepper << "/3" ;
491 }
492// cPulsetypes->cd(1);
493// gPad->SetGrid();
494// hSinglePeakOverlay->Draw("COLZ");
495//
496// cPulsetypes->cd(2);
497// gPad->SetGrid();
498// hDoublePeakOverlay->Draw("hLargePeakOverlayCOLZ");
499//
500// cPulsetypes->cd(3);
501// gPad->SetGrid();
502// hTripplePeakOverlay->Draw("COLZ");
503//
504// cPulsetypes->cd(4);
505// gPad->SetGrid();
506// hLargePeakOverlay->Draw("COLZ");
507 // cPulsetypes->cd(5);
508 // hSinglePeakMaximum->Draw("");
509 cout << "...done" << endl;
510 cPulsetypes->Modified();
511 cPulsetypes->Update();
512 }
513 }
514 if (breakout)
515 break;
516 }
517
518//-------------------------------------
519// end of loop over Events
520//-------------------------------------
521
522//-------------------------------------
523// Histogramm of Maximas in Overlay Spectra
524//-------------------------------------
525 cout << "...producing Maximum peaks:" << endl;
526 for (unsigned int cnt = 0 ; cnt < 4 ; cnt ++){
527 cout << "\t ...peak type " << cnt;
528 hOverlayTemp = hPeakOverlay[cnt];
529 hMaximumTemp = hPeakMaximum[cnt];
530// if (cnt == 1){
531// hMaximumTemp = hSinglePeakMaximum;
532// hOverlayTemp = hSinglePeakOverlay;
533// }
534// else if (cnt == 2){
535// hMaximumTemp = hDoublePeakMaximum;
536// hOverlayTemp = hDoublePeakOverlay;
537// }
538// else if (cnt == 3){
539// hMaximumTemp = hTripplePeakMaximum;
540// hOverlayTemp = hTripplePeakOverlay;
541// }
542// else if (cnt == 4){
543// hMaximumTemp = hLargePeakMaximum;
544// hOverlayTemp = hLargePeakOverlay;
545// }
546 // resize vector SingleMaximum to number of timeslices in Overlay Spectra
547 SingleMaximum.resize(hOverlayTemp->GetNbinsX());
548 cout << "\t ...peak type " << cnt;
549 //put maximumvalue of every Y-projection of every timeslice into vector
550 for (int TimeSlice = 0; TimeSlice <= hOverlayTemp->GetNbinsX(); TimeSlice++ ){
551 histotitle = "hproj_py";
552 histotitle += cnt ;
553 hProjPeak = hOverlayTemp->ProjectionY(histotitle, TimeSlice, TimeSlice, "");
554 SingleMaximum[ TimeSlice ].maxAmpl = (hProjPeak->GetMaximumBin() * 0.5) - 3.5;
555 SingleMaximum[ TimeSlice ].countOfMax = hProjPeak->GetBinContent( hProjPeak->GetMaximumBin() );
556 hMaximumTemp->SetBinContent(TimeSlice, SingleMaximum[ TimeSlice ].maxAmpl );
557 }
558 cout << "...done" << endl;
559 if (fitdata){
560 cout << "...calculating Landaufit" ;
561 hMaximumTemp->Fit("landau", "", "", -50, 250);
562 cout << "...done" << endl;
563 }
564 }
565//-------------------------------------
566// Draw Histograms
567//-------------------------------------
568 if (ProduceGraphic){
569 //OVERLAY PEAKS
570 cout << "...drawing overlay histograms" ;
571 for(int stepper = 0; stepper < 4; stepper++){
572 cPulsetypes->cd(stepper+1);
573 gPad->SetGrid();
574 hPeakOverlay[stepper]->ResetStats();
575 hPeakOverlay[stepper]->Draw("COLZ");
576 }
577 cout << "...done" << endl;
578 cout << "...drawing peak-shape histograms" ;
579 for(int stepper = 0; stepper < 4; stepper++){
580 cPulsetypes->cd(5+stepper);
581 gPad->SetGrid();
582 hPeakMaximum[stepper]->ResetStats();
583 hPeakMaximum[stepper]->Draw();
584 }
585 cout << "...done" << endl;
586
587//cout << "producing...So" << endl;
588// cPulsetypes->cd(1);
589// gPad->SetGrid();
590// hSinglePeakOverlay->ResetStats();
591// hSinglePeakOverlay->Draw("COLZ");
592//cout << "producing...DO" << endl;
593// cPulsetypes->cd(2);
594// gPad->SetGrid();
595// hDoublePeakOverlay->ResetStats();
596// hDoublePeakOverlay->Draw("COLZ");
597//cout << "producing...TO" << endl;
598// cPulsetypes->cd(3);
599// gPad->SetGrid();
600// hTripplePeakOverlay->ResetStats();
601// hTripplePeakOverlay->Draw("COLZ");
602//cout << "producing...LO" << endl;
603// cPulsetypes->cd(4);
604// gPad->SetGrid();
605// hLargePeakOverlay->ResetStats();
606// hLargePeakOverlay->Draw("COLZ");
607//cout << "producing...SM" << endl;
608// cPulsetypes->cd(5);
609// hSinglePeakMaximum->ResetStats();
610// hSinglePeakMaximum->Draw("");
611//cout << "producing...DM" << endl;
612// cPulsetypes->cd(6);
613// hDoublePeakMaximum->ResetStats();
614// hDoublePeakMaximum->Draw("");
615//cout << "producing...TM" << endl;
616// cPulsetypes->cd(7);
617// hTripplePeakMaximum->ResetStats();
618// hTripplePeakMaximum->Draw("");
619//cout << "producing...LM" << endl;
620// cPulsetypes->cd(8);
621// hLargePeakMaximum->ResetStats();
622// hLargePeakMaximum->Draw("");
623//cout << "producing...Done" << endl;
624
625
626
627
628 cPulsetypes->Modified();
629 cPulsetypes->Update();
630
631// cPixelPeakOverlay->cd();
632// hPixelPeakOverlay->Draw("COLZ");
633// cPixelPeakOverlay->Modified();
634// cPixelPeakOverlay->Update();
635
636 }
637 SaveHistograms( OutRootFileName );
638
639 return( 0 );
640}
641
642// booking and parameter settings for all histos
643void BookHistos( ){
644
645 debugHistos = new TH1F[ NumberOfDebugHistoTypes ];
646 for ( int type = 0; type < NumberOfDebugHistoTypes; type++){
647 debugHistos[ type ].SetBins(1024, 0, 1024);
648 debugHistos[ type ].SetLineColor(1);
649 debugHistos[ type ].SetLineWidth(2);
650
651 // set X axis paras
652 debugHistos[ type ].GetXaxis()->SetLabelSize(0.1);
653 debugHistos[ type ].GetXaxis()->SetTitleSize(0.1);
654 debugHistos[ type ].GetXaxis()->SetTitleOffset(1.2);
655 debugHistos[ type ].GetXaxis()->SetTitle(Form("Time slice (%.1f ns/slice)", 1./2.));
656
657 // set Y axis paras
658 debugHistos[ type ].GetYaxis()->SetLabelSize(0.1);
659 debugHistos[ type ].GetYaxis()->SetTitleSize(0.1);
660 debugHistos[ type ].GetYaxis()->SetTitleOffset(0.3);
661 debugHistos[ type ].GetYaxis()->SetTitle("Amplitude (a.u.)");
662 }
663/*
664 ProjHistos = new TH1D[ NumberOfTimeSlices ];
665 for ( int type = 0; type < NumberOfTimeSlices; type++){
666 ProjHistos[ type ].SetBins(1024, 0, 1024);
667 ProjHistos[ type ].SetLineColor(1);
668 ProjHistos[ type ].SetLineWidth(2);
669
670 // set X axis paras
671 ProjHistos[ type ].GetXaxis()->SetLabelSize(0.1);
672 ProjHistos[ type ].GetXaxis()->SetTitleSize(0.1);
673 ProjHistos[ type ].GetXaxis()->SetTitleOffset(1.2);
674 ProjHistos[ type ].GetXaxis()->SetTitle(Form("Time slice (%.1f ns/slice)", 1./2.));
675
676 // set Y axis paras
677 ProjHistos[ type ].GetYaxis()->SetLabelSize(0.1);
678 ProjHistos[ type ].GetYaxis()->SetTitleSize(0.1);
679 ProjHistos[ type ].GetYaxis()->SetTitleOffset(0.3);
680 ProjHistos[ type ].GetYaxis()->SetTitle("Amplitude (a.u.)");
681 }
682*/
683
684//-----------------------------------------------------------------------------
685
686 for (int stepper = 0; stepper < 4; stepper ++){
687 hPeakOverlay[stepper] = new TH2F(histSetting[stepper].name, histSetting[stepper].title,
688 gPeakOverlayXaxisLeft + hPeakOverlayXaxisRight , (-1*gPeakOverlayXaxisLeft)-0.5, hPeakOverlayXaxisRight-0.5 ,
689 512, -55.5, 300.5);
690 hPeakOverlay[stepper]->SetAxisRange(-5.5, histSetting[stepper].yMax+ 25.5, "Y");
691 hPeakOverlay[stepper]->GetXaxis()->SetTitle( "Timeslices" );
692 hPeakOverlay[stepper]->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
693 //hSinglePeakOverlay->SetBit(TH2F::kCanRebin);
694 hList.Add( hPeakOverlay[stepper] );
695
696 hPeakMaximum[stepper] = new TH1F(histSetting[stepper].Mname, histSetting[stepper].Mtitle,
697 gPeakOverlayXaxisLeft + hPeakOverlayXaxisRight ,
698 (-1*gPeakOverlayXaxisLeft)-0.5,
699 hPeakOverlayXaxisRight-0.5
700 );
701 hPeakMaximum[stepper]->SetAxisRange(-0.5, histSetting[stepper].yMax + 15.5, "Y");
702 hPeakMaximum[stepper]->GetXaxis()->SetTitle( "Timeslices" );
703 hPeakMaximum[stepper]->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
704 hList.Add( hPeakMaximum[stepper] );
705 }
706//-----------------------------------------------------------------------------
707
708
709
710//-----------------------------------------------------------------------------
711 hSinglePeakOverlay = new TH2F("hSinglePeakOverlay", "Overlay of detected Single Photon Peaks",
712 gPeakOverlayXaxisLeft + hPeakOverlayXaxisRight , (-1*gPeakOverlayXaxisLeft)-0.5, hPeakOverlayXaxisRight-0.5 ,
713 512, -55.5, 300.5);
714 hSinglePeakOverlay->SetAxisRange(-5.5, 35.5, "Y");
715 hSinglePeakOverlay->GetXaxis()->SetTitle( "Timeslices" );
716 hSinglePeakOverlay->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
717 //hSinglePeakOverlay->SetBit(TH2F::kCanRebin);
718 hList.Add( hSinglePeakOverlay );
719
720 hSinglePeakMaximum = new TH1F("hSinglePeakMaximum", "Single Peak derived with Maximum of above Spektrum",
721 gPeakOverlayXaxisLeft + hPeakOverlayXaxisRight ,
722 (-1*gPeakOverlayXaxisLeft)-0.5,
723 hPeakOverlayXaxisRight-0.5
724 );
725 hSinglePeakMaximum->SetAxisRange(-0.5, 25.5, "Y");
726 hSinglePeakMaximum->GetXaxis()->SetTitle( "Timeslices" );
727 hSinglePeakMaximum->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
728 hList.Add( hSinglePeakMaximum );
729//-----------------------------------------------------------------------------
730
731 hDoublePeakOverlay = new TH2F("hDoublePeakOverlay", "Overlay of detected Double Photon Peaks",
732 gPeakOverlayXaxisLeft + hPeakOverlayXaxisRight , (-1*gPeakOverlayXaxisLeft)-0.5, hPeakOverlayXaxisRight-0.5 ,
733 512, -55.5, 300.5);
734 hDoublePeakOverlay->SetAxisRange(-5.5, 35.5, "Y");
735 hDoublePeakOverlay->GetXaxis()->SetTitle( "Timeslices" );
736 hDoublePeakOverlay->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
737 //hSinglePeakOverlay->SetBit(TH2F::kCanRebin);
738 hList.Add( hDoublePeakOverlay );;
739
740 hDoublePeakMaximum = new TH1F("hSinglePeakMaximum", "Double Peak derived with Maximum of above Spektrum",
741 gPeakOverlayXaxisLeft + hPeakOverlayXaxisRight ,
742 (-1*gPeakOverlayXaxisLeft)-0.5,
743 hPeakOverlayXaxisRight-0.5
744 );
745 hDoublePeakMaximum->SetAxisRange(-0.5, 25.5, "Y");
746 hDoublePeakMaximum->GetXaxis()->SetTitle( "Timeslices" );
747 hDoublePeakMaximum->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
748 hList.Add( hDoublePeakMaximum );
749 //-----------------------------------------------------------------------------
750
751 hTripplePeakOverlay= new TH2F("hTripplePeakOverlay", "Overlay of detected Tripple Photon Peaks",
752 gPeakOverlayXaxisLeft + hPeakOverlayXaxisRight , (-1*gPeakOverlayXaxisLeft)-0.5, hPeakOverlayXaxisRight-0.5 ,
753 512, -55.5, 300.5);
754 hTripplePeakOverlay->SetAxisRange(-5.5, 35.5, "Y");
755 hTripplePeakOverlay->GetXaxis()->SetTitle( "Timeslices" );
756 hTripplePeakOverlay->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
757 //hTripplePeakOverlay->SetBit(TH2F::kCanRebin);
758 hList.Add( hTripplePeakOverlay );;
759
760 hTripplePeakMaximum = new TH1F("hSinglePeakMaximum", "Triple Peak derived with Maximum of above Spektrum",
761 gPeakOverlayXaxisLeft + hPeakOverlayXaxisRight ,
762 (-1*gPeakOverlayXaxisLeft)-0.5,
763 hPeakOverlayXaxisRight-0.5
764 );
765 hTripplePeakMaximum->SetAxisRange(-0.5, 25.5, "Y");
766 hTripplePeakMaximum->GetXaxis()->SetTitle( "Timeslices" );
767 hTripplePeakMaximum->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
768 hList.Add( hTripplePeakMaximum );
769 //-----------------------------------------------------------------------------
770
771 hLargePeakOverlay= new TH2F("hLargePeakOverlay", "Overlay of detected Multi Photon Peaks",
772 gPeakOverlayXaxisLeft + hPeakOverlayXaxisRight , (-1*gPeakOverlayXaxisLeft)-0.5, hPeakOverlayXaxisRight-0.5 ,
773 512, -5.5, 300.5 );
774 hLargePeakOverlay->SetAxisRange(-5.5, 200.5, "Y");
775 hLargePeakOverlay->GetXaxis()->SetTitle( "Timeslices" );
776 hLargePeakOverlay->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
777 //hLargePeakOverlay->SetBit(TH2F::kCanRebin);
778 hList.Add( hLargePeakOverlay );;
779
780 hLargePeakMaximum = new TH1F("hLargePeakMaximum", "Peak derived with Maximum of above Spektrum",
781 gPeakOverlayXaxisLeft + hPeakOverlayXaxisRight ,
782 (-1*gPeakOverlayXaxisLeft)-0.5,
783 hPeakOverlayXaxisRight-0.5
784 );
785 hLargePeakMaximum->SetAxisRange(-0.5, 50.5, "Y");
786 hLargePeakMaximum->GetXaxis()->SetTitle( "Timeslices" );
787 hLargePeakMaximum->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
788 hList.Add( hLargePeakMaximum );
789 //-----------------------------------------------------------------------------
790
791// hPixelPeakOverlay = new TH2F("hPixelPeakOverlay", "Maximum of Statistic of overlayed Peaks",
792// gPeakOverlayXaxisLeft + hPeakOverlayXaxisRight , (-1*gPeakOverlayXaxisLeft)-0.5, hPeakOverlayXaxisRight-0.5 ,
793// 512, -55.5, 200.5 );
794// hPixelPeakOverlay->GetXaxis()->SetTitle( "Timeslices" );
795// hPixelPeakOverlay->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
796// hList.Add( hPixelPeakOverlay );
797
798// hEventPeakOverlay = new TH2F("hEventPeakOverlay", "Overlay of detected Peaks of all Pixel of one Event",
799// gPeakOverlayXaxisLeft + hPeakOverlayXaxisRight , (-1*gPeakOverlayXaxisLeft)-0.5, hPeakOverlayXaxisRight-0.5 ,
800// 4096, -48.5, 200.5 );
801// hEventPeakOverlay->GetXaxis()->SetTitle( "Timeslices" );
802// hEventPeakOverlay->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
803// hList.Add( hEventPeakOverlay );
804
805}
806
807void SaveHistograms(const char * loc_fname ){
808 // create the histogram file (replace if already existing)
809 TFile tf( loc_fname, "RECREATE");
810
811 hList.Write(); // write the major histograms into the top level directory
812
813 tf.Close(); // close the file
814} // end of SaveHistograms( char * loc_fname )
815
Note: See TracBrowser for help on using the repository browser.