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

Last change on this file since 12556 was 12556, checked in by Jens Buss, 13 years ago
...now has most propable Single peak in one Histogram with Landau-fit
File size: 19.9 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
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{
129hPeakOverlayXaxisLeft = OverlayWindowLeft;
130hPeakOverlayXaxisRight = 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, 600);
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 = OpenDataFile( 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 FOpenCalibFile( drsfilename,
210 drs_basemean,
211 drs_gainmean,
212 drs_triggeroffsetmean,
213 drs_n);
214cout << "test4" << endl;
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
231 if (verbosityLevel > 0){
232 if (pix == firstpixel){
233 cout << "Processing Event: " << CurrentEventID << "/" << nevents << endl;
234 }
235 }
236
237 applyDrsCalibration( Ameas,pix,12,12,
238 drs_basemean, drs_gainmean, drs_triggeroffsetmean,
239 RegionOfInterest, AllPixelDataVector, StartCellVector);
240
241 // finds spikes in the raw data, and interpolates the value
242 // spikes are: 1 or 2 slice wide, positive non physical artifacts
243 removeSpikes (Ameas, Vcorr);
244
245 // filter Vcorr with sliding average using FIR filter function
246 sliding_avg(Vcorr, Vslide, avg1);
247 // filter Vslide with CFD using FIR filter function
248 factfir(b_cfd , a_cfd, k_cfd, Vslide, Vcfd);
249 // filter Vcfd with sliding average using FIR filter function
250 sliding_avg(Vcfd, Vcfd2, avg2);
251
252 // peaks in Ameas[] are found by searching for zero crossings
253 // in Vcfd2
254 // first Argument 1 means ... *rising* edge
255 // second Argument 1 means ... search with stepsize 1 ... 10 is okay as well
256 vector<Region> * zXings = zerosearch( Vcfd2 , 1 , 1);
257 // zXings means "zero cross ings"
258 EnlargeRegion(*zXings, 10, 10);
259 findAbsMaxInRegions(*zXings, Vslide);
260 removeMaximaBelow( *zXings, 3.0);
261 removeRegionWithMaxOnEdge( *zXings, 2);
262 removeRegionOnFallingEdge( *zXings, 100);
263
264 //following Code should produce the Overlay of peaks
265 vector<Region>::iterator it;
266
267 for (it = zXings->begin() ; it < zXings->end() ; it++){
268 int Left = it->maxPos - OverlayWindowLeft;
269 int Right = it->maxPos + OverlayWindowRight;
270 if (Left < 0)
271 Left =0;
272 if (Right > (int)Vcorr.size() )
273 Right=Vcorr.size() ;
274 for ( int pos = Left; pos < Right; pos++){
275
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
285
286 // TODO do this correct. The vectors should be the rigt ones... this is just luck
287 debugHistos[Ameas_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
288 debugHistos[Vcorr_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
289 debugHistos[Vslide_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
290 debugHistos[Vcfd_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
291 debugHistos[Vcfd2_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5);
292
293 for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
294 debugHistos[Ameas_].SetBinContent(sl, Ameas[sl]);
295 debugHistos[Vcorr_].SetBinContent(sl, Vcorr[sl]);
296 debugHistos[Vslide_].SetBinContent( sl, Vslide[sl] );
297 debugHistos[Vcfd_].SetBinContent( sl, Vcfd[sl] );
298 debugHistos[Vcfd2_].SetBinContent( sl, Vcfd2[sl] );
299 }
300
301
302 cFiltered->cd(1);
303 gPad->SetGrid();
304 debugHistos[Vcorr_].Draw();
305
306 cFiltered->cd(2);
307 gPad->SetGrid();
308 debugHistos[Vslide_].Draw();
309
310 TBox *OneBox;
311 vector<TBox*> MyBoxes;
312 for (unsigned int i=0; i<zXings->size(); i++){
313 OneBox = new TBox(
314 zXings->at(i).maxPos -10 ,
315 zXings->at(i).maxVal -0.5,
316 zXings->at(i).maxPos +10 ,
317 zXings->at(i).maxVal +0.5);
318 OneBox->SetLineColor(kBlue);
319 OneBox->SetLineWidth(1);
320 OneBox->SetFillStyle(0);
321 OneBox->SetFillColor(kRed);
322 MyBoxes.push_back(OneBox);
323 OneBox->Draw();
324 }
325
326 cFiltered->cd(3);
327 gPad->SetGrid();
328 debugHistos[Vcfd2_].Draw();
329 TLine *zeroline = new TLine(0, 0, 1024, 0);
330 zeroline->SetLineColor(kBlue);
331 zeroline->Draw();
332
333 cFiltered->Update();
334
335 //OVERLAY PEAKS
336 cPeakOverlay->cd(1);
337 hPeakOverlay->Draw("COLZ");
338 cPeakOverlay->Modified();
339 cPeakOverlay->Update();
340
341 //Process gui events asynchronously during input
342 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
343 timer.TurnOn();
344 TString input = Getline("Type 'q' to exit, <return> to go on: ");
345 timer.TurnOff();
346 if (input=="q\n") {
347 breakout=true;
348 }
349
350 //TODO!!!!!!!!!
351 // do some Garbage collection ...
352 // all the Objects on the heap should be deleted here.
353
354 }// end of if(spikeDebug)
355
356 delete zXings;
357 if (breakout)
358 break;
359
360 }
361 //-------------------------------------
362 // end of loop over pixels
363 //-------------------------------------
364 if (ProduceGraphic){
365 if (ev % 50 == 0){
366
367 //OVERLAY PEAKS
368 cPeakOverlay->cd(1);
369 hPeakOverlay->Draw("COLZ");
370 cPeakOverlay->cd(2);
371 // hMaxPeakOverlay->Draw("");
372 cPeakOverlay->Modified();
373 cPeakOverlay->Update();
374 }
375 }
376 if (breakout)
377 break;
378 }
379
380//-------------------------------------
381// end of loop over Events
382//-------------------------------------
383
384//-------------------------------------
385// Histogramm of Maximas in Overlay Spectra
386//-------------------------------------
387 // resize vector SingleMaximum to number of timeslices in Overlay Spectra
388 SingleMaximum.resize(hPeakOverlay->GetNbinsX());
389
390 //put maximumvalue of every Y-projection of every timeslice into vector
391 for (int TimeSlice = 0; TimeSlice <= hPeakOverlay->GetNbinsX(); TimeSlice++ ){
392 TH1D* hProjPeak = hPeakOverlay->ProjectionY("proj_py", TimeSlice, TimeSlice, "");
393 SingleMaximum[ TimeSlice ].maxAmpl = (hProjPeak->GetMaximumBin() * 0.5) - 54;
394 SingleMaximum[ TimeSlice ].countOfMax = hProjPeak->GetBinContent( hProjPeak->GetMaximumBin() );
395 hMaxPeakOverlay->SetBinContent(TimeSlice, SingleMaximum[ TimeSlice ].maxAmpl );
396 }
397 hMaxPeakOverlay->Fit("landau", "", "", -50, 250);
398
399//-------------------------------------
400// Draw Histograms
401//-------------------------------------
402 if (ProduceGraphic){
403 //OVERLAY PEAKS
404 cPeakOverlay->cd(1);
405 hPeakOverlay->Draw("COLZ");
406 cPeakOverlay->cd(2);
407 hMaxPeakOverlay->Draw("");
408 cPeakOverlay->Modified();
409 cPeakOverlay->Update();
410
411// cPixelPeakOverlay->cd();
412// hPixelPeakOverlay->Draw("COLZ");
413// cPixelPeakOverlay->Modified();
414// cPixelPeakOverlay->Update();
415
416 }
417 SaveHistograms( OutRootFileName );
418
419 return( 0 );
420}
421
422// booking and parameter settings for all histos
423void BookHistos( ){
424
425 debugHistos = new TH1F[ NumberOfDebugHistoTypes ];
426 for ( int type = 0; type < NumberOfDebugHistoTypes; type++){
427 debugHistos[ type ].SetBins(1024, 0, 1024);
428 debugHistos[ type ].SetLineColor(1);
429 debugHistos[ type ].SetLineWidth(2);
430
431 // set X axis paras
432 debugHistos[ type ].GetXaxis()->SetLabelSize(0.1);
433 debugHistos[ type ].GetXaxis()->SetTitleSize(0.1);
434 debugHistos[ type ].GetXaxis()->SetTitleOffset(1.2);
435 debugHistos[ type ].GetXaxis()->SetTitle(Form("Time slice (%.1f ns/slice)", 1./2.));
436
437 // set Y axis paras
438 debugHistos[ type ].GetYaxis()->SetLabelSize(0.1);
439 debugHistos[ type ].GetYaxis()->SetTitleSize(0.1);
440 debugHistos[ type ].GetYaxis()->SetTitleOffset(0.3);
441 debugHistos[ type ].GetYaxis()->SetTitle("Amplitude (a.u.)");
442 }
443/*
444 ProjHistos = new TH1D[ NumberOfTimeSlices ];
445 for ( int type = 0; type < NumberOfTimeSlices; type++){
446 ProjHistos[ type ].SetBins(1024, 0, 1024);
447 ProjHistos[ type ].SetLineColor(1);
448 ProjHistos[ type ].SetLineWidth(2);
449
450 // set X axis paras
451 ProjHistos[ type ].GetXaxis()->SetLabelSize(0.1);
452 ProjHistos[ type ].GetXaxis()->SetTitleSize(0.1);
453 ProjHistos[ type ].GetXaxis()->SetTitleOffset(1.2);
454 ProjHistos[ type ].GetXaxis()->SetTitle(Form("Time slice (%.1f ns/slice)", 1./2.));
455
456 // set Y axis paras
457 ProjHistos[ type ].GetYaxis()->SetLabelSize(0.1);
458 ProjHistos[ type ].GetYaxis()->SetTitleSize(0.1);
459 ProjHistos[ type ].GetYaxis()->SetTitleOffset(0.3);
460 ProjHistos[ type ].GetYaxis()->SetTitle("Amplitude (a.u.)");
461 }
462*/
463 hPeakOverlay = new TH2F("hPeakOverlay", "Overlay of detected Peaks",
464 hPeakOverlayXaxisLeft + hPeakOverlayXaxisRight , (-1*hPeakOverlayXaxisLeft)-0.5, hPeakOverlayXaxisRight-0.5 ,
465 512, -55.5, 200.5 );
466 hPeakOverlay->GetXaxis()->SetTitle( "Timeslices" );
467 hPeakOverlay->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
468 //hPeakOverlay->SetBit(TH2F::kCanRebin);
469 hList.Add( hPeakOverlay );
470
471 hMaxPeakOverlay = new TH1F("hMaxPeakOverlay", "Single Peak derived with Maximum of above Spektrum",
472 hPeakOverlayXaxisLeft + hPeakOverlayXaxisRight ,
473 (-1*hPeakOverlayXaxisLeft)-0.5,
474 hPeakOverlayXaxisRight-0.5
475 );
476 hMaxPeakOverlay->GetXaxis()->SetTitle( "Timeslices" );
477 hMaxPeakOverlay->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
478 hList.Add( hMaxPeakOverlay );
479
480// hPixelPeakOverlay = new TH2F("hPixelPeakOverlay", "Maximum of Statistic of overlayed Peaks",
481// hPeakOverlayXaxisLeft + hPeakOverlayXaxisRight , (-1*hPeakOverlayXaxisLeft)-0.5, hPeakOverlayXaxisRight-0.5 ,
482// 512, -55.5, 200.5 );
483// hPixelPeakOverlay->GetXaxis()->SetTitle( "Timeslices" );
484// hPixelPeakOverlay->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
485// hList.Add( hPixelPeakOverlay );
486
487// hEventPeakOverlay = new TH2F("hEventPeakOverlay", "Overlay of detected Peaks of all Pixel of one Event",
488// hPeakOverlayXaxisLeft + hPeakOverlayXaxisRight , (-1*hPeakOverlayXaxisLeft)-0.5, hPeakOverlayXaxisRight-0.5 ,
489// 4096, -48.5, 200.5 );
490// hEventPeakOverlay->GetXaxis()->SetTitle( "Timeslices" );
491// hEventPeakOverlay->GetYaxis()->SetTitle( "Amplitude [a.u.]" );
492// hList.Add( hEventPeakOverlay );
493
494}
495
496void SaveHistograms(const char * loc_fname ){
497 // create the histogram file (replace if already existing)
498 TFile tf( loc_fname, "RECREATE");
499
500 hList.Write(); // write the major histograms into the top level directory
501
502 tf.Close(); // close the file
503} // end of SaveHistograms( char * loc_fname )
504
Note: See TracBrowser for help on using the repository browser.