source: fact/tools/rootmacros/flightpulser.C

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