source: fact/tools/rootmacros/fpeak_cdf.C@ 12383

Last change on this file since 12383 was 12379, checked in by neise, 14 years ago
moved spike removal out of this file
  • Property svn:executable set to *
File size: 9.1 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 <TBox.h>
10#include <TMath.h>
11#include <TFile.h>
12#include <TStyle.h>
13
14#include <stdio.h>
15#include <stdint.h>
16#include <cstdio>
17#include <deque>
18
19#define NPIX 1440
20#define NCELL 1024
21#define FAD_MAX_SAMPLES 1024
22
23#define HAVE_ZLIB
24#include "fits.h"
25#include "FOpenCalibFile.c"
26
27#include "discriminator.h"
28#include "discriminator.C"
29#include "zerosearch.h"
30#include "zerosearch.C"
31#include "factfir.C"
32
33#include "FOpenDataFile.h"
34#include "FOpenDataFile.c"
35
36#include "DrsCalibration.C"
37#include "DrsCalibration.h"
38
39#include "SpikeRemoval.h"
40#include "SpikeRemoval.C"
41
42bool breakout=false;
43
44int NEvents;
45vector<int16_t> AllPixelDataVector;
46vector<int16_t> StartCellVector;
47unsigned int CurrentEventID;
48size_t PXLxROI;
49UInt_t RegionOfInterest;
50UInt_t NumberOfPixels;
51
52size_t drs_n;
53vector<float> drs_basemean;
54vector<float> drs_gainmean;
55vector<float> drs_triggeroffsetmean;
56
57vector<float> Ameas(FAD_MAX_SAMPLES); // copy of the data (measured amplitude
58vector<float> Vcorr(FAD_MAX_SAMPLES); // corrected Values
59vector<float> Vslide(FAD_MAX_SAMPLES); // sliding average result
60vector<float> Vcfd(FAD_MAX_SAMPLES); // CDF result
61vector<float> Vcfd2(FAD_MAX_SAMPLES); // CDF result + 2nd sliding average
62
63// histograms
64const int NumberOfDebugHistoTypes = 7;
65const unsigned int
66 Ameas_ = 0,
67 N1mean_ = 1,
68 Vcorr_ = 2,
69 Vtest_ = 3,
70 Vslide_ = 4,
71 Vcfd_ = 5,
72 Vcfd2_ = 6;
73
74TH1F* h;
75TH1F *debugHistos;
76TH2F* hStartCell; // id of the DRS physical pipeline cell where readout starts
77 // x = pixel id, y = DRS cell id
78
79TH2F * hAmplSpek_cfd;
80TObjArray hList;
81
82void BookHistos( );
83void SaveHistograms(const char * );
84
85///////////////////////////////////////////////////////////////////////////////
86///////////////////////////////////////////////////////////////////////////////
87///////////////////////////////////////////////////////////////////////////////
88// read FACT raw data
89// * remove spikes
90// * calculate baseline
91// * find peaks (CFD and normal discriminator)
92// * compute pulse height and pulse integral spektrum of the peaks
93int fpeak(
94 char *datafilename = "data/20111016_013.fits.gz",
95 const char *drsfilename = "../../20111016_011.drs.fits.gz",
96 const char *OutRootFileName = "../analysis/fpeak_cdf.Coutput.root",
97 int firstevent = 0,
98 int nevents = -1,
99 int firstpixel = 0,
100 int npixel = -1,
101 bool spikeDebug = false,
102 int avg1 = 14,
103 int avg2 = 8,
104 int verbosityLevel = 1, // different verbosity levels can be implemented here
105 bool ProduceGraphic = true
106 )
107{
108gStyle->SetPalette(1,0);
109gROOT->SetStyle("Plain");
110
111 TCanvas *cFiltered = NULL;
112 if (spikeDebug){
113 cFiltered = new TCanvas("cFiltered","filtered DRS Waveforms",410,410,400,400);
114 cFiltered->Divide(1, 3);
115 }
116
117
118
119// CFD filter settings
120 int k_cfd = 10;
121 vector<double> a_cfd(k_cfd, 0);
122 double b_cfd = 1.;
123 a_cfd[0]=-0.75;
124 a_cfd[k_cfd-1]=1.;
125
126 fits * datafile;
127 // Opens the raw data file and 'binds' the variables given as
128 // Parameters to the data file. So they are filled with
129 // raw data as soon as datafile->GetRow(int) is called.
130 NEvents = OpenDataFile( datafilename, &datafile,
131 AllPixelDataVector, StartCellVector, CurrentEventID,
132 RegionOfInterest, NumberOfPixels, PXLxROI, verbosityLevel);
133 if (NEvents == 0){
134 cout << "return code of FOpenDataFile:" << datafilename<< endl;
135 cout << "is zero -> aborting." << endl;
136 return 1;
137 }
138
139 if (verbosityLevel > 0)
140 cout <<"number of events in file: "<< NEvents << "\t";
141 if ( nevents == -1 || nevents > NEvents ) nevents = NEvents; // -1 means all!
142 if (verbosityLevel > 0)
143 cout <<"of, which "<< nevents << "will be processed"<< endl;
144
145 if (verbosityLevel > 0)
146 cout <<"Totel # of Pixel: "<< NumberOfPixels << "\t";
147 if ( npixel == -1 || npixel > (int)NumberOfPixels ) npixel = (int)NumberOfPixels; // -1 means all!
148 if (verbosityLevel > 0)
149 cout <<"of, which "<< nevents << "will be processed"<< endl;
150
151
152
153 //Get the DRS calibration
154 FOpenCalibFile( drsfilename,
155 drs_basemean,
156 drs_gainmean,
157 drs_triggeroffsetmean,
158 drs_n);
159
160 BookHistos( );
161
162 for ( int ev = firstevent; ev < firstevent + nevents; ev++) {
163 // Get an Event --> consists of 1440 Pixel ...erm....data
164 datafile->GetRow( ev );
165
166 for ( int pix = firstpixel; pix < firstpixel+npixel; pix++ ){
167 if (verbosityLevel > 0){
168 if (pix == firstpixel){
169 cout << "Processing Event: " << CurrentEventID << "/" << nevents << endl;
170 }
171 }
172
173 applyDrsCalibration( Ameas,pix,
174 drs_basemean, drs_gainmean, drs_triggeroffsetmean,
175 RegionOfInterest, AllPixelDataVector, StartCellVector);
176
177 // finds spikes in the raw data, and interpolates the value
178 // spikes are: 1 or 2 slice wide, positive non physical artifacts
179 removeSpikes (Ameas, Vcorr);
180
181 // filter Vcorr with sliding average using FIR filter function
182 sliding_avg(Vcorr, Vslide, avg1);
183 // filter Vslide with CFD using FIR filter function
184 factfir(b_cfd , a_cfd, k_cfd, Vslide, Vcfd);
185 // filter Vcfd with sliding average using FIR filter function
186 sliding_avg(Vcfd, Vcfd2, avg2);
187
188
189 // peaks in Ameas[] are found by searching for zero crossings
190 // in Vcfd2
191 // first Argument 1 means ... *rising* edge
192 // second Argument 1 means ... search with stepsize 1 ... 8 is okay as well
193 vector<Region> * zXings = zerosearch( Vcfd2 , 1 , 8);
194 // zXings means "zero cross ings"
195 EnlargeRegion( *zXings, 10, 10);
196 findAbsMaxInRegions( *zXings, Vslide);
197 removeMaximaBelow( *zXings, 3.0, 0);
198 removeRegionOnFallingEdge( *zXings, 100);
199
200 // fill maxima in Histogram
201 if (zXings->size() != 0 ){
202 for (unsigned int i=0; i<zXings->size(); i++){
203 if (verbosityLevel > 1){
204 cout << zXings->at(i).maxPos << ":\t"
205 << zXings->at(i).maxVal <<endl;
206 }
207 hAmplSpek_cfd->Fill(pix, zXings->at(i).maxVal);
208 }
209 }
210
211 if ( spikeDebug ){
212 for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
213 debugHistos[Ameas_].SetBinContent(sl, Ameas[sl]);
214 debugHistos[Vslide_].SetBinContent( sl, Vslide[sl] );
215 debugHistos[Vcfd_].SetBinContent( sl, Vcfd[sl] );
216 debugHistos[Vcfd2_].SetBinContent( sl, Vcfd2[sl] );
217 }
218
219 cFiltered->cd(1);
220 gPad->SetGrid();
221 debugHistos[Vcorr_].Draw();
222
223 cFiltered->cd(2);
224 gPad->SetGrid();
225 debugHistos[Vslide_].Draw();
226
227 TBox *OneBox;
228 vector<TBox*> MyBoxes;
229 for (unsigned int i=0; i<zXings->size(); i++){
230 OneBox = new TBox(
231 zXings->at(i).maxPos -10 ,
232 zXings->at(i).maxVal -0.5,
233 zXings->at(i).maxPos +10 ,
234 zXings->at(i).maxVal +0.5);
235 OneBox->SetLineColor(kBlue);
236 OneBox->SetLineWidth(1);
237 OneBox->SetFillStyle(0);
238 OneBox->SetFillColor(kRed);
239 MyBoxes.push_back(OneBox);
240 OneBox->Draw();
241 }
242
243 cFiltered->cd(3);
244 gPad->SetGrid();
245 debugHistos[Vcfd2_].Draw();
246 TLine *zeroline = new TLine(0, 0, 1024, 0);
247 zeroline->SetLineColor(kBlue);
248 zeroline->Draw();
249
250 cFiltered->Update();
251 //Process gui events asynchronously during input
252 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
253 timer.TurnOn();
254 TString input = Getline("Type 'q' to exit, <return> to go on: ");
255 timer.TurnOff();
256 if (input=="q\n") {
257 breakout=true;
258 }
259
260 //TODO!!!!!!!!!
261 // do some Garbage collection ...
262 }// end of if(spikeDebug)
263
264 delete zXings;
265 if (breakout)
266 break;
267
268 } // end of loop over pixels
269
270 if (breakout)
271 break;
272 } // end of loop over pixels
273if (ProduceGraphic){
274 TCanvas * cSpektrum;
275 cSpektrum = new TCanvas("cSpektrum","Amplitude Spektra of different discriminators",10,10,400,400);
276 cSpektrum->Divide(1,1);
277 cSpektrum->cd(1);
278 hAmplSpek_cfd->Draw("COLZ");
279 cSpektrum->Modified();
280 cSpektrum->Update();
281}
282 SaveHistograms( OutRootFileName );
283 return( 0 );
284}
285
286// booking and parameter settings for all histos
287void BookHistos( ){
288 // histograms for baseline extraction
289
290 TString name,title;
291 title = "all events all slices of pixel ";
292 name = "base";
293
294 hAmplSpek_cfd = new TH2F("hAmplSpek_cfd","amplitude spektrum - CFD",1440,-0.5,1439.5, 256, -27.5, 100.5);
295 hAmplSpek_cfd->GetXaxis()->SetTitle( "pixel" );
296 hAmplSpek_cfd->GetYaxis()->SetTitle( "amplitude in mV" );
297 hList.Add( hAmplSpek_cfd );
298
299
300 debugHistos = new TH1F[ NumberOfDebugHistoTypes ];
301 for ( int type = 0; type < NumberOfDebugHistoTypes; type++){
302 debugHistos[ type ].SetBins(1024, 0, 1024);
303 debugHistos[ type ].SetLineColor(1);
304 debugHistos[ type ].SetLineWidth(2);
305
306 // set X axis paras
307 debugHistos[ type ].GetXaxis()->SetLabelSize(0.1);
308 debugHistos[ type ].GetXaxis()->SetTitleSize(0.1);
309 debugHistos[ type ].GetXaxis()->SetTitleOffset(1.2);
310 debugHistos[ type ].GetXaxis()->SetTitle(Form("Time slice (%.1f ns/slice)", 1./2.));
311
312 // set Y axis paras
313 debugHistos[ type ].GetYaxis()->SetLabelSize(0.1);
314 debugHistos[ type ].GetYaxis()->SetTitleSize(0.1);
315 debugHistos[ type ].GetYaxis()->SetTitleOffset(0.3);
316 debugHistos[ type ].GetYaxis()->SetTitle("Amplitude (a.u.)");
317 }
318}
319
320void SaveHistograms( const char * loc_fname ){
321 // create the histogram file (replace if already existing)
322 TFile tf( loc_fname, "RECREATE");
323
324 hList.Write(); // write the major histograms into the top level directory
325
326 tf.Close(); // close the file
327}
Note: See TracBrowser for help on using the repository browser.