source: fact/tools/rootmacros/fana.C@ 18481

Last change on this file since 18481 was 12305, checked in by neise, 13 years ago
adjusted use of zerosearch() method ... this method doesn't return ints anymore, but struct Region see Region.h
File size: 11.8 KB
Line 
1
2#include <TROOT.h>
3#include <TCanvas.h>
4#include <TProfile.h>
5#include <TTimer.h>
6#include <TH1F.h>
7#include <TH2F.h>
8#include <Getline.h>
9#include <TLine.h>
10#include <TMath.h>
11
12#include <stdint.h>
13#include <cstdio>
14
15#define HAVE_ZLIB
16#include "fits.h"
17//#include "TPKplotevent.c"
18//#include "FOpenDataFile.c"
19#include "FOpenCalibFile.c"
20
21#include "Region.h"
22#include "zerosearch.h"
23#include "zerosearch.C"
24
25#define NPIX 1440
26#define NCELL 1024
27
28// data access and treatment
29#define FAD_MAX_SAMPLES 1024
30
31vector<int16_t> data;
32vector<int16_t> data_offset;
33
34unsigned int data_num;
35
36size_t data_n;
37UInt_t data_px;
38UInt_t data_roi;
39int NEvents;
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("PixelPedestal", "PixelPedestal", NCELL, 0., NCELL, 200, -50., 150.);
80
81void BookHistos( int );
82
83// Create a canvas
84TCanvas* CW;
85TCanvas* cFilter;
86
87int spikeDebug = 1;
88
89
90int fana(
91 char *datafilename = "../raw/20110916_025.fits",
92 const char *drsfilename = "../raw/20110916_024.drs.fits",
93 int pixelnr = 0,
94 int firstevent = 0,
95 int nevents = -1 ){
96// read and analyze FACT raw data
97
98// sliding window filter settings
99 int k_slide = 16;
100 vector<double> a_slide(k_slide, 1);
101 double b_slide = k_slide;
102
103// CFD filter settings
104 int k_cfd = 10;
105 vector<double> a_cfd(k_cfd, 0);
106 double b_cfd = 1.;
107 a_cfd[0]=0.75;
108 a_cfd[k_cfd-1]=-1.;
109
110// 2nd slinding window filter
111 //int ks2 = 16;
112 //vector<double> as2(ks2, 1);
113 //double bs2 = ks2;
114 gROOT->SetStyle("Plain");
115
116//-------------------------------------------
117// Open the file
118//-------------------------------------------
119 fits datafile( datafilename );
120 if (!datafile) {
121 printf( "Could not open the file: %s\n", datafilename );
122 return 1;
123 }
124
125 // access data
126 NEvents = FOpenDataFile( datafile );
127
128 printf("number of events in file: %d\n", NEvents);
129
130 // compare the number of events in the data file with the nevents the
131 // the user would like to read. nevents = -1 means: read all
132 if ( nevents == -1 || nevents > NEvents ) nevents = NEvents;
133
134
135//-------------------------------------------
136//Get the DRS calibration
137//-------------------------------------------
138
139 FOpenCalibFile(drsfilename, drs_basemean, drs_gainmean, drs_triggeroffsetmean, drs_n);
140
141//-------------------------------------------
142//Check the sizes of the data columns
143//-------------------------------------------
144 if(drs_n!=data_n)
145 {
146 cout << "Data and DRS file incompatible (Px*ROI disagree)" << endl;
147 return 1;
148 }
149// Book the histograms
150 BookHistos( data_roi );
151
152// Loop over events
153 cout << "--------------------- Data --------------------" << endl;
154
155 float value;
156 TH1F * sp = new TH1F("spektrum", "test of Stepktrum", 256, -0.5, 63.5);
157
158 for ( int ev = firstevent; ev < firstevent + nevents; ev++) {
159
160 datafile.GetRow( ev );
161
162 if (ev % 50 ==0){
163 cout << "Event number: " << data_num << endl;
164 }
165
166 for ( int pix = 0; pix < 1440; pix++ ){
167 hStartCell->Fill( pix, data_offset[pix] );
168 }
169
170 for ( unsigned int sl = 0; sl < data_roi; sl++){
171 value = getValue(sl, pixelnr);
172 //printf("value = %f\n", value);
173 Ameas[ sl ] = value;
174 h[tAmeas].SetBinContent(sl, value);
175
176 }
177
178 computeN1mean( data_roi );
179 removeSpikes( data_roi );
180
181 for ( unsigned int sl = 0; sl < data_roi; sl++){
182 hPixelCellData.Fill(sl, Vcorr[ sl ]);
183 }
184
185 // filter Vcorr with sliding average using FIR filter function
186 factfir(b_slide , a_slide, k_slide, Vcorr, Vslide);
187
188 // filter Vslide with CFD using FIR filter function
189 factfir(b_cfd , a_cfd, k_cfd, Vslide, Vcfd);
190 // filter Vcfd with sliding average using FIR filter function
191 factfir(b_slide , a_slide, k_slide, Vcfd, Vcfd2);
192
193 if ( spikeDebug ){
194 for ( unsigned int sl = 0; sl < data_roi; sl++){
195 h[tVslide].SetBinContent( sl, Vslide[sl] );
196 h[tVcfd].SetBinContent( sl, Vcfd[sl] );
197 h[tVcfd2].SetBinContent( sl, Vcfd2[sl] );
198 }
199 }
200
201 vector<Region> * zeros = zerosearch( Vcfd2, -1, 1, 0);
202 ShiftRegionBy(*zeros, -k_slide/2);
203 EnlargeRegion(*zeros, 10, 20);
204 findAbsMaxInRegions( *zeros, Vslide);
205 if (zeros->size() == 0 ){
206 continue;
207 }
208 // check value of Vside at zero position
209 for (unsigned int i=0; i<zeros->size(); i++){
210 cout << zeros->at(i).maxPos << ":\t" << zeros->at(i).maxVal <<endl;
211 sp->Fill( zeros->at(i).maxVal );
212 }
213
214 if ( spikeDebug ){
215
216 CW->cd( tAmeas + 1);
217 gPad->SetGrid();
218 h[tAmeas].Draw();
219
220 CW->cd( tN1mean + 1);
221 gPad->SetGrid();
222 h[tN1mean].Draw();
223
224 CW->cd( tVcorr + 1);
225 gPad->SetGrid();
226 h[tVcorr].Draw();
227
228 // CW->cd( tVtest + 1);
229 // gPad->SetGrid();
230 // h[tVtest].Draw();
231
232 cFilter->cd( Ntypes - tVslide );
233 cFilter->cd(1);
234 gPad->SetGrid();
235 h[tVslide].Draw();
236
237 cFilter->cd( Ntypes - tVcfd );
238 cFilter->cd(2);
239 gPad->SetGrid();
240 h[tVcfd].Draw();
241
242 TLine zeroline(0, 0, 1024, 0);
243 zeroline.SetLineColor(kBlue);
244 zeroline.Draw();
245
246 cFilter->cd( Ntypes - tVcfd2 );
247 cFilter->cd(3);
248 gPad->SetGrid();
249 h[tVcfd2].Draw();
250
251 zeroline.Draw();
252
253 CW->Update();
254 cFilter->Update();
255
256 //Process gui events asynchronously during input
257 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
258 timer.TurnOn();
259 TString input = Getline("Type 'q' to exit, <return> to go on: ");
260 timer.TurnOff();
261 if (input=="q\n") break;
262 }
263
264
265 }
266
267 TCanvas * cSpektrum = new TCanvas();
268 cSpektrum->cd();
269 sp->Draw();
270 TCanvas * cStartCell = new TCanvas();
271 cStartCell->cd();
272 hStartCell->Draw();
273 hPixelCellData.Draw();
274
275 delete cStartCell;
276
277 return( 0 );
278}
279
280void removeSpikes(int Samples){
281
282 const float fract = 0.8;
283 float x, xp, xpp, x3p;
284
285 // assume that there are no spikes
286 for ( int i = 0; i < Samples; i++) Vcorr[i] = Ameas[i];
287
288// find the spike and replace it by mean value of neighbors
289 for ( int i = 0; i < Samples; i++) {
290
291 // printf("Vcorr[%d] = %f, Ameas[%d] = %f\n", i, Vcorr[ i ], i, Ameas[ i ] );
292
293 x = Ameas[i] - N1mean[i];
294
295 if ( x < -5. ){ // a spike candidate
296 // check consistency with a single channel spike
297 xp = Ameas[i+1] - N1mean[i+1];
298 xpp = Ameas[i+2] - N1mean[i+2];
299 x3p = Ameas[i+3] - N1mean[i+3];
300
301 // printf("candidates x[%d] = %f; xp = %f; xpp = %f, x3p = %f\n", i, x, xp, xpp, x3p);
302
303 if ( Ameas[i+2] - ( Ameas[i] + Ameas[i+3] )/2. > 10. ){
304 // printf("double spike candidate\n");
305 Vcorr[i+1] = ( Ameas[i] + Ameas[i+3] )/2.;
306 Vcorr[i+2] = ( Ameas[i] + Ameas[i+3] )/2.;
307 // printf("Vcorr[%d] = %f %f %f %f\n", i, Vcorr[i], Vcorr[i+1], Vcorr[i+2], Vcorr[ i+3 ]);
308 // printf("Ameas[%d] = %f %f %f %f\n", i, Ameas[ i ], Ameas[ i+1 ], Ameas[ i+2 ], Ameas[i+3]);
309 i = i + 3;
310 }
311 else{
312
313 if ( ( xp > -2.*x*fract ) && ( xpp < -10. ) ){
314 Vcorr[i+1] = N1mean[i+1];
315 // printf("Vcorr[%d] = %f %f %f\n", i, Vcorr[i], Vcorr[i+1], Vcorr[i+2]);
316 // N1mean[i+1] = (Ameas[i] - Ameas[i+2] / 2.);
317 N1mean[i+2] = (Ameas[i+1] - Ameas[i+3] / 2.);
318 i = i + 2;//do not care about the next sample it was the spike
319 }
320 // treatment for the end of the pipeline must be added !!!
321 }
322 }
323 else{
324 // do nothing
325 }
326 } // end of spike search and correction
327 for ( int i = 0; i < Samples; i++ ) h[ tVcorr ].SetBinContent( i, Vcorr[i] );
328}
329
330void computeN1mean( int Samples ){
331
332// compute the mean of the left and right neighbors of a channel
333
334 for( int i = 0; i < Samples; i++){
335 if (i == 0){ // use right sample es mean
336 N1mean[i] = Ameas[i+1];
337 }
338 else if ( i == Samples-1 ){ //use left sample as mean
339 N1mean[i] = Ameas[i-1];
340 }
341 else{
342 N1mean[i] = ( Ameas[i-1] + Ameas[i+1] ) / 2.;
343 }
344 h[tN1mean].SetBinContent(i, Ameas[i] - N1mean[i]);
345 }
346} // end of computeN1mean computation
347
348float getValue( int slice, int pixel ){
349
350 const float dconv = 2000/4096.0;
351
352 float vraw, vcal;
353
354 unsigned int pixel_pt;
355 unsigned int slice_pt;
356 unsigned int cal_pt;
357 unsigned int drs_cal_offset;
358
359 // printf("pixel = %d, slice = %d\n", slice, pixel);
360
361 pixel_pt = pixel * data_roi;
362 slice_pt = pixel_pt + slice;
363 drs_cal_offset = ( slice + data_offset[ pixel ] )%data_roi;
364 cal_pt = pixel_pt + drs_cal_offset;
365
366 vraw = data[ slice_pt ] * dconv;
367 vcal = ( vraw - drs_basemean[ cal_pt ] - drs_triggeroffsetmean[ slice_pt ] ) / drs_gainmean[ cal_pt ]*1907.35;
368
369 return( vcal );
370}
371
372void BookHistos( int Samples ){
373// booking and parameter settings for all histos
374
375 h = new TH1F[ Ntypes ];
376
377 for ( int type = 0; type < Ntypes; type++){
378
379 h[ type ].SetBins(Samples, 0, Samples);
380 h[ type ].SetLineColor(1);
381 h[ type ].SetLineWidth(2);
382
383 // set X axis paras
384 h[ type ].GetXaxis()->SetLabelSize(0.1);
385 h[ type ].GetXaxis()->SetTitleSize(0.1);
386 h[ type ].GetXaxis()->SetTitleOffset(1.2);
387 h[ type ].GetXaxis()->SetTitle(Form("Time slice (%.1f ns/slice)", 1./2.));
388
389 // set Y axis paras
390 h[ type ].GetYaxis()->SetLabelSize(0.1);
391 h[ type ].GetYaxis()->SetTitleSize(0.1);
392 h[ type ].GetYaxis()->SetTitleOffset(0.3);
393 h[ type ].GetYaxis()->SetTitle("Amplitude (a.u.)");
394 }
395 CW = new TCanvas("CW","DRS Waveform",10,10,800,600);
396 CW->Divide(1, 3);
397 cFilter = new TCanvas("cFilter","filtered DRS Waveforms",10,10,800,600);
398 cFilter->Divide(1, 3);
399
400 hStartCell = new TH2F("StartCell", "StartCell", 1440, 0., 1440., 1024, 0., 1024);
401
402}
403int FOpenDataFile(fits &datafile){
404
405/* cout << "-------------------- Data Header -------------------" << endl;
406 datafile.PrintKeys();
407 cout << "------------------- Data Columns -------------------" << endl;
408 datafile.PrintColumns();
409 */
410
411 //Get the size of the data column
412 data_roi = datafile.GetUInt("NROI"); // Value from header
413 data_px = datafile.GetUInt("NPIX");
414 data_n = datafile.GetN("Data"); //Size of column "Data" = #Pixel x ROI
415
416 //Set the sizes of the data vectors
417 data.resize(data_n,0);
418 data_offset.resize(data_px,0);
419
420 //Link the data to variables
421 datafile.SetRefAddress("EventNum", data_num);
422 datafile.SetVecAddress("Data", data);
423 datafile.SetVecAddress("StartCellData", data_offset);
424 datafile.GetRow(0);
425
426 cout << "Opening data file successful..." << endl;
427
428 return datafile.GetNumRows() ;
429}
430
Note: See TracBrowser for help on using the repository browser.