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 |
|
---|
18 | #define NPIX 1440
|
---|
19 | #define NCELL 1024
|
---|
20 | #define FAD_MAX_SAMPLES 1024
|
---|
21 |
|
---|
22 | #define HAVE_ZLIB
|
---|
23 | #include "fits.h"
|
---|
24 | #include "FOpenCalibFile.c"
|
---|
25 |
|
---|
26 | #include "discriminator.h"
|
---|
27 | #include "discriminator.C"
|
---|
28 | #include "zerosearch.h"
|
---|
29 | #include "zerosearch.C"
|
---|
30 | #include "factfir.C"
|
---|
31 |
|
---|
32 |
|
---|
33 | vector<int16_t> AllPixelDataVector;
|
---|
34 | vector<int16_t> StartCellVector;
|
---|
35 |
|
---|
36 | unsigned int CurrentEventID;
|
---|
37 |
|
---|
38 | bool breakout=false;
|
---|
39 |
|
---|
40 | size_t ROIxNOP;
|
---|
41 | UInt_t NumberOfPixels;
|
---|
42 | UInt_t RegionOfInterest;
|
---|
43 | int NEvents;
|
---|
44 |
|
---|
45 | size_t drs_n;
|
---|
46 | vector<float> drs_basemean;
|
---|
47 | vector<float> drs_gainmean;
|
---|
48 | vector<float> drs_triggeroffsetmean;
|
---|
49 |
|
---|
50 | int FOpenDataFile( fits & );
|
---|
51 |
|
---|
52 |
|
---|
53 | vector<float> Ameas(FAD_MAX_SAMPLES); // copy of the data (measured amplitude
|
---|
54 | vector<float> N1mean(FAD_MAX_SAMPLES); // mean of the +1 -1 ch neighbors
|
---|
55 | vector<float> Vcorr(FAD_MAX_SAMPLES); // corrected Values
|
---|
56 | vector<float> Vdiff(FAD_MAX_SAMPLES); // numerical derivative
|
---|
57 |
|
---|
58 | vector<float> Vslide(FAD_MAX_SAMPLES); // sliding average result
|
---|
59 | vector<float> Vcfd(FAD_MAX_SAMPLES); // CDF result
|
---|
60 | vector<float> Vcfd2(FAD_MAX_SAMPLES); // CDF result + 2nd sliding average
|
---|
61 |
|
---|
62 | // not needed?
|
---|
63 | //vector<float> N2mean(FAD_MAX_SAMPLES); // mean of the +2 -2 ch neighbors
|
---|
64 |
|
---|
65 |
|
---|
66 | float getValue( int, int );
|
---|
67 | float correctDrsOffset( int slice, int pixel );
|
---|
68 | void computeN1mean( int );
|
---|
69 | void removeSpikes( int );
|
---|
70 |
|
---|
71 | // histograms
|
---|
72 | const int Ntypes = 7;
|
---|
73 | const unsigned int // arranged by Dominik
|
---|
74 | tAmeas = 0,
|
---|
75 | tN1mean = 1,
|
---|
76 | tVcorr = 2,
|
---|
77 | tVtest = 3,
|
---|
78 | tVslide = 4,
|
---|
79 | tVcfd = 5,
|
---|
80 | tVcfd2 = 6;
|
---|
81 |
|
---|
82 | TH1F* h;
|
---|
83 | TH1F *oldh;
|
---|
84 | TH2F* hStartCell; // id of the DRS physical pipeline cell where readout starts
|
---|
85 | // x = pixel id, y = DRS cell id
|
---|
86 | TH2F hPixelCellData(
|
---|
87 | "PixelPedestal", "PixelPedestal", NCELL, 0., NCELL, 200, -50., 150.);
|
---|
88 | TH1F *hBaseline[ NPIX ]; // histograms for baseline extraction
|
---|
89 | TH1F *hMeanBsl, *hpltMeanBsl;
|
---|
90 | TH1F *hRmsBsl, *hpltRmsBsl;
|
---|
91 | TH2F * hAmplSpek_cfd;
|
---|
92 | TH2F * hAmplSpek_discri;
|
---|
93 | TObjArray hList;
|
---|
94 | TObjArray hListBaseline;
|
---|
95 |
|
---|
96 | void BookHistos( );
|
---|
97 | void SaveHistograms( char * );
|
---|
98 |
|
---|
99 | // Create a canvas
|
---|
100 | TCanvas* CW;
|
---|
101 | TCanvas* cFilter;
|
---|
102 |
|
---|
103 |
|
---|
104 | int searchSinglesPeaks = 0;
|
---|
105 |
|
---|
106 | int fpeak_discri(
|
---|
107 | char *datafilename = "../../20111011_055.fits.gz",
|
---|
108 | const char *drsfilename = "../../20111011_054.drs.fits.gz",
|
---|
109 | int PixelID = -1,
|
---|
110 | int firstevent = 0,
|
---|
111 | int nevents = -1,
|
---|
112 | bool spikeDebug = false,
|
---|
113 | int verbosityLevel = 1 // different verbosity levels can be implemented here
|
---|
114 | )
|
---|
115 | {
|
---|
116 |
|
---|
117 | gStyle->SetPalette(1,0);
|
---|
118 | gROOT->SetStyle("Plain");
|
---|
119 | // read FACT raw data
|
---|
120 | // * remove spikes
|
---|
121 | // * calculate baseline
|
---|
122 | // * find peaks (CFD and normal discriminator)
|
---|
123 | // * compute pulse height and pulse integral spektrum of the peaks
|
---|
124 |
|
---|
125 | // sliding window filter settings
|
---|
126 | int k_slide = 16;
|
---|
127 | vector<double> a_slide(k_slide, 1);
|
---|
128 | double b_slide = k_slide;
|
---|
129 |
|
---|
130 | // CFD filter settings
|
---|
131 | int k_cfd = 10;
|
---|
132 | vector<double> a_cfd(k_cfd, 0);
|
---|
133 | double b_cfd = 1.;
|
---|
134 | a_cfd[0]=-0.75;
|
---|
135 | a_cfd[k_cfd-1]=1.;
|
---|
136 |
|
---|
137 | // 2nd slinding window filter
|
---|
138 | int ks2 = 16;
|
---|
139 | vector<double> as2(ks2, 1);
|
---|
140 | double bs2 = ks2;
|
---|
141 |
|
---|
142 | // Open the data file
|
---|
143 | fits *datafile = new fits( datafilename );
|
---|
144 | if (!datafile) {
|
---|
145 | printf( "Could not open the file: %s\n", datafilename );
|
---|
146 | return 1;
|
---|
147 | }
|
---|
148 |
|
---|
149 | // access data
|
---|
150 | NEvents = FOpenDataFile( *datafile );
|
---|
151 | printf("number of events in file: %d\n", NEvents);
|
---|
152 | if ( nevents == -1 || nevents > NEvents ) nevents = NEvents; // -1 means all!
|
---|
153 |
|
---|
154 | //Get the DRS calibration
|
---|
155 | FOpenCalibFile( drsfilename,
|
---|
156 | drs_basemean,
|
---|
157 | drs_gainmean,
|
---|
158 | drs_triggeroffsetmean,
|
---|
159 | drs_n);
|
---|
160 |
|
---|
161 | //Check the sizes of the data columns
|
---|
162 | if (drs_n != 1474560) {
|
---|
163 | cerr << "error: DRS calib file has wrong ...erm...size ... drs_n is: "
|
---|
164 | << drs_n << endl;
|
---|
165 | cerr << " Aborting." << endl;
|
---|
166 | return 1;
|
---|
167 | }
|
---|
168 |
|
---|
169 | if(ROIxNOP != 1474560)
|
---|
170 | {
|
---|
171 | cout << "warning: data_n should better be 1440x1024=1474560, but is "
|
---|
172 | << ROIxNOP << endl;
|
---|
173 | cout << "this script is not guaranteed to run under these "
|
---|
174 | <<" circumstances....any way ... it is never guaranteed." << endl;
|
---|
175 | }
|
---|
176 | // Book the histograms
|
---|
177 |
|
---|
178 | BookHistos( );
|
---|
179 | TCanvas * cSpektrum = new TCanvas();
|
---|
180 | cSpektrum->Divide(1,2);
|
---|
181 | cSpektrum->cd(1);
|
---|
182 | hAmplSpek_discri->Draw("COLZ");
|
---|
183 | cSpektrum->cd(2);
|
---|
184 | hAmplSpek_cfd->Draw("COLZ");
|
---|
185 | // TCanvas * cStartCell = new TCanvas();
|
---|
186 | // cStartCell->cd();
|
---|
187 | // hStartCell->Draw();
|
---|
188 | // hPixelCellData.Draw();
|
---|
189 |
|
---|
190 | float calibratedVoltage;
|
---|
191 | for ( int ev = firstevent; ev < firstevent + nevents; ev++) {
|
---|
192 | datafile->GetRow( ev );
|
---|
193 |
|
---|
194 |
|
---|
195 | for ( int pix = 0; pix < 1440; pix++ ){
|
---|
196 | // This means: GetEvent
|
---|
197 |
|
---|
198 | // this is a stupid hack ... these is more code at the
|
---|
199 | // end of this loop to complete this hack ...
|
---|
200 | if (PixelID != -1) {
|
---|
201 | pix = PixelID;
|
---|
202 | if (verbosityLevel > 0){
|
---|
203 | cout << "Processing Event number: " << CurrentEventID << "\t"
|
---|
204 | << "Pixel number: "<< pix << endl;
|
---|
205 |
|
---|
206 | }
|
---|
207 | }
|
---|
208 |
|
---|
209 | if (verbosityLevel > 0){
|
---|
210 | if (pix % 20 ==0){
|
---|
211 | cout << "Processing Event number: " << CurrentEventID << "\t"
|
---|
212 | << "Pixel number: "<< pix << endl;
|
---|
213 | }
|
---|
214 | }
|
---|
215 |
|
---|
216 | //hStartCell->Fill( pix, StartCellVector[pix] );
|
---|
217 |
|
---|
218 | for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
|
---|
219 | //calibratedVoltage = getValue( sl, pix);
|
---|
220 | calibratedVoltage = correctDrsOffset( sl, pix);
|
---|
221 | if (verbosityLevel > 10){
|
---|
222 | printf("calibratedVoltage = %f\n", calibratedVoltage);
|
---|
223 | }
|
---|
224 | Ameas[ sl ] = calibratedVoltage;
|
---|
225 | if (spikeDebug){
|
---|
226 | oldh[tAmeas].SetBinContent(sl, calibratedVoltage);
|
---|
227 | }
|
---|
228 | }
|
---|
229 |
|
---|
230 | computeN1mean( RegionOfInterest );
|
---|
231 | // operates on Ameas[] and writes to N1mean[]
|
---|
232 | removeSpikes( RegionOfInterest );
|
---|
233 | // operates on Ameas[] and N1mean[], then writes to Vcorr[]
|
---|
234 |
|
---|
235 | if (spikeDebug){
|
---|
236 | for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
|
---|
237 | hPixelCellData.Fill(sl, Vcorr[ sl ]);
|
---|
238 | }
|
---|
239 | }
|
---|
240 |
|
---|
241 | // filter Vcorr with sliding average using FIR filter function
|
---|
242 | factfir(b_slide , a_slide, k_slide, Vcorr, Vslide);
|
---|
243 | // filter Vslide with CFD using FIR filter function
|
---|
244 | factfir(b_cfd , a_cfd, k_cfd, Vslide, Vcfd);
|
---|
245 | // filter Vcfd with sliding average using FIR filter function
|
---|
246 | factfir(bs2 , as2, ks2, Vcfd, Vcfd2);
|
---|
247 |
|
---|
248 | float myTHR = 3.5;
|
---|
249 | vector<Region> *Regions = discriminator( Vslide, myTHR,true , 120);
|
---|
250 | for (unsigned int p=0; p<Regions->size(); p++ ){
|
---|
251 | hAmplSpek_discri->Fill(pix , Regions->at(p).maxVal);
|
---|
252 | }
|
---|
253 |
|
---|
254 | // peaks in Ameas[] are found by searching for zero crossings
|
---|
255 | // in Vcfd2
|
---|
256 | // first Argument 1 means ... *rising* edge
|
---|
257 | // second Argument 1 means ... search with stepsize 1 ... 10 is okay as well
|
---|
258 | vector<Region> * zXings = zerosearch( Vcfd2 , 1 , 1);
|
---|
259 |
|
---|
260 | // zXings means "zero cross ings"
|
---|
261 | ShiftRegionBy(*zXings, -ks2/2);
|
---|
262 | EnlargeRegion(*zXings, 10, 10);
|
---|
263 | findAbsMaxInRegions(*zXings, Vslide);
|
---|
264 |
|
---|
265 | if (zXings->size() != 0 ){
|
---|
266 | for (unsigned int i=0; i<zXings->size(); i++){
|
---|
267 | if (verbosityLevel > 1){
|
---|
268 | cout << zXings->at(i).maxPos << ":\t"
|
---|
269 | << zXings->at(i).maxVal <<endl;
|
---|
270 | }
|
---|
271 | hAmplSpek_cfd->Fill(pix, zXings->at(i).maxVal);
|
---|
272 | }
|
---|
273 | }
|
---|
274 |
|
---|
275 | if ( spikeDebug ){
|
---|
276 | for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
|
---|
277 | oldh[tVslide].SetBinContent( sl, Vslide[sl] );
|
---|
278 | oldh[tVcfd].SetBinContent( sl, Vcfd[sl] );
|
---|
279 | oldh[tVcfd2].SetBinContent( sl, Vcfd2[sl] );
|
---|
280 | }
|
---|
281 |
|
---|
282 | CW->cd( 1);
|
---|
283 | gPad->SetGrid();
|
---|
284 | oldh[tAmeas].Draw();
|
---|
285 |
|
---|
286 | // CW->cd( tN1mean + 1);
|
---|
287 | // gPad->SetGrid();
|
---|
288 | // oldh[tN1mean].Draw();
|
---|
289 |
|
---|
290 | CW->cd( 2);
|
---|
291 | gPad->SetGrid();
|
---|
292 | oldh[tVcorr].Draw();
|
---|
293 |
|
---|
294 | //cFilter->cd( Ntypes - tVslide );
|
---|
295 | cFilter->cd(1);
|
---|
296 | gPad->SetGrid();
|
---|
297 | oldh[tVslide].Draw();
|
---|
298 | TLine thrLine(0, myTHR, 1024, myTHR);
|
---|
299 | thrLine.SetLineColor(kRed);
|
---|
300 | thrLine.SetLineWidth(2);
|
---|
301 | thrLine.Draw();
|
---|
302 | TLine * OneLine;
|
---|
303 | vector<TLine*> MyLines;
|
---|
304 | for (unsigned int p=0; p<Regions->size(); p++ ){
|
---|
305 | OneLine = new TLine(
|
---|
306 | Regions->at(p).begin ,
|
---|
307 | Regions->at(p).maxVal,
|
---|
308 | Regions->at(p).end,
|
---|
309 | Regions->at(p).maxVal);
|
---|
310 | OneLine->SetLineColor(kRed);
|
---|
311 | OneLine->SetLineWidth(2);
|
---|
312 | MyLines.push_back(OneLine);
|
---|
313 | OneLine->Draw();
|
---|
314 | }
|
---|
315 | TBox *OneBox;
|
---|
316 | vector<TBox*> MyBoxes;
|
---|
317 | for (unsigned int i=0; i<zXings->size(); i++){
|
---|
318 | OneBox = new TBox(
|
---|
319 | zXings->at(i).maxPos -10 ,
|
---|
320 | zXings->at(i).maxVal -0.5,
|
---|
321 | zXings->at(i).maxPos +10 ,
|
---|
322 | zXings->at(i).maxVal +0.5);
|
---|
323 | OneBox->SetLineColor(kBlue);
|
---|
324 | OneBox->SetLineWidth(1);
|
---|
325 | OneBox->SetFillStyle(0);
|
---|
326 | OneBox->SetFillColor(kRed);
|
---|
327 | MyBoxes.push_back(OneBox);
|
---|
328 | OneBox->Draw();
|
---|
329 | }
|
---|
330 |
|
---|
331 | // cFilter->cd( Ntypes - tVcfd );
|
---|
332 | // cFilter->cd(2);
|
---|
333 | // gPad->SetGrid();
|
---|
334 | // oldh[tVcfd].Draw();
|
---|
335 | // TLine zeroline(0, 0, 1024, 0);
|
---|
336 | // zeroline.SetLineColor(kBlue);
|
---|
337 | // zeroline.Draw();
|
---|
338 |
|
---|
339 | //cFilter->cd( Ntypes - tVcfd2 );
|
---|
340 | cFilter->cd(2);
|
---|
341 | gPad->SetGrid();
|
---|
342 | oldh[tVcfd2].Draw();
|
---|
343 | TLine zeroline(0, 0, 1024, 0);
|
---|
344 | zeroline.SetLineColor(kBlue);
|
---|
345 | zeroline.Draw();
|
---|
346 |
|
---|
347 |
|
---|
348 |
|
---|
349 |
|
---|
350 |
|
---|
351 | CW->Update();
|
---|
352 | cFilter->Update();
|
---|
353 |
|
---|
354 | //Process gui events asynchronously during input
|
---|
355 | TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
|
---|
356 | timer.TurnOn();
|
---|
357 | TString input = Getline("Type 'q' to exit, <return> to go on: ");
|
---|
358 | timer.TurnOff();
|
---|
359 | if (input=="q\n") {
|
---|
360 | breakout=true;
|
---|
361 | delete OneLine;
|
---|
362 | // for (unsigned int i=0; i<MyLines.size(); i++){
|
---|
363 | // delete MyLines[i];
|
---|
364 | // }
|
---|
365 | // for (unsigned int i=0; i<MyBoxes.size(); i++){
|
---|
366 | // delete MyBoxes[i];
|
---|
367 | // }
|
---|
368 | }
|
---|
369 | }// end of if(spikeDebug)
|
---|
370 | delete Regions;
|
---|
371 | delete zXings;
|
---|
372 |
|
---|
373 | if (PixelID != -1){
|
---|
374 | pix = 2000;
|
---|
375 | }
|
---|
376 |
|
---|
377 | if (breakout)
|
---|
378 | break;
|
---|
379 |
|
---|
380 | } // end of loop over pixels
|
---|
381 |
|
---|
382 | if (ev % 10 ==0){
|
---|
383 | cSpektrum->Modified();
|
---|
384 | cSpektrum->Update();
|
---|
385 | // cStartCell->Modified();
|
---|
386 | // cStartCell->Update();
|
---|
387 | }
|
---|
388 | if (breakout)
|
---|
389 | break;
|
---|
390 | } // end of loop over pixels
|
---|
391 |
|
---|
392 |
|
---|
393 | // delete cStartCell;
|
---|
394 | return( 0 );
|
---|
395 | }
|
---|
396 |
|
---|
397 | void removeSpikes(int Samples){
|
---|
398 |
|
---|
399 | const float fract = 0.8;
|
---|
400 | float x, xp, xpp, x3p;
|
---|
401 |
|
---|
402 | // assume that there are no spikes
|
---|
403 | for ( int i = 0; i < Samples; i++) Vcorr[i] = Ameas[i];
|
---|
404 |
|
---|
405 | // find the spike and replace it by mean value of neighbors
|
---|
406 | for ( int i = 0; i < Samples; i++) {
|
---|
407 |
|
---|
408 | // printf("Vcorr[%d] = %f, Ameas[%d] = %f\n", i, Vcorr[ i ], i, Ameas[ i ] );
|
---|
409 |
|
---|
410 | x = Ameas[i] - N1mean[i];
|
---|
411 |
|
---|
412 | if ( x < -5. ){ // a spike candidate
|
---|
413 | // check consistency with a single channel spike
|
---|
414 | xp = Ameas[i+1] - N1mean[i+1];
|
---|
415 | xpp = Ameas[i+2] - N1mean[i+2];
|
---|
416 | x3p = Ameas[i+3] - N1mean[i+3];
|
---|
417 |
|
---|
418 | // printf("candidates x[%d] = %f; xp = %f; xpp = %f, x3p = %f\n", i, x, xp, xpp, x3p);
|
---|
419 |
|
---|
420 | if ( Ameas[i+2] - ( Ameas[i] + Ameas[i+3] )/2. > 10. ){
|
---|
421 | // printf("double spike candidate\n");
|
---|
422 | Vcorr[i+1] = ( Ameas[i] + Ameas[i+3] )/2.;
|
---|
423 | Vcorr[i+2] = ( Ameas[i] + Ameas[i+3] )/2.;
|
---|
424 | // printf("Vcorr[%d] = %f %f %f %f\n", i, Vcorr[i], Vcorr[i+1], Vcorr[i+2], Vcorr[ i+3 ]);
|
---|
425 | // printf("Ameas[%d] = %f %f %f %f\n", i, Ameas[ i ], Ameas[ i+1 ], Ameas[ i+2 ], Ameas[i+3]);
|
---|
426 | i = i + 3;
|
---|
427 | }
|
---|
428 | else{
|
---|
429 |
|
---|
430 | if ( ( xp > -2.*x*fract ) && ( xpp < -10. ) ){
|
---|
431 | Vcorr[i+1] = N1mean[i+1];
|
---|
432 | // printf("Vcorr[%d] = %f %f %f\n", i, Vcorr[i], Vcorr[i+1], Vcorr[i+2]);
|
---|
433 | // N1mean[i+1] = (Ameas[i] - Ameas[i+2] / 2.);
|
---|
434 | N1mean[i+2] = (Ameas[i+1] - Ameas[i+3] / 2.);
|
---|
435 | i = i + 2;//do not care about the next sample it was the spike
|
---|
436 | }
|
---|
437 | // treatment for the end of the pipeline must be added !!!
|
---|
438 | }
|
---|
439 | }
|
---|
440 | else{
|
---|
441 | // do nothing
|
---|
442 | }
|
---|
443 | } // end of spike search and correction
|
---|
444 | for ( int i = 0; i < Samples; i++ ) oldh[ tVcorr ].SetBinContent( i, Vcorr[i] );
|
---|
445 | }
|
---|
446 | /*
|
---|
447 | void computeN1mean( int Samples ){
|
---|
448 | cout << "In compute N1mean" << endl;
|
---|
449 | // compute the mean of the left and right neighbors of a channel
|
---|
450 |
|
---|
451 | for( int i = 0; i < Samples; i++){
|
---|
452 | if (i == 0){ // use right sample es mean
|
---|
453 | N1mean[i] = Ameas[i+1];
|
---|
454 | }
|
---|
455 | else if ( i == Samples-1 ){ //use left sample as mean
|
---|
456 | N1mean[i] = Ameas[i-1];
|
---|
457 | }
|
---|
458 | else{
|
---|
459 | N1mean[i] = ( Ameas[i-1] + Ameas[i+1] ) / 2.;
|
---|
460 | }
|
---|
461 | h[tN1mean].SetBinContent(i, Ameas[i] - N1mean[i]);
|
---|
462 | }
|
---|
463 | } // end of computeN1mean computation
|
---|
464 | */
|
---|
465 |
|
---|
466 | void computeN1mean( int Samples ){
|
---|
467 | // compute the mean of the left and right neighbors of a channel
|
---|
468 |
|
---|
469 | for( int i = 2; i < Samples - 2; i++){
|
---|
470 | /* if (i == 0){ // use right sample es mean
|
---|
471 | N1mean[i] = Ameas[i+1];
|
---|
472 | }
|
---|
473 | else if ( i == Samples-1 ){ //use left sample as mean
|
---|
474 | N1mean[i] = Ameas[i-1];
|
---|
475 | }
|
---|
476 | else{
|
---|
477 | N1mean[i] = ( Ameas[i-1] + Ameas[i+1] ) / 2.;
|
---|
478 | }
|
---|
479 | */
|
---|
480 | N1mean[i] = ( Ameas[i-1] + Ameas[i+1] ) / 2.;
|
---|
481 | }
|
---|
482 | } // end of computeN1mean computation
|
---|
483 |
|
---|
484 | float getValue( int slice, int pixel ){
|
---|
485 | const float dconv = 2000/4096.0;
|
---|
486 |
|
---|
487 | float vraw, vcal;
|
---|
488 |
|
---|
489 | unsigned int pixel_pt;
|
---|
490 | unsigned int slice_pt;
|
---|
491 | unsigned int cal_pt;
|
---|
492 | unsigned int drs_cal_offset;
|
---|
493 |
|
---|
494 | // printf("pixel = %d, slice = %d\n", slice, pixel);
|
---|
495 |
|
---|
496 | pixel_pt = pixel * RegionOfInterest;
|
---|
497 | slice_pt = pixel_pt + slice;
|
---|
498 | drs_cal_offset = ( slice + StartCellVector[ pixel ] )%RegionOfInterest;
|
---|
499 | cal_pt = pixel_pt + drs_cal_offset;
|
---|
500 |
|
---|
501 | vraw = AllPixelDataVector[ slice_pt ] * dconv;
|
---|
502 | vcal = ( vraw - drs_basemean[ cal_pt ] - drs_triggeroffsetmean[ slice_pt ] ) / drs_gainmean[ cal_pt ]*1907.35;
|
---|
503 |
|
---|
504 | return( vcal );
|
---|
505 | }
|
---|
506 | float correctDrsOffset( int slice, int pixel ){
|
---|
507 |
|
---|
508 | const float dconv = 2000/4096.0;
|
---|
509 |
|
---|
510 | // here 1024 is not the RegionOfInterest, but really the lenth of the pipeline
|
---|
511 | unsigned int physical_slice = ( slice + StartCellVector[ pixel ] ) % 1024;
|
---|
512 |
|
---|
513 | unsigned int slice_pt;
|
---|
514 | unsigned int physical_slice_pt;
|
---|
515 | slice_pt = pixel * RegionOfInterest + slice;
|
---|
516 | physical_slice_pt = pixel * RegionOfInterest + physical_slice;
|
---|
517 |
|
---|
518 | float vcal = AllPixelDataVector[ slice_pt ] *
|
---|
519 | dconv - drs_basemean[ physical_slice_pt ];
|
---|
520 | return( vcal );
|
---|
521 | }
|
---|
522 |
|
---|
523 | void BookHistos( ){
|
---|
524 | // booking and parameter settings for all histos
|
---|
525 |
|
---|
526 | // histograms for baseline extraction
|
---|
527 | char hName[500];
|
---|
528 | char hTitle[500];
|
---|
529 |
|
---|
530 | TH1F *h;
|
---|
531 |
|
---|
532 | for( int i = 0; i < NPIX; i++ ) {
|
---|
533 | sprintf(&hTitle[0],"all events all slices of pixel %d", i);
|
---|
534 | sprintf(&hName[0],"base%d", i);
|
---|
535 |
|
---|
536 | h = new TH1F( hName, hTitle, 400, -99.5 ,100.5 );
|
---|
537 |
|
---|
538 | h->GetXaxis()->SetTitle( "Sample value (mV)" );
|
---|
539 | h->GetYaxis()->SetTitle( "Entries / 0.5 mV" );
|
---|
540 | hListBaseline.Add( h );
|
---|
541 | hBaseline[i] = h;
|
---|
542 | }
|
---|
543 |
|
---|
544 | hMeanBsl = new TH1F("histo_mean","Value of maximal probability",400,-99.5,100.5);
|
---|
545 | hMeanBsl->GetXaxis()->SetTitle( "max value (mV)" );
|
---|
546 | hMeanBsl->GetYaxis()->SetTitle( "Entries / 0.5 mV" );
|
---|
547 | hList.Add( hMeanBsl );
|
---|
548 |
|
---|
549 | hpltMeanBsl = new TH1F("hplt_mean","Value of maximal probability",1440,-0.5,1439.5);
|
---|
550 | hpltMeanBsl->GetXaxis()->SetTitle( "pixel" );
|
---|
551 | hpltMeanBsl->GetYaxis()->SetTitle( "max value in mV" );
|
---|
552 | hList.Add( hpltMeanBsl );
|
---|
553 |
|
---|
554 | hRmsBsl = new TH1F("histo_rms","RMS in mV",2000,-99.5,100.5);
|
---|
555 | hRmsBsl->GetXaxis()->SetTitle( "RMS (mV)" );
|
---|
556 | hRmsBsl->GetYaxis()->SetTitle( "Entries / 0.5 mV" );
|
---|
557 | hList.Add( hRmsBsl );
|
---|
558 |
|
---|
559 | hpltRmsBsl = new TH1F("hplt_rms","Value of maximal probability",1440,-0.5,1439.5);
|
---|
560 | hpltRmsBsl->GetXaxis()->SetTitle( "pixel" );
|
---|
561 | hpltRmsBsl->GetYaxis()->SetTitle( "RMS in mV" );
|
---|
562 | hList.Add( hpltRmsBsl );
|
---|
563 |
|
---|
564 | hAmplSpek_cfd = new TH2F("hAmplSpek_cfd","amplitude spektrum - CFD",1440,-0.5,1439.5, 256, -27.5, 100.5);
|
---|
565 | hAmplSpek_cfd->GetXaxis()->SetTitle( "pixel" );
|
---|
566 | hAmplSpek_cfd->GetYaxis()->SetTitle( "amplitude in mV" );
|
---|
567 | hList.Add( hAmplSpek_cfd );
|
---|
568 |
|
---|
569 | hAmplSpek_discri = new TH2F("hAmplSpek_discri","amplitude spektrum - std discriminator",1440,-0.5,1439.5, 256, -27.5, 100.5);
|
---|
570 | hAmplSpek_discri->GetXaxis()->SetTitle( "pixel" );
|
---|
571 | hAmplSpek_discri->GetXaxis()->SetTitle( "amplitude in mV" );
|
---|
572 | hList.Add( hAmplSpek_discri );
|
---|
573 |
|
---|
574 | hStartCell = new TH2F("StartCell", "StartCell", 1440, 0., 1440., 1024, 0., 1024);
|
---|
575 | hStartCell->GetXaxis()->SetTitle( "pixel" );
|
---|
576 | hStartCell->GetXaxis()->SetTitle( "slice" );
|
---|
577 | hList.Add( hStartCell );
|
---|
578 |
|
---|
579 | oldh = new TH1F[ Ntypes ];
|
---|
580 |
|
---|
581 | for ( int type = 0; type < Ntypes; type++){
|
---|
582 |
|
---|
583 | oldh[ type ].SetBins(1024, 0, 1024);
|
---|
584 | oldh[ type ].SetLineColor(1);
|
---|
585 | oldh[ type ].SetLineWidth(2);
|
---|
586 |
|
---|
587 | // set X axis paras
|
---|
588 | oldh[ type ].GetXaxis()->SetLabelSize(0.1);
|
---|
589 | oldh[ type ].GetXaxis()->SetTitleSize(0.1);
|
---|
590 | oldh[ type ].GetXaxis()->SetTitleOffset(1.2);
|
---|
591 | oldh[ type ].GetXaxis()->SetTitle(Form("Time slice (%.1f ns/slice)", 1./2.));
|
---|
592 |
|
---|
593 | // set Y axis paras
|
---|
594 | oldh[ type ].GetYaxis()->SetLabelSize(0.1);
|
---|
595 | oldh[ type ].GetYaxis()->SetTitleSize(0.1);
|
---|
596 | oldh[ type ].GetYaxis()->SetTitleOffset(0.3);
|
---|
597 | oldh[ type ].GetYaxis()->SetTitle("Amplitude (a.u.)");
|
---|
598 | }
|
---|
599 |
|
---|
600 | CW = new TCanvas("CW","DRS Waveform",10,10,800,600);
|
---|
601 | CW->Divide(1, 2);
|
---|
602 | cFilter = new TCanvas("cFilter","filtered DRS Waveforms",10,10,800,600);
|
---|
603 | cFilter->Divide(1, 2);
|
---|
604 | }
|
---|
605 |
|
---|
606 | void SaveHistograms( char * loc_fname ){
|
---|
607 |
|
---|
608 | TString fName; // name of the histogram file
|
---|
609 |
|
---|
610 | // create the filename for the histogram file
|
---|
611 | fName = loc_fname; // use the name of the tree file
|
---|
612 | // TODO ... next statement doesn't work for ".fits.gz"
|
---|
613 | fName.Remove(fName.Length() - 5); // remove the extension .fits
|
---|
614 | fName += "_discri.root"; // add the new extension
|
---|
615 |
|
---|
616 | // create the histogram file (replace if already existing)
|
---|
617 | TFile tf( fName, "RECREATE");
|
---|
618 |
|
---|
619 | hList.Write(); // write the major histograms into the top level directory
|
---|
620 | tf.mkdir("BaselineHisto");
|
---|
621 | tf.cd("BaselineHisto"); // go to new subdirectory
|
---|
622 | hListBaseline.Write(); // write histos into subdirectory
|
---|
623 |
|
---|
624 | tf.Close(); // close the file
|
---|
625 | } // end of SaveHistograms( char * loc_fname )
|
---|
626 |
|
---|
627 | int FOpenDataFile(fits &datafile){
|
---|
628 |
|
---|
629 | // cout << "-------------------- Data Header -------------------" << endl;
|
---|
630 | // datafile.PrintKeys();
|
---|
631 | // cout << "------------------- Data Columns -------------------" << endl;
|
---|
632 | // datafile.PrintColumns();
|
---|
633 |
|
---|
634 | //Get the size of the data column
|
---|
635 | RegionOfInterest = datafile.GetUInt("NROI");
|
---|
636 | NumberOfPixels = datafile.GetUInt("NPIX");
|
---|
637 | //Size of column "Data" = #Pixel x ROI
|
---|
638 | ROIxNOP = datafile.GetN("Data");
|
---|
639 |
|
---|
640 | //Set the sizes of the data vectors
|
---|
641 | AllPixelDataVector.resize(ROIxNOP,0);
|
---|
642 | StartCellVector.resize(NumberOfPixels,0);
|
---|
643 |
|
---|
644 | //Link the data to variables
|
---|
645 | datafile.SetRefAddress("EventNum", CurrentEventID);
|
---|
646 | datafile.SetVecAddress("Data", AllPixelDataVector);
|
---|
647 | datafile.SetVecAddress("StartCellData", StartCellVector);
|
---|
648 | datafile.GetRow(0);
|
---|
649 |
|
---|
650 | return datafile.GetNumRows() ;
|
---|
651 | }
|
---|
652 |
|
---|
653 |
|
---|
654 | /////////////////////////////////////////////////////////////////////////////
|
---|
655 | /// old BookHistos
|
---|
656 | /*
|
---|
657 | void BookHistos( int Samples ){
|
---|
658 | // booking and parameter settings for all histos
|
---|
659 |
|
---|
660 | h = new TH1F[ Ntypes ];
|
---|
661 |
|
---|
662 | for ( int type = 0; type < Ntypes; type++){
|
---|
663 |
|
---|
664 | h[ type ].SetBins(Samples, 0, Samples);
|
---|
665 | h[ type ].SetLineColor(1);
|
---|
666 | h[ type ].SetLineWidth(2);
|
---|
667 |
|
---|
668 | // set X axis paras
|
---|
669 | h[ type ].GetXaxis()->SetLabelSize(0.1);
|
---|
670 | h[ type ].GetXaxis()->SetTitleSize(0.1);
|
---|
671 | h[ type ].GetXaxis()->SetTitleOffset(1.2);
|
---|
672 | h[ type ].GetXaxis()->SetTitle(Form("Time slice (%.1f ns/slice)", 1./2.));
|
---|
673 |
|
---|
674 | // set Y axis paras
|
---|
675 | h[ type ].GetYaxis()->SetLabelSize(0.1);
|
---|
676 | h[ type ].GetYaxis()->SetTitleSize(0.1);
|
---|
677 | h[ type ].GetYaxis()->SetTitleOffset(0.3);
|
---|
678 | h[ type ].GetYaxis()->SetTitle("Amplitude (a.u.)");
|
---|
679 | }
|
---|
680 | CW = new TCanvas("CW","DRS Waveform",10,10,800,600);
|
---|
681 | CW->Divide(1, 3);
|
---|
682 | cFilter = new TCanvas("cFilter","filtered DRS Waveforms",10,10,800,600);
|
---|
683 | cFilter->Divide(1, 3);
|
---|
684 |
|
---|
685 | hStartCell = new TH2F("StartCell", "StartCell", 1440, 0., 1440., 1024, 0., 1024);
|
---|
686 |
|
---|
687 | }
|
---|
688 | */
|
---|
689 |
|
---|