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

Last change on this file since 13616 was 12720, checked in by neise, 14 years ago
removed typo
File size: 19.2 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
37#include "openFits.h"
38#include "openFits.c"
39
40#include "discriminator.h"
41#include "discriminator.C"
42#include "zerosearch.h"
43#include "zerosearch.C"
44#include "factfir.C"
45
46#include "DrsCalibration.C"
47#include "DrsCalibration.h"
48
49#include "SpikeRemoval.h"
50#include "SpikeRemoval.C"
51
52//-----------------------------------------------------------------------------
53// Decleration of variables
54//-----------------------------------------------------------------------------
55bool breakout=false;
56
57int NEvents;
58vector<int16_t> AllPixelDataVector;
59vector<int16_t> StartCellVector;
60unsigned int CurrentEventID;
61size_t PXLxROI;
62UInt_t RegionOfInterest;
63UInt_t NumberOfPixels;
64
65size_t TriggerOffsetROI, RC;
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
77typedef struct{
78 float maxAmpl;
79 float countOfMax;
80 } OverlayMaximum;
81vector<OverlayMaximum> SingleMaximum;
82// histograms
83const int NumberOfDebugHistoTypes = 7;
84
85const unsigned int
86 Ameas_ = 0,
87 N1mean_ = 1,
88 Vcorr_ = 2,
89 Vtest_ = 3,
90 Vslide_ = 4,
91 Vcfd_ = 5,
92 Vcfd2_ = 6;
93
94TH1F* h;
95TH1F *debugHistos;
96TH2F *hPeakOverlay; //histogrammm for overlay of detected Peaks
97TH1F *hMaxPeakOverlay;
98//TH2F *hPixelPeakOverlay;
99//TH2F *hEventPeakOverlay;
100int hPeakOverlayXaxisLeft,hPeakOverlayXaxisRight;
101
102TObjArray hList;
103TObjArray hListBaseline;
104
105void BookHistos( );
106void SaveHistograms( const char * );
107
108//-----------------------------------------------------------------------------
109// Main
110//-----------------------------------------------------------------------------
111int tpeak(
112 char *datafilename = "data/20111016_013.fits.gz",
113 const char *drsfilename = "../../20111016_011.drs.fits.gz",
114 const char *OutRootFileName = "../analysis/tpeak_cdf.Coutput.root",
115 int firstevent = 0,
116 int nevents = -1,
117 int firstpixel = 0,
118 int npixel = -1,
119 bool spikeDebug = false,
120 int avg1 = 14,
121 int avg2 = 8,
122 int OverlayWindowLeft = 50,
123 int OverlayWindowRight = 150,
124 int verbosityLevel = 1, // different verbosity levels can be implemented here
125 bool ProduceGraphic = true
126 )
127
128{
129 hPeakOverlayXaxisLeft = OverlayWindowLeft;
130 hPeakOverlayXaxisRight = OverlayWindowRight;
131
132 gStyle->SetPalette(1,0);
133 gROOT->SetStyle("Plain");
134
135//-----------------------------------------------------------------------------
136// Create (pointer to) Canvases, which are used in every run,
137// also in 'non-debug' runs
138//-----------------------------------------------------------------------------
139 // Canvases only need if spike Debug, but I want to deklare
140 // the pointers anyway ...
141 TCanvas *cFiltered = NULL;
142 if (spikeDebug){
143 cFiltered = new TCanvas("cFiltered","filtered DRS Waveforms", 1,310,400,300);
144 cFiltered->Divide(1, 3);
145 }
146 // Canvases to show the peak template
147 TCanvas *cPeakOverlay = NULL;
148 cPeakOverlay = new TCanvas("cPeakOverlay", "Overlay of detected Peaks", 1, 1, 1200, 800);
149 cPeakOverlay->Divide(2,1);
150// // All peaks of one Pixel
151// TCanvas *cPixelPeakOverlay = NULL;
152// cPixelPeakOverlay = new TCanvas("cPixelPeakOverlay", "Overlay of detected Peaks of all Events of one Pixel", 401, 1, 400, 400);
153// // All Peaks of one Event
154// TCanvas *cYProjection = NULL;
155// cYProjection = new TCanvas("cYProjection", "Y Projection of peak overlay for each bin ", 802, 1, 400, 400);
156
157
158
159//-----------------------------------------------------------------------------
160// read FACT raw data
161// * remove spikes
162// * calculate baseline
163// * find peaks (CFD and normal discriminator)
164// * compute pulse height and pulse integral spektrum of the peaks
165//-----------------------------------------------------------------------------
166
167
168//-----------------------------------------------------------------------------
169// Filter-Settings
170//-----------------------------------------------------------------------------
171// CFD filter settings
172 int k_cfd = 10;
173 vector<double> a_cfd(k_cfd, 0);
174 double b_cfd = 1.;
175 a_cfd[0]=-0.75;
176 a_cfd[k_cfd-1]=1.;
177
178//-----------------------------------------------------------------------------
179// prepare datafile
180//-----------------------------------------------------------------------------
181
182// Open the data file
183
184 fits * datafile;
185 // Opens the raw data file and 'binds' the variables given as
186 // Parameters to the data file. So they are filled with
187 // raw data as soon as datafile->GetRow(int) is called.
188 NEvents = openDataFits( datafilename, &datafile,
189 AllPixelDataVector, StartCellVector, CurrentEventID,
190 RegionOfInterest, NumberOfPixels, PXLxROI, verbosityLevel);
191 if (NEvents == 0){
192 cout << "return code of OpenDataFile:" << datafilename<< endl;
193 cout << "is zero -> aborting." << endl;
194 return 1;
195 }
196
197 if (verbosityLevel > 0)
198 cout <<"number of events in file: "<< NEvents << "\t";
199 if ( nevents == -1 || nevents > NEvents ) nevents = NEvents; // -1 means all!
200 if (verbosityLevel > 0)
201 cout <<"of, which "<< nevents << "will be processed"<< endl;
202
203 if (verbosityLevel > 0)
204 cout <<"Total # of Pixel: "<< NumberOfPixels << "\t";
205 if ( npixel == -1 || npixel > (int)NumberOfPixels ) npixel = (int)NumberOfPixels; // -1 means all!
206 if (verbosityLevel > 0)
207 cout <<"of, which "<< npixel << "will be processed"<< endl;
208//Get the DRS calibration
209 RC = openCalibFits( drsfilename, drs_basemean, drs_gainmean, drs_triggeroffsetmean, TriggerOffsetROI);
210 if (RC == 0){
211 cout << "return code of openCalibFits:" << drsfilename << endl;
212 cout << "is zero -> aborting." << endl;
213 return 1;
214 }
215// Book the histograms
216 BookHistos( );
217
218//-----------------------------------------------------------------------------
219// Loops over Every Event and Pixel
220//-----------------------------------------------------------------------------
221 for ( int ev = firstevent; ev < firstevent + nevents; ev++) {
222 // Get an Event --> consists of 1440 Pixel ...erm....data
223 datafile->GetRow( ev );
224
225 //-------------------------------------
226 // Loop over every Pixel of Event
227 //-------------------------------------
228
229 for ( int pix = firstpixel; pix < firstpixel + npixel; pix++ ){
230 if (verbosityLevel > 0){
231 if (pix == firstpixel){
232 cout << "Processing Event: " << CurrentEventID << "/" << nevents << endl;
233 }
234 }
235
236 applyDrsCalibration( Ameas,pix,12,12,
237 drs_basemean, drs_gainmean, drs_triggeroffsetmean,
238 RegionOfInterest, AllPixelDataVector, StartCellVector);
239
240 // finds spikes in the raw data, and interpolates the value
241 // spikes are: 1 or 2 slice wide, positive non physical artifacts
242 removeSpikes (Ameas, Vcorr);
243
244 // filter Vcorr with sliding average using FIR filter function
245 sliding_avg(Vcorr, Vslide, avg1);
246 // filter Vslide with CFD using FIR filter function
247 factfir(b_cfd , a_cfd, k_cfd, Vslide, Vcfd);
248 // filter Vcfd with sliding average using FIR filter function
249 sliding_avg(Vcfd, Vcfd2, avg2);
250
251 // peaks in Ameas[] are found by searching for zero crossings
252 // in Vcfd2
253 // first Argument 1 means ... *rising* edge
254 // second Argument 1 means ... search with stepsize 1 ... 10 is okay as well
255 vector<Region> * zXings = zerosearch( Vcfd2 , 1 , 8);
256 // zXings means "zero cross ings"
257 EnlargeRegion(*zXings, 10, 10);
258 findAbsMaxInRegions(*zXings, Vslide);
259 removeMaximaBelow( *zXings, 3.0);
260 removeRegionWithMaxOnEdge( *zXings, 2);
261 removeRegionOnFallingEdge( *zXings, 100);
262
263 //following Code should produce the Overlay of peaks
264 vector<Region>::iterator it;
265 for (it = zXings->begin() ; it < zXings->end() ; it++){
266 if (it->maxPos < 12 || it->maxPos > RegionOfInterest-12)
267 continue;
268 //domstest->Fill(it->maxPos);
269 int Left = it->maxPos - OverlayWindowLeft;
270 int Right = it->maxPos + OverlayWindowRight;
271 if (Left < 0)
272 Left =0;
273 if (Right > (int)Vcorr.size() )
274 Right=Vcorr.size() ;
275 for ( int pos = Left; pos < Right; pos++){
276 hPeakOverlay->Fill( pos - (it->maxPos), Vslide[pos]) ;
277// hPixelPeakOverlay->Fill( pos - (it->maxPos), Vcorr[pos]) ;
278// hEventPeakOverlay->Fill( pos - (it->maxPos), 2*Vcorr[pos]) ;
279 }
280 }
281
282
283 if ( spikeDebug ){
284 // TODO do this correct. The vectors should be the rigt ones... this is just luck
285 debugHistos[Ameas_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
286 debugHistos[Vcorr_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
287 debugHistos[Vslide_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
288 debugHistos[Vcfd_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
289 debugHistos[Vcfd2_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
290
291 for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
292 debugHistos[Ameas_].SetBinContent(sl, Ameas[sl]);
293 debugHistos[Vcorr_].SetBinContent(sl, Vcorr[sl]);
294 debugHistos[Vslide_].SetBinContent( sl, Vslide[sl] );
295 debugHistos[Vcfd_].SetBinContent( sl, Vcfd[sl] );
296 debugHistos[Vcfd2_].SetBinContent( sl, Vcfd2[sl] );
297 }
298
299
300 cFiltered->cd(1);
301 gPad->SetGrid();
302 debugHistos[Vcorr_].Draw();
303
304 cFiltered->cd(2);
305 gPad->SetGrid();
306 debugHistos[Vslide_].Draw();
307
308 TBox *OneBox;
309 vector<TBox*> MyBoxes;
310 for (unsigned int i=0; i<zXings->size(); i++){
311 OneBox = new TBox(
312 zXings->at(i).maxPos -10 ,
313 zXings->at(i).maxVal -0.5,
314 zXings->at(i).maxPos +10 ,
315 zXings->at(i).maxVal +0.5);
316 OneBox->SetLineColor(kBlue);
317 OneBox->SetLineWidth(1);
318 OneBox->SetFillStyle(0);
319 OneBox->SetFillColor(kRed);
320 MyBoxes.push_back(OneBox);
321 OneBox->Draw();
322 }
323
324 cFiltered->cd(3);
325 gPad->SetGrid();
326 debugHistos[Vcfd2_].Draw();
327 TLine *zeroline = new TLine(0, 0, 1024, 0);
328 zeroline->SetLineColor(kBlue);
329 zeroline->Draw();
330
331 cFiltered->Update();
332
333 //OVERLAY PEAKS
334 cPeakOverlay->cd(1);
335 hPeakOverlay->Draw("COLZ");
336 cPeakOverlay->Modified();
337 cPeakOverlay->Update();
338
339 //Process gui events asynchronously during input
340 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
341 timer.TurnOn();
342 TString input = Getline("Type 'q' to exit, <return> to go on: ");
343 timer.TurnOff();
344 if (input=="q\n") {
345 breakout=true;
346 }
347
348 //TODO!!!!!!!!!
349 // do some Garbage collection ...
350 // all the Objects on the heap should be deleted here.
351
352 }// end of if(spikeDebug)
353
354 delete zXings;
355 if (breakout)
356 break;
357
358 }
359 //-------------------------------------
360 // end of loop over pixels
361 //-------------------------------------
362 if (ProduceGraphic){
363 if (ev % 50 == 0){
364
365 //OVERLAY PEAKS
366 cPeakOverlay->cd(1);
367 hPeakOverlay->Draw("COLZ");
368 cPeakOverlay->cd(2);
369 // hMaxPeakOverlay->Draw("");
370 cPeakOverlay->Modified();
371 cPeakOverlay->Update();
372 }
373 }
374 if (breakout)
375 break;
376 }
377
378//-------------------------------------
379// end of loop over Events
380//-------------------------------------
381
382//-------------------------------------
383// Histogramm of Maximas in Overlay Spectra
384//-------------------------------------
385 // resize vector SingleMaximum to number of timeslices in Overlay Spectra
386 SingleMaximum.resize(hPeakOverlay->GetNbinsX());
387
388 //put maximumvalue of every Y-projection of every timeslice into vector
389 for (int TimeSlice = 0; TimeSlice <= hPeakOverlay->GetNbinsX(); TimeSlice++ ){
390 TH1D* hProjPeak = hPeakOverlay->ProjectionY("proj_py", TimeSlice, TimeSlice, "");
391 SingleMaximum[ TimeSlice ].maxAmpl = (hProjPeak->GetMaximumBin() * 0.5) - 54;
392 SingleMaximum[ TimeSlice ].countOfMax = hProjPeak->GetBinContent( hProjPeak->GetMaximumBin() );
393 hMaxPeakOverlay->SetBinContent(TimeSlice, SingleMaximum[ TimeSlice ].maxAmpl );
394 }
395 hMaxPeakOverlay->Fit("landau", "", "", -50, 250);
396
397//-------------------------------------
398// Draw Histograms
399//-------------------------------------
400 if (ProduceGraphic){
401 //OVERLAY PEAKS
402 cPeakOverlay->cd(1);
403 hPeakOverlay->Draw("COLZ");
404 cPeakOverlay->cd(2);
405 hMaxPeakOverlay->Draw("");
406 cPeakOverlay->Modified();
407 cPeakOverlay->Update();
408
409// cPixelPeakOverlay->cd();
410// hPixelPeakOverlay->Draw("COLZ");
411// cPixelPeakOverlay->Modified();
412// cPixelPeakOverlay->Update();
413
414 }
415 SaveHistograms( OutRootFileName );
416
417 return( 0 );
418}
419
420// booking and parameter settings for all histos
421void BookHistos( ){
422
423 debugHistos = new TH1F[ NumberOfDebugHistoTypes ];
424 for ( int type = 0; type < NumberOfDebugHistoTypes; type++){
425 debugHistos[ type ].SetBins(1024, 0, 1024);
426 debugHistos[ type ].SetLineColor(1);
427 debugHistos[ type ].SetLineWidth(2);
428
429 // set X axis paras
430 debugHistos[ type ].GetXaxis()->SetLabelSize(0.1);
431 debugHistos[ type ].GetXaxis()->SetTitleSize(0.1);
432 debugHistos[ type ].GetXaxis()->SetTitleOffset(1.2);
433 debugHistos[ type ].GetXaxis()->SetTitle(Form("Time slice (%.1f ns/slice)", 1./2.));
434
435 // set Y axis paras
436 debugHistos[ type ].GetYaxis()->SetLabelSize(0.1);
437 debugHistos[ type ].GetYaxis()->SetTitleSize(0.1);
438 debugHistos[ type ].GetYaxis()->SetTitleOffset(0.3);
439 debugHistos[ type ].GetYaxis()->SetTitle("Amplitude (a.u.)");
440 }
441/*
442 ProjHistos = new TH1D[ NumberOfTimeSlices ];
443 for ( int type = 0; type < NumberOfTimeSlices; type++){
444 ProjHistos[ type ].SetBins(1024, 0, 1024);
445 ProjHistos[ type ].SetLineColor(1);
446 ProjHistos[ type ].SetLineWidth(2);
447
448 // set X axis paras
449 ProjHistos[ type ].GetXaxis()->SetLabelSize(0.1);
450 ProjHistos[ type ].GetXaxis()->SetTitleSize(0.1);
451 ProjHistos[ type ].GetXaxis()->SetTitleOffset(1.2);
452 ProjHistos[ type ].GetXaxis()->SetTitle(Form("Time slice (%.1f ns/slice)", 1./2.));
453
454 // set Y axis paras
455 ProjHistos[ type ].GetYaxis()->SetLabelSize(0.1);
456 ProjHistos[ type ].GetYaxis()->SetTitleSize(0.1);
457 ProjHistos[ type ].GetYaxis()->SetTitleOffset(0.3);
458 ProjHistos[ type ].GetYaxis()->SetTitle("Amplitude (a.u.)");
459 }
460*/
461 hPeakOverlay = new TH2F("hPeakOverlay", "Overlay of detected Peaks",
462 hPeakOverlayXaxisLeft + hPeakOverlayXaxisRight , (-1*hPeakOverlayXaxisLeft)-0.5, hPeakOverlayXaxisRight-0.5 ,
463 512, -55.5, 200.5 );
464 hPeakOverlay->GetXaxis()->SetTitle( "Timeslices" );
465 hPeakOverlay->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
466 //hPeakOverlay->SetBit(TH2F::kCanRebin);
467 hList.Add( hPeakOverlay );
468
469 hMaxPeakOverlay = new TH1F("hMaxPeakOverlay", "Single Peak derived with Maximum of above Spektrum",
470 hPeakOverlayXaxisLeft + hPeakOverlayXaxisRight ,
471 (-1*hPeakOverlayXaxisLeft)-0.5,
472 hPeakOverlayXaxisRight-0.5
473 );
474 hMaxPeakOverlay->GetXaxis()->SetTitle( "Timeslices" );
475 hMaxPeakOverlay->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
476 hList.Add( hMaxPeakOverlay );
477
478// hPixelPeakOverlay = new TH2F("hPixelPeakOverlay", "Maximum of Statistic of overlayed Peaks",
479// hPeakOverlayXaxisLeft + hPeakOverlayXaxisRight , (-1*hPeakOverlayXaxisLeft)-0.5, hPeakOverlayXaxisRight-0.5 ,
480// 512, -55.5, 200.5 );
481// hPixelPeakOverlay->GetXaxis()->SetTitle( "Timeslices" );
482// hPixelPeakOverlay->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
483// hList.Add( hPixelPeakOverlay );
484
485// hEventPeakOverlay = new TH2F("hEventPeakOverlay", "Overlay of detected Peaks of all Pixel of one Event",
486// hPeakOverlayXaxisLeft + hPeakOverlayXaxisRight , (-1*hPeakOverlayXaxisLeft)-0.5, hPeakOverlayXaxisRight-0.5 ,
487// 4096, -48.5, 200.5 );
488// hEventPeakOverlay->GetXaxis()->SetTitle( "Timeslices" );
489// hEventPeakOverlay->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
490// hList.Add( hEventPeakOverlay );
491
492}
493
494void SaveHistograms(const char * loc_fname ){
495 // create the histogram file (replace if already existing)
496 TFile tf( loc_fname, "RECREATE");
497
498 hList.Write(); // write the major histograms into the top level directory
499
500 tf.Close(); // close the file
501} // end of SaveHistograms( char * loc_fname )
502
Note: See TracBrowser for help on using the repository browser.