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

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