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

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