source: fact/tools/rootmacros/tpeak.C@ 12426

Last change on this file since 12426 was 12392, checked in by neise, 13 years ago
tidied up tpeak.C and put an entry into README.txt
File size: 15.5 KB
Line 
1/********************** Peak-Template ***********************
2
3this makro shell search for Peaks in data and form a template
4of these, so it can be used for detection of other peaks
5
6*************************************************************/
7
8//root libraries
9
10#include <TROOT.h>
11#include <TCanvas.h>
12#include <TProfile.h>
13#include <TTimer.h>
14#include <TH1F.h>
15#include <TH2F.h>
16#include <Getline.h>
17#include <TLine.h>
18#include <TBox.h>
19#include <TMath.h>
20#include <TFile.h>
21#include <TStyle.h>
22
23#include <stdio.h>
24#include <stdint.h>
25#include <cstdio>
26
27#define NPIX 1440
28#define NCELL 1024
29#define FAD_MAX_SAMPLES 1024
30
31#define HAVE_ZLIB
32
33//rootmacros
34
35#include "fits.h"
36#include "FOpenCalibFile.c"
37
38#include "discriminator.h"
39#include "discriminator.C"
40#include "zerosearch.h"
41#include "zerosearch.C"
42#include "factfir.C"
43
44#include "FOpenDataFile.h"
45#include "FOpenDataFile.c"
46
47#include "DrsCalibration.C"
48#include "DrsCalibration.h"
49
50#include "SpikeRemoval.h"
51#include "SpikeRemoval.C"
52
53//-----------------------------------------------------------------------------
54// Decleration of variables
55//-----------------------------------------------------------------------------
56bool breakout=false;
57
58int NEvents;
59vector<int16_t> AllPixelDataVector;
60vector<int16_t> StartCellVector;
61unsigned int CurrentEventID;
62size_t PXLxROI;
63UInt_t RegionOfInterest;
64UInt_t NumberOfPixels;
65
66size_t drs_n;
67vector<float> drs_basemean;
68vector<float> drs_gainmean;
69vector<float> drs_triggeroffsetmean;
70
71vector<float> Ameas(FAD_MAX_SAMPLES); // copy of the data (measured amplitude
72vector<float> Vcorr(FAD_MAX_SAMPLES); // corrected Values
73vector<float> Vslide(FAD_MAX_SAMPLES); // sliding average result
74vector<float> Vcfd(FAD_MAX_SAMPLES); // CDF result
75vector<float> Vcfd2(FAD_MAX_SAMPLES); // CDF result + 2nd sliding average
76
77// histograms
78const int NumberOfDebugHistoTypes = 7;
79const unsigned int
80 Ameas_ = 0,
81 N1mean_ = 1,
82 Vcorr_ = 2,
83 Vtest_ = 3,
84 Vslide_ = 4,
85 Vcfd_ = 5,
86 Vcfd2_ = 6;
87
88TH1F* h;
89TH1F *debugHistos;
90TH2F *hPeakOverlay; //histogrammm for overlay of detected Peaks
91int hPeakOverlayXaxisLeft,hPeakOverlayXaxisRight;
92
93TObjArray hList;
94TObjArray hListBaseline;
95
96void BookHistos( );
97void SaveHistograms( const char * );
98
99//-----------------------------------------------------------------------------
100// Main
101//-----------------------------------------------------------------------------
102int tpeak(
103 char *datafilename = "data/20111016_013.fits.gz",
104 const char *drsfilename = "../../20111016_011.drs.fits.gz",
105 const char *OutRootFileName = "../analysis/fpeak_cdf.Coutput.root",
106 int firstevent = 0,
107 int nevents = -1,
108 int firstpixel = 0,
109 int npixel = -1,
110 bool spikeDebug = false,
111 int avg1 = 14,
112 int avg2 = 8,
113 int OverlayWindowLeft = 50,
114 int OverlayWindowRight = 150,
115 int verbosityLevel = 1, // different verbosity levels can be implemented here
116 bool ProduceGraphic = true
117 )
118
119{
120hPeakOverlayXaxisLeft = OverlayWindowLeft;
121hPeakOverlayXaxisRight = OverlayWindowRight;
122
123 gStyle->SetPalette(1,0);
124 gROOT->SetStyle("Plain");
125
126//-----------------------------------------------------------------------------
127// Create (pointer to) Canvases, which are used in every run,
128// also in 'non-debug' runs
129//-----------------------------------------------------------------------------
130 // Canvases only need if spike Debug, but I want to deklare
131 // the pointers anyway ...
132 TCanvas *cFiltered = NULL;
133 if (spikeDebug){
134 cFiltered = new TCanvas("cFiltered","filtered DRS Waveforms",1,10,400,300);
135 cFiltered->Divide(1, 3);
136 }
137 // Canvases to show the peak template
138 TCanvas *cPeakOverlay = NULL;
139 cPeakOverlay = new TCanvas("cPeakOverlay", "Overlay of detected Peaks", 1, 310, 400, 400);
140
141
142//-----------------------------------------------------------------------------
143// read FACT raw data
144// * remove spikes
145// * calculate baseline
146// * find peaks (CFD and normal discriminator)
147// * compute pulse height and pulse integral spektrum of the peaks
148//-----------------------------------------------------------------------------
149
150
151//-----------------------------------------------------------------------------
152// Filter-Settings
153//-----------------------------------------------------------------------------
154// CFD filter settings
155 int k_cfd = 10;
156 vector<double> a_cfd(k_cfd, 0);
157 double b_cfd = 1.;
158 a_cfd[0]=-0.75;
159 a_cfd[k_cfd-1]=1.;
160
161//-----------------------------------------------------------------------------
162// prepare datafile
163//-----------------------------------------------------------------------------
164
165// Open the data file
166
167 fits * datafile;
168 // Opens the raw data file and 'binds' the variables given as
169 // Parameters to the data file. So they are filled with
170 // raw data as soon as datafile->GetRow(int) is called.
171 NEvents = OpenDataFile( datafilename, &datafile,
172 AllPixelDataVector, StartCellVector, CurrentEventID,
173 RegionOfInterest, NumberOfPixels, PXLxROI, verbosityLevel);
174 if (NEvents == 0){
175 cout << "return code of OpenDataFile:" << datafilename<< endl;
176 cout << "is zero -> aborting." << endl;
177 return 1;
178 }
179
180 if (verbosityLevel > 0)
181 cout <<"number of events in file: "<< NEvents << "\t";
182 if ( nevents == -1 || nevents > NEvents ) nevents = NEvents; // -1 means all!
183 if (verbosityLevel > 0)
184 cout <<"of, which "<< nevents << "will be processed"<< endl;
185
186 if (verbosityLevel > 0)
187 cout <<"Total # of Pixel: "<< NumberOfPixels << "\t";
188 if ( npixel == -1 || npixel > (int)NumberOfPixels ) npixel = (int)NumberOfPixels; // -1 means all!
189 if (verbosityLevel > 0)
190 cout <<"of, which "<< npixel << "will be processed"<< endl;
191
192//Get the DRS calibration
193 FOpenCalibFile( drsfilename,
194 drs_basemean,
195 drs_gainmean,
196 drs_triggeroffsetmean,
197 drs_n);
198
199// Book the histograms
200 BookHistos( );
201
202//-----------------------------------------------------------------------------
203// Loops over Every Event and Pixel
204//-----------------------------------------------------------------------------
205 for ( int ev = firstevent; ev < firstevent + nevents; ev++) {
206 // Get an Event --> consists of 1440 Pixel ...erm....data
207 datafile->GetRow( ev );
208
209 //-------------------------------------
210 // Loop over every Pixel of Event
211 //-------------------------------------
212
213 for ( int pix = firstpixel; pix < firstpixel + npixel; pix++ ){
214
215 if (verbosityLevel > 0){
216 if (pix == firstpixel){
217 cout << "Processing Event: " << CurrentEventID << "/" << nevents << endl;
218 }
219 }
220
221 applyDrsCalibration( Ameas,pix,12,12,
222 drs_basemean, drs_gainmean, drs_triggeroffsetmean,
223 RegionOfInterest, AllPixelDataVector, StartCellVector);
224
225 // finds spikes in the raw data, and interpolates the value
226 // spikes are: 1 or 2 slice wide, positive non physical artifacts
227 removeSpikes (Ameas, Vcorr);
228
229 // filter Vcorr with sliding average using FIR filter function
230 sliding_avg(Vcorr, Vslide, avg1);
231 // filter Vslide with CFD using FIR filter function
232 factfir(b_cfd , a_cfd, k_cfd, Vslide, Vcfd);
233 // filter Vcfd with sliding average using FIR filter function
234 sliding_avg(Vcfd, Vcfd2, avg2);
235
236 // peaks in Ameas[] are found by searching for zero crossings
237 // in Vcfd2
238 // first Argument 1 means ... *rising* edge
239 // second Argument 1 means ... search with stepsize 1 ... 10 is okay as well
240 vector<Region> * zXings = zerosearch( Vcfd2 , 1 , 1);
241 // zXings means "zero cross ings"
242 EnlargeRegion(*zXings, 10, 10);
243 findAbsMaxInRegions(*zXings, Vslide);
244 removeMaximaBelow( *zXings, 3.0);
245 removeRegionWithMaxOnEdge( *zXings, 2);
246 removeRegionOnFallingEdge( *zXings, 100);
247
248 //following Code should produce the Overlay of peaks
249 vector<Region>::iterator it;
250
251 for (it = zXings->begin() ; it < zXings->end() ; it++){
252 int Left = it->maxPos - OverlayWindowLeft;
253 int Right = it->maxPos + OverlayWindowRight;
254 if (Left < 0)
255 Left =0;
256 if (Right > (int)Vcorr.size() )
257 Right=Vcorr.size() ;
258 for ( int pos = Left; pos < Right; pos++){
259
260 hPeakOverlay->Fill( pos - (it->maxPos), Vcorr[pos]) ;
261 }
262 }
263
264
265 if ( spikeDebug ){
266
267
268 // TODO do this correct. The vectors should be the rigt ones... this is just luck
269 debugHistos[Ameas_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
270 debugHistos[Vcorr_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
271 debugHistos[Vslide_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
272 debugHistos[Vcfd_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
273 debugHistos[Vcfd2_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
274
275 for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
276 debugHistos[Ameas_].SetBinContent(sl, Ameas[sl]);
277 debugHistos[Vcorr_].SetBinContent(sl, Vcorr[sl]);
278 debugHistos[Vslide_].SetBinContent( sl, Vslide[sl] );
279 debugHistos[Vcfd_].SetBinContent( sl, Vcfd[sl] );
280 debugHistos[Vcfd2_].SetBinContent( sl, Vcfd2[sl] );
281 }
282
283
284 cFiltered->cd(1);
285 gPad->SetGrid();
286 debugHistos[Vcorr_].Draw();
287
288 cFiltered->cd(2);
289 gPad->SetGrid();
290 debugHistos[Vslide_].Draw();
291
292 TBox *OneBox;
293 vector<TBox*> MyBoxes;
294 for (unsigned int i=0; i<zXings->size(); i++){
295 OneBox = new TBox(
296 zXings->at(i).maxPos -10 ,
297 zXings->at(i).maxVal -0.5,
298 zXings->at(i).maxPos +10 ,
299 zXings->at(i).maxVal +0.5);
300 OneBox->SetLineColor(kBlue);
301 OneBox->SetLineWidth(1);
302 OneBox->SetFillStyle(0);
303 OneBox->SetFillColor(kRed);
304 MyBoxes.push_back(OneBox);
305 OneBox->Draw();
306 }
307
308 cFiltered->cd(3);
309 gPad->SetGrid();
310 debugHistos[Vcfd2_].Draw();
311 TLine *zeroline = new TLine(0, 0, 1024, 0);
312 zeroline->SetLineColor(kBlue);
313 zeroline->Draw();
314
315 cFiltered->Update();
316
317 //OVERLAY PEAKS
318 cPeakOverlay->cd();
319 hPeakOverlay->Draw("COLZ");
320 cPeakOverlay->Modified();
321 cPeakOverlay->Update();
322
323 //Process gui events asynchronously during input
324 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
325 timer.TurnOn();
326 TString input = Getline("Type 'q' to exit, <return> to go on: ");
327 timer.TurnOff();
328 if (input=="q\n") {
329 breakout=true;
330 }
331
332 //TODO!!!!!!!!!
333 // do some Garbage collection ...
334 // all the Objects on the heap should be deleted here.
335
336 }// end of if(spikeDebug)
337
338 delete zXings;
339 if (breakout)
340 break;
341
342 }
343 //-------------------------------------
344 // end of loop over pixels
345 //-------------------------------------
346 if (ProduceGraphic){
347 if (ev % 50 == 0){
348 //OVERLAY PEAKS
349 cPeakOverlay->cd();
350 hPeakOverlay->Draw("COLZ");
351 cPeakOverlay->Modified();
352 cPeakOverlay->Update();
353 }
354 }
355 if (breakout)
356 break;
357 }
358//-------------------------------------
359// end of loop over Events
360//-------------------------------------
361 if (ProduceGraphic){
362 //OVERLAY PEAKS
363 cPeakOverlay->cd();
364 hPeakOverlay->Draw("COLZ");
365 cPeakOverlay->Modified();
366 cPeakOverlay->Update();
367 }
368
369 SaveHistograms( OutRootFileName );
370
371 return( 0 );
372}
373
374// booking and parameter settings for all histos
375void BookHistos( ){
376
377 debugHistos = new TH1F[ NumberOfDebugHistoTypes ];
378 for ( int type = 0; type < NumberOfDebugHistoTypes; type++){
379 debugHistos[ type ].SetBins(1024, 0, 1024);
380 debugHistos[ type ].SetLineColor(1);
381 debugHistos[ type ].SetLineWidth(2);
382
383 // set X axis paras
384 debugHistos[ type ].GetXaxis()->SetLabelSize(0.1);
385 debugHistos[ type ].GetXaxis()->SetTitleSize(0.1);
386 debugHistos[ type ].GetXaxis()->SetTitleOffset(1.2);
387 debugHistos[ type ].GetXaxis()->SetTitle(Form("Time slice (%.1f ns/slice)", 1./2.));
388
389 // set Y axis paras
390 debugHistos[ type ].GetYaxis()->SetLabelSize(0.1);
391 debugHistos[ type ].GetYaxis()->SetTitleSize(0.1);
392 debugHistos[ type ].GetYaxis()->SetTitleOffset(0.3);
393 debugHistos[ type ].GetYaxis()->SetTitle("Amplitude (a.u.)");
394 }
395 hPeakOverlay = new TH2F("hPeakOverlay", "Overlay of detected Peaks",
396 hPeakOverlayXaxisLeft + hPeakOverlayXaxisRight , (-1*hPeakOverlayXaxisLeft)-0.5, hPeakOverlayXaxisRight-0.5 ,
397 4096, -48.5, 1999.5 );
398 hPeakOverlay->GetXaxis()->SetTitle( "Timeslices" );
399 hPeakOverlay->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
400 hList.Add( hPeakOverlay );
401
402}
403
404void SaveHistograms(const char * loc_fname ){
405 // create the histogram file (replace if already existing)
406 TFile tf( loc_fname, "RECREATE");
407
408 hList.Write(); // write the major histograms into the top level directory
409
410 tf.Close(); // close the file
411} // end of SaveHistograms( char * loc_fname )
412
Note: See TracBrowser for help on using the repository browser.