source: fact/tools/rootmacros/tpeak.C@ 12385

Last change on this file since 12385 was 12366, checked in by Jens Buss, 14 years ago
This a modified Version of dNeises fpeak, which produces a overlay-histogram of all darkcount-peaks
File size: 31.8 KB
Line 
1/********************** Peak-Template ***********************
2
3this makro shell search for Peaks in data and form a template
4of these, so it can be used for detection of other peaks
5
6*************************************************************/
7
8//root libraries
9
10#include <TROOT.h>
11#include <TCanvas.h>
12#include <TProfile.h>
13#include <TTimer.h>
14#include <TH1F.h>
15#include <TH2F.h>
16#include <Getline.h>
17#include <TLine.h>
18#include <TBox.h>
19#include <TMath.h>
20#include <TFile.h>
21#include <TStyle.h>
22
23#include <stdio.h>
24#include <stdint.h>
25#include <cstdio>
26
27#define NPIX 1440
28#define NCELL 1024
29#define FAD_MAX_SAMPLES 1024
30
31#define HAVE_ZLIB
32
33//rootmacros
34
35#include "fits.h"
36#include "FOpenCalibFile.c"
37
38#include "discriminator.h"
39#include "discriminator.C"
40#include "zerosearch.h"
41#include "zerosearch.C"
42#include "factfir.C"
43
44
45//-----------------------------------------------------------------------------
46// Decleration of variables
47//-----------------------------------------------------------------------------
48
49vector<int16_t> AllPixelDataVector;
50vector<int16_t> StartCellVector;
51
52unsigned int CurrentEventID;
53
54bool breakout=false;
55
56size_t ROIxNOP;
57UInt_t NumberOfPixels;
58UInt_t RegionOfInterest;
59int NEvents;
60
61size_t drs_n;
62vector<float> drs_basemean;
63vector<float> drs_gainmean;
64vector<float> drs_triggeroffsetmean;
65
66int FOpenDataFile( fits & );
67
68
69vector<float> Ameas(FAD_MAX_SAMPLES); // copy of the data (measured amplitude
70vector<float> N1mean(FAD_MAX_SAMPLES); // mean of the +1 -1 ch neighbors
71vector<float> Vcorr(FAD_MAX_SAMPLES); // corrected Values
72vector<float> Vdiff(FAD_MAX_SAMPLES); // numerical derivative
73
74vector<float> Vslide(FAD_MAX_SAMPLES); // sliding average result
75vector<float> Vcfd(FAD_MAX_SAMPLES); // CDF result
76vector<float> Vcfd2(FAD_MAX_SAMPLES); // CDF result + 2nd sliding average
77
78
79float getValue( int, int );
80float correctDrsOffset( int slice, int pixel );
81void computeN1mean( int );
82void removeSpikes( int );
83
84// histograms
85const int NumberOfDebugHistoTypes = 7;
86const unsigned int
87 Ameas_ = 0,
88 N1mean_ = 1,
89 Vcorr_ = 2,
90 Vtest_ = 3,
91 Vslide_ = 4,
92 Vcfd_ = 5,
93 Vcfd2_ = 6;
94
95TH1F* h;
96TH1F *debugHistos;
97TH2F* hStartCell; // id of the DRS physical pipeline cell where readout starts
98 // x = pixel id, y = DRS cell id
99
100TH1F *hBaseline[ NPIX ]; // histograms for baseline extraction
101TH1F *hMeanBsl, *hpltMeanBsl;
102TH1F *hRmsBsl, *hpltRmsBsl;
103TH2F * hAmplSpek_cfd;
104TH2F * hAmplSpek_discri;
105
106TH2F *hPeakOverlay; //histogrammm for overlay of detected Peaks
107
108TObjArray hList;
109TObjArray hListBaseline;
110
111void BookHistos( );
112void SaveHistograms( char * );
113
114int searchSinglesPeaks = 0;
115
116//-----------------------------------------------------------------------------
117// Main
118//-----------------------------------------------------------------------------
119
120int tpeak(
121 char *datafilename = "../../20111011_055.fits.gz",
122 const char *drsfilename = "../../20111011_054.drs.fits.gz",
123 int PixelID = -1,
124 int firstevent = 0,
125 int nevents = -1,
126 bool spikeDebug = false,
127 int verbosityLevel = 1 // different verbosity levels can be implemented here
128 )
129{
130
131//-----------------------------------------------------------------------------
132// Create (pointer to) Canvases, which are used in every run,
133// also in 'non-debug' runs
134//-----------------------------------------------------------------------------
135
136 TCanvas * cSpektrum;
137 TCanvas * cStartCell;
138 cSpektrum = new TCanvas("cSpektrum","Amplitude Spektra of different discriminators",10,10,400,400);
139 cSpektrum->Divide(1,2);
140 cStartCell = new TCanvas("cStartCell ", "The Startcells of this run", 10,410,400,400);
141
142 // Canvases only need if spike Debug, but I want to deklare
143 // the pointers anyway ...
144 TCanvas *cRawAndSpikeRemoval = NULL;
145 TCanvas *cFiltered = NULL;
146 if (spikeDebug){
147 cRawAndSpikeRemoval = new TCanvas("cRawAndSpikeRemoval","DRS Waveform",410,10,400,400);
148 cRawAndSpikeRemoval->Divide(1, 2);
149 cFiltered = new TCanvas("cFiltered","filtered DRS Waveforms",410,410,400,400);
150 cFiltered->Divide(1, 2);
151 }
152 // Canvases to show the peak template
153 TCanvas *cPeakOverlay = NULL;
154 cPeakOverlay = new TCanvas("cPeakOverlay", "Overlay of detected Peaks", 410, 10, 600, 600);
155
156gStyle->SetPalette(1,0);
157gROOT->SetStyle("Plain");
158
159//-----------------------------------------------------------------------------
160// read FACT raw data
161// * remove spikes
162// * calculate baseline
163// * find peaks (CFD and normal discriminator)
164// * compute pulse height and pulse integral spektrum of the peaks
165//-----------------------------------------------------------------------------
166
167
168//-----------------------------------------------------------------------------
169// Filter-Settings
170//-----------------------------------------------------------------------------
171
172// sliding window filter settings
173 int k_slide = 16;
174 vector<double> a_slide(k_slide, 1); //erzeugt vector der länge k_slide, und füllt die elemente mit 1
175 double b_slide = k_slide;
176
177// CFD filter settings
178 int k_cfd = 10;
179 vector<double> a_cfd(k_cfd, 0);
180 double b_cfd = 1.;
181 a_cfd[0]=-0.75;
182 a_cfd[k_cfd-1]=1.;
183
184// 2nd slinding window filter
185 int ks2 = 16;
186 vector<double> as2(ks2, 1);
187 double bs2 = ks2;
188
189//-----------------------------------------------------------------------------
190// prepare datafile
191//-----------------------------------------------------------------------------
192
193// Open the data file
194 fits *datafile = new fits( datafilename );
195 if (!datafile) {
196 printf( "Could not open the file: %s\n", datafilename );
197 return 1;
198 }
199
200// access data
201 NEvents = FOpenDataFile( *datafile );
202 printf("number of events in file: %d\n", NEvents);
203 if ( nevents == -1 || nevents > NEvents ) nevents = NEvents; // -1 means all!
204
205//Get the DRS calibration
206 FOpenCalibFile( drsfilename,
207 drs_basemean,
208 drs_gainmean,
209 drs_triggeroffsetmean,
210 drs_n);
211
212//Check the sizes of the data columns
213 if (drs_n != 1474560) {
214 cerr << "error: DRS calib file has wrong ...erm...size ... drs_n is: "
215 << drs_n << endl;
216 cerr << " Aborting." << endl;
217 return 1;
218 }
219
220 if(ROIxNOP != 1474560)
221 {
222 cout << "warning: data_n should better be 1440x1024=1474560, but is "
223 << ROIxNOP << endl;
224 cout << "this script is not guaranteed to run under these "
225 <<" circumstances....any way ... it is never guaranteed." << endl;
226 }
227
228// Book the histograms
229
230 BookHistos( );
231
232 float myTHR = 3.5;
233 float calibratedVoltage;
234//-----------------------------------------------------------------------------
235// Loops over Every Event and Pixel
236//-----------------------------------------------------------------------------
237 for ( int ev = firstevent; ev < firstevent + nevents; ev++) {
238 // Get an Event --> consists of 1440 Pixel ...erm....data
239 datafile->GetRow( ev );
240
241 //-------------------------------------
242 // Loop over every Pixel of Event
243 //-------------------------------------
244
245 for ( int pix = 0; pix < 1440; pix++ ){
246
247 hStartCell->Fill( pix, StartCellVector[pix] );
248
249 // this is a stupid hack ... there is more code at the
250 // end of this loop to complete this hack ...
251 // beginning with if (PixelID != -1) as well.
252 if (PixelID != -1) {
253 pix = PixelID;
254 if (verbosityLevel > 0){
255 cout << "Processing Event number: " << CurrentEventID << "\t"
256 << "Pixel number: "<< pix << endl;
257 }
258 }
259
260 if (verbosityLevel > 0){
261 if (pix % 20 ==0){
262 cout << "Processing Event number: " << CurrentEventID << "\t"
263 << "Pixel number: "<< pix << endl;
264 }
265 }
266
267 // compute the DRs calibrated values and put them into Ameas[]
268 for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
269 calibratedVoltage = getValue( sl, pix);
270 if (verbosityLevel > 10){
271 printf("calibratedVoltage = %f\n", calibratedVoltage);
272 }
273 Ameas[ sl ] = calibratedVoltage;
274
275 // in case we want to plot this ... we need to put it into a
276 // Histgramm
277 if (spikeDebug){
278 debugHistos[Ameas_].SetBinContent(sl, calibratedVoltage);
279 }
280 }
281 // operates on Ameas[] and writes to N1mean[]
282 computeN1mean( RegionOfInterest );
283
284 // operates on Ameas[] and N1mean[], then writes to Vcorr[]
285 removeSpikes( RegionOfInterest );
286
287 // filter Vcorr with sliding average using FIR filter function
288 factfir(b_slide , a_slide, k_slide, Vcorr, Vslide);
289 // filter Vslide with CFD using FIR filter function
290 factfir(b_cfd , a_cfd, k_cfd, Vslide, Vcfd);
291 // filter Vcfd with sliding average using FIR filter function
292 factfir(bs2 , as2, ks2, Vcfd, Vcfd2);
293
294
295 vector<Region> *Regions = discriminator( Vslide, myTHR,true , 120);
296 for (unsigned int p=0; p<Regions->size(); p++ ){
297 hAmplSpek_discri->Fill(pix , Regions->at(p).maxVal);
298 }
299
300 // peaks in Ameas[] are found by searching for zero crossings
301 // in Vcfd2
302 // first Argument 1 means ... *rising* edge
303 // second Argument 1 means ... search with stepsize 1 ... 10 is okay as well
304 vector<Region> * zXings = zerosearch( Vcfd2 , 1 , 1);
305 // zXings means "zero cross ings"
306 ShiftRegionBy(*zXings, -ks2/2);
307 EnlargeRegion(*zXings, 10, 10);
308 findAbsMaxInRegions(*zXings, Vslide);
309
310 for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
311 // hPixelCellData.Fill( sl, Vcorr[sl] );
312 hBaseline[pix]->Fill( Vslide[sl] ) ;
313 }
314
315 //following Code should produce the Overlay of peaks
316 vector<Region>::iterator it;
317
318 for (it = zXings->begin() ; it < zXings->end() ; it++){
319 int a = it->maxPos-100;
320 int b = it->maxPos+100;
321 if (a < 0)
322 a =0;
323 if (b>1023)
324 b=1023;
325 for ( int pos = a; pos <= b; pos++){
326
327 hPeakOverlay->Fill( pos - (it->maxPos), Vcorr[pos]) ;
328 }
329 }
330
331
332
333 if (zXings->size() != 0 ){
334 for (unsigned int i=0; i<zXings->size(); i++){
335 if (verbosityLevel > 1){
336 cout << zXings->at(i).maxPos << ":\t"
337 << zXings->at(i).maxVal <<endl;
338 }
339 hAmplSpek_cfd->Fill(pix, zXings->at(i).maxVal);
340 }
341 }
342
343 if ( spikeDebug ){
344 for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
345 debugHistos[Vslide_].SetBinContent( sl, Vslide[sl] );
346 debugHistos[Vcfd_].SetBinContent( sl, Vcfd[sl] );
347 debugHistos[Vcfd2_].SetBinContent( sl, Vcfd2[sl] );
348 }
349
350 cRawAndSpikeRemoval->cd( 1);
351 gPad->SetGrid();
352 debugHistos[Ameas_].Draw();
353
354 cRawAndSpikeRemoval->cd( 2);
355 gPad->SetGrid();
356 debugHistos[Vcorr_].Draw();
357
358 cFiltered->cd(1);
359 gPad->SetGrid();
360 debugHistos[Vslide_].Draw();
361 TLine * thrLine = new TLine(0, myTHR, 1024, myTHR);
362 thrLine->SetLineColor(kRed);
363 thrLine->SetLineWidth(2);
364 thrLine->Draw();
365 TLine * OneLine;
366 vector<TLine*> MyLines;
367 for (unsigned int p=0; p<Regions->size(); p++ ){
368 OneLine = new TLine(
369 Regions->at(p).begin ,
370 Regions->at(p).maxVal,
371 Regions->at(p).end,
372 Regions->at(p).maxVal);
373 OneLine->SetLineColor(kRed);
374 OneLine->SetLineWidth(2);
375 MyLines.push_back(OneLine);
376 OneLine->Draw();
377 }
378 TBox *OneBox;
379 vector<TBox*> MyBoxes;
380 for (unsigned int i=0; i<zXings->size(); i++){
381 OneBox = new TBox(
382 zXings->at(i).maxPos -10 ,
383 zXings->at(i).maxVal -0.5,
384 zXings->at(i).maxPos +10 ,
385 zXings->at(i).maxVal +0.5);
386 OneBox->SetLineColor(kBlue);
387 OneBox->SetLineWidth(1);
388 OneBox->SetFillStyle(0);
389 OneBox->SetFillColor(kRed);
390 MyBoxes.push_back(OneBox);
391 OneBox->Draw();
392 }
393
394 cFiltered->cd(2);
395 gPad->SetGrid();
396 debugHistos[Vcfd2_].Draw();
397 TLine *zeroline = new TLine(0, 0, 1024, 0);
398 zeroline->SetLineColor(kBlue);
399 zeroline->Draw();
400
401 cRawAndSpikeRemoval->Update();
402 cFiltered->Update();
403
404 //OVERLAY PEAKS
405 cPeakOverlay->cd();
406 hPeakOverlay->Draw("COLZ");
407 cPeakOverlay->Modified();
408 cPeakOverlay->Update();
409
410 //Process gui events asynchronously during input
411 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
412 timer.TurnOn();
413 TString input = Getline("Type 'q' to exit, <return> to go on: ");
414 timer.TurnOff();
415 if (input=="q\n") {
416 breakout=true;
417 }
418
419 //TODO!!!!!!!!!
420 // do some Garbage collection ...
421 // all the Objects on the heap should be deleted here.
422
423 }// end of if(spikeDebug)
424
425 delete Regions;
426 delete zXings;
427
428 // this is the 2nd part of the ugly hack in the beginning.
429 // this makes sure, that the loop ends
430 if (PixelID != -1){
431 pix = 2000;
432 }
433
434 if (breakout)
435 break;
436
437 }
438 //-------------------------------------
439 // end of loop over pixels
440 //-------------------------------------
441if (ev % 50 == 0){
442 cSpektrum->cd(1);
443 hAmplSpek_cfd->Draw("COLZ");
444 cSpektrum->cd(2);
445 hAmplSpek_discri->Draw("COLZ");
446 cSpektrum->Modified();
447 cSpektrum->Update();
448
449 // updating seems not to work ..
450 // debug cout
451 cStartCell->cd();
452 hStartCell->Draw();
453 cStartCell->Modified();
454 cStartCell->Update();
455
456 //OVERLAY PEAKS
457 cPeakOverlay->cd();
458 hPeakOverlay->Draw("COLZ");
459 cPeakOverlay->Modified();
460 cPeakOverlay->Update();
461
462}
463
464
465 if (breakout)
466 break;
467 }
468//-------------------------------------
469// end of loop over Events
470//-------------------------------------
471
472
473 for (int pix = 0; pix < NumberOfPixels; pix++) {
474 //printf("Maximum bin pix[%d] %f\n", pix ,hBaseline[pix]->GetMaximumBin() );
475
476 float vmaxprob = hBaseline[pix]->GetXaxis()->GetBinCenter(
477 hBaseline[pix]->GetMaximumBin() );
478
479 printf("%8.3f", vmaxprob );
480
481 hMeanBsl->Fill( vmaxprob );
482 hpltMeanBsl->SetBinContent(pix+1, vmaxprob );
483
484 hRmsBsl->Fill(hBaseline[pix]->GetRMS() );
485 hpltRmsBsl->SetBinContent( pix+1, hBaseline[pix]->GetRMS() );
486 }
487
488
489 SaveHistograms( datafilename );
490
491 return( 0 );
492}
493
494void removeSpikes(int Samples){
495
496 const float fract = 0.8;
497 float x, xp, xpp, x3p;
498
499 // assume that there are no spikes
500 for ( int i = 0; i < Samples; i++) Vcorr[i] = Ameas[i];
501
502// find the spike and replace it by mean value of neighbors
503 for ( int i = 0; i < Samples; i++) {
504
505 // printf("Vcorr[%d] = %f, Ameas[%d] = %f\n", i, Vcorr[ i ], i, Ameas[ i ] );
506
507 x = Ameas[i] - N1mean[i];
508
509 if ( x < -5. ){ // a spike candidate
510 // check consistency with a single channel spike
511 xp = Ameas[i+1] - N1mean[i+1];
512 xpp = Ameas[i+2] - N1mean[i+2];
513 x3p = Ameas[i+3] - N1mean[i+3];
514
515 // printf("candidates x[%d] = %f; xp = %f; xpp = %f, x3p = %f\n", i, x, xp, xpp, x3p);
516
517 if ( Ameas[i+2] - ( Ameas[i] + Ameas[i+3] )/2. > 10. ){
518 // printf("double spike candidate\n");
519 Vcorr[i+1] = ( Ameas[i] + Ameas[i+3] )/2.;
520 Vcorr[i+2] = ( Ameas[i] + Ameas[i+3] )/2.;
521 // printf("Vcorr[%d] = %f %f %f %f\n", i, Vcorr[i], Vcorr[i+1], Vcorr[i+2], Vcorr[ i+3 ]);
522 // printf("Ameas[%d] = %f %f %f %f\n", i, Ameas[ i ], Ameas[ i+1 ], Ameas[ i+2 ], Ameas[i+3]);
523 i = i + 3;
524 }
525 else{
526
527 if ( ( xp > -2.*x*fract ) && ( xpp < -10. ) ){
528 Vcorr[i+1] = N1mean[i+1];
529 // printf("Vcorr[%d] = %f %f %f\n", i, Vcorr[i], Vcorr[i+1], Vcorr[i+2]);
530 // N1mean[i+1] = (Ameas[i] - Ameas[i+2] / 2.);
531 N1mean[i+2] = (Ameas[i+1] - Ameas[i+3] / 2.);
532 i = i + 2;//do not care about the next sample it was the spike
533 }
534 // treatment for the end of the pipeline must be added !!!
535 }
536 }
537 else{
538 // do nothing
539 }
540 } // end of spike search and correction
541 for ( int i = 0; i < Samples; i++ ) debugHistos[ Vcorr_ ].SetBinContent( i, Vcorr[i] );
542}
543/*
544void computeN1mean( int Samples ){
545cout << "In compute N1mean" << endl;
546// compute the mean of the left and right neighbors of a channel
547
548 for( int i = 0; i < Samples; i++){
549 if (i == 0){ // use right sample es mean
550 N1mean[i] = Ameas[i+1];
551 }
552 else if ( i == Samples-1 ){ //use left sample as mean
553 N1mean[i] = Ameas[i-1];
554 }
555 else{
556 N1mean[i] = ( Ameas[i-1] + Ameas[i+1] ) / 2.;
557 }
558 h[tN1mean].SetBinContent(i, Ameas[i] - N1mean[i]);
559 }
560} // end of computeN1mean computation
561*/
562
563void computeN1mean( int Samples ){
564// compute the mean of the left and right neighbors of a channel
565
566 for( int i = 2; i < Samples - 2; i++){
567/* if (i == 0){ // use right sample es mean
568 N1mean[i] = Ameas[i+1];
569 }
570 else if ( i == Samples-1 ){ //use left sample as mean
571 N1mean[i] = Ameas[i-1];
572 }
573 else{
574 N1mean[i] = ( Ameas[i-1] + Ameas[i+1] ) / 2.;
575 }
576*/
577 N1mean[i] = ( Ameas[i-1] + Ameas[i+1] ) / 2.;
578 }
579} // end of computeN1mean computation
580
581
582
583/*
584// this function takes a vector<float> and tries to apply the DRS calibration to it.
585// therfor it uses a DRS calibration file, which must be given as a parameter.
586// in case the calibration file can not be opened the function generates a warning on cerr
587// and returns an empty vector
588//
589// the function assumes the first parameter conatins the data of all Pixel
590// the RegionOfInterest is deduced from its size().
591
592 //
593
594vector<float> * calibrateDrsData(
595 vector<float> & AllPixelData,
596 vector<float> & AllStartCells,
597 const char * DrsCalibFilePath,
598 bool ApplyCalibToSourceVector = false)
599{
600 vector<float> * CalibratedDRSData;
601
602 size_t DrsLengthOfVectors;
603 vector<float> DrsOffset;
604 vector<float> DrsGain;
605 vector<float> DrsOffsetSpecialTrigger;
606
607 // the Total Number of Pixel, we have data of
608 // is deduced from the size() of AllPixelStartCell
609 unsigned int NumberOfPixel = AllStartCells.size();
610 //sanity Check:
611 if (NumberOfPixel < 1) {
612 cerr << "calibrateDrsData - Error: NumberOfPixel=AllStartCells.size() ia:"
613 << AllStartCells.size() << " must be >= 1" << endl;
614 return CalibratedDRSData
615 }
616
617 // The RegionOfInterest is deduced from the size() of AllPixel Data
618 unsigned int RegionOfInterest = AllPixelData.size() / NumberOfPixel;
619 // sanity Check
620 if ( RegionOfInterest < 1){
621 cerr << "calibrateDrsData - Error: RegionOfInterest = AllPixelData.size() / NumberOfPixel is:"
622 << RegionOfInterest << " must be >= 1" << endl;
623 return CalibratedDRSData
624 }
625
626 // now lets see if we can find the drs calib file
627 FOpenCalibFile( DrsCalibFilePath,
628 DrsOffset,
629 DrsGain,
630 DrsOffsetSpecialTrigger,
631 DrsLengthOfVectors);
632 //sanity check is done in FOpenCalibFile, I guess ... check and maybe TODO
633 // maybe at least the return code of FOpenCalibFile should be checked here ... TODO
634
635
636
637
638}
639
640*/
641
642float getValue( int slice, int pixel ){
643 const float dconv = 2000/4096.0;
644
645 float vraw, vcal;
646
647 unsigned int pixel_pt;
648 unsigned int slice_pt;
649 unsigned int cal_pt;
650 unsigned int drs_cal_offset;
651
652 // printf("pixel = %d, slice = %d\n", slice, pixel);
653
654 pixel_pt = pixel * RegionOfInterest;
655 slice_pt = pixel_pt + slice;
656 drs_cal_offset = ( slice + StartCellVector[ pixel ] )%RegionOfInterest;
657 cal_pt = pixel_pt + drs_cal_offset;
658
659 vraw = AllPixelDataVector[ slice_pt ] * dconv;
660 vcal = ( vraw - drs_basemean[ cal_pt ] - drs_triggeroffsetmean[ slice_pt ] ) / drs_gainmean[ cal_pt ]*1907.35;
661
662 return( vcal );
663}
664
665// this is a faster version of getValue, that does not do a full calibration
666// was used for a performace test...
667float correctDrsOffset( int slice, int pixel ){
668
669 const float dconv = 2000/4096.0;
670
671 // here 1024 is not the RegionOfInterest, but really the lenth of the pipeline
672 unsigned int physical_slice = ( slice + StartCellVector[ pixel ] ) % 1024;
673
674 unsigned int slice_pt;
675 unsigned int physical_slice_pt;
676 slice_pt = pixel * RegionOfInterest + slice;
677 physical_slice_pt = pixel * RegionOfInterest + physical_slice;
678
679 float vcal = AllPixelDataVector[ slice_pt ] *
680 dconv - drs_basemean[ physical_slice_pt ];
681 return( vcal );
682}
683
684// booking and parameter settings for all histos
685void BookHistos( ){
686 // histograms for baseline extraction
687 char hName[500];
688 char hTitle[500];
689
690 TH1F *h;
691
692 for( int i = 0; i < NPIX; i++ ) {
693 sprintf(&hTitle[0],"all events all slices of pixel %d", i);
694 sprintf(&hName[0],"base%d", i);
695
696 h = new TH1F( hName, hTitle, 400, -99.5 ,100.5 );
697
698 h->GetXaxis()->SetTitle( "Sample value (mV)" );
699 h->GetYaxis()->SetTitle( "Entries / 0.5 mV" );
700 hListBaseline.Add( h );
701 hBaseline[i] = h;
702 }
703
704 hMeanBsl = new TH1F("histo_mean","Value of maximal probability",400,-99.5,100.5);
705 hMeanBsl->GetXaxis()->SetTitle( "max value (mV)" );
706 hMeanBsl->GetYaxis()->SetTitle( "Entries / 0.5 mV" );
707 hList.Add( hMeanBsl );
708
709 hpltMeanBsl = new TH1F("hplt_mean","Value of maximal probability",1440,-0.5,1439.5);
710 hpltMeanBsl->GetXaxis()->SetTitle( "pixel" );
711 hpltMeanBsl->GetYaxis()->SetTitle( "max value in mV" );
712 hList.Add( hpltMeanBsl );
713
714 hRmsBsl = new TH1F("histo_rms","RMS in mV",2000,-99.5,100.5);
715 hRmsBsl->GetXaxis()->SetTitle( "RMS (mV)" );
716 hRmsBsl->GetYaxis()->SetTitle( "Entries / 0.5 mV" );
717 hList.Add( hRmsBsl );
718
719 hpltRmsBsl = new TH1F("hplt_rms","Value of maximal probability",1440,-0.5,1439.5);
720 hpltRmsBsl->GetXaxis()->SetTitle( "pixel" );
721 hpltRmsBsl->GetYaxis()->SetTitle( "RMS in mV" );
722 hList.Add( hpltRmsBsl );
723
724 hAmplSpek_cfd = new TH2F("hAmplSpek_cfd","amplitude spektrum - CFD",1440,-0.5,1439.5, 256, -27.5, 100.5);
725 hAmplSpek_cfd->GetXaxis()->SetTitle( "pixel" );
726 hAmplSpek_cfd->GetYaxis()->SetTitle( "amplitude in mV" );
727 hList.Add( hAmplSpek_cfd );
728
729 hAmplSpek_discri = new TH2F("hAmplSpek_discri","amplitude spektrum - std discriminator",1440,-0.5,1439.5, 256, -27.5, 100.5);
730 hAmplSpek_discri->GetXaxis()->SetTitle( "pixel" );
731 hAmplSpek_discri->GetXaxis()->SetTitle( "amplitude in mV" );
732 hList.Add( hAmplSpek_discri );
733
734 hStartCell = new TH2F("StartCell", "StartCell", 1440, 0., 1440., 1024, 0., 1024);
735 hStartCell->GetXaxis()->SetTitle( "pixel" );
736 hStartCell->GetXaxis()->SetTitle( "slice" );
737 hList.Add( hStartCell );
738
739 debugHistos = new TH1F[ NumberOfDebugHistoTypes ];
740 for ( int type = 0; type < NumberOfDebugHistoTypes; type++){
741 debugHistos[ type ].SetBins(1024, 0, 1024);
742 debugHistos[ type ].SetLineColor(1);
743 debugHistos[ type ].SetLineWidth(2);
744
745 // set X axis paras
746 debugHistos[ type ].GetXaxis()->SetLabelSize(0.1);
747 debugHistos[ type ].GetXaxis()->SetTitleSize(0.1);
748 debugHistos[ type ].GetXaxis()->SetTitleOffset(1.2);
749 debugHistos[ type ].GetXaxis()->SetTitle(Form("Time slice (%.1f ns/slice)", 1./2.));
750
751 // set Y axis paras
752 debugHistos[ type ].GetYaxis()->SetLabelSize(0.1);
753 debugHistos[ type ].GetYaxis()->SetTitleSize(0.1);
754 debugHistos[ type ].GetYaxis()->SetTitleOffset(0.3);
755 debugHistos[ type ].GetYaxis()->SetTitle("Amplitude (a.u.)");
756 }
757 hPeakOverlay = new TH2F("hPeakOverlay", "Overlay of detected Peaks", 201, -100.5, 100.5, 1024, -99.5, 100.5);
758 hPeakOverlay->GetXaxis()->SetTitle( "Timeslices" );
759 hPeakOverlay->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
760 //hList.Add( hPeakOverlay h);
761
762}
763
764void SaveHistograms( char * loc_fname ){
765
766 TString fName; // name of the histogram file
767
768 // create the filename for the histogram file
769 fName = loc_fname; // use the name of the tree file
770 // TODO ... next statement doesn't work for ".fits.gz"
771 fName.Remove(fName.Length() - 5); // remove the extension .fits
772 fName += "_discri.root"; // add the new extension
773
774 // create the histogram file (replace if already existing)
775 TFile tf( fName, "RECREATE");
776
777 hList.Write(); // write the major histograms into the top level directory
778 tf.mkdir("BaselineHisto");
779 tf.cd("BaselineHisto"); // go to new subdirectory
780 hListBaseline.Write(); // write histos into subdirectory
781
782 tf.Close(); // close the file
783} // end of SaveHistograms( char * loc_fname )
784
785int FOpenDataFile(fits &datafile){
786
787// cout << "-------------------- Data Header -------------------" << endl;
788// datafile.PrintKeys();
789// cout << "------------------- Data Columns -------------------" << endl;
790// datafile.PrintColumns();
791
792 //Get the size of the data column
793 RegionOfInterest = datafile.GetUInt("NROI");
794 NumberOfPixels = datafile.GetUInt("NPIX");
795 //Size of column "Data" = #Pixel x ROI
796 ROIxNOP = datafile.GetN("Data");
797
798 //Set the sizes of the data vectors
799 AllPixelDataVector.resize(ROIxNOP,0);
800 StartCellVector.resize(NumberOfPixels,0);
801
802 //Link the data to variables
803 datafile.SetRefAddress("EventNum", CurrentEventID);
804 datafile.SetVecAddress("Data", AllPixelDataVector);
805 datafile.SetVecAddress("StartCellData", StartCellVector);
806 datafile.GetRow(0);
807
808 return datafile.GetNumRows() ;
809}
810
811
812/////////////////////////////////////////////////////////////////////////////
813/// old BookHistos
814/*
815void BookHistos( int Samples ){
816// booking and parameter settings for all histos
817
818 h = new TH1F[ Ntypes ];
819
820 for ( int type = 0; type < Ntypes; type++){
821
822 h[ type ].SetBins(Samples, 0, Samples);
823 h[ type ].SetLineColor(1);
824 h[ type ].SetLineWidth(2);
825
826 // set X axis paras
827 h[ type ].GetXaxis()->SetLabelSize(0.1);
828 h[ type ].GetXaxis()->SetTitleSize(0.1);
829 h[ type ].GetXaxis()->SetTitleOffset(1.2);
830 h[ type ].GetXaxis()->SetTitle(Form("Time slice (%.1f ns/slice)", 1./2.));
831
832 // set Y axis paras
833 h[ type ].GetYaxis()->SetLabelSize(0.1);
834 h[ type ].GetYaxis()->SetTitleSize(0.1);
835 h[ type ].GetYaxis()->SetTitleOffset(0.3);
836 h[ type ].GetYaxis()->SetTitle("Amplitude (a.u.)");
837 }
838 cRawAndSpikeRemoval = new TCanvas("cRawAndSpikeRemoval","DRS Waveform",10,10,800,600);
839 cRawAndSpikeRemoval->Divide(1, 3);
840 cFilter = new TCanvas("cFilter","filtered DRS Waveforms",10,10,800,600);
841 cFilter->Divide(1, 3);
842
843 hStartCell = new TH2F("StartCell", "StartCell", 1440, 0., 1440., 1024, 0., 1024);
844
845}
846*/
847
Note: See TracBrowser for help on using the repository browser.