source: fact/tools/rootmacros/flightpulser.C@ 12515

Last change on this file since 12515 was 12391, checked in by neise, 14 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.