source: fact/tools/rootmacros/fpeak.C@ 14535

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