source: fact/tools/rootmacros/fbsl.C@ 13521

Last change on this file since 13521 was 12527, checked in by neise, 14 years ago
faster abort in case of wrong ROI calib file
File size: 8.8 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 <TMath.h>
10#include <TFile.h>
11
12#include <stdio.h>
13#include <stdint.h>
14#include <cstdio>
15
16#define HAVE_ZLIB
17#include "fits.h"
18
19#include "openFits.h"
20#include "openFits.c"
21
22#include "DrsCalibration.h"
23#include "DrsCalibration.C"
24
25#include "SpikeRemoval.h"
26#include "SpikeRemoval.C"
27
28//#include "zerosearch.C"
29#include "factfir.C"
30
31#define NPIX 1440
32#define NCELL 1024
33
34// data access and treatment
35#define FAD_MAX_SAMPLES 1024
36
37int NEvents;
38vector<int16_t> Data; // vector, which will be filled with raw data
39vector<int16_t> StartCells; // vector, which will be filled with DRS start positions
40unsigned int EventID; // index of the current event
41UInt_t RegionOfInterest; // Width of the Region, read out of the DRS
42UInt_t NumberOfPixels; // Total number of pixel, read out of the camera
43size_t PXLxROI; // Size of column "Data" = #Pixel x ROI
44
45int NBSLeve = 1000;
46
47size_t TriggerOffsetROI, RC;
48vector<float> Offset, Gain, TriggerOffset;
49
50
51vector<float> Ameas(FAD_MAX_SAMPLES); // copy of the data (measured amplitude
52vector<float> N1mean(FAD_MAX_SAMPLES); // mean of the +1 -1 ch neighbors
53vector<float> N2mean(FAD_MAX_SAMPLES); // mean of the +2 -2 ch neighbors
54vector<float> Vcorr(FAD_MAX_SAMPLES); // corrected Values
55vector<float> Vdiff(FAD_MAX_SAMPLES); // numerical derivative
56
57vector<float> Vslide(FAD_MAX_SAMPLES); // sliding average result
58vector<float> Vcfd(FAD_MAX_SAMPLES); // CDF result
59vector<float> Vcfd2(FAD_MAX_SAMPLES); // CDF result + 2nd sliding average
60
61
62
63// histograms
64const int Ntypes = 7;
65const unsigned int // arranged by Dominik
66 tAmeas = 0,
67 tN1mean = 1,
68 tVcorr = 2,
69 tVtest = 3,
70 tVslide = 4,
71 tVcfd = 5,
72 tVcfd2 = 6;
73
74TH1F* h;
75TH2F hPixelCellData(
76 "PixelPedestal", "PixelPedestal", NCELL, 0., NCELL, 200, -50., 150.);
77TH1F *hBaseline[ NPIX ]; // histograms for baseline extraction
78TH1F *hMeanBsl, *hpltMeanBsl;
79TH1F *hRmsBsl, *hpltRmsBsl;
80TObjArray hList;
81TObjArray hListBaseline;
82
83void BookHistos( int , int);
84void SaveHistograms( const char * );
85
86// Create a canvas
87TCanvas* CW;
88TCanvas* cFilter;
89
90
91int fbsl(
92 const char *datafilename = "path-to-datafile.fits.gz",
93 const char *drsfilename = "path-to-calibfile.drs.fits.gz",
94 int firstevent = 0,
95 int nevents = -1,
96 int firstpixel = 0,
97 int npixel = -1,
98 const char *TextOutFileName = "",
99 const char *RootOutFileName = "",
100 bool produceGraphic = true
101){
102 fits * datafile = NULL;
103
104 NEvents = openDataFits(
105 datafilename,
106 &datafile,
107 Data,
108 StartCells,
109 EventID,
110 RegionOfInterest,
111 NumberOfPixels,
112 PXLxROI);
113
114 printf("number of events in file: %d\n", NEvents);
115 if (NEvents == 0){
116 cout << "return code of openDataFits:" << datafilename<< endl;
117 cout << "is zero -> aborting." << endl;
118 return 1;
119 }
120
121
122 // compare the number of events in the data file with the nevents the
123 // the user would like to read. nevents = -1 means: read all
124 if ( nevents == -1 || nevents > NEvents ) nevents = NEvents;
125 if ( npixel == -1 || npixel > (int)NumberOfPixels) npixel = NumberOfPixels;
126
127 RC = openCalibFits( drsfilename, Offset, Gain, TriggerOffset, TriggerOffsetROI);
128 if (RC == 0){
129 cout << "return code of openCalibFits:" << drsfilename << endl;
130 cout << "is zero -> aborting." << endl;
131 return 1;
132 }
133
134
135 BookHistos( RegionOfInterest, npixel );
136 size_t calib_RC = 1;
137
138 // loop over events
139 for ( int ev = firstevent; ev < firstevent + nevents; ev++) {
140
141 datafile->GetRow( ev );
142
143 if( ev % 100 == 0){
144 cout << "Event ID: " << EventID << endl;
145 }
146
147 // loop over pixel
148 for ( int pix = firstpixel ; pix < npixel+firstpixel ; pix++ ){
149
150
151 // get the data of this pixel from the Data vector
152 // apply the Drs Calibration and cut off 12 slices at the beginning
153 // and at the end.
154 calib_RC = applyDrsCalibration( Ameas,pix,12,12,
155 Offset, Gain, TriggerOffset,
156 RegionOfInterest, Data, StartCells);
157 if (calib_RC == 0){
158 break;
159 }
160
161 // finds spikes in the raw data, and interpolates the value
162 // spikes are: 1 or 2 slice wide, positive non physical artifacts
163 removeSpikes (Ameas, Vcorr);
164
165 // filter Vcorr with sliding average using FIR filter function
166 // 8 is here the HalfWidth of the sliding average window.
167 sliding_avg(Vcorr, Vslide, 8);
168
169 for ( unsigned int sl = 0; sl < Vslide.size() ; sl++){
170 // hPixelCellData.Fill( sl, Vcorr[sl] );
171 hBaseline[pix-firstpixel]->Fill( Vslide[sl] ) ;
172 }
173 }
174
175 if (calib_RC == 0){
176 break;
177 }
178
179 }
180
181
182 if (calib_RC == 0){
183 cout << "DRS Calibration didn't work ... aborting" << endl;
184 return 2;
185 }
186
187 FILE *fp;
188 TString fName;
189
190 fName = TextOutFileName;
191 if (fName.Length() != 0){
192 cout << "saving text to " << TextOutFileName << endl;
193
194
195 fp = fopen(fName, "a+");
196 fprintf( fp, "FILE: %s\n", datafilename );
197 fprintf( fp, "NEVENTS: %d\n", nevents);
198 NBSLeve = nevents; // this has to be changed when computation on a subset of a run is implemented
199 fprintf( fp, "NBSLEVE: %d\n", NBSLeve );
200
201 for (int pix = firstpixel; pix < firstpixel+npixel; pix++) {
202 float vmaxprob = hBaseline[pix-firstpixel]->GetXaxis()->GetBinCenter(
203 hBaseline[pix-firstpixel]->GetMaximumBin() );
204
205 fprintf( fp, "%8.3f", vmaxprob );
206
207 hMeanBsl->Fill( vmaxprob );
208 hpltMeanBsl->SetBinContent(pix+1, vmaxprob );
209
210 hRmsBsl->Fill(hBaseline[pix-firstpixel]->GetRMS() );
211 hpltRmsBsl->SetBinContent( pix+1, hBaseline[pix]->GetRMS() );
212 }
213 fprintf( fp, "\n" );
214
215 fclose( fp );
216 }
217 if (RootOutFileName != ""){
218 cout << "saving histograms to " << RootOutFileName << endl;
219 SaveHistograms( RootOutFileName );
220 }
221 if (produceGraphic){
222 TCanvas * cMeanBsl = new TCanvas();
223 cMeanBsl->cd();
224 hMeanBsl->Draw();
225 cMeanBsl->Update();
226
227 TCanvas * cRmsBsl = new TCanvas();
228 cRmsBsl->cd();
229 hRmsBsl->Draw();
230 cMeanBsl->Update();
231 }
232 return( 0 );
233}
234
235void BookHistos( int Samples , int NumberOfPixel){
236// booking and parameter settings for all histos
237
238 // histograms for baseline extraction
239 char hName[500];
240 char hTitle[500];
241
242 TH1F *h;
243
244 printf("inside BookHistos\n");
245
246 for( int i = 0; i < NumberOfPixel; i++ ) {
247
248 // printf("call sprintf [%d]\n", i );
249 sprintf(&hTitle[0],"all events all slices of pixel %d", i);
250 sprintf(&hName[0],"base%d", i);
251 // printf("call sprintf [%d] done\n", i );
252
253 h = new TH1F( hName, hTitle, 400, -99.5 ,100.5 );
254
255 // printf("histo booked\n");
256 h->GetXaxis()->SetTitle( "Sample value (mV)" );
257 h->GetYaxis()->SetTitle( "Entries / 0.5 mV" );
258 // printf("histo title set\n");
259 hListBaseline.Add( h );
260 // printf("histo added to List\n");
261 hBaseline[i] = h;
262 // printf("histo assigned to array\n");
263 }
264
265 printf("made HBaseline * 1440\n");
266
267 hMeanBsl = new TH1F("histo_mean","Value of maximal probability",400,-99.5,100.5);
268 hMeanBsl->GetXaxis()->SetTitle( "max value (mV)" );
269 hMeanBsl->GetYaxis()->SetTitle( "Entries / 0.5 mV" );
270 hList.Add( hMeanBsl );
271
272 hpltMeanBsl = new TH1F("hplt_mean","Value of maximal probability",1440,-0.5,1439.5);
273 hpltMeanBsl->GetXaxis()->SetTitle( "pixel" );
274 hpltMeanBsl->GetYaxis()->SetTitle( "max value in mV" );
275 hList.Add( hpltMeanBsl );
276
277 hRmsBsl = new TH1F("histo_rms","RMS in mV",2000,-99.5,100.5);
278 hRmsBsl->GetXaxis()->SetTitle( "RMS (mV)" );
279 hRmsBsl->GetYaxis()->SetTitle( "Entries / 0.5 mV" );
280 hList.Add( hRmsBsl );
281
282 hpltRmsBsl = new TH1F("hplt_rms","Value of maximal probability",1440,-0.5,1439.5);
283 hpltRmsBsl->GetXaxis()->SetTitle( "pixel" );
284 hpltRmsBsl->GetYaxis()->SetTitle( "RMS in mV" );
285 hList.Add( hpltRmsBsl );
286}
287
288
289void SaveHistograms( const char * loc_fname ){
290
291 TString fName; // name of the histogram file
292
293 /* create the filename for the histogram file */
294 fName = loc_fname; // use the name of the tree file
295 //fName.Remove(fName.Length() - 5); // remove the extension .root
296 //fName += "_histo.root"; // add the new extension
297 //fName += ".root";
298
299 TFile tf( fName, "RECREATE"); // create the histogram file (replace if already existing)
300
301 hList.Write(); // write the major histograms into the top level directory
302 tf.mkdir("BaselineHisto"); tf.cd("BaselineHisto"); // go to new subdirectory
303 hListBaseline.Write(); // write histos into subdirectory
304
305 tf.Close(); // close the file
306
307} // end of function: void ana::SaveHistograms( void )
308
Note: See TracBrowser for help on using the repository browser.