source: trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.cc@ 7859

Last change on this file since 7859 was 7844, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 41.9 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Hendrik Bartko, 09/2004 <mailto:hbartko@mppmu.mpg.de>
19! Author(s): Markus Gaug, 05/2004 <mailto:markus@ifae.es>
20! Author(s): Diego Tescaro, 05/2004 <mailto:tescaro@pd.infn.it>
21!
22! Copyright: MAGIC Software Development, 2000-2004
23!
24!
25\* ======================================================================== */
26
27//////////////////////////////////////////////////////////////////////////////
28//
29// MExtractTimeAndChargeDigitalFilter
30//
31// Hendrik has promised to write more documentation
32//
33// The following variables have to be set by the derived class and
34// do not have defaults:
35// - fNumHiGainSamples
36// - fNumLoGainSamples
37// - fSqrtHiGainSamples
38// - fSqrtLoGainSamples
39//
40// The reading of automatic weights files (color, type) can be switched
41// off using EnableAutomaticWeights(kFALSE).
42//
43// An empty name or "-" as the weights file name is a synonym for
44// setting all weights to 1
45//
46// Input Containers:
47// MRawEvtData
48// MRawRunHeader
49// MPedestalCam
50// [MCalibrationPattern]
51//
52// Output Containers:
53// MArrivalTimeCam
54// MExtractedSignalCam
55//
56//////////////////////////////////////////////////////////////////////////////
57#include "MExtractTimeAndChargeDigitalFilter.h"
58
59#include <errno.h>
60#include <fstream>
61
62#include <TH1.h>
63#include <TH2.h>
64#include <TMatrix.h>
65
66#include "MLog.h"
67#include "MLogManip.h"
68
69#include "MParList.h"
70
71#include "MRawRunHeader.h"
72#include "MCalibrationPattern.h"
73#include "MExtractedSignalCam.h"
74
75#include "MPedestalPix.h"
76
77ClassImp(MExtractTimeAndChargeDigitalFilter);
78
79using namespace std;
80
81const Byte_t MExtractTimeAndChargeDigitalFilter::fgHiGainFirst = 0;
82const Byte_t MExtractTimeAndChargeDigitalFilter::fgHiGainLast = 15;
83const Byte_t MExtractTimeAndChargeDigitalFilter::fgLoGainFirst = 1;
84const Byte_t MExtractTimeAndChargeDigitalFilter::fgLoGainLast = 14;
85const Int_t MExtractTimeAndChargeDigitalFilter::fgWindowSizeHiGain = 4;
86const Int_t MExtractTimeAndChargeDigitalFilter::fgWindowSizeLoGain = 6;
87const Int_t MExtractTimeAndChargeDigitalFilter::fgBinningResolutionHiGain = 10;
88const Int_t MExtractTimeAndChargeDigitalFilter::fgBinningResolutionLoGain = 10;
89const Int_t MExtractTimeAndChargeDigitalFilter::fgSignalStartBinHiGain = 4;
90const Int_t MExtractTimeAndChargeDigitalFilter::fgSignalStartBinLoGain = 4;
91const Float_t MExtractTimeAndChargeDigitalFilter::fgOffsetLoGain = 0.95;
92const Float_t MExtractTimeAndChargeDigitalFilter::fgLoGainStartShift = -1.8;
93
94// --------------------------------------------------------------------------
95//
96// Default constructor.
97//
98// Calls:
99// - SetWindowSize();
100// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast)
101// - SetBinningResolution();
102//
103// Sets all weights to 1.
104//
105MExtractTimeAndChargeDigitalFilter::MExtractTimeAndChargeDigitalFilter(const char *name, const char *title)
106 : fTimeShiftHiGain(0.), fTimeShiftLoGain(0.), fAutomaticWeights(kTRUE), fRandomIter(0)
107{
108 fName = name ? name : "MExtractTimeAndChargeDigitalFilter";
109 fTitle = title ? title : "Digital Filter";
110
111 SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
112 SetWindowSize();
113 SetBinningResolution();
114 SetSignalStartBin();
115
116 SetOffsetLoGain(fgOffsetLoGain); // Order between both
117 SetLoGainStartShift(fgLoGainStartShift); // is important
118}
119
120// ---------------------------------------------------------------------------------------
121//
122// Checks:
123// - if a window is bigger than the one defined by the ranges, set it to the available range
124//
125// Sets:
126// - fNumHiGainSamples to: (Float_t)fWindowSizeHiGain
127// - fNumLoGainSamples to: (Float_t)fWindowSizeLoGain
128//
129void MExtractTimeAndChargeDigitalFilter::SetWindowSize(Int_t windowh, Int_t windowl)
130{
131
132 if (windowh != fgWindowSizeHiGain)
133 *fLog << warn << GetDescriptor()
134 << ": ATTENTION!!! If you are not Hendrik Bartko, do NOT use a different window size than the default." << endl;
135 if (windowl != fgWindowSizeLoGain)
136 *fLog << warn << GetDescriptor()
137 << ": ATTENTION!!! If you are not Hendrik Bartko, do NOT use a different window size than the default" << endl;
138
139 fWindowSizeHiGain = windowh;
140 fWindowSizeLoGain = windowl;
141
142 const Int_t availhirange = (Int_t)(fHiGainLast-fHiGainFirst+1);
143
144 if (fWindowSizeHiGain > availhirange)
145 {
146 *fLog << warn << GetDescriptor() << ": Hi Gain window size: " << Form("%2i",fWindowSizeHiGain);
147 *fLog << " is bigger than available range: [" << Form("%2i", (int)fHiGainFirst);
148 *fLog << "," << Form("%21", (int)fHiGainLast) << "]" << endl;
149
150 fHiGainLast = fHiGainFirst + fWindowSizeHiGain;
151
152 *fLog << warn << GetDescriptor() << ": Will set the upper range to: " << (int)fHiGainLast << endl;
153 }
154
155 if (fWindowSizeHiGain < 2)
156 {
157 fWindowSizeHiGain = 2;
158 *fLog << warn << GetDescriptor() << ": High Gain window size set to two samples" << endl;
159 }
160
161 if (fLoGainLast != 0 && fWindowSizeLoGain != 0)
162 {
163 const Int_t availlorange = (Int_t)(fLoGainLast-fLoGainFirst+1);
164
165 if (fWindowSizeLoGain > availlorange)
166 {
167 *fLog << warn << GetDescriptor() << ": Lo Gain window size: " << Form("%2i",fWindowSizeLoGain);
168 *fLog << " is bigger than available range: [" << Form("%2i", (int)fLoGainFirst);
169 *fLog << "," << Form("%21", (int)fLoGainLast) << "]" << endl;
170
171 fLoGainLast = fLoGainFirst + fWindowSizeLoGain;
172
173 *fLog << warn << GetDescriptor() << ": Will set the upper range to: " << (int)fLoGainLast << endl;
174 }
175
176 if (fWindowSizeLoGain<2)
177 {
178 fWindowSizeLoGain = 2;
179 *fLog << warn << GetDescriptor() << ": Low Gain window size set to two samples" << endl;
180 }
181 }
182 //
183 // We need here the effective number of samples which is about 2.5 in the case of a window
184 // size of 6. The exact numbers have to be found still.
185 //
186 fNumHiGainSamples = (Float_t)fWindowSizeHiGain/2.4;
187 fNumLoGainSamples = (Float_t)fWindowSizeLoGain/2.4;
188 fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
189 fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
190
191}
192
193// --------------------------------------------------------------------------
194//
195// Executing MExtractTimeAndCharge::PreProcess and searching for
196// MCalibrationPattern
197//
198Int_t MExtractTimeAndChargeDigitalFilter::PreProcess(MParList *pList)
199{
200 if (!MExtractTimeAndCharge::PreProcess(pList))
201 return kFALSE;
202
203 fCalibPattern = (MCalibrationPattern*)pList->FindObject("MCalibrationPattern");
204 return kTRUE;
205}
206
207// --------------------------------------------------------------------------
208//
209// The weights are determined using GetAutimaticWeights().
210//
211// kFALSE is returned if it returned an error.
212// kTRUE is returned if no new weights were set.
213//
214// If new weights are set
215// fTimeShiftHiGain
216// fTimeShiftLoGain
217// fNumHiGainSamples
218// fNumLoGainSamples
219// fSqrtHiGainSamples
220// fSqrtLoGainSamples
221// and
222// fSignals->SetUsedFADCSlices(...)
223// is updated accordingly.
224//
225Bool_t MExtractTimeAndChargeDigitalFilter::GetWeights()
226{
227 switch (GetAutomaticWeights())
228 {
229 case kERROR: // An error occured
230 return kFALSE;
231 case kFALSE: // No new weights set
232 return kTRUE;
233 }
234
235 //
236 // shift the times back to the right reference (start counting from 0)
237 //
238 // The high-gain has the extraction offset (fHiGainFirst) already included
239 // here for speed reason. The low-gain has a changing extraction offset,
240 // so it will be added at every event (in FindTimeAndChargeLoGain)
241 fTimeShiftHiGain = 0.5 + 1./fBinningResolutionHiGain + fHiGainFirst;
242 fTimeShiftLoGain = 0.5 + 1./fBinningResolutionLoGain;
243
244 //
245 // We need here the effective number of samples which is about 2.5 in the case of a window
246 // size of 6. The exact numbers have to be found still.
247 //
248 fNumHiGainSamples = (Float_t)fWindowSizeHiGain/2.4;
249 fNumLoGainSamples = (Float_t)fWindowSizeLoGain/2.4;
250 fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
251 fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
252
253 // From MExtractTimeAndCharge::ReInit
254 if (fSignals)
255 fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast+fHiLoLast, fNumHiGainSamples,
256 fLoGainFirst, fLoGainLast, fNumLoGainSamples);
257 return kTRUE;
258}
259
260// --------------------------------------------------------------------------
261//
262// InitArrays
263//
264// Gets called in the ReInit() and initialized the arrays
265//
266Bool_t MExtractTimeAndChargeDigitalFilter::InitArrays()
267{
268 if (!fRunHeader)
269 return kFALSE;
270
271 const Int_t rangehi = (Int_t)(fHiGainLast - fHiGainFirst + 1 + fHiLoLast);
272 const Int_t rangelo = (Int_t)(fLoGainLast - fLoGainFirst + 1);
273
274 fHiGainSignal.Set(rangehi);
275 fLoGainSignal.Set(rangelo);
276
277 return GetWeights();
278}
279
280// --------------------------------------------------------------------------
281//
282// Check if reading a new weights file is necessary because the calibration
283// pattern has changed. (Cannot be done in ReInit, because at this time
284// the calibration pattern is not available.
285// Then process the event.
286//
287Int_t MExtractTimeAndChargeDigitalFilter::Process()
288{
289 // Change Weights if the calibration patter changes
290 if (!GetWeights())
291 return kERROR;
292
293 // Process event
294 return MExtractTimeAndCharge::Process();
295}
296
297void MExtractTimeAndChargeDigitalFilter::CalcBinningResArrays()
298{
299
300 fArrBinningResHiGain.Set(fWindowSizeHiGain);
301 fArrBinningResHalfHiGain.Set(fWindowSizeHiGain);
302
303 for (int i=0; i<fWindowSizeHiGain; i++)
304 {
305 fArrBinningResHiGain[i] = fBinningResolutionHiGain*i;
306 fArrBinningResHalfHiGain[i] = fArrBinningResHiGain[i] + fBinningResolutionHalfHiGain;
307 }
308
309 fArrBinningResLoGain.Set(fWindowSizeLoGain);
310 fArrBinningResHalfLoGain.Set(fWindowSizeLoGain);
311
312 for (int i=0; i<fWindowSizeLoGain; i++)
313 {
314 fArrBinningResLoGain[i] = fBinningResolutionLoGain*i;
315 fArrBinningResHalfLoGain[i] = fArrBinningResLoGain[i] + fBinningResolutionHalfLoGain;
316 }
317}
318
319// --------------------------------------------------------------------------
320//
321// Apply the digital filter algorithm to the high-gain slices.
322//
323void MExtractTimeAndChargeDigitalFilter::FindTimeAndChargeHiGain(Byte_t *ptr, Byte_t *logain, Float_t &sum, Float_t &dsum,
324 Float_t &time, Float_t &dtime,
325 Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag)
326{
327
328 Int_t range = fHiGainLast - fHiGainFirst + 1;
329
330 const Byte_t *end = ptr + range;
331 Byte_t *p = ptr;
332
333 //
334 // Preparations for the pedestal subtraction (with AB-noise correction)
335 //
336 const Float_t pedes = ped.GetPedestal();
337 const Float_t ABoffs = ped.GetPedestalABoffset();
338
339 const Float_t pedmean[2] = { pedes + ABoffs, pedes - ABoffs };
340
341 range += fHiLoLast;
342 fMaxBinContent = 0;
343 //
344 // Check for saturation in all other slices
345 //
346 Int_t ids = fHiGainFirst;
347 Float_t *sample = fHiGainSignal.GetArray();
348 while (p<end)
349 {
350 *sample++ = (Float_t)*p - pedmean[(ids++ + abflag) & 0x1];
351
352 if (*p > fMaxBinContent)
353 {
354 Byte_t maxpos = p-ptr;
355
356 // range-fWindowSizeHiGain+1 == fHiLoLast isn't it?
357 if (maxpos > 1 && maxpos < (range - fWindowSizeHiGain + 1))
358 fMaxBinContent = *p;
359 }
360
361 if (*p++ >= fSaturationLimit)
362 if (!sat)
363 sat = ids-fHiGainFirst;
364 }
365
366 //
367 // Slide with a window of size fWindowSizeHiGain over the sample
368 // and multiply the entries with the corresponding weights
369 //
370 if (IsNoiseCalculation())
371 {
372 if (fRandomIter == fBinningResolutionHiGain)
373 fRandomIter = 0;
374 for (Int_t ids=0; ids < fWindowSizeHiGain; ids++)
375 {
376 const Int_t idx = fArrBinningResHiGain[ids] + fRandomIter;
377 sum += fAmpWeightsHiGain [idx]*fHiGainSignal[ids];
378 }
379 fRandomIter++;
380 return;
381 }
382
383 if (fHiLoLast != 0)
384 {
385
386 end = logain + fHiLoLast;
387
388 while (logain<end)
389 {
390
391 *sample++ = (Float_t)*logain - pedmean[(ids++ + abflag) & 0x1];
392
393 if (*logain++ >= fSaturationLimit)
394 if (!sat)
395 sat = ids-fHiGainFirst;
396 }
397 }
398
399 //
400 // allow no saturated slice
401 //
402 // if (sat > 0)
403 // return;
404
405 Float_t time_sum = 0.;
406 Float_t fmax = -9999.;
407 Float_t ftime_max = 0.;
408 Int_t max_p = 0;
409
410 //
411 // Calculate the sum of the first fWindowSize slices
412 //
413 for (Int_t i=0;i<range-fWindowSizeHiGain+1;i++)
414 {
415 sum = 0.;
416 time_sum = 0.;
417
418 //
419 // Slide with a window of size fWindowSizeHiGain over the sample
420 // and multiply the entries with the corresponding weights
421 //
422 for (Int_t sample=0; sample < fWindowSizeHiGain; sample++)
423 {
424 const Int_t idx = fBinningResolutionHiGain*sample+fBinningResolutionHalfHiGain;
425 const Float_t pex = fHiGainSignal[sample+i];
426 sum += fAmpWeightsHiGain [idx]*pex;
427 time_sum += fTimeWeightsHiGain[idx]*pex;
428 }
429
430 if (sum>fmax)
431 {
432 fmax = sum;
433 ftime_max = time_sum;
434 max_p = i;
435 }
436 } /* for (Int_t i=0;i<range-fWindowSizeHiGain+1;i++) */
437
438 if (fmax==0)
439 return;
440
441 ftime_max /= fmax;
442 Int_t t_iter = Int_t(ftime_max*fBinningResolutionHiGain);
443 Int_t sample_iter = 0;
444
445 while ( t_iter > fBinningResolutionHalfHiGain-1 || t_iter < -fBinningResolutionHalfHiGain )
446 {
447 if (t_iter > fBinningResolutionHalfHiGain-1)
448 {
449 t_iter -= fBinningResolutionHiGain;
450 max_p--;
451 sample_iter--;
452 }
453 if (t_iter < -fBinningResolutionHalfHiGain)
454 {
455 t_iter += fBinningResolutionHiGain;
456 max_p++;
457 sample_iter++;
458 }
459 }
460
461 sum = 0.;
462 time_sum = 0.;
463 if (max_p < 0)
464 max_p = 0;
465 if (max_p+fWindowSizeHiGain > range)
466 max_p = range-fWindowSizeHiGain;
467 //
468 // Slide with a window of size fWindowSizeHiGain over the sample
469 // and multiply the entries with the corresponding weights
470 //
471 for (Int_t sample=0; sample < fWindowSizeHiGain; sample++)
472 {
473 const Int_t idx = fArrBinningResHalfHiGain[sample] + t_iter;
474 const Int_t ids = max_p + sample;
475 const Float_t pex = fHiGainSignal[ids];
476 sum += fAmpWeightsHiGain [idx]*pex;
477 time_sum += fTimeWeightsHiGain[idx]*pex;
478 }
479
480 if (sum == 0)
481 return;
482
483 // here, the first high-gain slice is already included in the fTimeShiftHiGain
484// time = fTimeShiftHiGain + max_p - Float_t(t_iter)/fBinningResolutionHiGain;
485 time = max_p + fTimeShiftHiGain + (Float_t)fHiGainFirst /* this shifts the time to the start of the rising edge */
486 - ((Float_t)t_iter)/fBinningResolutionHiGain;
487
488 const Float_t timefineadjust = time_sum/sum;
489
490 if (TMath::Abs(timefineadjust) < 4./fBinningResolutionHiGain)
491 time -= timefineadjust;
492
493}
494
495// --------------------------------------------------------------------------
496//
497// Apply the digital filter algorithm to the low-gain slices.
498//
499void MExtractTimeAndChargeDigitalFilter::FindTimeAndChargeLoGain(Byte_t *ptr, Float_t &sum, Float_t &dsum,
500 Float_t &time, Float_t &dtime,
501 Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag)
502{
503
504 const Int_t range = fLoGainLast - fLoGainFirst + 1;
505
506 const Byte_t *end = ptr + range;
507 Byte_t *p = ptr;
508 //
509 // Prepare the low-gain pedestal
510 //
511 const Float_t pedes = ped.GetPedestal();
512 const Float_t ABoffs = ped.GetPedestalABoffset();
513
514 const Float_t pedmean[2] = { pedes + ABoffs, pedes - ABoffs };
515
516 //
517 // Check for saturation in all other slices
518 //
519 Float_t *sample = fLoGainSignal.GetArray();
520 Int_t ids = fLoGainFirst;
521 while (p<end)
522 {
523 *sample++ = (Float_t)*p - pedmean[(ids++ + abflag) & 0x1];
524
525 if (*p++ >= fSaturationLimit)
526 sat++;
527 }
528
529 //
530 // Slide with a window of size fWindowSizeLoGain over the sample
531 // and multiply the entries with the corresponding weights
532 //
533 if (IsNoiseCalculation())
534 {
535 if (fRandomIter == fBinningResolutionLoGain)
536 fRandomIter = 0;
537 for (Int_t ids=0; ids < fWindowSizeLoGain; ids++)
538 {
539 const Int_t idx = fArrBinningResLoGain[ids] + fRandomIter;
540 sum += fAmpWeightsLoGain [idx]*fLoGainSignal[ids];
541 }
542 return;
543 }
544
545 Float_t time_sum = 0.;
546 Float_t fmax = 0.;
547 Float_t ftime_max = 0.;
548 Int_t max_p = 0;
549
550 //
551 // Calculate the sum of the first fWindowSize slices
552 //
553 for (Int_t i=0;i<range-fWindowSizeLoGain+1;i++)
554 {
555 sum = 0.;
556 time_sum = 0.;
557
558 //
559 // Slide with a window of size fWindowSizeLoGain over the sample
560 // and multiply the entries with the corresponding weights
561 //
562 for (Int_t sample=0; sample < fWindowSizeLoGain; sample++)
563 {
564 const Int_t idx = fArrBinningResHalfLoGain[sample];
565 const Float_t pex = fLoGainSignal[sample+i];
566 sum += fAmpWeightsLoGain [idx]*pex;
567 time_sum += fTimeWeightsLoGain[idx]*pex;
568 }
569
570 if (sum>fmax)
571 {
572 fmax = sum;
573 ftime_max = time_sum;
574 max_p = i;
575 }
576 } /* for (Int_t i=0;i<range-fWindowSizeLoGain+1;i++) */
577
578 time = 0;
579 if (fmax==0)
580 return;
581
582 ftime_max /= fmax;
583 Int_t t_iter = Int_t(ftime_max*fBinningResolutionLoGain);
584 Int_t sample_iter = 0;
585
586 while ( t_iter > fBinningResolutionHalfLoGain-1 || t_iter < -fBinningResolutionHalfLoGain )
587 {
588 if (t_iter > fBinningResolutionHalfLoGain-1)
589 {
590 t_iter -= fBinningResolutionLoGain;
591 max_p--;
592 sample_iter--;
593 }
594 if (t_iter < -fBinningResolutionHalfLoGain)
595 {
596 t_iter += fBinningResolutionLoGain;
597 max_p++;
598 sample_iter++;
599 }
600 }
601
602 sum = 0.;
603 time_sum = 0.;
604
605 //
606 // Slide with a window of size fWindowSizeLoGain over the sample
607 // and multiply the entries with the corresponding weights
608 //
609 for (Int_t sample=0; sample < fWindowSizeLoGain; sample++)
610 {
611 const Int_t idx = fArrBinningResHalfLoGain[sample] + t_iter;
612 const Int_t ids = max_p + sample;
613 const Float_t pex = ids < 0 ? 0. : ( ids >= range ? 0. : fLoGainSignal[ids]);
614 sum += fAmpWeightsLoGain [idx]*pex;
615 time_sum += fTimeWeightsLoGain[idx]*pex;
616 }
617
618 if (sum == 0)
619 return;
620
621 time = max_p + fTimeShiftLoGain + (Float_t)fLoGainFirst /* this shifts the time to the start of the rising edge */
622 - ((Float_t)t_iter)/fBinningResolutionLoGain;
623
624 const Float_t timefineadjust = time_sum/sum;
625
626 if (TMath::Abs(timefineadjust) < 4./fBinningResolutionLoGain)
627 time -= timefineadjust;
628
629}
630
631// --------------------------------------------------------------------------
632//
633// Read the setup from a TEnv, eg:
634// MJPedestal.MExtractor.WindowSizeHiGain: 6
635// MJPedestal.MExtractor.WindowSizeLoGain: 6
636// MJPedestal.MExtractor.BinningResolutionHiGain: 10
637// MJPedestal.MExtractor.BinningResolutionLoGain: 10
638// MJPedestal.MExtractor.WeightsFile: filename
639// MJPedestal.MExtractor.AutomaticWeights: off
640//
641Int_t MExtractTimeAndChargeDigitalFilter::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
642{
643
644 Byte_t hw = fWindowSizeHiGain;
645 Byte_t lw = fWindowSizeLoGain;
646 Bool_t rc = kFALSE;
647
648 if (IsEnvDefined(env, prefix, "AutomaticWeights", print))
649 {
650 EnableAutomaticWeights(GetEnvValue(env, prefix, "AutomaticWeights", fAutomaticWeights));
651 rc = kTRUE;
652 }
653
654 if (IsEnvDefined(env, prefix, "WindowSizeHiGain", print))
655 {
656 hw = GetEnvValue(env, prefix, "WindowSizeHiGain", hw);
657 rc = kTRUE;
658 }
659 if (IsEnvDefined(env, prefix, "WindowSizeLoGain", print))
660 {
661 lw = GetEnvValue(env, prefix, "WindowSizeLoGain", lw);
662 rc = kTRUE;
663 }
664
665 if (rc)
666 SetWindowSize(hw, lw);
667
668 Bool_t rc2 = kFALSE;
669 Int_t brh = fBinningResolutionHiGain;
670 Int_t brl = fBinningResolutionLoGain;
671
672 if (IsEnvDefined(env, prefix, "BinningResolutionHiGain", print))
673 {
674 brh = GetEnvValue(env, prefix, brh);
675 rc2 = kTRUE;
676 }
677 if (IsEnvDefined(env, prefix, "BinningResolutionLoGain", print))
678 {
679 brl = GetEnvValue(env, prefix, brl);
680 rc2 = kTRUE;
681 }
682
683 if (rc2)
684 {
685 SetBinningResolution(brh, brl);
686 rc = kTRUE;
687 }
688
689 if (IsEnvDefined(env, prefix, "WeightsFile", print))
690 {
691 SetNameWeightsFile(GetEnvValue(env, prefix, "WeightsFile", ""));
692 rc = kTRUE;
693 }
694
695 return MExtractTimeAndCharge::ReadEnv(env, prefix, print) ? kTRUE : rc;
696}
697
698//----------------------------------------------------------------------------
699//
700// If automatic weights are requested, no default weights (name.IsNull())
701// are requested, fRunHeader is available and fRunHeader->IsMonteCarloRun()
702// is true prepend "MC_" in from of the name.
703//
704// return poth+name;
705//
706TString MExtractTimeAndChargeDigitalFilter::CompileWeightFileName(TString path, const TString &name) const
707{
708 if (fAutomaticWeights && !name.IsNull() && fRunHeader && fRunHeader->IsMonteCarloRun())
709 path += "MC_";
710
711 path += name;
712
713 return path;
714}
715
716//----------------------------------------------------------------------------
717//
718// Read a pre-defined weights file into the class.
719// This is mandatory for the extraction
720//
721// If filenname is empty, then all weights will be set to 1.
722//
723// Returns:
724// kTRUE: new weights set
725// kFALSE: no weights set
726// kERROR: error
727//
728Int_t MExtractTimeAndChargeDigitalFilter::ReadWeightsFile(TString filename, TString path)
729{
730 if (filename.IsNull())
731 {
732 fAmpWeightsHiGain .Set(fBinningResolutionHiGain*fWindowSizeHiGain);
733 fAmpWeightsLoGain .Set(fBinningResolutionLoGain*fWindowSizeLoGain);
734 fTimeWeightsHiGain.Set(fBinningResolutionHiGain*fWindowSizeHiGain);
735 fTimeWeightsLoGain.Set(fBinningResolutionLoGain*fWindowSizeLoGain);
736
737 fAmpWeightsHiGain.Reset(1);
738 fTimeWeightsHiGain.Reset(1);
739 fAmpWeightsLoGain.Reset(1);
740 fTimeWeightsLoGain.Reset(1);
741 return kTRUE;
742 }
743
744 // Add "MC_" in front of the filename if necessary
745 filename = CompileWeightFileName(path, filename);
746
747 //filename = MJob::ExpandPath(filename);
748
749 if (fNameWeightsFileSet==filename)
750 return kFALSE; // No file read
751
752 ifstream fin(filename.Data());
753 if (!fin)
754 {
755 *fLog << err << GetDescriptor() << ": ERROR - Cannot open file " << filename << ": ";
756 *fLog << strerror(errno) << endl;
757 return kERROR;
758 }
759
760 *fLog << all << GetDescriptor() << ": Reading weights in " << filename << "..." << flush;
761
762 Int_t len = 0;
763 Int_t cnt = 0;
764 Int_t line = 0;
765 Bool_t hi = kFALSE;
766 Bool_t lo = kFALSE;
767
768 TString str;
769
770 while (1)
771 {
772 str.ReadLine(fin);
773 if (!fin)
774 break;
775
776 line++;
777
778 if (str.Contains("# High Gain Weights:"))
779 {
780 if (hi)
781 {
782 *fLog << err << "ERROR - 'High Gain Weights' found twice in line #" << line << "." << endl;
783 return kERROR;
784 }
785
786 if (2!=sscanf(str.Data(), "# High Gain Weights:%2i %2i", &fWindowSizeHiGain, &fBinningResolutionHiGain))
787 {
788 *fLog << err << "ERROR - Wrong number of arguments in line #" << line << ":" << endl;
789 *fLog << str << endl;
790 return kERROR;
791 }
792
793 len = fBinningResolutionHiGain*fWindowSizeHiGain;
794 fAmpWeightsHiGain .Set(len);
795 fTimeWeightsHiGain.Set(len);
796 hi = kTRUE;
797 continue;
798 }
799
800 if (str.Contains("# Low Gain Weights:"))
801 {
802 if (lo)
803 {
804 *fLog << err << "ERROR - 'Lo Gain Weights' found twice in line #" << line << "." << endl;
805 return kERROR;
806 }
807
808 if (2!=sscanf(str.Data(),"# Low Gain Weights:%2i %2i", &fWindowSizeLoGain, &fBinningResolutionLoGain))
809 {
810 *fLog << err << "ERROR - Wrong number of arguments in line #" << line << ":" << endl;
811 *fLog << str << endl;
812 return kERROR;
813 }
814
815 len = fBinningResolutionLoGain*fWindowSizeLoGain;
816 fAmpWeightsLoGain .Set(len);
817 fTimeWeightsLoGain.Set(len);
818 lo = kTRUE;
819 continue;
820 }
821
822 // Handle lines with comments
823 if (str.Contains("#"))
824 continue;
825
826 // Nothing found so far
827 if (len == 0)
828 continue;
829
830 if (2!=sscanf(str.Data(), "%f %f",
831 lo ? &fAmpWeightsLoGain [cnt] : &fAmpWeightsHiGain [cnt],
832 lo ? &fTimeWeightsLoGain[cnt] : &fTimeWeightsHiGain[cnt]))
833 {
834 *fLog << err << "ERROR - Wrong number of arguments in line #" << line << ":" << endl;
835 *fLog << str << endl;
836 return kERROR;
837 }
838
839 if (++cnt == len)
840 {
841 len = 0;
842 cnt = 0;
843 }
844 }
845
846 if (cnt != len)
847 {
848 *fLog << err << "ERROR - Size mismatch in weights file " << filename << endl;
849 return kERROR;
850 }
851
852 if (!hi)
853 {
854 *fLog << err << "ERROR - No correct header found in weights file " << filename << endl;
855 return kERROR;
856 }
857
858 *fLog << "done." << endl;
859
860 *fLog << inf << " File contains " << fWindowSizeHiGain << " hi-gain slices ";
861 *fLog << "with a resolution of " << fBinningResolutionHiGain << endl;
862
863 *fLog << inf << " File contains " << fWindowSizeLoGain << " lo-gain slices ";
864 *fLog << "with a resolution of " << fBinningResolutionLoGain << endl;
865
866 CalcBinningResArrays();
867
868 switch (fWindowSizeHiGain)
869 {
870 case 4:
871 SetResolutionPerPheHiGain(0.036);
872 break;
873 case 6:
874 SetResolutionPerPheHiGain(0.021);
875 break;
876 default:
877 *fLog << warn << "Could not set the high-gain extractor resolution per phe for window size "
878 << fWindowSizeHiGain << endl;
879 }
880
881 switch (fWindowSizeLoGain)
882 {
883 case 4:
884 SetResolutionPerPheLoGain(0.005);
885 break;
886 case 6:
887 SetResolutionPerPheLoGain(0.004);
888 break;
889 default:
890 *fLog << warn << "Could not set the low-gain extractor resolution per phe for window size "
891 << fWindowSizeLoGain << endl;
892 }
893
894 fNameWeightsFileSet = filename;
895
896 return kTRUE;
897}
898
899
900//----------------------------------------------------------------------------
901//
902// The default (+ prepending possible "MC_") is read for:
903//
904// - RunType: Pedestal (independant of fAutomaticWeights)
905// - fAutomaticWeights disabled
906//
907// if fAutomaticWeights enabled:
908// - fNameWeightsFile.IsNull()
909// - !fCalibPattern
910// - fCalibPattern->GetPulserColor()==MCalibrationCam::kNONE
911//
912// If automatic weights are enabled, the case above didn't take place and
913// fNameWeightsFile starts with "calibration_weights_"
914// - the color (blue, UV) is replaced by the appropriate one
915// taken from the calibration pattern
916//
917// In most cases a debug output is printed. Further output about the color
918// determination can be switched on with debug level > 5;
919//
920// Returns:
921// kFALSE: No new weights set
922// kTRUE: New weights set
923// kERROR: Error
924//
925Int_t MExtractTimeAndChargeDigitalFilter::GetAutomaticWeights()
926{
927 const Ssiz_t pos = fNameWeightsFile.Last('/')+1;
928 const Ssiz_t len = fNameWeightsFile.Length();
929
930 // Split file name in path and name
931 TString path = fNameWeightsFile(0, pos>=0?pos:len);
932 TString name = fNameWeightsFile(pos>=0?pos:0, len);
933
934 // Remove trailing "MC_" for automatic weights
935 if (fAutomaticWeights && name.BeginsWith("MC_"))
936 name.Remove(0, 3);
937
938 // In case of a pedetsal run no calibration pattern can be available
939 // the default weights are always used.
940 if (fRunHeader->GetRunType()==MRawRunHeader::kRTPedestal)
941 {
942 *fLog << dbg << "Pedestal file... using default weights: " << fNameWeightsFile << endl;
943 return ReadWeightsFile(name, path);
944 }
945
946 // If automatic weights are switched off use default weights
947 if (!fAutomaticWeights)
948 {
949 *fLog << dbg << "Automatic weights switched off... using default weights: " << fNameWeightsFile << endl;
950 return ReadWeightsFile(name, path);
951 }
952
953 // If automatic weights are switched on but no filename is given raise error
954 if (fNameWeightsFile.IsNull())
955 {
956 *fLog << err << "ERROR - Cannot get automatic weights without default filename." << endl;
957 return kERROR;
958 }
959
960 // If this is no pedestal run, automatic weights are requested and a
961 // filename for the weights file is given pedestal-extraction from
962 // cosmics data is assumed.
963 if (!fCalibPattern)
964 {
965 *fLog << dbg << "No decoded calibration pattern available... using default weights: " << fNameWeightsFile << endl;
966 return ReadWeightsFile(name, path);
967 }
968
969 const Bool_t debug = gLog.GetDebugLevel()>5;
970
971 // If no calibration pattern is available do not change the
972 // current weighs or current weights file name.
973 if (fCalibPattern->GetPulserColor()==MCalibrationCam::kNONE)
974 {
975 // If we are extracting data and the calibration pattern is kNONE
976 // we assume that it is a data file without interleaved events
977 // and calibration pattern information available.
978 if ((fRunHeader->GetRunType()!=MRawRunHeader::kRTData && !fRunHeader->IsMonteCarloRun()) || debug)
979 *fLog << dbg << "No calibration color set so far... guessing default: " << fNameWeightsFile << endl;
980
981 return ReadWeightsFile(name, path);
982 }
983
984 if (debug)
985 {
986 *fLog << dbg << endl;
987 *fLog << underline << GetDescriptor() << endl;
988 *fLog << " Trying to get automatic weight for " << fNameWeightsFile << endl;
989 *fLog << " Run type: ";
990 }
991
992 if (name.BeginsWith("calibration_weights_") && fCalibPattern)
993 {
994 if (debug)
995 *fLog << " Calibration with color " << fCalibPattern->GetPulserColorStr() << ", setting ";
996 switch (fCalibPattern->GetPulserColor())
997 {
998 case MCalibrationCam::kBLUE: // 2
999 case MCalibrationCam::kGREEN: // 1
1000 if (debug)
1001 *fLog << "blue/green, ";
1002 name.ReplaceAll("UV", "blue");
1003 break;
1004
1005 case MCalibrationCam::kUV: // 3
1006 case MCalibrationCam::kCT1: // 0
1007 if (debug)
1008 *fLog << "UV/CT1, ";
1009 name.ReplaceAll("blue", "UV");
1010 break;
1011 case MCalibrationCam::kNONE:
1012 break;
1013 default: // kNone + etc
1014 *fLog << err << "ERROR - Cannot get automatic weights for " << fCalibPattern->GetPulserColorStr() << endl;
1015 return kERROR;
1016 }
1017 }
1018
1019 return ReadWeightsFile(name, path);
1020}
1021
1022//----------------------------------------------------------------------------
1023//
1024// Create the weights file
1025// Beware that the shape-histogram has to contain the pulse starting at bin 1
1026//
1027Bool_t MExtractTimeAndChargeDigitalFilter::WriteWeightsFile(TString filename, TH1F *shapehi, TH2F *autocorrhi,
1028 TH1F *shapelo, TH2F *autocorrlo )
1029{
1030
1031 const Int_t nbinshi = shapehi->GetNbinsX();
1032 Float_t binwidth = shapehi->GetBinWidth(1);
1033
1034 TH1F *derivativehi = new TH1F(Form("%s%s",shapehi->GetName(),"_der"),
1035 Form("%s%s",shapehi->GetTitle()," derivative"),
1036 nbinshi,
1037 shapehi->GetBinLowEdge(1),
1038 shapehi->GetBinLowEdge(nbinshi)+binwidth);
1039
1040 //
1041 // Calculate the derivative of shapehi
1042 //
1043 for (Int_t i = 1; i<nbinshi+1;i++)
1044 {
1045 derivativehi->SetBinContent(i,
1046 ((shapehi->GetBinContent(i+1)-shapehi->GetBinContent(i-1))/2./binwidth));
1047 derivativehi->SetBinError(i,
1048 (sqrt(shapehi->GetBinError(i+1)*shapehi->GetBinError(i+1)
1049 +shapehi->GetBinError(i-1)*shapehi->GetBinError(i-1))/2./binwidth));
1050 }
1051
1052 //
1053 // normalize the shapehi, such that the integral for fWindowSize slices is one!
1054 //
1055 Float_t sum = 0;
1056 Int_t lasttemp = fBinningResolutionHiGain * (fSignalStartBinHiGain + fWindowSizeHiGain);
1057 lasttemp = lasttemp > nbinshi ? nbinshi : lasttemp;
1058
1059 for (Int_t i=fBinningResolutionHiGain*fSignalStartBinHiGain; i<lasttemp; i++) {
1060 sum += shapehi->GetBinContent(i);
1061 }
1062 sum /= fBinningResolutionHiGain;
1063
1064 shapehi->Scale(1./sum);
1065 derivativehi->Scale(1./sum);
1066
1067 //
1068 // read in the noise auto-correlation function:
1069 //
1070 TMatrix Bhi(fWindowSizeHiGain,fWindowSizeHiGain);
1071
1072 for (Int_t i=0; i<fWindowSizeHiGain; i++){
1073 for (Int_t j=0; j<fWindowSizeHiGain; j++){
1074 Bhi[i][j]=autocorrhi->GetBinContent(i+1,j+1); //+fSignalStartBinHiGain +fSignalStartBinHiGain
1075 }
1076 }
1077 Bhi.Invert();
1078
1079 const Int_t nsizehi = fWindowSizeHiGain*fBinningResolutionHiGain;
1080 fAmpWeightsHiGain.Set(nsizehi);
1081 fTimeWeightsHiGain.Set(nsizehi);
1082
1083 //
1084 // Loop over relative time in one BinningResolution interval
1085 //
1086 Int_t start = fBinningResolutionHiGain*(fSignalStartBinHiGain + 1);
1087
1088 for (Int_t i = -fBinningResolutionHalfHiGain+1; i<=fBinningResolutionHalfHiGain; i++)
1089 {
1090
1091 TMatrix g(fWindowSizeHiGain,1);
1092 TMatrix gT(1,fWindowSizeHiGain);
1093 TMatrix d(fWindowSizeHiGain,1);
1094 TMatrix dT(1,fWindowSizeHiGain);
1095
1096 for (Int_t count=0; count < fWindowSizeHiGain; count++){
1097
1098 g[count][0]=shapehi->GetBinContent(start
1099 +fBinningResolutionHiGain*count+i);
1100 gT[0][count]=shapehi->GetBinContent(start
1101 +fBinningResolutionHiGain*count+i);
1102 d[count][0]=derivativehi->GetBinContent(start
1103 +fBinningResolutionHiGain*count+i);
1104 dT[0][count]=derivativehi->GetBinContent(start
1105 +fBinningResolutionHiGain*count+i);
1106 }
1107
1108 TMatrix m_denom = (gT*(Bhi*g))*(dT*(Bhi*d)) - (dT*(Bhi*g))*(dT*(Bhi*g));
1109 Float_t denom = m_denom[0][0]; // ROOT thinks, m_denom is still a matrix
1110
1111 TMatrix m_first = dT*(Bhi*d); // ROOT thinks, m_first is still a matrix
1112 Float_t first = m_first[0][0]/denom;
1113
1114 TMatrix m_last = gT*(Bhi*d); // ROOT thinks, m_last is still a matrix
1115 Float_t last = m_last[0][0]/denom;
1116
1117 TMatrix m1 = gT*Bhi;
1118 m1 *= first;
1119
1120 TMatrix m2 = dT*Bhi;
1121 m2 *=last;
1122
1123 TMatrix w_amp = m1 - m2;
1124
1125 TMatrix m_first1 = gT*(Bhi*g);
1126 Float_t first1 = m_first1[0][0]/denom;
1127
1128 TMatrix m_last1 = gT*(Bhi*d);
1129 Float_t last1 = m_last1 [0][0]/denom;
1130
1131 TMatrix m11 = dT*Bhi;
1132 m11 *=first1;
1133
1134 TMatrix m21 = gT*Bhi;
1135 m21 *=last1;
1136
1137 TMatrix w_time= m11 - m21;
1138
1139 for (Int_t count=0; count < fWindowSizeHiGain; count++)
1140 {
1141 const Int_t idx = i+fBinningResolutionHalfHiGain+fBinningResolutionHiGain*count-1;
1142 fAmpWeightsHiGain [idx] = w_amp [0][count];
1143 fTimeWeightsHiGain[idx] = w_time[0][count];
1144 }
1145
1146 } // end loop over i
1147
1148 //
1149 // Low Gain histograms
1150 //
1151 TH1F *derivativelo = NULL;
1152 if (shapelo)
1153 {
1154 const Int_t nbinslo = shapelo->GetNbinsX();
1155 binwidth = shapelo->GetBinWidth(1);
1156
1157 derivativelo = new TH1F(Form("%s%s",shapelo->GetName(),"_der"),
1158 Form("%s%s",shapelo->GetTitle()," derivative"),
1159 nbinslo,
1160 shapelo->GetBinLowEdge(1),
1161 shapelo->GetBinLowEdge(nbinslo)+binwidth);
1162
1163 //
1164 // Calculate the derivative of shapelo
1165 //
1166 for (Int_t i = 1; i<nbinslo+1;i++)
1167 {
1168 derivativelo->SetBinContent(i,
1169 ((shapelo->GetBinContent(i+1)-shapelo->GetBinContent(i-1))/2./binwidth));
1170 derivativelo->SetBinError(i,
1171 (sqrt(shapelo->GetBinError(i+1)*shapelo->GetBinError(i+1)
1172 +shapelo->GetBinError(i-1)*shapelo->GetBinError(i-1))/2./binwidth));
1173 }
1174
1175 //
1176 // normalize the shapelo, such that the integral for fWindowSize slices is one!
1177 //
1178 sum = 0;
1179 lasttemp = fBinningResolutionLoGain * (fSignalStartBinLoGain + fWindowSizeLoGain);
1180 lasttemp = lasttemp > nbinslo ? nbinslo : lasttemp;
1181
1182 for (Int_t i=fBinningResolutionLoGain*fSignalStartBinLoGain; i<lasttemp; i++)
1183 sum += shapelo->GetBinContent(i);
1184
1185 sum /= fBinningResolutionLoGain;
1186
1187 shapelo->Scale(1./sum);
1188 derivativelo->Scale(1./sum);
1189
1190 //
1191 // read in the noise auto-correlation function:
1192 //
1193 TMatrix Blo(fWindowSizeLoGain,fWindowSizeLoGain);
1194
1195 for (Int_t i=0; i<fWindowSizeLoGain; i++){
1196 for (Int_t j=0; j<fWindowSizeLoGain; j++){
1197 Blo[i][j]=autocorrlo->GetBinContent(i+1+fSignalStartBinLoGain,j+1+fSignalStartBinLoGain);
1198 }
1199 }
1200 Blo.Invert();
1201
1202 const Int_t nsizelo = fWindowSizeLoGain*fBinningResolutionLoGain;
1203 fAmpWeightsLoGain.Set(nsizelo);
1204 fTimeWeightsLoGain.Set(nsizelo);
1205
1206 //
1207 // Loop over relative time in one BinningResolution interval
1208 //
1209 Int_t start = fBinningResolutionLoGain*fSignalStartBinLoGain + fBinningResolutionHalfLoGain;
1210
1211 for (Int_t i = -fBinningResolutionHalfLoGain+1; i<=fBinningResolutionHalfLoGain; i++)
1212 {
1213
1214 TMatrix g(fWindowSizeLoGain,1);
1215 TMatrix gT(1,fWindowSizeLoGain);
1216 TMatrix d(fWindowSizeLoGain,1);
1217 TMatrix dT(1,fWindowSizeLoGain);
1218
1219 for (Int_t count=0; count < fWindowSizeLoGain; count++){
1220
1221 g[count][0] = shapelo->GetBinContent(start
1222 +fBinningResolutionLoGain*count+i);
1223 gT[0][count]= shapelo->GetBinContent(start
1224 +fBinningResolutionLoGain*count+i);
1225 d[count][0] = derivativelo->GetBinContent(start
1226 +fBinningResolutionLoGain*count+i);
1227 dT[0][count]= derivativelo->GetBinContent(start
1228 +fBinningResolutionLoGain*count+i);
1229 }
1230
1231 TMatrix m_denom = (gT*(Blo*g))*(dT*(Blo*d)) - (dT*(Blo*g))*(dT*(Blo*g));
1232 Float_t denom = m_denom[0][0]; // ROOT thinks, m_denom is still a matrix
1233
1234 TMatrix m_first = dT*(Blo*d); // ROOT thinks, m_first is still a matrix
1235 Float_t first = m_first[0][0]/denom;
1236
1237 TMatrix m_last = gT*(Blo*d); // ROOT thinks, m_last is still a matrix
1238 Float_t last = m_last[0][0]/denom;
1239
1240 TMatrix m1 = gT*Blo;
1241 m1 *= first;
1242
1243 TMatrix m2 = dT*Blo;
1244 m2 *=last;
1245
1246 TMatrix w_amp = m1 - m2;
1247
1248 TMatrix m_first1 = gT*(Blo*g);
1249 Float_t first1 = m_first1[0][0]/denom;
1250
1251 TMatrix m_last1 = gT*(Blo*d);
1252 Float_t last1 = m_last1 [0][0]/denom;
1253
1254 TMatrix m11 = dT*Blo;
1255 m11 *=first1;
1256
1257 TMatrix m21 = gT*Blo;
1258 m21 *=last1;
1259
1260 TMatrix w_time= m11 - m21;
1261
1262 for (Int_t count=0; count < fWindowSizeLoGain; count++)
1263 {
1264 const Int_t idx = i+fBinningResolutionHalfLoGain+fBinningResolutionLoGain*count-1;
1265 fAmpWeightsLoGain [idx] = w_amp [0][count];
1266 fTimeWeightsLoGain[idx] = w_time[0][count];
1267 }
1268
1269 } // end loop over i
1270 }
1271
1272 ofstream fn(filename.Data());
1273
1274 fn << "# High Gain Weights: " << fWindowSizeHiGain << " " << fBinningResolutionHiGain << endl;
1275 fn << "# (Amplitude) (Time) " << endl;
1276
1277 for (Int_t i=0; i<nsizehi; i++)
1278 fn << "\t" << fAmpWeightsHiGain[i] << "\t" << fTimeWeightsHiGain[i] << endl;
1279
1280 fn << "# Low Gain Weights: " << fWindowSizeLoGain << " " << fBinningResolutionLoGain << endl;
1281 fn << "# (Amplitude) (Time) " << endl;
1282
1283 for (Int_t i=0; i<nsizehi; i++)
1284 fn << "\t" << fAmpWeightsLoGain[i] << "\t" << fTimeWeightsLoGain[i] << endl;
1285
1286 delete derivativehi;
1287 if (derivativelo)
1288 delete derivativelo;
1289
1290 return kTRUE;
1291}
1292
1293void MExtractTimeAndChargeDigitalFilter::Print(Option_t *o) const
1294{
1295 if (IsA()==Class())
1296 *fLog << GetDescriptor() << ":" << endl;
1297
1298 MExtractTimeAndCharge::Print(o);
1299 *fLog << " Time Shift HiGain: " << fTimeShiftHiGain << " LoGain: " << fTimeShiftLoGain << endl;
1300 *fLog << " Window Size HiGain: " << fWindowSizeHiGain << " LoGain: " << fWindowSizeLoGain << endl;
1301 *fLog << " Binning Res HiGain: " << fBinningResolutionHiGain << " LoGain: " << fBinningResolutionHiGain << endl;
1302 *fLog << " Weights File desired: " << (fNameWeightsFile.IsNull()?"-":fNameWeightsFile.Data()) << endl;
1303 if (!fNameWeightsFileSet.IsNull())
1304 *fLog << " Weights File set: " << fNameWeightsFileSet << endl;
1305
1306 TString opt(o);
1307 if (!opt.Contains("weights"))
1308 return;
1309
1310 *fLog << endl;
1311 *fLog << inf << "Using the following weights: " << endl;
1312 *fLog << "Hi-Gain:" << endl;
1313 for (Int_t i=0; i<fBinningResolutionHiGain*fWindowSizeHiGain; i++)
1314 *fLog << " " << fAmpWeightsHiGain[i] << " \t " << fTimeWeightsHiGain[i] << endl;
1315
1316 *fLog << "Lo-Gain:" << endl;
1317 for (Int_t i=0; i<fBinningResolutionLoGain*fWindowSizeLoGain; i++)
1318 *fLog << " " << fAmpWeightsLoGain[i] << " \t " << fTimeWeightsLoGain[i] << endl;
1319}
Note: See TracBrowser for help on using the repository browser.