source: fact/tools/rootmacros/fana2.C@ 13560

Last change on this file since 13560 was 12604, checked in by lusterma, 14 years ago
added fana2.C, a working version of fana.C
File size: 8.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#define HAVE_ZLIB
15#include "fits.h"
16
17#include "openFits.h"
18#include "openFits.c"
19
20#include "DrsCalibration.h"
21#include "DrsCalibration.C"
22
23#include "SpikeRemoval.h"
24#include "SpikeRemoval.C"
25
26#define NPIX 1440
27#define NCELL 1024
28
29// data access and treatment
30#define FAD_MAX_SAMPLES 1024
31
32int NEvents;
33vector<int16_t> Data; // vector, which will be filled with raw data
34vector<int16_t> StartCells; // vector, which will be filled with DRS start positions
35unsigned int EventID; // index of the current event
36UInt_t RegionOfInterest; // Width of the Region, read out of the DRS
37UInt_t NumberOfPixels; // Total number of pixel, read out of the camera
38size_t PXLxROI; // Size of column "Data" = #Pixel x ROI
39
40int NBSLeve = 1000;
41
42size_t TriggerOffsetROI, RC;
43vector<float> Offset, Gain, TriggerOffset;
44
45
46vector<float> Ameas(FAD_MAX_SAMPLES); // copy of the data (measured amplitude
47vector<float> N1mean(FAD_MAX_SAMPLES); // mean of the +1 -1 ch neighbors
48vector<float> N2mean(FAD_MAX_SAMPLES); // mean of the +2 -2 ch neighbors
49vector<float> SpikeEst(FAD_MAX_SAMPLES); // corrected Values
50vector<float> Vcorr(FAD_MAX_SAMPLES); // corrected Values
51vector<float> Vdiff(FAD_MAX_SAMPLES); // numerical derivative
52
53vector<float> Vslide(FAD_MAX_SAMPLES); // sliding average result
54vector<float> Vcfd(FAD_MAX_SAMPLES); // CDF result
55vector<float> Vcfd2(FAD_MAX_SAMPLES); // CDF result + 2nd sliding average
56
57#include "factfir.C"
58
59void computeSpikeEstimator( int );
60
61// histograms
62const int Ntypes = 7;
63const unsigned int // arranged by Dominik
64 tAmeas = 0,
65 tSpikeEst = 1,
66 tVcorr = 2,
67 tVtest = 3,
68 tVslide = 4,
69 tVcfd = 5,
70 tVcfd2 = 6;
71
72TH1F* h;
73//TH2F* 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 fana2(
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 = NULL;
116 NEvents = openDataFits(
117 datafilename,
118 &datafile,
119 Data,
120 StartCells,
121 EventID,
122 RegionOfInterest,
123 NumberOfPixels,
124 PXLxROI );
125-
126 printf( "number of events in file: %d\n", NEvents );
127 if (NEvents == 0){
128 cout << "return code of openDataFits:" << datafilename<< endl;
129 cout << "is zero -> aborting." << endl;
130 return 1;
131 }
132
133 // compare the number of events in the data file with the nevents the
134 // the user would like to read. nevents = -1 means: read all
135 if ( nevents == -1 || nevents > NEvents ) nevents = NEvents;
136
137
138//-------------------------------------------
139//Get the DRS calibration
140//-------------------------------------------
141
142 RC = openCalibFits( drsfilename, Offset, Gain, TriggerOffset, TriggerOffsetROI, 3);
143 if (RC == 0){
144 cout << "return code of openCalibFits:" << drsfilename << endl;
145 cout << "is zero -> aborting." << endl;
146 return 1;
147 }
148
149// Book the histograms
150 BookHistos( RegionOfInterest );
151
152// Loop over events
153 cout << "--------------------- Data --------------------" << endl;
154
155 // float value;
156
157 // TH1F * sp = new TH1F("spektrum", "test of Stepktrum", 256, -0.5, 63.5);
158
159 size_t calib_RC = 1;
160 for ( int ev = firstevent; ev < firstevent + nevents; ev++) {
161
162 datafile->GetRow( ev );
163 if (ev % 50 ==0){
164 cout << "Event number: " << EventID << endl;
165 }
166
167 // get the data of this pixel from the Data vector
168 // apply the Drs Calibration and cut off 12 slices at the beginning
169 // and at the end.
170 calib_RC = applyDrsCalibration( Ameas, pixelnr,12,12, Offset, Gain, TriggerOffset,
171 RegionOfInterest, Data, StartCells);
172 if (calib_RC == 0){
173 break;
174 }
175
176 computeSpikeEstimator( RegionOfInterest );
177 removeSpikes( Ameas, Vcorr );
178 sliding_avg( Vcorr, Vslide, 8 );
179
180 // for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
181 // hPixelCellData.Fill(sl, Vcorr[ sl ]);
182 //}
183
184 // filter Vcorr with sliding average using FIR filter function
185 factfir(b_slide , a_slide, k_slide, Vcorr, Vslide);
186
187 // filter Vslide with CFD using FIR filter function
188 factfir(b_cfd , a_cfd, k_cfd, Vslide, Vcfd);
189 // filter Vcfd with sliding average using FIR filter function
190 factfir(b_slide , a_slide, k_slide, Vcfd, Vcfd2);
191
192 if ( spikeDebug ){
193 for ( unsigned int sl = 0; sl < RegionOfInterest; sl++ ){
194
195 h[tAmeas].SetBinContent( sl, Ameas[ sl ] );
196 h[tVcorr].SetBinContent( sl, Vcorr[ sl ] );
197
198 h[tVslide].SetBinContent( sl, Vslide[ sl ] );
199 h[tVcfd].SetBinContent( sl, Vcfd[ sl ] );
200 h[tVcfd2].SetBinContent( sl, Vcfd2[ sl ] );
201 }
202 }
203
204 if ( spikeDebug ){
205 CW->cd( tAmeas + 1);
206 gPad->SetGrid();
207 h[tAmeas].Draw();
208
209 CW->cd( tSpikeEst + 1);
210 gPad->SetGrid();
211 h[tSpikeEst].Draw();
212
213 CW->cd( tVcorr + 1);
214 gPad->SetGrid();
215 h[tVcorr].Draw();
216
217 cFilter->cd( Ntypes - tVslide );
218 cFilter->cd(1);
219 gPad->SetGrid();
220 h[tVslide].Draw();
221 cFilter->cd( Ntypes - tVcfd );
222 cFilter->cd(2);
223 gPad->SetGrid();
224 h[tVcfd].Draw();
225 TLine zeroline(0, 0, 1024, 0);
226 zeroline.SetLineColor(kBlue);
227 zeroline.Draw();
228
229 cFilter->cd( Ntypes - tVcfd2 );
230 cFilter->cd(3);
231 gPad->SetGrid();
232 h[tVcfd2].Draw();
233
234 zeroline.Draw();
235
236 CW->Update();
237 cFilter->Update();
238
239 //Process gui events asynchronously during input
240 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
241 timer.TurnOn();
242 TString input = Getline("Type 'q' to exit, <return> to go on: ");
243 timer.TurnOff();
244 if (input=="q\n") break;
245 }
246
247
248 }
249
250 return( 0 );
251}
252
253void computeSpikeEstimator( int Samples ){
254
255// compute the mean of the left and right neighbors of a channel
256
257 for( int i = 1; i < Samples-1; i++){
258
259 N1mean[ i ] = ( Ameas[i-1] + Ameas[i+1] ) / 2.;
260 SpikeEst[ i ] = Ameas[ i ] - N1mean[ i ];
261 h[tSpikeEst].SetBinContent( i, SpikeEst[ i ] );
262 }
263} // end of computeN1mean computation
264
265
266void BookHistos( int Samples ){
267// booking and parameter settings for all histos
268
269 h = new TH1F[ Ntypes ];
270
271 for ( int type = 0; type < Ntypes; type++){
272
273 h[ type ].SetBins(Samples, 0, Samples);
274 h[ type ].SetLineColor(1);
275 h[ type ].SetLineWidth(2);
276
277 // set X axis paras
278 h[ type ].GetXaxis()->SetLabelSize(0.1);
279 h[ type ].GetXaxis()->SetTitleSize(0.1);
280 h[ type ].GetXaxis()->SetTitleOffset(1.2);
281 h[ type ].GetXaxis()->SetTitle(Form("Time slice (%.1f ns/slice)", 1./2.));
282
283 // set Y axis paras
284 h[ type ].GetYaxis()->SetLabelSize(0.1);
285 h[ type ].GetYaxis()->SetTitleSize(0.1);
286 h[ type ].GetYaxis()->SetTitleOffset(0.3);
287 h[ type ].GetYaxis()->SetTitle("Amplitude (a.u.)");
288 }
289 CW = new TCanvas("CW","DRS Waveform",10,10,800,600);
290 CW->Divide(1, 3);
291 cFilter = new TCanvas("cFilter","filtered DRS Waveforms",10,10,800,600);
292 cFilter->Divide(1, 3);
293
294 // hStartCell = new TH2F("StartCell", "StartCell", 1440, 0., 1440., 1024, 0., 1024);
295
296}
Note: See TracBrowser for help on using the repository browser.