source: fact/tools/rootmacros/fbsl.C@ 12352

Last change on this file since 12352 was 12307, checked in by neise, 13 years ago
added '#define HAVE_ZLIB'
File size: 12.7 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 <TMath.h>
10#include <TFile.h>
11
12#include <stdio.h>
13#include <stdint.h>
14#include <cstdio>
15
16#define HAVE_ZLIB
17#include "fits.h"
18//#include "TPKplotevent.c"
19//#include "FOpenDataFile.c"
20#include "FOpenCalibFile.c"
21
22#include "zerosearch.C"
23
24#define NPIX 1440
25#define NCELL 1024
26
27// data access and treatment
28#define FAD_MAX_SAMPLES 1024
29
30vector<int16_t> data;
31vector<int16_t> data_offset;
32
33unsigned int data_num;
34
35size_t data_n;
36UInt_t data_px;
37UInt_t data_roi;
38int NEvents;
39int NBSLeve = 1000;
40
41size_t drs_n;
42vector<float> drs_basemean;
43vector<float> drs_gainmean;
44vector<float> drs_triggeroffsetmean;
45
46int FOpenDataFile( fits & );
47
48
49vector<float> Ameas(FAD_MAX_SAMPLES); // copy of the data (measured amplitude
50vector<float> N1mean(FAD_MAX_SAMPLES); // mean of the +1 -1 ch neighbors
51vector<float> N2mean(FAD_MAX_SAMPLES); // mean of the +2 -2 ch neighbors
52vector<float> Vcorr(FAD_MAX_SAMPLES); // corrected Values
53vector<float> Vdiff(FAD_MAX_SAMPLES); // numerical derivative
54
55vector<float> Vslide(FAD_MAX_SAMPLES); // sliding average result
56vector<float> Vcfd(FAD_MAX_SAMPLES); // CDF result
57vector<float> Vcfd2(FAD_MAX_SAMPLES); // CDF result + 2nd sliding average
58
59#include "factfir.C"
60
61float getValue( int, int );
62void computeN1mean( int );
63void removeSpikes( int );
64
65// histograms
66const int Ntypes = 7;
67const unsigned int // arranged by Dominik
68 tAmeas = 0,
69 tN1mean = 1,
70 tVcorr = 2,
71 tVtest = 3,
72 tVslide = 4,
73 tVcfd = 5,
74 tVcfd2 = 6;
75
76TH1F* h;
77TH2F* hStartCell; // id of the DRS physical pipeline cell where readout starts
78 // x = pixel id, y = DRS cell id
79TH2F hPixelCellData(
80 "PixelPedestal", "PixelPedestal", NCELL, 0., NCELL, 200, -50., 150.);
81TH1F *hBaseline[ NPIX ]; // histograms for baseline extraction
82TH1F *hMeanBsl, *hpltMeanBsl;
83TH1F *hRmsBsl, *hpltRmsBsl;
84TObjArray hList;
85TObjArray hListBaseline;
86
87void BookHistos( int );
88void SaveHistograms( char * );
89
90// Create a canvas
91TCanvas* CW;
92TCanvas* cFilter;
93
94int spikeDebug = 0;
95int searchSinglesPeaks = 0;
96
97
98int fbsl(
99 char *datafilename = "../../20110804_024.fits.gz",
100 const char *drsfilename = "../../20110804_023.drs.fits.gz",
101 int pixelnr = 0,
102 int firstevent = 0,
103 int nevents = -1 ){
104// read and analyze FACT raw data
105
106// sliding window filter settings
107 int k_slide = 16;
108 vector<double> a_slide(k_slide, 1);
109 double b_slide = k_slide;
110
111// CFD filter settings
112 int k_cfd = 10;
113 vector<double> a_cfd(k_cfd, 0);
114 double b_cfd = 1.;
115 a_cfd[0]=0.75;
116 a_cfd[k_cfd-1]=-1.;
117
118// 2nd slinding window filter
119 int ks2 = 16;
120 vector<double> as2(ks2, 1);
121 double bs2 = ks2;
122 gROOT->SetStyle("Plain");
123
124//-------------------------------------------
125// Open the file
126//-------------------------------------------
127 fits datafile( datafilename );
128 if (!datafile) {
129 printf( "Could not open the file: %s\n", datafilename );
130 return 1;
131 }
132
133 // access data
134 NEvents = FOpenDataFile( datafile );
135
136 printf("number of events in file: %d\n", NEvents);
137
138 // compare the number of events in the data file with the nevents the
139 // the user would like to read. nevents = -1 means: read all
140 if ( nevents == -1 || nevents > NEvents ) nevents = NEvents;
141
142
143 // FILE *pedFile;
144 // pedFile = fopen("t.txt","u");
145 // fprintf(pedFile,"# Pedestal Data");
146 // fclose( pedFile );
147//-------------------------------------------
148//Get the DRS calibration
149//-------------------------------------------
150
151 FOpenCalibFile(drsfilename, drs_basemean, drs_gainmean, drs_triggeroffsetmean, drs_n);
152
153
154 //Check the sizes of the data columns
155 if(drs_n!=data_n)
156 {
157 cout << "Data and DRS file incompatible (Px*ROI disagree)" << endl;
158 return 1;
159 }
160 // Book the histograms
161 // printf("call BookHistos\n");
162 BookHistos( data_roi );
163 // printf("returned from BookHistos\n");
164 // Loop over events
165 cout << "--------------------- Data --------------------" << endl;
166
167 float value;
168
169 // loop over events
170// for ( int ev = firstevent; ev < nevents; ev++) {
171 for ( int ev = firstevent; ev < firstevent + nevents; ev++) {
172
173 datafile.GetRow( ev );
174
175 if ( ev % 100 == 0 ){
176 cout << "Event number: " << data_num << endl;
177 }
178
179 // loop over pixel
180 for ( int pix = 0; pix < 1440; pix++ ){
181
182 // hStartCell->Fill( pix, data_offset[pix] );
183
184 // loop over DRS slices
185 for ( unsigned int sl = 0; sl < data_roi; sl++){
186
187 Ameas[ sl ] = getValue(sl, pix);
188
189 }
190 // printf("Ameas[100] = %f\n", Ameas[100] );
191
192 computeN1mean( data_roi ); // prepare spike removal
193 removeSpikes( data_roi ); // output in Vcorr
194
195 // filter Vcorr with sliding average using FIR filter function
196 factfir(b_slide , a_slide, k_slide, Vcorr, Vslide);
197
198 for ( unsigned int sl = 0; sl < data_roi; sl++){
199 // hPixelCellData.Fill( sl, Vcorr[sl] );
200 hBaseline[pix]->Fill( Vslide[sl] ) ;
201 }
202 }
203 }
204
205 FILE *fp;
206 TString fName; // name of the histogram file
207 /* create the filename for the histogram file */
208 fName = datafilename; // use the name of the tree file
209 fName.Remove(fName.Length() - 5); // remove the extension .root
210 fName += "_bsl.txt"; // add the new extension
211
212 fp = fopen(fName, "a+");
213 fprintf( fp, "FILE: %s\n", datafilename );
214 fprintf( fp, "NEVENTS: %d\n", NEvents);
215 NBSLeve = NEvents; // this has to be changed when computation on a subset of a run is implemented
216 fprintf( fp, "NBSLEVE: %d\n", NBSLeve );
217
218 for (int pix = 0; pix < NPIX; pix++) {
219 //printf("Maximum bin pix[%d] %f\n", pix ,hBaseline[pix]->GetMaximumBin() );
220
221 float vmaxprob = hBaseline[pix]->GetXaxis()->GetBinCenter(
222 hBaseline[pix]->GetMaximumBin() );
223
224 fprintf( fp, "%8.3f", vmaxprob );
225
226 hMeanBsl->Fill( vmaxprob );
227 hpltMeanBsl->SetBinContent(pix+1, vmaxprob );
228
229 hRmsBsl->Fill(hBaseline[pix]->GetRMS() );
230 hpltRmsBsl->SetBinContent( pix+1, hBaseline[pix]->GetRMS() );
231 }
232
233 fclose( fp );
234
235 SaveHistograms( datafilename );
236
237 TCanvas * cMeanBsl = new TCanvas();
238 cMeanBsl->cd();
239 hMeanBsl->Draw();
240 //canv_mean->Modified();
241 cMeanBsl->Update();
242
243 TCanvas * cRmsBsl = new TCanvas();
244 cRmsBsl->cd();
245 hRmsBsl->Draw();
246 //canv_rms->Modified();
247 cMeanBsl->Update();
248
249 return( 0 );
250}
251
252void removeSpikes(int Samples){
253
254 const float fract = 0.8;
255 float x, xp, xpp, x3p;
256
257 // assume that there are no spikes
258 for ( int i = 0; i < Samples; i++) Vcorr[i] = Ameas[i];
259
260 // find the spike and replace it by mean value of neighbors
261 for ( int i = 2; i < Samples-2 ; i++) {
262
263 x = Ameas[i] - N1mean[i];
264
265 if ( x < -5. ){ // a spike candidate
266 // check consistency with a single channel spike
267 xp = Ameas[i+1] - N1mean[i+1];
268 xpp = Ameas[i+2] - N1mean[i+2];
269
270 if ( Ameas[i+2] - ( Ameas[i] + Ameas[i+3] )/2. > 10. ){
271 // printf("double spike candidate\n");
272 Vcorr[i+1] = ( Ameas[i] + Ameas[i+3] )/2.;
273 Vcorr[i+2] = ( Ameas[i] + Ameas[i+3] )/2.;
274 i = i + 3;
275 }
276 else{
277
278 if ( ( xp > -2.*x*fract ) && ( xpp < -10. ) ){
279 Vcorr[i+1] = N1mean[i+1];
280 // printf("Vcorr[%d] = %f %f %f\n", i, Vcorr[i], Vcorr[i+1], Vcorr[i+2]);
281 // N1mean[i+1] = (Ameas[i] - Ameas[i+2] / 2.);
282 N1mean[i+2] = (Ameas[i+1] - Ameas[i+3] / 2.);
283 i = i + 2;//do not care about the next sample it was the spike
284 }
285 // treatment for the end of the pipeline must be added !!!
286 }
287 }
288 else{
289 // do nothing
290 }
291 } // end of spike search and correction
292}
293
294void computeN1mean( int Samples ){
295// compute the mean of the left and right neighbors of a channel
296
297 for( int i = 2; i < Samples - 2; i++){
298/* if (i == 0){ // use right sample es mean
299 N1mean[i] = Ameas[i+1];
300 }
301 else if ( i == Samples-1 ){ //use left sample as mean
302 N1mean[i] = Ameas[i-1];
303 }
304 else{
305 N1mean[i] = ( Ameas[i-1] + Ameas[i+1] ) / 2.;
306 }
307*/
308 N1mean[i] = ( Ameas[i-1] + Ameas[i+1] ) / 2.;
309 }
310} // end of computeN1mean computation
311
312float getValue( int slice, int pixel ){
313
314 const float dconv = 2000/4096.0;
315
316 float vraw, vcal;
317
318 unsigned int pixel_pt;
319 unsigned int slice_pt;
320 unsigned int cal_pt;
321 unsigned int drs_cal_offset;
322
323 // printf("pixel = %d, slice = %d\n", slice, pixel);
324
325 pixel_pt = pixel * data_roi;
326 slice_pt = pixel_pt + slice;
327 drs_cal_offset = ( slice + data_offset[ pixel ] )%data_roi;
328 cal_pt = pixel_pt + drs_cal_offset;
329
330 vraw = data[ slice_pt ] * dconv;
331 vcal = ( vraw - drs_basemean[ cal_pt ] - drs_triggeroffsetmean[ slice_pt ] ) / drs_gainmean[ cal_pt ]*1907.35;
332
333 return( vcal );
334}
335
336void BookHistos( int Samples ){
337// booking and parameter settings for all histos
338
339 // histograms for baseline extraction
340 char hName[500];
341 char hTitle[500];
342
343 TH1F *h;
344
345 printf("inside BookHistos\n");
346
347 for( int i = 0; i < NPIX; i++ ) {
348
349 // printf("call sprintf [%d]\n", i );
350 sprintf(&hTitle[0],"all events all slices of pixel %d", i);
351 sprintf(&hName[0],"base%d", i);
352 // printf("call sprintf [%d] done\n", i );
353
354 h = new TH1F( hName, hTitle, 400, -99.5 ,100.5 );
355
356 // printf("histo booked\n");
357 h->GetXaxis()->SetTitle( "Sample value (mV)" );
358 h->GetYaxis()->SetTitle( "Entries / 0.5 mV" );
359 // printf("histo title set\n");
360 hListBaseline.Add( h );
361 // printf("histo added to List\n");
362 hBaseline[i] = h;
363 // printf("histo assigned to array\n");
364 }
365
366 printf("made HBaseline * 1440\n");
367
368 hMeanBsl = new TH1F("histo_mean","Value of maximal probability",400,-99.5,100.5);
369 hMeanBsl->GetXaxis()->SetTitle( "max value (mV)" );
370 hMeanBsl->GetYaxis()->SetTitle( "Entries / 0.5 mV" );
371 hList.Add( hMeanBsl );
372
373 hpltMeanBsl = new TH1F("hplt_mean","Value of maximal probability",1440,-0.5,1439.5);
374 hpltMeanBsl->GetXaxis()->SetTitle( "pixel" );
375 hpltMeanBsl->GetYaxis()->SetTitle( "max value in mV" );
376 hList.Add( hpltMeanBsl );
377
378 hRmsBsl = new TH1F("histo_rms","RMS in mV",2000,-99.5,100.5);
379 hRmsBsl->GetXaxis()->SetTitle( "RMS (mV)" );
380 hRmsBsl->GetYaxis()->SetTitle( "Entries / 0.5 mV" );
381 hList.Add( hRmsBsl );
382
383 hpltRmsBsl = new TH1F("hplt_rms","Value of maximal probability",1440,-0.5,1439.5);
384 hpltRmsBsl->GetXaxis()->SetTitle( "pixel" );
385 hpltRmsBsl->GetYaxis()->SetTitle( "RMS in mV" );
386 hList.Add( hpltRmsBsl );
387}
388
389
390void SaveHistograms( char * loc_fname ){
391
392 TString fName; // name of the histogram file
393
394 /* create the filename for the histogram file */
395 fName = loc_fname; // use the name of the tree file
396 fName.Remove(fName.Length() - 5); // remove the extension .root
397 fName += "_histo.root"; // add the new extension
398
399 TFile tf( fName, "RECREATE"); // create the histogram file (replace if already existing)
400
401 hList.Write(); // write the major histograms into the top level directory
402 tf.mkdir("BaselineHisto"); tf.cd("BaselineHisto"); // go to new subdirectory
403 hListBaseline.Write(); // write histos into subdirectory
404
405 tf.Close(); // close the file
406
407} // end of function: void ana::SaveHistograms( void )
408
409int FOpenDataFile(fits &datafile){
410
411/* cout << "-------------------- Data Header -------------------" << endl;
412 datafile.PrintKeys();
413 cout << "------------------- Data Columns -------------------" << endl;
414 datafile.PrintColumns();
415 */
416
417 //Get the size of the data column
418 data_roi = datafile.GetUInt("NROI"); // Value from header
419 data_px = datafile.GetUInt("NPIX");
420 data_n = datafile.GetN("Data"); //Size of column "Data" = #Pixel x ROI
421
422 //Set the sizes of the data vectors
423 data.resize(data_n,0);
424 data_offset.resize(data_px,0);
425
426 //Link the data to variables
427 datafile.SetRefAddress("EventNum", data_num);
428 datafile.SetVecAddress("Data", data);
429 datafile.SetVecAddress("StartCellData", data_offset);
430 datafile.GetRow(0);
431
432 cout << "Opening data file successful..." << endl;
433
434 return datafile.GetNumRows() ;
435}
436
Note: See TracBrowser for help on using the repository browser.