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

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