source: fact/tools/rootmacros/zippedfitstest.C@ 12196

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