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

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