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

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