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

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