source: fact/tools/rootmacros/fpeak_cfd.C@ 16858

Last change on this file since 16858 was 12673, checked in by neise, 13 years ago
treating every 9th channel differently - cutting off 60 slices in the end
  • Property svn:executable set to *
File size: 10.7 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
26#include "openFits.h"
27#include "openFits.c"
28
29#include "discriminator.h"
30#include "discriminator.C"
31#include "zerosearch.h"
32#include "zerosearch.C"
33#include "factfir.C"
34
35#include "FOpenDataFile.h"
36#include "FOpenDataFile.c"
37
38#include "DrsCalibration.C"
39#include "DrsCalibration.h"
40
41#include "SpikeRemoval.h"
42#include "SpikeRemoval.C"
43
44bool breakout=false;
45
46int NEvents;
47vector<int16_t> AllPixelDataVector;
48vector<int16_t> StartCellVector;
49unsigned int CurrentEventID;
50size_t PXLxROI;
51UInt_t RegionOfInterest;
52UInt_t NumberOfPixels;
53
54size_t TriggerOffsetROI, RC;
55vector<float> Offset, Gain, TriggerOffset;
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( TString &title );
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_cfd(
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 = 8,
103 int avg2 = 8,
104 float threshold = 5.0,
105 int verbosityLevel = 1, // different verbosity levels can be implemented here
106 bool ProduceGraphic = false
107 )
108{
109gStyle->SetPalette(1,0);
110gROOT->SetStyle("Plain");
111
112TString histogramtitle;
113histogramtitle += datafilename;
114histogramtitle += " ";
115histogramtitle += drsfilename ;
116histogramtitle += " ";
117
118histogramtitle += avg1;
119histogramtitle += " ";
120histogramtitle += avg2;
121histogramtitle += " ";
122histogramtitle += threshold ;
123histogramtitle += " 100";
124cout << histogramtitle.Data() << endl;
125
126
127 TCanvas *cFiltered = NULL;
128 if (spikeDebug){
129 cFiltered = new TCanvas("cFiltered","filtered DRS Waveforms",410,410,400,400);
130 cFiltered->Divide(1, 4);
131 }
132
133
134
135// CFD filter settings
136 int k_cfd = 10;
137 vector<double> a_cfd(k_cfd, 0);
138 double b_cfd = 1.;
139 a_cfd[0]=-0.75;
140 a_cfd[k_cfd-1]=1.;
141
142 fits * datafile;
143 // Opens the raw data file and 'binds' the variables given as
144 // Parameters to the data file. So they are filled with
145 // raw data as soon as datafile->GetRow(int) is called.
146 NEvents = openDataFits( datafilename, &datafile,
147 AllPixelDataVector, StartCellVector, CurrentEventID,
148 RegionOfInterest, NumberOfPixels, PXLxROI, verbosityLevel);
149 if (NEvents == 0){
150 cout << "return code of OpenDataFile:" << datafilename<< endl;
151 cout << "is zero -> aborting." << endl;
152 return 1;
153 }
154
155 if (verbosityLevel > 0)
156 cout <<"number of events in file: "<< NEvents << "\t";
157 if ( nevents == -1 || nevents > NEvents ) nevents = NEvents; // -1 means all!
158 if (verbosityLevel > 0)
159 cout <<"of, which "<< nevents << "will be processed"<< endl;
160
161 if (verbosityLevel > 0)
162 cout <<"Total # of Pixel: "<< NumberOfPixels << "\t";
163 if ( npixel == -1 || npixel > (int)NumberOfPixels ) npixel = (int)NumberOfPixels; // -1 means all!
164 if (verbosityLevel > 0)
165 cout <<"of, which "<< npixel << "will be processed"<< endl;
166
167 RC = openCalibFits( drsfilename, Offset, Gain, TriggerOffset, TriggerOffsetROI);
168 if (RC == 0){
169 cout << "return code of openCalibFits:" << drsfilename << endl;
170 cout << "is zero -> aborting." << endl;
171 return 1;
172 }
173
174 BookHistos( histogramtitle );
175//-----------------------------------------------------------------------------
176// Loops over Every Event and Pixel
177//-----------------------------------------------------------------------------
178 cout << "Processing Events in total: " << nevents << endl;
179
180 for ( int ev = firstevent; ev < firstevent + nevents; ev++) {
181 // Get an Event --> consists of 1440 Pixel ...erm....data
182 datafile->GetRow( ev );
183
184 for ( int pix = firstpixel; pix < firstpixel+npixel; pix++ ){
185 if (verbosityLevel > 0){
186 if ( pix == (firstpixel + npixel - 1) ) {
187 if (ev == firstevent) {
188 printf ("%06d done", CurrentEventID);
189 fflush(stdout);
190 } else {
191 printf ("\b\b\b\b\b\b\b\b\b\b\b%06d done", CurrentEventID);
192 fflush(stdout);
193 }
194 }
195 }
196
197 if ((pix+1)%9==0){
198 applyDrsCalibration( Ameas,pix,12,60,
199 Offset, Gain, TriggerOffset,
200 RegionOfInterest, AllPixelDataVector, StartCellVector);
201 }else{
202 applyDrsCalibration( Ameas,pix,12,12,
203 Offset, Gain, TriggerOffset,
204 RegionOfInterest, AllPixelDataVector, StartCellVector);
205 }
206
207 // finds spikes in the raw data, and interpolates the value
208 // spikes are: 1 or 2 slice wide, positive non physical artifacts
209 removeSpikes (Ameas, Vcorr);
210
211 // filter Vcorr with sliding average using FIR filter function
212 sliding_avg(Vcorr, Vslide, avg1);
213 // filter Vslide with CFD using FIR filter function
214 factfir(b_cfd , a_cfd, k_cfd, Vslide, Vcfd);
215 // filter Vcfd with sliding average using FIR filter function
216 sliding_avg(Vcfd, Vcfd2, avg2);
217
218
219 // peaks in Ameas[] are found by searching for zero crossings
220 // in Vcfd2
221 // first Argument 1 means ... *rising* edge
222 // second Argument 1 means ... search with stepsize 1 ... 8 is okay as well
223 vector<Region> * zXings = zerosearch( Vcfd2 , 1 , 8);
224 // zXings means "zero cross ings"
225 EnlargeRegion( *zXings, 10, 10);
226 findAbsMaxInRegions( *zXings, Vslide);
227 removeMaximaBelow( *zXings, threshold, 0);
228 removeRegionWithMaxOnEdge( *zXings, 2);
229 removeRegionOnFallingEdge( *zXings, 100);
230
231 // fill maxima in Histogram
232 if (zXings->size() != 0 ){
233 for (unsigned int i=0; i<zXings->size(); i++){
234 if (verbosityLevel > 1){
235 cout << zXings->at(i).maxPos << ":\t"
236 << zXings->at(i).maxVal <<endl;
237 }
238 hAmplSpek_cfd->Fill(pix, zXings->at(i).maxVal);
239 }
240 }
241
242 if ( spikeDebug ){
243
244 // TODO do this correct. The vectors should be the rigt ones... this is just luck
245 debugHistos[Ameas_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
246 debugHistos[Vcorr_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
247 debugHistos[Vslide_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
248 debugHistos[Vcfd_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
249 debugHistos[Vcfd2_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
250
251 for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
252 debugHistos[Ameas_].SetBinContent(sl, Ameas[sl]);
253 debugHistos[Vcorr_].SetBinContent(sl, Vcorr[sl]);
254 debugHistos[Vslide_].SetBinContent( sl, Vslide[sl] );
255 debugHistos[Vcfd_].SetBinContent( sl, Vcfd[sl] );
256 debugHistos[Vcfd2_].SetBinContent( sl, Vcfd2[sl] );
257 }
258
259 cFiltered->cd(1);
260 gPad->SetGrid();
261 debugHistos[Ameas_].Draw();
262
263 cFiltered->cd(2);
264 gPad->SetGrid();
265 debugHistos[Vcorr_].Draw();
266
267 cFiltered->cd(3);
268 gPad->SetGrid();
269 debugHistos[Vslide_].Draw();
270
271 TBox *OneBox;
272 vector<TBox*> MyBoxes;
273 for (unsigned int i=0; i<zXings->size(); i++){
274 OneBox = new TBox(
275 zXings->at(i).maxPos -10 ,
276 zXings->at(i).maxVal -0.5,
277 zXings->at(i).maxPos +10 ,
278 zXings->at(i).maxVal +0.5);
279 OneBox->SetLineColor(kBlue);
280 OneBox->SetLineWidth(1);
281 OneBox->SetFillStyle(0);
282 OneBox->SetFillColor(kRed);
283 MyBoxes.push_back(OneBox);
284 OneBox->Draw();
285 }
286
287 cFiltered->cd(4);
288 gPad->SetGrid();
289 debugHistos[Vcfd2_].Draw();
290 TLine *zeroline = new TLine(0, 0, 1024, 0);
291 zeroline->SetLineColor(kBlue);
292 zeroline->Draw();
293
294 cFiltered->Update();
295 //Process gui events asynchronously during input
296 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
297 timer.TurnOn();
298 TString input = Getline("Type 'q' to exit, <return> to go on: ");
299 timer.TurnOff();
300 if (input=="q\n") {
301 breakout=true;
302 }
303
304 //TODO!!!!!!!!!
305 // do some Garbage collection ...
306 }// end of if(spikeDebug)
307
308 delete zXings;
309 if (breakout)
310 break;
311
312 } // end of loop over pixels
313
314 if (breakout)
315 break;
316 } // end of loop over pixels
317
318 cout << "Event processing done." << endl;
319
320if (ProduceGraphic){
321 TCanvas * cSpektrum;
322 cSpektrum = new TCanvas("cSpektrum","Amplitude Spektrum",10,10,400,400);
323 cSpektrum->Divide(1,1);
324 cSpektrum->cd(1);
325 hAmplSpek_cfd->Draw("COLZ");
326 cSpektrum->Modified();
327 cSpektrum->Update();
328}
329 if ( strlen(OutRootFileName) != 0){
330 SaveHistograms( OutRootFileName );
331 }
332 return( 0 );
333}
334
335// booking and parameter settings for all histos
336void BookHistos( TString& title){
337 // histograms for baseline extraction
338
339
340
341 hAmplSpek_cfd = new TH2F("hAmplSpek_cfd", title,1440,-0.5,1439.5, 256, -27.5, 100.5);
342 hAmplSpek_cfd->GetXaxis()->SetTitle( "pixel" );
343 hAmplSpek_cfd->GetYaxis()->SetTitle( "amplitude in mV" );
344 hList.Add( hAmplSpek_cfd );
345
346
347 debugHistos = new TH1F[ NumberOfDebugHistoTypes ];
348 for ( int type = 0; type < NumberOfDebugHistoTypes; type++){
349 debugHistos[ type ].SetBins(1024, 0, 1024);
350 debugHistos[ type ].SetLineColor(1);
351 debugHistos[ type ].SetLineWidth(2);
352
353 // set X axis paras
354 debugHistos[ type ].GetXaxis()->SetLabelSize(0.1);
355 debugHistos[ type ].GetXaxis()->SetTitleSize(0.1);
356 debugHistos[ type ].GetXaxis()->SetTitleOffset(1.2);
357 debugHistos[ type ].GetXaxis()->SetTitle(Form("Time slice (%.1f ns/slice)", 1./2.));
358
359 // set Y axis paras
360 debugHistos[ type ].GetYaxis()->SetLabelSize(0.1);
361 debugHistos[ type ].GetYaxis()->SetTitleSize(0.1);
362 debugHistos[ type ].GetYaxis()->SetTitleOffset(0.3);
363 debugHistos[ type ].GetYaxis()->SetTitle("Amplitude (a.u.)");
364 }
365}
366
367void SaveHistograms( const char * loc_fname){
368 // create the histogram file (replace if already existing)
369 TFile tf( loc_fname, "RECREATE");
370
371 hList.Write(); // write the major histograms into the top level directory
372
373 tf.Close(); // close the file
374}
Note: See TracBrowser for help on using the repository browser.